Common hyperparameter for splines by categories

Interfaces brms
1

Hello!

In brms, when setting a cubic spline, for example

y ~ s (x, "cr", k =4, by = cat)

Weights of basis are normally distributed with estimated standard-deviation. When using a spline for each factor level, a standard-deviation is estimated for each factor level.

I’d like to estimate a common standard-deviation of basis for each factor level. Is this possible? I tried to use the “id” argument, but the code still produces a “sds” per category.

Thank you very much!
Lucas

2

Hello!

I’ll define my question more precisely. Here is my data generating process :

library(tidyverse)
library(splines)
library(brms)

D <- tibble(x = rep(seq(from=0, to=15, by=2),6),
            cat = rep(LETTERS[1:3],each = 16))# generating inputs

B <- t(bs(X, knots=seq(-5,5,1), degree=3, intercept = TRUE)) # creating the B-splines
num_data <- length(X); num_basis <- nrow(B)
a0 <- matrix(c(0.2,0.2,0.5,0.5,-0.2,-0.2), ncol = 6) # intercept
a_A1 <- rnorm(num_basis, 0, 1) # coefficients of B-splines
a_A2 <- rnorm(num_basis, a_A1, 0.2)
a_B1 <- rnorm(num_basis, 0, 1) # coefficients of B-splines
a_B2 <- rnorm(num_basis, a_B1, 0.2)
a_C1 <- rnorm(num_basis, 0, 1) # coefficients of B-splines
a_C2 <- rnorm(num_basis, a_C1, 0.2)
A <- matrix(c(a_A1, a_A2, a_B1, a_B2, a_C1, a_C2), ncol = 6)
Y_true <- as.vector(t(a0)%*%X + t(A)%*%B) # generating the output
D$y <- Y_true + rnorm(length(X),0,.2) # adding noise

D %>% ggplot(aes(x = x, y = y, color = cat)) + 
  geom_point() + 
  geom_smooth()

In math :

y_i \sim Normal(\mu_i, \sigma)\\ \mu_i = \beta_c x_i + f_c(x)\\ f_c(x) = \sum^K_{k=1} \alpha_{ck}b_k\\ \alpha_{ck} \sim normal(0, \tau)

Where c is the index of the category at which the observation i belongs. There is one standard deviation of weights common to each category (\tau).

However, when I fit this with brms using the following code, I get an sd estimated for each level of my category, such as :

\alpha_{ck} \sim normal(0, \tau_c) 
# Store formula
form_spline <- bf(y ~ cat + s(x, by = cat, k = 3))
make_stancode(form_spline, data = D)

Is it possible to estimate a unique standard-deviation, as in my data generating process? I am afraid my real data does not contain enough information to estimate C standard-deviation without causing divergences.

I tried the following approach, based upon the “id” argument from s(), which states

An id with a factor by variable causes the smooths at each factor level to have the same smoothing parameter.

Unfortunately, the following code still produces C standard-deviation to be estimated.

form_spline <- bf(y ~ cat + s(x, by = cat, k = 3, id = "a"))
make_stancode(form_spline, data = D)

Is someone aware of any solution?

Thank you very much!
Lucas

1 Like
3

Unfortunately, I don’t think this is currently possible from the brms side unless there is a functionality for this already built into mgcv.

4

Hi Paul!

The “id” argument of s () seems to be related, but I do not know how it is treated internally.

I thought afterward I could easily modify the brms produced code, using the “empty fit” functionality! I have not tested this approach yet, but I think it should work!

Alle the best,
Lucas

5

If I recall correctly: if the smoothing parameter is the same for all factors then it is equal to the random smooth setting. I will see if I find any support for this claim.

Regarding the id argument: it seems like there is a bug in the mgcv code for representing smooths as random effects? Since brms relies on that, that could be the origin of the discrepancy between brms’ behaviour and the stated behaviour in the mgcv help files.

6

Hello!

I do not really understand what you meant (because I think I do not know what is the random smooth setting).

I mark this post as solution because I simply modified the brms generated stan code. I retained only the first sds_1_1 and mutiply each un-centered basis by it, removing all the others sds. I sampled from the posterior with rstan and pushed the resulting fit in an empty brms object.

It worked like a charm! I get rid of divergent transitions with an adapt delta of 0.8, and with the number of basis I thought was straightforward.

Thanks everybody for help!
Lucas

7

The setting you are describing, i.e. the same smoothness for all factors, correspond as I understand it to model S in Pedersen et al. (2019) and could be fitted with “fs” as basis, with formulas corresponding to y∼s(x, fac, bs="fs") (for one-dimensional smoothers) or y∼t2(x1, x2, fac, bs=c("tp", "tp", "re") (for two-dimensional smoothers). The one-dimensional case would also be possible to specify as y∼s(x, fac, bs=c("tp","re")). You could use other bases instead of “tp” in this setting if you want cyclic smoothers that vary by some factor with the same smoothness, although thin plate splines and "re" correspond to "fs". The factor should be stored as a factor in R, but I guess you know that.

I highly recommend the paper a read, although it is focused only on mgcv.

Pedersen, E. J., Miller, D. L., Simpson, G. L., & Ross, N. (2019). Hierarchical generalized additive models in ecology: An introduction with mgcv. PeerJ, 7, e6876. https://doi.org/10.7717/peerj.6876

2 Likes
8

Hi Staffan,

Thank you very much for your explanation and the excellent reference. I saw it when it went out, but I completely forgot it!

I tried the y∼s(x, fac, bs=c("tp","re")) approach, but it still produces a sds per factor level in brms.

All the best,
Lucas

1 Like