Add support for rstanarm (generalized) models with random effects

[Selbosh]

Currently, crosspred has support out of box for rstanarm models, e.g. a stan_glm, since these model objects have coef and vcov methods defined that behave much in the same way as for an ordinary glm or gam.

However, the function does not handle Stan models with random effects terms in them (stan_glmer), the switch from coef to fixef is only activated for certain classes of models. I suggest a simple tweak, either by changing class(model) to the more robust form inherits(model, 'lmerMod') or otherwise.

I can make a quick MWE and pull request later.

[Selbosh]

Here is the problematic line. https://github.com/gasparrini/dlnm/blob/3068104d7ab137e1e2908f698ef560f00e4fe5bb/R/getcoef.R#L11-L13 A stan_glmer object (but not a stan_lmer or lmer) will pass this test and so coef(model) rather than fixef(model) will be called.

[gasparrini]

Hi David, thanks for spotting this issue. Can you please send me a simple reproducible example (maybe by email at antonio.gasparrini@lshtm.ac.uk), so I can have a look? I recognise this code could be improved, but it's a messy context, especially when dealing with lme4 stuff. I can then post a solution and possibly fixing the function. best -Antonio

[Selbosh]

Sure thing, I will post a MWE shortly.

In the meantime, the workaround for such Stan models is to pass args coef = fixef(model)[-1] and vcov = vcov(model)[-1, -1]; so that the (Intercept) term is omitted and the matrices conform properly.

[Selbosh]

library(dlnm) # 2.4.7
library(lme4) # 1.1.35.1
library(rstanarm) # 2.21.1

set.seed(2024)

data <- data.frame(
  x = rnorm(1000),
  y = rpois(1000, 1),
  subj = 1:10
)

cb <- crossbasis(data$x, lag = 5, argvar = list(fun = 'lin'), arglag = list(fun = 'bs', degree = 3))

freq_model <- glmer(y ~ cb + (1 | subj), data = data, family = poisson)
freq_pred <- crosspred(cb, freq_model)

bayes_model <- stan_glmer(y ~ cb + (1 | subj), data = data, family = poisson)
bayes_pred <- crosspred(cb, bayes_model)

The final line throws the error:

Error in crosspred(cb, bayes_model) : 
  coef/vcov not consistent with basis matrix. See help(crosspred)

But workaround is:

bayes_pred <- crosspred(cb, coef = fixef(bayes_model)[-1], vcov = vcov(bayes_model)[-1, -1])

which runs fine.

Then we can perform downstream tasks, e.g.

plot(freq_pred, 'slices', var=list(x=1))
plot(bayes_pred, 'slices', var=list(x=1))

[gasparrini]

Yes, this is correct. As stated in the help page of crosspred(), the function automatically works only with a selection of regression model objects, in addition to those with methods for coef()/vcov(). This is not the case for stan_glmer.

The proposed solution is to use the coef/vcov arguments to pass the parameters after manual extraction, which is what you suggest. Bear in mind to add the argument model.link="log" in the call to crosspred(), in order to exponentiate the results accordingly with the Poisson family.