1. sci666首页
  2. 实用技巧
  3. 生物信息学

R语言系列五:⑤R语言与多元回归

这一节里我们将要讨论包含多个预测变量的回归分析问题。不过模型设定和结果输出等内容与前面系列讲过的关于回归分析和方差分析的内容差别不大,链接:R语言系列第四期:②R语言多组样本方差分析与KW检验R语言系列第四期:④R语言简单相关与回归R语言系列第四期(番外篇):样本容量和把握度计算

因此我们这一节的内容新颖之处就在于模型探索方面。

R语言系列五:⑤R语言与多元回归
A. 多维数据绘图
R语言系列五:⑤R语言与多元回归

下面以Altman提到的一项关于囊胞性纤维症患者的肺功能的研究作为例子,数据集是ISwR包中的cystifibr

使用pairs()函数可以绘制数据集中任意两个变量之间的散点图:

> library(ISwR)

> par(mex=0.5)

>pairs(cystfibr,gap=0,cex.labels=0.9)

R语言系列五:⑤R语言与多元回归

#Tips:首先par()函数的一个参数mex用来调整图形边界的行间距。而pairs()的参数gap用来改变各个子图之间的图间距,cex.labels用来缩放图中的字号。当前这种情况下,就是把cystfibr数据集里的所有变量都进行两两的相关分析,依据位置查看结果。

上面的图形中各个子图虽然相对较小,但是这个图形提供了一种获知多维数据整体情况的有效方法。例如,从图中可以清晰地看到变量ageheightweight强相关关系。

为了便于直接引用cystfibr数据集中的变量,我们可以将该数据集加入到当前的搜索路径中:

> attach(cystfibr)

The following object is masked from package:ISwR:

    tlc

因为我们这个部分会用到很多诸如ageweightheight这样的常见变量名,所以最好确保当前工作空间中是否存在同名的变量。(无论大小写)

R语言系列五:⑤R语言与多元回归
B. 模型设定和模型输出
R语言系列五:⑤R语言与多元回归

多元回归分析的模型设定是通过在模型公式中的解释变量之间添加“+”来完成的:

lm(pemax~age+sex+height+weight+bmp+fev1+rv+frc+tlc)

上面的公式意味着变量pemax可由一个由变量agesex及其他变量组成的模型来描述(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外其他变量的显著性,显然是无统计学意义的。

R语言系列五:⑤R语言与多元回归
C. 模型筛选
R语言系列五:⑤R语言与多元回归

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

这些可能会帮助到你: 问答社区 | 共享百度SVIP | 留言建议

欢迎入群交流:生信分析群: 732179952 · Meta分析群: 797345521

发表评论

登录后才能评论