Verdächtig hohe Schrumpfung der logistischen Regression mit zufälligen Effekten

8

Betrachten Sie das folgende einfache Beispiel:

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 )

Das heißt, wir geben ein festes Effekt- und ein Zufallseffektmodell für dasselbe, sehr einfache Problem an (logistische Regression, nur Abfangen).

So sieht das Modell mit festen Effekten aus:

plot( summary( fitFE ) )

Modell mit festen Effekten

Und so wirken sich zufällige Effekte aus:

dotplot( ranef( fitRE, condVar = TRUE ) )

Modell mit zufälligen Effekten

Das Schrumpfen ist an sich nicht überraschend, aber sein Ausmaß ist es. Hier ist ein direkterer Vergleich:

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

Vorausgesagte Wahrscheinlichkeiten aus Modellen mit festen und zufälligen Effekten

Die Schätzungen für feste Effekte reichen von weniger als 3% bis mehr als 11%, die zufälligen Effekte liegen jedoch zwischen 7,5 und 9,5%. (Die Einbeziehung von Kovariaten macht dies noch extremer.)

Ich bin kein Experte für zufällige Effekte in der logistischen Regression, aber aufgrund der linearen Regression hatte ich den Eindruck, dass eine so erhebliche Schrumpfung nur bei sehr, sehr kleinen Gruppengrößen auftreten kann. Hier hat jedoch selbst die kleinste Gruppe fast hundert Beobachtungen, und die Stichprobengrößen liegen über 500.

Was ist der Grund dafür? Oder übersehe ich etwas ...?

BEARBEITEN (28. Juli 2017). Gemäß dem Vorschlag von @Ben Bolker habe ich versucht, was passiert, wenn die Antwort kontinuierlich ist (damit wir Probleme mit der effektiven Stichprobengröße beseitigen , die für Binomialdaten spezifisch ist).

Das Neue SimDataist also

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

und die neuen Modelle sind

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

Das Ergebnis mit

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

ist

Geben Sie hier die Bildbeschreibung ein

So weit, ist es gut!

Ich beschloss jedoch, eine weitere Überprüfung durchzuführen, um Bens Idee zu überprüfen, aber das Ergebnis erwies sich als ziemlich bizarr. Ich beschloss, die Theorie auf eine andere Weise zu überprüfen: Ich kehre zum binären Ergebnis zurück, erhöhe aber die Mittelwerte, damit die effektiven Stichproben größer werden. Ich bin einfach gelaufen params$means <- params$means + 0.5und habe dann das ursprüngliche Beispiel erneut versucht. Hier ist das Ergebnis:

Geben Sie hier die Bildbeschreibung ein

Trotz der minimalen (effektiven) Stichprobengröße steigt die ...

> 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 

... der Schrumpf hat tatsächlich zugenommen ! (Gesamt werden, wobei die Varianz Null geschätzt wird.)

Tamas Ferenci
quelle
Sie zeichnen das Quotenverhältnis im ersten Diagramm und das logarithmische Quotenverhältnis im zweiten Diagramm.
Douglas Bates
Ja, aber die dritte Darstellung, die sie tatsächlich vergleicht und das Problem dieser Frage zeigt, verwendet für beide den gleichen Maßstab! Genau wie mein verbaler Befehl unter der Handlung.
Tamas Ferenci

Antworten:

5

Ich vermute, dass die Antwort hier mit der Definition der "effektiven Stichprobengröße" zu tun hat. Eine Faustregel (aus Harrells Regression Modeling Strategies- Buch) lautet, dass die effektive Stichprobengröße für eine Bernoulli-Variable das Minimum der Anzahl von Erfolgen und Misserfolgen ist. Beispiel: Eine Stichprobe der Größe 10.000 mit nur 4 Erfolgen entspricht eher als . Die effektiven Stichprobengrößen sind hier nicht winzig, aber viel kleiner als die Anzahl der Beobachtungen.n = 10 4n=4n=104

Effektive Stichprobengrößen pro Gruppe:

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 

Stichprobengröße pro Gruppe:

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 

Eine Möglichkeit, diese Erklärung zu testen, wäre ein analoges Beispiel mit kontinuierlich variierenden (Gamma- oder Gaußschen) Antworten.

Ben Bolker
quelle
2
Wow, effektive Stichprobengröße, daran hätte ich nie gedacht. Vielen Dank! Mein Experiment mit der Gaußschen Antwort bestätigt Ihre Idee, eine Erhöhung der minimalen effektiven Stichprobengröße jedoch nicht. siehe meine Bearbeitung ...
Tamas Ferenci