3
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ I don't think so, that would be ~1|family/school. $\endgroup$ Nov 12, 2021 at 12:47
  • $\begingroup$ That does not run unfortunately, I get the following error: Error in MEEM(object, conLin, control$niterEM) : NAs in foreign function call (arg 2) In addition: Warning messages: 1: In ncols * isLast : longer object length is not a multiple of shorter object length 2: In ncols * c(rep(1, Q), 0, 0) : longer object length is not a multiple of shorter object length $\endgroup$
    – Abdel
    Nov 14, 2021 at 9:46

1 Answer 1

2
+50
$\begingroup$

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
$\endgroup$
5
  • $\begingroup$ In a previous version of this answer you see that the method did not work. The difference was in using the pdIdent function, which forces the random slope effects to be drawn from a distribution with the same variance. $\endgroup$ Nov 14, 2021 at 15:21
  • $\begingroup$ Thank you! So the current version of the answer does work, is this reliable? $\endgroup$
    – Abdel
    Nov 15, 2021 at 9:42
  • 1
    $\begingroup$ @Abdel, yes, it is tricky to get these crossed effects specified, but finally, it works. For this specific example, it is not so difficult to get to it (I picked it from somebody else) and you can check easily that it works (e.g. you can also check the estimated variances and see that they are the same). But the coding is still a bit like magic to me. I wonder if there is any good manual that explains the use of these positive definite matrix functions to specific the mixed models. $\endgroup$ Nov 15, 2021 at 9:43
  • $\begingroup$ Thank you very much! It works on the example data indeed, but on my real data (which is very similar, but much larger) I seem to be getting the following error: Error in MEEM(object, conLin, control[dollar_sign]niterEM) : NAs in foreign function call (arg 2) In addition: Warning message: In MEEM(object, conLin, control[dollar_sign]niterEM) : NAs introduced by coercion to integer range) $\endgroup$
    – Abdel
    Nov 15, 2021 at 13:51
  • $\begingroup$ @Abdel that error message is not so clear. I sometimes look up in the code where it comes from when I do not understand the error message and debug it by digging into it and go through everything step by step. Something must be causing those NAs. Possibly you could try to pinpoint it quickly by reducing your data and increase it until you find the point where it fails. Maybe there is some false input. $\endgroup$ Nov 20, 2021 at 13:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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