R语言系列第五期（番外篇）：R语言与线性模型相关问题

很多数据集本身非常复杂，按照标准的建模流程难以进行合适的处理，因此，需要构建特别的模型，线性模型提供了一个灵活的模型框架，在此框架内，我们得以对上述大部分复杂数据集拟合模型。

你可能已经注意到，lm()函数既可以应用到分组数据的情况，也可以应用到线性回归问题，详情点击：R语言系列第四期：④R语言简单相关与回归、R语言系列第四期：②R语言多组样本方差分析与KW检验、R语言系列五：①R语言与多元回归

但是，事实上，他们是同一个模型的特例而已。

这个部分包含一些复杂模型以及使用lm()构造模型的过程以及在这个过程中经常出现的问题的处理。

A. 多项式回归

在多元回归里有的时候不像它看起来那么简单，有时可以在多元回归分析中纳入变量的二次和高次幂，尽管这个看似是非线性关系的模型依旧属于线性模型的范畴，重点在于参数和预期的观测值是线性关系。

下面以之前使用的cystic fibrosis 数据集为例。

首先我们展示一下数据：

> cystfibr

   age sex height weight bmp fev1  rv frc tlc pemax

1    7   0    109   13.1  68   32 258 183 137    95

2    7   1    112   12.9  65   19 449 245 134    85

…

24  23   0    175   51.1  71   33 224 131 113    95

25  23   0    179   71.5  95   52 225 127 101   195

另外，我们展示一下之前在线性回归中，所给出的相关图：

展示一下变量pemax和变量height之间的关系。可以看到二者属于非线性关系。通过在模型公式中增加一个变量height的平方可对此进行改造。

> summary(lm(pemax~height+I(height^2)))

Call:

lm(formula = pemax ~ height + I(height^2))

Residuals:

    Min      1Q  Median      3Q     Max

-51.411 -14.932  -2.288  12.787  44.933

Coefficients:

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

(Intercept) 615.36248  240.95580   2.554   0.0181 *

height       -8.08324    3.32052  -2.434   0.0235 *

I(height^2)   0.03064    0.01126   2.721   0.0125 *

—

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

Residual standard error: 24.18 on 22 degrees of freedom

Multiple R-squared:  0.5205,    Adjusted R-squared:  0.4769

F-statistic: 11.94 on 2 and 22 DF,  p-value: 0.0003081

#Tips：需要注意的是，模型公式中的height^2需要用I()函数进行保护。这个技术常常用来防止模型公式中的操作符被特殊解释。这种解释不作用于函数命令内部，I()是反身函数，原封不动地返回自身的输入参数。

使用predict()函数可以绘制带预测值和置信带的拟合曲线。我们采用几个有序的的身高点来做图：

> pred.frame<-data.frame(height=seq(110,180,2))

> lm.pemax.hq<-lm(pemax~height+I(height^2))

> predict(lm.pemax.hq,interval=”pred”,newdata=pred.frame)

         fit      lwr      upr

1   96.90026 37.94461 155.8559

2   94.33611 36.82985 151.8424

…

35 147.21294 93.51117 200.9147

36 152.98174 98.36718 207.5963

> pp<-predict(lm.pemax.hq,newdata=pred.frame,interval=”pred”)

> pc<-predict(lm.pemax.hq,newdata=pred.frame,interval=”conf”)

> plot(height,pemax,ylim=c(0,200))

> matlines(pred.frame$height,pp,lty=c(1,2,2),col=”black”)

> matlines(pred.frame$height,pc,lty=c(1,3,3),col=”black”)

#Tips：我们可以从图里看出来，身高较小时，拟合曲线随着身高的增加而减少，而后随着身高的增加而快速增加。这是为拟合数据选取了二次多项式形式的模型导致的。而身高和pemax变量之间的关系也可以很好地使用二次项来解释。

B. 组间共线性

有时候，数据会根据某个连续尺度的分段进行分组，或者试验设计过程中令变量取几个特定的x值的集合。这两种情况都跟比较线性回归的结果和方差分析的结果相关。

对于同样的数据，我们有两种可供选择的数据模型。两者都属于线性模型的范畴，且都能通过lm()函数拟合。线性回归模型是单因素方差分析模型的子模型，因为前者可以通过向后者的参数添加约束来获得。

可以通过检验较大模型对应的残差方差来决定模型降价是否通过，这就是F检验。

我们以trypsin（胰蛋白酶）浓度的例子里，数据按照年龄进行分组，数据以6个组对应的均值和标准差的形式给出。R没法直接处理这种形式的数据，因此，有必要创造一个均值和标准差都相同的伪数据，对应的代码如下：

> attach(fake.trypsin)

方差分析的真是结果仅依赖于均值和标准差，对生成伪数据的过程感兴趣的可以看一下ISwR文件夹中的fake.trypsin.R文件。

数据框fake.trypsin共包含3个变量，可以运行下面代码查看：

> summary(fake.trypsin)

    trypsin            grp        grpf  

 Min.   :-39.96   Min.   :1.000   1: 32 

 1st Qu.:119.52   1st Qu.:2.000   2:137 

 Median :167.59   Median :2.000   3: 38 

 Mean   :168.68   Mean   :2.583   4: 44 

 3rd Qu.:213.98   3rd Qu.:3.000   5: 16 

 Max.   :390.13   Max.   :6.000   6:  4 

#Tips：数据框既包括数值型向量grp，也包括一个具有6个水平的因子型变量grpf。

对生成的数据执行单因素方差分析得到ANOVA表如下：

> anova(lm(trypsin~grpf))

Analysis of Variance Table

Response: trypsin

           Df Sum Sq Mean Sq F value    Pr(>F)   

grpf        5 224103   44821  13.508 9.592e-12 ***

Residuals 265 879272    3318                     

—

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

如果用变量grp代替模型公式中的grpfin变量，将得到一个关于分组号的线性回归模型。得到的ANOVA表格如下：

> anova(lm(trypsin~grp))

Analysis of Variance Table

Response: trypsin

           Df Sum Sq Mean Sq F value    Pr(>F)   

grp         1 206698  206698  62.009 8.451e-14 ***

Residuals 269 896677    3333                     

—

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

#Tips：残差均方的变化不大，这意味着这两个模型对数据的刻画效果基本上一样。

如果想做一个正规的检验来比较简单线性模型和各组具有独立均值的模型的话，可以直接运行下面代码：

> anova(lm(trypsin~grp+grpf))

Analysis of Variance Table

Response: trypsin

           Df Sum Sq Mean Sq F value    Pr(>F)   

grp         1 206698  206698 62.2959 7.833e-14 ***

grpf        4  17405    4351  1.3114    0.2661

Residuals 265 879272    3318                     

—

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

#Tips：从上面结果可以看出，grpf对模型的拟合没有grp好。而之前他们对整个模型的解释都是有意义的，而grpf在这个模型里的p值突然无意义了，说明了两个变量之间有明显的共线问题。

我们可以绘制带拟合线和经验均值的伪数据的散点图：

> xbar.trypsin<-tapply(trypsin,grpf,mean)

> stripchart(trypsin~grp,method=”jitter”,jitter=0.1,vertical=T,pch=20)

> lines(1:6,xbar.trypsin,type=”b”,pch=4,cex=2,lty=2)

> abline(lm(trypsin~grp))

#Tips：这个图的技巧本质上同之前的图相同，所以这里不再对绘图细节进行详细描述。图中有一个点表示胰岛素密度是负的，这个数据是伪造的，原始数据是没法看到。

C. 交互效应

多元回归模型的一个基本假设就是模型中的各变量对响应变量的影响具有叠加效应。然而，这并不是说线性模型模型没法刻画非叠加效应。我们可以通过添加特殊交互项来指定一个变量受另一个变量水平变动的影响程度。在R中的模型公式里，交互项可以使用“：”来生成，比如a：b。通常，我们还会在模型中包含a和b这两项，同时，R的模型里允许a*b或者a+b+a:b这种公式，这两个公式是等效的。

当然在模型建立的过程中还有很多需要注意很多事项，我们这里就不一一列举了。

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