1 Cargar/instalar librerías

pacman::p_load(lme4,foreign, stargazer, texreg, lattice, sjPlot, dplyr, ggeffects) # paquetes a cargar

2 Leer datos base High School and Beyond (HSB)

Como en sesión anterior

mlm = read.dta("http://www.stata-press.com/data/mlmus3/hsb.dta")
dim(mlm)
names(mlm)
attach(mlm)

#Agregacion de casos a escuelas
agg_mlm=mlm %>% group_by(schoolid) %>% summarise_all(funs(mean))

## Warning: funs() is soft deprecated as of dplyr 0.8.0
## please use list() instead
## 
##   # Before:
##   funs(name = f(.))
## 
##   # After: 
##   list(name = ~ f(.))
## This warning is displayed once per session.

3 Coeficientes aleatorios

Estos coeficientes se asocian a la posibilidad de estimar variabilidad de intercepto y de pendiente(s) a través de las unidades de nivel 2. Volviendo a la estimación de un modelo con intercepto y pendiente aleatoria:

yij=γ00+γ10Xij+μ0j+μ1jXij+rij

reg_mlm0=lmer(mathach ~ 1 + ( 1 | schoolid), data = mlm)
reg_mlm1=lmer(mathach ~ 1 + ses + sector +( 1 | schoolid), data = mlm)

reg_mlm2 = lmer(mathach ~ 1 + ses + (1 +ses | schoolid), data=mlm)
htmlreg(reg_mlm2, doctype = FALSE)

En este modelo vemos la estimación de efectos fijos y aleatorios: fijo para intercepto y pendiente, y aleatorios para intercepto y pendiente. Para poder ver el resultado de la estimación, se puede realizar la predicción de los efectos aleatorios y gráficos.

4 Empirical bayes / posterior means para null model

ˆβEB0j=γjˆβ0j+(1−γj)ˆγ00

Donde:

-ˆβEB0j: estimador empirical bayes del intercepto para el grupo j -γj es un ponderador que se define como la confiabilidad del promedio del grupo, y que equivale a

γj=τ00τ00+σ2/nj

-ˆβ0j: es el promedio del grupo -ˆγ00: gran promedio (efecto fijo intercepto)

Entonces, primero seleccionamos una escuela como ejemplo (id 1477)

attach(mlm)
mlm %>% group_by(schoolid) %>% summarise(mean(mathach),count=n()) %>% slice(1:20)

mlm %>% group_by(schoolid) %>% summarise(mean(mathach),count=n()) %>% filter(schoolid==1477)

id1477<- mlm[ which(mlm$schoolid==1477), ] # subset datos para id 1477
id1477
dim(id1477)

γj

varcomp=as.data.frame(VarCorr(reg_mlm0)) # entrega matriz componentes varianza como base de datos

#recordar reg_mlm0 es la estimación de arriba del modelo nulo

tau00=varcomp[1,4]
sigma2=varcomp[2,4]
gama_j=tau00/(tau00+(sigma2/62)) # 62 es el N de la escuela 1477

ˆβ0j

mean_j= mlm %>% filter(schoolid==1477) %>% summarise(mean(mathach))

mean_j

##   mean(mathach)
## 1      14.22847

ˆγ00

gama_00= reg_mlm0@beta # beta es el nombre del estimador gama00 en el objeto
gama_00

## [1] 12.63697

Estimación media posterior intercepto para id 1477

int_eb=gama_j*(mean_j)+((1-gama_j)*(reg_mlm0@beta))

Comparaciones (ver “shrinkage effect”)

int_eb

##   mean(mathach)
## 1      14.11978

mean_j

##   mean(mathach)
## 1      14.22847

mean(mathach)

## [1] 12.74785

Contrastar shrinkage para N más pequeño

gama_j2=tau00/(tau00+(sigma2/20))
int_eb2=gama_j2*(mean_j)+((1-gama_j2)*(reg_mlm0@beta))
int_eb2 # aumenta shrinkage a medida que disminuye N del grupo

##   mean(mathach)
## 1      13.93379

5 Predicción efectos aleatorios

Comandos coef y ranef

reg_mlm0=lmer(mathach ~ 1 + ( 1 | schoolid), data = mlm)
reg_mlm1=lmer(mathach ~ 1 + ses + sector +( 1 | schoolid), data = mlm)
head(ranef(reg_mlm0)) # U0j para cada escuela
head(coef(reg_mlm0)) # comparar, coef entrega intercepto para cada escuela
                    # = gamma_00 + U0j

Para el casode id 1477:

gama_00+1.4828 # = coef

## [1] 14.11977

Modelos con predictores (fijos)

reg_mlm0=lmer(mathach ~ 1 + ( 1 | schoolid), data = mlm)
reg_mlm1=lmer(mathach ~ 1 + ses + sector +( 1 | schoolid), data = mlm)

graf1=ggpredict(reg_mlm1, terms = c("ses","schoolid [sample=8]"), type="re")
plot(graf1)

Predictores aleatorios

reg_mlm2=lmer(mathach ~ 1 + ses + sector +( 1 + ses | schoolid), data = mlm)
graf2=ggpredict(reg_mlm2, terms = c("ses","schoolid [sample=9]"), type="re")
plot(graf2)