"""
statistical_testing.py
Complete statistical testing toolkit with automatic test selection.
"""
import warnings
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats import multitest
from statsmodels.stats.power import TTestIndPower, FTestAnovaPower
from typing import Union, Optional, Tuple, List, Dict, Any
def test_normality(
data: np.ndarray,
alpha: float = 0.05,
verbose: bool = True,
) -> Dict[str, Any]:
"""
Run Shapiro-Wilk, Kolmogorov-Smirnov, and D'Agostino-Pearson normality tests.
Parameters
----------
data : array-like
1-D numeric data.
alpha : float
Significance level.
verbose : bool
Print a summary.
Returns
-------
dict with keys: shapiro, ks, dagostino, is_normal
"""
data = np.asarray(data, dtype=float)
data = data[~np.isnan(data)]
n = len(data)
results: Dict[str, Any] = {}
if n < 5000:
stat_sw, p_sw = stats.shapiro(data)
results["shapiro"] = {"statistic": stat_sw, "p_value": p_sw, "normal": p_sw > alpha}
else:
results["shapiro"] = {"statistic": None, "p_value": None, "normal": None,
"note": "n>=5000; skipped Shapiro-Wilk"}
if n >= 4:
try:
from statsmodels.stats.diagnostic import lilliefors
stat_ks, p_ks = lilliefors(data, dist="norm")
results["ks"] = {"statistic": stat_ks, "p_value": p_ks, "normal": p_ks > alpha,
"note": "Lilliefors-corrected KS"}
except ImportError:
results["ks"] = {"statistic": None, "p_value": None, "normal": None,
"note": "statsmodels missing; skipped KS (do not substitute "
"a plain kstest with fitted parameters)"}
else:
results["ks"] = {"statistic": None, "p_value": None, "normal": None,
"note": "n<4; skipped KS"}
if n >= 20:
stat_da, p_da = stats.normaltest(data)
results["dagostino"] = {"statistic": stat_da, "p_value": p_da, "normal": p_da > alpha}
else:
results["dagostino"] = {"statistic": None, "p_value": None, "normal": None,
"note": "n<20; skipped D'Agostino"}
votes = [v["normal"] for v in results.values() if isinstance(v, dict) and v["normal"] is not None]
results["is_normal"] = sum(votes) >= (len(votes) / 2)
results["n"] = n
if verbose:
print(f"Normality tests (n={n}, alpha={alpha})")
for test_name, res in results.items():
if isinstance(res, dict) and "p_value" in res and res["p_value"] is not None:
verdict = "NORMAL" if res["normal"] else "NON-NORMAL"
print(f" {test_name:15s}: stat={res['statistic']:.4f}, p={res['p_value']:.4f} [{verdict}]")
print(f" Consensus: {'NORMAL' if results['is_normal'] else 'NON-NORMAL'}\n")
return results
def cohens_d(group1: np.ndarray, group2: np.ndarray) -> float:
"""Cohen's d for two independent groups."""
g1, g2 = np.asarray(group1), np.asarray(group2)
n1, n2 = len(g1), len(g2)
pooled_sd = np.sqrt(((n1 - 1) * np.var(g1, ddof=1) + (n2 - 1) * np.var(g2, ddof=1)) / (n1 + n2 - 2))
return (np.mean(g1) - np.mean(g2)) / pooled_sd
def eta_squared(f_statistic: float, df_between: int, df_within: int) -> float:
"""Eta-squared from F-statistic for one-way ANOVA."""
ss_between = f_statistic * df_between
ss_total = ss_between + df_within
return ss_between / ss_total
def cramers_v(contingency_table: np.ndarray) -> float:
"""Cramér's V for association in a contingency table."""
chi2 = stats.chi2_contingency(contingency_table, correction=False)[0]
n = contingency_table.sum()
min_dim = min(contingency_table.shape) - 1
return np.sqrt(chi2 / (n * min_dim))
def interpret_effect_size(d: float, measure: str = "cohens_d") -> str:
"""Return a qualitative label for an effect size."""
thresholds = {
"cohens_d": [(0.2, "small"), (0.5, "medium"), (0.8, "large")],
"eta_squared": [(0.01, "small"), (0.06, "medium"), (0.14, "large")],
"cramers_v": [(0.1, "small"), (0.3, "medium"), (0.5, "large")],
}
d = abs(d)
for threshold, label in thresholds.get(measure, []):
if d < threshold:
return label
return "large"
def run_comparison(
*groups: np.ndarray,
paired: bool = False,
alpha: float = 0.05,
verbose: bool = True,
) -> Dict[str, Any]:
"""
Automatically select and run the appropriate comparison test.
Parameters
----------
*groups : array-like
Two or more numeric arrays (one per group).
paired : bool
Whether observations are paired / repeated measures.
alpha : float
Significance level.
verbose : bool
Print result summary.
Returns
-------
dict with keys: test_name, statistic, p_value, effect_size, significant
"""
groups = [np.asarray(g, dtype=float) for g in groups]
k = len(groups)
if k < 2:
raise ValueError("Need at least 2 groups.")
min_n = min(len(g) for g in groups)
if min_n < 5:
note = (f"smallest group has n={min_n} < 5: statistical test not performed "
"(pack policy: no p-values below n=5 — report mean ± std or range, "
"not significance claims).")
result = {
"test_name": None,
"statistic": None,
"p_value": None,
"significant": None,
"alpha": alpha,
"effect_size": {},
"descriptives": [
{"n": int(len(g)),
"mean": round(float(np.mean(g)), 4),
"std": round(float(np.std(g, ddof=1)), 4) if len(g) > 1 else None,
"min": round(float(np.min(g)), 4),
"max": round(float(np.max(g)), 4)}
for g in groups
],
"note": note,
}
if verbose:
print(f"No test run: {note}")
for i, d in enumerate(result["descriptives"]):
print(f" group {i}: n={d['n']}, mean={d['mean']}, std={d['std']}, "
f"range=[{d['min']}, {d['max']}]")
print()
return result
all_normal = all(test_normality(g, alpha=alpha, verbose=False)["is_normal"] for g in groups)
result: Dict[str, Any] = {}
if k == 2:
g1, g2 = groups[0], groups[1]
if paired:
if all_normal:
stat, p = stats.ttest_rel(g1, g2)
result["test_name"] = "Paired t-test"
d = np.mean(g1 - g2) / np.std(g1 - g2, ddof=1)
result["effect_size"] = {"cohens_d": round(d, 4)}
else:
stat, p = stats.wilcoxon(g1, g2)
result["test_name"] = "Wilcoxon signed-rank"
diff = g1 - g2
diff = diff[diff != 0]
ranks = stats.rankdata(np.abs(diff))
t_plus = float(ranks[diff > 0].sum())
t_minus = float(ranks[diff < 0].sum())
rb = (t_plus - t_minus) / (t_plus + t_minus)
result["effect_size"] = {"rank_biserial": round(rb, 4)}
else:
if all_normal:
stat, p = stats.ttest_ind(g1, g2, equal_var=False)
result["test_name"] = "Welch's t-test"
d = cohens_d(g1, g2)
result["effect_size"] = {"cohens_d": round(d, 4),
"interpretation": interpret_effect_size(d)}
else:
stat, p = stats.mannwhitneyu(g1, g2, alternative="two-sided")
result["test_name"] = "Mann-Whitney U"
n1, n2 = len(g1), len(g2)
rb = 1 - (2 * stat) / (n1 * n2)
result["effect_size"] = {"rank_biserial_r": round(rb, 4)}
else:
if paired:
if all_normal:
try:
import pingouin as pg
df_long = pd.DataFrame({
"value": np.concatenate(groups),
"group": np.repeat(np.arange(k), [len(g) for g in groups]),
"subject": np.tile(np.arange(len(groups[0])), k),
})
aov = pg.rm_anova(data=df_long, dv="value", within="group", subject="subject")
stat = aov["F"].iloc[0]
p = aov["p-unc"].iloc[0]
result["test_name"] = "Repeated-measures ANOVA (pingouin)"
result["effect_size"] = {"eta_squared": round(aov["np2"].iloc[0], 4)}
except ImportError:
warnings.warn("pingouin not installed; falling back to Friedman test.")
stat, p = stats.friedmanchisquare(*groups)
result["test_name"] = "Friedman test (fallback)"
result["effect_size"] = {}
else:
stat, p = stats.friedmanchisquare(*groups)
result["test_name"] = "Friedman test"
result["effect_size"] = {}
else:
if all_normal:
stat, p = stats.f_oneway(*groups)
result["test_name"] = "One-way ANOVA"
total_n = sum(len(g) for g in groups)
df_between = k - 1
df_within = total_n - k
es = eta_squared(stat, df_between, df_within)
result["effect_size"] = {"eta_squared": round(es, 4),
"interpretation": interpret_effect_size(es, "eta_squared")}
else:
stat, p = stats.kruskal(*groups)
result["test_name"] = "Kruskal-Wallis H"
result["effect_size"] = {}
result["statistic"] = round(float(stat), 4)
result["p_value"] = round(float(p), 6)
result["significant"] = p < alpha
result["alpha"] = alpha
if verbose:
print(f"Test: {result['test_name']}")
print(f" Statistic = {result['statistic']}, p = {result['p_value']}")
print(f" Significant at alpha={alpha}: {result['significant']}")
if result.get("effect_size"):
print(f" Effect size: {result['effect_size']}")
print()
return result
def correct_pvalues(
p_values: List[float],
method: str = "fdr_bh",
alpha: float = 0.05,
) -> pd.DataFrame:
"""
Apply multiple comparison correction.
Parameters
----------
p_values : list of float
method : str
One of: bonferroni, holm, fdr_bh (Benjamini-Hochberg),
fdr_by (Benjamini-Yekutieli), sidak
alpha : float
Returns
-------
DataFrame with original p-values, corrected p-values, and rejection flags.
"""
reject, p_corrected, _, _ = multitest.multipletests(p_values, alpha=alpha, method=method)
return pd.DataFrame({
"p_original": p_values,
"p_corrected": p_corrected,
"reject_H0": reject,
"method": method,
})
def power_analysis_ttest(
effect_size: float = 0.5,
alpha: float = 0.05,
power: float = 0.80,
n_per_group: Optional[int] = None,
) -> Dict[str, float]:
"""
Two-sample t-test power analysis — for **a priori** planning only (choose n
before running the experiment, using an effect size from pilot data or the
literature). Do NOT compute "achieved/post-hoc power" from the observed effect
size of the same experiment: it is a monotone transform of the p-value and adds
no information. Note this solver is for *independent* two-sample designs; a
paired design has different power characteristics (use
statsmodels TTestPower on the difference scores instead).
Provide any three of (effect_size, alpha, power, n_per_group) to solve for the fourth.
"""
analysis = TTestIndPower()
if n_per_group is None:
n = analysis.solve_power(effect_size=effect_size, alpha=alpha, power=power, ratio=1.0)
return {"required_n_per_group": np.ceil(n), "effect_size": effect_size,
"alpha": alpha, "power": power}
else:
achieved_power = analysis.solve_power(effect_size=effect_size, alpha=alpha,
nobs1=n_per_group, ratio=1.0)
return {"n_per_group": n_per_group, "effect_size": effect_size,
"alpha": alpha, "achieved_power": round(achieved_power, 4)}
def chi_square_test_with_effect(
contingency: np.ndarray,
alpha: float = 0.05,
) -> Dict[str, Any]:
"""Chi-square test with Cramér's V effect size."""
chi2, p, dof, expected = stats.chi2_contingency(contingency)
v = cramers_v(contingency)
result = {
"test_name": "Chi-square" if expected.min() >= 5 else "Fisher exact (recommended)",
"chi2": round(chi2, 4),
"p_value": round(p, 6),
"dof": dof,
"cramers_v": round(v, 4),
"significant": p < alpha,
}
if expected.min() < 5:
if contingency.shape == (2, 2):
_, p_fisher = stats.fisher_exact(contingency)
result["fisher_p"] = round(p_fisher, 6)
return result
