在上一篇文章裡，我們給大家介紹了之前系列裡提及的線性回歸的擴充部分。
但醫學工作者最常接觸的結局預測變量多為二分類變量，比如陽性、陰性，病例、對照乃至生存、死亡這樣的變量。這樣我們就可以描述或推測在某些不同狀況下得某種疾病的風險或者說陽性時間發生的機率。這裡自然而然就引入我們今天的主題：
邏輯回歸模型——logistic regression model。關于邏輯回歸模型，需要注意的是，他與線性模型不同，沒有誤差項。我們是對一個事件發生的機率直接模組化，而二進制輸出的變異性将由此機率來确定。是以，與正态分布不同，這裡沒有方差這個參數。
我們這裡按照資料的原始類型分類來講解不同的原始資料應該怎樣通過R語言建立邏輯回歸模型。
A. 表格化資料的邏輯回歸這一個部分裡，我們采用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變量，重新進行模型的建立。
B. 原始資料的邏輯回歸我們同樣采用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
![r語言 回歸分析 分類變量_R語言 | 回歸分析（四）[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)