Nilai yang Anda dapatkan dari BLUP tidak diestimasi dengan cara yang sama seperti estimasi BLUE untuk efek tetap; menurut konvensi, BLUP disebut sebagai prediksi . Ketika Anda cocok dengan model efek campuran, yang diperkirakan awalnya adalah mean dan varians (dan mungkin kovarians) dari efek acak. Efek acak untuk unit studi yang diberikan (katakanlah seorang siswa) kemudian dihitung dari estimasi rata-rata dan varians, dan data. Dalam model linier sederhana, rata-rata diperkirakan (seperti juga varian residual), tetapi skor yang diamati dianggap terdiri dari keduanya dan kesalahan, yang merupakan variabel acak. Dalam model efek campuran, efek untuk unit yang diberikan juga merupakan variabel acak (meskipun dalam beberapa hal telah terealisasi).
Anda juga dapat memperlakukan unit tersebut sebagai efek tetap, jika diinginkan. Dalam hal ini, parameter untuk unit tersebut diperkirakan seperti biasa. Namun dalam kasus seperti itu, rata-rata (misalnya) dari populasi dari mana unit diambil tidak diperkirakan.
Selain itu, asumsi di balik efek acak adalah bahwa mereka disampel secara acak dari beberapa populasi, dan itu adalah populasi yang Anda sayangi. Asumsi yang mendasari efek tetap adalah bahwa Anda memilih unit-unit itu dengan sengaja karena itu adalah satu-satunya unit yang Anda pedulikan.
Jika Anda berbalik dan menyesuaikan model efek campuran dan memprediksi efek yang sama, mereka cenderung 'menyusut' terhadap rata-rata populasi relatif terhadap perkiraan efek tetap mereka. Anda dapat menganggap ini sebagai analog dengan analisis Bayesian di mana estimasi rata-rata dan varians menentukan sebelum normal dan BLUP seperti rata-rata posterior yang berasal dari penggabungan data secara optimal dengan sebelumnya.
Jumlah susut bervariasi berdasarkan beberapa faktor. Penentuan penting dari seberapa jauh prediksi efek acak akan berasal dari perkiraan efek tetap adalah rasio varians dari efek acak ke varians kesalahan. Berikut ini adalah R
demo cepat untuk kasus paling sederhana dengan unit 5 'level 2' dengan hanya sarana (penyadapan) yang pas. (Anda dapat menganggap ini sebagai nilai ujian untuk siswa di dalam kelas.)
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)
Jadi rasio varians dari efek acak ke varians kesalahan adalah 1/5 untuk model 1
, 5/5 untuk model 2
, dan 5/1 untuk model 3
. Perhatikan bahwa saya menggunakan level means coding untuk model efek tetap. Kita sekarang dapat memeriksa bagaimana perkiraan efek tetap dan efek acak yang diprediksi dibandingkan untuk ketiga skenario ini.
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
Cara lain untuk berakhir dengan prediksi efek acak yang lebih dekat dengan perkiraan efek tetap adalah ketika Anda memiliki lebih banyak data. Kita dapat membandingkan model 1
dari atas, dengan rasio rendah dari efek acak varians ke varians kesalahan, ke versi ( model 1b
) dengan rasio yang sama, tetapi lebih banyak data (perhatikan bahwa ni = 500
alih - alih 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)
Berikut efeknya:
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
Pada catatan yang agak terkait, Doug Bates (penulis paket R lme4) tidak menyukai istilah "BLUP" dan menggunakan "mode kondisional" sebagai gantinya (lihat hlm. 22-23 dari draft lme4 book pdf-nya ). Secara khusus, ia menunjukkan dalam bagian 1.6 bahwa "BLUP" hanya dapat secara bermakna digunakan untuk model efek campuran linier .