TopolearnDesigning an experiment with enough data, and testing many groups honestly, are core statistical skills.
ANOVA tables: sm.stats.anova_lm(res, typ=2) on a fitted OLS model gives an F-test per term. Type 2 tests each term after the others (respecting marginality) and is a good default without interactions; with interactions and unbalanced data, type 3 with sum-to-zero contrasts is common. anova_lm(m1, m2) compares nested models.
Multiple comparisons: after finding that groups differ, pairwise_tukeyhsd(endog, groups) from statsmodels.stats.multicomp tests all pairs while controlling the family-wise error rate. For an arbitrary list of p-values, multipletests(pvals, method="fdr_bh") from statsmodels.stats.multitest applies Benjamini-Hochberg, Holm, Bonferroni and others.
Power analysis answers "how many samples do I need?" before running an experiment. statsmodels.stats.power has classes such as TTestIndPower, NormalIndPower and FTestAnovaPower. Their solve_power method takes all but one of effect size, sample size, alpha and power, and solves for the missing one. Effect sizes are standardised: Cohen's d for means, proportion_effectsize(p1, p2) for proportions.
Underpowered studies are a major source of unreliable findings: they miss real effects, and the significant results they do produce overstate the effect size. Fix the sample size in advance, and don't stop early the moment p dips below 0.05.
Example
import numpy as np
import pandas as pd
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.power import NormalIndPower, TTestIndPower
from statsmodels.stats.proportion import proportion_effectsize
n = TTestIndPower().solve_power(effect_size=0.3, alpha=0.05, power=0.8)
print(round(n)) # 175 per group for d = 0.3
es = proportion_effectsize(0.12, 0.10) # conversion 10% -> 12%
print(round(NormalIndPower().solve_power(es, alpha=0.05, power=0.8))) # 3835 per group
rng = np.random.default_rng(12)
df = pd.DataFrame({"group": np.repeat(["A", "B", "C"], 30),
"y": np.r_[rng.normal(10, 2, 30), rng.normal(10.2, 2, 30), rng.normal(12, 2, 30)]})
tukey = pairwise_tukeyhsd(df["y"], df["group"])
print(tukey.reject.tolist()) # [False, True, False] (A-B, A-C, B-C)Remember
anova_lm(res, typ=2)tests each term; Tukey HSD compares all pairsmultipletestsadjusts lists of p-values (Holm, FDR)- Use
solve_powerto fix sample size before the experiment
Watch out for
- Peeking at results repeatedly and stopping as soon as p < 0.05.