Python 方差分析：单因素、双因素与事后检验

三组以上的均值比较，不能靠反复做 t 检验解决。三组两两比是 3 次检验，每次都用 α = 0.05，至少犯一次第一类错误的概率是 1 − 0.95³ ≈ 14.3%；五组要比较 10 次，这个概率涨到 1 − 0.95¹⁰ ≈ 40.1%。方差分析（ANOVA，analysis of variance）的做法是先做一次整体的 F 检验，把族错误率（family-wise error rate）控制在 5% 以内，确认存在差异后再做事后多重比较。

多重比较的问题有多大

上面的 14.3% 是三次检验相互独立时的上界。三组的两两比较共用同一批数据，检验之间正相关，实际错误率会低于这个上界，但离 5% 仍然很远。用模拟量一次：

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
reject = 0
for _ in range(10000):
    groups = [rng.normal(0, 1, 20) for _ in range(3)]
    p_values = [
        stats.ttest_ind(groups[i], groups[j]).pvalue
        for i, j in [(0, 1), (0, 2), (1, 2)]
    ]
    if min(p_values) < 0.05:
        reject += 1

print(reject / 10000)

0.1217

三组数据都抽自同一个正态分布，真实差异为零，但 10000 次实验里有 1217 次至少有一次 t 检验报出显著。实验整体的错误率是 12.2%，不是 5%。

单因素方差分析

下面用一组模拟数据演示单因素设计：三种处理（ctrl、trt1、trt2），每组 10 个观测，因变量是干重。数据由固定种子生成，结果可复现。

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

rng = np.random.default_rng(7)
spec = {"ctrl": (5.03, 0.58), "trt1": (4.66, 0.79), "trt2": (5.53, 0.44)}
rows = [
    {"group": g, "weight": round(float(v), 2)}
    for g, (m, s) in spec.items()
    for v in rng.normal(m, s, 10)
]
pg = pd.DataFrame(rows)

print(pg.groupby("group")["weight"].agg(["count", "mean", "std"]).round(3))

       count   mean    std
group
ctrl      10  4.911  0.396
trt1      10  4.320  0.705
trt2      10  5.262  0.419

statsmodels 用公式接口拟合模型，再用 anova_lm() 把拟合结果转成方差分析表。typ=2 指定平方和的类型，单因素设计下三种类型结果相同，双因素以上才有区别。

model = smf.ols("weight ~ C(group)", data=pg).fit()
table = sm.stats.anova_lm(model, typ=2)
table["mean_sq"] = table["sum_sq"] / table["df"]
print(table)

           sum_sq    df         F    PR(>F)   mean_sq
C(group)  4.53282   2.0  8.197275  0.001652  2.266410
Residual  7.46505  27.0       NaN       NaN  0.276483

C(group) 表示把 group 当作分类变量（categorical）处理。R 里写 aov(weight ~ group) 时，字符列自动按因子处理；Python 的 patsy 公式接口不会，字符串列必须用 C() 包起来，否则 group 会被当成连续变量，只占一个自由度。

anova_lm() 的输出顺序是 sum_sq、df、F、PR(>F)，不含均方（mean square），而均方是理解 F 值的关键，所以上面手动补了 mean_sq 一列。

读懂方差分析表

五个量的含义：

• sum_sq（平方和，sum of squares）：组间平方和 4.533 衡量各组均值围绕总均值的离散程度；残差平方和 7.465 衡量组内观测围绕各自组均值的离散程度。
• df（自由度，degrees of freedom）：组间自由度是组数减 1（3 − 1 = 2），残差自由度是总观测数减组数（30 − 3 = 27）。两者之和 29 是总自由度。
• mean_sq（均方）：sum_sq 除以对应的 df。组间均方 4.533 / 2 = 2.266，残差均方 7.465 / 27 = 0.276。残差均方就是组内方差，后面算标准误要用到它。
• F：组间均方除以残差均方，2.266 / 0.276 = 8.197。直观理解是「组间差异是组内噪声的多少倍」，接近 1 说明组间差异和随机波动差不多。
• PR(>F)：在「三组均值相等」的原假设下，出现这么大 F 值的概率。这里 p = 0.0017，拒绝原假设。

结论是「三组中至少有一组和其他组不同」，但没说具体是哪一组。这正是事后检验要回答的问题。

同一件事也可以用 SciPy 做，f_oneway() 接受若干个数组：

groups = [g["weight"].to_numpy() for _, g in pg.groupby("group")]
print(stats.f_oneway(*groups))

F_onewayResult(statistic=np.float64(8.197275302911567), pvalue=np.float64(0.0016520975278473656))

F 值和 p 值与上面的方差分析表完全一致，说明 f_oneway() 算的就是单因素方差分析。代价是它只返回这两个数，没有平方和与自由度，做报告时还得回到 statsmodels。

事后多重比较：Tukey HSD

Tukey HSD（honestly significant difference，诚实显著差异）在所有两两比较中控制整体错误率，是最常用的方法。statsmodels 的接口是 pairwise_tukeyhsd()，直接传因变量和分组变量：

from statsmodels.stats.multicomp import pairwise_tukeyhsd

print(pairwise_tukeyhsd(pg["weight"], pg["group"]))

Multiple Comparison of Means - Tukey HSD, FWER=0.05
==================================================
group1 group2 meandiff p-adj  lower  upper  reject
--------------------------------------------------
  ctrl   trt1   -0.591 0.0465 -1.174 -0.008   True
  ctrl   trt2    0.351 0.3102 -0.232  0.934  False
  trt1   trt2    0.942 0.0012  0.359  1.525   True
--------------------------------------------------

meandiff 是均值差，lower/upper 是 95% 置信区间，p-adj 是已经校正过的 p 值。报告结果时必须看 p-adj，不能自己用原始 p 值——Tukey 的校正已经包含在里面了。

三组比较的结论：trt1 比 ctrl 低 0.591（p = 0.0465，置信区间 [−1.174, −0.008] 不含 0），trt2 比 trt1 高 0.942（p = 0.0012），而 trt2 与 ctrl 的差异不显著（p = 0.31）。整体 F 检验显著、个别两两比较不显著，这种情况很常见，尤其组数多、每组样本少的时候。别因为整体显著就默认所有组间都有差异。

方差分析就是线性模型

smf.ols() 拟合的是普通线性模型，方差分析只是它的另一种输出格式。看系数表就清楚了：

print(model.params.round(4))
print(model.bse.round(4))

Intercept           4.911
C(group)[T.trt1]   -0.591
C(group)[T.trt2]    0.351
dtype: float64
Intercept           0.1663
C(group)[T.trt1]    0.2352
C(group)[T.trt2]    0.2352
dtype: float64

截距 4.911 是 ctrl 组的均值，后两个系数是相对 ctrl 的差值，和 Tukey 输出的 meandiff 列对得上。标准误 0.2352 来自残差均方：√(2 × 0.276 / 10) = 0.235。

三种输出回答的是三个不同的问题：

输出	比较方式	是否校正
系数表	每组对基准组	否，是普通 t 检验
方差分析表	整体 F 检验	不涉及多重比较
Tukey HSD	所有两两组合	是

双因素与交互作用

单因素设计只回答「组间有没有差异」。研究里更常见的是两个因素同时变化，这时要关心它们是否相互影响。下面构造一个 2×3 设计：两种处理方式（supp：OJ、VC）× 三个剂量（dose：0.5、1.0、2.0），每组 10 个观测。

rng = np.random.default_rng(0)
supp = np.repeat(["OJ", "VC"], 30)
dose = np.tile(np.repeat([0.5, 1.0, 2.0], 10), 2)
means = {
    ("OJ", 0.5): 13.23, ("VC", 0.5): 7.98,
    ("OJ", 1.0): 22.70, ("VC", 1.0): 16.77,
    ("OJ", 2.0): 26.06, ("VC", 2.0): 26.14,
}
length = np.round(
    [means[(s, d)] for s, d in zip(supp, dose)] + rng.normal(0, 3.0, 60), 1
)
tg = pd.DataFrame({"supp": supp, "dose": dose, "len": length})

print(tg.groupby(["supp", "dose"])["len"].mean().round(2))

supp  dose
OJ    0.5     13.47
      1.0     21.34
      2.0     26.07
VC    0.5      8.35
      1.0     19.43
      2.0     25.59
Name: len, dtype: float64

dose 在数据里是数值型，直接放进公式会被当成连续变量，只占一个自由度。要比较三个剂量水平之间的差异，必须用 C() 显式转成分类变量：

two_way = smf.ols("len ~ C(supp) * C(dose)", data=tg).fit()
table2 = sm.stats.anova_lm(two_way, typ=2)
table2["mean_sq"] = table2["sum_sq"] / table2["df"]
print(table2)

                      sum_sq    df           F        PR(>F)      mean_sq
C(supp)            94.000167   1.0   14.766099  3.224293e-04    94.000167
C(dose)          2280.200333   2.0  179.093641  1.469189e-24  1140.100167
C(supp):C(dose)    56.464333   2.0    4.434875  1.647042e-02    28.232167
Residual          343.761000  54.0         NaN           NaN     6.365944

C(supp) * C(dose) 展开成三项：两个主效应加一个交互项 C(supp):C(dose)。交互项的 p = 0.0165，说明处理方式的效果随剂量变化，不能只报一个「OJ 比 VC 平均高多少」的结论。

上面各组均值就说明了这件事：剂量 0.5 时 OJ 比 VC 高 5.12，剂量 1.0 时高 1.91，到剂量 2.0 时只差 0.48——剂量足够高之后，两种处理的差别基本消失。这就是交互作用的实际形态：一个因素的效应依赖于另一个因素的水平。报告这类结果应该分组呈现，而不是只给主效应。

如果建模时漏掉剂量会怎样：

one_way = smf.ols("len ~ C(supp)", data=tg).fit()
table3 = sm.stats.anova_lm(one_way, typ=2)
table3["mean_sq"] = table3["sum_sq"] / table3["df"]
print(table3)

               sum_sq    df         F    PR(>F)    mean_sq
C(supp)     94.000167   1.0  2.034009  0.159175  94.000167
Residual  2680.425667  58.0       NaN       NaN  46.214236

只看 supp 时 p = 0.159，会得出「处理方式没有影响」的结论。把剂量放进模型后，C(supp) 的 p 值降到 0.00032。原因是漏掉剂量这个重要因素时，它的变异全部混进了残差，残差均方从 6.37 膨胀到 46.21，F 检验自然无力。设计里记录了的因素，不要因为主效应不显著就把它从模型里删掉。

检查假设

方差分析有两条假设：残差独立且正态、各组方差齐性。正态性检验的对象是残差，不是原始因变量——分组之后各组内部的分布才和假设相关。

print(stats.shapiro(model.resid))
print(stats.levene(*groups))

ShapiroResult(statistic=np.float64(0.9821017726860367), pvalue=np.float64(0.8782576412293355))
LeveneResult(statistic=np.float64(2.959706204181805), pvalue=np.float64(0.06883476078535292))

两个 p 值都不显著，没有证据反对假设。这同样只是「证据不足」，不是「假设一定成立」——总共才 30 个观测，这类检验的功效本来就不高。

方差齐性一旦被拒绝，普通方差分析的标准误就不对了。iris 数据集的三个物种在花瓣长度上的方差差异很大：

from statsmodels.stats.oneway import anova_oneway
from sklearn.datasets import load_iris

iris = load_iris(as_frame=True).frame
iris.columns = ["sepal_len", "sepal_wid", "petal_len", "petal_wid", "species"]
petal = [g["petal_len"].to_numpy() for _, g in iris.groupby("species")]

print(stats.levene(*petal))
print(anova_oneway(petal, use_var="unequal"))

LeveneResult(statistic=np.float64(19.480338801923573), pvalue=np.float64(3.1287566394085397e-08))
AnovaResult(statistic=np.float64(1828.0919450856877), pvalue=np.float64(2.6933273587151527e-66), df=(2.0, np.float64(78.07295548394265)), df_num=2.0, df_denom=np.float64(78.07295548394265), nobs_total=np.float64(150.0), n_groups=3, means=array([1.462, 4.26 , 5.552]), nobs=array([50., 50., 50.]), vars_=array([0.03015918, 0.22081633, 0.30458776]), use_var='unequal', welch_correction=True, df2=None, df_num2=None, pvalue2=None)

Levene 检验 p = 3.1e-08，方差齐性被明确拒绝。anova_oneway(use_var="unequal") 使用 Welch 校正，分母自由度不再是整数（78.07），可以理解为方差分析版的 Welch t 检验。代价是它只给整体检验，pairwise_tukeyhsd() 不能用在它上面，事后比较要用 Games-Howell 方法（scikit_posthocs 包的 posthoc_gameshowell()）。

要不要在方差齐性被拒绝时就换方法，取决于方差差异有多大、各组样本量是否相等。样本量相等且接近时，普通方差分析对异方差还算稳健；样本量悬殊时必须换。

区组设计：把已知的干扰因素拿出来

研究里常有一些已知的系统性干扰，比如不同批次、不同地块、不同测量日。把这些因素作为区组（block）放进模型，能把它们的变异从残差里剥离出来，检验因此更有力。

下面是一个随机区组设计的模拟数据：6 个区组 × 3 种处理，区组之间有明显的基线差异。

rng = np.random.default_rng(1)
blocks = np.repeat([f"B{i}" for i in range(1, 7)], 3)
treat = np.tile(["ctrl", "low", "high"], 6)
block_effect = np.repeat([0, 2, 4, 1, 3, 5], 3)
treat_effect = np.tile([0, 3, 4.5], 6)
yield_kg = np.round(10 + block_effect + treat_effect + rng.normal(0, 1.0, 18), 2)
bd = pd.DataFrame({"block": blocks, "treat": treat, "yield_kg": yield_kg})

print(bd.groupby("treat")["yield_kg"].agg(["count", "mean"]).round(2))

       count   mean
treat
ctrl       6  12.28
high       6  17.15
low        6  15.87

把区组放进模型：

with_block = smf.ols("yield_kg ~ C(block) + C(treat)", data=bd).fit()
t4 = sm.stats.anova_lm(with_block, typ=2)
t4["mean_sq"] = t4["sum_sq"] / t4["df"]
print(t4)

             sum_sq    df           F        PR(>F)    mean_sq
C(block)  45.452467   5.0   28.367014  1.335667e-05   9.090493
C(treat)  76.681733   2.0  119.643221  1.038740e-07  38.340867
Residual   3.204600  10.0         NaN           NaN   0.320460

不进区组时：

without_block = smf.ols("yield_kg ~ C(treat)", data=bd).fit()
t5 = sm.stats.anova_lm(without_block, typ=2)
t5["mean_sq"] = t5["sum_sq"] / t5["df"]
print(t5)

             sum_sq    df          F    PR(>F)    mean_sq
C(treat)  76.681733   2.0  11.819722  0.000828  38.340867
Residual  48.657067  15.0        NaN       NaN   3.243804

区组效应显著（F = 28.37，p = 1.3e-05），把它放进模型的价值体现在残差均方上：不含区组时是 3.244，含区组后降到 0.320，缩小到原来的 1/10。残差变小意味着检验更有力——处理的 F 值从 11.82 升到 119.64，p 值从 0.00083 降到 1.0e-07。区组设计的意义就在这里：把已知的系统性干扰从残差里剥离出来。

注意处理的 sum_sq 在两张表里都是 76.68，没有变化。区组只改变残差，不改变处理效应本身的估计，这也是为什么均衡设计下缺失区组不至于让处理效应的点估计产生偏倚，只是让检验变钝。

下一步

方差分析的底层是线性模型，所以它处理不了的两件事——连续自变量、更复杂的交互结构——都能在回归里做，见 statsmodels 线性回归；因变量是二分类时要换成逻辑回归，见 Python 逻辑回归。如果数据是重复测量（同一个体测多次），那些观测不独立，普通的方差分析会低估标准误，需要用混合效应模型，R 语言的 /r/modeling/mixed-models/ 讲了这类模型的写法，Python 侧用 statsmodels 的 MixedLM 对应。

做方差分析前先用箱线图看一眼各组分布，能提前发现异常值和方差不齐，画法见 Seaborn 统计可视化。