引言：第8章（回归）和第9章（方差分析）中，我们探究了线性模型，它们可以通过一系列连续型和/或类别型预测变量来预测正态分布的响应变量。但在许多情况下，假设因变量为正态分布（甚至连续型变量）并不合理，比如结果变量可能是类别型的。二值变量（是/否、通过/未通过、活着/死亡）和多分类变量（比如差/良好/优秀）都显然不是正态分布。 或者是计数型的（比如，一周交通事故的数目，每日酒水消耗的数量）。这类变量都是非负的有限值，而且它们的均值和方差通常都是相关的（正态分布变量间不是如此，而是相互独立）。广义线性模型扩展了线性模型的框架，它包含了非正态因变量的分析。本章我们主要学习三部分，分别是建立广义线性模型；预测类别型变量和计数型数据建模。
13.1 广义线性模型和 glm() 函数
广义线性模型扩展了线性模型的框架，它包含了非正态的因变量分析：
广义线性模型拟合形式：
其中g(μY)是条件均值的函数（称为连接函数），广义线性模型是通过拟合响应变量的条件均值的一个函数（不是响应变量的条件均值），如果响应变量服从指数分布族中的某个分布（并不仅限于正态分布），极大地扩展了标准线性模型。模型参数估计的推导依据的是极大似然估计，而非最小二乘法。
13.1.1 glm() 函数
R中可通过glm() 函数（还可用其他专门的函数）拟合广义线性模型。形式主要是：
glm( f o r m u l a, family= f a m i l y(link= f u n c t i o n), data=)
表13-1列出了概率分布（ f am i l y）和相应默认的连接函数（ f u n c ti o n）
表13-1 g l m ( ) 的参数
glm()函数可以拟合许多流行的模型， 比如Logistic回归、泊松回归和生存分析（此处不考虑）
假设你有一个响应变量（ Y）、三个预测变量（ X1、 X2 、 X3 ）和一个包含数据的数据框（ mydata），下面对前两个模型进行阐述：
１、Logistic回归适用于二值响应变量（ 0和1 ），模型假设Y服从二项分布，
如下代码拟合Logistic回归模型：
glm(Y~X1+X2+X3, family=binomial(link="logit"), data=mydata)
２、泊松回归适用于在给定时间内响应变量为事件发生数目的情形。它假设Y服从泊松分布，
可用如下代码拟合泊松回归模型：
glm(Y~X1+X2+X3, family=poisson(link="log"), data=mydata)
３、标准线性模型也是广义线性模型的一个特例，如果令连接函数g(μY)=μY或恒
等函数，并设定概率分布为正态（高斯）分布，那么
glm(Y~X1+X2+X3, family=gaussian(link="identity"), data=mydata)二者生成的结果相同：
lm(Y~X1+X2+X3, data=mydata)
13.1.2 连用的函数
与分析标准线性模型时lm() 连用的许多函数在glm() 中都有对应的形式，如下表１３－２
表13-2 与 glm () 连用的函数
| 函数 | 描述 | |
| summary() | 展示拟合模型的细节 | |
| coefficients()/coef() | 列出拟合模型的参数（截距项和斜率） | |
| confint() | 给出模型参数的置信区间（默认为95%） | |
| residuals() | 列出拟合模型的残差值 | |
| anova() | 生成两个拟合模型的方差分析表 | |
| plot() | 生成评价拟合模型的诊断图 | |
| predict() | 用拟合模型对新数据集进行预测 | |
| deviance() | 拟合模型的偏差 | |
| df.residual() | 拟合模型的残差自由度 | |
13.2 Logistic 回归
当通过一系列连续型和/或类别型预测变量来预测二值型结果变量时， Logistic回归是一个非
常有用的工具。以AER包中的数据框Affairs为例：
> install.packages("AER") ＃下载和安装了ＡＥＲ包> data(Affairs, package="AER")> summary(Affairs)＃一些描述性的统计信息affairs gender age yearsmarried childrenMin. : 0.000 female:315 Min. :17.50 Min. : 0.125 no :1711st Qu. : 0.000 male :286 1st Qu.:27.00 1st Qu.: 4.000 yes:430Median : 0.000 Median :32.00 Median : 7.000Mean : 1.456 Mean :32.49 Mean : 8.1783rd Qu. : 0.000 3rd Qu.:37.00 3rd Qu.:15.000Max. :12.000 Max. :57.00 Max. :15.000Religiousness education occupation ratingMin. :1.000 Min. : 9.00 Min. :1.000 Min. :1.0001st Qu. :2.000 1st Qu.:14.00 1st Qu.: 3.000 1st Qu.:3.000Median :3.000 Median : 16.00 Median :5.000 Median :4.000Mean :3.116 Mean : 16.17 Mean :4.195 Mean :3.9323rd Qu. :4.000 3rd Qu.: 18.00 3rd Qu.:6.000 3rd Qu.:5.000Max. :5.000 Max. : 20.00 Max. :7.000 Max. :5.000> table(Affairs$affairs)0 1 2 3 7 12451 34 17 19 42 38> Affairs$ynaffair[Affairs$affairs > 0] <- 1> Affairs$ynaffair[Affairs$affairs == 0] <- 0> Affairs$ynaffair <- factor(Affairs$ynaffair,levels=c(0,1),labels=c("No","Yes"))> table(Affairs$ynaffair)＃将affairs转化为二值型因子ynaffairNo Yes451 150> fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation +rating, data=Affairs, family=binomial())> summary(fit.full)＃展示ｌｏｇｉｓｔｉｃ回归的结果变量Call:glm(formula = ynaffair ~ gender + age + yearsmarried + children +religiousness + education + occupation + rating, family = binomial(),data = Affairs)Deviance Residuals:Min 1Q Median 3Q Max-1.571 -0.750 -0.569 -0.254 2.519Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) 1.3773 0.8878 1.55 0.12081gendermale 0.2803 0.2391 1.17 0.24108age -0.0443 0.0182 -2.43 0.01530 *yearsmarried 0.0948 0.0322 2.94 0.00326 **childrenyes 0.3977 0.2915 1.36 0.17251religiousness -0.3247 0.0898 -3.62 0.00030 ***education 0.0211 0.0505 0.42 0.67685occupation 0.0309 0.0718 0.43 0.66663rating -0.4685 0.0909 -5.15 2.6e-07 ***---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(Dispersion parameter for binomial family taken to be 1)Null deviance: 675.38 on 600 degrees of freedomResidual deviance: 609.51 on 592 degrees of freedomAIC: 627.5Number of Fisher Scoring iterations: 4
从回归系数的p值（最后一栏）发现性别、是否有孩子、学历和职业对方程的贡献都不显著（你无法拒绝参数为0的假设）。去除这些变量重新拟合模型：
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data=Affairs, family=binomial())summary(fit.reduced)Call:= ynaffair ~ age + yearsmarried + religiousness + rating,family = binomial(), data = Affairs)Deviance Residuals:Min 1Q Median 3Q Max-0.755 -0.570 -0.262 2.400Coefficients:Estimate Std. Error z value Pr(>|z|)(Intercept) 1.9308 0.6103 3.16 0.00156 **age -0.0353 0.0174 -2.03 0.04213 *yearsmarried 0.1006 0.0292 3.44 0.00057 ***religiousness -0.3290 0.0895 -3.68 0.00023 ***rating -0.4614 0.0888 -5.19 2.1e-07 ***---codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1parameter for binomial family taken to be 1)Null deviance: 675.38 on 600 degrees of freedomResidual deviance: 615.36 on 596 degrees of freedomAIC: 625.4Number of Fisher Scoring iterations: 4
新模型的每个回归系数都非常显著（ p<0.05），使用anova() 函数对它们进行比较，对于广义线性回归，可用卡方检验：
anova(fit.reduced, fit.full, test="Chisq")Analysis of Deviance TableModel 1: ynaffair ~ age + yearsmarried + religiousness + ratingModel 2: ynaffair ~ gender + age + yearsmarried + children +religiousness + education + occupation + ratingResid. Df Resid. Dev Df Deviance P(>|Chi|)1 596 6152 592 610 4 5.85 0.21
结果的卡方值不显著（ p=0.21 ），表明四个预测变量的新模型与九个完整预测变量的模型拟
合程度一样好，因此我们可以选择更为简单的第二个模型fit.reduced。
13.2.1 解释模型参数
回归系数：当其他预测变量不变时，一单位预测变量的变化可引起的响应变量对数优势比的变化。
在Logistic回归中，响应变量是Y=1的对数优势比（ log）
接着求上面回归系数：
> coef(fit.reduced)(Intercept) age yearsmarried religiousness rating1.931 -0.035 0.101 -0.329 -0.461
由于对数优势比解释性差，你可对结果进行指数化：
> exp(coef(fit.reduced))(Intercept) age yearsmarried religiousness rating6.895 0.965 1.106 0.720 0.630
解读：婚龄增加一年，婚外情的优势比将乘以1.106；年龄增加一岁，婚外情的的优势比则乘以0.965；对于二值型Logistic回归，比起预测变量一单位的变化我们更关注某预测变量n单位的变化引起的较高值上优势比的变化为exp(βj)^n，比如，保持其他预测变量不变，婚龄增加一年，婚外情的优势比将乘以1.106，而如果婚龄增加10年，优势比将乘以1.106^10，即2.7。
13.2.2 评价预测变量对结果概率的影响
使用predict() 函数，可以观察某个预测变量在各个水平时对结果概率的影响。
testdata <- data.frame(rating=c(1, 2, 3, 4, 5), age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))#创建一个虚拟数据集，设定年龄、婚龄和宗教信仰为它们的均值，婚姻评分的范围为1～ 5testdatarating age yearsmarried religiousness1 1 32.5 8.18 3.122 2 32.5 8.18 3.123 3 32.5 8.18 3.123.125 5 32.5 8.18 3.12testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")#用测试数据集预测相应的概率> testdata #当婚姻评分从1（很不幸福）变为5（非常幸福）时，婚外情概率从0.53降低到了0.15（假定年龄、婚龄和宗教信仰不变）rating age yearsmarried religiousness prob8.18 3.12 0.5308.18 3.12 0.4168.18 3.12 0.3100.2203.12 0.151
用同样的方法看看年龄的影响：
> testdata <- data.frame(rating=mean(Affairs$rating),age=seq(17, 57, 10),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))> testdatarating age yearsmarried religiousness1 3.93 17 8.18 3.122 3.93 27 8.18 3.123 3.93 37 8.18 3.124 3.93 47 8.18 3.125 3.93 57 8.18 3.12> testdata$prob <- predict(fit.reduced, newdata=testdata, type="response")> testdata #当其他变量不变，年龄从17增加到57时，婚外情的概率将从0.34降低到0.11。rating age yearsmarried religiousness prob1 3.93 17 8.18 3.12 0.3352 3.93 27 8.18 3.12 0.2623 3.93 37 8.18 3.12 0.1994 3.93 47 8.18 3.12 0.1495 3.93 57 8.18 3.12 0.109
13.2.3 过度离势
过度离势，即观测到的响应变量的方差大于期望的二项分布的方差。过度离势会导致奇异的标准误检验和不精确的显著性检验。
1、当出现过度离势时，仍可使用glm() 函数拟合Logistic回归，但此时需要将二项分布改为类二项分布（ quasibinomial distribution）。
2、检测过度离势
方法一：比较二项分布模型的残差偏差与残差自由度a，如果比值比1大很多，你便可认为存在过度离势。
> deviance(fit.reduced)/df.residual(fit.reduced)#结果接近于1，表明没有过度离势[1] 1.032
方法二：分别拟合模型两次，family=binomial” 和family=”quasibinomial”，提供的p值即可对零假设H0：b=1与备择假设H1： b≠1进行检验。若p很小（小于0.05），你便可拒绝零假设。
pchisq(summary(fit.od)$dispersion * fit$df.residual, fit$df.residual, lower = F)应用到上面的例子中：
> fit <- glm(ynaffair ~ age + yearsmarried + religiousness +rating, family = binomial(), data = Affairs）#拟合模型第一次> fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +rating, family = quasibinomial(), data = Affairs)#拟合模型第二次> pchisq(summary(fit.od)$dispersion * fit$df.residual,fit$df.residual, lower = F)#此处p值（ 0.34）显然不显著（ p>0.05），接受原假设，即是认为不存在过度离势[1] 0.34
13.2.4 扩展
- 稳健Logistic回归 ：当拟合Logistic回归模型数据出现离群点和强影响点时
- 多项分布回归：若响应变量包含两个以上的无序类别（比如，已婚/寡居/离婚
），用mlogit包中的mlogit() 函数拟合多项Logistic回归 -
序数Logistic回归：若响应变量是一组有序的类别（比如，信用风险为差/良/好）
），用rms包中的lrm() 函数拟合序数Logistic回归。
13.3 泊松回归
当通过一系列连续型和/或类别型预测变量来预测计数型结果变量时，可使用泊松回归。
例子：robust包中的Breslow癫痫数据，研究受轻微或严重间歇性癫痫的病人在接受药物治疗前后的发作情况，以便研究药物治疗是否能减少癫痫发病数。响应变量为sumY（随机化后八周内癫痫发病数），预测变量为治疗条件（ Trt）、年龄（ Age）和前八周内的基础癫痫发病数（ Base）。
> data(breslow.dat, package="robust")#加载robust包> names(breslow.dat)[1] "ID" "Y1" "Y2" "Y3" "Y4" "Base" "Age" "Trt" "Ysum"[10] "sumY" "Age10" "Base4"> summary(breslow.dat[c(6,7,8,10)])#选取之前描述的四个变量Base Age Trt sumYMin. : 6.0 Min. :18.0 placebo : 28 Min. : 0.01st Qu.: 12.0 1st Qu.: 23.0 progabide:31 1st Qu.: 11.5Median : 22.0 Median : 28.0 Median : 16.0Mean : 31.2 Mean : 28.3 Mean : 33.13rd Qu.: 41.0 3rd Qu.: 32.0 3rd Qu.: 36.0Max. : 151.0 Max. : 42.0 Max. : 302.0
对以上的数据进行可视化，可以清楚地看到因变量的偏倚特性以及可能的离群点
如图13-1：
> opar <- par(no.readonly=TRUE)> par(mfrow=c(1,2))> attach(breslow.dat)> hist(sumY, breaks=20, xlab="Seizure Count",main="Distribution of Seizures")#建立sumY的直方图> boxplot(sumY ~ Trt, xlab="Treatment", main="Group Comparisons")#建立sumY关于Trt的箱式图> par(opar)
图13-1 随机化后的癫痫发病数的分布情况
接下来拟合泊松回归：
fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson())summary(fit)#输出结果列出了偏差、回归参数、标准误和参数为0的检验；而且此处预测变量在p<0.05的水平下都非常显著Call:= sumY ~ Base + Age + Trt, family = poisson(), data = breslow.dat)Deviance Residuals:Min 1Q Median 3Q Max-2.043 -0.940 0.793 11.006Coefficients:Estimate Std. Error z value Pr(>|z|)1.948826 0.135619 14.37 < 2e-16 ***Base 0.022652 0.000509 44.48 < 2e-16 ***Age 0.022740 0.004024 5.65 1.6e-08 ***Trtprogabide -0.152701 0.047805 -3.19 0.0014 **---codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1parameter for poisson family taken to be 1)Null deviance: 2122.73 on 58 degrees of freedomResidual deviance: 559.44 on 55 degrees of freedomAIC: 850.7Number of Fisher Scoring iterations: 5
13.3.1 解释模型参数
使 用 coef() 函 数 可 获 取 模 型 系 数 ， 或 者 调 用 summary() 函 数 的 输 出 结 果 中 的Coefficients表格：
> coef(fit)(Intercept) Base Age Trtprogabide1.9488 0.0227 0.0227 -0.1527
解读：在泊松回归中，因变量以条件均值的对数形式loge(λ)来建模。年龄的回归参数为0.0227，表示保持其他预测变量不变，年龄增加一岁，癫痫发病数的对数均值将相应增加0.03。通常在因变量的初始尺度（癫痫发病数，而非发病数的对数）上解释回归系数比较容易。为此，指数化系数：
> exp(coef(fit))(Intercept) Base Age Trtprogabide7.020 1.023 1.023 0.858
年龄增加一岁，期望的癫痫发病数将乘以 1.023；而对于期望的癫痫发病数将乘以0.86，也就是药物治疗后癫痫发病数降低了20%。注意，与Logistic回归中的指数化参数相似，泊松模型中的指数化参数对响应变量的影响都是成倍增加的，而不是线性相加。
13.3.2 过度离势
泊松分布的方差和均值相等。当响应变量观测的方差比依据泊松分布预测的方差大时，泊松回归可能发生过度离势，发生过度离势的原因主要有遗漏了某个重要的预测变量、事件相关、在纵向数据分析中，重复测量数据等等。
检验方法一：
> deviance(fit)/df.residual(fit)#[1] 10.17
检验方法二，qcc包的cc.overdispersion.test（）函数：
> library(qcc)> qcc.overdispersion.test(breslow.dat$sumY, type="poisson")#p值果然小于0.05，进一步表明确实存在过度离势Overdispersion test Obs.Var/Theor.Var Statistic p-valuepoisson data 62.9 3646 0
一旦存在过度离势，通过用family=”quasipoisson” 替换family=”poisson”仍可以用glm() 函数对该数据进行拟合。
fit.od <- glm(sumY ~ Base + Age + Trt, data=breslow.dat,family=quasipoisson())summary(fit.od)Call:= sumY ~ Base + Age + Trt, family = quasipoisson(),data = breslow.dat)Deviance Residuals:Min 1Q Median 3Q Max-2.043 -0.940 0.793 11.006Coefficients:Estimate Std. Error t value Pr(>|t|)1.94883 0.46509 4.19 0.00010 ***0.02265 0.00175 12.97 < 2e-16 ***Age 0.02274 0.01380 1.65 0.10509Trtprogabide -0.15270 0.16394 -0.93 0.35570---codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1parameter for quasipoisson family taken to be 11.8)Null deviance: 2122.73 on 58 degrees of freedomResidual deviance: 559.44 on 55 degrees of freedomAIC: NANumber of Fisher Scoring iterations: 5
注意：使用类泊松（ quasi-Poisson）方法所得的参数估计会使标准误变大。此处，标准误越大将会导致Trt（和Age）的p值越大于0.05。当考虑过度离势，并控制基础癫痫数和年龄时，并没有充足的证据表明药物治疗相对于使用安慰剂能显著降低癫痫发病次数。
13.3.3 扩展 （若有兴趣详细内容可进一步参考《R语言实战》）
1. 时间段变化的泊松回归
2. 零膨胀的泊松回归：pscl包中的zeroinfl() 函数
3. 稳健泊松回归：robust包中的glmRob() 函数
小结：本章使用广义线性模型扩展了经典方法的方法适用范围。具体来讲，该框架可以用来分析非正态的响应变量，包括分类数据和离散的计数型数据。我们重点学习了Logistic回归（分析二值型结果变量）和泊松回归（分析计数型或比率结果变量）并且讨论了过度离势问题，包括如何检测以及依据它进行调整等方法。最后，学习了它们在R中的一些可用扩展和变种。到目前为止，我们介绍的每种统计方法都是直接处理可观测、可记录的变量。在下一章中，我们将看看处理潜变量的统计模型，即那些你坚信存在并能解释可观测变量的、无法被观测到的、理论上的变量。具体而言，你将学习如何使用因子分析方法检测和检验这些无法被观测到的变量的假设。
