简介

LRA系列的第四篇，介绍了 (单因子) 协方差分析的内容，并通过结合SAS和R的结果，帮助读者深入地理解协方差分析的实质。

手敲\(\LaTeX\)难免出现纰漏，有任何疑似错误或者不清楚的地方请直接在下方评论区留言，谢谢各位读者。

这是在这学期修统计辅修的线性回归分析 (Linear Regression Analysis)后整理的笔记，主要是 (单因子) 协方差分析的内容，组织的方式是通过结合SAS和R的结果，深入地理解协方差分析的实质。

本文中用到的 SAS代码 R代码 数据文件 SAS结果 均可直接点击下载。其中SAS结果是利用SAS Studio 得到的。

单因子协方差分析的假设检验

我们假设读者已经对单变量线性回归模型

\[Y_i = \alpha + \beta X_i + \varepsilon_i\]

和单因子方差分析模型

\[Y_{ij} = \mu + \tau_i + \varepsilon_{ij}\]

有了基本的认识，那么单因子协方差分析就是二者的简单组合，即

\[Y_{ij}=\mu+\tau_i+\gamma X_{ij} + \varepsilon_{ij}\]

这一简单模型也是十分有用的，一个简单的例子就是研究植物在不同环境下的生长情况，那么\(X_{ij}\) 为第\(i\)组第\(j\)棵植物实验前的高度，\(Y_{ij}\) 为第\(i\)组第\(j\)棵植物实验后的高度，\(\tau_i\)则代表着组间的平均差距，而\(\mu\)通常为基准。值得注意的是，这样的模型相当于假设了不同组下的回归结果有着相同的斜率(均为\(\gamma\))而有着不同的截距(由\(\tau\)操控)。值得注意的是，如果有\(a\)个水平(即\(\tau_1,\ldots,\tau_a\))，其内部存在着约束。该约束可以定义为\(\sum_i\tau_i=0\)，因为若非零的话非零部分可被吸收进\(\mu\)；亦可定义为\(\tau_a=0\)，即选取第\(a\)组(当然你也可以选择其他组)为基准，此时\(\mu\)代表着这组的平均水平，而\(\tau\)代表着其他组的平均水平与该组的差别。

组间效应的检验

一个很自然的检验自然是，组间效应是否显著到我们需要将其囊括至模型中？事实上，这个检验和简单的单因子模型中的检验并无二致，简介如下。(出于减少记号带来的混乱，我们将因子变量写作A) 检验问题与备择假设

\[H_0: \sum_i \tau_i^2 = 0\quad \mathrm{v.s.}\quad H_1: \mathrm{o.w.} \]

原假设成立下的模型

\[Y_{ij} = \mu + \gamma X_{ij} + \varepsilon_{ij}\]

同样的，我们的检验量\(F^\star\) ，正如我们常使用的检验标准，为“两个模型\(SSE\)的差”在“两个模型自由度的差”上的平均，除以加入新变量后的模型的\(SSE\)在该模型自由度下的平均(即\(MSE\))。

\[F^\star=\frac{SSE(X)-SSE(A,X)}{(n-2)-(n-(a+1))}\div\frac{SSE(A,X)}{n-(a+1)}=\frac{SSR(A|X)}{a-1}\div\frac{SSE(A,X)}{n-a-1}\sim F_{a-1,n-a-1}\]

其中\(n\)为样本数，\(r\)为因子的水平数。线性回归那部分减少了两个自由度，\(r\)个水平的因子减少了\(r-1\)个自由度，因为这\(r\)个水平内部存在约束。

关于该统计量一个简单而直观的理解是，当我们新引入的变量对于减小\(MSE\)毫无作用时，原模型的\(SSE\)在两份自由度上的分配(新引入的变量的自由度\(r-1\)，和剩余的自由度\(n-r-1\))应当是是平均的。因此前者比后者接近于1。若该比值显著大于1，则说明原始模型的\(SSE\)更加集中在这\(r-1\)个自由度上，因此我们有必要将其囊括进模型而获得更好的性能。

截距效应的检验

分析过程与上述完全对偶，读者可以作为练习检验自己。

\[H_0: \beta = 0\quad \mathrm{v.s.}\quad H_1: \beta\neq 0 \]

\[F^\star=\frac{SSR(X|A)}{1}\div\frac{SSE(A,X)}{n-a-1}\sim F_{1,n-a-1}\]

是否同截距的检验

此时，扩充的模型如下，同样需要注意\(\beta_i\)内部存在约束。

\[Y_{ij} = \mu + \tau_i + \gamma x_{ij} + \beta_i x_{ij} + \varepsilon_{ij}\]

此时该模型即为对于不同组数据分别拟合不同的一元线性回归模型，然后使用因子调整截距和斜率将其“整合”在一起。而另一种理解是该模型引入了交叉项\(A*X\)。

检验问题与备择假设

\[H_0: \sum_i \beta_i^2 = 0\quad \mathrm{v.s.}\quad H_1: \mathrm{o.w.} \]

原假设成立下的模型

\[Y_{ij} = \mu + \tau_i + \gamma X_{ij} + \varepsilon_{ij}\]

\[F^{\ast}=\frac{SSR(A\ast X|A,X)}{a-1} \div \frac{SSE(A,X,A\ast X)}{n-2a}\sim F_{a-1,n-2a}\]

分析过程与上仍然一样，同样作为练习，注意两个模型的自由度的关系即可准确写出。

其他形式的检验

诸如\(H_0: \tau_1=\tau_2\)这样的假设同样也是可以被检验。利用contrast的方法会在下面介绍。当然直接通过对自变量做相应的变换，使得原假设的形式变为\(H_0: \tau^\prime_1=0\)也是可行的。

单因子协方差分析的代码

使用SAS可以用极少的代码获得相当多而全的统计结果，但我始终觉得为了知其然更知其所以然，了解SAS里变量和含义及其计算方法是检验自己理解，加深认识的一个十分有效的手段。在这里我们使用基础R尝试重复出SAS给出的结果。

回归结果

上图为SAS回归的结果，从TRT5那一行可以看出SAS取其为基准，也就是\(tau_5=0\)，我们在R里重复这一过程。

1
2
3
4
5
6
7
8
9
10
11
12
oyster <- read.table('oyster.txt')
names(oyster) <- c('TRT','IDX','INIT','NOW')
oyster$TRT <- factor(oyster$TRT)
summary(lm(NOW~INIT+TRT, oyster))
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.25040    1.44308   1.559 0.141205    
## INIT         1.08318    0.04762  22.746 1.87e-12 ***
## TRT2        -0.03581    0.40723  -0.088 0.931169    
## TRT3         1.89922    0.45802   4.147 0.000988 ***
## TRT4         1.35157    0.41937   3.223 0.006135 ** 
## TRT5         0.24446    0.57658   0.424 0.678022

可以发现得到的结果并不一致，但INIT的结果一样还是给了我们一剂强心剂。自己看看R的结果，我们不难猜出R是选取了TRT1作为基准，因此得到的截距以及TRT的系数都不一样。是不是这样呢？我们可以手动算一下TRT2的组内平均水平。在SAS中约为\(2.495-0.280=2.215\)，R中则为\(2.250-0.036=2.214\)，可见结果基本一致。

事实上，我们只需要给TRT重新排序即可得到一模一样的结果。

1
2
3
4
5
6
7
8
9
10
11
oyster$TRT <- factor(oyster$TRT, levels=c(5,1,2,3,4))
m.lm <- lm(FINAL~TRT+INIT, oyster)
summary(m.lm)
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.49486    1.02786   2.427  0.02930 *  
## INIT         1.08318    0.04762  22.746 1.87e-12 ***
## TRT1        -0.24446    0.57658  -0.424  0.67802    
## TRT2        -0.28027    0.49291  -0.569  0.57863    
## TRT3         1.65476    0.42943   3.853  0.00176 ** 
## TRT4         1.10711    0.47175   2.347  0.03417 *

顺便一提，在R中使用anova和drop1求取I型误差和III型误差（事实上drop1并不是用来求III型误差的，当drop1不能给出所想要的III型误差时，使用car::Anova函数可以求多种SS。

1
2
3
4
5
6
7
8
9
10
11
12
13
anova(m.lm)
## Response: FINAL
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## TRT        4 198.407  49.602  164.47 1.340e-11 ***
## INIT       1 156.040 156.040  517.38 1.867e-12 ***
## Residuals 14   4.222   0.302 
drop1(m.lm)
## Model:
## FINAL ~ TRT + INIT
##        Df Sum of Sq     RSS     AIC
## <none>                4.222 -19.107
## TRT     4    12.089  16.312  -0.077
## INIT    1   156.040 160.262  51.622

更深入的理解

事情好像结束了？且慢，这张表里的均值和方差又是如何得到的呢？

我们知道，方差分析和多元线性回归有着千丝万缕的联系。事实上这些系数都是通过多元线性回归的算法(其实就是矩阵运算)得到的。具体说来，要复现出上面的结果，我们需要构造一组输入并拟合多元线性回归模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
X<-cbind(rep(c(0,1,0),c(0,4,16)),rep(c(0,1,0),c(4,4,12)),
         rep(c(0,1,0),c(8,4,8)), rep(c(0,1,0),c(12,4,4)))
X 
##       [,1] [,2] [,3] [,4]
##  [1,]    1    0    0    0
##  [2,]    1    0    0    0
##  [3,]    1    0    0    0
##  [4,]    1    0    0    0
##  [5,]    0    1    0    0
##  [6,]    0    1    0    0
##  [7,]    0    1    0    0
##  [8,]    0    1    0    0
##  [9,]    0    0    1    0
## [10,]    0    0    1    0
## [11,]    0    0    1    0
## [12,]    0    0    1    0
## [13,]    0    0    0    1
## [14,]    0    0    0    1
## [15,]    0    0    0    1
## [16,]    0    0    0    1
## [17,]    0    0    0    0
## [18,]    0    0    0    0
## [19,]    0    0    0    0
## [20,]    0    0    0    0
m.lm <- lm(oyster$FINAL~oyster$INIT+X)
summary(m.lm)
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.49486    1.02786   2.427  0.02930 *  
## oyster$INIT  1.08318    0.04762  22.746 1.87e-12 ***
## X1          -0.24446    0.57658  -0.424  0.67802    
## X2          -0.28027    0.49291  -0.569  0.57863    
## X3           1.65476    0.42943   3.853  0.00176 ** 
## X4           1.10711    0.47175   2.347  0.03417 *

由于X的后4行的所有列都为零，因此这个模型的截距也就是TRT5的水平均值。而X的1-4列的系数则分别代表着对应组的水平均值相对于TRT5的差距(因为截距是共享的)。因此得到一模一样的结果是毫不奇怪的。

是的，我们“手算”出了系数的均值和方差(这里调用的lm和你直接用矩阵乘法并无不同)。然而值得注意的是，这里\(X_1,\ldots,X_4\)的标准差，也就是\(\tau_1,\ldots,\tau_4\)的标准差，也就是截距间的差值的标准差\(s\{b_1-b_5\},\ldots,s\{b_2-b_5\}\)。 那我们能否得到\(s\{b_1\},\ldots,s\{b_5\}\)的值呢？答案当然是可以的，我们同样给TRT5附上自变量即可，但要注意，此时不再需要截距项了(想想为什么？)。

1
2
3
4
5
6
7
8
9
10
11
12
XX<-cbind(X,rep(c(0,1),c(16,4)))
tmpX<-cbind(XX,oyster$INIT)
colnames(tmpX) <- c(paste0('TRT', c(1,2,3,4,5)), 'INIT')
summary(lm(oyster$FINAL~tmpX-1))
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)    
## tmpXTRT1  2.25040    1.44308   1.559  0.14120    
## tmpXTRT2  2.21459    1.32290   1.674  0.11631    
## tmpXTRT3  4.14962    1.20553   3.442  0.00397 ** 
## tmpXTRT4  3.60197    1.28798   2.797  0.01428 *  
## tmpXTRT5  2.49486    1.02786   2.427  0.02930 *  
## tmpXINIT  1.08318    0.04762  22.746 1.87e-12 ***

从均值可以看到我们确实消去了截距项(而将截距项加入每个TRT中)，此时的标准差即为我们感兴趣的\(s\{b_1\},...,s\{b_5\}\)。为了保险起见我们可以验算一下，注意\(s^2\{b_i-b_j\}=s^2\{b_i\}-2s^2\{b_i,b_j\}+s^2\{b_j\}\)。但从哪里获得\(s^2\{b_i,b_j\}\)呢？最简单的方式自然是矩阵了。

1
2
3
4
5
6
SSE <- sum(m.lm$residuals^2)
MSE <- SSE/14
B <- solve(t(tmpX)%*%tmpX)%*%t(tmpX)%*%oyster$FINAL
COV <- MSE*solve(t(tmpX)%*%(tmpX))
sqrt(COV[1,1]+COV[5,5]-2*COV[1,5])
## [1] 0.576582

得到的正是SAS中TRT1的标准差。

预测其他变量

有了上述的数据，算其线性组合的均值和标准差岂非小菜一碟？

1
2
3
4
5
6
7
8
9
TRT1.coef <- c(4, -1, -1, -1, -1) / 5
E.TRT1 <- sum(TRT1.coef * B[1:5])
S.TRT1 <- sqrt(sum(TRT1.coef %*% t(TRT1.coef) * COV[1:5,1:5]))
E.TRT1
## [1] -0.6918875
S.TRT1
## [1] 0.3105175
E.TRT1 / S.TRT1
## -2.228175

Contrast

对于TRT1来说，一种做法就是利用上一小节的结果。对于自由度为1的变量，其t检验和F检验是等价的，有关系式\(F=t^2\)(为什么？想想t分布的F分布的定义，想想多元回归分析的假设检验)。因此上式的 \((-2.228175)^2\) 正为SAS contrast 结果中的F值。我希望读者对这一结果并不感到意外。

当然，使用上述的检验方差分配的方式也能获得同样的结果。对于\((4, -1, -1, -1, -1)\)这个contrast，直观的理解就是检验水平1是否和其他四个水平的均值相同。因此对于其他四个水平，有其各自的均值，而水平1的均值则建立在其他四个水平的均值上。反应在自变量上则应该如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
TRT1.X <- XX[,2:5]
TRT1.X[1:4,]=1/4
TRT1.X
##      [,1] [,2] [,3] [,4]
##  [1,] 0.25 0.25 0.25 0.25
##  [2,] 0.25 0.25 0.25 0.25
##  [3,] 0.25 0.25 0.25 0.25
##  [4,] 0.25 0.25 0.25 0.25
##  [5,] 1.00 0.00 0.00 0.00
##  [6,] 1.00 0.00 0.00 0.00
##  [7,] 1.00 0.00 0.00 0.00
##  [8,] 1.00 0.00 0.00 0.00
##  [9,] 0.00 1.00 0.00 0.00
## [10,] 0.00 1.00 0.00 0.00
## [11,] 0.00 1.00 0.00 0.00
## [12,] 0.00 1.00 0.00 0.00
## [13,] 0.00 0.00 1.00 0.00
## [14,] 0.00 0.00 1.00 0.00
## [15,] 0.00 0.00 1.00 0.00
## [16,] 0.00 0.00 1.00 0.00
## [17,] 0.00 0.00 0.00 1.00
## [18,] 0.00 0.00 0.00 1.00
## [19,] 0.00 0.00 0.00 1.00
## [20,] 0.00 0.00 0.00 1.00

对于这个模型，我们拟合后计算其\(SSE\)，并与之前的做差，即能得到SAS中的“对比组平方和”

1
2
3
SSE.1 <- sum(lm(oyster$FINAL~oyster$INIT+TRT1.X-1)$residuals^2)
SSE.1 - SSE
## [1] 1.497346

2LVL检验了什么呢？根据其constrast的形式不难发现其所检验的问题为 \[H_0: \tau_1-\tau_2=0\quad\mathrm{and}\quad\tau_1-\tau_5=0\quad\mathrm{and}\quad\tau_3-\tau_4=0 \quad \mathrm{v.s.} \quad H_1: \mathrm{o.w.}\] 也就是是否可以将其视为只有两个水平。 当明确了这个问题后，列出缩减后模型的形式，再计算。其代码如下。但我希望读者在这里先听一下，思考下现在的自变量\(X\)应该是什么形式

1
2
3
4
LVL.X <- cbind(rep(c(0,1,0),c(8,8,4)),rep(c(1,0,1),c(8,8,4)))
SSE.2 <- sum(lm(oyster$FINAL~oyster$INIT+LVL.X-1)$residuals^2)
SSE.2 - SSE
## [1] 0.6219307

而关于SAS的最后一页（有交叉项的模型），分析过程与上述大同小异，就不再赘述了。

总结

本文通过在R里复现SAS的结果，并尝试用多元回归的方式求解相应的统计量，旨在加深读者对方差分析的理解和对方差分析和回归分析内部本质的认识，希望读者有所收获。