¿Por qué los valores estimados de un mejor predictor lineal insesgado (BLUP) difieren de un mejor estimador lineal insesgado (BLUE)?

Question

Entiendo que la diferencia entre ellas está relacionada con si la variable de agrupación en el modelo se estima como efecto fijo o aleatorio, pero no me queda claro por qué no son iguales (si es que no lo son).

Estoy específicamente interesado en cómo funciona esto cuando se utiliza la estimación de áreas pequeñas, si eso es relevante, pero sospecho que la pregunta es relevante para cualquier aplicación de efectos fijos y aleatorios.

Answer 1

Los valores que se obtienen de los BLUP no se estiman de la misma manera que las estimaciones BLUE de los efectos fijos; por convención, los BLUP se denominan predicciones . Cuando se ajusta un modelo de efectos mixtos, lo que se estima inicialmente es la media y la varianza (y posiblemente la covarianza) de los efectos aleatorios. El efecto aleatorio para una unidad de estudio determinada (por ejemplo, un estudiante) se calcula posteriormente a partir de la media y la varianza estimadas, y de los datos. En un modelo lineal simple, la media se estima (al igual que la varianza residual), pero las puntuaciones observadas se consideran compuestas por ella y por el error, que es una variable aleatoria. En un modelo de efectos mixtos, el efecto para una unidad dada es igualmente una variable aleatoria (aunque en cierto sentido ya se ha realizado).

También puede tratar dichas unidades como efectos fijos, si lo desea. En ese caso, los parámetros de esa unidad se estiman como de costumbre. Sin embargo, en este caso no se estima la media (por ejemplo) de la población de la que se extrajeron las unidades.

Además, la suposición que subyace a los efectos aleatorios es que se tomaron muestras al azar de alguna población, y es la población la que te interesa. La suposición que subyace a los efectos fijos es que se seleccionaron esas unidades a propósito porque son las únicas que le interesan.

Si se da la vuelta y se ajusta un modelo de efectos mixtos y se predicen esos mismos efectos, tienden a "reducirse" hacia la media de la población en relación con sus estimaciones de efectos fijos. Se puede pensar en esto como algo análogo a un análisis bayesiano en el que la media y la varianza estimadas especifican una previa normal y el BLUP es como la media de la posterior que viene de combinar óptimamente los datos con la previa.

La cantidad de contracción varía en función de varios factores. Un determinante importante de lo lejos que estarán las predicciones de los efectos aleatorios de las estimaciones de los efectos fijos es la relación entre la varianza de los efectos aleatorios y la varianza del error. A continuación se presenta una rápida R demostración para el caso más simple con 5 unidades de "nivel 2" con sólo medias (interceptos) ajustadas. (Puede pensar en esto como en las puntuaciones de los exámenes de los estudiantes dentro de las clases).

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)

Por lo tanto, la relación entre la varianza de los efectos aleatorios y la varianza del error es de 1/5 para model 1 5/5 para model 2 y 5/1 para model 3 . Obsérvese que he utilizado la codificación de medias de nivel para los modelos de efectos fijos. Ahora podemos examinar cómo se comparan los efectos fijos estimados y los efectos aleatorios predichos para estos tres escenarios.

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

Otra forma de acabar con predicciones de efectos aleatorios que se acercan más a las estimaciones de efectos fijos es cuando se tienen más datos. Podemos comparar model 1 de arriba, con su baja relación entre la varianza de los efectos aleatorios y la varianza del error, a una versión ( model 1b ) con la misma proporción, pero con muchos más datos (observe que ni = 500 en lugar 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)

Aquí están los efectos:

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

---

En una nota algo relacionada, a Doug Bates (el autor del paquete R lme4) no le gusta el término "BLUP" y utiliza en su lugar "modo condicional" (ver las páginas 22-23 de su borrador del libro lme4 pdf ). En particular, señala en la sección 1.6 que el "BLUP" sólo puede utilizarse con sentido para lineal modelos de efectos mixtos.