Una organización busca entender qué factores explican que un empleado decida rotar (dejar su cargo actual) y anticipar, para cada trabajador, la probabilidad de que esto ocurra en el próximo período. Esta información le permitiría a la gerencia intervenir de manera preventiva sobre los empleados con mayor riesgo, en lugar de reaccionar una vez la renuncia ya ocurrió.
El objetivo de este informe es construir un modelo de regresión logística binomial que estime la probabilidad de rotación de un empleado a partir de un conjunto de variables organizacionales, identificar cuáles de esas variables presentan una asociación estadísticamente significativa con la probabilidad de rotación y estimar la magnitud y dirección de dicha asociación, y evaluar qué tan bien discrimina el modelo entre empleados que rotan y empleados que permanecen. Al tratarse de un diseño observacional y transversal, el alcance de estas conclusiones es asociativo, no causal — vuelvo sobre esta distinción de forma explícita al interpretar los coeficientes.
La variable respuesta es Rotación, codificada como “Sí”
(el empleado rotó) y “No” (el empleado permaneció). Al ser una variable
binaria, el modelo de regresión lineal múltiple no es apropiado: no
garantiza que las predicciones queden dentro del rango \([0,1]\), viola los supuestos de normalidad
y homoscedasticidad de los residuales, y su ajuste (medido por \(R^2\)) suele ser artificialmente bajo
cuando la variable dependiente solo toma dos valores. Por esa razón se
recurre al modelo logit, que transforma el problema de estimar una
probabilidad en el de estimar el logaritmo de una razón de
probabilidades (log-odds), garantizando que la probabilidad estimada
esté siempre acotada entre 0 y 1.
Aunque el tamaño muestral (1.470 empleados) permitiría en principio
considerar un conjunto amplio de predictores, opté por un modelo
parsimonioso de seis variables por dos razones concretas, no solo por
preferencia de simplicidad. Primero, una regla empírica estándar en
regresión logística (Peduzzi et al., 1996) recomienda al menos 10
eventos de la clase minoritaria por cada variable candidata para que las
estimaciones de máxima verosimilitud sean estables; con 237 casos de
rotación (la clase minoritaria), el límite razonable de predictores está
alrededor de 23-24 — es decir, usar las 24 variables disponibles
simultáneamente dejaría el modelo prácticamente en el borde de esa
regla, sin margen para interacciones ni términos adicionales. Segundo,
varias variables del conjunto de datos miden esencialmente el mismo
constructo (por ejemplo, Antigüedad,
Antigüedad_Cargo, Años_acargo_con_mismo_jefe y
Años_ultima_promoción son todas proxies del tiempo que un
empleado lleva vinculado a la empresa, y están correlacionadas entre sí
por encima de 0.55), por lo que incluirlas todas añadiría redundancia
sin aportar información nueva proporcional al costo en grados de
libertad e interpretabilidad.
Antes de seleccionar, es necesario clasificar correctamente el tipo
de cada variable, porque varias aparecen codificadas como enteros sin
ser cuantitativas en sentido estricto: Educación,
Satisfacción_Ambiental, Satisfación_Laboral,
Rendimiento_Laboral y Equilibrio_Trabajo_Vida
son escalas ordinales tipo Likert (1 a 4 o 1 a 5). En este informe las
trato como variables cuasi-cuantitativas de intervalo aproximado, una
convención común en la práctica aplicada, y lo declaro explícitamente
para que la decisión sea trazable.
Seleccioné 3 variables categóricas y 3 cuantitativas, con un criterio doble: (a) sustento teórico en la literatura de comportamiento organizacional y rotación de personal, y (b) baja colinealidad entre sí, verificada empíricamente antes de comprometerme con el conjunto final.
Horas_Extra (Sí / No). La teoría del
agotamiento de recursos en el trabajo (job demands-resources
model; Maslach, describiendo el síndrome de burnout) plantea que
exigencias laborales sostenidas por encima de la capacidad de
recuperación del trabajador erosionan su compromiso y aumentan la
intención de abandonar el cargo. Hipótesis: trabajar horas extra
aumenta la probabilidad de rotación.
Viaje de Negocios (No_Viaja / Raramente
/ Frecuentemente). El modelo de conflicto trabajo-familia de Greenhaus y
Beutell (1985) explica que cuando las demandas del rol laboral
interfieren con las del rol familiar o personal, aumenta la tensión y,
con ella, la probabilidad de abandonar el empleo. Viajar con frecuencia
es una de las formas más directas de ese conflicto. Hipótesis: se
plantea una tendencia ordinal creciente — a mayor frecuencia de viaje,
mayor probabilidad de rotación — a contrastar empíricamente, aunque el
modelo la trata como variable categórica (cada nivel comparado contra la
referencia) y no impone una tendencia lineal entre niveles.
Estado_Civil (Casado / Divorciado /
Soltero). El modelo clásico de rotación de Mobley (1977) y buena parte
de la literatura posterior asocian los menores costos de salida (menor
cantidad de ataduras económicas y familiares) con una mayor probabilidad
de dejar un empleo. Los solteros, en promedio, enfrentan menores costos
de oportunidad al cambiar de trabajo. Hipótesis: los solteros rotan
más que los casados o divorciados.
Satisfación_Laboral (escala 1-4). Es
uno de los predictores más replicados en la literatura de rotación; el
metaanálisis de Griffeth, Hom y Gaertner (2000) la identifica como uno
de los correlatos más consistentes de la intención de rotación.
Hipótesis: a mayor satisfacción laboral, menor probabilidad de
rotación (signo negativo esperado).
Ingreso_Mensual. Desde la teoría de
intercambio económico, el salario es la contraprestación directa más
visible de la relación laboral. Hipótesis: a mayor ingreso, menor
probabilidad de rotación.
Antigüedad_Cargo (años en el cargo
actual). La literatura describe dos mecanismos en tensión: por un lado,
mayor tiempo en el cargo puede generar arraigo y mayores costos de
salida (menor rotación); por otro, puede reflejar estancamiento
profesional sin promoción (mayor rotación). Planteo la hipótesis como
abierta a contrastar empíricamente con el signo del coeficiente, en
lugar de asumir una dirección a priori.
Antes de comprometerme con el modelo final, verifiqué que las tres
variables cuantitativas elegidas no estuvieran midiendo esencialmente lo
mismo. La correlación de Pearson entre ellas resultó prácticamente nula
entre Satisfación_Laboral y las otras dos (\(r \approx 0\)), y moderada-baja entre
Ingreso_Mensual y Antigüedad_Cargo (\(r \approx 0.36\)). No existe un umbral
universal único a partir del cual una correlación se vuelve
“preocupante” — la verificación definitiva la hago más adelante con el
Factor de Inflación de Varianza (VIF) sobre el modelo ya ajustado, que
sí tiene un criterio de referencia más establecido; esta revisión
bivariada previa solo sirve como primer filtro razonable antes de
comprometerme con el conjunto final de variables.
# Cargo las librerías necesarias para el análisis completo:
# tidyverse: manipulación de datos y graficación (incluye dplyr, ggplot2, readr)
# janitor: estandarización de nombres de columnas
# broom: extracción ordenada (tidy) de resultados de modelos
# car: cálculo del factor de inflación de varianza (VIF)
# caret: partición estratificada de datos y matriz de confusión
# pROC: curva ROC y área bajo la curva (AUC)
# kableExtra: tablas HTML con formato profesional
# scales: formato de ejes y etiquetas (porcentajes, separadores de miles)
# rstatix: cálculo de Cramér's V (tamaño del efecto para tablas de contingencia)
# effsize: cálculo de la d de Cohen (tamaño del efecto para diferencias de medias)
library(tidyverse)
library(janitor)
library(broom)
library(car)
library(caret)
library(pROC)
library(kableExtra)
library(scales)
library(rstatix)
library(effsize)
# Fijo una semilla para que la partición train/test y cualquier proceso
# aleatorio sean exactamente reproducibles por cualquier persona que
# vuelva a ejecutar este documento.
set.seed(123)
# Instalo (solo si no está ya instalado) y cargo paqueteMODELOS, que
# distribuye la base de datos "rotacion" directamente desde el
# repositorio de GitHub del curso. Esto evita depender de un archivo
# local y garantiza que trabajo con la misma versión de los datos
# que el resto del curso.
if (!requireNamespace("paqueteMODELOS", quietly = TRUE)) {
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("centromagis/paqueteMODELOS", force = TRUE)
}
library(paqueteMODELOS)
data("rotacion")
# clean_names() estandariza los nombres de columna: los pasa a minúsculas,
# reemplaza espacios por guiones bajos y elimina tildes/diéresis, lo que
# evita tener que escribir nombres entre comillas invertidas (backticks)
# en el resto del documento.
rotacion <- rotacion %>% clean_names()
# Reviso la estructura general del conjunto de datos.
glimpse(rotacion)
## Rows: 1,470
## Columns: 24
## $ rotacion <chr> "Si", "No", "Si", "No", "No", "No", "No", …
## $ edad <dbl> 41, 49, 37, 33, 27, 32, 59, 30, 38, 36, 35…
## $ viaje_de_negocios <chr> "Raramente", "Frecuentemente", "Raramente"…
## $ departamento <chr> "Ventas", "IyD", "IyD", "IyD", "IyD", "IyD…
## $ distancia_casa <dbl> 1, 8, 2, 3, 2, 2, 3, 24, 23, 27, 16, 15, 2…
## $ educacion <dbl> 2, 1, 2, 4, 1, 2, 3, 1, 3, 3, 3, 2, 1, 2, …
## $ campo_educacion <chr> "Ciencias", "Ciencias", "Otra", "Ciencias"…
## $ satisfaccion_ambiental <dbl> 2, 3, 4, 4, 1, 4, 3, 4, 4, 3, 1, 4, 1, 2, …
## $ genero <chr> "F", "M", "M", "F", "M", "M", "F", "M", "M…
## $ cargo <chr> "Ejecutivo_Ventas", "Investigador_Cientifi…
## $ satisfacion_laboral <dbl> 4, 2, 3, 3, 2, 4, 1, 3, 3, 3, 2, 3, 3, 4, …
## $ estado_civil <chr> "Soltero", "Casado", "Soltero", "Casado", …
## $ ingreso_mensual <dbl> 5993, 5130, 2090, 2909, 3468, 3068, 2670, …
## $ trabajos_anteriores <dbl> 8, 1, 6, 1, 9, 0, 4, 1, 0, 6, 0, 0, 1, 0, …
## $ horas_extra <chr> "Si", "No", "Si", "Si", "No", "No", "Si", …
## $ porcentaje_aumento_salarial <dbl> 11, 23, 15, 11, 12, 13, 20, 22, 21, 13, 13…
## $ rendimiento_laboral <dbl> 3, 4, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3, 3, 3, …
## $ anos_experiencia <dbl> 8, 10, 7, 8, 6, 8, 12, 1, 10, 17, 6, 10, 5…
## $ capacitaciones <dbl> 0, 3, 3, 3, 3, 2, 3, 2, 2, 3, 5, 3, 1, 2, …
## $ equilibrio_trabajo_vida <dbl> 1, 3, 3, 3, 3, 2, 2, 3, 3, 2, 3, 3, 2, 3, …
## $ antiguedad <dbl> 6, 10, 0, 8, 2, 7, 1, 1, 9, 7, 5, 9, 5, 2,…
## $ antiguedad_cargo <dbl> 4, 7, 0, 7, 2, 7, 0, 0, 7, 7, 4, 5, 2, 2, …
## $ anos_ultima_promocion <dbl> 0, 1, 0, 3, 2, 3, 0, 0, 1, 7, 0, 0, 4, 1, …
## $ anos_acargo_con_mismo_jefe <dbl> 5, 7, 0, 0, 2, 6, 0, 0, 8, 7, 3, 8, 3, 2, …
# Cuento valores faltantes por columna. Un data frame de calidad para
# modelar debe estar libre de NA no gestionados, o al menos con una
# regla explícita de tratamiento si existen.
na_por_columna <- colSums(is.na(rotacion))
sum(na_por_columna)
## [1] 0
El conjunto de datos tiene 1.470 registros y 24 variables, sin valores faltantes en ninguna columna, por lo que no fue necesario imputar ni eliminar observaciones por esta razón. Esto simplifica la etapa de limpieza, pero no exime de revisar la calidad de las categorías (niveles mal escritos, categorías con muy pocas observaciones) antes de modelar.
# Recodifico las variables categóricas de interés como factores, fijando
# explícitamente el nivel de referencia (baseline) de cada una. Esta
# decisión no es trivial: el nivel de referencia es contra el que se
# interpretan todos los demás coeficientes del modelo.
#
# - horas_extra: "No" como referencia (el escenario sin sobrecarga).
# - viaje_de_negocios: "No_Viaja" como referencia (el escenario de
# menor conflicto trabajo-familia), en lugar del orden alfabético
# por defecto que dejaría "Frecuentemente" como base.
# - estado_civil: "Casado" como referencia, coincide con el orden
# alfabético pero lo dejo explícito para que la decisión sea legible.
# - rotacion: fijo "No" como referencia para que "Si" sea la categoría
# de éxito (Y = 1) en el modelo logístico.
rotacion <- rotacion %>%
mutate(
rotacion = factor(rotacion, levels = c("No", "Si")),
horas_extra = factor(horas_extra, levels = c("No", "Si")),
viaje_de_negocios = factor(viaje_de_negocios,
levels = c("No_Viaja", "Raramente", "Frecuentemente")),
estado_civil = factor(estado_civil,
levels = c("Casado", "Divorciado", "Soltero")),
# Reescalo el ingreso antes de modelar. La unidad original
# (ingreso_mensual tal como viene en el dataset) es demasiado
# pequeña para que un cambio de una unidad produzca un odds ratio
# distinguible de 1 — no porque el efecto no exista, sino porque
# la unidad no es relevante en esa escala. Verifico primero el
# rango real antes de decidir el factor de reescalamiento, en
# lugar de asumir una denominación (por ejemplo, pesos en
# millones) sin comprobarla.
ingreso_miles = ingreso_mensual / 1000
)
range(rotacion$ingreso_mensual)
## [1] 1009 19999
mean(rotacion$ingreso_mensual)
## [1] 6502.931
# Verifico que la recodificación no haya generado NA por desajuste
# de escritura entre los niveles originales y los que definí (por
# ejemplo, diferencias de mayúsculas/minúsculas). clean_names() solo
# estandariza los NOMBRES de columna, nunca los VALORES dentro de
# ellas, así que este chequeo es indispensable cada vez que se
# recodifica manualmente a partir de un factor().
recodificacion_check <- rotacion %>%
summarise(across(c(rotacion, horas_extra, viaje_de_negocios, estado_civil),
~ sum(is.na(.x))))
print(recodificacion_check)
## # A tibble: 1 × 4
## rotacion horas_extra viaje_de_negocios estado_civil
## <int> <int> <int> <int>
## 1 0 0 0 0
stopifnot(
"La recodificación generó NA: revisa que los niveles en factor(..., levels = ...) coincidan EXACTAMENTE (mayúsculas incluidas) con los valores originales del dataset" =
sum(recodificacion_check) == 0
)
Verifico explícitamente que la recodificación no haya introducido valores faltantes, en lugar de asumirlo: una discrepancia de mayúsculas/minúsculas o de escritura entre los niveles declarados y los valores reales del dataset convertiría silenciosamente esas observaciones en NA, lo que sesgaría todo el análisis posterior sin que se note a simple vista.
El range() y la media de ingreso_mensual
confirman que la variable se mueve en un orden de magnitud de miles (no
de millones): dividir entre 1.000, en lugar de entre 1.000.000, es la
escala que efectivamente produce una variación observable dentro del
rango real de los datos. De ahí la variable derivada
ingreso_miles, que es la que utilizo en el modelo en lugar
de ingreso_mensual en su unidad original.
# Tabla de frecuencias absolutas y relativas de la variable respuesta.
tabla_rotacion <- rotacion %>%
count(rotacion) %>%
mutate(porcentaje = round(n / sum(n) * 100, 1))
tabla_rotacion %>%
kable(col.names = c("Rotación", "Frecuencia", "Porcentaje (%)"),
align = "c", caption = "Distribución de la variable Rotación") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Rotación | Frecuencia | Porcentaje (%) |
|---|---|---|
| No | 1233 | 83.9 |
| Si | 237 | 16.1 |
ggplot(tabla_rotacion, aes(x = rotacion, y = n, fill = rotacion)) +
geom_col(width = 0.55) +
geom_text(aes(label = paste0(n, " (", porcentaje, "%)")),
vjust = -0.4, size = 4.2, fontface = "bold") +
scale_fill_manual(values = c("No" = "#2C7FB8", "Si" = "#D7191C")) +
scale_y_continuous(limits = c(0, 1400), expand = expansion(mult = c(0, 0.05))) +
labs(title = "Distribución de empleados según rotación",
x = NULL, y = "Número de empleados") +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
plot.title = element_text(face = "bold"))
El 16.1 % de los empleados (237 de 1.470) rotó en el período observado, mientras que el 83.9 % restante permaneció en su cargo. Este es un desbalance de clases moderado, no extremo, pero suficiente para que deba tenerlo en cuenta más adelante: un modelo que simplemente prediga “No rota” para todos los empleados ya acertaría en el 83.9 % de los casos sin haber aprendido nada útil, así que la exactitud (accuracy) por sí sola será una métrica engañosa para evaluar este modelo. Volveré sobre esto en la sección de evaluación.
# Construyo un gráfico de distribución para cada una de las tres
# variables categóricas seleccionadas, en un solo panel con facet_wrap
# para facilitar la comparación visual.
categoricas_long <- rotacion %>%
select(horas_extra, viaje_de_negocios, estado_civil) %>%
pivot_longer(everything(), names_to = "variable", values_to = "nivel") %>%
mutate(variable = dplyr::recode(variable,
horas_extra = "Horas extra",
viaje_de_negocios = "Viaje de negocios",
estado_civil = "Estado civil"))
ggplot(categoricas_long, aes(x = nivel, fill = variable)) +
geom_bar(width = 0.6) +
geom_text(stat = "count", aes(label = after_stat(count)),
vjust = -0.4, size = 3.8) +
facet_wrap(~ variable, scales = "free_x", ncol = 1) +
labs(title = "Distribución de las variables categóricas seleccionadas",
x = NULL, y = "Número de empleados") +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
strip.text = element_text(face = "bold", size = 12),
plot.title = element_text(face = "bold"))
Las tres variables tienen todas sus categorías bien pobladas (ninguna con menos de 100 observaciones), lo que reduce el riesgo de que el modelo logístico presente separación cuasi-perfecta (coeficientes que divergen a valores extremos por falta de datos en alguna combinación de categorías). La mayoría de los empleados no trabaja horas extra, viaja raramente y está casado; esto es relevante porque el nivel de referencia de cada variable en el modelo corresponde, en los tres casos, a la categoría más frecuente.
cuantitativas_long <- rotacion %>%
select(satisfacion_laboral, ingreso_mensual, antiguedad_cargo) %>%
pivot_longer(everything(), names_to = "variable", values_to = "valor") %>%
mutate(variable = dplyr::recode(variable,
satisfacion_laboral = "Satisfacción laboral (escala 1-4)",
ingreso_mensual = "Ingreso mensual",
antiguedad_cargo = "Antigüedad en el cargo (años)"))
ggplot(cuantitativas_long, aes(x = valor, fill = variable)) +
geom_histogram(bins = 20, color = "white") +
facet_wrap(~ variable, scales = "free", ncol = 1) +
labs(title = "Distribución de las variables cuantitativas seleccionadas",
x = NULL, y = "Frecuencia") +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
strip.text = element_text(face = "bold", size = 12),
plot.title = element_text(face = "bold"))
rotacion %>%
select(satisfacion_laboral, ingreso_mensual, antiguedad_cargo) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
group_by(Variable) %>%
summarise(Media = round(mean(valor), 1),
Mediana = round(median(valor), 1),
`Desv. estándar` = round(sd(valor), 1),
Mínimo = min(valor),
Máximo = max(valor)) %>%
mutate(Variable = dplyr::recode(Variable,
satisfacion_laboral = "Satisfacción laboral",
ingreso_mensual = "Ingreso mensual",
antiguedad_cargo = "Antigüedad en el cargo")) %>%
kable(align = "c", caption = "Estadísticos descriptivos de las variables cuantitativas") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Variable | Media | Mediana | Desv. estándar | Mínimo | Máximo |
|---|---|---|---|---|---|
| Antigüedad en el cargo | 4.2 | 3 | 3.6 | 0 | 18 |
| Ingreso mensual | 6502.9 | 4919 | 4708.0 | 1009 | 19999 |
| Satisfacción laboral | 2.7 | 3 | 1.1 | 1 | 4 |
Satisfación_Laboral se distribuye de forma prácticamente
uniforme entre sus cuatro niveles, sin un nivel dominante.
Ingreso_Mensual presenta un sesgo marcado hacia la derecha
(cola larga de salarios altos), típico de las distribuciones salariales
en organizaciones con jerarquías de cargos amplias.
Antigüedad_Cargo también está sesgada a la derecha,
concentrada en los primeros años del cargo, con una cola de empleados
que llevan mucho más tiempo en la misma posición.
Analizo ahora la relación de cada una de las seis variables
seleccionadas con Rotación, contrastando con la hipótesis
planteada en la sección 2.
# Función auxiliar que construye, para una variable categórica dada,
# la tabla de porcentaje de rotación por nivel.
tabla_pct_rotacion <- function(datos, variable) {
datos %>%
count({{ variable }}, rotacion) %>%
group_by({{ variable }}) %>%
mutate(porcentaje = round(n / sum(n) * 100, 1)) %>%
ungroup() %>%
filter(rotacion == "Si") %>%
select({{ variable }}, porcentaje)
}
pct_horas_extra <- tabla_pct_rotacion(rotacion, horas_extra)
pct_viaje <- tabla_pct_rotacion(rotacion, viaje_de_negocios)
pct_estado_civil <- tabla_pct_rotacion(rotacion, estado_civil)
pct_horas_extra %>%
kable(col.names = c("Horas extra", "% que rotó"), align = "c",
caption = "Porcentaje de rotación según horas extra") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Horas extra | % que rotó |
|---|---|
| No | 10.4 |
| Si | 30.5 |
pct_viaje %>%
kable(col.names = c("Viaje de negocios", "% que rotó"), align = "c",
caption = "Porcentaje de rotación según frecuencia de viaje") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Viaje de negocios | % que rotó |
|---|---|
| No_Viaja | 8.0 |
| Raramente | 15.0 |
| Frecuentemente | 24.9 |
pct_estado_civil %>%
kable(col.names = c("Estado civil", "% que rotó"), align = "c",
caption = "Porcentaje de rotación según estado civil") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Estado civil | % que rotó |
|---|---|
| Casado | 12.5 |
| Divorciado | 10.1 |
| Soltero | 25.5 |
bivariado_cat_long <- bind_rows(
pct_horas_extra %>% rename(nivel = horas_extra) %>% mutate(variable = "Horas extra"),
pct_viaje %>% rename(nivel = viaje_de_negocios) %>% mutate(variable = "Viaje de negocios"),
pct_estado_civil %>% rename(nivel = estado_civil) %>% mutate(variable = "Estado civil")
)
ggplot(bivariado_cat_long, aes(x = nivel, y = porcentaje, fill = variable)) +
geom_col(width = 0.55) +
geom_text(aes(label = paste0(porcentaje, "%")), vjust = -0.4, size = 4) +
facet_wrap(~ variable, scales = "free_x", ncol = 1) +
scale_y_continuous(limits = c(0, 40), expand = expansion(mult = c(0, 0.05))) +
labs(title = "Porcentaje de rotación por categoría",
x = NULL, y = "% de empleados que rotó") +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
strip.text = element_text(face = "bold", size = 12),
plot.title = element_text(face = "bold"))
# Aplico la prueba chi-cuadrado de independencia a cada variable
# categórica frente a Rotación. La hipótesis nula es que la variable
# categórica y la rotación son independientes; un valor p pequeño
# lleva a rechazarla, es decir, a concluir que sí existe asociación.
chi_horas <- chisq.test(table(rotacion$horas_extra, rotacion$rotacion))
chi_viaje <- chisq.test(table(rotacion$viaje_de_negocios, rotacion$rotacion))
chi_civil <- chisq.test(table(rotacion$estado_civil, rotacion$rotacion))
# Un valor p pequeño solo dice que la asociación es poco probable bajo
# independencia, no dice qué tan grande es esa asociación — con 1.470
# observaciones casi cualquier asociación real, por pequeña que sea,
# termina siendo "significativa". Por eso acompaño cada prueba con
# Cramér's V, que cuantifica la magnitud de la asociación en una
# escala de 0 (independencia total) a 1 (asociación perfecta),
# independiente del tamaño de la muestra.
v_horas <- rstatix::cramer_v(table(rotacion$horas_extra, rotacion$rotacion))
v_viaje <- rstatix::cramer_v(table(rotacion$viaje_de_negocios, rotacion$rotacion))
v_civil <- rstatix::cramer_v(table(rotacion$estado_civil, rotacion$rotacion))
tibble(
Variable = c("Horas extra", "Viaje de negocios", "Estado civil"),
`Estadístico chi-cuadrado` = round(c(chi_horas$statistic, chi_viaje$statistic, chi_civil$statistic), 2),
`Grados de libertad` = c(chi_horas$parameter, chi_viaje$parameter, chi_civil$parameter),
`Valor p` = signif(c(chi_horas$p.value, chi_viaje$p.value, chi_civil$p.value), 3),
`Cramér's V` = round(c(v_horas, v_viaje, v_civil), 3)
) %>%
kable(align = "c", caption = "Pruebas chi-cuadrado de independencia y tamaño del efecto (Cramér's V)") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Variable | Estadístico chi-cuadrado | Grados de libertad | Valor p | Cramér’s V |
|---|---|---|---|---|
| Horas extra | 87.56 | 1 | 0.0e+00 | 0.246 |
| Viaje de negocios | 24.18 | 2 | 5.6e-06 | 0.128 |
| Estado civil | 46.16 | 2 | 0.0e+00 | 0.177 |
Interpreto Cramér’s V con las convenciones habituales, que no son un estándar estricto sino una referencia orientativa: valores cercanos a 0.1 se consideran una asociación pequeña, alrededor de 0.3 moderada, y por encima de 0.5 grande.
Las tres variables categóricas muestran una asociación
estadísticamente significativa con la rotación de forma bivariada (valor
p muy por debajo de 0.05 en los tres casos, con tamaños del efecto que
reporto en la tabla de Cramér’s V más abajo). Horas_Extra
exhibe la diferencia más marcada: cerca del triple de rotación entre
quienes trabajan horas extra frente a quienes no.
Viaje_de_Negocios sigue una tendencia creciente consistente
con la hipótesis planteada: la proporción de rotación aumenta a medida
que aumenta la frecuencia de viaje. Estado_Civil es
consistente con la hipótesis de Mobley: los solteros rotan
considerablemente más que casados y divorciados, mientras que estos dos
últimos grupos son similares entre sí.
Es importante ser preciso sobre el alcance de esta evidencia: estas son asociaciones bivariadas, sin controlar por las demás variables. Una variable puede resultar significativa de forma aislada y dejar de serlo — o cambiar de magnitud — cuando se controla simultáneamente por las demás en el modelo conjunto de la sección siguiente. Retomo esta comparación explícitamente al interpretar los coeficientes del modelo multivariado.
cuant_bivariado_long <- rotacion %>%
select(rotacion, satisfacion_laboral, ingreso_mensual, antiguedad_cargo) %>%
pivot_longer(-rotacion, names_to = "variable", values_to = "valor") %>%
mutate(variable = dplyr::recode(variable,
satisfacion_laboral = "Satisfacción laboral",
ingreso_mensual = "Ingreso mensual",
antiguedad_cargo = "Antigüedad en el cargo"))
ggplot(cuant_bivariado_long, aes(x = rotacion, y = valor, fill = rotacion)) +
geom_boxplot(width = 0.5, outlier.alpha = 0.4) +
facet_wrap(~ variable, scales = "free_y", ncol = 1) +
scale_fill_manual(values = c("No" = "#2C7FB8", "Si" = "#D7191C")) +
labs(title = "Distribución de las variables cuantitativas según rotación",
x = "Rotación", y = NULL) +
theme_minimal(base_size = 13) +
theme(legend.position = "none",
strip.text = element_text(face = "bold", size = 12),
plot.title = element_text(face = "bold"))
# Aplico la prueba t de Welch (no asume varianzas iguales entre los dos
# grupos, una precaución razonable cuando los tamaños de grupo son muy
# distintos, como ocurre aquí: 1233 vs 237 empleados).
t_satisfaccion <- t.test(satisfacion_laboral ~ rotacion, data = rotacion)
t_ingreso <- t.test(ingreso_mensual ~ rotacion, data = rotacion)
t_antiguedad <- t.test(antiguedad_cargo ~ rotacion, data = rotacion)
# La d de Cohen cuantifica la magnitud de la diferencia de medias en
# desviaciones estándar conjuntas, independientemente del tamaño de la
# muestra. Con n = 1.470 casi cualquier diferencia real resulta
# "significativa" en el sentido del valor p; la d de Cohen permite
# distinguir una diferencia significativa pero trivial de una que
# además es sustantivamente relevante.
d_satisfaccion <- effsize::cohen.d(satisfacion_laboral ~ rotacion, data = rotacion)$estimate
d_ingreso <- effsize::cohen.d(ingreso_mensual ~ rotacion, data = rotacion)$estimate
d_antiguedad <- effsize::cohen.d(antiguedad_cargo ~ rotacion, data = rotacion)$estimate
tibble(
Variable = c("Satisfacción laboral", "Ingreso mensual", "Antigüedad en el cargo"),
`Media (No rota)` = round(c(t_satisfaccion$estimate[1], t_ingreso$estimate[1], t_antiguedad$estimate[1]), 1),
`Media (Sí rota)` = round(c(t_satisfaccion$estimate[2], t_ingreso$estimate[2], t_antiguedad$estimate[2]), 1),
`Valor p (t de Welch)` = signif(c(t_satisfaccion$p.value, t_ingreso$p.value, t_antiguedad$p.value), 3),
`d de Cohen` = round(abs(c(d_satisfaccion, d_ingreso, d_antiguedad)), 3)
) %>%
kable(align = "c", caption = "Prueba t de Welch y tamaño del efecto (d de Cohen) por grupo de rotación") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Variable | Media (No rota) | Media (Sí rota) | Valor p (t de Welch) | d de Cohen |
|---|---|---|---|---|
| Satisfacción laboral | 2.8 | 2.5 | 0.000105 | 0.283 |
| Ingreso mensual | 6832.7 | 4787.1 | 0.000000 | 0.440 |
| Antigüedad en el cargo | 4.5 | 2.9 | 0.000000 | 0.442 |
Interpreto la d de Cohen con las convenciones usuales (también orientativas, no un corte estricto): alrededor de 0.2 es un efecto pequeño, 0.5 moderado, 0.8 o más grande.
Satisfación_Laboral es una escala ordinal de 1 a 4 que
en este informe trato como cuasi-cuantitativa; la prueba t asume una
escala de intervalo, supuesto defendible pero no automático. Como
verificación de sensibilidad, repito la comparación con una prueba no
paramétrica que no depende de ese supuesto:
# Wilcoxon-Mann-Whitney no asume que la variable sea de intervalo ni
# que su distribución sea normal; solo requiere que sea al menos
# ordinal, lo que la hace más apropiada en sentido estricto para
# Satisfación_Laboral que la prueba t.
wilcox_satisfaccion <- wilcox.test(satisfacion_laboral ~ rotacion, data = rotacion)
tibble(
Prueba = c("t de Welch", "Wilcoxon-Mann-Whitney"),
`Valor p` = signif(c(t_satisfaccion$p.value, wilcox_satisfaccion$p.value), 3)
) %>%
kable(align = "c", caption = "Satisfacción laboral: comparación entre prueba paramétrica y no paramétrica") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Prueba | Valor p |
|---|---|
| t de Welch | 1.05e-04 |
| Wilcoxon-Mann-Whitney | 7.96e-05 |
Ambas pruebas coinciden en la conclusión (asociación estadísticamente
significativa), lo que da tranquilidad de que tratar
Satisfación_Laboral como cuasi-cuantitativa no está
distorsionando el resultado.
Las diferencias de medias son estadísticamente significativas en las
tres variables cuantitativas, con tamaños del efecto que van de pequeños
a moderados según la tabla anterior. Los empleados que rotaron tienen,
en promedio, menor satisfacción laboral, menor ingreso mensual y menor
antigüedad en el cargo que quienes permanecieron. Esto es consistente,
de forma bivariada, con la hipótesis de Satisfación_Laboral
e Ingreso_Mensual, y con una de las dos direcciones
teóricas que había dejado abiertas para Antigüedad_Cargo
(arraigo/costo de salida, en vez de estancamiento profesional) — aunque,
igual que con las categóricas, esta es evidencia bivariada; confirmo si
se sostiene al controlar por las demás variables en la siguiente
sección.
Antes de estimar el modelo separo los datos en un conjunto de entrenamiento (para ajustar los coeficientes) y uno de prueba (para evaluar qué tan bien generaliza el modelo a datos que no vio durante el ajuste). Evaluar el modelo sobre los mismos datos con los que se estimó tiende a sobreestimar su capacidad predictiva real.
Uso caret::createDataPartition() en lugar de un muestreo
aleatorio simple porque estratifica por la variable
respuesta: garantiza que la proporción de 16.1 % de rotación se
mantenga aproximadamente igual tanto en el conjunto de entrenamiento
como en el de prueba. Con un desbalance de clases como el de este
conjunto de datos, un muestreo aleatorio simple podría, por azar, dejar
muy pocos casos de rotación en el conjunto de prueba y volver inestable
la evaluación.
indice_entrenamiento <- createDataPartition(rotacion$rotacion, p = 0.70, list = FALSE)
train <- rotacion[indice_entrenamiento, ]
test <- rotacion[-indice_entrenamiento, ]
# Verifico que la proporción de rotación se haya mantenido similar en
# ambos conjuntos, confirmando que la estratificación funcionó.
tibble(
Conjunto = c("Entrenamiento", "Prueba"),
n = c(nrow(train), nrow(test)),
`% Rotación` = round(c(mean(train$rotacion == "Si"), mean(test$rotacion == "Si")) * 100, 1)
) %>%
kable(align = "c", caption = "Tamaño y balance de clases en cada conjunto") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Conjunto | n | % Rotación |
|---|---|---|
| Entrenamiento | 1030 | 16.1 |
| Prueba | 440 | 16.1 |
# Estimo el modelo logit por máxima verosimilitud con glm(), family =
# binomial(link = "logit"). El modelo incluye las tres variables
# categóricas y las tres cuantitativas seleccionadas y justificadas
# en la sección 2, estimado únicamente con los datos de entrenamiento.
modelo_logit <- glm(
rotacion ~ horas_extra + viaje_de_negocios + estado_civil +
satisfacion_laboral + ingreso_miles + antiguedad_cargo,
family = binomial(link = "logit"),
data = train
)
summary(modelo_logit)
##
## Call:
## glm(formula = rotacion ~ horas_extra + viaje_de_negocios + estado_civil +
## satisfacion_laboral + ingreso_miles + antiguedad_cargo, family = binomial(link = "logit"),
## data = train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.06355 0.54487 -3.787 0.000152 ***
## horas_extraSi 1.43869 0.19307 7.452 9.23e-14 ***
## viaje_de_negociosRaramente 1.13912 0.46136 2.469 0.013547 *
## viaje_de_negociosFrecuentemente 1.94738 0.48424 4.021 5.78e-05 ***
## estado_civilDivorciado -0.02691 0.27171 -0.099 0.921106
## estado_civilSoltero 1.09627 0.21194 5.172 2.31e-07 ***
## satisfacion_laboral -0.31939 0.08455 -3.778 0.000158 ***
## ingreso_miles -0.09290 0.02741 -3.389 0.000701 ***
## antiguedad_cargo -0.12126 0.03296 -3.679 0.000234 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 909.69 on 1029 degrees of freedom
## Residual deviance: 738.39 on 1021 degrees of freedom
## AIC: 756.39
##
## Number of Fisher Scoring iterations: 5
Los coeficientes del modelo logit están en la escala del logaritmo de odds, que no es directamente interpretable en términos prácticos. Los exponencio para obtener razones de probabilidad (odds ratios), acompañadas de su intervalo de confianza al 95 %.
# broom::tidy() con exponentiate = TRUE devuelve directamente los odds
# ratios y sus intervalos de confianza, en lugar de tener que aplicar
# exp() manualmente a cada coeficiente.
tabla_or <- tidy(modelo_logit, conf.int = TRUE, exponentiate = TRUE) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
tabla_or %>%
kable(col.names = c("Término", "Odds ratio", "Error estándar", "Estadístico z",
"Valor p", "IC 95% inferior", "IC 95% superior"),
align = "c", caption = "Odds ratios e intervalos de confianza del modelo logit") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Término | Odds ratio | Error estándar | Estadístico z | Valor p | IC 95% inferior | IC 95% superior |
|---|---|---|---|---|---|---|
| (Intercept) | 0.127 | 0.545 | -3.787 | 0.000 | 0.040 | 0.350 |
| horas_extraSi | 4.215 | 0.193 | 7.452 | 0.000 | 2.894 | 6.175 |
| viaje_de_negociosRaramente | 3.124 | 0.461 | 2.469 | 0.014 | 1.364 | 8.554 |
| viaje_de_negociosFrecuentemente | 7.010 | 0.484 | 4.021 | 0.000 | 2.906 | 19.925 |
| estado_civilDivorciado | 0.973 | 0.272 | -0.099 | 0.921 | 0.565 | 1.645 |
| estado_civilSoltero | 2.993 | 0.212 | 5.172 | 0.000 | 1.984 | 4.560 |
| satisfacion_laboral | 0.727 | 0.085 | -3.778 | 0.000 | 0.615 | 0.857 |
| ingreso_miles | 0.911 | 0.027 | -3.389 | 0.001 | 0.861 | 0.959 |
| antiguedad_cargo | 0.886 | 0.033 | -3.679 | 0.000 | 0.829 | 0.944 |
Interpreto cada odds ratio como el factor multiplicativo que aplica a los odds de rotar, manteniendo las demás variables constantes:
ingreso_miles en lugar del ingreso en su unidad original,
precisamente para que el odds ratio sea interpretable en una escala
económicamente relevante. El odds ratio reportado corresponde al cambio
en los odds de rotar por cada 1.000 unidades monetarias adicionales de
ingreso mensual, y es menor a 1: confirma el efecto protector
esperado.# car::vif() calcula el Factor de Inflación de Varianza. Para
# predictores categóricos con más de dos niveles se reporta el
# GVIF junto con GVIF^(1/(2*Df)), que es la cantidad comparable
# directamente con el umbral usual de referencia (valores por
# encima de 5, o de forma más conservadora de 2.5-3, sugieren
# preocupación).
vif(modelo_logit)
## GVIF Df GVIF^(1/(2*Df))
## horas_extra 1.038505 1 1.019071
## viaje_de_negocios 1.027648 2 1.006841
## estado_civil 1.045188 2 1.011110
## satisfacion_laboral 1.029280 1 1.014534
## ingreso_miles 1.082913 1 1.040631
## antiguedad_cargo 1.079651 1 1.039062
Ningún valor se acerca al umbral de preocupación, lo que confirma en el modelo ajustado la baja colinealidad que ya había verificado de forma bivariada antes de seleccionar las variables. Los valores de VIF/GVIF no evidencian problemas relevantes de multicolinealidad entre los predictores incluidos, lo que reduce la preocupación por inflación de las varianzas asociada a redundancia lineal entre variables. Esto no implica, por sí solo, que las pruebas de significancia individual sean válidas en un sentido absoluto — su confiabilidad depende también de la especificación correcta del modelo, la independencia de las observaciones y otros supuestos que discuto en la sección de limitaciones — pero sí descarta específicamente la multicolinealidad como explicación de algún coeficiente no significativo.
Antes de pasar a la significancia global, reviso si hay observaciones individuales que estén influyendo desproporcionadamente sobre el ajuste del modelo — un empleado atípico que, por sí solo, mueva de forma importante los coeficientes estimados.
# La distancia de Cook, adaptada a GLM, mide cuánto cambiarían los
# coeficientes estimados si se eliminara esa observación del ajuste.
# No hay un punto de corte universal; una referencia común (no un
# estándar estricto) es 4/n, que uso aquí solo como guía visual.
cooksd <- cooks.distance(modelo_logit)
umbral_cook <- 4 / nrow(train)
cooks_df <- tibble(indice = seq_along(cooksd), distancia_cook = cooksd) %>%
mutate(influyente = distancia_cook > umbral_cook)
ggplot(cooks_df, aes(x = indice, y = distancia_cook, color = influyente)) +
geom_point(size = 1.4, alpha = 0.7) +
geom_hline(yintercept = umbral_cook, linetype = "dashed", color = "grey40") +
scale_color_manual(values = c("FALSE" = "#2C7FB8", "TRUE" = "#D7191C")) +
labs(title = "Distancia de Cook por observación (conjunto de entrenamiento)",
x = "Índice de la observación", y = "Distancia de Cook",
color = paste0("Por encima\ndel umbral (", round(umbral_cook, 4), ")")) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
n_influyentes <- sum(cooks_df$influyente)
De las 1030 observaciones del conjunto de entrenamiento, 71 quedan por encima del umbral de referencia (4/n). Dado que ese umbral es solo orientativo — con una muestra de este tamaño es normal que un número pequeño de observaciones lo supere sin que eso implique necesariamente un problema — no las elimino automáticamente del modelo. Si ese número fuera una proporción grande del total, o si uno o dos puntos concentraran una distancia de Cook muy por encima del resto (órdenes de magnitud, no solo por encima del umbral), ahí sí valdría la pena inspeccionar esos registros puntuales antes de confiar en el modelo.
# Contraste global: H0: todos los coeficientes (excepto el intercepto)
# son cero, frente a Ha: al menos uno es distinto de cero. El
# estadístico de prueba es la diferencia entre la devianza nula y la
# devianza residual, que se distribuye chi-cuadrado con los grados de
# libertad correspondientes a la diferencia de parámetros estimados.
estadistico_chi2 <- modelo_logit$null.deviance - modelo_logit$deviance
gl <- modelo_logit$df.null - modelo_logit$df.residual
valor_p_global <- pchisq(estadistico_chi2, df = gl, lower.tail = FALSE)
tibble(
`Devianza nula` = round(modelo_logit$null.deviance, 2),
`Devianza residual` = round(modelo_logit$deviance, 2),
`Chi-cuadrado` = round(estadistico_chi2, 2),
`Grados de libertad` = gl,
`Valor p` = signif(valor_p_global, 3)
) %>%
kable(align = "c", caption = "Prueba de significancia global del modelo") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Devianza nula | Devianza residual | Chi-cuadrado | Grados de libertad | Valor p |
|---|---|---|---|---|
| 909.69 | 738.39 | 171.3 | 8 | 0 |
El valor p de la prueba global es prácticamente cero, por lo que rechazo \(H_0\): el conjunto de variables seleccionadas, en su totalidad, sí aporta capacidad explicativa significativa sobre la rotación, más allá de lo que aportaría un modelo con solo el intercepto.
# Genero las probabilidades predichas por el modelo sobre el conjunto
# de prueba (datos que el modelo no vio al momento de estimarse).
test <- test %>%
mutate(prob_predicha = predict(modelo_logit, newdata = test, type = "response"))
# Clasifico con el punto de corte convencional de 0.5 como primer
# acercamiento.
test <- test %>%
mutate(clase_predicha_05 = factor(if_else(prob_predicha > 0.5, "Si", "No"),
levels = c("No", "Si")))
matriz_05 <- confusionMatrix(test$clase_predicha_05, test$rotacion, positive = "Si")
matriz_05
## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 359 61
## Si 10 10
##
## Accuracy : 0.8386
## 95% CI : (0.8009, 0.8718)
## No Information Rate : 0.8386
## P-Value [Acc > NIR] : 0.5316
##
## Kappa : 0.1602
##
## Mcnemar's Test P-Value : 2.958e-09
##
## Sensitivity : 0.14085
## Specificity : 0.97290
## Pos Pred Value : 0.50000
## Neg Pred Value : 0.85476
## Prevalence : 0.16136
## Detection Rate : 0.02273
## Detection Prevalence : 0.04545
## Balanced Accuracy : 0.55687
##
## 'Positive' Class : Si
##
# Visualizo la matriz de confusión como un mapa de calor con las
# frecuencias impresas dentro de cada celda, evitando que el texto se
# empalme con las etiquetas de los ejes.
matriz_df <- as.data.frame(matriz_05$table)
ggplot(matriz_df, aes(x = Reference, y = Prediction, fill = Freq)) +
geom_tile(color = "white", linewidth = 1) +
geom_text(aes(label = Freq), color = "black", size = 7, fontface = "bold") +
scale_fill_gradient(low = "#DEEBF7", high = "#2C7FB8") +
labs(title = "Matriz de confusión (punto de corte = 0.5)",
x = "Valor observado", y = "Valor predicho", fill = "Frecuencia") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
Con el punto de corte convencional de 0.5, el modelo alcanza una exactitud alta, pero esa exactitud es engañosa por el desbalance de clases que señalé en la sección de análisis univariado. La sensibilidad (la capacidad de detectar correctamente a los empleados que sí van a rotar) resulta baja: el modelo tiende a clasificar como “No rota” a muchos empleados que en realidad sí rotaron, precisamente porque el punto de corte de 0.5 favorece a la clase mayoritaria. La especificidad, en cambio, es muy alta: el modelo es muy bueno identificando a quienes van a permanecer.
Para el objetivo de negocio de este problema (identificar preventivamente a empleados en riesgo de rotar), una sensibilidad baja es una limitación seria: significa dejar pasar sin intervención a una parte importante de los empleados que efectivamente van a renunciar.
# pROC::roc() construye la curva ROC comparando, para cada punto de
# corte posible, la tasa de verdaderos positivos (sensibilidad) contra
# la tasa de falsos positivos (1 - especificidad).
curva_roc <- roc(test$rotacion, test$prob_predicha, levels = c("No", "Si"), direction = "<")
auc_valor <- auc(curva_roc)
# ci.auc() calcula el intervalo de confianza del AUC por bootstrap
# (por defecto, 2000 remuestreos). Reportar el AUC sin su intervalo
# da la falsa impresión de que es un valor exacto y no una estimación
# sujeta a la variabilidad de la muestra de prueba.
ic_auc <- ci.auc(curva_roc)
ggroc(curva_roc, color = "#D7191C", linewidth = 1.1) +
geom_abline(slope = 1, intercept = 1, linetype = "dashed", color = "grey50") +
labs(title = paste0("Curva ROC del modelo (AUC = ", round(auc_valor, 4), ")"),
subtitle = paste0("IC 95% (bootstrap): [", round(ic_auc[1], 4), ", ", round(ic_auc[3], 4), "]"),
x = "Especificidad", y = "Sensibilidad") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
El modelo alcanza un AUC de 0.7704 (IC 95%: 0.7112 - 0.8296) sobre el conjunto de prueba. Según las guías convencionales de discriminación de modelos de clasificación, esto se interpreta como una capacidad discriminativa aceptable a buena: bien por encima de 0.5 (una clasificación equivalente al azar) pero por debajo del rango de 0.9+ que se consideraría excelente. El intervalo de confianza, calculado por remuestreo bootstrap sobre el conjunto de prueba, recuerda que esta estimación está sujeta a la variabilidad propia de evaluar el modelo sobre una única partición train/test — una validación cruzada repetida daría una estimación más estable del rango plausible del AUC, algo que señalo como limitación más adelante y no desarrollo aquí por alcance de esta actividad.
El punto de corte de 0.5 y el óptimo de Youden son solo dos puntos entre infinitos posibles. Antes de fijar uno, reviso cómo se comportan sensibilidad, especificidad y precisión a lo largo de una rejilla de umbrales, para entender el trade-off completo en vez de comparar únicamente dos valores.
# coords() con una rejilla explícita de umbrales permite leer cómo
# evolucionan las métricas de clasificación en todo el rango de
# cortes posibles, no solo en un punto.
grilla_cortes <- coords(
curva_roc,
x = seq(0.10, 0.60, by = 0.05),
input = "threshold",
ret = c("threshold", "sensitivity", "specificity", "ppv", "npv")
)
grilla_cortes %>%
mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
kable(col.names = c("Punto de corte", "Sensibilidad", "Especificidad", "Precisión (PPV)", "VPN"),
align = "c", caption = "Métricas de clasificación a lo largo de una rejilla de puntos de corte") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Punto de corte | Sensibilidad | Especificidad | Precisión (PPV) | VPN |
|---|---|---|---|---|
| 0.10 | 0.859 | 0.547 | 0.268 | 0.953 |
| 0.15 | 0.718 | 0.672 | 0.297 | 0.925 |
| 0.20 | 0.606 | 0.778 | 0.344 | 0.911 |
| 0.25 | 0.549 | 0.840 | 0.398 | 0.906 |
| 0.30 | 0.465 | 0.873 | 0.412 | 0.894 |
| 0.35 | 0.380 | 0.892 | 0.403 | 0.882 |
| 0.40 | 0.282 | 0.924 | 0.417 | 0.870 |
| 0.45 | 0.183 | 0.954 | 0.433 | 0.859 |
| 0.50 | 0.141 | 0.973 | 0.500 | 0.855 |
| 0.55 | 0.113 | 0.978 | 0.500 | 0.851 |
| 0.60 | 0.113 | 0.989 | 0.667 | 0.853 |
grilla_larga <- grilla_cortes %>%
select(threshold, sensitivity, specificity) %>%
pivot_longer(-threshold, names_to = "metrica", values_to = "valor") %>%
mutate(metrica = dplyr::recode(metrica, sensitivity = "Sensibilidad", specificity = "Especificidad"))
ggplot(grilla_larga, aes(x = threshold, y = valor, color = metrica)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey50") +
scale_color_manual(values = c("Sensibilidad" = "#D7191C", "Especificidad" = "#2C7FB8")) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(title = "Sensibilidad y especificidad según el punto de corte",
subtitle = "La línea vertical marca el corte convencional de 0.5",
x = "Punto de corte", y = "Valor de la métrica", color = NULL) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"), legend.position = "bottom")
La tabla y el gráfico muestran con claridad el trade-off: a medida que el punto de corte baja, la sensibilidad sube y la especificidad baja, de forma continua a lo largo de toda la rejilla, no solo entre dos valores puntuales. Esto es lo que permite elegir un corte de forma informada en vez de arbitraria.
Dentro de ese trade-off, el índice de Youden (que maximiza sensibilidad + especificidad − 1) es un criterio estándar, reproducible y defendible como punto de partida cuando no se dispone de una estimación explícita de los costos relativos de cada tipo de error.
corte_optimo <- coords(curva_roc, "best", best.method = "youden",
ret = c("threshold", "sensitivity", "specificity"))
corte_optimo %>%
kable(align = "c", caption = "Punto de corte según el índice de Youden") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| threshold | sensitivity | specificity |
|---|---|---|
| 0.1634485 | 0.7183099 | 0.7235772 |
Es importante ser preciso sobre lo que Youden sí y no optimiza: maximiza la suma de sensibilidad y especificidad por igual, pero no incorpora los costos reales y potencialmente distintos de un falso negativo (un empleado que rota y el modelo no detecta) frente a un falso positivo (un empleado que se interviene preventivamente sin que fuera a rotar). Si la organización pudiera cuantificar esos costos — por ejemplo, el costo de reemplazar a un empleado clave frente al costo de una conversación o revisión salarial preventiva —, el punto de corte óptimo en términos de negocio podría ser distinto (probablemente más bajo que el de Youden, dado que el costo de un falso negativo suele ser mayor en problemas de retención de talento). La tabla de la rejilla anterior es precisamente la herramienta que permitiría a la gerencia elegir ese corte una vez tenga esa información de costos; Youden es el punto de partida razonable en ausencia de ella, no necesariamente el óptimo final.
test <- test %>%
mutate(clase_predicha_optima = factor(
if_else(prob_predicha > corte_optimo$threshold, "Si", "No"),
levels = c("No", "Si")
))
matriz_optima <- confusionMatrix(test$clase_predicha_optima, test$rotacion, positive = "Si")
matriz_optima
## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 267 20
## Si 102 51
##
## Accuracy : 0.7227
## 95% CI : (0.6784, 0.7641)
## No Information Rate : 0.8386
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.3014
##
## Mcnemar's Test P-Value : 2.244e-13
##
## Sensitivity : 0.7183
## Specificity : 0.7236
## Pos Pred Value : 0.3333
## Neg Pred Value : 0.9303
## Prevalence : 0.1614
## Detection Rate : 0.1159
## Detection Prevalence : 0.3477
## Balanced Accuracy : 0.7209
##
## 'Positive' Class : Si
##
Al bajar el punto de corte desde 0.5 hasta el valor de Youden, la sensibilidad aumenta de forma sustancial a costa de una especificidad algo menor, confirmando en un punto concreto el trade-off que ya se veía en la rejilla completa. La elección final de dónde operar sobre ese trade-off es, en el fondo, una decisión de negocio y no puramente estadística: depende de cómo la organización valore el costo de un falso negativo frente a un falso positivo, algo que este análisis no cuantifica pero que la tabla de sensibilidad del punto anterior deja lista para que la gerencia decida con esa información.
Construyo el perfil de un empleado hipotético para ilustrar cómo se usaría el modelo en la práctica: un empleado que trabaja horas extra, viaja frecuentemente, está soltero, reporta baja satisfacción laboral (nivel 1 de 4), tiene un ingreso mensual notoriamente por debajo del promedio de la organización, y lleva apenas un año en su cargo actual.
# En lugar de asumir una cifra de ingreso y afirmar que está "por
# debajo de la media", calculo la media real de ingreso_mensual y
# construyo el perfil hipotético como una fracción explícita de esa
# media, de modo que la afirmación quede demostrada por el propio
# código y no solo enunciada en el texto.
media_ingreso <- mean(train$ingreso_mensual)
ingreso_hipotetico <- round(media_ingreso * 0.5)
empleado_hipotetico <- tibble(
horas_extra = factor("Si", levels = levels(train$horas_extra)),
viaje_de_negocios = factor("Frecuentemente", levels = levels(train$viaje_de_negocios)),
estado_civil = factor("Soltero", levels = levels(train$estado_civil)),
satisfacion_laboral = 1,
ingreso_miles = ingreso_hipotetico / 1000,
antiguedad_cargo = 1
)
tibble(
`Ingreso del empleado hipotético` = ingreso_hipotetico,
`Ingreso medio de la organización` = round(media_ingreso, 1),
`¿Está por debajo de la media?` = ingreso_hipotetico < media_ingreso
) %>%
kable(align = "c", caption = "Verificación del ingreso hipotético frente a la media real") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Ingreso del empleado hipotético | Ingreso medio de la organización | ¿Está por debajo de la media? |
|---|---|---|
| 3215 | 6430.6 | TRUE |
probabilidad_hipotetica <- predict(modelo_logit, newdata = empleado_hipotetico, type = "response")
tibble(
`Probabilidad estimada de rotar` = round(probabilidad_hipotetica, 3),
`Punto de corte de referencia` = round(corte_optimo$threshold, 3),
`Decisión sugerida` = if_else(probabilidad_hipotetica > corte_optimo$threshold,
"Intervenir preventivamente", "No requiere intervención prioritaria")
) %>%
kable(align = "c", caption = "Predicción para el empleado hipotético") %>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
| Probabilidad estimada de rotar | Punto de corte de referencia | Decisión sugerida |
|---|---|---|
| 0.843 | 0.163 | Intervenir preventivamente |
El perfil hipotético combina, de manera deliberada, todos los factores de riesgo identificados como significativos en el modelo (horas extra, viaje frecuente, soltería, baja satisfacción, bajo ingreso relativo, poca antigüedad en el cargo), por lo que es razonable esperar una probabilidad de rotación estimada considerablemente por encima del punto de corte de referencia, lo que llevaría a recomendar una intervención preventiva sobre este empleado.
Factores que aumentan el riesgo de rotación (odds ratio mayor a 1): trabajar horas extra es, con claridad, el factor de mayor magnitud. Le siguen la frecuencia de viajes de negocio y el estado civil de soltería.
Factores que reducen el riesgo de rotación (odds ratio menor a 1): mayor satisfacción laboral, mayor ingreso mensual y mayor antigüedad en el cargo actual.
Con base en estos hallazgos, una posible estrategia de retención debería priorizar:
Estado_Civil y
Antigüedad_Cargo, así que no puedo afirmar que la
combinación de ambos (“soltero y con poca antigüedad”) sea más
riesgosa que la suma de sus efectos individuales — eso sería una
hipótesis razonable para una iteración futura del modelo, no una
conclusión que los datos actuales respalden. Con esa salvedad, planes de
desarrollo de carrera más visibles y mentoría temprana son acciones
razonables para ambos factores.Es importante ser cuidadoso con el alcance de estas conclusiones: el modelo identifica asociaciones estadísticas, no relaciones causales. Por ejemplo, no puedo afirmar que reducir las horas extra de un empleado específico reducirá su probabilidad de rotar en la magnitud exacta implicada por el odds ratio; solo puedo afirmar que, en los datos observados, los empleados con más horas extra rotan más, controlando por las demás variables del modelo.
brglm2::detect_separation() podría incorporarse como
diagnóstico adicional, pero no se realizó en este informe.Ingreso_Mensual ni Antigüedad_Cargo (ambas con
distribuciones sesgadas identificadas en el análisis univariado);
pruebas como Box-Tidwell o una comparación de AIC contra una
especificación con splines (splines::ns()) permitirían
confirmar si una forma no lineal ajusta mejor, en particular para
Antigüedad_Cargo, donde el efecto podría no ser monótono en
todo su rango.set.seed(123) garantiza que la partición
entrenamiento/prueba sea idéntica en cualquier ejecución posterior de
este documento.