4
$\begingroup$

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)

hsb <- read.csv('https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')

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

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

1 Answer 1

3
$\begingroup$

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 were not aware of (see the paper by Bates 2015 Parsimoneous Mixed Models for more details).

$\endgroup$
10
  • $\begingroup$ 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? $\endgroup$
    – rnorouzian
    Sep 26, 2020 at 15:11
  • $\begingroup$ 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/… $\endgroup$ Sep 26, 2020 at 15:19
  • $\begingroup$ 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." $\endgroup$ Sep 26, 2020 at 15:43
  • $\begingroup$ 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. $\endgroup$ Sep 26, 2020 at 15:58
  • 1
    $\begingroup$ Rob, I think this question is right up your alley. $\endgroup$
    – rnorouzian
    Oct 7, 2020 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.