Dlaczego nie mogę dopasować wyjścia glmer (rodzina = dwumianowy) do ręcznej implementacji algorytmu Gaussa-Newtona?

Chciałbym dopasować wyjścia lmera (naprawdę blasku) do dwumianowego przykładu zabawki. Przeczytałem winiety i wierzę, że rozumiem, co się dzieje.

Ale najwyraźniej nie. Po utknięciu, naprawiłem „prawdę” w kategoriach efektów losowych i poszedłem po ocenie samych ustalonych efektów. Podaję ten kod poniżej. Aby zobaczyć, że jest to uzasadnione, możesz skomentować + Z %*% b.ki będzie pasowało do wyników zwykłego glm. Mam nadzieję, że pożyczę trochę mocy mózgu, aby dowiedzieć się, dlaczego nie jestem w stanie dopasować mocy lmera, jeśli uwzględnione zostaną efekty losowe.

# 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")
}

Jeśli zmienisz polecenie dopasowania modelu na następujące, podejście dopasowywania działa:

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

Kluczową zmianą jest ta nAGQ = 0, która odpowiada twojemu podejściu, podczas gdy default ( nAGQ = 1) nie. nAGQoznacza „liczbę adaptacyjnych punktów kwadraturowych Gaussa-Hermity” i określa, w jaki sposób glmerzintegrować efekty losowe podczas dopasowywania modelu mieszanego. Kiedy nAGQjest większa niż 1, wówczas do nAGQpunktów używana jest kwadratura adaptacyjna . Kiedy nAGQ = 1używane jest przybliżenie Laplace'a, a kiedy nAGQ = 0całka jest „ignorowana”. Nie będąc zbyt szczegółowym (a zatem być może zbyt technicznym), nAGQ = 0oznacza, że ​​efekty losowe wpływają tylko na oszacowania efektów stałych poprzez ich szacowane tryby warunkowe - dlatego teżnAGQ = 0nie uwzględnia w pełni losowości efektów losowych. Aby w pełni uwzględnić efekty losowe, należy je zintegrować. Jednak, jak odkryłeś, różnica między nAGQ = 0i nAGQ = 1często może być dość mała.

Twoje podejście do dopasowywania nie będzie działać nAGQ > 0. Wynika to z faktu, że w tych przypadkach optymalizacja składa się z trzech etapów: (1) ukarana iteracyjnie ponownie ważona metodą najmniejszych kwadratów (PIRLS) w celu oszacowania trybów warunkowych efektów losowych, (2) (w przybliżeniu) zintegrowanie efektów losowych dotyczących ich trybów warunkowych oraz (3) nieliniowa optymalizacja funkcji celu (tj. wynik całkowania). Kroki te są powtarzane aż do konwergencji. Po prostu wykonujesz iteracyjnie przeważoną serię najmniejszych kwadratów (IRLS), która zakłada, że bjest znana i wprowadza Z%*%btermin przesunięcia. Twoje podejście okazuje się być równoważne PIRLS, ale ta równoważność obowiązuje tylko dlatego, że używasz glmerszacunkowych trybów warunkowych (których inaczej byś nie wiedział).

Przepraszamy, jeśli nie jest to dobrze wyjaśnione, ale nie jest to temat, który nadaje się do szybkiego opisu. Może się przydać https://github.com/lme4/lme4pureR , który jest (niepełną) implementacją lme4podejścia w czystym kodzie R. lme4pureRjest zaprojektowany tak, aby był bardziej czytelny niż lme4sam (choć znacznie wolniej).