brms非线性模型中预测变量N遇公式错误，求解决方法

原模型与需求 #

我用brms拟合了如下非线性模型，y、N、X、A和B均为实测连续变量，模型拟合效果良好：

mod <- bf(y~(I*N*X^b)/(A+k*B),
          I ~ 1,
          b ~ 1,
          k~1+(1|group),
          nl = T,
          family = gaussian())

现在需要利用已有的y测量值和模型估计的参数，反向预测N。尝试用mi()函数修改模型时出现错误，具体如下：

尝试的错误方法 #

修改后的模型代码：

mod <- bf(y~(I*N*X^b)/(A+k*B),
          y~(I*mi(N)*X^b)/(A+k*B),
          I ~ 1,
          b ~ 1,
          k~1+(1|group),
          nl = T,
          family = gaussian())

运行get_prior(mod)或brm()时触发错误：

Error in terms.formula(formula, ...) : 
incorrect power in the formula

数据示例 #

structure(list(y = c(592200, 551335.652173913, 1720408.69551196, 
5135100, 3710068.69230388, 1904819.95681897), N = c(145L, 41L, 
72L, 3173L, 2966L, 1262L), X = c(404.822115384615, 398.5, 514.76, 
184.786096256684, 184.460601961447, 245.710784313725), A = c(10662371.6311457, 
1924044.03258699, 8519963.12725198, 44606197.3266835, 7148806.05247308, 
40475049.9899619), B = c(107809.545032839, 107809.545032839, 
107809.545032839, 319346.983077122, 319346.983077122, 319346.983077122
), group = c("A", "A", "A", "B", "B", "B")), row.names = c(NA, 
6L), class = "data.frame")

可行解决方案 #

方法1：代数变形后直接计算后验预测 #

你的模型公式可以直接代数解出N的表达式：
[ N = \frac{y \times (A + k \times B)}{I \times X^b} ]
步骤如下：

1. 从拟合好的mod中提取参数后验样本：

   post_samples <- posterior_samples(mod)
   

2. 结合原始数据，对每个后验样本计算N的预测值：

   library(dplyr)
   # 注意：参数名需与posterior_samples输出的列名匹配，brms非线性模型参数默认前缀为b_，随机效应列名如r_group[group,k]
   pred_N <- bind_cols(your_data, post_samples) %>%
     rowwise() %>%
     mutate(N_pred = (y * (A + b_k + r_group[group, k])) / (b_I * X^b_b)) %>%
     ungroup()
   

3. 计算每个观测的N预测统计量（均值、95%置信区间）：

   N_summary <- pred_N %>%
     group_by(y, X, A, B, group) %>%
     summarise(
       N_mean = mean(N_pred),
       N_lower = quantile(N_pred, 0.025),
       N_upper = quantile(N_pred, 0.975),
       .groups = "drop"
     )
   

方法2：重新构建以N为响应变量的模型 #

如果不想手动代数求解，可以直接把N设为响应变量，y作为预测变量，复用原模型的参数估计作为先验加速拟合：

# 构建逆预测模型
mod_N <- bf(N ~ (y * (A + k*B))/(I * X^b),
            I ~ 1,
            b ~ 1,
            k ~ 1 + (1|group),
            nl = TRUE,
            family = gaussian())

# 提取原模型参数作为先验
prior_N <- prior(normal(b_I, 0.1), class = b, coef = I) +
  prior(normal(b_b, 0.1), class = b, coef = b) +
  prior(normal(b_k, 0.1), class = b, coef = k) +
  prior(normal(0, 1), class = sd, coef = k, group = group)

# 拟合模型
fit_N <- brm(mod_N, data = your_data, prior = prior_N, chains = 4, cores = 4)

为什么mi()方法失效？ #

mi()函数主要用于缺失数据的多重插补，在非线性模型公式中直接嵌套会破坏语法解析逻辑（比如X^b的解析出现冲突），并不适用于这种逆预测场景。

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