-
Notifications
You must be signed in to change notification settings - Fork 65
Open
Description
I'd like to use spread_draws(), or some other function if more appropriate, to pull parameters in a coef() metric (i.e., population mean + random effect). I can do it, but the code seems clunky and I'm wondering if there's a more elegant way. Here's an example of what I mean.
# Load the packages
library(tidyverse)
library(brms)
library(tidybayes)
# Load the data
data(ChickWeight)
# Adjust the data
cw <- ChickWeight |>
data.frame() |>
rename_all(.funs = tolower) |>
mutate(time_z = (time - mean(time)) / sd(time),
weight_z = (weight - mean(weight)) / sd(weight))
# Fit a model with several random effects
my_fit <- brm(
data = cw,
weight_z ~ 1 + time_z + I(time_z^2) + (1 + time_z + I(time_z^2) | chick),
prior = prior(normal(0, 1), class = Intercept) +
prior(normal(0, 1), class = b) +
prior(exponential(1), class = sd) +
prior(exponential(1), class = sigma) +
prior(lkj(2), class = cor),
cores = 4, seed = 1
)
# Extract the posteriors, wrangle, and plot
beta_levels <- c("Intercept", "time_z", "Itime_zE2")
beta_labels <- str_c("beta[", 0:2, "]")
my_fit |>
spread_draws(b_Intercept, b_time_z, b_Itime_zE2, r_chick[chick, ranef] | ranef) |>
# 50 chicks is overkill
filter(chick <= 10) |>
pivot_longer(cols = c(b_Intercept:b_Itime_zE2, Intercept:time_z), names_to = "beta") |>
mutate(type = ifelse(str_detect(beta, "b_"), "fixef", "ranef"), # This line is not needed, but can bring conceptual clarity
beta = str_remove(beta, "b_"),
chick = factor(chick)) |>
group_by(.draw, chick, beta) |>
summarise(value = sum(value)) |>
mutate(beta = factor(beta, levels = beta_levels, labels = beta_labels)) |>
ggplot(aes(x = value, y = chick)) +
stat_pointinterval(.width = 0.95,
linewidth = 1/4, size = 1) +
facet_wrap(~ beta, labeller = label_parsed, nrow = 1)
Any thoughts on better versions of the spread_draws() part of the code for combining the fixed and random effects?
Reactions are currently unavailable
Metadata
Metadata
Assignees
Labels
No labels