Schätzungen zufälliger Effekte im Binomialmodell (lme4)

Ich simuliere Bernoulli-Versuche mit einem zufälligen zwischen Gruppen und passe dann das entsprechende Modell mit an das Paket:$\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)

Dann vergleiche ich die Anteile der Erfolge nach Gruppen mit ihren Schätzungen und erhalte immer ein solches Ergebnis:

Mit " immer " meine ich, dass die glmer-Schätzungen immer höher sind als die empirischen Anteile für kleine Anteile und immer niedriger für hohe Anteile. Die Glmer-Schätzung liegt nahe am empirischen Anteil für Werte um den Gesamtanteil ( in meinem Beispiel). Nach dem Erhöhen von der Unterschied zwischen Schätzungen und Anteilen vernachlässigbar, aber man bekommt immer dieses Bild. Ist es eine bekannte Tatsache und warum gilt es? Ich erwartete Schätzungen, die sich auf die empirischen Proportionen konzentrierten.$0.7$$J$

Was Sie sehen, ist ein Phänomen namens Schrumpfung , das eine grundlegende Eigenschaft gemischter Modelle ist. Einzelne Gruppenschätzungen werden in Abhängigkeit von der relativen Varianz jeder Schätzung in Richtung des Gesamtmittelwerts "geschrumpft". (Während das Schrumpfen in verschiedenen Antworten auf CrossValidated diskutiert wird, beziehen sich die meisten auf Techniken wie Lasso- oder Ridge-Regression. Antworten auf diese Frage stellen Verbindungen zwischen gemischten Modellen und anderen Ansichten des Schrumpfens her.)

Schrumpfen ist wohl wünschenswert; es wird manchmal als Kreditstärke bezeichnet . Insbesondere wenn wir nur wenige Stichproben pro Gruppe haben, werden die separaten Schätzungen für jede Gruppe weniger genau sein als Schätzungen, die eine gewisse Zusammenfassung aus jeder Population nutzen. In einem Bayes'schen oder empirischen Bayes'schen Rahmen können wir uns die Verteilung auf Bevölkerungsebene als Prior für die Schätzungen auf Gruppenebene vorstellen. Schrumpfungsschätzungen sind besonders nützlich / leistungsfähig, wenn (wie in diesem Beispiel nicht der Fall ist) die Informationsmenge pro Gruppe (Stichprobengröße / -genauigkeit) stark variiert, z. B. in einem räumlichen epidemiologischen Modell, in dem es Regionen mit sehr kleinen und sehr großen Populationen gibt .

Die Schrumpfeigenschaft sollte sowohl für Bayes'sche als auch für frequentistische Anpassungsansätze gelten - die tatsächlichen Unterschiede zwischen den Ansätzen liegen auf der obersten Ebene (die "bestrafte gewichtete Restquadratsumme" des Frequentisten ist die log-posteriore Abweichung des Bayesian auf Gruppenebene ... ) Der Hauptunterschied im Bild unten, der zeigt lme4und MCMCglmmresultiert, besteht darin, dass sich die Schätzungen für verschiedene Gruppen mit den gleichen beobachteten Anteilen geringfügig unterscheiden, da MCMCglmm einen stochastischen Algorithmus verwendet.

Mit etwas mehr Arbeit könnten wir den genauen Grad der erwarteten Schrumpfung herausfinden, indem wir die Binomialvarianzen für die Gruppen und den Gesamtdatensatz vergleichen. In der Zwischenzeit ist hier eine Demonstration (die Tatsache, dass der Fall J = 10 weniger aussieht geschrumpft als J = 20 ist nur eine Stichprobenvariation, denke ich). (Ich habe versehentlich die Simulationsparameter auf Mittelwert = 0,5, RE-Standardabweichung = 0,7 (auf der Logit-Skala) geändert ...)

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)