广义线性回归模型
广义线性回归模型
例题1 R.Norell实验
为研究高压电线对牲畜的影响，R.Norell研究小的电流对农场动物的影响。他在实验中，选择了7头，6种电击强度， 0,1,2,3,4,5毫安，每头牛被电击30下，每种强度5下，按随机的次序进行，然后重复整个实验，每头牛总共被电击60下。对每次电击，相应变量——嘴巴运动，或者出现，或者未出现。下表中的数据给出每种电击强度70次试验中响应的总次数。试分析电击对牛的影响
| 电流（毫安） | 试验次数 | 响应次数 | 响应的比例 |
|---|---|---|---|
| 0 | 70 | 0 | 0.000 |
| 1 | 70 | 9 | 0.129 |
| 2 | 70 | 21 | 0.300 |
| 3 | 70 | 47 | 0.671 |
| 4 | 70 | 60 | 0.857 |
| 5 | 70 | 63 | 0.900 |
解：用数据框形式输入数据，再构造矩阵，一列是成功（响应）的次数，另一列是失败（不响应）的次数，然后再作logistic回归。其程序如下（程序名：exam0619.R）
> norell<-data.frame(
+ x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63)
+ )
> norell$Ymat<-cbind(norell$success, norell$n-norell$success)
> glm.sol<-glm(Ymat~x, family=binomial, data=norell)
> 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
即
。并且回归方程通过了检验，因此回归模型为 
其中
是电流强度（单位：毫安）
与线性回归模型相同，在得到回归模型后，可以作预测，例如，当电流强度为3.5毫安时，有响应的牛的概率为多少？
> pre<-predict(glm.sol,data.frame(x=3.5))
> p<-exp(pre)/(1+exp(pre));p
1
0.742642
即74.26%
可以作控制，如有50%的牛有响应，其电流强度是多少？当P=0.5时，
，所以
> X<--glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]]
> X
[1] 2.649439
即2.65毫安的电流强度，可以使50%的牛有响应
最后画出响应的比例与logistic回归曲线。
> d<-seq(0,5,len=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)
广义线性模型建模函数qlm()
用法如下
> norell<-data.frame(
+ x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63)
+ )
> norell$Ymat<-cbind(norell$success, norell$n-norell$success)
> glm.sol<-glm(Ymat~x, family=binomial, data=norell)
> summary(glm.sol)