4. 单因素方差分析及事后多重比较实现

任务对象与数据形态

本案例考察三种培养基配方（A、B、C）对某酶活指标（U/mL）的影响。与均衡设计不同，实际实验中因污染、样本损耗等原因，各组有效重复数往往不等：A 组 6 个、B 组 5 个、C 组 7 个。不均衡样本量会直接影响平方和分解方式与事后比较的置信区间宽度，因此必须使用Type II 平方和（对单因素而言与 Type III 等价）并采用 Tukey-Kramer 校正，而非仅适用于等样本量的 Tukey HSD 原版公式。

示例数据（教学假设值，非实测）：

方差分析的计算路径

总样本量 N = 18，组数 k = 3。组间平方和 SSB = Σ nᵢ(x̄ᵢ − x̄..)²，组内平方和 SSW = Σ(nᵢ − 1)sᵢ²，自由度分别为 k−1 = 2 与 N−k = 15。代入数据得 x̄.. ≈ 12.006，SSB ≈ 19.63，SSW ≈ 2.72，F = (19.63/2)/(2.72/15) ≈ 54.1，远大于 F₀.₉₅(2,15) ≈ 3.68，整体组间差异显著。

整体显著只说明“至少有一组不同”，具体哪两组不同需事后比较。Tukey-Kramer 的检验统计量为 q = |x̄ᵢ − x̄ⱼ| / √(MSE/2 · (1/nᵢ + 1/nⱼ))，其中 MSE = SSW/(N−k) ≈ 0.1813。由于 nᵢ 不等，每对比较的标准误不同，这是与等样本量情形的关键区别。

Python 实现（含不均衡样本量处理）

下面代码自包含数据、方差分析、Tukey-Kramer 事后比较与可视化，可直接运行。关键点在于使用 statsmodels 的 anova_lm 指定 Type II 平方和，并用 pairwise_tukeyhsd 输出校正后的比较矩阵。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 1. 构造不均衡数据
data = {
    "A": [12.1, 11.4, 12.8, 11.9, 12.5, 11.7],
    "B": [10.3, 10.9, 10.1, 10.6, 10.4],
    "C": [13.2, 12.6, 13.8, 13.1, 12.9, 13.5, 13.0],
}
df = pd.DataFrame(
    [(g, v) for g, vals in data.items() for v in vals],
    columns=["group", "activity"]
)

# 2. 单因素方差分析（Type II 平方和）
model = ols("activity ~ C(group)", data=df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

# 3. Tukey-Kramer 事后比较（自动处理不等 n）
tukey = pairwise_tukeyhsd(
    endog=df["activity"], groups=df["group"], alpha=0.05
)
print(tukey)

# 4. 提取显著性矩阵
pivot = tukey.summary().data[1:]  # 每行: group1, group2, meandiff, p-adj, lower, upper, reject
matrix = pd.DataFrame(
    [(r[0], r[1], float(r[3]), bool(r[6])) for r in pivot],
    columns=["g1", "g2", "p_adj", "reject"]
)
print(matrix)

# 5. 均值与同时置信区间可视化
fig, ax = plt.subplots(figsize=(6, 4))
means = df.groupby("group")["activity"].mean()
sems = df.groupby("group")["activity"].sem()
ax.errorbar(means.index, means.values, yerr=1.96 * sems.values,
            fmt="o", capsize=6, color="#2f6fa8")
ax.set_ylabel("酶活 (U/mL)")
ax.set_title("各组均值与 95% 置信区间")
plt.tight_layout()
plt.savefig("tukey_result.png", dpi=150)

结果判读与不均衡情形的注意事项

• 整体检验：F = 54.1，p < 0.001，拒绝各组均值相等的原假设。
• 事后比较：三对比较的校正 p 值均小于 0.05，说明 C 组显著高于 A 组，A 组显著高于 B 组，排序为 C > A > B。
• 不均衡的影响：B 组仅 5 个重复，其与 C 组比较的标准误最大，置信区间最宽；若 B 组样本量进一步减少，可能出现“整体显著但某对比较不显著”的情形，此时应报告效应量与置信区间，而非仅看 p 值。
• 前提核查：Tukey-Kramer 仍要求方差齐性与近似正态。若 Levene 检验提示方差不齐，应改用 Games-Howell 事后比较。