¿Por qué no puedo igualar la salida de glmer (familia=binomial) con la implementación manual del algoritmo Gauss-Newton?

Question

Me gustaría igualar los resultados de lmer (realmente glmer) con un ejemplo binomial de juguete. He leído los documentos y creo que entiendo lo que está sucediendo.

Pero aparentemente no es así. Después de quedarme atascado, arreglé la "verdad" en términos de los efectos aleatorios y fui tras la estimación de los efectos fijos solamente. Incluyo el código a continuación. Para ver que es legítimo, puedes comentar + Z %*% b.k y coincidirá con los resultados de un glm regular. Espero poder contar con ayuda para averiguar por qué no puedo igualar la salida de lmer cuando se incluyen los efectos aleatorios.

# Configuración - codificación manual de un conjunto de datos sencillo
df <- data.frame(x1 = rep(c(1:5), 3), subject = sort(rep(c(1:3), 5)))
df$subject <- factor(df$subject)

# Valores verdaderos de los coeficientes
beta <- matrix(c(-3.3, 1), ncol = 1) # Intercepto y pendiente, respectivamente 
u <- matrix(c(-.5, .6, .9), ncol = 1) # efectos aleatorios para los 3 sujetos 

# Matrices de diseño Z (efectos aleatorios) y X (efectos fijos)
Z <- model.matrix(~ 0 + factor(subject), data = df)
X <- model.matrix(~ 1 + x1, data = df)

# Respuesta  
df$y <- c(1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1)
y <- df$y

### Objetivo: igualar las estimaciones de la siguiente salida de lmer! 
library(lme4)
my.lmer <- lmer( y ~ x1 + (1 | subject), data = df, family = binomial)
summary(my.lmer)
ranef(my.lmer)

### Comienza el esfuerzo para igualar 

beta.k <- matrix(c(-3, 1.5), ncol = 1) # Valores iniciales (cercanos a la verdad)
b.k <- matrix(c(1.82478, -1.53618, -.5139356), ncol = 1) # efectos aleatorios de lmer

# Algoritmo iterativo de Gauss-Newton
for (iter in 1:6) {
  lin.pred <- as.numeric(X %*% beta.k +  Z %*% b.k)
  mu.k <- plogis(lin.pred)
  variances <- mu.k * (1 - mu.k)
  W.k <- diag(1/variances)

  y.star <- W.k^(.5) %*% (y - mu.k)
  X.star <- W.k^(.5) %*% (variances * X)
  delta.k <- solve(t(X.star) %*% X.star) %*% t(X.star) %*% y.star

  # Actualización de Gauss-Newton 
  beta.k <- beta.k + delta.k
  cat(iter, "Efectos Fijos: ", beta.k, "\n")
}

Answer 1

Si cambias tu comando de ajuste de modelo a lo siguiente, tu enfoque de emparejamiento funciona:

my.lmer <- glmer(y ~ x1 + (1 | subject), data = df, family = binomial, nAGQ = 0)

El cambio clave es el nAGQ = 0, que coincide con tu enfoque, mientras que el valor por defecto (nAGQ = 1) no lo hace. nAGQ significa 'número de puntos de cuadratura adaptativa de Gauss-Hermite', y establece cómo glmer integrará los efectos aleatorios al ajustar el modelo mixto. Cuando nAGQ es mayor que 1, se utiliza cuadratura adaptativa con nAGQ puntos. Cuando nAGQ = 1, se utiliza la aproximación de Laplace, y cuando nAGQ = 0, la integral se 'ignora'. Sin ser demasiado específico (y por lo tanto tal vez demasiado técnico), nAGQ = 0 significa que los efectos aleatorios solo influyen en las estimaciones de los efectos fijos a través de sus modos condicionales estimados; por lo tanto, nAGQ = 0 no tiene en cuenta completamente la aleatoriedad de los efectos aleatorios. Para tener en cuenta completamente los efectos aleatorios, estos deben integrarse. Sin embargo, como descubriste, esta diferencia entre nAGQ = 0 y nAGQ = 1 a menudo es bastante pequeña.

Tu enfoque de emparejamiento no funcionará con nAGQ > 0. Esto se debe a que en estos casos hay tres pasos para la optimización: (1) mínimos cuadrados reponderados iterativamente con penalización (PIRLS) para estimar los modos condicionales de los efectos aleatorios, (2) integrar (aproximadamente) los efectos aleatorios alrededor de sus modos condicionales, y (3) optimización no lineal de la función objetivo (es decir, el resultado de la integración). Estos pasos se iteran hasta la convergencia. Simplemente estás ejecutando un ajuste de mínimos cuadrados reponderados iterativamente (IRLS), que asume que b es conocido y colocando Z%*%b en un término de desplazamiento. Tu enfoque resulta ser equivalente a PIRLS, pero esta equivalencia solo se mantiene porque usas glmer para obtener los modos condicionales estimados (que de otra manera no conocerías).

Disculpas si esto no está bien explicado, pero no es un tema que se preste bien a una descripción rápida. Puede que encuentres útil https://github.com/lme4/lme4pureR, que es una implementación (incompleta) del enfoque de lme4 en código R puro. lme4pureR está diseñado para ser más legible que lme4 en sí mismo (aunque mucho más lento).