Why is this linear mixed model singular?

Question

I'm trying to understand why I get a singular fit when a linear mixed-effect model is fitted to the data below.

I used R lme4::lmer and the model is very simple having only the intercept as fixed effect and a factor variable as random.

Here's the dataset (can be copy & pasted to R)

data <- read.table(text= "
    group_id     y
           1  6.38
           1 10.83
           1 13.25
           1  2.96
           1 11.29
           1 11.52
           1  8.28
           1  8.36
           1  8.31
           1  7.33
           2  8.57
           2  7.00
           2  7.67
           2 10.19
           2 12.88
           2  9.67
           2  8.47
           2  7.27
           2  7.49
           2 17.25
           3 10.40
           3  8.53
           3  8.68
           3 11.38
           3  7.92
           3  5.66
           3 11.72
           3  6.93
           3  9.95
           3  7.19
           4 13.31
           4  8.57
           4  7.87
           4  8.50
           4  5.11
           4  6.50
           4  3.46
           4  5.98
           4  9.12
           4  8.60
           5 14.35
           5  6.79
           5  7.43
           5  9.16
           5  7.02
           5  7.09
           5  6.68
           5  6.24
           5  8.43
           5  8.51", 
    header= TRUE, colClasses= c('factor', 'numeric'))

This is the fitted model:

library(lme4)

fit <- lmer(data= data, y ~ 1 + (1|group_id))

boundary (singular) fit: see ?isSingular      <<<<<< 


summary(fit)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ 1 + (1 | group_id)
   Data: data

REML criterion at convergence: 239

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.139 -0.604 -0.093  0.467  3.242 

Random effects:
 Groups   Name        Variance Std.Dev.
 group_id (Intercept) 0.00     0.00    
 Residual             7.05     2.66    
Number of obs: 50, groups:  group_id, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)    8.641      0.376      23
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular

Help for isSingular says that some "dimensions" of the variance-covariance matrix have been estimated as exactly zero and I'd like to see in the data why this is happening.

Answer 1

As you have discovered, this happens when one of the variance components is estimated as zero. This typically has one of two explanations:

• the random effects structure is over-fitted - usually because of too many random slopes

• one of more variance components is actually very close to zero and there is insufficient data to estimate it above zero.

Obviously the first scenario is not the case with your data, since you only have random intercepts.

So it is likely that the actual variation of the random intercepts for group_id is very close to zero. If it is, then with only 5 groups, the software may not be able to estimate a variance above zero.

A good place to start if by plotting the data:

We can already see that the variation in the means of the groups is small compared to the variation within the groups.

We can investigate this more formally three ways (at least):

First, let us look at the means of the data in each group:

library(tidyverse)
data %>% group_by(group_id) %>% summarize(mean = mean(y))

## 1         8.85
## 2         9.65
## 3         8.84
## 4         7.70
## 5         8.17

Note that there is fairly small variation among all the groups, but note that the mean of group 1 and 3 is almost identical. Let us remove group 1 and see what happens:

data %>% filter(group_id != 1) %>% lmer(y ~ 1 + (1|group_id), data = .) %>% summary()

## Random effects:
##  Groups   Name        Variance Std.Dev.
##  group_id (Intercept) 0.03789  0.1947  
##  Residual             6.75636  2.5993  
## Number of obs: 40, groups:  group_id, 4

## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   8.5885     0.4224   20.34

So the model converges without singularity, but the variance component for group_id is very small, as we suspected.

Next, we can add some additional variance to the group_id component. The problem with doing this is that with only 5 groups, if we were to sample 5 observations from, say rnorm(5, 0, 1) (with a standard deviation of 1) the sample standard deviation is likely to be not close to 1 and the mean is likely to be not close to zero. A good approach to solve this is to use monte-carlo simulation (basically just do it many times and take averages).

Here we will do 100 simulations:

n.sim <- 100
simvec_rint <- numeric(n.sim)  # vector to hold the random intercepts variances
simvec_fint <- numeric(n.sim)  # vector to hold the fixed intercepts

for (i in 1:n.sim) {
  set.seed(i)
  data$y1 = data$y + rep(rnorm(5, 0, 1), each = 10)
  m0 <- lmer(y1 ~ 1 + (1|group_id), data = data)

  if (!isSingular(m0)) {
    # If the model is not singular then extract the random and fixed effects
    VarCorr(m0) %>% as.data.frame() %>% pull(vcov) %>% nth(1) -> simvec_rint[i]
    summary(m0) %>% coef() %>% as.vector() %>% nth(1) -> simvec_fint[i]
  } else {
    simvec_rint[i] <- simvec_fint[i] <- NA
  }
}

So we have added random noise to the groups with a variance of 1. The monte carlo estimates are:

> mean(simvec_rint, na.rm = TRUE)
[1] 1.116416
> mean(simvec_fint, na.rm = TRUE)
[1] 8.637063

Note that:

• The variance of the random intercepts is 1.12. However we have added variance to the groups equal to 1, so this implies that the variance of random intercepts in the original data is close to zero, as we suspected.

• The fixed intercept is 8.64 which is basically the same as the model fitted to the original data.

Lastly, let us look at a model without random effects, which will obviously simply be an ANOVA:

> aov(y ~ group_id, data = data) %>% summary()
            Df Sum Sq Mean Sq F value Pr(>F)
group_id     4   22.0   5.489   0.763  0.555
Residuals   45  323.6   7.190      

So there is very little evidence that the means of the 5 groups are different from each other. Another way to look at this is with:

> lm(y ~ group_id, data = data) %>% summary()

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.8510     0.8479  10.438 1.33e-13 ***
group_id2     0.7950     1.1992   0.663    0.511    
group_id3    -0.0150     1.1992  -0.013    0.990    
group_id4    -1.1490     1.1992  -0.958    0.343    
group_id5    -0.6810     1.1992  -0.568    0.573    

So there is also very little evidence that groups 2, 3, 4, and 5 have different means from group 1. Both these models are consistent with there being very small variation of the random intercepts in the mixed model.

So, to sum up, in this case we can conclude that due to a combination of the small number of groups and the estimated variation between groups being small, the software is unable to estimate the random intercepts variation above zero, and hence the model is singular, although the model estimates seem to be reliable.