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),
m1 <- lme(score~
random= list(int = pdIdent(form = ~ 0 + as.factor(school)) ,family = ~1),
control = lmeControl(opt = "optim"),
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
m2 <- lme4::lmer(score ~ training_1 + training_2 + (1|family) + (1|school), data = data)
With the data for your example, the random effects have little influence. Below is a data set that makes the comparison easier.
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:
with conditional variances for “family” “school”
as.factor(school)1 as.factor(school)2 as.factor(school)3
1 0.3500093 -0.7326827 0.3826734