Brms GAM intercept too low

Modeling
,
1

Hi!

I’m using brms to fit a GAM with 8 covariates, only one smooth term, and all seems well and good in terms of model performance. But when I plot the resulting smooth term together with my raw data the intercept seems to sometimes (not always!) be too low; see the attached image where the purple lines (model) are very low relative to the points which are binned averages of raw data.

tmp
tmp733×462 5.27 KB

It’s a mixed effects model where there’s a lot of dispersion, and it’s all positive data that is being modelled with a lognormal family.

brms call, simplified and in generic terms:
tmp.mod ← brm(formula = value ~ s(x1, bs = ‘tp’) + x2 + x3 + (1 | Area),
data = dat, family = lognormal(),
warmup = 1000, iter = 4000, chains = 4, cores = 4,
prior = tmp.priors,
control = list(adapt_delta = 0.99, max_treedepth = 12))

priors (not sure if they matter in this case):
tmp.priors ← c(set_prior(‘normal(0, 0.5)’, class = ‘b’),
set_prior(‘normal(-3, 2.5)’, class = ‘Intercept’),
set_prior(‘normal(0.5, 0.15)’, class = ‘sigma’))

It doesn’t matter whether calling conditional_effects or extracting the outcome using epred_draws, when holding x2 and x3 at their mean, and re_formula = NA. I’ve tried using the formula above as well as including “0 + Intercept” before the other terms.

Is there anyone knows why this is, and whether I’m simply ignorant of something that makes this behaviour explainable/logical or whether there’s something going wrong?

Let me know if I can provide additional information (I’m not allowed to share/upload the data)!

R version 4.2.2
brms 2.20.4
rstan 2.32.3

2

can you show us your plotting code?

3

Hi and thanks for the quick reply! Did not expect such a quick one so haven’t checked on the topic until today, but will check in more often from now :)

Yes, of course. In the image in the original post, each colored line is a separate model. So to plot the smooth of (for example) the purple line model, the simplest is just:
conditional_smooths(tmp.mod, effects = ‘x1’)

With the epred_draws function I’ve tried putting the co-variates at both their mean and median level and that only makes a minor difference, i.e. does not solve the problem.

tmp.dat <- bind_cols(lapply(tmp.dat[, c('x2', 'x3')], mean))
tmp.dat$Area <- dat$Area[1]
tmp.dat <- expand_grid(x1 = 1:55, tmp.dat)
tmp.dat <- epred_draws(tmp.mod, newdata = tmp.dat, re_formula = NA, scale = 'response')
tmp.dat <- bind_rows(lapply(1:55, function(x){
    quantile(as.vector(unlist(tmp.dat[tmp.dat$x1 == x, '.epred'])), probs = c(0.025, 0.05, 0.075, 0.5, 0.925, 0.95, 0.975))
  }))
colnames(tmp.dat) <- c('lwr95', 'lwr90', 'lwr85', 'avg', 'upr85', 'upr90', 'upr95')

tmp.dat %>% ggplot(aes(x = tSinceCC, y = avg)) + 
     geom_line() + 
     geom_ribbon(fill = NA, linetype = 'dashed', aes(ymin = lwr90, ymax = upr90))

Thankful for any suggestions/ideas! :)

4

Apparently I cannot edit once more, but I’ve noticed a typo in the code I pasted here. The first row should be:

tmp.dat <- bind_cols(lapply(dat[, c('x2', 'x3')], mean))

It was a typo in my post, not in the original code (I’ve double checked!).

5

Hi! There are still a few things I don’t understand here. What are the different colors on the plot and where do they come from?

That said, one possible source of the discrepancy is that there is no guarantee that the grand mean of the random effect will correspond to the mean of the data. This could happen either if there are imbalanced group sizes, or simply because the mean of the linear predictor on the link scale does not correspond to the mean of the linear predictor on the response scale. That is, exp(mean(x)) is not the same as mean(exp(x)).

6

Hi again jsocoloar!

The different colors were mentioned in my first response (“In the image in the original post, each colored line is a separate model”), but I’ll do my best to clarify: each of the different colors represent a distinct model, and different data feed into the parameterization of each of these models. You can imagine these different colors being different species, countries, or whatever grouping variable that needs its own model parameterizaton because they have unique respones not only to x1, but also to x2 and x3. I could of course use random slopes for all these terms, but that’s a different story and perhaps best to not discuss that here. The solid line is the 50th percentile of the epred_draws, and the dashed lines the 5th and 95th percentiles of the epred_draws (credible interval), as indicated from the code I posted on 20th of May.

Ok. Would checking that be as simple as using the ranef function and extacting the mean of the “Estimate.Intercept”, and then comparing that to the “Intercept” given from the summary function? They are indeed very different, one even being positive (0.001) and one being negative (-1.01), but I’m not sure if they’re on the same scale.
From what I can see the default link for lognormal is identity, i.e. the predictors should be centered, but otherwise on the same scale as when I send them to the model. Or am I missing something?

I’m sorry for the confusion, I really feel like I’m entering a new, to me unknown, world!

7

Thanks; I just didn’t understand the correspondence between the plot and the code posted. The upshot is that I can ignore everything that’s not purple on the plot?

I was careless in my last post: brms does use an identity link for the lognormal family, but the parameter predicted on the identity scale by the linear predictor is not the expectation of the response but the expectation of the distribution’s logarithm (NB: which is different from the logarithm of the expectation of the response, so this isn’t quite a log link).

8

@questionmark. Why use:

conditional_smooths(tmp.mod, effects = ‘x1’)

instead of:

conditional_effects(tmp.mod, effects = ‘x1’)

Just a thought.

9

You can indeed ignore the not-purple ones. I just found it strange that the model was giving me outcomes that I expected in some cases, but not in others: orange, green and blue lines/models are giving epred_draws that are closer to the data compared to the purple, where the intercept seems to be offset.

Okay, so the intercept that is being shown in both the summary function and the ranef function are intercepts for “the expectation of the distribution’s logarithm” when all covariates are 0 (when not using the 0 + Intercept syntax)? That would mean that the output from the summary and ranef functions are one the same scale?

In any case, is there any way to handle what you describe as “the mean of the linear predictor on the link scale does not correspond to the mean of the linear predictor on the response scale”?
There might be a slight imbalance in groups, but the non-purple models in my original post would suffer the same imbalance as the other color models.

@MilaniC I could do that to get the functional form/shape of the smooth, but that plot does not come out on the response scale. From other posts I’ve read this plot shows the smooth under the “sum-to-zero” constraint, and without taking any intercepts into account. But if possible to adjust the plot from conditional_smooths so that it comes out on the response scale then that would of course of be good.

10

@questionmark.

conditional_effects(tmp.mod, effects = ‘x1’) is on the response scale — else what scale do you think it is on?

11

@MilaniC
I know that, and in my original post I mention that I’ve tried both conditional_effects and epred_draws:
“It doesn’t matter whether calling conditional_effects or extracting the outcome using epred_draws, when holding x2 and x3 at their mean, and re_formula = NA. I’ve tried using the formula above as well as including “0 + Intercept” before the other terms.”
But I thought you were suggesting to use conditional_smooths (which I have not said that I’m using), not conditional_effects. But it seems you suggested I try something I’ve already tried (conditional_effects).