均值多重比较的校正 - Python 中的 Tukey HSD答案

【问题标题】：Correction for Multiple Comparison of Means - Tukey HSD in Python均值多重比较的校正 - Python 中的 Tukey HSD
【发布时间】：2020-11-10 13:18:07
【问题描述】：

我有一个包含 4 个条件（A、B、C、D）的日期集。我在运行单向方差分析时观察到，在 4 种条件下，我的因变量（反应时间，RT）呈线性增加。

我想运行一个事后测试，看看 RT 从 A 到 B、从 B 到 C 和 C 到 D 的增加对于 Tukey HSD 事后测试是否显着。

要在 Python 中运行测试，我使用以下代码：

#Multiple Comparison of Means - Tukey HSD
from statsmodels.stats.multicomp import pairwise_tukeyhsd
print(pairwise_tukeyhsd(df["RT"], df['Cond']))

我面临的问题是，这里假设我对所有可能的比较（A vs B，A vs C，A vs D，B vs C，B vs D，C vs D）感兴趣。因此，所应用的校正基于 6 次测试。但是，我只对 3 次比较（A vs B、B vs C、C vs D）做出假设。

我如何告知事后测试我感兴趣的比较数量/类型？

【问题讨论】：

标签： python python-3.x statsmodels posthoc tukey

【解决方案1】：

很遗憾，你不能。 Tukey HSD 不像您对原始 p 值进行多重比较调整的成对 t 检验。您看到的 p 值基于the studentized range (q) distribution。

您可以做到这一点的一种方法是拟合一个线性模型，就像您的方差分析一样，您对系数进行成对 t 检验，并在您需要的那些上进行子集。

为了说明这一点，我使用了一些模拟数据，这就是 TukeyHSD 的样子：

import pandas as pd
import numpy as np
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multitest import multipletests

np.random.seed(123)

df = pd.DataFrame({'RT':np.random.randn(100),'Cond':np.random.choice(['A','B','C','D'],100)})

hs_res=pairwise_tukeyhsd(df["RT"], df['Cond'])
print(hs_res)

Multiple Comparison of Means - Tukey HSD, FWER=0.05
===================================================
group1 group2 meandiff p-adj   lower  upper  reject
---------------------------------------------------
     A      B  -0.6598 0.2428 -1.5767 0.2571  False
     A      C  -0.3832 0.6946 -1.3334  0.567  False
     A      D   -0.634 0.2663 -1.5402 0.2723  False
     B      C   0.2766 0.7861 -0.5358 1.0891  False
     B      D   0.0258    0.9 -0.7347 0.7864  False
     C      D  -0.2508 0.8257 -1.0513 0.5497  False
---------------------------------------------------

现在我们做 ols，你可以看到它非常相似：

res = ols("RT ~ Cond", df).fit()
pw = res.t_test_pairwise("Cond",method="sh")
pw.result_frame

    coef    std err t   P>|t|   Conf. Int. Low  Conf. Int. Upp. pvalue-sh   reject-sh
B-A -0.659798   0.350649    -1.881645   0.062914    -1.355831   0.036236    0.352497    False
C-A -0.383176   0.363404    -1.054407   0.294343    -1.104528   0.338176    0.829463    False
D-A -0.633950   0.346604    -1.829032   0.070499    -1.321954   0.054054    0.352497    False
C-B 0.276622    0.310713    0.890281    0.375541    -0.340138   0.893382    0.829463    False
D-B 0.025847    0.290885    0.088858    0.929380    -0.551555   0.603250    0.929380    False
D-C -0.250774   0.306140    -0.819147   0.414731    -0.858458   0.356910    0.829463    False

然后我们选择correction的子集和方法，下面我像上面一样使用simes-hochberg：

subdf = pw.result_frame.loc[['B-A','C-B','D-C']]
subdf['adj_p'] = multipletests(subdf['P>|t|'].values,method='sh')[1]
subdf

    coef    std err t   P>|t|   Conf. Int. Low  Conf. Int. Upp. pvalue-sh   reject-sh   adj_p
B-A -0.659798   0.350649    -1.881645   0.062914    -1.355831   0.036236    0.352497    False   0.188742
C-B 0.276622    0.310713    0.890281    0.375541    -0.340138   0.893382    0.829463    False   0.414731
D-C -0.250774   0.306140    -0.819147   0.414731    -0.858458   0.356910    0.829463    False   0.414731

作为评论，如果您看到一种趋势，可能会有其他模型对其进行建模，而不是依赖事后测试。此外，对您需要的测试进行子集化并进行校正可以被认为是某种类型的樱桃采摘。如果比较的次数（如您的示例 6），我建议您使用 Tukey。这是您可以在交叉验证上发布的另一个讨论。

【讨论】：