如何将列标题插入到贝叶斯回归模型 (brm) 的 for 循环中？

Question

我正在尝试在成对的列上运行 for 循环，以实现零一膨胀的 beta 回归。我已经对其他模型（例如 glm 或 betareg）运行了类似的循环，但在这种情况下它无法识别列标题。

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' 

如果只有一组列，我会运行类似下面的命令来生成一个列表，但我不知道如何在 ? 中引入第二组列？位置。

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
    )))
  }))

我只想运行 a1 和 b1、a2 和 b2 等对，因此运行一个 double 来对每个 a 执行所有 bs 会非常低效。

任何循环无法运行的想法或其他方法的建议都很棒。

Answer 1

一般规则是您应该使用

as.formula()

或

formula()

，或者在本例中使用

reformulate()

（这并没有什么不同，但我发现更美观）来构建您的公式。

首先是一个辅助函数，因为我们将重复相同的过程，但有细微的变化：

ffun <- function(i, response) {
    reformulate(c("x", paste0("b", i)), response = response)
}

在您的上下文中使用它：

library(brms)
for (i in 1:3) {
    model <- brm(
        formula = bf(
            reformulate("x", response = sprintf("y | weights(a%d)", i)),
            ffun(i, "phi"),
            ffun(i, "zoi"),
            ffun(i, "coi"),
            family = zero_one_inflated_beta()
        ),
        data = dat
    )
}