这一节里我们将要讨论包含多个预测变量的回归分析问题。不过模型设定和结果输出等内容与前面系列讲过的关于回归分析和方差分析的内容差别不大，链接：R语言系列第四期：②R语言多组样本方差分析与KW检验、R语言系列第四期：④R语言简单相关与回归。
因此我们这一节的内容新颖之处就在于模型探索方面。
A. 多维数据绘图
下面以Altman提到的一项关于囊胞性纤维症患者的肺功能的研究作为例子，数据集是ISwR包中的cystifibr。
使用pairs()函数可以绘制数据集中任意两个变量之间的散点图：
> library(ISwR)
> par(mex=0.5)
>pairs(cystfibr,gap=0,cex.labels=0.9)
#Tips：首先par()函数的一个参数mex用来调整图形边界的行间距。而pairs()的参数gap用来改变各个子图之间的图间距，cex.labels用来缩放图中的字号。当前这种情况下，就是把cystfibr数据集里的所有变量都进行两两的相关分析，依据位置查看结果。
上面的图形中各个子图虽然相对较小，但是这个图形提供了一种获知多维数据整体情况的有效方法。例如，从图中可以清晰地看到变量age、height和weight强相关关系。
为了便于直接引用cystfibr数据集中的变量，我们可以将该数据集加入到当前的搜索路径中：
> attach(cystfibr)
The following object is masked from package:ISwR:
tlc
因为我们这个部分会用到很多诸如age、weight、height这样的常见变量名，所以最好确保当前工作空间中是否存在同名的变量。（无论大小写）
B. 模型设定和模型输出
多元回归分析的模型设定是通过在模型公式中的解释变量之间添加“ ”来完成的：
lm(pemax~age sex height weight bmp fev1 rv frc tlc)
上面的公式意味着变量pemax可由一个由变量age、sex及其他变量组成的模型来描述（pemax是指患者的最大呼气压力，数据集cystfibr中其他变量的解释可以参考R中的数据集解释）
与之前谈到简单回归一样，lm函数返回的结果有限。我们可以借助summary()函数可以得到更多有趣的输出结果：
> summary(lm(pemax~age sex height weight bmp fev1 rv frc tlc))
Call:
lm(formula = pemax ~ age sex height weight bmp fev1
rv frc tlc)
Residuals:
Min 1Q Median 3Q Max
-37.338 -11.532 1.081 13.386 33.405
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 176.0582 225.8912 0.779 0.448
age -2.5420 4.8017 -0.529 0.604
sex -3.7368 15.4598 -0.242 0.812
height -0.4463 0.9034 -0.494 0.628
weight 2.9928 2.0080 1.490 0.157
bmp -1.7449 1.1552 -1.510 0.152
fev1 1.0807 1.0809 1.000 0.333
rv 0.1970 0.1962 1.004 0.331
frc -0.3084 0.4924 -0.626 0.540
tlc 0.1886 0.4997 0.377 0.711
Residual standard error: 25.47 on 15 degrees of freedom
Multiple R-squared: 0.6373, Adjusted R-squared: 0.4197
F-statistic: 2.929 on 9 and 15 DF, p-value: 0.03195
#Tips：注意，上面结果表明所有变量对应的t值都不显著，但是，联合F检验的结果却是显著的，原因在于t检验说明的仅仅是当从模型中删除某个变量而保留其他变量时模型的变化结果；对于变量在简化模型中是否统计显著，则没有做出说明；t检验认为没有一个变量是不能从模型中删除的。
通过Anova函数可以得到多元回归分析对应的方差分析表，该表给出的结果就跟上面的结果截然不同：
> anova(lm(pemax~age sex height weight bmp fev1 rv frc tlc))
Analysis of Variance Table
Response: pemax
Df Sum Sq Mean Sq F value Pr(>F)
age 1 10098.5 10098.5 15.5661 0.001296 **
sex 1 955.4 955.4 1.4727 0.243680
height 1 155.0 155.0 0.2389 0.632089
weight 1 632.3 632.3 0.9747 0.339170
bmp 1 2862.2 2862.2 4.4119 0.053010 .
fev1 1 1549.1 1549.1 2.3878 0.143120
rv 1 561.9 561.9 0.8662 0.366757
frc 1 194.6 194.6 0.2999 0.592007
tlc 1 92.4 92.4 0.1424 0.711160
Residuals 15 9731.2 648.7
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：除了最后一行（对应变量tlc）之外，这里的F检验结果和summary函数输出的t检验结果几乎是完全相悖的。Age变量的检验结果变得显著了，导致这种结果的原因在于这里的检验过程是逐步进行的。
Anova表的输出结果表明在模型中已包含age变量的情况下，再添加其他变量，模型准确度并未得到显著提高。可以进行联合检验，看看是否可以将age之外的变量全部去掉，做法是求贡献值的平方和，再对总和进行F检验，对应的程序如下：
> m1<-lm(pemax~age sex height weight bmp fev1 rv frc tlc)
> m2<-lm(pemax~age)
> anova(m1,m2)
Analysis of Variance Table
Model 1: pemax ~ age sex height weight bmp fev1 rv frc
tlc
Model 2: pemax ~ age
Res.Df RSS Df Sum of Sq F Pr(>F)
1 15 9731.2
2 23 16734.2 -8 -7002.9 1.3493 0.2936
上述结果中给出的是近似的F检验。
#Tips：上述两个模型应该是嵌套关系。从方差分析表可以看出，删除“age”外的其他变量是合理的，这里的p值0.2936指代的是模型除age外其他变量的显著性，显然是无统计学意义的。
C. 模型筛选
R中有一个按照赤池信息准则（Akaike Information Criterion）进行模型筛选的函数step()。不过我们也可以人工向后消元：
> summary(lm(pemax~age sex height weight bmp fev1 rv frc tlc))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 176.0582 225.8912 0.779 0.448
age -2.5420 4.8017 -0.529 0.604
sex -3.7368 15.4598 -0.242 0.812
height -0.4463 0.9034 -0.494 0.628
weight 2.9928 2.0080 1.490 0.157
bmp -1.7449 1.1552 -1.510 0.152
fev1 1.0807 1.0809 1.000 0.333
rv 0.1970 0.1962 1.004 0.331
frc -0.3084 0.4924 -0.626 0.540
tlc 0.1886 0.4997 0.377 0.711
人工进行模型降维的优点在于可以在该过程自由判断变量的取舍。显然我们会考虑把最后一个变量tlc优先剔除。
> summary(lm(pemax~age sex height weight bmp fev1 rv frc))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 221.8055 185.4350 1.196 0.2491
age -3.1346 4.4144 -0.710 0.4879
sex -4.6933 14.8363 -0.316 0.7558
height -0.5428 0.8428 -0.644 0.5286
weight 3.3157 1.7672 1.876 0.0790 .
bmp -1.9403 1.0047 -1.931 0.0714 .
fev1 1.0183 1.0392 0.980 0.3417
rv 0.1857 0.1887 0.984 0.3396
frc -0.2605 0.4628 -0.563 0.5813
---
> summary(lm(pemax~age sex height weight bmp fev1 rv))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 166.71822 154.31294 1.080 0.2951
age -1.81783 3.66773 -0.496 0.6265
sex 0.10239 11.89990 0.009 0.9932
height -0.40981 0.79257 -0.517 0.6118
weight 2.87386 1.55120 1.853 0.0814 .
bmp -1.94971 0.98415 -1.981 0.0640 .
fev1 1.41526 0.74788 1.892 0.0756 .
rv 0.09567 0.09798 0.976 0.3425
---
> summary(lm(pemax~age sex height weight bmp))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 280.4482 124.9556 2.244 0.0369 *
age -3.0750 3.6352 -0.846 0.4081
sex -11.5281 10.3720 -1.111 0.2802
height -0.6853 0.7962 -0.861 0.4001
weight 3.5546 1.5281 2.326 0.0312 *
bmp -1.9613 0.9263 -2.117 0.0476 *
---
最后我们根据我们的标准0.05，选择了保留一个变量weight，这个变量最终的p值是0.000646，其实这个方法并不是标准的，存在偶然性，况且这个方法需要很多次的计算。而step()的方法可以这样：
> a1<-lm(pemax~age sex height weight bmp fev1 rv frc tlc))
> final=step(a1,direction="both")
Start: AIC=169.11
pemax ~ age sex height weight bmp fev1 rv frc
tlc
Df Sum of Sq RSS AIC
- sex 1 37.90 9769.2 167.20
- tlc 1 92.40 9823.7 167.34
- height 1 158.32 9889.6 167.51
- age 1 181.81 9913.1 167.57
- frc 1 254.55 9985.8 167.75
- fev1 1 648.45 10379.7 168.72
- rv 1 653.78 10385.0 168.73
<none> 9731.2 169.11
- weight 1 1441.21 11172.5 170.56
- bmp 1 1480.12 11211.4 170.65
Step: AIC=167.2
pemax ~ age height weight bmp fev1 rv frc tlc
Df Sum of Sq RSS AIC
- tlc 1 115.94 9885.1 165.50
- height 1 131.21 9900.4 165.54
- age 1 145.55 9914.7 165.57
- frc 1 221.50 9990.7 165.76
- rv 1 636.15 10405.3 166.78
<none> 9769.2 167.20
- weight 1 1446.23 11215.4 168.65
- bmp 1 1474.72 11243.9 168.72
sex 1 37.90 9731.2 169.11
- fev1 1 1770.40 11539.6 169.37
Step: AIC=165.5
pemax ~ age height weight bmp fev1 rv frc
Df Sum of Sq RSS AIC
- frc 1 133.16 10018.3 163.83
- height 1 215.79 10100.9 164.04
- age 1 252.20 10137.3 164.13
- rv 1 543.55 10428.6 164.84
<none> 9885.1 165.50
tlc 1 115.94 9769.2 167.20
sex 1 61.44 9823.7 167.34
- fev1 1 1727.40 11612.5 167.52
- weight 1 2132.53 12017.6 168.38
- bmp 1 2354.31 12239.4 168.84
Step: AIC=163.83
pemax ~ age height weight bmp fev1 rv
Df Sum of Sq RSS AIC
- age 1 145.34 10163.6 162.19
- height 1 158.20 10176.5 162.22
- rv 1 568.07 10586.3 163.21
<none> 10018.3 163.83
frc 1 133.16 9885.1 165.50
tlc 1 27.59 9990.7 165.76
sex 1 0.04 10018.2 165.83
- weight 1 2027.22 12045.5 166.44
- bmp 1 2324.06 12342.3 167.05
- fev1 1 2851.24 12869.5 168.09
Step: AIC=162.19
pemax ~ height weight bmp fev1 rv
Df Sum of Sq RSS AIC
- height 1 191.0 10355 160.66
- rv 1 829.0 10993 162.15
<none> 10164 162.19
age 1 145.3 10018 163.83
tlc 1 102.3 10061 163.94
frc 1 26.3 10137 164.13
sex 1 0.6 10163 164.19
- weight 1 2603.5 12767 165.89
- bmp 1 2743.5 12907 166.17
- fev1 1 3210.9 13374 167.06
Step: AIC=160.66
pemax ~ weight bmp fev1 rv
Df Sum of Sq RSS AIC
<none> 10355 160.66
- rv 1 1183.6 11538 161.36
tlc 1 197.1 10158 162.18
height 1 191.0 10164 162.19
age 1 178.1 10176 162.22
frc 1 3.4 10351 162.65
sex 1 2.4 10352 162.65
- bmp 1 3072.6 13427 165.15
- fev1 1 3717.1 14072 166.33
- weight 1 10930.2 21285 176.67
> summary(final)
Call:
lm(formula = pemax ~ weight bmp fev1 rv)
Residuals:
Min 1Q Median 3Q Max
-39.77 -11.74 4.33 15.66 35.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 63.94669 53.27673 1.200 0.244057
weight 1.74891 0.38063 4.595 0.000175 ***
bmp -1.37724 0.56534 -2.436 0.024322 *
fev1 1.54770 0.57761 2.679 0.014410 *
rv 0.12572 0.08315 1.512 0.146178
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 22.75 on 20 degrees of freedom
Multiple R-squared: 0.6141, Adjusted R-squared: 0.5369
F-statistic: 7.957 on 4 and 20 DF, p-value: 0.000523
这个方法就是按照AIC的值进行选择的。Direction有三个选择“both”“forward”和“backward”分别代表逐步法、向前法和向后法，这里就是我们统计上采用的逐步法筛选变量。
当然，筛选变量还需要考虑变量的实际意义，统计的结果不能作为唯一的标准。另外，我们在平常使用线性模型中也经遇到一些问题，比如共线性，交互效应等问题，我们会在这个系列的番外——R语言系列5番外为大家介绍。
好了，这部分的内容就先介绍到这里，我们下期再见。
参考资料： 1.《R语言统计入门（第二版）》人民邮电出版社 Peter Dalgaard著 2.《R语言初学者指南》人民邮电出版社 Brian Dennis著 3.Vicky的小笔记本《blooming for you》by Vicky
