均值多重比较校正 - Python 中的 Tukey HSD
Correction for Multiple Comparison of Means - Tukey HSD in Python
我有一个包含 4 个条件(A、B、C、D)的日期集。我观察到 运行单向方差分析是我的因变量(反应时间,RT)在 4 种情况下线性增加。
我想 运行 进行 post-hoc 测试,看看 RT 从 A 到 B、从 B 到 C 和 C 到 D 的增加对于 Tukey HSD 是否显着post-临时测试。
为了运行在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)做出假设。
如何通知 post-hoc 测试我感兴趣的 number/type 比较?
很遗憾,你不能。 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
作为评论,如果您看到趋势,可能有其他模型可以对其建模,而不是依赖 posthoc 测试。此外,对您需要的测试进行子集化并执行更正可以被认为是某种类型的樱桃采摘。如果比较的数量(如您的示例 6),我建议您使用 Tukey。这是您可以 post 进行交叉验证的另一个讨论。
我有一个包含 4 个条件(A、B、C、D)的日期集。我观察到 运行单向方差分析是我的因变量(反应时间,RT)在 4 种情况下线性增加。
我想 运行 进行 post-hoc 测试,看看 RT 从 A 到 B、从 B 到 C 和 C 到 D 的增加对于 Tukey HSD 是否显着post-临时测试。
为了运行在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)做出假设。
如何通知 post-hoc 测试我感兴趣的 number/type 比较?
很遗憾,你不能。 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
作为评论,如果您看到趋势,可能有其他模型可以对其建模,而不是依赖 posthoc 测试。此外,对您需要的测试进行子集化并执行更正可以被认为是某种类型的樱桃采摘。如果比较的数量(如您的示例 6),我建议您使用 Tukey。这是您可以 post 进行交叉验证的另一个讨论。