R方差分析完全指南：aov 单因素、双因素与 TukeyHSD 事后检验
三组以上的均值比较，不能靠反复做 t 检验解决。三组两两比是 3 次检验，每次都用 α = 0.05，至少犯一次第一类错误的概率是 1 − 0.95³ ≈ 14%；五组要比较 10 次，这个概率涨到 40%。方差分析（ANOVA，analysis of variance）的做法是先做一次整体的 F 检验，把「族错误率」（family-wise error rate）控制在 5% 以内，确认存在差异后再做事后多重比较。
多重比较的问题有多大
Section titled “多重比较的问题有多大”上面的 14% 是三次检验相互独立时的上界。三组的两两比较共用同一批数据，检验之间正相关，实际错误率会低于这个上界，但离 5% 仍然很远。用模拟量一次：
set.seed(42)
reject <- replicate(10000, { g <- list(rnorm(20), rnorm(20), rnorm(20)) p <- c( t.test(g[[1]], g[[2]])$p.value, t.test(g[[1]], g[[3]])$p.value, t.test(g[[2]], g[[3]])$p.value ) min(p) < 0.05})
mean(reject)[1] 0.1234三组数据都抽自同一个标准正态分布，真实差异为零，但 10000 次实验里有 1234 次至少有一次 t 检验报出显著。实验整体的错误率是 12.3%，不是 5%。这就是「先做整体 F 检验，再做事后比较」的必要性。
单因素方差分析
Section titled “单因素方差分析”PlantGrowth 是标准的单因素设计：三组处理（ctrl、trt1、trt2），每组 10 株，因变量是干重。
fit <- aov(weight ~ group, data = PlantGrowth)summary(fit) 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---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1五个列的含义：
- Df（自由度，degrees of freedom）：组间自由度是组数减 1（3 − 1 = 2），残差自由度是总观测数减组数（30 − 3 = 27）。总自由度 29 = 2 + 27。
- Sum Sq（平方和，sum of squares）：组间平方和 3.766 衡量各组均值围绕总均值的离散程度；残差平方和 10.492 衡量组内观测围绕各自组均值的离散程度。
- Mean Sq（均方，mean square）：Sum Sq 除以对应的 Df。组间均方 3.766 / 2 = 1.883，残差均方 10.492 / 27 = 0.3886。残差均方就是组内方差，也是后面算标准误的基础。
- F value：组间均方除以残差均方，1.8832 / 0.3886 = 4.846。直观理解是「组间差异是组内噪声的多少倍」，接近 1 说明组间差异和随机波动差不多。
- Pr(>F)：在「所有组均值相等」的原假设下，出现这么大 F 值的概率。这里 p = 0.0159，拒绝原假设。
结论是「三组中至少有一组和其他组不同」，但没说具体是哪一组。这正是事后检验要回答的问题。
效应量：差异有多大
Section titled “效应量：差异有多大”p 值回答「有没有差异」，不回答「差异有多大」。30 个观测下 p = 0.016 可能对应很大的组间差异，3000 个观测下同一个 p 值对应的差异可能小到没有实际意义。效应量补的就是这一半信息。单因素设计最常用的是 η²（eta squared）：
ss <- summary(fit)[[1]][["Sum Sq"]]eta2 <- ss[1] / sum(ss)eta2[1] 0.2641483η² = 组间平方和 / 总平方和 = 3.766 / (3.766 + 10.492) = 0.264，读作「因变量的变异中有 26.4% 由分组解释」。经验分档是 0.01 小、0.06 中、0.14 大，0.264 属于大效应。
这个数不是新算出来的：前面 lm() 输出的 Multiple R-squared: 0.2641 就是它。单因素方差分析里 R² 和 η² 是同一个量，只是前者从回归角度命名、后者从方差分解角度命名。报告时两者写哪个都可以，但同一篇文章里保持统一。
η² 有一个和 R² 相同的缺点：它随自变量个数增加而膨胀，所以双因素设计里更严谨的做法是报偏 η²（partial eta squared）。car::Anova() 的输出里有各项的 Sum Sq，偏 η² 按 SS_项 / (SS_项 + SS_残差) 算。以 ToothGrowth 的双因素模型为例（结果见下文「双因素与交互作用」一节），supp 的偏 η² 是 205.35 / (205.35 + 712.11) = 0.224，dose 是 2426.43 / (2426.43 + 712.11) = 0.773。
修正 η² 过度乐观的还有一个 ω²（omega squared），样本量小时比 η² 更保守，手算公式是 (SS_组间 − (k−1) × MS_残差) / (SS_总 + MS_残差)。代进 PlantGrowth：(3.766 − 2 × 0.3886) / (14.258 + 0.3886) = 0.204，比 η² 的 0.264 小一截。样本量越小，ω² 与 η² 的差距越大，所以小样本研究报 ω² 更稳妥。
事后多重比较：TukeyHSD
Section titled “事后多重比较：TukeyHSD”Tukey HSD（honestly significant difference，诚实显著差异）在所有两两比较中控制整体错误率，是最常用的方法。
TukeyHSD(fit) Tukey multiple comparisons of means 95% family-wise confidence level
Fit: aov(formula = weight ~ group, data = PlantGrowth)
$group diff lwr upr p adjtrt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960trt2-trt1 0.865 0.1737839 1.5562161 0.0120064diff 是均值差，lwr/upr 是 95% 置信区间，p adj 是已经校正过的 p 值。报告结果时必须看 p adj，不能自己用原始 p 值——Tukey 的校正已经包含在里面了。
三组比较的结论很清晰：trt2 比 trt1 高 0.865（置信区间 [0.174, 1.556]，不含 0，p = 0.012），而 trt1 和 trt2 相对对照组 ctrl 的差异都不显著（p 分别是 0.39 和 0.20）。注意「整体 F 检验显著、但没有任何两两比较显著」是完全可能出现的，尤其组数多、每组样本少的时候，这时候别急着下「有差异」的结论。
aov() 就是 lm()
Section titled “aov() 就是 lm()”aov() 不是独立的建模系统，它拟合的就是普通线性模型，只是换了输出格式。用 lm() 拟合同一个模型会得到系数表：
summary.lm(lm(weight ~ group, data = PlantGrowth))Coefficients: Estimate Std. Error t value Pr(>|t|)(Intercept) 5.0320 0.1971 25.527 <2e-16 ***grouptrt1 -0.3710 0.2788 -1.331 0.1944grouptrt2 0.4940 0.2788 1.772 0.0877 .---Residual standard error: 0.6234 on 27 degrees of freedomMultiple R-squared: 0.2641, Adjusted R-squared: 0.2096F-statistic: 4.846 on 2 and 27 DF, p-value: 0.01591截距 5.032 是 ctrl 组的均值，后两个系数是相对 ctrl 的差值，和 Tukey 输出的 diff 列完全对得上。F-statistic: 4.846 也和方差分析表的 F 值一致。两种写法等价，区别在于：
- 系数表给的是每组对基准组的比较（默认基准是第一个水平），是未校正的 t 检验；
- 方差分析表给的是整体的 F 检验；
TukeyHSD()给的是所有两两组合且已校正的结果。
想对 lm() 对象出一张和 aov() 一样的表，用 anova(lm(weight ~ group, data = PlantGrowth))。
双因素与交互作用
Section titled “双因素与交互作用”ToothGrowth 是 2×3 设计：两种维生素 C 补充方式（supp：OJ 橙汁、VC 抗坏血酸）× 三个剂量（dose：0.5、1、2 mg）。每组 10 只豚鼠，因变量是牙齿长度。
第一个坑：dose 在数据里是数值型（num），直接放进公式会被当成连续变量，只有一个自由度：
class(ToothGrowth$dose)[1] "numeric"要分析三个剂量水平之间的差异，必须显式转成因子：
summary(aov(len ~ supp * factor(dose), data = ToothGrowth)) Df Sum Sq Mean Sq F value Pr(>F)supp 1 205.4 205.4 15.572 0.000231 ***factor(dose) 2 2426.4 1213.2 92.000 < 2e-16 ***supp:factor(dose) 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 ‘ ’ 1supp * factor(dose) 展开成三项：两个主效应加一个交互项。交互项的 p = 0.0219，说明补充方式的效果随剂量变化，不能只报一个「OJ 比 VC 平均高多少」的结论。
看各组均值就明白了：
aggregate(len ~ supp + dose, data = ToothGrowth, FUN = mean) supp dose len1 OJ 0.5 13.232 VC 0.5 7.983 OJ 1.0 22.704 VC 1.0 16.775 OJ 2.0 26.066 VC 2.0 26.14剂量 0.5 时 OJ 比 VC 高 5.25，剂量 1.0 时高 5.93，到剂量 2.0 时两者只差 0.08——剂量足够高之后，补充方式的差别消失了。这就是交互作用的实际形态：一个因素的效应依赖于另一个因素的水平。报告这种结果时应该分组呈现，而不是只给主效应，必要时用 interaction.plot(ToothGrowth$dose, ToothGrowth$supp, ToothGrowth$len) 画一张交互作用图。
顺带一提，如果不控制剂量会怎样：
summary(aov(len ~ supp, data = ToothGrowth)) Df Sum Sq Mean Sq F value Pr(>F)supp 1 205 205.35 3.668 0.0604 .Residuals 58 3247 55.98只看 supp 时 p = 0.060，几乎要得出「补充方式没影响」的结论。把剂量放进模型后 p 值降到 0.00023——因为漏掉剂量这个重要因素时，它的变异全部混进了残差，残差均方从 13.2 膨胀到 55.98，F 检验自然无力。设计里记录了的因素，不要因为「主效应不显著」就把它从模型里删掉。
方差分析有两条假设：残差独立且正态、各组方差齐性。正态性看残差的 QQ 图或检验：
shapiro.test(residuals(fit)) Shapiro-Wilk normality test
data: residuals(fit)W = 0.96607, p-value = 0.4379p = 0.44，没有证据反对正态性。注意是检验残差，不是检验原始因变量——分组之后各组内部的分布才和假设相关。
方差齐性用 Bartlett 检验或 Levene 检验：
library(car) # leveneTest() 来自 car 包
bartlett.test(weight ~ group, data = PlantGrowth)leveneTest(weight ~ group, data = PlantGrowth) Bartlett test of homogeneity of variances
data: weight by groupBartlett's K-squared = 2.8786, df = 2, p-value = 0.2371
Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F)group 2 1.1192 0.3412 27两种检验都不显著，但这同样只是「证据不足」，不是「方差一定相等」——总共才 30 个观测，这类检验的功效本来就不高。这里还有个取舍：Bartlett 检验假设正态，对偏离正态敏感；Levene 检验（默认用中位数中心化）更稳健，数据有偏态时优先用它。如果方差确实不齐，别硬做 ANOVA，改用 oneway.test()：
oneway.test(weight ~ group, data = PlantGrowth) One-way analysis of means (not assuming equal variances)
data: weight and groupF = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739它不假设等方差，自由度也不再是整数（17.128），相当于方差分析版的 Welch t 检验。代价是它只给整体检验，TukeyHSD() 不能用在它上面，事后比较要用 Games-Howell 方法（PMCMRplus 包或 rstatix 包的 games_howell_test()）。
oneway.test() 有一个 var.equal 参数，默认 FALSE。设成 TRUE 它就退化成普通方差分析，输出里能直接对照 F 值：
oneway.test(weight ~ group, data = PlantGrowth, var.equal = TRUE) One-way analysis of means
data: weight and groupF = 4.8461, num df = 2, denom df = 27, p-value = 0.01591F = 4.8461、分母自由度 27，和开头 aov() 的方差分析表完全一致，var.equal = TRUE 时它算的就是同一个 F 检验。反过来也能看出 Welch 校正改了什么：自由度从整数 27 变成 17.128，F 从 4.846 变成 5.181。自由度变小是校正的代价，换来的是方差不齐时更可信的 p 值。
上面这个例子里方差本来是齐的，所以两种算法差别不大。换一组方差确实不齐的数据，差距就明显了。iris 的三个物种在花瓣长度上的离散程度差得很远：
leveneTest(Petal.Length ~ Species, data = iris)Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F)group 2 19.48 3.129e-08 *** 147---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1p = 3.1e-08，方差齐性被明确拒绝。这时候两种算法的结果分道扬镳：
oneway.test(Petal.Length ~ Species, data = iris)summary(aov(Petal.Length ~ Species, data = iris)) One-way analysis of means (not assuming equal variances)
data: Petal.Length and SpeciesF = 1828.1, num df = 2.000, denom df = 78.073, p-value < 2.2e-16
Df Sum Sq Mean Sq F value Pr(>F)Species 2 437.1 218.55 1180 <2e-16 ***Residuals 147 27.2 0.19---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1F 值从 1180 变成 1828，差了 55%。这个例子里两个 p 值都远小于 0.05，结论不受影响，但换成效应没那么强的研究，55% 的差距足以让一个边缘结果跨过或退离显著性门槛。
要不要一遇到方差不齐就换方法，取决于方差差异有多大、各组样本量是否相等。三组样本量相等（各 50）时，普通方差分析对异方差还算稳健，上面也就是 F 值变了、结论没变；各组样本量悬殊时必须换，那时普通方差分析的 I 类错误率会明显偏离设定的 α。判断样本量是否悬殊比判断方差是否齐性更重要，也更常被忽略。
区组设计：把已知的干扰因素拿出来
Section titled “区组设计：把已知的干扰因素拿出来”npk 是 6 个区组的析因试验，研究氮、磷、钾三种肥料对产量的影响。先把区组（block）放进模型：
summary(aov(yield ~ block + N * P * K, data = npk)) Df Sum Sq Mean Sq F value Pr(>F)block 5 343.3 68.66 4.447 0.01594 *N 1 189.3 189.28 12.259 0.00437 **P 1 8.4 8.40 0.544 0.47490K 1 95.2 95.20 6.166 0.02880 *N:P 1 21.3 21.28 1.378 0.26317N:K 1 33.1 33.14 2.146 0.16865P:K 1 0.5 0.48 0.031 0.86275Residuals 12 185.3 15.44区组效应显著（F = 4.447，p = 0.0159）。把它放进模型的价值在残差均方上。同样的数据不放进区组再拟合一次：
summary(aov(yield ~ N * P * K, data = npk)) Df Sum Sq Mean Sq F value Pr(>F)N 1 189.3 189.28 6.161 0.0245 *P 1 8.4 8.40 0.273 0.6082K 1 95.2 95.20 3.099 0.0975 .N:P 1 21.3 21.28 0.693 0.4175N:K 1 33.1 33.14 1.078 0.3145P:K 1 0.5 0.48 0.016 0.9019N:P:K 1 37.0 37.00 1.204 0.2887Residuals 16 491.6 30.72---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1两处对比值得看，方向正好相反。
残差均方几乎砍半：30.72 降到 15.44。残差变小意味着检验更有力——氮肥的 p 值从 0.0245 降到 0.00437，钾肥从 0.0975 降到 0.0288，两个本来边缘的结果变得明确了。区组设计的意义就在这里：把已知的系统性干扰从残差里剥离出来。
但效应的平方和一点没变：N 是 189.3、P 是 8.4、K 是 95.2，两张表里完全相同。区组只改变残差，不改变各处理效应的点估计。这也是为什么均衡设计里即使区组没测或没记录，处理效应的估计仍然是无偏的，代价只是检验变钝（p 值偏大）。反过来说，漏掉区组不会让结论朝着错误的方向走，只会让你更容易错过真实存在的差异。
再看一张表的差异：含区组的模型里 N:P:K 三阶交互项没有出现，不含区组的模型里它出现了（Df = 1，F = 1.204）。npk 是部分析因（fractional factorial）设计，数据集文档写明「confounding the NPK interaction」，三阶交互与区组效应混杂（confounded）——两者用的是同一个自由度，模型里只能留一个。留了区组就估不出三阶交互，丢了区组就能估出来，但那时区组的系统性差异会混进残差。这个实验的设计意图是估区组和处理主效应，三阶交互从一开始就不在估计范围内。这不是报错，但你得知道模型里少了什么——写方法部分时要说明这一点。
方差分析的底层是线性模型，所以 aov() 能做的事 lm() 基本都能做，而且 lm() 还能扩展到连续自变量和交互项，见 /r/modeling/linear-regression；如果数据是重复测量（同一个体测多次），那些观测不独立，普通的 aov() 会低估标准误，要用混合效应模型，见 /r/modeling/mixed-models。做方差分析前用箱线图看一眼各组分布，能提前发现异常值和方差不齐，画法见 /python/visualization/seaborn。