Edit (to address the point raised in the first comment below):
This random effects structure:
(1 | eth1/eth2/eth3)
Expands to:
(1 | eth1) + (1 | eth1:eth2) + (1 | eth1:eth2:eth3)
However, provided that the nesting is explicit in the coding of factors, this is also equivalent to :
(1 | eth1) + (1 | eth2) + (1 | eth3)
A simple simulation will demonstrate this:
> set.seed(15)
> dtA <- expand.grid(eth1 = c(1,2,3,4,5), eth2 = c(1,2,3,4,5), eth3 = c(1,2,3,4,5,6), measure = c(1,2))
> dtA$y <- dtA$eth1 + dtA$eth2 + dtA$eth3 + dtA$measure + rnorm(nrow(dtA),0,1)
> fm01 <- lmer(y ~ 1 + (1|eth1/eth2/eth3), data = dtA)
> summary(fm01)
Random effects:
Groups Name Variance Std.Dev.
eth3:(eth2:eth1) (Intercept) 3.038 1.743
eth2:eth1 (Intercept) 1.718 1.311
eth1 (Intercept) 1.948 1.396
Residual 1.377 1.173
Number of obs: 300, groups: eth3:(eth2:eth1), 150; eth2:eth1, 25; eth1, 5
However, the "nesting" here is not explicit. Each level of eth2
occurs in every level of eth1
:
> xtabs(~ eth1 + eth2, dtA)
eth2
eth1 1 2 3 4 5
1 12 12 12 12 12
2 12 12 12 12 12
3 12 12 12 12 12
4 12 12 12 12 12
5 12 12 12 12 12
Thus it is ambiguous whether these data are crossed or nested. Nesting is a property of the experimental design. So, if the data are nested, then eth2
= 1 in eth1
= 1 is not the same unit of measurement as eth2
= 1 in eth1
= 2. If the data are nested and the factors are not coded uniquely between clusters, then it is necessary to write the random effect structure as (1|eth1/eth2/eth3)
or (1|eth1) + (1|eth1:eth2) + (1|eth1:eth2:eth3)
which tells lmer
that the data are nested.
If we fit this model:
> fm02 <- lmer(y ~ 1 + (1|eth1) + (1|eth2) + (1|eth3), data = dtA)
...then we obtain different estimates for the variance components, because unless we specify the nesting explicitly, it means that the random effects are crossed, and not nested:
> summary(fm02)
Random effects:
Groups Name Variance Std.Dev.
eth3 (Intercept) 3.132 1.770
eth2 (Intercept) 2.232 1.494
eth1 (Intercept) 2.395 1.547
Residual 1.291 1.136
Number of obs: 300, groups: eth3, 6; eth2, 5; eth1, 5
Now, let us encode the factors uniquely:
> dtA$eth2.u <- paste(dtA$eth1, dtA$eth2, sep=".")
> dtA$eth3.u <- paste(dtA$eth1, dtA$eth2, dtA$eth3, sep=".")
> xtabs(~ eth1 + eth2.u, dtA)
eth2.u
eth1 1.1 1.2 1.3 1.4 1.5 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 4.1 4.2 4.3 4.4 4.5 5.1 5.2 5.3 5.4 5.5
1 12 12 12 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 0 0 0 0 0 12 12 12 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0 0 0 12 12 12 12 12 0 0 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 12 12 12 12 0 0 0 0 0
5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 12 12 12 12
and now we see that the nesting is explicit.
So to complete the demonstration, we fit:
> fm03 <- lmer(y ~ 1 + (1|eth1/eth2.u/eth3.u), data = dtA)
or equivalently:
> fm04 <- lmer(y ~ 1 + (1|eth1) + (1|eth2.u:eth1) + (1|eth3.u:eth2.u:eth1), data = dtA)
then we obtain:
> summary(fm03)
Random effects:
Groups Name Variance Std.Dev.
eth3.u:(eth2.u:eth1) (Intercept) 3.038 1.743
eth2.u:eth1 (Intercept) 1.718 1.311
eth1 (Intercept) 1.948 1.396
Residual 1.377 1.173
Number of obs: 300, groups: eth3.u:(eth2.u:eth1), 150; eth2.u:eth1, 25; eth1, 5
And finally, if we fit:
> fm05 <- lmer(y ~ 1 + (1|eth1) + (1|eth2.u) + (1|eth3.u), data = dtA)
we obtain:
> summary(fm03)
Random effects:
Groups Name Variance Std.Dev.
eth3.u (Intercept) 3.038 1.743
eth2.u (Intercept) 1.718 1.311
eth1 (Intercept) 1.948 1.396
Residual 1.377 1.173
Number of obs: 300, groups: eth3.u, 150; eth2.u, 25; eth1, 5
..which indeed is exactly the same as fm01
, fm03
and fm04
.
The main take-home point here is that when factors are coded uniquely between clusters, lmer
will take care of the nesting, and the model can be specified in several ways, depending on what is most convenient, but if the factors are not unique then the nesting must either be made explicit by changing the coding of the factors, as we did above, or the nesting must be specified using /
or :
when writing the grouping variables in the call to lmer
.