Por que não consigo combinar a saída glmer (family = binomial) com a implementação manual do algoritmo de Gauss-Newton?

Gostaria de combinar as saídas do lmer (realmente glmer) com um exemplo binomial de brinquedo. Eu li as vinhetas e acredito que entendo o que está acontecendo.

Mas aparentemente eu não. Depois de ficar preso, consertei a "verdade" em termos dos efeitos aleatórios e fui atrás da estimativa dos efeitos fixos. Estou incluindo este código abaixo. Para ver que é legítimo, você pode comentar + Z %*% b.ke ele corresponderá aos resultados de um glm regular. Espero emprestar algum poder intelectual para descobrir por que não sou capaz de igualar a produção de Imer quando os efeitos aleatórios são incluídos.

# Setup - hard coding simple data set 
df <- data.frame(x1 = rep(c(1:5), 3), subject = sort(rep(c(1:3), 5)))
df$subject <- factor(df$subject)

# True coefficient values  
beta <- matrix(c(-3.3, 1), ncol = 1) # Intercept and slope, respectively 
u <- matrix(c(-.5, .6, .9), ncol = 1) # random effects for the 3 subjects 

# Design matrices Z (random effects) and X (fixed effects)
Z <- model.matrix(~ 0 + factor(subject), data = df)
X <- model.matrix(~ 1 + x1, data = df)

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

### Goal: match estimates from the following lmer output! 
library(lme4)
my.lmer <- lmer( y ~ x1 + (1 | subject), data = df, family = binomial)
summary(my.lmer)
ranef(my.lmer)

### Matching effort STARTS HERE 

beta.k <- matrix(c(-3, 1.5), ncol = 1) # Initial values (close to truth)
b.k <- matrix(c(1.82478, -1.53618, -.5139356), ncol = 1) # lmer's random effects

# Iterative Gauss-Newton algorithm
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

  # Gauss-Newton Update 
  beta.k <- beta.k + delta.k
  cat(iter, "Fixed Effects: ", beta.k, "\n")
}

Se você alterar o comando de ajuste de modelo para o seguinte, sua abordagem de correspondência funcionará:

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

A mudança de chave é o nAGQ = 0, que corresponde à sua abordagem, enquanto o padrão ( nAGQ = 1) não. nAGQsignifica "número de pontos de quadratura adaptativos de Gauss-Hermite" e define comoglmer integrará os efeitos aleatórios ao ajustar o modelo misto. Quando nAGQé maior que 1, a quadratura adaptativa é usada com nAGQpontos. Quando nAGQ = 1, a aproximação de Laplace é usada e quando nAGQ = 0, a integral é 'ignorada'. Sem ser muito específico (e, portanto, talvez muito técnico), nAGQ = 0significa que os efeitos aleatórios influenciam apenas as estimativas dos efeitos fixos por meio de seus modos condicionais estimados - portanto,nAGQ = 0não explica completamente a aleatoriedade dos efeitos aleatórios. Para explicar completamente os efeitos aleatórios, eles precisam ser integrados. No entanto, como você descobriu essa diferença entre nAGQ = 0e nAGQ = 1geralmente pode ser bastante pequena.

Sua abordagem de correspondência não funcionará nAGQ > 0. Isso ocorre porque nesses casos há três etapas para a otimização: (1) mínimos quadrados com ponderação iterativa e reponderada (PIRLS) para estimar os modos condicionais dos efeitos aleatórios, (2) (aproximadamente) integrar os efeitos aleatórios sobre seus modos condicionais e (3) otimização não linear da função objetivo (ou seja, o resultado da integração). Essas etapas são iteradas até a convergência. Você está simplesmente executando uma execução IRLS (IQM) de mínimos quadrados iterativamenteb é conhecida e está colocando Z%*%bum termo de compensação. Sua abordagem acaba sendo equivalente a PIRLS, mas essa equivalência é válida apenas porque você usa glmerpara obter modos condicionais estimados (que você não conheceria de outra forma).

Desculpas se isso não for bem explicado, mas não é um tópico que se presta bem a uma descrição rápida. Você pode achar https://github.com/lme4/lme4pureR útil, que é uma implementação (incompleta) da lme4abordagem em código R puro. lme4pureRfoi projetado para ser mais legível que lme4ele (embora muito mais lento).