Divergent transitions fitting smooths in brms

Modeling
,
1

I’m trying to get my head around using smoothing splines with brms.

I generated a simple dataset (sin wave plus normal noise):

dataset1
dataset12100×1499 111 KB

And fit the following brms model, for B-splines with 10 knots:

fit_brms <- brm(
  y ~ s(x, bs="bs", k=10),
  data = df, 
  family = gaussian,
  backend = "cmdstanr",
  refresh = 0
)

This fits quickly and the posterior looks reasonable (see plots below), but I ended up with 27 of 4000 divergent transitions. I would like to understand and eliminate these transitions!

In order to understand what was going on I decided to code the model myself in stan:

data {
    int n;
    int len_a;
    vector[n] y;
    matrix[n, len_a] X;
    matrix[len_a, len_a] S;
}

parameters {
    vector[len_a] a_tilde;
    real<lower=0> sigma_s;
    real<lower=0> sigma;
}

transformed parameters {
    vector[len_a] a = a_tilde * sigma_s;
    vector[n] mu = X * a;
}

model {
    target+= - quad_form(S, a_tilde);
    target+= normal_lpdf(y | mu, sigma);
}

I’m including the penalisation via the prior: p(a|\sigma_s) = exp(- \frac{a^T S a}{2 \sigma_s^2}). The a_tilde term is to avoid divergent transitions by using a non-centred parametrisation.

I get the model matrix, X, and the penalisation matrix, S, using mgcv’s function smoothCon in R:

sm_spec <- s(x, bs = "bs", k = 10)

sm_obj <- smoothCon(sm_spec, data = df, absorb.cons = TRUE)[[1]]

X <- sm_obj$X
S <- sm_obj$S[[1]]

Fitting this stan model, again the posterior looks reasonable and now I only get 1 to 2 divergent transitions in 4000 which I am more comfortable with getting rid of by increasing adapt_delta.

I know of two big differences between my parametrisation and brms’s. First, I have flat priors on all parameters, whereas brms has a set of default priors (I’m not sure how to replicate equivalent priors to the brms ones in mine).

Second, I understand (from this blog by Tristan Mahr: Random effects and penalized splines are the same thing - Higher Order Functions ), that under the hood brms uses mgcv to convert the penalized spline problem to a mixed effects model.

My questions are as follows:

  1. How should I go about trying to reduce divergent transitions for brms in this example?
  2. brms have chosen a sensible parametrisation for computational reasons, if I used my stan code for more complex problems in the future, would I run into problems?

Fit using my stan code:

stan_gam1
stan_gam12100×1499 109 KB

Fit using brms:

brms_gam1
brms_gam12100×1499 114 KB

Fit using mgcv (REML):

mgcv_gam1
mgcv_gam11054×870 38.5 KB

2

Quick update to this, I forgot to include an intercept term, new model is:

transformed parameters {
    vector[len_a] a = a_tilde * sigma_s;
    vector[n] mu = a_0 + X * a;
}

model {
    target+= cauchy_lpdf(sigma_s |0, 2.5);
    target+= - quad_form(S, a_tilde);
    target+= normal_lpdf(y | mu, sigma);
}

I had to include a prior on sigma_s to avoid divergences but otherwise this didn’t change the previous comparison.

3

I’ve had similar issues implementing various types of smooths – a few divergent transitions but the overall fit looks fine. I suspect it has something to do with either potential collinearity between the basis columns and/or scaling of the design matrix. Can you post a pairs plot of your posterior draws? It might be also worth looking into doing a QR decomposition of X, but in that case I’m not sure how the penalty matrix S would be derived.

4

Thanks for the reply. Here is the pairs plot for a selection of brms parameters. It looks OK to me.

pairs
pairs646×470 122 KB

Do you mean for the version I wrote or brms?

5

I’m not sure what brms does under the hood. I meant for the Stan implementation – if you’re doing an unpenalized spline then you can just swap the design matrix into something like Regression Models, which could help if collinearity between the basis functions is an issue. But since you’re then doing inference on the transformed parameters, I don’t know if the penalty matrix formulation of the likelihood would still behave as intended.

6

Smooths tend to have the problem that you get funnels with both centered and non-centered parameterizations as illustrated by @Niko in his blog post Divergent transitions in Hilbert Space Gaussian process posteriors and how to avoid them. Although the blog post uses HSGP as the example, the same problem happens with the usual spline implementations. Sometimes increasing adapt_delta a little bit can help, but sometimes it may require switching to adaptive parameterization, but as it is not automatically used by brms, it may require some effort.

7

Thanks for the answer and the link to the blog post, that’s very helpful for my understanding.