lmer Dropping levels from summary information in models with main effect and interactions

Question

First off, here are my models:

  lm1 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + (1|x10/x11), data=df, na.action=na.exclude)
  lm2 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x12 + (1|x10/x11), data=df, na.action=na.exclude)
  lm3 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12:x6 + x12:x7 + x12:x8 + x12:x9 + (1|x10/x11), data=df, na.action=na.exclude)
  lm4 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x7 + x12*x8 + x12*x9 + (1|x10/x11), data=df, na.action=na.exclude)
  lm5 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12:x6 + x12:x9:x7 + x12:x9:x8 + (1|x10/x11), data=df, na.action=na.exclude)
  lm6 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x9*x7 + x12*x9*x8 + (1|x10/x11), data=df, na.action=na.exclude)

Dataset is quite large so it's hosted here (csv saved as txt): https://1drv.ms/t/s!AvjuPq_dbORpaySu39F0vPXbv5Q?e=ah7rsi

I am struggling to understand why the addition of main effect terms in a model results in multiple levels being dropped from interaction terms as observed through model summary info.

I am aware that x12 and x9 are nominal (i.e., 2 levels and 9 levels respectively), and that some of interactions involving both variables are rank deficient. This is not the source of the problem, interactions are being dropped from the summary information which are not rank deficient.

lm4 and lm6 are variants of lm3 and lm5 respectively where main effect terms have been included along with interactions (i.e., lm3 and lm5 contain no main effect terms).

The problem is that I'm finding that some of the levels of the interaction terms, which are key in my application, are dropped from the summary information where the main effects are included (i.e., in lm4 and lm6). This problem is particularly bad in lm4 (see below). I can't figure out why inclusion of the main effect terms would warrant removal of interaction levels from the summary information.

NOTE: I use the lmerTest package to bolster the summary info. in my analysis.

  > summary(lm4)
    Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
    Formula: log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x7 + x12*x8 + x12*x9 + (1|x10/x11)
       Data: df
    
    REML criterion at convergence: 229807.1
    
    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -6.4111 -0.4604  0.0917  0.5612  5.0303 
    
    Random effects:
     Groups   Name        Variance Std.Dev.
     x10:x11  (Intercept) 0.2678   0.5175  
     x10      (Intercept) 0.2735   0.5230  
     Residual             0.2036   0.4512  
    Number of obs: 132083, groups:  PTR:PLT, 39282; PLT, 526
    
    Fixed effects:
                  Estimate Std. Error         df t value Pr(>|t|)    
    (Intercept) -5.772e-01  2.775e-01  8.191e+02  -2.080  0.03785 *  
    x1           1.089e-01  1.125e-03  7.422e+04  96.801  < 2e-16 ***
    x2          -1.881e-02  3.498e-04  8.952e+04 -53.775  < 2e-16 ***
    x3           1.498e-02  1.624e-04  9.646e+04  92.204  < 2e-16 ***
    x4          -1.305e-02  1.144e-03  4.606e+03 -11.407  < 2e-16 ***
    x5           1.888e-05  2.370e-05  4.689e+02   0.797  0.42604    
    x12NF       -4.808e-01  1.789e-01  3.366e+03  -2.687  0.00724 ** 
    x6          7.401e-03  7.154e-03  4.768e+02   1.035  0.30140    
    x7          1.139e-04  4.641e-05  9.878e+04   2.455  0.01411 *  
    x8          -8.096e-04  1.180e-04  6.469e+04  -6.858 7.03e-12 ***
    x9BW        -2.252e-01  4.507e-02  4.245e+04  -4.997 5.84e-07 ***
    x9BY        2.444e-01  1.479e-01  3.051e+04   1.653  0.09838 .  
    x9LT        -2.574e-01  1.174e-01  4.461e+04  -2.192  0.02837 *  
    x9MR        -4.036e-01  1.536e-01  3.875e+04  -2.627  0.00862 ** 
    x9PW         1.225e-01  1.421e-01  3.486e+04   0.862  0.38874    
    x9RT        -1.846e-01  1.064e-01  3.850e+04  -1.735  0.08269 .  
    x9SB        -2.126e-01  2.381e-02  4.841e+04  -8.930  < 2e-16 ***
    x9SW        -5.054e-01  9.772e-02  5.010e+04  -5.172 2.33e-07 ***
    x12NF:x6     2.415e-03  8.312e-03  4.736e+02   0.291  0.77150    
    x12NF:x7    -2.817e-04  4.660e-05  9.853e+04  -6.045 1.50e-09 ***
    x12NF:x8     8.554e-04  1.191e-04  6.624e+04   7.183 6.91e-13 ***
    x12NF:x9BW   9.307e-03  4.866e-02  4.112e+04   0.191  0.84834    
    x12NF:x9LT   2.630e-01  1.289e-01  4.242e+04   2.041  0.04128 *  
    x12NF:x9MR   1.404e-01  1.639e-01  3.842e+04   0.856  0.39179    
    x12NF:x9RT   7.611e-02  1.204e-01  3.696e+04   0.632  0.52719    
    x12NF:x9SB   3.738e-01  2.644e-02  4.469e+04  14.136  < 2e-16 ***
    x12NF:x9SW   5.014e-01  1.027e-01  4.789e+04   4.881 1.06e-06 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Correlation matrix not shown by default, as p = 27 > 12.
    Use print(x, correlation=TRUE)  or
        vcov(x)        if you need it
    
    fit warnings:
    fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients

Answer 1

Update:

As OP pointed out I missunderstood the Question:

[...] it is not the rank deficiency I am concerned with. The reasons for those interactions not being in the summary information is well understood. There are many interactions missing from the summary information I provide above which are not rank deficient. For instance, why are there no interactions between the x12LD level and any levels in the x9 variable despite many of those interactions being possible?

@Roland detected the right reason for this behavior:

That's a direct consequence of the default treatment contrasts.

I share his Opinon and found an excellent Answer to a simular Question on Stackoverflow which can be read here.

---

While dealing with rank deficient you can look:

• at first at your data and check for missing values with e.g. df[!complete.cases(df), ], and perhaps remove them an example of this Problem can be found here
  note that the term na.action = na.exclude in your models should be na.action = na.omit since this is the default command of lmer, you can read about it here
• however your data has no incomplet cases so this doese not cause the warning
  fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients

Checking your model summary e.g. summary(lm3) and anova(lm3) the anova gives you the reason for dropping 2 columns/coefficients:

> anova(lm3)
Missing cells for: x12LD:x9BY, x12LD:x9PW.  
Interpret type III hypotheses with care.
Type III Analysis of Variance Table with Satterthwaite's method
        Sum Sq Mean Sq NumDF DenDF   F value    Pr(>F)    
x1     1907.61 1907.61     1 73643 9367.3611 < 2.2e-16 ***
x2      588.91  588.91     1 89254 2891.8584 < 2.2e-16 ***
x3     1732.84 1732.84     1 95991 8509.1166 < 2.2e-16 ***
x4       26.37   26.37     1  4554  129.4769 < 2.2e-16 ***
x5        0.15    0.15     1   466    0.7380   0.39075    
x12:x6    1.66    0.83     2   567    4.0723   0.01754 *  
x12:x7   53.71   26.86     2 41303  131.8788 < 2.2e-16 ***
x12:x8   13.56    6.78     2  4899   33.2939 4.345e-15 ***
x12:x9  118.36    8.45    14 37055   41.5141 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Your dataset is missing cases with an interaction between x12LD:x9BY and x12LD:x9PW this can be checked and shown with:

df[df$x12 == "LD" & df$x9 == "BY", ] #checking for existing cases -> none
 [1] X   x1  x2  x3  x4  x5  x12 x6  x7  x8  x9  x10 x11 y  
<0 Rows> (or row.names with Length 0)

df[df$x12 == "LD" & df$x9 == "FB", ] #checking for existing cases -> plenty
        X    x1   x2  x3       x4       x5 x12    x6   x7      x8 x9
73    125 19.60 61.2  89 23.00000  9031.35  LD  5.96 1181  862.35 FB
84    137 26.89 38.0  93  0.00670  8742.60  LD  1.55 1248 1004.33 FB
126   179 37.70 52.7  98  3.00000 10412.22  LD 16.56  824  995.43 FB
[ reached 'max' / getOption("max.print") -- omitted 3298 rows ]

Since there is no data avaiable for these interactions it is clear that lmer can not estimate any coefficients (this is also the case with your other models lm4, lm5)

It is greatly explained in this thread: What is rank deficiency, and how to deal with it?

Rank deficiency in this context says there is insufficient information contained in your data to estimate the model you desire. [...] The deficiency may stem from simply too little data. In general, you cannot uniquely estimate n parameters with less than n data points.

---

Two possible solutions can be found in this Post from @Ben Bolker:

The warning (not an error) occurs because you don't have all combinations of timing, intensity, and year in your data set, but you are asking R to estimate parameters for every combination. A few reasonable choices are:

• ignore the warning (you'll probably get reasonable answers anyway when comparing the overall effect of each factor)
• reduce the complexity of the model

---

The code I used:

#load data
df <- read.table("df.txt", 
                 sep = ",", 
                 header = T)

#check data
str(df)
summary(df)
head(df)
length(df) == length(df[complete.cases(df), ]) #checking if there are any "NA's" in your dataframe

#specifing models
library(lmerTest)
lm1 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + (1|x10/x11), data=df, na.action=na.exclude)
lm2 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x12 + (1|x10/x11), data=df, na.action=na.exclude)
lm3 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12:x6 + x12:x7 + x12:x8 + x12:x9 + (1|x10/x11), data=df, na.action=na.omit)
lm4 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x7 + x12*x8 + x12*x9 + (1|x10/x11), data=df, na.action=na.omit)
lm5 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12:x6 + x12:x9:x7 + x12:x9:x8 + (1|x10/x11), data=df, na.action=na.exclude)
lm6 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x9*x7 + x12*x9*x8 + (1|x10/x11), data=df, na.action=na.exclude)

#checking models
summary(lm3)
anova(lm3) #Note the Warning "Missing cells for: x12LD:x9BY, x12LD:x9PW."

df[df$x12 == "LD" & df$x9 == "BY", ] #checking for existing cases -> none
df[df$x12 == "LD" & df$x9 == "FB", ] #checking for existing cases -> plenty