为什么最佳线性无偏预测器 (BLUP) 的估计值与最佳线性无偏估计器 (BLUE) 不同?

机器算法验证 混合模式 蓝色的 虚张声势 小面积估计
2022-01-29 19:39:52

我知道它们之间的差异与模型中的分组变量是被估计为固定效应还是随机效应有关,但我不清楚为什么它们不一样(如果它们不一样)。

如果相关的话,我对使用小面积估计时它的工作原理特别感兴趣,但我怀疑这个问题与固定和随机效应的任何应用有关。

1个回答

您从 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”只能有意义地用于线性混合效应模型。