statistics

name	statistics
description	Reference for applied statistics on tabular data: descriptives, confidence intervals, one/two-sample and paired tests (t-test, Mann-Whitney, Wilcoxon, chi-square, Fisher exact, ANOVA, Kruskal-Wallis), assumption checks (Shapiro, Levene, Bartlett), correlation (Pearson, Spearman, Kendall), effect sizes (Cohen's d, Cliff's delta, eta-squared, phi), OLS with residual diagnostics, multiple-comparison correction (Bonferroni, BH), and a single power-analysis example via statsmodels.
metadata	{"dependencies":["scipy","statsmodels","numpy","pandas"]}

Statistics Reference

Pick the test by answering three questions: what is the unit of comparison (one sample, two samples, k samples, paired, or association), is the variable continuous or categorical, and do parametric assumptions hold. The rest is mechanical: correct call, correct effect size, correct correction.

1. Question to method

Question	Continuous, parametric OK	Continuous, non-parametric	Categorical
One sample vs known value	scipy.stats.ttest_1samp	scipy.stats.wilcoxon (vs median)	scipy.stats.binomtest
Two independent samples	scipy.stats.ttest_ind (Welch if var unequal)	scipy.stats.mannwhitneyu	scipy.stats.chi2_contingency / fisher_exact
Paired samples	scipy.stats.ttest_rel	scipy.stats.wilcoxon	McNemar via statsmodels.stats.contingency_tables.mcnemar
3+ independent groups	scipy.stats.f_oneway (one-way ANOVA)	scipy.stats.kruskal	chi2_contingency on full table
Association, two continuous	scipy.stats.pearsonr	spearmanr, kendalltau	—
Continuous outcome from predictors	statsmodels.api.OLS	rank-transform OLS	logistic via statsmodels.api.Logit
CI for mean	scipy.stats.t.interval or bootstrap	scipy.stats.bootstrap	statsmodels.stats.proportion.proportion_confint

scipy.stats.ttest_ind(..., equal_var=False) is Welch's t-test — prefer it unless variances are known equal. chi2_contingency requires expected counts

= 5 in each cell; otherwise use fisher_exact (2x2) or scipy.stats.fisher_exact extended via permutation for larger tables.

2. Descriptives and confidence intervals

Quantity	Call	Notes
Location	np.mean, np.median, scipy.stats.trim_mean(x, 0.1)	Median for skewed / outlier-heavy
Spread	np.std(ddof=1), scipy.stats.iqr, scipy.stats.median_abs_deviation	ddof=1 for sample SD
Shape	scipy.stats.skew, scipy.stats.kurtosis(fisher=True)	Fisher kurtosis: 0 = normal
95% CI for mean	scipy.stats.t.interval(0.95, n-1, loc=m, scale=se)	se = std / sqrt(n)
95% CI, non-parametric	scipy.stats.bootstrap((x,), np.mean, confidence_level=0.95)	Default 9999 resamples
95% CI for proportion	statsmodels.stats.proportion.proportion_confint(k, n, method='wilson')	Wilson > normal-approx for small n

3. Assumption checks before parametric tests

Test	Required assumption	Check	Library call
One/two-sample t-test, paired t-test	Approximate normality of group(s) or differences	Shapiro-Wilk (n <= 5000)	scipy.stats.shapiro(x)
Two-sample t-test (pooled)	Equal variances	Levene (robust) or Bartlett (normal)	scipy.stats.levene(a, b) / bartlett(a, b)
One-way ANOVA	Normal residuals + equal variances	Shapiro on residuals + Levene across groups	as above
Pearson correlation	Bivariate normal, linear	Shapiro on each + scatter inspection	as above
OLS regression	Linearity, normal residuals, homoscedastic, independent	residual plots + Breusch-Pagan + Durbin-Watson	see Section 6

Rule of thumb: with n >= 30 per group, the CLT covers mild non-normality for the t-test. Severe skew or n < 15 — switch to the non-parametric column.

4. Effect sizes

Always report alongside p-values. Magnitude, not just existence.

Test family	Effect size	Formula / call	Small / Medium / Large
Two-sample t	Cohen's d	(m1 - m2) / pooled_sd	0.2 / 0.5 / 0.8
Mann-Whitney	Cliff's delta	2 * U / (n1 * n2) - 1 from mannwhitneyu	0.15 / 0.33 / 0.47
Paired t	Cohen's d_z	mean(diff) / sd(diff)	0.2 / 0.5 / 0.8
One-way ANOVA	eta-squared	SS_between / SS_total	0.01 / 0.06 / 0.14
Pearson r	r itself	pearsonr(x, y).statistic	0.1 / 0.3 / 0.5
Chi-square 2x2	phi	sqrt(chi2 / n)	0.1 / 0.3 / 0.5
Chi-square RxC	Cramer's V	sqrt(chi2 / (n * min(r-1, c-1)))	0.1 / 0.3 / 0.5
OLS	R-squared, adjusted R-squared	model.rsquared, model.rsquared_adj	context-dependent

5. Test call signatures and what to pull from the result

import numpy as np
from scipy import stats

# Two-sample t-test (Welch by default — safe choice)
res = stats.ttest_ind(a, b, equal_var=False)
t_stat, p = res.statistic, res.pvalue
# Cohen's d (pooled SD)
sd_p = np.sqrt(((len(a)-1)*a.var(ddof=1) + (len(b)-1)*b.var(ddof=1))
               / (len(a) + len(b) - 2))
d = (a.mean() - b.mean()) / sd_p

# Mann-Whitney U (non-parametric two-sample)
u, p = stats.mannwhitneyu(a, b, alternative="two-sided")
cliffs_delta = 2 * u / (len(a) * len(b)) - 1

# Paired test
t_stat, p = stats.ttest_rel(before, after)         # parametric
w, p     = stats.wilcoxon(before, after)            # non-parametric

# One-way ANOVA + Kruskal-Wallis
f, p = stats.f_oneway(g1, g2, g3)
h, p = stats.kruskal(g1, g2, g3)

# Categorical
chi2, p, dof, expected = stats.chi2_contingency(table)   # >= 2x2
odds, p = stats.fisher_exact(table_2x2)                   # 2x2, small n

# Correlations
r,   p = stats.pearsonr(x, y)
rho, p = stats.spearmanr(x, y)
tau, p = stats.kendalltau(x, y)

Every scipy.stats test returns a named tuple with .statistic and .pvalue. Effect size is rarely returned — compute it.

6. OLS regression with residual diagnostics

import numpy as np
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson

X = sm.add_constant(predictors_df)            # adds intercept column
model = sm.OLS(y, X).fit()
print(model.summary())                         # coefs, SE, t, p, R2, F

resid  = model.resid
fitted = model.fittedvalues

# Diagnostics
shapiro_p = stats.shapiro(resid).pvalue                    # normality
bp_lm, bp_p, _, _ = het_breuschpagan(resid, model.model.exog)  # homoscedast.
dw = durbin_watson(resid)                                  # ~2 = independent
# Multicollinearity (per predictor)
from statsmodels.stats.outliers_influence import variance_inflation_factor
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

Red flags: shapiro_p < 0.05 (non-normal residuals — consider transform), bp_p < 0.05 (heteroscedastic — use model.get_robustcov_results('HC3')), dw < 1.5 or > 2.5 (autocorrelation), VIF > 10 (multicollinearity).

7. Multiple-comparison correction

When running k tests, control family-wise error or false-discovery rate:

from statsmodels.stats.multitest import multipletests

pvals = [0.001, 0.01, 0.03, 0.04, 0.20]
reject_bonf, p_bonf, _, _ = multipletests(pvals, alpha=0.05, method="bonferroni")
reject_bh,   p_bh,   _, _ = multipletests(pvals, alpha=0.05, method="fdr_bh")

Use Bonferroni when the k tests are few and you need strict FWER. Use Benjamini-Hochberg (fdr_bh) for screening many hypotheses (e.g. column-wise group differences) where some discoveries-per-list is acceptable.

8. A priori power and sample size

from statsmodels.stats.power import TTestIndPower

# How many per group to detect d=0.5 at alpha=0.05, power=0.80?
analysis = TTestIndPower()
n = analysis.solve_power(effect_size=0.5, alpha=0.05, power=0.80, ratio=1.0)
# What effect could we detect with n=50 per group?
d_min = analysis.solve_power(effect_size=None, nobs1=50, alpha=0.05, power=0.80)

Same pattern for FTestAnovaPower, NormalIndPower. Pre-study planning only; post-hoc "observed power" from a non-significant result is uninformative.

Pitfalls

• Reporting p without an effect size — significance does not imply magnitude.
• Using Student's t-test (equal_var=True) without a Levene check.
• Running pairwise t-tests across k groups without correction.
• Pearson r on non-linear or outlier-driven data — use Spearman / Kendall.
• Chi-square on a sparse contingency table (expected < 5) — use Fisher exact.
• 3-sigma outlier rules on heavy-tailed data — outliers inflate the SD they are flagged against.
• Bootstrapping a paired statistic by resampling the two arrays independently — resample paired indices instead.
• Treating Shapiro-Wilk failures as a fatal blocker at large n; the test is over-powered for n > ~5000 and rejects trivial deviations.