singular fit issue for simple mixed model

Question

I'm trying to run a simple MLM, but I'm bumping into singular fit warnings:

    y <- Reduce(c,replicate(10, scale(rnorm(700, 0, 1))))
    x <- rep(c("A","B"), each=700, times=10)
    g <- rep(c("g1", "g2", "g3", "g4", "g5", "g6", "g7", "g8", 
              "g9", "g10"), each=7000)
    
    df <- data.frame(y=y, x=x, g=g)
    
    m <- lmer(y ~ x + (1|g), data=df)
    boundary (singular) fit: see ?isSingular

I believe the reason for this warning is that there seems to be no variation of the random intercept as ranef(m) shows.

I read several answers here in CV suggesting to reduce the complexity of the model, often by removing random slopes. However, this model is already as simple as it can get, and there are no random slopes. I also have over 200 groups in my actual data set so it cannot be that I have too few groups.

What are my options? Should I remove the random intercept entirely?

Answer 1

Your example model samples fine, if a little slowly (~45 minutes on my laptop), using rstanarm. So one solution is to use a Bayesian MLM/LMM:

library(rstanarm)
sm <- stan_lmer(y ~ x + (1|g), data=df,
            chains= 4,
            cores = 4)

which recovers the expected values:

                                   mean   sd   10%   50%   90%
(Intercept)                       0.0    0.0  0.0   0.0   0.0 
xB                                0.0    0.0 -0.1   0.0   0.0 
b[(Intercept) g:g1]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g10]              0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g2]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g3]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g4]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g5]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g6]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g7]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g8]               0.0    0.0  0.0   0.0   0.0 
b[(Intercept) g:g9]               0.0    0.0  0.0   0.0   0.0 
sigma                             1.0    0.0  1.0   1.0   1.0 

This is just using the default priors, which you should tune to your specific question. My experience has been that mildly regularizing priors allow these otherwise singular models in lmer() to sample adequately.