文|何品言

邏輯回歸，也稱之為邏輯模型，用于預(yù)測(cè)二分結(jié)果變量。在邏輯模型當(dāng)中，輸出結(jié)果所占的比率就是預(yù)測(cè)變量的線性組合。

這篇文章將要使用下面這幾個(gè)包，而且你們需要保證在運(yùn)行我所舉的例子的時(shí)候，你已經(jīng)把這些包都裝好了。如果你還沒裝好這些包，那么，運(yùn)行install.packages(“R包名稱”)這個(gè)操作，或者你可能需要更新版本，運(yùn)行update.packages()。

library(aod)

library(ggplot2)

library(Rcpp)

版本說明：本文基于R3.2.3進(jìn)行測(cè)試的，一下是包的版本：

Rcpp：0.12.3 ggplot2：2.0.0 knitr：1.12.3

記住：本文意在要你知道如何用相關(guān)的指令進(jìn)行邏輯回歸分析，而并沒有涵蓋所有研究人員可能會(huì)做的事情，尤其是數(shù)據(jù)是沒有進(jìn)行清洗和查閱的，而且假設(shè)并非最嚴(yán)謹(jǐn)，其它方面也不會(huì)相當(dāng)?shù)臉?biāo)準(zhǔn)。

案例

案例1： 假如我們對(duì)這些因素感興趣，它們表示政治候選人是否贏得選舉的因子，其中，我們把結(jié)果變量表示為0或1，也可以表達(dá)成贏或者輸。而預(yù)測(cè)變量的利益可由一場(chǎng)運(yùn)動(dòng)中所投入的金額表示，或者是選舉人所花的時(shí)間，再或者在職人員是否獲得足夠的支持。

案例2：一個(gè)研究者可能會(huì)對(duì)GRE（研究生入學(xué)考試成績(jī)）、GPA（大學(xué)平均績(jī)點(diǎn)），以及研究生學(xué)院的名譽(yù)感興趣，因?yàn)樗鼈冇绊憣W(xué)校的招生問題。這里，我們用允許/不允許這個(gè)二進(jìn)制結(jié)果表示其因變量。

數(shù)據(jù)的描述

對(duì)于我們接下來要進(jìn)行的數(shù)據(jù)分析來說，我們要對(duì)案例2的入學(xué)問題進(jìn)行深入的探討。我們有了通常情況下假設(shè)所產(chǎn)生的數(shù)據(jù)，而它們可從R的相關(guān)網(wǎng)站得到。記住，當(dāng)我們要調(diào)用相關(guān)函數(shù)導(dǎo)入數(shù)據(jù)的時(shí)候，如果我們要具體表示一個(gè)硬盤驅(qū)動(dòng)器的文件，我們需要打上斜杠（/），而不是反斜杠（\）。

1.mydata <- 
read.csv("//www.ats.ucla.edu/stat/data/binary.csv")## view the first few rows of the datahead(mydata)
2. ##admit gre gpa rank  
3.
4.##1 0 380 3.61 3 
5.## 2 1 660 3.67 3 
6.## 3 1 800 4.00 1 
7.## 4 1 640 3.19 4 
8.## 5 0 520 2.93 4 
9.## 6 1 760 3.00 2 

這個(gè)數(shù)據(jù)集有二進(jìn)制的結(jié)果（輸出值，依賴），它表示允許。這里有3個(gè)預(yù)測(cè)變量：gre、gpa以及rank。我們把gre和gpa看作是連續(xù)變量。rank表示有4個(gè)值為1。這里，為0的那所學(xué)校聲望最高，其它的這4所高校聲望最低。這時(shí)，我們可以用summary()函數(shù)來匯總一下這個(gè)數(shù)據(jù)集的情況，而且，如果想要計(jì)算里面的標(biāo)準(zhǔn)差，我們可以使用sapply()函數(shù)，并在里面寫上sd來獲取其標(biāo)準(zhǔn)差。

1.summary(mydata) 

2.## admit gre gpa rank 

3.## Min. :0.000 Min. :220 Min. :2.26 Min. :1.00 

4.## 1st Qu.:0.000 1st Qu.:520 1st Qu.:3.13 1st Qu.:2.00 

5.## Median :0.000 Median :580 Median :3.40 Median :2.00 

6.## Mean :0.318 Mean :588 Mean :3.39 Mean :2.48 

7.## 3rd Qu.:1.000 3rd Qu.:660 3rd Qu.:3.67 3rd Qu.:3.00 

8.## Max. :1.000 Max. :800 Max. :4.00 Max. :4.00 

9.sapply(mydata, sd) 

10.## admit gre gpa rank 

11.## 0.466 115.517 0.381 0.944 

12.## two-way contingency table of categorical outcome and predictors## we want to make sure there are not 0 cellsxtabs(~ admit + rank, data = mydata) 

13.## rank 

14.## admit 1 2 3 4 

15.## 0 28 97 93 55 

16.## 1 33 54 28 12 

你可能會(huì)考慮到的分析方法

接下來，我會(huì)列舉一些你可能會(huì)用到的方法。這里所列舉的一些方法相對(duì)來說比較合理，畢竟有些其他方法不能執(zhí)行或者尤其局限性。

1.邏輯回歸，也是本文的重點(diǎn)

2.概率回歸。概率分析所產(chǎn)生的結(jié)果類似于邏輯回歸。我們可以依據(jù)我們的需要進(jìn)行有選擇性的進(jìn)行概率回歸或邏輯回歸。

3.最小二乘法。當(dāng)你使用二進(jìn)制輸出變量的時(shí)候，這種模型就是常用于線性回歸分析的，而且也可以用它來進(jìn)行條件概率的運(yùn)算。然而，這里的誤差（例如殘差）來自于線性概率模型，而且會(huì)影響到異方差性最小二乘法回歸的正態(tài)誤差檢驗(yàn)，它也影響無效標(biāo)準(zhǔn)誤差和假設(shè)性分析。想要了解更多關(guān)于線性概率模型的相關(guān)問題，可以查閱Long（1997，p.38-40）。

4.雙組判別函數(shù)分析。這是一個(gè)多元的方法處理輸出變量。

5.Hotelling 的T2。0/1的輸出結(jié)果轉(zhuǎn)換到分組變量中，而前一個(gè)預(yù)測(cè)值則成為輸出變量。這樣可以進(jìn)行一個(gè)顯著性測(cè)試，然而并沒有給沒個(gè)變量分別給一個(gè)系數(shù)，而且這對(duì)于某個(gè)變量根據(jù)影響情況調(diào)整到其它變量的程度是不清楚的。

使用邏輯模型

接下來的代碼，通過使用glm()函數(shù)（廣義線性模型）進(jìn)行相關(guān)評(píng)估。首先，我們要把rank（秩）轉(zhuǎn)換成因子，并預(yù)示著rank在這里被視為分類變量。

1.mydata$rank <- factor(mydata$rank)mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial") 

由于我們得到了模型的名字（mylogit)，而R并不會(huì)從我們的回歸中產(chǎn)生任何輸出結(jié)果。為了得到結(jié)果，我們使用summary()函數(shù)進(jìn)行提取：

1.summary(mylogit) 

2.## 

3.## Call: 

4.## glm(formula = admit ~ gre + gpa + rank, family = "binomial", 

5.## data = mydata) 

6.## 

7.## Deviance Residuals: 

8.## Min 1Q Median 3Q Max 

9.## -1.627 -0.866 -0.639 1.149 2.079 

10.## 

11.## Coefficients: 

12.## Estimate Std. Error z value Pr(>|z|) 

13.## (Intercept) -3.98998 1.13995 -3.50 0.00047 *** 

14.## gre 0.00226 0.00109 2.07 0.03847 * 

15.## gpa 0.80404 0.33182 2.42 0.01539 * 

16.## rank2 -0.67544 0.31649 -2.13 0.03283 * 

17.## rank3 -1.34020 0.34531 -3.88 0.00010 *** 

18.## rank4 -1.55146 0.41783 -3.71 0.00020 *** 

19.## --- 

20.## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

21.## 

22.## (Dispersion parameter for binomial family taken to be 1) 

23.## 

24.## Null deviance: 499.98 on 399 degrees of freedom 

25.## Residual deviance: 458.52 on 394 degrees of freedom 

26.## AIC: 470.5 

27.## 

28.## Number of Fisher Scoring iterations: 4

1.在上面的結(jié)果中，我們首先看到的就是call，這提示我們這時(shí)的R在運(yùn)行什么東西，我們?cè)O(shè)定了什么選項(xiàng)，等等。

2.接下來，我們看到了偏差殘差，用于測(cè)量模型的擬合度。這部分的結(jié)果顯示了這個(gè)分布的偏差殘差，而針對(duì)每個(gè)使用在模型里的個(gè)案。接下來，我們討論一下如何匯總偏差估計(jì)，從而知道模型的擬合結(jié)果。

3.下一部分的結(jié)果顯示了相關(guān)系數(shù)，標(biāo)準(zhǔn)誤差，z統(tǒng)計(jì)量（有時(shí)也叫Wald Z統(tǒng)計(jì)量），以及它的相關(guān)結(jié)果。gre和gpa都是同等重要的統(tǒng)計(jì)量，并作為rank的3個(gè)變量。邏輯回歸系數(shù)給改變了在預(yù)測(cè)變量中增加一個(gè)單位的輸出結(jié)果誤差。

對(duì)于gre，每改變一個(gè)單位，輸出結(jié)果的允許誤差（相對(duì)于不允許來說）增加0.002。

對(duì)于gpa，每改變一個(gè)單位，入學(xué)的研究生的允許誤差增加0.804。

ranK的指示變量的觀察方法就有點(diǎn)不同了。例如，層在本科學(xué)習(xí)過的在rank的值是2，其它為1，使之允許誤差變?yōu)?0.675。

1.4.下面的圖表里的系數(shù)表示不錯(cuò)的擬合效果，包含了空值、偏差殘差以及AIC。后面，我們就會(huì)學(xué)習(xí)如何用這些信息來判斷這個(gè)模型是否擬合。 

我們可以使用confint（）函數(shù)來獲取相關(guān)的區(qū)間的預(yù)測(cè)信息。對(duì)于邏輯回歸模型來說，其置信區(qū)間時(shí)基于異形對(duì)數(shù)似然函數(shù)求出來的。同樣，我們得到的CL是基于默認(rèn)方法求出的標(biāo)準(zhǔn)差算的。

1.## CIs using profiled log-likelihoodconfint(mylogit) 

2.## Waiting for profiling to be done... 

3.## 2.5 % 97.5 % 

4.## (Intercept) -6.271620 -1.79255 

5.## gre 0.000138 0.00444 

6.## gpa 0.160296 1.46414 

7.## rank2 -1.300889 -0.05675 

8.## rank3 -2.027671 -0.67037 

9.## rank4 -2.400027 -0.75354 

10.## CIs using standard errorsconfint.default(mylogit) 

11.## 2.5 % 97.5 % 

12.## (Intercept) -6.22424 -1.75572 

13.## gre 0.00012 0.00441 

14.## gpa 0.15368 1.45439 

15.## rank2 -1.29575 -0.05513 

16.## rank3 -2.01699 -0.66342 

17.## rank4 -2.37040 -0.73253 

我們可以調(diào)用aod包里的wald.test()函數(shù)來測(cè)試rank的所有影響。圖表里的系數(shù)的順序和模型里的項(xiàng)順序是一樣的。這樣很重要，因?yàn)閣ald.test()函數(shù)就是基于這些模型的項(xiàng)順序進(jìn)行測(cè)試的。b提供了系數(shù)，而Sigma提供了誤差項(xiàng)的方差協(xié)方差矩陣，而且最終，這些項(xiàng)告訴R哪些項(xiàng)用來進(jìn)行測(cè)試，而在這種情況下，第4、5、6這三項(xiàng)作為rank的層次進(jìn)行測(cè)試。

1.wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6) 

2.## Wald test: 

3.## ---------- 

4.## 

5.## Chi-squared test: 

6.## X2 = 20.9, df = 3, P(> X2) = 0.00011 

卡方檢驗(yàn)算出來的值是20.9，這里涉及到3個(gè)自由度，p值算出來是0.00011，這預(yù)示著我們所假設(shè)的這些項(xiàng)之間具有顯著的影響效果。

我們也可以測(cè)試而外的假設(shè)，這些假設(shè)包含rank里不同層次的差異。下面，我們通過測(cè)試得知rank=2的測(cè)試結(jié)果和rank=3的時(shí)候一樣。下面的第一行代碼創(chuàng)造的向量是1，這定義了測(cè)試?yán)镂覀円獔?zhí)行的內(nèi)容。在這種情況下，我們想要測(cè)試rank=2和rank=3（即模型的第4項(xiàng)、第5項(xiàng)）這兩個(gè)不同的項(xiàng)所產(chǎn)生的差異（減）。為了對(duì)比這2項(xiàng)，我們對(duì)其中一項(xiàng)乘以1，另一項(xiàng)乘以-1。其它項(xiàng)不包含在測(cè)試?yán)铮运鼈兘y(tǒng)一都乘上0。第二行代碼寫上L=1，這告訴R基于向量為1時(shí)執(zhí)行這次測(cè)試（并不是我們之前所選擇的進(jìn)行測(cè)試）。

1.l <- cbind="" 0="" 0="" 0="" 1="" -1="" 0="" wald="" test="" b="coef(mylogit)," sigma="vcov(mylogit)," l="l)" wald="" test:="" ----------="" chi-squared="" test:="" x2="5.5," df="1," p=""> X2) = 0.019 

卡方檢驗(yàn)所在自由度為1的情況下算出來的結(jié)果為5.5，并得出相關(guān)的p值為0.019，這預(yù)示著rank=2和rank=3之間存在著顯著的差異。

你一可以把系數(shù)指數(shù)化，并從誤差率進(jìn)行解讀，而R會(huì)自動(dòng)幫你算出來。為了要得到指數(shù)化系數(shù)，你可以告訴R你想要指數(shù)化（exp），而R也會(huì)按照你的要求把它們指數(shù)化，它屬于mylogit（coef（mylogit））的一部分。我們可以使用相同的誤差率以及它們的置信區(qū)間，并可先前就把置信區(qū)間指數(shù)化。把它們都放到其中一個(gè)圖表，我們可以使用cbind系數(shù)和系數(shù)置信區(qū)間這些列整合起來。

1.## odds ratios onlyexp(coef(mylogit)) 

2.## (Intercept) gre gpa rank2 rank3 rank4 

3.## 0.0185 1.0023 2.2345 0.5089 0.2618 0.2119 

4.## odds ratios and 95% CIexp(cbind(OR = coef(mylogit), confint(mylogit))) 

5.## Waiting for profiling to be done... 

6.## OR 2.5 % 97.5 % 

7.## (Intercept) 0.0185 0.00189 0.167 

8.## gre 1.0023 1.00014 1.004 

9.## gpa 2.2345 1.17386 4.324 

10.## rank2 0.5089 0.27229 0.945 

11.## rank3 0.2618 0.13164 0.512 

12.## rank4 0.2119 0.09072 0.471 

現(xiàn)在，我們可以說gpa增加了一個(gè)單位，而研究生入學(xué)（反之就是沒有入學(xué)的）的誤差則在因子上增加2.23。對(duì)于更多關(guān)于誤差率信息的解讀，查看我們的FAQ頁 How do I interpret odds ratios in logistic regression? 。注意，當(dāng)R產(chǎn)生了這個(gè)結(jié)果時(shí)，關(guān)于誤差率的攔截一般都不被解讀。

你也可以使用預(yù)測(cè)概率來幫助你解讀這個(gè)模型。預(yù)測(cè)概率可以均可由分類預(yù)測(cè)變量或連續(xù)預(yù)測(cè)變量計(jì)算出來。為了算出預(yù)測(cè)概率，我們首先要?jiǎng)?chuàng)建含有我們需要的獨(dú)立變量來創(chuàng)建新的數(shù)據(jù)框。

我們將要開始計(jì)算每個(gè)rank值的預(yù)測(cè)概率的允許值，并計(jì)算gre和gpa的平均值。首先，我們要?jiǎng)?chuàng)建新的數(shù)據(jù)框：

1.newdata1 <- with(mydata, 

2.data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4))) 

3.## view data framenewdata1 

4.## gre gpa rank 

5.## 1 588 3.39 1 

6.## 2 588 3.39 2 

7.## 3 588 3.39 3 

8.## 4 588 3.39 4 

這些值必須含有和你之前所創(chuàng)建的邏輯回歸分析相同的名字（例如，在這個(gè)例子中，gre的就必須命名為gre）。既然，我們現(xiàn)在已經(jīng)創(chuàng)建好了我們需要進(jìn)行運(yùn)算的數(shù)據(jù)框，那么我們可以告訴R根據(jù)這個(gè)來創(chuàng)建預(yù)測(cè)概率。第一行代碼經(jīng)過了壓縮，我們現(xiàn)在就把它分開來，討論這些值是怎樣執(zhí)行的。Newdata$rankP告訴R我們要根據(jù)數(shù)據(jù)集（數(shù)據(jù)框）的newdata1里的rankP創(chuàng)建一個(gè)新的變量剩余的指令告訴R這些rankP值應(yīng)當(dāng)使用prediction()函數(shù)進(jìn)行預(yù)測(cè)。圓括號(hào)里的選項(xiàng)告訴R這些預(yù)測(cè)值應(yīng)該基于分析mylogit進(jìn)行預(yù)測(cè)，mylogit的值源自newdata1以及它的預(yù)測(cè)值類型就是預(yù)測(cè)概率（type=”response”)。第二行代碼列舉了數(shù)據(jù)框newdata1的值，盡管它不是十分理想，而這就是圖表的預(yù)測(cè)概率。

1.newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response")newdata1 

2.## gre gpa rank rankP 

3.## 1 588 3.39 1 0.517 

4.## 2 588 3.39 2 0.352 

5.## 3 588 3.39 3 0.219 

6.## 4 588 3.39 4 0.185 

在上面的預(yù)測(cè)結(jié)果中，我們看到來自最好的名校（rank=1)并被接收到研究生的本科生預(yù)測(cè)概率是0.52，而0.18的學(xué)生來自最低檔次的學(xué)校（rank=4），以gre和gpa作為平均值。我們可以做相似的事情來創(chuàng)建一個(gè)針對(duì)不斷變化的變量gre和gpa的預(yù)測(cè)概率圖表。我們可以基于此作圖，所以我們可以在200到800之間創(chuàng)建100個(gè)gre值，基于它的rank（1,2,3,4）。

newdata2 <- with(mydata,

data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 4),

gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))

這些代碼所產(chǎn)生的預(yù)測(cè)概率（下面第一行）和之前算的一樣，除非我們還想要對(duì)標(biāo)準(zhǔn)差進(jìn)行要求，否則我們可以對(duì)置信區(qū)間作圖。我們可以對(duì)關(guān)聯(lián)規(guī)模進(jìn)行預(yù)測(cè)，同時(shí)反向變換預(yù)測(cè)值和置信區(qū)間的臨近值到概率中。

1.newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type="link", se=TRUE))newdata3 <- within(newdata3, { 

2.PredictedProb <- plogis(fit) 

3.LL <- plogis(fit - (1.96 * se.fit)) 

4.UL <- plogis(fit + (1.96 * se.fit))}) 

5.## view first few rows of final datasethead(newdata3) 

6.## gre gpa rank fit se.fit residual.scale UL LL PredictedProb 

7.## 1 200 3.39 1 -0.811 0.515 1 0.549 0.139 0.308 

8.## 2 206 3.39 1 -0.798 0.509 1 0.550 0.142 0.311 

9.## 3 212 3.39 1 -0.784 0.503 1 0.551 0.145 0.313 

10.## 4 218 3.39 1 -0.770 0.498 1 0.551 0.149 0.316 

11.## 5 224 3.39 1 -0.757 0.492 1 0.552 0.152 0.319 

12.## 6 230 3.39 1 -0.743 0.487 1 0.553 0.155 0.322 

當(dāng)然，使用圖像描繪預(yù)測(cè)概率來解讀和展示模型也是相當(dāng)有用的。我們會(huì)使用ggplot2包來作圖。下面我們作圖描繪預(yù)測(cè)概率，和95%置信區(qū)間。

1.ggplot(newdata3, aes(x = gre, y = PredictedProb)) + 

2.geom_ribbon(aes(ymin = LL, ymax = UL, fill = rank), alpha = .2) + 

3.geom_line(aes(colour = rank), size=1) 

我們也許很想看到這個(gè)模型的擬合效果怎么樣，而剛剛那樣做是非常有用的，尤其是對(duì)比這些競(jìng)爭(zhēng)的模型。Summary(mylogit)所產(chǎn)生的結(jié)果包含了擬合系數(shù)（下面展示了其系數(shù)），包含了空值、偏差殘差和AIC。模型擬合度的一個(gè)衡量標(biāo)準(zhǔn)就是整個(gè)模型的顯著程度。這個(gè)測(cè)試問了我們使用了預(yù)測(cè)值的模型是否比僅僅含有截距的模型（即，空模型）更加顯著，而這個(gè)測(cè)試?yán)锏慕y(tǒng)計(jì)量是含有預(yù)測(cè)值的預(yù)測(cè)擬合指數(shù)和空模型的差，而且這個(gè)模型的統(tǒng)計(jì)量也是含有其自由度等于現(xiàn)有的模型（即，模型里的預(yù)測(cè)變量個(gè)數(shù)）和空模型的自由度之差的卡方分布。為了要找到這兩個(gè)模型的偏差（即，這個(gè)測(cè)試的統(tǒng)計(jì)量），我們可以使用下面的指令：

1.with(mylogit, null.deviance - deviance) 

2.## [1] 41.5 

這兩個(gè)模型的自由度之差等于這個(gè)模型里預(yù)測(cè)變量的個(gè)數(shù)，而且我們可以按照下面的方法獲取它：

1.with(mylogit, df.null - df.residual) 

2.## [1] 5 

最后，我們提取一下p值：

1.with(mylogit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE)) 

2.## [1] 7.58e-08 

在自由度為5的情況下算出來的卡方是41.46，而相關(guān)的p值則小于0.001。這告訴我們，此模型的擬合效果比一個(gè)空模型所產(chǎn)生的擬合效果更加顯著。這個(gè)，我們有時(shí)稱它為似然比測(cè)試（偏差殘差為-2*log似然值）。要看到它的對(duì)數(shù)似然，我們可以這樣寫：

1.logLik(mylogit) 

2.

3.## 'log Lik.' -229 (df=6)