library(tidyverse)
library(dplyr)
library(tidyr)
library(tidymodels)
library(themis)
library(yardstick)
library(learningtower)

Introducción

Pruebas PISA: origen y definición

Desde el año 2000, la OCDE (Organización para la Cooperación y el Desarrollo Económico) lleva adelante la prueba PISA (Programme for International Student Assessment), que consiste en un conjunto de evaluaciones a estudiantes de 15 años, de diversos países, para medir su desempeño en Matemática, Lectura y Ciencia y poder evaluar la calidad de los sistemas educativos. Estas pruebas se hacen cada tres años, a excepción de períodos extraordinarios como la época de pandemia por Covid-19 que obligó a reprogramarlas. Argentina ha participado en todas las ediciones salvo en el año 2003, por elección, y en el 2015 que pese a haber realizado la prueba quedó fuera del resultado por un problema de muestreo. Para que el país pueda participar se le exige una rigurosa muestra de, al menos, 5000 estudiantes seleccionados teniendo en cuenta diversas zonas y tipos de escuela. Desde el año 2015 las pruebas comenzaron a tomarse mediadas por computadoras y no solo exigen a los estudiantes la resolución de tareas y ejercicios de áreas temáticas sino también que completen ciertos cuestionarios acerca de su entorno social y de vida, que sirve para contextualizar los resultados. Escalas de puntuación

Resulta cuanto menos complejo estandarizar procesos cognitivos, cuyos resultados se hacen evidentes a largo plazo, como la calidad educativa. Aun así, resulta necesario poder contar con ciertos indicadores de calidad educativa que permitan el diseño y la ejecución de políticas públicas acordes a las problemáticas de la coyuntura local. Las pruebas internacionales estandarizadas, como las PISA, permite no solo que los países puedan conocer cuál es el estado de la educación de sus poblaciones sino que los vuelve comparables, permitiendo tener un marco de análisis y poder diseñar políticas articuladas sin perder de vista los contextos locales. Las pruebas pisa tienen un puntaje final obtenido de la suma de las tres áreas, en una escala del 0 al 1000. La media de puntajes en las pruebas de 2022 fue de 476 (472 en Matemática, 476 en Lectura y 485 en Ciencias). Entre los países más destacados, que superaron esta media, Se encuentra Singapur con resultados de 543 puntos para Lectura y 561 para Ciencias. Por el otro lado, los páises con peor desempeño en la última edición son Rumania dentro de la Unión Europea, con una media de 428, y República Dominicana que se posicionó 80 de 81 participantes. Argentina, por su parte, no logró llegar a la media de OCDE y obtuvo una media de 378 puntos en Matemática, 401 en Lectura y 406 en Ciencias (Argentina, 2023).

Para el análisis comparativo los resultados obtenidos por los distintos países pueden ser agrupados en niveles, siendo el nivel I y II los dos mas bajos. El nivel II, alcanzado por el 70/80% de los estudiantes que participan es considerado el umbral de competencias mínimas. Por otro lado, los niveles más altos son el v y VI. En el nivel V se esperan resultados de Matemática entre 607 y 669 puntos; mientras que por encima de 669 ya sería un nivel VI. En Lectura, un nivel V espera resultados de entre 626 y 698 puntos. Por encima de estas cifras se entraría en nivel VI. Y, por último, en Ciencias lo esperable es una suma de entre 688 y 708 puntos en nivel V y mayor a 708 en nivel VI. Teniendo en cuenta los resultados de los países comentados previamente se observa que ninguno logra llegar a estos puntajes finales. En Estados Unidos tan solo el 7% de los estudiantes logran llegar a estos niveles (OECD, 2023).

Base de datos

Para poder analizar resultados de las pruebas PISA y otros indicadores asociados se usará el paquete learningtower desarrollado por Wang, Yacobellis, Siregar, Romanes, Fitter y Dalla Riva G en 2024. Puede accederse a información adicional sobre el mismo en https://github.com/kevinwang09/learningtower. Learningtower permite un acceso sencillo y rápido a un subconjunto de variables de datos PISA recogidos de la OCDE.

Contiene principalmente tres conjuntos de datos: student, school y countrycode. Cada uno de ellos arroja datos sobre datos y calificaciones de los estudiantes, datos de las escuelas y datos del país. Es una visión ordenada de los datos de la OCDE que también podrían ser descargados y analizados desde su propio sitio web.

Para los fines de este reporte se seleccionó solo el conjunto de datos student, que se explorará en profundidad a continuación mediante funciones como str(), summary() o head(). La función load_student(), según indica en el post de github de los desarrolladores, sirve para descargar los datos de interés.

students_2022 <- load_student(year = 2022)

En primer lugar se trata de una base que contiene 22 variables y que agrupa 613744 observaciones. A continuación una descripción de las variables:

La decisión de comenzar a indagar sobre variables socioeconómicas o demográficas responde a la necesidad de analizar el fenómeno del rendimiento educativo desde diferentes aristas, comprendiendo su complejidad y su naturaleza multifactorial. En ese marco, importa saber ¿cuánta incidencia tiene el contexto en el resultado de las prubeas PISA?. Para ello, y continuando con el trabajo realizado en Estadística computacional [https://rpubs.com/ldidier/1347662] acerca de la relación entre el índice de estatus económico y los resultados de las pruebas a nivel de cada estudiante, este análisis buscará, mediante el uso de machine learning predecir si los resultados de las pruebas de cada estudiante estarán por debajo o por encima de la media pero teniendo en cuenta solo las variables vinculadas al contexto de vida.

Preprocesamiento de datos

En primer lugar, a partir de la exploración de la base, observamos que la variable wealth para el año de estudio tomado (2022) no presenta información. Por ello se toma la decisión de sacarla del data set, para que los NA no alteren los demás procesos.

Además, la variable dishwasher es engañosa, puesto que en algunos países o regiones es de acceso común tener un lavavajillas pero en Argentina y otras regiones no, por lo que no es representativa de toda la heterogeneidad muestral.

Al comienzo se habían incluído las variables school_id y student_id; sin embargo, dieron problemas en el análisis ya que arrojaban un error de falta de datos. Es decir, el data de test tenía datos que no estaban en el data de train y por ello R no supo cómo interpretarlos. Si bien en otros trabajos es relevante el ID, en este caso no era relevante para lo que se está queriendo indagar y se procedió a quitarlos del data set. De este modo, el data set quedó conformado por 613744 observaciones y 18 variables.

students_2022 <- students_2022 %>% 
  select(-wealth, -dishwasher, -school_id, -student_id)

Los resultados de PISA, como se dijo previamente, son datos numéricos contínuos provenientes de las calificaciones en tres áreas: Matemática, Lectura y Ciencias. Para poder analizar si los estudiantes obtienen resultados por encima o por debajo de la media, en primer lugar fue necesario crear una nueva variable denominada total_score que sumara las tres áreas.

students_2022 <- students_2022 %>% 
  mutate(total_score = math+read+science)

Sobre ese total de calificación obtenido se calculó el promedio y se le llamó media_OCDE.

media_OCDE <- mean(students_2022$total_score, na.rm = TRUE)

La media_OCDE arrojó un resultado de 1330 puntos; es decir: los estudiantes cuyo total_score < 1330 estarán Por debajo de la media OCDE mientras que los estudiantes cuyo total_score > 1330 estarán en el grupo Por encima de la media. Ese clasificación de grupo es la variable target y es lo que se le pedirá al modelo de ML. Por eso, en primer lugar se debe crear la columna para la variable de clasificación.

students_2022 <- students_2022 %>%
 group_by(country) %>%
  mutate(grupo = if_else(total_score < media_OCDE,
                    "Por debajo de la media",
                    "Por encima de la media")) %>%
ungroup()

Modelo

Para poder predecir si los resultados obtenidos por cada estudiante estarán por debajo o por encima de la media de total_score se entrenará un modelo de decision tree para clasificación que tomará los datos de todas las variables y hará una predicción de ese puntaje para poder ubicar a cada estudiante en su grupo correspondiente. Al intentar correr el modelo por primer vez daba un resultado perfecto: la matriz de confusión no arrojaba valores falsos o erróneos, el accucary rondaba los 0.92 y la curva ROC estaba casi pegada al eje superior. Esto no podía ser normal, había una falla. Tras intentar diferentes cosas se llegó a la conclusión de que el error estaba en que se estaba entrenando un modelo para predecir pero se le estaban dando las respuestas en el entrenamiento. Al incorporar las variables math, read, science y total_score en la receta el modelo sabía exactamente qué puntaje había obtenido cada estudiante: no estaba haciendo predicciones sino solo ordenando los datos, algo que se podría hacer simplemente corriendo un summarise y que de hecho había sido realizado en el trabajo mencionado de Estadística Computacional. Eso se alejaba del objetivo del análisis que es predecir los resultados, y la clasificación en los grupos, a partir solo de las variables de contexto de vida. Por ello, en un segundo intento se quitaron de la receta esas variables que daban respuestas de puntajes.

División de datos

Una vez seleccionado el modelo de acuerdo a los intereses es necesario dividir los datos, para poder entrenar el modelo con un subconjunto y testearlo con otro. Por eso se dividieron en tres: inicial, entrenamiento y testeo, con una proporción se 0.75; es decir, 75% para entrenamiento y 25% para testeo. Además, se le pidió que la división fuera por grupo y, previamente, se relaizó el seteo de la semilla.

set.seed(123)
inicial_data <- initial_split(students_2022, prop=0.75, strata = grupo)
train_data <-  training(inicial_data)
test_data <- testing(inicial_data)

Una vez obtenidos los tres grupos se factorizó la variable grupo en los subconjuntos train y test.

train_data$grupo <- as.factor(train_data$grupo)
test_data$grupo <- as.factor(test_data$grupo)

Una vez que los datos han sido preprocesados, que se crearon las variables necesarias y se dividieron los datos se está en condiciones de plasmar la receta del modelo.

Receta

Como se comentó previamente no puede entrenarse al modelo con las respuestas, por ello, en la receta se coloca primero que la variable target será grupo y que lo hará a partir del subconjunto de datos train_data. Pero, además, se agrega la función step_rm() para indicarle que tome todas las variables del data salvo las indicadas entre paréntesis, a saber: math, read, science y total_score.

receta <- recipe(grupo ~ .,  data = train_data) %>% 
   step_rm(total_score, math, read, science)

Una vez armada la receta se procede a seleccionar el modelo de ML que será utilizado. Con la función decision_tree() se le indica que se trata de un árol de decisión y se le agrega, mediante set_engine() y set_mode() el paquete desde donde se toma el árbol y si será usado en modo clasificación o regresión.

arbol <- decision_tree() %>% 
  set_engine("rpart") %>% 
  set_mode("classification")

Workflow

Para no tener que repetir grandes líneas de código que entorpezcan el análisis se realiza un workflow donde se agregan la receta y el modelo de ML escogido, y todo queda condensado allí.

wf <- workflow() %>% 
  add_recipe(receta) %>% 
  add_model(arbol)

Hatsa ahora no se realizó el entrenamiento del modelo ni hay predicciones, solo se realizaron los ajustes necesarios sobre el data set para poder correr los modelos escogidos.

Entrenamiento

Con la función fit() se dará inicio al trabajo específico de ML y se comenzará a entrenar el modelo. Como se observa en el siguiente chunk, primero se setea la semilla, luego se crea un objeto de “árbol entrenado” que opera sobre el workflow configurado previamente y al cual le indico que el conjunto de datos sobre el que debe trabajar es el de train_data. Ya no es necesario aclarar cuál es la variable target o qué variables debe utilizar para predecir puesto que eso quedó almacenado en la receta que, a su vez, quedó almacenada en el workflow. Este paso es el que más demoras arroja puesto que está entrenando el modelo con el conjunto de datos provisto para realizar las predicciones que utilizaremos a continuación.

set.seed(123)
arbol_fit <- wf %>% 
  fit(data = train_data)
arbol_fit
## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_rm()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## n= 460307 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 460307 221100 Por debajo de la media (0.5196684 0.4803316)  
##    2) country=ALB,ARE,ARG,BGR,BRA,BRN,COL,CRI,DOM,GEO,GTM,IDN,JAM,JOR,KAZ,KHM,KSV,MAR,MDA,MEX,MKD,MNE,MNG,MYS,PAN,PER,PHL,PRY,PSE,QAT,QAZ,SAU,SLV,THA,URY,UZB 212359  55129 Por debajo de la media (0.7403972 0.2596028) *
##    3) country=AUS,AUT,BEL,CAN,CHE,CHL,CZE,DEU,DNK,ESP,EST,FIN,FRA,GBR,GRC,HKG,HRV,HUN,IRL,ISL,ISR,ITA,JPN,KOR,LTU,LVA,MAC,MLT,NLD,NOR,NZL,POL,PRT,QUR,ROU,SGP,SRB,SVK,SVN,SWE,TAP,TUR,USA,VNM 247948  81977 Por encima de la media (0.3306217 0.6693783)  
##      6) book=0-10,11-25,26-100 86434  41809 Por debajo de la media (0.5162899 0.4837101)  
##       12) book=0-10,11-25 45408  18148 Por debajo de la media (0.6003347 0.3996653) *
##       13) book=26-100 41026  17365 Por encima de la media (0.4232682 0.5767318) *
##      7) book=101-200,201-500,more than 500 161514  37352 Por encima de la media (0.2312617 0.7687383) *

Predicciones

En primer lugar se generan las clases predichas por el modelo. Para cada observación del conjunto de test, que el modelo hasta ahora no había visto, se obtiene una clase predicha. Es decir, se ordenan los datos en la etiqueta ganadora en cada caso y el resultado son muchos datos en cada una de las dos etiquetas posibles de la variable target con la que se está trabajando (Por encima o por debajo de la media).

Sobre esas etiquetas que asignó calcula las probabilidades de que la clase sea la correcta. El modelo no solo indica a qué clase puede pertenecer cada observación sino también con qué nivel de confianza estadística lo hace. Estas probabilidades, que permiten calcular las posibilidades de error, serán utilizadas más adelante para la Curva ROC.

clase_pre <- predict(arbol_fit, new_data = test_data, type = "class")
clase_prob <- predict(arbol_fit, new_data = test_data, type = "prob")

Resultados

resultados_test <- test_data %>% 
  select(grupo) %>%
  bind_cols(clase_pre) %>%
  bind_cols(clase_prob)

Una vez obtenidas esas predicciones, se las asigna a un objeto llamado resultados_test mediante la función bind_cols().

Métricas

Para saber si el modelo funciona, si sus resultados son confiables, existen métricas a las cuales se les debe prestar atención, tales como: - Accuracy: mide la proporción de predicciones correctas. - Sensivity: Mide la cantidad de casos positivos detectados correctamente. - Specificity: Mide la cantidad de casos negativos detectados correctamente. - Precision: Mide cuántos casos de los que predijo como positivo efectivamente lo eran. - F_meas: Balance entre sensivity y precision.

Mediante la función metric_set() se arma la lista de las métricas a solicitar y, luego, se las aplica a los resultados del modelo para evaluar su rendimiento.

metricas <- metric_set(accuracy, sensitivity, specificity, precision, f_meas)

metricas(resultados_test,
                 truth   = grupo,
                 estimate = .pred_class)
## # A tibble: 5 × 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 accuracy    binary         0.722
## 2 sensitivity binary         0.772
## 3 specificity binary         0.669
## 4 precision   binary         0.716
## 5 f_meas      binary         0.743

¿Qué indican los valores obtenidos en cada métrica? En primer lugar, un accuracy de 0.72 indica que el modelo logra clasificar correctamente el 72% de los casos dados. Previamente se habló de un valor anormal de 0.92 que se había obtenido al incluir todas las variables en la receta. Como se le daban las respuestas al modelo de antemano el margen de error fue muy pequeño. Sin embargo, cuando se dejaron solo las variables vinculadas a contexto de vida y se quitaron las que otorgaban algún indicio sobre puntaje, este valor bajó. Si bien es un valor aceptable hay que tener cuidado con el desbalance de las clases. Para seguir corroborando el rendimiento tenemos sensivity y specificity. La primera significa que cuando la clase positiva (Por encima de la media) se presenta el modelo es capaz de acertar 77% de los casos. Por otro lado, el 0.67 de specificity indica que el modelo es capaz de reconocer 67 de cada 100 casos Por debajo de la media, lo que evita los falsos negativos. La precision indica que de todas las clasificaciones que el modelo predice hay un 72% que son correctas. Por último, f_meas sirve para los casos en los cuales las clases están desbalanceadas (hay muchos más casos en una clase que en otra) y arroja un balance entre sensivity y precision. Que su valor sea de 0.74 indica que ambas métricas tenían un valor alto y equilibrado. Es útil sobre todo para detectar correctamente los positivos. En este estudio, si se dice que las con determinadas variables socioeconómicas el puntaje puede estar por encima de la media es necesario que este dato sea lo más fidedigno posible.

Matriz de confusión

Para poder visualizar de manera más clara cuántos valores verdaderamente positivos y verdaderamente negativos, tanto como los falsos de ambas clases, quedaron a partir del modelo, se realiza una matriz de confusión y luego se le agega un formato cuadrícula para una interpretación más clara, mediante las funciones conf_mat() y autoplot().

conf_mat(resultados_test,
         truth   = grupo,
         estimate = .pred_class)
##                         Truth
## Prediction               Por debajo de la media Por encima de la media
##   Por debajo de la media                  61549                  24415
##   Por encima de la media                  18187                  49286
conf_mat(resultados_test,
         truth   = grupo,
         estimate = .pred_class) %>%  
autoplot(type = "heatmap")

Curva ROC

Por último, se grafica la curva ROC, que muestra que tan bien separa las clases el modelo evaluando todos los posibles cortos, en un continuo. Lo que marca la curva es cómo cambia la detección de positivos y de falsos positivos a distintos umbrales de clasificación y lo hace colocando en el eje Y sensibility y en el eje X specificity. Antes de graficar, se calcula roc_auc que es el área que queda bajo la curva; es decir, la probabilidad de que el modelo asigne un puntaje más alto a un caso positivo que a un caso negativo tomado al azar. Si, por ejemplo, se toma un caso “Por encima” y otro “Por debajo” hay un 74% de probabilidades de que se asigne correctamente la categoría o etiqueta positiva. Es interesante detenerse aquí porque, inicialmente, no se le indicó a roc_auc la categoría “correcta” y todos los valores daban alarmantemente mal. El roc_auc era menos de 0.5 es decir que el modelo estaba ordenando los casos peor que el azar y la curva estaba graficada con la panza hacia abajo de la línea gris. Esto llamo la atneicón puesto que las métricas analizadas previamente indicaban que el rendimiento del modelo era bueno. Parecía que el modelo etiquetaba bien pero ordenaba mal las probabilidades. Al ver eso se indicó que roc_auc tomara la etiqueta “Por debajo de la media” puesto que previamente estaba calculando utilizando como clase positiva la probabilidad asociada a la categoría equivocada.

Al modificar la etiqueta el roc_auc alcanzó un valor de 0.74, super aceptable (los rangos son de 0a 1) y la curva “se dio vuelta” mostrando una curvatura hacia arriba de la línea de puntos. Si solo se hubiera analizado el roc_auc o se hubiera graficado la curva se podría haber llegado a la conclusión, errónea, de que el modelo no funcionaba. Sin embargo, al contar también con las otras métricas se constató que el modelo tenía un buen rendimiento pero que al asignarle incorrectamente la etiqueta se estaban ordenando mal los datos. Esto destaca la importancia de realizar todos los pasos para poder tener la mayor cantidad de certezas.

roc_auc(resultados_test, truth = grupo, `.pred_Por debajo de la media`)
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc binary         0.745
#roc_auc(resultados_test, truth = grupo, `.pred_Por encima de la media`)
resultados_test %>%
  roc_curve(truth = grupo, `.pred_Por debajo de la media`) %>%
  ggplot(aes(x = 1 - specificity, y = sensitivity)) +
  geom_line(size = 1.5, color = "midnightblue") +
  geom_abline(lty = 2, alpha = 0.5, color = "black", size = 1.2) +
  theme_minimal()+
  coord_equal()

Conclusión

Een trabajos previos de Estadística computacional interesaba conocer si el Índice ESCS tenía incidencia sobre los resultados obtenidos por los estudiantes que rindieron la prueba PISA en 2022. Tras diveros cálculos, fue posible rechazar la hipótesis de que no hay diferencia entre aquellos estudiantes con bajo nivel de índice ESCS, sin computadora y sin internet y aquellos con alto índice ESCS, internet y computadora. Es decir, a mayor índice ESCS, teniendo computadora y teniendo conexión a internet el puntaje total mejora.

Para este trabajo, y partiendo de esos resultados, interesaba saber si a partir de las variables vinculadas solo a contexto de vida el modelo de ML era capaz de predecir si el puntaje total obtenido estaría por encima o por debajo de la media. ¿Para qué serviría esto? en Sociología de la Educación y Políticas Educativas es siempre un tema de discusión cuánta incidencia tiene el contexto de vida sobre el rendimiento educativo. Factores como nivel socioeconómico, educación de madre y padre, condiciones habitacionales e inclusión digital, están sobre la mesa cuando se discute el rendimiento. Conocer si el árbol de decisión es capaz de solo con estas variables predecir categorías podría arrojar valiosa información sobre cuánta incidencia tiene el contexto y focalizar hacia allí el diseño y la implementación de políticas públicas, basadas en datos, para la mejora de los resultados.

En este caso, el modelo obtuvo un buen rendimiento. Predijo de manera correcta el 72% de los casos positivos. Esto no quiere decir que no haya errores ni quiere decir que el contexto sea determinante en el resultado de PISA, pero si es un valioso puntapié inicial para incorporar el aprendizaje maquínico a fenómenos multivariables como el rendimiento educativo.

Para seguir complejizando el análisis se podría realizar un tunning de hiperparámetros o correr modelos de ML más sofisticados y complejos como random forest o boosting.