R語言統(tǒng)計系列第11篇-Logistic回歸

今天是各類統(tǒng)計方法R語言實現的第11期，我們主要介紹Logistic回歸。Logistic回歸屬于廣義線性回歸，因此我們從廣義線性回歸講起。

廣義線性回歸

線性回歸模型要求因變量服從正態(tài)分布，但是當結果變量是分類型（有無，患病與否等，二分類常用Logistic回歸）、計數型（某地區(qū)某年發(fā)生腫瘤患者的人數等，常用泊松回歸）或者臨床上經常使用的無復發(fā)生存期數據等（常用Cox回歸），因變量不符合正態(tài)分布，無法直接使用線性回歸。而廣義線性模型擴展了線性模型的框架，可以進行非正態(tài)因變量的分析，在R語言中可以通過glm()函數實現。

glm()函數的參數

分布族 默認的連接函數
binomial (link = “l(fā)ogit”)
gaussian (link = “identity”)
gamma (link = “inverse”)
inverse.gaussian (link = “1/mu^2”)
poisson (link = “l(fā)og”)
quasi (link = “identity”, variance = “constant”)
quasibinomial (link = “l(fā)ogit”)
quasipoisson (link = “l(fā)og”)

連用的函數

函數 描述
summary() 展示擬合模型的細節(jié)
coefficients(), coef() 列出擬合模型的參數（截距項和斜率）
confint() 給出模型參數的置信區(qū)間（默認為95%）
residuals() 列出擬合模型的殘差值
anova() 生成兩個擬合模型的方差分析表
plot() 生成評價擬合模型的診斷圖
predict() 用擬合模型對新數據集進行預測

Logistic回歸

二分類因變量常用Logistic回歸，假設因變量Y服從二項分布，查表得(link = “l(fā)ogit”)

此處是一份婚外情數據，我們用性別、年齡等因素預測參與者是否發(fā)生婚外情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
# 創(chuàng)建二分類因變量，1表示發(fā)生婚外情，0表示不發(fā)生婚外情
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 = “l(fā)ogit”) 
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

從結果中，我們發(fā)現性別gendermale、是否有孩子childrenyes、教育水平education和職業(yè)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

發(fā)現最終p=0.21，表明兩個模型沒有顯著差異，因此用age + yearsmarried + religiousness + rating四個變量就能很好地預測是否發(fā)生婚外情。

# 輸出回歸系數，解釋回歸系數
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
#置信區(qū)間
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(回歸系數)可以求得自變量引起因變量變化的優(yōu)勢比OR，對于發(fā)病率較低的慢性疾病，OR可作為相對危險度RR的估計。

OR=1，自變量X對于應變量發(fā)生與否不起作用，OR>1是一個危險因素，OR<1是一個保護因素。（當1表示發(fā)生，0表示不發(fā)生的情況）

此處結婚年齡yearsmarried是發(fā)生婚外情的危險因素，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()等，之前推文均已講過，此處不再贅述。

還要注意年齡等連續(xù)變量每增加一個變量對于二分類結果影響不大，經常會分組為有序多分類變量。有序多分類變量按照各個分類與因變量是否線性變化決定是否啞變量化。無序多分類自變量需要啞變量化。R中使用factor函數即可實現啞變量化。

好了，今天的R語言實現統(tǒng)計方法系列推文暫時告一段落，我們下次再見吧！ 小伙伴們如果有什么統(tǒng)計上的問題，或者如果想要學習什么方面的生物信息內容，可以在微信群或者知識星球提問，沒準哪天的推文就是專門解答你的問題哦！

?著作權歸作者所有,轉載或內容合作請聯系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。
禁止轉載，如需轉載請通過簡信或評論聯系作者。

友情鏈接更多精彩內容