(This is a fairly long answer, there is a summary at the end)
You are not wrong in your understanding of what nested and crossed random effects are in the scenario that you describe. However, your definition of crossed random effects is a little narrow. A more general definition of crossed random effects is simply: not nested. We will look at this at the end of this answer, but the bulk of the answer will focus on the scenario you presented, of classrooms within schools.
First note that:
Nesting is a property of the data, or rather the experimental design, not the model.
Also,
Nested data can be encoded in at least 2 different ways, and this is at the heart of the issue you found.
The dataset in your example is rather large, so I will use another schools example from the internet to explain the issues. But first, consider the following over-simplified example:
Here we have classes nested in schools, which is a familiar scenario. The important point here is that, between each school, the classes have the same identifier, even though they are distinct if they are nested. Class1 appears in School1, School2 and School3. However if the data are nested then Class1 in School1 is not the same unit of measurement as Class1 in School2 and School3. If they were the same, then we would have this situation:
which means that every class belongs to every school. The former is a nested design, and the latter is a crossed design (some might also call it multiple membership. Edit: For a discussion of the differences between multiple membership and crossed random effects, see here ), and we would formulate these in lme4 using:
(1|School/Class) or equivalently (1|School) + (1|Class:School)
and
(1|School) + (1|Class)
respectively. Due to the ambiguity of whether there is nesting or crossing of random effects, it is very important to specify the model correctly as these models will produce different results, as we shall show below. Moreover, it is not possible to know, just by inspecting the data, whether we have nested or crossed random effects. This can only be determined with knowledge of the data and the experimental design.
But first let us consider a case where the Class variable is coded uniquely across schools:
There is no longer any ambiguity concerning nesting or crossing. The nesting is explicit. Let us now see this with an example in R, where we have 6 schools (labelled I-VI) and 4 classes within each school (labelled a to d):
> dt <- read.table("http://bayes.acs.unt.edu:8083/BayesContent/class/Jon/R_SC/Module9/lmm.data.txt",
header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)
> # data was previously publicly available from
> # http://researchsupport.unt.edu/class/Jon/R_SC/Module9/lmm.data.txt
> # but the link is now broken
> xtabs(~ school + class, dt)
class
school a b c d
I 50 50 50 50
II 50 50 50 50
III 50 50 50 50
IV 50 50 50 50
V 50 50 50 50
VI 50 50 50 50
We can see from this cross tabulation that every class ID appears in every school, which satisfies your definition of crossed random effects (in this case we have fully, as opposed to partially, crossed random effects, because every class occurs in every school). So this is the same situation that we had in the first figure above. However, if the data are really nested and not crossed, then we need to explicitly tell lme4:
> m0 <- lmer(extro ~ open + agree + social + (1 | school/class), data = dt)
> summary(m0)
Random effects:
Groups Name Variance Std.Dev.
class:school (Intercept) 8.2043 2.8643
school (Intercept) 93.8421 9.6872
Residual 0.9684 0.9841
Number of obs: 1200, groups: class:school, 24; school, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 60.2378227 4.0117909 15.015
open 0.0061065 0.0049636 1.230
agree -0.0076659 0.0056986 -1.345
social 0.0005404 0.0018524 0.292
> m1 <- lmer(extro ~ open + agree + social + (1 | school) + (1 |class), data = dt)
summary(m1)
Random effects:
Groups Name Variance Std.Dev.
school (Intercept) 95.887 9.792
class (Intercept) 5.790 2.406
Residual 2.787 1.669
Number of obs: 1200, groups: school, 6; class, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 60.198841 4.212974 14.289
open 0.010834 0.008349 1.298
agree -0.005420 0.009605 -0.564
social -0.001762 0.003107 -0.567
As expected, the results differ because m0 is a nested model while m1 is a crossed model.
Now, if we introduce a new variable for the class identifier:
> dt$classID <- paste(dt$school, dt$class, sep=".")
> xtabs(~ school + classID, dt)
classID
school I.a I.b I.c I.d II.a II.b II.c II.d III.a III.b III.c III.d IV.a IV.b
I 50 50 50 50 0 0 0 0 0 0 0 0 0 0
II 0 0 0 0 50 50 50 50 0 0 0 0 0 0
III 0 0 0 0 0 0 0 0 50 50 50 50 0 0
IV 0 0 0 0 0 0 0 0 0 0 0 0 50 50
V 0 0 0 0 0 0 0 0 0 0 0 0 0 0
VI 0 0 0 0 0 0 0 0 0 0 0 0 0 0
classID
school IV.c IV.d V.a V.b V.c V.d VI.a VI.b VI.c VI.d
I 0 0 0 0 0 0 0 0 0 0
II 0 0 0 0 0 0 0 0 0 0
III 0 0 0 0 0 0 0 0 0 0
IV 50 50 0 0 0 0 0 0 0 0
V 0 0 50 50 50 50 0 0 0 0
VI 0 0 0 0 0 0 50 50 50 50
The cross tabulation shows that each level of class occurs only in one level of school, as per your definition of nesting. This is also the case with your data, however it is difficult to show that with your data because it is very sparse. Both model formulations will now produce the same output (that of the nested model m0 above):
> m2 <- lmer(extro ~ open + agree + social + (1 | school/classID), data = dt)
> summary(m2)
Random effects:
Groups Name Variance Std.Dev.
classID:school (Intercept) 8.2043 2.8643
school (Intercept) 93.8419 9.6872
Residual 0.9684 0.9841
Number of obs: 1200, groups: classID:school, 24; school, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 60.2378227 4.0117882 15.015
open 0.0061065 0.0049636 1.230
agree -0.0076659 0.0056986 -1.345
social 0.0005404 0.0018524 0.292
> m3 <- lmer(extro ~ open + agree + social + (1 | school) + (1 |classID), data = dt)
> summary(m3)
Random effects:
Groups Name Variance Std.Dev.
classID (Intercept) 8.2043 2.8643
school (Intercept) 93.8419 9.6872
Residual 0.9684 0.9841
Number of obs: 1200, groups: classID, 24; school, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 60.2378227 4.0117882 15.015
open 0.0061065 0.0049636 1.230
agree -0.0076659 0.0056986 -1.345
social 0.0005404 0.0018524 0.292
It is worth noting that crossed random effects do not have to occur within the same factor - in the above the crossing was completely within school. However, this does not have to be the case, and very often it is not. For example, sticking with a school scenario, if instead of classes within schools we have pupils within schools, and we were also interested in the doctors that the pupils were registered with, then we would also have nesting of pupils within doctors. There is no nesting of schools within doctors, or vice versa, so this is also an example of crossed random effects, and we say that schools and doctors are crossed. A similar scenario where crossed random effects occur is when individual observations are nested within two factors simultaneously, which commonly occurs with so-called repeated measures subject-item data. Typically each subject is measured/tested multiple times with/on different items and these same items are measured/tested by different subjects. Thus, observations are clustered within subjects and within items, but items are not nested within subjects or vice-versa. Again, we say that subjects and items are crossed.
Summary: TL;DR
The difference between crossed and nested random effects is that nested random effects occur when one factor (grouping variable) appears only within a particular level of another factor (grouping variable). This is specified in lme4 with:
(1|group1/group2)
where group2 is nested within group1.
Crossed random effects are simply: not nested. This can occur with three or more grouping variables (factors) where one factor is separately nested in both of the others, or with two or more factors where individual observations are nested separately within the two factors. These are specified in lme4 with:
(1|group1) + (1|group2)
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.
Best Answer
It might be that
lmeis not build for crossed random effects. I am not sure about this, but at least, it is very difficult to find examples that use it in such a way.However, there is a working method from this message by Peter Dalgaard
Here you are tricking
lmeby creating a group with a single level. The model is still nested but now it is the single level group that is part of the nesting which is no problem.The crossed random effect is incorporated by treating the random effect for the one group as a parameter for a slope instead of an intercept.
You can compare it with the same model in
lme4::lmerWith the data for your example, the random effects have little influence. Below is a data set that makes the comparison easier.
The models give the same results. For instance, see the random effects: