rstatisticsp-valueemmeanstukey

How does emmeans adjust the p-values when using "Tukey" as adjustment method?


I want to calculate a critical difference for my pairwise comparisons. For this I need to understand how the p-values are adjusted, since I want to have the tukey method active.

I was able to reconstruct the T-Statistics calculation and unadjusted p-values. But using both AI and the internet, I couldn't come up with a formula for the adjusted p-values that gave me the same result as emmeans does.

This is the example-case I want to understand:

              Estimate  2.5_ci  97.5_ci     SE       DF  ci_upper  ci_lower
algorithm                                                                  
PriorBand+BO     2.556   2.504    2.608  0.027  12000.0     2.504     2.608
BOHB             2.625   2.573    2.677  0.027  12000.0     2.573     2.677
RS+Prior         3.223   3.171    3.275  0.027  12000.0     3.171     3.275
PiBO             3.223   3.171    3.275  0.027  12000.0     3.171     3.275
BO               3.373   3.321    3.425  0.027  12000.0     3.321     3.425

       algorithm_1     algorithm_2  Estimate  2.5_ci  97.5_ci     SE       DF  T-stat  P-val  Sig
7   (PriorBand+BO)            BOHB    -0.069  -0.172    0.033  0.038  12000.0  -1.843  0.349     

So 0.349 is my "target".

The t_stat is -1.8157894736842106, and the unadjusted p-value is 0.06942762000543438

AI proposed this formula, which calculates a q value from the studentized range and from this a correction factor:

q=studentized_range.ppf(1 - 0.05, k=k,df=df)
factor=q/np.sqrt(k)
adjusted_p=p*factor

But the result does not fit my target:

q=3.8582424830229995
factor=1.725458493143401
adjusted_p=0.11979447659710946

Solution

  • As mentioned by @Matt Haberland, a similar question has been asked here. So to calculate the p-value of my example, you use

    1 - ptukey(t_stat * sqrt(2), k, df) (in R) or

    1 - scipy.stats.studentized_range.cdf(t_stat*np.sqrt(2), k=k, df=df) (in Python)

    where t_stat=abs(mean_difference/SE), k=5 is the number of groups and df=12000. And to calculate the critical difference HSD (of this specific pair):

    HSD=scipy.stats.studentized_range.ppf(1-0.05, k=k,df=df)/np.sqrt(2)*SE