Contracción sospechosamente alta en regresión logística de efectos aleatorios

Considere el siguiente ejemplo sencillo:

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 )

Es decir, estamos dando efectos fijos y un modelo de efectos aleatorios para el mismo problema muy simple (regresión logística, intercepción solamente).

Así es como se ve el modelo de efectos fijos:

plot( summary( fitFE ) )

Y así es como los efectos aleatorios:

dotplot( ranef( fitRE, condVar = TRUE ) )

La contracción no es sorprendente en sí misma, pero sí lo es. Aquí hay una comparación más directa:

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

Las estimaciones de efectos fijos oscilan entre menos del 3% y más de 11, sin embargo, los efectos aleatorios están entre el 7,5 y el 9,5%. (La inclusión de covariables hace que esto sea aún más extremo).

No soy un experto en efectos aleatorios en la regresión logística, pero a partir de la regresión lineal, tenía la impresión de que una reducción tan sustancial solo puede ocurrir a partir de grupos muy pequeños. Aquí, sin embargo, incluso el grupo más pequeño tiene casi un centenar de observaciones, y los tamaños de muestra superan los 500.

¿Cuál es la razón para esto? ¿O estoy pasando por alto algo ...?

EDITAR (28 de julio de 2017). Según la sugerencia de @Ben Bolker, probé lo que sucede si la respuesta es continua (de modo que eliminemos los problemas sobre el tamaño de muestra efectivo , que es específico para los datos binomiales).

Lo nuevo SimDataes por lo tanto

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

y los nuevos modelos son

fitFE2 <- ols( Res2 ~ ID, data = SimData )
fitRE2 <- lmer( Res2 ~ ( 1|ID ), data = SimData )

El resultado con

xyplot( fe~re, data = data.frame( re = coef( fitRE2 )$ID[ , 1 ],
                       fe = c( 0, coef( fitFE2 )[ -1  ] )+coef( fitFE2 )[ 1 ] ),
    abline = c( 0, 1 ) )

es

¡Hasta aquí todo bien!

Sin embargo, decidí realizar otra verificación para verificar la idea de Ben, pero el resultado resultó ser bastante extraño. Decidí verificar la teoría de otra manera: vuelvo al resultado binario, pero aumento los medios para que los tamaños de muestra efectivos crezcan. Simplemente ejecuté params$means <- params$means + 0.5y luego volví a intentar el ejemplo original, aquí está el resultado:

A pesar del tamaño mínimo (efectivo) de la muestra de hecho drásticamente el aumento de ...

> 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 

... la contracción en realidad aumentó ! (Se convierte en total, con cero varianza estimada).

Sospecho que la respuesta aquí tiene que ver con la definición de "tamaño de muestra efectivo". Una regla general (del libro de Estrategias de modelado de regresión de Harrell ) es que el tamaño de muestra efectivo para una variable de Bernoulli es el mínimo del número de éxitos y fracasos; Por ejemplo, una muestra de tamaño 10,000 con solo 4 éxitos es más como tener$n=4$ que $n=10^4$. Los tamaños de muestra efectivos aquí no son pequeños, pero son mucho más pequeños que el número de observaciones.

Tamaños de muestra efectivos por grupo:

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 

Tamaños de muestra por grupo:

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 

Una forma de probar esta explicación sería hacer un ejemplo análogo con respuestas continuamente variables (gamma o gaussianas).