使用statsmodels做重复测量ANOVA的事后HSD检验是否正确？

问题描述 #

我尝试为重复测量数据执行事后HSD检验，目前使用statsmodels.stats.multicomp.MultiComparison，但不确定该方法能否处理重复测量问题。以下是我的代码：

from statsmodels.stats.anova import AnovaRM
import statsmodels.stats.multicomp as mc

aovrm2way = AnovaRM(df, amp, 'subject', within=['cond', 'type'])
res2way = aovrm2way.fit()

print(res2way)


comp = mc.MultiComparison(df[amp], df['cond'])
post_hoc_res = comp.tukeyhsd()
post_hoc_res.summary()
print(post_hoc_res.summary())

请问这种做法对于重复测量数据是否正确？若不正确，是否有其他库可用于重复测量的事后检验？

回答 #

你的做法是否正确？ #

不正确。MultiComparison搭配TukeyHSD的默认实现是针对独立样本设计的，它没有考虑重复测量数据中同一受试者多次观测的相关性，会错误地估计标准误，导致统计检验的结果不可靠（比如I类错误率偏高）。

替代方案 #

• 使用pingouin库（最简便）
  pingouin专门封装了重复测量设计的事后检验接口，直接支持TukeyHSD，无需手动处理相关性：

  import pingouin as pg
  
  # 执行重复测量的TukeyHSD检验
  post_hoc = pg.pairwise_tukey(data=df, dv=amp, within="cond", subject="subject")
  print(post_hoc)
  

  如果需要检验两因素交互效应的事后比较（比如按type分组后比较cond的差异），可以添加groupby参数：

  post_hoc_interaction = pg.pairwise_tukey(data=df, dv=amp, within="cond", subject="subject", groupby="type")
  print(post_hoc_interaction)
  

• 使用emmeans库（灵活适配复杂模型）
  emmeans可以基于拟合的重复测量模型计算边际均值，并执行校正后的两两比较，支持任意复杂的实验设计：

  from statsmodels.formula.api import mixedlm
  from emmeans import emmeans
  
  # 先拟合混合效应模型（替代AnovaRM的另一种方式）
  model = mixedlm(f"{amp} ~ cond * type", df, groups=df["subject"])
  result = model.fit()
  
  # 针对cond因子计算边际均值并执行TukeyHSD
  emm = emmeans(result, "cond")
  tukey_results = emm.tukey()
  print(tukey_results)
  
  # 针对交互效应的事后检验：按type分组比较cond
  emm_interaction = emmeans(result, ["cond", "type"])
  tukey_interaction = emm_interaction.tukey()
  print(tukey_interaction)
  

• statsmodels结合边际均值（原生实现）
  如果你不想引入额外库，可以基于statsmodels的线性混合模型，手动提取边际均值后执行多重比较校正，但步骤相对繁琐：

  import statsmodels.api as sm
  from statsmodels.formula.api import mixedlm
  from statsmodels.stats.multitest import multipletests
  
  # 拟合模型
  model = mixedlm(f"{amp} ~ cond * type", df, groups=df["subject"])
  result = model.fit()
  
  # 生成各cond水平的预测值（近似边际均值）
  cond_levels = df["cond"].unique()
  mean_estimates = []
  for cond in cond_levels:
      temp_df = df[df["cond"] == cond].copy()
      temp_df["cond"] = cond
      pred = result.predict(temp_df)
      mean_estimates.append(pred.mean())
  
  # 手动计算两两比较的p值并校正（Bonferroni或Tukey）
  # 此方法需自行计算标准误，推荐优先使用前两种方案
  

内容的提问来源于stack exchange，提问作者Noa Guttman