Modelos lineales mixtos para explicar los puntajes de las pruebas ICFES en el año 2018
Santiago Restrepo Álvarez
library(lme4)
library(nlme)
library(ggpubr)
library(ggplot2)
library(knitr)
punt<-read.csv("clasifcol.csv",header=T,stringsAsFactors = TRUE)Los datos del siguiente análisis fueron extraidos de la página de datos abiertos de la alcaldía de Medellín (http://medata.gov.co/dataset/clasificación-establecimientos-prueba-saber-11)
kable(head(punt))| establecimiento_educativo | comuna_establecimiento | sector_educativo | numero_matriculados | total |
|---|---|---|---|---|
| inst educ jose maria velaz | 2 | oficial | 105 | 0.7134 |
| inst educ barrio santa cruz | 2 | oficial | 143 | 0.6298 |
| inst educ cefa | 10 | oficial | 2807 | 0.7491 |
| inst educ jose maria bernal | 16 | oficial | 193 | 0.7630 |
| inst educ camilo torres restrepo | 4 | oficial | 88 | 0.6670 |
| col la pastora-en administracion (ac) | 1 | no oficial | 46 | 0.6394 |
Nota: Se muestran solamente las variables utilizadas para evitar una tabla muy grande. Para más información de las variables, remitirse al link anteriormente mencionado.
Variables
Estos datos tienen 317 observaciones con las siguientes variables:
año: El año de realización de la prueba.
establecimiento_educativo: Nombre de la institución.
prestacion_servicio: Tipo de prestación de servicio con 3 niveles(Oficial, no oficial y contratación).
comuna_establecimiento: Comuna donde está ubicada el establecimiento educativo.
sector_educativo: Variable que dice si el colegio es público o privado(oficial o no oficial)
clasificacion: Clasificación a partir del puntaje logrado en la prueba.
numero_matriculados: Número total de matriculados.
matematicas: Puntaje en matemáticas de la prueba.
ciencias_naturales: Puntaje en ciencias naturales de la prueba.
sociales_ciudadanas: Puntaje en ciencias sociales de la prueba.
lectura_escritura: Puntaje en letura y escritura de la prueba.
ingles: Puntaje en inglés de la prueba.
total: Puntaje ponderado total de la prueba.
Limpieza y formato de datos.
punt$comuna_establecimiento<-as.factor(punt$comuna_establecimiento)
#Formato del nombre a chr
punt$establecimiento_educativo<-as.character(punt$establecimiento_educativo)Análisis exploratorio.
#Exploratorio
#Gráfico por comunas
ggplot(data = punt, aes(x = numero_matriculados, y = total, color = comuna_establecimiento)) +
geom_point() +
theme_bw() +
facet_wrap(~ comuna_establecimiento) +
theme(legend.position = "none")
Notamos que las comunas que parecen tener mejor puntaje son las 14(Poblado), la 12(Américas) y la 11(Laureles), parece que la comuna es candidata a ser un efecto aleatorio para el modelamiento.
#Gráfico por público o privado
ggplot(punt, aes(numero_matriculados, total, color=sector_educativo ) ) +
geom_point() +
labs(colour="Tipo de colegio") + theme_bw()
Se nota que hay algunos outliers a la derecha del eje x del gráfico, por consiguiente se analizarán solamente los establecimientos educativos con menos de 1000 estudiantes.
punt<-punt[which(punt$numero_matriculados<1000) ,]
ggplot(punt, aes(numero_matriculados, total, color=sector_educativo ) ) +
geom_point() +
labs(colour="Tipo de colegio") + theme_bw()
Aquí se puede notar que el tipo de colegio puede servir como efecto aleatorio igual que la comuna.
Etapa de modelamiento
Modelo lineal clásico
\[\begin{align} Total&\sim N(\mu,\sigma^2)\\ \mu&=\beta_0+\beta_1\text{(Número_matriculados)}_{ij}+\beta_2(\text{prestacion_servicio})_{ij} \end{align}\]
mod1<-gls(total~numero_matriculados+sector_educativo,method='REML',data=punt)
summary(mod1)## Generalized least squares fit by REML
## Model: total ~ numero_matriculados + sector_educativo
## Data: punt
## AIC BIC logLik
## -1039.753 -1024.426 523.8766
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.7297757 0.005537226 131.79446 0
## numero_matriculados 0.0001593 0.000023312 6.83521 0
## sector_educativooficial -0.0710233 0.006002064 -11.83314 0
##
## Correlation:
## (Intr) nmr_mt
## numero_matriculados -0.535
## sector_educativooficial -0.499 -0.297
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.76996465 -0.64702484 0.02982624 0.60891742 2.69363728
##
## Residual standard error: 0.05015373
## Degrees of freedom: 344 total; 341 residualModelo de intercepto aleatorio y pendiente fija
Para este modelo se usará un modelo lineal mixto con la comuna como intercepto aleatorio.
\[\begin{align} Total_{ij}|b_0&\sim N(\mu_{ij},\sigma^2_{Total})\\ \mu_{ij}&=\beta_0+\beta_1(\text{Número_matriculados})_{ij}+\beta_2(\text{Prestación_servicio})_{ij}+b_{0i}\\ b_0&\sim N(0,\sigma^2_{b0}) \end{align}\]
mod2<-lme(total~numero_matriculados+sector_educativo, random=~1|comuna_establecimiento,method="REML",data=punt)
summary(mod2)## Linear mixed-effects model fit by REML
## Data: punt
## AIC BIC logLik
## -1089.746 -1070.586 549.8729
##
## Random effects:
## Formula: ~1 | comuna_establecimiento
## (Intercept) Residual
## StdDev: 0.02694406 0.04408794
##
## Fixed effects: total ~ numero_matriculados + sector_educativo
## Value Std.Error DF t-value p-value
## (Intercept) 0.7233937 0.007993147 321 90.50174 0
## numero_matriculados 0.0001388 0.000021162 321 6.56086 0
## sector_educativooficial -0.0553242 0.005862775 321 -9.43652 0
## Correlation:
## (Intr) nmr_mt
## numero_matriculados -0.307
## sector_educativooficial -0.353 -0.316
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.310832521 -0.599721865 0.007062578 0.704593392 2.895077895
##
## Number of Observations: 344
## Number of Groups: 21Modelo de intercepto y pendiente aleatoria
Para este modelo se usará un intercepto aleatorio a partir de la comuna y una pendiente aleatoria a partir del tipo de colegio(público o privado).
\[\begin{align} Total_{ij}|b_0,b_1 &\sim N(\mu_{ij},\sigma^2_{Total}) \\ \mu_{ij}&=\beta_0+\beta_1(\text{Número_matriculados)}_{ij}+\beta_2(\text{Prestación_servicio})_{ij}+b_{0i}+b_{1i}(\text{Prestación_servicio})\\ \begin{pmatrix} b_{0} \\ b_{1} \end{pmatrix}&\sim N\bigg(\begin{bmatrix} 0 \\ 0 \end{bmatrix},\begin{bmatrix} \sigma^2_{b0} & \sigma_{b10}\\ \sigma_{b01} & \sigma^2_{b1} \end{bmatrix}\bigg)\\ \end{align}\]
mod3<-lme(total~numero_matriculados+sector_educativo, random=~1+sector_educativo|comuna_establecimiento,method="REML",data=punt,control=lmeControl(msMaxIter = 1000, msMaxEval = 1000))
summary(mod3)## Linear mixed-effects model fit by REML
## Data: punt
## AIC BIC logLik
## -1099.824 -1073.001 556.9121
##
## Random effects:
## Formula: ~1 + sector_educativo | comuna_establecimiento
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.04116649 (Intr)
## sector_educativooficial 0.02830715 -0.945
## Residual 0.04290702
##
## Fixed effects: total ~ numero_matriculados + sector_educativo
## Value Std.Error DF t-value p-value
## (Intercept) 0.7193851 0.010979719 321 65.51945 0
## numero_matriculados 0.0001420 0.000020594 321 6.89773 0
## sector_educativooficial -0.0546316 0.009045563 321 -6.03960 0
## Correlation:
## (Intr) nmr_mt
## numero_matriculados -0.210
## sector_educativooficial -0.809 -0.210
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.66937414 -0.64044719 -0.03531933 0.67943191 3.42683146
##
## Number of Observations: 344
## Number of Groups: 21Comparación de los modelos
Se usará una tabla ANOVA para comparar los 3 modelos.
Correlación entre valores ajustados y reales.
corre <- data.frame(cor_mod1 = cor(fitted(mod1),
punt$total),
cor_mod2 = cor(fitted(mod2),punt$total),
cor_mod3 = cor(fitted(mod3),punt$total))
kable(corre)| cor_mod1 | cor_mod2 | cor_mod3 |
|---|---|---|
| 0.5554582 | 0.7013152 | 0.7245113 |
A partir de la correlación de las predicciones de los 3 modelos y los datos reales, el modelo 3 es el mejor de todos.
kable(anova(mod1,mod2,mod3))| call | Model | df | AIC | BIC | logLik | Test | L.Ratio | p-value | |
|---|---|---|---|---|---|---|---|---|---|
| mod1 | gls(model = total ~ numero_matriculados + sector_educativo, data = punt, method = “REML”) | 1 | 4 | -1039.753 | -1024.426 | 523.8766 | NA | NA | |
| mod2 | lme.formula(fixed = total ~ numero_matriculados + sector_educativo, data = punt, random = ~1 | comuna_establecimiento, method = “REML”) | 2 | 5 | -1089.746 | -1070.586 | 549.8729 | 1 vs 2 | 51.99268 | 0.0000000 |
| mod3 | lme.formula(fixed = total ~ numero_matriculados + sector_educativo, data = punt, random = ~1 + sector_educativo | comuna_establecimiento, method = “REML”, control = lmeControl(msMaxIter = 1000, msMaxEval = 1000)) | 3 | 7 | -1099.824 | -1073.001 | 556.9121 | 2 vs 3 | 14.07828 | 0.0008769 |
A partir de la tabla ANOVA, podemos inferir que el modelo 3(intercepto y pendientes aleatorias) es el mejor de los 3 modelos debido a tener menor AIC, BIC y la prueba para likelihood ratio.
Validación de supuestos del modelo escogido.
Normalidad.
ggplot(punt, aes(x=total)) +
geom_histogram(aes(y=..density..), colour='black', fill='white') +
geom_density(alpha=0.4, fill="green") +
labs(x='Puntaje ponderado de la prueba', y = 'Frecuencia', title= 'Histograma ICFES')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La variable del puntaje parece ser normal aunque tiene una cola derecha un poco pesada.
Análisis de residuales.
plot(mod3)
Podemos ver que el supuesto de que la media de los residuales es igual a 0 parece cumplirse.
Añadiendo un gráfico de cuantiles para los residuales, obtenemos:
residuales <- residuals(mod3)
ajustados <- fitted(mod3)
resid <- data.frame(residuales = residuales, fitt = ajustados)
ggplot(resid, aes(sample = residuales)) +
stat_qq() +
stat_qq_line(color = "red") +
theme_bw()
La cola superior del gráfico de cuantiles es un poco más pesada de lo que debería ser pero en su mayoría se debería verificar a partir de una prueba de hipótesis el supuesto de normalidad en los residuales.
Conclusión.
El mejor modelo es el modelo 3 con la siguiente ecuación:
\[\begin{align} Total_{ij}|b_0,b_1 &\sim N(\mu_{ij},\sigma^2_{Total}\approx 0.0429) \\ \mu_{ij}&=0.7193851+0.000142(\text{Número_matriculados})_{ij}-0.0546316(\text{Prestación_servicio})_\text{oficialj}+b_{0i}+b_{1i}(\text{Prestación_servicio})\\ \begin{pmatrix} b_{0} \\ b_{1} \end{pmatrix}&\sim N\bigg(\begin{bmatrix} 0 \\ 0 \end{bmatrix},\begin{bmatrix} \sigma^2_{b0}=0.04116649 & \sigma_{b10}=-0.000109046\\ \sigma_{b01}=-0.000109046 & \sigma^2_{b1}=0.02830715 \end{bmatrix}\bigg)\\ \end{align}\]
Despejando desde la definición de correlación de Pearson:
\(\sigma_{b10}=(0.04116649)(0.028030715)(-0.0945)=-0.000109046\)
Y tenemos entonces que el mejor modelo de los propuestos es el que tiene la comuna del establecimiento educativo como intercepto aleatorio y el tipo de establecimiento(público o privado) como pendiente aleatoria.