1

Background

When fitting generalised additive models (GAMs) in the past, I have used methods for the mgcv package outlined in

Pedersen EJ, Miller DL, Simpson GL, Ross N (2019) Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7: e6876.

I am currently transitioning to Bayesian alternatives and preferentially use the rethinking package, but have stumbled into some problems when introducing splines to Stan and adding categorical predictor variables.

MRE

Below I follow the code outlined in this book:

McElreath R (2019) Statistical rethinking2: A Bayesian course with examples in R and Stan. Second Edition. Chapman and Hall/CRC

and this Bookdown:

Kurz S (2023) Statistical rethinking with brms, ggplot2, and the tidyverse: Second edition. https://bookdown.org/content/4857

# generate some sinusoidal data
d <- data.frame(x = rep(seq(0, 5 * pi, length.out = 100), 2),
                group = rep(c("a", "b"), each = 100))
d$y[1:100] <- 2 * sin(2 * d$x[1:100]) + rnorm(100, 0, 1)
d$y[101:200] <- 5 * sin(3 * d$x[101:200]) + rnorm(100, 0, 0.5)

# generate b splines
knots <- 15
knot_list <- quantile(d$x , probs = seq(0, 1, length.out = knots))
require(splines)
B <- bs(d$x,
        knots = knot_list[-c(1, knots)],
        degree = 3, intercept = TRUE)

# fit quap model according to McElreath (2019)
require(rethinking)
r.mod <- quap(
  alist(
    y ~ dnorm(mu, sigma),
    mu <- a + B %*% w,
    a ~ dnorm(0, 1),
    w ~ dnorm(0, 1),
    sigma ~ dexp(1)
  ), data = list(y = d$y, B = B),
     start = list(w = rep(0, ncol(B)))
)
precis(r.mod) # model produces parameter estimates

# attempt to fit identical ulam (Stan) model
r.mod.stan <- ulam( # note ulam instead of quap
  alist(
    y ~ dnorm(mu, sigma),
    mu <- a + B %*% w,
    a ~ dnorm(0, 1),
    w ~ dnorm(0, 1),
    sigma ~ dexp(1)
  ), data = list(y = d$y, B = B),
     start = list(w = rep(0, ncol(B))),
     chains = 8, cores = 8, iter = 1e4
)

# model fails due to apparently semantic error:
# Ill-typed arguments supplied to assignment operator =: lhs has type real and rhs has type row_vector

# fit brms model according to Kurz (2023)
require(brms)
d$B <- B
b.mod.stan <- brm(data = d,
              family = gaussian,
              y ~ 1 + B,
              prior = c(prior(normal(0, 1), class = Intercept),
                        prior(normal(0, 1), class = b),
                        prior(exponential(1), class = sigma)),
              iter = 1e4, warmup = 5e3, chains = 4, cores = 4,
              seed = 4)
summary(b.mod.stan) # model produces same estimates as quap model

Questions

  1. Why does the rethinking Stan GAM not run while the brms GAM does?
  2. How can I incorporate the grouping factor group?

I would appreciate a worked example in rethinking.

Shawn Hemelstrand
  • 2,676
  • 4
  • 17
  • 30
Luka Seamus Wright
  • 75
  • 6
  • 1
    While you can drive Stan via Richard's rethinking package, you have to remember that it is motivated as a tool for learning about Bayesian modelling and not necessarily something you should use as a daily driver for fitting models. *brms* is a bit different in that regard. I wouldn't fit this model in exactly the same way as Richard does (I've never understood why the model and the basis both include an intercept term for example) as this isn't penalising the coefficients of the smooth in anything like how *mgcv* does it. Solomon's translation of the rethiking code to brms is just that 1/2 – Gavin Simpson May 08 '23 at 11:27
  • 1
    ...a direct translation, so again it isn't something I would recommend for a production model, largely because the spline is not penalised (beyond the prior on the coefs coercing the coefs towards 0, slightly.) Instead, I would use `brm(y ~ group + s(x, k = 15, bs = "cr"), ...)` to include a penalised spline in the model, and then you'd need to specify the other bits as you had them, but also look at the priors for terms like this; there is another variance term you need a prior on, which controls the wiggliness of the fitted spline). 2/2 – Gavin Simpson May 08 '23 at 11:30
  • Your Q is off topic here because it's ostensibly about coaxing software to do what you want and not a statistical problem associated with software. – Gavin Simpson May 08 '23 at 11:33
  • @GavinSimpson, thanks for your response. I understand that `ulam` is also a pedagogical tool but it has been applied in the scientific literature and I intend to do so too. I like that its syntax is closer to that of Stan, which I would eventually like to learn. I appreciate the line of `brm` and it's similarity to `gam`. However, while I can see the underlying Stan code using `stancode`, I do not have the necessary knowledge of Stan to be able to recreate your `brm` model in `ulam`. I have flagged the question to be moved to StackOverflow where I'll hopefully get a worked example in `ulam`. – Luka Seamus Wright May 10 '23 at 05:57

0 Answers0