I’m trying to get my head around using smoothing splines with brms.
I generated a simple dataset (sin wave plus normal noise):
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:
- How should I go about trying to reduce divergent transitions for brms in this example?
- 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:
Fit using brms:
Fit using mgcv (REML):
