Estimates of random effects in binomial model (lme4)

Question

I'm simulating Bernoulli trials with a random $\text{logit}\, \theta \sim {\cal N}(\text{logit}\, \theta_0, 1^2)$ between groups and then I fit the corresponding model with the lme4 package:

library(lme4)
library(data.table)
I <- 30 # number of groups
J <- 10 # number of Bernoulli trials within each group
logit <- function(p) log(p)-log(1-p)
expit <- function(x) exp(x)/(1+exp(x))
theta0 <- 0.7
ddd <- data.table(subject=factor(1:I),logittheta=rnorm(I, logit(theta0)))[, list(result=rbinom(J, 1, expit(logittheta))), by=subject]
fit <- glmer(result~(1|subject), data=ddd, family="binomial")
props <- ddd[, list(p=mean(result)), by=subject]$p
estims <- expit(coef(fit)$subject[,1])
par(pty="s")
plot(props, estims, asp=1, xlim=c(0,1), ylim=c(0,1),
     xlab="proportion", ylab="glmer estimate")
abline(0,1)

Then I compare the proportions of successes by group with their estimates and I always get such a result:

By "always", I mean that the glmer estimates are always higher than the empirical proportions for small proportions, and always lower for high proportions. The glmer estimate is close to the empirical proportion for values around the overall proportion ($0.7$ in my example). After increasing $J$ the difference between estimates and proportions becomes negligible but one always get this picture. Is it a known fact and why does it hold ? I expected to get estimates centered around the empirical proportions.

Answer 1

What you're seeing is a phenomenon called shrinkage, which is a fundamental property of mixed models; individual group estimates are "shrunk" toward the overall mean as a function of the relative variance of each estimate. (While shrinkage is discussed in various answers on CrossValidated, most refer to techniques such as lasso or ridge regression; answers to this question provide connections between mixed models and other views of shrinkage.)

Shrinkage is arguably desirable; it is sometimes referred to as borrowing strength. Especially when we have few samples per group, the separate estimates for each group are going to be less precise than estimates that take advantage of some pooling from each population. In a Bayesian or empirical Bayesian framework, we can think of the population-level distribution as acting as a prior for the group-level estimates. Shrinkage estimates are especially useful/powerful when (as is not the case in this example) the amount of information per group (sample size/precision) varies widely, e.g. in a spatial epidemiological model where there are regions with very small and very large populations.

The shrinkage property should apply to both Bayesian and frequentist fitting approaches -- the real differences between the approaches lie at the top level (the frequentist's "penalized weighted residual sum of squares" is the Bayesian's log-posterior deviance at the group level ...) The main difference in the picture below, which shows lme4 and MCMCglmm results, is that because MCMCglmm uses a stochastic algorithm the estimates for different groups with the same observed proportions differ slightly.

With a little more work, I think we could figure out the precise degree of shrinkage expected by comparing the binomial variances for the groups and the overall data set, but in the meantime here's a demonstration (the fact that the J=10 case looks less shrunk than J=20 is just sampling variation, I think). (I accidentally changed the simulation parameters to mean=0.5, RE standard deviation=0.7 (on logit scale) ...)

library("lme4")
library("MCMCglmm")
##' @param I number of groups
##' @param J number of Bernoulli trials within each group
##' @param theta random effects standard deviation (logit scale)
##' @param beta intercept (logit scale)
simfun <- function(I=30,J=10,theta=0.7,beta=0,seed=NULL) {
    if (!is.null(seed)) set.seed(seed)
    ddd <- expand.grid(subject=factor(1:I),rep=1:J)
    ddd <- transform(ddd,
                     result=suppressMessages(simulate(~1+(1|subject),
                     family=binomial,
                     newdata=ddd,
                     newparams=list(theta=theta,beta=beta))[[1]]))
}
sumfun <- function(ddd) {
    fit <- glmer(result~(1|subject), data=ddd, family="binomial")
    fit2 <- MCMCglmm(result~1,random=~subject, data=ddd,
                    family="categorical",verbose=FALSE,
                    pr=TRUE)
    res <- data.frame(
        props=with(ddd,tapply(result,list(subject),mean)),
        lme4=plogis(coef(fit)$subject[,1]),
        MCMCglmm=plogis(colMeans(fit2$Sol[,-1])))
    return(res)
}
set.seed(101)
res <- do.call(rbind,
        lapply(c(10,20,50,100,500),
               function(J) {
                   data.frame(J=J,sumfun(simfun(J=J)))
               }))
library("reshape2")
m <- melt(res,id.vars=c("J","props"))
library("ggplot2"); theme_set(theme_bw())
ggplot(m,aes(props,value))+
    geom_point(aes(colour=factor(J),shape=variable))+
    geom_abline(intercept=0,slope=1,colour="gray")+
      labs(x="observed proportion",y="estimate")
ggsave("shrinkage.png",width=5,height=5)