Mengapa saya tidak dapat mencocokkan output glmer (keluarga = binomial) dengan penerapan algoritma Gauss-Newton secara manual?

Saya ingin mencocokkan output dari lmer (benar-benar glmer) dengan contoh binomial mainan. Saya sudah membaca sketsa dan yakin saya mengerti apa yang sedang terjadi.

Tapi ternyata saya tidak. Setelah macet, saya memperbaiki "kebenaran" dalam hal efek acak dan pergi setelah memperkirakan efek tetap saja. Saya memasukkan kode ini di bawah. Untuk melihat bahwa itu sah, Anda dapat berkomentar + Z %*% b.kdan itu akan cocok dengan hasil GLM biasa. Saya berharap untuk meminjam kekuatan otak untuk mencari tahu mengapa saya tidak dapat mencocokkan output lmer ketika efek acak dimasukkan.

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

Jika Anda mengubah perintah pemasangan model Anda ke yang berikut, pendekatan pencocokan Anda berfungsi:

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

Perubahan kuncinya adalah nAGQ = 0, yang cocok dengan pendekatan Anda, sedangkan default ( nAGQ = 1) tidak. nAGQberarti 'jumlah titik quadrature Gauss-Hermite adaptif', dan menetapkan bagaimana glmerakan mengintegrasikan efek acak ketika menyesuaikan model campuran. Ketika nAGQlebih besar dari 1, maka kuadratur adaptif digunakan dengan nAGQpoin. Kapan nAGQ = 1, pendekatan Laplace digunakan, dan kapan nAGQ = 0, integralnya 'diabaikan'. Tanpa terlalu spesifik (dan karena itu mungkin terlalu teknis), nAGQ = 0berarti bahwa efek acak hanya memengaruhi estimasi efek tetap melalui estimasi mode kondisional mereka - oleh karena itu,nAGQ = 0tidak sepenuhnya menjelaskan keacakan dari efek acak. Untuk sepenuhnya memperhitungkan efek acak, mereka perlu diintegrasikan. Namun, ketika Anda menemukan perbedaan ini nAGQ = 0dan nAGQ = 1sering kali bisa sangat kecil.

Pendekatan pencocokan Anda tidak akan bekerja dengan nAGQ > 0. Ini karena dalam kasus-kasus ini ada tiga langkah untuk optimasi: (1) dihukum kuadrat terkecil berulang berulang (PIRLS) untuk memperkirakan mode kondisional dari efek acak, (2) (kurang-lebih) mengintegrasikan efek acak tentang mode kondisional mereka , dan (3) optimasi nonlinier dari fungsi tujuan (yaitu hasil integrasi). Langkah-langkah ini sendiri diulang sampai konvergensi. Anda hanya melakukan menjalankan kuadrat berulang berulang berulang (IRLS), yang mengasumsikan bdiketahui dan memasukkan Z%*%bistilah offset. Pendekatan Anda ternyata setara dengan PIRLS, tetapi kesetaraan ini hanya berlaku karena Anda gunakan glmeruntuk mendapatkan estimasi mode bersyarat (yang Anda tidak akan tahu).

Permintaan maaf jika ini tidak dijelaskan dengan baik, tetapi itu bukan topik yang cocok dengan deskripsi cepat. Anda mungkin menemukan https://github.com/lme4/lme4pureR berguna, yang merupakan implementasi (tidak lengkap) dari lme4pendekatan dalam kode R murni. lme4pureRdirancang agar lebih mudah dibaca daripada lme4dirinya sendiri (meskipun jauh lebih lambat).