How to run Tukey pairwise comparisons after ANOVA in Statsmodels (with Code)

Editorial · Published 1 November 2026

The Short Answer (TL;DR)

Use pairwise_tukeyhsd with the numeric outcomes, their group labels and use_var="equal" to compare all pairs of means with Tukey-adjusted inference. Report each directed mean difference, adjusted p-value and simultaneous confidence interval while checking independent observations and the common-variance model.

import numpy as np; from statsmodels.stats.multicomp import pairwise_tukeyhsd
y = np.array([8,9,10,11,9,10,8,9,12,13,14,12,15,13,14,13,15,17,16,18,17,16,19,18]); group = np.repeat(["A","B","C"],8)
print(pairwise_tukeyhsd(y,group,alpha=.05,use_var="equal").summary())

Detailed Implementation

This example uses synthetic data and does not need a network connection. Install the recorded package versions in a separate Python environment before running it.

python -m pip install numpy==2.5.4 scipy==1.18.1 pandas==3.0.6 statsmodels==0.15.0

from pathlib import Path
from itertools import combinations
import numpy as np
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd

# Synthetic independent groups, measured in the same response units.
samples={
    "A":np.array([8,9,10,11,9,10,8,9],float),
    "B":np.array([12,13,14,12,15,13,14,13],float),
    "C":np.array([15,17,16,18,17,16,19,18],float),
}
if len(samples) < 2 or any(x.ndim != 1 or x.size < 2 or not np.isfinite(x).all() for x in samples.values()):
    raise ValueError("Use independent finite one-dimensional groups with at least two observations.")
y=np.concatenate(list(samples.values()))
group=np.concatenate([np.repeat(label,x.size) for label,x in samples.items()])
omnibus=stats.f_oneway(*samples.values(),equal_var=True)
result=pairwise_tukeyhsd(y,group,alpha=.05,use_var="equal")
labels=result.groupsunique
pairs=list(combinations(range(len(labels)),2))
lines=[f"Overall ANOVA F={omnibus.statistic:.9g}; p={omnibus.pvalue:.9g}",
       "Tukey family: all three pairwise means; FWER target=.05",
       "comparison group2-minus-group1 adjusted_p simultaneous_95%_CI reject"]
for index,(i,j) in enumerate(pairs):
    low,high=result.confint[index]
    lines.append(f"{labels[i]}->{labels[j]} {result.meandiffs[index]:.6f} "
                 f"{result.pvalues[index]:.9g} [{low:.6f},{high:.6f}] {bool(result.reject[index])}")
report="\n".join(lines)+"\n"
Path("tukey-pairwise.txt").write_text(report,encoding="utf-8")
print(report)

The differences B-A, C-A and C-B are 4.000000, 7.750000 and 3.750000 response units. All three Tukey-adjusted comparisons reject at the prespecified .05 family-wise level, with strictly positive p-values and simultaneous intervals above zero.

What does each comparison mean?

The table uses group2 minus group1, so A->B is B minus A, not the reverse. The three observed differences are 4.00, 7.75 and 3.75 response units. Their simultaneous 95% intervals are approximately [2.570971,5.429029], [6.320971,9.179029] and [2.320971,5.179029]. All three exclude zero in this synthetic example. The full report uses scientific notation for small adjusted p-values; a rounded summary displaying 0.0 does not mean the probability is exactly zero.

How do ANOVA and Tukey answer different questions?

The overall ANOVA tests equality of all group means, whereas Tukey identifies differences within the planned family of all three pairwise comparisons. The code reports both and does not change the family after seeing which differences look large. Tukey-adjusted p-values and simultaneous intervals account for this family; repeatedly running unadjusted two-sample tests would answer a different multiplicity problem. A significant omnibus test is not a universal prerequisite for a prespecified Tukey procedure.

Which assumptions support the stated family-wise inference?

The equal-variance procedure assumes independent normal observations within groups with a common within-group population variance. Balanced synthetic groups illustrate the computation rather than prove those assumptions for real data. Dependence, extreme outliers and unequal variances require an analysis suited to the design. In Statsmodels 0.15, use_var="unequal" selects Games-Howell, a different procedure; it must not be described as the same Tukey calculation.

What do the independent checks reconstruct?

Exact rational arithmetic gives pooled within-group sum of squares 27 and residual degrees of freedom 21, yielding variance 9/7. Each observed mean difference and its standard error is then combined with the studentized range reference distribution to reconstruct the adjusted p-value and simultaneous interval. Reordering observations preserves the results, while changing response units by a positive factor scales differences and interval endpoints without changing the p-values. Rejection is evidence about the specified mean contrasts under the model, not proof of a causal effect or its practical importance.

Mathematical Context

For k independent groups, estimate the common error variance as MSE=sum_i sum_j(y_ij-mean_i)^2/(N-k). For the balanced example k=3, N=24 and MSE=27/21=9/7. A pairwise difference mean_j-mean_i uses s_ij=sqrt((MSE/2)(1/n_i+1/n_j)); its studentized range statistic is q=|difference|/s_ij. Simultaneous limits are difference plus or minus q_(1-alpha;k,N-k) times s_ij. The response differences and limits retain the original measurement units, while q is dimensionless. Under the stated normal common-variance model, the all-pairs family has the intended simultaneous coverage; this is a family-wise criterion rather than a separate unadjusted 95% assertion for each selected comparison.

\[ q_{ij}=\frac{|\bar y_j-\bar y_i|}{\sqrt{\frac{\mathrm{MSE}}{2}(1/n_i+1/n_j)}} \]

Reproducibility

Executed with Python 3.12.14; numpy 2.5.4, scipy 1.18.1, pandas 3.0.6, statsmodels 0.15.0. Automated checks confirm the stated output properties.