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).
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, x, 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 ) )
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.)