Why can't I match glmer (family=binomial) output with manual implementation of Gauss-Newton algorithm? I'd like to match the outputs of lmer (really glmer) with a toy binomial example. I've read the vignettes and believe I understand what's going on.
But apparently I do not. After getting stuck, I fixed the "truth" in terms of the random effects and went after estimation of the fixed effects alone. I'm including this code below. To see that it's legit, you can comment out + Z %*% b.k and it will match the results of a regular glm. I'm hoping to borrow some brainpower to figure out why I'm not able to match lmer's output when the random effects are included.
# 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")
}

 A: If you change your model fitting command to the following, your matching
approach works:
my.lmer <- glmer(y ~ x1 + (1 | subject), data = df, family = binomial, nAGQ = 0)

The key change is the nAGQ = 0, which matches your approach, whereas the
default (nAGQ = 1) does not.  nAGQ means 'number of adaptive
Gauss-Hermite quadrature points', and sets how glmer will integrate out
the random effects when fitting the mixed model.  When nAGQ is greater
than 1, then adaptive quadrature is used with nAGQ points.  When nAGQ =
1, the Laplace approximation is used, and when nAGQ = 0, the integral is
'ignored'.  Without being too specific (and therefore perhaps too
technical), nAGQ = 0 means that the random effects only influence the
estimates of the fixed effects through their estimated conditional modes
-- therefore, nAGQ = 0 does not completely account for the randomness of
the random effects.  To fully account for the random effects, they need
to be integrated out.  However, as you discovered this difference
between nAGQ = 0 and nAGQ = 1 can often be fairly small.
Your matching approach will not work with nAGQ > 0.  This is because in
these cases there are three steps to the optimization:  (1) penalized
iteratively reweighted least squares (PIRLS) to estimate the conditional
modes of the random effects, (2) (approximately) integrate out the
random effects about their conditional modes, and (3) nonlinear
optimization of the objective function (i.e. the result of the
integration).  These steps are themselves iterated until convergence.
You are simply doing an iteratively reweighted least squares (IRLS) run,
which assumes b is known and putting Z%*%b in an offset term.  Your
approach turns out to be equivalent to PIRLS, but this equivalence only
holds because you use glmer to get estimated conditional modes (which
you wouldn't otherwise know).
Apologies if this isn't well explained, but it isn't a topic that lends
itself well to a quick description. You might find
https://github.com/lme4/lme4pureR useful, which is an (incomplete)
implementation of the lme4 approach in pure R code.  lme4pureR is
designed to be more readable than lme4 itself (although much slower).
