[R-sig-ME] (lme4) Calculating variance component and random-effect shrinkage in a logistic random-intercept model

[Juho Kristian Ruohonen]

Hi,

I'm working with a simple, pedagogical dataset, trying to teach myself how
the random-effect SD is estimated and how the individual random effects are
shrunk towards zero.

The dataset is based on Agresti 2018 (*Introduction to Categorical Data
Analysis, 3rd ed*. pp. 278-281), but I've modified it to get rid of all 0%
and 100% proportions, since those would yield infinite logits and
complicate the calculus. There are just 20 binomial observations,
representing clusters with sample sizes varying from 5 to 20. And there are
only two parameters: a fixed intercept and a random intercept:



*set.seed(2019)d <- data.frame(id = letters[1:20], TrueLogit = rnorm(20), y
= NA, n = round(runif(20, min = 5, max = 20)))*
*d$y <- round(plogis(d$TrueLogit) * d$n) *#Create cluster-specific success
counts based on the "true" logits

Each cluster's sample proportion (on the logit scale) should be helpful in
estimating its random intercept:
*d$EmpiricalLogit <- qlogis(d$y / d$n) *

The fixed intercept should likewise be easy to calculate -- it is the
overall sample proportion on the logit scale, no? i.e:
*d$CommonLogit <- qlogis(sum(d$y) / sum(d$n)) *

The empirical, unshrunken "random effects" should be easy to calculate as
differences between the overall sample proportion and the cluster-specific
ones, all on the logit scale:
*d$Difference <- d$EmpiricalLogit - d$CommonLogit*

The empirical SD of the (as yet) unshrunken random effects should be easy
to calculate as the SD of the cluster-specific differences from the common
intercept, each weighted by the cluster size:
*sd(rep(d$Difference, times = d$n)) *#1.028238

This SD equals 1.03 logits. However, lme4 estimates it at only 0.77:
*require(lme4)*
*(mix <- glmer(y/n ~ (1|id) , weights = n, family = binomial, data = d))*

I thought shrinkage was only applied to the individual random effects, not
to their estimated variance. But the output of lme4 suggests that even the
estimated variance undergoes some shrinkage. How?

Secondly, Agresti (2018: 279), whose examples indeed use lme4, states that
in an example like ours, the fitted value of each cluster is "a weighted
average of the sample proportion for that area (cluster) alone and the
sample proportion after pooling all n (20) samples." (parentheticals mine).
If this is so, then at least the fitted probabilities should be easy to
calculate manually by using the aforementioned weighted-average approach:



*(d$fitted.manual <- sapply(1:nrow(d), function(i){  weighted.mean(c(d$y[i]
/ d$n[i], sum(d$y) / sum(d$n)), w = c(d$n[i], sum(d$n)))  }))*

Yet this misses the mark as well. The manually fitted values are way more
shrunken than the lme4-fitted values:


*d$fitted.lme4 <- fitted(mix)cbind(ManualFitted = d$fitted.manual,
SampleProp = d$y/d$n, lme4Fitted = d$fitted.lme4)*

Trying to calculate the random effects manually, as weighted averages
between 0 and the cluster-specific Difference from the fixed intercept,
fails equally. The manual calculation results in much more shrinkage than
that seen in the output of *ranef() *(not shown).

So, even though the principle behind the calculation of the variance
component and the individual intercepts in a GLMM sounds straightforward, I
cannot replicate even one step of the procedure, not even with a maximally
simple example dataset. This is demoralizing to someone who is trying to
gain a solid understanding of the method. Any pointers?

Best,

Juho

	[[alternative HTML version deleted]]