Incluyendo variables en el Nivel 1 y en el Nivel 2
Gonzalo J. Muñoz
2019-07-03
Comparación de modelos bottom-up
Cuando incluimos variables de Nivel 1 y Nivel 2 hay que tener un sistema para determinar qué efectos son relevantes para el modelo definitivo. Hox (2010) recomienda seguir cinco pasos para evaluar un modelo multinivel. Este es un proceso bottom-up que parte con un modelo donde sólo se estima la varianza del intercepto sin ningún predictor y avanza progresivamente hasta un modelo con una interacción entre una variable de Nivel 2 y una variable de Nivel 1.
Paso 1: Intercepto Paso 2: Factores fijos de Nivel 1 Paso 3: Factores fijos de Nivel 1 + Factores Fijos de Nivel 2 Paso 4: Factores aleatorios factores de Nivel 1 significativos Paso 5: Interacción entre factores de Nivel 1 y Nivel 2
Estadísticos descriptivos
En general es una buena idea mirar los estadísticos descriptivos de la data que estamos analizando. Hay que tener en consideración que las variables de Nivel 1 están anidadas en el Nivel 2, por lo que sólo en el caso del Nivel 1 podemos usar el promedio de todas las observaciones. Las variables de Nivel 2 deben ser agregadas de acuerdo a su variable de agrupación para poder describirlas correctamente (i.e, grados de libertad correctos). Por ahora, es importante consignar que la variable ses tiene un promedio de 0, por lo que se puede inferir que está centrada (más sobre esto más adelante)1 En la versión agregada de la data el tamaño (size) se refiere al tamaño del colegio, no de la muestra tomada del colegio para el estudio. Lo mismo ocurre con meanses
hsb<-readRDS(file = "hsb.rds")
require(psych)
## Loading required package: psych
describe(hsb[,2:5])
## vars n mean sd median trimmed mad min max range skew
## minority 1 7185 0.27 0.45 0.00 0.22 0.00 0.00 1.00 1.00 1.01
## female 2 7185 0.53 0.50 1.00 0.54 0.00 0.00 1.00 1.00 -0.11
## ses 3 7185 0.00 0.78 0.00 0.02 0.85 -3.76 2.69 6.45 -0.23
## mathach 4 7185 12.75 6.88 13.13 12.91 8.12 -2.83 24.99 27.82 -0.18
## kurtosis se
## minority -0.98 0.01
## female -1.99 0.01
## ses -0.38 0.01
## mathach -0.92 0.08
aggdata<-aggregate(hsb[,6:8],by=list(hsb$schid),FUN=mean)
#El `-1` es para simplemente para evitar que aparezca la media de de la variable de agrupación
describe(aggdata[,-1])
## vars n mean sd median trimmed mad min max
## size 1 160 1097.83 629.51 1061.00 1058.24 695.34 100.00 2713.00
## schtype 2 160 0.44 0.50 0.00 0.42 0.00 0.00 1.00
## meanses 3 160 0.00 0.41 0.03 0.01 0.47 -1.19 0.83
## range skew kurtosis se
## size 2613.00 0.46 -0.61 49.77
## schtype 1.00 0.25 -1.95 0.04
## meanses 2.02 -0.28 -0.45 0.03
Paso 1: Intercepto
require(lme4)
## Loading required package: lme4
## Loading required package: Matrix
require(lmerTest)
## Loading required package: lmerTest
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | schid)
## Data: hsb
##
## REML criterion at convergence: 47116.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0631 -0.7539 0.0267 0.7606 2.7426
##
## Random effects:
## Groups Name Variance Std.Dev.
## schid (Intercept) 8.614 2.935
## Residual 39.148 6.257
## Number of obs: 7185, groups: schid, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6370 0.2444 156.6473 51.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La correlación intraclase en este caso sería: \[\rho=\frac{\tau_{00}}{\tau_{00}+\sigma^2_e}=\frac{8.61}{8.61+39.15}=0.18\]
Paso 2: Factores fijos de Nivel 1
model.2<-lmer(mathach ~ ses + (1|schid),REML=TRUE,verbose=FALSE, data=hsb)
summary(model.2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ ses + (1 | schid)
## Data: hsb
##
## REML criterion at convergence: 46645.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.12607 -0.72720 0.02188 0.75772 2.91912
##
## Random effects:
## Groups Name Variance Std.Dev.
## schid (Intercept) 4.768 2.184
## Residual 37.034 6.086
## Number of obs: 7185, groups: schid, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6575 0.1880 148.3022 67.33 <2e-16 ***
## ses 2.3902 0.1057 6838.0776 22.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## ses 0.003
Una manera de comparar los modelos 1 y 2 es determinar la proporción de varianza explicada por el modelo 2. Hox (2010) propone comparar los residuos de cada nivel por separado.2 b=baseline;m=comparison model Para los residuos de nivel 1 tenemos:
Mientras que para los residuos de nivel 2 tenemos:
\[R^2_2=(\frac{\sigma^2_{u0|b}-\sigma^2_{u0|m}}{\sigma^2_{u0|b}})\]\[R^2_1=(\frac{8.614-4.768}{8.614})=0.44548\approx0.45\] Es interesante notar que la variable de Nivel 1 reduce la varianza entre grupos o clases. Esto es así porque la composición de los grupos varía en SES. Si algunos grupos tienen estudiantes mayoritariamente con alto SES y otros grupos mayoritariamente estudiantes con bajo SES la variable individual puede influier en el promedio de los grupos. Si, en cambio, añadimos una variable de Nivel 2 (e.g., meanses) ese factor no debiese tener un efecto en reducir la varianza residual del nivel 1, o si ocurre sería ruido derivado del proceso de estimación.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ meanses + (1 | schid)
## Data: hsb
##
## REML criterion at convergence: 46961.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.13480 -0.75256 0.02409 0.76773 2.78501
##
## Random effects:
## Groups Name Variance Std.Dev.
## schid (Intercept) 2.639 1.624
## Residual 39.157 6.258
## Number of obs: 7185, groups: schid, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6494 0.1493 153.7425 84.74 <2e-16 ***
## meanses 5.8635 0.3615 153.4067 16.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## meanses -0.004
En cualquier caso, al comparar modelos con el método de comparación de varianzas puede darse el caso que al incluir una variable explicatoria aument la varianza no explicada lo que conduciría a valores de r cuadrado negativos, lo que obviamente es un problema.
Otro método para comparar modelos es comparar sus respectivas desviaciones (deviances). Cuando se usa REML este valor aparece etiquetado como REML criterion at convergence: y se llama así porque es -2 \(\times\) ln el valor de la función likelihood al momento de converger. La desviación del modelo 1 es 47116.8 mientras que en el modelo 2 46645.2lo que sugiere que el ajuste del modelo 2 es mejor que el ajuste del modelo 1. La diferencia entre ambos modelos se puede evaluar mediante una prueba \(\chi^2\) con grados de libertad igual a la diferencia de parámetros estimados en los dos modelos.
El modelo 1 se obtiene removiendo la variable ses del modelo 2. La diferencia entre la desviación del modelo 1 y la desviación del modelo 2 refleja la pérdida de ajuste del modelo al no incluir este predictor.
47116.8-46645.2
## [1] 471.6
El valor crítico de \(\chi^2\) para df=1 se puede obtener mediante
qchisq(p=.95,df=1)
## [1] 3.841459
Es decir, si la diferencia entre el modelo 1 y 2 es igual o mayor que 3.84 sería estadísticamente significativa. Para obtener el valor exacto de p para el valor de 471.6 podemos usar:
pchisq(471.6, df=1,lower.tail = FALSE)
## [1] 1.437468e-104
Una función muy util para comparar modelos es anova. Es importante notar que el modelo es reajustado pero usando ahora maximum likelihood en vez de restricted maximum likelihood, por lo que los valores de las desviaciones del modelo 1 y 2 son algo distintos de los que obtuvimos en el paso anterior.
Paso 4: Factores aleatorios factores de Nivel 1 significativos
Si tuviésemos más de un factor fijo de Nivel 1, es recomendado evaluarlos 1 a 1 y dejar en el modelo final sólo aquellos que tengan efectos estadísticamente significativos
model.4<-lmer(mathach ~ ses + schtype + (ses|schid),REML=TRUE,verbose=FALSE, data=hsb)
summary(model.4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ ses + schtype + (ses | schid)
## Data: hsb
##
## REML criterion at convergence: 46601.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1373 -0.7296 0.0225 0.7568 2.8920
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schid (Intercept) 3.9646 1.991
## ses 0.4343 0.659 0.55
## Residual 36.8008 6.066
## Number of obs: 7185, groups: schid, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.4729 0.2315 153.7963 49.568 < 2e-16 ***
## ses 2.3854 0.1179 157.8423 20.238 < 2e-16 ***
## schtype 2.5408 0.3445 151.3087 7.375 1.01e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ses
## ses 0.228
## schtype -0.655 -0.079
Dado que en este caso la variable de Nivel 2 (schtype) sólo modera el intercepto, la pendiente para ambos tipos de colegio se asume que es igual (en promedio). En este caso tenemos los siguientes valores predichos: b0_publico = 11.47 b0_catolico = 11.47 + 2.54 = 14.01 b1 = 2.39
La matriz de varianza-covarianza de los efectos del modelo se puede escribir como:
donde \(\tau{00}\) y \(\tau{11}\) corresponden a la varianza del intercepto y la pendiente, respectivamente; mientras que \(\tau{10}=\tau{01}\) es la covarianza entre el intercepto y la pendiente.
Podemos extraer los valores para la matriz de varianza-covarianza de los efectos del modelo mediante la función VarCorr
Visualizar modelo 4 con correlación entre intercepto y pendiente
Como se puede apreciar más abajo, el modelo permite que cada colegio tenga su propio intecepto y pendiente. Esto quiere decir que nuestro modelo incluye un efecto de interacción inter-nivel (i.e., cross-level effect) que se observa en que las tanto las pendientes de cada modelo se cruzan. ¿Depende la relación entre SES (individual) y matemáticas (individual) del tipo de colegio (grupal)? Mirando la figura de más abajo, la respuesta parece se que sí.
Veamos si el tipo de colegio (público vs católico) puede servir para explicar esta diferencia. Primero, veamos gráficamente cómo se podría ver este modelo.3 Considere que el tipo de colegio es 0=Público y 1=CatólicoAparentemente la relación entre ses y mathach es más fuerte en el caso de los colegios públicos
Paso 5: Interacción entre factores de Nivel 1 y Nivel 2
Ecuaciones
Modelo Nivel 1: \[Y_{ij}=\beta_{0_j}+\beta_{1j}SES_{ij}+e_{ij}\] Modelos Nivel 2: \[\beta_{0_j}=\gamma_{00}+\gamma_{01}schtype_j+u_{0j}\]\[\beta_{1_j}=\gamma_{10}+\gamma_{11}schtype_j+u_{1j}\]
En el siguiente modelo (model.5) podemos evaluar si la relación entre SES y mathach varía por tipo colegio (público vs católico). Para incluir la interacción podemos usar ses + schtype + ses:schtype o, en su forma abreviada ses*schtype4 A veces se agrega un 1 delante de ses para señalar el intercepto (i.e., 1 + ses*schtype + (1 + ses|schid)), pero esto no es necesario
Si bien el modelo 5 es superior al modelo 4, hay un warning en el modelo 5 sobre la convergencia del modelo. Una forma de resolver este problema (y otros problemas de interpretación) es centrar las variables.
Bauer y Curran (2005) recomiendan centrar los predictores de Nivel 1 en relación al grupo (group-mean centering) y centrar los predictores de Nivel 2 en relación a la gran media. Usando group-mean centering vemos que ahora no hay problemas de convergencia.5https://cran.r-project.org/web/packages/lme4/vignettes/lmerperf.html
Lo que ocurre porque para poder hacer este test tenemos que tener los mismos predictores en cada modelo. Agregemos al modelo 4 la variable ses centrada
#Ojo que al centrar schtype hay que usar la data agregada
hsb$gm.cen.schtype<-hsb$schtype-mean(aggdata$schtype)
model.4C<-lmer(mathach ~ grp.cen.ses+gm.cen.schtype + (grp.cen.ses|schid),REML=TRUE,data=hsb)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00414013
## (tol = 0.002, component 1)
summary(model.4C)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ grp.cen.ses + gm.cen.schtype + (grp.cen.ses | schid)
## Data: hsb
##
## REML criterion at convergence: 46667.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.12077 -0.73655 0.01832 0.75294 2.85781
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## schid (Intercept) 6.8585 2.6189
## grp.cen.ses 0.6971 0.8349 0.58
## Residual 36.7063 6.0586
## Number of obs: 7185, groups: schid, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6349 0.2200 150.3171 57.440 < 2e-16 ***
## grp.cen.ses 2.2147 0.1280 154.0569 17.306 < 2e-16 ***
## gm.cen.schtype 3.4604 0.4239 151.7308 8.163 1.17e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grp.c.
## grp.cen.ses 0.287
## gm.cn.schty -0.014 0.002
## convergence code: 0
## Model failed to converge with max|grad| = 0.00414013 (tol = 0.002, component 1)
