How to translate (relate) **mgcv** `gam()` and **brms** `brm()` results

Interfaces brms
1

I am keen to get a detailed grasp of the relationship of some non-linear models analysed with mgcv vs. brms. The syntax is virtually the same, apart from the functions themselves gam() vs. brm(). But the summary tables differ. To illustrate this, I am using the first toy-example from mgcv:

require(mgcv)
require(brms)

# example from mgcv
set.seed(2)
dat <- gamSim(1, n=400, dist='normal', scale=2)

# gam from mgcv
summary(gam1 <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data=dat))
# Parametric coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  7.83328    0.09878    79.3   <2e-16

# Approximate significance of smooth terms:
#         edf Ref.df      F  p-value
# s(x0) 2.500  3.115  6.921 0.000131
# s(x1) 2.401  2.984 81.858  < 2e-16
# s(x2) 7.698  8.564 88.158  < 2e-16
# s(x3) 1.000  1.000  4.343 0.037818

# brm from brms
summary(brm1 <- brm(y ~ s(x0) + s(x1) + s(x2) + s(x3), data=dat))
# Smooth Terms: 
#            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sds(sx0_1)     2.18      1.19     0.74     5.11 1.00     1856     2533
# sds(sx1_1)     2.05      1.11     0.70     4.88 1.00     1928     2416
# sds(sx2_1)    19.74      5.03    12.16    32.31 1.00     1142     1741
# sds(sx3_1)     0.90      0.91     0.03     3.23 1.00     1468     1731

# Population-Level Effects: 
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept     7.83      0.10     7.64     8.03 1.00     5138     3185
# sx0_1        -1.06      3.77    -8.79     6.21 1.00     2113     2201
# sx1_1         9.31      3.74     1.90    17.20 1.00     2019     2222
# sx2_1        35.47     11.85    11.58    58.22 1.00     2581     2655
# sx3_1        -0.61      2.44    -4.89     5.45 1.00     1881     1648

# Family Specific Parameters: 
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     1.98      0.07     1.84     2.13 1.00     4643     3062

So, all four smooths are significant according to the gam. Equally, by brm, CIs do not contain zeros. Nicely, x3 is flagged weak by both, so to say.

However, how should I interpret Population-Level Effects from brm? Are these simple comparisons with the Intercept? Is there more to it? Interestingly enough, besides x3 x0 also does not differ from the Intercept.

  • Operating System: macOS Monterey (12.5.1)
  • brms Version: 2.17.0

Thanks!

2

@striatum hello, I think the answer is to stop trying to interpret the coefficientsā€¦ If you wanted to compare the predictions from each of these models then Iā€™d examine each of the smooth functions directly, thereā€™s a clean shortcut for this from brms with conditional_smooths() and of course you can show any predictions you want with fitted() or predict().

I believe that both of these models are equivalently thin-plate regression splines using s(), by default. The default priors in brms are not all flat though, so you wouldnā€™t expect the estimates to be exactly the same, even for the same model.

1 Like
3

Hey @AWoodward, thanks so much for this! It makes perfect sense. What I like to know/have is a sort of initial ā€œtriageā€ of the smooths worth further inspection with conditional_smooths(). For that, I suppose, CIs should suffice. (By the way, very often, one needs to make reviewers happy rather than present/summarise what is actually meaningful. šŸ™„ šŸ˜‰)

4

See this post and the article by Tristan Mahr linked there. As @AWoodward said, the population-level effects are not very interpretable.
Also, brms treats the smooths as varying effects (aka random effects) and the hyperparameter that controls the wiggliness is like the sd in a varying effects model (that is what those sds() are). So this is more like gamm4.

1 Like