Почему я не могу сопоставить вывод glmer (family = binomial) с ручной реализацией алгоритма Гаусса-Ньютона?

Я хотел бы сравнить выходные данные lmer (действительно glmer) с примером игрушечного бинома. Я прочитал виньетки и, кажется, понимаю, что происходит.

Но, видимо, я не. Застряв, я исправил «правду» в терминах случайных эффектов и пошел оценивать только фиксированные эффекты. Я включаю этот код ниже. Чтобы увидеть, что это законно, вы можете закомментировать, + Z %*% b.kи это будет соответствовать результатам обычного GLM. Я надеюсь позаимствовать некоторые умственные способности, чтобы понять, почему я не могу сравниться с выходом Лмера, когда включены случайные эффекты.

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

Если вы измените команду подбора модели на следующую, ваш подход подбора будет работать:

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

Ключевое изменение - это то nAGQ = 0, что соответствует вашему подходу, тогда как default ( nAGQ = 1) - нет. nAGQозначает «число адаптивных квадратурных точек Гаусса-Эрмита» и задает, как glmerбудут интегрироваться случайные эффекты при подборе смешанной модели. Когда nAGQбольше 1, то адаптивная квадратура используется с nAGQточками. Когда nAGQ = 1используется приближение Лапласа, а когда nAGQ = 0интеграл «игнорируется». Не будучи слишком конкретным (и, следовательно, возможно, слишком техническим), nAGQ = 0означает, что случайные эффекты влияют только на оценки фиксированных эффектов через их оцененные условные режимы - поэтомуnAGQ = 0не полностью учитывает случайность случайных эффектов. Чтобы полностью учесть случайные эффекты, их необходимо интегрировать. Тем не менее, как вы обнаружили, эта разница nAGQ = 0и nAGQ = 1часто может быть довольно небольшой.

Ваш подход соответствия не будет работать с nAGQ > 0. Это происходит потому, что в этих случаях есть три шага к оптимизации: (1) штрафовать итеративно взвешенные наименьшие квадраты (PIRLS) для оценки условных мод случайных эффектов, (2) (приблизительно) интегрировать случайные эффекты относительно их условных мод и (3) нелинейная оптимизация целевой функции (т.е. результат интегрирования). Эти шаги сами повторяются до сходимости. Вы просто выполняете повторяющееся повторение взвешенных наименьших квадратов (IRLS), в котором предполагается, что bон известен, и вводите Z%*%bсмещение. Ваш подход оказывается эквивалентным PIRLS, но эта эквивалентность имеет место только потому, что вы используете glmerдля получения оценочных условных режимов (которые вы не знали бы иначе).

Извинения, если это не очень хорошо объяснено, но это не тема, которая поддается краткому описанию. Вам может пригодиться https://github.com/lme4/lme4pureR , который является (неполной) реализацией lme4подхода в чистом R-коде. lme4pureRразработан, чтобы быть более читабельным, чем он lme4сам (хотя и намного медленнее).