快放假了，这应该是年前最后一篇文章，提前祝大家新年快乐！
接着上次的课程，本周刷完了逻辑回归课程。选取一则医疗相关的案例进行分析，通过医疗保险数据来评估医疗护理的好坏。有兴趣的可以在公众号后台回复“quality”获取数据下载链接，一起参与分析。
还是先加载数据，看看数据的基本特征。
> quality<-read.csv("quality.csv")
> str(quality)
可以看出有14个变量和131名观测者。第一个变量MemberID表示观测者编号，从InpatientDays到AcuteDrugGapSmall均为自变量，分别为：InpatientDays（住院天数）、ERVisits（病人进急诊室的次数）、OfficeVisits（诊室看病的病人数量）、Narcotics（病人使用麻醉药品处方的数量）、DaysSinceLastERVisit（病人上一次到急诊室的时间至该研究结束的时间）、Pain（病人中抱怨疼痛的数量）、TotalVisits（病人来访任何医护人员的总次数）、ProviderCount（服务病人的服务者数量）、MedicalClaims（病人申请医疗理赔的天数）、ClaimLines（病人医疗理赔的所有次数）、StartedOnCombination（病人是否开始联合使用药物来治疗糖尿病）、AcuteDrugGapSmall（处方用完后迅速补充的急性药物的部分）。最后一个PoorCare为因变量（0表示病人受到了高质量的护理-goodcare；1表示病人受到了低质量的护理-PoorCare）。
我们使用OfficeVisits（来访病人的数量）和Narcotics（麻醉药品处方的数量）两个变量对因变量进行预测。
> table(quality$PoorCare)
0 1
98 33
先看看baseline model，有33人是PoorCare，98人是goodcare。所以，baseline model的准确率为98/(98+33)=0.748，用来和逻辑回归的模型进行对比。
由于这节课只给到一个数据集，故需要使用random.split（）函数将数据分成training set 和testing set两部分。注意：此处需要安装caTools包。
> install.packages(“caTools”)
> library(caTools)
而random.split（ ）会将数据进行随机分配，为了便于教学，让所有的分析者得到同样的数据，需要设置相同的seed，即所有的分析者在set.seed（ ）函数中选取一个相同的数字（课程中填了88）。
> set.seed(88)
> split<-sample.split(quality$PoorCare,SplitRatio = 0.75)# 比例的设置范围一般在0.5-0.8，数据量多，则可以把Training set少分配，而Testing set多分配，以增加预测信度。
> split
返回TRUE表示抽取观测值放在Training set中，FALSE表示抽取观测值放在Testing set中。故需要创建两个数据集，然后建立回归模型。
> qualityTrain <- subset(quality, split == TRUE)
> qualityTest <- subset(quality, split == FALSE)
> qualityLog<-glm(PoorCare~OfficeVisits+Narcotics, data=qualityTrain,
family=binomial)
> summary(qualityLog)
上图的模型结果可以看出，和线性回归很相似。两个变量前的系数均为正数，且变量显著，表明这两个变量增加会对因变量产生正的效果，即更好的预测接收到PoorCare的情况。除此之外，最下面是AIC，主要用来将模型进行对比，从而检测模型的拟合度。
> predictTrain<-predict(qualityLog,type = "response")#告知预测函数给出概率
> summary(predictTrain)
> tapply(predictTrain, qualityTrain$PoorCare, mean)#对预测值按照训练数据集PoorCare的类别进行分类计算平均数
0 1
0.1894512 0.4392246
对预测值按照训练数据集PoorCare的类别进行分类，并计算不同类别的概率平均数（1表示训练集是PoorCare，预测值也是PoorCare；0表示训练集是GoodCare，预测值是PoorCare）。预测PoorCare的概率要大于实际training set中PoorCare概率的值，表示是好事，可以提前进行干预。
逻辑回归的结果是概率值，一般情况下，我们想要二分变量的预测，比如这个病人接受到高质量的护理，还是低质量的护理，这时需要用到阈值t：
如果P(PoorCare)>t，预测为poor quality care；
如果P(PoorCare)
那么如何选取t值呢？
t值的选取和一类错误、二类错误有关，即不同的t值会对sensitivity和specificity造成不同的影响，所以选取t值需要根据具体情况而定。比如在这个案例中，主要目的是想提高对PoorCare的预测，那么可以适当降低t值，获得较高的sensitivity（True positive rate）。反之，如果想提高对goodcare的预测，可以适当提高t值，获得较高的specificity（True negative rate）。
（此处实际上就是信号检测论的基本原理，如果不是很明白，强烈建议观看该链接中的视频：https://www.youtube.com/watch?v=vtYDyGGeQyo）
> table(qualityTrain$PoorCare,predictTrain>0.2)#第二个参数表示预测值大于阈值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会相应降低，而specificity会相应提高。
那么，阈值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进行预测。
> 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
计算出模型预测testing set样本的准确率为（19+6）/（19+5+2+6）= 0.78，预测的质量AUC=0.77。该模型可以识别大多数接受低质量护理的病人。
上一篇学的是线性模型，可以通过一系列连续变量来预测正态分布的结果变量。但实际情况中，很多结果变量并不是正态分布，比如：通过/失败，活着/死亡。这就需要用逻辑回归来做预测。
说实话，这个unit我看到更多的是理论，对于实际来说，逻辑回归的模型并不是很容易解释，比如模型的两个自变量相关系数都为正数，但如何解释对PoorCare和goodcare的影响？这难道就是“听过很多道理，依然过不好这一生”？
简单看了下Unit3，表示Trees可以更好的解释模型，那就拭目以待吧。
相关文章
