Warum unterscheiden sich die geschätzten Werte eines Best Linear Unbias Predictor (BLUP) von einem Best Linear Unbias Estimator (BLUE)?

20

Ich verstehe, dass der Unterschied zwischen ihnen damit zusammenhängt, ob die Gruppierungsvariable im Modell als fester oder zufälliger Effekt geschätzt wird, aber mir ist nicht klar, warum sie nicht gleich sind (wenn sie nicht gleich sind).

Ich bin speziell daran interessiert, wie dies bei der Verwendung der Schätzung kleiner Flächen funktioniert, wenn dies relevant ist, aber ich vermute, dass die Frage für jede Anwendung fester und zufälliger Effekte relevant ist.

Jeremy Miles
quelle

Antworten:

26

Die Werte, die Sie von BLUPs erhalten, werden nicht auf die gleiche Weise geschätzt wie die BLUE-Schätzungen fester Effekte. Konventionell werden BLUPs als Vorhersagen bezeichnet . Wenn Sie ein Modell mit gemischten Effekten anpassen, werden zunächst der Mittelwert und die Varianz (und möglicherweise die Kovarianz) der zufälligen Effekte geschätzt. Der zufällige Effekt für eine bestimmte Lerneinheit (z. B. einen Schüler) wird anschließend aus dem geschätzten Mittelwert und der geschätzten Varianz sowie den Daten berechnet. In einem einfachen linearen Modell wird der Mittelwert geschätzt (ebenso wie die Restvarianz), aber die beobachteten Punktzahlen setzen sich sowohl aus diesem als auch dem Fehler zusammen, der eine Zufallsvariable ist. In einem Modell mit gemischten Effekten ist der Effekt für eine bestimmte Einheit ebenfalls eine Zufallsvariable (obwohl er in gewissem Sinne bereits realisiert wurde).

Sie können solche Einheiten auch als feste Effekte behandeln, wenn Sie möchten. In diesem Fall werden die Parameter für diese Einheit wie gewohnt geschätzt. In einem solchen Fall wird jedoch der Mittelwert (zum Beispiel) der Bevölkerung, aus der die Einheiten stammen, nicht geschätzt.

Die Annahme hinter zufälligen Effekten ist außerdem, dass sie zufällig aus einer bestimmten Population entnommen wurden, und es ist die Population, die Sie interessiert. Die Annahme, die festen Effekten zugrunde liegt, ist, dass Sie diese Einheiten gezielt ausgewählt haben, da dies die einzigen Einheiten sind, die Sie interessieren.

Wenn Sie sich umdrehen und ein Modell mit gemischten Effekten anpassen und dieselben Effekte vorhersagen, werden sie im Verhältnis zu ihren Schätzungen der festen Effekte tendenziell in Richtung des Bevölkerungsmittelwerts „geschrumpft“. Sie können sich dies analog zu einer Bayes'schen Analyse vorstellen, bei der der geschätzte Mittelwert und die geschätzte Varianz einen normalen Prior angeben und der BLUP dem Mittelwert des Seitenzahns entspricht, der aus der optimalen Kombination der Daten mit dem Prior resultiert.

Das Ausmaß der Schrumpfung hängt von mehreren Faktoren ab. Eine wichtige Bestimmung, wie weit die Vorhersagen für zufällige Effekte von den Schätzungen für feste Effekte entfernt sein werden, ist das Verhältnis der Varianz der zufälligen Effekte zur Fehlervarianz. Hier ist eine kurze RDemo für den einfachsten Fall mit 5 'Level 2' Einheiten mit nur Mittelwerten (Abschnitten). (Sie können sich dies als Testergebnisse für Schüler in Klassen vorstellen.)

library(lme4)   # we'll need to use this package
set.seed(1673)  # this makes the example exactly reproducible
nj = 5;    ni = 5;    g = as.factor(rep(c(1:nj), each=ni))

##### model 1
pop.mean = 16;    sigma.g = 1;    sigma.e = 5
r.eff1   = rnorm(nj,    mean=0, sd=sigma.g)
error    = rnorm(nj*ni, mean=0, sd=sigma.e)
y        = pop.mean + rep(r.eff1, each=ni) + error

re.mod1  = lmer(y~(1|g))
fe.mod1  = lm(y~0+g)
df1      = data.frame(fe1=coef(fe.mod1), re1=coef(re.mod1)$g)

##### model 2
pop.mean = 16;    sigma.g = 5;    sigma.e = 5
r.eff2   = rnorm(nj,    mean=0, sd=sigma.g)
error    = rnorm(nj*ni, mean=0, sd=sigma.e)
y        = pop.mean + rep(r.eff2, each=ni) + error

re.mod2  = lmer(y~(1|g))
fe.mod2  = lm(y~0+g)
df2      = data.frame(fe2=coef(fe.mod2), re2=coef(re.mod2)$g)

##### model 3
pop.mean = 16;    sigma.g = 5;    sigma.e = 1
r.eff3   = rnorm(nj,    mean=0, sd=sigma.g)
error    = rnorm(nj*ni, mean=0, sd=sigma.e)
y        = pop.mean + rep(r.eff3, each=ni) + error

re.mod3  = lmer(y~(1|g))
fe.mod3  = lm(y~0+g)
df3      = data.frame(fe3=coef(fe.mod3), re3=coef(re.mod3)$g)

Das Verhältnis der Varianz der Zufallseffekte zur Fehlervarianz beträgt also 1/5 für model 1, 5/5 für model 2und 5/1 für model 3. Beachten Sie, dass ich Level verwendet habe, um die Modelle mit festen Effekten zu codieren. Wir können nun untersuchen, wie die geschätzten festen Effekte und die vorhergesagten zufälligen Effekte für diese drei Szenarien verglichen werden.

df1
#         fe1     re1
# g1 17.88528 15.9897
# g2 18.38737 15.9897
# g3 14.85108 15.9897
# g4 14.92801 15.9897
# g5 13.89675 15.9897

df2
#          fe2      re2
# g1 10.979130 11.32997
# g2 13.002723 13.14321
# g3 26.118189 24.89537
# g4 12.109896 12.34319
# g5  9.561495 10.05969

df3
#         fe3      re3
# g1 13.08629 13.19965
# g2 16.36932 16.31164
# g3 17.60149 17.47962
# g4 15.51098 15.49802
# g5 13.74309 13.82224

Eine andere Möglichkeit, zufällige Effektvorhersagen zu erhalten, die näher an den Schätzungen für feste Effekte liegen, besteht darin, dass Sie über mehr Daten verfügen. Wir können model 1von oben mit seinem geringen Verhältnis von Zufallseffektvarianz zu Fehlervarianz mit einer Version ( model 1b) mit demselben Verhältnis, aber viel mehr Daten vergleichen (beachten Sie, dass ni = 500anstelle von ni = 5).

##### model 1b
nj = 5;    ni = 500;    g = as.factor(rep(c(1:nj), each=ni))
pop.mean = 16;    sigma.g = 1;    sigma.e = 5
r.eff1b  = rnorm(nj,    mean=0, sd=sigma.g)
error    = rnorm(nj*ni, mean=0, sd=sigma.e)
y        = pop.mean + rep(r.eff1b, each=ni) + error

re.mod1b = lmer(y~(1|g))
fe.mod1b = lm(y~0+g)
df1b     = data.frame(fe1b=coef(fe.mod1b), re1b=coef(re.mod1b)$g)

Hier sind die Effekte:

df1
#         fe1     re1
# g1 17.88528 15.9897
# g2 18.38737 15.9897
# g3 14.85108 15.9897
# g4 14.92801 15.9897
# g5 13.89675 15.9897

df1b
#        fe1b     re1b
# g1 15.29064 15.29543
# g2 14.05557 14.08403
# g3 13.97053 14.00061
# g4 16.94697 16.92004
# g5 17.44085 17.40445

Aus einem ähnlichen Grund mag Doug Bates (der Autor des R-Pakets lme4) den Begriff "BLUP" nicht und verwendet stattdessen den "bedingten Modus" (siehe S. 22-23 seines Entwurfs lme4 Buch pdf ). Insbesondere weist er in Abschnitt 1.6 darauf hin, dass "BLUP" nur für lineare Mixed-Effects-Modelle sinnvoll verwendet werden kann .

gung - Wiedereinsetzung von Monica
quelle
3
+1. Ich bin mir jedoch nicht sicher, ob ich die terminologische Unterscheidung zwischen "Vorhersagen" und "Schätzen" richtig einschätze. Ein Verteilungsparameter ist also "geschätzt", aber eine latente Variable kann nur "vorhergesagt" werden? Verstehe ich dann richtig, dass zB Faktorladungen in der Faktoranalyse "geschätzt" werden, aber Faktorbewertungen "vorhergesagt" werden? Abgesehen davon finde ich es bemerkenswert verwirrend, dass etwas, das als "bester linearer unverzerrter Prädiktor" bezeichnet wird, tatsächlich ein voreingenommener Schätzer ist (weil es eine Schrumpfung implementiert und daher voreingenommen sein muss), wenn man es als einen "Schätzer" der festen Effekte behandeln würde. ..
Amöbe sagt Reinstate Monica
@amoeba, was bedeutet "am besten" überhaupt? Am besten was? Ist es die beste Schätzung des Datenmittelwerts oder die beste Kombination der in den Daten enthaltenen Informationen und der vorherigen? Hilft Ihnen die Bayes'sche Analogie?
gung - Wiedereinsetzung von Monica
2
Zumindest ist klar, was "linear" bedeutet :-) Im Ernst, ich habe diese sehr hilfreiche Antwort von @whuber auf den terminologischen Unterschied zwischen "Vorhersage" und "Schätzung" gefunden. Ich denke, es hat mir die Terminologie klargestellt, aber es hat mein Gefühl noch verstärkt, dass BLUP trotz seines Namens eher ein Schätzer ist. [Forts.]
Amöbe sagt Reinstate Monica
2
@amoeba, ja das ist alles vernünftig. Aber ich würde nicht den gleichen Namen für beide verwenden wollen, weil Sie etwas anderes tun (dh die Gleichungen unterscheiden sich) und es nützlich ist, wenn die Namen unterschiedlich sind.
gung - Wiedereinsetzung von Monica
1
@amoeba, ich habe die Formulierung im ersten Absatz angepasst, um diese Begriffe zu entwerten, um "Vorhersage" nicht zu verdrängen, sondern um die Unterscheidung beizubehalten. Sehen Sie, ob Sie glauben, ich hätte die Nadel eingefädelt, oder ob dies weiter geklärt werden sollte.
gung - Wiedereinsetzung von Monica