R Shiny中向brms传递用户选择系数变量的参数解析问题

问题场景 #

开发R Shiny Web应用时，支持用户上传数据集、选择分析变量，并用自定义先验值拟合贝叶斯ANCOVA模型。点击触发按钮input$activate后，出现如下错误：

Warning: Error in : The following priors do not correspond to any model parameter:
b_input$baseline ~ normal(baseline_mean, baseline_sd)
b_input$group ~ normal(treatment_mean, treatment_sd)
Function 'default_prior' might be helpful to you.

错误核心：brms将coef = input$baseline中的input$baseline当作字面量字符串处理，而非读取它存储的实际系数名（如measurement_1）。

相关代码片段 #

服务端核心逻辑 #

bfit <- eventReactive(
  input$activate,
  ignoreNULL = F,
  {
    req(lmmformula, anovaformula, ancovaformula, input$interceptMean,
        input$interceptVariance, input$baselineMean, input$baselineVariance,
        input$treatmentMean, input$treatmentVariance, input$sigmaScale,
        input$seed)
    
    stanparams <- stanvar(input$interceptMean, "intercept_mean") +
      stanvar(sqrt(input$interceptVariance), "intercept_sd") +
      stanvar(input$baselineMean, "baseline_mean") + 
      stanvar(sqrt(input$baselineVariance), "baseline_sd") + 
      stanvar(input$treatmentMean, "treatment_mean") + 
      stanvar(sqrt(input$treatmentVariance), "treatment_sd") + 
      stanvar(input$sigmaScale, "sigma_scale")
    
    priors <- c(prior(normal(intercept_mean,
                             intercept_sd),
                      class = "Intercept"),
                prior(normal(baseline_mean,
                             baseline_sd),
                      class = "b",
                      coef = input$baseline),
                prior(normal(treatment_mean,
                             treatment_sd),
                      class = "b",
                      coef = input$group),
                prior(cauchy(0, sigma_scale),
                      class = "sigma"))
    
    list(
      bayes = brm(ancovaformula(), data = sel_data(), silent = 1,
                  thin = 1, iter = input$iter, warmup = (input$iter/2),
                  prior = priors, stanvars = stanparams, seed = input$seed),
      
      bayesx2 = brm(ancovaformula(), data = sel_data(), silent = 1,
                    thin = 1, iter = (input$iter)*2, warmup = input$iter,
                    prior = priors, stanvars = stanparams, seed = input$seed)
    )
  }
)

observeEvent(input$activate, {
  bfit()
})

ANCOVA公式定义 #

ancovaformula <- reactive({
    req(sel_data(), longdata())
    ancovaformula = as.formula(paste(input$fwup, " ~ ", input$baseline,
                                     " + ", input$group))
  })

示例：当输入值为post、pre、Tx时，生成公式为post ~ pre + Tx

解决方案 #

要让brms正确解析变量存储的系数名，有两种可行方式：

方式1：先赋值为本地变量 #

将input$baseline和input$group的值先存入本地字符串变量，再传递给coef参数：

# 先缓存变量值
baseline_coef <- input$baseline
group_coef <- input$group

# 重新定义先验
priors <- c(prior(normal(intercept_mean, intercept_sd), class = "Intercept"),
            prior(normal(baseline_mean, baseline_sd), class = "b", coef = baseline_coef),
            prior(normal(treatment_mean, treatment_sd), class = "b", coef = group_coef),
            prior(cauchy(0, sigma_scale), class = "sigma"))

方式2：使用!!强制解析变量 #

利用R的bang-bang操作符，强制brms解析变量的实际值而非字面量：

priors <- c(prior(normal(intercept_mean, intercept_sd), class = "Intercept"),
            prior(normal(baseline_mean, baseline_sd), class = "b", coef = !!input$baseline),
            prior(normal(treatment_mean, treatment_sd), class = "b", coef = !!input$group),
            prior(cauchy(0, sigma_scale), class = "sigma"))

原理：brms的prior()函数默认将coef参数视为字面量表达式，通过上述两种方式，可让函数读取变量存储的实际系数名，匹配模型中的参数。

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