如何在R的betaBayes包中为贝叶斯Beta回归添加马蹄形收缩先验？

betaBayes确实支持通过自定义先验实现马蹄形收缩，核心是利用包的custom先验模式，结合Stan代码片段定义马蹄形先验的结构。以下是具体实现步骤和示例代码：

1. 准备数据与加载包 #

首先加载包并模拟示例数据（可替换为你的真实数据）：

library(betaBayes)
set.seed(123)

# 模拟数据：前2个特征有真实效应，后3个无效应
n <- 100
p <- 5
X <- matrix(rnorm(n*p), ncol = p)
mu <- plogis(X %*% c(1, -0.5, 0, 0, 0))
phi <- 10
y <- rbeta(n, shape1 = mu*phi, shape2 = (1-mu)*phi)
data_df <- data.frame(y, X)

2. 定义马蹄形先验的Stan代码片段 #

马蹄形先验的核心结构是：回归系数$\beta_j \sim \mathcal{N}(0, \lambda_j^2 \tau^2)$，其中$\lambda_j \sim \text{Half-Cauchy}(0,1)$，$\tau \sim \text{Half-Cauchy}(0,1)$。将这个结构写成Stan代码片段：

horseshoe_prior <- "
// 定义马蹄形先验的超参数
real<lower=0> tau;
vector[p] lambda;

// 超参数先验（半柯西分布）
tau ~ cauchy(0, 1);
lambda ~ cauchy(0, 1);

// 回归系数的马蹄形先验
beta ~ normal(0, lambda * tau);
"

注：Stan中通过real<lower=0>约束变量为正，从而实现半柯西分布。

3. 运行带马蹄形先验的Beta回归 #

调用beta_reg()函数，指定prior = "custom"并传入自定义先验代码：

hs_model <- beta_reg(
  formula = y ~ .,
  data = data_df,
  prior = "custom",
  prior_code = horseshoe_prior,
  phi_prior = "gamma(1, 1)",  # 精度参数phi的先验，可按需调整
  chains = 4,
  iter = 2000,
  warmup = 1000,
  refresh = 0  # 关闭迭代过程输出
)

4. 基于后验结果做特征选择 #

查看系数后验统计 #

用summary()查看系数的后验均值、置信区间，若95%置信区间包含0，说明该特征的效应不显著：

summary(hs_model, pars = "beta")

计算后验包含概率（PIP） #

PIP是系数后验绝对值大于阈值的概率，值越高说明特征越重要：

# 提取系数后验样本
beta_samples <- rstan::extract(hs_model$stanfit, pars = "beta")$beta

# 计算每个特征的PIP（阈值设为0.1，可调整）
pip_vals <- apply(beta_samples, 2, function(x) mean(abs(x) > 0.1))
names(pip_vals) <- colnames(X)
print(pip_vals)

注意事项 #

• 超参数尺度可调整：如果数据特征尺度差异大，可修改半柯西分布的尺度参数（比如把cauchy(0,1)改成cauchy(0,2.5)）。
• 收敛检查：务必检查模型的Rhat值（应接近1）和有效样本量（ESS），确保后验样本可靠。
• 多场景适配：如果用beta_binomial_reg()处理计数型Beta数据，只需调整公式并保持先验代码逻辑一致即可。

内容的提问来源于stack exchange，提问作者te time