How do I specify crossed random intercepts in lme?

Question

I want to run a random intercept mixed-effect model, with two random intercepts. I made some example data below, which consists of 10 subjects from 3 families that go to 3 different schools. They get two types of training, that I want to include as fixed effects.

Mock data:

data <- data.frame(child=rep(1:10), 
                   school=c(2,1,1,2,3,1,2,3,2,1), 
                   family=c(1,3,2,1,2,2,3,2,1,3), 
                   training_1=c(5,4,6,5,6,7,8,6,5,6), 
                   training_2=c(4,7,4,5,5,8,6,5,8,5), 
                   score=c(46,56,60,68,70,55,67,60,59,47))

I run the following lme command with training_1 & training_2 as fixed effects, and family and school as random intercepts:

m1 <- lme(score~
            training_1+
            training_2,
          random=list(~1|family, ~1|school),
          method="ML", 
          data=data)
summary(m1)

If I'm not mistaken, family and school are now analyzed as "nested", which means that a certain family appears only within a particular school. In my data however, different family members from the same family can end up in different schools, which means that I should analyze them as "crossed", right? How do I adjust my lme command to do this?

Answer 1

It might be that lme is 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

data <- data.frame(child=rep(1:10), 
                   int = rep(1,10),
                   school=c(2,1,1,2,3,1,2,3,2,1), 
                   family=c(1,3,2,1,2,2,3,2,1,3), 
                   training_1=c(5,4,6,5,6,7,8,6,5,6), 
                   training_2=c(4,7,4,5,5,8,6,5,8,5), 
                   score=c(46,56,60,68,70,55,67,60,59,47))

m1 <- lme(score~ 
                training_1+
                training_2,
          random= list(int = pdIdent(form = ~ 0 + as.factor(school)) ,family = ~1),
          control = lmeControl(opt = "optim"),
          data=data)
summary(m1)

Here you are tricking lme by 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::lmer

m2 <- lme4::lmer(score ~ training_1 + training_2 + (1|family) + (1|school), data = data)
summary(m2)

With the data for your example, the random effects have little influence. Below is a data set that makes the comparison easier.

set.seed(1)
n = 50
child = 1:n
int = rep(1,n)
school = sample(1:3,n, replace = 1)
family = sample(1:3,n, replace = 1)
schoolrand = rnorm(3)
familyrand = rnorm(3)
epsilonrand = rnorm(n)
training_1 = sample(4:9, n , replace = 1)
training_2 = sample(4:9, n , replace = 1)
score = 30+training_1*3+training_2*6+schoolrand[school]+familyrand[family]+epsilonrand[child]

data <- data.frame(child = child, int = int, school = school, family = family,
                   training_1 = training_1, training_2 = training_2, score = score)

The models give the same results. For instance, see the random effects:

> lme4::ranef(m2)
$family
  (Intercept)
1  -1.4667324
2   0.5262105
3   0.9405218

$school
  (Intercept)
1   0.3500098
2  -0.7326836
3   0.3826738

with conditional variances for “family” “school” 
> m1$coefficients$random
$int
  as.factor(school)1 as.factor(school)2 as.factor(school)3
1          0.3500093         -0.7326827          0.3826734

$family
    (Intercept)
1/1  -1.4667334
1/2   0.5262108
1/3   0.9405226