Por que os valores estimados de um melhor preditor linear e não tendencioso (BLUP) diferem de um melhor estimador linear e imparcial (AZUL)?

20

Entendo que a diferença entre eles está relacionada ao fato de a variável de agrupamento no modelo ser estimada como um efeito fixo ou aleatório, mas não está claro para mim por que eles não são os mesmos (se não são os mesmos).

Estou especificamente interessado em como isso funciona ao usar a estimativa de área pequena, se isso for relevante, mas suspeito que a questão seja relevante para qualquer aplicação de efeitos fixos e aleatórios.

Jeremy Miles
fonte

Respostas:

26

Os valores que você obtém dos BLUPs não são estimados da mesma maneira que as estimativas BLUE de efeitos fixos; por convenção, os BLUPs são referidos como previsões . Quando você ajusta um modelo de efeitos mistos, o que é estimado inicialmente são a média e a variação (e possivelmente a covariância) dos efeitos aleatórios. O efeito aleatório de uma determinada unidade de estudo (digamos, um aluno) é subsequentemente calculado a partir da média e variância estimadas e dos dados. Em um modelo linear simples, a média é estimada (assim como a variância residual), mas as pontuações observadas são consideradas compostas por isso e pelo erro, que é uma variável aleatória. Em um modelo de efeitos mistos, o efeito para uma determinada unidade também é uma variável aleatória (embora, em certo sentido, ela já tenha sido realizada).

Você também pode tratar essas unidades como efeitos fixos, se quiser. Nesse caso, os parâmetros para essa unidade são estimados como de costume. Nesse caso, no entanto, a média (por exemplo) da população da qual as unidades foram extraídas não é estimada.

Além disso, a suposição por trás dos efeitos aleatórios é que eles foram amostrados aleatoriamente em alguma população, e é a população com a qual você se importa. A suposição subjacente aos efeitos fixos é que você selecionou essas unidades intencionalmente porque essas são as únicas unidades com as quais você se importa.

Se você se virar e ajustar um modelo de efeitos mistos e prever esses mesmos efeitos, eles tendem a ser "encolhidos" em relação à média da população em relação às estimativas de efeitos fixos. Você pode pensar nisso como análogo a uma análise bayesiana em que a média estimada e a variância especificam um prior normal e o BLUP é como a média do posterior resultante da combinação otimizada dos dados com o prior.

A quantidade de encolhimento varia com base em vários fatores. Uma determinação importante de quão longe as previsões de efeitos aleatórios estarão das estimativas de efeitos fixos é a razão entre a variação dos efeitos aleatórios e a variação do erro. Aqui está uma Rdemonstração rápida do caso mais simples, com unidades de 5 'nível 2', apenas com ajuste de meios (interceptações). (Você pode pensar nisso como pontuações de teste para alunos das turmas.)

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)

Portanto, as proporções da variação dos efeitos aleatórios para a variação do erro são 1/5 para model 1, 5/5 para model 2e 5/1 para model 3. Note que usei level significa codificar para os modelos de efeitos fixos. Agora podemos examinar como os efeitos fixos estimados e os efeitos aleatórios previstos se comparam para esses três cenários.

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

Outra maneira de terminar com previsões de efeitos aleatórios mais próximas das estimativas de efeitos fixos é quando você tem mais dados. Podemos comparar a model 1partir de cima, com sua baixa proporção de efeitos aleatórios e variação de erro, com uma versão ( model 1b) com a mesma proporção, mas com muito mais dados (observe isso em ni = 500vez de 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)

Aqui estão os efeitos:

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

Em uma nota algo relacionado, Doug Bates (o autor do lme4 pacote R) não gosta do termo "BLUP" e usa "modo condicional" em vez (ver pp. 22-23 do seu projecto de lme4 livro pdf ). Em particular, ele aponta na seção 1.6 que "BLUP" pode ser usado apenas de maneira significativa para modelos lineares de efeitos mistos.

- Reinstate Monica
fonte
3
+1. Mas não tenho certeza se aprecio totalmente a distinção terminológica entre "prever" e "estimar". Portanto, um parâmetro de distribuição é "estimado", mas uma variável latente só pode ser "prevista"? Entendo então corretamente que, por exemplo, as cargas fatoriais na análise fatorial são "estimadas", mas as pontuações fatoriais são "previstas"? Além disso, acho extraordinariamente confuso que algo chamado "melhor preditor imparcial linear" seja realmente um estimador tendencioso (porque implementa encolhimento e, portanto, deve ser tendencioso), se alguém o tratar como um "estimador" dos efeitos fixos. ..
ameba diz Restabelecer Monica
@amoeba, o que significa "melhor", afinal? Melhor o que? É a melhor estimativa da média dos dados ou a melhor combinação das informações contidas nos dados e nas anteriores? A analogia bayesiana ajuda você?
gung - Restabelece Monica
2
Pelo menos, está claro o que "linear" significa :-) Sério, eu encontrei esta resposta muito útil do @whuber na diferença terminológica entre "previsão" e "estimativa". Acho que me esclareceu a terminologia, mas até reforçou meu sentimento de que o BLUP é um estimador, apesar do nome. [cont.]
ameba diz Restabelecer Monica
2
@amoeba, sim, tudo isso é razoável. Mas eu não gostaria de usar o mesmo nome para ambos, porque você está fazendo algo diferente (ou seja, as equações diferem) e é útil que os nomes sejam distintos.
gung - Restabelece Monica
1
@amoeba, alterei o fraseado no primeiro parágrafo para enfatizar esses termos, de modo a não reificar a "previsão", mas para manter a distinção. Veja se você acha que eu passei a agulha ou se deve ser esclarecido mais.
gung - Restabelece Monica