面向研究数据的引导式统计分析——检验选择、假设检查、效应量、功效分析、贝叶斯替代方法和 APA 格式报告。适用于用户想要比较组、检验假设、分析实验或调查数据、检
进行假设检验(t 检验、方差分析、卡方检验)、回归、相关和贝叶斯分析,并系统地进行假设检查、效应量计算和 APA 风格报告是一项面向实际任务的技能,主要用于目标是产生一份审稿人无法挑剔的分析:正确的检验、经过验证的假设、诚实的效应量和完整的报告;何时使用此技能;
运行贝叶斯统计分析;检查统计假设和诊断。它将相关步骤、工具调用和结果整理方式集中到统一流程中,帮助使用者更快完成目标并减少重复操作。从功能定位来看,该技能强调把分散的操作要求整理成清晰、可复用的处理流程,使用户能够围绕既定目标快速准备输入、选择执行方式并获得结构化结果。
实际使用前应先确认任务范围、数据来源、运行环境、必要权限和关键参数,再依据技能说明逐步执行;若输入条件不完整,应先补齐信息或采用保守配置,避免因错误假设导致结果偏离需求。执行过程中需要关注工具调用是否成功、接口或依赖是否可用、输出格式是否符合预期,并对异常提示、缺失字段和边界情况进行处理;涉及批量任务时,还应保存进度,避免中断后重复操作。
进行假设检验(t 检验、方差分析、卡方检验)、回归、相关和贝叶斯分析,并系统地进行假设检查、效应量计算和 APA 风格报告。目标是产生一份审稿人无法挑剔的分析:正确的检验、经过验证的假设、诚实的效应量和完整的报告。
在以下情况下使用此技能:
使用 uv 安装此技能中使用的库。生产环境固定版本;非固定安装适合探索。
# 核心频率统计栈(Python 3.10+;推荐 3.12+ 以获得最新的 SciPy/ArviZ) uv pip install "pingouin>=0.6" "scipy>=1.11" "statsmodels>=0.14.6" pandas matplotlib seaborn # 贝叶斯建模(PyMC 5 + ArviZ) uv pip install "pymc>=5.0" "arviz>=1.0"
兼容性说明(已验证 pingouin 0.6.1、statsmodels 0.14.6、arviz 1.2,2026 年):
p_val、cohen_d、CI95、p_unc(之前为 0.5.x 中的 p-val、cohen-d、CI95%、p-unc)。以下示例使用当前名称;如果停留在 0.5.x,请使用连字符形式。statsmodels>=0.14.6 和 scipy>=1.11 以避免 SciPy 1.16+ 上的 _lazywhere 导入错误。az.summary() 现在默认为 89% 区间(eti89 列),宽度参数为 ci_prob(而非 hdi_prob)。要报告常规的 95% 可信区间,请传递 az.summary(trace, ci_prob=0.95)。pg.ttest(..., alternative='greater') 静默丢弃 BF10 列,pg.bayesfactor_ttest 在单侧替代上引发错误。对于单侧贝叶斯检验,直接使用 PyMC(计算方向性假设的后验概率)或 JASP/R 的 BayesFactor。对于模型特定 API(OLS、GLM、ARIMA),参见 statsmodels 技能。对于 PyMC 工作流,参见 pymc 技能。
每项合理的分析都遵循相同的流程。跳过步骤是分析最终被撤稿的原因,因此请按顺序进行,并在每一步说明你所做的工作。
references/test_selection_guide.md。scripts/assumption_checks.py 检查假设。 如果假设不满足,切换到补救性检验(见下表),并同时报告原计划和变更。如果用户只需要一个步骤(例如"我需要多少参与者?"),直接跳到该部分——但仍然确认计算所依赖的设计假设。
参见 references/test_selection_guide.md 获取全面指导(计数、生存、可靠性、因子设计)。快速参考:
比较两组:
比较 3 组及以上:
关系:
贝叶斯替代方法:
所有检验都有贝叶斯版本,提供关于假设的直接概率陈述、量化证据的贝叶斯因子以及支持零假设的能力。参见 references/bayesian_statistics.md。
在解释检验结果之前始终检查假设,并报告检查结果——审稿人会关注这些。
使用捆绑的 scripts/assumption_checks.py 模块。从技能目录(skills/statistical-analysis/)运行 Python,或将 scripts/ 添加到 sys.path:
from assumption_checks import comprehensive_assumption_check
# 异常值 + 正态性(每组)+ 方差齐性,含图表
results = comprehensive_assumption_check(
data=df,
value_col='score',
group_col='group', # 可选:用于组间比较
alpha=0.05
)
针对特定检查,导入单独的函数:
from assumption_checks import (
check_normality, # Shapiro-Wilk + Q-Q 图 + 直方图
check_normality_per_group,
check_homogeneity_of_variance, # Levene 检验 + 箱线图
check_linearity, # 简单回归的散点图 + 残差图
check_regression_diagnostics, # 完整 OLS 诊断(见下面的回归)
detect_outliers # IQR 或 z-score 方法
)
result = check_normality(data=df['score'], name='Test Score', alpha=0.05, plot=True)
print(result['interpretation'])
print(result['recommendation'])
正态性违反:
方差齐性违反:
pg.ttest 在 correction='auto' 时自动应用)pg.welch_anova)或 Brown-Forsythe线性违反(回归):
形式化检验随着 n 增大而变得过于敏感:对于 n ≥ 100,Q-Q 图的权重应高于 Shapiro-Wilk p 值。参见 references/assumptions_and_diagnostics.md 获取全面指导。
主要库:
import pingouin as pg
# correction='auto' 在方差不相等时应用 Welch 校正
result = pg.ttest(group_a, group_b, correction='auto')
# Pingouin >= 0.6 列名
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
print(f"t({df:.0f}) = {t_stat:.2f}, p = {p_value:.3f}, d = {cohens_d:.2f}")
import pingouin as pg
aov = pg.anova(dv='score', between='group', data=df, detailed=True)
print(aov)
# 效应量:偏 eta 平方
eta_p2 = aov['np2'].values[0]
# 如果显著,进行事后检验(Tukey HSD 控制族系误差)
if aov['p_unc'].values[0] < 0.05:
posthoc = pg.pairwise_tukey(dv='score', between='group', data=df)
print(posthoc) # 包含每对的 Hedges' g
import statsmodels.api as sm
from assumption_checks import check_regression_diagnostics
X = sm.add_constant(X_predictors) # 添加截距
model = sm.OLS(y, X).fit()
print(model.summary())
# 4 面板残差图 + Shapiro-Wilk、Breusch-Pagan、Durbin-Watson、VIF
diag = check_regression_diagnostics(model)
print(diag['interpretation'])
print(diag['vif'])
# 如果检测到异方差性,改用稳健标准误
robust = model.get_robustcov_results('HC3')
import pymc as pm
import arviz as az
import numpy as np
with pm.Model() as model:
# 先验
mu1 = pm.Normal('mu_group1', mu=0, sigma=10)
mu2 = pm.Normal('mu_group2', mu=0, sigma=10)
sigma = pm.HalfNormal('sigma', sigma=10)
# 似然
y1 = pm.Normal('y1', mu=mu1, sigma=sigma, observed=group_a)
y2 = pm.Normal('y2', mu=mu2, sigma=sigma, observed=group_b)
# 导出量
diff = pm.Deterministic('difference', mu1 - mu2)
trace = pm.sample(2000, tune=1000)
# ArviZ 1.x 默认为 89% 区间;为报告显式请求 95%
print(az.summary(trace, var_names=['difference'], ci_prob=0.95))
# 直接概率陈述(这就是单侧问题的转化形式)
prob_greater = np.mean(trace.posterior['difference'].values > 0)
print(f"P(mu1 > mu2 | data) = {prob_greater:.3f}")
# ArviZ 1.x 移除了 az.plot_posterior;使用 plot_dist(在 0.x 上,plot_posterior 仍然可用)
az.plot_dist(trace, var_names=['difference'], ci_prob=0.95)
根据数据缩放先验(例如,sigma=10 适用于 SD 接近 10 的结果;以观察到的 SD 为参考),并在报告中说明先验。
效应量量化幅度;p 值仅指示存在性。 为每个检验报告效应量。参见 references/effect_sizes_and_power.md 获取完整指南。
| 检验 | 效应量 | 小 | 中 | 大 |
|---|---|---|---|---|
| T 检验 | Cohen's d | 0.20 | 0.50 | 0.80 |
| 方差分析 | η²_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 检验:检测 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"每组所需样本量: {n_required:.0f}")
# 单因素方差分析:检测 Cohen's f = 0.25 需要多少样本?
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"所需总样本量: {math.ceil(n_total)}(每组 {math.ceil(n_total / 3)})")
确定研究能够检测到的效应量:
# 每组 n=50,在 80% 功效下能检测到什么效应?
detectable_d = tt_ind_solve_power(
effect_size=None, # 求解此项
nobs1=50,
alpha=0.05,
power=0.80,
ratio=1.0,
alternative='two-sided'
)
print(f"研究能检测到 d >= {detectable_d:.2f}")
注意:事后"观察功效"(从观察到的效应计算功效)是循环论证且具有误导性——它是 p 值的确定性函数。如果研究已经完成并且有人询问功效,请改为进行敏感性分析。
参见 references/effect_sizes_and_power.md 获取详细指南。
遵循 references/reporting_standards.md 的 APA 风格。每份报告需要:
p = .034,而非 p < .05;仅在低于 .001 时使用 p < .001)A 组(n = 48, M = 75.2, SD = 8.5)的得分显著高于 B 组(n = 52, M = 68.3, SD = 9.2),t(98) = 3.82, p < .001, d = 0.77, 95% CI [0.36, 1.18],双尾。正态性假设(Shapiro-Wilk: A 组 W = 0.97, p = .18;B 组 W = 0.96, p = .12)和方差 齐性假设(Levene's F(1, 98) = 1.23, p = .27)均满足。
单因素方差分析显示处理条件对测试分数有显著的主效应, F(2, 147) = 8.45, p < .001, η²_p = .10。使用 Tukey's HSD 的 事后比较表明,条件 A(M = 78.2, SD = 7.3)显著高于 条件 B(M = 71.5, SD = 8.1, p = .002, d = 0.87)和条件 C (M = 70.1, SD = 7.9, p < .001, d = 1.07)。条件 B 与 C 之间 无显著差异(p = .52, d = 0.18)。
进行了多元线性回归,以根据学习时间、先前 GPA 和出勤率 预测考试成绩。整体模型显著,F(3, 146) = 45.2, p < .001, R² = .48, 调整后 R² = .47。学习时间(B = 1.80, SE = 0.31, β = .35, t = 5.78, p < .001, 95% CI [1.18, 2.42])和 先前 GPA(B = 8.52, SE = 1.95, β = .28, t = 4.37, p < .001, 95% CI [4.66, 12.38])是显著预测因子,而出勤率则不是 (B = 0.15, SE = 0.12, β = .08, t = 1.25, p = .21, 95% CI [-0.09, 0.39])。 多重共线性不是问题(所有 VIF < 1.5)。
使用弱信息先验(组均值的 Normal(0, 10))进行了贝叶斯 独立样本 t 检验。后验分布表明 A 组得分高于 B 组 (M_diff = 6.8, 95% 可信区间 [3.2, 10.4]),A 组均值 超过 B 组均值的后验概率为 99.8%。收敛诊断令人满意 (所有 R-hat < 1.01, ESS > 1000)。
如果使用了非参数检验,报告中位数而非均值、U/W/H 统计量以及基于秩的效应量(例如,pg.mwu 以 RBC 返回的秩-双列相关)。
在以下情况下考虑贝叶斯方法:
参见 references/bayesian_statistics.md 了解先验设定、贝叶斯因子、可信区间、分层模型和收敛检查(R-hat < 1.01、足够的 ESS、后验预测检查)。
references/)scripts/)comprehensive_assumption_check():一次调用完成异常值 + 正态性 + 方差齐性check_normality()、check_normality_per_group():带 Q-Q 图的 Shapiro-Wilkcheck_homogeneity_of_variance():带箱线图的 Levene 检验check_regression_diagnostics():4 面板残差图 + 已拟合 OLS 模型的 Shapiro-Wilk、Breusch-Pagan、Durbin-Watson、VIFcheck_linearity()、detect_outliers()这些是使分析站得住脚的做法。它们很重要,因为最常见的统计失败不是计算错误——而是静默的灵活性(一直测试直到有结果)和选择性报告。
相关专题
热门下载
相关下载
精品课程
共6课时 | 54.7万人学习
共89课时 | 134.1万人学习
共49课时 | 82.6万人学习
最新文章