2025-05-16
Zamora Thowinsson Jesús David1:
info.thowinsson@gmail.com
Barranquilla, Colombia.
Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses y schtyp como variables independientes. Repita todos los análisis realizados en este documento.
Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses y honors como variables independientes. Repita todos los análisis realizados en este documento.
Considere los datos hbsdemo, tomando a prog como variable dependiente y gender, ses y awards como variables independientes. Repita todos los análisis realizados en este documento.
library(repmis) # Permite importar archivos desde internet, como Rdata o .csv
source_data("https://github.com/hllinas/DatosPublicos/blob/main/hsbdemo.Rdata?raw=TRUE") [1] "hsbdemo"#source_data del paquete repmis, que Descarga y carga directamente un archivo .Rdata desde una URL.
Datos <- hsbdemo
attach(Datos)library(knitr)
library(kableExtra)
library(dplyr)
nombres <- names(Datos) # variables en el orden del dataset
enumerados <- paste(seq_along(nombres), nombres)
mitad <- ceiling(length(enumerados) / 2)
col1 <- enumerados[1:mitad]
col2 <- enumerados[(mitad + 1):length(enumerados)]
if (length(col2) < length(col1)) {
col2 <- c(col2, rep("", length(col1) - length(col2)))
}
tabla <- data.frame(col1, col2)
kbl(tabla, col.names = c("", "")) %>%
add_header_above(c("Nombre de las Variables" = 2), bold = TRUE) %>%
kable_classic(full_width = FALSE, html_font = "Arial") %>%
column_spec(1, extra_css = "text-align: center;") %>%
column_spec(2, extra_css = "text-align: center;") %>%
row_spec(0, bold = TRUE, background = "#D3D3D3")| 1 Student | 10 science |
| 2 id | 11 socst |
| 3 gender | 12 honors |
| 4 ses | 13 awards |
| 5 schtyp | 14 cid |
| 6 prog | 15 prog0 |
| 7 read | 16 prog1 |
| 8 write | 17 prog2 |
| 9 math |
summary(Datos) Student id gender ses
Min. : 1.00 Min. : 1.00 Length:200 Length:200
1st Qu.: 50.75 1st Qu.: 50.75 Class :character Class :character
Median :100.50 Median :100.50 Mode :character Mode :character
Mean :100.50 Mean :100.50
3rd Qu.:150.25 3rd Qu.:150.25
Max. :200.00 Max. :200.00
schtyp prog read write
Length:200 Length:200 Min. :28.00 Min. :31.00
Class :character Class :character 1st Qu.:44.00 1st Qu.:45.75
Mode :character Mode :character Median :50.00 Median :54.00
Mean :52.23 Mean :52.77
3rd Qu.:60.00 3rd Qu.:60.00
Max. :76.00 Max. :67.00
math science socst honors
Min. :33.00 Min. :26.00 Min. :26.0 Length:200
1st Qu.:45.00 1st Qu.:44.00 1st Qu.:46.0 Class :character
Median :52.00 Median :53.00 Median :52.0 Mode :character
Mean :52.65 Mean :51.85 Mean :52.4
3rd Qu.:59.00 3rd Qu.:58.00 3rd Qu.:61.0
Max. :75.00 Max. :74.00 Max. :71.0
awards cid prog0 prog1 prog2
Min. :0.00 Min. : 1.00 Min. :0.000 Min. :0.00 Min. :0.000
1st Qu.:0.00 1st Qu.: 5.00 1st Qu.:0.000 1st Qu.:0.00 1st Qu.:0.000
Median :1.00 Median :10.50 Median :0.000 Median :0.00 Median :1.000
Mean :1.67 Mean :10.43 Mean :0.225 Mean :0.25 Mean :0.525
3rd Qu.:2.00 3rd Qu.:15.00 3rd Qu.:0.000 3rd Qu.:0.25 3rd Qu.:1.000
Max. :7.00 Max. :20.00 Max. :1.000 Max. :1.00 Max. :1.000 Organizando las Variables en Tablas de Mapeo de Variablesx
| Variable | female | male |
|---|---|---|
| gender | 1 | 0 |
| 0 | 1 |
| Variable | Enrolled | Not.Enrolled |
|---|---|---|
| Honors | 0 | 1 |
| 1 | 0 |
| Variable | high | low | middle |
|---|---|---|---|
| ses | 0 | 0 | 0 |
| 0 | 1 | 0 | |
| 0 | 0 | 1 |
| Variable | academic | general | vocation |
|---|---|---|---|
| prog | 0 | 0 | 0 |
| 0 | 1 | 0 | |
| 0 | 0 | 1 |
Los vectores de parámetros:
\[ \boldsymbol{\alpha} = \underbrace{(\delta_0, \beta_{0ses=low}, \beta_{0ses=middle},\beta_{0gender=male}, \beta_{0honors=notenrolled}}_{\text{Primer conjunto de parámetros}}, \underbrace{\delta_1, \beta_{1ses=low}, \beta_{1ses=middle},\beta_{1gender=male}, \beta_{1honors=noenrolled})^T}_{\text{Segundo conjunto de parámetros}} \]
y sus estimadores:
\[ \hat{\boldsymbol{\alpha}} = \underbrace{(\hat{\delta_0}, \hat{\beta}_{0ses=low}, \hat{\beta}_{0ses=middle},\hat{\beta}_{0gender=male}, \hat{\beta}_{0honors=enrolled}}_{\text{Primer conjunto de parámetros estimados}}, \underbrace{\hat{\delta_1}, \hat{\beta}_{1ses=low}, \hat{\beta}_{1ses=middle},\hat{\beta}_{1gender=male}, \hat{\beta}_{1honors=enrolled})^T}_{\text{Segundo conjunto de parámetros estimados}} \]
\[ \hat{p}_{rj} = P(\text{prog} = r \mid \text{ses} = s_j, \text{gender} = g_j, \text{honors} = h_j, \text{schtyp} = sch_j) \]
Variables:
Las variables independientes que usamos para predecir el programa educativo son ses, gender, honors, schtyp.
ses = \(s_j\): Esto indica el nivel socioeconómico del individuo \(j\). \(s_j\) tomará uno de los valores categóricos definidos para ses (e.g., “low”, “middle”, “high”).
gender = \(g_j\): Esto indica el género del individuo \(j\). \(g_j\) tomará uno de los valores categóricos definidos para gender (e.g., “male”, “female”).
honors = \(hr_j\): Esto indica si el individuo \(j\) ha sido inscrito en honores. \(h_j\) tomará uno de los valores categóricos definidos para honors (e.g., “no”, “yes”).
schtyp = \(sch_j\): Esto indica el tipo de escuela del individuo \(j\). Es una variable categórica con dos niveles (e.g., “private”, “public”)
La Probabilidad de Pertenecer a una Categoría Específica
La probabilidad de que un individuo pertenezca a un programa educativo específico \(r\), dadas sus características (\(ses_j\), \(gender_j\), \(honors_j\) \(schtyp\)), se puede expresar formalmente como:
\[ g_{r_j} := \delta_r + \beta_{r_1} {s_j} + \beta_{r_2} {g_j} + \beta_{r_3} {hr_j} + \beta_{r_4} {sch}_j \]
Donde:
Sabiendo que \(\hat{p}_{r_j}\) es como en el inciso anterior, el modelo estimado se puede escribir teniendo en cuenta la ecuación (10.1) o la (11.1):
\[ \text{Logit}(\hat{p}_{r_j}) := \ln \left( \frac{\hat{p}_{r_j}}{1 - \hat{p}_{r_j}} \right) = \hat{g}_{r_j} \] y
\[ \hat{p}_{r_j} = \frac{e^{\hat{g}_{r_j}}}{1 + e^{\hat{g}_{0_j}} + e^{\hat{g}_{1_j}}} \]
library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + gender + honors, data=Datos)# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergedsummary(modelo)Call:
multinom(formula = prog ~ ses + gender + honors, data = Datos)
Coefficients:
(Intercept) seslow sesmiddle gendermale honorsnot enrolled
general -2.237413 1.231003 0.5148617 0.07539753 1.079877
vocation -2.574235 1.152161 1.1710126 -0.04301746 1.277882
Std. Errors:
(Intercept) seslow sesmiddle gendermale honorsnot enrolled
general 0.5200380 0.5185866 0.4707092 0.3815807 0.4780577
vocation 0.5682089 0.5679726 0.4851326 0.3700698 0.4971312
Residual Deviance: 380.4986
AIC: 400.4986 library(VGAM)
modelo <- vglm(prog ~ ses + gender + honors, multinomial(refLevel = 1), data=Datos)
summary(modelo)
Call:
vglm(formula = prog ~ ses + gender + honors, family = multinomial(refLevel = 1),
data = Datos)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -2.23736 0.52003 -4.302 1.69e-05 ***
(Intercept):2 -2.57432 0.56822 -4.530 5.89e-06 ***
seslow:1 1.23101 0.51858 2.374 0.0176 *
seslow:2 1.15210 0.56798 2.028 0.0425 *
sesmiddle:1 0.51486 0.47071 1.094 0.2740
sesmiddle:2 1.17100 0.48513 2.414 0.0158 *
gendermale:1 0.07536 0.38158 0.197 0.8434
gendermale:2 -0.04306 0.37007 -0.116 0.9074
honorsnot enrolled:1 1.07981 0.47805 2.259 0.0239 *
honorsnot enrolled:2 1.27800 0.49715 2.571 0.0102 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
Residual deviance: 380.4986 on 390 degrees of freedom
Log-likelihood: -190.2493 on 390 degrees of freedom
Number of Fisher scoring iterations: 4
No Hauck-Donner effect found in any of the estimates
Reference group is level 1 of the responseSabemos que:
\[ g_{r_j} := \delta_r + \beta_{r_1} s_j + \beta_{r_2} g_j + \beta_{r_3} hr_j \]
donde \(s_j\) es un posible valor de la variable \(\text{ses}\), \(g_j\) es un posible valor de la variable \(\text{gender}\) y \(hr_j\) es un posiblevalor de la variable \(\text{honors}\). Sabemos que \(\hat{p}_{r_j}\) es como en el inciso anterior, el modelo estimado se puede escribir utilizando la ecuación (7.1) o la (8.1):
\[ \text{Logit}(\hat{p}_{r_j}) := \ln \left( \frac{\hat{p}_{r_j}}{1 - \hat{p}_{r_j}} \right) = \hat{g}_{r_j} \]
y
\[ \hat{p}_{r_j} = \frac{e^{\hat{g}_{r_j}}}{1 + e^{\hat{g}_{0_j}} + e^{\hat{g}_{1_j}} + e^{\hat{g}_{2_j}}} \]
El modelo estimado se puede escribir utilizando una de las dos expresiones anteriores. Por ejemplo, para ses=low:
\[ general: \hat{g_{0_j}} = \hat{\delta_{0_j}} + \hat{\beta_{0_1}} sl_j + \hat{\beta_{0_1}} sm_j + \hat{\beta_{0_2}} g_j + \hat{\beta_{0_3}} hr_j = -2.237413 + 1.231003(1) + 0.5148617(0) + 0.07539753(1) + 1.079877(1) = 0.14886453 \]
\[ vocation: \hat{g_{1_j}} = \hat{\delta_{1_j}} + \hat{\beta_{1_1}} sl_j + \hat{\beta_{1_1}} sm_j + \hat{\beta_{1_2}} g_j + \hat{\beta_{1_3}} hr_j = -2.574235 + 1.152161(1) + 1.1710126(0) -0.04301746 + 1.277882 = - 0.18720946 \]
library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ 1, data=Datos)# weights: 6 (2 variable)
initial value 219.722458
final value 204.096674
convergedsummary(modelo)Call:
multinom(formula = prog ~ 1, data = Datos)
Coefficients:
(Intercept)
general -0.8472980
vocation -0.7419374
Std. Errors:
(Intercept)
general 0.1781742
vocation 0.1718249
Residual Deviance: 408.1933
AIC: 412.1933 c <- summary(modelo)$coefficients[1]
d <- summary(modelo)$coefficients[2]
p0 <- exp(c)/(1+ exp(c)+ exp(d))
p0[1] 0.225p1 <- exp(d)/(1+ exp(c)+ exp(d))
p1[1] 0.25álculo de las probabilidades
Por lo tanto, las ecuaciones correspondientes para las probabilidades estimadas son:
\[ \hat{p}_{0_j} = \frac{e^{\hat{g}_{0_j}}}{1 + e^{\hat{g}_{0_j}} + e^{\hat{g}_{1_j}}} \]
y
\[ \hat{p}_{1_j} = \frac{e^{\hat{g}_{1_j}}}{1 + e^{\hat{g}_{0_j}} + e^{\hat{g}_{1_j}}} \]
library(nnet)
library(ggplot2)
# Ajuste del modelo multinomial
prog <- relevel(as.factor(prog), ref = "academic") # Relevel para tener "academic" como referencia
modelo <- multinom(prog ~ ses + gender + honors, data = Datos) # Ajuste del modelo# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergeds <- summary(modelo)
# Crear un conjunto de datos para realizar predicciones
gender_values <- c("male", "female") # Valores para gender
honors_values <- c("not enrolled", "enrolled") # Valores para honors
ses_values <- c("low", "middle", "high") # Posibles valores para ses
# Inicialización de los resultados de las probabilidades
probabilidades <- data.frame(gender = character(0), honors = character(0), ses = character(0), p_0 = numeric(0), p_1 = numeric(0), p_2 = numeric(0))
# Predecir para cada combinación de gender, honors y ses
for (g in gender_values) {
for (s_type in honors_values) {
for (ses_value in ses_values) { # Loop adicional para ses
# Codificar el valor de ses_value como un factor
ses_factor <- factor(ses_value, levels = c("low", "middle", "high"))
# Cálculos para las predicciones
g_0 <- s$coefficients[1, 1] + s$coefficients[1, 2] * (ses_value == "low") + s$coefficients[1, 3] * (ses_value == "middle") + s$coefficients[1, 4] * (g == "male") + s$coefficients[1, 5] * (s_type == "enrolled")
g_1 <- s$coefficients[2, 1] + s$coefficients[2, 2] * (ses_value == "low") + s$coefficients[2, 3] * (ses_value == "middle") + s$coefficients[2, 4] * (g == "male") + s$coefficients[2, 5] * (s_type == "enrolled")
# Probabilidades de pertenecer a los programas
p_0 <- exp(g_0) / (1 + exp(g_0) + exp(g_1)) # Probabilidad para el programa "general"
p_1 <- exp(g_1) / (1 + exp(g_0) + exp(g_1)) # Probabilidad para el programa "vocational"
p_2 <- 1 - (p_0 + p_1) # Probabilidad para el programa "academic"
# Guardar los resultados en el data frame
probabilidades <- rbind(probabilidades, data.frame(gender = g, honors = s_type, ses = ses_value, p_0 = p_0, p_1 = p_1, p_2 = p_2))
}
}
}
# Graficar los resultados
ggplot(probabilidades) +
geom_point(mapping = aes(x = ses, y = p_0, color = "p_0"), size = 1.5) +
geom_point(mapping = aes(x = ses, y = p_1, color = "p_1"), size = 1.5) +
geom_point(mapping = aes(x = ses, y = p_2, color = "p_2"), size = 1.5) +
labs(x = "Nivel Socioeconómico (ses)", y = "Probabilidad de exito") +
facet_wrap(gender ~ honors + ses) + # Facet por gender, honors y ses
theme_bw(base_size = 12) +
scale_color_discrete(name = "Tipo de programa", labels = c("(a) General", "(b) Vocational", "(c) Academic"))
Análisis de la gráfica
Estructura general de la gráfica
La gráfica está organizada en 12 paneles, cada uno representando una combinación de:
female o male.enrolled (inscrito) o not enrolled (no inscrito).high, middle, o low.Cada panel muestra:
ses).Los puntos representan las probabilidades de éxito para tres tipos de programas:
Variables clave
a) Género
female (mujeres)male (hombres)b) Estado de inscripción
enrolled: El individuo está inscrito en algún programa.not enrolled: El individuo no está inscrito en ningún programa.c) Nivel socioeconómico (SES)
high: Alto nivel socioeconómico.middle: Medio nivel socioeconómico.low: Bajo nivel socioeconómico.d) Tipo de programa
3. Interpretación de los paneles
Cada panel representa una combinación específica de las variables mencionadas. Por ejemplo:
Panel superior izquierdo
femaleenrolledhighObservaciones: - Las probabilidades de éxito varían según el tipo de programa: - Programas académicos (azul) tienen probabilidades más altas. - Programas generales (rojo) tienen probabilidades más bajas. - Para mujeres con alto SES inscritas, los programas académicos parecen ser más exitosos.
Panel inferior derecho
malenot enrolledmiddleObservaciones: - Las probabilidades de éxito son más bajas en general. - Los programas académicos (azul) aún muestran mejores resultados, aunque todos los valores son bajos.
4. Tendencias globales
Al observar todos los paneles:
Inscripción vs. No inscripción
enrolled) que para personas no inscritas (not enrolled).Nivel socioeconómico
high) y más bajas para personas con bajo SES (low).Tipo de programa
azul) suelen tener las mayores probabilidades de éxito.rojo) tienden a tener las menores probabilidades de éxito.verde) están en un punto intermedio.Diferencias por género
5. Conclusiones principales
Factores clave para el éxito
Recomendaciones
Conclusión general:
El nivel socioeconómico (ses) tiene un gran impacto en las probabilidades de asignación a diferentes tipos de programas educativos, con los individuos de niveles socioeconómicos bajos mostrando una mayor probabilidad de estar en programas “General”.
El género también tiene una ligera influencia, pero las diferencias no son tan notorias como en el caso del nivel socioeconómico.
El tipo de escuela (not enrolled vs enrolled) parece tener una influencia moderada en las probabilidades, especialmente en el programa “Academic” (azul), donde las escuelas privadas tienden a tener mayores probabilidades de asignación.
En resumen, las variables nivel socioeconómico y tipo de escuela son las que más influyen en la asignación a los programas educativos, con el género mostrando menos variabilidad en las probabilidades de éxito para los diferentes tipos de programas.
library(nnet)
library(ggplot2)
library(dplyr)
library(tidyr)
# Asegurar que todas las variables relevantes son factores
Datos <- Datos %>%
mutate(across(c(prog, ses, gender, honors), as.factor))
# Ajustar modelo multinomial
modelo <- multinom(prog ~ ses + gender + honors, data = Datos)# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
converged# Crear combinaciones posibles
nuevos_datos <- expand.grid(
ses = levels(Datos$ses),
gender = levels(Datos$gender),
honors = levels(Datos$honors)
)
# Convertir columnas a factores con los niveles del dataset original
nuevos_datos <- nuevos_datos %>%
mutate(across(everything(), ~ factor(.x, levels = levels(Datos[[cur_column()]]))))
# Predecir probabilidades
probabilidades <- predict(modelo, newdata = nuevos_datos, type = "probs")
# Unir predicciones con combinaciones y mantener columnas limpias
df_probs <- bind_cols(nuevos_datos, as.data.frame(probabilidades))
# Pivotear solo las columnas de las probabilidades
df_probs_largo <- df_probs %>%
pivot_longer(
cols = colnames(probabilidades), # Solo columnas de predicción
names_to = "Programa",
values_to = "Probabilidad"
)
# Crear una etiqueta de combinación de factores para el eje X
df_probs_largo <- df_probs_largo %>%
unite("Combinacion", ses, gender, honors, sep = "_")
# Graficar
ggplot(df_probs_largo, aes(x = Combinacion, y = Probabilidad, fill = Programa)) +
geom_bar(stat = "identity", position = "dodge") +
labs(x = "Combinación SES-Gender-Honors", y = "Probabilidad estimada",
title = "Predicción de éxito por tipo de programa y categorías") +
theme_bw(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + gender + honors, data = Datos)# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergedlibrary(VGAM)
modelo <- vglm(prog ~ ses + gender + honors, multinomial(refLevel = 1), data=Datos)
exp(coef(modelo)) (Intercept):1 (Intercept):2 seslow:1
0.10674006 0.07620594 3.42470121
seslow:2 sesmiddle:1 sesmiddle:2
3.16483231 1.67340902 3.22523080
gendermale:1 gendermale:2 honorsnot enrolled:1
1.07827231 0.95785809 2.94410867
honorsnot enrolled:2
3.58944029 Conclusión general: Nivel socioeconómico tiene un efecto significativo en las probabilidades de pertenecer a las categorías de programas educativos, especialmente para los niveles bajo y medio. Los individuos con un nivel socioeconómico más bajo o medio tienen una probabilidad más alta de estar en categorías no académicas, como “vocational” o “general”.
Género también juega un papel, pero el efecto no es tan pronunciado, ya que el coeficiente gendermale:1 no es muy grande, lo que sugiere que el género masculino tiene una probabilidad ligeramente mayor de estar en las categorías no académicas en comparación con las mujeres.
Tipo de escuela (pública vs. privada) también tiene un impacto significativo en las probabilidades de pertenecer a los diferentes programas. Las escuelas públicas parecen estar asociadas con una probabilidad más alta de que los individuos pertenezcan a categorías diferentes de “academic”.
En resumen, el nivel socioeconómico es el factor que parece tener el mayor impacto en la asignación a las categorías de prog, seguido por el género y el tipo de escuela.
head(fitted(modelo)) academic general vocation
1 0.3399135 0.3658242 0.2942623
2 0.4145803 0.2350826 0.3503371
3 0.6246638 0.2116684 0.1636677
4 0.3344838 0.3881571 0.2773591
5 0.4145803 0.2350826 0.3503371
6 0.6298058 0.1979192 0.1722750# Codificar 'gender' como variable numérica
Datos$gender_numeric <- ifelse(Datos$gender == "female", 0, 1)
# Calcular el promedio de 'gender' codificado
mean_gender <- mean(Datos$gender_numeric, na.rm = TRUE)
# Crear el dataframe 'Nuevo' con el promedio de 'gender' y otros valores
Nuevo <- data.frame(
ses = c("low", "middle", "high"),
gender = mean_gender
)
# Ver el resultado
print(Nuevo) ses gender
1 low 0.455
2 middle 0.455
3 high 0.455gender fue codificado con 0 para “female” y 1 para “male”. Un valor de 0.455 indica que, aproximadamente, el 45.5% de los individuos en cada nivel socioeconómico son hombres (gender = 1), y el 54.5% son mujeres (gender = 0).
El valor 0.455 es el promedio de las codificaciones de género para cada categoría de ses, lo que sugiere que la proporción de hombres y mujeres es bastante similar en los tres niveles socioeconómicos (low, middle, y high).
En cada nivel socioeconómico, “low”, “middle”, y “high”, la proporción de hombres y mujeres es casi igual, con un valor promedio de 0.455 para gender en todos los niveles de ses.
El promedio 0.455 es el resultado de codificar “female” como 0 y “male” como 1, y muestra una ligera inclinación hacia “female”, ya que el valor está más cerca de 0 que de 1.
library(nnet)
library(ggplot2)
# Ajuste del modelo multinomial
prog <- relevel(as.factor(prog), ref = "academic") # Relevel para tener "academic" como referencia
modelo <- multinom(prog ~ ses + gender + honors, data = Datos) # Ajuste del modelo# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergeds <- summary(modelo)
# Crear un conjunto de datos para realizar predicciones
gender_values <- c("male", "female") # Valores para gender
honors_values <- c("not enrolled", "enrolled") # Valores para honors
ses_values <- c("low", "middle", "high") # Posibles valores para ses
# Inicialización de los resultados de las probabilidades
probabilidades <- data.frame(gender = character(0), honors = character(0), ses = character(0), p_0 = numeric(0), p_1 = numeric(0), p_2 = numeric(0))
# Predecir para cada combinación de gender, honors y ses
for (g in gender_values) {
for (s_type in honors_values) {
for (ses_value in ses_values) { # Loop adicional para ses
# Codificar el valor de ses_value como un factor
ses_factor <- factor(ses_value, levels = c("low", "middle", "high"))
# Cálculos para las predicciones
g_0 <- s$coefficients[1, 1] + s$coefficients[1, 2] * (ses_value == "low") + s$coefficients[1, 3] * (ses_value == "middle") + s$coefficients[1, 4] * (g == "male") + s$coefficients[1, 5] * (s_type == "enrolled")
g_1 <- s$coefficients[2, 1] + s$coefficients[2, 2] * (ses_value == "low") + s$coefficients[2, 3] * (ses_value == "middle") + s$coefficients[2, 4] * (g == "male") + s$coefficients[2, 5] * (s_type == "enrolled")
# Probabilidades de pertenecer a los programas
p_0 <- exp(g_0) / (1 + exp(g_0) + exp(g_1)) # Probabilidad para el programa "general"
p_1 <- exp(g_1) / (1 + exp(g_0) + exp(g_1)) # Probabilidad para el programa "vocational"
p_2 <- 1 - (p_0 + p_1) # Probabilidad para el programa "academic"
# Guardar los resultados en el data frame
probabilidades <- rbind(probabilidades, data.frame(gender = g, honors = s_type, ses = ses_value, p_0 = p_0, p_1 = p_1, p_2 = p_2))
}
}
}
predichos <- c(p_0,p_1,p_2)
predichos[1] 0.1979231 0.1722697 0.6298072Conclusión: El modelo ha convergido correctamente y ha ajustado los coeficientes,Las probabilidades predichas para cada tipo de programa (general, vocational, academic) son: General: 15.03% Vocational: 14.01% Academic: 70.96%
Esto indica que, manteniendo gender fijo en su media, la probabilidad más alta es que el individuo pertenezca al programa “academic”.
str(hsbdemo$gender) chr [1:200] "female" "male" "male" "male" "male" "female" "male" "male" ...str(hsbdemo$honors) chr [1:200] "not enrolled" "not enrolled" "not enrolled" "not enrolled" ...# Convertir 'gender' y 'honors' a factores si no lo son
hsbdemo$gender <- factor(hsbdemo$gender, levels = c("female", "male"))
hsbdemo$honors <- factor(hsbdemo$honors, levels = c("enrolled", "not enrolled"))# Paso 1: Crear el vector para 'ses' con 200 elementos
ses_values <- rep(c("low", "middle", "high"), length.out = 200) # Esto asegura 200 elementos
# Paso 2: Crear el vector para 'gender' con 200 elementos (ya lo tienes como secuencia continua o categórica)
gender_values <- rep(c("male", "female"), length.out = 200) # Esto asegura 200 elementos para 'gender'
# Paso 3: Crear el data frame 'Nuevo' con 200 filas
Nuevo <- data.frame(
ses = ses_values, # Valores de 'ses' con 200 elementos
gender = gender_values, # Valores de 'gender' con 200 elementos
honors = rep(c("enrolled", "not enrolled"), length.out = 200) # Combinación de 'enrolled' y 'not enrolled' para 200 filas
)
# Verifica la estructura del nuevo data frame 'Nuevo'
str(Nuevo)'data.frame': 200 obs. of 3 variables:
$ ses : chr "low" "middle" "high" "low" ...
$ gender: chr "male" "female" "male" "female" ...
$ honors: chr "enrolled" "not enrolled" "enrolled" "not enrolled" ...# Paso 1: Crear un vector para 'ses' con 200 elementos
# Usamos 'length.out = 200' para asegurar que tenga exactamente 200 filas
ses_values <- rep(c("low", "middle", "high"), length.out = 200) # Repetir 'ses' para 200 filas
# Paso 2: Crear un vector para 'gender' con 200 elementos
# Usamos 'length.out = 200' para asegurar que tenga exactamente 200 filas
gender_values <- rep(c("male", "female"), length.out = 200) # Repetir 'gender' para 200 filas
# Paso 3: Crear un vector para 'honors' con 200 elementos
# Usamos 'length.out = 200' para asegurar que tenga exactamente 200 filas
honors_values <- rep(c("enrolled", "not enrolled"), length.out = 200) # Repetir 'honors' para 200 filas
# Crear el data frame 'Nuevo' con las combinaciones de 'ses', 'gender' y 'honors'
Nuevo <- data.frame(
ses = ses_values, # Valores de 'ses' con 200 elementos
gender = gender_values, # Valores de 'gender' con 200 elementos
honors = honors_values # Valores de 'honors' con 200 elementos
)
# Verifica la estructura del nuevo data frame 'Nuevo'
str(Nuevo)'data.frame': 200 obs. of 3 variables:
$ ses : chr "low" "middle" "high" "low" ...
$ gender: chr "male" "female" "male" "female" ...
$ honors: chr "enrolled" "not enrolled" "enrolled" "not enrolled" ...# Paso 2: Realizar las predicciones para los nuevos datos
# Realizamos las predicciones utilizando el modelo ajustado
Predichos <- cbind(Nuevo, predict(modelo, newdata = Nuevo, type = "probs"))
# Ver las primeras predicciones
head(Predichos)# Paso 3: Agrupar por 'ses' y calcular la media de las probabilidades para cada tipo de programa
# Usamos dplyr para agrupar y calcular las medias de las probabilidades
library(dplyr)
Resumen <- Predichos %>%
group_by(ses) %>% # Agrupar por 'ses' (low, middle, high)
summarise(
P0_general = mean(general), # Promedio de la probabilidad de "general"
P1_vocation = mean(vocation), # Promedio de la probabilidad de "vocational"
P2_academic = mean(academic) # Promedio de la probabilidad de "academic"
)
# Ver el resumen de las probabilidades predichas
print(Resumen)# A tibble: 3 × 4
ses P0_general P1_vocation P2_academic
<chr> <dbl> <dbl> <dbl>
1 high 0.147 0.117 0.736
2 low 0.303 0.217 0.480
3 middle 0.177 0.267 0.556El nivel socioeconómico parece tener una influencia importante en las probabilidades de que un individuo participe en programas académicos versus generales o vocacionales. Los individuos con alto nivel socioeconómico tienen una probabilidad mucho mayor de estar en un programa académico, mientras que aquellos con bajo nivel socioeconómico tienden a estar más distribuidos en programas generales o vocacionales.
Esta interpretación puede variar dependiendo de otros factores que no están reflejados en esta tabla, pero proporciona un buen punto de partida para analizar la relación entre el nivel socioeconómico y la elección del programa educativo.
library(ggplot2)
library(reshape2)
# Suponiendo que 'Predichos' es tu data frame original
Tabla <- melt(Predichos, id.vars = c("ses", "gender"), value.name = "probability")
ggplot(Tabla, aes(x = gender, y = probability, group = ses, colour = ses)) +
geom_line(size = 1) +
geom_point(size = 3) +
facet_wrap(~ variable, scales = "free_y") +
theme_minimal() +
labs(title = "Probabilidades por Genero y Nivel Socioeconómico",
x = "Género", y = "Probabilidad", colour = "Nivel SES") +
theme(strip.text = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
legend.position = "bottom")
Conclusión general: El nivel socioeconómico influye fuertemente en el tipo de programa educativo al que acceden los estudiantes. El género también influye, especialmente en el tipo de escuela (privada vs. pública). Las mujeres de nivel alto tienen mayores oportunidades de acceder a programas académicos, mientras que las personas de nivel bajo, sin importar el género, tienden a estar en programas generales.
library(nnet)
modelo <- multinom(prog ~ ses + gender + honors, data = hsbdemo)# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergedsummary_modelo <- summary(modelo)
std_errors <- summary_modelo$standard.errors
print(std_errors) (Intercept) seslow sesmiddle gendermale honorsnot enrolled
general 0.5200380 0.5185866 0.4707092 0.3815807 0.4780577
vocation 0.5682089 0.5679726 0.4851326 0.3700698 0.4971312library(VGAM)
modelo <- vglm(prog ~ ses + gender + honors, family = multinomial(refLevel = 1), data = Datos)
# Obtener el resumen del modelo
summary_modelo <- summary(modelo)
summary_modelo
Call:
vglm(formula = prog ~ ses + gender + honors, family = multinomial(refLevel = 1),
data = Datos)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -2.23736 0.52003 -4.302 1.69e-05 ***
(Intercept):2 -2.57432 0.56822 -4.530 5.89e-06 ***
seslow:1 1.23101 0.51858 2.374 0.0176 *
seslow:2 1.15210 0.56798 2.028 0.0425 *
sesmiddle:1 0.51486 0.47071 1.094 0.2740
sesmiddle:2 1.17100 0.48513 2.414 0.0158 *
gendermale:1 0.07536 0.38158 0.197 0.8434
gendermale:2 -0.04306 0.37007 -0.116 0.9074
honorsnot enrolled:1 1.07981 0.47805 2.259 0.0239 *
honorsnot enrolled:2 1.27800 0.49715 2.571 0.0102 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
Residual deviance: 380.4986 on 390 degrees of freedom
Log-likelihood: -190.2493 on 390 degrees of freedom
Number of Fisher scoring iterations: 4
No Hauck-Donner effect found in any of the estimates
Reference group is level 1 of the responseLa estimación del logaritmo de la función de máxima verosimilitud en el modelo logístico se puede obtener de varias maneras. Una es reemplazando los valores correspondientes en la ecuación (8.2):
\[ 𝓛(\hat{\alpha}) = \sum_{j=1}^J \left[z_{0_j} g^{0_j} + (n_j - z_{0_j} - z_{1_j}) g^{1_j} - n_j \ln(1 + e^{g^{0_j}} + e^{g^{1_j}})\right] = -190.1483 \]
donde,
\[ g^{r_j} := \delta^r + \beta^r_1 x_{j1} + \cdots + \beta^r_K x_{jK} \]
Las otras maneras para calcular \(𝓛(\hat{\alpha})\) son con las funciones multinom::nnet o vglm::VGAM. Con cualquier camino encontramos que
\[ 𝓛(\hat{\alpha}) = -190.1483 \]
library(nnet)
prog <- relevel(as.factor(prog), ref = "academic")
modelo <- multinom(prog ~ ses + gender + honors, data=Datos)# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergedsummary(modelo)Call:
multinom(formula = prog ~ ses + gender + honors, data = Datos)
Coefficients:
(Intercept) seslow sesmiddle gendermale honorsnot enrolled
general -2.237413 1.231003 0.5148617 0.07539753 1.079877
vocation -2.574235 1.152161 1.1710126 -0.04301746 1.277882
Std. Errors:
(Intercept) seslow sesmiddle gendermale honorsnot enrolled
general 0.5200380 0.5185866 0.4707092 0.3815807 0.4780577
vocation 0.5682089 0.5679726 0.4851326 0.3700698 0.4971312
Residual Deviance: 380.4986
AIC: 400.4986 library(VGAM)
modelo <- vglm(prog ~ ses + gender + honors, multinomial(refLevel = 1), data=Datos)
summary(modelo)
Call:
vglm(formula = prog ~ ses + gender + honors, family = multinomial(refLevel = 1),
data = Datos)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -2.23736 0.52003 -4.302 1.69e-05 ***
(Intercept):2 -2.57432 0.56822 -4.530 5.89e-06 ***
seslow:1 1.23101 0.51858 2.374 0.0176 *
seslow:2 1.15210 0.56798 2.028 0.0425 *
sesmiddle:1 0.51486 0.47071 1.094 0.2740
sesmiddle:2 1.17100 0.48513 2.414 0.0158 *
gendermale:1 0.07536 0.38158 0.197 0.8434
gendermale:2 -0.04306 0.37007 -0.116 0.9074
honorsnot enrolled:1 1.07981 0.47805 2.259 0.0239 *
honorsnot enrolled:2 1.27800 0.49715 2.571 0.0102 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
Residual deviance: 380.4986 on 390 degrees of freedom
Log-likelihood: -190.2493 on 390 degrees of freedom
Number of Fisher scoring iterations: 4
No Hauck-Donner effect found in any of the estimates
Reference group is level 1 of the responseCalcular Intervalos confianza sea \(\text{logit}(p_j)\) o esccala de probabilidad
Cálculo de la Varianza para el Estimador del Logit
La varianza de \(\text{logit}(p_j)\) se puede calcular usando la matriz de covarianzas de los coeficientes \(\hat{\beta}\). Dado que tenemos dos variables independientes, la matriz de covarianzas tendrá la forma:
\[ \hat{\text{Cov}}(\hat{\beta}) = \begin{pmatrix} \hat{\text{Var}}(\hat{\beta_0}) & \hat{\text{Cov}}(\hat{\beta_0}, \hat{\beta_1}) & \hat{\text{Cov}}(\hat{\beta_0}, \hat{\beta_2}) \\ \hat{\text{Cov}}(\hat{\beta_1}, \hat{\beta_0}) & \hat{\text{Var}}(\hat{\beta_1}) & \hat{\text{Cov}}(\hat{\beta_1}, \hat{\beta_2}) \\ \hat{\text{Cov}}(\hat{\beta_2}, \hat{\beta_0}) & \hat{\text{Cov}}(\hat{\beta_2}, \hat{\beta_1}) & \hat{\text{Var}}(\hat{\beta_2}) \end{pmatrix} \]
Para simplificar, podemos asumir que ya tenemos la matriz de covarianzas \(\text{Cov}(\hat{\beta})\). Luego, la varianza del logit para la categoría de café (por ejemplo) se calcularía de la siguiente manera:
\[ \text{Var}(\text{logit}(p_{\text{café}})) = x_1^2 \cdot \text{Var}(\hat{\beta_1}) + x_2^2 \cdot \text{Var}(\hat{\beta_2}) + 2 \cdot x_1 \cdot x_2 \cdot \text{Cov}(\hat{\beta_1}, \hat{\beta_2}) \]
Y para la categoría de jugo:
\[ \text{Var}(\text{logit}(p_{\text{jugo}})) = x_1^2 \cdot \text{Var}(\hat{\beta_1}) + x_2^2 \cdot \text{Var}(\hat{\beta_2}) + 2 \cdot x_1 \cdot x_2 \cdot \text{Cov}(\hat{\beta_1}, \hat{\beta_2}) \]
modelo <- multinom(prog ~ ses + gender + honors, data=Datos)# weights: 18 (10 variable)
initial value 219.722458
iter 10 value 190.375915
final value 190.249312
convergedsummary(modelo)Call:
multinom(formula = prog ~ ses + gender + honors, data = Datos)
Coefficients:
(Intercept) seslow sesmiddle gendermale honorsnot enrolled
general -2.237413 1.231003 0.5148617 0.07539753 1.079877
vocation -2.574235 1.152161 1.1710126 -0.04301746 1.277882
Std. Errors:
(Intercept) seslow sesmiddle gendermale honorsnot enrolled
general 0.5200380 0.5185866 0.4707092 0.3815807 0.4780577
vocation 0.5682089 0.5679726 0.4851326 0.3700698 0.4971312
Residual Deviance: 380.4986
AIC: 400.4986 | Pruebas de Hipotesis | Vector de parámetros | No de parámetros | Logaritmo de la función de verosimilitud | Estadístico de prueba (𝐷) | Distribución asintótica de 𝐷. | Grado de libertad (𝜈) | P-valor 𝑃(𝜒2𝜈≥𝐷) |
|---|---|---|---|---|---|---|---|
| Modelo logístico | 𝛼̂ | 2(1+𝐾) | \(\mathcal{L}(\hat{\alpha})\) | ||||
| \(𝐻_0\): Logit vs \(H_1\): Completo | \(\hat{p} = u\) | \(2n\) | \[\mathcal{L}(\hat{\alpha}) = 0\] | \[-2\mathcal{L}(\hat{\alpha})=\] | \(\chi^2\) | \(2[n - (1 + K)] =\) | |
| Pruebas de Hipótesis | Vector de parámetros | Nº de parámetros | Logaritmo de la función de verosimilitud | Estadístico de prueba (\(D\)) | Distribución asintótica de \(D\) | Grados de libertad (\(\nu\)) | P-valor \(P(\chi^2_{\nu} \geq D)\) |
|---|---|---|---|---|---|---|---|
| Modelo logístico | \(\hat{\boldsymbol{\alpha}}\) | \(2(1+K) = 10\) | \(\mathcal{L}(\hat{\boldsymbol{\alpha}}) = -190.2493\) | ||||
| \(H_0\): Logit vs \(H_1\): Completo | \(\hat{p} = u\) | \(2n=400\) | \(\mathcal{L}(\hat{\boldsymbol{\alpha}}) = 0\) | \(-2 \log \mathcal{L}(\hat{\boldsymbol{\alpha}}) = 380.4986\) | \(\chi^2\) | \(2[n - (1 + K)] = 390\) | \(P(\chi^2_{390} \geq 380.4986) = 0.638\) |
pchisq(380.4986, df = 390, lower.tail = FALSE)[1] 0.625169Una vez que tienes la varianza del logit, el intervalo de confianza para \(\text{logit}(p_{\text{café}})\) se calcula usando una distribución normal (aproximación de máxima verosimilitud) con un valor \(z\) asociado al nivel de confianza deseado (por ejemplo, \(z = 1.96\) para un intervalo de confianza del 95%).
El intervalo de confianza para \(\text{logit}(p_{\text{café}})\) es:
\[ \hat{\text{logit}}(p_{\text{café}}) \pm z \cdot \sqrt{\text{Var}(\text{logit}(p_{\text{café}}))} \]
Estimación del logit: La estimación de \(\text{logit}(p_{\text{café}})\) para un consumidor con valores conocidos de las variables independientes (edad \(x_1\) e ingresos \(x_2\)) es:
\[ \hat{\text{logit}}(p_{\text{café}}) = \hat{\beta}_{\text{café0}} + \hat{\beta}_{\text{café1}} \cdot x_1 + \hat{\beta}_{\text{café2}} \cdot x_2 \]
Cálculo de la varianza: La varianza de \(\text{logit}(p_{\text{café}})\) es:
\[ \text{Var}(\text{logit}(p_{\text{café}})) = x_1^2 \cdot \text{Var}(\hat{\beta}_{\text{café1}}) + x_2^2 \cdot \text{Var}(\hat{\beta}_{\text{café2}}) + 2 \cdot x_1 \cdot x_2 \cdot \text{Cov}(\hat{\beta}_{\text{café1}}, \hat{\beta}_{\text{café2}}) \]
Cálculo del intervalo de confianza: Con el valor \(z = 1.96\) para un intervalo de confianza del 95%, el intervalo de confianza para \(\text{logit}(p_{\text{café}})\) es:
\[ \hat{\text{logit}}(p_{\text{café}}) \pm 1.96 \cdot \sqrt{\text{Var}(\text{logit}(p_{\text{café}}))} \]
Transformación a la escala de probabilidad: Para obtener el intervalo de confianza para la probabilidad \(p_{\text{café}}\), debemos aplicar la función logística inversa:
\[ p_{\text{café}} = \frac{1}{1 + e^{-\hat{\text{logit}}(p_{\text{café}})}} \]
OJO SE REEMPLAZA LOS VALORES DE LOS INTERVALOS EN LA FUNCION Y transformarlo a la escala de probabilidad en un modelo multinomial. lo que quiere decir va a estar entre 20% y 40%
EJEMPLO
\[ IC_{\text{logit}} = [-872.23, 1173.23] \]
\[ p_{\text{café}} = \frac{1}{1 + e^{-(-872.23)}} \quad \text{y} \quad p_{\text{café}} = \frac{1}{1 + e^{-1173.23}} \]
Jesús David Zamora Thowinsson. Economista y Administrador Público, Especialista en Estadística Aplicada, Candidato a Magíster Scientiarum en Gerencia Empresarial, Magíster en Estadística Aplicada.↩︎