Encolhimento suspeito de regressão logística de efeitos aleatórios

8

Considere o seguinte exemplo simples:

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 )

Ou seja, estamos fornecendo efeitos fixos e um modelo de efeitos aleatórios para o mesmo problema muito simples (regressão logística, somente interceptação).

É assim que o modelo de efeitos fixos se parece:

plot( summary( fitFE ) )

Modelo de efeitos fixos

E é assim que os efeitos aleatórios:

dotplot( ranef( fitRE, condVar = TRUE ) )

Modelo de efeitos aleatórios

O encolhimento não é surpreendente, mas sua extensão é. Aqui está uma comparação mais direta:

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

Probabilidades previstas a partir de modelos de efeitos fixos e aleatórios

As estimativas de efeitos fixos variam de menos de 3% a mais de 11, mas os efeitos aleatórios estão entre 7,5 e 9,5%. (A inclusão de covariáveis ​​torna isso ainda mais extremo.)

Não sou especialista em efeitos aleatórios em regressão logística, mas, a partir da regressão linear, tive a impressão de que uma retração substancial pode ocorrer apenas em grupos muito pequenos. Aqui, no entanto, até o menor grupo tem quase uma centena de observações e o tamanho da amostra ultrapassa 500.

Qual é a razão para isto? Ou estou com vista para algo ...?

EDIT (28 de julho de 2017). Conforme a sugestão de @Ben Bolker, tentei o que acontece se a resposta for contínua (para removermos problemas sobre o tamanho efetivo da amostra, específico para dados binomiais).

O novo SimDataé, portanto,

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

e os novos modelos são

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

O resultado com

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

é

insira a descrição da imagem aqui

Por enquanto, tudo bem!

No entanto, decidi fazer outra verificação para verificar a ideia de Ben, mas o resultado acabou sendo bastante bizarro. Decidi verificar a teoria de outra maneira: volto ao resultado binário, mas aumento os meios para que os tamanhos efetivos das amostras aumentem. Simplesmente corri params$means <- params$means + 0.5e tentei novamente o exemplo original, eis o resultado:

insira a descrição da imagem aqui

Apesar do tamanho mínimo (efetivo) da amostra aumentar drasticamente ...

> 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 

... o encolhimento realmente aumentou ! (Tornando-se total, com variação zero estimada.)

Tamas Ferenci
fonte
Você está plotando o odds ratio no primeiro gráfico e o log odds ratio no segundo gráfico.
Douglas Bates
Sim, mas o terceiro gráfico, que realmente os compara, e mostra o problema dessa pergunta, usa a mesma escala para ambos! Assim como meu comando verbal abaixo da trama.
Tamas Ferenci

Respostas:

5

n=4n=104

Tamanhos efetivos de amostra 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 

Tamanhos de amostra 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 

Uma maneira de testar essa explicação seria fazer um exemplo análogo com respostas continuamente variáveis ​​(gama ou gaussianas).

Ben Bolker
fonte
2
uau, tamanho efetivo da amostra, eu nunca teria pensado nisso. Obrigado! Meu experimento com a resposta gaussiana confirma sua ideia, mas aumentar o tamanho mínimo efetivo da amostra não; veja minha edição ...
Tamas Ferenci