R语言统计系列第11篇-Logistic回归
今天是各类统计方法R语言实现的第11期，我们主要介绍Logistic回归。Logistic回归属于广义线性回归，因此我们从广义线性回归讲起。
广义线性回归
线性回归模型要求因变量服从正态分布，但是当结果变量是分类型（有无，患病与否等，二分类常用Logistic回归）、计数型（某地区某年发生肿瘤患者的人数等，常用泊松回归）或者临床上经常使用的无复发生存期数据等（常用Cox回归），因变量不符合正态分布，无法直接使用线性回归。而广义线性模型扩展了线性模型的框架，可以进行非正态因变量的分析，在R语言中可以通过glm()函数实现。
glm()函数的参数
| 分布族 | 默认的连接函数 |
|---|---|
| binomial | (link = “logit”) |
| gaussian | (link = “identity”) |
| gamma | (link = “inverse”) |
| inverse.gaussian | (link = “1/mu^2”) |
| poisson | (link = “log”) |
| quasi | (link = “identity”, variance = “constant”) |
| quasibinomial | (link = “logit”) |
| quasipoisson | (link = “log”) |
连用的函数
| 函数 | 描述 |
|---|---|
| summary() | 展示拟合模型的细节 |
| coefficients(), coef() | 列出拟合模型的参数（截距项和斜率） |
| confint() | 给出模型参数的置信区间（默认为95%） |
| residuals() | 列出拟合模型的残差值 |
| anova() | 生成两个拟合模型的方差分析表 |
| plot() | 生成评价拟合模型的诊断图 |
| predict() | 用拟合模型对新数据集进行预测 |
Logistic回归
二分类因变量常用Logistic回归，假设因变量Y服从二项分布，查表得(link = “logit”)
此处是一份婚外情数据，我们用性别、年龄等因素预测参与者是否发生婚外情affairs。
# 载入数据
data(Affairs, package = "AER")
summary(Affairs)
## affairs gender age yearsmarried children
## Min. : 0.000 female:315 Min. :17.50 Min. : 0.125 no :171
## 1st Qu.: 0.000 male :286 1st Qu.:27.00 1st Qu.: 4.000 yes:430
## Median : 0.000 Median :32.00 Median : 7.000
## Mean : 1.456 Mean :32.49 Mean : 8.178
## 3rd Qu.: 0.000 3rd Qu.:37.00 3rd Qu.:15.000
## Max. :12.000 Max. :57.00 Max. :15.000
## religiousness education occupation rating
## Min. :1.000 Min. : 9.00 Min. :1.000 Min. :1.000
## 1st Qu.:2.000 1st Qu.:14.00 1st Qu.:3.000 1st Qu.:3.000
## Median :3.000 Median :16.00 Median :5.000 Median :4.000
## Mean :3.116 Mean :16.17 Mean :4.195 Mean :3.932
## 3rd Qu.:4.000 3rd Qu.:18.00 3rd Qu.:6.000 3rd Qu.:5.000
## Max. :5.000 Max. :20.00 Max. :7.000 Max. :5.000
table(Affairs$affairs)
##
## 0 1 2 3 7 12
## 451 34 17 19 42 38
# 创建二分类因变量，1表示发生婚外情，0表示不发生婚外情
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)
##
## No Yes
## 451 150
# 拟合模型(link = “logit”)
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.5713 -0.7499 -0.5690 -0.2539 2.5191
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.37726 0.88776 1.551 0.120807
## gendermale 0.28029 0.23909 1.172 0.241083
## age -0.04426 0.01825 -2.425 0.015301 *
## yearsmarried 0.09477 0.03221 2.942 0.003262 **
## childrenyes 0.39767 0.29151 1.364 0.172508
## religiousness -0.32472 0.08975 -3.618 0.000297 ***
## education 0.02105 0.05051 0.417 0.676851
## occupation 0.03092 0.07178 0.431 0.666630
## rating -0.46845 0.09091 -5.153 2.56e-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 freedom
## Residual deviance: 609.51 on 592 degrees of freedom
## AIC: 627.51
##
## Number of Fisher Scoring iterations: 4
从结果中，我们发现性别gendermale、是否有孩子childrenyes、教育水平education和职业occupation对于模型贡献不显著，因此去除这些变量重新拟合模型。
# 重新拟合模型
fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data = Affairs, family = binomial())
summary(fit.reduced)
##
## Call:
## glm(formula = ynaffair ~ age + yearsmarried + religiousness +
## rating, family = binomial(), data = Affairs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6278 -0.7550 -0.5701 -0.2624 2.3998
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.93083 0.61032 3.164 0.001558 **
## age -0.03527 0.01736 -2.032 0.042127 *
## yearsmarried 0.10062 0.02921 3.445 0.000571 ***
## religiousness -0.32902 0.08945 -3.678 0.000235 ***
## rating -0.46136 0.08884 -5.193 2.06e-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 freedom
## Residual deviance: 615.36 on 596 degrees of freedom
## AIC: 625.36
##
## Number of Fisher Scoring iterations: 4
此时，每一个变量对模型贡献都非常显著，我们可以使用卡方检验比较两个模型。
# 比较模型
anova(fit.reduced, fit.full, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: ynaffair ~ age + yearsmarried + religiousness + rating
## Model 2: ynaffair ~ gender + age + yearsmarried + children + religiousness +
## education + occupation + rating
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 596 615.36
## 2 592 609.51 4 5.8474 0.2108
发现最终p=0.21，表明两个模型没有显著差异，因此用age + yearsmarried + religiousness + rating四个变量就能很好地预测是否发生婚外情。
# 输出回归系数，解释回归系数
coef(fit.reduced)
## (Intercept) age yearsmarried religiousness rating
## 1.93083017 -0.03527112 0.10062274 -0.32902386 -0.46136144
exp(coef(fit.reduced))
## (Intercept) age yearsmarried religiousness rating
## 6.8952321 0.9653437 1.1058594 0.7196258 0.6304248
#置信区间
exp (confint(fit.reduced))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 2.1255764 23.3506030
## age 0.9323342 0.9981470
## yearsmarried 1.0448584 1.1718250
## religiousness 0.6026782 0.8562807
## rating 0.5286586 0.7493370
回归系数可用coef()获取，exp(回归系数)可以求得自变量引起因变量变化的优势比OR，对于发病率较低的慢性疾病，OR可作为相对危险度RR的估计。
OR=1，自变量X对于应变量发生与否不起作用，OR>1是一个危险因素，OR<1是一个保护因素。（当1表示发生，0表示不发生的情况）
此处结婚年龄yearsmarried是发生婚外情的危险因素，age、religiousness、rating是保护因素。
评价预测变量对结果概率的影响
我们可以假定其他因素不变，仅改变其中一个变量，从而评价这个变量对于结果概率的影响。
# 评价婚姻评分rating
testdata <- data.frame(rating = c(1, 2, 3, 4, 5), age = mean(Affairs$age), yearsmarried = mean(Affairs$yearsmarried), religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata = testdata, type = "response")
testdata
## rating age yearsmarried religiousness prob
## 1 1 32.48752 8.177696 3.116473 0.5302296
## 2 2 32.48752 8.177696 3.116473 0.4157377
## 3 3 32.48752 8.177696 3.116473 0.3096712
## 4 4 32.48752 8.177696 3.116473 0.2204547
## 5 5 32.48752 8.177696 3.116473 0.1513079
婚姻评分从1到5，婚外情概率从0.53降到0.15
# 评价年龄age
testdata <- data.frame(rating = mean(Affairs$rating), age = seq(17, 57, 10), yearsmarried = mean(Affairs$yearsmarried), religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata = testdata, type = "response")
testdata
## rating age yearsmarried religiousness prob
## 1 3.93178 17 8.177696 3.116473 0.3350834
## 2 3.93178 27 8.177696 3.116473 0.2615373
## 3 3.93178 37 8.177696 3.116473 0.1992953
## 4 3.93178 47 8.177696 3.116473 0.1488796
## 5 3.93178 57 8.177696 3.116473 0.1094738
年龄从17到57，婚外情概率从0.34降到0.11
过度离势
过度离势是指观测到的响应变量的方差大于期望的二项分布的方差，过度离势会导致奇异的标准误检验和不精确的显著性检验。
deviance(fit.reduced)/df.residual(fit.reduced)
## [1] 1.03248
结果非常接近1，表示没有过度离势。
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)
## [1] 0.340122
p=0.34，二者之间没有显著差异，表明没有过度离势。
如果存在过度离势，可使用类二项分布 family = quasibinomial()。
条件logistic回归
使用survival包中的clogit()，用于分析配对数据
library(survival)
## Warning: package 'survival' was built under R version 3.6.3
data(logan)
summary(logan)
## occupation focc education race
## farm : 19 farm : 92 Min. : 2.00 non-black:764
## operatives :217 operatives :235 1st Qu.:12.00 black : 74
## craftsmen :202 craftsmen :232 Median :13.00
## sales :105 sales : 82 Mean :13.58
## professional:295 professional:197 3rd Qu.:16.00
## Max. :20.00
#整理数据
resp <- levels(logan$occupation)
n <- nrow(logan)
indx <- rep(1:n, length(resp))
logan2 <- data.frame(logan[indx,],
id = indx,
tocc = factor(rep(resp, each=n)))
logan2$case <- (logan2$occupation == logan2$tocc)
# strata(id)表示配对样本的编号，其余与之前一致，不过此处分析了交互作用
model <-clogit(case ~ tocc + tocc:education + strata(id), logan2)
summary(model)
## Call:
## coxph(formula = Surv(rep(1, 4190L), case) ~ tocc + tocc:education +
## strata(id), data = logan2, method = "exact")
##
## n= 4190, number of events= 838
##
## coef exp(coef) se(coef) z Pr(>|z|)
## toccfarm -1.8964629 0.1500986 1.3807822 -1.373 0.16961
## toccoperatives 1.1667502 3.2115388 0.5656465 2.063 0.03914
## toccprofessional -8.1005492 0.0003034 0.6987244 -11.593 < 2e-16
## toccsales -5.0292297 0.0065438 0.7700862 -6.531 6.54e-11
## tocccraftsmen:education -0.3322842 0.7172835 0.0568682 -5.843 5.13e-09
## toccfarm:education -0.3702858 0.6905370 0.1164100 -3.181 0.00147
## toccoperatives:education -0.4222188 0.6555906 0.0584328 -7.226 4.98e-13
## toccprofessional:education 0.2782469 1.3208122 0.0510212 5.454 4.94e-08
## toccsales:education NA NA 0.0000000 NA NA
##
## toccfarm
## toccoperatives *
## toccprofessional ***
## toccsales ***
## tocccraftsmen:education ***
## toccfarm:education **
## toccoperatives:education ***
## toccprofessional:education ***
## toccsales:education
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## toccfarm 0.1500986 6.6623 1.002e-02 2.247505
## toccoperatives 3.2115388 0.3114 1.060e+00 9.731781
## toccprofessional 0.0003034 3296.2778 7.713e-05 0.001193
## toccsales 0.0065438 152.8152 1.447e-03 0.029603
## tocccraftsmen:education 0.7172835 1.3941 6.416e-01 0.801857
## toccfarm:education 0.6905370 1.4481 5.497e-01 0.867512
## toccoperatives:education 0.6555906 1.5253 5.846e-01 0.735141
## toccprofessional:education 1.3208122 0.7571 1.195e+00 1.459723
## toccsales:education NA NA NA NA
##
## Concordance= 0.766 (se = 0.012 )
## Likelihood ratio test= 665.5 on 8 df, p=<2e-16
## Wald test = 413.5 on 8 df, p=<2e-16
## Score (logrank) test = 682.1 on 8 df, p=<2e-16
另外，有时我们还需要分析交互作用，使用逐步回归法step()等，之前推文均已讲过，此处不再赘述。
还要注意年龄等连续变量每增加一个变量对于二分类结果影响不大，经常会分组为有序多分类变量。有序多分类变量按照各个分类与因变量是否线性变化决定是否哑变量化。无序多分类自变量需要哑变量化。R中使用factor函数即可实现哑变量化。
好了，今天的R语言实现统计方法系列推文暂时告一段落，我们下次再见吧！小伙伴们如果有什么统计上的问题，或者如果想要学习什么方面的生物信息内容，可以在微信群或者知识星球提问，没准哪天的推文就是专门解答你的问题哦！