| name | statistical-testing |
| description | Select and run the correct hypothesis test based on data properties. Covers parametric/non-parametric tests, effect sizes, and multiple comparison correction.
|
| tags | ["statistics","hypothesis-testing","scipy","effect-size","multiple-comparisons","power-analysis"] |
| version | 1.0.0 |
| authors | [{"name":"awesome-rosetta-skills contributors","github":"@awesome-rosetta-skills"}] |
| license | MIT |
| platforms | ["claude-code","codex","gemini-cli","cursor"] |
| dependencies | {"python":["scipy>=1.10.0","numpy>=1.24.0","pandas>=2.0.0","statsmodels>=0.14.0","pingouin>=0.5.3"]} |
| last_updated | 2026-03-17 |
Statistical Testing
本整合包适配（ML/时序实验小样本规则）：
- n=3 种子不要报 p 值——检验功效不足，任何"显著"都不可信；报均值±标准差
（
mean ± std，不要写"方差"除非真的报的是 variance）或极差。这是硬策略：
run_comparison 在最小组 n<5 时直接返回描述统计（p_value / significant 置
None），不会返回可被误写进论文的 p 值。
- ≥5 种子才考虑显著性检验；同种子/同数据集天然配对，用配对检验
（小样本首选 Wilcoxon signed-rank；分布近似正态可用 paired t-test）。
- 时序预测特有注意：不同预测长度（96/192/336/720）与不同数据集的误差不可
合并成一个总体做检验（分布不同质）；按（数据集, horizon）分层比较。滚动
评测窗口之间高度自相关，不满足独立样本假设，不要把窗口当独立观测灌进 t 检验。
- 多重比较：数据集 × horizon × baseline × 指标的全组合检验必须做
Holm 或 Benjamini-Hochberg FDR 校正，否则只对预先指定的主对比做检验。
- 实验级比较台账（comparison ledger）：校正必须覆盖整个实验中做过的全部检验，
不是每次调用
correct_pvalues 时手头恰好有的那几个 p 值。在
research_run/<slug>/stage4_experiments/comparison_ledger.md 维护一张累计表
（比较对象 / 数据集 / horizon / 指标 / 原始 p 值 / 所属比较族），投稿前对每个
比较族整体重跑一次校正；表里没登记的检验不得写进论文。
- 两组独立比较默认 Welch：不要用"先 Levene 检验再选 equal_var"的两步法
（膨胀错误率）；
run_comparison 已默认 equal_var=False。
Choose and execute the right statistical test for your data. This skill covers
normality checks, parametric and non-parametric tests, effect size computation,
and multiple comparison correction.
Decision Tree
| Scenario | Groups | Normal? | Recommended Test |
|---|
| Compare means | 2 independent | Yes | Independent t-test |
| Compare means | 2 independent | No | Mann-Whitney U |
| Compare means | 2 paired | Yes | Paired t-test |
| Compare means | 2 paired | No | Wilcoxon signed-rank |
| Compare means | ≥3 independent | Yes | One-way ANOVA |
| Compare means | ≥3 independent | No | Kruskal-Wallis |
| Compare means | ≥3 repeated | Yes | Repeated-measures ANOVA |
| Compare means | ≥3 repeated | No | Friedman test |
| Association | 2 categorical | Expected ≥5 | Chi-square |
| Association | 2 categorical | Expected <5 | Fisher exact |
| Correlation | 2 continuous | Both normal | Pearson r |
| Correlation | 2 continuous | Non-normal | Spearman ρ |
Setup
pip install scipy numpy pandas statsmodels pingouin
Core Implementation
"""
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
Example 1 — Comparing Three Treatment Groups
import numpy as np
from statistical_testing import run_comparison, correct_pvalues, power_analysis_ttest
rng = np.random.default_rng(42)
control = rng.exponential(scale=5, size=40)
treatment1 = rng.exponential(scale=7, size=38)
treatment2 = rng.exponential(scale=9, size=42)
result = run_comparison(control, treatment1, treatment2, alpha=0.05)
from scipy.stats import mannwhitneyu
pairs = [
("control vs t1", mannwhitneyu(control, treatment1, alternative="two-sided").pvalue),
("control vs t2", mannwhitneyu(control, treatment2, alternative="two-sided").pvalue),
("t1 vs t2", mannwhitneyu(treatment1, treatment2, alternative="two-sided").pvalue),
]
labels, raw_p = zip(*pairs)
correction_df = correct_pvalues(list(raw_p), method="fdr_bh")
correction_df.index = labels
print(correction_df)
Example 2 — Paired Comparison + A-Priori Power Planning
import numpy as np
from statistical_testing import test_normality, run_comparison, power_analysis_ttest
rng = np.random.default_rng(0)
pre = rng.normal(loc=100, scale=15, size=30)
post = pre + rng.normal(loc=8, scale=10, size=30)
diff = post - pre
norm_result = test_normality(diff)
result = run_comparison(pre, post, paired=True)
d_planned = 0.5
needed = power_analysis_ttest(effect_size=d_planned, alpha=0.05, power=0.90)
print(f"N needed for 90% power at d={d_planned}: {int(needed['required_n_per_group'])}")
Categorical Data Example
import numpy as np
from statistical_testing import chi_square_test_with_effect, correct_pvalues
table = np.array([[45, 15],
[30, 30]])
result = chi_square_test_with_effect(table)
print(result)
multi_table = np.array([[50, 20, 10],
[40, 25, 15],
[35, 30, 20]])
result3 = chi_square_test_with_effect(multi_table)
print(result3)
Quick Reference: Effect Size Benchmarks
| Measure | Small | Medium | Large |
|---|
| Cohen's d | 0.2 | 0.5 | 0.8 |
| η² (eta-squared) | 0.01 | 0.06 | 0.14 |
| Cramér's V (2×2) | 0.1 | 0.3 | 0.5 |
| Pearson r | 0.1 | 0.3 | 0.5 |
Multiple Comparison Methods
| Method | Controls | Best For |
|---|
| Bonferroni | Family-wise error rate (FWER) | Few comparisons, strict control |
| Holm | FWER (less conservative) | General use |
| Benjamini-Hochberg (FDR) | False discovery rate | Many comparisons (genomics, etc.) |
| Benjamini-Yekutieli | FDR under dependence | Correlated tests |
Common Pitfalls
- Assuming normality without testing: Always run normality checks, especially for n < 30.
- Ignoring equal-variance assumption: Run Levene's test before independent t-test.
- Multiple comparisons inflation: Any time you run ≥3 tests on the same dataset, apply correction.
- Over-relying on p < 0.05: Always report effect sizes alongside p-values.
- Wrong test for paired data: Pre/post measurements are paired; use paired tests.
Environment Variables
No API keys required. All computation is local.
export STATS_RANDOM_SEED=42
