35 R方差分析
35.1 单因素方差分析
单因素方差分析可以看成基础统计中两样本t检验的一个推广, 要比较试验观测值的某个因变量(称为“指标”)按照一个分组变量(称为“因素”)分组后, 各组的因变量均值有无显著差异。
设因素\(A\)将所有观测分为\(m\)个组, 每组对因变量进行\(r\)次观测, 且各次观测相互独立, 模型为 \[\begin{aligned} y_{ij} =& \mu_i + e_{ij}, i=1,2,\dots, m, \ j=1,2, \dots, r. \\ e_{ij} & \text{ iid N}(0, \sigma^2) . \end{aligned}\]
希望研究如下问题: \[\begin{aligned} H_0:\ \mu_1 = \mu_2 = \dots = \mu_m . \end{aligned}\]
如果拒绝\(H_0\), 希望找出哪些组的均值两两之间是有显著差异的。
这样的问题称为单因素方差分析问题, “单因素”是指按照一种分组方式分组。 主要模型中假定各个观测独立, 误差\(e_{ij}\)服从正态分布, 且误差的方差相等。
方差分析在最初是采用平方和分解的方式进行检验, 但是随着模型变得更复杂, 现在处理方差分析问题一般都采用线性模型的一般理论进行处理。
单因素方差分析问题可以推广到不同组的重复观测次数不同的情形。
35.1.1 单因素方差分析示例
例子:考虑如下的方差分析问题。 某工厂要比较三种不同组装工艺(A, B, C)的工作效率, 将15名工人随机分为3个组, 每个组采用一种工艺, 一周后各组的成品数如下:
A | B | C |
---|---|---|
58 | 58 | 48 |
64 | 69 | 57 |
55 | 71 | 59 |
66 | 64 | 47 |
67 | 68 | 49 |
我们需要长表格式,进行转换:
<- d.manu.w |>
d.manu pivot_longer(A:C,
names_to = "grp",
values_to = "y") |>
mutate(grp = factor(grp, levels=c("A", "B", "C"))) |>
arrange(grp)
::kable(d.manu) knitr
grp | y |
---|---|
A | 58 |
A | 64 |
A | 55 |
A | 66 |
A | 67 |
B | 58 |
B | 69 |
B | 71 |
B | 64 |
B | 68 |
C | 48 |
C | 57 |
C | 59 |
C | 47 |
C | 49 |
用aov()
进行方差分析:
<- aov(y ~ grp, data=d.manu)
aov.manu summary(aov.manu)
## Df Sum Sq Mean Sq F value Pr(>F)
## grp 2 520 260.00 9.176 0.00382 **
## Residuals 12 340 28.33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
主效应grp
(分组)的F检验的p值为0.004,
若检验水平为0.05则分组效应显著,
各组之间有显著差异。
可以用并列盒形图比较各组的取值, 并检查分布对称性和方差齐性:
ggplot(data = d.manu, mapping = aes(
x = grp, y = y)) +
geom_boxplot()
从图形看,C组的工作效率显著低于A组和B组, A、B两组之间差异不明显。
为了检查各组方差是否相等,可以进行Bartlett检验:
bartlett.test(y ~ grp, data = d.manu)
##
## Bartlett test of homogeneity of variances
##
## data: y by grp
## Bartlett's K-squared = 0.024476, df = 2, p-value = 0.9878
结果表明可以认为各组的方差相等。
除了各组并列盒形图, 也可以做点图并将各组的均值、均值置信区间作图,如:
<- d.manu |>
d2.manu group_by(grp) |>
summarise(ym = mean(y),
ystd = sd(y),
yrad = 1.96 * ystd / sqrt(n()))
ggplot(data=d2.manu) +
geom_errorbar(data = d2.manu,
mapping = aes(
x = grp, ymin = ym - yrad, ymax = ym + yrad )) +
geom_point(mapping = aes(
x = grp, y = ym ), size=3, color="turquoise4") +
geom_jitter(data=d.manu, mapping=aes(
x = grp, y = y ), width=0.1)
可以看出C的误差条(均值置信区间)低于A和B的误差条, 而A和B的误差条有交集。
35.1.2 多重比较
为了找到各组两两之间是否有显著差异, 可以进行两两的独立两样本t检验, 但这样不能利用共同的模型参数, 进行多次重复检验也会使得总第一类错误概率变得比较高, 发生过度拟合。
进行多个假设检验(如均值比较)的操作称为“多重比较”(multiple comparison, 或multiple testing), 多次检验会使得总第一类错误概率增大。 以两两比较为例, 当比较次数很多时, 即使所有的组之间都没有真实的差异, 很多个两两比较也会找到差异最大的一对, 使得发现的的显著差异有很大可能性不是真实存在的。
为此,可以进行一些调整,
使得报告的检验p值能够控制总第一类错误概率。
multcomp包的glht()
函数可以对方差分析结果进行多重比较并控制总错误率,
一种方法是利用Tukey的HSD(Honest Significant Difference)方法,
程序如下:
library(multcomp, quietly=TRUE)
glht(aov.manu, linfct = mcp(grp = "Tukey")) |>
summary()
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = y ~ grp, data = d.manu)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 4.000 3.367 1.188 0.48201
## C - A == 0 -10.000 3.367 -2.970 0.02931 *
## C - B == 0 -14.000 3.367 -4.159 0.00373 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Tukey HSD检验的结果是, 在0.05水平下, A和B没有显著差异, C与A、B均有显著差异。
35.1.3 方差不相等情形
方差分析模型要求误差项独立同正态分布N(0, \(\sigma^2\)),
这意味着各组的因变量方差相等。
实际中不同组的因变量可能有不同的方差。
R提供了oneway.test()
函数作为独立两样本t检验的Welch方法的推广,
可以不要求方差相等。
如
oneway.test(y ~ grp, data=d.manu)
##
## One-way analysis of means (not assuming equal variances)
##
## data: y and grp
## F = 8.1976, num df = 2.0000, denom df = 7.9914, p-value = 0.01159
p值为0.011,在0.05水平下显著, 说明三种生产工艺的时间有显著差异。
为了进行多重比较, 可以进行两两t检验并不使用合并的标准差估计, 使用Holm方法进行p值调整以控制总错误率:
with(d.manu,
pairwise.t.test(y, grp, pool.sd=FALSE,
p.adjust.method="holm"))
##
## Pairwise comparisons using t tests with non-pooled SD
##
## data: y and grp
##
## A B
## B 0.258 -
## C 0.039 0.010
##
## P value adjustment method: holm
在0.05水平下, A和B没有显著差异,C与A和B都有显著差异。
35.1.4 非参数方差分析
如果各组的因变量(指标)分布严重偏离正态, 则单因素方差分析所依据的F检验会有很大的误差。 这时, 可以使用非参数方法, Kruskal-Wallis检验是独立两样本比较的Wilcoxon秩和检验的推广。 如:
kruskal.test(y ~ grp, data=d.manu)
##
## Kruskal-Wallis rank sum test
##
## data: y by grp
## Kruskal-Wallis chi-squared = 7.7677, df = 2, p-value = 0.02057
检验p值为0.02,在0.05水平下拒绝零假设, 认为各组之间有显著差异。
35.2 协方差分析
在单因素方差分析问题中, 有可能存在一个连续型的解释变量, 对因变量(指标)有影响, 希望在排除此连续型解释变量影响后对各组进行比较。 这相当于如下的模型: \[\begin{aligned} y_{ij} = \mu_i + \beta x_{ij} + e_{ij}, \end{aligned}\] 其中\(x_{ij}\)为已知值的解释变量。
例子 考虑multcomp包的litter数据。
为研究某种药物对出生小数体重的影响,
将若干怀孕小鼠随机分配到4个不同的处理组,
变量dose
表示这4个组,
其取值是母鼠服药的剂量。
研究的因变量(指标)是母鼠所生产的幼鼠的平均体重(weight),
要排除协变量gesttime
(怀孕时间)的影响。
这个数据的各组试验数不相同,频数统计为:
data(litter, package="multcomp")
|>
litter count(dose) |>
::kable() knitr
dose | n |
---|---|
0 | 20 |
5 | 19 |
50 | 18 |
500 | 17 |
进行协方差分析:
<- aov(weight ~ gesttime + dose,
aov.lit1 data = litter)
summary(aov.lit1)
## Df Sum Sq Mean Sq F value Pr(>F)
## gesttime 1 134.3 134.30 8.049 0.00597 **
## dose 3 137.1 45.71 2.739 0.04988 *
## Residuals 69 1151.3 16.69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
从结果中dose
的检验结果看,
排除协变量gesttime
影响后,
在0.05水平下,
四个处理组之间有显著差异。
将服药的三个组的均值与不服药(dose=0)的组比较的检验如下:
<- rbind("服药与不服药比较" = c(3, -1, -1, -1))
contr glht(aov.lit1, linfct = mcp(dose = contr)) |>
summary()
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = weight ~ gesttime + dose, data = litter)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 服药与不服药比较 == 0 8.284 3.209 2.581 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
检验的零假设为\(H_0: \alpha_1 = (\alpha_2 + \alpha_3 + \alpha_4)/3\),
其中\(\alpha_j\)表示第\(j\)个组的均值,\(\alpha_1\)对应于dose=0
的组,
即不服药的组。
检验结果表明在0.05水平下,
服药的和不服药的指标均值有显著差异。
35.3 两因素方差分析
两因素方差分析, 最简单的情形是有两个分组变量(称为因素)A和B, \(A\)有\(s\)个分类, \(B\)有\(t\)个分类, 对\(A\)和\(B\)的每一种搭配\((i,j)\), 重复试验\(r\)次得到观测值\(y_{ijk}\), 设各次观测值相互独立。 模型为 \[ y_{ijk} = \mu_{ij} + e_{ijk}, \ i=1,\dots,s, \ j=1,\dots,t, \ k=1,\dots, r, \] 误差项\(e_{ijk}\)独立同正态分布N(0,\(\sigma^2\))。
一般将上述模型用另外的参数表示成: \[\begin{aligned} y_{ijk} =& \mu + \alpha_i + \beta_j + \gamma_{ij} + e_{ijk}, \end{aligned}\] 其中\(\alpha_i\)表示因素A的“主效应”, \(\beta_j\)表示因素\(B\)的主效应, \(\gamma_{ij}\)表示因素A和因素B的交互作用效应。 这样模型中的参数是冗余的, 一般会加限制, 比如设\(\alpha_1 = 0\), \(\beta_1 = 0\), \(\gamma_{i1} = \gamma_{1j} = 0\)。
二元方差分析的问题是分别检验两个因素的主效应是否显著(不等于零), 交互作用效应是否存在(不等于零)。
35.3.1 两因素方差分析计算示例
以R的datasets包的ToothGrowth数据为例。
考虑维生素C对豚鼠的成牙质细胞生长的影响,
因变量(指标)len
是某牙齿长度测量值,
考虑两种不同的维生素C类型(变量supp
,取值为OJ或VC),
以及三种剂量(变量dose
, 取值为0.5, 1.0, 2.0)。
需要将这两个因素都转换成R的因子:
data(ToothGrowth, package="datasets")
<- ToothGrowth |>
d.tooth mutate(supp = factor(supp, levels=c("OJ", "VC")),
dose = factor(dose, levels=c(0.5, 1.0, 2.0)))
|>
d.tooth count(supp, dose)
## supp dose n
## 1 OJ 0.5 10
## 2 OJ 1 10
## 3 OJ 2 10
## 4 VC 0.5 10
## 5 VC 1 10
## 6 VC 2 10
可见这个试验是均衡有重复完全试验设计。
进行二元方差分析:
<- aov(len ~ dose + supp + dose:supp,
aov.to1 data = d.tooth)
summary(aov.to1)
## Df Sum Sq Mean Sq F value Pr(>F)
## dose 2 2426.4 1213.2 92.000 < 2e-16 ***
## supp 1 205.4 205.4 15.572 0.000231 ***
## dose:supp 2 108.3 54.2 4.107 0.021860 *
## Residuals 54 712.1 13.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
在0.05水平下, 两个主效应和交互作用效应都显著。
35.3.2 交互作用图
R函数interaction.plot()
可以绘制表示交互作用的图形。
此图形以因变量为y坐标,
以第一个因素的水平为横坐标作折线图,
但根据第二个因素的不同水平分组作图。
如果不同组的折线图基本平行则没有交互作用效应,
否则提示有交互作用效应。
程序如:
with(d.tooth,
interaction.plot(dose, supp, len) )
图形结果提示有交互作用效应。