R數(shù)據(jù)分析案例：邏輯回歸

邏輯回歸，也稱之為邏輯模型，用于預(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

記?。罕疚囊庠谝阒廊绾斡孟嚓P(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("http://www.ats.ucla.edu/stat/data/binary.csv")##view the first few rows of the datahead(mydata)
2. ##admit gre gpa rank
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)?，所以它們統(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ì)量），我們可以使用下面的指令：文章來源cda數(shù)據(jù)分析培訓(xùn)

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)
3. ## 'log Lik.' -229 (df=6)