I'm trying to plot the density curves from a brms model and add a line showing the mean predicted value. I've calculated the mean in two different ways, and the yeild slightly different results. Am I missing something? How do I know which is correct?
library(tidyverse)
library(brms)
library(tidybayes)
set.seed(42)
df <- data_frame(site = rep(LETTERS[1:4], each = 25),
sex = rep(c("f", "m"), 50),
length = rnorm(100, 20, 2),
mass = rnorm(100, 50, 5))
m <- brm(mass ~ length + site + sex,
data = df,
prior = c(prior(normal(0, 10), class = Intercept),
prior(normal(0, 1), class = b),
prior(normal(0, 1), class = sigma)))
# obtain predicted values from posterior distribution draws
pred <- data_frame(site = rep(LETTERS[1:4], each = 2),
sex = rep(c("f", "m"), 4),
length = 20) %>%
add_predicted_draws(m)
# APPROACH 1: calculate mean value from predicted values
predM <- group_by(pred, site, sex) %>% summarise(mP = mean(.prediction))
# APPROACH 2: calculate mean value from posterior summary
postSum <- posterior_summary(m)
calcM <- data_frame(site = rep(LETTERS[1:4], each = 2),
sex = rep(c("f", "m"), 4),
length = 20) %>%
mutate(mC = postSum[1] + # intercept
postSum[2] * length + # slope
case_when(site == "B" ~ postSum[3], # site-specific intercept adjustments
site == "C" ~ postSum[4],
site == "D" ~ postSum[5],
TRUE ~ 0) +
if_else(sex == "m", postSum[6], 0)) # sex-specific intercept adjustment
ggplot(data = pred, aes(x = .prediction, col = sex)) + geom_density() +
geom_vline(data = predM, aes(xintercept = mP, col = sex), linetype = "dashed") +
geom_vline(data = calcM, aes(xintercept = mC, col = sex), linetype = "solid") +
facet_wrap(~site)
