matrices de covarianza singulares = estimaciones de la varianza cero de los efectos aleatorios en un modelo con interceptos aleatorios

Question

Me preguntaba por qué para los modelos de efectos mixtos que tienen matrices de covarianza degeneradas o singulares (se estima que algunas combinaciones lineales de los efectos aleatorios no tienen variabilidad) se dice que corresponden:

Estimaciones de la varianza de efectos aleatorios cero en un modelo con interceptos aleatorios solamente

Para concretar esta discusión, utilicemos el modelo ajustado a continuación (en 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
```

Answer 1

No estoy seguro de dónde está el allFit de la función, pero su código produce errores en mi sistema.

De todos modos, no estoy muy de acuerdo con la premisa de que un modelo con un ajuste singular se diga que corresponde a estimaciones de varianza de efectos aleatorios cero en un modelo con interceptos aleatorios solamente.

Esa podría ser una razón, pero en mi experiencia es más frecuente que los modelos singulares se produzcan cuando los efectos aleatorios están sobreajustados -muy a menudo demasiadas pendientes aleatorias- en relación con los datos. Otra razón común es cuando hay muy poca variación en las pendientes dentro de algún factor de agrupación, de manera que no se necesitan pendientes aleatorias para la variable en cuestión. Es fácil construir un ejemplo sencillo:

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

que resulta en un ajuste singular (correlación estimada en 1) donde el verdadero problema es que hay muy poca varianza entre las pendientes aleatorias.

En mi opinión, la razón por la que vemos esto tan a menudo se debe al infame artículo de Barr et al (2013) "Keep it Maximal", donde el consejo es básicamente ajustar pendientes aleatorias, con correlaciones, para todos los efectos fijos. Incluso uno de los modelos de su artículo tiene un ajuste singular, del que no se dieron cuenta.