A random effect nested in another random effect in R with mgcv package

Question

I am interested in examining the relationship between yearly income and student success, taking into account the hierarchical structure of the data. The data includes schools as random effects, with classrooms nested within schools, and there may be correlations among observations within the same school and classroom.

To address this, my objective is to develop a generalized additive model that considers student success as the response variable, with yearly income as a fixed effect, and includes a random effect for schools and also for classrooms nested within schools.

Here is the proposed code for fitting the model using the mgcv package in R:

library(mgcv)

gam_model <- mgcv::gam(success ~ s(incomeperyear) + s(classroom, bs = "re", by = school),
                       data = my_data,
                       method = "REML",
                       family = gaussian())

summary(gam_model)

Furthermore, considering that schools and classrooms have distinct numeric values, I am uncertain whether it is required to convert them into factors prior to executing the provided code. I am unsure whether a random effect for a continuous variable is statistically applicable. Hence, I am unclear whether the mgcv package will automatically treat each numeric level as a distinct category if I do not specify the classroom and school variables as factors.

I would appreciate any insights or recommendations regarding the accuracy of the provided code and the appropriate handling of numeric variables as factors in this context.

Answer 1

If you have many schools, the factor by smooths may be appropriate, but you might get better/different results using a simple nested random effect.

From a technical point of view the factor by smooths are not themselves shrunk towards zero (as a normal random effect would). Hence you model would need to be modified to:

success ~ s(school, bs = "re") + s(incomeperyear) +
  s(classroom, bs = "re", by = school)

where we include the school random effect (first smooth) to model the school means/intercepts.

Then this model would allow for different levels of shrinkage of the classrooms towards their respective school means.

If schools and classrooms are coded uniquely, then you can also do:

success ~ s(classroom, bs = "re") + s(school, bs = "re") +
  s(school, classroom, bs = "re") + s(incomeperyear)

but this doesn't have different amounts of shrinkage per school.

The key point for both these models is to have classrooms coded uniquely (A1, A2, ..., B1, B2, ... Z1, ... Z10), where A, B, ... Z code for the schools and 1, 2, etc code for classrooms with schools.

Both school and classroom must be coded as a factor for either of these models to work, regardless of whether they are numeric or character. If you don't you'll get a varying coefficient model (for the by term) and a single linear random slope. Which is not what you want.

Finally, how is success measured? Can it be negative? Your family = gaussian() would allow that. If success is something else (a positive score, a category, a percentage, a probability, etc?) then it is unlikely that the Gaussian family is appropriate for such data.