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")
}
Best Answer
If you change your model fitting command to the following, your matching approach works:
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 howglmer
will integrate out the random effects when fitting the mixed model. WhennAGQ
is greater than 1, then adaptive quadrature is used withnAGQ
points. WhennAGQ = 1
, the Laplace approximation is used, and whennAGQ = 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 betweennAGQ = 0
andnAGQ = 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 assumesb
is known and puttingZ%*%b
in an offset term. Your approach turns out to be equivalent to PIRLS, but this equivalence only holds because you useglmer
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 thanlme4
itself (although much slower).