Rétrécissement étrangement élevé dans la régression logistique à effets aléatoires

Prenons l'exemple simple suivant:

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 )

C'est-à-dire que nous donnons un effet fixe et un modèle à effets aléatoires pour le même problème très simple (régression logistique, interception uniquement).

Voici à quoi ressemble le modèle à effets fixes:

plot( summary( fitFE ) )

Et voici comment les effets aléatoires:

dotplot( ranef( fitRE, condVar = TRUE ) )

Le retrait n'est pas surprenant en soi, mais son ampleur l'est. Voici une comparaison plus directe:

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

Les estimations des effets fixes varient de moins de 3% à plus de 11, mais les effets aléatoires se situent entre 7,5 et 9,5%. (L'inclusion de covariables rend cela encore plus extrême.)

Je ne suis pas un expert des effets aléatoires dans la régression logistique, mais à partir de la régression linéaire, j'avais l'impression qu'un rétrécissement si important ne peut se produire qu'à partir de très très petites tailles de groupe. Ici, cependant, même le plus petit groupe compte près d'une centaine d'observations et la taille des échantillons dépasse 500.

Quelle est la raison pour ça? Ou est-ce que j'oublie quelque chose ...?

EDIT (28 juil.2017). Conformément à la suggestion de @Ben Bolker, j'ai essayé ce qui se passe si la réponse est continue (afin d'éliminer les problèmes de taille d'échantillon efficace , qui est spécifique aux données binomiales).

Le nouveau SimDataest donc

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

et les nouveaux modèles sont

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

Le résultat avec

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

est

Jusqu'ici tout va bien!

Cependant, j'ai décidé d'effectuer une autre vérification pour vérifier l'idée de Ben, mais le résultat s'est avéré assez bizarre. J'ai décidé de vérifier la théorie d'une autre manière: je reviens au résultat binaire, mais j'augmente les moyens pour que les tailles d'échantillon efficaces deviennent plus grandes. J'ai simplement couru params$means <- params$means + 0.5puis réessayé l'exemple d'origine, voici le résultat:

Malgré la taille minimale (effective) de l'échantillon, elle augmente en effet considérablement ...

> 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 

... le retrait a en fait augmenté ! (Devenir total, avec une variance nulle estimée.)

Je soupçonne que la réponse ici a à voir avec la définition de "taille d'échantillon effective". Une règle d'or (tirée du livre Harrell's Regression Modeling Strategies ) est que la taille réelle de l'échantillon d'une variable de Bernoulli est le minimum du nombre de réussites et d'échecs; Par exemple, un échantillon de taille 10 000 avec seulement 4 succès ressemble plus à qu'à . Les tailles d'échantillon efficaces ici ne sont pas minuscules, mais elles sont beaucoup plus petites que le nombre d'observations.$n=4$$n=10^4$

Tailles d'échantillon efficaces par groupe:

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 

Taille des échantillons par groupe:

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 

Une façon de tester cette explication serait de faire un exemple analogue avec des réponses variant en continu (gamma ou gaussiennes).