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
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.
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).
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.
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).
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.
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()
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.
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
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)}\]
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")
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")
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.
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
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.
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.
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.
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
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.
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.