Les valeurs que vous obtenez des BLUP ne sont pas estimées de la même manière que les estimations BLEUES des effets fixes; par convention, les BLUP sont appelés prédictions . Lorsque vous ajustez un modèle à effets mixtes, ce qui est estimé initialement est la moyenne et la variance (et éventuellement la covariance) des effets aléatoires. L'effet aléatoire pour une unité d'étude donnée (disons un étudiant) est ensuite calculé à partir de la moyenne et de la variance estimées, et des données. Dans un modèle linéaire simple, la moyenne est estimée (tout comme la variance résiduelle), mais les scores observés sont considérés comme composés à la fois de cela et de l'erreur, qui est une variable aléatoire. Dans un modèle à effets mixtes, l'effet pour une unité donnée est également une variable aléatoire (bien que dans un certain sens, il ait déjà été réalisé).
Vous pouvez également traiter ces unités comme des effets fixes, si vous le souhaitez. Dans ce cas, les paramètres de cette unité sont estimés comme d'habitude. Dans un tel cas cependant, la moyenne (par exemple) de la population dont les unités ont été tirées n'est pas estimée.
De plus, l'hypothèse derrière les effets aléatoires est qu'ils ont été échantillonnés au hasard dans une certaine population, et c'est la population qui vous intéresse. L'hypothèse sous-jacente aux effets fixes est que vous avez choisi ces unités délibérément parce que ce sont les seules unités qui vous intéressent.
Si vous vous retournez et ajustez un modèle à effets mixtes et que vous prédisez ces mêmes effets, ils ont tendance à être «réduits» vers la moyenne de la population par rapport à leurs estimations d'effets fixes. Vous pouvez considérer cela comme analogue à une analyse bayésienne où la moyenne et la variance estimées spécifient un a priori normal et le BLUP est comme la moyenne du postérieur qui résulte de la combinaison optimale des données avec l'a priori.
La quantité de retrait varie en fonction de plusieurs facteurs. Le rapport entre la variance des effets aléatoires et la variance d'erreur est un facteur déterminant important de la distance entre les prévisions des effets aléatoires et les estimations des effets fixes. Voici une Rdémo rapide pour le cas le plus simple avec des unités de 5 'niveau 2' avec seulement des moyens (interceptions) en forme. (Vous pouvez considérer cela comme des résultats aux tests pour les élèves des classes.)
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)
Ainsi, les rapports de la variance des effets aléatoires à la variance d'erreur sont de 1/5 pour model 1, 5/5 pour model 2et 5/1 pour model 3. Notez que j'ai utilisé le codage de niveau pour les modèles d'effets fixes. Nous pouvons maintenant examiner comment les effets fixes estimés et les effets aléatoires prévus se comparent pour ces trois scénarios.
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
Une autre façon de se retrouver avec des prédictions d'effets aléatoires plus proches des estimations d'effets fixes consiste à disposer de plus de données. Nous pouvons comparer model 1ci-dessus, avec son faible rapport de la variance des effets aléatoires à la variance de l'erreur, à une version ( model 1b) avec le même rapport, mais beaucoup plus de données (notez qu'au ni = 500lieu 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)
Voici les effets:
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
Sur une note quelque peu connexe, Doug Bates (l'auteur du package R lme4) n'aime pas le terme "BLUP" et utilise plutôt le "mode conditionnel" (voir pp. 22-23 de son projet de livre lme4 pdf ). En particulier, il souligne dans la section 1.6 que "BLUP" ne peut être utilisé de manière significative que pour les modèles linéaires à effets mixtes.
