如何在数据帧上执行双向方差分析后的 Sidak 测试（多计算）？

Posted 2023-03-12

【中文标题】如何在数据帧上执行双向方差分析后的 Sidak 测试（多计算）？

【英文标题】：How to perform Sidak's test (multi comp) following two-way anova on a dataframe?

【发布时间】：2022-01-13 03:26:16

【问题描述】：

我有一个包含以下列的数据框：整整 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

然后我使用以下代码执行双向 anova 以比较药物在每个时间点对流动性的影响：

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

然后我想做一个多重比较测试（事后），特别是 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

任何帮助都会有很大帮助

【参考方案1】：

如果您只对组内比较感兴趣，那么您的系数中已经包含它们。

使用示例数据集：

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+A和 Time[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

【讨论】：

非常感谢，这很有意义，并且能够将其应用于我的实际数据框。只是为了更好地理解它，T. 是什么意思？例如在 Drug[T.B]?

是的，这意味着药物 B 为真，就像 T.19 一样，T.20 意味着时间 == 19，时间==20 的真

感谢您的澄清