这一讲为实践内容，我们将学习如何使用R语言建立和拟合逻辑回归模型，并对模型和参数进行假设检验。

本讲默认读者已经了解R语言的基本语法，能够在自己的电脑中成功编写和运行R脚本。

一、单个0-1响应的逻辑回归模型

1.1 数据集介绍与模型建立

我们先从最简单的情况开始。这里我们使用faraway包中的wcgs数据集，该数据集的收集过程是：在1960年挑选了约3000个年龄在39到59岁的健康男性作为受访者，随访调查八年半，其间记录每个人每天的抽烟量，最后检查他们是否患上了心脏病。

这个数据集涉及的变量很多，这里我们只看这几列：chd是二分变量，表示是否患病；height是身高（英寸），cigs是每日吸烟量（支）：

> data(wcgs, package = "faraway")
> summary(wcgs[,c("chd","height","cigs")])
  chd           height           cigs     
 no :2897   Min.   :60.00   Min.   : 0.0  
 yes: 257   1st Qu.:68.00   1st Qu.: 0.0  
            Median :70.00   Median : 0.0  
            Mean   :69.78   Mean   :11.6  
            3rd Qu.:72.00   3rd Qu.:20.0  
            Max.   :78.00   Max.   :99.0

我们进行简单的可视化，看看在患病和不患病的受试者中，身高和吸烟量的分布情况：

library(ggplot2)
library(cowplot)
p1 <- ggplot(wcgs, aes(x=height, color=chd)) + 
	geom_histogram(position = "dodge", binwidth = 1)
p2 <- ggplot(wcgs, aes(x=cigs, color=chd)) +
	geom_histogram(position = "dodge", binwidth = 5, 
					aes(y=..density..))
plot_grid(p1, p2, ncol=2)

wcgs数据集可视化

从左图中可以看到，患病与不患病人群在身高上的分布较为类似；从右图中可以看到，吸烟量的增加使得患病人数反超未患病人数。但是只从图上观察是不准确的，我们需要建立模型来确认患病与否跟身高、吸烟量的关系。

接下来我们开始建立逻辑回归模型。逻辑回归模型是广义线性回归模型，所以类似于线性回归模型的lm()语句，广义线性回归也有内置的glm()语句。由于我们要建立逻辑回归模型，响应变量是二分变量，所以只需额外指定family参数为binomial即可：

> lmod <- glm(chd ~ height + cigs, family = binomial, wcgs)
> summary(lmod)
...
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.50161    1.84186  -2.444   0.0145 *  
height       0.02521    0.02633   0.957   0.3383    
cigs         0.02313    0.00404   5.724 1.04e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1781.2  on 3153  degrees of freedom
Residual deviance: 1749.0  on 3151  degrees of freedom
AIC: 1755

Number of Fisher Scoring iterations: 5

如果你学习过之前的所有内容，那么这个输出就很容易看懂了。Coefficients部分给出了回归系数的估计值；Deviance值有两个，一个是零模型/最小模型的deviance，一个是残差/回归模型deviance；AIC是信息准则；最后给出了Fisher Scoring算法的迭代次数，通常是4~8次。

对于逻辑回归模型来说，只知道回归系数是不够的，不具有实际含义。由于回归系数代表了对 $\log(\pi/1-\pi)$ 的贡献，所以我们把回归系数转成几率来进行解读：

> (beta <- coef(lmod))
(Intercept)      height        cigs 
-4.50161397  0.02520779  0.02312740 
> exp(beta)
(Intercept)      height        cigs 
 0.01109108  1.02552819  1.02339691 
> exp(beta[3]*20)
    cigs 
1.588115

可以看到，身高和吸烟量的几率分别为1.026和1.0233，都是大于1的，所以如果两个系数的检验都通过的话，我们就可以说，身高增加1英寸，患病几率增加2.6%；或者每天多吸1支烟，患病几率增加2.3%。另外，如果每天吸了20支烟，患病优势增加59%。

1.2 模型的检验

拟合完一个模型之后，最重要的一件事就是做检验。我们首先考虑每一个回归系数的检验。

你可以看到输出结果中Coefficients部分有一个z value列，而在线性回归中这一列是t value。在理论部分我们推导过：

$b\sim N(\beta,\mathfrak I^{-1})$

由于我们检验的目标是 $b_i=0$ ，所以进行标准化：

$z_i=\frac{b-0}{\sqrt{diag(\mathfrak I^{-1})_i}}\sim N(0,1)$

这就是z value列的由来，你仍然可以尝试把Estimate列和Std. Error列相除，就可以得到z value的结果。

为什么在逻辑回归中要使用z值（标准正态分布）而不是t值（t分布）？

这是因为， $\mathfrak I$ 是确定的，无需估计的。

如果你还记得理论知识： $\mathfrak I=X^TWX$ ，其中 $W$ 就是二项分布的方差，你就会发现方差是来自于结构（二项分布）而非样本，在大样本下，上述标准化就是服从正态分布的。

如果你回忆线性回归的回归系数的分布： $\hat\beta\sim N(\beta,\sigma^2(X^TX)^-)$ ，你就会发现 $\sigma$ 是未知的，于是我们使用了 $\hat\sigma$ 来代替了 $\sigma$ ，即用 $(n-r)\hat\sigma^2/n$ 估计 $\sigma$ 。既然我们对方差进行了估计，那么这种估计肯定是有误差的，这种误差会导致统计量存在拖尾的情况，因此使用t分布更为准确。

事实上，广义线性回归模型几乎都是使用z值来检验回归系数的。

接下来我们考虑模型的变量是否显著。我们之前学习的方法是，拟合不同的模型，然后调用anova()函数。这个方法在这里依然成立，你可以进行尝试，不过这里我们推荐一个更好的函数——drop1()，该函数可以自动地依次去掉一个自变量，然后与原模型做似然比检验，这样你就不需要手动搭建模型了：

> drop1(lmod, test="Chi")
Single term deletions

Model:
chd ~ height + cigs
       Df Deviance    AIC     LRT Pr(>Chi)    
<none>      1749.0 1755.0                     
height  1   1750.0 1754.0  0.9202   0.3374    
cigs    1   1780.1 1784.1 31.0695 2.49e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到，当去掉cigs项时似然比检验显著，说明cigs确实是必不可少的项；而height在去掉前后对模型没有显著的影响，所以我们可以认为身高和患病之间不存在相关性。

最后提一下模型结构检验。在理论部分我们学过，要检验模型结构是否存在就可以使用Deviance检验，也就是分别求全模型和拟合模型的对数似然值，然后做似然比检验。但是，在R实践当中，你不能这样做，你不能真的这样写：

glm(chd ~ factor(1:n), family=binomial, data=wcgs)

这样做是完全过拟合，往往不收敛甚至发生拟合错误。就算可以计算全模型的对数似然值，也无法直接检验，因为全模型的自由度为0。所以你会看到，summary(lmod)输出了Residual deviance: 1749.0 on 3151 degrees of freedom的字样，这里的Deviance值表示拟合模型和全模型的对数似然值的差，但是并不能做检验。

因此，一个更优的做法就是拟合最小模型/零模型，然后调用anova()函数：

> lmod_null <- glm(chd ~ 1, family=binomial, data=wcgs)
> anova(lmod_null, lmod)
Analysis of Deviance Table

Model 1: chd ~ 1
Model 2: chd ~ height + cigs
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1      3153     1781.2                          
2      3151     1749.0  2   32.195 1.021e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p值显著，说明模型结构是显著的。

1.3 模型的诊断

1.2小节完成了显著性的检验，现在我们还需要进行模型的诊断。

首先是拟合优度检验。你会发现上面1.1小节中的summary(lmod)里没有输出决定系数，这是因为线性回归中的决定系数无法直接套用在逻辑回归中。在实际应用当中，我们一般使用HL检验和伪决定系数。

要实现HL检验，可以直接调用ResourceSelection包中的hoslem.test()函数：

#install.packages("ResourceSelection")
> library(ResourceSelection)
> hoslem.test(lmod$y, fitted(lmod), g=10)
	Hosmer and Lemeshow goodness of fit (GOF) test

data:  lmod$y, fitted(lmod)
X-squared = 11.378, df = 8, p-value = 0.1812

我们把数据分成了10组，即g=10，结果是p值不显著，说明拟合优度偏差未达到显著性水平。

要计算伪决定系数，我们可以使用DescTools中的PseudoR2()函数来实现：

#install.packages("DescTools")
> library(DescTools)
> PseudoR2(lmod, which = "McFadden")
  McFadden 
0.01807417

这里展示了理论知识中提到的McFadden版本的伪决定系数，你还可以选择其他版本的伪决定系数，详情可以通过?PseudoR2来查询。

接下来看看残差。如果要获取每个样本的残差，可以使用residuals(lmod)。这里的残差默认是Deviance残差，如果你要使用Pearson残差，那就调用residuals(lmod, type="pearson")；如果你想调用的是观测值（0/1）和拟合值（ $\pi$ ）之间的“残差”，那就调用residuals(lmod, tupe="response")。获得残差以后，就可以做各种残差诊断，如正态性、等方差性、独立性等等，和线性回归没有区别。

这里还没提到过度分散问题，我们在第二节进行展开。

二、多个0-1响应的逻辑回归模型

上一节我们介绍了单个0-1响应的逻辑回归模型，接下来我们要拓展到多个0-1响应。多个0-1响应本质上就是理论学习部分的“分组处理”，即在某个 $\pi_i$ 水平下有多条数据，而不是以每条数据单独分组。这种处理类似于“批次”的概念，我们讨论的是每个批次的成功概率，而不是单个样本的成功与否。

2.1 数据集介绍与模型建立

1986年，航天飞机“挑战者号”在发射升空后爆炸。经过调查，工程师们把目光聚焦在火箭助推器的橡胶O型密封圈上。他们发现，在低温环境下，橡胶会变脆导致密封性变差。而在回顾之前已经完成的23次火箭发射任务，他们的确也发现了O型密封圈的损坏。那么是否可以把温度和密封圈的损坏联系起来呢？

faraway包中的orings数据集描述了每一个火箭的O型密封圈所处的温度和损坏的个数。每一个火箭有2个助推器，每个助推器有3个O型密封圈，也就是说，一个火箭有6个密封圈。damage列记录了每个火箭O型密封圈的损坏个数，temp列记录了温度。我们可以绘制O型密封圈损坏比例和温度的关系图：

> data(orings, package = "faraway")
> plot(damage/6 ~ temp, orings, xlim=c(25,85), ylim=c(0,1),
+      xlab="Temperature", ylab = "Prob of damage")

orings数据集可视化

在这个数据中，我们讨论的是每个火箭的O型密封圈损坏的比例，而不是每个密封圈是否损坏；另外，处于同一个火箭的密封圈的温度应该是一致的，所以没必要重复地使用同一水平的自变量。因此，这里我们采取n次二项分布来描述响应变量的分布。

对于多个0-1响应的逻辑回归模型，在搭建模型的步骤方面和单个0-1响应的模型是一致的，但是存在一点不同：对于多个0-1响应来说，公式部分的y不能只放damage，而是应该放一个两列的矩阵，第1列是事件成功次数（这里是发生损坏的次数），第2列是事件失败次数（这里是6-损坏次数），这样才能正确建模：

> lmod <- glm(cbind(damage,6-damage) ~ temp, family = binomial, orings)
> summary(lmod)
...
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 11.66299    3.29626   3.538 0.000403 ***
temp        -0.21623    0.05318  -4.066 4.78e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 38.898  on 22  degrees of freedom
Residual deviance: 16.912  on 21  degrees of freedom
AIC: 33.675

Number of Fisher Scoring iterations: 6

后续的分析步骤就和第一节完全一样，这里不再赘述，下面我们着重探讨过度分散的问题。

2.2 过度分散及其解决

在理论知识部分我们提到，逻辑回归有可能会出现过度分散问题，也就是说观测值的方差大于模型估计方差，使得模型不能完全解释数据中的方差波动。解决方法是引入一个分散参数 $\phi$ ，使得 $Var(Y_i)=\phi \pi_i(1-\pi_i)$ 。这个分散参数是可以手动估计的：

$\phi=\frac{\chi^2}{N-p}$

其中 $\chi^2$ 是模型拟合优度的Pearson卡方值。

在 $R$ 语言中，我们可以这样写：

sum(residuals(model, type="pearson")^2/(N-p))

上面是伪代码，在实际使用时需要替换model, N, p。如果你发现分散参数大于1，那说明模型存在过度分散问题。

我们使用2.1小节的例子，估计分散参数：

> sum(residuals(lmod, type = "pearson")^2/21)
[1] 1.336542

可以看到，上述搭建的模型是存在过度分散的。其实也不奇怪，从2.1小节的图中就可以看到，>0.5的数据和<0.5数据的量完全不对等。

要解决过度分散问题，我们可以使用拟二项分布（quasi-binomial）。要实现这一点很简单，只需要指定family参数为quasibinomial即可：

> lmod <- glm(cbind(damage,6-damage) ~ temp, family = quasibinomial, orings)
> summary(lmod)
...
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 11.66299    3.81077   3.061  0.00594 **
temp        -0.21623    0.06148  -3.517  0.00205 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 1.336542)

    Null deviance: 38.898  on 22  degrees of freedom
Residual deviance: 16.912  on 21  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

拟二项分布不会改变参数估计的结果，但是会改变方差，从而对显著性p值进行了修正。

三、广义逻辑回归模型

最后我们介绍广义逻辑回归模型，这里我们只考虑模型建立的问题，后续的检验和诊断就不再反复赘述。

3.1 多分类无序逻辑回归模型

多分类无序逻辑回归是指，联系函数仍然为逻辑函数，但是响应变量服从多项分布，且类与类之间不存在先后关系。这类模型之所以是广义逻辑回归模型，是因为它仍然使用了逻辑函数，并且我们在建模时，选取某一类为参考类，剩下的每一类都分别与参考类构成二项分布的关系。

我们以Rosenstone发表的关于1996年美国大选研究的数据为例，展示多分类无序逻辑回归模型的建立过程。为了简化，我们只考虑受访者的年龄、教育水平、经济收入水平，并探讨这些变量和他们加入的党派之间的关系。我们只考虑三类党派：民主党（Democrat）、独立派（Independent）、共和党（Republican），因此我们需要进行一定的数据重排：

> data(nes96, package = "faraway")
> party <- nes96$PID
> levels(party) <- c("Democrat", "Democrat", "Independent", "Independent",
+ "Independent", "Republican", "Republican")
> inca <- c(1.5,4,6,8,9.5,10.5,11.5,12.5,13.5,14.5,16,18.5,21,23.5,
+ 27.5,32.5,37.5,42.5,47.5,55,67.5,82.5,97.5,115)
> income <- inca[unclass(nes96$income)]
> rnes96 <- data.frame(party, income, education=nes96$educ, age=nes96$age)
> summary(rnes96)
         party         income        education        age       
 Democrat   :380   Min.   :  1.50   MS    : 13   Min.   :19.00  
 Independent:239   1st Qu.: 23.50   HSdrop: 52   1st Qu.:34.00  
 Republican :325   Median : 37.50   HS    :248   Median :44.00  
                   Mean   : 46.58   Coll  :187   Mean   :47.04  
                   3rd Qu.: 67.50   CCdeg : 90   3rd Qu.:58.00  
                   Max.   :115.00   BAdeg :227   Max.   :91.00  
                                    MAdeg :127

接下来我们来建立模型。要建立多分类无序逻辑回归，我们可以使用nnet包中的multinom()函数。nnet包是一个神经网络包，但是我们只需要使用神经网络训练器里面的优化方法来计算极大似然估计，不涉及到更深的网络。这样做的原因是，用这个优化方法会更快：

> library(nnet)
> mmod <- multinom(party ~ age + education + income, rnes96)
# weights:  30 (18 variable)
initial  value 1037.090001 
iter  10 value 990.568608
iter  20 value 984.319052
final  value 984.166272 
converged
> mmodi <- step(mmod)
...
> summary(mmodi)
Call:
multinom(formula = party ~ income, data = rnes96)

Coefficients:
            (Intercept)     income
Independent  -1.1749331 0.01608683
Republican   -0.9503591 0.01766457

Std. Errors:
            (Intercept)      income
Independent   0.1536103 0.002849738
Republican    0.1416859 0.002652532

Residual Deviance: 1985.424 
AIC: 1993.424

这种方法训练速度很快。我们通过step()进行了简单的筛选，最终只保留了income这个变量。

或者，如果你想要类似于glm()的写法，一个更通用的包是VGAM，使用其中的vglm()函数，指定family=multinomial即可：

#install.packages("VGAM")
library(VGAM)
mmod <- vglm(party ~ age + education + income, family = multinomial, rnes96)

3.2 多分类有序逻辑回归模型

如果响应变量不仅服从多项分布，每个类别之间还有顺序关系，那么我们就需要使用多分类有序逻辑回归模型了。在理论部分我们提到，多分类有序逻辑回归模型大致有三类：

比例几率模型（Proportional Odds Model, PO）
相邻类别模型（Adjacent Categories Model, AC）
连续比例模型（Continuation Ratio Model, CR）

在实际情况下，我们需要根据研究内容的特点来选择合适的模型。这里我们以PO模型为例。

我们仍然以3.1小节中构造的数据集为例，我们假设三个党派是有顺序关系的（虽然事实上并没有，这里只是为了简便），例如，我们假设独立派在民主党与共和党之间。要建立PO模型，我们使用MASS包中的polr()函数来完成：

> library(MASS)
> pomod <- polr(party ~ age + education + income, rnes96)
> summary(pomod)
Re-fitting to get Hessian
...
Coefficients:
                Value Std. Error  t value
age          0.005775   0.003887  1.48581
education.L  0.724087   0.384388  1.88374
education.Q -0.781361   0.351173 -2.22500
education.C  0.040168   0.291762  0.13767
education^4 -0.019925   0.232429 -0.08573
education^5 -0.079413   0.191533 -0.41462
education^6 -0.061104   0.157747 -0.38735
income       0.012739   0.002140  5.95187

Intercepts:
                       Value   Std. Error t value
Democrat|Independent    0.6449  0.2435     2.6479
Independent|Republican  1.7374  0.2493     6.9694

Residual Deviance: 1984.211 
AIC: 2004.211

如果你对模型使用step()，你会发现还是只跟income有关。

当然，VGAM包更全面，你仍然只需要使用vglm()函数，并且指定family = cumulative(link="logit")即可。另外，对于AC模型，可以参考VGAM包中的acat()；对于CR模型，可以参考VGAM包中的cratio()。所以说，VGAM包算是更为通用的包。

3.3 二项响应的不同联系函数

最后我们来看二项响应的不同联系函数。虽然使用别的联系函数已经不能称之为”逻辑回归“了，但是为了完整性，我们还是要介绍如何使用其他联系函数。不过请注意，这里所谓的”其他联系函数“要和响应变量的二项分布相适应，不是随便用一个联系函数就可以。

如果响应变量服从二项分布，那么除了逻辑函数，还有常见的比如Probit模型（标准正态分布连接）、log-log模型（两层对数连接）、cauchit模型（Cauchy分布连接）等。在glm()的family参数处，binomial默认是logit连接，也就是逻辑回归。要使用其他联系函数只需要在这个基础上进行修改即可。

我们以bliss数据集为例，该数据集记录了不同杀虫剂浓度对昆虫死亡比例的影响：

> data(bliss, package = "faraway")
> bliss
  dead alive conc
1    2    28    0
2    8    22    1
3   15    15    2
4   23     7    3
5   27     3    4

> mlogit <- glm(cbind(dead,alive)~conc, family=binomial, data=bliss)
> mprobit <- glm(cbind(dead,alive)~conc, family=binomial(link=probit), data=bliss)
> mcloglog <- glm(cbind(dead,alive)~conc, family=binomial(link=cloglog), data=bliss)
> mcauchit <- glm(cbind(dead,alive)~conc, family=binomial(link=cauchit), data=bliss)

我们来对比一下四个模型的拟合值：

> predval <- sapply(list(mlogit,mprobit,mcloglog,mcauchit), fitted)
> dimnames(predval) <- list(0:4,c("logit","probit","cloglog","cauchit"))
> round(predval,3)
  logit probit cloglog cauchit
0 0.089  0.084   0.127   0.119
1 0.238  0.245   0.250   0.213
2 0.500  0.498   0.455   0.506
3 0.762  0.752   0.722   0.791
4 0.911  0.914   0.933   0.882

看起来它们的区别是有但不多的。如何选择这几个联系函数不能单单靠数据本身，有时候还需要实际经验、领域知识以及便利性。比如在生物医学和临床队列方面，可能就用逻辑回归模型是最佳的，因为它数学表达简单、容易解释OR值、适合前瞻性和回顾性研究；对于一些金融、博弈方面的模型，或许用probit或log-log模型更好。