Análisis Multinomial del Tipo de Programa Escolar según Desempeño Académico
MJ
2025-11-26
Introducción
El conjunto de datos hsbdemo describe el rendimiento acádemico de estudiantes de secundaria en diferentes áreas (lectura, escritura, matemáticas, ciencias y estudios sociales). Además, cada estudiante está clasificado en uno de los 3 tipos de programas: - General - Academic - Volcation El objetivo de este análisis es estudar cómo los puntajes acádemicos influyen en la probabilidad de pertenecer a cada programa, utilizando un modelo multinomial con categoría base tipo logit, donde la categoría de referencia es General.
Data
Se cargan los paquetes necesarios para manejar datos, ajustar el modelo multinomial, graficar y realizar diagnósticos.
library(vcdExtra)## Loading required package: vcd## Loading required package: grid## Loading required package: gnmlibrary(nnet)
library(haven)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:vcdExtra':
##
## summarise## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(tidyr)
setwd("C:/Users/mjsal.DESKTOP-AEUGMMS/Downloads/")
hsb <- read_dta("hsbdemo.dta")
head(hsb, 3)## # A tibble: 3 × 13
## id female ses schtyp prog read write math science socst honors
## <dbl> <dbl+lb> <dbl+l> <dbl+l> <dbl+l> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl+l>
## 1 45 1 [fema… 1 [low] 1 [pub… 3 [voc… 34 35 41 29 26 0 [not…
## 2 108 0 [male] 2 [mid… 1 [pub… 1 [gen… 34 33 41 36 36 0 [not…
## 3 15 0 [male] 3 [hig… 1 [pub… 3 [voc… 39 39 44 26 42 0 [not…
## # ℹ 2 more variables: awards <dbl>, cid <dbl>Preparación de los Datos
La variable prog viene codificada como 1,2,3. Se convierte a nombres descriptivos y a factor para usarla en el modelo.
hsb <- hsb %>%
mutate(
progr = case_when(
prog == 1 ~ "General",
prog == 2 ~ "Academic",
prog == 3 ~ "Vocation"
)
)
# Convieerto a factor y mantengo el orden
hsb$progr <- factor(hsb$progr,
levels = c("General", "Academic", "Vocation"))
table(hsb$prog, hsb$progr)##
## General Academic Vocation
## 1 45 0 0
## 2 0 105 0
## 3 0 0 50Modelo Logit Multinomial
Ajustaré el modelo usando cinco puntajes académicos como predictores: lectura, escritura, matemáticas, ciencias y estudios sociales.
mo <- multinom(
progr ~ read + write + math + science + socst,
data = hsb
)## # weights: 21 (12 variable)
## initial value 219.722458
## iter 10 value 183.838126
## iter 20 value 164.975570
## final value 164.975567
## convergedsummary(mo)## Call:
## multinom(formula = progr ~ read + write + math + science + socst,
## data = hsb)
##
## Coefficients:
## (Intercept) read write math science
## Academic -4.804712 0.045583547 0.02299571 0.09837651 -0.09359544
## Vocation 4.601228 0.009746056 -0.01443702 -0.01742029 -0.03457200
## socst
## Academic 0.03207963
## Vocation -0.03543572
##
## Std. Errors:
## (Intercept) read write math science socst
## Academic 1.498370 0.02958171 0.02958053 0.03329577 0.02985787 0.02567747
## Vocation 1.589176 0.03223947 0.03010542 0.03585979 0.02993956 0.02548155
##
## Residual Deviance: 329.9511
## AIC: 353.9511El modelo estima la probabilidad de pertenecer a cada programa comparado con la categoría base (General).
Los coeficientes positivos indican mayor probabilidad de estar en ese programa con respecto al programa General, a medida que aumenta el puntaje.
Probabilidades Ajustadas
Se varía el puntaje de lectura manteniendo los demás puntajes en su media para observar el efecto en la probabilidad de pertenecer a un programa.
read_seq <- seq(min(hsb$read), max(hsb$read), length.out = 100)
newdat <- data.frame(
read = read_seq,
write = mean(hsb$write),
math = mean(hsb$math),
science = mean(hsb$science),
socst = mean(hsb$socst)
)
probs <- as.data.frame(predict(mo, newdata = newdat, type = "probs"))Transformamos los datos para graficar:
plotdat <- bind_cols(newdat, probs) |>
pivot_longer(cols = levels(hsb$progr),
names_to = "prog",
values_to = "prob")Gráfico de probabilidades
ggplot(plotdat, aes(x = read, y = prob, colour = prog)) +
geom_line(linewidth = 1) +
labs(x = "Reading score",
y = "Fitted probability",
colour = "Program",
title = "Probabilidades Ajustadas vs Reading Score") +
theme_minimal()
A mayor puntaje en lectura, es más probable que el estudiante pertenezca al programa Academic.
Los estudiantes con baja lectura se asocian más con General.
Diagósticos Globales del Modelo
Se evalúa el ajuste mediante log-likelihood, AIC, BIC, pseudo-R² y test de razón de verosimilitudes.
logLik_fit <- logLik(mo)
AIC_fit <- AIC(mo)
BIC_fit <- BIC(mo)
monull <- multinom(progr ~ 1, data = hsb, trace = FALSE)
R2_McFadden <- 1 - (as.numeric(logLik_fit) / as.numeric(logLik(monull)))
LR <- 2 * (as.numeric(logLik_fit) - as.numeric(logLik(monull)))
df_LR <- attr(logLik(mo), "df") - attr(logLik(monull), "df")
p_LR <- pchisq(LR, df = df_LR, lower.tail = FALSE)
list(
logLik_full = as.numeric(logLik_fit),
logLik_null = as.numeric(logLik(monull)),
AIC = AIC_fit,
BIC = BIC_fit,
R2_McFadden = R2_McFadden,
LR_stat = LR,
LR_df = df_LR,
LR_p_value = p_LR
)## $logLik_full
## [1] -164.9756
##
## $logLik_null
## [1] -204.0967
##
## $AIC
## [1] 353.9511
##
## $BIC
## [1] 393.5309
##
## $R2_McFadden
## [1] 0.1916793
##
## $LR_stat
## [1] 78.24222
##
## $LR_df
## [1] 10
##
## $LR_p_value
## [1] 1.108843e-12- El pseudo-R² indica que los puntajes explican la variación entre programas.
- El test LR confirma que el modelo con covariables es significativamente mejor que el modelo nulo.
Tabla de Clasificación
Se comparan las categorías predichas con las reales.
hsb$prog_hat <- predict(mo, type = "class")
tab <- table(Observed = hsb$progr, Predicted = hsb$prog_hat)
accuracy <- mean(hsb$prog_hat == hsb$progr)
class_accuracy <- prop.table(tab, 1)
list(
confusion_matrix = tab,
overall_accuracy = accuracy,
class_accuracy = class_accuracy
)## $confusion_matrix
## Predicted
## Observed General Academic Vocation
## General 4 27 14
## Academic 4 92 9
## Vocation 6 17 27
##
## $overall_accuracy
## [1] 0.615
##
## $class_accuracy
## Predicted
## Observed General Academic Vocation
## General 0.08888889 0.60000000 0.31111111
## Academic 0.03809524 0.87619048 0.08571429
## Vocation 0.12000000 0.34000000 0.54000000- El modelo clasifica bien en general, con mejor desempeño en la categoría Academic.
Análisis de Residuales
Identificamos observaciones difíciles de predecir.
probs <- as.data.frame(predict(mo, type = "probs"))
colnames(probs) <- levels(hsb$progr)
Y <- model.matrix(~ progr - 1, data = hsb)
resid_mat <- Y - as.matrix(probs)
resid_L2 <- sqrt(rowSums(resid_mat^2))
hsb_diag <- hsb |>
mutate(resid_L2 = resid_L2)
hsb_diag |>
arrange(desc(resid_L2)) |>
select(progr, read, write, math, science, socst, resid_L2) |>
head(10)## # A tibble: 10 × 7
## progr read write math science socst resid_L2
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Vocation 63 63 75 72 66 1.33
## 2 Vocation 68 62 56 50 51 1.24
## 3 General 55 59 52 42 56 1.20
## 4 Vocation 50 67 66 66 56 1.20
## 5 Academic 39 33 38 47 41 1.18
## 6 Vocation 46 52 55 44 56 1.14
## 7 General 60 65 58 61 66 1.14
## 8 Academic 47 40 43 45 31 1.13
## 9 Vocation 68 59 53 63 61 1.13
## 10 General 57 62 56 58 66 1.11- Las observaciones con mayores residuos son estudiantes cuyo rendimiento no coincide con el patrón típico del programa.
por lo tanto…
- El rendimiento académico influye significativamente en el tipo de programa.
- Lectura es la variable más determinante.
- El modelo posee buen ajuste y una capacidad aceptable de clasificación.
- Existen casos atípicos que pueden representar perfiles educativos específicos