library(tidyverse)
library(dplyr)
library(tidyr)
library(pwr)
library(data.table)
library(gt)
library(learningtower)
library(weights)
library(lme4)
library(nortest)
library(lmtest)
library(sandwich)
library(patchwork)  #Interesante paquete que encontré gracias a Estación R y que permite hacer un dashboard con gráficos

Acceso Digital y pruebas PISA: ¿la brecha digital es también educativa?

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).

Learningtower: explorando PISA

students_2022 <- load_student(year = 2022)

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 mi interés.

Student_subset

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:

year: variable entera del año en que se tomó la prueba PISA.
country: Variable de tipo factor que asigna un código de tres letras a cada país. Ejemplo: Islandia = ISL a Argentina= ARG.
school_id y student_id: indicadores únicos para escuela y estudiante. Como estoy utilizando el conjunto de datos de estudiantes todos los que pertenezcan a la misma escuela tendrán el mismo school_id. Son variables numéricas (enteras).
mother_educ y father_educ: refiere al nivel educativo de padre y madre (más alto alcanzado). Son variables de tipo factor cuyo nivel mínimo es ISCED1 y el máximo es ISCED 3A.
gender: variable de facotr que determina si el género es masculino o femenino. No se contemplan en esta base de datos otras opciones.
computer e internet también son variables de tipo factor que determinan si el estudiante tiene acceso, o no, a una computadora o a internet. Las dos opciones posibles son yes o no.
math, read y science: son los puntajes obtenidos en cada área por cada estudiantes. Se trata de variables numéricas.
stu_wgt: también es una variable numérica pero del peso muestral del estudiante, utilizado para garantizar que los análisis sean representativos del total de estudiantes de 15 años en un país determinado.
desk, room y diswasher: se trata de ariables tipo factor, con solo dos respuestas posibles yes o no acerca de la poseción de escritorio, habitación privada y lavavajillas. Sobre el lavavajillas me pregunto si no será una variable proxy puesto que, por ejemplo en Argentina, las clases medias y altas solamente tienen acceso a un lavavajillas.
television, computer_n, car y book: variables de tipo facotr, con respuestas posibles de 1, 2, 3 o +3 sobre la tenencia de televisión, computadoras, libros y autos.
wealth: variable numérica que indica la riqueza familiar material. s un índice estandarizado (con media 0 y desviación estándar 1 a nivel OCDE). Valores positivos → hogares con más recursos materiales que el promedio OCDE. Valores negativos → hogares con menos recursos materiales que el promedio. Esta variable está ausente en el dataset seleccionado, es decir, no contiene observaciones (todas NA).
escs: índice de estatus económico, social y cultural. Variable numérica. Es una variable combinada, construida por la OCDE. La misma combina: 1) Nivel educativo más alto de los padres/tutores (años de estudio). 2) Ocupación más alta de los padres/tutores (clasificada en un índice internacional de estatus socioeconómico) y 3) Bienes y recursos en el hogar, que incluyen: Riqueza material (índice wealth); Recursos educativos en casa (ej. cantidad de libros, escritorio, material de estudio) y Bienes culturales (ej. obras de arte, literatura clásica).

Brecha digital y calidad educativa ¿relacionadas?

Mucho se ha hablado ya de las posibles relaciones entre tener acceso a las tecnologías digitales y posibles mejoras en los índices de calidad educativa. Políticas públicas de gran envergadura como el Programa Conectar Igualdad en Argentina, inspirado en el One laptop per child de Nicholas Negroponte, apuntaba justamente a ello: reducir la brecha digital como una apuesta por la inclusión social que podría relacionarse con una mejora en la equidad y calidad educativa.

La escuela, desde fines de silo pasado cuando irrumpe la revolución digital, oscila entre mirar a la tecnología con buenos ojos y ‘dejarla entrar’ y normativas regulacionistas o prohibitistas como las que están teniendo auge en la actualidad, vinculadas a la prohibición de que los estudiantes asistan a clase con sus dispositivos móviles.

A los fines de este análisis no entraremos en debate sobrecuál deberia ser la postura de la escuela sino que interesa Analizar si existe una relación entre el acceso digital (tener o no computadora o internet) y el índice de estatus económico, social y cultural en relación a los resultados finales de las pruebas PISA. Interesa ver estas relaciones pero estableciendo una comparativa entre los países cuyas notas finales superan la media de OCDE y aquellos que no alcanzan esta puntuación, para ver si el acceso digital y el índice escs resultan más significativo en uno u otro caso. Para llevar a cabo este análisis se trabajará, en primera instancia, con las variables total_score, como variable dependiente (y), la cual voy a crear nueva en mi dataset a partir de la sumatoria entre math, read y science, grupo que es una variable de tipo factor con dos niveles: por debajo o por encima de la media de OCDE. Por otro lado, las variables independientes, predictoras o explicativas son: computer e internet, variables también de tipo factor que en el dataset original aparecen como yes/no; y, escs variable numérica contínua creada por OCDE. En los análisis voy a pedirle a R que me incluya la variable stu_wgt que es el peso muestral y que puede ayudarme a evitar sesgos en mi análisis.

La hipótesis planteada en este trabajo es que: si los estudiantes cuentan con computadora e internet y tienen un elevado ESCS su puntuación final en la prueba será mayor que la de aquellos estudiantes sin estos accesos.

Análisis de variables

Como no tengo una columna específica que sume las tres áreas de evaluación de las pruebas la voy a crear.

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

Luego, voy a calcular la media de esa nueva columna llamada total_score:

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

Ahora agrupo los países por si sus total_score están por encima o por debajo de la media_OCDE.

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()

Países por debajo de la media

En este primer momento del reporte solo se trabajará con la muestra recortada que corresponde a los países cuyo total_score están por debajo de la media de OCDE. Interesa comenzar trabajando por este grupo particular porque partimos del supuesto de que la educación en general es una variable que puede estar influenciada por múltiples factores individuales y contextuales, tales como: edad, género, condiciones material, nivel de alfabetización de padres, etcétera. En este sentido, interesa ver cómo se comportan estas relaciones en países que tienen promedios por debajo de los esperados y sobre los cuáles se podrían diseñar programas o políticas educativas específicas.

paises_inf <- students_2022 %>% 
  filter(grupo== "Por debajo de la media")

Para comenzar a explorar este nuevo dataset veamos primero a qué países incluye para ver si ‘tiene sentido’ en relación con datos que conocemos de antemano y que están en la introducción de este reporte. Tal como se comentó previamente Rumania, República Dominicana y Argentina deberían haber quedado en este grupo, ya que vimos que sus promedios están por debajo de la media. Pero, por el contrario, Singapur debería quedar en el otro grupo siendo el país con mejor desempeño.

unique(paises_inf$country)
##  [1] ALB QAZ ARG AUS AUT BEL BRA BRN BGR KHM CAN CHL TAP COL CRI HRV CZE DNK DOM
## [20] SLV EST FIN FRA GEO PSE DEU GRC GTM HKG HUN ISL IDN IRL ISR ITA KSV JAM JPN
## [39] KAZ JOR KOR LVA LTU MAC MYS MLT MEX MNG MDA MNE MAR NLD NZL NOR PAN PRY PER
## [58] PHL POL PRT QAT ROU SAU SRB SGP SVK VNM SVN ESP SWE CHE THA ARE TUR QUR MKD
## [77] GBR USA URY UZB
## 80 Levels: ALB ARE ARG AUS AUT BEL BGR BRA BRN CAN CHE CHL COL CRI CZE ... VNM

Tal como suponíamos vemos los códigos DOM (República Dominicana), ROU (Rumania) y ARG (Argentina. Y no vemos SGP (Singapur). En principio este recorte del dataset parecería estar concordando con los datos que se conocían previamente. Profundizemos un poco:

summary(paises_inf)
##       year         country        school_id           student_id      
##  Min.   :2022   ARE    : 13918   Length:318943      Min.   :  800001  
##  1st Qu.:2022   KAZ    : 11820   Class :character   1st Qu.:21405336  
##  Median :2022   IDN    : 11011   Mode  :character   Median :39816862  
##  Mean   :2022   ESP    :  9431                      Mean   :42889129  
##  3rd Qu.:2022   ARG    :  8312                      3rd Qu.:64201872  
##  Max.   :2022   BRA    :  7651                      Max.   :86007492  
##                 (Other):256800                                        
##            mother_educ               father_educ        gender      
##  ISCED 1         : 25789   ISCED 1         : 25855   female:154249  
##  ISCED 2         : 57215   ISCED 2         : 58044   male  :164678  
##  ISCED 3A        :145741   ISCED 3A        :133870   NA's  :    16  
##  ISCED 3B, C     : 49173   ISCED 3B, C     : 53183                  
##  less than ISCED1: 16796   less than ISCED1: 16238                  
##  NA's            : 24229   NA's            : 31753                  
##                                                                     
##  computer      internet           math            read          science     
##  no  : 84639   no  : 48737   Min.   :  0.0   Min.   :  0.0   Min.   :  0.0  
##  yes :213010   yes :247998   1st Qu.:326.6   1st Qu.:311.6   1st Qu.:330.7  
##  NA's: 21294   NA's: 22208   Median :367.2   Median :361.0   Median :374.8  
##                              Mean   :364.3   Mean   :356.7   Mean   :370.7  
##                              3rd Qu.:404.9   3rd Qu.:405.6   3rd Qu.:414.5  
##                              Max.   :589.3   Max.   :593.4   Max.   :635.8  
##                                                                             
##     stu_wgt          desk          room        dishwasher    television   
##  Min.   :   1.00   NA's:318943   no  : 74257   NA's:318943   0   : 11248  
##  1st Qu.:   5.10                 yes :223285                 1   :187495  
##  Median :  13.32                 NA's: 21401                 2   : 72440  
##  Mean   :  50.94                                             3+  : 15350  
##  3rd Qu.:  48.96                                             NA's: 32410  
##  Max.   :2371.65                                                          
##                                                                           
##  computer_n      car                    book           wealth      
##  0   :133389   0   : 98359   0-10         :39554   Min.   : NA     
##  1   :119406   1   :101303   101-200      :60235   1st Qu.: NA     
##  2   : 17928   2   : 59346   11-25        :98094   Median : NA     
##  3+  :  5067   3+  : 40352   201-500      :20790   Mean   :NaN     
##  NA's: 43153   NA's: 19583   26-100       :69607   3rd Qu.: NA     
##                              more than 500:10858   Max.   : NA     
##                              NA's         :19805   NA's   :318943  
##       escs          total_score        grupo          
##  Min.   :-6.8407   Min.   : 212.4   Length:318943     
##  1st Qu.:-1.4888   1st Qu.: 983.5   Class :character  
##  Median :-0.6891   Median :1109.5   Mode  :character  
##  Mean   :-0.7464   Mean   :1091.6                     
##  3rd Qu.: 0.1125   3rd Qu.:1220.1                     
##  Max.   : 7.3800   Max.   :1329.6                     
##  NA's   :16613

Se trata de un dataset de 318943 observaciones en 24 variables, de las cuales total_score y grupo fueron creadas para este reporte.

Un dato interesante es que la muestra registra valores menores para mujeres (154249) como para varones (164678) cuyos padres y madres habrían alcanzado estudios secundarios (ISCED 3A) pero donde muy pocos han llegado a niveles de educación terciaria o universitaria. En relación al índice de estatus (escs), y directamente relacionado con el nivel educativo familiar, hay un promedio y una mediana negativos, lo cual indicaría que los casos registrados de estos países tienen un bajo nivel de estatus.

En relación a los puntajes obtenidos em las pruebas PISA de 2022, el promedio de este grupo de países es de 1091 puntos, y la mediana arroja un valor de 1109. Mientras que el valor usado de referencia (media de OCDE) es de 1329 puntos. Sin embargo, estos promedios entre los países del grupo inferior oscilan entre los 212 puntos y los 1329. Veamos algunos de estos datos de manera gráfica:

#Filtrar los NA
#Sacar proporciones
compu_prop <- paises_inf %>%
  filter(!is.na(computer)) %>%
  count(computer) %>%
  mutate(prop = n / sum(n))
#Graficar
graf_compu <- ggplot(compu_prop, aes(x = computer, y = prop, fill = computer)) +
  geom_col(width = 0.4) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("#8B3E2F", "darkolivegreen4", "gray")) +
  labs(title = "Estudiantes con o sin computadora",
       subtitle = "Pruebas PISA, 2022",
       caption = "Figura 1: Elaboración propia con datos de OCDE (2023)",
       x = "Acceso a computadora",
       y = "Porcentaje") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5),
        axis.text.x = element_blank())
graf_compu

#Filtrar los NA
#Sacar proporciones
inet_prop <- paises_inf %>%
  filter(!is.na(internet)) %>%
  count(internet) %>%
  mutate(prop = n / sum(n))
#Graficar
graf_inet <- ggplot(inet_prop, aes(x = internet, y = prop, fill = internet)) +
  geom_col(width = 0.4) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("#8B3E2F", "darkolivegreen4", "gray")) +
  labs(title = "Estudiantes con o sin internet",
       subtitle = "Pruebas PISA, 2022",
       caption = "Figura 2: Elaboración propia con datos de OCDE (2023)",
       x = "Acceso a internet",
       y = "Porcentaje") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5),
        axis.text.x = element_blank())
graf_inet

graf_compu+graf_inet

Tal como puede observarse en los gráficos previos (Figuras 1 y 2) es mayor la cantidad de estudiantes que poseen una computadora e internet que aquellos que no poseen.

Análisis de correlación

Lo primero que interesa conocer, para responder el objetivo propuesto, es si hay relación entre las variables seleccionadas. Para ello, se utilizará la función cor() para calcular el coeficiente de correlación. Pero, primero, debemos transformar las respuestas de computer y de internet en variables numéricas, por lo que en vez de yes o no serán 1=yes y 0=no. Además, voy a eliminar las filas que arrojan valores nulos para las variables que debo analizar porque si no R no podrá procesarlo:

paises_inf <- paises_inf %>% 
  mutate(computer_num = ifelse(computer== "yes", 1, 0))
paises_inf <- paises_inf %>% 
  mutate(internet_num = ifelse(internet== "yes", 1, 0))
paises_inf <- paises_inf %>% 
  drop_na(stu_wgt, computer_num, internet_num, escs)

Tal como se observa el dataset de nuestra muestro sumó dos nuevas columnas, que son computer_num e internet_num y que corresponden a variables de tipo factor de dos niveles. Ahora si pueden hacerse los tres cálculos de correlación.

cor_escs <- cor(paises_inf$total_score, paises_inf$escs)  #Correlaciones sin peso muestral
cor_comp <- cor(paises_inf$total_score, paises_inf$computer_num)
cor_inet <- cor(paises_inf$total_score, paises_inf$internet_num)

Los coeficientes de correlación calculados previamente no consideran el peso muestral, que sirve para ayudar a evitar sesgos derivados de problemas de elegibilidad y representatividad. Por eso, a continuación se calculan también con el paquete weiths que sirve para indicarle a R que tenga en cuenta este peso en el cálculo.

cor_escs_wgt <- wtd.cor(x= paises_inf$escs, y= paises_inf$total_score, weight = paises_inf$stu_wgt)
cor_comp_wgt <- wtd.cor(x= paises_inf$computer_num, y= paises_inf$total_score, weight = paises_inf$stu_wgt)
cor_inet_wgt <- wtd.cor(x= paises_inf$internet_num, y= paises_inf$total_score, weight = paises_inf$stu_wgt)

A continuación se analizan los valores de estos coeficientes. En primer lugar se observa que los coeficientes de correlación ponderados, es decir con el peso muestral agregado, son menores que aquellos calculados sin peso muestral. En el caso de ESCS con peso arroja un valor de 0.19 mientras que sin peso es 0.20. En el caso de computer, con peso es 0.16 y sin peso 0.14 y en el caso de internet con peso es 0.17 y sin peso 0.17. En primer lugar estamos hablando de que si existe una relación entre las variables, puesto que ningún coeficiente dio cero. Además, todas las relaciones son positivas (si una variable aumenta la otra podría aumentar también) ya que ningún valor es negativo. El hecho de que los coeficientes que involucran peso muestral sean levemente más bajos se debe a que esta fórmula ‘ajusta’ los datos en el caso de que algunos estudiantes puedan estar sobre o subrepresentados. Que los valores ponderados sean menores también indica que las correlaciones son más débiles de lo que parece ser cuando no se incluyen los pesos muestrales. Las variables que más relación parecen tener con el puntaje promedio son el índice ESCS e internet, mientras que la computadora está un poco por debajo. Estos valores, aunque bajos, demuestran que existe una relación moderada. En estudios vinculados a la educación, como PISA, esto puede deberse a la multiplicidad de factores que podrían incidir en el desempeño de los estudiantes. Más adelante se hará un Análisis de Regresión Múltiple para ver si las variables que por si solas resultan en una relación moderada o débil combinadas con otras permiten mejores explicaciones de la variabilidad del puntaje.

Para tener una primera aproximación de cómo se comportan nuestras variables, de manera aislada, previo a un modelo de regresión, acontinuación están graficadas la variable dependiente total_score y la variable independiente (x1) escs. Con estos gráficos puede evaluarse la distribución de los datos:

histo_total_score <- ggplot(data = paises_inf, aes(x= total_score))+
  geom_histogram(aes(y = after_stat(density)),color= "#8B3E2F", fill= "lightsalmon2")+
  geom_density(color = "gray35", linewidth = 1.5, alpha = 0.3)+
   labs(
    title = "Histograma de total_score",
    subtitle = "Variable dependiente",
    caption = "Figura 3: elaboración propia con datos de OCDE (2023)")+
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))
histo_total_score

histo_escs <-  ggplot(data = paises_inf, aes(x= escs))+
  geom_histogram(aes(y = after_stat(density)),color= "#8B3E2F", fill= "lightsalmon2")+
  geom_density(color = "gray35", linewidth = 1.5, alpha = 0.3)+
   labs(
    title = "Histograma de ESCS",
    subtitle = "Variable independiente (X1)",
    caption = "Figura 4: Elaboración propia con datos de OCDE (2023)")+
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))
histo_escs

histo_total_score+histo_escs

Total score, computadora e internet: ¿relacionadas?

Para poder verificar la hipótesis planteada en un inicio, que indica que los estudiantes que poseen computadora, conexión a internet y mayor índice de ESCS obtienen puntajes más altos en las prubeas PISA, y debido a que las variables computer e internet son de tipo factor con dos niveles, se procederá a realizar un Test de medias. Este test sirve para saber si la media del grupo 0 es igual que la media del grupo 1 y poder verificar relaciones entre variables. En este caso, si la media de estudiantes que NO poseen computadora es igual a la media del grupo que SI posee, en cuanto a sus puntajes obtenidos.

\[\mu_{(con\;compu)} = \mu_{(sin\;compu)}\]

\[\mu_{(con\;compu)} \neq \mu_{(sin\;compu)}\]

Promedio relacionado con tener o no computadora

t_test_comp <- t.test(total_score ~ computer_num, data = paises_inf)
t_test_comp
## 
##  Welch Two Sample t-test
## 
## data:  total_score by computer_num
## t = -91.212, df = 149839, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -57.94381 -55.50598
## sample estimates:
## mean in group 0 mean in group 1 
##        1057.735        1114.460

El Test de medias realizado demuestra que se rechaza la hipótesis nula (que el puntaje de los estudiantes no varía en relación a tener o no computadora) ya que el p-value cercano a cero nos indica que hay relación entre tener o no computadora y el puntaje total final obtenido en las pruebas PISA. Pero, además, pueden observarse otros datos: en primer lugar hay un t-estadístico cuyo valor es de -91 que está en negativo puesto que quienes no tienen computadora obtienen un puntaje menor con respecto a aquellos que si tienen. Por otro lado, el intervalo de confianza -58/-55 indica que en el 95% de los estudios que se realicen, el estimador de las medias estará entre esos valores. Por último, el promedio del grupo 0, es decir quienes NO tienen computadora, es de 1058 puntos en el examen mientras que el promedio del grupo 1, quienes SI tienen computadora, es, en promedio, de 1114 puntos, aún por debajo claramente de la media OCDE. Observemos gráficamente esta distribución:

ggplot(paises_inf, aes(x = factor(computer_num), y = total_score, fill = factor(computer_num))) +
  geom_boxplot(width = 0.5, outlier.alpha = 0.3) +  
  scale_x_discrete(labels = c("No tiene", "Sí tiene")) +
  scale_fill_manual(values= c("#BF3EFF", "darkorchid4"))+
  labs(title = "Distribución  de puntaje según tenga computadora",
       subtitle = "Pruebas PISA 2022",
       x = "Computadora",
    y = "Total score",
    caption = "Figura 5: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5),
        legend.position = "none")

Puntajes relacionados con tener o no conexión a internet

t_test_inet <- t.test(total_score ~ internet_num, data = paises_inf)
t_test_inet
## 
##  Welch Two Sample t-test
## 
## data:  total_score by internet_num
## t = -91.667, df = 66772, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -71.71054 -68.70816
## sample estimates:
## mean in group 0 mean in group 1 
##        1039.612        1109.822

Al igual que en el caso de la computadora, el hecho de tener internet se relaciona con el puntaje total. Se rechaza la hipótesis alternativa de que la diferencia entre medias es igual a cero. El p-value por debajo del 0.05 indica que hay una relación entre las variables y los intervalos de confianza permiten afirmar que, en un 95% de casos, los estimadores de las medias podrían estar entre los -72 y -69 puntos. Por último, el puntaje de los estudiantes sin internet podría rondar los 1039 puntos mientras que en el caso de aquellos que si tienen podría ascender a 1110 puntos. Si bien las diferencias no son tan grandes son significativas y permiten rechazar H0.

Veamos también gráficamente estas diferencias:

ggplot(paises_inf, aes(x = factor(internet_num), y= total_score, fill = factor(internet_num)))+
  geom_boxplot(width = 0.5, outlier.alpha = 0.3) +  
  scale_x_discrete(labels = c("No tiene", "Sí tiene")) +
  scale_fill_manual(values= c("#1E90FF", "#104E8B"))+
  labs(title = "Distribución de puntaje según tenga internet",
       subtitle = "Pruebas PISA 2022",
       x = "Internet",
    y = "Total Score",
    fill = "Internet",
    caption = "Figura 6: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5),
        legend.position = "none")

Análisis de regresión lineal

Se sabe que la brecha digital está directamente vinculada con la brecha económica y social. La desigualdad estructural entre ciertos sectores se traduce también, en la actualidad, en una desigualdad digital. Diversos programas y políticas, tanto de organismos internacionales como la ONU y UNICEF, como locales, por ejemplo el Programa Conectar Igualdad o la Red Federal de Fibra Óptica en Argentina, han buscado disminuir esta brecha, reduciendo las barreras de acceso tanto a los equipamientos como a la conectividad.

Teniendo en cuenta esto a continuación se realizará un Análisis de regresión lineal entre el promedio OCDE del puntaje obtenido por cada estudiante en las pruebas del 2022, como variable dependiente, y el índice de estatus social, económico y cultural definido por OCDE, como variable explicativa.

options(scipen = 999)
modelo_1 <- lm(formula = total_score ~ escs, 
               data = paises_inf,
               weights = stu_wgt)
summary(modelo_1)
## 
## Call:
## lm(formula = total_score ~ escs, data = paises_inf, weights = stu_wgt)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -20202.0   -331.9     61.7    423.2  10572.8 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 1114.7887     0.3921  2843.2 <0.0000000000000002 ***
## escs          25.2911     0.2397   105.5 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1086 on 290324 degrees of freedom
## Multiple R-squared:  0.03694,    Adjusted R-squared:  0.03694 
## F-statistic: 1.114e+04 on 1 and 290324 DF,  p-value: < 0.00000000000000022

A través del análisis se podría rechazar la hipótesis nula de que no hay diferencia en los estudiantes con mayor o menor índice de estatus en relación a sus puntajes totales, puesto que el valor del p-value es prácticamente cero. Cabe destacar que el índice ESCS creado por OCDE resulta de una combinación de otras variables, como fuera explicado anteriormente, por lo cual no puede simplemente medirse en unidades sino en aumento de desviaciones estándar; por ejemplo: el aumento de 1 unidad de escs significaría el aumento de 1 desviación estándar. Aclarado esto, Alpha muestra el puntaje promedio de los estudiantes cuando el índice es cercano a cero (puesto que sería un sinsentido pensar que alguien tenga ‘cero estatus’), que arroja un valor de 1119. Por su parte, Beta indica cuál es la pendiente, es decir cuánto podría aumentar el puntaje total de los estudiantes si aumentara una unidad (una desviación) el índice de escs, arrojando un valor de 25.29. Es decir, los estudiantes podrían sacarse, en promedio, 25.29 puntos más en su examen si aumenta una unidad su escs. Se observa un error estándar residual de 1086 que se verá con mayor detalle en el gráfico de residuos.

IC_escs_inf <- 25.29-1.96*0.23
IC_escs_sup <- 25.29+1.96*0.23

Puede afirmarse, con un 95% de confianza, que los estudiantes cuyo escs aumente en una unidad podrían tener un aumento en sus puntajes totales de entre 25 a 26 puntos. La reducida amplitud entre estos valores podría indicar que la estimación realizada resulta precisa.

Por último, el valor de R2 de 0.03 indica que el modelo utilizado permite explicar el 3% de la relación, es decir: la variable escs explica la variación de los puntajes en un 3%. Esto indica a simple vista que NO nos alcanza solo con analizar la relación entre estas dos variables, puesto que el porcentaje que este modelo logra explicar acerca de cómo varían los promedios de puntaje es BAJO. Por eso, debemos seguir sumando variables, siempre teniendo en cuenta cuánto nos está ‘penalizando’ el modelo con el R2 adjusted. Pero, antes de eso, grafiquemos esta información:

ggplot(paises_inf, aes(x = escs, y = total_score)) +
  geom_point(color = "#8B3A3A") +
  geom_smooth(method = "lm", color = "gray2", se = TRUE) +
  labs(
    title = "Relación entre ESCS y  puntajes obtenidos",
    subtitle = "Pruebas PISA 2022",
    caption = "Figura 7: Elaboración propia con datos de OCDE (2023)",
    x = "Índice ESCS",
    y = "Total score") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

En el gráfico 7 pueden observarse datos agrupados alrededor de la línea ‘negra’ de ajuste, esto podría indicar que hay una relación entre ESCS y Total_score pero que el índice por si solo no estaría logrando explicar todos los puntajes: seguramente hay diferentes puntajes con el mismo índice de ESCE ¿entonces qué los explica? lo podremos ver al correr los modelos siguientes.

residuos <- data.frame(residuos= modelo_1$residuals, ajustados = modelo_1$fitted.values)
ggplot(data = residuos, aes(x= ajustados, y= residuos))+
  geom_point(color= "gray")+
  geom_smooth(se = FALSE, color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(
    title = "Residuos del modelo",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 8: elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

Tal como podría esperarse con los resultados del summary del modelo el mismo presenta problemas de linealidad: la línea roja que representa a los valores ajustados se curva hacia abajo, razón por la cual podría no estar cumpliéndose el supuesto de que los datos responden a un modelo lineal.

Previamente se dijo que este modelo presenta un R2 muy bajo, o sea que la variable independiente escogida resulta poco explicativa de la variable dependiente. Por este motivo vamos a correr un modelo de regresión múltiple, que combine otras variables independientes y algunas variables de control.

Análisis de regresión lineal múltiple

Se explicó previamente que variables vinculadas a la educación, como el promedio de puntaje obtenido por cada estudiante en las pruebas PISA 2022, pueden estar vinculadas a una multiplicidad de factores, tanto individuales como contextuales. Al respecto, se correrá un modelo múltiple que al modelo original (puntaje total en relación al índice ESCS) le sume las variables independientes computer_num e internet_num y las variables de control gender y country.

modelo_2 <- lm(formula = total_score ~ escs+ computer_num + internet_num + gender+ country,
               data = paises_inf,
               weights = stu_wgt)   
summary(modelo_2)
## 
## Call:
## lm(formula = total_score ~ escs + computer_num + internet_num + 
##     gender + country, data = paises_inf, weights = stu_wgt)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -21102.4   -317.3     51.8    388.4   9681.4 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  1048.7242     7.5639 138.649 < 0.0000000000000002 ***
## escs           15.9119     0.2850  55.832 < 0.0000000000000002 ***
## computer_num    6.5574     0.6371  10.292 < 0.0000000000000002 ***
## internet_num   24.7294     0.7137  34.651 < 0.0000000000000002 ***
## gendermale    -30.0894     0.5346 -56.284 < 0.0000000000000002 ***
## countryARE     -2.9489     9.4985  -0.310             0.756210    
## countryARG     41.7151     7.6900   5.425 0.000000058122094248 ***
## countryAUS     81.0785     8.4242   9.625 < 0.0000000000000002 ***
## countryAUT     86.5022    10.1100   8.556 < 0.0000000000000002 ***
## countryBEL     78.3905     9.1588   8.559 < 0.0000000000000002 ***
## countryBGR     26.5230     9.6028   2.762             0.005745 ** 
## countryBRA     38.5154     7.5608   5.094 0.000000350642461890 ***
## countryBRN     69.4015    20.2630   3.425             0.000615 ***
## countryCAN     99.1830     8.3478  11.881 < 0.0000000000000002 ***
## countryCHE    100.0827    10.3840   9.638 < 0.0000000000000002 ***
## countryCHL     75.5864     8.1612   9.262 < 0.0000000000000002 ***
## countryCOL     56.3520     7.6970   7.321 0.000000000000246296 ***
## countryCZE    107.6470     9.5690  11.250 < 0.0000000000000002 ***
## countryDEU     79.4727     7.8827  10.082 < 0.0000000000000002 ***
## countryDNK    108.3866    11.2028   9.675 < 0.0000000000000002 ***
## countryDOM    -14.2834     8.1886  -1.744             0.081107 .  
## countryESP     97.6131     7.9848  12.225 < 0.0000000000000002 ***
## countryEST    139.2545    21.9316   6.349 0.000000000216357414 ***
## countryFIN     91.3602    11.1381   8.202 0.000000000000000236 ***
## countryFRA     77.5567     7.7821   9.966 < 0.0000000000000002 ***
## countryGBR     84.2802     7.9586  10.590 < 0.0000000000000002 ***
## countryGEO     21.4744     9.7127   2.211             0.027040 *  
## countryGRC     68.7547     8.8139   7.801 0.000000000000006175 ***
## countryGTM     33.4540     8.0174   4.173 0.000030116474300198 ***
## countryHKG    114.3336    13.1615   8.687 < 0.0000000000000002 ***
## countryHRV    108.7410    11.9158   9.126 < 0.0000000000000002 ***
## countryHUN     79.1681     9.6615   8.194 0.000000000000000253 ***
## countryIDN     42.8672     7.5396   5.686 0.000000013047124963 ***
## countryIRL    125.9981    11.5193  10.938 < 0.0000000000000002 ***
## countryISL     52.5663    24.3039   2.163             0.030551 *  
## countryISR     31.4521     8.9879   3.499             0.000466 ***
## countryITA    102.0035     7.9330  12.858 < 0.0000000000000002 ***
## countryJAM     28.0339    11.1860   2.506             0.012206 *  
## countryJOR     -1.7383     8.0778  -0.215             0.829618    
## countryJPN    125.6815     7.9506  15.808 < 0.0000000000000002 ***
## countryKAZ     76.1755     7.8834   9.663 < 0.0000000000000002 ***
## countryKHM      3.0384     8.1557   0.373             0.709482    
## countryKOR     89.9031     8.3246  10.800 < 0.0000000000000002 ***
## countryKSV    -21.9260    10.8498  -2.021             0.043294 *  
## countryLTU    104.4262    13.7466   7.596 0.000000000000030516 ***
## countryLVA    120.9583    16.6396   7.269 0.000000000000362280 ***
## countryMAC    137.9697    42.4335   3.251             0.001148 ** 
## countryMAR     11.3994     7.6954   1.481             0.138521    
## countryMDA     57.5260    10.7841   5.334 0.000000095954786597 ***
## countryMEX     72.2804     7.5860   9.528 < 0.0000000000000002 ***
## countryMKD      5.2291    12.0996   0.432             0.665619    
## countryMLT     41.5268    28.1624   1.475             0.140335    
## countryMNE     33.9651    17.9443   1.893             0.058384 .  
## countryMNG     69.0575     9.7441   7.087 0.000000000001372397 ***
## countryMYS     64.4709     7.7833   8.283 < 0.0000000000000002 ***
## countryNLD     49.3906     8.7215   5.663 0.000000014881941641 ***
## countryNOR     60.5020    10.5025   5.761 0.000000008383629712 ***
## countryNZL     84.5631    11.3355   7.460 0.000000000000086765 ***
## countryPAN     36.1170     9.7436   3.707             0.000210 ***
## countryPER     61.8651     7.7279   8.005 0.000000000000001195 ***
## countryPHL    -28.7100     7.5592  -3.798             0.000146 ***
## countryPOL    100.6285     8.2246  12.235 < 0.0000000000000002 ***
## countryPRT    104.0156     9.4802  10.972 < 0.0000000000000002 ***
## countryPRY     10.1953     8.6149   1.183             0.236636    
## countryPSE      2.7429     8.4355   0.325             0.745055    
## countryQAT     36.4019    12.8035   2.843             0.004468 ** 
## countryQAZ     18.7763    10.2765   1.827             0.067685 .  
## countryQUR     70.6603     8.3645   8.448 < 0.0000000000000002 ***
## countryROU     38.5131     8.2943   4.643 0.000003429737969328 ***
## countrySAU     38.7713     7.8000   4.971 0.000000667565313598 ***
## countrySGP    118.6565    16.1195   7.361 0.000000000000182944 ***
## countrySLV     21.2061     8.7791   2.416             0.015713 *  
## countrySRB     74.4684     9.6582   7.710 0.000000000000012580 ***
## countrySVK     65.8878    10.6475   6.188 0.000000000609677141 ***
## countrySVN    101.4805    15.4687   6.560 0.000000000053765143 ***
## countrySWE     72.1870     9.4408   7.646 0.000000000000020753 ***
## countryTAP    105.1281     9.2755  11.334 < 0.0000000000000002 ***
## countryTHA     57.6347     7.6692   7.515 0.000000000000057033 ***
## countryTUR    119.6011     7.6909  15.551 < 0.0000000000000002 ***
## countryURY     66.1627    10.2651   6.445 0.000000000115457127 ***
## countryUSA     79.1439     7.5736  10.450 < 0.0000000000000002 ***
## countryUZB     -3.2188     7.6859  -0.419             0.675366    
## countryVNM    145.4914     7.7362  18.807 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1037 on 290231 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.1219, Adjusted R-squared:  0.1217 
## F-statistic: 491.4 on 82 and 290231 DF,  p-value: < 0.00000000000000022

Análisis de regresión múltiple

A primera vista todas las variables introducidas son significativas pero es al graficar (más adelante está el gráfico) que nos deparamos con el primer problema: la variable de control country explica toda la variabilidad del promedio, lógicamente porque hay promedios diferentes en cada país, entonces pareciera ‘tapar’ a las demás variables explicativas. ¿Aún así es lógico un R2 de 0.12? ¿Cuáles podrían ser entonces las dos alternativas? en primer lugar voy a correr un modelo múltiple tal como se viene trabajando, con interacciones entre variables que se solapan (computer_num e internet_num) y realizaré sobre este los test correspondientes para ver si cumple con los supuestos mínimos de un modelo de regresión. Por otro lado, correré un modelo multinivel, con el paquete lme4 para calcular los coeficientes de las variables explicativas dentro de cada país y que la variable country no ‘tape’ todo el análisis.

residuos_2 <- data.frame(Residuos= modelo_2$residuals, Ajustados= modelo_2$fitted.values)
ggplot(data = residuos_2, aes(x= Ajustados, y= Residuos))+
  geom_point(color= "gray")+
   geom_smooth(se = FALSE, color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed")+
   labs(
    title = "Residuos del modelo",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 9: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

Si bien el gráfico de residuos (Figura 9) muestra “menos” problemas que el anterior, puesto que la curvatura no es tan pronunciada, aún se observa un problema de linealidad. Pero, además, el summary del modelo resulta muy poco intuitivo para las interpretaciones y el listado de coeficientes de cada país acaba por confundir más de lo que aclara.

Análisis de regresión múltiple con interacciones

Algo que también podría estar sucediendo es que tengamos variables solapadas, es decir, variables que se explican por o con otras y que generan no necesariamente un problema de correlación pero si de multicolinearidad. Estas variables podrían ser computer_num e internet_num. Podría pensarse que quien tenga computadora tendrá internet, puesto que la computadora es, entre otras cosas, un puerto de acceso a la red. Aún así, algunos de los estudiantes que tengan internet podrían tener computadora, pero también podrían conectarse a través de otros dispositivos (smartphone, por ejemplo). Por ello, vamos a correr nuevamente el modelo de regresión múltiple pero contemplaando las interacciones posibles. Además vamos a sacar país que genera el conflicto analizado previamente y solo se dejará como variable de control a gender y los pesos muestrales.

modelo_int <- lm(total_score ~ computer_num * internet_num + escs + gender, 
                 data = paises_inf,
                 weights = stu_wgt)
summary(modelo_int)
## 
## Call:
## lm(formula = total_score ~ computer_num * internet_num + escs + 
##     gender, data = paises_inf, weights = stu_wgt)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -20311.9   -335.9     50.2    404.8  10358.9 
## 
## Coefficients:
##                            Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)               1094.5726     0.9339 1172.08 <0.0000000000000002 ***
## computer_num               -12.0428     1.1463  -10.51 <0.0000000000000002 ***
## internet_num                20.0279     0.8829   22.68 <0.0000000000000002 ***
## escs                        17.5124     0.2803   62.48 <0.0000000000000002 ***
## gendermale                 -29.4749     0.5523  -53.37 <0.0000000000000002 ***
## computer_num:internet_num   37.0021     1.3240   27.95 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1072 on 290308 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.06122,    Adjusted R-squared:  0.06121 
## F-statistic:  3786 on 5 and 290308 DF,  p-value: < 0.00000000000000022

En el modelo corrido previamente le pedimos a R que analice la relación entre total_score (y) y las variables explicativas escs, computer_num e internet_num, teniendo como variable de control gender y ponderando los pesos muestrales. Pero, además, le pedimos que genere una interacción entre computer_num e internet_num para analizar la relación entre el puntaje y la computadora, teniendo en cuenta el efecto de tener o no internet, y viceversa. Esto último para ver qué ocurre con estas variables que podrían estar explicándose mutuamente. Vamos a ver qué cambió con respecto al modelo sin interacciones.

Con respecto a la media de los residuos, que podrían indicar cómo se distribuyen los datos, en el modelo previo era de 51.8 mientras que en el modelo con interacción y sin country es de 50. Mejora levemente pero sigue habiendo valores extremos que podrían indicar la presencia de outliers. Por otro lado, el R2 que en el modelo anterior era de 0.12 y que indicaba que el modelo permitía explicar el 12% de la variabilidad del puntaje desciende a 0.06 en el modelo nuevo, por lo que este solo explicaría el 6% de la variabilidad del puntaje. Tiene sentido sobre todo recordando que quitamos la variable country que era por demás explicativa.

Alpha, que es el valor de puntaje promedio de los estudiantes cuando NO tienen computadora ni internet y su índice de ESCS es bajo, arroja un valor de 1094, por encima del valor de Alpha en el modelo 2 de este análisis que era de 1048. El índice de ESCS se mantiene con valores ‘relativamente’ similares: en el modelo con interacción es 17.5 (Beta3) y en el modelo 2 era de 15.9 (Beta1). Esto significa que al aumentar una unidad de ESCS el puntaje podría aumentar, en promedio, 17.9 o 15.9 dependiendo el modelo. Los valores de Beta para computer_num y para internet_num deben ser analizados a partir de la interacción que arroja un valor de 37 lo cual significa que: el efecto de ambas variables combinadas es relevante ya que el hecho de tener computadora y de tener internet podría aumentar el puntaje, en promedio, 37 puntos. Veamos el gráfico de residuos:

residuos_int <- data.frame(Residuos= modelo_int$residuals, Ajustados= modelo_int$fitted.values)
ggplot(data = residuos_int, aes(x= Ajustados, y= Residuos))+
  geom_point(color= "gray")+
   geom_smooth(se = FALSE, color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed")+
   labs(
    title = "Residuos del modelo",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 10: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

El gráfico 10 sigue demostrando lo que venimos suponiendo: no hay un problema grave del modelo pero parecería haber un problema de linealidad puesto que la línea de los valores ajustados se curva hacia abajo. Antes de los test para evaluar si esto efectivamente es así veamos que pasa si corremos un análisis multinivel.

Análisis de regresión múltiple multinivel

modelo_ml <- lmer(total_score ~ escs + computer_num + internet_num + gender + (1 | country),
                  data = paises_inf)
summary(modelo_ml)
## Linear mixed model fit by REML ['lmerMod']
## Formula: total_score ~ escs + computer_num + internet_num + gender + (1 |  
##     country)
##    Data: paises_inf
## 
## REML criterion at convergence: 3700317
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2962 -0.6686  0.1089  0.7691  2.6906 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  country  (Intercept)  1629     40.36  
##  Residual             20062    141.64  
## Number of obs: 290314, groups:  country, 79
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)  1106.0158     4.6533  237.69
## escs           20.0294     0.2887   69.37
## computer_num   16.7169     0.6913   24.18
## internet_num   28.4090     0.8237   34.49
## gendermale    -32.2908     0.5278  -61.18
## 
## Correlation of Fixed Effects:
##             (Intr) escs   cmptr_ intrn_
## escs         0.111                     
## computer_nm -0.104 -0.299              
## internet_nm -0.147 -0.217 -0.151       
## gendermale  -0.066 -0.057  0.035  0.007

A primera vista puede advertirse que el summary() de este modelo difiere del realizado previamente con los otros modelos, esto es debido a la forma en que lme4 analiza las variables e incluye el “efecto país”. Recordemos que la idea era ver cuánto de los puntajes totales podrían estar explicándose por este efecto país lo que podría llevar a que los R2 de los modelos que no lo incluyen como variable de control diera tan bajo. por ello, se corrió este modelo multiniveles para ver cómo se relacionan las variables entre sí incluyendo, o teniendo en cuenta, este efecto.

En primer lugar, hay dos clases de elementos que pueden observarse: random effects y fixed effects. En relación a los random effects como primera medida pueden analizarse los desvíos estandar de los países que nos permiten afirmar que la variabilidad de los puntajes total podría cambiar entre países, en promedio, ± 40 puntos. Sin embargo, esa variabilidad aumenta dentro de los mismos países, ascendiendo a ± 142 puntos. Esto tiene sentido si se piensa que todos los países de este grupo tienen total_score de sus estudiantes por debajo de la media de OCDE pero, si recordamos, aún así las diferencias entre los puntajes mínimos y máximos era amplia dentro de cada país. Si se hiciera la comparación con la muestra de países que están por encima de la media podríamos, quizás, ver otros resultados.

En relación a los fixed effects podemos analizar, tal como se hizo en los demás modelos, los coeficientes. En primer lugar, Alpha que será el puntaje para los estudiantes que no tengan computadora, ni internet y que tengan un bajo índice de ESCS, es de 1106, levemente por debajo de Alpha en los otros modelos. Beta1 que marca aumento promedio en el puntaje de los estudiantes cuando aumenta una unidad de ESCS es de 20 puntos, mientras que el hecho de tener computadora marca una diferencia positiva en el puntaje de, en promedio, 18 puntos (Beta2) y, tener internet, marca una diferencia positiva de puntaje de, en promedio, 28.4 puntos (Beta3) siendo la covariable más significativa si se tiene en cuenta el efecto país.

Por otro lado, puede advertirse que en este modelo la media de los residuos es menor, un valor de 0.11 y no se observan valores extremos.

df_modelo_ml <- paises_inf[rownames(modelo_ml@frame), ]

residuos_ml <- df_modelo_ml %>%
  mutate(
    resid = resid(modelo_ml, type = "pearson"),
    fitted = fitted(modelo_ml))
ggplot(data = residuos_ml, aes(x= fitted, y= resid))+
  geom_point(color= "gray")+
   geom_smooth(se = FALSE, color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed")+
   labs(
    title = "Residuos del modelo",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 11: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

Conclusión luego de correr tres tipos de modelos de regresión lineal diferentes(múltiple, múltiple con interacciones y multinivel) hay que revisar si se cumplen los supuestos porque podría estar habiendo problemas en, al menos, la linealidad y la homocedasticidad.

Supuestos del modelo

Para chequear las sospechas que tenemos en cuanto a los modelos corridos vamos a realizar, en primer lugar, los test correspondientes para chequear el Modelo con interacciones. Luego, le aplicaremos las transformaciones que fueran necesarias. En segundo lugar, realizaremos también los test para el Modelo Multinivel y se procederá también a efectuar transformaciones requeridas.

Test al modelo con interacciones

test_norm <- ad.test(modelo_int$residuals)
test_norm
## 
##  Anderson-Darling normality test
## 
## data:  modelo_int$residuals
## A = 2010.8, p-value < 0.00000000000000022

Tal como puede observarse en el resultado del Test, donde el P-Value arroja un valor por debajo del 0.05 HABRIA problemas de distribución de los datos. Pero revisemos los gráficos de los residuos para poder constatarlo de forma gráfica:

qqnorm(modelo_int$residuals, col= "azure4")          
qqline(modelo_int$residuals, col="red")  

El gráfico previo demuestra que hay una separación de la línea recta hacia el extremo, lo cual podría indicar la presencia de outliers o valores muy extremos (algo de esto ya habíamos chequeado en el summary) o bien que no hay una linealidad sino otra relación entre los datos. Chequeemos ahora un historgrama de los residuos:

media_res_int <- mean(modelo_int$residuals)
histo_mod_int <- ggplot(data = modelo_int, aes(x=modelo_int$residuals))+
  geom_histogram(color= "azure4")+
  geom_vline(xintercept = media_res_int, color="red", linetype = "dashed", size = 1)+
   labs(
    title = "Histograma de residuos",
    subtitle = "Modelo con interacciones",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 13: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))
histo_mod_int

El histograma de la Figura 13 que solo muestra cómo están distribuídos los residuos en general arroja un resultado bastante “similar” a una campana, lo cual podría indicar que el desvío del QQPlot (Figura 12) podría no ser de gravedad en nuestro modelo y que, probablemente, pueda corregirse a través de una transformación logarítmica.

library(lmtest)
reset_result <- resettest(modelo_int)
reset_result
## 
##  RESET test
## 
## data:  modelo_int
## RESET = 83.23, df1 = 2, df2 = 290306, p-value < 0.00000000000000022

El test de normalidad, cuyo valor de P-value es menor a 0.05 indicaría que, además de tener problemas de distribución normal de los valores también podría tratarse de problemas de linealidad, que podrían advertir que no siempre que se mueva X se moverá Y. Todo parece indicar que si realizamos una prueba para ver si nuestros datos tienen heterocedasticidad seguramente arroje un resultado positivo, veamos:

homoces <- bptest(modelo_int)
homoces
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_int
## BP = 15113085, df = 5, p-value < 0.00000000000000022

Conclusión: chequeados normalidad, linealidad y homocedasticidad puede concluirse que el modelo de interacciones presenta una distribución NO normal de residuos, problemas en cuanto a la linealidad y una fuerte heterocedasticidad. ¿Qué hacer? en primer lugar vamos a ver si podemos robustecer los errores estánbdar para que no de estimaciones ineficientes y haga pasar por significativo algo que en verdad no lo es.

summary(modelo_int)
## 
## Call:
## lm(formula = total_score ~ computer_num * internet_num + escs + 
##     gender, data = paises_inf, weights = stu_wgt)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -20311.9   -335.9     50.2    404.8  10358.9 
## 
## Coefficients:
##                            Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)               1094.5726     0.9339 1172.08 <0.0000000000000002 ***
## computer_num               -12.0428     1.1463  -10.51 <0.0000000000000002 ***
## internet_num                20.0279     0.8829   22.68 <0.0000000000000002 ***
## escs                        17.5124     0.2803   62.48 <0.0000000000000002 ***
## gendermale                 -29.4749     0.5523  -53.37 <0.0000000000000002 ***
## computer_num:internet_num   37.0021     1.3240   27.95 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1072 on 290308 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.06122,    Adjusted R-squared:  0.06121 
## F-statistic:  3786 on 5 and 290308 DF,  p-value: < 0.00000000000000022
coeftest(modelo_int, vcov. = vcovHC(modelo_int, type = "HC1"))
## 
## t test of coefficients:
## 
##                             Estimate Std. Error  t value              Pr(>|t|)
## (Intercept)               1094.57257    2.30758 474.3385 < 0.00000000000000022
## computer_num               -12.04281    3.10171  -3.8826             0.0001034
## internet_num                20.02792    2.13997   9.3590 < 0.00000000000000022
## escs                        17.51243    0.66134  26.4802 < 0.00000000000000022
## gendermale                 -29.47485    1.31352 -22.4395 < 0.00000000000000022
## computer_num:internet_num   37.00210    3.42881  10.7915 < 0.00000000000000022
##                              
## (Intercept)               ***
## computer_num              ***
## internet_num              ***
## escs                      ***
## gendermale                ***
## computer_num:internet_num ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ese cambio “robustece” los errores de estándar y eso hace que nos arroje valores más eficientes o reales de p-value, pero no modifica los coeficientes. Además de eso, y considerando los problemas de normalidad, linealidad y la posible presencia de outliers, voy a transformar también el modelo mediante una Transformación logarítmica aplicada sobre total_score, la variable dependiente. Elegí aplicar la transformación sobre esta variable porque 1) el puntaje total es el que arroja los valores extremos que pueden estar distorsionando el modelo 2) las variables computer_num e internet_num son factoriales, por lo que considero que no tendría sentido hacerlo. ¿Cómo sé que mi variable dependiente, total_score, es la que arroja valores extremos? mediante la observación de sus valores en el summary de la base de datos; pero, para ayudar a la comprensión, también graficaré con un histograma estos valores:

ggplot(paises_inf, aes(x = total_score)) +
      geom_histogram(aes(y = after_stat(density)), fill = "skyblue3", color= "skyblue3") +
      geom_density(color = "dodgerblue4", linewidth = 1.5, fill = "dodgerblue4", alpha = 0.3)+
  labs(caption = "Figura 14: Elaboración propia con datos de OCDE (2023)")+
  theme_minimal()+
   theme (plot.caption = element_text(hjust = 0.5))

Como puede observarse en el histograma de total_score (Figura 14) hay una cola larga hacia la izquierda, y casi nada hacia la derecha, por lo que no hay una distribución normal de los valores.

paises_inf$log_total_score <- log(paises_inf$total_score)
modelo_logy <- lm(log_total_score ~ computer_num * internet_num + escs + gender,
                  data = paises_inf, weights = stu_wgt)
residuos_log <- data.frame(Residuos= modelo_logy$residuals, Ajustados= modelo_logy$fitted.values)
ggplot(data = residuos_log, aes(x= Ajustados, y= Residuos))+
  geom_point(color= "gray")+
   geom_smooth(se = FALSE, color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed")+
   labs(
    title = "Residuos del modelo",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 15: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

media_res_log <- mean(modelo_logy$residuals)
histo_mod_log <- ggplot(data = modelo_logy, aes(x=modelo_logy$residuals))+
  geom_histogram(color= "azure4")+
  geom_vline(xintercept = media_res_log, color="red", linetype = "dashed", size = 1)+
   labs(
    title = "Histograma de residuos",
    subtitle = "Modelo con transformación algorítmica",
    x = "Valores ajustados",
    y = "Residuos",
    caption = "Figura 16: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal()+
   theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

Se realiza nuevamente en Test de Anderson Darling para chequear normalidad pero, comp puede verse a continuación, y comop podía visualizarse en la figura 16, el modelo sigue teniendo problemas de distribución de los datos:

ad.test(modelo_logy$residuals)
## 
##  Anderson-Darling normality test
## 
## data:  modelo_logy$residuals
## A = 3665.3, p-value < 0.00000000000000022
resettest(modelo_logy)
## 
##  RESET test
## 
## data:  modelo_logy
## RESET = 72.286, df1 = 2, df2 = 290306, p-value < 0.00000000000000022
histo_mod_int+histo_mod_log

Test de hipótesis

A medida que el tiempo va transcurriendo los Estados intentar diseñar e implementar políticas educativas vinculadas a la mejora en la calidad, lo cual debería poder verse reflejado en resultados de pruebas internacionales como, por ejemplo, las pruebas PISA. Partiendo de esta premisa se plantea que los estudiantes que realizaron la prueba en 2022, pertenecientes a los países de mejores desempeños, alcanzaron resultados más altos que los estudiantes de los mismos países en la prueba de 2018. Para comprobar esto se hará una Test de Hipótesis. En primer lugar se agrupan los países que realizaron la prueba en 2022 y se calcula el promedio de los puntajes obtenidos por los estudiantes para ordenarlos y visualizar cuáles son aquellos con mejor desempeño. Luego, se comparan las poblaciones de 2018 y 2022 para ver cuáles de estos países con alto desempeño realizaron la prueba en las dos ediciones. Posteriormente se agrupan estos datos en un solo dataframe que permita realizar el test. Por último, se realiza el test para comprobar o rechazar la hipótesis.

promedios_paises <- students_2022 %>%
  group_by(country) %>%
  summarise(promedio = mean(total_score, na.rm = TRUE)) %>%
  arrange(promedio)

Una vez agrupados y ordenados según resultados del 2022, se analizan cuáles de estos países tuvieron también participación en la pueba PISA inmediatamente anterior: 2018.

paises_2018 <- unique(student_subset_2018$country)
paises_2022 <- unique(student_subset_2022$country)
paises_comunes <- intersect(paises_2018, paises_2022)

paises_comunes
##  [1] AUS AUT BEL CAN CHE CHL COL CRI CZE DEU DNK ESP EST FIN FRA GBR GRC HUN IRL
## [20] ISL ISR ITA JPN KOR LTU LVA MEX NLD NOR NZL POL PRT SVK SVN SWE TUR USA
## 89 Levels: ALB ARE ARG AUS AUT BEL BGR BIH BLR BRA BRN CAN CHE CHL COL ... UZB

Entre los países con mejores resultados en las pruebas PISA de 2022, que realizaron también la prueba en 2018, se encuentran los siguientes países: JPN: Japón; KOR: Corea; EST: Estonia; IRL: Irlanda y, CZE: República Checa. Con esta información se armarán dos nuevas muestras: paises_2018 y países_2022, cada una con los resultados filtrados tan solo a estos países.

paises_2018 <- student_subset_2018 %>% 
  filter(country %in% c("JPN", "KOR", "EST", "IRL", "CZE"))
paises_2022 <- student_subset_2022 %>% 
  filter(country %in% c("JPN", "KOR", "EST", "IRL", "CZE"))
paises_2018 <- paises_2018 %>% 
  mutate(year = 2018)  #Agregar columna de año n cada muestra para que no se me mezclen al unirlos
paises_2022 <- paises_2022 %>% 
  mutate(year = 2022)
paises_2018 <- paises_2018 %>% 
  mutate(total_score = math+read+science)
paises_2022 <- paises_2022 %>% 
  mutate(total_score = math+read+science)
muestra_testhip <- bind_rows(paises_2018, paises_2022)  #Unirlos
test_results <- t.test(total_score ~ year, data = muestra_testhip)
test_results
## 
##  Welch Two Sample t-test
## 
## data:  total_score by year
## t = -2.532, df = 497.91, p-value = 0.01165
## alternative hypothesis: true difference in means between group 2018 and group 2022 is not equal to 0
## 95 percent confidence interval:
##  -104.38225  -13.16782
## sample estimates:
## mean in group 2018 mean in group 2022 
##           1511.255           1570.030

El resultado del Test, con un valor de p-value de 0.011 demuestra que puede rechazarse la hipótesis nula o H0: hay diferencia significativa entre las medias de total_score obtenidas por los estudiantes en 2018 y en 2022 y los resultados demuestran un mejor desempeño en este último período. Además del p-value, el test también indica que el puntaje obtenido, en promedio, en 2018 fue de 1511 mientras que en 2022 ascendió hasta los 1570, con 59 puntos de diferencia. Además, puede afirmarse que los estimadores de las medias se encontrarán entre los -104 y -13 puntos. Cabe destacar que estos valores son negativos puesto que los estudiantes de 2018 obtuvieron, en promedio, menor puntaje que los del 2022. Visualicemos gráficamente estas diferencias:

set.seed(123) 
ggplot(muestra_testhip, aes(x = factor(year), y = total_score, fill = factor(year))) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 1, color= "black") +
  scale_fill_manual(values = c("2018"= "#104E8B", "2022"= "#68228B")) +
  stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "white") +
  labs(
    x = "Año",
    y = "Puntaje Total",
    title = "Distribución de puntajes por año (2018 vs 2022)",
    caption = "Figura 17: Elaboración propia con datos de OCDE (2023)") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5))

En la figura 17 puede observarse como la media de los puntajes totales de los estudiantes para la prueba 2022 está por encima de la media de 2018, al igual que el extremo de abajo de la figura del violín, reafirmando lo que se concluyó al interpretar el test de hipótesis previamente: los estudiantes de países con alto desempeño en las pruebas PISA obtuvieron, en promedio, mejores resultados que en la prueba de 2018.

Conclusión

Luego de correr diversos modelos de regresión y de hacer cálculos de coeficientes de correlación y test de medias, en la primera parte de este reporte, puede concluirse que si bien los modelos presentan algunas dificultades en cuanto a los supuesto de normalidad y linealidad, es 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. ]A mayor índice ESCS, teniendo computadora y teniendo conexión a internet el puntaje total mejora. Además, se concluye que, para variables complejas como el puntaje de pruebas educativas y, teniendo en cuenta la cantidad de datos de la muestra, sería necesario incluir más variables explicativas que den cuenta de la multiplicidad de factores que indicen en el desempeño de los estudiantes. En la segunda parte de este reporte se buscó indagar si el resultado de 2018 y el del 2022 eran diferentes en los estudiantes pertenecientes a países con alto nivel de desempeño en PISA, con la hipótesis de que los resultados de 2022 serían mayores. A través del test de hipótesis se concluyó que SI hay difernecia entre ambos resultados y que, tal como indicaba la hipótesis, los resultados de 2022 fueron mejores que los del período anterior.