Mixed-Model

為什麼最佳線性無偏預測器 (BLUP) 的估計值與最佳線性無偏估計器 (BLUE) 不同?

  • October 31, 2014

我知道它們之間的差異與模型中的分組變量是估計為固定效應還是隨機效應有關,但我不清楚為什麼它們不一樣(如果它們不一樣)。

如果相關的話,我對使用小面積估計時它的工作原理特別感興趣,但我懷疑這個問題與固定和隨機效應的任何應用有關。

您從 BLUP 獲得的值的估計方式與 BLUE 對固定效應的估計不同;按照慣例,BLUP 被稱為預測。當您擬合混合效應模型時,最初估計的是隨機效應的均值和方差(可能還有協方差)。隨後根據估計的均值和方差以及數據計算給定學習單元(例如學生)的隨機效應。在一個簡單的線性模型中,均值是估計的(殘差也是如此),但觀察到的分數被認為是由該值和誤差組成,誤差是一個隨機變量。在混合效應模型中,給定單位的效應同樣是一個隨機變量(儘管在某種意義上它已經實現)。

如果您願意,您也可以將這些單位視為固定效果。在這種情況下,該單元的參數將照常估計。然而,在這種情況下,不會估計從中提取單位的總體的平均值(例如)。

此外,隨機效應背後的假設是它們是從某些人群中隨機抽樣的,而您關心的是人群。固定效應的假設是您有目的地選擇這些單位,因為這些是您唯一關心的單位。

如果您轉身擬合混合效應模型並預測這些相同的效應,則相對於其固定效應估計值,它們往往會向總體均值“縮小”。您可以將其視為類似於貝葉斯分析,其中估計的均值和方差指定了一個正常的先驗,而 BLUP 類似於將數據與先驗進行最佳組合的後驗均值。

收縮量因幾個因素而異。隨機效應預測與固定效應估計之間的差距的一個重要決定因素是隨機效應的方差與誤差方差的比率。這是最簡單情況的快速R演示,其中包含 5 個“2 級”單元,只有均值(截距)適合。(您可以將其視為班級內學生的考試成績。)

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)

因此,隨機效應的方差與誤差方差之比為 1/5 model 1、 5/5model 2和 5/1 model 3。請注意,我對固定效應模型使用了級別表示編碼。我們現在可以檢查估計的固定效應和預測的隨機效應如何比較這三種情景。

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

獲得更接近固定效應估計的隨機效應預測的另一種方法是當您擁有更多數據時。我們可以model 1從上面比較,其隨機效應方差與誤差方差的比率較低,model 1b與具有相同比率但數據更多的版本 ( ) 相比(請注意,ni = 500而不是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)

以下是效果:

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


在一些相關的說明中,Doug Bates(R 包 lme4 的作者)不喜歡術語“BLUP”,而是使用“條件模式”(參見他的草案 lme4 書pdf的第 22-23 頁)。特別是,他在第 1.6 節中指出,“BLUP”只能有意義地用於線性混合效應模型。

引用自:https://stats.stackexchange.com/questions/122218

comments powered by Disqus