Hey @ucfagls, thank you for the helpful answer. I apologize for the response lag. After putting this on the shelf for a while, and with the help from Steve Wild, it’s now clear how to make this work with brms. For posterity, I’ll work it out in two ways:
The single-level option
You’re right to bring up the identifiability issue. Happily, even modest priors solve that problem within the Bayesian context.
Load the focal packages and the data.
library(tidyverse)
library(brms)
library(splines)
data(cherry_blossoms, package = "rethinking")
d <- cherry_blossoms
rm(cherry_blossoms)
# drop missing cases
d2 <-
d %>%
drop_na(doy)
Define the knots and the b-splines.
num_knots <- 15
knot_list <- quantile(d2$year, probs = seq(from = 0, to = 1, length.out = num_knots))
B <- bs(d2$year,
knots = knot_list[-c(1, num_knots)],
degree = 3,
intercept = TRUE)
Make a new data set which includes B as a matrix column.
d3 <-
d2 %>%
mutate(B = B)
Fit the model.
b4.8 <-
brm(data = d3,
family = gaussian,
doy ~ 1 + B,
prior = c(prior(normal(100, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(exponential(1), class = sigma)),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
seed = 4)
Here’s the model summary.
print(b4.8)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: doy ~ 1 + B
Data: d3 (Number of observations: 827)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 103.59 2.49 98.58 108.45 1.01 761 1071
B1 -3.19 3.88 -10.80 4.21 1.00 1591 2435
B2 -1.11 3.92 -8.82 6.62 1.00 1463 1859
B3 -1.27 3.68 -8.58 6.07 1.00 1332 1833
B4 4.56 2.96 -1.13 10.54 1.00 1046 1648
B5 -1.08 2.99 -6.89 4.78 1.00 954 1399
B6 4.03 3.02 -1.86 10.07 1.00 1056 1515
B7 -5.55 2.91 -11.23 0.30 1.00 981 1672
B8 7.57 2.92 1.92 13.29 1.00 1005 1583
B9 -1.22 2.97 -6.90 4.58 1.01 1001 1458
B10 2.79 3.03 -3.17 9.02 1.00 1019 1527
B11 4.42 3.00 -1.49 10.27 1.00 1056 1846
B12 -0.37 2.99 -5.98 5.70 1.00 1013 1551
B13 5.28 2.99 -0.52 11.07 1.00 1024 1768
B14 0.46 3.09 -5.44 6.74 1.00 1036 1712
B15 -1.05 3.37 -7.77 5.59 1.00 1215 2147
B16 -7.25 3.44 -14.01 -0.49 1.00 1287 1936
B17 -7.88 3.26 -14.24 -1.29 1.00 1246 1824
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 5.94 0.14 5.67 6.23 1.00 4173 2921
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Make the bottom panel of McElreath’s Figure 4.13 (p. 118).
fitted(b4.8) %>%
data.frame() %>%
bind_cols(d3) %>%
ggplot(aes(x = year, y = doy, ymin = Q2.5, ymax = Q97.5)) +
geom_hline(yintercept = fixef(b4.8)[1, 1], linetype = 2) +
geom_point(color = "steelblue") +
geom_ribbon(alpha = 2/3) +
labs(x = "year",
y = "day in year") +
theme_classic()
The multilevel option
It seems like we can get a pretty close s()-based analogue with a little tweeking. Here’s my attempt.
b4.11 <-
brm(data = d2,
family = gaussian,
doy ~ 1 + s(year, bs = "bs", k = 19),
prior = c(prior(normal(100, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(student_t(3, 0, 5.9), class = sds),
prior(exponential(1), class = sigma)),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
seed = 4,
control = list(adapt_delta = .99))
Parameter summary:
posterior_summary(b4.11)[-22, ] %>% round(2)
Estimate Est.Error Q2.5 Q97.5
b_Intercept 104.54 0.20 104.14 104.95
bs_syear_1 -0.09 0.33 -0.75 0.57
sds_syear_1 1.31 0.59 0.54 2.80
sigma 5.98 0.15 5.70 6.28
s_syear_1[1] -1.37 1.67 -5.50 1.02
s_syear_1[2] -0.38 0.10 -0.60 -0.19
s_syear_1[3] 0.40 0.24 -0.09 0.88
s_syear_1[4] 1.19 0.42 0.40 2.08
s_syear_1[5] 0.78 0.61 -0.43 2.04
s_syear_1[6] 0.11 0.78 -1.38 1.79
s_syear_1[7] 0.48 0.93 -1.36 2.39
s_syear_1[8] 0.33 1.01 -1.59 2.44
s_syear_1[9] 0.95 1.00 -0.96 3.09
s_syear_1[10] -0.80 0.95 -2.82 1.02
s_syear_1[11] -0.43 0.92 -2.38 1.33
s_syear_1[12] -0.20 0.93 -2.08 1.69
s_syear_1[13] 0.57 1.06 -1.33 2.88
s_syear_1[14] 1.04 1.21 -0.92 3.82
s_syear_1[15] -0.03 1.10 -2.33 2.20
s_syear_1[16] 1.47 1.46 -0.82 4.86
s_syear_1[17] 1.54 1.54 -0.77 5.38
Make the s()-based alternative to the bbottom panel of McElreath’s Figure 4.13.
fitted(b4.11) %>%
data.frame() %>%
bind_cols(d3) %>%
ggplot(aes(x = year, y = doy, ymin = Q2.5, ymax = Q97.5)) +
geom_hline(yintercept = fixef(b4.11)[1, 1], linetype = 2) +
geom_point(color = "steelblue") +
geom_ribbon(alpha = 2/3) +
labs(x = "year",
y = "day in year") +
theme_classic()
Here we might compare the bias weights from the two models with a coefficient plot.
bind_rows(
# single level
fixef(b4.8)[-1, -2] %>%
data.frame() %>%
mutate(number = 1:n(),
fit = "single level") %>%
select(fit, number, everything()),
# multilevel (with `s()`)
posterior_samples(b4.11) %>%
select(bs_syear_1, contains("[")) %>%
pivot_longer(-bs_syear_1, names_to = "number") %>%
mutate(bias = bs_syear_1 + value,
number = str_remove(number, "s_syear_1") %>% str_extract(., "\\d+") %>% as.integer()) %>%
group_by(number) %>%
summarise(Estimate = mean(bias),
Q2.5 = quantile(bias, probs = 0.025),
Q97.5 = quantile(bias, probs = 0.975)) %>%
mutate(fit = "multilevel") %>%
select(fit, number, everything())
) %>%
mutate(number = if_else(fit == "multilevel", number + 0.1, number - 0.1)) %>%
# plot
ggplot(aes(x = Estimate, xmin = Q2.5, xmax = Q97.5, y = number, color = fit)) +
geom_pointrange(fatten = 1) +
scale_color_viridis_d(option = "A", end = 2/3, direction = -1) +
theme_classic()
I’ll be uploading a more detailed walk-out of this workflow in my ebook, soon. In the meantime, the updated code lives on GitHub.
