Estimations des effets aléatoires dans le modèle binomial (lme4)

Je simule des essais de Bernoulli avec un aléatoire entre les groupes, puis j'adapte le modèle correspondant avec le paquet:$\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)

Ensuite, je compare les proportions de succès par groupe avec leurs estimations et j'obtiens toujours un tel résultat:

Par " toujours ", je veux dire que les estimations glmer sont toujours plus élevées que les proportions empiriques pour les petites proportions, et toujours plus faibles pour les proportions élevées. L'estimation glmer est proche de la proportion empirique pour les valeurs autour de la proportion globale ( dans mon exemple). Après avoir augmenté la différence entre les estimations et les proportions devient négligeable, mais on obtient toujours cette image. Est-ce un fait connu et pourquoi tient-il? Je m'attendais à obtenir des estimations centrées sur les proportions empiriques.$0.7$$J$

Ce que vous voyez est un phénomène appelé rétrécissement , qui est une propriété fondamentale des modèles mixtes; les estimations de chaque groupe sont «rétrécies» vers la moyenne globale en fonction de la variance relative de chaque estimation. (Bien que le retrait soit discuté dans diverses réponses sur CrossValidated, la plupart se réfèrent à des techniques telles que le lasso ou la régression de crête; les réponses à cette question fournissent des connexions entre des modèles mixtes et d'autres vues du retrait.)

Le retrait est sans doute souhaitable; il est parfois appelé force d'emprunt . Surtout lorsque nous avons peu d'échantillons par groupe, les estimations distinctes pour chaque groupe vont être moins précises que les estimations qui tirent parti d'une mise en commun de chaque population. Dans un cadre bayésien ou empirique bayésien, nous pouvons penser que la distribution au niveau de la population agit comme un a priori pour les estimations au niveau du groupe. Les estimations du retrait sont particulièrement utiles / puissantes lorsque (comme ce n'est pas le cas dans cet exemple) la quantité d'informations par groupe (taille / précision de l'échantillon) varie considérablement, par exemple dans un modèle épidémiologique spatial où il y a des régions avec des populations très petites et très grandes. .

La propriété de retrait devrait s'appliquer aux approches d'ajustement bayésien et fréquentiste - les vraies différences entre les approches se situent au niveau supérieur (la "somme résiduelle pondérée des carrés pénalisée" du fréquentiste est la déviance log-postérieure du bayésien au niveau du groupe ... ) La principale différence dans l'image ci-dessous, qui montre lme4et MCMCglmmdonne des résultats, est que parce que MCMCglmm utilise un algorithme stochastique, les estimations pour différents groupes avec les mêmes proportions observées diffèrent légèrement.

Avec un peu plus de travail, je pense que nous pourrions déterminer le degré précis de retrait attendu en comparant les variances binomiales pour les groupes et l'ensemble de données global, mais en attendant, voici une démonstration (le fait que le cas J = 10 semble moins rétréci que J = 20 est juste une variation d'échantillonnage, je pense). (J'ai accidentellement changé les paramètres de simulation pour signifier = 0,5, écart type RE = 0,7 (sur l'échelle 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)