4
$\begingroup$

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?

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

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.