Hello,
I am fitting splines in brms, and I am wondering how I can plot the smooth of the difference between two individual smooths in a binary by-factor smooth (the difference smooth). This is very easy to do using mgcv, but I’m not sure how to get the same results (and plots) with brms. I’m also wondering how to interpret the brms output with a by-factor smooth when I use an ordered binary factor.
Here’s an example of what I am trying (using the mtcars dataset)
mtcars$vs ← as.ordered(mtcars$vs)
brms_mod ← brm(mpg ~ s(hp) + vs + s(hp, by=vs), data = mtcars)
gam_mod ← gam(mpg ~ s(hp) + vs + s(hp, by=vs), data = mtcars)
As I understand it in mgcv, if I use an ordered factor in a by-factor smooth, the model will produce the difference in the smooth between two categories of the factor (if it’s binary). I assume this would be the same in brms, but the output looks a bit different, so I’m not 100% sure.
Here is the output from the brms model above (brms_mod):
And here is the screenshot from the mgcv model (gam_mod):
From what I understand for splines in brms models, the sds() parameters for the smooth terms are the estimates of the “wiggliness”. With an ordered factor in brms, would the term “sds(shpvs1_1)” represent the wiggliness of the difference smooth? In the mgcv model, I would interpret the term s(hp):vs1 as the estimate of the difference smooth. This term in the mgcv model has an EDF of 1, so there’s no wiggliness in the difference smooth.
I would also like to be able to plot the difference between two smooths. With an ordered binary factor, mgcv produces this automatically. For the mgcv model above, I can obtain this plot from the mgcv model above with plot(gam_mod):
When I use conditional_effects(brms_mod, “hp:vs”) to obtain the smooth plots from the brms model , I obtain each smooth within each factor level on the same plot. Here is what I get:
This plot is very helpful for visualizing the results, but I would also like to see the difference smooth (with intervals) from brms. Is there a way to obtain that with conditional_effects()?
Thank you!