# singular covariance matrices = estimates of zero random-effects variance in a model with random-intercepts

I was wondering why for mixed-effects models that have degenerate or singular covariance matrices (some linear combinations of the random effects are estimated to having no variability) are said to correspond to:

Estimates of zero random-effects variance in a model with random-intercepts only

To make this discussion concrete, please let's use the model fit below (in R).

library(lme4)

fit <- lmer(math ~ ses*sector + (ses|sch.id), data = hsb)

lme4::allFit(fit) # You'll see several boundary (singular) fit messages under different optimizations
$$$$


I'm not sure where the allFit function comes from, but your code produces errors on my system.

Anyway, I don't really agree with the premise that a model with a singular fit are said to correspond to estimates of zero random-effects variance in a model with random-intercepts only.

That could be one reason, but in my experience it is more often the case that singular models occur when the random effects are overfitted - very often too many random slopes - relative to the data. Another common reason is where there is very little variation in the slopes within some grouping factor, such that random slopes for the variable in question are not needed. It is easy to construct a simple example:

set.seed(15)
n.school <- 20

dt <- expand.grid(sch.id = 1:n.school, ses = 1:10)
dt$Y <- 1 X <- model.matrix( ~ ses, data = dt) myFormula <- "Y ~ ses + (ses | sch.id)" foo <- lFormula(eval(myFormula), dt) Z <- t(as.matrix(foo$$reTrms$$Zt)) betas <- c(1, 2) s1 <- 3 # SD of random intercepts s2 <- 0.1 # SD of random slopes rho <- 0.4 # correlation between intercepts and slopes cormat <- matrix(c(s1, rho, rho, s2), 2, 2) # correlation matrix covmat <- lme4::sdcor2cov(cormat) # covariance matrix (needed for mvrnorm) umat <- MASS::mvrnorm(n.school, c(0, 0), covmat, empirical = TRUE) # simulate the random effects u <- c(rbind(umat[, 1], umat[, 2])) # lme4 needs the random effects in this order interleaved) e <- rnorm(nrow(dt), 0, 2) # residual error dt$Y <- X %*% betas + Z %*% u + e

lmer(myFormula, dt) %>% summary()


which resuls in a singular fit (correlation estimated at 1) where the real problem is that there is too little variance among the random slopes.

In my opinion, the reason we see this so often is due to the infamous paper by Barr et al (2013) "Keep it Maximal" where the advice is basically to fit random slopes, with correlations, for all fixed effects. Even one of the models in their paper has a singular fit, which they didn't realise.

• Thanks Robert! The above quote was from Douglas Bates' et al where he criticizes keeping it maximal. Also, if we do lme4::rePCA(fit) does that reveal the problem better? – rnorouzian Sep 26 '20 at 15:11
• rePCA` does a principal components analysis on the VCV matrix which is a good way to identify where the problem in the random effects structure is. I've written an answer about this previously: stats.stackexchange.com/questions/449095/… – Robert Long Sep 26 '20 at 15:19
• Note that the quote you took from Bates' paper is missing the end of the sentence: "This corresponds to estimates of zero random-effects variance in a model with random-intercepts only or a correlation of ±1 in a model with correlated random intercepts and slopes." – Robert Long Sep 26 '20 at 15:43
• Ahh, well, it means that if you had very small variance of random intercepts, even with meaningful variation in random slopes, you could get a singular fit. If you change the seed to 5, s1 to 0.1 and s2 to 3 in the above code you will get a singular fit for that reason. – Robert Long Sep 26 '20 at 15:58
• Rob, I think this question is right up your alley. – rnorouzian Oct 7 '20 at 17:09