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.
See here for more info about coding crossed and nested factors.
