statistical-analysis

统计分析

概述

进行假设检验（t 检验、方差分析、卡方检验）、回归、相关性和贝叶斯分析，并系统地检查假设、计算效应量以及按 APA 格式报告结果。目标是完成一项经得起审稿人严格审查的分析：选对检验、验证假设、如实报告效应量，并提供完整的说明。

何时使用此技能

在以下情况下使用此技能：

• 进行统计假设检验（t 检验、方差分析、卡方检验、非参数检验）
• 进行回归或相关性分析
• 进行贝叶斯统计分析
• 检查统计假设和诊断指标
• 计算效应量并进行功效分析
• 按 APA 格式报告统计结果
• 分析用于研究的实验数据或观察数据

---

安装

使用 uv 安装此技能使用的库。在生产环境中固定版本；探索时可以不固定版本。

# Core frequentist stack (Python 3.10+; 3.12+ recommended for latest SciPy/ArviZ)
uv pip install "pingouin>=0.6" "scipy>=1.11" "statsmodels>=0.14.6" pandas matplotlib seaborn

# Bayesian modeling (PyMC 5 + ArviZ)
uv pip install "pymc>=5.0" "arviz>=1.0"

兼容性说明（已针对 pingouin 0.6.1、statsmodels 0.14.6、arviz 1.2 和 2026 年进行验证）：

• Pingouin 0.6.0 重命名了输出列，移除了特殊字符：p_val、cohen_d、CI95、p_unc（在 0.5.x 中分别为 p-val、cohen-d、CI95%、p-unc）。以下示例使用当前名称；如果仍使用 0.5.x，请使用带连字符的形式。
• statsmodels + SciPy：使用 statsmodels>=0.14.6 和 scipy>=1.11，以避免在 SciPy 1.16+ 中出现 _lazywhere 导入错误。
• ArviZ 1.x：az.summary() 现在默认使用 89% 区间（eti89 列），区间宽度参数为 ci_prob（而非 hdi_prob）。若要报告常规的 95% 可信区间，请传入 az.summary(trace, ci_prob=0.95)。
• Pingouin 不再支持单侧贝叶斯因子：pg.ttest(..., alternative='greater') 会静默地移除 BF10 列，而 pg.bayesfactor_ttest 在使用单侧备择假设时会报错。进行单侧贝叶斯检验时，请直接使用 PyMC（计算方向性假设的后验概率），或使用 JASP/R 的 BayesFactor。

关于特定模型的 API（OLS、GLM、ARIMA），请参阅 statsmodels 技能。关于 PyMC 工作流，请参阅 pymc 技能。

---

分析工作流

每项可靠的分析都遵循相同的流程。跳过步骤会导致分析最终被撤回，因此请按顺序完成每一步，并说明你在每一步做了什么。

1. 先明确问题，再接触数据。 说明假设、结果变量和预测变量，以及研究设计（独立样本还是配对样本、组数）。现在就确定计划使用的检验——看过结果后再选择检验属于 p 值操纵，即便并非有意为之。
2. 检查数据。 按组检查：样本量 n、均值、标准差、中央値、缺失值。进行任何检验前，先绘制原始数据图（直方图或箱线图）。组间样本量不等、缺失值、地板／天花板效应以及离群值都会影响适用的检验——应向用户指出这些情况，而不是悄悄绕开。
3. 选择检验，可使用下方的快速参考，或查阅 references/test_selection_guide.md，了解基础设计以外的情况（计数、事件发生时间、信度、析因设计）。
4. 检查假设，使用 scripts/assumption_checks.py。如果某项假设不成立，请改用补救检验（见下表），并报告原计划和变更。
5. 运行检验，并始终同时计算效应量——p 值说明效应是否存在；效应量则说明效应是否值得关注。
6. 报告结果，使用下方的 APA 模板，其中包括描述性统计、精确统计量、带置信区间的效应量，以及进行过的假设检查。

如果用户只需要某一个步骤（例如“我需要多少名参与者？”），可直接跳到相应部分，但仍需确认计算所依据的设计假设。

---

检验选择指南

快速参考：选择合适的检验

请参阅 references/test_selection_guide.md，了解全面指南（计数、生存分析、可靠性、因子设计）。快速参考：

比较两组：

• 独立、连续、正态 → 独立样本 t 检验
• 独立、连续、非正态 → Mann-Whitney U 检验
• 配对、连续、正态 → 配对样本 t 检验
• 配对、连续、非正态 → Wilcoxon 符号秩检验
• 二元结果 → 卡方检验或 Fisher 精确检验

比较 3 组或以上：

• 独立、连续、正态 → 单因素 ANOVA
• 独立、连续、非正态 → Kruskal-Wallis 检验
• 配对、连续、正态 → 重复测量 ANOVA
• 配对、连续、非正态 → Friedman 检验

关系：

• 两个连续变量 → Pearson（正态）或 Spearman 相关（非正态）
• 连续结果变量与预测变量 → 线性回归
• 二元结果变量与预测变量 → Logistic 回归

贝叶斯替代方法： 所有检验都有贝叶斯版本，可直接对假设给出概率陈述、用贝叶斯因子量化证据，并支持零假设。请参阅 references/bayesian_statistics.md。

---

假设检查

解释检验结果前，务必检查假设，并报告检查结果——审稿人会关注这一点。

请使用捆绑的 scripts/assumption_checks.py 模块。在技能目录（skills/statistical-analysis/）中运行 Python，或将 scripts/ 添加到 sys.path：

from assumption_checks import comprehensive_assumption_check

# Outliers + normality (per group) + homogeneity of variance, with plots
results = comprehensive_assumption_check(
    data=df,
    value_col='score',
    group_col='group',  # Optional: for group comparisons
    alpha=0.05
)

如需有针对性地检查，请导入单个函数：

from assumption_checks import (
    check_normality,                # Shapiro-Wilk + Q-Q plot + histogram
    check_normality_per_group,
    check_homogeneity_of_variance,  # Levene's test + box plots
    check_linearity,                # scatter + residual plot for simple regression
    check_regression_diagnostics,   # full OLS diagnostics (see Regression below)
    detect_outliers                 # IQR or z-score methods
)

result = check_normality(data=df['score'], name='Test Score', alpha=0.05, plot=True)
print(result['interpretation'])
print(result['recommendation'])

假设不满足时的处理方法

正态性不满足：

• 轻微违反，且每组 n > 30 → 使用参数检验（稳健）
• 中度违反 → 使用非参数替代方法
• 严重违反 → 转换数据或使用非参数检验

方差齐性不满足：

• 对于 t 检验 → 使用 Welch's t 检验（pg.ttest 会在 correction='auto' 时自动应用）
• 对于 ANOVA → 使用 Welch's ANOVA（pg.welch_anova）或 Brown-Forsythe 检验
• 对于回归 → 使用稳健标准误或加权最小二乘法

线性假设不满足（回归）：

• 添加多项式项、转换变量，或使用非线性模型 / GAM

随着 n 增大，正式检验会变得过于敏感：当 n ≥ 100 时，应更多依赖 Q-Q 图，而不是 Shapiro-Wilk p 值。请参阅 references/assumptions_and_diagnostics.md，了解全面指南。

---

运行统计检验

主要库：

• pingouin：易用的检验，默认返回效应量——标准检验优先使用它
• scipy.stats：核心统计检验
• statsmodels：回归、诊断、功效分析
• pymc + arviz：贝叶斯建模与诊断

完整报告 t 检验结果

import pingouin as pg

# correction='auto' applies Welch's correction when variances are unequal
result = pg.ttest(group_a, group_b, correction='auto')

# Pingouin >= 0.6 column names
t_stat = result['T'].values[0]
df = result['dof'].values[0]
p_value = result['p_val'].values[0]
cohens_d = result['cohen_d'].values[0]
ci_lower, ci_upper = result['CI95'].values[0]  # CI for the mean difference

print(f"t({df:.0f}) = {t_stat:.2f}, p = {p_value:.3f}, d = {cohens_d:.2f}")

使用事后检验的 ANOVA

import pingouin as pg

aov = pg.anova(dv='score', between='group', data=df, detailed=True)
print(aov)

# Effect size: partial eta-squared
eta_p2 = aov['np2'].values[0]

# If significant, conduct post-hoc tests (Tukey HSD controls family-wise error)
if aov['p_unc'].values[0] < 0.05:
    posthoc = pg.pairwise_tukey(dv='score', between='group', data=df)
    print(posthoc)  # includes Hedges' g per pair

带诊断的线性回归

import statsmodels.api as sm
from assumption_checks import check_regression_diagnostics

X = sm.add_constant(X_predictors)  # Add intercept
model = sm.OLS(y, X).fit()
print(model.summary())

# 4-panel residual plot + Shapiro-Wilk, Breusch-Pagan, Durbin-Watson, VIF
diag = check_regression_diagnostics(model)
print(diag['interpretation'])
print(diag['vif'])

# If heteroscedasticity was flagged, report robust standard errors instead
robust = model.get_robustcov_results('HC3')

贝叶斯 t 检验

import pymc as pm
import arviz as az
import numpy as np

with pm.Model() as model:
    # Priors
    mu1 = pm.Normal('mu_group1', mu=0, sigma=10)
    mu2 = pm.Normal('mu_group2', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=10)

    # Likelihood
    y1 = pm.Normal('y1', mu=mu1, sigma=sigma, observed=group_a)
    y2 = pm.Normal('y2', mu=mu2, sigma=sigma, observed=group_b)

    # Derived quantity
    diff = pm.Deterministic('difference', mu1 - mu2)

    trace = pm.sample(2000, tune=1000)

# ArviZ 1.x defaults to 89% intervals; request 95% explicitly for reporting
print(az.summary(trace, var_names=['difference'], ci_prob=0.95))

# Direct probability statement (this is what one-sided questions become)
prob_greater = np.mean(trace.posterior['difference'].values > 0)
print(f"P(mu1 > mu2 | data) = {prob_greater:.3f}")

# ArviZ 1.x removed az.plot_posterior; use plot_dist (on 0.x, plot_posterior still works)
az.plot_dist(trace, var_names=['difference'], ci_prob=0.95)

先验应根据数据进行缩放（例如，结果变量的 SD 接近 10 时，sigma=10 较合适；可参考观测到的 SD），并在报告中说明先验。

---

效应量

效应量用于量化大小；p 值仅表示是否存在效应。 每项检验都应报告一个效应量。完整指南请参阅 references/effect_sizes_and_power.md。

快速参考：常见效应量

检验	效应量	小	中	大
T 检验	Cohen's d	0.20	0.50	0.80
ANOVA	η²_p	0.01	0.06	0.14
相关	r	0.10	0.30	0.50
回归	R²	0.02	0.13	0.26
卡方检验	Cramér's V	0.07	0.21	0.35

这些基准是惯例，并非定律——“小”效应可能极为重要（如药物副作用），“大”效应也可能微不足道。请结合具体情境进行解释。

计算效应量

Pingouin 会在检验结果中返回效应量（pg.ttest 返回 cohen_d，pg.anova 返回 np2，pg.pairwise_tukey 返回 hedges；pg.corr 中的 r 本身就是效应量）。

效应量的置信区间

报告效应量的置信区间（CI），以展示其精确度。使用 pg.compute_esci（注意：pg.compute_effsize_from_t 仅返回点估计，不会返回 CI）：

import pingouin as pg

d = pg.compute_effsize(group_a, group_b, eftype='cohen')
ci_lower, ci_upper = pg.compute_esci(stat=d, nx=len(group_a), ny=len(group_b),
                                     eftype='cohen', confidence=0.95)
print(f"d = {d:.2f}, 95% CI [{ci_lower:.2f}, {ci_upper:.2f}]")

---

功效分析

先验功效分析（研究规划）

在收集数据前确定所需的样本量：

from statsmodels.stats.power import tt_ind_solve_power, FTestAnovaPower

# T-test: What n per group is needed to detect d = 0.5?
n_required = tt_ind_solve_power(
    effect_size=0.5,
    alpha=0.05,
    power=0.80,
    ratio=1.0,
    alternative='two-sided'
)
print(f"Required n per group: {n_required:.0f}")

# One-way ANOVA: What n is needed to detect Cohen's f = 0.25?
# Notes: the parameter is k_groups; effect_size is Cohen's f (f = sqrt(eta2/(1-eta2)));
# and solve_power returns the TOTAL sample size, not n per group.
import math
anova_power = FTestAnovaPower()
n_total = anova_power.solve_power(
    effect_size=0.25,
    k_groups=3,
    alpha=0.05,
    power=0.80
)
print(f"Required total N: {math.ceil(n_total)} ({math.ceil(n_total / 3)} per group)")

敏感性分析（研究后）

确定研究能够检测到多大的效应量：

# With n=50 per group, what effect could we detect at 80% power?
detectable_d = tt_ind_solve_power(
    effect_size=None,  # Solve for this
    nobs1=50,
    alpha=0.05,
    power=0.80,
    ratio=1.0,
    alternative='two-sided'
)
print(f"Study could detect d >= {detectable_d:.2f}")

注意：事后“观察功效”（根据观察到的效应计算功效）是循环论证且具有误导性——它是 p 值的确定性函数。如果研究已经完成，而有人询问功效，请改为进行敏感性分析。

详细指南请参阅 references/effect_sizes_and_power.md。

---

报告结果

请遵循 references/reporting_standards.md 中的 APA 格式。每份报告都需要包含：

1. 描述性统计：所有组别／变量的 M、SD、n
2. 检验统计量：检验名称、统计量、df、精确 p 值（p = .034，而非 p < .05；仅当 p < .001 时使用 p < .001）
3. 效应量：附置信区间
4. 假设检验：进行了哪些检验、结果以及采取的措施
5. 所有计划的分析：包括不显著的结果——省略这些结果属于挑选性报告

报告示例模板

独立样本 T 检验

Group A (n = 48, M = 75.2, SD = 8.5) scored significantly higher than
Group B (n = 52, M = 68.3, SD = 9.2), t(98) = 3.82, p < .001, d = 0.77,
95% CI [0.36, 1.18], two-tailed. Assumptions of normality (Shapiro-Wilk:
Group A W = 0.97, p = .18; Group B W = 0.96, p = .12) and homogeneity
of variance (Levene's F(1, 98) = 1.23, p = .27) were satisfied.

单因素方差分析

A one-way ANOVA revealed a significant main effect of treatment condition
on test scores, F(2, 147) = 8.45, p < .001, η²_p = .10. Post hoc
comparisons using Tukey's HSD indicated that Condition A (M = 78.2,
SD = 7.3) scored significantly higher than Condition B (M = 71.5,
SD = 8.1, p = .002, d = 0.87) and Condition C (M = 70.1, SD = 7.9,
p < .001, d = 1.07). Conditions B and C did not differ significantly
(p = .52, d = 0.18).

多元回归

Multiple linear regression was conducted to predict exam scores from
study hours, prior GPA, and attendance. The overall model was significant,
F(3, 146) = 45.2, p < .001, R² = .48, adjusted R² = .47. Study hours
(B = 1.80, SE = 0.31, β = .35, t = 5.78, p < .001, 95% CI [1.18, 2.42])
and prior GPA (B = 8.52, SE = 1.95, β = .28, t = 4.37, p < .001,
95% CI [4.66, 12.38]) were significant predictors, while attendance was
not (B = 0.15, SE = 0.12, β = .08, t = 1.25, p = .21, 95% CI [-0.09, 0.39]).
Multicollinearity was not a concern (all VIF < 1.5).

贝叶斯分析

A Bayesian independent samples t-test was conducted using weakly
informative priors (Normal(0, 10) for group means). The posterior
distribution indicated that Group A scored higher than Group B
(M_diff = 6.8, 95% credible interval [3.2, 10.4]), with a 99.8%
posterior probability that Group A's mean exceeded Group B's mean.
Convergence diagnostics were satisfactory (all R-hat < 1.01, ESS > 1000).

如果使用了非参数检验，请报告中位数而非均值、U／W／H 统计量，以及基于秩的效应量（例如秩双列相关；pg.mwu 以 RBC 的形式返回该值）。

---

贝叶斯统计

在以下情况下，可考虑使用贝叶斯方法：

• 有需要纳入的先验信息
• 希望直接陈述假设的概率（“效应有 95% 的概率落在此区间内”）
• 样本量较小或采用序贯数据收集（无需对可选停止进行校正）
• 需要量化支持零假设的证据
• 模型较复杂（具有层级结构、存在缺失数据）

有关先验设定、贝叶斯因子、可信区间、层级模型和收敛性检查（R-hat < 1.01、ESS 充足、后验预测检查），请参阅 references/bayesian_statistics.md。

---

附带资源

参考资料（references/）

• test_selection_guide.md：决策树，涵盖组间比较、关系、计数、生存时间、一致性／信度和分类分析
• assumptions_and_diagnostics.md：检查和处理假设违反情况的详细指南
• effect_sizes_and_power.md：效应量的计算、解释和报告；功效分析
• bayesian_statistics.md：先验、贝叶斯因子、可信区间、层级模型、诊断
• reporting_standards.md：APA 格式报告指南及示例

脚本（scripts/）

• assumption_checks.py：自动检查假设并生成可视化图表
  • comprehensive_assumption_check()：一次调用即可检查离群值、正态性和方差齐性
  • check_normality()、check_normality_per_group()：Shapiro-Wilk 检验及 Q-Q 图
  • check_homogeneity_of_variance()：Levene 检验及箱线图
  • check_regression_diagnostics()：对拟合的 OLS 模型生成 4 面板残差图，并进行 Shapiro-Wilk、Breusch-Pagan、Durbin-Watson 和 VIF 检验
  • check_linearity()、detect_outliers()

---

统计完整性

以下做法能确保分析经得起检验。它们之所以重要，是因为最常见的统计问题并非计算错误，而是隐性的灵活操作（反复检验，直到得到理想结果）和选择性报告。

1. 区分验证性分析与探索性分析。 运行分析前说明计划采用的分析方法；将过程中发现的任何内容标记为探索性发现。
2. 不要为了显著性反复尝试。 如果计划采用的检验结果不显著，这就是结果。反复尝试不同的检验、子群或异常值剔除方案，直到 p < .05，会使 p 值失去效力。
3. 进行多重比较校正，适用于成组检验（事后 ANOVA 使用 Tukey HSD；其他检验组使用 Holm 或 Benjamini-Hochberg FDR），并说明采用了哪种校正方法。
4. 结果不显著并不代表没有效应。 如果样本量 n 较小，研究可能只是检验效能不足——进行敏感性分析，或采用贝叶斯分析／等效性检验，实际量化对零假设的支持程度。
5. 统计显著性不等于实际重要性。 样本量 n 较大时，微不足道的效应也可能达到 p < .001。解读时应首先说明效应量。
6. 删除行之前先了解缺失数据。 只有数据完全随机缺失时，列表删除才安全；否则应考虑多重插补，并说明所采取的做法。
7. 确保分析可复现。 设置随机种子，报告基于模拟的方法所用的库版本，并将分析保存在可运行的脚本中。

相关技能

• statistical-power：用于研究规划的样本量和统计检验效能计算。
• statsmodels：Python 统计模型库。
• bgpt-paper-search：通过 BGPT MCP 服务搜索科学论文，并获取从全文研究中提取的结构化实验数据。