专注系列化、高质量的R语言教程

推文索引 | 联系小编 | 付费合集

当存在两个及以上的分组变量时，可以使用多因素方差分析（N-way  ANOVA、Multifactor ANOVA）检验各组的样本均值是否存在显著差异。本篇主要以双因素方差分析（Two-way  ANOVA）为例介绍相关内容。

本篇的目录如下：

• 2 多因素方差分析

  • 2.1 示例数据

  • 2.2 平衡试验设计

  • 2.3 I型方差分析

  • 2.4 交互效应

  • 2.5 多因素方差分析

  • 未完待续

2 多因素方差分析

2.1 示例数据

本篇使用两个示例数据。

第一个示例数据是来自基础包datasets的npk：

npk
##    block N P K yield
## 1      1 0 1 1  49.5
## 2      1 1 1 0  62.8
## 3      1 0 0 0  46.8
## 4      1 1 0 1  57.0
## 5      2 1 0 0  59.8
## 6      2 1 1 1  58.5
## 7      2 0 0 1  55.5
## ...

npk共有5个变量：yield表示产量，为观测值；block为土地编号，取值为1-6，因子类型；N、P、K分别表示氮、磷、钾的使用情况，均为二分变量。

第二个示例数据是由mtcars数据集的四个变量组成，记为mpg：

library(dplyr)
mpg <- mtcars %>%select(mpg, cyl, gear, carb) %>%mutate(cyl = factor(cyl),gear = factor(gear),carb = factor(carb))
mpg
##                      mpg cyl gear carb
## Mazda RX4           21.0   6    4    4
## Mazda RX4 Wag       21.0   6    4    4
## Datsun 710          22.8   4    4    1
## Hornet 4 Drive      21.4   6    3    1
## Hornet Sportabout   18.7   8    3    2
## Valiant             18.1   6    3    1
## Duster 360          14.3   8    3    4
## ...

• 以mpg为观测变量，cyl、gear、gear为分组变量，它们分别将样本分为3、3、6组。

• 当分组数量在三个及以上时，分组变量应为因子类型。

2.2 平衡试验设计

参照单因素方差分析的模型形式，可以很自然地推测双因素方差分析的模型形式。

例2.1：观察并比较如下两组模型的结果

## npk
aov(yield ~ N + P, npk)
##                        N        P Residuals
## Sum of Squares  189.2817   8.4017  678.6817
## Deg. of Freedom        1        1        21
aov(yield ~ P + N, npk)
## Terms:
##                        P        N Residuals
## Sum of Squares    8.4017 189.2817  678.6817
## Deg. of Freedom        1        1        21## mpg
aov(mpg ~ cyl + gear, mpg)
## Terms:
##                      cyl     gear Residuals
## Sum of Squares  824.7846   8.2519  293.0107
## Deg. of Freedom        2        2        27
aov(mpg ~ gear + cyl, mpg)
## Terms:
##                     gear      cyl Residuals
## Sum of Squares  483.2432 349.7933  293.0107
## Deg. of Freedom        2        2        27

上述两组模型各包含两个方差分析，模型形式的区别仅在于分组变量的顺序不同：

• 对于第一组模型（npk）来说，两个方差分析的结果完全相同；

• 对于第二组模型（mpg）来说，cyl和gear变量所对应的离差平方和在两个方差分析中是不同的，但二者之和仍然相同（824.7846 + 8.2519 = 483.2432 + 349.7933 = 833.0365），自由度和组内离差平方和（Residuals）也相同。

以上结果意味着多因素方差分析可能与分组变量顺序有关。实际上，这与样本分组是否符合平衡试验设计（balanced design）有关。平衡设计是指，由多个分组变量形成的交叉组所对应的样本数量完全一致。

对于平衡设计来说，多因素方差分析的结果与变量顺序无关；而对于非平衡设计来说，存在三种类型的多因素方差分析：I型、II型和III型，而aov()函数使用的是I型方差分析，它的结果与变量顺序有关。

在aov()函数的帮助文档中作者写道：aov()函数是为平衡设计准备的，对于非平衡设计它的结果很难解释。

统计两个示例数据的交叉组样本数量：

with(npk, {tapply(yield, list(N, P), length)})
##   0 1
## 0 6 6
## 1 6 6with(mpg, {tapply(mpg, list(cyl, gear), length)})
##    3  4 5
## 4  1  8 2
## 6  2  4 1
## 8 12 NA 2

可以看出，npk数据的4个交叉组的样本均为6，属于平衡设计；而mpg数据的9个交叉组的样本并不完全一致，并且有一个交叉组没有样本（NA），属于非平衡设计。

2.3 I型方差分析

对于平衡设计来说，三种类型的方差分析的结果一致，因此下文仅针对非平衡设计的方差分析，示例数据为mpg。

例2.2：aov()函数执行的是I型方差分析，使用它运行如下模型

aov(mpg ~ cyl, mpg)
aov(mpg ~ gear, mpg)
aov(mpg ~ cyl + gear, mpg)
aov(mpg ~ gear + cyl, mpg)

为了方便比较，我们将模型的离差平方和汇总如下：

SS表示离差平方和（Sum of Squares）。

通过比较可以得到如下结论：

• 所有模型的离差平方和之和（即总离差平方和）是相等的；

• 在双因素方差分析中，分组变量的离差平方和与它的顺序有关，其中第一个分组变量的离差平方和与它对应的单因素方差分析的结果相同；

• 在双因素方差分析中，cyl和gear的离差平方和与顺序无关。

总离差平方和（）：

只与观察变量有关，而与分组变量无关，其计算方法和单因素方差分析完全一致。使用表示样本标识，表示全体样本的均值，有

组间离差平方和（和）：

这里组间离差平方和由两部分组成，分别对应两个分组变量的离差平方和，记为、。使用、表示两个分组变量的水平数，、表示分组标识。

假设对应的是第一个分组变量的离差平方和，我们已经知道它和单因素方差分析的结果一样：

其中，表示第组的样本均值。

不易直接计算，但是我们知道它与之和为定值：

因此只要计算出即可反推出。

组内离差平方和（）：

相当于“剩余”部分，在输出结果中也使用Residuals标识。方差分析与线性回归存在内在联系联系，实际就等于线性回归的残差平方和。

fit.lm <- lm(mpg ~ cyl + gear, mpg)
sum(residuals(fit.lm)^2)
## [1] 293.0107

自由度和F统计量：

的自由度为（为样本总量），、的自由度分别为和，的自由度为。

2.4 交互效应

分组变量之间可能存在交互效应。

例2.3：比较不含和含有交互项的双因素方差分析

aov(mpg ~ cyl + gear, mpg)
## Terms:
##                      cyl     gear Residuals
## Sum of Squares  824.7846   8.2519  293.0107
## Deg. of Freedom        2        2        27aov(mpg ~ cyl*gear, mpg)
## Terms:
##                      cyl     gear cyl:gear Residuals
## Sum of Squares  824.7846   8.2519  23.8907  269.1200
## Deg. of Freedom        2        2        3        24

在有交互项的情况下，组间离差平方和由三部分组成：、、，其中表示交互项的离差平方和。

可以看出，交互项的出现不影响、、的值。那的值和自由度应如何确定呢？

我们依然可以先使用线性回归求出再反推：

fit2.lm <- lm(mpg ~ cyl*gear, mpg)
sum(residuals(fit2.lm)^2)
## [1] 269.12

不过也可以直接求出、、的和：

其中，表示交叉组的样本均值。

将交叉组看作是由单个变量分组形成的，也可以使用类似单因素方差分析的方法直接计算：

假设交叉组的个数为，则、、的自由度之和为，因此的自由度为。

显然，交叉组的最大数目为，在此情况下：

前文已经提到，mpg数据中有一个交叉组不存在，因此，而对于平衡设计来说，总是成立的。

例2.4：手动计算和

data <- mpg %>%mutate(mean = mean(mpg)) %>%group_by(cyl, gear) %>%mutate(mean2 = mean(mpg)) %>%ungroup() %>%mutate(SSAB = (mean2 - mean)^2,SSE = (mpg - mean2)^2)sum(data$SSAB) - 824.7846 - 8.2519
## [1] 23.89069
sum(data$SSE)
## [1] 269.12

2.5 多因素方差分析

从前文可以看出，即使只有两个分组变量，I型方差分析的都很难直接计算。那对于多因素方差分析来说该如何计算各个离差平方和呢？下面通过两个例子进行说明。

例2.5：运行形式最简单的三因素方差分析

aov(mpg ~ cyl + gear + carb, mpg)
## Terms:
##                      cyl     gear     carb Residuals
## Sum of Squares  824.7846   8.2519  88.5199  204.4908
## Deg. of Freedom        2        2        5        22

cyl的离差平方和：

fit1 <- lm(mpg ~ cyl, mpg)
sum((fitted(fit1) - mean(mpg$mpg))^2)
## [1] 824.7846

gear的离差平方和：

fit2 <- lm(mpg ~ cyl + gear, mpg)
sum((fitted(fit2) - mean(mpg$mpg))^2) - sum((fitted(fit1) - mean(mpg$mpg))^2)
## [1] 8.251855

carb的离差平方和：

fit3 <- lm(mpg ~ cyl + gear + carb, mpg)
sum((fitted(fit3) - mean(mpg$mpg))^2) - sum((fitted(fit2) - mean(mpg$mpg))^2)
## [1] 88.5199

Residuals的离差平方和：

sum(residuals(fit3)^2)
## [1] 204.4908

因此，I型方差分析可以使用逐步回归的方法计算离差平方和。I型平方和又称为顺序平方和。

例2.6：运行含有交互项的三因素方差分析

aov(mpg ~ cyl*gear + carb, mpg)
## Terms:
##                      cyl     gear     carb cyl:gear Residuals
## Sum of Squares  824.7846   8.2519  88.5199  25.3878  179.1030
## Deg. of Freedom        2        2        5        2        20aov(terms(mpg ~ cyl*gear + carb, keep.order = T), mpg)
## Terms:
##                      cyl     gear cyl:gear     carb Residuals
## Sum of Squares  824.7846   8.2519  23.8907  90.0170  179.1030
## Deg. of Freedom        2        2        3        4        20

默认情况下，交互项会排在所有主效应变量之后；而使用terms()函数可以更改设置。位置不同会导致离差平方和发生变化。

两种情况下cyl:gear的离差平方和：

## 默认情况
fit11 <- lm(mpg ~ cyl + gear + carb, mpg)
fit22 <- lm(mpg ~ cyl*gear + carb, mpg)
sum((fitted(fit22) - mean(mpg$mpg))^2) - sum((fitted(fit11) - mean(mpg$mpg))^2)
## [1] 25.38784## 更改后的情况
fit33 <- lm(mpg ~ cyl + gear, mpg)
fit44 <- lm(mpg ~ cyl*gear, mpg)
sum((fitted(fit44) - mean(mpg$mpg))^2) - sum((fitted(fit33) - mean(mpg$mpg))^2)
## [1] 23.89074

正如aov()函数的帮助文档所言，非平衡设计的I型方差分析的结果很难解释。学堂君注意到<医学方>公众号的一篇推文对其做了非常直观的解释，有兴趣的读者可以阅读：

关于方差分析，您必须知道的四种类型平方和