Los valores que obtiene de BLUP no se estiman de la misma manera que las estimaciones AZUL de efectos fijos; por convención, los BLUP se denominan predicciones . Cuando 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 dada (digamos un estudiante) se calcula posteriormente a partir de la media estimada y la varianza, y los datos. En un modelo lineal simple, la media se estima (al igual que la varianza residual), pero las puntuaciones observadas se consideran compuestas tanto por eso como por el error, que es una variable aleatoria. En un modelo de efectos mixtos, el efecto para una unidad dada es también una variable aleatoria (aunque en cierto sentido ya se ha realizado).
También puede tratar tales unidades como efectos fijos, si lo desea. En ese caso, los parámetros para esa unidad se estiman como de costumbre. En tal caso, sin embargo, no se estima la media (por ejemplo) de la población de la que se extrajeron las unidades.
Además, la suposición detrás de los efectos aleatorios es que se tomaron muestras al azar de alguna población, y es la población que le interesa. La suposición subyacente de los efectos fijos es que seleccionó esas unidades a propósito porque esas son las únicas unidades que le interesan.
Si se da vuelta y ajusta un modelo de efectos mixtos y predice esos mismos efectos, tienden a "reducirse" hacia la media de la población en relación con sus estimaciones de efectos fijos. Puede pensar en esto como análogo a un análisis bayesiano donde la media y la varianza estimadas especifican un previo normal y el BLUP es como el promedio del posterior que proviene de combinar de manera óptima los datos con el previo.
La cantidad de contracción varía en función de varios factores. Un determinante importante de cuán lejos estarán las predicciones de efectos aleatorios de las estimaciones de efectos fijos es la relación entre la varianza de los efectos aleatorios y la varianza del error. Aquí hay una Rdemostración rápida para el caso más simple con unidades de 5 'nivel 2' con solo ajuste de medios (intercepciones). (Puede pensar en esto como resultados de exámenes para 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)
Entonces, las proporciones de la varianza de los efectos aleatorios a la varianza del error es 1/5 para model 1, 5/5 para model 2y 5/1 para model 3. Tenga en cuenta que usé nivel significa codificación 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 terminar con predicciones de efectos aleatorios que están más cerca de las estimaciones de efectos fijos es cuando tiene más datos. Podemos comparar model 1desde arriba, con su baja proporción de varianza de efectos aleatorios a varianza de error, a una versión ( model 1b) con la misma proporción, pero muchos más datos (observe eso en ni = 500lugar 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, Doug Bates (el autor del paquete R lme4) no le gusta el término "BLUP" y utiliza "modo condicional" en su lugar (ver páginas 22-23 de su borrador del libro lme4 pdf ). En particular, señala en la sección 1.6 que "BLUP" solo puede usarse significativamente para modelos lineales de efectos mixtos.
