R语言系列第五期：②R语言与逻辑回归建立
但医学工作者最常接触的结局预测变量多为二分类变量，比如阳性、阴性，病例、对照乃至生存、死亡这样的变量。这样我们就可以描述或推测在某些不同状况下得某种疾病的风险或者说阳性时间发生的概率。这里自然而然就引入我们今天的主题：逻辑回归模型——logistic regression model。
关于逻辑回归模型，需要注意的是，他与线性模型不同，没有误差项。我们是对一个事件发生的概率直接建模，而二元输出的变异性将由此概率来确定。因此，与正态分布不同，这里没有方差这个参数。
我们这里按照数据的原始类型分类来讲解不同的原始数据应该怎样通过R语言建立逻辑回归模型。
这一个部分里，我们采用Altman收集的高血压数据：
> no.yes<-c("No","Yes")
> smoking<-gl(2,1,8,no.yes)
> obesity<-gl(2,2,8,no.yes)
> snoring<-gl(2,4,8,no.yes)
> n.tot<-c(60,17,8,2,187,85,51,23)
> n.hyp<-c(5,2,1,0,35,13,15,8)
> data.frame(smoking,obesity,snoring,n.tot,n.hyp)
smoking obesity snoring n.tot n.hyp
1 No No No 60 5
2 Yes No No 17 2
3 No Yes No 8 1
4 Yes Yes No 2 0
5 No No Yes 187 35
6 Yes No Yes 85 13
7 No Yes Yes 51 15
8 Yes Yes Yes 23 8
#Tips：gl()函数用来“生成水平”，专门为平衡的实验设计而出现，它有三个参数：水平的数目，每块长度（每一水平需要重复多少次），以及结果的总长度。，第四个参数用来指定所生成的因子的水平名称。而把这些变量放到一个数据框中，输出更加直观好看。
对于表格化的数据进行逻辑回归分析，在R中有两种途径。你需要将数据表示成一个矩阵，其中一列是“患病”的个数，一列是“健康”的个数（或者“成功”、“失败”，基于自己的场景而定。）
> hyp.tbl<-cbind(n.hyp,n.con=n.tot-n.hyp)
> hyp.tbl
n.hyp n.con
[1,] 5 55
[2,] 2 15
[3,] 1 7
[4,] 0 2
[5,] 35 152
[6,] 13 72
[7,] 15 36
[8,] 8 15
#Tips：cbind()函数将变量按列合并在一起，形成一个矩阵，而我们将发生高血压的患者和正常人分列成两列，才能完成我们的逻辑回归（一定要注意，不能用合计的样本数代替对照组，不然会发生错误）。
随后，我们就可以建立逻辑回归模型了：
> glm(hyp.tbl~smoking+obesity+snoring,family=binomial("logit"))
Call: glm(formula = hyp.tbl ~ smoking + obesity + snoring, family = binomial("logit"))
Coefficients:
(Intercept) smokingYes obesityYes snoringYes
-2.37766 -0.06777 0.69531 0.87194
Degrees of Freedom: 7 Total (i.e. Null); 4 Residual
Null Deviance: 14.13
Residual Deviance: 1.618 AIC: 34.54
或者：
> glm(hyp.tbl~smoking+obesity+snoring,binomial)
另外一种建立逻辑回归模型的方法是给出每个水平组合中得病数的占比以及当前水平组合的总数：
> prop.hyp=n.hyp/n.tot
> glm.hyp<-glm(prop.hyp~smoking+obesity+snoring,binomial,weights=n.tot)
#Tips：结果是跟上面的输出结果完全一致。注意这里的weights参数是必须的，因为R无法识别这个占比所基于的基数是多少。其实这两种方法都是一样的，主要是看你有什么样子的数据。另外glm()是建立广义线性模型的函数。逻辑回归要采用的就是广义线性模型。当然这个glm()函数不止能建立逻辑回归模型，其他的模型我们这里不做详细介绍。
#Tips：如果仅仅从当前结果看的话，我们可以得出smoking的系数为-0.0677，而obesity和snoring的系数都是大于0的，说明吸烟可能是高血压的保护因素，后两个是高血压的危险因素，但是我们并不知道它们具不具有显著性。
同样我们可以使用summary()函数来提取必要的信息：
> summary(glm.hyp)
Call:
glm(formula = prop.hyp ~ smoking + obesity + snoring, family = binomial,
weights = n.tot)
Deviance Residuals:
1 2 3 4 5 6 7 8
-0.04344 0.54145 -0.25476 -0.80051 0.19759 -0.46602 -0.21262 0.56231
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.37766 0.38018 -6.254 4e-10 ***
smokingYes -0.06777 0.27812 -0.244 0.8075
obesityYes 0.69531 0.28509 2.439 0.0147 *
snoringYes 0.87194 0.39757 2.193 0.0283 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 14.1259 on 7 degrees of freedom
Residual deviance: 1.6184 on 4 degrees of freedom
AIC: 34.537
Number of Fisher Scoring iterations: 4
#Tips：划线的部分是我们最感兴趣的表，这张表中给出了回归系数的估计、标准差以及每一个系数显著性的假设检验结果。这里的结果告诉我们，尽管smoking的系数是负的，但是不具有统计学意义，而后两个变量的P值是小于0.05的，所以可以大胆地说肥胖和打鼾与高血压是正相关的。当然，这种情况下，我们会去掉smoking变量，重新进行模型的建立。
我们同样采用juul数据集，首先我们要把这个数据集里的分类变量转化成因子以便后续计算：
> library(ISwR)
> juul$menarche<-factor(juul$menarche,labels=c("No","Yes"))
> juul$tanner<-factor(juul$tanner)
接下来我们看看作为应变量的“menarche”。这个变量指示每个女孩是否已经经历了月经初潮。其被编码为1和2，分别表示“没有”和“有”，查看8—20岁女孩的数据很方便，可以通过如下代码实现：
> juul.girl<-subset(juul,age>8 & age<20 & complete.cases(menarche))
> attach(juul.girl)
#Tips：显然男孩子不会有这个生理现象，因此没必要对性别进行区分，我们使用complete.cases()函数和相应的年龄区间条件选择把menarche的有效观测给输出到juul.girl的数据集里。
然后，我们可以分析初潮与年龄的关系，代码如下：
> summary(glm(menarche~age,binomial))
Call:
glm(formula = menarche ~ age, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.32759 -0.18998 0.01253 0.12132 2.45922
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -20.0132 2.0284 -9.867 <2e-16 ***
age 1.5173 0.1544 9.829 <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: 719.39 on 518 degrees of freedom
Residual deviance: 200.66 on 517 degrees of freedom
AIC: 204.66
Number of Fisher Scoring iterations: 7
#Tips：应变量“menarche”是一个两水平的因子，第二个水平表示事件发生，当然如果变量被编码成0和1也是可以的。而R做的就是以小的数字做参照，来计算大的数字发生的概率（有参数可以设置那个值作为参照）。我们计算一下这个群体月经初潮年龄的预期中位数（P=0.5），其实就是logit P=0的年龄。大概是13.19岁（1.5173*age-20.0132=0）
再复杂一点，我们可以引入青春期分期变量tanner变量，tanner变量是一个分类变量，这件事我们之前已经告诉过R，所以R将它进行哑变量化处理：
> summary(glm(menarche~age+tanner,binomial))
Call:
glm(formula = menarche ~ age + tanner, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.56180 -0.12461 0.02475 0.08055 2.86120
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -13.7758 2.7630 -4.986 6.17e-07 ***
age 0.8603 0.2311 3.723 0.000197 ***
tanner2 -0.5211 1.4846 -0.351 0.725609
tanner3 0.8264 1.2377 0.668 0.504313
tanner4 2.5645 1.2172 2.107 0.035132 *
tanner5 5.1897 1.4140 3.670 0.000242 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 604.19 on 435 degrees of freedom
Residual deviance: 106.60 on 430 degrees of freedom
(83 observations deleted due to missingness)
AIC: 118.6
Number of Fisher Scoring iterations: 8
#Tips：这里的tanner变量以1作为参照，有的哑变量是有统计学意义的，有的则没有，但是这些变量必须同进同出，如果你想查看一下这个变量总体是否有统计学意义，可以做一下联合检验：
> drop1(glm(menarche~age+tanner,binomial),test="Chisq")
Single term deletions
Model:
menarche ~ age + tanner
Df Deviance AIC LRT Pr(>Chi)
<none> 106.60 118.60
age 1 124.50 134.50 17.901 2.327e-05 ***
tanner 4 161.88 165.88 55.282 2.835e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：这里drop1()函数是剔除掉某个变量之后做设定的检验。英文解释：Compute all the single terms in the scope argument that can be dropped from the model, fit those models and compute a table of the changes in fit.显然这两个变量都是有统计学意义的变量。
关于逻辑回归模型建立的部分我们已经介绍完了，根据我们数据类型分为表格类型数据和原始数据，两种数据的输入方式是不同，下面一个部分会为大家介绍逻辑回归模型的预测和检验。敬请期待。
参考资料：
1.《R语言统计入门（第二版）》人民邮电出版社 Peter Dalgaard著
2.《R语言初学者指南》人民邮电出版社 Brian Dennis著
3.Vicky的小笔记本《blooming for you》by Vicky
生信发文助手
多点好看，少点脱发