9.8　用回归来做ANOVA

在9.2节中，我们提到ANOVA和回归都是广义线性模型的特例。因此，本章所有的设计都可以用lm()函数来分析。但是，为了更好地理解输出结果，需要弄明白在拟合模型时，R是如何处理类别型变量的。

以9.3节的单因素ANOVA问题为例，即比较五种降低胆固醇药物疗法（trt）的影响。

  1. > library(multcomp)
  2. > levels(cholesterol$trt)
  4. [1] "1time"  "2times" "4times" "drugD" "drugE"

首先，用aov()函数拟合模型：

  1. > fit.aov <- aov(response ~ trt, data=cholesterol)
  2. > summary(fit.aov)
  4.             Df  Sum Sq Mean Sq F value    Pr(>F)
  5. trt          4 1351.37  337.84  32.433 9.819e-13 ***
  6. Residuals   45  468.75   10.42

现在，用lm()函数拟合同样的模型。结果见代码清单9-11。

代码清单9-11　解决9.3节ANOVA问题的回归方法

  1. > fit.lm <- lm(response ~ trt, data=cholesterol)
  2. > summary(fit.lm)
  4. Coefficients:
  5.             Estimate Std. Error t value Pr(>|t|)
  6. (Intercept)    5.782      1.021   5.665 9.78e-07 ***
  7. trt2times      3.443      1.443   2.385   0.0213 *
  8. trt4times      6.593      1.443   4.568 3.82e-05 ***
  9. trtdrugD       9.579      1.443   6.637 3.53e-08 ***
  10. trtdrugE      15.166      1.443  10.507 1.08e-13 ***
  12. Residual standard error: 3.227 on 45 degrees of freedom
  13. Multiple R-squared: 0.7425,     Adjusted R-squared: 0.7196
  14. F-statistic: 32.43 on 4 and 45 DF,  p-value: 9.819e-13

我们能发现什么？因为线性模型要求预测变量是数值型，当lm()函数碰到因子时，它会用一系列与因子水平相对应的数值型对照变量来代替因子。如果因子有k个水平，将会创建k?1个对照变量。R提供了五种创建对照变量的内置方法（见表9-6），你也可以自己重新创建（此处不做介绍）。默认情况下，对照处理用于无序因子，正交多项式用于有序因子。 表9-6　内置对照

对照变量创建方法 描　　述
contr.helmert 第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对照前三个的均值，以此类推
contr.poly 基于正交多项式的对照，用于趋势分析（线性、二次、三次等）和等距水平的有序因子
contr.sum 对照变量之和限制为0。也称作偏差找对，对各水平的均值与所有水平的均值进行比较
contr.treatment 各水平对照基线水平（默认第一个水平）。也称作虚拟编码
contr.SAS 类似于contr.treatment，只是基线水平变成了最后一个水平。生成的系数类似于大部分SAS过程中使用的对照变量

以对照（treatment contrast）为例，因子的第一个水平变成了参考组，随后的变量都以它为标准。可以通过contrasts()函数查看它的编码过程。

  1. > contrasts(cholesterol$trt)
  2.        2times 4times drugD drugE
  3. 1time       0      0     0     0
  4. 2times      1      0     0     0
  5. 4times      0      1     0     0
  6. drugD       0      0     1     0
  7. drugE       0      0     0     1

若患者处于drugD条件下，变量drugD等于1，其他变量2times、4times和drugE都等于0。无需列出第一组的变量值，因为其他四个变量都为0，这已经说明患者处于1time条件。

在代码清单9-11中，变量trt2times表示水平1time和2times的一个对照。类似地，trt4times是1time和4times的一个对照，其余以此类推。从输出的概率值来看，各药物条件与第一组（1time）显著不同。

通过设定contrasts选项，可以修改lm()中默认的对照方法。例如，若使用Helmert对照：

  1. fit.lm <- lm(response ~ trt, data=cholesterol, contrasts="contr.helmert")

还能通过options()函数修改R会话中的默认对照方法，例如，

  1. options(contrasts = c("contr.SAS", "contr.helmert"))

将设定无序因子默认对比方法为contr.SAS，有序因子的为contr.helmert。虽然我们一直都是在线性模型范围中讨论对照方法的使用，但是在R中，你完全可以将其应用到其他模型中，包括第13章将会介绍的广义线性模型。