Wartości, które otrzymujesz z BLUP nie są szacowane w taki sam sposób, jak NIEBIESKIE oszacowania ustalonych efektów; zgodnie z konwencją BLUP są nazywane prognozami . Po dopasowaniu modelu efektów mieszanych początkowo szacuje się średnią i wariancję (i ewentualnie kowariancję) efektów losowych. Efekt losowy dla danej jednostki badawczej (powiedzmy ucznia) jest następnie obliczany na podstawie oszacowanej średniej i wariancji oraz danych. W prostym modelu liniowym szacuje się średnią (podobnie jak wariancję resztkową), ale obserwowane wyniki uważa się za złożone zarówno z tego, jak i z błędu, który jest zmienną losową. W modelu efektów mieszanych efekt dla danej jednostki jest również zmienną losową (chociaż w pewnym sensie została już zrealizowana).
Możesz także traktować takie jednostki jako efekty stałe, jeśli chcesz. W takim przypadku parametry dla tej jednostki są szacowane jak zwykle. W takim przypadku jednak średnia (na przykład) populacji, z której zostały pobrane jednostki, nie jest szacowana.
Co więcej, założeniem za efektami losowymi jest to, że próbki zostały pobrane losowo z pewnej populacji, i jest to populacja, na której Ci zależy. Założeniem leżącym u podstaw stałych efektów jest to, że wybrałeś te jednostki celowo, ponieważ są to jedyne jednostki, na których ci zależy.
Jeśli odwrócisz się i dopasujesz model efektów mieszanych i przewidzisz te same efekty, będą one „kurczyć się” w stosunku do średniej populacji w stosunku do ich stałych oszacowań efektów. Można to traktować jako analogię do analizy bayesowskiej, w której szacowana średnia i wariancja określają normalny wcześniejszy, a BLUP jest jak średnia z tyłu, która pochodzi z optymalnego połączenia danych z wcześniejszym.
Wielkość skurczu zależy od kilku czynników. Ważnym wyznacznikiem odległości przewidywań efektów losowych od oszacowań efektów stałych jest stosunek wariancji efektów losowych do wariancji błędu. Oto szybkie Rdemo dla najprostszego przypadku z jednostkami 5 'poziomu 2' z dopasowaniem tylko środków (przechwytuje). (Można to traktować jako wyniki testów dla uczniów w ramach klas).
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)
Tak więc stosunek wariancji efektów losowych do wariancji błędu wynosi 1/5 dla model 1, 5/5 dla model 2i 5/1 dla model 3. Zauważ, że użyłem poziomu oznacza kodowanie modeli efektów stałych. Możemy teraz zbadać, jak szacowane efekty stałe i przewidywane efekty losowe są porównywane dla tych trzech scenariuszy.
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
Innym sposobem na uzyskanie losowych prognoz efektów, które są bliższe oszacowaniom efektów stałych, jest posiadanie większej ilości danych. Możemy porównać model 1z góry, z niskim współczynnikiem wariancji efektów losowych do wariancji błędów, do wersji ( model 1b) o tym samym współczynniku, ale o wiele więcej danych (zauważ, że ni = 500zamiast 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)
Oto efekty:
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
W dość pokrewnej notatce Doug Bates (autor pakietu R lme4) nie lubi terminu „BLUP” i zamiast tego używa „trybu warunkowego” (patrz str. 22-23 jego wersji pdf książki Lme4 ). W szczególności zwraca uwagę w sekcji 1.6, że „BLUP” może być sensownie używany tylko w liniowych modelach efektów mieszanych.
