线性回归实例
junjun
2016年2月7日
实例一、一元线性回归
测试沸点与气压的关系:Forbes数据共有四列，分别表示 沸点（F）、气压（英寸汞柱）、log（气压）、100*log（气压）
#1、读取数据
X <- matrix(c(
194.5, 20.79, 1.3179, 131.79,
194.3, 20.79, 1.3179, 131.79,
197.9, 22.40, 1.3502, 135.02,
198.4, 22.67, 1.3555, 135.55,
199.4, 23.15, 1.3646, 136.46,
199.9, 23.35, 1.3683, 136.83,
200.9, 23.89, 1.3782, 137.82,
201.1, 23.99, 1.3800, 138.00,
201.4, 24.02, 1.3806, 138.06,
201.3, 24.01, 1.3805, 138.05,
203.6, 25.14, 1.4004, 140.04,
204.6, 26.57, 1.4244, 142.44,
209.5, 28.49, 1.4547, 145.47,
208.6, 27.76, 1.4434, 144.34,
210.7, 29.04, 1.4630, 146.30,
211.9, 29.88, 1.4754, 147.54,
212.2, 30.06, 1.4780, 147.80),
ncol=4, byrow=T,
dimnames = list(1:17, c("F", "h", "log", "log100")))
#2、数据预处理
#1）将数据转化为数据框格式
forbes <- as.data.frame(X)
#2）画图分析因变量与自变量之间的关系
plot(forbes$F, forbes$log100)#3、建模：从散点图上可以看出，这些点基本上落在一条一直上，所以做回归分析
model_lm <- lm(log100~F, data=forbes)
#4、模型评估
#1)分析summary函数结果
summary(model_lm)##
## Call:
## lm(formula = log100 ~ F, data = forbes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32261 -0.14530 -0.06750 0.02111 1.35924
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -42.13087 3.33895 -12.62 2.17e-09 ***
## F 0.89546 0.01645 54.45 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3789 on 15 degrees of freedom
## Multiple R-squared: 0.995, Adjusted R-squared: 0.9946
## F-statistic: 2965 on 1 and 15 DF, p-value: < 2.2e-16#从上面可以看出：对应于两个系数的P值均<2.17*10-9，是非常显著的；关于方程的检验，残差的标准差为0.3789，相关系数（R2）为0.995，关于F分布的P值<2.2*10-16也是分成显著的。所以，该模型能通过t检验和F检验。因此，回归方程为 Y-42.13087+0.89546X
#2）分析残差：函数residuals函数计算回归方程的残差。计算残差并画出关于残差的散点图
y.residuals <- residuals(model_lm)
plot(y.residuals)#从图中可以看出第12个样本可能会有问题，比其他的样本点的残差大得多。因此，这个点可能不正确，或者模型的差假设不正确，或者是方差不是常数等。需要对这个问题进行分析。
#5、模型优化
#1）把每个点都标记出来
a <- as.data.frame(y.residuals)
a$id <- 1:nrow(a)
names(a) <- list("value", "id")
b <- a[, c("id", "value")]
plot(value~id, b, xlab=b$id)
for(i in 1:nrow(b))
text(b[i, 1], b[i, 2], labels = i, adj=1.2)#2）去掉第12号样本点
i <- 1:17
forbes <- as.data.frame(X[i!=12, ])
model_lm12 <- lm(log100~F, data=forbes)
summary(model_lm12)##
## Call:
## lm(formula = log100 ~ F, data = forbes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21175 -0.06194 0.01590 0.09077 0.13042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -41.30180 1.00038 -41.29 5.01e-16 ***
## F 0.89096 0.00493 180.73 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1133 on 14 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9995
## F-statistic: 3.266e+04 on 1 and 14 DF, p-value: < 2.2e-16#去掉第12号样本后，回归方程的系数没有太大的变化，但系数的标准差和残差的标准差有很大的变化，小了3倍左右，相关系数R2也有提高。
#6、预测
x <- data.frame(F=195)
pred <- predict(model_lm12, x)
pred## 1
## 132.4347
实例二、多元线性回归
试根据提供的数据建立一个数学模型，分析牙膏销售量与其他因素的关系，为定制加个人策略提供数量依据。 牙膏销量为Y，价格差为X1，公司的广告费为X2，假设基本模型为线性模型 Y=b0+b1X1+b2X2+e 输入数据，调用R软件中的lm（）函数求解，并用summary（）显示计算结果
#1、加载数据
toothpaste<-data.frame(
X1=c(-0.05, 0.25,0.60,0, 0.25,0.20, 0.15,0.05,-0.15, 0.15, 0.20, 0.10,0.40,0.45,0.35,0.30, 0.50,0.50, 0.40,-0.05, -0.05,-0.10,0.20,0.10,0.50,0.60,-0.05,0, 0.05, 0.55),
X2=c( 5.50,6.75,7.25,5.50,7.00,6.50,6.75,5.25,5.25,6.00, 6.50,6.25,7.00,6.90,6.80,6.80,7.10,7.00,6.80,6.50, 6.25,6.00,6.50,7.00,6.80,6.80,6.50,5.75,5.80,6.80),
Y =c( 7.38,8.51,9.52,7.50,9.33,8.28,8.75,7.87,7.10,8.00, 7.89,8.15,9.10,8.86,8.90,8.87,9.26,9.00,8.75,7.95, 7.65,7.27,8.00,8.50,8.75,9.21,8.27,7.67,7.93,9.26)
)
#2、建立模型
model_lm <- lm(Y~X1+X2, data=toothpaste)
#3、模型评估
#1）查看summary函数的结果
summary(model_lm)##
## Call:
## lm(formula = Y ~ X1 + X2, data = toothpaste)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49779 -0.12031 -0.00867 0.11084 0.58106
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4075 0.7223 6.102 1.62e-06 ***
## X1 1.5883 0.2994 5.304 1.35e-05 ***
## X2 0.5635 0.1191 4.733 6.25e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2383 on 27 degrees of freedom
## Multiple R-squared: 0.886, Adjusted R-squared: 0.8776
## F-statistic: 105 on 2 and 27 DF, p-value: 1.845e-13#由此得到销售量与价格差与广告费之间的关系为 Y=4.4075+1.5883X1+0.5635X2
#2）分别做y与x1和y与x2的散点图，可以看出y与x1的散点图拟合为一条直线，y与x2的散点图拟合为一条曲线。所以，修改模型为：
attach(toothpaste)
plot(Y~X1)plot(Y~X2)model_lm2 <- update(model_lm, .~.+I(X2^2))
summary(model_lm2)##
## Call:
## lm(formula = Y ~ X1 + X2 + I(X2^2), data = toothpaste)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40330 -0.14509 -0.03035 0.15488 0.46602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.3244 5.6415 3.071 0.00495 **
## X1 1.3070 0.3036 4.305 0.00021 ***
## X2 -3.6956 1.8503 -1.997 0.05635 .
## I(X2^2) 0.3486 0.1512 2.306 0.02934 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2213 on 26 degrees of freedom
## Multiple R-squared: 0.9054, Adjusted R-squared: 0.8945
## F-statistic: 82.94 on 3 and 26 DF, p-value: 1.944e-13#此时，模型残差的标准差有所下降，相关系数的平方有所上升，这说明修正是合理的。但是，对应于b2的P值>0.05。进一步分析，作b的区间估计
#3)编写函数：求拟合函数fm的区间估计（a=0.05）
beta.int <- function(fm, alpha=0.05) {
A <- summary(fm)$coefficients
df <- fm$df.residual
left <- A[,1]-A[,2]*qt(1-alpha/2, df)
right <- A[,1]+A[,2]*qt(1-alpha/2, df)
rowname <- dimnames(A)[[1]]
colname <- c("Estimate", "Left", "Right")
matrix(c(A[,1], left, right), ncol=3, dimnames=list(rowname, colname))
}
beta.int(model_lm2)## Estimate Left Right
## (Intercept) 17.3243685 5.72818421 28.9205529
## X1 1.3069887 0.68290927 1.9310682
## X2 -3.6955867 -7.49886317 0.1076898
## I(X2^2) 0.3486117 0.03786354 0.6593598#由此可知，b2的区间估计是[-7.49886317, 0.1076898]，包含了0，也就是说b2的值可能会使0。所以，去掉X2的一次项，再进行分析
model_lm3 <- update(model_lm2, .~.-X2)
summary(model_lm3)##
## Call:
## lm(formula = Y ~ X1 + I(X2^2), data = toothpaste)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4859 -0.1141 -0.0046 0.1053 0.5592
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.07667 0.35531 17.102 5.17e-16 ***
## X1 1.52498 0.29859 5.107 2.28e-05 ***
## I(X2^2) 0.04720 0.00952 4.958 3.41e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2332 on 27 degrees of freedom
## Multiple R-squared: 0.8909, Adjusted R-squared: 0.8828
## F-statistic: 110.2 on 2 and 27 DF, p-value: 1.028e-13#由此可知，此模型虽然通过了F检验和T检验，但是与上一模型对比来看，残差的标准差上升了，相关系数下降了，所以模型有不足之处，需要进一步修正，考虑x1与x2的交互作用
#4)添加X1与X2的相互作用
model_lm4 <- update(model_lm3, .~.+X1*X2)
summary(model_lm4)##
## Call:
## lm(formula = Y ~ X1 + I(X2^2) + X2 + X1:X2, data = toothpaste)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43725 -0.11754 0.00489 0.12263 0.38410
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.1133 7.4832 3.890 0.000656 ***
## X1 11.1342 4.4459 2.504 0.019153 *
## I(X2^2) 0.6712 0.2027 3.312 0.002824 **
## X2 -7.6080 2.4691 -3.081 0.004963 **
## X1:X2 -1.4777 0.6672 -2.215 0.036105 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2063 on 25 degrees of freedom
## Multiple R-squared: 0.9209, Adjusted R-squared: 0.9083
## F-statistic: 72.78 on 4 and 25 DF, p-value: 2.107e-13#由上可知：模型通过t检验和F检验，并且残差的标准差减少了，相关系数增加了。因此，模型最终为Y=29.1133+11.1342X2-7.6080X2+0.6712X2^2-1.4777X1X2+e
实例三、逐步回归
某种水泥在凝固时放出的热量Y（卡/克）与水泥中四种化学成分X1、X2、X3、X4有关，现测得13组数据，选出主要的变量，建立Y关于他们的线性回归方程
#1、加载数据
cement<-data.frame(
X1=c( 7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10),
X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68),
X3=c( 6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8),
X4=c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12),
Y =c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1,115.9, 83.8, 113.3, 109.4)
)
#2、建模
model_lm <- lm(Y~., data=cement)
#3、模型评估
summary(model_lm)##
## Call:
## lm(formula = Y ~ ., data = cement)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1750 -1.6709 0.2508 1.3783 3.9254
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.4054 70.0710 0.891 0.3991
## X1 1.5511 0.7448 2.083 0.0708 .
## X2 0.5102 0.7238 0.705 0.5009
## X3 0.1019 0.7547 0.135 0.8959
## X4 -0.1441 0.7091 -0.203 0.8441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.446 on 8 degrees of freedom
## Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736
## F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07#从上可以看到，如果选择全部变量作回归方程，效果是不好的。因为，回归方程的系数全部都没有通过检验。
#4、用函数step（）作逐步回归
model_step <- step(model_lm)## Start: AIC=26.94
## Y ~ X1 + X2 + X3 + X4
##
## Df Sum of Sq RSS AIC
## - X3 1 0.1091 47.973 24.974
## - X4 1 0.2470 48.111 25.011
## - X2 1 2.9725 50.836 25.728
## <none> 47.864 26.944
## - X1 1 25.9509 73.815 30.576
##
## Step: AIC=24.97
## Y ~ X1 + X2 + X4
##
## Df Sum of Sq RSS AIC
## <none> 47.97 24.974
## - X4 1 9.93 57.90 25.420
## - X2 1 26.79 74.76 28.742
## - X1 1 820.91 868.88 60.629##这里去掉x3后，再去掉任意变量，都会使AIC值增大。
#分析计算结果
model_update <- update(model_lm, .~.-X3)
summary(model_update)##
## Call:
## lm(formula = Y ~ X1 + X2 + X4, data = cement)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0919 -1.8016 0.2562 1.2818 3.8982
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 71.6483 14.1424 5.066 0.000675 ***
## X1 1.4519 0.1170 12.410 5.78e-07 ***
## X2 0.4161 0.1856 2.242 0.051687 .
## X4 -0.2365 0.1733 -1.365 0.205395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.309 on 9 degrees of freedom
## Multiple R-squared: 0.9823, Adjusted R-squared: 0.9764
## F-statistic: 166.8 on 3 and 9 DF, p-value: 3.323e-08#由上可知：回归系数检验的显著性水平有很大提高，但变量X2、X4系数检验的显著水平仍不理想。在R软件中，还有add1()和drop1()方法，尝试去掉一个变量。
#5、下面，用drop1()计算：
drop1(model_step)## Single term deletions
##
## Model:
## Y ~ X1 + X2 + X4
## Df Sum of Sq RSS AIC
## <none> 47.97 24.974
## X1 1 820.91 868.88 60.629
## X2 1 26.79 74.76 28.742
## X4 1 9.93 57.90 25.420#去掉变量x4残差平方和增大9.93，AIC增长也是最小的，所以去掉x4
#6、建立最终的模型
model <- lm(Y~X1+X2, data=cement)
summary(model)##
## Call:
## lm(formula = Y ~ X1 + X2, data = cement)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.893 -1.574 -1.302 1.363 4.048
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.57735 2.28617 23.00 5.46e-10 ***
## X1 1.46831 0.12130 12.11 2.69e-07 ***
## X2 0.66225 0.04585 14.44 5.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.406 on 10 degrees of freedom
## Multiple R-squared: 0.9787, Adjusted R-squared: 0.9744
## F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09#从上可知：所有的检验均是显著的，最后得到“最优”的回归方程为 Y=52.57735+1.46831X1+0.66225X2
实例四、回归诊断
在建立模型前后，需要做回归诊断：异常样品的存在往往会给回归模型带来不稳定，所以提出了回归诊断的问题，主要内容有
（1）关于误差项是否满足：①独立性；②等方差性；③正态性 （2）选择线性模型是否合适？ （3）是否存在异常样本？ （4）回归分析的结果是否对某些样本的依赖过重？即回归模型是否具备稳定性？ （5）自变量之间是否存在高度相关？即是否有多重共线性问题存在？
#1、残差：判断模型的参数是否符合正态性分布
#residuals()或resid()函数提供了模型残差的计算，得到残差后，可以对残差进行检验，如正态性检验等。shapiro.test()函数为正态性检验函数。
#实例：对Forbes数据得到回归模型，得到的残差做W正态性检验
X <- matrix(c(
194.5, 20.79, 1.3179, 131.79,
194.3, 20.79, 1.3179, 131.79,
197.9, 22.40, 1.3502, 135.02,
198.4, 22.67, 1.3555, 135.55,
199.4, 23.15, 1.3646, 136.46,
199.9, 23.35, 1.3683, 136.83,
200.9, 23.89, 1.3782, 137.82,
201.1, 23.99, 1.3800, 138.00,
201.4, 24.02, 1.3806, 138.06,
201.3, 24.01, 1.3805, 138.05,
203.6, 25.14, 1.4004, 140.04,
204.6, 26.57, 1.4244, 142.44,
209.5, 28.49, 1.4547, 145.47,
208.6, 27.76, 1.4434, 144.34,
210.7, 29.04, 1.4630, 146.30,
211.9, 29.88, 1.4754, 147.54,
212.2, 30.06, 1.4780, 147.80),
ncol=4, byrow=T,
dimnames = list(1:17, c("F", "h", "log", "log100")))
forbes <- as.data.frame(X)
model_lm <- lm(log100~F, data=forbes)
y.residuals <- residuals(model_lm)
shapiro.test(y.residuals)##
## Shapiro-Wilk normality test
##
## data: y.residuals
## W = 0.54654, p-value = 3.302e-06#此时，P值小于0.05拒绝原假设，即该数据不符合正态分布。因此，在这个例子中我们也可发现较大的W统计量的情况下，我们依然可能拒绝其符合正态分布。
i <- 1:17
forbes <- as.data.frame(X[i!=12, ])
model_lm2 <- lm(log100~F, data=forbes)
y.residuals <- residuals(model_lm2)
shapiro.test(y.residuals)##
## Shapiro-Wilk normality test
##
## data: y.residuals
## W = 0.92215, p-value = 0.1827#此时，其W统计量接近1，p值显著大于0.05，所以我们没办法拒绝其符合正态分布.所以，通过正态性检验，去掉 第12号样本点是合理的。
#2、残差图：rstandard()计算回归模型的标准化残差
#利用上面的血压与身高的数据：
blood<-data.frame(
X1=c(76.0, 91.5, 85.5, 82.5, 79.0, 80.5, 74.5, 79.0, 85.0, 76.5, 82.0, 95.0, 92.5),
X2=c(50, 20, 20, 30, 30, 50, 60, 50, 40, 55, 40, 40, 20),
Y= c(120, 141, 124, 126, 117, 125, 123, 125, 132, 123, 132, 155, 147)
)
#画残差图：
model_lm <- lm(Y~X1+X2, data=blood)
y.residuals <- resid(model_lm)
y.fit <- predict(model_lm)
plot(y.residuals~y.fit)#画标准化残差图：
dev.new()
y.rst <- rstandard(model_lm)
plot(y.rst~y.fit)
#从上可知：当残差服从正态分布的假设成立时，标准化残差应近似的服从标准正态分布。对于标准化残差，应该有95%的样本点落在区间[-2, 2]中，通过标准化残差图，更容易诊断出回归模型是否出现问题。
#3、残差的Q-Q图：正态分布的检验方法——QQ图，可以检验残差的正态性。
a <- c(1:100)
qqnorm(a)
qqline(a)
#qqplot（）函数提供了更为精确的正态假设检验方法，画出了n-p-1个自由度的t分布下的学生化残差图形，n为样本大小，p是回归参数的数目（包括截距项） 用法：qqplot(lm(),data=))
#注意：car包中的qqPlot(model, labels, id.method, simulate)函数：分位比较图。labels为文本字符串向量；id.method为点识别方法；simulate为T,计算置信区间。用于检测正态性。
#4、多重共线性：一个变量可以由其他变量求出，例如，学生的总成绩可以由各科成绩求出。度量多重共线性严重程度的一个重要指标是矩阵的条件数，可以由函数kappa()求出。在R中，函数kappa()计算矩阵的条件数，其使用方法为：kappa(z, exact=F, …) 其中，z是矩阵（相关系数的矩阵,可由cor函数求得），exact是逻辑变量，当exact=T时，精确计算条件数，否则近似计算条件数。
#注意：一般条件数K<100,则认为多重共线性的程度很小；若100<=K<=1000则认为存在中等程度或较强的多重共线性；若K>1000则认为存在严重的多重共线性。
collinear<-data.frame(
Y=c(10.006, 9.737, 15.087, 8.422, 8.625, 16.289, 5.958, 9.313, 12.960, 5.541, 8.756, 10.937),
X1=rep(c(8, 0, 2, 0), c(3, 3, 3, 3)),
X2=rep(c(1, 0, 7, 0), c(3, 3, 3, 3)),
X3=rep(c(1, 9, 0), c(3, 3, 6)),
X4=rep(c(1, 0, 1, 10), c(1, 2, 6, 3)),
X5=c(0.541, 0.130, 2.116, -2.397, -0.046, 0.365, 1.996, 0.228, 1.38, -0.798, 0.257, 0.440),
X6=c(-0.099, 0.070, 0.115, 0.252, 0.017, 1.504, -0.865, -0.055, 0.502, -0.399, 0.101, 0.432) )
XX <- cor(collinear[2:7])
kappa(XX, exact = T)## [1] 2195.908#由上可知：得到条件数 K=2195.908>1000，认为有严重的多重共线性。
#找出哪些变量是多重共线性的，计算矩阵的特征值和相应的特征向量：
eigen(XX)## $values
## [1] 2.428787365 1.546152096 0.922077664 0.793984690 0.307892134 0.001106051
##
## $vectors
## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.3907189 0.33968212 0.67980398 -0.07990398 0.2510370
## [2,] -0.4556030 0.05392140 -0.70012501 -0.05768633 0.3444655
## [3,] 0.4826405 0.45332584 -0.16077736 -0.19102517 -0.4536372
## [4,] 0.1876590 -0.73546592 0.13587323 0.27645223 -0.0152087
## [5,] -0.4977330 0.09713874 -0.03185053 0.56356440 -0.6512834
## [6,] 0.3519499 0.35476494 -0.04864335 0.74817535 0.4337463
## [,6]
## [1,] -0.447679719
## [2,] -0.421140280
## [3,] -0.541689124
## [4,] -0.573371872
## [5,] -0.006052127
## [6,] -0.002166594
实例五、广义线性回归模型
Norell实验，高压电线对牲畜的影响
#1、加载数据
norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
#2、数据预处理
norell$Ymat<-cbind(norell$success, norell$n-norell$success)
#3、建模
glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
#4、模型评估
summary(glm.sol)##
## Call:
## glm(formula = Ymat ~ x, family = binomial, data = norell)
##
## Deviance Residuals:
## 1 2 3 4 5 6
## -2.2507 0.3892 -0.1466 1.1080 0.3234 -1.6679
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.3010 0.3238 -10.20 <2e-16 ***
## x 1.2459 0.1119 11.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 250.4866 on 5 degrees of freedom
## Residual deviance: 9.3526 on 4 degrees of freedom
## AIC: 34.093
##
## Number of Fisher Scoring iterations: 4#5、预测：与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
pre <- predict(glm.sol, data.frame(x=3.5))
p <- exp(pre)/(1+exp(pre))
#6、求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
glm.sol$coefficients## (Intercept) x
## -3.301035 1.245937X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]]
#7、画出响应比例与logistic回归曲线
d <- seq(0, 5, length=100)
pre <- predict(glm.sol, data.frame(x=d))
p <- exp(pre)/(1+exp(pre))
norell$y <- norell$success/norell$n
plot(norell$x, norell$y)
lines(d, p)#其中，d是给出曲线横坐标的点，pre是计算预测值，p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
