Nested effects in lmer model

Question

I have the following data structure:

  X       R          Y       Z         ID    YZ          error
  <fct>   <fct>      <fct>   <fct>     <fct> <fct>       <dbl>
1 G       1          18      6        3      18,6        20.1
2 G       1          36      12       3      36,12       35.9 
3 G       2          18      6        3      18,6        4.0
4 G       2          36      12       3      36,12       6.0
5 M       1          18      6        3      18,6        27.5 
6 M       1          36      12       3      36,12       10.3 
7 M       2          18      6        3      18,6        1.5 
8 M       2          36      12       3      36,12       8.5 

I need to determine the main effects of XYZ, using an ANOVA. From my understanding, Z is nested in Y (or vice versa), because the combinations 36:6 and 18:12 do not exist. ID is the subject identifier and R is the measurement repetition in a repeated measures design.

So when I run the lmer model:

m = lmer(error ~ (X*Y*Z) + (1|ID/R), data = data_transform)

I get the warning: "fixed-effect model matrix is rank deficient", which makes sense because not all combinations of Y and Z exist, so lmer drops those combinations. So far, this not a problem, however, when running post-hoc tests:

emm <- emmeans(m, ~ X | Z)
post_hoc_tests <- summary(pairs(emm, adjust = "bonferroni"))

all the estimates are NA. It must have something to do with the nesting, since

m = lmer(error ~ (X*Z) + (1|ID/R), data = data_transform)

and running emmeans does produce results. Maybe it could be solved by creating the variable YZ, which contains all combinations but then I cannot figure out how to test the effects of Y and Z separately.

My questions are:

1. How can I specify the nesting in the model correctly to get results from emmeans?
2. Does it make sense to specify the random effects like this in a repeated measures design?

Update: While trying to construct a compact minimum example of the data, I failed to mention that Y and Z actually have more than 2 levels each. More specifically, Y has 6 levels and Z has 11 levels, leading to a total of 66 possible combinations, where only 18 of them exist in the data. So consequently, YZ has 18 distinct levels. Here is an example (which does not include all combinations):

    X       R          Y       Z         ID    YZ          error
      <fct>   <fct>      <fct>   <fct>     <fct> <fct>       <dbl>
    1  G       1          18      6        3      18,6        20.1
    2  G       1          36      12       3      36,12       35.9 
    3  G       1          27      10       3      27,10       5.4
    4  G       2          18      6        3      18,6        4.0
    5  G       2          36      12       3      36,12       6.0
    6  G       2          27      10       3      27,10       8.8
    7  M       1          18      6        3      18,6        27.5 
    8  M       1          36      12       3      36,12       10.3
    9  M       1          27      10       3      27,10       1.1
    10 M       2          18      6        3      18,6        1.5 
    11 M       2          36      12       3      36,12       8.5
    12 M       2          27      10       3      27,10       3.3

Update:

I found that "telling" emmeans about the nesting works:

    grid = ref_grid(m, nesting = "Y %in% Z*X")
    post_hoc_tests = summary(pairs(emmeans(grid, ~ X | Z), adjust = "bonferroni"))

produces results.

Answer 1

Based on your data structure, $Y=18$ whenever $Z = 6$ and $Y = 36$ whenever $Z = 12$. In this case, these effects are not nested but instead are collinear. There is not enough variation for us to infer a nesting structure. To determine which factor is nested in the other, we need at least one more factor level; for example if $Y = 18$ was associated with both $Z = 6$ and $Z = 7$ (and these values never occurred when $Y = 36$).

The warning message you get has to do with this collinearity. You can visualize the problem yourself if you are interested. In order to fit a factor variable in a linear model, we need to construct indicator variables. So we add an intercept to the model, which is a column of all ones. Then we include an indicator variable $I_y = I(Y = 36)$ and an indicator variable $I_z = I(Z = 12)$. Neglecting the other variables for a moment (because they do not change what we are about to see), you can see that $I_y = I_z$. That means these two variables are perfectly collinear and the effects are not identifiable because the model cannot tell which variable is causing the variation in the outcome.

Now, if we instead construct the indicator $I_{z} = I(Z = 6)$, you can see that $I_y + I_z = \langle 1, 1, 1, \ldots, 1 \rangle$. This is collinear with the intercept, so again, we lose identifiability in the model. If we take out the intercept, this would be estimable but does not disentangle the effects of $Y$ and $Z$. Instead, it would give us separate effects of the first level and the second level, without having a reference level.

Think about it this way: if we know the value of $Y$ we automatically know the value of $Z$. There is no variability for the model to detect, so it can only detect the combined effect of $Y$ and $Z$ at the same time.