Estimativas de efeitos aleatórios no modelo binomial (lme4)

Estou simulando ensaios de Bernoulli com um aleatório entre grupos e, em seguida, encaixo o modelo correspondente em o pacote:$\text{logit}\, \theta \sim {\cal N}(\text{logit}\, \theta_0, 1^2)$lme4

library(lme4)
library(data.table)
I <- 30 # number of groups
J <- 10 # number of Bernoulli trials within each group
logit <- function(p) log(p)-log(1-p)
expit <- function(x) exp(x)/(1+exp(x))
theta0 <- 0.7
ddd <- data.table(subject=factor(1:I),logittheta=rnorm(I, logit(theta0)))[, list(result=rbinom(J, 1, expit(logittheta))), by=subject]
fit <- glmer(result~(1|subject), data=ddd, family="binomial")
props <- ddd[, list(p=mean(result)), by=subject]$p
estims <- expit(coef(fit)$subject[,1])
par(pty="s")
plot(props, estims, asp=1, xlim=c(0,1), ylim=c(0,1),
     xlab="proportion", ylab="glmer estimate")
abline(0,1)

Depois, comparo as proporções de sucessos por grupo com suas estimativas e sempre obtenho esse resultado:

Por " sempre ", quero dizer que as estimativas do glmer são sempre mais altas que as proporções empíricas para pequenas proporções e sempre mais baixas para altas proporções. A estimativa de glmer está próxima da proporção empírica para valores em torno da proporção geral ( no meu exemplo). Depois de aumentar a diferença entre estimativas e proporções se torna insignificante, mas sempre se obtém essa imagem. É um fato conhecido e por que isso ocorre? Eu esperava obter estimativas centradas em torno das proporções empíricas.$0.7$$J$

O que você está vendo é um fenômeno chamado encolhimento , que é uma propriedade fundamental dos modelos mistos; as estimativas de grupos individuais são "reduzidas" em relação à média geral em função da variação relativa de cada estimativa. (Embora o encolhimento seja discutido em várias respostas no CrossValidated, a maioria se refere a técnicas como regressão de laço ou crista; as respostas a essa pergunta fornecem conexões entre modelos mistos e outras visualizações de encolhimento.)

O encolhimento é indiscutivelmente desejável; às vezes é chamado de força de empréstimo . Especialmente quando temos poucas amostras por grupo, as estimativas separadas para cada grupo serão menos precisas do que as estimativas que tiram vantagem de alguns agrupamentos de cada população. Em uma estrutura bayesiana ou bayesiana empírica, podemos pensar na distribuição no nível da população como agindo como um prior para as estimativas no nível do grupo. As estimativas de retração são especialmente úteis / poderosas quando (como não é o caso neste exemplo) a quantidade de informações por grupo (tamanho / precisão da amostra) varia muito, por exemplo, em um modelo epidemiológico espacial em que existem regiões com populações muito pequenas e muito grandes .

A propriedade de encolhimento deve se aplicar às abordagens de ajuste bayesiana e freqüentista - as diferenças reais entre as abordagens estão no nível superior (a "soma residual quadrada ponderada penalizada" do frequentista é o desvio log-posterior do bayesiano no nível do grupo ... ) A principal diferença na figura abaixo, que mostra lme4e MCMCglmmresulta, é que, como o MCMCglmm usa um algoritmo estocástico, as estimativas para diferentes grupos com as mesmas proporções observadas diferem ligeiramente.

Com um pouco mais de trabalho, acho que poderíamos descobrir o grau exato de encolhimento esperado comparando as variações binomiais dos grupos e o conjunto de dados geral, mas, entretanto, aqui está uma demonstração (o fato de o caso J = 10 parecer menos encolhido que J = 20 é apenas variação de amostragem, eu acho). (Alterei acidentalmente os parâmetros da simulação para média = 0,5, desvio padrão de ER = 0,7 (na escala logit) ...)

library("lme4")
library("MCMCglmm")
##' @param I number of groups
##' @param J number of Bernoulli trials within each group
##' @param theta random effects standard deviation (logit scale)
##' @param beta intercept (logit scale)
simfun <- function(I=30,J=10,theta=0.7,beta=0,seed=NULL) {
    if (!is.null(seed)) set.seed(seed)
    ddd <- expand.grid(subject=factor(1:I),rep=1:J)
    ddd <- transform(ddd,
                     result=suppressMessages(simulate(~1+(1|subject),
                     family=binomial,
                     newdata=ddd,
                     newparams=list(theta=theta,beta=beta))[[1]]))
}
sumfun <- function(ddd) {
    fit <- glmer(result~(1|subject), data=ddd, family="binomial")
    fit2 <- MCMCglmm(result~1,random=~subject, data=ddd,
                    family="categorical",verbose=FALSE,
                    pr=TRUE)
    res <- data.frame(
        props=with(ddd,tapply(result,list(subject),mean)),
        lme4=plogis(coef(fit)$subject[,1]),
        MCMCglmm=plogis(colMeans(fit2$Sol[,-1])))
    return(res)
}
set.seed(101)
res <- do.call(rbind,
        lapply(c(10,20,50,100,500),
               function(J) {
                   data.frame(J=J,sumfun(simfun(J=J)))
               }))
library("reshape2")
m <- melt(res,id.vars=c("J","props"))
library("ggplot2"); theme_set(theme_bw())
ggplot(m,aes(props,value))+
    geom_point(aes(colour=factor(J),shape=variable))+
    geom_abline(intercept=0,slope=1,colour="gray")+
      labs(x="observed proportion",y="estimate")
ggsave("shrinkage.png",width=5,height=5)