R: Fit Linear Mixed-Effects Model with GAM (mgcv)

Question

Imagine a linear mixed effects model with one random intercept:

library(lme4)
LMM1 <- lmer(response ~ experience + (1|subject), REML=FALSE, data=train)

I'd like to fit a conceptual equivalent using the mgcv package, that is:

library(mgcv)
LMM2 <- bam(response ~ experience, s(subject, bs='re'), method='ML', data=train)

Note that I use the function bam which is optimized for big datasets. The summaries of the two models yield the same results w.r.t. the estimated intercept, and effect size, standard deviation and t-value of the fixed effect. However, AIC(LMM1) and AIC(LMM2) yield different results. Why?

In addition, I'm wondering whether the below listed pairs of models are conceptual equal:

Random intercept:
lmer(response ~ experience + (1|subject), data=train)
bam(response ~ experience + s(subject, bs=’re’), data=train)

Random slope:
lmer(response ~ experience + (0+experience|subject), data=train)
bam(response ~ experience + s(experience, subject, bs=’re’), data=train)

Random intercept and random slope:
lmer(response ~ experience + (1+experience|subject), data=train)
bam(response ~ experience + s(subject, bs='re') + s(experience, subject, bs=’re’), data=train)

If not, please put me right.

Answer 1

Those two model definitions should be equivalent. However, you have to be careful comparing likelihood's between software implementations as different constants can often be ignored in the computation; not sure that's the case here but it is something to be aware of. What may be causing differences here is that mgcv is adjusting the AIC to account for the selection of the "smoothness parameter" for the random effects terms; here this basically just means choosing the amount of shrinkage, but mgcv uses a correction to the AIC for smoothness selection and that could be the cause of the difference.

Of the other model forms, these are all equivalent, except for the random slopes and intercepts model. mgcv can't fit correlated random effect models and the version of the random slopes and intercepts models you showed has correlated random effects. The lmer() version of the model has an additional covariance parameter which will be missing from the bam() fit.

You can fit and compare uncorrelated random effects models using the alternate form for lmer():

lmer(response ~ experience + (1 | subject) + (0 + experience | subject), 
     data = train)

bam(response ~ experience + s(subject, bs = 're') +
      s(experience, subject, bs = ’re’),
    data = train)