如何在brm零一膨胀Beta回归循环中正确引用列名？

问题背景 #

针对成对列（a1&b1、a2&b2、a3&b3）运行零一膨胀Beta回归循环，此前在glm或betareg等模型中类似循环正常，但brm模型无法识别列名，仅需成对组合，无需全组合。

测试数据 #

dat <- data.frame(y =  rbeta(10, 0.7, 1.5),
                  x = sample(c(0,1), 10, replace=TRUE),
                  a1 = abs(rnorm(10)),
                  a2 = abs(rnorm(10)),
                  a3 = abs(rnorm(10)),
                  b1 = abs(rnorm(10)),
                  b2 = abs(rnorm(10)),
                  b3 = abs(rnorm(10)))

尝试的循环代码 #

library(brms)

for (i in 1:3){
model <- brm(
  formula = bf(
    y | weights(paste0("a",i)) ~ x,
    phi ~ x + paste0("b",i),
    zoi ~ x + paste0("b",i),
    coi ~ x + paste0("b",i), 
    family = zero_one_inflated_beta()
  ),
  data = dat
)
mod <- summary(model)
e[i] <- mod$fixed["x1", "Estimate"]
}

报错信息 #

Error: The following variables can neither be found in 'data' nor in 'data2':'i'

尝试的reframe方法（未完成） #

data %>% reframe(across(a1:a3,\(a) {
    list(summary(brm(
      formula = bf(
        y | weights(a) ~ x,
        phi ~ x + ?,
        zoi ~ x + ?,
        coi ~ x + ?, 
        family = zero_one_inflated_beta()
      ),
      data = dat
    )))
  }))

报错原因及解决方法 #

报错原因 #

brms的bf()公式接口不会自动解析字符串拼接的变量名，而是直接把paste0("a",i)当成变量名去数据框中查找，而数据框里没有名为paste0("a",i)或i的列，因此触发报错。

解决方法1：动态构建公式字符串并转换为brmsformula #

通过字符串拼接构建完整公式，再转换为brms可识别的格式：

library(brms)

# 初始化结果向量
e <- numeric(3)

for (i in 1:3) {
  # 构建各部分公式字符串
  resp_form <- paste0("y | weights(a", i, ") ~ x")
  phi_form <- paste0("phi ~ x + b", i)
  zoi_form <- paste0("zoi ~ x + b", i)
  coi_form <- paste0("coi ~ x + b", i)
  
  # 合并为完整的brms公式
  full_formula <- bf(
    as.formula(resp_form),
    as.formula(phi_form),
    as.formula(zoi_form),
    as.formula(coi_form),
    family = zero_one_inflated_beta()
  )
  
  # 拟合模型
  model <- brm(full_formula, data = dat)
  
  # 提取结果
  mod_summary <- summary(model)
  e[i] <- mod_summary$fixed["x1", "Estimate"]
}

解决方法2：使用get()函数引用变量 #

在bf()中通过get()直接提取数据框对应列，需指定环境确保变量可被找到：

library(brms)

e <- numeric(3)

for (i in 1:3) {
  model <- brm(
    formula = bf(
      y | weights(get(paste0("a", i))) ~ x,
      phi ~ x + get(paste0("b", i)),
      zoi ~ x + get(paste0("b", i)),
      coi ~ x + get(paste0("b", i)), 
      family = zero_one_inflated_beta()
    ),
    data = dat,
    envir = environment()
  )
  mod_summary <- summary(model)
  e[i] <- mod_summary$fixed["x1", "Estimate"]
}

解决方法3：重构数据为长格式（更简洁） #

将成对的a、b列转换为长格式，按分组拟合模型，避免动态列名问题：

library(brms)
library(dplyr)
library(tidyr)

# 转换为长格式
dat_long <- dat %>%
  pivot_longer(cols = starts_with(c("a", "b")),
               names_to = c(".value", "group"),
               names_pattern = "(a|b)(\\d)") %>%
  mutate(group = as.integer(group))

# 按group分组拟合模型
models <- dat_long %>%
  group_by(group) %>%
  summarise(
    model = list(brm(
      bf(
        y | weights(a) ~ x,
        phi ~ x + b,
        zoi ~ x + b,
        coi ~ x + b,
        family = zero_one_inflated_beta()
      ),
      data = cur_data()
    )),
    estimate = list(summary(model[[1]])$fixed["x1", "Estimate"])
  )

# 提取结果
e <- pull(models, estimate) %>% unlist()

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