使用brms绘制二项式模型预测后验分布遇报错求助

算力阿峰

2026-7-5

关于brms绘制多变量二项式模型预测图的问题

有人用过brms绘制多变量二项式模型的预测结果吗？我只有单响应变量的代码，运行绘图代码后出现了如下警告：

Warning message:
Computation failed in stat_eye()
Caused by error in bw.SJ():
! sample is too sparse to find TD

我的代码如下：

My_model <- brm(FutureReproduction~ Reproductive_success+ Age 
     + Precipitation+ Scaled_Leg_index 
     + Reproductive_success*Age
     + (1 | ID),
     family = bernoulli()) 

condition <-expand.grid("Age" = (c(3, 4, 5, 6, 7,8, 9, 10, 11)),
                        Year = 2011,
                        ID = 3,
                        Precipitation= 0,
                        Leg_index= 0,
                        Reproductive_success = c(0, 1))


Model_prediction <- condition %>% 
  add_predicted_rvars(My_model, allow_new_levels = TRUE, re_formula = NA, newdata = .) %>% 
  mutate(Reproductive_success = factor(ifelse(Reproductive_success == 1, "Yes", "No"), levels = c("No", "Yes")))

Model_prediction的输出：

# A tibble: 18 × 7
     Age  Year    ID       Precipitation      Leg_index  Reproductive_success  prediction[,"Jumpheight"]  [,"Sprintspeed"] [,"Future Reproduction"] [,"Survival"]
   <dbl> <dbl> <dbl>              <dbl>            <dbl> <fct>                        <rvar[,1]>    <rvar[,1]>          <rvar[,1]>      <rvar[,1]>
 1     3  2011     3                  0                0 No                          21.33 ± 5.9    2.41 ± 1.8        0.261 ± 0.44     0.97 ± 0.18
 2     4  2011     3                  0                0 No                          10.63 ± 6.0    1.57 ± 1.8        0.230 ± 0.42     0.97 ± 0.17
 3     5  2011     3                  0                0 No                           4.60 ± 5.9    1.03 ± 1.8        0.202 ± 0.40     0.97 ± 0.18
 4     6  2011     3                  0                0 No                           2.49 ± 5.9    0.74 ± 1.8        0.188 ± 0.39     0.96 ± 0.19
 5     7  2011     3                  0                0 No                           2.63 ± 5.9    0.67 ± 1.8        0.168 ± 0.37     0.96 ± 0.19
 6     8  2011     3                  0                0 No                           3.59 ± 5.9    0.64 ± 1.8        0.153 ± 0.36     0.95 ± 0.22
 7     9  2011     3                  0                0 No                           5.12 ± 5.9    0.66 ± 1.8        0.148 ± 0.36     0.93 ± 0.25
 8    10  2011     3                  0                0 No                           6.88 ± 6.0    0.69 ± 1.8        0.136 ± 0.34     0.90 ± 0.29
 9    11  2011     3                  0                0 No                           8.97 ± 6.4    0.79 ± 1.9        0.134 ± 0.34     0.85 ± 0.35
10     3  2011     3                  0                0 Yes                          9.54 ± 6.6   -0.66 ± 2.0        0.029 ± 0.17     0.91 ± 0.28
11     4  2011     3                  0                0 Yes                          5.97 ± 5.9   -0.41 ± 1.8        0.025 ± 0.16     0.93 ± 0.25
12     5  2011     3                  0                0 Yes                          2.89 ± 6.1   -0.71 ± 1.8        0.022 ± 0.15     0.94 ± 0.23
13     6  2011     3                  0                0 Yes                          0.46 ± 6.0   -1.37 ± 1.8        0.018 ± 0.13     0.95 ± 0.23
14     7  2011     3                  0                0 Yes                         -1.18 ± 6.0   -2.11 ± 1.8        0.016 ± 0.13     0.95 ± 0.23
15     8  2011     3                  0                0 Yes                         -2.21 ± 6.1   -2.56 ± 1.8        0.017 ± 0.13     0.94 ± 0.24
16     9  2011     3                  0                0 Yes                         -2.53 ± 6.1   -2.63 ± 1.8        0.018 ± 0.13     0.92 ± 0.27
17    10  2011     3                  0                0 Yes                         -2.30 ± 6.2   -2.54 ± 1.8        0.017 ± 0.13     0.89 ± 0.32
18    11  2011     3                  0                0 Yes                         -1.54 ± 6.4   -2.22 ± 1.9        0.019 ± 0.14     0.86 ± 0.35

绘图代码：

Model_prediction %>% 
  ggplot(aes(x = Reproductive_success , ydist = .prediction[,"Future Reproduction"]) )+
  stat_eye()

运行后出现上述警告，但用同样方法绘制高斯分布（连续响应变量）的图时完全正常。

解决方法

这个警告的核心原因是二项式模型的预测值（概率）过于集中，导致stat_eye()使用的SJ带宽估计方法无法计算有效带宽——SJ方法需要数据有足够的分散度，而你的部分预测概率（比如Reproductive_success=Yes组）几乎都集中在0附近，数据稀疏度太高。

可以用以下几种方法解决：

手动指定带宽：给stat_eye()加上bw参数，手动设置一个合适的带宽值，比如：
```
Model_prediction %>% 
  ggplot(aes(x = Reproductive_success , ydist = .prediction[,"Future Reproduction"]) )+
  stat_eye(bw = 0.05)
```
可以根据数据分布调整bw的数值，比如0.02、0.1等，直到图形正常显示。

改用更适合离散/集中数据的可视化方式：二项式模型的预测是概率分布，用stat_halfeye()或者直接展示中位数和置信区间可能更合适，比如：

# 用stat_halfeye，默认带宽方法更鲁棒
Model_prediction %>% 
  ggplot(aes(x = Reproductive_success , ydist = .prediction[,"Future Reproduction"]) )+
  stat_halfeye()

# 直接展示中位数和95%置信区间
Model_prediction %>% 
  mutate(
    median_pred = median(.prediction[,"Future Reproduction"]),
    lower_pred = quantile(.prediction[,"Future Reproduction"], 0.025),
    upper_pred = quantile(.prediction[,"Future Reproduction"], 0.975)
  ) %>% 
  ggplot(aes(x = Reproductive_success, y = median_pred)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower_pred, ymax = upper_pred), width = 0.2)

调整预测分布的采样：如果模型的预测采样数太少，也可能导致数据稀疏，可以在add_predicted_rvars()中增加ndraws参数，让采样更密集：

Model_prediction <- condition %>% 
  add_predicted_rvars(My_model, allow_new_levels = TRUE, re_formula = NA, newdata = ., ndraws = 1000) %>% 
  mutate(Reproductive_success = factor(ifelse(Reproductive_success == 1, "Yes", "No"), levels = c("No", "Yes")))

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