# Suspiciously high shrinkage in random effects logistic regression

Consider the following simple example:

library( rms )
library( lme4 )
params <- structure(list(Ns = c(181L, 191L, 147L, 190L, 243L, 164L, 83L,
383L, 134L, 238L, 528L, 288L, 214L, 502L, 307L, 302L, 199L, 156L,
183L), means = c(0.09, 0.05, 0.03, 0.06, 0.07, 0.07, 0.1, 0.1,
0.06, 0.11, 0.1, 0.11, 0.07, 0.11, 0.1, 0.09, 0.1, 0.09, 0.08
)), .Names = c("Ns", "means"), row.names = c(NA, -19L), class = "data.frame")
SimData <- data.frame( ID = as.factor( rep( 1:nrow( params ), params$Ns ) ), Res = do.call( c, apply( params, 1, function( x ) c( rep( 0, x[ 1 ]-round( x[ 1 ]*x[ 2 ] ) ), rep( 1, round( x[ 1 ]*x[ 2 ] ) ) ) ) ) ) tapply( SimData$Res, SimData$ID, mean ) dd <- datadist( SimData ) options( datadist = "dd" ) fitFE <- lrm( Res ~ ID, data = SimData ) fitRE <- glmer( Res ~ ( 1|ID ), data = SimData, family = binomial( link = logit ), nAGQ = 50 )  I.e. we are giving a fixed effects and a random effects model for the the same, very simple problem (logistic regression, intercept only). This is how the fixed effects model looks like: plot( summary( fitFE ) )  And this is how random effects: dotplot( ranef( fitRE, condVar = TRUE ) )  The shrinkage is not surprising itself, but its extent is. Here is a more direct comparison: xyplot( plogis(fe)~plogis(re), data = data.frame( re = coef( fitRE )$ID[ , 1 ],
fe = c( 0, coef( fitFE )[ -1  ] )+coef( fitFE )[ 1 ] ),
abline = c( 0, 1 ) )


The fixed effects estimates range from less then 3% to more than 11, the random effects, however, are between 7.5 and 9.5%. (Inclusion of covariates makes this even more extreme.)

I'm no expert in random effects in logistic regression, but from linear regression, I was under the impression that so substantial shrinkage can occur only from very-very small group sizes. Here, however, even the smallest group has almost a hundred observation, and sample sizes go above 500.

What is the reason for this? Or am I overlooking something...?

EDIT (Jul 28, 2017). As per @Ben Bolker's suggestion, I tried what happens if the response is continuous (so that we remove problems about effective sample size, which is specific for binomial data).

The new SimData is therefore

SimData <- data.frame( ID = as.factor( rep( 1:nrow( params ), params$Ns ) ), Res = do.call( c, apply( params, 1, function( x ) c( rep( 0, x[ 1 ]-round( x[ 1 ]*x[ 2 ] ) ), rep( 1, round( x[ 1 ]*x[ 2 ] ) ) ) ) ), Res2 = do.call( c, apply( params, 1, function( x ) rnorm( x[1], x[2], 0.1 ) ) ) ) data.frame( params, Res = tapply( SimData$Res, SimData$ID, mean ), Res2 = tapply( SimData$Res2, SimData$ID, mean ) )  and the new models are fitFE2 <- ols( Res2 ~ ID, data = SimData ) fitRE2 <- lmer( Res2 ~ ( 1|ID ), data = SimData )  The result with xyplot( fe~re, data = data.frame( re = coef( fitRE2 )$ID[ , 1 ],
fe = c( 0, coef( fitFE2 )[ -1  ] )+coef( fitFE2 )[ 1 ] ),
abline = c( 0, 1 ) )


is

So far so good!

However, I decided to perform another check to verify Ben's idea, but the outcome turned out to be pretty bizarre. I decided to check the theory another way: I return to binary outcome, but increase the means so that effective sample sizes get bigger. I simply ran params$means <- params$means + 0.5 and then retried the original example, here is the result:

Despite the minimum (effective) sample size indeed drastically increasing...

> summary(with(SimData,tapply(Res,list(ID),
+                             function(x) min(sum(x==0),sum(x==1)))))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
33.0    72.5    86.0   100.3   117.5   211.0


...the shrinkage actually increased! (Becoming total, with zero variance estimated.)

• You are plotting the odds ratio in the first plot and the log odds ratio in the second plot. – Douglas Bates Jul 25 '17 at 16:29
• Yes, but the third plot, which actually compares them, and shows the problem of this question, uses the same scale for both! Just as my verbal command below the plot. – Tamas Ferenci Jul 25 '17 at 16:47

I suspect that the answer here has to do with the definition of "effective sample size". A rule of thumb (from Harrell's Regression Modeling Strategies book) is that effective sample size for a Bernoulli variable is the minimum of the number of successes and failures; e.g. a sample of size 10,000 with only 4 successes is more like having $n=4$ than $n=10^4$. The effective sample sizes here are not tiny, but they're a lot smaller than the number of observations.

Effective sample sizes per group:

summary(with(SimData,tapply(Res,list(ID),
function(x) min(sum(x==0),sum(x==1)))))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.00   11.00   16.00   21.63   29.00   55.00


Sample sizes per group:

summary(c(table(SimData\$ID)))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
83.0   172.5   199.0   243.8   295.0   528.0


One way to test this explanation would be to do an analogous example with continuously varying (Gamma or Gaussian) responses.

• wow, effective sample size, I'd have never thought of this. Thanks! My experiment with Gaussian response confirms your idea, but increasing the minimum effective sample size does not; see my edit... – Tamas Ferenci Jul 28 '17 at 17:06