Fixing Random Effects Issues in Multilevel Modeling with mgcv: A Simple Solution

The problem with the code is that it’s not properly modeling the random effects. The bs = "re" argument in the smooth function implies that it’s a random effect model, but the predict function doesn’t understand this and instead treats it as if it were a fixed effect.

To fix this, you need to exclude the terms you consider ‘random’ from the prediction using the exclude argument in the predict function. This will ensure that only the fixed effects are included in the predictions.

Here’s the corrected code:

library(mgcv)
bam_ex.2 <- bam(y ~ 
                 factor5 + 
                 factor4 + 
                 s(x, by = interaction(factor5, factor4), k= 3) + # one smooth per each level of the two factors
                 s(factor1, by = factor4, bs = "re") + # factor1 is nested in factor4, factor 1 is random
                 s(factor2, by = factor3, bs = "re") + # factor2 is nested in factor3,  factor 2 is random       
                 s(factor3, bs = "re"), # factor 3 is random
               data = data_exp,
               discrete = TRUE)

pred_ex.2 <- data_exp %>% 
  mutate(y.p = predict(bam_ex.2,
                        exclude=c("s(factor1):factor4 DE",
                                  "s(factor1):factor4 PE",
                                  "s(factor2):factor3 A15",
                                  "s(factor2):factor3 A16",
                                  "s(factor3)"),
                        type="response"))

By excluding the terms with random effects, you’re ensuring that only the fixed effects are included in the predictions, which should fix the issue.

Last modified on 2024-07-16