R语言系列第五期：③R语言逻辑回归预测和检验

在上一篇文章里，无论原始数据是表格式的还是罗列式的，我们都可以建立起相应的逻辑回归模型。详情点击：R语言系列五：②R语言与逻辑回归建立

但是模型建立起来之后，是用来做什么的？我们当然需要利用模型来解释变量，但是我们也可以利用模型来预测结局，我们建立起来模型之后，可以通过个人的数据来计算这个人发生阳性事件的概率大小，从而最终给出结局分类，并且做出相应的对策。

我们首先考虑之前的高血压的例子，这个例子中共有8个分类组合的水平，我们为了方便后续的操作，我们把上一节的表列在这里：

smoking obesity snoring n.tot n.hyp

1 No No No 60 5

2 Yes No No 17 2

3 No Yes No 8 1

4 Yes Yes No 2 0

5 No No Yes 187 35

6 Yes No Yes 85 13

7 No Yes Yes 51 15

8 Yes Yes Yes 23 8

我们接着上一部分的内容，把smoking变量给删除掉重新建立模型，而我们做预测的函数就是predict(),我们得到的预测结果是以列表的形式给出：

> glm.hyp=glm(hyp.tbl~obesity+snoring,family=binomial("logit"))

Call: glm(formula = hyp.tbl ~ obesity + snoring, family = binomial("logit"))

Coefficients:

(Intercept) obesityYes snoringYes

-2.3921 0.6954 0.8655

Degrees of Freedom: 7 Total (i.e. Null); 5 Residual

Null Deviance: 14.13

Residual Deviance: 1.678 AIC: 32.6

> summary(glm.hyp)

Call:

glm(formula = hyp.tbl ~ obesity + snoring, family = binomial("logit"))

Deviance Residuals:

1 2 3 4 5 6 7 8

-0.01247 0.47756 -0.24050 -0.82050 0.30794 -0.62742 -0.14449 0.45770

Coefficients:

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

(Intercept) -2.3921 0.3757 -6.366 1.94e-10 ***

obesityYes 0.6954 0.2851 2.440 0.0147 *

snoringYes 0.8655 0.3967 2.182 0.0291 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 14.1259 on 7 degrees of freedom

Residual deviance: 1.6781 on 5 degrees of freedom

AIC: 32.597

Number of Fisher Scoring iterations: 4

> predict(glm.hyp)

1 2 3 4

-2.3920763 -2.3920763 -1.6966575 -1.6966575

5 6 7 8

-1.5266180 -1.5266180 -0.8311991 -0.8311991

#Tips：当前代码是在上一节的内容的基础上进行的！

#Tips：首先，我们删掉了smoking变量之后，剩下的两个变量都是有统计学意义的变量。而删掉了smoking之后，预测的结果就是两两相等的。另外这些值都是在logit刻度下的，所以我们再去计算这个公式里P值解之后就可以得到阳性率。比如第一个数值，-2.3920763，我们由下面公式计算得到P值是0.08378.

#Tips：另一件事是，你会发现2.392-1.697=1.527-0.831=0.695（四舍五入的误差不考虑），这个值是obesity的回归系数，而snoring的回归系数则为2.392-1.527=1.697-0.831=0.866。

当然，我们也可以直接显示概率刻度下的预测值，需要在predict函数中设定参数type=“response”：

> predict(glm.hyp,type="response")

1 2 3 4

0.08377892 0.08377892 0.15490233 0.15490233

5 6 7 8

0.17848906 0.17848906 0.30339158 0.30339158

另一种情况，连续型自变量里，比如menarche的分析中，我们可以给每个点都计算概率，但是对于我们而言没有什么意义，我们的主要兴趣点在于看到年龄和期望概率的关系图。我们可以用如下代码画一下：

> glm.menarche<-glm(menarche~age,binomial)

> Age<-seq(8,20,0.1)

> newages<-data.frame(age=Age)

> predicted.probability<-predict(glm.menarche,newages,type="resp")

> plot(predicted.probability~Age,type="l")

#Tips：Age变量是用来做横轴的点，seq()函数生成等距元素的向量，这里年龄是从8-20岁，间隔为0.1，所以点连起来会很光滑。predicted.probability可以通过我们第一行代码建立的回归模型，计算Age变量所对应的每个点的概率值，这样我们每个x，y都计算出来之后，使用plot()函数就可以把他们绘制出来了。

B. 模型检查

当然，我们建立了模型之后，肯定要利用模型说明问题，但是我们建立的模型到底好不好，我们又必须给出适当的判断。

对于表格式的数据，很明显，我们应该去比较观测和拟合出来的值的占比。在前面高血压的例子中，我们可以计算各组水平概率（下面的是实际概率）：

> fitted(glm.hyp)

1 2 3 4

0.08377892 0.08377892 0.15490233 0.15490233

5 6 7 8

0.17848906 0.17848906 0.30339158 0.30339158

#Tips：fitted()函数是用来计算利用原建模数据和已建成模型所得最终预测值的函数。fitted()函数和predict()函数很类似，但是它不可以利用外部数据计算新的概率。下面是实际的对应阳性事件发生率。

> prop.hyp

[1] 0.08333333 0.11764706 0.12500000 0.00000000

[5] 0.18716578 0.15294118 0.29411765 0.34782609

问题在于，你不知道相对频率是如何计算得到的，如果能看一看观测到的频数和，估计得到的数量会更好一些：

> fitted(glm.hyp)*n.tot

1 2 3 4

5.0267351 1.4242416 1.2392186 0.3098047

5 6 7 8

33.3774535 15.1715698 15.4729705 6.9780063

我们可以把比较的内容以更好的方式打印出来：

> data.frame(fit=fitted(glm.hyp)*n.tot,n.hyp,n.tot)

fit n.hyp n.tot

1 5.0267351 5 60

2 1.4242416 2 17

3 1.2392186 1 8

4 0.3098047 0 2

5 33.3774535 35 187

6 15.1715698 13 85

7 15.4729705 15 51

8 6.9780063 8 23

#Tips：第四个观测，出现15%的期望和0%的观测这样的差别，虽然差别大，但是模型预测里只有0.3个高血压，而实际是0个，所以要注意相对数使用时容易出现这样的问题。

对于连续多变量的复杂模型，对其进行充分的检查是蛮困难的。当观测数据只有两个不同的取值时，就没有方法利用残差图进行检验了。

我们利用menarche的例子。我们试着将x轴划分为几个区间，然后看看每个区间里的点的数量占比与估计的概率之间是否相符：

> age.group<-cut(age,c(8,10,12,13,14,15,16,18,20))

> tb<-table(age.group,menarche)

> tb

menarche

age.group No Yes

(8,10] 100 0

(10,12] 97 4

(12,13] 32 21

(13,14] 22 20

(14,15] 5 36

(15,16] 0 31

(16,18] 0 105

(18,20] 0 46

> rel.freq<-prop.table(tb,1)[,2]

> rel.freq

(8,10] (10,12] (12,13] (13,14]

0.00000000 0.03960396 0.39622642 0.47619048

(14,15] (15,16] (16,18] (18,20]

0.87804878 1.00000000 1.00000000 1.00000000

> points(rel.freq~c(9,11,12.5,13.5,14.5,15.5,17,19),pch=5)

#Tips：我们需要对上面所用到的方法进行解释，首先，cut()函数用于定义因子age.group，它把原来的连续性数据划分成多个组，有点像factor()函数，但是它会给出明确的分割点，而factor()是原先就分好的类。cut()的第二个参数就是分割点；然后table()函数利用age.table变量和menarche变量进行制表。使用prop.table()函数，我们之前提过，它会计算tb表格中每行行内数据构成比（1表示行，2表示列），随后[,2]表示只保留第二列，即yes的那一列；最后，绘制关于期望概率的图，与观测占比的图叠加起来。我们是在原先的图的基础上把点给叠加上去的。pch参数表示使用什么符号，5代表菱形。

整体来看，这个图还是有意义的，尽管12-13岁年龄段和13-14年龄段原始数据和预测数据略有差池。

但是这样的偏差是否有统计学意义呢？其实这个方法并不能非常好的给出结果，我们只能直观的感觉预测效果如何而异。

这里推荐使用一个检验模型很好的工具ROC曲线，我们可以一步一步告诉你ROC曲线是如何画出来的：

> glm.menarche<-glm(menarche~age,binomial)

# 原始模型输入到glm.menarche变量里

> pre<-predict(glm.menarche,type="response")

# 模型对应的每一个点预测值

> pre.obs<-data.frame(prob=pre,obs=menarche)

# 我们把模型的预测值和其对应的原始观测值合并到一个数据框里

> pre.obs<-pre.obs[order(pre.obs$prob),]

# 我们把所有预测和实际值对子按照预测概率的大小由小到大排序

> n<-nrow(pre.obs)

# 计算总行数，就是总的观测数

> tpr<-fpr<-rep(0,n)

# 先将真阳性率和真阴性率设为对应观测数目的0

> for(i in 1:n){

+ threshold<-pre.obs$prob[i]

+ tp<-sum(pre.obs$prob>threshold & pre.obs$obs=="Yes")

+ fp<-sum(pre.obs$prob>threshold & pre.obs$obs=="No")

+ tn<-sum(pre.obs$prob<threshold & pre.obs$obs=="No")

+ fn<-sum(pre.obs$prob<threshold & pre.obs$obs=="Yes")

+ tpr[i]<-tp/(tp+fn)

+ fpr[i]<-fp/(tn+fp)

+ }

# 计算以预测的每个概率为截断概率所对应的实际真阳性率和真阴性率

> plot(fpr,tpr,type="l")

# 把对应的真阳性率（灵敏度）和真阴性率（1-特异度）画点并连线

> abline(a=0,b=1)

# 添加一条过原点和（1,1）的线

#Tips：ROC曲线是诊断逻辑回归模型很好的工具，如果曲线下的面积越大，说明模型的预测能力越强，而这个模型显然是一个很不错的模型。

另外，R里也有专门的ROC绘制包，比如ROCR包（不做详细解释，直接套用就可以了，这个包会直接给出来AUC值，即曲线下面积的大小。另外这部分的数据是在上一个方法所求的部分变量的基础上）：

> library(ROCR)

> pred <- prediction(pre,menarche)

> performance(pred,'auc')@y.values

[[1]]

[1] 0.9768076

> perf <- performance(pred,'tpr','fpr')

> plot(perf)

但是ROCR包画图函数功能比较单一，笔者比较偏好使用功能更强大的pROC包。它可以方便比较两个分类器，还能自动标注出最优的临界点，图看起来也比较漂亮。（同样套用即可。只需要更改的是第二行代码，roc()函数里改成自己的y，pre同样是第一个方法的计算结果。）

> library(pROC)

> modelroc <- roc(menarche,pre)

> plot(modelroc, print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2),

+ grid.col=c("green", "red"),

+ max.auc.polygon=TRUE,auc.polygon.col="skyblue", print.thres=TRUE)

#Tips：这种方法做出来的图形还给出了最优截断点（cut-off点）在左上角，图形的各种参数都可以通过plot()函数里的内容更改，比如颜色等。

关于逻辑回归的建立，预测和诊断我们都已经介绍完了，下面就是关于生存分析的部分了。敬请期待。

参考资料： 1.《R语言统计入门（第二版）》人民邮电出版社 Peter Dalgaard著 2.《R语言初学者指南》人民邮电出版社 Brian Dennis著 3.Vicky的小笔记本《blooming for you》by Vicky