快放假了，這應(yīng)該是年前最后一篇文章，提前祝大家新年快樂！
接著上次的課程，本周刷完了邏輯回歸課程。選取一則醫(yī)療相關(guān)的案例進(jìn)行分析，通過醫(yī)療保險數(shù)據(jù)來評估醫(yī)療護(hù)理的好壞。有興趣的可以在公眾號后臺回復(fù)“quality”獲取數(shù)據(jù)下載鏈接，一起參與分析。
還是先加載數(shù)據(jù)，看看數(shù)據(jù)的基本特征。
> quality<-read.csv("quality.csv")
> str(quality)
可以看出有14個變量和131名觀測者。第一個變量MemberID表示觀測者編號，從InpatientDays到AcuteDrugGapSmall均為自變量，分別為：InpatientDays（住院天數(shù)）、ERVisits（病人進(jìn)急診室的次數(shù)）、OfficeVisits（診室看病的病人數(shù)量）、Narcotics（病人使用麻醉藥品處方的數(shù)量）、DaysSinceLastERVisit（病人上一次到急診室的時間至該研究結(jié)束的時間）、Pain（病人中抱怨疼痛的數(shù)量）、TotalVisits（病人來訪任何醫(yī)護(hù)人員的總次數(shù)）、ProviderCount（服務(wù)病人的服務(wù)者數(shù)量）、MedicalClaims（病人申請醫(yī)療理賠的天數(shù)）、ClaimLines（病人醫(yī)療理賠的所有次數(shù)）、StartedOnCombination（病人是否開始聯(lián)合使用藥物來治療糖尿?。?、AcuteDrugGapSmall（處方用完后迅速補(bǔ)充的急性藥物的部分）。最后一個PoorCare為因變量（0表示病人受到了高質(zhì)量的護(hù)理-goodcare；1表示病人受到了低質(zhì)量的護(hù)理-PoorCare）。
我們使用OfficeVisits（來訪病人的數(shù)量）和Narcotics（麻醉藥品處方的數(shù)量）兩個變量對因變量進(jìn)行預(yù)測。
> table(quality$PoorCare)
0? 1
98 33
先看看baseline model，有33人是PoorCare，98人是goodcare。所以，baseline model的準(zhǔn)確率為98/(98+33)=0.748，用來和邏輯回歸的模型進(jìn)行對比。
由于這節(jié)課只給到一個數(shù)據(jù)集，故需要使用random.split（）函數(shù)將數(shù)據(jù)分成training set 和testing set兩部分。注意：此處需要安裝caTools包。
> install.packages(“caTools”)
> library(caTools)
而random.split（ ）會將數(shù)據(jù)進(jìn)行隨機(jī)分配，為了便于教學(xué)，讓所有的分析者得到同樣的數(shù)據(jù)，需要設(shè)置相同的seed，即所有的分析者在set.seed（ ）函數(shù)中選取一個相同的數(shù)字（課程中填了88）。
> set.seed(88)
> split<-sample.split(quality$PoorCare,SplitRatio = 0.75)# 比例的設(shè)置范圍一般在0.5-0.8，數(shù)據(jù)量多，則可以把Training set少分配，而Testing set多分配，以增加預(yù)測信度。
> split
返回TRUE表示抽取觀測值放在Training set中，F(xiàn)ALSE表示抽取觀測值放在Testing set中。故需要創(chuàng)建兩個數(shù)據(jù)集，然后建立回歸模型。
> qualityTrain <- subset(quality, split == TRUE)
> qualityTest <- subset(quality, split == FALSE)
> qualityLog<-glm(PoorCare~OfficeVisits+Narcotics, data=qualityTrain,
family=binomial)
> summary(qualityLog)
上圖的模型結(jié)果可以看出，和線性回歸很相似。兩個變量前的系數(shù)均為正數(shù)，且變量顯著，表明這兩個變量增加會對因變量產(chǎn)生正的效果，即更好的預(yù)測接收到PoorCare的情況。除此之外，最下面是AIC，主要用來將模型進(jìn)行對比，從而檢測模型的擬合度。
> predictTrain<-predict(qualityLog,type = "response")#告知預(yù)測函數(shù)給出概率
> summary(predictTrain)
> tapply(predictTrain, qualityTrain$PoorCare, mean)#對預(yù)測值按照訓(xùn)練數(shù)據(jù)集PoorCare的類別進(jìn)行分類計算平均數(shù)
??????? 0???????? 1
0.1894512 ???0.4392246
對預(yù)測值按照訓(xùn)練數(shù)據(jù)集PoorCare的類別進(jìn)行分類，并計算不同類別的概率平均數(shù)（1表示訓(xùn)練集是PoorCare，預(yù)測值也是PoorCare；0表示訓(xùn)練集是GoodCare，預(yù)測值是PoorCare）。預(yù)測PoorCare的概率要大于實際training set中PoorCare概率的值，表示是好事，可以提前進(jìn)行干預(yù)。
邏輯回歸的結(jié)果是概率值，一般情況下，我們想要二分變量的預(yù)測，比如這個病人接受到高質(zhì)量的護(hù)理，還是低質(zhì)量的護(hù)理，這時需要用到閾值t：
如果P(PoorCare)>t，預(yù)測為poor quality care；
如果P(PoorCare)
那么如何選取t值呢？
t值的選取和一類錯誤、二類錯誤有關(guān)，即不同的t值會對sensitivity和specificity造成不同的影響，所以選取t值需要根據(jù)具體情況而定。比如在這個案例中，主要目的是想提高對PoorCare的預(yù)測，那么可以適當(dāng)降低t值，獲得較高的sensitivity（True positive rate）。反之，如果想提高對goodcare的預(yù)測，可以適當(dāng)提高t值，獲得較高的specificity（True negative rate）。
（此處實際上就是信號檢測論的基本原理，如果不是很明白，強(qiáng)烈建議觀看該鏈接中的視頻：https://www.youtube.com/watch?v=vtYDyGGeQyo）
> table(qualityTrain$PoorCare,predictTrain>0.2)#第二個參數(shù)表示預(yù)測值大于閾值0.2，則返回TRUE，反之，返回FALSE
?? FALSE TRUE
? 0??? 54?? 20
? 1???? 9?? 16
因此，sensitivity1=16/(16+9)=0.64，specificity1=54/(54+20)=0.73
> table(qualityTrain$PoorCare,predictTrain>0.5)
??? FALSE TRUE
? 0??? 70??? 4
? 1??? 15?? 10
sensitivity2=10/(10+15)=0.4，specificity2=70/(70+4)=0.99
我們可以看到，隨著閾值t的增加，sensitivity會相應(yīng)降低，而specificity會相應(yīng)提高。
那么，閾值t要選取多少呢？ 可以用ROC（接受者操作特性曲線）來綜合判斷。
> install.packages(“ROCR”)
> library(ROCR)
> ROCRperd<-prediction(predictTrain,qualityTrain$PoorCare)
> ROCRperf<-performance(ROCRperd,"tpr", "fpr")
>?plot(ROCRperf,colorize=TRUE,print.cutoffs.at=seq(0,1,by=0.1),text.adj=c(-0.2,0.5)）
如果想要較低的虛報率（圖中橫軸，表示1-specificity），那么就最大化命中率（圖中的縱軸，表示sensitivity）的同時，保持較低的虛報率，比如使用（0.1,0.5）的閾值，接近0.3。反之亦然。
假如我們使用閾值0.3，此時對testing set進(jìn)行預(yù)測。
> predictTest<-predict(QualityLog,newdata = qualityTest,type = "response")
> table(qualityTest$PoorCare,predictTest>0.3)
? FALSE TRUE
? 0??? 19??? 5
? 1???? 2??? 6
> as.numeric(performance(ROCRperd,"auc")@y.values)?#計算AUC
[1] 0.7745946
計算出模型預(yù)測testing set樣本的準(zhǔn)確率為（19+6）/（19+5+2+6）= 0.78，預(yù)測的質(zhì)量AUC=0.77。該模型可以識別大多數(shù)接受低質(zhì)量護(hù)理的病人。
上一篇學(xué)的是線性模型，可以通過一系列連續(xù)變量來預(yù)測正態(tài)分布的結(jié)果變量。但實際情況中，很多結(jié)果變量并不是正態(tài)分布，比如：通過/失敗，活著/死亡。這就需要用邏輯回歸來做預(yù)測。
說實話，這個unit我看到更多的是理論，對于實際來說，邏輯回歸的模型并不是很容易解釋，比如模型的兩個自變量相關(guān)系數(shù)都為正數(shù)，但如何解釋對PoorCare和goodcare的影響？這難道就是“聽過很多道理，依然過不好這一生”？
簡單看了下Unit3，表示Trees可以更好的解釋模型，那就拭目以待吧。
相關(guān)文章
