R数据分析案例：逻辑回归

逻辑回归，也称之为逻辑模型，用于预测二分结果变量。在逻辑模型当中，输出结果所占的比率就是预测变量的线性组合。

这篇文章将要使用下面这几个包，而且你们需要保证在运行我所举的例子的时候，你已经把这些包都装好了。如果你还没装好这些包，那么，运行install.packages(“R包名称”)这个操作，或者你可能需要更新版本，运行update.packages()。

library(aod)

library(ggplot2)

library(Rcpp)

版本说明：本文基于R3.2.3进行测试的，一下是包的版本：

Rcpp：0.12.3 ggplot2：2.0.0 knitr：1.12.3

记住：本文意在要你知道如何用相关的指令进行逻辑回归分析，而并没有涵盖所有研究人员可能会做的事情，尤其是数据是没有进行清洗和查阅的，而且假设并非最严谨，其它方面也不会相当的标准。

案例

案例1： 假如我们对这些因素感兴趣，它们表示政治候选人是否赢得选举的因子，其中，我们把结果变量表示为0或1，也可以表达成赢或者输。而预测变量的利益可由一场运动中所投入的金额表示，或者是选举人所花的时间，再或者在职人员是否获得足够的支持。

案例2：一个研究者可能会对GRE（研究生入学考试成绩）、GPA（大学平均绩点），以及研究生学院的名誉感兴趣，因为它们影响学校的招生问题。这里，我们用允许/不允许这个二进制结果表示其因变量。

数据的描述

对于我们接下来要进行的数据分析来说，我们要对案例2的入学问题进行深入的探讨。我们有了通常情况下假设所产生的数据，而它们可从R的相关网站得到。记住，当我们要调用相关函数导入数据的时候，如果我们要具体表示一个硬盘驱动器的文件，我们需要打上斜杠（/），而不是反斜杠（\）。

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

这个数据集有二进制的结果（输出值，依赖），它表示允许。这里有3个预测变量：gre、gpa以及rank。我们把gre和gpa看作是连续变量。rank表示有4个值为1。这里，为0的那所学校声望最高，其它的这4所高校声望最低。这时，我们可以用summary()函数来汇总一下这个数据集的情况，而且，如果想要计算里面的标准差，我们可以使用sapply()函数，并在里面写上sd来获取其标准差。

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

你可能会考虑到的分析方法

接下来，我会列举一些你可能会用到的方法。这里所列举的一些方法相对来说比较合理，毕竟有些其他方法不能执行或者尤其局限性。

1.逻辑回归，也是本文的重点

2.概率回归。概率分析所产生的结果类似于逻辑回归。我们可以依据我们的需要进行有选择性的进行概率回归或逻辑回归。

3.最小二乘法。当你使用二进制输出变量的时候，这种模型就是常用于线性回归分析的，而且也可以用它来进行条件概率的运算。然而，这里的误差（例如残差）来自于线性概率模型，而且会影响到异方差性最小二乘法回归的正态误差检验，它也影响无效标准误差和假设性分析。想要了解更多关于线性概率模型的相关问题，可以查阅Long（1997，p.38-40）。

4.双组判别函数分析。这是一个多元的方法处理输出变量。

5.Hotelling 的T2。0/1的输出结果转换到分组变量中，而前一个预测值则成为输出变量。这样可以进行一个显著性测试，然而并没有给没个变量分别给一个系数，而且这对于某个变量根据影响情况调整到其它变量的程度是不清楚的。

使用逻辑模型

接下来的代码，通过使用glm()函数（广义线性模型）进行相关评估。首先，我们要把rank（秩）转换成因子，并预示着rank在这里被视为分类变量。

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

由于我们得到了模型的名字（mylogit)，而R并不会从我们的回归中产生任何输出结果。为了得到结果，我们使用summary()函数进行提取：

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.在上面的结果中，我们首先看到的就是call，这提示我们这时的R在运行什么东西，我们设定了什么选项，等等。

2.接下来，我们看到了偏差残差，用于测量模型的拟合度。这部分的结果显示了这个分布的偏差残差，而针对每个使用在模型里的个案。接下来，我们讨论一下如何汇总偏差估计，从而知道模型的拟合结果。

3.下一部分的结果显示了相关系数，标准误差，z统计量（有时也叫Wald Z统计量），以及它的相关结果。gre和gpa都是同等重要的统计量，并作为rank的3个变量。逻辑回归系数给改变了在预测变量中增加一个单位的输出结果误差。

对于gre，每改变一个单位，输出结果的允许误差（相对于不允许来说）增加0.002。

对于gpa，每改变一个单位，入学的研究生的允许误差增加0.804。

ranK的指示变量的观察方法就有点不同了。例如，层在本科学习过的在rank的值是2，其它为1，使之允许误差变为-0.675。

1. 4.下面的图表里的系数表示不错的拟合效果，包含了空值、偏差残差以及AIC。后面，我们就会学习如何用这些信息来判断这个模型是否拟合。

我们可以使用confint（）函数来获取相关的区间的预测信息。对于逻辑回归模型来说，其置信区间时基于异形对数似然函数求出来的。同样，我们得到的CL是基于默认方法求出的标准差算的。

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

我们可以调用aod包里的wald.test()函数来测试rank的所有影响。图表里的系数的顺序和模型里的项顺序是一样的。这样很重要，因为wald.test()函数就是基于这些模型的项顺序进行测试的。b提供了系数，而Sigma提供了误差项的方差协方差矩阵，而且最终，这些项告诉R哪些项用来进行测试，而在这种情况下，第4、5、6这三项作为rank的层次进行测试。

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

卡方检验算出来的值是20.9，这里涉及到3个自由度，p值算出来是0.00011，这预示着我们所假设的这些项之间具有显著的影响效果。

我们也可以测试而外的假设，这些假设包含rank里不同层次的差异。下面，我们通过测试得知rank=2的测试结果和rank=3的时候一样。下面的第一行代码创造的向量是1，这定义了测试里我们要执行的内容。在这种情况下，我们想要测试rank=2和rank=3（即模型的第4项、第5项）这两个不同的项所产生的差异（减）。为了对比这2项，我们对其中一项乘以1，另一项乘以-1。其它项不包含在测试里，所以它们统一都乘上0。第二行代码写上L=1，这告诉R基于向量为1时执行这次测试（并不是我们之前所选择的进行测试）。

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

卡方检验所在自由度为1的情况下算出来的结果为5.5，并得出相关的p值为0.019，这预示着rank=2和rank=3之间存在着显著的差异。

你一可以把系数指数化，并从误差率进行解读，而R会自动帮你算出来。为了要得到指数化系数，你可以告诉R你想要指数化（exp），而R也会按照你的要求把它们指数化，它属于mylogit（coef（mylogit））的一部分。我们可以使用相同的误差率以及它们的置信区间，并可先前就把置信区间指数化。把它们都放到其中一个图表，我们可以使用cbind系数和系数置信区间这些列整合起来。

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

现在，我们可以说gpa增加了一个单位，而研究生入学（反之就是没有入学的）的误差则在因子上增加2.23。对于更多关于误差率信息的解读，查看我们的FAQ页 How do I interpret odds ratios in logistic regression? 。注意，当R产生了这个结果时，关于误差率的拦截一般都不被解读。

你也可以使用预测概率来帮助你解读这个模型。预测概率可以均可由分类预测变量或连续预测变量计算出来。为了算出预测概率，我们首先要创建含有我们需要的独立变量来创建新的数据框。

我们将要开始计算每个rank值的预测概率的允许值，并计算gre和gpa的平均值。首先，我们要创建新的数据框：

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

这些值必须含有和你之前所创建的逻辑回归分析相同的名字（例如，在这个例子中，gre的就必须命名为gre）。既然，我们现在已经创建好了我们需要进行运算的数据框，那么我们可以告诉R根据这个来创建预测概率。第一行代码经过了压缩，我们现在就把它分开来，讨论这些值是怎样执行的。Newdata$rankP告诉R我们要根据数据集（数据框）的newdata1里的rankP创建一个新的变量剩余的指令告诉R这些rankP值应当使用prediction()函数进行预测。圆括号里的选项告诉R这些预测值应该基于分析mylogit进行预测，mylogit的值源自newdata1以及它的预测值类型就是预测概率（type=”response”)。第二行代码列举了数据框newdata1的值，尽管它不是十分理想，而这就是图表的预测概率。

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

在上面的预测结果中，我们看到来自最好的名校（rank=1)并被接收到研究生的本科生预测概率是0.52，而0.18的学生来自最低档次的学校（rank=4），以gre和gpa作为平均值。我们可以做相似的事情来创建一个针对不断变化的变量gre和gpa的预测概率图表。我们可以基于此作图，所以我们可以在200到800之间创建100个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))))
这些代码所产生的预测概率（下面第一行）和之前算的一样，除非我们还想要对标准差进行要求，否则我们可以对置信区间作图。我们可以对关联规模进行预测，同时反向变换预测值和置信区间的临近值到概率中。

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

当然，使用图像描绘预测概率来解读和展示模型也是相当有用的。我们会使用ggplot2包来作图。下面我们作图描绘预测概率，和95%置信区间。

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)

我们也许很想看到这个模型的拟合效果怎么样，而刚刚那样做是非常有用的，尤其是对比这些竞争的模型。Summary(mylogit)所产生的结果包含了拟合系数（下面展示了其系数），包含了空值、偏差残差和AIC。模型拟合度的一个衡量标准就是整个模型的显著程度。这个测试问了我们使用了预测值的模型是否比仅仅含有截距的模型（即，空模型）更加显著，而这个测试里的统计量是含有预测值的预测拟合指数和空模型的差，而且这个模型的统计量也是含有其自由度等于现有的模型（即，模型里的预测变量个数）和空模型的自由度之差的卡方分布。为了要找到这两个模型的偏差（即，这个测试的统计量），我们可以使用下面的指令：文章来源cda数据分析培训

1. with(mylogit,null.deviance-deviance)
2. ## [1] 41.5

这两个模型的自由度之差等于这个模型里预测变量的个数，而且我们可以按照下面的方法获取它：

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，而相关的p值则小于0.001。这告诉我们，此模型的拟合效果比一个空模型所产生的拟合效果更加显著。这个，我们有时称它为似然比测试（偏差残差为-2*log似然值）。要看到它的对数似然，我们可以这样写：

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