# Multilevel model with 4 levels?

My dataset is a pretty typical educational dataset: we have data about students, courses, faculty, schools, and I plan to include partially-crossed random effects because students could enroll in multiple courses from the same faculty member or in courses from multiple faculty members.

I have predictors at all levels except for school (school n = 22). Is a four-level model ideal? Or should I run a three-level model and include the school data as fixed dummy variables?

My understanding is that I can't use lme4 for four-level models--so if I should be modeling the data with four levels, what are packages that might help me do so?

• Why do you think lme4 cannot fit a 4-level model ? Mar 5, 2019 at 19:21
• I've heard that from people with much more lme4 experience than me, but are we wrong? Mar 6, 2019 at 0:12
• Yes, I believe so Mar 6, 2019 at 13:51

There is no restriction to the number of "levels" in lme4. The package will be able to fit any number of levels provided that the data supports such a random effects structure.

We can demonstrate with the following simulation of a 4-level dataset similar to that as described in the OP:

> set.seed(15)
> library(lme4)
> dt1 <- data.frame(expand.grid(SchoolID = LETTERS[1:6], FacultyID = LETTERS[1:6], CourseID = LETTERS[1:10], StudentID = 1:100, Score = c(NA, NA, NA)))
> dt1$$Score <- as.numeric(dt1$$SchoolID) + as.numeric(dt1$$FacultyID) + as.numeric(dt1$$CourseID) + as.numeric(dt1$StudentID) + rnorm(nrow(dt1), 0,5) > lmm1 <- lmer(Score ~ 1 + (1 | SchoolID/FacultyID/CourseID/StudentID), data = dt1) > summary(lmm1) Random effects: Groups Name Variance Std.Dev. StudentID:(CourseID:(FacultyID:SchoolID)) (Intercept) 841.6574 29.0113 CourseID:(FacultyID:SchoolID) (Intercept) 0.8581 0.9263 FacultyID:SchoolID (Intercept) 2.5579 1.5993 SchoolID (Intercept) 2.8880 1.6994 Residual 24.9743 4.9974 Number of obs: 108000, groups: StudentID:(CourseID:(FacultyID:SchoolID)), 36000; CourseID:(FacultyID:SchoolID), 360; FacultyID:SchoolID, 36; SchoolID, 6  We could also fit a 5-level model if we wished: > dt2 <- data.frame(expand.grid(CityID = LETTERS[1:6], SchoolID = LETTERS[1:6], FacultyID = LETTERS[1:6], CourseID = LETTERS[1:10], StudentID = 1:20, Score = c(NA, NA, NA))) > dt2$$Score <- as.numeric(dt2$$CityID) + as.numeric(dt2$$SchoolID) + as.numeric(dt2$$FacultyID) + as.numeric(dt2$$CourseID) + as.numeric(dt2$$StudentID) + rnorm(nrow(dt2), 0, 5) > lmm2 <- lmer(Score ~ 1 + (1 | CityID/SchoolID/FacultyID/CourseID/StudentID), data = dt2) > summary(lmm2) Random effects: Groups Name Variance Std.Dev. StudentID:(CourseID:(FacultyID:(SchoolID:CityID))) (Intercept) 34.778 5.897 CourseID:(FacultyID:(SchoolID:CityID)) (Intercept) 7.418 2.724 FacultyID:(SchoolID:CityID) (Intercept) 2.516 1.586 SchoolID:CityID (Intercept) 2.873 1.695 CityID (Intercept) 2.922 1.709 Residual 24.940 4.994 Number of obs: 129600, groups: StudentID:(CourseID:(FacultyID:(SchoolID:CityID))), 43200; CourseID:(FacultyID:(SchoolID:CityID)), 2160; FacultyID:(SchoolID:CityID), 216; SchoolID:CityID, 36; CityID, 6  [ Note that this 2nd model may take a while to fit ! ] The partially crossed structure would be best represented by ensuring that the factors in each clusters are coded uniquely and lme4 should then be able to handle the partially crossed / partially nested structure simply by specifying the random effects as (1 | SchoolID) + (1 | FacultyID) + (1 | CourseID) + (1 | StudentID)  This means that, for example, if you have StudentID 1 in Faculty A and Student 1 in Faculty B and these are different (ie, these 2 students are nested in their respective Faculties), then they should be coded as something like StudentID 1A and StudentID 1B respectively. We can demonstrate this with the dt1 dataset above, by re-coding the factors as follows: > dt1.1 <- dt1 > dt1.1$$FacultyID <- paste(dt1$$SchoolID, dt1$$FacultyID, sep = ".") > dt1.1$$CourseID <- paste(dt1.1$$FacultyID, dt1$$CourseID, sep = ".") > dt1.1$$StudentID <- paste(dt1.1$$CourseID, dt1$StudentID, sep = ".")
> lmm1.1 <- lmer(Score ~ 1 + (1 | SchoolID) + (1 | FacultyID) + (1 | CourseID) + (1 | StudentID), data = dt1.1)
> summary(lmm1.1)

Random effects:
Groups    Name        Variance Std.Dev.
StudentID (Intercept) 841.6568 29.0113
CourseID  (Intercept)   0.8584  0.9265
FacultyID (Intercept)   2.5585  1.5995
SchoolID  (Intercept)   2.8893  1.6998
Residual               24.9743  4.9974
Number of obs: 108000, groups:  StudentID, 36000; CourseID, 360; FacultyID, 36; SchoolID, 6


Note that the model output is the same as for lmm1 above, although presented slightly differently.

So far the data are fully nested. That is, each student is enrolled on 1 and only 1 course, one course "belongs" to one and only 1 Faculty etc. To simulate a crossed factor, for example a student that is enrolled on 2 courses, we can simply combine the relevant student IDs: First we identify the student IDs that we want to combine:

> dt1.1[dt1.1$$StudentID == "A.A.A.31" | dt1.1$$StudentID == "A.A.B.31", ]
SchoolID FacultyID CourseID StudentID    Score
10801        A       A.A    A.A.A  A.A.A.31 33.00600
10837        A       A.A    A.A.B  A.A.B.31 33.69633
46801        A       A.A    A.A.A  A.A.A.31 33.03089
46837        A       A.A    A.A.B  A.A.B.31 33.00802
82801        A       A.A    A.A.A  A.A.A.31 41.68804
82837        A       A.A    A.A.B  A.A.B.31 31.26155


and gives them the same (unique) ID:

> dt1.1[dt1.1$$StudentID == "A.A.A.31" | dt1.1$$StudentID == "A.A.B.31", ]\$StudentID   <- "CCCC"


And then we can fit the model with the same same call:

lmm1.1 <- lmer(Score ~ 1 + (1 | SchoolID) + (1 | FacultyID) + (1 | CourseID) + (1 | StudentID), data = dt1.1)
> summary(lmm1.1)
Random effects:
Groups    Name        Variance Std.Dev.
StudentID (Intercept) 841.6867 29.0118
CourseID  (Intercept)   0.8312  0.9117
FacultyID (Intercept)   2.5570  1.5991
SchoolID  (Intercept)   2.8851  1.6986
Residual               24.9742  4.9974
Number of obs: 108000, groups:  StudentID, 35999; CourseID, 360; FacultyID, 36; SchoolID, 6


Note that we now have 35,999 StudentIDs, rather than 36,000.