Skip to main content
Addressing the point in the first comment
Source Link
Robert Long
  • 65.9k
  • 11
  • 133
  • 248

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.


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.

Bounty Ended with 200 reputation awarded by Lepidopterist
Source Link
Robert Long
  • 65.9k
  • 11
  • 133
  • 248

My question is, can one usefully specify a multilevel-model with this "partially-nested, partially non-nested" structure? This need not be able to be fit using lmer().

Yes, you can specify a model with a nested / partially nested structure. Provided that you encode the nested/crossed/partially-crossed factors uniquely, then lmer will handle the structure correctly if you specify the grouping variables as

(1 | eth1) + (1 | eth2) + (1 | eth3) + (1|eth3:age) + (1|age) 

Of course you can add random slopes for x-centered as necessary (It may not be necessary or warranted to fit random slopes to all the factors of the nesting hierarchy).

A nice little trick, if you want to estimate the same model in a Bayesian framework, and/or you are having problems with convergence, is simply to replace the call to lmer with a call to stan_lmer from the rstanarm package, which was developed by Andrew Gelman and colleagues. This is designed as a plug-and-play replacement for lmer which fits a Bayesian model with default priors and estimates it using Hamilton Monte Carlo. If you seek to set you own priors, then this can also be done quite easily. There is a lot of documentation at https://mc-stan.org