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
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