R语言统计入门-第十一章
第十一章 多元回归
终于到了在统计中最常见的多元回归了~这本书之前都在讲什么
首先介绍下多元回归的基本模型： \[ y=\beta_0+\beta_1x_1+\cdot\cdot\cdot+\beta_kx_k+\epsilon \] 其中，\(x_1,...,x_k\)等是解释变量（也叫作预测变量）；模型参数\(\beta_1,...,\beta_k\)可通过最小二乘法估计得出（具体参见6.1节）。
11.1 多维数据绘图
下面以Altman（1991）提到的一个关于囊胞性纤维症患者肺功能的研究为例子。数据包含在ISwR包内的cystfibr数据框中。
使用pairs函数可以绘制数据集中任意两个变量间的散点图：
library(ISwR)
par(mex = 0.5)
pairs(cystfibr, gap = 0, cex.labels = 0.9)
在这个代码中，参数gap和cex.lables用来控制图形的外观。前者用来移除各个子图之间的图间距，后者用来缩放图中的字号，而绘图命令mex则用来减少图形边界的行间距。
用plot命令也能做出类似的图：
plot(cystfibr)
用pairs命令得到的那个图，各个子图相对较小，不适合直接放在论文中。不过，这种图形可以清晰的看到多维数据的整体情况，比如，就可以看到age、height和weight具有强相关关系。
为了便于直接引用cystfibr数据集中的变量。可以将该数据集加入到当前的搜索路径中：
attach(cystfibr)## The following object is masked from package:ISwR:
##
## tlc11.2 模型设定和模型输出
多元回归分析的模型是通过在模型公式中的解释变量（~应该就是自变量）之间添加+号来实现：
#注意一点，在之前分工作目录中可能会存在名称相同的变量，比如age这个变量，因此在运行这个命令前，需要清空工作目录。
lm(pemax~age+sex+height+weight+bmp+fev1+rv+frc+tlc)##
## Call:
## lm(formula = pemax ~ age + sex + height + weight + bmp + fev1 +
## rv + frc + tlc)
##
## Coefficients:
## (Intercept) age sex height weight
## 176.0582 -2.5420 -3.7368 -0.4463 2.9928
## bmp fev1 rv frc tlc
## -1.7449 1.0807 0.1970 -0.3084 0.1886这个公式意思就是，pemax这个变量可由一个包含变量age、sex等组成的模型来描述（pemax是患者的最大呼气压力）。
与前面一样，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注意，这个结果表明所有变量对应的 t 值都不显著，但是，联合 F 检验的结果却是显著的，这一定是有原因的。
这个原因就在于，t 检验说明的是仅仅是，当从模型中删除某个变量而保留其他变量时模型的变化结果；对于变量在简化模型中是否统计显著，则没有做出说明；值得注意的是，t 检验认为没有一个变量是不能从模型中删除的。
注意，输出结果中未调整R2（Multiple R-squared）和调整后R2（Adjusted R-squared）有较大差异，这归咎于模型中较多的变量个数，而这个与方差的自由度密切相关。前者表示的是与空模型相对的残差平方和的变化，后者对应的是残差方差的类似变化：
1-25.5^2/var(pemax)## [1] 0.4183949其中，25.5这个数字取自summary函数输出结果中的残差标准误（Residual standard error）。通过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注意，除了最后一行（对应于变量tlc）之外，这里的 F 检验结果与summary函数输出的 t 检验结果几乎完全相悖。这里，age变量的检验结果变得显著了，导致这种结果的原因在于这里的检验过程是逐步进行的；具体而言，对应于（从下至上）将变量逐个从模型中移除，直至剩下age变量。在该过程中，变量bmp的检验结果一度接近5%的临界点，但考虑到检验的个数，这一结果几乎不显著。
在8次独立的检验中，结果给出小于等于0.053的 p 值的概率仅仅略高于35%。虽然ANOVA表中的检验并非完全独立，但是其近似结果还是不错的。
ANOVA表的输出结果表明在模型已包含age变量的情况下，再添加其他变量，模型准确度并未得到显著的提高。可以进行联合检验，看看是否可以将age以外的变量全部去掉，做法是求贡献值的平方和加，再对总和进行 F 检验：
955.4+155.0+632.3+2862.2+1549.1+561.9+194.6+92.4## [1] 7002.97002.9/8## [1] 875.3625875.36/648.7## [1] 1.3494071-pf(1.349407, 8,15)## [1] 0.2935148对应于去掉边框线的表格，它看起来是这样的：
| Df | Sum Sq | Mean Sq | F | Pr(>F) | |
|---|---|---|---|---|---|
| age | 1 | 10098.5 | 10098.5 | 15.566 | 0.00130 |
| others | 8 | 7002.9 | 875.4 | 1.349 | 0.29351 |
| Residual | 15 | 9731.2 | 648.7 |
这个表格是自己根据数据整理的~不是系统的输出。
如要直接得到上述结果，可运行：
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.293611.3 模型筛选
R有一个按照赤池信息准则（Akaike Information Criterion）进行模型筛选的函数step()。而这个本书不讲~~
本书只使用一种较为简单的人工向后消元法。
下面是一个模型降阶的例子，注意，为减少输出结果占用的空间，对输出信息进行了编辑：
#书上是只保留了Coefficients这一项
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人工进行模型降阶的优点在于该模型引入逻辑结构。在本利中，很自然地会先想到去掉肺功能的指标。
summary(lm(pemax~age+sex+height+weight+bmp+fev1+rv+frc))##
## Call:
## lm(formula = pemax ~ age + sex + height + weight + bmp + fev1 +
## rv + frc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.072 -10.067 0.113 13.526 36.990
##
## Coefficients:
## 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
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.78 on 16 degrees of freedom
## Multiple R-squared: 0.6339, Adjusted R-squared: 0.4508
## F-statistic: 3.463 on 8 and 16 DF, p-value: 0.01649summary(lm(pemax~age+sex+height+weight+bmp+fev1+rv))##
## Call:
## lm(formula = pemax ~ age + sex + height + weight + bmp + fev1 +
## rv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.425 -12.391 3.834 14.797 36.693
##
## Coefficients:
## 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
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.28 on 17 degrees of freedom
## Multiple R-squared: 0.6266, Adjusted R-squared: 0.4729
## F-statistic: 4.076 on 7 and 17 DF, p-value: 0.008452summary(lm(pemax~age+sex+height+weight+bmp+fev1))##
## Call:
## lm(formula = pemax ~ age + sex + height + weight + bmp + fev1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.238 -7.403 -0.081 15.534 36.028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 260.6313 120.5215 2.163 0.0443 *
## age -2.9062 3.4898 -0.833 0.4159
## sex -1.2115 11.8083 -0.103 0.9194
## height -0.6067 0.7655 -0.793 0.4384
## weight 3.3463 1.4719 2.273 0.0355 *
## bmp -2.3042 0.9136 -2.522 0.0213 *
## fev1 1.0274 0.6329 1.623 0.1219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.24 on 18 degrees of freedom
## Multiple R-squared: 0.6057, Adjusted R-squared: 0.4743
## F-statistic: 4.608 on 6 and 18 DF, p-value: 0.00529summary(lm(pemax~age+sex+height+weight+bmp))##
## Call:
## lm(formula = pemax ~ age + sex + height + weight + bmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.194 -9.412 -2.425 9.157 40.112
##
## Coefficients:
## 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 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.27 on 19 degrees of freedom
## Multiple R-squared: 0.548, Adjusted R-squared: 0.429
## F-statistic: 4.606 on 5 and 19 DF, p-value: 0.006388上上述结果看，去掉4个肺功能相关的变量没什么不妥，接下来尝试删除那些描述病人身体发育状态或尺寸信息的变量。在开始时，尽量避免删除weight和bmp变量，因为它们对应的 p 值很接近5%的显著性界限。
summary(lm(pemax~age+height+weight+bmp))##
## Call:
## lm(formula = pemax ~ age + height + weight + bmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.501 -15.460 -2.838 11.082 42.991
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 274.5307 125.5745 2.186 0.0409 *
## age -3.0832 3.6566 -0.843 0.4091
## height -0.6985 0.8008 -0.872 0.3934
## weight 3.6338 1.5354 2.367 0.0282 *
## bmp -1.9621 0.9317 -2.106 0.0480 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.41 on 20 degrees of freedom
## Multiple R-squared: 0.5186, Adjusted R-squared: 0.4223
## F-statistic: 5.386 on 4 and 20 DF, p-value: 0.004137summary(lm(pemax~height+weight+bmp))##
## Call:
## lm(formula = pemax ~ height + weight + bmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.794 -11.764 -1.218 13.202 43.631
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 245.3936 119.8927 2.047 0.0534 .
## height -0.8264 0.7808 -1.058 0.3019
## weight 2.7717 1.1377 2.436 0.0238 *
## bmp -1.4876 0.7375 -2.017 0.0566 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.24 on 21 degrees of freedom
## Multiple R-squared: 0.5015, Adjusted R-squared: 0.4302
## F-statistic: 7.041 on 3 and 21 DF, p-value: 0.00187summary(lm(pemax~weight+bmp))##
## Call:
## lm(formula = pemax ~ weight + bmp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.924 -13.399 4.361 16.642 48.404
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 124.8297 37.4786 3.331 0.003033 **
## weight 1.6403 0.3900 4.206 0.000365 ***
## bmp -1.0054 0.5814 -1.729 0.097797 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.31 on 22 degrees of freedom
## Multiple R-squared: 0.4749, Adjusted R-squared: 0.4271
## F-statistic: 9.947 on 2 and 22 DF, p-value: 0.0008374summary(lm(pemax~weight))##
## Call:
## lm(formula = pemax ~ weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.30 -22.69 2.23 15.91 48.41
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.5456 12.7016 5.003 4.63e-05 ***
## weight 1.1867 0.3009 3.944 0.000646 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.38 on 23 degrees of freedom
## Multiple R-squared: 0.4035, Adjusted R-squared: 0.3776
## F-statistic: 15.56 on 1 and 23 DF, p-value: 0.0006457注意，一旦删除age和height变量，变量bmp就不再显著了。在原文献（Altman，1991）中，变量weight、feval和bmp在最终结果中对应的 p值都低于5%。然而，并非所有模型降阶过程都是如此。
特别关注变量age、weight和height是个不错的想法，因为在处理儿童和青少年对应的数据时，这些变量表现出很强的相关性。
summary(lm(pemax~age+height+weight))##
## Call:
## lm(formula = pemax ~ age + height + weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.675 -21.566 3.229 16.274 48.068
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.65555 82.40935 0.785 0.441
## age 1.56755 3.14363 0.499 0.623
## height -0.07608 0.80278 -0.095 0.925
## weight 0.86949 0.85922 1.012 0.323
##
## Residual standard error: 27.41 on 21 degrees of freedom
## Multiple R-squared: 0.4118, Adjusted R-squared: 0.3278
## F-statistic: 4.901 on 3 and 21 DF, p-value: 0.009776summary(lm(pemax~age+height))##
## Call:
## lm(formula = pemax ~ age + height)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.817 -17.883 3.815 18.275 53.824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.8600 68.2493 0.262 0.796
## age 2.7178 2.9325 0.927 0.364
## height 0.3397 0.6900 0.492 0.627
##
## Residual standard error: 27.43 on 22 degrees of freedom
## Multiple R-squared: 0.3831, Adjusted R-squared: 0.3271
## F-statistic: 6.832 on 2 and 22 DF, p-value: 0.00492summary(lm(pemax~age))##
## Call:
## lm(formula = pemax ~ age)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.666 -17.174 6.209 16.209 51.334
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.408 16.657 3.026 0.00601 **
## age 4.055 1.088 3.726 0.00111 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.97 on 23 degrees of freedom
## Multiple R-squared: 0.3764, Adjusted R-squared: 0.3492
## F-statistic: 13.88 on 1 and 23 DF, p-value: 0.001109summary(lm(pemax~height))##
## Call:
## lm(formula = pemax ~ height)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.876 -19.306 1.787 18.170 61.464
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.2757 40.0445 -0.831 0.41453
## height 0.9319 0.2596 3.590 0.00155 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.34 on 23 degrees of freedom
## Multiple R-squared: 0.3591, Adjusted R-squared: 0.3312
## F-statistic: 12.89 on 1 and 23 DF, p-value: 0.001549从上面的结果可以看出，对于变量weight、height和age，没有证据表明哪个比另外两个好。上面所用的消元方法之所以在最后仅留下了weight作为自变量，完全是出于偶然。
练习题
- 在数据集secher中，对出生体重、腰围和二项骨的直径变量进行对数变换后可得到很好的数据分析结果。请拟合出生体重的预测表达式。在模型中同时纳入腹部直径和二项骨直径时，模型结果如何，模型中两个回归系数之和约为3，如何对其进行解释。
- 数据集tlc有一个同名变量tlc，这不是一个很好的命名方式，请解释原因。用数据集中的其余变量来解释tlc变量，并对模型的有效性进行解释。
- 数据集cystfibr的分析过程设计sex变量，它是一个二元变量，如何解释回归结果中对应的系数。
- 考虑juul2数据集，并筛选出该数据集中年龄超过25岁的子集，用age变量对\(\sqrt{igfl}\)变量进行回归分析。在扩展模型中加入变量height和变量weight，计算扩展模型对应的方差分析表，有没有意想不到的结果出现，为什么会这样呢。
- 使用多远回归模型，分析冰结师kfmdata数据集中各个解释变量对牛奶摄入量的影响。注意，这里的sex变量是因子型变量，这对分析过程有什么影响？
