Singularity issue in lmer (mixed model)

Question

I was trying to test state-level variation in a mixed model to predict attitudes using lmer() function.

My original code was

summary(out3 <- lmer(attitude ~age+ gender+race+education.z+
                       PID.z+c.ex+l.ex+
                       a.z+b.z+
                       cases.z+state.bi+
                       (1|state), long))

all variables are individual-level variables, except cases.z and state.bi. However, this model causes singularity issue-- state-level variance is zero and warning message.

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.7255 -0.6600 -0.0543  0.6332  4.8199 

Random effects:
 Groups   Name        Variance Std.Dev.
 state    (Intercept) 0.0000   0.0000  
 Residual             0.5719   0.7563  
Number of obs: 8280, groups:  state, 51

Fixed effects:
                  Estimate   Std. Error           df t value             Pr(>|t|)    
(Intercept)      2.7188699    0.0418414 8268.9999998  64.980 < 0.0000000000000002 ***
age              0.0005328    0.0005156 8268.9999998   1.033              0.30149    
race            -0.0426514    0.0199646 8268.9999998  -2.136              0.03268 *  
education.z     -0.0257635    0.0086102 8268.9999998  -2.992              0.00278 ** 
PID.z            0.4058194    0.0138666 8268.9999998  29.266 < 0.0000000000000002 ***
c.ex             0.0520329    0.0067580 8268.9999998   7.699   0.0000000000000152 ***
l.ex            -0.0413710    0.0040821 8268.9999999 -10.135 < 0.0000000000000002 ***
a.z             -0.3638447    0.0127821 8268.9999998 -28.465 < 0.0000000000000002 ***
b.z              0.5385977    0.0134168 8268.9999998  40.143 < 0.0000000000000002 ***
cases.z          0.0147368    0.0089858 8268.9999997   1.640              0.10104    
state.bi        -0.0052052    0.0179407 8268.9999998  -0.290              0.77172    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

> isSingular(out3)
[1] TRUE

So I tried a number of different sets of models, by removing variables.

Then I found when I remove a.z and b.z in the model, the singularity issue seems to be gone.

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3838 -0.6590 -0.0610  0.6936  4.1528 

Random effects:
 Groups   Name        Variance Std.Dev.
 state    (Intercept) 0.001854 0.04305 
 Residual             0.804934 0.89718 
Number of obs: 8280, groups:  state, 51

Fixed effects:
                  Estimate   Std. Error           df t value             Pr(>|t|)    
(Intercept)      2.7217213    0.0504252  927.3686109  53.975 < 0.0000000000000002 ***
age             -0.0015578    0.0005987 8269.3682573  -2.602              0.00929 ** 
race             0.0358835    0.0238024 7210.4125421   1.508              0.13171    
education.z     -0.0810985    0.0101691 8220.5507442  -7.975  0.00000000000000173 ***
PID.z            0.9792668    0.0115351 8246.0744340  84.894 < 0.0000000000000002 ***
c.ex             0.1366403    0.0078149 8255.4395441  17.485 < 0.0000000000000002 ***
l.ex            -0.0867261    0.0047532 8094.0071038 -18.246 < 0.0000000000000002 ***
cases.z          0.0295554    0.0129040   33.9892491   2.290              0.02833 *  
state.bi        -0.0112124    0.0264556   27.1321735  -0.424              0.67504    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I wonder I should remove a.z and b.z in the model. But I don't want to do it, since these variables are fundamental for my hypothesis. Do you have any solutions?

Also, my second question is -- why do they drop gender?

Answer 1

To awnser your last question first: lme4::lmer reports "fixed-effect model matrix is rank deficient", do I need a fix and how to?

And the second question : When you obtain a singular fit, this is often indicating that the model is overfitted – that is, the random effects structure is too complex to be supported by the data, which naturally leads to the advice to remove the most complex part of the random effects structure (usually random slopes). The benefit of this approach is that it leads to a more parsimonious model that is not over-fitted.

See: https://stats.stackexchange.com/questions/378939/dealing-with-singular-fit-in-mixed-models

Let me know if these links helped :)