skillhub-086-statistical-analysis/references/bayesian_statistics.md

name, description, metadata

name	description	metadata
bayesian-statistical-analysis	贝叶斯统计分析指南——涵盖贝叶斯 vs 频率学派哲学、贝叶斯定理、先验分布、贝叶斯假设检验、贝叶斯估计、常见贝叶斯分析、层次模型、模型比较、模型检查及报告规范	type reference

贝叶斯统计分析

本文档为进行和解读贝叶斯统计分析提供指导，贝叶斯统计为频率学派（经典）统计提供了一种替代框架。

贝叶斯 vs. 频率学派哲学

根本差异

方面	频率学派	贝叶斯学派
概率解释	事件的长期发生频率	信念/不确定性的程度
参数	固定但未知	具有分布的随机变量
推断	基于抽样分布	基于后验分布
主要输出	p 值、置信区间	后验概率、可信区间
先验信息	不正式纳入	通过先验显式纳入
假设检验	拒绝/不拒绝原假设	给定数据下假设的概率
样本量	通常需要最小样本量	可处理任意样本量
解释	间接（给定 H₀ 下数据的概率）	直接（给定数据下假设的概率）

关键问题差异

频率学派：「如果原假设为真，观察到如此极端或更极端数据的概率是多少？」

贝叶斯学派：「给定观测数据，假设为真的概率是多少？」

贝叶斯问题更加直观，直接回答了研究者想要知道的问题。

---

贝叶斯定理

公式：

P(θ|D) = P(D|θ) × P(θ) / P(D)

用文字表述：

后验 = 似然 × 先验 / 证据

其中：

• θ（theta）：感兴趣的参数（例如均值差、相关系数）
• D：观测数据
• P(θ|D): 后验分布（看到数据后对 θ 的信念）
• P(D|θ): 似然（给定 θ 下数据的概率）
• P(θ): 先验分布（看到数据前对 θ 的信念）
• P(D): 边际似然/证据（归一化常数）

---

先验分布

先验的类型

1. 信息性先验

何时使用：当你拥有大量先验知识时，来源包括：

• 先前的研究
• 专家知识
• 理论
• 预试验数据

示例：元分析显示效应量 d ≈ 0.40，SD = 0.15

• 先验：Normal(0.40, 0.15)

优点：

• 纳入已有知识
• 更高效（所需样本量更小）
• 可在小样本情况下稳定估计

缺点：

• 具有主观性（但主观性也可成为优势）
• 必须被合理证明并保持透明
• 若强先验与数据冲突，可能引发争议

---

2. 弱信息性先验

何时使用：大多数应用场景的默认选择

特征：

• 正则化估计（防止极端值）
• 在中等样本量下对后验影响极小
• 防止计算问题

示例先验：

• 效应量：Normal(0, 1) 或 Cauchy(0, 0.707)
• 方差：Half-Cauchy(0, 1)
• 相关系数：Uniform(-1, 1) 或 Beta(2, 2)

优点：

• 在客观性与正则化之间取得平衡
• 计算稳定
• 广泛可接受

---

3. 无信息先验（平坦/均匀先验）

何时使用：试图保持「客观」时

示例：Uniform(-∞, ∞) 对任意值

⚠️ 注意：

• 可能导致非正常后验
• 可能产生不合理的结果
• 并非真正的「无信息」（仍然做了假设）
• 在现代贝叶斯实践中通常不推荐

更好的替代方案：使用弱信息性先验

---

先验敏感性分析

始终进行：检验结果如何随不同先验变化

流程：

1. 用默认/计划先验拟合模型
2. 用更分散的先验拟合模型
3. 用更集中的先验拟合模型
4. 比较后验分布

报告：

• 若结果相似：证据稳健
• 若结果差异显著：数据不足以压倒先验

Python 示例：

import pymc as pm

prior_specs = [
    ('weakly_informative', 0, 1),
    ('diffuse', 0, 10),
    ('informative', 0.5, 0.3),
]

results = {}
for name, mu_prior, sigma_prior in prior_specs:
    with pm.Model() as model:
        effect = pm.Normal('effect', mu=mu_prior, sigma=sigma_prior)
        # ... 似然函数和观测数据
        trace = pm.sample(2000, tune=1000, return_inferencedata=True)
        results[name] = trace

---

贝叶斯假设检验

贝叶斯因子（BF）

定义：两个竞争假设的证据比率

公式：

BF₁₀ = P(D|H₁) / P(D|H₀)

解释：

BF₁₀	证据强度
>100	支持 H₁ 的决定性证据
30-100	支持 H₁ 的极强证据
10-30	支持 H₁ 的强证据
3-10	支持 H₁ 的中等证据
1-3	支持 H₁ 的微弱证据
1	无证据
1/3-1	支持 H₀ 的微弱证据
1/10-1/3	支持 H₀ 的中等证据
1/30-1/10	支持 H₀ 的强证据
1/100-1/30	支持 H₀ 的极强证据
<1/100	支持 H₀ 的决定性证据

相对于 p 值的优势：

1. 可提供支持原假设的证据
2. 不依赖于抽样意图（无「偷看」问题）
3. 直接量化证据
4. 可随更多数据更新

Python 计算：

# Pingouin 0.5+：用于独立双边 t 检验的 BF10；单边 BF 已被移除。
import pingouin as pg

result = pg.ttest(group1, group2, correction=False)
bf10 = result['BF10'].values[0]

# 严格贝叶斯因子：BayesFactor（R）、JASP 或 PyMC 模型比较（参见 pymc 技能）

---

实际等价区间（ROPE）

目的：定义可忽略效应量的范围

流程：

1. 定义 ROPE（例如 d ∈ [-0.1, 0.1] 为可忽略效应）
2. 计算后验落在 ROPE 内的百分比
3. 做出判断：

   • 95% 在 ROPE 内：接受实际等价

   • 95% 在 ROPE 外：拒绝等价

   • 其他情况：无定论

优势：直接检验实际显著性

Python 示例：

# 定义 ROPE
rope_lower, rope_upper = -0.1, 0.1

# 计算后验在 ROPE 内的百分比
in_rope = np.mean((posterior_samples > rope_lower) &
                  (posterior_samples < rope_upper))

print(f"{in_rope*100:.1f}% 的后验在 ROPE 内")

---

贝叶斯估计

可信区间

定义：以 X% 的概率包含参数的区间

95% 可信区间的解释：

「真实参数有 95% 的概率落在此区间内。」

这正是人们以为置信区间所表示的含义（但在频率学派框架中并非如此）

类型：

等尾区间（ETI）

• 第 2.5 至第 97.5 百分位
• 计算简单
• 对于偏态分布可能不包含众数

最高密度区间（HDI）

• 包含分布 95% 的最窄区间
• 始终包含众数
• 更适合偏态分布

Python 计算：

import arviz as az

# 等尾区间
eti = np.percentile(posterior_samples, [2.5, 97.5])

# HDI
hdi = az.hdi(posterior_samples, hdi_prob=0.95)

---

后验分布

解读后验分布：

1. 集中趋势：

   • 均值：后验的平均值
   • 中位数：第 50 百分位
   • 众数：最可能的值（MAP——最大后验估计）

2. 不确定性：

   • 标准差：后验的离散程度
   • 可信区间：量化不确定性

3. 形状：

   • 对称：近似正态
   • 偏态：非对称不确定性
   • 多峰：存在多个合理值

可视化：

import matplotlib.pyplot as plt
import arviz as az

# 带 HDI 的后验图
az.plot_posterior(trace, hdi_prob=0.95)

# 迹线图（检查收敛性）
az.plot_trace(trace)

# 森林图（多个参数）
az.plot_forest(trace)

---

常见贝叶斯分析

贝叶斯 T 检验

目的：比较两组（t 检验的贝叶斯替代方法）

输出：

1. 均值差的后验分布
2. 95% 可信区间
3. 贝叶斯因子（BF₁₀）
4. 方向性假设的概率（例如 P(μ₁ > μ₂)）

Python 实现：

import pymc as pm
import arviz as az

# 贝叶斯独立样本 t 检验
with pm.Model() as model:
    # 组均值的先验
    mu1 = pm.Normal('mu1', mu=0, sigma=10)
    mu2 = pm.Normal('mu2', mu=0, sigma=10)

    # 合并标准差的前沿
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 似然函数
    y1 = pm.Normal('y1', mu=mu1, sigma=sigma, observed=group1)
    y2 = pm.Normal('y2', mu=mu2, sigma=sigma, observed=group2)

    # 衍生量：均值差
    diff = pm.Deterministic('diff', mu1 - mu2)

    # 后验抽样
    trace = pm.sample(2000, tune=1000, return_inferencedata=True)

# 分析结果
print(az.summary(trace, var_names=['mu1', 'mu2', 'diff']))

# group1 > group2 的概率
prob_greater = np.mean(trace.posterior['diff'].values > 0)
print(f"P(μ₁ > μ₂) = {prob_greater:.3f}")

# 绘制后验图
az.plot_posterior(trace, var_names=['diff'], ref_val=0)

---

贝叶斯方差分析（ANOVA）

目的：比较三个或更多组

模型：

import pymc as pm

with pm.Model() as anova_model:
    # 超先验
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_between = pm.HalfNormal('sigma_between', sigma=5)
    sigma_within = pm.HalfNormal('sigma_within', sigma=5)

    # 组均值（层次结构）
    group_means = pm.Normal('group_means',
                            mu=mu_global,
                            sigma=sigma_between,
                            shape=n_groups)

    # 似然函数
    y = pm.Normal('y',
                  mu=group_means[group_idx],
                  sigma=sigma_within,
                  observed=data)

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

# 后验对比
contrast_1_2 = trace.posterior['group_means'][:,:,0] - trace.posterior['group_means'][:,:,1]

---

贝叶斯相关分析

目的：估计两个变量之间的相关系数

优势：提供相关系数值的分布

Python 实现：

import pymc as pm

with pm.Model() as corr_model:
    # 相关系数的先验
    rho = pm.Uniform('rho', lower=-1, upper=1)

    # 转换为协方差矩阵
    cov_matrix = pm.math.stack([[1, rho],
                                [rho, 1]])

    # 似然函数（二元正态分布）
    obs = pm.MvNormal('obs',
                     mu=[0, 0],
                     cov=cov_matrix,
                     observed=np.column_stack([x, y]))

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

# 相关系数汇总
print(az.summary(trace, var_names=['rho']))

# 相关系数为正的概率
prob_positive = np.mean(trace.posterior['rho'].values > 0)

---

贝叶斯线性回归

目的：建模预测变量与结果变量之间的关系

优势：

• 所有参数均带有不确定性
• 自然正则化（通过先验）
• 可纳入先验知识
• 预测的可信区间

Python 实现：

import pymc as pm

with pm.Model() as regression_model:
    # 系数的先验
    alpha = pm.Normal('alpha', mu=0, sigma=10)  # 截距
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 期望值
    mu = alpha + pm.math.dot(X, beta)

    # 似然函数
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

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

# 后验预测检验
with regression_model:
    ppc = pm.sample_posterior_predictive(trace)

az.plot_ppc(ppc)

# 带不确定性的预测
with regression_model:
    pm.set_data({'X': X_new})
    posterior_pred = pm.sample_posterior_predictive(trace)

---

层次（多水平）模型

何时使用：

• 嵌套/聚类数据（学生嵌套于学校）
• 重复测量
• 元分析
• 跨组的变异性效应

关键概念：部分池化

• 完全池化：忽略组别（有偏）
• 无池化：分别分析各组（高方差）
• 部分池化：跨组借用信息（贝叶斯）

示例：变截距模型：

with pm.Model() as hierarchical_model:
    # 超先验
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_between = pm.HalfNormal('sigma_between', sigma=5)
    sigma_within = pm.HalfNormal('sigma_within', sigma=5)

    # 组级截距
    alpha = pm.Normal('alpha',
                     mu=mu_global,
                     sigma=sigma_between,
                     shape=n_groups)

    # 似然函数
    y_obs = pm.Normal('y_obs',
                     mu=alpha[group_idx],
                     sigma=sigma_within,
                     observed=y)

    trace = pm.sample()

---

模型比较

方法

1. 贝叶斯因子

• 直接比较模型证据
• 对先验设定敏感
• 计算量可能较大

2. 信息准则

WAIC（广泛适用信息准则）：

• AIC 的贝叶斯类比
• 越小越好
• 考虑了参数的有效数量

LOO（留一法交叉验证）：

• 估计样本外预测误差
• 越小越好
• 比 WAIC 更稳健

Python 计算：

import arviz as az

# 计算 WAIC 和 LOO
waic = az.waic(trace)
loo = az.loo(trace)

print(f"WAIC: {waic.elpd_waic:.2f}")
print(f"LOO: {loo.elpd_loo:.2f}")

# 比较多个模型
comparison = az.compare({
    'model1': trace1,
    'model2': trace2,
    'model3': trace3
})
print(comparison)

---

检查贝叶斯模型

1. 收敛诊断

R-hat（Gelman-Rubin 统计量）：

• 比较链内方差与链间方差
• 接近 1.0 的值表明收敛
• R-hat < 1.01：良好
• R-hat > 1.05：收敛不佳

有效样本量（ESS）：

• 独立样本的数量
• 越高越好
• 建议每链 ESS > 400

迹线图：

• 应呈现「毛虫状」
• 无趋势，无卡顿链

Python 检查：

# 带诊断的自动汇总
print(az.summary(trace, var_names=['parameter']))

# 可视化诊断
az.plot_trace(trace)
az.plot_rank(trace)  # 秩图

---

2. 后验预测检验

目的：模型生成的数据是否与观测数据相似？

流程：

1. 从后验生成预测
2. 与实际数据比较
3. 寻找系统性偏差

Python 实现：

with model:
    ppc = pm.sample_posterior_predictive(trace)

# 可视化检验
az.plot_ppc(ppc, num_pp_samples=100)

# 定量检验
obs_mean = np.mean(observed_data)
pred_means = [np.mean(sample) for sample in ppc.posterior_predictive['y_obs']]
p_value = np.mean(pred_means >= obs_mean)  # 贝叶斯 p 值

---

报告贝叶斯结果

T 检验报告示例

「采用贝叶斯独立样本 t 检验比较 A 组和 B 组。使用了弱信息性先验：均值差采用 Normal(0, 1)，合并标准差采用 Half-Cauchy(0, 1)。均值差的后验分布均值为 5.2（95% CI [2.3, 8.1]），表明 A 组得分高于 B 组。贝叶斯因子 BF₁₀ = 23.5 为组间差异提供了强证据，A 组均值超过 B 组均值的概率为 99.7%。」

回归报告示例

「采用弱信息性先验（系数使用 Normal(0, 10)，残差标准差使用 Half-Cauchy(0, 5)）拟合了贝叶斯线性回归模型。模型解释了可观的方差（R² = 0.47，95% CI [0.38, 0.55]）。学习时长（β = 0.52，95% CI [0.38, 0.66]）和先前的 GPA（β = 0.31，95% CI [0.17, 0.45]）是可信的预测变量（95% CI 排除零）。后验预测检验显示模型拟合良好。收敛诊断令人满意（所有 R-hat < 1.01，ESS > 1000）。」

---

优势与局限

优势

1. 直观的解释：关于参数的直接概率陈述
2. 纳入先验知识：利用所有可用信息
3. 灵活：轻松处理复杂模型
4. 无 p 值操纵：可随时查看数据
5. 量化不确定性：完整的后验分布
6. 小样本：可处理任意样本量

局限

1. 计算量大：需要 MCMC 抽样（可能较慢）
2. 先验设定：需要思考和论证
3. 复杂性：学习曲线较陡
4. 软件工具：相比频率学派方法工具较少
5. 沟通成本：可能需要向审稿人/读者做解释

---

关键 Python 包

使用 uv 安装（参见 SKILL.md）。ArviZ 0.23+ 需要 Python 3.12+。

• PyMC (pymc>=5)：完整的贝叶斯建模框架
• ArviZ (arviz>=0.17)：可视化与诊断工具（文档）
• Bambi：回归模型的高级接口（uv pip install bambi）
• PyStan：Stan 的 Python 接口
• TensorFlow Probability：使用 TensorFlow 进行贝叶斯推断

---

何时使用贝叶斯方法

使用贝叶斯方法当：

• 你有先验信息需要纳入
• 你想要直接的概率陈述
• 样本量较小
• 模型较复杂（层次模型、缺失数据等）
• 你希望随着数据到来不断更新分析

频率学派方法可能足够当：

• 标准分析且样本量大
• 无先验信息
• 计算资源有限
• 审稿人不熟悉贝叶斯方法