R 语言方差分析：aov()、ANOVA 表与 Tukey HSD

为什么用方差分析，而不是一堆 t 检验

t 检验比较的是两个均值。有三组或更多组时，最容易想到的做法是对每一对都做 t 检验，而这恰恰是个陷阱。每次检验在 0.05 水平下都有 5% 的假阳性风险，这些风险会不断累积：4 组需要做 6 次两两检验，至少出现一次虚假“显著”结果的概率约为 26%；5 组（10 次检验）则约为 40%。你会凭空从纯粹的噪声中“发现”结论。

方差分析（ANOVA，ANalysis Of VAriance）的解决办法是：用一个 p 值回答一个问题：所有组的均值都相同，还是至少有一组不同？别看名字叫方差分析，它比较的其实是均值，只不过是通过分析方差来比较：如果各组平均值之间的分散程度，超出了组内噪声所能解释的范围，那就说明确实存在差异。

用 aov() 做单因素方差分析

PlantGrowth 正是为此准备的数据集：植物在一个对照条件和两个处理条件下的重量。用 aov() 拟合，然后（这一步很重要）用 summary() 打印表格：

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

把公式读作“weight 是否取决于 group？”。分组变量必须是因子：PlantGrowth$group 本身就是因子，但如果你的组是用数字编码的（剂量、批次 ID），就要把它包起来：aov(y ~ factor(dose), ...)。否则 aov() 会悄悄地对组编码拟合一条回归直线，而不是比较各组均值：模型错了，却没有任何报错。

读懂方差分析表

表格有两行（因子和残差）、五列。逐格来看：

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df：自由度。因子行是“组数 − 1”（3 组 → 2）；残差行是“观测数 − 组数”（30 − 3 = 27）。可以快速确认 R 理解的实验设计和你想的一致。
• Sum Sq：变异，分成两部分。group 行是组间部分：各组均值离总均值有多远。Residuals 是组内部分：植物围绕各自组均值的波动有多大。两者相加就是数据中的总变异。
• Mean Sq：每个 Sum Sq 除以对应的 Df，把变异总量变成可以相互比较的“每个自由度的变异”。残差的 Mean Sq（0.389）就是噪声水平。
• F value：比值，即 group 的 Mean Sq 除以 Residuals 的 Mean Sq（1.8832 / 0.3886 ≈ 4.85）。如果各组均值真的全都相等，这个比值会在 1 附近徘徊。F 越大，说明组间差异越超出噪声所能解释的范围。
• Pr(>F)：p 值，即假如三组的真实均值完全相同，得到这么大的 F 的概率。这里是 0.016，在常用的 0.05 水平下足够小，可以拒绝“全部相等”。和往常一样，它不是原假设为真的概率，也不会告诉你是哪几组不同、差多少。

最后这一点正是方差分析的关键局限：F 显著只说明“某处存在某种差异”。仅此而已。

哪些组不同？TukeyHSD()

接下来的问题需要事后检验。Tukey 的 HSD（Honest Significant Difference）检验会对每一对组做检验，同时把所有比较的族系错误率整体控制在 5%：

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

每一对组占一行，每行四个数：

• diff：估计的均值差（后一个组减去前一个组）。
• lwr, upr：这个差值的 95% 族系置信区间。
• p adj：已经针对三次比较做过校正的 p 值。

对于 PlantGrowth：trt2-trt1 的差值约为 0.87，p adj ≈ 0.012，区间不包含零：处理 2 的植物长得比处理 1 好。两个处理组对比对照组的行（trt1-ctrl、trt2-ctrl）区间都跨过了零，p adj 远高于 0.05：根据这份样本，两种处理都无法与对照组区分开。这条经验法则和置信区间的用法完全一致：区间不包含零 ⇔ p adj 低于 0.05。

注意，整体 F 检验说的是“存在差异”，而 Tukey 检验把差异精确定位到了某一对组上。这种两步结构正是整套方法的精髓：先做一次可靠的整体检验，再在正确校正的前提下逐对排查。

双因素方差分析：两个因子及其交互作用

有两个分组变量时，一次调用就能同时检验两者，还能检验它们之间有没有交互作用。ToothGrowth 把补充剂类型（supp）和剂量（0.5、1、2 mg，是数值，所以需要 factor()）交叉组合：

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) 会展开成三个效应，每个效应在表中占一行：

• supp：对各剂量取平均后，补充剂类型有没有影响？（有：p ≈ 0.0002。）
• factor(dose)：对各补充剂取平均后，剂量有没有影响？（非常显著：p 极小。）
• supp:factor(dose)：交互作用，即补充剂的效果是否随剂量而变？这里 p ≈ 0.022，是的。低剂量时橙汁比抗坏血酸效果好，但到了 2 mg 剂量，差距就消失了。

显著的交互作用相当于给主效应贴上了一张警示标签：“补充剂类型有影响”只在平均意义上成立，而这个平均值掩盖了随剂量变化的情况。当交互作用显著时，应该描述各个组合的情况（用 TukeyHSD(fit2, "supp:factor(dose)")，或者通过分组汇总得到各组均值表），而不是只报告主效应。只有当你有意要拟合一个没有交互作用的模型时，才用 + 代替 *。

前提假设，以及基于秩的替代方法

方差分析的数学基础依赖三个假设，按从最不能通融到最可以通融排列：

• 独立性：观测值之间互不影响。一旦违反，事后没有任何办法补救；这是研究设计本身的属性。
• 各组方差相等：残差的 Mean Sq 是一个合并后的噪声估计，所以各组的噪声大小应该差不多。一行代码就能检查：bartlett.test(weight ~ group, data = PlantGrowth)（p 值大说明没有证据表明方差不等）。
• 残差近似正态：在各组样本量均衡且足够大时影响最小；用 residuals(fit) 取出残差，画个直方图或箱线图看一看。

当数据严重偏斜、是有序数据或者有很多离群值时，单因素方差分析基于秩的替代方法是 kruskal.test(weight ~ group, data = PlantGrowth)，它相当于 wilcox.test() 的多组版本。只需一行代码，用一些检验效能换来了稳健性。

要点回顾

• 方差分析用一个 p 值回答关于 3 组及以上均值的一个问题，解决了多次 t 检验的多重比较问题。
• fit <- aov(y ~ group, data = df); summary(fit)，而且分组变量必须是因子。
• F 值是组间信号与组内噪声之比；Pr(>F) 小意味着“某处存在某种差异”，仅此而已。
• TukeyHSD(fit) 找出具体是哪几对组不同，并且内置了多重比较校正。
• a * b 拟合两个因子及其交互作用；交互作用显著时，单看主效应不能说明全部情况。
• 检查方差齐性（bartlett.test）和残差正态性；kruskal.test() 是基于秩的替代方法。

下一篇：从比较组均值转向为变量关系建模，用 lm() 做线性回归。