La data que analizaremos a continuación está disponible en el paquete merTools y es un ejemplo que se usa bastante para ilutrar el uso de RCM. HSB es una muestra de 7.185 observationes y contiene las siguientes variables:
schid = School ID (160 colegios)
mathach = resultados test de matemáticas
female = sexo (0=male; 1=female)
ses = socio-economic estatus
minority = raza (0=white; 1=non-white students)
schtype= School type (0=public; 1=private/catholic)
meanses= ses promedio para cada colegio en la data
size=número de estudiantes por colegio
Las preguntas que intentaremos responder via RCM son las siguientes:
¿Cuánto varían los colegios en su promedio de matemáticas?
¿Tiene mejores puntajes los colegios con SES promedio más alto?
En promedio, ¿cuál es la relación entre el SES del alumnos con sus puntajes de matemáticas? ¿Es esta asociación similar entre colegios?
¿Es la relación entre SES y matemáticas más fuerte en algunos colegios (e.g., católico, o con alto SES promedio) que en otros (e.g., Público, Bajo SES promedio)?
Random Intercept Model
¿Cuánto varían los colegios en su promedio de matemáticas?
A continuación probamos un modelo inicial (model.0). Cuando se usa la opción REML=TRUE el valor de REML criterion at convergence: corresponde al valor de deviance1 En SPSS este valor se denomina -2 Restricted Log Likelihood cuando se usa la opción REML=FALSE. Este valor es un indicador del ajuste del modelo (más sobre esto más adelante). La opción verbose=TRUE muestra los resultados del proceso de iteración (no es muy relevante)
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
En la Figura de más abajo se observa que los que existe los colegios que tienen un promedio de SES más alto tienen mejores puntajes en matemáticas. Lo importante es darse cuenta que tenemos variables de nivel 1 y nivel 2 combinadas en un mismo modelo.
## Loading required package: ggplot2
Obtengamos formalmente los valores para los factores fijos-
## 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.4066 16.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## meanses -0.004
El Intercepto para todos los colegios es 12.6494 y representa la media de los colegios cuando meanses = 0. El coeficiente para meanses de 5.86 es la pendiente del modelo y quiere decir que por cada punto de meanses los colegios suben en promedio 5.86 puntos en matemáticas. Además vemos que este coeficiente es estadísticamente significativo. Nótese que la varianza residual del modelo es casi la misma que en el modelo sin meanses, pero la varianza entre colegios baja de 8.614 a 2.639. Esto quiere decir que controlando por meanses, la varianza entre colegios es menor.
Llamémosle modelo incondicional al modelo sin meanses, y condicional al modelo CON meanses. Podemos calcular la varianza explicada por meanses comparando la varianza de la variable de Nivel dos entre ambos modelos (8.614 vs 2.639). En este caso, la varianza explicada por meanses es .69 o 69%
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mathach ~ schtype + (1 | schid)
## Data: hsb
##
## REML criterion at convergence: 47080.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0130 -0.7523 0.0253 0.7602 2.7472
##
## Random effects:
## Groups Name Variance Std.Dev.
## schid (Intercept) 6.677 2.584
## Residual 39.151 6.257
## Number of obs: 7185, groups: schid, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.3930 0.2928 158.5372 38.907 < 2e-16 ***
## schtype 2.8049 0.4391 153.5006 6.388 1.9e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## schtype -0.667
De nuevo, nótese la diferencia entre RCM y un modelo de regresión que no considera la variable de agrupación. El coeficiente es aproximadamente el mismo (2.8062 vs 2.8049) pero el error estándar es casi tres veces más pequeño (0.1589 vs 0.4391)
##
## Call:
## lm(formula = mathach ~ schtype, data = hsb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.0023 -5.2163 0.2167 5.4169 13.6289
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.3641 0.1116 101.84 <2e-16 ***
## schtype 2.8062 0.1589 17.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.734 on 7183 degrees of freedom
## Multiple R-squared: 0.04161, Adjusted R-squared: 0.04148
## F-statistic: 311.9 on 1 and 7183 DF, p-value: < 2.2e-16
El error estándar del coeficiente estimado se define como.
\[SE{\hat\gamma_{01} = \sqrt{\frac{\tau_{00}+\frac{\sigma^2_e}{n}}{J\sigma^2}}}\]\[SE{\hat\gamma_{01}}=\sqrt\frac{6.677+\frac{39.151}{45}}{160(.25)}=.434\] El tamaño aproximado de los colegios es 45 y la varianza de la variable colegio es .25. Los valores exactos son:
var(hsb$schtype)
## [1] 0.2499873
mean(table(hsb$size))
## [1] 48.22148
En definitiva, usar RCM no va a incrementar el poder estadístico ya que al considerar la variable de agrupación se REDUCE el tamaño del error estándar. No obstante para INCREMENTAR el poder para detectar un \(\gamma_{01}\) significativo habría que (1) Incrementar el tamaño de los grupos o (2) incrementar el número de grupos.
