Alto índice sospechoso de contracción en regresión logística de efectos aleatorios

Question

Considere el siguiente ejemplo simple:

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 proporcionando un modelo de efectos fijos y un modelo de efectos aleatorios para el mismo problema, muy simple (regresión logística, solo intercepción).

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

plot( summary( fitFE ) )

Y así es como se ve el modelo de efectos aleatorios:

dotplot( ranef( fitRE, condVar = TRUE ) )

La contracción no es sorprendente en sí misma, pero su extensión 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 van desde menos del 3% hasta más del 11, sin embargo, los efectos aleatorios están entre 7.5 y 9.5%. (La inclusión de covariables hace que esto sea aún más extremo.)

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

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

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

El nuevo SimData es 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 comprobación para verificar la idea de Ben, pero el resultado resultó ser bastante extraño. Decidí verificar la teoría de otra manera: volví al resultado binario, pero aumenté las medias para que los tamaños de muestra efectivos fueran más grandes. Simplemente ejecuté params$means <- params$means + 0.5 y luego volví a intentar el ejemplo original, aquí está el resultado:

A pesar de que el tamaño de muestra (efectivo) mínimo aumentó drásticamente...

> 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ó! (Convirtiéndose en total, con una varianza estimada de cero.)

Answer 1

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

Tamaños efectivos de muestra 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 de variación continua (Gamma o Gaussiana).