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

**Answer**

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 =`

, the Laplace approximation is used, and when

1`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).

**Attribution***Source : Link , Question Author : Ben Ogorek , Answer Author : Steve Walker*