R中rvar数据叠加原始点到stat_halfeye图及logistic回归绘图问题

问题概述 #

1. 用brms构建Grasshopper_model得到按年龄划分的rvar格式预测分布，希望在stat_halfeye绘制的分布上叠加原始数据点，尝试将数值转为sd=0的rvar失败。
2. 绘制logistic回归预测时出现错误：Computation failed in stat_slabinterval()，根源是bw.SJ()中样本过于稀疏无法找到TD。

解决方案及修改代码 #

问题1：叠加原始数据点 #

不需要将原始数值转为rvar，直接用geom_point即可，只需确保数据的分组（Age、Hopped）与绘图的分面、x轴映射对齐。

问题2：解决稀疏样本带宽错误 #

当样本稀疏时，SJ带宽估计器会失效，可通过手动指定带宽（如bw="nrd0"）、调整带宽缩放系数（adjust=2），或改用直方图作为密度估计方式来解决。

以下是修改后的完整代码：

num_ids <- 10
df <- expand.grid(
    "Age" = c(3, 4, 5, 6, 7, 8, 9, 10, 11),
    "Year" = 2011,
    "ID" = 1:num_ids,  
    "Hopped" = c(0, 1)
)
condition <- expand.grid(
    "Age" = c(3, 4, 5, 6, 7,8, 9, 10, 11),
    Year = 2011,
    ID = 1,
    Hopped = c(0, 1)
) 
df <- df %>%
    group_by(ID) %>% 
    mutate(kJs = ifelse(Hopped == 0, rnorm(1, mean = 3, sd = 2), rnorm(1, mean = 15, sd = 2)))

Filtered_data <- df%>% 
    filter(!is.na(kJs),
           !is.na(Age),
           !is.na(Hopped),
           Age>2,
           Age < 12) %>% 
    mutate(Hopped= ifelse(Hopped== 1, "Yes", "No"))

Grasshopper_model <- brm(kJs ~ Hopped + Age + Hopped*Age + (1|ID),
                         family = gaussian(),
                         data = df)

# 自定义调色板（示例用默认色调板，可替换为自己的配色）
custom_palette <- scales::hue_pal()(length(unique(df$Age)))
strip_colours <- ggplot2::strip_themed()

condition %>% 
    add_predicted_rvars(Grasshopper_model, allow_new_levels = TRUE, re_formula = NA, newdata = .) %>% 
    mutate(Hopped = factor(ifelse(Hopped == 1, "Yes", "No"), levels = c("No", "Yes"))) %>% 
    ggplot(aes(ydist = .prediction[,"kJs"], x = Hopped)) +
    # 指定稳健带宽估计解决稀疏样本问题
    stat_halfeye(aes(fill = factor(Age)), bw = "nrd0") +
    # 叠加原始数据点，自动匹配分面的Age分组
    geom_point(data = Filtered_data , 
               aes(x = Hopped, y = kJs), alpha = 0.1) +
    facet_wrap2(~ Age, strip = strip_colours, scale = "free_x") +
    labs(color = "Age", fill = "Age") +
    ylab("Predicted energy expenditure (kJ)") +
    xlab("Hopped") +
    theme(strip.background=element_blank(),
          panel.background = element_rect(fill = "white"),  
          axis.line = element_line(color = "black"),
          axis.text = element_text(face = "bold", size = 11),
          legend.text = element_text(size = 14), 
          legend.title = element_text(size = 14),
          strip.text = element_text(size = 14),
          axis.title = element_text(face = "bold", size = 15),
          legend.position = "none") +
    scale_color_manual(values = custom_palette) +
    scale_fill_manual(values = custom_palette) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .25)

关键修改说明 #

1. 叠加原始数据点：取消geom_point的注释，将y映射改为原始数据的kJs字段，Filtered_data的Hopped因子水平与绘图数据一致，分面会自动匹配Age分组完成对齐。
2. 解决带宽错误：在stat_halfeye中添加bw = "nrd0"参数，使用更稳健的带宽估计方法；若仍有问题，可尝试adjust=2（放大带宽）或density_estimator="histogram"改用直方图展示分布。

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