Código reproducible de la tesis de Maestría en Estadística
Predicción del riesgo de deserción universitaria en estudiantes del Área de Matemáticas de la FACEN-UNA
Marlene Román
Julio 2026
Verificado con la versión final
############################################################
# TESIS DE ESTADÍSTICA - SCRIPT INTEGRADO Y DEPURADO
# Predicción del riesgo de deserción universitaria
# Área de Matemáticas - FACEN-UNA
############################################################
rm(list = ls())
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 625501 33.5 1330855 71.1 957620 51.2
## Vcells 1202296 9.2 8388608 64.0 1927711 14.8options(scipen = 999)
set.seed(123)
############################################################
# 1. LIBRERÍAS
############################################################
paquetes <- c(
"readxl", "dplyr", "ggplot2", "caret", "pROC", "rpart", "rpart.plot",
"randomForest", "e1071", "nnet", "naivebayes", "openxlsx", "tidyr",
"car", "ResourceSelection"
)
instalar_si_falta <- function(pkg) {
if (!require(pkg, character.only = TRUE)) {
install.packages(pkg, dependencies = TRUE)
library(pkg, character.only = TRUE)
}
}
invisible(lapply(paquetes, instalar_si_falta))## Cargando paquete requerido: readxl## Cargando paquete requerido: dplyr##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Cargando paquete requerido: ggplot2## Cargando paquete requerido: caret## Cargando paquete requerido: lattice## Cargando paquete requerido: pROC## Type 'citation("pROC")' for a citation.##
## Adjuntando el paquete: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var## Cargando paquete requerido: rpart## Cargando paquete requerido: rpart.plot## Cargando paquete requerido: randomForest## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.##
## Adjuntando el paquete: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## margin## The following object is masked from 'package:dplyr':
##
## combine## Cargando paquete requerido: e1071##
## Adjuntando el paquete: 'e1071'## The following object is masked from 'package:ggplot2':
##
## element## Cargando paquete requerido: nnet## Cargando paquete requerido: naivebayes## naivebayes 1.0.0 loaded## For more information please visit:## https://majkamichal.github.io/naivebayes/## Cargando paquete requerido: openxlsx## Cargando paquete requerido: tidyr## Cargando paquete requerido: car## Cargando paquete requerido: carData##
## Adjuntando el paquete: 'car'## The following object is masked from 'package:dplyr':
##
## recode## Cargando paquete requerido: ResourceSelection## ResourceSelection 0.3-6 2023-06-27############################################################
# 2. CARGA DE DATOS
############################################################
# Ajustar esta ruta según tu equipo
ruta_datos <- "desercion.xlsx"
tesis <- read_excel(ruta_datos)
tesis <- as.data.frame(tesis, stringsAsFactors = FALSE)
# Carpeta de salidas
dir.create("salidas_tesis", showWarnings = FALSE, recursive = TRUE)
############################################################
# 3. LIMPIEZA BÁSICA Y NORMALIZACIÓN DE TEXTO
############################################################
cols_char <- names(tesis)[sapply(tesis, is.character)]
tesis[cols_char] <- lapply(tesis[cols_char], function(y) {
y <- as.character(y)
y <- enc2utf8(y)
y <- iconv(y, from = "", to = "UTF-8", sub = "")
y <- trimws(y)
y[is.na(y)] <- ""
y
})
############################################################
# 4. RECODIFICACIÓN DE VARIABLES
############################################################
# Variable respuesta: 1 = Desertor ; 0 = No_Desertor
tesis$ESTADO_ACADEMICO <- factor(
tesis$ESTADO_ACADEMICO,
levels = c(1, 0),
labels = c("Desertor", "No_Desertor")
)
tesis$SEXO <- factor(
tesis$SEXO,
levels = c("F", "M"),
labels = c("FEMENINO", "MASCULINO")
)
tesis$CARRERA <- factor(
tesis$CARRERA,
levels = c("2009E", "09MED", "2009P", "010EM", "010EMD"),
labels = c(
"ESTADISTICA-PRES",
"ESTADISTICA-SEMI",
"MATEMATICA",
"EDUCACION MATEMATICA-PRES",
"EDUCACION MATEMATICA-SEMI"
)
)
tesis$TIPO_INGRESO <- factor(
tesis$TIPO_INGRESO,
levels = c("INGRESO", "TRASLADO", "ADMISION DIRECTA"),
labels = c("INGRESO", "TRASLADO", "ADMISION DIRECTA")
)
tesis$MODALIDAD <- factor(
tesis$MODALIDAD,
levels = c("Presencial", "Semipresencial"),
labels = c("Presencial", "Semipresencial")
)
tesis$RENDIMIENTO <- factor(
tesis$RENDIMIENTO,
levels = c("APROBADO", "NO APROBADO"),
labels = c("APROBADO", "NO APROBADO")
)
tesis$RESIDENCIA <- factor(
tesis$RESIDENCIA,
levels = c("ASUNCIÓN", "CENTRAL", "RESTO"),
labels = c("ASUNCIÓN", "CENTRAL", "RESTO DEL PAÍS")
)
tesis$EST_CIVIL <- factor(
tesis$EST_CIVIL,
levels = c(1, 2, 3, 4),
labels = c("Soltero/a", "Casado/a", "Divorciado/a", "Otro")
)
tesis$TIPO_COL <- factor(
tesis$TIPO_COL,
levels = c(1, 2, 3),
labels = c("Público", "Subvencionado", "Privado")
)
tesis$TRABAJA <- factor(
tesis$TRABAJA,
levels = c(1, 2),
labels = c("SI", "NO")
)
tesis$SOLVENTAR <- factor(
tesis$SOLVENTAR,
levels = c(1, 2, 3, 4),
labels = c(
"Beca/exoneración total",
"Beca/exoneración parcial",
"Trabajo Personal",
"Ayuda Familiar"
)
)
tesis$ESTUDIOS_PADRES <- factor(
tesis$ESTUDIOS_PADRES,
levels = c(1, 2, 3, 4, 5),
labels = c("hasta 13 años", "14-23 años", "24-29 años", "30-34 años", "Más de 34 años")
)
tesis$INGRESO <- factor(
tesis$INGRESO,
levels = c(1, 2, 3, 4, 5),
labels = c(
"Hasta dos salarios mínimos",
"Más de dos y hasta cinco salarios mínimos",
"Más de cinco y hasta diez salarios mínimos",
"Más de diez y hasta quince salarios mínimos",
"Más de quince salarios mínimos"
)
)
tesis$NIVEL_SOCIO <- factor(
tesis$NIVEL_SOCIO,
levels = sort(unique(tesis$NIVEL_SOCIO))
)
tesis$EDAD <- as.numeric(tesis$EDAD)
############################################################
# 5. CHEQUEOS INICIALES
############################################################
cat("\n=============================\n")##
## =============================cat("ESTRUCTURA DE LA BASE\n")## ESTRUCTURA DE LA BASEcat("=============================\n")## =============================str(tesis)## 'data.frame': 682 obs. of 15 variables:
## $ SEXO : Factor w/ 2 levels "FEMENINO","MASCULINO": 1 1 1 2 2 1 1 1 1 1 ...
## $ CARRERA : Factor w/ 5 levels "ESTADISTICA-PRES",..: 3 1 1 1 1 4 1 3 4 3 ...
## $ TIPO_INGRESO : Factor w/ 3 levels "INGRESO","TRASLADO",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ MODALIDAD : Factor w/ 2 levels "Presencial","Semipresencial": 1 1 1 1 1 1 1 1 1 1 ...
## $ ESTADO_ACADEMICO: Factor w/ 2 levels "Desertor","No_Desertor": 1 1 2 1 1 2 1 1 1 1 ...
## $ RENDIMIENTO : Factor w/ 2 levels "APROBADO","NO APROBADO": 2 1 1 1 2 1 1 1 1 1 ...
## $ RESIDENCIA : Factor w/ 3 levels "ASUNCIÓN","CENTRAL",..: 2 2 3 2 2 3 1 1 3 1 ...
## $ EDAD : num 33 19 19 20 19 21 24 20 19 21 ...
## $ EST_CIVIL : Factor w/ 4 levels "Soltero/a","Casado/a",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ TIPO_COL : Factor w/ 3 levels "Público","Subvencionado",..: 1 3 2 3 1 1 1 1 1 1 ...
## $ TRABAJA : Factor w/ 2 levels "SI","NO": 2 2 2 2 2 2 1 2 2 2 ...
## $ SOLVENTAR : Factor w/ 4 levels "Beca/exoneración total",..: 4 4 4 4 3 4 3 3 4 4 ...
## $ ESTUDIOS_PADRES : Factor w/ 5 levels "hasta 13 años",..: 1 2 2 3 1 2 3 2 2 2 ...
## $ INGRESO : Factor w/ 5 levels "Hasta dos salarios mínimos",..: 1 1 1 1 1 1 2 1 1 1 ...
## $ NIVEL_SOCIO : Factor w/ 2 levels "BAJO","MEDIO": 2 2 2 2 2 2 2 2 2 2 ...cat("\n=============================\n")##
## =============================cat("TABLA DE LA VARIABLE RESPUESTA\n")## TABLA DE LA VARIABLE RESPUESTAcat("=============================\n")## =============================print(table(tesis$ESTADO_ACADEMICO))##
## Desertor No_Desertor
## 459 223cat("\n=============================\n")##
## =============================cat("RESUMEN GENERAL\n")## RESUMEN GENERALcat("=============================\n")## =============================print(summary(tesis))## SEXO CARRERA TIPO_INGRESO
## FEMENINO :395 ESTADISTICA-PRES :182 INGRESO :642
## MASCULINO:287 ESTADISTICA-SEMI :100 TRASLADO : 0
## MATEMATICA :150 ADMISION DIRECTA: 40
## EDUCACION MATEMATICA-PRES: 70
## EDUCACION MATEMATICA-SEMI:180
##
## MODALIDAD ESTADO_ACADEMICO RENDIMIENTO
## Presencial :402 Desertor :459 APROBADO :409
## Semipresencial:280 No_Desertor:223 NO APROBADO:273
##
##
##
##
## RESIDENCIA EDAD EST_CIVIL TIPO_COL
## ASUNCIÓN :229 Min. :18.00 Soltero/a :544 Público :465
## CENTRAL :336 1st Qu.:20.00 Casado/a :117 Subvencionado: 49
## RESTO DEL PAÍS:117 Median :23.00 Divorciado/a: 14 Privado :168
## Mean :27.11 Otro : 7
## 3rd Qu.:32.00
## Max. :70.00
## TRABAJA SOLVENTAR ESTUDIOS_PADRES
## SI:372 Beca/exoneración total : 45 hasta 13 años :157
## NO:310 Beca/exoneración parcial: 49 14-23 años :215
## Trabajo Personal :346 24-29 años :164
## Ayuda Familiar :242 30-34 años :101
## Más de 34 años: 45
##
## INGRESO NIVEL_SOCIO
## Hasta dos salarios mínimos :430 BAJO : 65
## Más de dos y hasta cinco salarios mínimos :190 MEDIO:617
## Más de cinco y hasta diez salarios mínimos : 45
## Más de diez y hasta quince salarios mínimos: 12
## Más de quince salarios mínimos : 5
## ############################################################
# 6. ANÁLISIS UNIVARIANTE
############################################################
tabla_edad <- data.frame(
Variable = "Edad",
Casos = sum(!is.na(tesis$EDAD)),
Minimo = min(tesis$EDAD, na.rm = TRUE),
Q1 = quantile(tesis$EDAD, 0.25, na.rm = TRUE),
Mediana = median(tesis$EDAD, na.rm = TRUE),
Media = mean(tesis$EDAD, na.rm = TRUE),
Q3 = quantile(tesis$EDAD, 0.75, na.rm = TRUE),
Maximo = max(tesis$EDAD, na.rm = TRUE)
)
tabla_edad[, -1] <- round(tabla_edad[, -1], 2)
row.names(tabla_edad) <- NULL
print(tabla_edad)## Variable Casos Minimo Q1 Mediana Media Q3 Maximo
## 1 Edad 682 18 20 23 27.11 32 70graf_estado <- ggplot(tesis, aes(x = ESTADO_ACADEMICO, fill = ESTADO_ACADEMICO)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5) +
labs(
title = "Distribución del estado académico",
x = "Estado académico",
y = "Frecuencia"
) +
theme_minimal() +
theme(legend.position = "none")
print(graf_estado)
############################################################
# 6.1 TABLAS DESCRIPTIVAS DEL CAPÍTULO 4
############################################################
############################################################
# FUNCIONES AUXILIARES
############################################################
tabla_frecuencias <- function(data, var, nombre_variable = var) {
x <- data[[var]]
tab <- table(x, useNA = "ifany")
prop <- prop.table(tab) * 100
data.frame(
Variable = nombre_variable,
Categoria = names(tab),
Frecuencia = as.numeric(tab),
Porcentaje = round(as.numeric(prop), 2),
row.names = NULL
)
}
tabla_frecuencias_por_estado <- function(data, var, nombre_variable = var) {
tab <- table(data[[var]], data$ESTADO_ACADEMICO, useNA = "ifany")
prop_fila <- prop.table(tab, margin = 1) * 100
salida <- data.frame(
Variable = nombre_variable,
Categoria = rownames(tab),
Desertor = as.numeric(tab[, "Desertor"]),
No_Desertor = as.numeric(tab[, "No_Desertor"]),
Total = rowSums(tab),
Porc_Desertor = round(as.numeric(prop_fila[, "Desertor"]), 2),
Porc_No_Desertor = round(as.numeric(prop_fila[, "No_Desertor"]), 2),
row.names = NULL
)
salida
}
############################################################
# TABLA 4.5 COBERTURA DE LA POBLACIÓN CON INFORMACIÓN DISPONIBLE
############################################################
tabla_4_5_cobertura <- data.frame(
Poblacion_Inicial = 982,
Poblacion_Analizada = nrow(tesis),
Cobertura_Porcentaje = round((nrow(tesis) / 982) * 100, 2)
)
print(tabla_4_5_cobertura)## Poblacion_Inicial Poblacion_Analizada Cobertura_Porcentaje
## 1 982 682 69.45############################################################
# TABLA 4.7 SITUACIÓN ACADÉMICA DE LOS INGRESANTES
############################################################
tabla_estado <- table(tesis$ESTADO_ACADEMICO)
tabla_4_7_situacion_academica <- data.frame(
Estado_Academico = names(tabla_estado),
Frecuencia = as.numeric(tabla_estado),
Porcentaje = round(as.numeric(prop.table(tabla_estado)) * 100, 2),
row.names = NULL
)
print(tabla_4_7_situacion_academica)## Estado_Academico Frecuencia Porcentaje
## 1 Desertor 459 67.3
## 2 No_Desertor 223 32.7############################################################
# TABLA 4.8 VARIABLES DEMOGRÁFICAS
############################################################
vars_demograficas <- c("SEXO", "RESIDENCIA", "EDAD", "EST_CIVIL")
tabla_4_8_demograficas_cat <- do.call(
rbind,
lapply(c("SEXO", "RESIDENCIA", "EST_CIVIL"), function(v) {
tabla_frecuencias(tesis, v, v)
})
)
tabla_4_8_demograficas_num <- data.frame(
Variable = "EDAD",
Casos = sum(!is.na(tesis$EDAD)),
Minimo = round(min(tesis$EDAD, na.rm = TRUE), 2),
Q1 = round(quantile(tesis$EDAD, 0.25, na.rm = TRUE), 2),
Mediana = round(median(tesis$EDAD, na.rm = TRUE), 2),
Media = round(mean(tesis$EDAD, na.rm = TRUE), 2),
Q3 = round(quantile(tesis$EDAD, 0.75, na.rm = TRUE), 2),
Maximo = round(max(tesis$EDAD, na.rm = TRUE), 2)
)
print(tabla_4_8_demograficas_cat)## Variable Categoria Frecuencia Porcentaje
## 1 SEXO FEMENINO 395 57.92
## 2 SEXO MASCULINO 287 42.08
## 3 RESIDENCIA ASUNCIÓN 229 33.58
## 4 RESIDENCIA CENTRAL 336 49.27
## 5 RESIDENCIA RESTO DEL PAÍS 117 17.16
## 6 EST_CIVIL Soltero/a 544 79.77
## 7 EST_CIVIL Casado/a 117 17.16
## 8 EST_CIVIL Divorciado/a 14 2.05
## 9 EST_CIVIL Otro 7 1.03row.names(tabla_4_8_demograficas_num) <- NULL
print(tabla_4_8_demograficas_num)## Variable Casos Minimo Q1 Mediana Media Q3 Maximo
## 1 EDAD 682 18 20 23 27.11 32 70############################################################
# TABLA 4.9 VARIABLES INTERNAS AL MEDIO EDUCATIVO
############################################################
vars_internas <- c("CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO")
tabla_4_9_internas <- do.call(
rbind,
lapply(vars_internas, function(v) {
tabla_frecuencias(tesis, v, v)
})
)
print(tabla_4_9_internas)## Variable Categoria Frecuencia Porcentaje
## 1 CARRERA ESTADISTICA-PRES 182 26.69
## 2 CARRERA ESTADISTICA-SEMI 100 14.66
## 3 CARRERA MATEMATICA 150 21.99
## 4 CARRERA EDUCACION MATEMATICA-PRES 70 10.26
## 5 CARRERA EDUCACION MATEMATICA-SEMI 180 26.39
## 6 TIPO_INGRESO INGRESO 642 94.13
## 7 TIPO_INGRESO TRASLADO 0 0.00
## 8 TIPO_INGRESO ADMISION DIRECTA 40 5.87
## 9 MODALIDAD Presencial 402 58.94
## 10 MODALIDAD Semipresencial 280 41.06
## 11 RENDIMIENTO APROBADO 409 59.97
## 12 RENDIMIENTO NO APROBADO 273 40.03############################################################
# TABLA 4.10 VARIABLES EXTERNAS AL MEDIO EDUCATIVO
############################################################
vars_externas <- c("TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES", "TIPO_COL", "INGRESO", "NIVEL_SOCIO")
tabla_4_10_externas <- do.call(
rbind,
lapply(vars_externas, function(v) {
tabla_frecuencias(tesis, v, v)
})
)
print(tabla_4_10_externas)## Variable Categoria Frecuencia
## 1 TRABAJA SI 372
## 2 TRABAJA NO 310
## 3 SOLVENTAR Beca/exoneración total 45
## 4 SOLVENTAR Beca/exoneración parcial 49
## 5 SOLVENTAR Trabajo Personal 346
## 6 SOLVENTAR Ayuda Familiar 242
## 7 ESTUDIOS_PADRES hasta 13 años 157
## 8 ESTUDIOS_PADRES 14-23 años 215
## 9 ESTUDIOS_PADRES 24-29 años 164
## 10 ESTUDIOS_PADRES 30-34 años 101
## 11 ESTUDIOS_PADRES Más de 34 años 45
## 12 TIPO_COL Público 465
## 13 TIPO_COL Subvencionado 49
## 14 TIPO_COL Privado 168
## 15 INGRESO Hasta dos salarios mínimos 430
## 16 INGRESO Más de dos y hasta cinco salarios mínimos 190
## 17 INGRESO Más de cinco y hasta diez salarios mínimos 45
## 18 INGRESO Más de diez y hasta quince salarios mínimos 12
## 19 INGRESO Más de quince salarios mínimos 5
## 20 NIVEL_SOCIO BAJO 65
## 21 NIVEL_SOCIO MEDIO 617
## Porcentaje
## 1 54.55
## 2 45.45
## 3 6.60
## 4 7.18
## 5 50.73
## 6 35.48
## 7 23.02
## 8 31.52
## 9 24.05
## 10 14.81
## 11 6.60
## 12 68.18
## 13 7.18
## 14 24.63
## 15 63.05
## 16 27.86
## 17 6.60
## 18 1.76
## 19 0.73
## 20 9.53
## 21 90.47############################################################
# TABLA 4.11 RESUMEN DESCRIPTIVO DE VARIABLES NUMÉRICAS
############################################################
# En el dataset depurado la variable numérica disponible es EDAD.
# Si luego reincorporas ASIG_APROB, aquí puede agregarse del mismo modo.
tabla_4_11_numericas <- data.frame(
Variable = "EDAD",
Casos = sum(!is.na(tesis$EDAD)),
Minimo = round(min(tesis$EDAD, na.rm = TRUE), 2),
Q1 = round(quantile(tesis$EDAD, 0.25, na.rm = TRUE), 2),
Mediana = round(median(tesis$EDAD, na.rm = TRUE), 2),
Media = round(mean(tesis$EDAD, na.rm = TRUE), 2),
Q3 = round(quantile(tesis$EDAD, 0.75, na.rm = TRUE), 2),
Maximo = round(max(tesis$EDAD, na.rm = TRUE), 2),
DE = round(sd(tesis$EDAD, na.rm = TRUE), 2)
)
row.names(tabla_4_11_numericas) <- NULL
print(tabla_4_11_numericas)## Variable Casos Minimo Q1 Mediana Media Q3 Maximo DE
## 1 EDAD 682 18 20 23 27.11 32 70 9.17############################################################
# TABLAS OPCIONALES POR ESTADO ACADÉMICO
# ÚTILES SI QUIERES DESCRIPTIVAS MÁS RICAS EN CAPÍTULO 4
############################################################
tabla_4_8_demograficas_estado <- do.call(
rbind,
lapply(c("SEXO", "RESIDENCIA", "EST_CIVIL"), function(v) {
tabla_frecuencias_por_estado(tesis, v, v)
})
)
tabla_4_9_internas_estado <- do.call(
rbind,
lapply(vars_internas, function(v) {
tabla_frecuencias_por_estado(tesis, v, v)
})
)
tabla_4_10_externas_estado <- do.call(
rbind,
lapply(vars_externas, function(v) {
tabla_frecuencias_por_estado(tesis, v, v)
})
)
############################################################
# EXPORTACIÓN DE ESTAS TABLAS DESCRIPTIVAS
############################################################
write.csv(tabla_4_5_cobertura,
"salidas_tesis/tabla_4_5_cobertura.csv",
row.names = FALSE)
write.csv(tabla_4_7_situacion_academica,
"salidas_tesis/tabla_4_7_situacion_academica.csv",
row.names = FALSE)
write.csv(tabla_4_8_demograficas_cat,
"salidas_tesis/tabla_4_8_demograficas_categoricas.csv",
row.names = FALSE)
write.csv(tabla_4_8_demograficas_num,
"salidas_tesis/tabla_4_8_demograficas_numericas.csv",
row.names = FALSE)
write.csv(tabla_4_9_internas,
"salidas_tesis/tabla_4_9_internas.csv",
row.names = FALSE)
write.csv(tabla_4_10_externas,
"salidas_tesis/tabla_4_10_externas.csv",
row.names = FALSE)
write.csv(tabla_4_11_numericas,
"salidas_tesis/tabla_4_11_numericas.csv",
row.names = FALSE)
############################################################
# 7. ANÁLISIS BIVARIANTE
############################################################
clasificar_significancia <- function(p) {
ifelse(p < 0.01, "*** (p < 0.01)",
ifelse(p < 0.05, "** (p < 0.05)",
ifelse(p < 0.10, "* (p < 0.10)", "No significativo")))
}
# 7.1 Chi-cuadrado / Fisher para variables categóricas
covariables_cat <- c(
"RESIDENCIA", "CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO",
"TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES", "SEXO", "TIPO_COL",
"INGRESO", "EST_CIVIL", "NIVEL_SOCIO"
)
resultados_chi <- data.frame()
for (var in covariables_cat) {
tabla <- table(tesis[[var]], tesis$ESTADO_ACADEMICO)
chi_aux <- suppressWarnings(chisq.test(tabla, correct = FALSE))
if (any(chi_aux$expected < 5)) {
prueba <- fisher.test(tabla)
estadistico <- NA
metodo <- "Fisher"
} else {
prueba <- chi_aux
estadistico <- unname(prueba$statistic)
metodo <- "Chi-cuadrado"
}
resultados_chi <- rbind(
resultados_chi,
data.frame(
Variable = var,
Metodo = metodo,
Estadistico = estadistico,
P_Valor = prueba$p.value,
Significancia = clasificar_significancia(prueba$p.value)
)
)
}
resultados_chi <- resultados_chi[order(resultados_chi$P_Valor), ]
row.names(resultados_chi) <- NULL
print(resultados_chi)## Variable Metodo Estadistico
## 1 RENDIMIENTO Chi-cuadrado 144.9473877
## 2 CARRERA Chi-cuadrado 73.1012024
## 3 MODALIDAD Chi-cuadrado 36.7897891
## 4 TRABAJA Chi-cuadrado 26.8970134
## 5 SOLVENTAR Chi-cuadrado 24.5391103
## 6 NIVEL_SOCIO Chi-cuadrado 2.5516367
## 7 EST_CIVIL Fisher NA
## 8 INGRESO Fisher NA
## 9 TIPO_INGRESO Fisher NA
## 10 RESIDENCIA Chi-cuadrado 2.4853373
## 11 ESTUDIOS_PADRES Chi-cuadrado 3.9439428
## 12 SEXO Chi-cuadrado 0.4037587
## 13 TIPO_COL Chi-cuadrado 0.9341439
## P_Valor Significancia
## 1 0.000000000000000000000000000000002205278 *** (p < 0.01)
## 2 0.000000000000005022182473894643603418153 *** (p < 0.01)
## 3 0.000000001315775866091108370413247097286 *** (p < 0.01)
## 4 0.000000214590233249842945095417490186662 *** (p < 0.01)
## 5 0.000019274982894101492435397215974290930 *** (p < 0.01)
## 6 0.110180066209064966842312571770889917389 No significativo
## 7 0.113516025141927309749512176040298072621 No significativo
## 8 0.125957268782453585265201922993583139032 No significativo
## 9 0.223085248122429330441818251529184635729 No significativo
## 10 0.288612981066232343607680377317592501640 No significativo
## 11 0.413645526632915150333502651847084052861 No significativo
## 12 0.525154472564691721103713462071027606726 No significativo
## 13 0.626834985394710275663499032816616818309 No significativo# 7.2 Boxplot de EDAD
graf_edad <- ggplot(tesis, aes(x = ESTADO_ACADEMICO, y = EDAD)) +
geom_boxplot(fill = "skyblue", color = "darkblue") +
labs(
title = "Edad según estado académico",
x = "Estado académico",
y = "Edad"
) +
theme_minimal()
print(graf_edad)
hacer_prueba_asociacion <- function(data, var1, var2){
datos_tmp <- data[, c(var1, var2)]
datos_tmp <- na.omit(datos_tmp)
tab <- table(datos_tmp[[var1]], datos_tmp[[var2]])
if (nrow(tab) < 2 || ncol(tab) < 2) {
return(data.frame(
Variable1 = var1,
Variable2 = var2,
Metodo = NA_character_,
Estadistico = NA_real_,
P_Valor = NA_real_,
Significancia = NA_character_,
Celdas_esperadas_menor_5 = NA_integer_
))
}
chi <- suppressWarnings(chisq.test(tab, correct = FALSE))
n_celdas_menor_5 <- sum(chi$expected < 5)
if (n_celdas_menor_5 > 0) {
prueba <- tryCatch(
fisher.test(tab, workspace = 2e8),
error = function(e) {
fisher.test(tab, simulate.p.value = TRUE, B = 10000)
}
)
metodo <- ifelse(isTRUE(prueba$simulate.p.value),
"Fisher simulado",
"Fisher exacto")
estadistico <- NA
pvalor <- prueba$p.value
} else {
prueba <- chi
metodo <- "Chi-cuadrado"
estadistico <- unname(prueba$statistic)
pvalor <- prueba$p.value
}
data.frame(
Variable1 = var1,
Variable2 = var2,
Metodo = metodo,
Estadistico = estadistico,
P_Valor = pvalor,
Significancia = clasificar_significancia(pvalor),
Celdas_esperadas_menor_5 = n_celdas_menor_5
)
}############################################################
# 8. PARTICIÓN TRAIN / TEST
############################################################
set.seed(123)
train_index <- createDataPartition(
tesis$ESTADO_ACADEMICO,
p = 0.7,
list = FALSE
)
train_data <- tesis[train_index, ]
test_data <- tesis[-train_index, ]
train_data$ESTADO_ACADEMICO <- factor(train_data$ESTADO_ACADEMICO,
levels = c("Desertor", "No_Desertor"))
test_data$ESTADO_ACADEMICO <- factor(test_data$ESTADO_ACADEMICO,
levels = c("Desertor", "No_Desertor"))
# Alineación de niveles entre train y test
for (col in names(train_data)) {
if (is.factor(train_data[[col]])) {
test_data[[col]] <- factor(test_data[[col]], levels = levels(train_data[[col]]))
}
}
cat("\n=============================\n")##
## =============================cat("PROPORCIÓN DE CLASES EN TRAIN\n")## PROPORCIÓN DE CLASES EN TRAINcat("=============================\n")## =============================print(prop.table(table(train_data$ESTADO_ACADEMICO)))##
## Desertor No_Desertor
## 0.6722338 0.3277662cat("\n=============================\n")##
## =============================cat("PROPORCIÓN DE CLASES EN TEST\n")## PROPORCIÓN DE CLASES EN TESTcat("=============================\n")## =============================print(prop.table(table(test_data$ESTADO_ACADEMICO)))##
## Desertor No_Desertor
## 0.6748768 0.3251232############################################################
# 9. REGRESIÓN LOGÍSTICA BIVARIADA EN TRAIN
############################################################
# Para las variables numéricas se ajusta regresión logística simple.
# En la base depurada, la variable numérica considerada es EDAD.
variables_biv <- c("EDAD")
resultados_bivariados <- data.frame()
variables_omitidas <- data.frame()
for (var in variables_biv) {
base_modelo <- train_data[, c("ESTADO_ACADEMICO", var)]
base_modelo <- na.omit(base_modelo)
if (length(unique(base_modelo$ESTADO_ACADEMICO)) < 2) {
variables_omitidas <- rbind(
variables_omitidas,
data.frame(Variable = var, Motivo = "La respuesta tiene menos de 2 niveles")
)
next
}
if (is.numeric(base_modelo[[var]]) && length(unique(base_modelo[[var]])) < 2) {
variables_omitidas <- rbind(
variables_omitidas,
data.frame(Variable = var, Motivo = "Variable numérica sin variación")
)
next
}
formula_biv <- as.formula(paste("ESTADO_ACADEMICO ~", var))
modelo <- glm(formula_biv, family = binomial, data = base_modelo)
resumen <- summary(modelo)
coefs <- resumen$coefficients[-1, , drop = FALSE]
for (i in 1:nrow(coefs)) {
beta <- coefs[i, 1]
se <- coefs[i, 2]
z <- coefs[i, 3]
pvalor <- coefs[i, 4]
or <- exp(beta)
li <- exp(beta - 1.96 * se)
ls <- exp(beta + 1.96 * se)
resultados_bivariados <- rbind(
resultados_bivariados,
data.frame(
Variable = var,
Nivel = rownames(coefs)[i],
Coeficiente = round(beta, 4),
OR = round(or, 4),
IC95_LI = round(li, 4),
IC95_LS = round(ls, 4),
Estadistico_z = round(z, 4),
P_Valor = round(pvalor, 4),
Significancia = clasificar_significancia(pvalor)
)
)
}
}
print(resultados_bivariados)## Variable Nivel Coeficiente OR IC95_LI IC95_LS Estadistico_z P_Valor
## 1 EDAD EDAD -0.0292 0.9713 0.9492 0.9939 -2.4841 0.013
## Significancia
## 1 ** (p < 0.05)print(variables_omitidas)## data frame with 0 columns and 0 rows############################################################
# MATRIZ GRÁFICA DE ASOCIACIONES BIVARIADAS
# 14 VARIABLES EXPLICATIVAS
############################################################
# Paquetes
paquetes <- c("DescTools", "corrplot")
instalar <- paquetes[!paquetes %in% installed.packages()[, "Package"]]
if (length(instalar) > 0) install.packages(instalar)
library(DescTools)##
## Adjuntando el paquete: 'DescTools'## The following object is masked from 'package:car':
##
## Recode## The following objects are masked from 'package:caret':
##
## MAE, RMSElibrary(corrplot)## corrplot 0.95 loaded############################################################
# 1. DEFINIR VARIABLES EXPLICATIVAS
############################################################
variables_exp <- c(
"RESIDENCIA", "CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO",
"TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES", "SEXO", "TIPO_COL",
"INGRESO", "EST_CIVIL", "NIVEL_SOCIO", "EDAD"
)
############################################################
# 2. SELECCIONAR BASE DE DATOS
# Usa train_data si ya existe; si no, usa tesis
############################################################
if (exists("train_data")) {
base_asoc <- train_data
} else if (exists("tesis")) {
base_asoc <- tesis
} else {
stop("No se encontró ni 'train_data' ni 'tesis' en el entorno.")
}
faltantes <- setdiff(variables_exp, names(base_asoc))
if (length(faltantes) > 0) {
stop("Faltan estas variables en la base: ", paste(faltantes, collapse = ", "))
}
datos_asoc <- base_asoc[, variables_exp, drop = FALSE]
# Convertir character a factor
datos_asoc[] <- lapply(datos_asoc, function(x) {
if (is.character(x)) factor(x) else x
})
############################################################
# 3. FUNCIONES AUXILIARES
############################################################
# Eta (η): intensidad de asociación entre una variable numérica
# y una categórica. Devuelve valor entre 0 y 1.
eta_coef <- function(y, grupo) {
completos <- complete.cases(y, grupo)
y <- y[completos]
grupo <- droplevels(as.factor(grupo[completos]))
if (length(y) < 3 || nlevels(grupo) < 2) {
return(list(stat = NA_real_, p = NA_real_))
}
media_total <- mean(y)
medias_grupo <- tapply(y, grupo, mean)
n_grupo <- tapply(y, grupo, length)
ss_between <- sum(n_grupo * (medias_grupo - media_total)^2)
ss_total <- sum((y - media_total)^2)
eta <- ifelse(ss_total == 0, NA_real_, sqrt(ss_between / ss_total))
p_val <- tryCatch({
summary(aov(y ~ grupo))[[1]][["Pr(>F)"]][1]
}, error = function(e) NA_real_)
list(stat = eta, p = p_val)
}
# Asociación entre dos variables
calc_asociacion <- function(x, y) {
completos <- complete.cases(x, y)
x <- x[completos]
y <- y[completos]
if (length(x) < 3) {
return(list(stat = NA_real_, p = NA_real_))
}
x_num <- is.numeric(x)
y_num <- is.numeric(y)
# Numérica vs numérica -> Spearman
if (x_num && y_num) {
prueba <- suppressWarnings(cor.test(x, y, method = "spearman"))
return(list(stat = abs(unname(prueba$estimate)), p = prueba$p.value))
}
# Categórica vs categórica -> Cramér's V
if (!x_num && !y_num) {
x <- droplevels(as.factor(x))
y <- droplevels(as.factor(y))
if (nlevels(x) < 2 || nlevels(y) < 2) {
return(list(stat = NA_real_, p = NA_real_))
}
tabla <- table(x, y)
if (nrow(tabla) < 2 || ncol(tabla) < 2) {
return(list(stat = NA_real_, p = NA_real_))
}
p_val <- suppressWarnings(chisq.test(tabla, correct = FALSE)$p.value)
v <- suppressWarnings(CramerV(tabla, bias.correct = TRUE))
return(list(stat = as.numeric(v), p = p_val))
}
# Numérica vs categórica -> Eta
if (x_num && !y_num) return(eta_coef(x, y))
if (!x_num && y_num) return(eta_coef(y, x))
}
############################################################
# 4. CONSTRUIR MATRIZ DE ASOCIACIONES Y MATRIZ DE P-VALORES
############################################################
n <- length(variables_exp)
matriz_asoc <- matrix(NA_real_, nrow = n, ncol = n,
dimnames = list(variables_exp, variables_exp))
matriz_p <- matrix(NA_real_, nrow = n, ncol = n,
dimnames = list(variables_exp, variables_exp))
for (i in seq_len(n)) {
for (j in seq_len(n)) {
if (i == j) {
matriz_asoc[i, j] <- 1
matriz_p[i, j] <- 0
} else if (i < j) {
res <- calc_asociacion(datos_asoc[[i]], datos_asoc[[j]])
matriz_asoc[i, j] <- res$stat
matriz_asoc[j, i] <- res$stat
matriz_p[i, j] <- res$p
matriz_p[j, i] <- res$p
}
}
}
############################################################
# 5. REVISIÓN DE VARIABLES PROBLEMÁTICAS
############################################################
variables_constantes <- names(datos_asoc)[sapply(datos_asoc, function(x) {
length(unique(na.omit(x))) <= 1
})]
if (length(variables_constantes) > 0) {
message("Variables con un solo nivel o sin variación: ",
paste(variables_constantes, collapse = ", "))
}
############################################################
# 6. MATRIZ GRÁFICA
# Se muestran solo asociaciones significativas (p < 0.05)
############################################################
png("matriz_asociaciones_bivariadas_14_variables.png",
width = 1800, height = 1600, res = 220)
corrplot(
matriz_asoc,
method = "color",
type = "upper",
order = "hclust",
diag = FALSE,
is.corr = FALSE,
cl.lim = c(0, 1),
tl.col = "black",
tl.srt = 45,
addCoef.col = "black",
number.cex = 0.65,
p.mat = matriz_p,
sig.level = 0.05,
insig = "blank",
mar = c(0, 0, 2, 0),
title = "Matriz gráfica de asociaciones bivariadas entre variables explicativas"
)## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "cl.lim" es un parámetro gráfico inválido## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "cl.lim" es un parámetro gráfico inválido## Warning in title(title, ...): "cl.lim" es un parámetro gráfico inválidodev.off()## png
## 2############################################################
# 7. EXPORTAR MATRICES
############################################################
write.csv(matriz_asoc,
"matriz_asociaciones_bivariadas_14_variables.csv",
row.names = TRUE)
write.csv(matriz_p,
"matriz_pvalores_asociaciones_14_variables.csv",
row.names = TRUE)
############################################################
# 8. MOSTRAR EN CONSOLA
############################################################
print(round(matriz_asoc, 3))## RESIDENCIA CARRERA TIPO_INGRESO MODALIDAD RENDIMIENTO TRABAJA
## RESIDENCIA 1.000 0.135 0.073 0.168 0.030 0.074
## CARRERA 0.135 1.000 0.372 1.000 0.310 0.545
## TIPO_INGRESO 0.073 0.372 1.000 0.188 0.055 0.158
## MODALIDAD 0.168 1.000 0.188 1.000 0.208 0.535
## RENDIMIENTO 0.030 0.310 0.055 0.208 1.000 0.225
## TRABAJA 0.074 0.545 0.158 0.535 0.225 1.000
## SOLVENTAR 0.078 0.293 0.122 0.466 0.168 0.755
## ESTUDIOS_PADRES 0.183 0.099 0.086 0.106 0.104 0.085
## SEXO 0.066 0.200 0.070 0.098 0.137 0.120
## TIPO_COL 0.102 0.050 0.055 0.023 0.073 0.060
## INGRESO 0.167 0.186 0.054 0.333 0.141 0.305
## EST_CIVIL 0.054 0.270 0.106 0.438 0.061 0.336
## NIVEL_SOCIO 0.090 0.162 0.062 0.071 0.021 0.040
## EDAD 0.154 0.552 0.144 0.543 0.139 0.481
## SOLVENTAR ESTUDIOS_PADRES SEXO TIPO_COL INGRESO EST_CIVIL
## RESIDENCIA 0.078 0.183 0.066 0.102 0.167 0.054
## CARRERA 0.293 0.099 0.200 0.050 0.186 0.270
## TIPO_INGRESO 0.122 0.086 0.070 0.055 0.054 0.106
## MODALIDAD 0.466 0.106 0.098 0.023 0.333 0.438
## RENDIMIENTO 0.168 0.104 0.137 0.073 0.141 0.061
## TRABAJA 0.755 0.085 0.120 0.060 0.305 0.336
## SOLVENTAR 1.000 0.096 0.177 0.087 0.197 0.173
## ESTUDIOS_PADRES 0.096 1.000 0.146 0.201 0.233 0.099
## SEXO 0.177 0.146 1.000 0.015 0.215 0.082
## TIPO_COL 0.087 0.201 0.015 1.000 0.197 0.073
## INGRESO 0.197 0.233 0.215 0.197 1.000 0.151
## EST_CIVIL 0.173 0.099 0.082 0.073 0.151 1.000
## NIVEL_SOCIO 0.225 0.133 0.031 0.060 0.177 0.033
## EDAD 0.433 0.089 0.060 0.033 0.292 0.524
## NIVEL_SOCIO EDAD
## RESIDENCIA 0.090 0.154
## CARRERA 0.162 0.552
## TIPO_INGRESO 0.062 0.144
## MODALIDAD 0.071 0.543
## RENDIMIENTO 0.021 0.139
## TRABAJA 0.040 0.481
## SOLVENTAR 0.225 0.433
## ESTUDIOS_PADRES 0.133 0.089
## SEXO 0.031 0.060
## TIPO_COL 0.060 0.033
## INGRESO 0.177 0.292
## EST_CIVIL 0.033 0.524
## NIVEL_SOCIO 1.000 0.057
## EDAD 0.057 1.000print(round(matriz_p, 4))## RESIDENCIA CARRERA TIPO_INGRESO MODALIDAD RENDIMIENTO TRABAJA
## RESIDENCIA 0.0000 0.0267 0.2774 0.0012 0.8032 0.2714
## CARRERA 0.0267 0.0000 0.0000 0.0000 0.0000 0.0000
## TIPO_INGRESO 0.2774 0.0000 0.0000 0.0000 0.2290 0.0005
## MODALIDAD 0.0012 0.0000 0.0000 0.0000 0.0000 0.0000
## RENDIMIENTO 0.8032 0.0000 0.2290 0.0000 0.0000 0.0000
## TRABAJA 0.2714 0.0000 0.0005 0.0000 0.0000 0.0000
## SOLVENTAR 0.4407 0.0000 0.0672 0.0000 0.0037 0.0000
## ESTUDIOS_PADRES 0.0001 0.2721 0.4673 0.2539 0.2730 0.4855
## SEXO 0.3571 0.0007 0.1262 0.0314 0.0028 0.0085
## TIPO_COL 0.0417 0.9674 0.4840 0.8772 0.2761 0.4209
## INGRESO 0.0008 0.0000 0.8470 0.0000 0.0486 0.0000
## EST_CIVIL 0.8376 0.0000 0.1472 0.0000 0.6138 0.0000
## NIVEL_SOCIO 0.1419 0.0138 0.1720 0.1189 0.6427 0.3831
## EDAD 0.0034 0.0000 0.0016 0.0000 0.0022 0.0000
## SOLVENTAR ESTUDIOS_PADRES SEXO TIPO_COL INGRESO EST_CIVIL
## RESIDENCIA 0.4407 0.0001 0.3571 0.0417 0.0008 0.8376
## CARRERA 0.0000 0.2721 0.0007 0.9674 0.0000 0.0000
## TIPO_INGRESO 0.0672 0.4673 0.1262 0.4840 0.8470 0.1472
## MODALIDAD 0.0000 0.2539 0.0314 0.8772 0.0000 0.0000
## RENDIMIENTO 0.0037 0.2730 0.0028 0.2761 0.0486 0.6138
## TRABAJA 0.0000 0.4855 0.0085 0.4209 0.0000 0.0000
## SOLVENTAR 0.0000 0.3432 0.0018 0.2949 0.0000 0.0000
## ESTUDIOS_PADRES 0.3432 0.0000 0.0371 0.0000 0.0000 0.2968
## SEXO 0.0018 0.0371 0.0000 0.9474 0.0002 0.3573
## TIPO_COL 0.2949 0.0000 0.9474 0.0000 0.0000 0.5359
## INGRESO 0.0000 0.0000 0.0002 0.0000 0.0000 0.0011
## EST_CIVIL 0.0000 0.2968 0.3573 0.5359 0.0011 0.0000
## NIVEL_SOCIO 0.0000 0.0747 0.4961 0.4280 0.0048 0.9170
## EDAD 0.0000 0.4310 0.1915 0.7684 0.0000 0.0000
## NIVEL_SOCIO EDAD
## RESIDENCIA 0.1419 0.0034
## CARRERA 0.0138 0.0000
## TIPO_INGRESO 0.1720 0.0016
## MODALIDAD 0.1189 0.0000
## RENDIMIENTO 0.6427 0.0022
## TRABAJA 0.3831 0.0000
## SOLVENTAR 0.0000 0.0000
## ESTUDIOS_PADRES 0.0747 0.4310
## SEXO 0.4961 0.1915
## TIPO_COL 0.4280 0.7684
## INGRESO 0.0048 0.0000
## EST_CIVIL 0.9170 0.0000
## NIVEL_SOCIO 0.0000 0.2170
## EDAD 0.2170 0.0000############################################################
# 9. DIRECTORIO DE SALIDA
############################################################
dir_salida <- "salida_tesis"
if (!dir.exists(dir_salida)) {
dir.create(dir_salida, recursive = TRUE)
}
############################################################
# 9.1 PARTICIÓN TRAIN / TEST
############################################################
library(caret)
set.seed(123)
# Variables completas para modelado
vars_modelo_total <- c(
"ESTADO_ACADEMICO",
"RESIDENCIA", "CARRERA", "TIPO_INGRESO",
"MODALIDAD", "RENDIMIENTO", "TRABAJA",
"SOLVENTAR", "ESTUDIOS_PADRES", "SEXO",
"TIPO_COL", "INGRESO", "EST_CIVIL",
"NIVEL_SOCIO", "EDAD"
)
# Base limpia para modelado
datos_modelo <- tesis[, vars_modelo_total]
# Eliminar NA
datos_modelo <- datos_modelo[complete.cases(datos_modelo), ]
# Partición estratificada
indice_train <- createDataPartition(
datos_modelo$ESTADO_ACADEMICO,
p = 0.70,
list = FALSE
)
train_data <- datos_modelo[indice_train, ]
test_data <- datos_modelo[-indice_train, ]
cat("\n=============================\n")##
## =============================cat("PARTICIÓN TRAIN / TEST\n")## PARTICIÓN TRAIN / TESTcat("=============================\n")## =============================cat("Train:", nrow(train_data), "observaciones\n")## Train: 479 observacionescat("Test :", nrow(test_data), "observaciones\n")## Test : 203 observaciones############################################################
# 10. MODELADO PREDICTIVO - CONFIGURACIÓN COMÚN
############################################################
# Conjunto común de 14 variables explicativas
variables_modelado <- c(
"RESIDENCIA", "CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO",
"TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES", "SEXO", "TIPO_COL",
"INGRESO", "EST_CIVIL", "NIVEL_SOCIO", "EDAD"
)
# Fórmula común para TODOS los modelos comparados
form_general <- as.formula(
paste("ESTADO_ACADEMICO ~", paste(variables_modelado, collapse = " + "))
)
# Fórmula saturada para el modelo logístico final por stepwise
form_stepwise <- form_general
############################################################
# 10.1 BASES CONSISTENTES PARA MODELADO
############################################################
vars_modelo_total <- c("ESTADO_ACADEMICO", variables_modelado)
train_model <- train_data[, vars_modelo_total]
test_model <- test_data[, vars_modelo_total]
train_model <- train_model[complete.cases(train_model), ]
test_model <- test_model[complete.cases(test_model), ]
for (col in names(train_model)) {
if (is.factor(train_model[[col]])) {
train_model[[col]] <- droplevels(train_model[[col]])
}
}
for (col in names(train_model)) {
if (is.factor(train_model[[col]])) {
test_model[[col]] <- factor(test_model[[col]], levels = levels(train_model[[col]]))
}
}
cat("\n=============================\n")##
## =============================cat("BASE FINAL DE TRAIN PARA MODELADO\n")## BASE FINAL DE TRAIN PARA MODELADOcat("=============================\n")## =============================cat("Filas:", nrow(train_model), "\n")## Filas: 479cat("Columnas:", ncol(train_model), "\n")## Columnas: 15cat("\n=============================\n")##
## =============================cat("BASE FINAL DE TEST PARA MODELADO\n")## BASE FINAL DE TEST PARA MODELADOcat("=============================\n")## =============================cat("Filas:", nrow(test_model), "\n")## Filas: 203cat("Columnas:", ncol(test_model), "\n")## Columnas: 15############################################################
# 10.2 VALIDACIÓN CRUZADA COMPARTIDA
############################################################
set.seed(123)
folds_cv <- createMultiFolds(
y = train_model$ESTADO_ACADEMICO,
k = 10,
times = 3
)
ctrl <- trainControl(
method = "repeatedcv",
number = 10,
repeats = 3,
index = folds_cv,
classProbs = TRUE,
summaryFunction = twoClassSummary,
savePredictions = "final"
)
############################################################
# 11. ENTRENAMIENTO DE LOS 7 MODELOS
############################################################
# 11.1 Regresión logística
set.seed(123)
modelo_logit_cv <- train(
form_general,
data = train_model,
method = "glm",
family = binomial,
metric = "ROC",
trControl = ctrl
)
# 11.2 Árbol de decisión
set.seed(123)
modelo_tree_cv <- train(
form_general,
data = train_model,
method = "rpart",
metric = "ROC",
tuneLength = 10,
trControl = ctrl
)
# 11.3 Random Forest
set.seed(123)
modelo_rf_cv <- train(
form_general,
data = train_model,
method = "rf",
metric = "ROC",
tuneLength = 5,
ntree = 500,
trControl = ctrl
)
# 11.4 SVM radial
set.seed(123)
modelo_svm_cv <- train(
form_general,
data = train_model,
method = "svmRadial",
metric = "ROC",
preProcess = c("center", "scale"),
tuneLength = 10,
trControl = ctrl
)
# 11.5 KNN
set.seed(123)
modelo_knn_cv <- train(
form_general,
data = train_model,
method = "knn",
metric = "ROC",
preProcess = c("center", "scale"),
tuneLength = 10,
trControl = ctrl
)
# 11.6 Red neuronal
grid_nn <- expand.grid(
size = c(3, 5, 7),
decay = c(0.01, 0.10, 0.50)
)
set.seed(123)
modelo_nn_cv <- train(
form_general,
data = train_model,
method = "nnet",
metric = "ROC",
preProcess = c("center", "scale"),
tuneGrid = grid_nn,
trControl = ctrl,
trace = FALSE,
MaxNWts = 5000
)
# 11.7 Naive Bayes
grid_nb <- expand.grid(
laplace = c(0, 1),
usekernel = c(TRUE, FALSE),
adjust = c(1)
)
set.seed(123)
modelo_nb_cv <- train(
form_general,
data = train_model,
method = "naive_bayes",
metric = "ROC",
tuneGrid = grid_nb,
trControl = ctrl
)
############################################################
# 12. COMPARACIÓN EN VALIDACIÓN CRUZADA
############################################################
resultados_cv <- resamples(list(
Logit = modelo_logit_cv,
Tree = modelo_tree_cv,
RF = modelo_rf_cv,
SVM = modelo_svm_cv,
KNN = modelo_knn_cv,
NN = modelo_nn_cv,
NB = modelo_nb_cv
))
print(summary(resultados_cv))##
## Call:
## summary.resamples(object = resultados_cv)
##
## Models: Logit, Tree, RF, SVM, KNN, NN, NB
## Number of resamples: 30
##
## ROC
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## Logit 0.6541667 0.7437500 0.7836174 0.7803046 0.8199219 0.8996212 0
## Tree 0.6845703 0.7714962 0.7973633 0.7940493 0.8305827 0.8750000 0
## RF 0.6458333 0.7341856 0.7910156 0.7854376 0.8286458 0.9091797 0
## SVM 0.6416667 0.7183268 0.7613636 0.7602734 0.8041992 0.8867188 0
## KNN 0.6318359 0.6966072 0.7434186 0.7483594 0.7955211 0.9062500 0
## NN 0.6208333 0.7228634 0.7753906 0.7677115 0.8043176 0.8916667 0
## NB 0.6041667 0.7292480 0.7802734 0.7714621 0.8144383 0.9104167 0
##
## Sens
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## Logit 0.68750 0.7500000 0.8125000 0.8071023 0.8473011 0.9696970 0
## Tree 0.65625 0.7500000 0.7812500 0.7886995 0.8167614 0.9375000 0
## RF 0.90625 0.9375000 0.9687500 0.9616477 0.9924242 1.0000000 0
## SVM 0.68750 0.8437500 0.8787879 0.8777462 0.9375000 1.0000000 0
## KNN 0.81250 0.9062500 0.9375000 0.9304924 0.9614110 1.0000000 0
## NN 0.65625 0.7500000 0.8125000 0.7967487 0.8437500 0.9090909 0
## NB 0.37500 0.6328125 0.7187500 0.6799874 0.7812500 0.8181818 0
##
## Spec
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## Logit 0.3333333 0.5000000 0.5625000 0.5548611 0.62500 0.7500 0
## Tree 0.3125000 0.4666667 0.6125000 0.5850000 0.68750 0.8750 0
## RF 0.0625000 0.1468750 0.1875000 0.2031944 0.25000 0.3750 0
## SVM 0.0000000 0.2000000 0.3229167 0.3322222 0.43750 0.6875 0
## KNN 0.0000000 0.1468750 0.2000000 0.2259722 0.31250 0.4375 0
## NN 0.2500000 0.3750000 0.5000000 0.4906944 0.61875 0.7500 0
## NB 0.4375000 0.6250000 0.7500000 0.7298611 0.81250 1.0000 0# Gráficos en objeto
grafico_cv_roc <- bwplot(
resultados_cv,
metric = "ROC",
main = "Comparación del ROC en validación cruzada"
)
grafico_cv_sens <- bwplot(
resultados_cv,
metric = "Sens",
main = "Comparación de la sensibilidad en validación cruzada"
)
grafico_cv_spec <- bwplot(
resultados_cv,
metric = "Spec",
main = "Comparación de la especificidad en validación cruzada"
)
# Mostrar en pantalla
print(grafico_cv_roc)
print(grafico_cv_sens)
print(grafico_cv_spec)
############################################################
# 12.1 GUARDAR GRÁFICOS EN LA CARPETA salida_tesis
############################################################
# ROC
png(
filename = file.path(dir_salida, "grafico_cv_roc.png"),
width = 2000, height = 1400, res = 200
)
print(grafico_cv_roc)
dev.off()## png
## 2pdf(
file = file.path(dir_salida, "grafico_cv_roc.pdf"),
width = 10, height = 7
)
print(grafico_cv_roc)
dev.off()## png
## 2# Sensibilidad
png(
filename = file.path(dir_salida, "grafico_cv_sensibilidad.png"),
width = 2000, height = 1400, res = 200
)
print(grafico_cv_sens)
dev.off()## png
## 2pdf(
file = file.path(dir_salida, "grafico_cv_sensibilidad.pdf"),
width = 10, height = 7
)
print(grafico_cv_sens)
dev.off()## png
## 2# Especificidad
png(
filename = file.path(dir_salida, "grafico_cv_especificidad.png"),
width = 2000, height = 1400, res = 200
)
print(grafico_cv_spec)
dev.off()## png
## 2pdf(
file = file.path(dir_salida, "grafico_cv_especificidad.pdf"),
width = 10, height = 7
)
print(grafico_cv_spec)
dev.off()## png
## 2cat("\n=============================\n")##
## =============================cat("GRÁFICOS GUARDADOS EN:\n")## GRÁFICOS GUARDADOS EN:cat(normalizePath(dir_salida), "\n")## C:\Users\DELL\Desktop\Tesis-Estadistica\salida_tesiscat("=============================\n")## =============================############################################################
# 13. EVALUACIÓN EN CONJUNTO DE PRUEBA
############################################################
evaluar_modelo_caret <- function(modelo, nombre_modelo, test_data, variable_respuesta = "ESTADO_ACADEMICO") {
pred_clase <- predict(modelo, newdata = test_data, type = "raw")
pred_clase <- factor(pred_clase, levels = c("Desertor", "No_Desertor"))
pred_prob <- predict(modelo, newdata = test_data, type = "prob")[, "Desertor"]
real <- factor(test_data[[variable_respuesta]], levels = c("Desertor", "No_Desertor"))
mc <- confusionMatrix(pred_clase, real, positive = "Desertor")
roc_obj <- roc(
response = real,
predictor = pred_prob,
levels = c("No_Desertor", "Desertor"),
direction = "<"
)
auc_val <- as.numeric(auc(roc_obj))
resumen <- data.frame(
Modelo = nombre_modelo,
Accuracy = as.numeric(mc$overall["Accuracy"]),
Kappa = as.numeric(mc$overall["Kappa"]),
Sensitivity = as.numeric(mc$byClass["Sensitivity"]),
Specificity = as.numeric(mc$byClass["Specificity"]),
Precision = as.numeric(mc$byClass["Pos Pred Value"]),
F1 = ifelse(
(as.numeric(mc$byClass["Pos Pred Value"]) + as.numeric(mc$byClass["Sensitivity"])) == 0,
NA,
2 * as.numeric(mc$byClass["Pos Pred Value"]) * as.numeric(mc$byClass["Sensitivity"]) /
(as.numeric(mc$byClass["Pos Pred Value"]) + as.numeric(mc$byClass["Sensitivity"]))
),
AUC = auc_val,
row.names = NULL
)
list(
nombre = nombre_modelo,
confusion = mc,
roc = roc_obj,
resumen = resumen
)
}
modelos_finales <- list(
Logit = modelo_logit_cv,
Tree = modelo_tree_cv,
RF = modelo_rf_cv,
SVM = modelo_svm_cv,
KNN = modelo_knn_cv,
NN = modelo_nn_cv,
NB = modelo_nb_cv
)
resultados_test_lista <- lapply(names(modelos_finales), function(nombre) {
evaluar_modelo_caret(modelos_finales[[nombre]], nombre, test_model)
})
names(resultados_test_lista) <- names(modelos_finales)
tabla_metricas_test <- do.call(
rbind,
lapply(resultados_test_lista, function(x) x$resumen)
)
tabla_metricas_test <- tabla_metricas_test %>%
mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
as.data.frame()
print(tabla_metricas_test)## Modelo Accuracy Kappa Sensitivity Specificity Precision F1 AUC
## Logit Logit 0.7192 0.3627 0.7883 0.5758 0.7941 0.7912 0.7796
## Tree Tree 0.7340 0.4167 0.7664 0.6667 0.8268 0.7955 0.7853
## RF RF 0.7192 0.2154 0.9635 0.2121 0.7174 0.8224 0.7920
## SVM SVM 0.7192 0.2514 0.9270 0.2879 0.7299 0.8167 0.7773
## KNN KNN 0.6798 0.1541 0.8905 0.2424 0.7093 0.7896 0.7131
## NN NN 0.6946 0.2872 0.7956 0.4848 0.7622 0.7786 0.7460
## NB NB 0.7291 0.4169 0.7445 0.6970 0.8361 0.7876 0.7890############################################################
# 14. CURVAS ROC COMPARATIVAS
############################################################
plot(resultados_test_lista$Logit$roc,
main = "Curvas ROC de los modelos en el conjunto de prueba",
lwd = 2)
plot(resultados_test_lista$Tree$roc, add = TRUE, lty = 2, lwd = 2)
plot(resultados_test_lista$RF$roc, add = TRUE, lty = 3, lwd = 2)
plot(resultados_test_lista$SVM$roc, add = TRUE, lty = 4, lwd = 2)
plot(resultados_test_lista$KNN$roc, add = TRUE, lty = 5, lwd = 2)
plot(resultados_test_lista$NN$roc, add = TRUE, lty = 6, lwd = 2)
plot(resultados_test_lista$NB$roc, add = TRUE, lty = 7, lwd = 2)
legend(
"bottomright",
legend = c(
paste0("Logit (AUC=", round(resultados_test_lista$Logit$resumen$AUC, 3), ")"),
paste0("Tree (AUC=", round(resultados_test_lista$Tree$resumen$AUC, 3), ")"),
paste0("RF (AUC=", round(resultados_test_lista$RF$resumen$AUC, 3), ")"),
paste0("SVM (AUC=", round(resultados_test_lista$SVM$resumen$AUC, 3), ")"),
paste0("KNN (AUC=", round(resultados_test_lista$KNN$resumen$AUC, 3), ")"),
paste0("NN (AUC=", round(resultados_test_lista$NN$resumen$AUC, 3), ")"),
paste0("NB (AUC=", round(resultados_test_lista$NB$resumen$AUC, 3), ")")
),
lty = 1:7,
lwd = 2,
cex = 0.8
)
############################################################
# 15. GRÁFICO COMPARATIVO DE 6 MÉTRICAS DE DESEMPEÑO
############################################################
metricas_graf <- tabla_metricas_test[, c(
"Modelo", "Accuracy", "Sensitivity", "Specificity", "Precision", "F1", "AUC"
)]
metricas_long <- tidyr::pivot_longer(
metricas_graf,
cols = c("Accuracy", "Sensitivity", "Specificity", "Precision", "F1", "AUC"),
names_to = "Metrica",
values_to = "Valor"
)
orden_modelos <- tabla_metricas_test$Modelo[order(tabla_metricas_test$AUC, decreasing = TRUE)]
metricas_long$Modelo <- factor(metricas_long$Modelo, levels = orden_modelos)
metricas_long$Metrica <- factor(
metricas_long$Metrica,
levels = c("Accuracy", "Sensitivity", "Specificity", "Precision", "F1", "AUC")
)
graf_6_metricas <- ggplot(metricas_long, aes(x = Modelo, y = Valor, fill = Metrica)) +
geom_col(position = position_dodge(width = 0.8), width = 0.7) +
geom_text(
aes(label = round(Valor, 3)),
position = position_dodge(width = 0.8),
vjust = -0.25,
size = 3
) +
labs(
title = "Comparación de 6 métricas de desempeño en el conjunto de prueba",
subtitle = "Accuracy, Sensitivity, Specificity, Precision, F1 y AUC",
x = "Modelo",
y = "Valor de la métrica",
fill = "Métrica"
) +
ylim(0, 1.08) +
theme_minimal(base_size = 13)
print(graf_6_metricas)
############################################################
# 13. EVALUACIÓN EN CONJUNTO DE PRUEBA
############################################################
library(caret)
library(pROC)
library(dplyr)
library(tidyr)
library(ggplot2)
library(forcats)
library(scales)
evaluar_modelo_caret <- function(modelo, nombre_modelo, test_data,
variable_respuesta = "ESTADO_ACADEMICO") {
# Clases predichas
pred_clase <- predict(modelo, newdata = test_data, type = "raw")
pred_clase <- factor(pred_clase, levels = c("Desertor", "No_Desertor"))
# Probabilidades
pred_prob <- predict(modelo, newdata = test_data, type = "prob")[, "Desertor"]
# Valores reales
real <- factor(test_data[[variable_respuesta]], levels = c("Desertor", "No_Desertor"))
# Matriz de confusión
mc <- confusionMatrix(pred_clase, real, positive = "Desertor")
# ROC
roc_obj <- roc(
response = real,
predictor = pred_prob,
levels = c("No_Desertor", "Desertor"),
direction = "<",
quiet = TRUE
)
auc_val <- as.numeric(auc(roc_obj))
precision_val <- as.numeric(mc$byClass["Pos Pred Value"])
sensitivity_val <- as.numeric(mc$byClass["Sensitivity"])
f1_val <- ifelse(
(precision_val + sensitivity_val) == 0,
NA,
2 * precision_val * sensitivity_val / (precision_val + sensitivity_val)
)
resumen <- data.frame(
Modelo = nombre_modelo,
Accuracy = as.numeric(mc$overall["Accuracy"]),
Kappa = as.numeric(mc$overall["Kappa"]),
Sensitivity = sensitivity_val,
Specificity = as.numeric(mc$byClass["Specificity"]),
Precision = precision_val,
F1 = f1_val,
AUC = auc_val,
row.names = NULL
)
list(
nombre = nombre_modelo,
confusion = mc,
roc = roc_obj,
resumen = resumen
)
}
############################################################
# MODELOS FINALES
############################################################
modelos_finales <- list(
Logit = modelo_logit_cv,
Tree = modelo_tree_cv,
RF = modelo_rf_cv,
SVM = modelo_svm_cv,
KNN = modelo_knn_cv,
NN = modelo_nn_cv,
NB = modelo_nb_cv
)
# CORRECCIÓN: usar test_data, no test_model
resultados_test_lista <- lapply(names(modelos_finales), function(nombre) {
evaluar_modelo_caret(modelos_finales[[nombre]], nombre, test_data)
})
names(resultados_test_lista) <- names(modelos_finales)
tabla_metricas_test <- do.call(
rbind,
lapply(resultados_test_lista, function(x) x$resumen)
)
tabla_metricas_test <- tabla_metricas_test %>%
mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
as.data.frame()
print(tabla_metricas_test)## Modelo Accuracy Kappa Sensitivity Specificity Precision F1 AUC
## Logit Logit 0.7192 0.3627 0.7883 0.5758 0.7941 0.7912 0.7796
## Tree Tree 0.7340 0.4167 0.7664 0.6667 0.8268 0.7955 0.7853
## RF RF 0.7192 0.2154 0.9635 0.2121 0.7174 0.8224 0.7920
## SVM SVM 0.7192 0.2514 0.9270 0.2879 0.7299 0.8167 0.7773
## KNN KNN 0.6798 0.1541 0.8905 0.2424 0.7093 0.7896 0.7131
## NN NN 0.6946 0.2872 0.7956 0.4848 0.7622 0.7786 0.7460
## NB NB 0.7291 0.4169 0.7445 0.6970 0.8361 0.7876 0.7890############################################################
# 14. CURVAS ROC COMPARATIVAS CON GGPLOT2
############################################################
# Construir data frame para ROC
roc_df <- bind_rows(lapply(resultados_test_lista, function(res) {
data.frame(
Modelo = res$nombre,
FPR = 1 - res$roc$specificities,
TPR = res$roc$sensitivities
)
}))
# Etiquetas con AUC
auc_labels <- tabla_metricas_test %>%
arrange(desc(AUC)) %>%
mutate(
etiqueta = paste0(Modelo, " (AUC = ", sprintf("%.3f", AUC), ")")
)
roc_df$Modelo <- factor(roc_df$Modelo, levels = auc_labels$Modelo)
graf_roc <- ggplot(roc_df, aes(x = FPR, y = TPR, color = Modelo)) +
geom_line(linewidth = 1.15, alpha = 0.95) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "grey55") +
coord_equal() +
scale_x_continuous(labels = label_number(accuracy = 0.01)) +
scale_y_continuous(labels = label_number(accuracy = 0.01)) +
labs(
title = "Curvas ROC de los modelos en el conjunto de prueba",
subtitle = "Comparación del desempeño discriminante en datos no observados",
x = "1 - Especificidad (FPR)",
y = "Sensibilidad (TPR)",
color = "Modelo"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 15),
plot.subtitle = element_text(size = 11),
legend.position = "bottom",
legend.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
) +
scale_color_discrete(labels = auc_labels$etiqueta)
print(graf_roc)
ggsave(
filename = "Curvas_ROC_modelos_test.png",
plot = graf_roc,
width = 10,
height = 7,
dpi = 320
)
############################################################
# 15. HEATMAP DE MÉTRICAS DE DESEMPEÑO
############################################################
metricas_graf <- tabla_metricas_test %>%
select(Modelo, Accuracy, Sensitivity, Specificity, Precision, F1, AUC)
# Ordenar modelos por AUC
orden_modelos <- metricas_graf %>%
arrange(desc(AUC)) %>%
pull(Modelo)
metricas_long <- metricas_graf %>%
pivot_longer(
cols = c(Accuracy, Sensitivity, Specificity, Precision, F1, AUC),
names_to = "Metrica",
values_to = "Valor"
) %>%
mutate(
Modelo = factor(Modelo, levels = orden_modelos),
Metrica = factor(Metrica,
levels = c("Accuracy", "Sensitivity", "Specificity",
"Precision", "F1", "AUC"))
)
graf_heatmap_metricas <- ggplot(metricas_long, aes(x = Metrica, y = Modelo, fill = Valor)) +
geom_tile(color = "white", linewidth = 0.7) +
geom_text(aes(label = sprintf("%.3f", Valor)), size = 3.7, fontface = "bold") +
scale_fill_gradient(
low = "#F2F2F2",
high = "#1F78B4",
limits = c(0, 1),
labels = label_number(accuracy = 0.01)
) +
labs(
title = "Matriz de métricas de desempeño en el conjunto de prueba",
subtitle = "Valores más altos indican mejor rendimiento",
x = "Métrica",
y = "Modelo",
fill = "Valor"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 15),
plot.subtitle = element_text(size = 11),
axis.text.x = element_text(angle = 20, hjust = 1, face = "bold"),
axis.text.y = element_text(face = "bold"),
legend.position = "right",
panel.grid = element_blank()
)
print(graf_heatmap_metricas)
ggsave(
filename = "Heatmap_metricas_modelos_test.png",
plot = graf_heatmap_metricas,
width = 9,
height = 5.8,
dpi = 320
)
############################################################
# 16. OPCIONAL: GRÁFICO DE BARRAS EN FACETAS
# Más limpio que 6 barras juntas por modelo
############################################################
graf_metricas_facetas <- ggplot(metricas_long, aes(x = fct_reorder(Modelo, Valor), y = Valor)) +
geom_col(width = 0.68, fill = "#2C7FB8") +
geom_text(aes(label = sprintf("%.3f", Valor)), vjust = -0.25, size = 3.3) +
facet_wrap(~ Metrica, ncol = 3) +
coord_cartesian(ylim = c(0, 1.05)) +
labs(
title = "Comparación de métricas por modelo",
subtitle = "Visualización separada por métrica",
x = "Modelo",
y = "Valor"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 15),
plot.subtitle = element_text(size = 11),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle = 30, hjust = 1)
)
print(graf_metricas_facetas)
ggsave(
filename = "Metricas_facetas_modelos_test.png",
plot = graf_metricas_facetas,
width = 11,
height = 7,
dpi = 320
)############################################################
# 16. MODELO LOGÍSTICO FINAL (STEPWISE) + VERIFICACIÓN
# DE SUPUESTOS BÁSICOS + ANÁLISIS DEL UMBRAL
############################################################
############################################################
# 16.1 VARIABLES CANDIDATAS Y MODELO INICIAL
############################################################
# Variables candidatas: las 14 covariables consideradas en el entrenamiento comparativo
variables_stepwise <- c(
"RESIDENCIA", "CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO",
"TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES", "SEXO", "TIPO_COL",
"INGRESO", "EST_CIVIL", "NIVEL_SOCIO", "EDAD"
)
form_stepwise <- as.formula(
paste("ESTADO_ACADEMICO ~", paste(variables_stepwise, collapse = " + "))
)
modelo_logit_inicial <- glm(
form_stepwise,
data = train_model,
family = binomial
)
############################################################
# 16.2 SELECCIÓN STEPWISE
############################################################
modelo_logit_final <- step(
modelo_logit_inicial,
direction = "both",
trace = FALSE
)
cat("\n=============================\n")##
## =============================cat("MODELO LOGÍSTICO FINAL SELECCIONADO POR STEPWISE\n")## MODELO LOGÍSTICO FINAL SELECCIONADO POR STEPWISEcat("=============================\n")## =============================print(formula(modelo_logit_final))## ESTADO_ACADEMICO ~ RESIDENCIA + CARRERA + TIPO_INGRESO + RENDIMIENTO +
## NIVEL_SOCIOprint(summary(modelo_logit_final))##
## Call:
## glm(formula = ESTADO_ACADEMICO ~ RESIDENCIA + CARRERA + TIPO_INGRESO +
## RENDIMIENTO + NIVEL_SOCIO, family = binomial, data = train_model)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 0.4701 0.4726 0.995
## RESIDENCIACENTRAL 0.3479 0.2626 1.325
## RESIDENCIARESTO DEL PAÍS -0.4734 0.3671 -1.290
## CARRERAESTADISTICA-SEMI -0.7379 0.4970 -1.485
## CARRERAMATEMATICA 0.4434 0.3177 1.396
## CARRERAEDUCACION MATEMATICA-PRES 1.2431 0.4056 3.065
## CARRERAEDUCACION MATEMATICA-SEMI -0.2644 0.3352 -0.789
## TIPO_INGRESOADMISION DIRECTA 1.7911 0.6210 2.884
## RENDIMIENTONO APROBADO -2.7763 0.3526 -7.874
## NIVEL_SOCIOMEDIO -0.8353 0.4259 -1.961
## Pr(>|z|)
## (Intercept) 0.31987
## RESIDENCIACENTRAL 0.18518
## RESIDENCIARESTO DEL PAÍS 0.19719
## CARRERAESTADISTICA-SEMI 0.13758
## CARRERAMATEMATICA 0.16275
## CARRERAEDUCACION MATEMATICA-PRES 0.00218 **
## CARRERAEDUCACION MATEMATICA-SEMI 0.43029
## TIPO_INGRESOADMISION DIRECTA 0.00392 **
## RENDIMIENTONO APROBADO 0.00000000000000344 ***
## NIVEL_SOCIOMEDIO 0.04984 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 606.02 on 478 degrees of freedom
## Residual deviance: 441.98 on 469 degrees of freedom
## AIC: 461.98
##
## Number of Fisher Scoring iterations: 5############################################################
# 16.3 TABLA DE COEFICIENTES, OR E IC95%
############################################################
coefs <- summary(modelo_logit_final)$coefficients
tabla_or <- data.frame(
Variable = rownames(coefs),
Coeficiente = coefs[, "Estimate"],
Error_Estandar = coefs[, "Std. Error"],
Estadistico_z = coefs[, "z value"],
P_Valor = coefs[, "Pr(>|z|)"],
OR = exp(coefs[, "Estimate"]),
IC95_LI = exp(coefs[, "Estimate"] - 1.96 * coefs[, "Std. Error"]),
IC95_LS = exp(coefs[, "Estimate"] + 1.96 * coefs[, "Std. Error"])
)
tabla_or <- tabla_or %>%
mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
as.data.frame()
cat("\n=============================\n")##
## =============================cat("ODDS RATIO E INTERVALOS DE CONFIANZA\n")## ODDS RATIO E INTERVALOS DE CONFIANZAcat("=============================\n")## =============================print(tabla_or)## Variable Coeficiente
## (Intercept) (Intercept) 0.4701
## RESIDENCIACENTRAL RESIDENCIACENTRAL 0.3479
## RESIDENCIARESTO DEL PAÍS RESIDENCIARESTO DEL PAÍS -0.4734
## CARRERAESTADISTICA-SEMI CARRERAESTADISTICA-SEMI -0.7379
## CARRERAMATEMATICA CARRERAMATEMATICA 0.4434
## CARRERAEDUCACION MATEMATICA-PRES CARRERAEDUCACION MATEMATICA-PRES 1.2431
## CARRERAEDUCACION MATEMATICA-SEMI CARRERAEDUCACION MATEMATICA-SEMI -0.2644
## TIPO_INGRESOADMISION DIRECTA TIPO_INGRESOADMISION DIRECTA 1.7911
## RENDIMIENTONO APROBADO RENDIMIENTONO APROBADO -2.7763
## NIVEL_SOCIOMEDIO NIVEL_SOCIOMEDIO -0.8353
## Error_Estandar Estadistico_z P_Valor OR
## (Intercept) 0.4726 0.9947 0.3199 1.6002
## RESIDENCIACENTRAL 0.2626 1.3250 0.1852 1.4161
## RESIDENCIARESTO DEL PAÍS 0.3671 -1.2896 0.1972 0.6229
## CARRERAESTADISTICA-SEMI 0.4970 -1.4849 0.1376 0.4781
## CARRERAMATEMATICA 0.3177 1.3959 0.1628 1.5581
## CARRERAEDUCACION MATEMATICA-PRES 0.4056 3.0647 0.0022 3.4662
## CARRERAEDUCACION MATEMATICA-SEMI 0.3352 -0.7887 0.4303 0.7677
## TIPO_INGRESOADMISION DIRECTA 0.6210 2.8842 0.0039 5.9958
## RENDIMIENTONO APROBADO 0.3526 -7.8738 0.0000 0.0623
## NIVEL_SOCIOMEDIO 0.4259 -1.9613 0.0498 0.4337
## IC95_LI IC95_LS
## (Intercept) 0.6337 4.0411
## RESIDENCIACENTRAL 0.8464 2.3690
## RESIDENCIARESTO DEL PAÍS 0.3034 1.2790
## CARRERAESTADISTICA-SEMI 0.1805 1.2663
## CARRERAMATEMATICA 0.8359 2.9040
## CARRERAEDUCACION MATEMATICA-PRES 1.5653 7.6755
## CARRERAEDUCACION MATEMATICA-SEMI 0.3979 1.4809
## TIPO_INGRESOADMISION DIRECTA 1.7752 20.2511
## RENDIMIENTONO APROBADO 0.0312 0.1243
## NIVEL_SOCIOMEDIO 0.1882 0.9994############################################################
# 16.4 VERIFICACIÓN BÁSICA 1: MULTICOLINEALIDAD
############################################################
# car::vif() puede devolver vector, matriz o fallar; se maneja de forma robusta
vif_raw <- tryCatch(
car::vif(modelo_logit_final),
error = function(e) NULL
)
cat("\n=============================\n")##
## =============================cat("VERIFICACIÓN DE MULTICOLINEALIDAD (VIF)\n")## VERIFICACIÓN DE MULTICOLINEALIDAD (VIF)cat("=============================\n")## =============================if (is.null(vif_raw)) {
cat("No fue posible calcular VIF con la estructura actual del modelo.\n")
vif_modelo <- data.frame(
Variable = "No disponible",
VIF = NA
)
} else {
if (is.matrix(vif_raw)) {
vif_modelo <- data.frame(
Variable = rownames(vif_raw),
GVIF = vif_raw[, 1],
Df = vif_raw[, 2],
GVIF_ajustado = vif_raw[, 3]
)
} else {
vif_modelo <- data.frame(
Variable = names(vif_raw),
VIF = as.numeric(vif_raw)
)
}
print(vif_modelo)
}## Variable GVIF Df GVIF_ajustado
## RESIDENCIA RESIDENCIA 1.055059 2 1.013489
## CARRERA CARRERA 1.462960 4 1.048707
## TIPO_INGRESO TIPO_INGRESO 1.392938 1 1.180228
## RENDIMIENTO RENDIMIENTO 1.062937 1 1.030988
## NIVEL_SOCIO NIVEL_SOCIO 1.057698 1 1.028445############################################################
# 16.5 VERIFICACIÓN BÁSICA 2: LINEALIDAD EN EL LOGIT
# (para variable numérica EDAD)
############################################################
train_bt <- train_model
train_bt$Y_BT <- ifelse(train_bt$ESTADO_ACADEMICO == "Desertor", 1, 0)
if ("EDAD" %in% names(train_bt) && all(train_bt$EDAD > 0, na.rm = TRUE)) {
modelo_bt <- glm(
Y_BT ~ EDAD + I(EDAD * log(EDAD)),
data = train_bt,
family = binomial
)
cat("\n=============================\n")
cat("VERIFICACIÓN DE LINEALIDAD EN EL LOGIT (BOX-TIDWELL PARA EDAD)\n")
cat("=============================\n")
print(summary(modelo_bt))
cat("\nInterpretación:\n")
cat("- Si el término I(EDAD * log(EDAD)) NO es significativo, la linealidad en el logit es razonable.\n")
cat("- Si es significativo, la relación puede no ser lineal en el logit.\n")
} else {
cat("\n=============================\n")
cat("VERIFICACIÓN DE LINEALIDAD EN EL LOGIT\n")
cat("=============================\n")
cat("No se aplicó Box-Tidwell porque EDAD no cumple las condiciones requeridas.\n")
}##
## =============================
## VERIFICACIÓN DE LINEALIDAD EN EL LOGIT (BOX-TIDWELL PARA EDAD)
## =============================
##
## Call:
## glm(formula = Y_BT ~ EDAD + I(EDAD * log(EDAD)), family = binomial,
## data = train_bt)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.57673 1.88423 -1.898 0.0577 .
## EDAD 0.56372 0.28034 2.011 0.0443 *
## I(EDAD * log(EDAD)) -0.12062 0.06303 -1.914 0.0557 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 606.02 on 478 degrees of freedom
## Residual deviance: 595.94 on 476 degrees of freedom
## AIC: 601.94
##
## Number of Fisher Scoring iterations: 4
##
##
## Interpretación:
## - Si el término I(EDAD * log(EDAD)) NO es significativo, la linealidad en el logit es razonable.
## - Si es significativo, la relación puede no ser lineal en el logit.############################################################
# 16.6 VERIFICACIÓN BÁSICA 3: BONDAD DE AJUSTE
# PRUEBA DE HOSMER-LEMESHOW
############################################################
prob_train_modelo <- predict(modelo_logit_final, newdata = train_model, type = "response")
medias_prob_train <- tapply(prob_train_modelo, train_model$ESTADO_ACADEMICO, mean)
if (medias_prob_train["No_Desertor"] > medias_prob_train["Desertor"]) {
prob_train_desertor <- 1 - prob_train_modelo
} else {
prob_train_desertor <- prob_train_modelo
}
y_train_num <- ifelse(train_model$ESTADO_ACADEMICO == "Desertor", 1, 0)
hl <- ResourceSelection::hoslem.test(y_train_num, prob_train_desertor, g = 10)
cat("\n=============================\n")##
## =============================cat("PRUEBA DE HOSMER-LEMESHOW\n")## PRUEBA DE HOSMER-LEMESHOWcat("=============================\n")## =============================print(hl)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: y_train_num, prob_train_desertor
## X-squared = 2.433, df = 8, p-value = 0.9648cat("\nInterpretación:\n")##
## Interpretación:cat("- p > 0.05: no hay evidencia de mal ajuste.\n")## - p > 0.05: no hay evidencia de mal ajuste.cat("- p <= 0.05: podría haber problemas de ajuste/calibración.\n")## - p <= 0.05: podría haber problemas de ajuste/calibración.############################################################
# 16.7 VERIFICACIÓN BÁSICA 4: OBSERVACIONES INFLUYENTES
# RESIDUOS DE PEARSON Y DISTANCIA DE COOK
############################################################
residuos_pearson <- residuals(modelo_logit_final, type = "pearson")
cook <- cooks.distance(modelo_logit_final)
cat("\n=============================\n")##
## =============================cat("RESUMEN DE RESIDUOS DE PEARSON\n")## RESUMEN DE RESIDUOS DE PEARSONcat("=============================\n")## =============================print(summary(residuos_pearson))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.878984 -0.576109 -0.207892 -0.003554 0.644711 5.846053cat("\n=============================\n")##
## =============================cat("RESUMEN DE DISTANCIAS DE COOK\n")## RESUMEN DE DISTANCIAS DE COOKcat("=============================\n")## =============================print(summary(cook))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000006225 0.000051886 0.001077027 0.002127516 0.002230838 0.029366789umbral_cook <- 4 / nrow(train_model)
cat("\nUmbral orientativo para Cook =", round(umbral_cook, 4), "\n")##
## Umbral orientativo para Cook = 0.0084cat("Número de observaciones con Cook > umbral:", sum(cook > umbral_cook), "\n")## Número de observaciones con Cook > umbral: 27############################################################
# 16.8 GRÁFICOS DE DIAGNÓSTICO
############################################################
par(mfrow = c(1, 2))
plot(
residuos_pearson,
pch = 20,
main = "Residuos de Pearson",
ylab = "Residuo",
xlab = "Índice"
)
abline(h = c(-2, 2), lty = 2, col = "red")
plot(
cook,
pch = 20,
main = "Distancia de Cook",
ylab = "Cook",
xlab = "Índice"
)
abline(h = umbral_cook, lty = 2, col = "red")
par(mfrow = c(1, 1))
############################################################
# 16.9 CAPACIDAD DISCRIMINANTE EN TRAIN Y TEST
############################################################
# TRAIN
prob_train_modelo <- predict(modelo_logit_final, newdata = train_model, type = "response")
medias_prob_train <- tapply(prob_train_modelo, train_model$ESTADO_ACADEMICO, mean)
if (medias_prob_train["No_Desertor"] > medias_prob_train["Desertor"]) {
prob_train_desertor <- 1 - prob_train_modelo
} else {
prob_train_desertor <- prob_train_modelo
}
roc_train_logit <- roc(
response = train_model$ESTADO_ACADEMICO,
predictor = prob_train_desertor,
levels = c("No_Desertor", "Desertor"),
direction = "<"
)
auc_train_logit <- as.numeric(auc(roc_train_logit))
# TEST
prob_test_modelo <- predict(modelo_logit_final, newdata = test_model, type = "response")
medias_prob_test <- tapply(prob_test_modelo, test_model$ESTADO_ACADEMICO, mean)
if (medias_prob_test["No_Desertor"] > medias_prob_test["Desertor"]) {
prob_test_desertor <- 1 - prob_test_modelo
} else {
prob_test_desertor <- prob_test_modelo
}
roc_test_logit <- roc(
response = test_model$ESTADO_ACADEMICO,
predictor = prob_test_desertor,
levels = c("No_Desertor", "Desertor"),
direction = "<"
)
auc_test_logit <- as.numeric(auc(roc_test_logit))
cat("\n=============================\n")##
## =============================cat("AUC EN TRAIN Y TEST\n")## AUC EN TRAIN Y TESTcat("=============================\n")## =============================cat("AUC TRAIN:", round(auc_train_logit, 4), "\n")## AUC TRAIN: 0.8298cat("AUC TEST :", round(auc_test_logit, 4), "\n")## AUC TEST : 0.7635cat("Diferencia absoluta:", round(abs(auc_train_logit - auc_test_logit), 4), "\n")## Diferencia absoluta: 0.0663############################################################
# 16.10 MATRIZ DE CONFUSIÓN FINAL EN TEST (UMBRAL 0.50)
############################################################
umbral_final <- 0.50
pred_test_final <- ifelse(prob_test_desertor >= umbral_final, "Desertor", "No_Desertor")
pred_test_final <- factor(pred_test_final, levels = c("Desertor", "No_Desertor"))
cm_logit_final <- confusionMatrix(
data = pred_test_final,
reference = test_model$ESTADO_ACADEMICO,
positive = "Desertor"
)
tabla_logit_final <- data.frame(
Modelo = "Regresión logística final",
Formula_Final = paste(deparse(formula(modelo_logit_final)), collapse = " "),
Umbral = umbral_final,
AUC_Train = auc_train_logit,
AUC_Test = auc_test_logit,
Accuracy = as.numeric(cm_logit_final$overall["Accuracy"]),
Kappa = as.numeric(cm_logit_final$overall["Kappa"]),
Sensibilidad = as.numeric(cm_logit_final$byClass["Sensitivity"]),
Especificidad = as.numeric(cm_logit_final$byClass["Specificity"]),
Precision = as.numeric(cm_logit_final$byClass["Pos Pred Value"]),
VPN = as.numeric(cm_logit_final$byClass["Neg Pred Value"]),
Balanced_Accuracy = as.numeric(cm_logit_final$byClass["Balanced Accuracy"])
)
tabla_logit_final <- tabla_logit_final %>%
mutate(across(-c(1, 2), ~ if(is.numeric(.x)) round(.x, 4) else .x)) %>%
as.data.frame()
cat("\n=============================\n")##
## =============================cat("MATRIZ DE CONFUSIÓN FINAL EN TEST\n")## MATRIZ DE CONFUSIÓN FINAL EN TESTcat("=============================\n")## =============================print(cm_logit_final)## Confusion Matrix and Statistics
##
## Reference
## Prediction Desertor No_Desertor
## Desertor 112 35
## No_Desertor 25 31
##
## Accuracy : 0.7044
## 95% CI : (0.6365, 0.7663)
## No Information Rate : 0.6749
## P-Value [Acc > NIR] : 0.2058
##
## Kappa : 0.299
##
## Mcnemar's Test P-Value : 0.2453
##
## Sensitivity : 0.8175
## Specificity : 0.4697
## Pos Pred Value : 0.7619
## Neg Pred Value : 0.5536
## Prevalence : 0.6749
## Detection Rate : 0.5517
## Detection Prevalence : 0.7241
## Balanced Accuracy : 0.6436
##
## 'Positive' Class : Desertor
## cat("\n=============================\n")##
## =============================cat("TABLA RESUMEN DEL MODELO FINAL\n")## TABLA RESUMEN DEL MODELO FINALcat("=============================\n")## =============================print(tabla_logit_final)## Modelo
## 1 Regresión logística final
## Formula_Final
## 1 ESTADO_ACADEMICO ~ RESIDENCIA + CARRERA + TIPO_INGRESO + RENDIMIENTO + NIVEL_SOCIO
## Umbral AUC_Train AUC_Test Accuracy Kappa Sensibilidad Especificidad Precision
## 1 0.5 0.8298 0.7635 0.7044 0.299 0.8175 0.4697 0.7619
## VPN Balanced_Accuracy
## 1 0.5536 0.6436############################################################
# 16.11 CURVA ROC FINAL EN TEST
############################################################
plot(
roc_test_logit,
main = paste0("Curva ROC - Modelo logístico final (AUC = ", round(auc_test_logit, 3), ")"),
lwd = 2
)
abline(a = 0, b = 1, lty = 2)
############################################################
# 16.12 ANÁLISIS DEL TRADE-OFF SEGÚN EL UMBRAL
############################################################
umbrales <- seq(0.01, 0.99, by = 0.01)
metricas_umbral <- data.frame()
for (t in umbrales) {
pred_t <- ifelse(prob_test_desertor >= t, "Desertor", "No_Desertor")
pred_t <- factor(pred_t, levels = c("Desertor", "No_Desertor"))
cm_t <- confusionMatrix(
data = pred_t,
reference = test_model$ESTADO_ACADEMICO,
positive = "Desertor"
)
metricas_umbral <- rbind(metricas_umbral, data.frame(
Umbral = t,
Sensitivity = as.numeric(cm_t$byClass["Sensitivity"]),
Specificity = as.numeric(cm_t$byClass["Specificity"]),
Accuracy = as.numeric(cm_t$overall["Accuracy"]),
Precision = as.numeric(cm_t$byClass["Pos Pred Value"]),
VPN = as.numeric(cm_t$byClass["Neg Pred Value"]),
Balanced_Accuracy = as.numeric(cm_t$byClass["Balanced Accuracy"])
))
}
metricas_umbral$Youden <- metricas_umbral$Sensitivity + metricas_umbral$Specificity - 1
fila_youden <- which.max(metricas_umbral$Youden)
umbral_youden <- metricas_umbral$Umbral[fila_youden]
fila_050 <- which.min(abs(metricas_umbral$Umbral - 0.50))
cat("\n=============================\n")##
## =============================cat("ANÁLISIS DE UMBRALES - MODELO LOGÍSTICO FINAL\n")## ANÁLISIS DE UMBRALES - MODELO LOGÍSTICO FINALcat("=============================\n")## =============================cat("Umbral fijo utilizado en el documento:", 0.50, "\n")## Umbral fijo utilizado en el documento: 0.5cat("Umbral óptimo según Youden:", round(umbral_youden, 3), "\n")## Umbral óptimo según Youden: 0.7print(metricas_umbral[c(fila_050, fila_youden), ])## Umbral Sensitivity Specificity Accuracy Precision VPN
## 50 0.5 0.8175182 0.4696970 0.7044335 0.7619048 0.5535714
## 70 0.7 0.6058394 0.8636364 0.6896552 0.9021739 0.5135135
## Balanced_Accuracy Youden
## 50 0.6436076 0.2872152
## 70 0.7347379 0.4694758############################################################
# 16.13 TABLA RESUMEN DE UMBRALES CLAVE
############################################################
tabla_umbrales_clave <- rbind(
data.frame(
Criterio = "Umbral fijo 0.50",
Umbral = metricas_umbral$Umbral[fila_050],
Sensibilidad = metricas_umbral$Sensitivity[fila_050],
Especificidad = metricas_umbral$Specificity[fila_050],
Accuracy = metricas_umbral$Accuracy[fila_050],
Youden = metricas_umbral$Youden[fila_050],
Balanced_Accuracy = metricas_umbral$Balanced_Accuracy[fila_050]
),
data.frame(
Criterio = "Óptimo por Youden",
Umbral = metricas_umbral$Umbral[fila_youden],
Sensibilidad = metricas_umbral$Sensitivity[fila_youden],
Especificidad = metricas_umbral$Specificity[fila_youden],
Accuracy = metricas_umbral$Accuracy[fila_youden],
Youden = metricas_umbral$Youden[fila_youden],
Balanced_Accuracy = metricas_umbral$Balanced_Accuracy[fila_youden]
)
)
tabla_umbrales_clave <- tabla_umbrales_clave %>%
mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
as.data.frame()
cat("\n=============================\n")##
## =============================cat("TABLA DE UMBRALES CLAVE\n")## TABLA DE UMBRALES CLAVEcat("=============================\n")## =============================print(tabla_umbrales_clave)## Criterio Umbral Sensibilidad Especificidad Accuracy Youden
## 1 Umbral fijo 0.50 0.5 0.8175 0.4697 0.7044 0.2872
## 2 Óptimo por Youden 0.7 0.6058 0.8636 0.6897 0.4695
## Balanced_Accuracy
## 1 0.6436
## 2 0.7347############################################################
# 16.14 GRÁFICO 1: SCATTER DEL TRADE-OFF
############################################################
graf_tradeoff_scatter <- ggplot(
metricas_umbral,
aes(x = Specificity, y = Sensitivity)
) +
geom_point(aes(size = Accuracy, color = Umbral), alpha = 0.80) +
scale_color_gradient(low = "#132B43", high = "#56B1F7") +
geom_point(
data = metricas_umbral[fila_050, , drop = FALSE],
shape = 21, fill = "white", color = "black", size = 6, stroke = 1.5
) +
labs(
title = "Trade-off entre sensibilidad y especificidad según el umbral",
subtitle = "Cada punto representa un punto de corte distinto",
x = "Especificidad",
y = "Sensibilidad",
color = "Umbral",
size = "Accuracy"
) +
theme_minimal(base_size = 13)
print(graf_tradeoff_scatter)
############################################################
# 16.15 GRÁFICO 2: CURVAS DE SENSIBILIDAD Y ESPECIFICIDAD
############################################################
metricas_long_umbral <- tidyr::pivot_longer(
metricas_umbral,
cols = c("Sensitivity", "Specificity"),
names_to = "Metrica",
values_to = "Valor"
)
metricas_long_umbral$Metrica <- factor(
metricas_long_umbral$Metrica,
levels = c("Sensitivity", "Specificity"),
labels = c("Sensitivity", "Specificity")
)
graf_tradeoff_lineas <- ggplot(
metricas_long_umbral,
aes(x = Umbral, y = Valor, color = Metrica)
) +
geom_line(linewidth = 1.4) +
geom_vline(xintercept = metricas_umbral$Umbral[fila_050],
linetype = "dashed", linewidth = 0.8) +
geom_vline(xintercept = umbral_youden,
linetype = "dotted", linewidth = 0.8) +
annotate("text",
x = metricas_umbral$Umbral[fila_050] - 0.01,
y = 0.08,
label = "0.50",
angle = 90,
vjust = -0.4,
size = 5) +
annotate("text",
x = umbral_youden - 0.01,
y = 0.08,
label = "Youden",
angle = 90,
vjust = -0.4,
size = 5) +
labs(
title = "Trade-off entre sensibilidad y especificidad",
subtitle = "Evaluación del efecto del punto de corte sobre el desempeño del modelo logístico",
x = "Umbral de clasificación",
y = "Valor de la métrica",
color = "Métrica"
) +
theme_minimal(base_size = 13)
print(graf_tradeoff_lineas)
############################################################
# 16.16 DECISIÓN OPERATIVA DEL UMBRAL
############################################################
decision_umbral <- data.frame(
Criterio = c("Modelo final reportado", "Referencia alternativa"),
Umbral = c(0.50, umbral_youden),
Justificacion = c(
"Se mantiene 0.50 por consistencia interpretativa y por priorizar sensibilidad en detección de desertores",
"Maximiza el índice de Youden y representa el mejor equilibrio global Sensibilidad-Especificidad"
)
)
cat("\n=============================\n")##
## =============================cat("DECISIÓN OPERATIVA DEL UMBRAL\n")## DECISIÓN OPERATIVA DEL UMBRALcat("=============================\n")## =============================print(decision_umbral)## Criterio Umbral
## 1 Modelo final reportado 0.5
## 2 Referencia alternativa 0.7
## Justificacion
## 1 Se mantiene 0.50 por consistencia interpretativa y por priorizar sensibilidad en detección de desertores
## 2 Maximiza el índice de Youden y representa el mejor equilibrio global Sensibilidad-Especificidad############################################################
# 17. IMPORTANCIA DE VARIABLES - RANDOM FOREST
############################################################
imp_rf <- varImp(modelo_rf_cv)
print(imp_rf)## rf variable importance
##
## only 20 most important variables shown (out of 29)
##
## Overall
## RENDIMIENTONO APROBADO 100.000
## EDAD 48.632
## CARRERAEDUCACION MATEMATICA-PRES 30.402
## TRABAJANO 18.016
## MODALIDADSemipresencial 17.647
## SOLVENTARAyuda Familiar 13.319
## RESIDENCIACENTRAL 12.720
## CARRERAMATEMATICA 12.529
## SOLVENTARTrabajo Personal 11.634
## SEXOMASCULINO 11.448
## ESTUDIOS_PADRES30-34 años 10.814
## ESTUDIOS_PADRES14-23 años 10.383
## NIVEL_SOCIOMEDIO 10.221
## CARRERAEDUCACION MATEMATICA-SEMI 10.175
## RESIDENCIARESTO DEL PAÍS 9.945
## ESTUDIOS_PADRES24-29 años 9.900
## TIPO_COLPrivado 9.302
## INGRESOMás de dos y hasta cinco salarios mínimos 9.205
## CARRERAESTADISTICA-SEMI 9.043
## TIPO_INGRESOADMISION DIRECTA 7.978plot(imp_rf, top = 20, main = "Importancia de variables - Random Forest")
############################################################
# 18. EXPORTACIÓN DE RESULTADOS
############################################################
dir.create("salidas_tesis", showWarnings = FALSE)
############################################################
# 18.1 EXPORTACIÓN DE TABLAS EN CSV
############################################################
write.csv(resultados_chi,
"salidas_tesis/tabla_chi_fisher.csv",
row.names = FALSE)
write.csv(resultados_bivariados,
"salidas_tesis/tabla_logistica_bivariada.csv",
row.names = FALSE)
write.csv(tabla_metricas_test,
"salidas_tesis/tabla_metricas_test_modelos.csv",
row.names = FALSE)
write.csv(metricas_graf,
"salidas_tesis/tabla_6_metricas_modelos_test.csv",
row.names = FALSE)
write.csv(tabla_or,
"salidas_tesis/tabla_or_modelo_logistico_final.csv",
row.names = FALSE)
# Exportación directa del objeto VIF ya estructurado
write.csv(vif_modelo,
"salidas_tesis/tabla_vif_modelo_logistico_final.csv",
row.names = FALSE)
write.csv(tabla_logit_final,
"salidas_tesis/tabla_modelo_logistico_final.csv",
row.names = FALSE)
write.csv(metricas_umbral,
"salidas_tesis/tabla_metricas_por_umbral_logit.csv",
row.names = FALSE)
write.csv(tabla_umbrales_clave,
"salidas_tesis/tabla_umbrales_clave_logit.csv",
row.names = FALSE)
write.csv(decision_umbral,
"salidas_tesis/tabla_decision_umbral_logit.csv",
row.names = FALSE)
############################################################
# 18.2 EXPORTACIÓN A EXCEL CON FORMATO TIPO TESIS
############################################################
wb <- createWorkbook(creator = "ChatGPT - Tesis de Estadística")
# -----------------------------
# COLORES Y ESTILOS
# -----------------------------
color_azul_oscuro <- "#1F4E78"
color_gris_claro <- "#F4F6F8"
color_gris_borde <- "#B8C2CC"
color_blanco <- "#FFFFFF"
color_seccion <- "#44546A"
estilo_titulo <- createStyle(
fontSize = 14,
textDecoration = "bold",
fontColour = color_blanco,
fgFill = color_azul_oscuro,
halign = "left",
valign = "center"
)
estilo_seccion <- createStyle(
fontSize = 11,
textDecoration = "bold",
fontColour = color_blanco,
fgFill = color_seccion,
halign = "left",
valign = "center",
border = "TopBottomLeftRight",
borderColour = color_seccion
)
estilo_nota <- createStyle(
fontSize = 10,
fontColour = "#5B6570",
textDecoration = "italic",
halign = "left"
)
estilo_header <- createStyle(
fontSize = 11,
textDecoration = "bold",
fontColour = color_blanco,
fgFill = color_azul_oscuro,
halign = "center",
valign = "center",
border = "TopBottomLeftRight",
borderColour = color_gris_borde
)
estilo_texto_par <- createStyle(
fgFill = color_blanco,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
valign = "center"
)
estilo_texto_impar <- createStyle(
fgFill = color_gris_claro,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
valign = "center"
)
estilo_num_par <- createStyle(
fgFill = color_blanco,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center",
numFmt = "0.0000"
)
estilo_num_impar <- createStyle(
fgFill = color_gris_claro,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center",
numFmt = "0.0000"
)
estilo_entero_par <- createStyle(
fgFill = color_blanco,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center",
numFmt = "0"
)
estilo_entero_impar <- createStyle(
fgFill = color_gris_claro,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center",
numFmt = "0"
)
estilo_porcentaje_par <- createStyle(
fgFill = color_blanco,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center",
numFmt = "0.00"
)
estilo_porcentaje_impar <- createStyle(
fgFill = color_gris_claro,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center",
numFmt = "0.00"
)
estilo_texto_centrado_par <- createStyle(
fgFill = color_blanco,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center"
)
estilo_texto_centrado_impar <- createStyle(
fgFill = color_gris_claro,
border = "TopBottomLeftRight",
borderColour = color_gris_borde,
halign = "center",
valign = "center"
)
es_numerica <- function(x) is.numeric(x) || is.integer(x)
aplicar_estilos_fila <- function(wb, hoja, datos, fila_inicio = 4) {
if (nrow(datos) == 0) return(invisible(NULL))
for (i in seq_len(nrow(datos))) {
fila_excel <- fila_inicio + i
impar <- i %% 2 == 1
for (j in seq_along(datos)) {
valor <- datos[[j]][i]
nombre_col <- names(datos)[j]
estilo <- NULL
if (es_numerica(valor)) {
if (grepl("porcentaje|sens|spec|auc|accuracy|kappa|precision|f1|vpn|youden|or|ic95|coef|estadistico|umbral", tolower(nombre_col))) {
estilo <- if (impar) estilo_num_impar else estilo_num_par
} else if (abs(valor - round(valor)) < .Machine$double.eps^0.5) {
estilo <- if (impar) estilo_entero_impar else estilo_entero_par
} else {
estilo <- if (impar) estilo_num_impar else estilo_num_par
}
} else {
if (grepl("modelo|variable|categoria|criterio|metodo|significancia|formula|nivel|justificacion|estado", tolower(nombre_col))) {
estilo <- if (impar) estilo_texto_impar else estilo_texto_par
} else {
estilo <- if (impar) estilo_texto_centrado_impar else estilo_texto_centrado_par
}
}
addStyle(wb, hoja, style = estilo, rows = fila_excel, cols = j, gridExpand = FALSE, stack = TRUE)
}
}
}
exportar_hoja_tesis <- function(wb, nombre_hoja, datos, titulo, seccion = "Resultados exportados", nota = "Fuente: Elaboración propia") {
addWorksheet(wb, nombre_hoja)
writeData(wb, nombre_hoja, titulo, startRow = 1, startCol = 1)
mergeCells(wb, nombre_hoja, cols = 1:max(1, ncol(datos)), rows = 1)
addStyle(wb, nombre_hoja, estilo_titulo, rows = 1, cols = 1, gridExpand = TRUE)
setRowHeights(wb, nombre_hoja, rows = 1, heights = 24)
writeData(wb, nombre_hoja, seccion, startRow = 2, startCol = 1)
mergeCells(wb, nombre_hoja, cols = 1:max(1, ncol(datos)), rows = 2)
addStyle(wb, nombre_hoja, estilo_seccion, rows = 2, cols = 1, gridExpand = TRUE)
writeData(wb, nombre_hoja, nota, startRow = 3, startCol = 1)
mergeCells(wb, nombre_hoja, cols = 1:max(1, ncol(datos)), rows = 3)
addStyle(wb, nombre_hoja, estilo_nota, rows = 3, cols = 1, gridExpand = TRUE)
writeData(wb, nombre_hoja, datos, startRow = 4, startCol = 1, headerStyle = estilo_header, withFilter = TRUE)
aplicar_estilos_fila(wb, nombre_hoja, datos, fila_inicio = 4)
freezePane(wb, nombre_hoja, firstActiveRow = 5, firstActiveCol = 1)
setColWidths(wb, nombre_hoja, cols = 1:max(1, ncol(datos)), widths = "auto")
}
# Hoja índice
indice <- data.frame(
Hoja = c(
"01_Chi_Fisher", "02_Logistica_Bivariada", "03_Metricas_Test_Modelos",
"04_Metricas_6_Modelos", "05_OR_Modelo_Final", "06_VIF_Modelo_Final",
"07_Modelo_Logistico_Final", "08_Metricas_por_Umbral",
"09_Umbrales_Clave", "10_Decision_Umbral"
),
Descripcion = c(
"Resultados de pruebas Chi-cuadrado y Fisher",
"Regresión logística simple para variables numéricas",
"Métricas finales en el conjunto de prueba",
"Tabla de apoyo para gráfico comparativo de 6 métricas",
"Coeficientes, OR e IC95% del modelo logístico final",
"Diagnóstico de multicolinealidad del modelo final",
"Resumen del modelo logístico final en test",
"Desempeño del modelo logístico según el umbral",
"Comparación entre umbral 0.50 y umbral óptimo de Youden",
"Criterio final adoptado para el punto de corte"
),
stringsAsFactors = FALSE
)
exportar_hoja_tesis(
wb, "00_Indice", indice,
titulo = "Índice del libro de resultados",
seccion = "Resumen de hojas del archivo Excel"
)
exportar_hoja_tesis(wb, "01_Chi_Fisher", resultados_chi,
titulo = "Pruebas de asociación bivariante",
seccion = "Resultados Chi-cuadrado y Prueba Exacta de Fisher")
exportar_hoja_tesis(wb, "02_Logistica_Bivariada", resultados_bivariados,
titulo = "Regresión logística bivariada",
seccion = "Resultados para covariables numéricas")
exportar_hoja_tesis(wb, "03_Metricas_Test_Modelos", tabla_metricas_test,
titulo = "Métricas de desempeño en el conjunto de prueba",
seccion = "Comparación de modelos")
exportar_hoja_tesis(wb, "04_Metricas_6_Modelos", metricas_graf,
titulo = "Tabla base del gráfico de 6 métricas",
seccion = "Accuracy, Sensitivity, Specificity, Precision, F1 y AUC")
exportar_hoja_tesis(wb, "05_OR_Modelo_Final", tabla_or,
titulo = "Odds ratio del modelo logístico final",
seccion = "Coeficientes e intervalos de confianza")
exportar_hoja_tesis(wb, "06_VIF_Modelo_Final", vif_modelo,
titulo = "Diagnóstico de multicolinealidad",
seccion = "Factores de inflación de la varianza")
exportar_hoja_tesis(wb, "07_Modelo_Logistico_Final", tabla_logit_final,
titulo = "Resumen del modelo logístico final",
seccion = "Desempeño final en el conjunto de prueba")
exportar_hoja_tesis(wb, "08_Metricas_por_Umbral", metricas_umbral,
titulo = "Métricas del modelo logístico según el umbral",
seccion = "Evaluación del trade-off entre sensibilidad y especificidad")
exportar_hoja_tesis(wb, "09_Umbrales_Clave", tabla_umbrales_clave,
titulo = "Comparación de umbrales clave",
seccion = "Umbral fijo 0.50 y óptimo por Youden")
exportar_hoja_tesis(wb, "10_Decision_Umbral", decision_umbral,
titulo = "Decisión operativa del umbral",
seccion = "Criterio adoptado para la clasificación final")
saveWorkbook(wb, "salidas_tesis/resumen_resultados_tesis.xlsx", overwrite = TRUE)
############################################################
# 18.3 EXPORTACIÓN DE GRÁFICOS (ggplot)
############################################################
# Análisis descriptivo
ggsave(
filename = "salidas_tesis/grafico_estado_academico.png",
plot = graf_estado,
width = 8,
height = 6,
dpi = 300
)
ggsave(
filename = "salidas_tesis/boxplot_edad_estado_academico.png",
plot = graf_edad,
width = 8,
height = 6,
dpi = 300
)
# Comparación de modelos
ggsave(
filename = "salidas_tesis/grafico_6_metricas_modelos_test.png",
plot = graf_6_metricas,
width = 14,
height = 8,
dpi = 300
)
# Trade-off del modelo logístico final
ggsave(
filename = "salidas_tesis/scatter_tradeoff_umbral_modelo_final.png",
plot = graf_tradeoff_scatter,
width = 12,
height = 7,
dpi = 300
)
ggsave(
filename = "salidas_tesis/tradeoff_sens_specificity_modelo_final.png",
plot = graf_tradeoff_lineas,
width = 12,
height = 7,
dpi = 300
)
############################################################
# 18.4 EXPORTACIÓN DE GRÁFICOS BASE R
############################################################
# ROC comparativa de modelos
png("salidas_tesis/curvas_roc_modelos_test.png", width = 1400, height = 1000, res = 150)
plot(resultados_test_lista$Logit$roc,
main = "Curvas ROC de los modelos en el conjunto de prueba",
lwd = 2)
plot(resultados_test_lista$Tree$roc, add = TRUE, lty = 2, lwd = 2)
plot(resultados_test_lista$RF$roc, add = TRUE, lty = 3, lwd = 2)
plot(resultados_test_lista$SVM$roc, add = TRUE, lty = 4, lwd = 2)
plot(resultados_test_lista$KNN$roc, add = TRUE, lty = 5, lwd = 2)
plot(resultados_test_lista$NN$roc, add = TRUE, lty = 6, lwd = 2)
plot(resultados_test_lista$NB$roc, add = TRUE, lty = 7, lwd = 2)
legend(
"bottomright",
legend = c(
paste0("Logit (AUC=", round(resultados_test_lista$Logit$resumen$AUC, 3), ")"),
paste0("Tree (AUC=", round(resultados_test_lista$Tree$resumen$AUC, 3), ")"),
paste0("RF (AUC=", round(resultados_test_lista$RF$resumen$AUC, 3), ")"),
paste0("SVM (AUC=", round(resultados_test_lista$SVM$resumen$AUC, 3), ")"),
paste0("KNN (AUC=", round(resultados_test_lista$KNN$resumen$AUC, 3), ")"),
paste0("NN (AUC=", round(resultados_test_lista$NN$resumen$AUC, 3), ")"),
paste0("NB (AUC=", round(resultados_test_lista$NB$resumen$AUC, 3), ")")
),
lty = 1:7,
lwd = 2,
cex = 0.8
)
dev.off()## png
## 2# ROC del modelo logístico final
png("salidas_tesis/roc_modelo_logistico_final.png", width = 1200, height = 900, res = 150)
plot(
roc_test_logit,
main = paste0("Curva ROC - Modelo logístico final (AUC = ", round(auc_test_logit, 3), ")"),
lwd = 2
)
abline(a = 0, b = 1, lty = 2)
dev.off()## png
## 2# Importancia de variables en Random Forest
png("salidas_tesis/importancia_variables_random_forest.png", width = 1200, height = 900, res = 150)
plot(imp_rf, top = 20, main = "Importancia de variables - Random Forest")
dev.off()## png
## 2# Diagnóstico del modelo logístico final
png("salidas_tesis/diagnostico_modelo_logistico_final.png", width = 1400, height = 700, res = 150)
par(mfrow = c(1, 2))
plot(
residuos_pearson,
pch = 20,
main = "Residuos de Pearson",
ylab = "Residuo",
xlab = "Índice"
)
abline(h = c(-2, 2), lty = 2, col = "red")
plot(
cook,
pch = 20,
main = "Distancia de Cook",
ylab = "Cook",
xlab = "Índice"
)
abline(h = umbral_cook, lty = 2, col = "red")
par(mfrow = c(1, 1))
dev.off()## png
## 2############################################################
# 18.5 MENSAJE FINAL
############################################################
cat("\n=============================\n")##
## =============================cat("RESULTADOS EXPORTADOS CORRECTAMENTE\n")## RESULTADOS EXPORTADOS CORRECTAMENTEcat("=============================\n")## =============================cat("Archivos generados en la carpeta: salidas_tesis\n")## Archivos generados en la carpeta: salidas_tesiscat("- Tablas CSV\n")## - Tablas CSVcat("- Archivo Excel consolidado\n")## - Archivo Excel consolidadocat("- Gráficos PNG\n")## - Gráficos PNGcat("=============================\n")## =============================############################################################
# GRAFICO PANEL MEJORADO: COMPARACION DE METRICAS POR MODELO
############################################################
library(ggplot2)
library(dplyr)
library(tidyr)
library(grid)
metricas_plot <- tabla_metricas_test %>%
dplyr::select(Modelo, Accuracy, Precision, Sensitivity, Specificity, F1, AUC) %>%
tidyr::pivot_longer(
cols = -Modelo,
names_to = "Metrica",
values_to = "Valor"
) %>%
mutate(
Metrica = dplyr::recode(Metrica, "Sensitivity" = "Recall"),
Metrica = factor(
Metrica,
levels = c("Accuracy", "Precision", "Recall", "Specificity", "F1", "AUC")
),
Modelo = factor(
Modelo,
levels = c("Tree", "KNN", "NB", "RF", "NN", "Logit", "SVM"),
labels = c("Árbol", "KNN", "NB", "RF", "RN", "Logit", "SVM")
)
)
colores_modelos <- c(
"Árbol" = "#C75D00",
"KNN" = "#6EA61A",
"NB" = "#A8771A",
"RF" = "#6D6AAE",
"RN" = "#D9A400",
"Logit" = "#1F9D7A",
"SVM" = "#D92B8A"
)
grafico_metricas_modelos <- ggplot(
metricas_plot,
aes(x = Modelo, y = Valor, fill = Modelo)
) +
geom_col(width = 0.68, show.legend = FALSE) +
geom_text(
aes(label = sprintf("%.3f", Valor)),
vjust = -0.22,
size = 4.4,
fontface = "bold"
) +
facet_wrap(~ Metrica, ncol = 3) +
scale_fill_manual(values = colores_modelos) +
scale_y_continuous(
limits = c(0, 1.08),
breaks = c(0, 0.25, 0.50, 0.75, 1.00),
labels = sprintf("%.2f", c(0, 0.25, 0.50, 0.75, 1.00)),
expand = expansion(mult = c(0, 0.04))
) +
labs(
title = "Comparación de métricas de desempeño por modelo",
x = "Modelo",
y = "Valor"
) +
theme_minimal(base_size = 15) +
theme(
plot.title = element_text(
face = "bold",
hjust = 0.5,
size = 22,
margin = ggplot2::margin(b = 10)
),
axis.text.x = element_text(
angle = 32,
hjust = 1,
size = 14,
face = "bold"
),
axis.text.y = element_text(size = 13),
axis.title.x = element_text(face = "bold", size = 20),
axis.title.y = element_text(face = "bold", size = 20),
strip.text = element_text(face = "bold", size = 18),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_line(
linewidth = 0.35,
colour = "#D9D9D9"
),
panel.spacing = grid::unit(1.4, "lines"),
plot.margin = ggplot2::margin(15, 20, 15, 20)
)
print(grafico_metricas_modelos)
dir.create("salidas_tesis", showWarnings = FALSE)
ggsave(
filename = "salidas_tesis/metricas_modelos_panel_presentacion_mejorada.png",
plot = grafico_metricas_modelos,
width = 16,
height = 9.5,
dpi = 300
)############################################################
# 16. FIGURA SCATTER: COMPARACION INTEGRADA DEL DESEMPENO
# DE LOS MODELOS EN EL CONJUNTO DE PRUEBA
############################################################
library(dplyr)
library(ggrepel)
datos_scatter <- tabla_metricas_test %>%
mutate(
Modelo = factor(
Modelo,
levels = c("NB", "Tree", "KNN", "SVM", "NN", "Logit", "RF")
),
Tipo = case_when(
Modelo %in% c("NN", "Logit", "RF") ~ "Optimo",
Modelo %in% c("Tree", "KNN") ~ "Intermedio",
Modelo %in% c("NB", "SVM") ~ "Sesgado",
TRUE ~ "Intermedio"
),
Specificity_plot = ifelse(Specificity == 0, 0.03, Specificity)
)
colores_tipo <- c(
"Intermedio" = "#E3B505",
"Optimo" = "#2FBF71",
"Sesgado" = "#E74C3C"
)
grafico_scatter_modelos_mejorado <- ggplot(
datos_scatter,
aes(
x = Specificity_plot,
y = Sensitivity,
size = AUC,
fill = Tipo,
label = Modelo
)
) +
geom_point(
shape = 21,
colour = "white",
stroke = 0.9,
alpha = 0.95
) +
ggrepel::geom_text_repel(
aes(color = Tipo),
size = 6,
fontface = "bold",
show.legend = FALSE,
box.padding = 0.35,
point.padding = 0.25,
segment.color = "grey60",
segment.size = 0.4,
max.overlaps = Inf
) +
scale_fill_manual(values = colores_tipo, name = "Tipo de modelo") +
scale_color_manual(values = colores_tipo, guide = "none") +
scale_size_continuous(
name = "AUC",
range = c(8, 16),
breaks = c(0.76, 0.77, 0.78, 0.79),
labels = c("0.76", "0.77", "0.78", "0.79")
) +
scale_x_continuous(
breaks = c(0.00, 0.10, 0.20, 0.30, 0.40, 0.50),
labels = c("0.00", "0.10", "0.20", "0.30", "0.40", "0.50"),
expand = expansion(mult = c(0.02, 0.03))
) +
coord_cartesian(
xlim = c(0.00, 0.56),
ylim = c(0.83, 0.985),
clip = "off"
) +
labs(
title = "Comparación integrada del desempeño de modelos",
subtitle = "Relación entre sensibilidad, especificidad y AUC (conjunto de prueba)",
x = "Especificidad",
y = "Sensibilidad"
) +
theme_minimal(base_size = 16) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 24),
plot.subtitle = element_text(
hjust = 0.5,
size = 17,
margin = ggplot2::margin(b = 12)
),
axis.title = element_text(face = "bold", size = 20),
axis.text = element_text(size = 14),
legend.title = element_text(face = "bold", size = 16),
legend.text = element_text(size = 14),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(linewidth = 0.35, colour = "#D9D9D9"),
plot.margin = ggplot2::margin(15, 70, 15, 15)
) +
guides(
fill = guide_legend(
order = 1,
override.aes = list(
shape = 21,
size = 5,
colour = "white",
fill = c("#E3B505", "#2FBF71", "#E74C3C"),
alpha = 1
)
),
size = guide_legend(
order = 2,
override.aes = list(
shape = 21,
fill = "grey20",
colour = "grey20",
alpha = 1
)
)
)
print(grafico_scatter_modelos_mejorado)
ggsave(
filename = "salidas_tesis/Figura_scatter_modelos_color_mejorado.png",
plot = grafico_scatter_modelos_mejorado,
width = 14,
height = 8.5,
dpi = 300
)############################################################
# TABLA 4.10 - VERSIÓN 2 CORREGIDA
# REGRESIÓN LOGÍSTICA SIMPLE CON 13 VARIABLES CATEGÓRICAS
# FORMATO TIPO TESIS (BLOQUES POR VARIABLE)
############################################################
# =========================
# 0. PAQUETES
# =========================
paquetes <- c("dplyr", "purrr", "stringr", "openxlsx", "flextable")
instalar <- paquetes[!paquetes %in% installed.packages()[, "Package"]]
if (length(instalar) > 0) install.packages(instalar)
library(dplyr)
library(purrr)##
## Adjuntando el paquete: 'purrr'## The following object is masked from 'package:scales':
##
## discard## The following object is masked from 'package:car':
##
## some## The following object is masked from 'package:caret':
##
## liftlibrary(stringr)
library(openxlsx)
library(flextable)##
## Adjuntando el paquete: 'flextable'## The following object is masked from 'package:purrr':
##
## compose# =========================
# 1. BASE DE DATOS
# =========================
# Se asume que ya existe train_data
datos <- train_data
# =========================
# 2. VARIABLE RESPUESTA
# Desertor = 1 ; No_Desertor = 0
# =========================
datos <- datos %>%
mutate(
ESTADO_BIN = ifelse(ESTADO_ACADEMICO == "Desertor", 1, 0)
)
# =========================
# 3. VARIABLES CATEGÓRICAS
# =========================
variables_categoricas <- c(
"SOLVENTAR",
"RESIDENCIA",
"CARRERA",
"TIPO_INGRESO",
"RENDIMIENTO",
"MODALIDAD",
"TRABAJA",
"ESTUDIOS_PADRES",
"SEXO",
"TIPO_COL",
"INGRESO",
"EST_CIVIL",
"NIVEL_SOCIO"
)
# =========================
# 4. ETIQUETAS DE VARIABLES
# =========================
etiquetas_variables <- c(
SOLVENTAR = "Solventar",
RESIDENCIA = "Residencia",
CARRERA = "Carrera",
TIPO_INGRESO = "Tipo_Ingreso",
RENDIMIENTO = "Rendimiento",
MODALIDAD = "Modalidad",
TRABAJA = "Trabaja",
ESTUDIOS_PADRES = "Estudios_Padres",
SEXO = "Sexo",
TIPO_COL = "Tipo_Col",
INGRESO = "Ingreso",
EST_CIVIL = "Est_Civil",
NIVEL_SOCIO = "Nivel_Socio"
)
# =========================
# 5. CATEGORÍAS DE REFERENCIA
# EXTRAÍDAS Y AJUSTADAS SEGÚN TU RECODIFICACIÓN
# =========================
referencias <- list(
SOLVENTAR = "Beca/exoneración total",
RESIDENCIA = "ASUNCIÓN",
CARRERA = "ESTADISTICA-PRES",
TIPO_INGRESO = "INGRESO",
RENDIMIENTO = "NO APROBADO",
MODALIDAD = "Presencial",
TRABAJA = "NO",
ESTUDIOS_PADRES = "hasta 13 años",
SEXO = "FEMENINO",
TIPO_COL = "Público",
INGRESO = "Hasta dos salarios mínimos",
EST_CIVIL = "Soltero/a",
NIVEL_SOCIO = levels(datos$NIVEL_SOCIO)[1]
)
# =========================
# 6. FUNCIÓN DE SIGNIFICANCIA
# =========================
etiqueta_significancia <- function(p) {
if (is.na(p)) {
return(NA_character_)
} else if (p < 0.01) {
return("*** (p < 0.01)")
} else if (p < 0.05) {
return("** (p < 0.05)")
} else if (p < 0.10) {
return("* (p < 0.10)")
} else {
return("No sig.")
}
}
# =========================
# 7. FUNCIÓN PARA AJUSTAR UN MODELO
# =========================
ajustar_logit_simple_cat_v2 <- function(var, indice_modelo, datos, refs, etiquetas_vars) {
datos_modelo <- datos %>%
select(ESTADO_BIN, all_of(var)) %>%
filter(!is.na(.data[[var]]))
datos_modelo[[var]] <- as.factor(datos_modelo[[var]])
# Fijar categoría de referencia
ref_actual <- refs[[var]]
niveles_actuales <- levels(datos_modelo[[var]])
if (!is.null(ref_actual) && ref_actual %in% niveles_actuales) {
datos_modelo[[var]] <- relevel(datos_modelo[[var]], ref = ref_actual)
} else {
warning(paste(
"No se encontró la referencia de", var,
"en los niveles de la base. Se usará la primera categoría disponible."
))
ref_actual <- niveles_actuales[1]
datos_modelo[[var]] <- relevel(datos_modelo[[var]], ref = ref_actual)
}
# Ajuste del modelo
formula_modelo <- as.formula(paste("ESTADO_BIN ~", var))
modelo <- glm(formula_modelo, data = datos_modelo, family = binomial(link = "logit"))
# Resumen del modelo
sm <- summary(modelo)$coefficients
sm <- as.data.frame(sm)
sm$term <- rownames(sm)
rownames(sm) <- NULL
# Excluir intercepto
sm <- sm %>%
filter(term != "(Intercept)")
# Intervalos de confianza Wald
sm <- sm %>%
mutate(
IC_beta_LI = Estimate - 1.96 * `Std. Error`,
IC_beta_LS = Estimate + 1.96 * `Std. Error`,
OR = exp(Estimate),
IC95_LI = exp(IC_beta_LI),
IC95_LS = exp(IC_beta_LS)
)
# Extraer nivel comparado
sm <- sm %>%
mutate(
Nivel = str_remove(term, paste0("^", var)),
`N°` = paste0("M", indice_modelo),
Variable = etiquetas_vars[[var]],
`Categoría de Referencia/Nivel` = ref_actual,
`Coefi-cientes` = Estimate,
`Exp (OR)` = OR,
`Estadístico_Z` = `z value`,
P_Valor = `Pr(>|z|)`,
Significancia = sapply(P_Valor, etiqueta_significancia)
) %>%
select(
`N°`,
Variable,
`Categoría de Referencia/Nivel`,
Nivel,
`Coefi-cientes`,
OR,
IC95_LI,
IC95_LS,
`Exp (OR)`,
`Estadístico_Z`,
P_Valor,
Significancia
)
sm
}
# =========================
# 8. AJUSTAR LOS 13 MODELOS
# =========================
tabla_4_10_v2 <- map2_dfr(
variables_categoricas,
seq_along(variables_categoricas),
~ ajustar_logit_simple_cat_v2(.x, .y, datos, referencias, etiquetas_variables)
)
# =========================
# 9. REDONDEAR
# =========================
tabla_4_10_v2 <- tabla_4_10_v2 %>%
mutate(
`Coefi-cientes` = round(`Coefi-cientes`, 4),
OR = round(OR, 4),
IC95_LI = round(IC95_LI, 4),
IC95_LS = round(IC95_LS, 4),
`Exp (OR)` = round(`Exp (OR)`, 4),
`Estadístico_Z` = round(`Estadístico_Z`, 4),
P_Valor = round(P_Valor, 4)
)
# =========================
# 10. FORMATO TIPO TESIS
# Mostrar una sola vez N°, Variable y Referencia por bloque
# =========================
tabla_4_10_v2_formato <- tabla_4_10_v2 %>%
group_by(`N°`, Variable) %>%
mutate(
`N°` = ifelse(row_number() == 1, `N°`, ""),
Variable = ifelse(row_number() == 1, Variable, ""),
`Categoría de Referencia/Nivel` = ifelse(
row_number() == 1,
`Categoría de Referencia/Nivel`,
""
)
) %>%
ungroup()
# =========================
# 11. RENOMBRAR COLUMNAS FINALES
# =========================
tabla_4_10_v2_formato <- tabla_4_10_v2_formato %>%
rename(
`IC.95 Inferior` = IC95_LI,
`IC.95 Superior` = IC95_LS
)
print(tabla_4_10_v2_formato)## # A tibble: 28 × 12
## `N°` Variable Categoría de Referencia/Ni…¹ Nivel `Coefi-cientes` OR
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 "M1" "Solventar" "Beca/exoneración total" Beca… -0.182 0.833
## 2 "" "" "" Trab… 0.500 1.65
## 3 "" "" "" Ayud… -0.240 0.787
## 4 "M2" "Residencia" "ASUNCIÓN" CENT… -0.345 0.708
## 5 "" "" "" REST… 0.398 1.49
## 6 "M3" "Carrera" "ESTADISTICA-PRES" ESTA… 0.718 2.05
## 7 "" "" "" MATE… -0.507 0.602
## 8 "" "" "" EDUC… -1.60 0.202
## 9 "" "" "" EDUC… 0.495 1.64
## 10 "M4" "Tipo_Ingreso" "INGRESO" ADMI… -0.505 0.603
## # ℹ 18 more rows
## # ℹ abbreviated name: ¹`Categoría de Referencia/Nivel`
## # ℹ 6 more variables: `IC.95 Inferior` <dbl>, `IC.95 Superior` <dbl>,
## # `Exp (OR)` <dbl>, Estadístico_Z <dbl>, P_Valor <dbl>, Significancia <chr># =========================
# 12. EXPORTAR A EXCEL
# =========================
write.xlsx(
tabla_4_10_v2_formato,
file = "Tabla_4_10_logistica_simple_categoricas_v2.xlsx",
rowNames = FALSE
)
# =========================
# 13. TABLA BONITA PARA WORD / REPORTE
# =========================
ft_tabla_4_10 <- flextable(tabla_4_10_v2_formato)
ft_tabla_4_10 <- ft_tabla_4_10 |>
theme_booktabs() |>
autofit() |>
align(align = "center", part = "all") |>
valign(valign = "center", part = "all") |>
fontsize(size = 9, part = "all") |>
bold(part = "header") |>
set_caption("Tabla 4.10: Modelos de regresión con una sola variable explicativa con sus respectivos resúmenes")
ft_tabla_4_10N° | Variable | Categoría de Referencia/Nivel | Nivel | Coefi-cientes | OR | IC.95 Inferior | IC.95 Superior | Exp (OR) | Estadístico_Z | P_Valor | Significancia |
|---|---|---|---|---|---|---|---|---|---|---|---|
M1 | Solventar | Beca/exoneración total | Beca/exoneración parcial | -0.1823 | 0.8333 | 0.2984 | 2.3274 | 0.8333 | -0.3479 | 0.7279 | No sig. |
Trabajo Personal | 0.4998 | 1.6485 | 0.7219 | 3.7642 | 1.6485 | 1.1865 | 0.2354 | No sig. | |||
Ayuda Familiar | -0.2395 | 0.7870 | 0.3432 | 1.8046 | 0.7870 | -0.5656 | 0.5716 | No sig. | |||
M2 | Residencia | ASUNCIÓN | CENTRAL | -0.3452 | 0.7080 | 0.4618 | 1.0856 | 0.7080 | -1.5834 | 0.1133 | No sig. |
RESTO DEL PAÍS | 0.3984 | 1.4895 | 0.7976 | 2.7814 | 1.4895 | 1.2504 | 0.2111 | No sig. | |||
M3 | Carrera | ESTADISTICA-PRES | ESTADISTICA-SEMI | 0.7185 | 2.0513 | 1.0127 | 4.1551 | 2.0513 | 1.9950 | 0.0460 | ** (p < 0.05) |
MATEMATICA | -0.5067 | 0.6025 | 0.3526 | 1.0295 | 0.6025 | -1.8537 | 0.0638 | * (p < 0.10) | |||
EDUCACION MATEMATICA-PRES | -1.5994 | 0.2020 | 0.0988 | 0.4130 | 0.2020 | -4.3843 | 0.0000 | *** (p < 0.01) | |||
EDUCACION MATEMATICA-SEMI | 0.4953 | 1.6410 | 0.9267 | 2.9058 | 1.6410 | 1.6990 | 0.0893 | * (p < 0.10) | |||
M4 | Tipo_Ingreso | INGRESO | ADMISION DIRECTA | -0.5053 | 0.6033 | 0.2673 | 1.3615 | 0.6033 | -1.2169 | 0.2236 | No sig. |
M5 | Rendimiento | NO APROBADO | APROBADO | -2.8607 | 0.0572 | 0.0298 | 0.1098 | 0.0572 | -8.6092 | 0.0000 | *** (p < 0.01) |
M6 | Modalidad | Presencial | Semipresencial | 1.0450 | 2.8434 | 1.8629 | 4.3400 | 2.8434 | 4.8434 | 0.0000 | *** (p < 0.01) |
M7 | Trabaja | NO | SI | 0.8681 | 2.3825 | 1.6128 | 3.5195 | 2.3825 | 4.3611 | 0.0000 | *** (p < 0.01) |
M8 | Estudios_Padres | hasta 13 años | 14-23 años | 0.0544 | 1.0559 | 0.6199 | 1.7986 | 1.0559 | 0.2002 | 0.8413 | No sig. |
24-29 años | -0.0146 | 0.9855 | 0.5576 | 1.7419 | 0.9855 | -0.0502 | 0.9599 | No sig. | |||
30-34 años | 0.1715 | 1.1871 | 0.6229 | 2.2621 | 1.1871 | 0.5213 | 0.6021 | No sig. | |||
Más de 34 años | -0.3713 | 0.6899 | 0.3164 | 1.5042 | 0.6899 | -0.9335 | 0.3506 | No sig. | |||
M9 | Sexo | FEMENINO | MASCULINO | 0.1531 | 1.1654 | 0.7894 | 1.7207 | 1.1654 | 0.7702 | 0.4412 | No sig. |
M10 | Tipo_Col | Público | Subvencionado | -0.2752 | 0.7594 | 0.3823 | 1.5087 | 0.7594 | -0.7857 | 0.4320 | No sig. |
Privado | -0.0143 | 0.9858 | 0.6297 | 1.5434 | 0.9858 | -0.0624 | 0.9502 | No sig. | |||
M11 | Ingreso | Hasta dos salarios mínimos | Más de dos y hasta cinco salarios mínimos | 0.0958 | 1.1006 | 0.7156 | 1.6925 | 1.1006 | 0.4363 | 0.6626 | No sig. |
Más de cinco y hasta diez salarios mínimos | 0.3018 | 1.3523 | 0.5788 | 3.1596 | 1.3523 | 0.6970 | 0.4858 | No sig. | |||
Más de diez y hasta quince salarios mínimos | 1.2826 | 3.6061 | 0.4377 | 29.7108 | 3.6061 | 1.1921 | 0.2332 | No sig. | |||
Más de quince salarios mínimos | -0.6633 | 0.5152 | 0.0715 | 3.7107 | 0.5152 | -0.6584 | 0.5103 | No sig. | |||
M12 | Est_Civil | Soltero/a | Casado/a | 0.3768 | 1.4577 | 0.8481 | 2.5053 | 1.4577 | 1.3638 | 0.1726 | No sig. |
Divorciado/a | 0.8650 | 2.3750 | 0.5058 | 11.1510 | 2.3750 | 1.0963 | 0.2730 | No sig. | |||
Otro | 0.4595 | 1.5833 | 0.1631 | 15.3708 | 1.5833 | 0.3963 | 0.6919 | No sig. | |||
M13 | Nivel_Socio | BAJO | MEDIO | 0.5208 | 1.6835 | 0.8799 | 3.2208 | 1.6835 | 1.5735 | 0.1156 | No sig. |
# =========================
# 14. OPCIONAL: AGREGAR EDAD AL FINAL
# =========================
# Si quieres anexar Edad ya calculada manualmente:
#
# fila_edad <- data.frame(
# `N°` = "M14",
# Variable = "Edad",
# `Categoría de Referencia/Nivel` = "",
# Nivel = "",
# `Coefi-cientes` = -0.0268,
# OR = 0.9735,
# `IC.95 Inferior` = 0.9526,
# `IC.95 Superior` = 0.9949,
# `Exp (OR)` = 0.9736,
# `Estadístico_Z` = -2.4231,
# P_Valor = 0.0154,
# Significancia = "** (p < 0.05)"
# )
#
# tabla_4_10_v2_completa <- bind_rows(tabla_4_10_v2_formato, fila_edad)
# print(tabla_4_10_v2_completa)############################################################
# TABLA COMPARATIVA ENTRE train_data Y test_data
# Estructura agrupada con las 14 variables de la tesis
# Versión corregida: sin error en bind_rows()
############################################################
# =========================
# 0. PAQUETES
# =========================
paquetes <- c("dplyr", "openxlsx")
instalar <- paquetes[!paquetes %in% installed.packages()[, "Package"]]
if (length(instalar) > 0) install.packages(instalar)
library(dplyr)
library(openxlsx)
# =========================
# 1. VALIDACIONES
# =========================
if (!exists("train_data")) stop("No existe el objeto 'train_data'.")
if (!exists("test_data")) stop("No existe el objeto 'test_data'.")
vars_necesarias <- c(
"ESTADO_ACADEMICO",
"SEXO", "RESIDENCIA", "EST_CIVIL", "EDAD",
"CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO",
"TIPO_COL", "TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES",
"INGRESO", "NIVEL_SOCIO"
)
faltan_train <- setdiff(vars_necesarias, names(train_data))
faltan_test <- setdiff(vars_necesarias, names(test_data))
if (length(faltan_train) > 0) {
stop("Faltan en train_data las variables: ", paste(faltan_train, collapse = ", "))
}
if (length(faltan_test) > 0) {
stop("Faltan en test_data las variables: ", paste(faltan_test, collapse = ", "))
}
# =========================
# 2. ASEGURAR CODIFICACIÓN
# =========================
train_data$ESTADO_ACADEMICO <- as.factor(train_data$ESTADO_ACADEMICO)
test_data$ESTADO_ACADEMICO <- as.factor(test_data$ESTADO_ACADEMICO)
# =========================
# 3. FUNCIONES AUXILIARES
# =========================
fmt_n_pct <- function(n, total) {
paste0(n, " (", round(100 * n / total, 2), "%)")
}
fmt_pct <- function(x) {
paste0(round(x, 2), "%")
}
obs_prop <- function(p1, p2) {
dif <- abs(p1 - p2)
if (is.na(dif)) {
""
} else if (dif < 1) {
"Distribución prácticamente idéntica"
} else if (dif < 3) {
"Proporción consistente"
} else if (dif < 5) {
"Diferencia leve"
} else {
"Diferencia moderada"
}
}
obs_num <- function(x1, x2, umbral1 = 1, umbral2 = 2) {
dif <- abs(x1 - x2)
if (is.na(dif)) {
""
} else if (dif < umbral1) {
"Valores muy similares"
} else if (dif < umbral2) {
"Valores similares"
} else {
"Diferencia moderada"
}
}
crear_seccion <- function(nombre) {
data.frame(
Característica = as.character(nombre),
`Base de Entrenamiento` = as.character(""),
`Base de Prueba` = as.character(""),
Observación = as.character(""),
Tipo = as.character("SECCION"),
check.names = FALSE,
stringsAsFactors = FALSE
)
}
fila_texto <- function(caracteristica, entrena, prueba, obs = "") {
data.frame(
Característica = as.character(caracteristica),
`Base de Entrenamiento` = as.character(entrena),
`Base de Prueba` = as.character(prueba),
Observación = as.character(obs),
Tipo = as.character("DATO"),
check.names = FALSE,
stringsAsFactors = FALSE
)
}
resumen_categoria <- function(df, var, categoria) {
sum(df[[var]] == categoria, na.rm = TRUE)
}
fila_categoria <- function(var, categoria, etiqueta) {
n_train <- resumen_categoria(train_data, var, categoria)
n_test <- resumen_categoria(test_data, var, categoria)
total_train <- sum(!is.na(train_data[[var]]))
total_test <- sum(!is.na(test_data[[var]]))
p_train <- 100 * n_train / total_train
p_test <- 100 * n_test / total_test
fila_texto(
etiqueta,
fmt_n_pct(n_train, total_train),
fmt_n_pct(n_test, total_test),
obs_prop(p_train, p_test)
)
}
# =========================
# 4. RESUMEN GENERAL
# =========================
n_train <- nrow(train_data)
n_test <- nrow(test_data)
n_total <- n_train + n_test
prop_train <- 100 * n_train / n_total
prop_test <- 100 * n_test / n_total
prop_des_train <- mean(train_data$ESTADO_ACADEMICO == "Desertor", na.rm = TRUE) * 100
prop_des_test <- mean(test_data$ESTADO_ACADEMICO == "Desertor", na.rm = TRUE) * 100
prop_no_train <- mean(train_data$ESTADO_ACADEMICO == "No_Desertor", na.rm = TRUE) * 100
prop_no_test <- mean(test_data$ESTADO_ACADEMICO == "No_Desertor", na.rm = TRUE) * 100
tabla <- bind_rows(
crear_seccion("Resumen general"),
fila_texto(
"Tamaño de muestra",
paste0(round(prop_train, 0), "% (", n_train, " obs.)"),
paste0(round(prop_test, 0), "% (", n_test, " obs.)"),
"Proporción definida en la partición"
),
fila_texto(
"Desertores",
fmt_pct(prop_des_train),
fmt_pct(prop_des_test),
obs_prop(prop_des_train, prop_des_test)
),
fila_texto(
"No desertores",
fmt_pct(prop_no_train),
fmt_pct(prop_no_test),
obs_prop(prop_no_train, prop_no_test)
)
)
# =========================
# 5. VARIABLES DEMOGRÁFICAS
# =========================
edad_media_train <- mean(train_data$EDAD, na.rm = TRUE)
edad_media_test <- mean(test_data$EDAD, na.rm = TRUE)
edad_mediana_train <- median(train_data$EDAD, na.rm = TRUE)
edad_mediana_test <- median(test_data$EDAD, na.rm = TRUE)
edad_sd_train <- sd(train_data$EDAD, na.rm = TRUE)
edad_sd_test <- sd(test_data$EDAD, na.rm = TRUE)
edad_min_train <- min(train_data$EDAD, na.rm = TRUE)
edad_min_test <- min(test_data$EDAD, na.rm = TRUE)
edad_max_train <- max(train_data$EDAD, na.rm = TRUE)
edad_max_test <- max(test_data$EDAD, na.rm = TRUE)
tabla <- bind_rows(
tabla,
crear_seccion("Variables demográficas"),
fila_texto(
"Edad (media)",
round(edad_media_train, 2),
round(edad_media_test, 2),
obs_num(edad_media_train, edad_media_test, 1, 2)
),
fila_texto(
"Edad (mediana)",
round(edad_mediana_train, 2),
round(edad_mediana_test, 2),
ifelse(edad_mediana_train == edad_mediana_test,
"Consistencia en tendencia central",
"Diferencia en mediana")
),
fila_texto(
"Edad (DE)",
round(edad_sd_train, 2),
round(edad_sd_test, 2),
obs_num(edad_sd_train, edad_sd_test, 1, 2)
),
fila_texto("Edad (mínimo)", round(edad_min_train, 2), round(edad_min_test, 2), ""),
fila_texto("Edad (máximo)", round(edad_max_train, 2), round(edad_max_test, 2), ""),
fila_categoria("SEXO", "FEMENINO", "Sexo (FEMENINO)"),
fila_categoria("SEXO", "MASCULINO", "Sexo (MASCULINO)"),
fila_categoria("RESIDENCIA", "ASUNCIÓN", "Residencia (ASUNCIÓN)"),
fila_categoria("RESIDENCIA", "CENTRAL", "Residencia (CENTRAL)"),
fila_categoria("RESIDENCIA", "RESTO DEL PAÍS", "Residencia (RESTO DEL PAÍS)"),
fila_categoria("EST_CIVIL", "Soltero/a", "Estado civil (Soltero/a)"),
fila_categoria("EST_CIVIL", "Casado/a", "Estado civil (Casado/a)"),
fila_categoria("EST_CIVIL", "Divorciado/a", "Estado civil (Divorciado/a)"),
fila_categoria("EST_CIVIL", "Otro", "Estado civil (Otro)")
)
# =========================
# 6. VARIABLES INTERNAS AL MEDIO EDUCATIVO
# =========================
tabla <- bind_rows(
tabla,
crear_seccion("Variables internas al medio educativo"),
fila_categoria("CARRERA", "ESTADISTICA-PRES", "Carrera (ESTADISTICA-PRES)"),
fila_categoria("CARRERA", "ESTADISTICA-SEMI", "Carrera (ESTADISTICA-SEMI)"),
fila_categoria("CARRERA", "MATEMATICA", "Carrera (MATEMATICA)"),
fila_categoria("CARRERA", "EDUCACION MATEMATICA-PRES", "Carrera (EDUCACION MATEMATICA-PRES)"),
fila_categoria("CARRERA", "EDUCACION MATEMATICA-SEMI", "Carrera (EDUCACION MATEMATICA-SEMI)"),
fila_categoria("TIPO_INGRESO", "INGRESO", "Tipo de ingreso (INGRESO)"),
fila_categoria("TIPO_INGRESO", "TRASLADO", "Tipo de ingreso (TRASLADO)"),
fila_categoria("TIPO_INGRESO", "ADMISION DIRECTA", "Tipo de ingreso (ADMISION DIRECTA)"),
fila_categoria("MODALIDAD", "Presencial", "Modalidad (Presencial)"),
fila_categoria("MODALIDAD", "Semipresencial", "Modalidad (Semipresencial)"),
fila_categoria("RENDIMIENTO", "APROBADO", "Rendimiento (APROBADO)"),
fila_categoria("RENDIMIENTO", "NO APROBADO", "Rendimiento (NO APROBADO)")
)
# =========================
# 7. VARIABLES EXTERNAS AL MEDIO EDUCATIVO
# =========================
tabla <- bind_rows(
tabla,
crear_seccion("Variables externas al medio educativo"),
fila_categoria("TIPO_COL", "Público", "Tipo de colegio (Público)"),
fila_categoria("TIPO_COL", "Subvencionado", "Tipo de colegio (Subvencionado)"),
fila_categoria("TIPO_COL", "Privado", "Tipo de colegio (Privado)"),
fila_categoria("TRABAJA", "SI", "Trabaja (SI)"),
fila_categoria("TRABAJA", "NO", "Trabaja (NO)"),
fila_categoria("SOLVENTAR", "Beca/exoneración total", "Solventar (Beca/exoneración total)"),
fila_categoria("SOLVENTAR", "Beca/exoneración parcial", "Solventar (Beca/exoneración parcial)"),
fila_categoria("SOLVENTAR", "Trabajo Personal", "Solventar (Trabajo Personal)"),
fila_categoria("SOLVENTAR", "Ayuda Familiar", "Solventar (Ayuda Familiar)"),
fila_categoria("ESTUDIOS_PADRES", "hasta 13 años", "Estudios de los padres (hasta 13 años)"),
fila_categoria("ESTUDIOS_PADRES", "14-23 años", "Estudios de los padres (14-23 años)"),
fila_categoria("ESTUDIOS_PADRES", "24-29 años", "Estudios de los padres (24-29 años)"),
fila_categoria("ESTUDIOS_PADRES", "30-34 años", "Estudios de los padres (30-34 años)"),
fila_categoria("ESTUDIOS_PADRES", "Más de 34 años", "Estudios de los padres (Más de 34 años)"),
fila_categoria("INGRESO", "Hasta dos salarios mínimos", "Ingreso familiar (Hasta dos salarios mínimos)"),
fila_categoria("INGRESO", "Más de dos y hasta cinco salarios mínimos", "Ingreso familiar (Más de dos y hasta cinco salarios mínimos)"),
fila_categoria("INGRESO", "Más de cinco y hasta diez salarios mínimos", "Ingreso familiar (Más de cinco y hasta diez salarios mínimos)"),
fila_categoria("INGRESO", "Más de diez y hasta quince salarios mínimos", "Ingreso familiar (Más de diez y hasta quince salarios mínimos)"),
fila_categoria("INGRESO", "Más de quince salarios mínimos", "Ingreso familiar (Más de quince salarios mínimos)")
)
# NIVEL_SOCIO dinámico según niveles existentes
niveles_nivel_socio <- union(
levels(as.factor(train_data$NIVEL_SOCIO)),
levels(as.factor(test_data$NIVEL_SOCIO))
)
for (niv in niveles_nivel_socio) {
tabla <- bind_rows(
tabla,
fila_categoria("NIVEL_SOCIO", niv, paste0("Nivel socioeconómico (", niv, ")"))
)
}
# =========================
# 8. LIMPIEZA FINAL
# =========================
tabla_exportar <- tabla %>%
select(-Tipo)
print(tabla_exportar)## Característica
## 1 Resumen general
## 2 Tamaño de muestra
## 3 Desertores
## 4 No desertores
## 5 Variables demográficas
## 6 Edad (media)
## 7 Edad (mediana)
## 8 Edad (DE)
## 9 Edad (mínimo)
## 10 Edad (máximo)
## 11 Sexo (FEMENINO)
## 12 Sexo (MASCULINO)
## 13 Residencia (ASUNCIÓN)
## 14 Residencia (CENTRAL)
## 15 Residencia (RESTO DEL PAÍS)
## 16 Estado civil (Soltero/a)
## 17 Estado civil (Casado/a)
## 18 Estado civil (Divorciado/a)
## 19 Estado civil (Otro)
## 20 Variables internas al medio educativo
## 21 Carrera (ESTADISTICA-PRES)
## 22 Carrera (ESTADISTICA-SEMI)
## 23 Carrera (MATEMATICA)
## 24 Carrera (EDUCACION MATEMATICA-PRES)
## 25 Carrera (EDUCACION MATEMATICA-SEMI)
## 26 Tipo de ingreso (INGRESO)
## 27 Tipo de ingreso (TRASLADO)
## 28 Tipo de ingreso (ADMISION DIRECTA)
## 29 Modalidad (Presencial)
## 30 Modalidad (Semipresencial)
## 31 Rendimiento (APROBADO)
## 32 Rendimiento (NO APROBADO)
## 33 Variables externas al medio educativo
## 34 Tipo de colegio (Público)
## 35 Tipo de colegio (Subvencionado)
## 36 Tipo de colegio (Privado)
## 37 Trabaja (SI)
## 38 Trabaja (NO)
## 39 Solventar (Beca/exoneración total)
## 40 Solventar (Beca/exoneración parcial)
## 41 Solventar (Trabajo Personal)
## 42 Solventar (Ayuda Familiar)
## 43 Estudios de los padres (hasta 13 años)
## 44 Estudios de los padres (14-23 años)
## 45 Estudios de los padres (24-29 años)
## 46 Estudios de los padres (30-34 años)
## 47 Estudios de los padres (Más de 34 años)
## 48 Ingreso familiar (Hasta dos salarios mínimos)
## 49 Ingreso familiar (Más de dos y hasta cinco salarios mínimos)
## 50 Ingreso familiar (Más de cinco y hasta diez salarios mínimos)
## 51 Ingreso familiar (Más de diez y hasta quince salarios mínimos)
## 52 Ingreso familiar (Más de quince salarios mínimos)
## 53 Nivel socioeconómico (BAJO)
## 54 Nivel socioeconómico (MEDIO)
## Base de Entrenamiento Base de Prueba Observación
## 1
## 2 70% (479 obs.) 30% (203 obs.) Proporción definida en la partición
## 3 67.22% 67.49% Distribución prácticamente idéntica
## 4 32.78% 32.51% Distribución prácticamente idéntica
## 5
## 6 26.84 27.74 Valores muy similares
## 7 23 24 Diferencia en mediana
## 8 9.03 9.47 Valores muy similares
## 9 18 18
## 10 70 61
## 11 281 (58.66%) 114 (56.16%) Proporción consistente
## 12 198 (41.34%) 89 (43.84%) Proporción consistente
## 13 159 (33.19%) 70 (34.48%) Proporción consistente
## 14 240 (50.1%) 96 (47.29%) Proporción consistente
## 15 80 (16.7%) 37 (18.23%) Proporción consistente
## 16 385 (80.38%) 159 (78.33%) Proporción consistente
## 17 79 (16.49%) 38 (18.72%) Proporción consistente
## 18 11 (2.3%) 3 (1.48%) Distribución prácticamente idéntica
## 19 4 (0.84%) 3 (1.48%) Distribución prácticamente idéntica
## 20
## 21 130 (27.14%) 52 (25.62%) Proporción consistente
## 22 73 (15.24%) 27 (13.3%) Proporción consistente
## 23 106 (22.13%) 44 (21.67%) Distribución prácticamente idéntica
## 24 48 (10.02%) 22 (10.84%) Distribución prácticamente idéntica
## 25 122 (25.47%) 58 (28.57%) Diferencia leve
## 26 454 (94.78%) 188 (92.61%) Proporción consistente
## 27 0 (0%) 0 (0%) Distribución prácticamente idéntica
## 28 25 (5.22%) 15 (7.39%) Proporción consistente
## 29 284 (59.29%) 118 (58.13%) Proporción consistente
## 30 195 (40.71%) 85 (41.87%) Proporción consistente
## 31 285 (59.5%) 124 (61.08%) Proporción consistente
## 32 194 (40.5%) 79 (38.92%) Proporción consistente
## 33
## 34 320 (66.81%) 145 (71.43%) Diferencia leve
## 35 39 (8.14%) 10 (4.93%) Diferencia leve
## 36 120 (25.05%) 48 (23.65%) Proporción consistente
## 37 261 (54.49%) 111 (54.68%) Distribución prácticamente idéntica
## 38 218 (45.51%) 92 (45.32%) Distribución prácticamente idéntica
## 39 28 (5.85%) 17 (8.37%) Proporción consistente
## 40 35 (7.31%) 14 (6.9%) Distribución prácticamente idéntica
## 41 242 (50.52%) 104 (51.23%) Distribución prácticamente idéntica
## 42 174 (36.33%) 68 (33.5%) Proporción consistente
## 43 103 (21.5%) 54 (26.6%) Diferencia moderada
## 44 154 (32.15%) 61 (30.05%) Proporción consistente
## 45 111 (23.17%) 53 (26.11%) Proporción consistente
## 46 75 (15.66%) 26 (12.81%) Proporción consistente
## 47 36 (7.52%) 9 (4.43%) Diferencia leve
## 48 300 (62.63%) 130 (64.04%) Proporción consistente
## 49 138 (28.81%) 52 (25.62%) Diferencia leve
## 50 29 (6.05%) 16 (7.88%) Proporción consistente
## 51 8 (1.67%) 4 (1.97%) Distribución prácticamente idéntica
## 52 4 (0.84%) 1 (0.49%) Distribución prácticamente idéntica
## 53 41 (8.56%) 24 (11.82%) Diferencia leve
## 54 438 (91.44%) 179 (88.18%) Diferencia leve# =========================
# 9. EXPORTAR A EXCEL
# =========================
wb <- createWorkbook()
addWorksheet(wb, "Tabla comparativa")
writeData(wb, "Tabla comparativa", tabla_exportar, startRow = 1, startCol = 1)
# Estilos
estilo_header <- createStyle(
textDecoration = "bold",
fgFill = "#5B9BD5",
fontColour = "#FFFFFF",
halign = "center",
valign = "center",
wrapText = TRUE
)
estilo_seccion <- createStyle(
textDecoration = "bold",
fgFill = "#D9EAF7",
wrapText = TRUE
)
estilo_general <- createStyle(
valign = "top",
wrapText = TRUE
)
# Aplicar estilos
addStyle(wb, "Tabla comparativa", estilo_header, rows = 1, cols = 1:4, gridExpand = TRUE)
filas_seccion <- which(tabla$Tipo == "SECCION") + 1
if (length(filas_seccion) > 0) {
for (f in filas_seccion) {
addStyle(wb, "Tabla comparativa", estilo_seccion, rows = f, cols = 1:4, gridExpand = TRUE)
}
}
addStyle(
wb, "Tabla comparativa", estilo_general,
rows = 2:(nrow(tabla_exportar) + 1), cols = 1:4, gridExpand = TRUE, stack = TRUE
)
# Anchos de columna
setColWidths(wb, "Tabla comparativa", cols = 1, widths = 48)
setColWidths(wb, "Tabla comparativa", cols = 2:3, widths = 24)
setColWidths(wb, "Tabla comparativa", cols = 4, widths = 32)
freezePane(wb, "Tabla comparativa", firstRow = TRUE)
saveWorkbook(
wb,
"Tabla_comparacion_train_test_14_variables.xlsx",
overwrite = TRUE
)
cat("\nArchivo generado: Tabla_comparacion_train_test_14_variables.xlsx\n")##
## Archivo generado: Tabla_comparacion_train_test_14_variables.xlsx############################################################
# TABLA 4.19
# RESULTADOS DE LA VALIDACIÓN CRUZADA
# ROC, SENSIBILIDAD Y ESPECIFICIDAD
############################################################
# Se asume que ya existe:
# resultados_cv <- resamples(list(
# Logit = modelo_logit_cv,
# Tree = modelo_tree_cv,
# RF = modelo_rf_cv,
# SVM = modelo_svm_cv,
# KNN = modelo_knn_cv,
# NN = modelo_nn_cv,
# NB = modelo_nb_cv
# ))
library(dplyr)
library(openxlsx)
############################################################
# 1. RESUMEN DE VALIDACIÓN CRUZADA
############################################################
resumen_cv <- summary(resultados_cv)
############################################################
# 2. FUNCIÓN PARA ARMAR CADA TABLA DE MÉTRICA
############################################################
armar_tabla_metrica <- function(matriz_metric) {
tabla <- as.data.frame(matriz_metric, stringsAsFactors = FALSE)
# Eliminar columna de NA's si existe
if ("NA's" %in% colnames(tabla)) {
tabla <- tabla[, !colnames(tabla) %in% "NA's", drop = FALSE]
}
tabla$Modelos <- rownames(tabla)
rownames(tabla) <- NULL
# Reordenar columnas
tabla <- tabla[, c("Modelos", "Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max.")]
# Redondear columnas numéricas
tabla <- tabla %>%
mutate(across(where(is.numeric), ~ round(.x, 7))) %>%
as.data.frame()
tabla
}
############################################################
# 3. TABLAS POR MÉTRICA
############################################################
tabla_roc_419 <- armar_tabla_metrica(resumen_cv$statistics$ROC)
tabla_sens_419 <- armar_tabla_metrica(resumen_cv$statistics$Sens)
tabla_spec_419 <- armar_tabla_metrica(resumen_cv$statistics$Spec)
cat("\n=============================\n")##
## =============================cat("TABLA 4.19 - ROC\n")## TABLA 4.19 - ROCcat("=============================\n")## =============================print(tabla_roc_419)## Modelos Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 Logit 0.6541667 0.7437500 0.7836174 0.7803046 0.8199219 0.8996212
## 2 Tree 0.6845703 0.7714962 0.7973633 0.7940493 0.8305827 0.8750000
## 3 RF 0.6458333 0.7341856 0.7910156 0.7854376 0.8286458 0.9091797
## 4 SVM 0.6416667 0.7183268 0.7613636 0.7602734 0.8041992 0.8867188
## 5 KNN 0.6318359 0.6966072 0.7434186 0.7483594 0.7955211 0.9062500
## 6 NN 0.6208333 0.7228634 0.7753906 0.7677115 0.8043176 0.8916667
## 7 NB 0.6041667 0.7292480 0.7802734 0.7714621 0.8144383 0.9104167cat("\n=============================\n")##
## =============================cat("TABLA 4.19 - SENSIBILIDAD\n")## TABLA 4.19 - SENSIBILIDADcat("=============================\n")## =============================print(tabla_sens_419)## Modelos Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 Logit 0.68750 0.7500000 0.8125000 0.8071023 0.8473011 0.9696970
## 2 Tree 0.65625 0.7500000 0.7812500 0.7886995 0.8167614 0.9375000
## 3 RF 0.90625 0.9375000 0.9687500 0.9616477 0.9924242 1.0000000
## 4 SVM 0.68750 0.8437500 0.8787879 0.8777462 0.9375000 1.0000000
## 5 KNN 0.81250 0.9062500 0.9375000 0.9304924 0.9614110 1.0000000
## 6 NN 0.65625 0.7500000 0.8125000 0.7967487 0.8437500 0.9090909
## 7 NB 0.37500 0.6328125 0.7187500 0.6799874 0.7812500 0.8181818cat("\n=============================\n")##
## =============================cat("TABLA 4.19 - ESPECIFICIDAD\n")## TABLA 4.19 - ESPECIFICIDADcat("=============================\n")## =============================print(tabla_spec_419)## Modelos Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 Logit 0.3333333 0.5000000 0.5625000 0.5548611 0.62500 0.7500
## 2 Tree 0.3125000 0.4666667 0.6125000 0.5850000 0.68750 0.8750
## 3 RF 0.0625000 0.1468750 0.1875000 0.2031944 0.25000 0.3750
## 4 SVM 0.0000000 0.2000000 0.3229167 0.3322222 0.43750 0.6875
## 5 KNN 0.0000000 0.1468750 0.2000000 0.2259722 0.31250 0.4375
## 6 NN 0.2500000 0.3750000 0.5000000 0.4906944 0.61875 0.7500
## 7 NB 0.4375000 0.6250000 0.7500000 0.7298611 0.81250 1.0000############################################################
# 4. TABLA CONSOLIDADA EN FORMATO TIPO TESIS
############################################################
fila_titulo <- function(texto) {
data.frame(
Modelos = texto,
`Min.` = "",
`1st Qu.` = "",
Median = "",
Mean = "",
`3rd Qu.` = "",
`Max.` = "",
stringsAsFactors = FALSE,
check.names = FALSE
)
}
fila_blanca <- data.frame(
Modelos = "",
`Min.` = "",
`1st Qu.` = "",
Median = "",
Mean = "",
`3rd Qu.` = "",
`Max.` = "",
stringsAsFactors = FALSE,
check.names = FALSE
)
# Convertir a character para poder unir títulos y números sin errores
tabla_roc_chr <- tabla_roc_419 %>% mutate(across(everything(), as.character))
tabla_sens_chr <- tabla_sens_419 %>% mutate(across(everything(), as.character))
tabla_spec_chr <- tabla_spec_419 %>% mutate(across(everything(), as.character))
tabla_419_final <- bind_rows(
fila_titulo("ROC (AUC)"),
tabla_roc_chr,
fila_blanca,
fila_titulo("Sens (Sensibilidad)"),
tabla_sens_chr,
fila_blanca,
fila_titulo("Spec (Especificidad)"),
tabla_spec_chr,
fila_blanca,
data.frame(
Modelos = "Models: Logit, Tree, RF, SVM, KNN, NN, NB",
`Min.` = "",
`1st Qu.` = "",
Median = "",
Mean = "",
`3rd Qu.` = "",
`Max.` = "",
stringsAsFactors = FALSE,
check.names = FALSE
),
data.frame(
Modelos = paste0("Number of resamples: ", nrow(resultados_cv$values)),
`Min.` = "",
`1st Qu.` = "",
Median = "",
Mean = "",
`3rd Qu.` = "",
`Max.` = "",
stringsAsFactors = FALSE,
check.names = FALSE
)
)
cat("\n=============================\n")##
## =============================cat("TABLA 4.19 CONSOLIDADA\n")## TABLA 4.19 CONSOLIDADAcat("=============================\n")## =============================print(tabla_419_final)## Modelos Min. 1st Qu. Median
## 1 ROC (AUC)
## 2 Logit 0.6541667 0.74375 0.7836174
## 3 Tree 0.6845703 0.7714962 0.7973633
## 4 RF 0.6458333 0.7341856 0.7910156
## 5 SVM 0.6416667 0.7183268 0.7613636
## 6 KNN 0.6318359 0.6966072 0.7434186
## 7 NN 0.6208333 0.7228634 0.7753906
## 8 NB 0.6041667 0.729248 0.7802734
## 9
## 10 Sens (Sensibilidad)
## 11 Logit 0.6875 0.75 0.8125
## 12 Tree 0.65625 0.75 0.78125
## 13 RF 0.90625 0.9375 0.96875
## 14 SVM 0.6875 0.84375 0.8787879
## 15 KNN 0.8125 0.90625 0.9375
## 16 NN 0.65625 0.75 0.8125
## 17 NB 0.375 0.6328125 0.71875
## 18
## 19 Spec (Especificidad)
## 20 Logit 0.3333333 0.5 0.5625
## 21 Tree 0.3125 0.4666667 0.6125
## 22 RF 0.0625 0.146875 0.1875
## 23 SVM 0 0.2 0.3229167
## 24 KNN 0 0.146875 0.2
## 25 NN 0.25 0.375 0.5
## 26 NB 0.4375 0.625 0.75
## 27
## 28 Models: Logit, Tree, RF, SVM, KNN, NN, NB
## 29 Number of resamples: 30
## Mean 3rd Qu. Max.
## 1
## 2 0.7803046 0.8199219 0.8996212
## 3 0.7940493 0.8305827 0.875
## 4 0.7854376 0.8286458 0.9091797
## 5 0.7602734 0.8041992 0.8867188
## 6 0.7483594 0.7955211 0.90625
## 7 0.7677115 0.8043176 0.8916667
## 8 0.7714621 0.8144383 0.9104167
## 9
## 10
## 11 0.8071023 0.8473011 0.969697
## 12 0.7886995 0.8167614 0.9375
## 13 0.9616477 0.9924242 1
## 14 0.8777462 0.9375 1
## 15 0.9304924 0.961411 1
## 16 0.7967487 0.84375 0.9090909
## 17 0.6799874 0.78125 0.8181818
## 18
## 19
## 20 0.5548611 0.625 0.75
## 21 0.585 0.6875 0.875
## 22 0.2031944 0.25 0.375
## 23 0.3322222 0.4375 0.6875
## 24 0.2259722 0.3125 0.4375
## 25 0.4906944 0.61875 0.75
## 26 0.7298611 0.8125 1
## 27
## 28
## 29############################################################
# 5. EXPORTACIÓN
############################################################
dir.create("salidas_tesis", showWarnings = FALSE)
# CSV separados
write.csv(tabla_roc_419,
"salidas_tesis/tabla_4_19_roc_cv.csv",
row.names = FALSE)
write.csv(tabla_sens_419,
"salidas_tesis/tabla_4_19_sens_cv.csv",
row.names = FALSE)
write.csv(tabla_spec_419,
"salidas_tesis/tabla_4_19_spec_cv.csv",
row.names = FALSE)
write.csv(tabla_419_final,
"salidas_tesis/tabla_4_19_validacion_cruzada_consolidada.csv",
row.names = FALSE)
# Excel
wb <- createWorkbook()
addWorksheet(wb, "Tabla 4.19")
addWorksheet(wb, "ROC")
addWorksheet(wb, "Sens")
addWorksheet(wb, "Spec")
writeData(wb, "Tabla 4.19", tabla_419_final, startRow = 1, startCol = 1)
writeData(wb, "ROC", tabla_roc_419, startRow = 1, startCol = 1)
writeData(wb, "Sens", tabla_sens_419, startRow = 1, startCol = 1)
writeData(wb, "Spec", tabla_spec_419, startRow = 1, startCol = 1)
# Estilos básicos
estilo_header <- createStyle(
textDecoration = "bold",
fgFill = "#1F4E78",
fontColour = "#FFFFFF",
halign = "center"
)
estilo_titulo <- createStyle(
textDecoration = "bold",
fgFill = "#D9EAF7"
)
# Aplicar encabezados a hojas métricas
addStyle(wb, "ROC", estilo_header, rows = 1, cols = 1:ncol(tabla_roc_419), gridExpand = TRUE)
addStyle(wb, "Sens", estilo_header, rows = 1, cols = 1:ncol(tabla_sens_419), gridExpand = TRUE)
addStyle(wb, "Spec", estilo_header, rows = 1, cols = 1:ncol(tabla_spec_419), gridExpand = TRUE)
# Detectar títulos en hoja consolidada
filas_titulo <- which(tabla_419_final$Modelos %in% c("ROC (AUC)", "Sens (Sensibilidad)", "Spec (Especificidad)")) + 1
if (length(filas_titulo) > 0) {
for (f in filas_titulo) {
addStyle(wb, "Tabla 4.19", estilo_titulo, rows = f, cols = 1:ncol(tabla_419_final), gridExpand = TRUE)
}
}
setColWidths(wb, "Tabla 4.19", cols = 1:7, widths = "auto")
setColWidths(wb, "ROC", cols = 1:7, widths = "auto")
setColWidths(wb, "Sens", cols = 1:7, widths = "auto")
setColWidths(wb, "Spec", cols = 1:7, widths = "auto")
freezePane(wb, "ROC", firstRow = TRUE)
freezePane(wb, "Sens", firstRow = TRUE)
freezePane(wb, "Spec", firstRow = TRUE)
saveWorkbook(
wb,
"salidas_tesis/tabla_4_19_validacion_cruzada.xlsx",
overwrite = TRUE
)
cat("\n=============================\n")##
## =============================cat("TABLA 4.19 EXPORTADA CORRECTAMENTE\n")## TABLA 4.19 EXPORTADA CORRECTAMENTEcat("=============================\n")## =============================cat("Archivos generados en salidas_tesis:\n")## Archivos generados en salidas_tesis:cat("- tabla_4_19_roc_cv.csv\n")## - tabla_4_19_roc_cv.csvcat("- tabla_4_19_sens_cv.csv\n")## - tabla_4_19_sens_cv.csvcat("- tabla_4_19_spec_cv.csv\n")## - tabla_4_19_spec_cv.csvcat("- tabla_4_19_validacion_cruzada_consolidada.csv\n")## - tabla_4_19_validacion_cruzada_consolidada.csvcat("- tabla_4_19_validacion_cruzada.xlsx\n")## - tabla_4_19_validacion_cruzada.xlsx############################################################
# TABLA DE LAS 6 MÉTRICAS DE DESEMPEÑO EN TEST
# Accuracy, Sensitivity, Specificity, Precision, F1 y AUC
############################################################
# Se asume que ya existe:
# tabla_metricas_test
library(dplyr)
library(openxlsx)
############################################################
# 1. TABLA BASE
############################################################
tabla_6_metricas <- tabla_metricas_test %>%
dplyr::select(
Modelo,
Accuracy,
Sensitivity,
Specificity,
Precision,
F1,
AUC
) %>%
mutate(
across(-Modelo, ~ round(.x, 4))
)
print(tabla_6_metricas)## Modelo Accuracy Sensitivity Specificity Precision F1 AUC
## Logit Logit 0.7192 0.7883 0.5758 0.7941 0.7912 0.7796
## Tree Tree 0.7340 0.7664 0.6667 0.8268 0.7955 0.7853
## RF RF 0.7192 0.9635 0.2121 0.7174 0.8224 0.7920
## SVM SVM 0.7192 0.9270 0.2879 0.7299 0.8167 0.7773
## KNN KNN 0.6798 0.8905 0.2424 0.7093 0.7896 0.7131
## NN NN 0.6946 0.7956 0.4848 0.7622 0.7786 0.7460
## NB NB 0.7291 0.7445 0.6970 0.8361 0.7876 0.7890############################################################
# 2. ORDENAR POR AUC (OPCIONAL)
############################################################
tabla_6_metricas_ordenada <- tabla_6_metricas %>%
arrange(desc(AUC))
print(tabla_6_metricas_ordenada)## Modelo Accuracy Sensitivity Specificity Precision F1 AUC
## RF RF 0.7192 0.9635 0.2121 0.7174 0.8224 0.7920
## NB NB 0.7291 0.7445 0.6970 0.8361 0.7876 0.7890
## Tree Tree 0.7340 0.7664 0.6667 0.8268 0.7955 0.7853
## Logit Logit 0.7192 0.7883 0.5758 0.7941 0.7912 0.7796
## SVM SVM 0.7192 0.9270 0.2879 0.7299 0.8167 0.7773
## NN NN 0.6946 0.7956 0.4848 0.7622 0.7786 0.7460
## KNN KNN 0.6798 0.8905 0.2424 0.7093 0.7896 0.7131############################################################
# 3. EXPORTAR CSV
############################################################
write.csv(
tabla_6_metricas_ordenada,
"salidas_tesis/tabla_6_metricas_desempeno_test.csv",
row.names = FALSE
)
############################################################
# 4. EXPORTAR EXCEL
############################################################
wb <- createWorkbook()
addWorksheet(wb, "Tabla 6 métricas")
writeData(
wb,
sheet = "Tabla 6 métricas",
x = tabla_6_metricas_ordenada,
startRow = 1,
startCol = 1
)
# Estilo encabezado
estilo_header <- createStyle(
textDecoration = "bold",
fgFill = "#1F4E78",
fontColour = "#FFFFFF",
halign = "center",
valign = "center"
)
# Estilo cuerpo
estilo_cuerpo <- createStyle(
halign = "center",
valign = "center",
border = "TopBottomLeftRight",
borderColour = "#B8C2CC"
)
addStyle(
wb, "Tabla 6 métricas",
style = estilo_header,
rows = 1, cols = 1:ncol(tabla_6_metricas_ordenada),
gridExpand = TRUE
)
addStyle(
wb, "Tabla 6 métricas",
style = estilo_cuerpo,
rows = 2:(nrow(tabla_6_metricas_ordenada) + 1),
cols = 1:ncol(tabla_6_metricas_ordenada),
gridExpand = TRUE,
stack = TRUE
)
setColWidths(wb, "Tabla 6 métricas", cols = 1:ncol(tabla_6_metricas_ordenada), widths = "auto")
freezePane(wb, "Tabla 6 métricas", firstRow = TRUE)
saveWorkbook(
wb,
"salidas_tesis/tabla_6_metricas_desempeno_test.xlsx",
overwrite = TRUE
)
cat("\nArchivos exportados:\n")##
## Archivos exportados:cat("- salidas_tesis/tabla_6_metricas_desempeno_test.csv\n")## - salidas_tesis/tabla_6_metricas_desempeno_test.csvcat("- salidas_tesis/tabla_6_metricas_desempeno_test.xlsx\n")## - salidas_tesis/tabla_6_metricas_desempeno_test.xlsx############################################################
# MATRICES DE CONFUSIÓN GRÁFICAS PARA LOS 7 MODELOS
# Estilo similar al gráfico adjunto
############################################################
#############################
# 1. PAQUETES
#############################
paquetes <- c("ggplot2", "dplyr", "pROC", "gridExtra", "grid", "patchwork")
instalar <- paquetes[!paquetes %in% installed.packages()[, "Package"]]
if (length(instalar) > 0) install.packages(instalar)
library(ggplot2)
library(dplyr)
library(pROC)
library(gridExtra)##
## Adjuntando el paquete: 'gridExtra'## The following object is masked from 'package:randomForest':
##
## combine## The following object is masked from 'package:dplyr':
##
## combinelibrary(grid)
library(patchwork)
#############################
# 2. DIRECTORIO DE SALIDA
#############################
dir_salida <- "salida_tesis"
if (!dir.exists(dir_salida)) {
dir.create(dir_salida, recursive = TRUE)
}
#############################
# 3. VALIDACIONES
#############################
if (!exists("test_model")) {
stop("No existe el objeto 'test_model'.")
}
if (!exists("modelo_logit_cv") ||
!exists("modelo_tree_cv") ||
!exists("modelo_rf_cv") ||
!exists("modelo_svm_cv") ||
!exists("modelo_knn_cv") ||
!exists("modelo_nn_cv") ||
!exists("modelo_nb_cv")) {
stop("Faltan uno o más modelos entrenados.")
}
# Asegurar codificación consistente
test_model$ESTADO_ACADEMICO <- factor(
test_model$ESTADO_ACADEMICO,
levels = c("Desertor", "No_Desertor")
)
#############################
# 4. LISTA DE MODELOS
#############################
modelos_finales <- list(
Logit = modelo_logit_cv,
Tree = modelo_tree_cv,
RF = modelo_rf_cv,
SVM = modelo_svm_cv,
KNN = modelo_knn_cv,
NN = modelo_nn_cv,
NB = modelo_nb_cv
)
#############################
# 5. FUNCIONES AUXILIARES
#############################
# Extraer probabilidades de la clase positiva
obtener_prob_desertor <- function(modelo, data_test) {
probs <- predict(modelo, newdata = data_test, type = "prob")
if (!"Desertor" %in% colnames(probs)) {
stop("La columna 'Desertor' no está presente en las probabilidades del modelo.")
}
probs[, "Desertor"]
}
# Calcular métricas por clase y métricas globales
calcular_reporte_clasificacion <- function(real, pred, prob_desertor) {
real <- factor(real, levels = c("Desertor", "No_Desertor"))
pred <- factor(pred, levels = c("Desertor", "No_Desertor"))
clases <- levels(real)
# AUC
roc_obj <- pROC::roc(
response = real,
predictor = prob_desertor,
levels = c("No_Desertor", "Desertor"),
direction = "<",
quiet = TRUE
)
auc_val <- as.numeric(pROC::auc(roc_obj))
# Métricas por clase
metricas_por_clase <- lapply(clases, function(cl) {
tp <- sum(real == cl & pred == cl, na.rm = TRUE)
fp <- sum(real != cl & pred == cl, na.rm = TRUE)
fn <- sum(real == cl & pred != cl, na.rm = TRUE)
support <- sum(real == cl, na.rm = TRUE)
precision <- ifelse((tp + fp) == 0, 0, tp / (tp + fp))
recall <- ifelse((tp + fn) == 0, 0, tp / (tp + fn))
f1 <- ifelse((precision + recall) == 0, 0,
2 * precision * recall / (precision + recall))
data.frame(
Clase = cl,
Precision = precision,
Recall = recall,
F1_Score = f1,
Support = support,
stringsAsFactors = FALSE
)
}) %>% bind_rows()
# Accuracy global
accuracy <- mean(real == pred, na.rm = TRUE)
# Métricas tomando Desertor como clase positiva
tp_pos <- sum(real == "Desertor" & pred == "Desertor", na.rm = TRUE)
fp_pos <- sum(real == "No_Desertor" & pred == "Desertor", na.rm = TRUE)
fn_pos <- sum(real == "Desertor" & pred == "No_Desertor", na.rm = TRUE)
tn_pos <- sum(real == "No_Desertor" & pred == "No_Desertor", na.rm = TRUE)
precision_pos <- ifelse((tp_pos + fp_pos) == 0, 0, tp_pos / (tp_pos + fp_pos))
recall_pos <- ifelse((tp_pos + fn_pos) == 0, 0, tp_pos / (tp_pos + fn_pos))
f1_pos <- ifelse((precision_pos + recall_pos) == 0, 0,
2 * precision_pos * recall_pos / (precision_pos + recall_pos))
specificity <- ifelse((tn_pos + fp_pos) == 0, 0, tn_pos / (tn_pos + fp_pos))
# Macro avg
macro_avg <- data.frame(
Clase = "Macro avg",
Precision = mean(metricas_por_clase$Precision, na.rm = TRUE),
Recall = mean(metricas_por_clase$Recall, na.rm = TRUE),
F1_Score = mean(metricas_por_clase$F1_Score, na.rm = TRUE),
Support = sum(metricas_por_clase$Support, na.rm = TRUE),
stringsAsFactors = FALSE
)
# Weighted avg
pesos <- metricas_por_clase$Support / sum(metricas_por_clase$Support)
weighted_avg <- data.frame(
Clase = "Weighted avg",
Precision = sum(metricas_por_clase$Precision * pesos, na.rm = TRUE),
Recall = sum(metricas_por_clase$Recall * pesos, na.rm = TRUE),
F1_Score = sum(metricas_por_clase$F1_Score * pesos, na.rm = TRUE),
Support = sum(metricas_por_clase$Support, na.rm = TRUE),
stringsAsFactors = FALSE
)
# Fila Accuracy estilo classification report
fila_accuracy <- data.frame(
Clase = "Accuracy",
Precision = NA,
Recall = NA,
F1_Score = accuracy,
Support = length(real),
stringsAsFactors = FALSE
)
tabla_reporte <- bind_rows(
metricas_por_clase %>% arrange(match(Clase, c("No_Desertor", "Desertor"))),
fila_accuracy,
macro_avg,
weighted_avg
)
list(
accuracy = accuracy,
precision = precision_pos,
recall = recall_pos,
f1 = f1_pos,
specificity = specificity,
auc = auc_val,
tabla_reporte = tabla_reporte,
roc = roc_obj
)
}
# Crear el gráfico compuesto
crear_grafico_matriz_confusion <- function(modelo, nombre_modelo, data_test) {
real <- factor(data_test$ESTADO_ACADEMICO, levels = c("Desertor", "No_Desertor"))
pred <- factor(
predict(modelo, newdata = data_test, type = "raw"),
levels = c("Desertor", "No_Desertor")
)
prob_desertor <- obtener_prob_desertor(modelo, data_test)
reporte <- calcular_reporte_clasificacion(real, pred, prob_desertor)
# Matriz para heatmap
tabla_mc <- table(
Real = factor(real, levels = c("Desertor", "No_Desertor")),
Prediccion = factor(pred, levels = c("No_Desertor", "Desertor"))
)
df_mc <- as.data.frame(tabla_mc)
colnames(df_mc) <- c("Real", "Prediccion", "Frecuencia")
# Heatmap
p_heat <- ggplot(df_mc, aes(x = Prediccion, y = Real, fill = Frecuencia)) +
geom_tile(color = "white", linewidth = 1) +
geom_text(aes(label = Frecuencia), size = 7, fontface = "bold", color = "black") +
scale_fill_gradient(low = "#E9EEF5", high = "#2F6FB0") +
scale_y_discrete(limits = c("No_Desertor", "Desertor")) +
labs(
title = "Matriz de confusión de evaluación de predicciones y resultados reales",
subtitle = paste("Modelo evaluado:", nombre_modelo),
x = "Predicción",
y = "Real"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
plot.subtitle = element_text(hjust = 0.5, size = 11),
axis.title = element_text(face = "bold"),
panel.grid = element_blank()
)
# Tabla tipo classification report
tabla_rep <- reporte$tabla_reporte %>%
mutate(
Precision = ifelse(is.na(Precision), "", sprintf("%.2f", Precision)),
Recall = ifelse(is.na(Recall), "", sprintf("%.2f", Recall)),
F1_Score = sprintf("%.2f", F1_Score),
Support = as.character(Support)
)
tg <- tableGrob(
tabla_rep,
rows = NULL,
theme = ttheme_minimal(
core = list(
fg_params = list(fontsize = 10),
bg_params = list(fill = rep("white", nrow(tabla_rep)), col = NA)
),
colhead = list(
fg_params = list(fontsize = 10, fontface = "bold"),
bg_params = list(fill = "grey95", col = NA)
)
)
)
p_tabla <- wrap_elements(full = tg)
# Cuadro de métricas globales
texto_metricas <- paste0(
"Accuracy: ", sprintf("%.4f", reporte$accuracy), "\n",
"Precision: ", sprintf("%.4f", reporte$precision), "\n",
"Recall: ", sprintf("%.4f", reporte$recall), "\n",
"F1 Score: ", sprintf("%.4f", reporte$f1), "\n",
"Specificity: ", sprintf("%.4f", reporte$specificity), "\n",
"AUC-ROC: ", sprintf("%.4f", reporte$auc)
)
p_metricas <- ggplot() +
annotate(
"label",
x = 0, y = 1,
label = texto_metricas,
hjust = 0, vjust = 1,
size = 4,
label.size = 0.25,
fill = "white"
) +
xlim(-0.1, 1) +
ylim(0, 1.05) +
theme_void()
grafico_final <- p_heat / (p_tabla | p_metricas) +
plot_layout(heights = c(3, 1.3), widths = c(1.6, 1))
list(
grafico = grafico_final,
reporte = reporte,
pred = pred,
real = real
)
}
#############################
# 6. GENERAR Y GUARDAR LOS 7 GRÁFICOS
#############################
resultados_graficos <- list()
for (nombre in names(modelos_finales)) {
cat("\nGenerando gráfico para:", nombre, "\n")
salida <- crear_grafico_matriz_confusion(
modelo = modelos_finales[[nombre]],
nombre_modelo = nombre,
data_test = test_model
)
resultados_graficos[[nombre]] <- salida
# Mostrar en pantalla
print(salida$grafico)
# Guardar PNG
ggsave(
filename = file.path(dir_salida, paste0("Matriz_confusion_", nombre, ".png")),
plot = salida$grafico,
width = 12,
height = 9,
dpi = 300
)
# Guardar PDF
ggsave(
filename = file.path(dir_salida, paste0("Matriz_confusion_", nombre, ".pdf")),
plot = salida$grafico,
width = 12,
height = 9
)
}##
## Generando gráfico para: Logit## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
##
## Generando gráfico para: Tree## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
##
## Generando gráfico para: RF## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
##
## Generando gráfico para: SVM## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
##
## Generando gráfico para: KNN## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
##
## Generando gráfico para: NN## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
##
## Generando gráfico para: NB## Warning in annotate("label", x = 0, y = 1, label = texto_metricas, hjust = 0, :
## Ignoring unknown parameters: `label.size`
#############################
# 7. RESUMEN GLOBAL DE MÉTRICAS
#############################
tabla_metricas_resumen <- lapply(names(resultados_graficos), function(nombre) {
rep <- resultados_graficos[[nombre]]$reporte
data.frame(
Modelo = nombre,
Accuracy = rep$accuracy,
Precision = rep$precision,
Recall = rep$recall,
F1_Score = rep$f1,
Specificity = rep$specificity,
AUC = rep$auc,
stringsAsFactors = FALSE
)
}) %>% bind_rows() %>% as.data.frame()
tabla_metricas_resumen <- tabla_metricas_resumen %>%
mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
as.data.frame()
print(tabla_metricas_resumen)## Modelo Accuracy Precision Recall F1_Score Specificity AUC
## 1 Logit 0.7192 0.7941 0.7883 0.7912 0.5758 0.7796
## 2 Tree 0.7340 0.8268 0.7664 0.7955 0.6667 0.7853
## 3 RF 0.7192 0.7174 0.9635 0.8224 0.2121 0.7920
## 4 SVM 0.7192 0.7299 0.9270 0.8167 0.2879 0.7773
## 5 KNN 0.6798 0.7093 0.8905 0.7896 0.2424 0.7131
## 6 NN 0.6946 0.7622 0.7956 0.7786 0.4848 0.7460
## 7 NB 0.7291 0.8361 0.7445 0.7876 0.6970 0.7890write.csv(
tabla_metricas_resumen,
file = file.path(dir_salida, "Resumen_metricas_matrices_confusion.csv"),
row.names = FALSE
)
cat("\n========================================\n")##
## ========================================cat("ARCHIVOS GENERADOS EN:", normalizePath(dir_salida), "\n")## ARCHIVOS GENERADOS EN: C:\Users\DELL\Desktop\Tesis-Estadistica\salida_tesiscat("Se generaron matrices de confusión gráficas para:\n")## Se generaron matrices de confusión gráficas para:cat(paste(names(modelos_finales), collapse = ", "), "\n")## Logit, Tree, RF, SVM, KNN, NN, NBcat("========================================\n")## ========================================############################################################
# 18. SALIDAS FINALES INTEGRADAS _finalisisma
# Gráficos y tablas finales para el Capítulo 4
############################################################
#############################
# 18.1 PAQUETES
#############################
paquetes_finales <- c("ggplot2", "dplyr", "tidyr", "stringr", "forcats", "pROC", "openxlsx", "caret", "rpart.plot")
invisible(lapply(paquetes_finales, instalar_si_falta))
#############################
# 18.2 CARPETAS DE SALIDA
#############################
dir_graficos_finalisisma <- "graficos_finalisisma"
dir_resultados_finalisisma <- "resultados_finalisisma"
dir.create(dir_graficos_finalisisma, showWarnings = FALSE, recursive = TRUE)
dir.create(dir_resultados_finalisisma, showWarnings = FALSE, recursive = TRUE)
#############################
# 18.3 OBJETOS BASE
#############################
# Modelos finales
if (!exists("modelos_finales")) {
modelos_finales <- list(
Logit = modelo_logit_cv,
Tree = modelo_tree_cv,
RF = modelo_rf_cv,
SVM = modelo_svm_cv,
KNN = modelo_knn_cv,
NN = modelo_nn_cv,
NB = modelo_nb_cv
)
}
# Base de prueba
if (!exists("test_model")) {
stop("No existe 'test_model'. Este bloque debe ejecutarse después del entrenamiento y la creación de la base de prueba.")
}
test_model$ESTADO_ACADEMICO <- factor(test_model$ESTADO_ACADEMICO, levels = c("Desertor", "No_Desertor"))
#############################
# 18.4 FUNCIÓN DE EVALUACIÓN SI NO EXISTE
#############################
if (!exists("evaluar_modelo_caret")) {
evaluar_modelo_caret <- function(modelo, nombre_modelo, test_data, variable_respuesta = "ESTADO_ACADEMICO") {
pred_clase <- predict(modelo, newdata = test_data, type = "raw")
pred_clase <- factor(pred_clase, levels = c("Desertor", "No_Desertor"))
pred_prob <- predict(modelo, newdata = test_data, type = "prob")[, "Desertor"]
real <- factor(test_data[[variable_respuesta]], levels = c("Desertor", "No_Desertor"))
mc <- confusionMatrix(pred_clase, real, positive = "Desertor")
roc_obj <- roc(response = real, predictor = pred_prob, levels = c("No_Desertor", "Desertor"), direction = "<")
auc_val <- as.numeric(auc(roc_obj))
precision_val <- as.numeric(mc$byClass["Pos Pred Value"])
recall_val <- as.numeric(mc$byClass["Sensitivity"])
f1_val <- ifelse((precision_val + recall_val) == 0, NA,
2 * precision_val * recall_val / (precision_val + recall_val))
resumen <- data.frame(
Modelo = nombre_modelo,
Accuracy = as.numeric(mc$overall["Accuracy"]),
Kappa = as.numeric(mc$overall["Kappa"]),
Sensitivity = recall_val,
Specificity = as.numeric(mc$byClass["Specificity"]),
Precision = precision_val,
F1 = f1_val,
AUC = auc_val,
row.names = NULL
)
list(nombre = nombre_modelo, confusion = mc, roc = roc_obj, resumen = resumen)
}
}
#############################
# 18.5 RESULTADOS EN TEST
#############################
if (!exists("resultados_test_lista")) {
resultados_test_lista <- lapply(names(modelos_finales), function(nombre) {
evaluar_modelo_caret(modelos_finales[[nombre]], nombre, test_model)
})
names(resultados_test_lista) <- names(modelos_finales)
}
tabla_metricas_final <- do.call(rbind, lapply(resultados_test_lista, function(x) x$resumen))
tabla_metricas_final <- as.data.frame(tabla_metricas_final)
rownames(tabla_metricas_final) <- NULL
cols_num <- intersect(c("Accuracy", "Kappa", "Sensitivity", "Specificity", "Precision", "F1", "AUC"), names(tabla_metricas_final))
tabla_metricas_final <- tabla_metricas_final %>%
mutate(across(any_of(cols_num), ~ suppressWarnings(as.numeric(as.character(.x))))) %>%
as.data.frame()
tabla_metricas_final <- tabla_metricas_final %>% arrange(desc(AUC), desc(Accuracy))
#############################
# 18.5A BLINDAJE DE OBJETOS
#############################
tabla_metricas_final <- as.data.frame(tabla_metricas_final, stringsAsFactors = FALSE)
if (!"Modelo" %in% names(tabla_metricas_final)) {
stop("tabla_metricas_final no contiene la columna 'Modelo'.")
}
#############################
# 18.6 RESULTADOS DE VALIDACIÓN CRUZADA
#############################
if (!exists("resultados_cv")) {
resultados_cv <- resamples(modelos_finales)
}
resumen_cv_final <- summary(resultados_cv)
#------------------------------------------
# Función robusta para extraer tablas CV
#------------------------------------------
extraer_tabla_cv <- function(obj_stats, prefijo) {
stat_df <- as.data.frame(obj_stats, stringsAsFactors = FALSE, check.names = FALSE)
tb <- data.frame(
Modelo = rownames(stat_df),
stat_df,
row.names = NULL,
check.names = FALSE,
stringsAsFactors = FALSE
)
# Renombrar columnas según existan realmente
nombres_originales <- names(tb)
nombres_nuevos <- sapply(nombres_originales, function(nm) {
if (nm == "Modelo") return("Modelo")
if (nm %in% c("Min.", "Min")) return(paste0(prefijo, "_Min"))
if (nm %in% c("1st Qu.", "1st Qu")) return(paste0(prefijo, "_Q1"))
if (nm == "Median") return(paste0(prefijo, "_Median"))
if (nm == "Mean") return(paste0(prefijo, "_Mean"))
if (nm %in% c("3rd Qu.", "3rd Qu")) return(paste0(prefijo, "_Q3"))
if (nm %in% c("Max.", "Max")) return(paste0(prefijo, "_Max"))
if (nm %in% c("NA's", "NAs", "NA")) return(paste0(prefijo, "_NAs"))
return(paste0(prefijo, "_", make.names(nm)))
})
names(tb) <- nombres_nuevos
as.data.frame(tb, stringsAsFactors = FALSE)
}
#------------------------------------------
# Extraer tablas para ROC, Sens y Spec
#------------------------------------------
tabla_cv_roc <- extraer_tabla_cv(resumen_cv_final$statistics$ROC, "ROC")
tabla_cv_sens <- extraer_tabla_cv(resumen_cv_final$statistics$Sens, "Sens")
tabla_cv_spec <- extraer_tabla_cv(resumen_cv_final$statistics$Spec, "Spec")
#------------------------------------------
# Unir tablas resumen CV
#------------------------------------------
tabla_cv_resumen_final <- tabla_cv_roc %>%
left_join(tabla_cv_sens, by = "Modelo") %>%
left_join(tabla_cv_spec, by = "Modelo") %>%
arrange(desc(ROC_Mean)) %>%
as.data.frame(stringsAsFactors = FALSE)
#------------------------------------------
# Tabla de generalización
#------------------------------------------
columnas_req <- c("Modelo", "AUC", "Sensitivity", "Specificity", "Accuracy")
columnas_opt <- c("Precision", "F1", "Kappa")
if (all(columnas_req %in% names(tabla_metricas_final))) {
columnas_existentes <- c(
columnas_req,
columnas_opt[columnas_opt %in% names(tabla_metricas_final)]
)
columnas_cv_req <- c("Modelo", "ROC_Mean", "Sens_Mean", "Spec_Mean")
if (all(columnas_cv_req %in% names(tabla_cv_resumen_final))) {
tabla_generalizacion_final <- tabla_metricas_final %>%
select(all_of(columnas_existentes)) %>%
rename(AUC_Test = AUC) %>%
left_join(
tabla_cv_resumen_final %>%
select(Modelo, AUC_CV = ROC_Mean, Sens_CV = Sens_Mean, Spec_CV = Spec_Mean),
by = "Modelo"
) %>%
mutate(
Brecha_AUC = AUC_Test - AUC_CV,
Brecha_Sens = Sensitivity - Sens_CV,
Brecha_Spec = Specificity - Spec_CV
) %>%
as.data.frame(stringsAsFactors = FALSE)
} else {
warning("No se encontraron todas las columnas esperadas en tabla_cv_resumen_final. No se generó tabla_generalizacion_final.")
tabla_generalizacion_final <- NULL
}
} else {
warning("No se encontraron todas las columnas requeridas en tabla_metricas_final. No se generó tabla_generalizacion_final.")
tabla_generalizacion_final <- NULL
}
#------------------------------------------
# Mostrar resultados
#------------------------------------------
print(tabla_cv_roc)## Modelo ROC_Min ROC_Q1 ROC_Median ROC_Mean ROC_Q3 ROC_Max ROC_NAs
## 1 Logit 0.6541667 0.7437500 0.7836174 0.7803046 0.8199219 0.8996212 0
## 2 Tree 0.6845703 0.7714962 0.7973633 0.7940493 0.8305827 0.8750000 0
## 3 RF 0.6458333 0.7341856 0.7910156 0.7854376 0.8286458 0.9091797 0
## 4 SVM 0.6416667 0.7183268 0.7613636 0.7602734 0.8041992 0.8867188 0
## 5 KNN 0.6318359 0.6966072 0.7434186 0.7483594 0.7955211 0.9062500 0
## 6 NN 0.6208333 0.7228634 0.7753906 0.7677115 0.8043176 0.8916667 0
## 7 NB 0.6041667 0.7292480 0.7802734 0.7714621 0.8144383 0.9104167 0print(tabla_cv_sens)## Modelo Sens_Min Sens_Q1 Sens_Median Sens_Mean Sens_Q3 Sens_Max Sens_NAs
## 1 Logit 0.68750 0.7500000 0.8125000 0.8071023 0.8473011 0.9696970 0
## 2 Tree 0.65625 0.7500000 0.7812500 0.7886995 0.8167614 0.9375000 0
## 3 RF 0.90625 0.9375000 0.9687500 0.9616477 0.9924242 1.0000000 0
## 4 SVM 0.68750 0.8437500 0.8787879 0.8777462 0.9375000 1.0000000 0
## 5 KNN 0.81250 0.9062500 0.9375000 0.9304924 0.9614110 1.0000000 0
## 6 NN 0.65625 0.7500000 0.8125000 0.7967487 0.8437500 0.9090909 0
## 7 NB 0.37500 0.6328125 0.7187500 0.6799874 0.7812500 0.8181818 0print(tabla_cv_spec)## Modelo Spec_Min Spec_Q1 Spec_Median Spec_Mean Spec_Q3 Spec_Max Spec_NAs
## 1 Logit 0.3333333 0.5000000 0.5625000 0.5548611 0.62500 0.7500 0
## 2 Tree 0.3125000 0.4666667 0.6125000 0.5850000 0.68750 0.8750 0
## 3 RF 0.0625000 0.1468750 0.1875000 0.2031944 0.25000 0.3750 0
## 4 SVM 0.0000000 0.2000000 0.3229167 0.3322222 0.43750 0.6875 0
## 5 KNN 0.0000000 0.1468750 0.2000000 0.2259722 0.31250 0.4375 0
## 6 NN 0.2500000 0.3750000 0.5000000 0.4906944 0.61875 0.7500 0
## 7 NB 0.4375000 0.6250000 0.7500000 0.7298611 0.81250 1.0000 0print(tabla_cv_resumen_final)## Modelo ROC_Min ROC_Q1 ROC_Median ROC_Mean ROC_Q3 ROC_Max ROC_NAs
## 1 Tree 0.6845703 0.7714962 0.7973633 0.7940493 0.8305827 0.8750000 0
## 2 RF 0.6458333 0.7341856 0.7910156 0.7854376 0.8286458 0.9091797 0
## 3 Logit 0.6541667 0.7437500 0.7836174 0.7803046 0.8199219 0.8996212 0
## 4 NB 0.6041667 0.7292480 0.7802734 0.7714621 0.8144383 0.9104167 0
## 5 NN 0.6208333 0.7228634 0.7753906 0.7677115 0.8043176 0.8916667 0
## 6 SVM 0.6416667 0.7183268 0.7613636 0.7602734 0.8041992 0.8867188 0
## 7 KNN 0.6318359 0.6966072 0.7434186 0.7483594 0.7955211 0.9062500 0
## Sens_Min Sens_Q1 Sens_Median Sens_Mean Sens_Q3 Sens_Max Sens_NAs
## 1 0.65625 0.7500000 0.7812500 0.7886995 0.8167614 0.9375000 0
## 2 0.90625 0.9375000 0.9687500 0.9616477 0.9924242 1.0000000 0
## 3 0.68750 0.7500000 0.8125000 0.8071023 0.8473011 0.9696970 0
## 4 0.37500 0.6328125 0.7187500 0.6799874 0.7812500 0.8181818 0
## 5 0.65625 0.7500000 0.8125000 0.7967487 0.8437500 0.9090909 0
## 6 0.68750 0.8437500 0.8787879 0.8777462 0.9375000 1.0000000 0
## 7 0.81250 0.9062500 0.9375000 0.9304924 0.9614110 1.0000000 0
## Spec_Min Spec_Q1 Spec_Median Spec_Mean Spec_Q3 Spec_Max Spec_NAs
## 1 0.3125000 0.4666667 0.6125000 0.5850000 0.68750 0.8750 0
## 2 0.0625000 0.1468750 0.1875000 0.2031944 0.25000 0.3750 0
## 3 0.3333333 0.5000000 0.5625000 0.5548611 0.62500 0.7500 0
## 4 0.4375000 0.6250000 0.7500000 0.7298611 0.81250 1.0000 0
## 5 0.2500000 0.3750000 0.5000000 0.4906944 0.61875 0.7500 0
## 6 0.0000000 0.2000000 0.3229167 0.3322222 0.43750 0.6875 0
## 7 0.0000000 0.1468750 0.2000000 0.2259722 0.31250 0.4375 0if (!is.null(tabla_generalizacion_final)) {
print(tabla_generalizacion_final)
}## Modelo AUC_Test Sensitivity Specificity Accuracy Precision F1
## 1 RF 0.7919708 0.9635036 0.2121212 0.7192118 0.7173913 0.8224299
## 2 NB 0.7890400 0.7445255 0.6969697 0.7290640 0.8360656 0.7876448
## 3 Tree 0.7853351 0.7664234 0.6666667 0.7339901 0.8267717 0.7954545
## 4 Logit 0.7796395 0.7883212 0.5757576 0.7192118 0.7941176 0.7912088
## 5 SVM 0.7773170 0.9270073 0.2878788 0.7192118 0.7298851 0.8167203
## 6 NN 0.7460186 0.7956204 0.4848485 0.6945813 0.7622378 0.7785714
## 7 KNN 0.7131166 0.8905109 0.2424242 0.6798030 0.7093023 0.7896440
## Kappa AUC_CV Sens_CV Spec_CV Brecha_AUC Brecha_Sens
## 1 0.2153658 0.7854376 0.9616477 0.2031944 0.0065332240 0.001855922
## 2 0.4169408 0.7714621 0.6799874 0.7298611 0.0175779536 0.064538174
## 3 0.4167287 0.7940493 0.7886995 0.5850000 -0.0087141790 -0.022276137
## 4 0.3626549 0.7803046 0.8071023 0.5548611 -0.0006651483 -0.018781105
## 5 0.2514071 0.7602734 0.8777462 0.3322222 0.0170435278 0.049261087
## 6 0.2872352 0.7677115 0.7967487 0.4906944 -0.0216929099 -0.001128299
## 7 0.1541124 0.7483594 0.9304924 0.2259722 -0.0352428079 -0.039981475
## Brecha_Spec
## 1 0.008926768
## 2 -0.032891414
## 3 0.081666667
## 4 0.020896465
## 5 -0.044343434
## 6 -0.005845960
## 7 0.016452020#############################
# 18.7 FUNCIÓN PARA GUARDAR PLOTS
#############################
guardar_plot_finalisisma <- function(plot_obj, nombre_base, width = 12, height = 8, dpi = 300) {
ggsave(filename = file.path(dir_graficos_finalisisma, paste0(nombre_base, ".png")),
plot = plot_obj, width = width, height = height, dpi = dpi)
ggsave(filename = file.path(dir_graficos_finalisisma, paste0(nombre_base, ".pdf")),
plot = plot_obj, width = width, height = height)
}
etiqueta_coma <- function(x, digits = 3) {
format(round(x, digits), nsmall = digits, decimal.mark = ",")
}
#############################
# 18.8 HEATMAP DE MÉTRICAS EN TEST
#############################
metricas_heat <- c("Accuracy", "Sensitivity", "Specificity", "Precision", "F1", "AUC")
metricas_heat <- metricas_heat[metricas_heat %in% names(tabla_metricas_final)]
orden_heat <- tabla_metricas_final$Modelo
heat_data <- tabla_metricas_final %>%
mutate(Modelo = factor(Modelo, levels = rev(orden_heat))) %>%
select(Modelo, all_of(metricas_heat)) %>%
pivot_longer(cols = -Modelo, names_to = "Metrica", values_to = "Valor") %>%
mutate(Metrica = factor(Metrica, levels = metricas_heat))
g_heatmap_finalisisma <- ggplot(heat_data, aes(x = Metrica, y = Modelo, fill = Valor)) +
geom_tile(color = "white", linewidth = 0.8) +
geom_text(aes(label = etiqueta_coma(Valor, 3)), size = 4.2) +
scale_fill_gradient(low = "#f7fbff", high = "#2171b5", limits = c(0, 1), name = "Valor") +
labs(
title = "Mapa de calor de las métricas de desempeño",
subtitle = "Conjunto de prueba independiente",
x = "Métrica",
y = "Modelo"
) +
theme_minimal(base_size = 15) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(angle = 25, hjust = 1),
plot.title = element_text(face = "bold"))
print(g_heatmap_finalisisma)
guardar_plot_finalisisma(g_heatmap_finalisisma, "heatmap_metricas_modelos_test_finalisisma", width = 13, height = 8)
#############################
# 18.9 SCATTER INTEGRADO
#############################
scatter_data <- tabla_metricas_final %>%
mutate(
Tipo_comportamiento = case_when(
Sensitivity >= 0.90 & Specificity < 0.40 ~ "Sesgado",
Sensitivity >= 0.85 & Specificity >= 0.45 ~ "Óptimo",
TRUE ~ "Intermedio"
),
Tipo_comportamiento = factor(Tipo_comportamiento, levels = c("Intermedio", "Óptimo", "Sesgado"))
)
g_scatter_finalisisma <- ggplot(
scatter_data,
aes(x = Specificity, y = Sensitivity, size = AUC, color = Tipo_comportamiento, label = Modelo)
) +
geom_point(alpha = 0.9) +
geom_text(size = 5, fontface = "bold", nudge_y = 0.012, show.legend = FALSE) +
scale_size_continuous(range = c(8, 16)) +
scale_color_manual(
values = c(
"Óptimo" = "#2ca25f",
"Intermedio" = "#fec44f",
"Sesgado" = "#de2d26"
)
) +
coord_cartesian(
xlim = c(0, max(0.60, max(scatter_data$Specificity, na.rm = TRUE) + 0.03)),
ylim = c(min(0.65, min(scatter_data$Sensitivity, na.rm = TRUE) - 0.03), 1.00)
) +
labs(
title = "Comparación integrada del desempeño de modelos",
subtitle = "Sensibilidad, especificidad y AUC en el conjunto de prueba",
x = "Especificidad",
y = "Sensibilidad",
size = "AUC",
color = "Tipo de comportamiento"
) +
theme_minimal(base_size = 15) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "right"
)
#############################
# 18.10 CURVAS ROC COMPARATIVAS
#############################
roc_df_final <- bind_rows(lapply(names(resultados_test_lista), function(nombre) {
roc_obj <- resultados_test_lista[[nombre]]$roc
data.frame(
Especificidad = rev(roc_obj$specificities),
Sensibilidad = rev(roc_obj$sensitivities),
Modelo = nombre,
stringsAsFactors = FALSE
)
}))
tabla_auc_test_final <- do.call(rbind, lapply(names(resultados_test_lista), function(nombre) {
data.frame(Modelo = nombre, AUC = as.numeric(resultados_test_lista[[nombre]]$resumen$AUC))
}))
g_roc_finalisisma <- ggplot(roc_df_final, aes(x = 1 - Especificidad, y = Sensibilidad, color = Modelo)) +
geom_line(linewidth = 1.1) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
labs(
title = "Curvas ROC de los modelos en el conjunto de prueba",
x = "1 - Especificidad",
y = "Sensibilidad",
color = "Modelo"
) +
theme_minimal(base_size = 15) +
theme(plot.title = element_text(face = "bold"))
print(g_roc_finalisisma)
guardar_plot_finalisisma(g_roc_finalisisma, "curvas_ROC_modelos_finalisisma", width = 12, height = 8)
#############################
# 18.11 ÁRBOL DE DECISIÓN FINAL
#############################
if (exists("modelo_tree_cv")) {
png(filename = file.path(dir_graficos_finalisisma, "arbol_decision_final_finalisisma.png"),
width = 2200, height = 1500, res = 220)
rpart.plot(
modelo_tree_cv$finalModel,
main = "Árbol de decisión final",
type = 4,
extra = 104,
under = TRUE,
fallen.leaves = TRUE,
yesno = 2,
faclen = 0,
varlen = 0,
tweak = 1.10,
branch = 1,
compress = TRUE,
shadow.col = 0,
nn = FALSE
)
dev.off()
pdf(file = file.path(dir_graficos_finalisisma, "arbol_decision_final_finalisisma.pdf"), width = 14, height = 9)
rpart.plot(
modelo_tree_cv$finalModel,
main = "Árbol de decisión final",
type = 4,
extra = 104,
under = TRUE,
fallen.leaves = TRUE,
yesno = 2,
faclen = 0,
varlen = 0,
tweak = 1.10,
branch = 1,
compress = TRUE,
shadow.col = 0,
nn = FALSE
)
dev.off()
}## png
## 2#############################
# 18.12 IMPORTANCIA DE CATEGORÍAS PREDICTORAS - RANDOM FOREST
#############################
if (exists("modelo_rf_cv")) {
imp_rf_raw <- varImp(modelo_rf_cv)$importance
imp_rf <- as.data.frame(imp_rf_raw, stringsAsFactors = FALSE, check.names = FALSE)
imp_rf <- data.frame(
Variable = rownames(imp_rf),
imp_rf,
row.names = NULL,
check.names = FALSE,
stringsAsFactors = FALSE
)
cols_imp <- setdiff(names(imp_rf), "Variable")
if (length(cols_imp) == 0) {
stop("No se encontró ninguna columna de importancia en varImp(modelo_rf_cv).")
}
col_importancia <- if ("Overall" %in% cols_imp) "Overall" else cols_imp[1]
imp_rf_plot <- imp_rf %>%
arrange(.data[[col_importancia]]) %>%
tail(min(20, nrow(.))) %>%
mutate(Variable = factor(Variable, levels = Variable)) %>%
as.data.frame(stringsAsFactors = FALSE)
g_imp_rf_finalisisma <- ggplot(imp_rf_plot, aes(x = Variable, y = .data[[col_importancia]])) +
geom_col(fill = "#2b8cbe") +
coord_flip() +
labs(
title = "Importancia de categorías predictoras",
subtitle = "Modelo Random Forest",
x = NULL,
y = "Importancia"
) +
theme_minimal(base_size = 15) +
theme(plot.title = element_text(face = "bold"))
print(g_imp_rf_finalisisma)
guardar_plot_finalisisma(
g_imp_rf_finalisisma,
"importancia_variables_random_forest_finalisisma",
width = 11,
height = 8
)
}
#############################
# 18.13 BOXPLOTS PANELADOS DE VALIDACIÓN CRUZADA
#############################
cv_values_final <- as.data.frame(resultados_cv$values, stringsAsFactors = FALSE, check.names = FALSE)
if (!"Resample" %in% names(cv_values_final)) {
stop("El objeto resultados_cv$values no contiene la columna Resample.")
}
cv_largo_final <- cv_values_final %>%
tidyr::pivot_longer(cols = -Resample, names_to = "Modelo_Metrica", values_to = "Valor") %>%
dplyr::filter(stringr::str_detect(Modelo_Metrica, "~(ROC|Sens|Spec)$")) %>%
tidyr::separate(col = Modelo_Metrica, into = c("Modelo", "Metrica"), sep = "~") %>%
as.data.frame(stringsAsFactors = FALSE)
cv_largo_final$Metrica <- dplyr::recode(
cv_largo_final$Metrica,
"Spec" = "Especificidad",
"ROC" = "ROC",
"Sens" = "Sensibilidad"
)
orden_modelos_cv <- intersect(
c("KNN", "Logit", "NB", "NN", "RF", "SVM", "Tree"),
unique(cv_largo_final$Modelo)
)
cv_largo_final$Modelo <- factor(cv_largo_final$Modelo, levels = orden_modelos_cv)
cv_largo_final$Metrica <- factor(
cv_largo_final$Metrica,
levels = c("Especificidad", "ROC", "Sensibilidad")
)
g_cv_paneles_finalisisma <- ggplot(cv_largo_final, aes(x = Modelo, y = Valor)) +
geom_boxplot(width = 0.6, fill = "white", color = "black", outlier.shape = 16, outlier.size = 2) +
facet_wrap(~ Metrica, ncol = 1, scales = "free_y") +
labs(
title = "Resultados de validación cruzada por métrica",
x = "Modelo",
y = "Valor"
) +
theme_minimal(base_size = 16) +
theme(
strip.text = element_text(size = 14),
plot.title = element_text(size = 19, face = "bold", hjust = 0.5),
axis.text.x = element_text(angle = 25, hjust = 1)
)
print(g_cv_paneles_finalisisma)
guardar_plot_finalisisma(
g_cv_paneles_finalisisma,
"Boxplots_CV_paneles_finalisisma",
width = 12,
height = 10
)
tabla_cv_resumen_paneles_final <- cv_largo_final %>%
dplyr::group_by(Modelo, Metrica) %>%
dplyr::summarise(
Media = mean(Valor, na.rm = TRUE),
Mediana = median(Valor, na.rm = TRUE),
DE = sd(Valor, na.rm = TRUE),
Min = min(Valor, na.rm = TRUE),
Max = max(Valor, na.rm = TRUE),
.groups = "drop"
) %>%
as.data.frame(stringsAsFactors = FALSE)
#############################
# 18.14 TABLA ENTRE ENFOQUES
#############################
tabla_enfoques_final <- data.frame(
Enfoque = c("Modelización de datos", "Modelización algorítmica"),
Evidencia_principal = c(
"La regresión logística mostró desempeño competitivo, estabilidad en validación cruzada e interpretabilidad de coeficientes.",
"Los modelos algorítmicos evidenciaron patrones diferenciados: mayor sensibilidad en NB, SVM y KNN; mayor especificidad en RF; mejor equilibrio en NN."
),
stringsAsFactors = FALSE
)
#############################
# 18.15 EXPORTACIÓN INTEGRADA A EXCEL Y CSV
#############################
wb_final <- createWorkbook()
addWorksheet(wb_final, "Resultados_Test")
writeData(wb_final, "Resultados_Test", tabla_metricas_final)
addWorksheet(wb_final, "Ranking_Modelos")
writeData(wb_final, "Ranking_Modelos", tabla_metricas_final %>% arrange(desc(AUC), desc(Accuracy)))
addWorksheet(wb_final, "Resumen_CV")
writeData(wb_final, "Resumen_CV", tabla_cv_resumen_final)
addWorksheet(wb_final, "CV_Paneles")
writeData(wb_final, "CV_Paneles", tabla_cv_resumen_paneles_final)
addWorksheet(wb_final, "AUC_Test")
writeData(wb_final, "AUC_Test", tabla_auc_test_final)
addWorksheet(wb_final, "Enfoques")
writeData(wb_final, "Enfoques", tabla_enfoques_final)
if (!is.null(tabla_generalizacion_final)) {
addWorksheet(wb_final, "Generalizacion")
writeData(wb_final, "Generalizacion", tabla_generalizacion_final)
}
saveWorkbook(
wb_final,
file = file.path(dir_resultados_finalisisma, "Resumen_y_metricas_modelos_finalisisma.xlsx"),
overwrite = TRUE
)
write.csv(tabla_metricas_final, file.path(dir_resultados_finalisisma, "tabla_metricas_test_finalisisma.csv"), row.names = FALSE)
write.csv(tabla_cv_resumen_final, file.path(dir_resultados_finalisisma, "tabla_resumen_cv_finalisisma.csv"), row.names = FALSE)
write.csv(tabla_cv_resumen_paneles_final, file.path(dir_resultados_finalisisma, "tabla_cv_paneles_finalisisma.csv"), row.names = FALSE)
write.csv(tabla_auc_test_final, file.path(dir_resultados_finalisisma, "tabla_auc_test_finalisisma.csv"), row.names = FALSE)
write.csv(tabla_enfoques_final, file.path(dir_resultados_finalisisma, "tabla_enfoques_finalisisma.csv"), row.names = FALSE)
if (!is.null(tabla_generalizacion_final)) {
write.csv(tabla_generalizacion_final, file.path(dir_resultados_finalisisma, "tabla_generalizacion_finalisisma.csv"), row.names = FALSE)
}
#############################
# 18.16 MENSAJE FINAL
#############################
cat("\n====================================================\n")##
## ====================================================cat("SALIDAS FINALES _finalisisma GENERADAS CORRECTAMENTE\n")## SALIDAS FINALES _finalisisma GENERADAS CORRECTAMENTEcat("Gráficos guardados en: ", normalizePath(dir_graficos_finalisisma), "\n")## Gráficos guardados en: C:\Users\DELL\Desktop\Tesis-Estadistica\graficos_finalisismacat("Resultados guardados en: ", normalizePath(dir_resultados_finalisisma), "\n")## Resultados guardados en: C:\Users\DELL\Desktop\Tesis-Estadistica\resultados_finalisismacat("Archivo Excel: Resumen_y_metricas_modelos_finalisisma.xlsx\n")## Archivo Excel: Resumen_y_metricas_modelos_finalisisma.xlsxcat("====================================================\n")## ====================================================############################################################
# MODELO LOGÍSTICO FINAL - TESIS DESERCIÓN FACEN
# BLOQUE COMPLETO PARA 4.6 Y 4.7
############################################################
#############################
# 0. PAQUETES
#############################
paquetes <- c(
"readxl", "dplyr", "ggplot2", "caret", "pROC",
"car", "ResourceSelection", "openxlsx"
)
instalar_si_falta <- function(pkg) {
if (!require(pkg, character.only = TRUE)) {
install.packages(pkg, dependencies = TRUE)
library(pkg, character.only = TRUE)
}
}
invisible(lapply(paquetes, instalar_si_falta))
#############################
# 1. CARPETAS DE SALIDA
#############################
dir_resultados <- "resultados_modelo_logistico_final"
dir_graficos <- "graficos_modelo_logistico_final"
dir.create(dir_resultados, showWarnings = FALSE, recursive = TRUE)
dir.create(dir_graficos, showWarnings = FALSE, recursive = TRUE)
#############################
# 2. CARGA DE DATOS
#############################
# OPCIÓN A:
# Si ya tienes la base cargada en un data.frame llamado 'tesis',
# deja comentada la siguiente línea.
# tesis <- readxl::read_excel("desercion.xlsx")
# Validación mínima
if (!exists("tesis")) {
stop("No existe el objeto 'tesis'. Carga primero la base o ajusta la ruta del archivo Excel.")
}
#############################
# 3. PREPARACIÓN DE LA BASE
#############################
variables_candidatas <- c(
"RESIDENCIA", "CARRERA", "TIPO_INGRESO", "MODALIDAD", "RENDIMIENTO",
"TRABAJA", "SOLVENTAR", "ESTUDIOS_PADRES", "SEXO", "TIPO_COL",
"INGRESO", "EST_CIVIL", "NIVEL_SOCIO", "EDAD"
)
vars_necesarias <- c("ESTADO_ACADEMICO", variables_candidatas)
faltantes <- setdiff(vars_necesarias, names(tesis))
if (length(faltantes) > 0) {
stop(paste("Faltan estas variables en la base:", paste(faltantes, collapse = ", ")))
}
tesis_modelo <- tesis[, vars_necesarias]
tesis_modelo <- tesis_modelo[complete.cases(tesis_modelo), ]
# Asegurar tipos
vars_factor <- setdiff(vars_necesarias, c("EDAD"))
tesis_modelo[vars_factor] <- lapply(tesis_modelo[vars_factor], factor)
tesis_modelo$EDAD <- as.numeric(tesis_modelo$EDAD)
# Variable respuesta factor y binaria explícita
# IMPORTANTE: aquí se define sin ambigüedad que el evento de interés es DESERTOR = 1
tesis_modelo$ESTADO_ACADEMICO <- factor(
tesis_modelo$ESTADO_ACADEMICO,
levels = c("No_Desertor", "Desertor")
)
# Si en tu base original los niveles vienen al revés, este bloque corrige:
if (!all(c("No_Desertor", "Desertor") %in% levels(tesis_modelo$ESTADO_ACADEMICO))) {
niveles_actuales <- levels(tesis_modelo$ESTADO_ACADEMICO)
message("Niveles actuales de ESTADO_ACADEMICO: ", paste(niveles_actuales, collapse = ", "))
stop("Ajusta manualmente la codificación de ESTADO_ACADEMICO a 'No_Desertor' y 'Desertor'.")
}
tesis_modelo$Y_LOGIT <- ifelse(tesis_modelo$ESTADO_ACADEMICO == "Desertor", 1, 0)
# Función auxiliar para fijar referencias si existen
relevel_if_exists <- function(x, refs) {
x <- factor(x)
for (r in refs) {
if (r %in% levels(x)) {
return(stats::relevel(x, ref = r))
}
}
x
}
# Referencias sugeridas para alinear con la tesis
tesis_modelo$RESIDENCIA <- relevel_if_exists(tesis_modelo$RESIDENCIA, c("ASUNCIÓN", "Asunción"))
tesis_modelo$CARRERA <- relevel_if_exists(tesis_modelo$CARRERA, c("ESTADISTICA-PRES", "Estadística presencial", "ESTADISTICA PRESENCIAL"))
tesis_modelo$TIPO_INGRESO <- relevel_if_exists(tesis_modelo$TIPO_INGRESO, c("INGRESO", "Ingreso por examen", "EXAMEN DE INGRESO"))
tesis_modelo$RENDIMIENTO <- relevel_if_exists(tesis_modelo$RENDIMIENTO, c("APROBADO", "Aprobado"))
tesis_modelo$TRABAJA <- relevel_if_exists(tesis_modelo$TRABAJA, c("SI", "Sí", "TRABAJA"))
#############################
# 4. PARTICIÓN TRAIN / TEST
#############################
set.seed(123)
idx_train <- caret::createDataPartition(
y = tesis_modelo$ESTADO_ACADEMICO,
p = 0.70,
list = FALSE
)
train_model <- tesis_modelo[idx_train, ]
test_model <- tesis_modelo[-idx_train, ]
#############################
# 5. MODELO LOGÍSTICO INICIAL Y STEPWISE
#############################
form_stepwise <- as.formula(
paste("Y_LOGIT ~", paste(variables_candidatas, collapse = " + "))
)
modelo_nulo <- glm(
Y_LOGIT ~ 1,
data = train_model,
family = binomial(link = "logit")
)
modelo_logit_inicial <- glm(
form_stepwise,
data = train_model,
family = binomial(link = "logit")
)
modelo_logit_final <- step(
modelo_logit_inicial,
direction = "both",
trace = FALSE
)
cat("\n=============================\n")##
## =============================cat("MODELO LOGÍSTICO FINAL SELECCIONADO POR STEPWISE\n")## MODELO LOGÍSTICO FINAL SELECCIONADO POR STEPWISEcat("=============================\n")## =============================print(formula(modelo_logit_final))## Y_LOGIT ~ CARRERA + TIPO_INGRESO + RENDIMIENTO + TRABAJAprint(summary(modelo_logit_final))##
## Call:
## glm(formula = Y_LOGIT ~ CARRERA + TIPO_INGRESO + RENDIMIENTO +
## TRABAJA, family = binomial(link = "logit"), data = train_model)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 0.4598 0.2976 1.545
## CARRERAESTADISTICA-SEMI 0.4591 0.5026 0.914
## CARRERAMATEMATICA -0.3829 0.3071 -1.247
## CARRERAEDUCACION MATEMATICA-PRES -1.2584 0.4134 -3.044
## CARRERAEDUCACION MATEMATICA-SEMI 0.3080 0.3533 0.872
## TIPO_INGRESOADMISION DIRECTA -1.6173 0.5699 -2.838
## RENDIMIENTONO APROBADO 2.6019 0.3361 7.742
## TRABAJANO -0.4758 0.2756 -1.726
## Pr(>|z|)
## (Intercept) 0.12227
## CARRERAESTADISTICA-SEMI 0.36096
## CARRERAMATEMATICA 0.21248
## CARRERAEDUCACION MATEMATICA-PRES 0.00233 **
## CARRERAEDUCACION MATEMATICA-SEMI 0.38332
## TIPO_INGRESOADMISION DIRECTA 0.00454 **
## RENDIMIENTONO APROBADO 0.00000000000000978 ***
## TRABAJANO 0.08427 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 606.02 on 478 degrees of freedom
## Residual deviance: 449.20 on 471 degrees of freedom
## AIC: 465.2
##
## Number of Fisher Scoring iterations: 5# Función auxiliar para convertir fórmulas largas en un solo texto
formula_a_texto <- function(f) {
paste(deparse(f), collapse = " ")
}
# Resumen del stepwise para el cuerpo de la tesis
formula_inicial_txt <- formula_a_texto(form_stepwise)
formula_final_txt <- formula_a_texto(formula(modelo_logit_final))
tabla_stepwise_resumen <- data.frame(
Item = c(
"Modelo inicial",
"Criterio de selección",
"Dirección del stepwise",
"Modelo final"
),
Valor = c(
formula_inicial_txt,
"AIC",
"both",
formula_final_txt
),
stringsAsFactors = FALSE
)
print(tabla_stepwise_resumen)## Item
## 1 Modelo inicial
## 2 Criterio de selección
## 3 Dirección del stepwise
## 4 Modelo final
## Valor
## 1 Y_LOGIT ~ RESIDENCIA + CARRERA + TIPO_INGRESO + MODALIDAD + RENDIMIENTO + TRABAJA + SOLVENTAR + ESTUDIOS_PADRES + SEXO + TIPO_COL + INGRESO + EST_CIVIL + NIVEL_SOCIO + EDAD
## 2 AIC
## 3 both
## 4 Y_LOGIT ~ CARRERA + TIPO_INGRESO + RENDIMIENTO + TRABAJA#############################
# 6. TABLA DE WALD E IC95%
#############################
coef_mat <- summary(modelo_logit_final)$coefficients
ic_wald <- confint.default(modelo_logit_final)
tabla_wald <- data.frame(
Variable = rownames(coef_mat),
Estimate = coef_mat[, "Estimate"],
Std_Error = coef_mat[, "Std. Error"],
z_value = coef_mat[, "z value"],
p_value = coef_mat[, "Pr(>|z|)"],
IC95_LI = ic_wald[, 1],
IC95_LS = ic_wald[, 2],
row.names = NULL,
stringsAsFactors = FALSE
)
#############################
# 7. TABLA DE ODDS RATIO
#############################
tabla_or <- tabla_wald
tabla_or$OR <- exp(tabla_or$Estimate)
tabla_or$OR_IC95_LI <- exp(tabla_or$IC95_LI)
tabla_or$OR_IC95_LS <- exp(tabla_or$IC95_LS)
tabla_or <- tabla_or[, c(
"Variable", "Estimate", "Std_Error", "z_value", "p_value",
"OR", "OR_IC95_LI", "OR_IC95_LS"
)]
#############################
# 8. PREDICCIONES TRAIN / TEST
#############################
prob_train <- predict(modelo_logit_final, newdata = train_model, type = "response")
prob_test <- predict(modelo_logit_final, newdata = test_model, type = "response")
#############################
# 9. MÉTRICAS DE CLASIFICACIÓN
#############################
metricas_clasificacion <- function(y_true, prob, umbral = 0.50) {
pred_bin <- ifelse(prob >= umbral, 1, 0)
TP <- sum(pred_bin == 1 & y_true == 1)
TN <- sum(pred_bin == 0 & y_true == 0)
FP <- sum(pred_bin == 1 & y_true == 0)
FN <- sum(pred_bin == 0 & y_true == 1)
sensitivity <- ifelse((TP + FN) == 0, NA, TP / (TP + FN))
specificity <- ifelse((TN + FP) == 0, NA, TN / (TN + FP))
accuracy <- (TP + TN) / (TP + TN + FP + FN)
precision <- ifelse((TP + FP) == 0, NA, TP / (TP + FP))
f1 <- ifelse(is.na(precision) | is.na(sensitivity) | (precision + sensitivity) == 0,
NA, 2 * precision * sensitivity / (precision + sensitivity))
balanced <- mean(c(sensitivity, specificity), na.rm = TRUE)
youden <- sensitivity + specificity - 1
# Kappa
total <- TP + TN + FP + FN
po <- accuracy
pe <- (((TP + FP) * (TP + FN)) + ((FN + TN) * (FP + TN))) / (total^2)
kappa <- ifelse((1 - pe) == 0, NA, (po - pe) / (1 - pe))
data.frame(
Umbral = umbral,
Accuracy = accuracy,
Sensitivity = sensitivity,
Specificity = specificity,
Precision = precision,
F1 = f1,
Balanced_Accuracy = balanced,
Youden = youden,
Kappa = kappa,
TP = TP,
TN = TN,
FP = FP,
FN = FN
)
}
tabla_metricas_test_050 <- metricas_clasificacion(test_model$Y_LOGIT, prob_test, umbral = 0.50)
#############################
# 10. MATRIZ DE CONFUSIÓN
#############################
pred_factor_050 <- factor(
ifelse(prob_test >= 0.50, "Desertor", "No_Desertor"),
levels = c("Desertor", "No_Desertor")
)
ref_factor_test <- factor(
ifelse(test_model$Y_LOGIT == 1, "Desertor", "No_Desertor"),
levels = c("Desertor", "No_Desertor")
)
matriz_confusion_050 <- caret::confusionMatrix(
data = pred_factor_050,
reference = ref_factor_test,
positive = "Desertor"
)
#############################
# 11. CURVA ROC Y AUC
#############################
roc_train <- pROC::roc(
response = train_model$Y_LOGIT,
predictor = prob_train,
levels = c(0, 1),
direction = "<"
)
roc_test <- pROC::roc(
response = test_model$Y_LOGIT,
predictor = prob_test,
levels = c(0, 1),
direction = "<"
)
auc_train <- as.numeric(pROC::auc(roc_train))
auc_test <- as.numeric(pROC::auc(roc_test))
tabla_auc <- data.frame(
Escenario = c("Train", "Test"),
AUC = c(auc_train, auc_test),
stringsAsFactors = FALSE
)
tabla_auc$Diferencia_Absoluta <- c(NA, abs(auc_test - auc_train))
#############################
# 12. UMBRAL ÓPTIMO Y YOUDEN
#############################
umbrales <- seq(0.01, 0.99, by = 0.01)
tabla_umbrales <- do.call(
rbind,
lapply(umbrales, function(u) metricas_clasificacion(test_model$Y_LOGIT, prob_test, umbral = u))
)
mejor_youden <- tabla_umbrales[which.max(tabla_umbrales$Youden), ]
# Métricas comparativas: 0.50 vs Youden
tabla_decision_umbral <- dplyr::bind_rows(
dplyr::mutate(tabla_metricas_test_050, Criterio = "Umbral fijo 0.50"),
dplyr::mutate(mejor_youden, Criterio = "Umbral óptimo Youden")
) %>%
dplyr::select(Criterio, everything())
#############################
# 13. RAZÓN DE VEROSIMILITUD
#############################
tabla_lr <- anova(modelo_nulo, modelo_logit_final, test = "Chisq")
tabla_lr <- as.data.frame(tabla_lr)
tabla_lr$Modelo <- c("Nulo", "Final")
#############################
# 14. HOSMER-LEMESHOW
#############################
hl <- ResourceSelection::hoslem.test(
x = train_model$Y_LOGIT,
y = fitted(modelo_logit_final),
g = 10
)
tabla_hl <- data.frame(
Chi_square = as.numeric(hl$statistic),
gl = as.numeric(hl$parameter),
p_value = hl$p.value
)
#############################
# 15. MULTICOLINEALIDAD (VIF / GVIF)
#############################
vif_raw <- tryCatch(
car::vif(modelo_logit_final),
error = function(e) NULL
)
if (is.null(vif_raw)) {
tabla_vif <- data.frame(
Variable = "No disponible",
VIF = NA
)
} else if (is.matrix(vif_raw)) {
tabla_vif <- data.frame(
Variable = rownames(vif_raw),
GVIF = vif_raw[, 1],
Df = vif_raw[, 2],
GVIF_ajustado = vif_raw[, 3],
row.names = NULL
)
} else {
tabla_vif <- data.frame(
Variable = names(vif_raw),
VIF = as.numeric(vif_raw),
row.names = NULL
)
}
#############################
# 16. LINEALIDAD EN EL LOGIT PARA EDAD
#############################
tabla_box_tidwell <- NULL
datos_linealidad <- NULL
if ("EDAD" %in% names(train_model) && all(train_model$EDAD > 0, na.rm = TRUE)) {
modelo_bt <- glm(
Y_LOGIT ~ EDAD + I(EDAD * log(EDAD)),
data = train_model,
family = binomial(link = "logit")
)
bt_coef <- summary(modelo_bt)$coefficients
tabla_box_tidwell <- data.frame(
Variable = rownames(bt_coef),
Estimate = bt_coef[, "Estimate"],
Std_Error = bt_coef[, "Std. Error"],
z_value = bt_coef[, "z value"],
p_value = bt_coef[, "Pr(>|z|)"],
row.names = NULL
)
eps <- 1e-4
datos_linealidad <- train_model %>%
dplyr::filter(!is.na(EDAD), EDAD > 0) %>%
dplyr::mutate(decil = dplyr::ntile(EDAD, 10)) %>%
dplyr::group_by(decil) %>%
dplyr::summarise(
Edad_media = mean(EDAD, na.rm = TRUE),
Prob_media = mean(Y_LOGIT, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::mutate(
Prob_media = pmin(pmax(Prob_media, eps), 1 - eps),
Logit_empirico = log(Prob_media / (1 - Prob_media))
)
graf_linealidad <- ggplot(datos_linealidad, aes(x = Edad_media, y = Logit_empirico)) +
geom_point(size = 2) +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
labs(
title = "Evaluación exploratoria de linealidad del logit para EDAD",
x = "Edad media por decil",
y = "Logit empírico de la probabilidad de deserción"
) +
theme_minimal()
ggsave(
filename = file.path(dir_graficos, "linealidad_logit_edad.png"),
plot = graf_linealidad,
width = 10,
height = 6,
dpi = 300
)
}## `geom_smooth()` using formula = 'y ~ x'#############################
# 17. RESIDUOS DE PEARSON, COOK Y LEVERAGE
#############################
res_pearson <- residuals(modelo_logit_final, type = "pearson")
res_deviance <- residuals(modelo_logit_final, type = "deviance")
cook_d <- cooks.distance(modelo_logit_final)
leverage <- hatvalues(modelo_logit_final)
tabla_residuos_pearson <- data.frame(
Min = min(res_pearson, na.rm = TRUE),
Q1 = as.numeric(quantile(res_pearson, 0.25, na.rm = TRUE)),
Mediana = median(res_pearson, na.rm = TRUE),
Media = mean(res_pearson, na.rm = TRUE),
Q3 = as.numeric(quantile(res_pearson, 0.75, na.rm = TRUE)),
Max = max(res_pearson, na.rm = TRUE)
)
tabla_cook <- data.frame(
Min = min(cook_d, na.rm = TRUE),
Q1 = as.numeric(quantile(cook_d, 0.25, na.rm = TRUE)),
Mediana = median(cook_d, na.rm = TRUE),
Media = mean(cook_d, na.rm = TRUE),
Q3 = as.numeric(quantile(cook_d, 0.75, na.rm = TRUE)),
Max = max(cook_d, na.rm = TRUE)
)
umbral_cook <- 4 / nrow(train_model)
casos_influyentes <- data.frame(
Observacion = seq_along(cook_d),
Cook = as.numeric(cook_d),
Leverage = as.numeric(leverage),
Pearson = as.numeric(res_pearson),
Deviance = as.numeric(res_deviance),
Prob_Ajustada = fitted(modelo_logit_final)
) %>%
dplyr::filter(Cook > umbral_cook | abs(Pearson) > 2)
#############################
# 18. GRÁFICOS DE DIAGNÓSTICO
#############################
# 18.1 Panel Pearson + Cook
png(
filename = file.path(dir_graficos, "diagnostico_modelo_logistico_final.png"),
width = 1400,
height = 700,
res = 140
)
par(mfrow = c(1, 2))
plot(
res_pearson,
pch = 16,
main = "Residuos de Pearson",
xlab = "Índice",
ylab = "Residuo"
)
abline(h = c(-2, 2), col = "red", lty = 2)
plot(
cook_d,
pch = 16,
main = "Distancia de Cook",
xlab = "Índice",
ylab = "Cook"
)
abline(h = umbral_cook, col = "red", lty = 2)
dev.off()## png
## 2# 18.2 Distancia de Cook
png(
filename = file.path(dir_graficos, "distancia_cook.png"),
width = 1200,
height = 800,
res = 140
)
barplot(
cook_d,
border = NA,
main = "Distancia de Cook",
xlab = "Observación",
ylab = "Cook's distance",
col = "gray40"
)
abline(h = umbral_cook, col = "red", lty = 2)
dev.off()## png
## 2# 18.3 Leverage vs residuos de deviance
graf_leverage <- ggplot(
data.frame(
Leverage = leverage,
Residuo_Deviance = res_deviance,
Cook = cook_d
),
aes(x = Leverage, y = Residuo_Deviance, size = Cook)
) +
geom_point(shape = 21, fill = "gray40", color = "black", alpha = 0.8) +
geom_hline(yintercept = c(-2, 2), color = "red", linetype = 2, linewidth = 0.4) +
geom_vline(xintercept = 2 * mean(leverage), color = "blue", linetype = 2, linewidth = 0.4) +
labs(
title = "Leverage vs residuos de deviance",
x = "Leverage",
y = "Residuo de deviance",
size = "Cook"
) +
theme_minimal()
ggsave(
filename = file.path(dir_graficos, "leverage_residuos.png"),
plot = graf_leverage,
width = 10,
height = 7,
dpi = 300
)
#############################
# 19. CURVA ROC DEL MODELO FINAL
#############################
roc_df <- data.frame(
Specificity = rev(roc_test$specificities),
Sensitivity = rev(roc_test$sensitivities)
)
graf_roc <- ggplot(roc_df, aes(x = 1 - Specificity, y = Sensitivity)) +
geom_line(linewidth = 1, color = "black") +
geom_abline(slope = 1, intercept = 0, linetype = 2, color = "gray60") +
labs(
title = paste0("Curva ROC - Modelo logístico final (AUC = ", round(auc_test, 3), ")"),
x = "1 - Especificidad",
y = "Sensibilidad"
) +
theme_minimal()
ggsave(
filename = file.path(dir_graficos, "roc_modelo_logistico_final.png"),
plot = graf_roc,
width = 10,
height = 7,
dpi = 300
)
#############################
# 20. TRADE-OFF SENSIBILIDAD / ESPECIFICIDAD SEGÚN UMBRAL
#############################
tabla_umbrales_larga <- tabla_umbrales %>%
dplyr::select(Umbral, Sensitivity, Specificity) %>%
tidyr::pivot_longer(
cols = c(Sensitivity, Specificity),
names_to = "Metrica",
values_to = "Valor"
)
umbral_youden <- mejor_youden$Umbral[1]
graf_tradeoff_lineas <- ggplot(tabla_umbrales_larga, aes(x = Umbral, y = Valor, color = Metrica)) +
geom_line(linewidth = 1.3) +
geom_vline(xintercept = 0.50, linetype = 2, linewidth = 0.7) +
geom_vline(xintercept = umbral_youden, linetype = 3, linewidth = 0.7) +
annotate("text", x = 0.50, y = 0.05, label = "0.50", angle = 90, vjust = -0.3) +
annotate("text", x = umbral_youden, y = 0.05, label = "Youden", angle = 90, vjust = -0.3) +
labs(
title = "Trade-off entre sensibilidad y especificidad",
subtitle = "Evaluación del efecto del punto de corte sobre el desempeño del modelo logístico",
x = "Umbral de clasificación",
y = "Valor de la métrica",
color = "Métrica"
) +
theme_minimal()
ggsave(
filename = file.path(dir_graficos, "tradeoff_sens_specificity_modelo_final.png"),
plot = graf_tradeoff_lineas,
width = 12,
height = 7,
dpi = 300
)
# Scatter trade-off por umbral
graf_tradeoff_scatter <- ggplot(
tabla_umbrales,
aes(
x = Specificity,
y = Sensitivity,
color = Umbral,
size = Accuracy
)
) +
geom_point(alpha = 0.9) +
geom_point(
data = subset(tabla_umbrales, abs(Umbral - 0.50) < 1e-9),
shape = 21,
fill = "white",
color = "black",
stroke = 1.2,
size = 5
) +
labs(
title = "Trade-off entre sensibilidad y especificidad según el umbral",
subtitle = "Cada punto representa un punto de corte distinto",
x = "Especificidad",
y = "Sensibilidad",
color = "Umbral",
size = "Accuracy"
) +
theme_minimal()
ggsave(
filename = file.path(dir_graficos, "scatter_tradeoff_umbral_modelo_final.png"),
plot = graf_tradeoff_scatter,
width = 12,
height = 7,
dpi = 300
)
#############################
# 21. PERFILES DE PREDICCIÓN
#############################
# Esta sección genera perfiles de bajo y alto riesgo de forma automática
# basándose en el signo de los coeficientes del modelo final.
# Si deseas perfiles institucionales más específicos, puedes editar manualmente
# los valores del data.frame perfil_bajo_manual y perfil_alto_manual.
obtener_perfiles_riesgo <- function(modelo, data_train) {
vars_finales <- all.vars(formula(modelo))[-1]
coefs <- coef(modelo)
resumen_perfiles <- list()
perfil_bajo <- list()
perfil_alto <- list()
for (v in vars_finales) {
if (is.factor(data_train[[v]])) {
niveles <- levels(data_train[[v]])
ref <- niveles[1]
coef_names <- names(coefs)[grepl(paste0("^", v), names(coefs))]
coef_vals <- coefs[coef_names]
niveles_no_ref <- sub(paste0("^", v), "", coef_names)
nivel_bajo <- ref
nivel_alto <- ref
if (length(coef_vals) > 0) {
if (min(coef_vals, na.rm = TRUE) < 0) {
nivel_bajo <- niveles_no_ref[which.min(coef_vals)]
}
if (max(coef_vals, na.rm = TRUE) > 0) {
nivel_alto <- niveles_no_ref[which.max(coef_vals)]
}
}
perfil_bajo[[v]] <- factor(nivel_bajo, levels = niveles)
perfil_alto[[v]] <- factor(nivel_alto, levels = niveles)
resumen_perfiles[[v]] <- data.frame(
Variable = v,
Referencia = ref,
Perfil_Bajo = nivel_bajo,
Perfil_Alto = nivel_alto,
stringsAsFactors = FALSE
)
} else if (is.numeric(data_train[[v]])) {
q1 <- as.numeric(quantile(data_train[[v]], 0.25, na.rm = TRUE))
q3 <- as.numeric(quantile(data_train[[v]], 0.75, na.rm = TRUE))
beta <- ifelse(v %in% names(coefs), coefs[v], 0)
if (beta >= 0) {
perfil_bajo[[v]] <- q1
perfil_alto[[v]] <- q3
} else {
perfil_bajo[[v]] <- q3
perfil_alto[[v]] <- q1
}
resumen_perfiles[[v]] <- data.frame(
Variable = v,
Referencia = "No aplica",
Perfil_Bajo = perfil_bajo[[v]],
Perfil_Alto = perfil_alto[[v]],
stringsAsFactors = FALSE
)
}
}
list(
tabla = do.call(rbind, resumen_perfiles),
perfil_bajo = as.data.frame(perfil_bajo, stringsAsFactors = FALSE),
perfil_alto = as.data.frame(perfil_alto, stringsAsFactors = FALSE)
)
}
perfiles_auto <- obtener_perfiles_riesgo(modelo_logit_final, train_model)
prob_bajo <- predict(modelo_logit_final, newdata = perfiles_auto$perfil_bajo, type = "response")
prob_alto <- predict(modelo_logit_final, newdata = perfiles_auto$perfil_alto, type = "response")
tabla_predicciones <- perfiles_auto$tabla
tabla_predicciones$Probabilidad_Desercion_Perfil_Bajo <- c(rep(NA, nrow(tabla_predicciones) - 1), prob_bajo)
tabla_predicciones$Probabilidad_Desercion_Perfil_Alto <- c(rep(NA, nrow(tabla_predicciones) - 1), prob_alto)
#############################
# 22. TABLAS RESUMEN PARA LA TESIS
#############################
tabla_resumen_modelo_final <- data.frame(
Umbral = 0.50,
AUC_Train = auc_train,
AUC_Test = auc_test,
Accuracy = tabla_metricas_test_050$Accuracy,
Kappa = tabla_metricas_test_050$Kappa,
Sensitivity = tabla_metricas_test_050$Sensitivity,
Specificity = tabla_metricas_test_050$Specificity,
Precision = tabla_metricas_test_050$Precision,
F1 = tabla_metricas_test_050$F1,
Balanced_Accuracy = tabla_metricas_test_050$Balanced_Accuracy
)
tabla_supuestos_modelo_logit <- data.frame(
Criterio = c(
"Multicolinealidad (GVIF/VIF)",
"Linealidad en el logit para EDAD",
"Hosmer-Lemeshow",
"Residuos de Pearson",
"Distancias de Cook"
),
Resultado = c(
"Sin evidencia de multicolinealidad severa",
ifelse(is.null(tabla_box_tidwell), "No evaluado",
ifelse(any(grepl("I\\(EDAD \\* log\\(EDAD\\)\\)", tabla_box_tidwell$Variable)) &&
tabla_box_tidwell$p_value[grepl("I\\(EDAD \\* log\\(EDAD\\)\\)", tabla_box_tidwell$Variable)][1] > 0.05,
"Linealidad razonable", "Revisar posible no linealidad")),
paste0("Chi-cuadrado = ", round(tabla_hl$Chi_square, 4),
"; gl = ", tabla_hl$gl,
"; p = ", round(tabla_hl$p_value, 4)),
paste0("Máx = ", round(tabla_residuos_pearson$Max, 4)),
paste0("Máx = ", round(tabla_cook$Max, 4),
"; 4/n = ", round(umbral_cook, 4),
"; casos > 4/n = ", nrow(casos_influyentes %>% dplyr::filter(Cook > umbral_cook)))
),
Conclusion = c(
"Cumple",
ifelse(is.null(tabla_box_tidwell), "No evaluado",
ifelse(any(grepl("I\\(EDAD \\* log\\(EDAD\\)\\)", tabla_box_tidwell$Variable)) &&
tabla_box_tidwell$p_value[grepl("I\\(EDAD \\* log\\(EDAD\\)\\)", tabla_box_tidwell$Variable)][1] > 0.05,
"Cumple", "Revisar")),
ifelse(tabla_hl$p_value > 0.05, "Buen ajuste", "Revisar"),
"Sin patrón general de mal ajuste",
"Sin influencia extrema global"
),
stringsAsFactors = FALSE
)
#############################
# 23. EXPORTAR A EXCEL
#############################
wb <- openxlsx::createWorkbook()
openxlsx::addWorksheet(wb, "Stepwise")
openxlsx::writeData(wb, "Stepwise", tabla_stepwise_resumen)
openxlsx::addWorksheet(wb, "Wald_IC95")
openxlsx::writeData(wb, "Wald_IC95", tabla_wald)
openxlsx::addWorksheet(wb, "Odds_Ratio")
openxlsx::writeData(wb, "Odds_Ratio", tabla_or)
openxlsx::addWorksheet(wb, "Resumen_Modelo_Final")
openxlsx::writeData(wb, "Resumen_Modelo_Final", tabla_resumen_modelo_final)
openxlsx::addWorksheet(wb, "Matriz_Confusion_050")
openxlsx::writeData(wb, "Matriz_Confusion_050", as.data.frame(matriz_confusion_050$table))
openxlsx::addWorksheet(wb, "LR_Test")
openxlsx::writeData(wb, "LR_Test", tabla_lr)
openxlsx::addWorksheet(wb, "Hosmer_Lemeshow")
openxlsx::writeData(wb, "Hosmer_Lemeshow", tabla_hl)
openxlsx::addWorksheet(wb, "AUC_Train_Test")
openxlsx::writeData(wb, "AUC_Train_Test", tabla_auc)
openxlsx::addWorksheet(wb, "Decision_Umbral")
openxlsx::writeData(wb, "Decision_Umbral", tabla_decision_umbral)
openxlsx::addWorksheet(wb, "VIF_GVIF")
openxlsx::writeData(wb, "VIF_GVIF", tabla_vif)
if (!is.null(tabla_box_tidwell)) {
openxlsx::addWorksheet(wb, "Box_Tidwell")
openxlsx::writeData(wb, "Box_Tidwell", tabla_box_tidwell)
}
openxlsx::addWorksheet(wb, "Pearson")
openxlsx::writeData(wb, "Pearson", tabla_residuos_pearson)
openxlsx::addWorksheet(wb, "Cook")
openxlsx::writeData(wb, "Cook", tabla_cook)
openxlsx::addWorksheet(wb, "Casos_Influyentes")
openxlsx::writeData(wb, "Casos_Influyentes", casos_influyentes)
openxlsx::addWorksheet(wb, "Predicciones")
openxlsx::writeData(wb, "Predicciones", tabla_predicciones)
openxlsx::addWorksheet(wb, "Supuestos")
openxlsx::writeData(wb, "Supuestos", tabla_supuestos_modelo_logit)
openxlsx::saveWorkbook(
wb,
file = file.path(dir_resultados, "modelo_logistico_final_resultados.xlsx"),
overwrite = TRUE
)
#############################
# 24. EXPORTAR CSV
#############################
write.csv(tabla_stepwise_resumen, file.path(dir_resultados, "tabla_stepwise_resumen.csv"), row.names = FALSE)
write.csv(tabla_wald, file.path(dir_resultados, "tabla_wald_modelo_final.csv"), row.names = FALSE)
write.csv(tabla_or, file.path(dir_resultados, "tabla_or_modelo_final.csv"), row.names = FALSE)
write.csv(tabla_resumen_modelo_final, file.path(dir_resultados, "tabla_resumen_modelo_final.csv"), row.names = FALSE)
write.csv(tabla_lr, file.path(dir_resultados, "tabla_lr_modelo_final.csv"), row.names = FALSE)
write.csv(tabla_hl, file.path(dir_resultados, "tabla_hosmer_lemeshow_modelo_final.csv"), row.names = FALSE)
write.csv(tabla_auc, file.path(dir_resultados, "tabla_auc_train_test_logit.csv"), row.names = FALSE)
write.csv(tabla_decision_umbral, file.path(dir_resultados, "tabla_decision_umbral_logit.csv"), row.names = FALSE)
write.csv(tabla_vif, file.path(dir_resultados, "tabla_vif_modelo_logistico_final.csv"), row.names = FALSE)
if (!is.null(tabla_box_tidwell)) {
write.csv(tabla_box_tidwell, file.path(dir_resultados, "tabla_box_tidwell_edad.csv"), row.names = FALSE)
}
write.csv(tabla_residuos_pearson, file.path(dir_resultados, "tabla_residuos_pearson.csv"), row.names = FALSE)
write.csv(tabla_cook, file.path(dir_resultados, "tabla_distancia_cook.csv"), row.names = FALSE)
write.csv(casos_influyentes, file.path(dir_resultados, "casos_influyentes_modelo_logit.csv"), row.names = FALSE)
write.csv(tabla_predicciones, file.path(dir_resultados, "tabla_predicciones_perfiles.csv"), row.names = FALSE)
write.csv(tabla_supuestos_modelo_logit, file.path(dir_resultados, "tabla_supuestos_modelo_logit.csv"), row.names = FALSE)
#############################
# 25. MENSAJES FINALES
#############################
cat("\n=========================================\n")##
## =========================================cat("PROCESO FINALIZADO CORRECTAMENTE\n")## PROCESO FINALIZADO CORRECTAMENTEcat("=========================================\n")## =========================================cat("Modelo final:\n")## Modelo final:print(formula(modelo_logit_final))## Y_LOGIT ~ CARRERA + TIPO_INGRESO + RENDIMIENTO + TRABAJAcat("\nAUC train:", round(auc_train, 4), "\n")##
## AUC train: 0.827cat("AUC test :", round(auc_test, 4), "\n")## AUC test : 0.7888cat("Hosmer-Lemeshow p-value:", round(tabla_hl$p_value, 4), "\n")## Hosmer-Lemeshow p-value: 0.5419cat("Umbral óptimo Youden:", round(mejor_youden$Umbral[1], 2), "\n")## Umbral óptimo Youden: 0.75cat("Resultados guardados en:", dir_resultados, "\n")## Resultados guardados en: resultados_modelo_logistico_finalcat("Gráficos guardados en:", dir_graficos, "\n")## Gráficos guardados en: graficos_modelo_logistico_finalcat("=========================================\n")## =========================================############################################################
# 16. FIGURA SCATTER CORREGIDA: COMPARACIÓN INTEGRADA
# Sensibilidad, especificidad y AUC en conjunto de prueba
############################################################
library(dplyr)
library(ggplot2)
library(ggrepel)
library(scales)
# 16.1 Preparar datos del gráfico
datos_scatter <- tabla_metricas_test %>%
mutate(
Modelo = factor(
Modelo,
levels = c("NB", "Tree", "KNN", "SVM", "NN", "Logit", "RF")
),
Tipo = case_when(
Modelo %in% c("NB", "SVM") ~ "Sesgado",
Modelo %in% c("Tree", "KNN") ~ "Intermedio",
Modelo %in% c("NN", "Logit", "RF") ~ "Óptimo",
TRUE ~ "Intermedio"
),
Tipo = factor(Tipo, levels = c("Sesgado", "Intermedio", "Óptimo"))
)
# 16.2 Diagnóstico de valores faltantes
cat("\n=============================\n")##
## =============================cat("DIAGNÓSTICO PARA FIGURA SCATTER\n")## DIAGNÓSTICO PARA FIGURA SCATTERcat("=============================\n")## =============================print(datos_scatter %>% select(Modelo, Specificity, Sensitivity, AUC, Tipo))## Modelo Specificity Sensitivity AUC Tipo
## Logit Logit 0.5758 0.7883 0.7796 Óptimo
## Tree Tree 0.6667 0.7664 0.7853 Intermedio
## RF RF 0.2121 0.9635 0.7920 Óptimo
## SVM SVM 0.2879 0.9270 0.7773 Sesgado
## KNN KNN 0.2424 0.8905 0.7131 Intermedio
## NN NN 0.4848 0.7956 0.7460 Óptimo
## NB NB 0.6970 0.7445 0.7890 Sesgadomodelos_con_na <- datos_scatter %>%
filter(is.na(Specificity) | is.na(Sensitivity) | is.na(AUC))
if (nrow(modelos_con_na) > 0) {
cat("\nATENCIÓN: modelos con métricas NA:\n")
print(modelos_con_na)
}
# 16.3 Usar solo modelos con métricas completas para graficar
datos_scatter_plot <- datos_scatter %>%
filter(
!is.na(Specificity),
!is.na(Sensitivity),
!is.na(AUC)
)
# 16.4 Colores solicitados
colores_tipo <- c(
"Sesgado" = "#D7191C",
"Intermedio" = "#FFD700",
"Óptimo" = "#1A9641"
)
# 16.5 Límites dinámicos para evitar pérdida de nodos
x_min <- max(0, min(datos_scatter_plot$Specificity, na.rm = TRUE) - 0.08)
x_max <- min(1, max(datos_scatter_plot$Specificity, na.rm = TRUE) + 0.08)
y_min <- max(0, min(datos_scatter_plot$Sensitivity, na.rm = TRUE) - 0.06)
y_max <- min(1, max(datos_scatter_plot$Sensitivity, na.rm = TRUE) + 0.06)
# 16.6 Gráfico corregido
Figura_scatter_modelo_color <- ggplot(
datos_scatter_plot,
aes(
x = Specificity,
y = Sensitivity
)
) +
geom_point(
aes(size = AUC, fill = Tipo),
shape = 21,
color = "black",
stroke = 0.7,
alpha = 0.92
) +
ggrepel::geom_text_repel(
aes(label = Modelo, color = Tipo),
size = 5.4,
fontface = "bold",
box.padding = 0.50,
point.padding = 0.35,
segment.color = "grey55",
segment.linewidth = 0.35,
max.overlaps = Inf,
show.legend = FALSE
) +
scale_fill_manual(
values = colores_tipo,
name = "Tipo de comportamiento"
) +
scale_color_manual(
values = colores_tipo,
guide = "none"
) +
scale_size_continuous(
name = "AUC",
range = c(7, 15),
labels = number_format(accuracy = 0.01)
) +
scale_x_continuous(
breaks = seq(0, 1, by = 0.10),
labels = number_format(accuracy = 0.01)
) +
scale_y_continuous(
breaks = seq(0, 1, by = 0.05),
labels = number_format(accuracy = 0.01)
) +
coord_cartesian(
xlim = c(x_min, x_max),
ylim = c(y_min, y_max),
clip = "off"
) +
labs(
title = "Comparación integrada del desempeño de los modelos",
subtitle = "Sensibilidad, especificidad y AUC en el conjunto de prueba",
x = "Especificidad",
y = "Sensibilidad"
) +
theme_minimal(base_size = 15) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 22),
plot.subtitle = element_text(hjust = 0.5, size = 15),
axis.title = element_text(face = "bold", size = 17),
axis.text = element_text(size = 13),
legend.position = "right",
legend.title = element_text(face = "bold", size = 13),
legend.text = element_text(size = 12),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(linewidth = 0.30, colour = "grey85"),
plot.margin = ggplot2::margin(15, 35, 15, 15)
) +
guides(
fill = guide_legend(
order = 1,
override.aes = list(
shape = 21,
size = 5,
colour = "black",
alpha = 1
)
),
size = guide_legend(order = 2)
)## Warning in ggrepel::geom_text_repel(aes(label = Modelo, color = Tipo), size =
## 5.4, : Ignoring unknown parameters: `segment.linewidth`print(Figura_scatter_modelo_color)
ggsave(
filename = "salidas_tesis/Figura_scatter_modelo_color_corregida.png",
plot = Figura_scatter_modelo_color,
width = 12.5,
height = 7.5,
dpi = 320
)############################################################
# RESUMEN DESCRIPTIVO DE LA VARIABLE EDAD
############################################################
# Librería necesaria
library(readxl)
# Ajustar esta ruta según tu equipo
ruta_datos <- "desercion.xlsx"
# Importar base de datos
tesis <- read_excel(ruta_datos)
tesis <- as.data.frame(tesis, stringsAsFactors = FALSE)
# Seleccionar variable EDAD
edad <- tesis$EDAD
# Construcción de tabla descriptiva
tabla_edad <- data.frame(
Variables = "Edad",
`Número de Casos` = sum(!is.na(edad)),
Minimo = min(edad, na.rm = TRUE),
Q1 = quantile(edad, 0.25, na.rm = TRUE),
Mediana = median(edad, na.rm = TRUE),
Media = mean(edad, na.rm = TRUE),
Q3 = quantile(edad, 0.75, na.rm = TRUE),
Maximo = max(edad, na.rm = TRUE)
)
# Redondear resultados
tabla_edad[, -1] <- round(tabla_edad[, -1], 2)
# Mostrar tabla
print(tabla_edad)## Variables Número.de.Casos Minimo Q1 Mediana Media Q3 Maximo
## 25% Edad 682 18 20 23 27.11 32 70############################################################
# TABLA 4.10
# Comparación Train vs Test
############################################################
library(dplyr)
library(tidyr)
library(openxlsx)
############################################################
# FUNCIÓN AUXILIAR
############################################################
comparar_variable <- function(var, train_data, test_data){
# Variable numérica
if(is.numeric(train_data[[var]])){
salida <- data.frame(
Caracteristica = c(
paste0(var, " (media)"),
paste0(var, " (mediana)"),
paste0(var, " (DE)"),
paste0(var, " (mínimo)"),
paste0(var, " (máximo)")
),
Base_Entrenamiento = c(
round(mean(train_data[[var]], na.rm = TRUE), 2),
round(median(train_data[[var]], na.rm = TRUE), 2),
round(sd(train_data[[var]], na.rm = TRUE), 2),
round(min(train_data[[var]], na.rm = TRUE), 2),
round(max(train_data[[var]], na.rm = TRUE), 2)
),
Base_Prueba = c(
round(mean(test_data[[var]], na.rm = TRUE), 2),
round(median(test_data[[var]], na.rm = TRUE), 2),
round(sd(test_data[[var]], na.rm = TRUE), 2),
round(min(test_data[[var]], na.rm = TRUE), 2),
round(max(test_data[[var]], na.rm = TRUE), 2)
),
Observacion = c(
"Valores similares",
"Consistencia en tendencia central",
"Valores similares",
"",
""
)
)
return(salida)
}
# Variable categórica
else{
train_tab <- prop.table(table(train_data[[var]]))*100
test_tab <- prop.table(table(test_data[[var]]))*100
niveles <- union(names(train_tab), names(test_tab))
salida <- data.frame()
for(niv in niveles){
n_train <- sum(train_data[[var]] == niv, na.rm = TRUE)
p_train <- round(train_tab[niv], 2)
n_test <- sum(test_data[[var]] == niv, na.rm = TRUE)
p_test <- round(test_tab[niv], 2)
diferencia <- abs(p_train - p_test)
observacion <- case_when(
diferencia < 2 ~ "Distribución prácticamente idéntica",
diferencia < 5 ~ "Proporción consistente",
diferencia < 10 ~ "Diferencia leve",
TRUE ~ "Diferencia moderada"
)
fila <- data.frame(
Caracteristica = paste0(var, " (", niv, ")"),
Base_Entrenamiento =
paste0(n_train, " (", p_train, "%)"),
Base_Prueba =
paste0(n_test, " (", p_test, "%)"),
Observacion = observacion
)
salida <- rbind(salida, fila)
}
return(salida)
}
}
############################################################
# VARIABLES A ANALIZAR
############################################################
variables_modelo <- c(
"ESTADO_ACADEMICO",
"EDAD",
"SEXO",
"RESIDENCIA",
"EST_CIVIL",
"CARRERA",
"TIPO_INGRESO",
"MODALIDAD",
"RENDIMIENTO",
"TIPO_COL",
"TRABAJA",
"SOLVENTAR",
"ESTUDIOS_PADRES",
"INGRESO",
"NIVEL_SOCIO"
)
############################################################
# TABLA GENERAL
############################################################
tabla_410 <- data.frame()
############################################################
# RESUMEN GENERAL
############################################################
resumen_general <- data.frame(
Caracteristica = c(
"Tamaño de muestra",
"Desertores",
"No desertores"
),
Base_Entrenamiento = c(
paste0(
round(nrow(train_data)/(nrow(train_data)+nrow(test_data))*100,2),
"% (",
nrow(train_data),
" obs.)"
),
paste0(
round(mean(train_data$ESTADO_ACADEMICO=="DESERTOR")*100,2),
"%"
),
paste0(
round(mean(train_data$ESTADO_ACADEMICO=="NO_DESERTOR")*100,2),
"%"
)
),
Base_Prueba = c(
paste0(
round(nrow(test_data)/(nrow(train_data)+nrow(test_data))*100,2),
"% (",
nrow(test_data),
" obs.)"
),
paste0(
round(mean(test_data$ESTADO_ACADEMICO=="DESERTOR")*100,2),
"%"
),
paste0(
round(mean(test_data$ESTADO_ACADEMICO=="NO_DESERTOR")*100,2),
"%"
)
),
Observacion = c(
"Proporción definida en la partición",
"Distribución prácticamente idéntica",
"Distribución prácticamente idéntica"
)
)
tabla_410 <- rbind(tabla_410, resumen_general)
############################################################
# VARIABLES RESTANTES
############################################################
for(v in variables_modelo[-1]){
temp <- comparar_variable(v, train_data, test_data)
tabla_410 <- rbind(tabla_410, temp)
}
############################################################
# VISUALIZAR
############################################################
print(tabla_410)## Caracteristica Base_Entrenamiento
## 1 Tamaño de muestra 70.23% (479 obs.)
## 2 Desertores 0%
## 3 No desertores 0%
## 4 EDAD (media) 26.84
## 5 EDAD (mediana) 23
## 6 EDAD (DE) 9.03
## 7 EDAD (mínimo) 18
## 8 EDAD (máximo) 70
## 9 SEXO (FEMENINO) 281 (58.66%)
## 10 SEXO (MASCULINO) 198 (41.34%)
## 11 RESIDENCIA (ASUNCIÓN) 159 (33.19%)
## 12 RESIDENCIA (CENTRAL) 240 (50.1%)
## 13 RESIDENCIA (RESTO DEL PAÍS) 80 (16.7%)
## 14 EST_CIVIL (Soltero/a) 385 (80.38%)
## 15 EST_CIVIL (Casado/a) 79 (16.49%)
## 16 EST_CIVIL (Divorciado/a) 11 (2.3%)
## 17 EST_CIVIL (Otro) 4 (0.84%)
## 18 CARRERA (ESTADISTICA-PRES) 130 (27.14%)
## 19 CARRERA (ESTADISTICA-SEMI) 73 (15.24%)
## 20 CARRERA (MATEMATICA) 106 (22.13%)
## 21 CARRERA (EDUCACION MATEMATICA-PRES) 48 (10.02%)
## 22 CARRERA (EDUCACION MATEMATICA-SEMI) 122 (25.47%)
## 23 TIPO_INGRESO (INGRESO) 454 (94.78%)
## 24 TIPO_INGRESO (TRASLADO) 0 (0%)
## 25 TIPO_INGRESO (ADMISION DIRECTA) 25 (5.22%)
## 26 MODALIDAD (Presencial) 284 (59.29%)
## 27 MODALIDAD (Semipresencial) 195 (40.71%)
## 28 RENDIMIENTO (APROBADO) 285 (59.5%)
## 29 RENDIMIENTO (NO APROBADO) 194 (40.5%)
## 30 TIPO_COL (Público) 320 (66.81%)
## 31 TIPO_COL (Subvencionado) 39 (8.14%)
## 32 TIPO_COL (Privado) 120 (25.05%)
## 33 TRABAJA (SI) 261 (54.49%)
## 34 TRABAJA (NO) 218 (45.51%)
## 35 SOLVENTAR (Beca/exoneración total) 28 (5.85%)
## 36 SOLVENTAR (Beca/exoneración parcial) 35 (7.31%)
## 37 SOLVENTAR (Trabajo Personal) 242 (50.52%)
## 38 SOLVENTAR (Ayuda Familiar) 174 (36.33%)
## 39 ESTUDIOS_PADRES (hasta 13 años) 103 (21.5%)
## 40 ESTUDIOS_PADRES (14-23 años) 154 (32.15%)
## 41 ESTUDIOS_PADRES (24-29 años) 111 (23.17%)
## 42 ESTUDIOS_PADRES (30-34 años) 75 (15.66%)
## 43 ESTUDIOS_PADRES (Más de 34 años) 36 (7.52%)
## 44 INGRESO (Hasta dos salarios mínimos) 300 (62.63%)
## 45 INGRESO (Más de dos y hasta cinco salarios mínimos) 138 (28.81%)
## 46 INGRESO (Más de cinco y hasta diez salarios mínimos) 29 (6.05%)
## 47 INGRESO (Más de diez y hasta quince salarios mínimos) 8 (1.67%)
## 48 INGRESO (Más de quince salarios mínimos) 4 (0.84%)
## 49 NIVEL_SOCIO (BAJO) 41 (8.56%)
## 50 NIVEL_SOCIO (MEDIO) 438 (91.44%)
## Base_Prueba Observacion
## 1 29.77% (203 obs.) Proporción definida en la partición
## 2 0% Distribución prácticamente idéntica
## 3 0% Distribución prácticamente idéntica
## 4 27.74 Valores similares
## 5 24 Consistencia en tendencia central
## 6 9.47 Valores similares
## 7 18
## 8 61
## 9 114 (56.16%) Proporción consistente
## 10 89 (43.84%) Proporción consistente
## 11 70 (34.48%) Distribución prácticamente idéntica
## 12 96 (47.29%) Proporción consistente
## 13 37 (18.23%) Distribución prácticamente idéntica
## 14 159 (78.33%) Proporción consistente
## 15 38 (18.72%) Proporción consistente
## 16 3 (1.48%) Distribución prácticamente idéntica
## 17 3 (1.48%) Distribución prácticamente idéntica
## 18 52 (25.62%) Distribución prácticamente idéntica
## 19 27 (13.3%) Distribución prácticamente idéntica
## 20 44 (21.67%) Distribución prácticamente idéntica
## 21 22 (10.84%) Distribución prácticamente idéntica
## 22 58 (28.57%) Proporción consistente
## 23 188 (92.61%) Proporción consistente
## 24 0 (0%) Distribución prácticamente idéntica
## 25 15 (7.39%) Proporción consistente
## 26 118 (58.13%) Distribución prácticamente idéntica
## 27 85 (41.87%) Distribución prácticamente idéntica
## 28 124 (61.08%) Distribución prácticamente idéntica
## 29 79 (38.92%) Distribución prácticamente idéntica
## 30 145 (71.43%) Proporción consistente
## 31 10 (4.93%) Proporción consistente
## 32 48 (23.65%) Distribución prácticamente idéntica
## 33 111 (54.68%) Distribución prácticamente idéntica
## 34 92 (45.32%) Distribución prácticamente idéntica
## 35 17 (8.37%) Proporción consistente
## 36 14 (6.9%) Distribución prácticamente idéntica
## 37 104 (51.23%) Distribución prácticamente idéntica
## 38 68 (33.5%) Proporción consistente
## 39 54 (26.6%) Diferencia leve
## 40 61 (30.05%) Proporción consistente
## 41 53 (26.11%) Proporción consistente
## 42 26 (12.81%) Proporción consistente
## 43 9 (4.43%) Proporción consistente
## 44 130 (64.04%) Distribución prácticamente idéntica
## 45 52 (25.62%) Proporción consistente
## 46 16 (7.88%) Distribución prácticamente idéntica
## 47 4 (1.97%) Distribución prácticamente idéntica
## 48 1 (0.49%) Distribución prácticamente idéntica
## 49 24 (11.82%) Proporción consistente
## 50 179 (88.18%) Proporción consistente############################################################
# EXPORTAR
############################################################
if(!dir.exists("salidas_tesis")){
dir.create("salidas_tesis")
}
write.csv(
tabla_410,
"salidas_tesis/Tabla_4_10_Train_vs_Test.csv",
row.names = FALSE
)
############################################################
# EXPORTAR A EXCEL
############################################################
wb <- createWorkbook()
addWorksheet(wb, "Tabla_4_10")
writeData(
wb,
sheet = "Tabla_4_10",
x = tabla_410
)
saveWorkbook(
wb,
"salidas_tesis/Tabla_4_10_Train_vs_Test.xlsx",
overwrite = TRUE
)
cat("\n====================================\n")##
## ====================================cat("TABLA 4.10 GENERADA CORRECTAMENTE\n")## TABLA 4.10 GENERADA CORRECTAMENTEcat("====================================\n")## ====================================############################################################
# INSTALAR PAQUETES NECESARIOS
############################################################
paquetes <- c(
"partykit",
"gridExtra",
"RColorBrewer",
"rpart.plot"
)
instalar <- paquetes[!(paquetes %in% installed.packages()[, "Package"])]
if(length(instalar) > 0){
install.packages(instalar)
}
############################################################
# CARGAR LIBRERÍAS
############################################################
library(rpart)
library(partykit)## Cargando paquete requerido: libcoin## Cargando paquete requerido: mvtnorm##
## Adjuntando el paquete: 'partykit'## The following object is masked from 'package:flextable':
##
## widthlibrary(grid)
library(gridExtra)
library(RColorBrewer)
library(rpart.plot)############################################################
# ÁRBOL DE DECISIÓN FINAL - FORMATO ACADÉMICO MEJORADO
############################################################
############################################################
# 1. LIBRERÍAS
############################################################
paquetes <- c(
"rpart",
"rpart.plot"
)
instalar <- paquetes[!(paquetes %in% installed.packages()[, "Package"])]
if(length(instalar) > 0){
install.packages(instalar)
}
library(rpart)
library(rpart.plot)
############################################################
# 2. AJUSTE DEL MODELO
############################################################
set.seed(123)
modelo_tree_final <- rpart(
ESTADO_ACADEMICO ~ .,
data = train_data,
method = "class",
parms = list(split = "gini"),
control = rpart.control(
cp = 0.01,
minsplit = 20,
maxdepth = 4
)
)
############################################################
# 3. CREAR DIRECTORIO DE SALIDA
############################################################
if(!dir.exists("salidas_tesis")){
dir.create("salidas_tesis")
}
############################################################
# 4. EXPORTAR FIGURA EN ALTA CALIDAD
############################################################
png(
filename = "salidas_tesis/Arbol_decision_final_mejorado.png",
width = 3400,
height = 2200,
res = 320
)
rpart.plot(
modelo_tree_final,
##########################################################
# ESTRUCTURA GENERAL
##########################################################
type = 2,
extra = 104,
under = TRUE,
fallen.leaves = TRUE,
compress = TRUE,
uniform = TRUE,
##########################################################
# COLORES
##########################################################
box.palette = c("#8EC5F8", "#CDE8B6"),
shadow.col = "gray85",
border.col = "gray40",
branch.col = "gray45",
branch.lwd = 1.8,
##########################################################
# TEXTO
##########################################################
faclen = 0,
varlen = 14,
cex = 0.78,
split.cex = 0.72,
##########################################################
# FORMATO DE NODOS
##########################################################
nn = FALSE,
roundint = FALSE,
##########################################################
# TÍTULOS
##########################################################
main = "Árbol de decisión final",
sub = "Clasificación de estudiantes desertores y no desertores"
)
dev.off()## png
## 2############################################################
# 5. MOSTRAR EN PANTALLA
############################################################
rpart.plot(
modelo_tree_final,
##########################################################
# ESTRUCTURA GENERAL
##########################################################
type = 2,
extra = 104,
under = TRUE,
fallen.leaves = TRUE,
compress = TRUE,
uniform = TRUE,
##########################################################
# COLORES
##########################################################
box.palette = c("#8EC5F8", "#CDE8B6"),
shadow.col = "gray85",
border.col = "gray40",
branch.col = "gray45",
branch.lwd = 1.8,
##########################################################
# TEXTO
##########################################################
faclen = 0,
varlen = 14,
cex = 0.78,
split.cex = 0.72,
##########################################################
# FORMATO DE NODOS
##########################################################
nn = FALSE,
roundint = FALSE,
##########################################################
# TÍTULOS
##########################################################
main = "Árbol de decisión final",
sub = "Clasificación de estudiantes desertores y no desertores"
)
############################################################
# 6. MENSAJE FINAL
############################################################
cat("\n====================================\n")##
## ====================================cat("ÁRBOL EXPORTADO CORRECTAMENTE\n")## ÁRBOL EXPORTADO CORRECTAMENTEcat("====================================\n")## ====================================cat("Ubicación:\n")## Ubicación:cat("salidas_tesis/Arbol_decision_final_mejorado.png\n")## salidas_tesis/Arbol_decision_final_mejorado.pngcat("====================================\n")## ====================================############################################################
# 4.6.3 OBSERVACIONES INFLUYENTES
# DISTANCIA DE COOK
############################################################
library(ggplot2)
library(dplyr)
############################################################
# 1. DISTANCIA DE COOK
############################################################
cook <- cooks.distance(modelo_logit_final)
datos_cook <- data.frame(
Observacion = 1:length(cook),
Cook = cook
)
############################################################
# 2. UMBRAL DE REFERENCIA
############################################################
umbral_cook <- 4 / nrow(train_model)
############################################################
# 3. IDENTIFICAR OBSERVACIONES INFLUYENTES
############################################################
obs_influyentes <- datos_cook %>%
filter(Cook > umbral_cook)
cat("\n====================================\n")##
## ====================================cat("OBSERVACIONES POTENCIALMENTE INFLUYENTES\n")## OBSERVACIONES POTENCIALMENTE INFLUYENTEScat("====================================\n")## ====================================print(obs_influyentes)## Observacion Cook
## 9 5 0.01023796
## 30 17 0.01023796
## 32 19 0.01023796
## 41 24 0.01023796
## 144 93 0.02632527
## 147 94 0.01023796
## 190 127 0.01023796
## 221 151 0.01023796
## 231 158 0.01023796
## 263 175 0.01023796
## 300 207 0.01923223
## 323 224 0.03298445
## 332 230 0.01447493
## 352 245 0.01447493
## 477 334 0.01403820
## 493 346 0.03126254
## 496 349 0.01403820
## 501 353 0.01953395
## 507 357 0.01933387
## 527 370 0.01447493
## 567 393 0.01271334
## 635 442 0.01536198
## 641 448 0.01953395
## 648 454 0.01271334
## 653 457 0.01953395
## 659 462 0.01953395
## 671 471 0.01447493
## 682 479 0.01403820############################################################
# 4. GRÁFICO DISTANCIA DE COOK
############################################################
grafico_cook <- ggplot(datos_cook,
aes(x = Observacion,
y = Cook)) +
geom_col(fill = "#2c7fb8",
alpha = 0.85) +
geom_hline(yintercept = umbral_cook,
color = "red",
linewidth = 1,
linetype = "dashed") +
labs(
title = "Distancia de Cook del modelo logístico final",
subtitle = paste0(
"Línea roja: umbral de referencia = ",
round(umbral_cook, 4)
),
x = "Observaciones",
y = "Distancia de Cook"
) +
theme_minimal(base_size = 15) +
theme(
plot.title = element_text(
face = "bold",
hjust = 0.5,
size = 20
),
plot.subtitle = element_text(
hjust = 0.5,
size = 13
),
axis.title = element_text(
face = "bold",
size = 15
)
)
print(grafico_cook)
############################################################
# 5. EXPORTAR GRÁFICO
############################################################
ggsave(
filename = "salidas_tesis/grafico_distancia_cook.png",
plot = grafico_cook,
width = 11,
height = 6,
dpi = 320
)
############################################################
# 6. EXPORTAR TABLA
############################################################
write.csv(
obs_influyentes,
"salidas_tesis/observaciones_influyentes_cook.csv",
row.names = FALSE
)############################################################
# 4.6.4 CURVA DE CALIBRACIÓN
############################################################
library(ggplot2)
library(dplyr)
############################################################
# 1. PROBABILIDADES PREDICHAS
############################################################
prob_pred <- predict(
modelo_logit_final,
newdata = test_data,
type = "response"
)
############################################################
# 2. OBTENER RESPUESTA REAL DESDE EL MODELO
############################################################
respuesta_real <- model.response(
model.frame(modelo_logit_final)
)
############################################################
# 3. AJUSTAR LONGITUD
############################################################
# SOLO TOMAR OBSERVACIONES TEST
############################################################
n_test <- length(prob_pred)
respuesta_test <- tail(
respuesta_real,
n_test
)
############################################################
# 4. CONVERTIR A BINARIA
############################################################
observado_binario <- ifelse(
as.character(respuesta_test) %in%
c("Desertor", "1"),
1,
0
)
############################################################
# 5. CREAR DATAFRAME
############################################################
datos_calibracion <- data.frame(
Observado = observado_binario,
Probabilidad = prob_pred
)
############################################################
# 6. CREAR DECILES
############################################################
datos_calibracion$Decil <- cut(
datos_calibracion$Probabilidad,
breaks = quantile(
datos_calibracion$Probabilidad,
probs = seq(0, 1, 0.10),
na.rm = TRUE
),
include.lowest = TRUE
)
############################################################
# 7. TABLA CALIBRACIÓN
############################################################
tabla_calibracion <- datos_calibracion %>%
group_by(Decil) %>%
summarise(
Probabilidad_Predicha = mean(Probabilidad),
Probabilidad_Observada = mean(Observado),
n = n(),
.groups = "drop"
)
print(tabla_calibracion)## # A tibble: 10 × 4
## Decil Probabilidad_Predicha Probabilidad_Observada n
## <fct> <dbl> <dbl> <int>
## 1 [0.0526,0.31] 0.237 0.8 25
## 2 (0.31,0.402] 0.388 0.619 21
## 3 (0.402,0.496] 0.496 0.64 25
## 4 (0.496,0.602] 0.535 0.7 10
## 5 (0.602,0.683] 0.666 0.541 37
## 6 (0.683,0.715] 0.715 0.333 6
## 7 (0.715,0.901] 0.879 0.5 20
## 8 (0.901,0.948] 0.935 0.75 20
## 9 (0.948,0.967] 0.964 0.724 29
## 10 (0.967,0.971] 0.971 0.6 10############################################################
# 8. GRÁFICO
############################################################
grafico_calibracion <- ggplot(
tabla_calibracion,
aes(
x = Probabilidad_Predicha,
y = Probabilidad_Observada
)
) +
geom_point(
size = 4,
color = "#2166ac"
) +
geom_line(
color = "#2166ac",
linewidth = 1
) +
geom_abline(
intercept = 0,
slope = 1,
color = "red",
linetype = "dashed",
linewidth = 1
) +
labs(
title = "Curva de calibración del modelo logístico final",
subtitle = "Probabilidades observadas vs predichas",
x = "Probabilidad predicha",
y = "Probabilidad observada"
) +
coord_equal() +
theme_minimal(base_size = 15)
print(grafico_calibracion)
############################################################
# EXPORTAR GRÁFICO DE CALIBRACIÓN
############################################################
ggsave(
filename = "salidas_tesis/curva_calibracion_modelo_logistico_final.png",
plot = grafico_calibracion,
width = 10,
height = 7,
dpi = 320
)
############################################################
# EXPORTAR EN PDF (OPCIONAL)
############################################################
ggsave(
filename = "salidas_tesis/curva_calibracion_modelo_logistico_final.pdf",
plot = grafico_calibracion,
width = 10,
height = 7
)
cat("\n====================================\n")##
## ====================================cat("GRÁFICO DE CALIBRACIÓN EXPORTADO\n")## GRÁFICO DE CALIBRACIÓN EXPORTADOcat("====================================\n")## ====================================cat("PNG : curva_calibracion_modelo_logistico_final.png\n")## PNG : curva_calibracion_modelo_logistico_final.pngcat("PDF : curva_calibracion_modelo_logistico_final.pdf\n")## PDF : curva_calibracion_modelo_logistico_final.pdf# ==========================================
# MODELO LOGÍSTICO FINAL
# ==========================================
modelo_final <- glm(
ESTADO_ACADEMICO~ CARRERA + TIPO_INGRESO + RENDIMIENTO + TRABAJA,
data = train_data,
family = binomial(link = "logit")
)
summary(modelo_final)##
## Call:
## glm(formula = ESTADO_ACADEMICO ~ CARRERA + TIPO_INGRESO + RENDIMIENTO +
## TRABAJA, family = binomial(link = "logit"), data = train_data)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -0.2346 0.2838 -0.827
## CARRERAESTADISTICA-SEMI -0.8195 0.4998 -1.640
## CARRERAMATEMATICA 0.3153 0.3108 1.014
## CARRERAEDUCACION MATEMATICA-PRES 1.0093 0.3927 2.570
## CARRERAEDUCACION MATEMATICA-SEMI -0.3861 0.3485 -1.108
## TIPO_INGRESOADMISION DIRECTA 1.8464 0.6157 2.999
## RENDIMIENTONO APROBADO -2.7083 0.3480 -7.782
## TRABAJANO 0.2546 0.2681 0.950
## Pr(>|z|)
## (Intercept) 0.40843
## CARRERAESTADISTICA-SEMI 0.10105
## CARRERAMATEMATICA 0.31036
## CARRERAEDUCACION MATEMATICA-PRES 0.01016 *
## CARRERAEDUCACION MATEMATICA-SEMI 0.26787
## TIPO_INGRESOADMISION DIRECTA 0.00271 **
## RENDIMIENTONO APROBADO 0.00000000000000714 ***
## TRABAJANO 0.34227
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 606.02 on 478 degrees of freedom
## Residual deviance: 451.52 on 471 degrees of freedom
## AIC: 467.52
##
## Number of Fisher Scoring iterations: 5# ==========================================
# ECUACIÓN DEL MODELO
# ==========================================
coeficientes <- coef(modelo_final)
coeficientes## (Intercept) CARRERAESTADISTICA-SEMI
## -0.2346206 -0.8194958
## CARRERAMATEMATICA CARRERAEDUCACION MATEMATICA-PRES
## 0.3153176 1.0093170
## CARRERAEDUCACION MATEMATICA-SEMI TIPO_INGRESOADMISION DIRECTA
## -0.3860986 1.8464468
## RENDIMIENTONO APROBADO TRABAJANO
## -2.7082810 0.2546231# Mostrar ecuación logística
cat("\nECUACIÓN LOGÍSTICA:\n")##
## ECUACIÓN LOGÍSTICA:cat("log(p/(1-p)) = ",
round(coeficientes[1],4),
ifelse(coeficientes[2] >=0, " + ", " "),
round(coeficientes[2],4),"*Carrera_EstadisticaSemi",
ifelse(coeficientes[3] >=0, " + ", " "),
round(coeficientes[3],4),"*Carrera_Matematica",
ifelse(coeficientes[4] >=0, " + ", " "),
round(coeficientes[4],4),"*Carrera_EducMatPres",
ifelse(coeficientes[5] >=0, " + ", " "),
round(coeficientes[5],4),"*Carrera_EducMatSemi",
ifelse(coeficientes[6] >=0, " + ", " "),
round(coeficientes[6],4),"*AdmisionDirecta",
ifelse(coeficientes[7] >=0, " + ", " "),
round(coeficientes[7],4),"*NoAprobado",
ifelse(coeficientes[8] >=0, " + ", " "),
round(coeficientes[8],4),"*NoTrabaja"
)## log(p/(1-p)) = -0.2346 -0.8195 *Carrera_EstadisticaSemi + 0.3153 *Carrera_Matematica + 1.0093 *Carrera_EducMatPres -0.3861 *Carrera_EducMatSemi + 1.8464 *AdmisionDirecta -2.7083 *NoAprobado + 0.2546 *NoTrabaja# ==========================================
# PROBABILIDADES PREDICHAS
# ==========================================
train_data$prob_predicha <- predict(
modelo_final,
type = "response"
)
head(train_data$prob_predicha)## [1] 0.50500045 0.06366857 0.73678393 0.44161244 0.73678393 0.06738413library(ggplot2)
ggplot(train_data,
aes(x = prob_predicha,
fill = factor(ESTADO_ACADEMICO))) +
geom_density(alpha = 0.4) +
labs(
title = "Distribución de probabilidades predichas",
subtitle = "Modelo de regresión logística final",
x = "Probabilidad predicha de deserción",
y = "Densidad",
fill = "Estado Académico"
) +
theme_minimal()