Una comparación bayesiana de los resultados de los estudiantes que abandonan la escuela con R y brms
Comparación bayesiana de resultados de estudiantes que abandonan escuela con R y brms
ANOVA — Estilo Bayesiano
Se le da mucha importancia a lo que queremos hacer cuando dejamos la escuela. Se nos pregunta desde pequeños qué queremos hacer cuando crezcamos, y luego pasamos 13 años en educación pre-terciaria. En política pública, se hace mucho hincapié en las diferencias entre los sistemas escolares del Gobierno, Católico e Independiente, especialmente en lo que respecta a la financiación pública de las escuelas no gubernamentales y cómo se deben asignar los recursos.
¿Hay alguna diferencia real entre los resultados de los estudiantes según el sector escolar?
En Victoria, Australia, el Gobierno Estatal realiza una encuesta anual para evaluar los resultados posteriores a la escuela secundaria (Encuesta On Track). Este conjunto de datos (disponible bajo la licencia CC BY 4.0) está disponible para varios años, y solo analizaremos el más reciente disponible al momento de la redacción de este artículo, 2021.
Este artículo buscará responder a estas preguntas aplicando metodologías de análisis bayesiano, utilizando el paquete R brms.
Cargar Librerías y Conjunto de Datos
A continuación, cargamos los paquetes necesarios. El conjunto de datos se presenta en un formato de informe amplio, con muchas celdas fusionadas. R no lo soporta, por lo que debemos hacer un poco de codificación para cambiar los nombres de los vectores y crear un marco de datos ordenado.
- Optimización de viaje por Europa Algoritmos Genéticos y la API de G...
- Procesamiento del Lenguaje Natural para Principiantes Absolutos
- Arquitecturas Multi-Tarea Una Guía Completa
library(tidyverse) #Paquete meta de Tidyverselibrary(brms) #Paquete de Modelado Bayesiano library(tidybayes) #Funciones Auxiliares Ordenadas y Geoms de Visualización para Modelos Bayesianoslibrary(readxl)df <- read_excel("DestinationData2022.xlsx", sheet = "SCHOOL PUBLICATION TABLE 2022", skip = 3)colnames(df) <- c('vcaa_code', 'school_name', 'sector', 'locality', 'total_completed_year12', 'on_track_consenters', 'respondants', 'bachelors', 'deferred', 'tafe', 'apprentice_trainee', 'employed', 'looking_for_work', 'other') df <- drop_na(df)Aquí se muestra una muestra de cómo se ve nuestro conjunto de datos.
El conjunto de datos tiene 14 campos.
- Código VCAA: ID administrativo para cada código
- Nombre de la Escuela
- Sector: G = Gobierno, C = Católico e I = Independiente
- Localidad o Suburbio
- Total de Estudiantes que Completan el Año 12
- Consentimiento de Seguimiento
- Encuestados
- El Resto de las Columnas Representan el Porcentaje de Encuestados para Cada Resultado
Para los propósitos de nuestro estudio transversal, nos interesa las proporciones de resultados por sector, por lo que debemos convertir este marco de datos amplio en un formato largo.
df_long <- df |> mutate(across(5:14, ~ as.numeric(.x)), #convertir todos los campos numéricos capturados como caracteres across(8:14, ~ .x/100 * respondants), #calcular recuentos a partir de proporciones across(8:14, ~ round(.x, 0)), #redondear a números enteros respondants = bachelors + deferred + tafe + apprentice_trainee + employed + looking_for_work + other) |> #recalcular el total de encuestados |> filter(sector != 'A') |> #Volumen bajo select(sector, school_name, 7:14) |> pivot_longer(c(-1, -2, -3), names_to = 'outcome', values_to = 'proportion') |> mutate(proportion = proportion/respondants)
Análisis Exploratorio de Datos
Vamos a visualizar y resumir brevemente nuestro conjunto de datos. El sector gubernamental tiene 157 escuelas, el sector independiente tiene 57 escuelas y el sector católico tiene 50 escuelas en este conjunto de datos.
df |> mutate(sector = fct_infreq(sector)) |> ggplot(aes(sector)) + geom_bar(aes(fill = sector)) + geom_label(aes(label = after_stat(count)), stat = 'count', vjust = 1.5) + labs(x = 'Sector', y = 'Cantidad', title = 'Cantidad de Escuelas por Sector', fill = 'Sector') + scale_fill_viridis_d(begin = 0.2, end = 0.8) + theme_ggdist()
df_long |> ggplot(aes(sector, proportion, fill = outcome)) + geom_boxplot(alpha = 0.8) + geom_point(position = position_dodge(width=0.75), alpha = 0.1, color = 'grey', aes(group = outcome)) + labs(x = 'Sector', y = 'Proporción', fill = 'Resultado', title = 'Diagrama de Caja de las Proporciones de los Encuestados por Sector y Resultado') + scale_fill_viridis_d(begin = 0.2, end = 0.8) + theme_ggdist()
Las distribuciones cuentan una historia interesante. El resultado de Bachelors muestra la mayor variabilidad en todos los sectores, con las escuelas independientes teniendo la mediana más alta de proporción de estudiantes que eligen seguir estudios universitarios. Es interesante notar que las escuelas gubernamentales tienen la mediana más alta de proporción de estudiantes que optan por el empleo después de terminar la secundaria. Todos los demás resultados no varían mucho, lo veremos más adelante.
Modelado Estadístico – Regresión de Probabilidad Beta
Nos centramos en las proporciones de estudiantes por escuela y sus resultados después de graduarse de la secundaria. Una probabilidad beta es nuestra mejor opción en estos casos. brms es un paquete brillante desarrollado por Buerkner para desarrollar modelos bayesianos. Nuestro objetivo es modelar la distribución de las proporciones por resultado y sector.
La regresión beta modela dos parámetros, μ y φ. μ representa la proporción media, y φ es una medida de precisión (o varianza).
Nuestro primer modelo se representa a continuación, notemos que ya estamos comenzando con una interacción entre el Sector y el Resultado porque esta es la pregunta que queremos que el modelo responda, y por lo tanto este es un modelo de ANOVA.
Estamos pidiendo al modelo que cree términos Beta individuales para cada combinación de Sector y Resultado, con un término φ agrupado, o que modele diferentes medias de proporción con la misma varianza. Esperamos que el 50% de estas proporciones se encuentren entre (-3, 1) en la escala logit o (0.01, 0.73) como proporciones. Esto es lo suficientemente amplio pero con una estimación previa informada. La estimación global de Phi es un número positivo, por lo que usamos una distribución gamma que es lo suficientemente amplia.
# Modelando la Proporción con el Término de Interacción Sector:Resultado utilizando Beta - Término de Phi Agrupadom3a <- brm( proportion ~ sector:outcome + 0, family = Beta, data = df_long |> mutate(proportion = proportion + 0.01), # La regresión beta no puede tomar resultados que sean 0, así que aquí agregamos un incremento pequeño prior = c(prior(normal(-1, 2), class = 'b'), prior(gamma(6, 1), class = 'phi')), iter = 8000, warmup = 1000, cores = 4, chains = 4, seed = 246, save_pars = save_pars(all = T), control = list(adapt_delta = 0.99, max_treedepth = 15)) |> add_criterion(c('waic', 'loo'), moment_match = T)summary(m3a)
Tenga en cuenta que en la configuración del modelo hemos agregado un incremento del 1% a todas las proporciones — esto se debe a que la regresión beta no maneja valores cero. Intentamos modelar esto con una beta inflada de cero, pero tardó mucho más en converger.
De manera similar, podemos modelar sin un phi agrupado, esto intuitivamente tiene sentido dado lo que observamos en las distribuciones anteriores, hay variación entre cada una de las combinaciones de sector y resultado. El modelo se define a continuación.
m3b <- brm( bf(proportion ~ sector:outcome + 0, phi ~ sector:outcome + 0), family = Beta, data = df_long |> mutate(proportion = proportion + 0.01), prior = c(prior(normal(-1, 2), class = 'b')), iter = 8000, warmup = 1000, cores = 4, chains = 4, seed = 246, save_pars = save_pars(all = T), control = list(adapt_delta = 0.99)) |> add_criterion(c('waic', 'loo'), moment_match = T)summary(m3b)
Diagnóstico y Comparación del Modelo
Con dos modelos en mano, ahora compararemos su precisión predictiva fuera de muestra utilizando la estimación LOO bayesiana de la densidad predictiva logarítmica puntual esperada (elpd_loo).
comp <- loo_compare(m3a, m3b)print(comp, simplify = F)
En pocas palabras, cuanto mayor sea el valor esperado del logaritmo puntual al dejar uno fuera, mayor será la precisión predictiva en datos no vistos. Esto nos da una buena medida relativa de la precisión del modelo entre modelos. Podemos verificar esto aún más mediante una verificación predictiva posterior, una comparación de valores observados y simulados. En nuestro caso, el modelo m3b es el mejor modelo para los datos observados.
alt_df <- df_long |> select(sector, outcome, proportion) |> rename(value = proportion) |> mutate(y = 'y', draw = 99) |> select(sector, outcome, draw, value, y)sim_df <- expand_grid(sector = c('C', 'I', 'G'), outcome = unique(df_long$outcome)) |> add_predicted_draws(m3b, ndraws = 1200) |> rename(value = .prediction) |> ungroup() |> mutate(y = 'y_rep', draw = rep(seq(from = 1, to = 50, length.out = 50), times = 504)) |> select(sector, outcome, draw, value, y) |> bind_rows(alt_df)sim_df |> ggplot(aes(value, group = draw)) + geom_density(aes(color = y)) + facet_grid(outcome ~ sector, scales = 'free_y') + scale_color_manual(name = '', values = c('y' = "darkblue", 'y_rep' = "grey")) + theme_ggdist() + labs(y = 'Densidad', x = 'y', title = 'Distribución de Proporciones Observadas y Replicadas por Sector y Resultado')
El modelo m3b, dado su término de varianza no agrupada o phi, es capaz de capturar mejor la variación en las distribuciones de proporciones por sector y resultado.
ANOVA — Estilo Bayesiano
Recuerde, nuestra pregunta de investigación es si hay diferencias en los resultados entre los sectores y en qué medida. En estadística frecuentista, podríamos usar ANOVA, un enfoque de diferencia de medias entre grupos. La debilidad de este enfoque es que los resultados proporcionan una estimación y un intervalo de confianza, sin sentido de incertidumbre acerca de estas estimaciones y un valor p perverso que indica si la diferencia entre las medias es estadísticamente significativa o no. No, gracias.
A continuación, generamos un conjunto de valores esperados para cada combinación de interacción entre sector y resultado, y luego utilizamos la excelente función tidybayes::compare_levels() que realiza el trabajo pesado. Calcula la diferencia en las medias posteriores entre los sectores para cada resultado. Excluimos el resultado ‘otro’ por motivos de brevedad.
new_df <- expand_grid(sector = c('I', 'G', 'C'), outcome = c('apprentice_trainee', 'bachelors', 'deferred', 'employed', 'looking_for_work', 'tafe'))epred_draws(m3b, newdata = new_df) |> compare_levels(.epred, by = sector, comparison = rlang::exprs(C - G, I - G, C - I)) |> mutate(sector = fct_inorder(sector), sector = fct_shift(sector, -1), sector = fct_rev(sector)) |> ggplot(aes(x = .epred, y = sector, fill = sector)) + stat_halfeye() + geom_vline(xintercept = 0, lty = 2) + facet_wrap(~ outcome, scales = 'free_x') + theme_ggdist() + theme(legend.position = 'bottom') + scale_fill_viridis_d(begin = 0.4, end = 0.8) + labs(x = 'Diferencia Proporcional', y = 'Sector', title = 'Diferencias en las Medias Posteriores por Sector y Resultado', fill = 'Sector')
Alternativamente, podemos resumir todas estas distribuciones con una tabla ordenada para una interpretación más fácil y un intervalo de credibilidad del 89%.
marg_gt <- epred_draws(m3b, newdata = new_df) |> compare_levels(.epred, by = sector, comparison = rlang::exprs(C - G, I - G, C - I)) |> median_qi(.width = .89) |> mutate(across(where(is.numeric), ~round(.x, 3))) |> select(-c(.point, .interval, .width)) |> arrange(outcome, sector) |> rename(diff_in_means = .epred, Q5.5 = .lower, Q94.5 = .upper) |> group_by(outcome) |> gt() |> tab_header(title = 'Distribuciones Marginales de los Sectores por Resultado') |> #tab_stubhead(label = 'Comparación de Sectores') |> fmt_percent() |> gtExtras::gt_theme_pff()
Por ejemplo, si tuviéramos que resumir una comparación en un informe formal, podríamos escribir lo siguiente.
Los estudiantes de las escuelas gubernamentales tienen menos probabilidades de comenzar estudios universitarios que sus contrapartes en escuelas católicas e independientes después de completar el VCE.
En promedio, el 42.5% (entre el 39.5% y el 45.6%) de los estudiantes de escuelas gubernamentales, el 53.2% (entre el 47.7% y el 58.4%) de los estudiantes de escuelas católicas y el 65% (entre el 60.1% y el 69.7%) de los estudiantes de escuelas independientes comienzan su educación universitaria después de completar el año 12.
Hay una probabilidad posterior del 89% de que la diferencia entre la inscripción universitaria de estudiantes católicos y gubernamentales esté entre el 5.6% y el 15.7% con una media del 10.7%. Además, la diferencia entre la inscripción universitaria de estudiantes independientes y gubernamentales está entre el 17.8% y el 27% con una media del 22.5%.
Estas diferencias son sustanciales y hay un 100% de probabilidad de que estas diferencias no sean cero.
Resumen y Conclusión
En este artículo hemos demostrado cómo modelar datos proporcionales utilizando una función de verosimilitud Beta mediante el modelado bayesiano, y luego realizar un ANOVA bayesiano para explorar las diferencias en los resultados proporcionales entre sectores.
No hemos buscado crear una comprensión causal de estas diferencias. Se puede imaginar que hay varios factores que influyen en el rendimiento de los estudiantes individuales, antecedentes socioeconómicos, logros educativos de los padres, además de impactos a nivel escolar, recursos, etc. Es un área increíblemente compleja de la política pública que a menudo se enreda en retórica de suma cero.
Personalmente, soy la primera persona en mi familia extendida en asistir y completar educación terciaria. Asistí a una escuela secundaria pública de nivel medio y obtuve calificaciones razonablemente buenas para ingresar a mi primera preferencia. Mis padres me animaron a seguir estudiando, ambos optaron por dejar la escuela a los 16 años. Si bien este artículo proporciona evidencia de que la diferencia entre las escuelas gubernamentales y no gubernamentales es real, es puramente descriptivo en su naturaleza.
