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