如何用spread_rvars提取brms模型中嵌套随机效应的后验分布？

数据结构 #

'data.frame':   2492 obs. of  5 variables:
 $ ALive        : int  6 9 9 10 10 1 0 0 8 5 ...
 $ Total        : int  10 10 10 10 10 10 10 10 10 10 ...
 $ Cl_act       : num  1804 1523 1188 296 18 ...
 $ Lakename     : chr  "15mile" "15mile" "15mile" "15mile" ...
 $ CloneNo      : chr  "Mile_1" "Mile_1" "Mile_1" "Mile_1" ...

brms模型定义 #

bf(Response | trials(Total) ~ 1/(1+exp(b*(log(cloride_act)-log(LC50)))),
            LC50 ~ 1+(1|Lakename/CloneNo),
            b ~ 1+(1|Lakename/CloneNo),
            family = binomial(link = "identity"),
            nl=T)

问题详情 #

使用get_variable(fit) %>% dput()已获取模型完整变量列表，但调用spread_rvars()提取后验分布时，常规参数可正常提取，嵌套随机效应提取失败：

• 直接使用变量名时，报错NA/NaN argument
• 用反引号包裹变量名时，报错“变量未找到”
  推测是语法匹配问题，无法自行修正，寻求解决方法。

解决方案 #

brms非线性模型的嵌套随机效应变量有特定命名规则，必须严格匹配get_variables(fit)输出的变量名才能正确提取。

操作步骤 #

1. 确认准确变量名
   先执行get_variables(fit)查看变量列表，非线性模型的嵌套随机效应变量名通常类似：

   get_variables(fit)
   # 示例输出："b" "LC50" "r_Lakename[LC50]" "r_Lakename__CloneNo[LC50]" "r_Lakename[b]" "r_Lakename__CloneNo[b]"
   

   其中r_Lakename__CloneNo[LC50]就是CloneNo嵌套在Lakename下的LC50随机效应，r_Lakename[LC50]是Lakename层级的随机效应。

2. 正确调用spread_rvars()
   变量名包含__、[]这类特殊字符，必须用反引号包裹，示例代码：

   library(tidybayes)
   
   # 提取嵌套随机效应及固定效应
   posterior_draws <- fit %>%
     spread_rvars(
       `r_Lakename__CloneNo[LC50]`,
       `r_Lakename__CloneNo[b]`,
       b,
       LC50
     )
   

3. 结合数据计算预测值
   如果需要基于原始数据计算预测概率，可合并数据后按模型公式计算：

   # 合并原始数据与后验分布
   predicted <- posterior_draws %>%
     left_join(your_data, by = c("Lakename", "CloneNo")) %>%
     # 按模型公式计算预测存活概率
     mutate(
       pred_prob = 1 / (1 + exp( (b + `r_Lakename__CloneNo[b]`) * (log(Cl_act) - log(LC50 + `r_Lakename__CloneNo[LC50]`)) ))
     )
   

关键注意事项 #

• 绝对不要自行猜测变量名，必须以get_variables(fit)的输出为准，brms对非线性模型的随机效应命名逻辑和线性模型不同
• 变量名中的所有特殊符号（包括[]、__）都是命名的一部分，必须完整包含在反引号内
• 若仍报错，检查是否存在拼写错误，或是否在spread_rvars()中混入了不存在的变量名

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