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”分别代表逐步法、向前法和向后法，这里就是我们统计上采用的逐步法筛选变量。

当然，筛选变量还需要考虑变量的实际意义，统计的结果不能作为唯一的标准。

参考资料： 1.《R语言统计入门（第二版）》人民邮电出版社 Peter Dalgaard著 2.《R语言初学者指南》人民邮电出版社 Brian Dennis著