2

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
```
Robert Long
  • 53,316
  • 10
  • 84
  • 148
rnorouzian
  • 3,056
  • 2
  • 16
  • 40

1 Answers1

2

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).

Robert Long
  • 53,316
  • 10
  • 84
  • 148
  • 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: https://stats.stackexchange.com/questions/449095/how-to-simplify-a-singular-random-structure-when-reported-correlations-are-not-n/449122#449122 – 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
  • I see, lastly, the singularity could also occur if no correlation exists between the slopes and intercepts as well? Does `rePCA()` reveal whether that could be the case? – rnorouzian Sep 26 '20 at 16:06
  • No, I don't see why a correlation of zero would cause a singular fit. Go ahead and run the code above in a loop with different seeds, with s1=2 and s2=1 and rho=0. I bet you won't find many cases of singular fit – Robert Long Sep 26 '20 at 16:12
  • But won't the model be overparameterized and not converged when we ask for a correlation that in realty (or at least in the current dataset) doesn't exist? – rnorouzian Sep 26 '20 at 16:14
  • No, a correlation of zero is a perfectly valid correlation. – Robert Long Sep 26 '20 at 16:32
  • 1
    Rob, I think [this question](https://stats.stackexchange.com/questions/490891/specifying-the-random-part-of-a-mixed-model-in-r) is right up your alley. – rnorouzian Oct 07 '20 at 17:09
  • Will take a look soon – Robert Long Oct 07 '20 at 17:22