Hi,
First of all, thanks a lot for the implementation of this package (@paul.buerkner) and all incredible books and posts (@Solomon ) that are available for newbies like me in the field of probabilistic modelling, this is amazing!
I would like to reproduce Figure 4.10 from Solomon Kurz’s book: Statistical rethinking with brms, ggplot2, and the tidyverse (original Fig 4.8 of Statistical rethinking). The idea is to plot the joint consequence of μ AND σ with predict(), but in this case, for a multilevel model in which time is a categorical variable for seven time points. I found very useful the examples of the book with the quadratic and cubic model version and apply them to my multilevel data:
Example data:
t ← (rep(0:6,10))
ab ← ((runif(1, 10, 20) * t) / (runif(1, 0, 10) + t)) + rnorm(10, 0, 1)
plot(t, ab,pch=19,las = 1)
X ← as.data.frame(cbind(t,ab));
X$t_fact ← as.factor(t)
X$id ← paste(“id”,rep(0:6,each=10),sep="_")
Cubic model
fitCub ←
brm(data = X,
family = gaussian,
ab ~ 0 + Intercept + t + I(t^2) + I(t^3) + (1 | id),
chains = 4, cores = 4, iter = 2000, warmup = 1000,
seed = 9,
file = “Ex_cubic”)
Data to evaluate predictions:
time_seq ←
tibble(t = seq(from = 0, to = 6, length.out = 40)) %>%
mutate(ItE2 = t^2) %>%
mutate(ItE3 = t^3)
id<-unique(X$id)
time_seq ← time_seq %>%
expand(time_seq,id)
Fitted values
fitd_cub ← fitted(fitCub, newdata = time_seq) %>%
as_tibble() %>%
bind_cols(time_seq)
Predicted responses
pred_cub ← predict(fitCub, newdata = time_seq) %>%
as_tibble() %>%
bind_cols(time_seq)
Plot the data
ggplot(data = X, aes(x = t)) +
geom_ribbon(data = pred_cub, aes(ymin = Q2.5, ymax = Q97.5),fill = “#cfd8e6”) +
geom_smooth(data = fitd_cub, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = “identity”, fill = “#a6a3ff”, color = “darkblue”, alpha = .5) +
geom_point(aes(y = ab),color = “navyblue”, shape = 19, size = 1.5, alpha = 1/3) +
coord_cartesian(xlim = range(X$t),ylim = range(X$ab)) +
theme_bw()
However, the final aim is to obtain the predicted values for the next model, in which the time is treated as a categorical variable, with a coefficient per time point:
Model
fitGT ←
brm(data =X,
family = gaussian,
formula = ab ~ 0 + Intercept + t_fact + (1| id),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
control = list(adapt_delta = 0.9),
seed = 9,
file = “Ex_GT”)
fitted
fitd_GT ← fitted(fitGT) %>%
as_tibble() %>%
bind_cols(X)
predict
pred_GT ← predict(fitGT) %>%
as_tibble() %>%
bind_cols(X)
Plot the data
ggplot(data = X, aes(x = t)) +
geom_ribbon(data = pred_GT, aes(ymin = Q2.5, ymax = Q97.5),fill = “#cfd8e6”) +
geom_smooth(data = fitd_GT, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = “identity”, fill = “#a6a3ff”, color = “darkblue”, alpha = .5) +
geom_point(aes(y = ab),color = “navyblue”, shape = 19, size = 1.5, alpha = 1/3) +
coord_cartesian(xlim = range(X$t),ylim = c(-10,10)) +
theme_bw()
In this case, I wonder if, given the nature of time as a factor, it is only possible to get the values for the time points already defined in the model, as is done in the example above. I really would like to obtain “smoother” intervals. Is there any way to get that? Is the code properly implemented?
Any help/guidance much appreciated!
- Operating System: ubuntu
- brms Version: brms_2.14.4