如何在数据帧上执行双向方差分析后的 Sidak 测试(多比较)?

How to perform Sidak's test (multi comp) following two-way anova on a dataframe?

我有一个包含以下列的数据框:整整 24 小时的时间、药物和流动性。

数据帧的快照

Time Drug Mobility
18   A    1.2 
19   A    1.3
20   A    1.3
21   A    1.2
18   B    3.2
19   B    3.2
20   B    3.3
21   B    3.3

然后我使用以下代码执行双向方差分析来比较药物在每个时间点对流动性的影响:

mod = ols('Mobility~Time+Drug+Time*Drug', data = fdf).fit()
aov = anova_lm(mod, type=2) 

然后我想做一个多重比较测试(post-hoc),特别是 sidak 但不知道如何做。努力在网上找到任何资源来学习。

我知道我可以使用 tukey 的测试,但它比较了一切,我也只对同一时间点的药物效果感兴趣:

18+A - 18+B
19+A - 19+B
20+A - 20+B

不是:

18+A - 19+B
18+A - 20+B
20+A - 18+A

任何帮助都会有很大帮助

如果您只对组内比较感兴趣,那么您的系数中已经有了它们。

使用示例数据集:

import statsmodels.formula.api as sm
import pandas as pd
import numpy as np
from statsmodels.stats.anova import anova_lm
from statsmodels.stats.multitest import multipletests

fdf = pd.DataFrame({'Time':np.random.choice([18,19,20],50),
'Drug':np.random.choice(['A','B'],50), 
'Mobility':np.random.uniform(0,1,50)})

fdf['Time'] = pd.Categorical(fdf['Time'])

mod = sm.ols('Mobility~Time+Drug+Time:Drug', data = fdf).fit()
aov = anova_lm(mod, type=2)

结果如下所示:

mod.summary()

                            OLS Regression Results                            
==============================================================================
Dep. Variable:               Mobility   R-squared:                       0.083
Model:                            OLS   Adj. R-squared:                 -0.021
Method:                 Least Squares   F-statistic:                    0.7994
Date:                Wed, 08 Dec 2021   Prob (F-statistic):              0.556
Time:                        07:13:14   Log-Likelihood:                -4.4485
No. Observations:                  50   AIC:                             20.90
Df Residuals:                      44   BIC:                             32.37
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
========================================================================================
                           coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------------
Intercept                0.4063      0.094      4.323      0.000       0.217       0.596
Time[T.19]               0.1622      0.149      1.091      0.281      -0.137       0.462
Time[T.20]              -0.0854      0.133     -0.643      0.524      -0.353       0.182
Drug[T.B]                0.0046      0.149      0.031      0.975      -0.295       0.304
Time[T.19]:Drug[T.B]    -0.1479      0.206     -0.717      0.477      -0.564       0.268
Time[T.20]:Drug[T.B]     0.2049      0.199      1.028      0.310      -0.197       0.607

在这种情况下,Time=18 是参考,所以 Drug[T.B] 将是 B 对 A 的影响,在 Time=18,这是 18+B - 18+A 的结果,Time[T.19]:Drug[T.B] 将是 19+B - 19+ATime[T.20]:Drug[T.B] 将是 20+B - 20+A .

对于多重比较调整,您可以简单地提取这些结果并计算更正后的 pvalues:

res = pd.concat([mod.params,mod.pvalues],axis=1)
res.columns=['coefficient','pvalues']
res = res[res.index.str.contains('Drug')]
res['corrected_p'] = multipletests(res['pvalues'],method="sidak")[1]
res['comparison'] = ['18+B - 18+A','19+B - 19+A','20+B - 20+A']

                      coefficient   pvalues  corrected_p   comparison
Drug[T.B]                0.004630  0.975284     0.999985  18+B - 18+A
Time[T.19]:Drug[T.B]    -0.147928  0.477114     0.857038  19+B - 19+A
Time[T.20]:Drug[T.B]     0.204925  0.309616     0.670942  20+B - 20+A

另见 help page for multipletests