Factores asociados a la rotación de empleados: un modelo de regresión logística

Actividad 3 — Modelos Estadísticos para la Toma de Decisiones (20262-A)

# La instalación del entorno se documenta en el Anexo A.
library(paqueteMODELOS)   # base de datos `rotacion`
library(dplyr)            # manipulación de datos
library(tidyr)            # reestructuración
library(ggplot2)          # gráficos
library(kableExtra)       # formato de tablas
library(car)              # Anova tipo II (LR) y GVIF
library(pROC)             # ROC, AUC e IC de DeLong
library(rms)              # bootstrap de optimismo (validate)
library(e1071)            # asimetría y curtosis
library(nortest)          # Anderson-Darling

# ---- Paleta base del informe
COL_JAV <- "#003087"; COL_AZUL <- "#0073B1"; COL_GRIS <- "#D4D4D4"
COL_NAR <- "#E8833A"; COL_TXT  <- "#4D4D4D"

# ---- Paleta ampliada.
# Se construye sobre el eje azul-naranja de la identidad visual, que es una
# pareja complementaria (matices ~205 y ~27 grados). Se añaden los matices
# análogos de cada extremo: hacia el azul, el verde azulado y el violeta;
# hacia el naranja, la terracota y el ámbar. Todos se fijan a una saturación
# y una luminosidad comparables, para que ningún color domine por brillo.
COL_TEAL <- "#1B8A8F"   # análogo frío del azul (verde azulado)
COL_VIO  <- "#6A5D9E"   # análogo frío del azul (violeta)
COL_TERR <- "#B4552F"   # análogo cálido del naranja (terracota)
COL_AMB  <- "#D9A22B"   # análogo cálido del naranja (ámbar)
PALETA   <- c(COL_AZUL, COL_NAR, COL_TEAL, COL_VIO, COL_TERR, COL_AMB, COL_JAV)

# Familias temáticas: el color codifica el contenido de la variable, no la
# decora. La misma familia conserva su matiz en todas las figuras.
FAMILIA <- c(
  Edad = "Perfil personal", Distancia_Casa = "Perfil personal",
  Genero = "Perfil personal", Estado_Civil = "Perfil personal",
  Ingreso_Mensual = "Compensación", Porcentaje_aumento_salarial = "Compensación",
  Antigüedad = "Tiempo y trayectoria", Antigüedad_Cargo = "Tiempo y trayectoria",
  Años_acargo_con_mismo_jefe = "Tiempo y trayectoria",
  Años_ultima_promoción = "Tiempo y trayectoria",
  Años_Experiencia = "Tiempo y trayectoria",
  Trabajos_Anteriores = "Formación y desarrollo",
  Capacitaciones = "Formación y desarrollo",
  Educación = "Formación y desarrollo", Campo_Educación = "Formación y desarrollo",
  Departamento = "Estructura de la empresa", Cargo = "Estructura de la empresa",
  Horas_Extra = "Condiciones de trabajo",
  `Viaje de Negocios` = "Condiciones de trabajo",
  Viaje_Negocios = "Condiciones de trabajo",
  Satisfacción_Ambiental = "Percepción del empleado",
  Satisfación_Laboral = "Percepción del empleado",
  Equilibrio_Trabajo_Vida = "Percepción del empleado",
  Rendimiento_Laboral = "Percepción del empleado",
  Rotación = "Respuesta")
COL_FAMILIA <- c(`Perfil personal` = COL_VIO, `Compensación` = COL_NAR,
                 `Tiempo y trayectoria` = COL_AZUL, `Formación y desarrollo` = COL_TEAL,
                 `Estructura de la empresa` = COL_JAV, `Condiciones de trabajo` = COL_TERR,
                 `Percepción del empleado` = COL_AMB, Respuesta = COL_NAR)
# Degradado de un mismo matiz, para los niveles ordenados de una ordinal
degradado <- function(col, k) colorRampPalette(c(colorspace_claro(col), col))(k)
colorspace_claro <- function(col, f = 0.72) {
  rgb1 <- grDevices::col2rgb(col)[, 1]
  grDevices::rgb(t(rgb1 + (255 - rgb1) * f), maxColorValue = 255)
}

theme_set(
  theme_minimal(base_size = 11) +
    theme(plot.title       = element_text(face = "bold", colour = COL_JAV),
          plot.subtitle    = element_text(colour = "grey45"),
          axis.title       = element_text(colour = "grey40"),
          axis.text        = element_text(colour = "grey40"),
          panel.grid.minor = element_blank(),
          panel.grid.major.x = element_blank(),
          legend.position  = "none")
)

# Etiquetas cortas de niveles para los ejes de los gráficos
ETQ_NIVEL <- function(x) dplyr::recode(x, No_Viaja = "No viaja", Frecuentemente = "Frecuente",
                                       Raramente = "Rara vez", Si = "Sí")

# ---- Números: coma decimal y punto de millar
num <- function(x, d = 2) {
  ifelse(is.na(x), "—",
         formatC(round(x, d), format = "f", digits = d,
                 big.mark = ".", decimal.mark = ","))
}
ent <- function(x) formatC(x, format = "d", big.mark = ".", decimal.mark = ",")
pct <- function(x, d = 1) paste0(num(100 * x, d), "\u00a0%")

# ---- Valores p: nunca "0"; bajo la precisión de máquina, < 2e-16
pval <- function(p, d = 3) {
  sapply(p, function(pi) {
    if (is.na(pi)) return("—")
    if (pi < 2.2e-16) return("$<2\\times10^{-16}$")
    if (pi < 10^(-d)) {
      e <- floor(log10(pi)); m <- pi / 10^e
      return(sprintf("$%s\\times10^{%d}$", num(m, 2), e))
    }
    num(pi, d)
  })
}

# ---- Celdas HTML: se escapan &, < y >, y el modo matemático $...$ pasa a \(...\),
# que MathJax procesa también dentro de las tablas
html_math <- function(x) {
  x <- as.character(x)
  # dentro de $...$ los signos < y > pasan a \\lt y \\gt, que MathJax interpreta sin
  # depender del escape HTML
  x <- vapply(x, function(si) {
    if (is.na(si)) return(si)
    m <- gregexpr("\\$[^$]*\\$", si)
    regmatches(si, m) <- lapply(regmatches(si, m), function(z)
      gsub(">", "\\\\gt ", gsub("<", "\\\\lt ", z)))
    si
  }, character(1), USE.NAMES = FALSE)
  x <- gsub("&", "&amp;", x, fixed = TRUE)
  x <- gsub("<", "&lt;", x, fixed = TRUE)
  x <- gsub(">", "&gt;", x, fixed = TRUE)
  gsub("\\$([^$]*)\\$", "\\\\(\\1\\\\)", x)
}

# Versión para texto en línea, dentro de $...$: incluye el signo (= o <)
pvm <- function(p, d = 3) {
  x <- pval(p, d)
  ifelse(grepl("^\\$<", x), gsub("\\$", "", x), paste0("=", gsub("\\$", "", x)))
}

# Decisión de un contraste con alfa = 0,05
decide <- function(p, a = 0.05) ifelse(p < a, "Se rechaza $H_0$", "No se rechaza $H_0$")

# ---- Envoltura única para todas las tablas (mismo formato de las Actividades 1 y 2)
tabla <- function(df, caption, digits = 2, align = NULL, ...) {
  df <- as.data.frame(df, check.names = FALSE)
  digits <- rep_len(digits, ncol(df))
  for (j in seq_along(df)) {
    v <- df[[j]]
    if (is.numeric(v)) {
      df[[j]] <- if (all(is.na(v) | v == round(v))) ifelse(is.na(v), "—", ent(v))
                 else num(v, digits[j])
    } else df[[j]] <- ifelse(is.na(v), "—", as.character(v))
  }
  if (is.null(align)) align <- ifelse(sapply(df, function(v)
    all(grepl("^[-$0-9.,<>\\\\{}times% —]+$", v))), "r", "l")
  for (j in seq_along(df)) df[[j]] <- html_math(df[[j]])
  kbl(df, format = "html", escape = FALSE, caption = html_math(caption),
      align = align, col.names = html_math(names(df)), row.names = FALSE) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = FALSE, position = "center", font_size = 13) %>%
    row_spec(0, bold = TRUE, background = COL_JAV, color = "white")
}

# ---- Etiquetas de coeficientes: se identifican por nombre, nunca por posición
ETIQUETAS <- c(
  "(Intercept)"                  = "Intercepto",
  "Horas_ExtraSi"                = "Horas extra: Sí vs. No",
  "Estado_CivilDivorciado"       = "Estado civil: Divorciado vs. Casado",
  "Estado_CivilSoltero"          = "Estado civil: Soltero vs. Casado",
  "Viaje_NegociosRaramente"      = "Viaje: Raramente vs. No viaja",
  "Viaje_NegociosFrecuentemente" = "Viaje: Frecuentemente vs. No viaja",
  "Edad"                         = "Edad (por año)",
  "log(Ingreso_Mensual)"         = "ln(Ingreso mensual)",
  "Antigüedad_Cargo"             = "Antigüedad en el cargo (por año)"
)
etiqueta <- function(n) {
  faltan <- setdiff(n, names(ETIQUETAS))
  if (length(faltan)) stop("Coeficiente sin etiqueta: ", paste(faltan, collapse = ", "))
  unname(ETIQUETAS[n])
}

data("rotacion")
crudo <- rotacion

# Tipología estadística de las 24 variables
NOMINALES  <- c("Rotación", "Departamento", "Campo_Educación", "Genero", "Cargo",
                "Estado_Civil", "Horas_Extra")
ORDINALES  <- c("Viaje de Negocios", "Educación", "Satisfacción_Ambiental",
                "Satisfación_Laboral", "Rendimiento_Laboral", "Equilibrio_Trabajo_Vida")
CONTINUAS  <- c("Edad", "Ingreso_Mensual")
DISCRETAS  <- c("Distancia_Casa", "Trabajos_Anteriores", "Porcentaje_aumento_salarial",
                "Años_Experiencia", "Capacitaciones", "Antigüedad", "Antigüedad_Cargo",
                "Años_ultima_promoción", "Años_acargo_con_mismo_jefe")
CUANTITATIVAS <- c(CONTINUAS, DISCRETAS)
stopifnot(setequal(c(NOMINALES, ORDINALES, CUANTITATIVAS), names(crudo)))

# Base de trabajo: factores con referencias explícitas y niveles ordinales en su orden
datos <- crudo %>%
  rename(Viaje_Negocios = `Viaje de Negocios`) %>%
  mutate(
    y              = as.integer(Rotación == "Si"),
    Rotación       = factor(Rotación, levels = c("No", "Si")),
    Horas_Extra    = relevel(factor(Horas_Extra), ref = "No"),
    Estado_Civil   = relevel(factor(Estado_Civil), ref = "Casado"),
    Viaje_Negocios = factor(Viaje_Negocios,
                            levels = c("No_Viaja", "Raramente", "Frecuentemente")),
    across(c(Departamento, Campo_Educación, Genero, Cargo), factor)
  )
stopifnot(!anyNA(datos$Viaje_Negocios), !anyNA(datos$y))
N  <- nrow(datos)
N1 <- sum(datos$y); N0 <- N - N1

# Etiquetas de las escalas ordinales según el diccionario de datos del curso
SATISF  <- c("1" = "Muy insatisfecho", "2" = "Insatisfecho", "3" = "Satisfecho",
             "4" = "Muy satisfecho")
ETQ_ORD <- list(
  Educación = c("1" = "Primaria", "2" = "Secundaria", "3" = "Técnico/tecnólogo",
                "4" = "Pregrado", "5" = "Posgrado"),
  Satisfacción_Ambiental  = SATISF,
  Satisfación_Laboral     = SATISF,
  Rendimiento_Laboral     = c("1" = "Bajo", "2" = "Medio", "3" = "Alto", "4" = "Muy alto"),
  Equilibrio_Trabajo_Vida = c("1" = "Muy bajo", "2" = "Bajo", "3" = "Medio", "4" = "Alto"))
etq_ord <- function(v, x) {
  e <- ETQ_ORD[[v]]
  if (is.null(e)) as.character(x) else paste0(x, " · ", e[as.character(x)])
}

# Núcleo de resultados que alimenta el resumen ejecutivo. Cada sección
# posterior reutiliza estos objetos y los desarrolla en detalle.
FORMULA <- y ~ Horas_Extra + Estado_Civil + Viaje_Negocios +
               Edad + log(Ingreso_Mensual) + Antigüedad_Cargo
modelo  <- glm(FORMULA, family = binomial(link = "logit"), data = datos)
k_param <- length(coef(modelo)) - 1

# Partición estratificada 70/30
set.seed(2026)
idx_tr <- sort(unlist(lapply(split(seq_len(N), datos$y),
                             function(i) sample(i, round(0.7 * length(i))))))
entren <- datos[idx_tr, ]; prueba <- datos[-idx_tr, ]
modelo_tr <- glm(FORMULA, family = binomial, data = entren)
p_tr <- predict(modelo_tr, entren, type = "response")
p_te <- predict(modelo_tr, prueba, type = "response")
roc_tr <- roc(entren$y, p_tr, levels = c(0, 1), direction = "<", quiet = TRUE)
roc_te <- roc(prueba$y, p_te, levels = c(0, 1), direction = "<", quiet = TRUE)
ic_te  <- ci.auc(roc_te, method = "delong")

# Métricas de clasificación para un corte c
metricas <- function(y, p, c) {
  pred <- as.integer(p >= c)
  VP <- sum(pred == 1 & y == 1); FP <- sum(pred == 1 & y == 0)
  VN <- sum(pred == 0 & y == 0); FN <- sum(pred == 0 & y == 1)
  sens <- VP / (VP + FN); esp <- VN / (VN + FP)
  vpp  <- if (VP + FP > 0) VP / (VP + FP) else NA
  vpn  <- if (VN + FN > 0) VN / (VN + FN) else NA
  data.frame(corte = c, Sens = sens, Esp = esp, VPP = vpp, VPN = vpn,
             Exactitud = (VP + VN) / length(y),
             F1 = if (is.na(vpp) || vpp + sens == 0) NA else 2 * vpp * sens / (vpp + sens),
             Youden = sens + esp - 1)
}
REJILLA <- seq(0.05, 0.95, by = 0.01)
rej_tr <- bind_rows(lapply(REJILLA, function(c) metricas(entren$y, p_tr, c)))
rej_te <- bind_rows(lapply(REJILLA, function(c) metricas(prueba$y, p_te, c)))
c_youden <- rej_tr$corte[which.max(rej_tr$Youden)]
m_youden_te <- metricas(prueba$y, p_te, c_youden)

# Matriz de confusión y métricas completas para un corte c (positivo: Y = 1)
cm_metricas <- function(y, p, c) {
  pred <- as.integer(p >= c); n <- length(y)
  VP <- sum(pred == 1 & y == 1); FP <- sum(pred == 1 & y == 0)
  VN <- sum(pred == 0 & y == 0); FN <- sum(pred == 0 & y == 1)
  wil  <- function(k, m) prop.test(k, m, correct = FALSE)$conf.int[1:2]   # IC de Wilson
  acc  <- (VP + VN) / n
  nir  <- max(mean(y), 1 - mean(y))           # tasa de no información
  sens <- VP / (VP + FN); esp <- VN / (VN + FP)
  vpp  <- VP / (VP + FP); vpn <- VN / (VN + FN)
  pe   <- ((VP + FP) * (VP + FN) + (VN + FN) * (VN + FP)) / n^2   # acuerdo por azar
  list(
    corte = c, cm = c(VP = VP, FP = FP, VN = VN, FN = FN),
    acc = acc, acc_ic = binom.test(VP + VN, n)$conf.int[1:2],     # IC de Clopper-Pearson
    nir = nir, p_nir = binom.test(VP + VN, n, p = nir, alternative = "greater")$p.value,
    kappa = (acc - pe) / (1 - pe),
    sens = sens, sens_ic = wil(VP, VP + FN), esp = esp, esp_ic = wil(VN, VN + FP),
    vpp = vpp, vpp_ic = wil(VP, VP + FP), vpn = vpn, vpn_ic = wil(VN, VN + FN),
    f1 = 2 * vpp * sens / (vpp + sens), bal = (sens + esp) / 2, J = sens + esp - 1,
    mcc = (VP * VN - FP * FN) / sqrt(as.numeric(VP + FP) * (VP + FN) * (VN + FP) * (VN + FN)),
    fpr = 1 - esp, fnr = 1 - sens, lrp = sens / (1 - esp), lrn = (1 - sens) / esp,
    prev = mean(y), detec = VP / n, senal = (VP + FP) / n,
    p_mcnemar = mcnemar.test(matrix(c(VP, FN, FP, VN), 2))$p.value)
}
met_05  <- cm_metricas(prueba$y, p_te, 0.5)
met_yt  <- cm_metricas(prueba$y, p_te, c_youden)
met_ytr <- cm_metricas(entren$y, p_tr, c_youden)

# Empleado hipotético
PERFIL <- data.frame(Edad = 32, Ingreso_Mensual = 2500, Antigüedad_Cargo = 1,
                     Horas_Extra = factor("No", levels = levels(datos$Horas_Extra)),
                     Estado_Civil = factor("Soltero", levels = levels(datos$Estado_Civil)),
                     Viaje_Negocios = factor("Raramente", levels = levels(datos$Viaje_Negocios)))
pr_eta <- predict(modelo, PERFIL, type = "link", se.fit = TRUE)
pi_hat <- plogis(pr_eta$fit)
pi_ic  <- plogis(pr_eta$fit + c(-1, 1) * qnorm(0.975) * pr_eta$se.fit)

OR_he <- exp(coef(modelo)["Horas_ExtraSi"])

Resumen ejecutivo

Con la base rotacion (1.470 empleados, 24 variables) se estimó un modelo de regresión logística para la probabilidad de rotación. En la base rotaron 237 empleados (16,1 %). Se seleccionaron seis covariables con hipótesis previas: horas extra, estado civil, frecuencia de viajes de negocios, edad, ingreso mensual y antigüedad en el cargo. Las seis resultaron asociadas con la rotación en el análisis bivariado y conservaron su significancia en el modelo ajustado (pruebas de razón de verosimilitudes, \(\alpha=0{,}05\)). Todas las direcciones coinciden con las hipótesis.

El factor de mayor peso es trabajar horas extra: manteniendo constantes las demás covariables, multiplica los odds de rotación por 4,34. Le siguen el estado civil (los solteros rotan más que los casados), los viajes frecuentes y el ingreso; la edad y la antigüedad en el cargo tienen efectos protectores de menor magnitud. En una muestra de prueba independiente, el modelo alcanza un AUC de 0,740 (IC 95 % de DeLong: 0,671–0,810). El corte que maximiza el índice de Youden en entrenamiento es 0,25; en prueba da una sensibilidad (recall) de 50,7 %, una especificidad de 85,1 %, una precisión de 39,6 % y un \(F_1\) de 0,444. La exactitud (79,6 %) no supera la de clasificar a todos como “no rota” (83,9 %), algo esperable con una prevalencia del 16 %; por eso el desempeño se juzga con métricas que no dependen de la clase mayoritaria (exactitud balanceada 0,679, \(\kappa=0,322\), MCC \(=0,325\)). Para un empleado hipotético de 32 años, soltero, sin horas extra, que viaja raramente, con ingreso de 2.500 y un año en el cargo, la probabilidad estimada de rotar es 0,254. La estrategia de retención prioriza la redistribución de las horas extra, la política de viajes y la revisión salarial.

1 Contexto y objetivos

Una organización busca comprender y prever los factores asociados a la rotación de sus empleados, con el fin de tomar medidas proactivas de retención. La gerencia propone un modelo de regresión logística que estime la probabilidad de que un empleado rote y que identifique los factores de mayor incidencia.

Objetivo general. Estimar la probabilidad de rotación de un empleado en función de características laborales y sociodemográficas, e identificar los factores asociados de mayor peso.

Objetivos específicos.

1. Seleccionar tres covariables categóricas y tres cuantitativas y formular, para cada una, una hipótesis con mecanismo y dirección esperada.
2. Caracterizar la base de datos con indicadores y gráficos adecuados al tipo de cada variable.
3. Contrastar la asociación de cada covariable con la rotación mediante la prueba apropiada a su tipo.
4. Estimar e interpretar un modelo de regresión logística con las seis covariables.
5. Evaluar el poder predictivo del modelo con la curva ROC y el AUC en datos no usados para estimarlo.
6. Predecir la probabilidad de rotación de un empleado hipotético y analizar la sensibilidad de la decisión de intervenir al punto de corte.
7. Proponer una estrategia de retención basada en los factores significativos.

Nota

Inconsistencias del enunciado y cómo se resuelven. (i) El punto 3 remite a “la hipótesis planteada en el punto 2”, pero las hipótesis se formulan en el punto 1; se contrasta contra las del punto 1. (ii) El punto 7 pide la estrategia con base en las variables significativas “en el punto 3”; se reportan las significativas del bivariado y del modelo ajustado, y la estrategia prioriza las que lo son en ambos. (iii) El problema habla de “cambio de cargo”, pero Rotación no distingue entre cambio interno y retiro; la interpretación se ciñe a lo que la variable registra. (iv) Algunos nombres de la base tienen erratas (Satisfación_Laboral); se usan tal cual.

2 Datos

2.1 Origen y codificación de la respuesta

Los datos provienen de la base rotacion del paquete paqueteMODELOS, adaptada de Weiers (2006). La respuesta se codifica como \(Y=1\) si el empleado rotó (Rotación = "Si") y \(Y=0\) si no rotó (Rotación = "No"). La Tabla 2.1 verifica la recodificación.

tabla(as.data.frame.matrix(table(`Rotación (original)` = crudo$Rotación,
                                 `Y (recodificada)` = datos$y)) %>%
        tibble::rownames_to_column("Rotación (original)") %>%
        rename(`Y = 0` = `0`, `Y = 1` = `1`),
      "Verificación de la recodificación de la variable respuesta.")

Tabla 2.1: Verificación de la recodificación de la variable respuesta.

Rotación (original)	Y = 0	Y = 1
No	1.233	0
Si	0	237

Decisión metodológica

La variable Viaje de Negocios se renombra Viaje_Negocios en la base de trabajo, porque los espacios en el nombre obligan a usar comillas invertidas en las fórmulas. No se modifica ningún valor.

2.2 Diccionario y tipología de las variables

DESCRIPCION <- c(
  "Rotación" = "El empleado rotó (Si/No)",
  "Edad" = "Edad en años",
  "Viaje de Negocios" = "Frecuencia de viajes de negocios",
  "Departamento" = "Departamento de la empresa",
  "Distancia_Casa" = "Distancia desde la casa (km)",
  "Educación" = "Nivel educativo: 1 = primaria, 2 = secundaria, 3 = técnico/tecnólogo, 4 = pregrado, 5 = posgrado",
  "Campo_Educación" = "Área de formación",
  "Satisfacción_Ambiental" = "Satisfacción ambiental: 1 = muy insatisfecho, 2 = insatisfecho, 3 = satisfecho, 4 = muy satisfecho",
  "Genero" = "Género (F/M)",
  "Cargo" = "Cargo actual",
  "Satisfación_Laboral" = "Satisfacción laboral: 1 = muy insatisfecho, 2 = insatisfecho, 3 = satisfecho, 4 = muy satisfecho",
  "Estado_Civil" = "Estado civil",
  "Ingreso_Mensual" = "Ingreso mensual",
  "Trabajos_Anteriores" = "Cantidad de trabajos antes de ingresar a la empresa",
  "Horas_Extra" = "Trabaja horas extra (Si/No)",
  "Porcentaje_aumento_salarial" = "Último aumento salarial (%)",
  "Rendimiento_Laboral" = "Rendimiento laboral: 1 = bajo, 2 = medio, 3 = alto, 4 = muy alto",
  "Años_Experiencia" = "Años de experiencia laboral total",
  "Capacitaciones" = "Capacitaciones en el último año",
  "Equilibrio_Trabajo_Vida" = "Equilibrio trabajo-vida: 1 = muy bajo, 2 = bajo, 3 = medio, 4 = alto",
  "Antigüedad" = "Años en la empresa",
  "Antigüedad_Cargo" = "Años en el cargo actual",
  "Años_ultima_promoción" = "Años desde la última promoción",
  "Años_acargo_con_mismo_jefe" = "Años con el mismo jefe")
stopifnot(setequal(names(DESCRIPCION), names(crudo)))

tipo_est <- function(v) {
  if (v %in% NOMINALES) "Cualitativa nominal"
  else if (v %in% ORDINALES) "Cualitativa ordinal"
  else if (v %in% CONTINUAS) "Cuantitativa continua"
  else "Cuantitativa discreta"
}
dominio <- function(v) {
  x <- crudo[[v]]
  if (v == "Viaje de Negocios") "No_Viaja, Raramente, Frecuentemente"
  else if (is.character(x)) paste(sort(unique(x)), collapse = ", ")
  else if (v %in% ORDINALES) paste(sort(unique(x)), collapse = ", ")
  else paste0("[", num(min(x), 0), "; ", num(max(x), 0), "]")
}
tipologia <- data.frame(
  Variable = names(crudo),
  Descripción = unname(DESCRIPCION[names(crudo)]),
  `Tipo en R` = sapply(crudo, function(x) class(x)[1]),
  `Tipo estadístico` = sapply(names(crudo), tipo_est),
  `Valores o rango` = sapply(names(crudo), dominio),
  check.names = FALSE, row.names = NULL)
tabla(tipologia, "Diccionario y tipología estadística de las 24 variables.", font = 7.5) %>%
  column_spec(2, width = "4.2cm") %>% column_spec(5, width = "4.6cm")

Tabla 2.2: Diccionario y tipología estadística de las 24 variables.

Variable	Descripción	Tipo en R	Tipo estadístico	Valores o rango
Rotación	El empleado rotó (Si/No)	character	Cualitativa nominal	No, Si
Edad	Edad en años	numeric	Cuantitativa continua	[18; 60]
Viaje de Negocios	Frecuencia de viajes de negocios	character	Cualitativa ordinal	No_Viaja, Raramente, Frecuentemente
Departamento	Departamento de la empresa	character	Cualitativa nominal	IyD, RH, Ventas
Distancia_Casa	Distancia desde la casa (km)	numeric	Cuantitativa discreta	[1; 29]
Educación	Nivel educativo: 1 = primaria, 2 = secundaria, 3 = técnico/tecnólogo, 4 = pregrado, 5 = posgrado	numeric	Cualitativa ordinal	1, 2, 3, 4, 5
Campo_Educación	Área de formación	character	Cualitativa nominal	Ciencias, Humanidades, Mercadeo, Otra, Salud, Tecnicos
Satisfacción_Ambiental	Satisfacción ambiental: 1 = muy insatisfecho, 2 = insatisfecho, 3 = satisfecho, 4 = muy satisfecho	numeric	Cualitativa ordinal	1, 2, 3, 4
Genero	Género (F/M)	character	Cualitativa nominal	F, M
Cargo	Cargo actual	character	Cualitativa nominal	Director_Investigación, Director_Manofactura, Ejecutivo_Ventas, Gerente, Investigador_Cientifico, Recursos_Humanos, Representante_Salud, Representante_Ventas, Tecnico_Laboratorio
Satisfación_Laboral	Satisfacción laboral: 1 = muy insatisfecho, 2 = insatisfecho, 3 = satisfecho, 4 = muy satisfecho	numeric	Cualitativa ordinal	1, 2, 3, 4
Estado_Civil	Estado civil	character	Cualitativa nominal	Casado, Divorciado, Soltero
Ingreso_Mensual	Ingreso mensual	numeric	Cuantitativa continua	[1.009; 19.999]
Trabajos_Anteriores	Cantidad de trabajos antes de ingresar a la empresa	numeric	Cuantitativa discreta	[0; 9]
Horas_Extra	Trabaja horas extra (Si/No)	character	Cualitativa nominal	No, Si
Porcentaje_aumento_salarial	Último aumento salarial (%)	numeric	Cuantitativa discreta	[11; 25]
Rendimiento_Laboral	Rendimiento laboral: 1 = bajo, 2 = medio, 3 = alto, 4 = muy alto	numeric	Cualitativa ordinal	3, 4
Años_Experiencia	Años de experiencia laboral total	numeric	Cuantitativa discreta	[0; 40]
Capacitaciones	Capacitaciones en el último año	numeric	Cuantitativa discreta	[0; 6]
Equilibrio_Trabajo_Vida	Equilibrio trabajo-vida: 1 = muy bajo, 2 = bajo, 3 = medio, 4 = alto	numeric	Cualitativa ordinal	1, 2, 3, 4
Antigüedad	Años en la empresa	numeric	Cuantitativa discreta	[0; 40]
Antigüedad_Cargo	Años en el cargo actual	numeric	Cuantitativa discreta	[0; 18]
Años_ultima_promoción	Años desde la última promoción	numeric	Cuantitativa discreta	[0; 15]
Años_acargo_con_mismo_jefe	Años con el mismo jefe	numeric	Cuantitativa discreta	[0; 17]

La base tiene 7 variables nominales, 6 ordinales, 2 continuas y 9 discretas (Tabla 2.2). Las escalas codificadas con números (Educación, las satisfacciones, Rendimiento_Laboral, Equilibrio_Trabajo_Vida) se tratan como ordinales: sus valores indican orden, pero la distancia entre categorías no tiene significado métrico. Viaje de Negocios es ordinal por la frecuencia creciente de sus niveles.

escalas <- bind_rows(lapply(names(ETQ_ORD), function(v)
  data.frame(Variable = c(v, rep("", length(ETQ_ORD[[v]]) - 1)),
             Código = names(ETQ_ORD[[v]]), Categoría = unname(ETQ_ORD[[v]]),
             `n observado` = as.integer(table(factor(crudo[[v]], levels = names(ETQ_ORD[[v]])))),
             check.names = FALSE)))
tabla(escalas, "Codificación de las escalas ordinales según el diccionario de datos del curso.",
      align = c("l", "c", "l", "r"))

Tabla 2.3: Codificación de las escalas ordinales según el diccionario de datos del curso.

Variable	Código	Categoría	n observado
Educación	1	Primaria	170
	2	Secundaria	282
	3	Técnico/tecnólogo	572
	4	Pregrado	398
	5	Posgrado	48
Satisfacción_Ambiental	1	Muy insatisfecho	284
	2	Insatisfecho	287
	3	Satisfecho	453
	4	Muy satisfecho	446
Satisfación_Laboral	1	Muy insatisfecho	289
	2	Insatisfecho	280
	3	Satisfecho	442
	4	Muy satisfecho	459
Rendimiento_Laboral	1	Bajo	0
	2	Medio	0
	3	Alto	1.244
	4	Muy alto	226
Equilibrio_Trabajo_Vida	1	Muy bajo	80
	2	Bajo	344
	3	Medio	893
	4	Alto	153

Nota sobre el diccionario de datos

Las etiquetas de las escalas ordinales, la unidad de Distancia_Casa (kilómetros) y la definición de Trabajos_Anteriores (trabajos previos al ingreso a la empresa) provienen del diccionario de datos que la profesora del curso publicó para esta actividad (Tabla 2.3). En ese diccionario las dos escalas de satisfacción asignan “Muy insatisfecho” tanto al código 1 como al 4; por el orden creciente de la escala (1 = muy insatisfecho, 2 = insatisfecho, 3 = satisfecho), el código 4 se interpreta como “Muy satisfecho”. La escala de Equilibrio_Trabajo_Vida termina en “Alto”, como la define el diccionario. Rendimiento_Laboral solo registra los códigos 3 (alto) y 4 (muy alto): ningún empleado tiene una calificación baja o media.

3 Análisis exploratorio y calidad de los datos

Esta etapa se ejecuta sobre las 24 variables antes de seleccionar covariables. Ningún análisis posterior usa la base sin haber pasado por ella.

3.1 Calidad de los datos

n_dup <- sum(duplicated(crudo))
n_na  <- sum(is.na(crudo))
n_completos <- sum(complete.cases(crudo))
ESCALAS_1_4 <- c("Satisfacción_Ambiental", "Satisfación_Laboral", "Equilibrio_Trabajo_Vida")
CHK_DOM <- list(
  "Cuantitativas sin valores negativos" = apply(crudo[CUANTITATIVAS] >= 0, 1, all),
  "Edad en el rango laboral [18; 65]" = crudo$Edad >= 18 & crudo$Edad <= 65,
  "Educación entera entre 1 y 5" = crudo$Educación %in% 1:5,
  "Satisfacciones y equilibrio enteros entre 1 y 4" =
    apply(crudo[ESCALAS_1_4], 1, function(r) all(r %in% 1:4)),
  "Rendimiento laboral entero entre 1 y 4" = crudo$Rendimiento_Laboral %in% 1:4,
  "Porcentaje de aumento en [0; 100]" =
    crudo$Porcentaje_aumento_salarial >= 0 & crudo$Porcentaje_aumento_salarial <= 100,
  "Etiquetas cualitativas sin espacios sobrantes" =
    apply(crudo[sapply(crudo, is.character)], 1, function(r) all(r == trimws(r))))
viol_dom <- sapply(CHK_DOM, function(r) !r)
REGLAS <- list(
  "Antigüedad en el cargo $\\le$ antigüedad en la empresa" =
    crudo$Antigüedad_Cargo <= crudo$Antigüedad,
  "Años con el mismo jefe $\\le$ antigüedad en la empresa" =
    crudo$Años_acargo_con_mismo_jefe <= crudo$Antigüedad,
  "Años desde la última promoción $\\le$ antigüedad en la empresa" =
    crudo$Años_ultima_promoción <= crudo$Antigüedad,
  "Antigüedad en la empresa $\\le$ experiencia total" =
    crudo$Antigüedad <= crudo$Años_Experiencia,
  "Edad $-$ experiencia $\\ge 14$ años" =
    crudo$Edad - crudo$Años_Experiencia >= 14)
viol <- sapply(REGLAS, function(r) !r)
n_viol_unicos <- sum(rowSums(viol) > 0)
edad_inicio_min <- min(crudo$Edad - crudo$Años_Experiencia)
# Varianza casi nula: razón de frecuencias > 95/5 y menos del 10 % de valores únicos
nzv <- bind_rows(lapply(names(crudo), function(v) {
  f <- sort(table(crudo[[v]]), decreasing = TRUE)
  data.frame(razon = if (length(f) > 1) f[1] / f[2] else NA,
             unicos = 100 * length(f) / N,
             fmin = if (v %in% c(NOMINALES, ORDINALES)) as.numeric(min(f)) else NA)
}))
n_nzv <- sum(nzv$razon > 95 / 5 & nzv$unicos < 10, na.rm = TRUE)
min_nivel <- min(nzv$fmin, na.rm = TRUE)
# Cada registro excluido se cuenta una sola vez, aunque falle varias reglas
excluir <- duplicated(crudo) | rowSums(viol_dom) > 0 | rowSums(viol) > 0
n_excl  <- sum(excluir)
stopifnot(nrow(datos) == nrow(crudo) - n_excl)
calidad <- data.frame(
  Bloque = c("Integridad", "", "Dominio", rep("", length(CHK_DOM) - 1),
             "Consistencia", rep("", length(REGLAS) - 1), "Variabilidad", "", "Depuración"),
  Verificación = c("Filas duplicadas exactas", "Celdas con dato faltante",
                   names(CHK_DOM), names(REGLAS),
                   "Variables con varianza casi nula",
                   "Casos en el nivel cualitativo menos frecuente",
                   "Registros excluidos / registros finales"),
  Resultado = c(ent(n_dup), ent(n_na), ent(colSums(viol_dom)), ent(colSums(viol)),
                ent(n_nzv), ent(min_nivel), paste0(ent(n_excl), " / ", ent(N))),
  check.names = FALSE)
tabla(calidad, "Verificaciones de calidad de la base (violaciones o conteos).",
      align = c("l", "l", "r"), font = 8.5)

Tabla 3.1: Verificaciones de calidad de la base (violaciones o conteos).

Bloque	Verificación	Resultado
Integridad	Filas duplicadas exactas	0
	Celdas con dato faltante	0
Dominio	Cuantitativas sin valores negativos	0
	Edad en el rango laboral [18; 65]	0
	Educación entera entre 1 y 5	0
	Satisfacciones y equilibrio enteros entre 1 y 4	0
	Rendimiento laboral entero entre 1 y 4	0
	Porcentaje de aumento en [0; 100]	0
	Etiquetas cualitativas sin espacios sobrantes	0
Consistencia	Antigüedad en el cargo \(\le\) antigüedad en la empresa	0
	Años con el mismo jefe \(\le\) antigüedad en la empresa	0
	Años desde la última promoción \(\le\) antigüedad en la empresa	0
	Antigüedad en la empresa \(\le\) experiencia total	0
	Edad \(-\) experiencia \(\ge 14\) años	0
Variabilidad	Variables con varianza casi nula	0
	Casos en el nivel cualitativo menos frecuente	27
Depuración	Registros excluidos / registros finales	0 / 1.470

La Tabla 3.1 resume las verificaciones de calidad. La base no tiene filas duplicadas (0) ni datos faltantes: los 1.470 registros son casos completos, así que no aplica la prueba MCAR de Little ni la imputación. Todas las variables están dentro de su dominio y ningún registro viola las restricciones lógicas; la diferencia mínima entre edad y experiencia es de 18 años. Con el criterio usual (razón entre la frecuencia del valor más común y la del segundo mayor que \(95/5\), y menos del 10 % de valores únicos) hay 0 variables con varianza casi nula, y el nivel cualitativo menos frecuente tiene 27 casos, así que no hay niveles escasos que agrupar. Rendimiento_Laboral solo toma los valores 3 y 4 de su escala (alto y muy alto). La base depurada conserva los 1.470 registros, y todos los análisis siguientes la usan.

3.2 Valores atípicos

Se aplican dos criterios a las 11 variables cuantitativas: el de Tukey (1977), que marca los valores fuera de \([Q_1-1{,}5\,\mathrm{RIC};\;Q_3+1{,}5\,\mathrm{RIC}]\), y el \(z\) robusto, que marca los valores con \(|x-\tilde x|/(1{,}4826\,\mathrm{MAD})>3{,}5\), donde \(\tilde x\) es la mediana.

atip <- bind_rows(lapply(CUANTITATIVAS, function(v) {
  x <- crudo[[v]]; q <- quantile(x, c(.25, .75)); r <- diff(q)
  tk <- sum(x < q[1] - 1.5 * r | x > q[2] + 1.5 * r)
  s  <- 1.4826 * mad(x, constant = 1)
  zr <- if (s > 0) sum(abs(x - median(x)) / s > 3.5) else NA
  data.frame(Variable = v, `Tukey (n)` = tk, `Tukey (%)` = 100 * tk / length(x),
             `z robusto (n)` = zr, `z robusto (%)` = 100 * zr / length(x),
             check.names = FALSE)
}))
tabla(atip, "Valores atípicos en las variables cuantitativas según dos criterios.",
      digits = c(0, 0, 1, 0, 1))

Tabla 3.2: Valores atípicos en las variables cuantitativas según dos criterios.

Variable	Tukey (n)	Tukey (%)	z robusto (n)	z robusto (%)
Edad	0	0,0	0	0,0
Ingreso_Mensual	114	7,8	118	8,0
Distancia_Casa	0	0,0	0	0,0
Trabajos_Anteriores	52	3,5	101	6,9
Porcentaje_aumento_salarial	0	0,0	18	1,2
Años_Experiencia	63	4,3	46	3,1
Capacitaciones	238	16,2	0	0,0
Antigüedad	104	7,1	66	4,5
Antigüedad_Cargo	21	1,4	0	0,0
Años_ultima_promoción	107	7,3	183	12,4
Años_acargo_con_mismo_jefe	14	1,0	0	0,0

largo_q <- crudo %>% select(all_of(CUANTITATIVAS)) %>%
  pivot_longer(everything(), names_to = "var", values_to = "x") %>%
  mutate(var = factor(var, levels = CUANTITATIVAS))
largo_q <- largo_q %>% mutate(familia = factor(FAMILIA[as.character(var)],
                                               levels = names(COL_FAMILIA)))
ggplot(largo_q, aes(x = 0, y = x, fill = familia)) +
  geom_boxplot(width = 0.5, colour = COL_TXT, alpha = 0.85,
               outlier.colour = COL_TXT, outlier.size = 0.8, outlier.alpha = 0.5) +
  scale_fill_manual(values = COL_FAMILIA, drop = TRUE, name = NULL) +
  facet_wrap(~ var, scales = "free_y", ncol = 4) +
  labs(title = "Los atípicos están en la cola superior, salvo en capacitaciones",
       subtitle = "El color identifica la familia temática de la variable", x = NULL, y = NULL) +
  theme(axis.text.x = element_blank(), strip.text = element_text(size = 7.5),
        legend.position = "bottom", legend.text = element_text(size = 7),
        legend.key.size = unit(0.35, "cm"))

Figura 3.1: Diagramas de caja de las 11 variables cuantitativas. Los puntos marcan los valores fuera de las cercas de Tukey.

colas <- sapply(CUANTITATIVAS, function(v) {
  x <- crudo[[v]]; q <- quantile(x, c(.25, .75)); r <- diff(q)
  c(bajo = sum(x < q[1] - 1.5 * r), alto = sum(x > q[2] + 1.5 * r))
})

La Figura 3.1 muestra dónde están los atípicos de la Tabla 3.2. Salvo en Capacitaciones, todos están en la cola superior, como corresponde a variables acotadas en cero y asimétricas a la derecha. En Capacitaciones la mitad central de los datos está en 2 y 3 capacitaciones, así que las cercas de Tukey dejan fuera tanto el 0 (54 casos) como el 5 y el 6 (184 casos): son valores posibles, no errores. Edad, Distancia_Casa y Porcentaje_aumento_salarial no tienen atípicos según Tukey.

Decisión metodológica

No se elimina ninguna observación por un criterio univariado. Un valor extremo pero posible (un ingreso alto, muchos años en la empresa) es información, no un error; solo se excluirían los valores imposibles, y los pasos anteriores no detectaron ninguno. En la regresión logística un atípico en \(x\) importa solo si es influyente, lo que se evalúa en la Sección 7.5 con la distancia de Cook, el apalancamiento y los DFBETAS. Para Ingreso_Mensual, que es asimétrico, se evalúa la escala logarítmica en la Sección 6.4.

3.3 Distribución de las variables cuantitativas

ancho <- function(x) {
  if (all(x == round(x)) && diff(range(x)) <= 40) 1
  else 2 * IQR(x) / length(x)^(1/3)          # regla de Freedman-Diaconis
}
capas <- lapply(CUANTITATIVAS, function(v)
  geom_histogram(data = filter(largo_q, var == v), aes(x),
                 binwidth = ancho(crudo[[v]]), boundary = if (ancho(crudo[[v]]) == 1) 0.5 else NULL,
                 fill = COL_FAMILIA[FAMILIA[v]], colour = "white", linewidth = 0.2,
                 alpha = 0.9))
ggplot() + capas +
  geom_vline(data = largo_q %>% group_by(var) %>% summarise(m = median(x), .groups = "drop"),
             aes(xintercept = m), colour = COL_TXT, linewidth = 0.5, linetype = "22") +
  facet_wrap(~ var, scales = "free", ncol = 4) +
  scale_x_continuous(n.breaks = 4, labels = function(x) num(x, 0)) +
  labs(title = "Casi todas las cuantitativas son asimétricas a la derecha",
       subtitle = "La línea punteada marca la mediana; el color, la familia temática",
       x = NULL, y = "Frecuencia") +
  theme(strip.text = element_text(size = 7.5))

Figura 3.2: Distribución de las 11 variables cuantitativas. Las variables discretas con pocos valores se grafican con un ancho de clase de 1; las demás, con el ancho de Freedman-Diaconis.

La Figura 3.2 anticipa lo que cuantifica la Tabla 5.2 (Sección 5.2): la edad es casi simétrica y el resto tiene cola derecha. Tres rasgos de forma son relevantes para el análisis. Ingreso_Mensual tiene una masa principal por debajo de 10.000 y un grupo reducido de ingresos altos. Antigüedad_Cargo y Años_acargo_con_mismo_jefe son multimodales, con picos en 0, 2 y 7 años; Años_ultima_promoción se concentra en 0 y 1 (el 63,8 % de los empleados). Porcentaje_aumento_salarial empieza en 11 %: no hay aumentos menores.

3.4 Distribución de las variables cualitativas

CUALI_EDA <- c(NOMINALES, ORDINALES)

bar_eda <- bind_rows(lapply(CUALI_EDA, function(v) {
  x <- as.character(crudo[[v]])
  niv <- if (v == "Viaje de Negocios") c("No_Viaja", "Raramente", "Frecuentemente")
         else if (v %in% ORDINALES) as.character(sort(unique(crudo[[v]])))
         else names(sort(table(x)))                       # de menor a mayor: arriba el más frecuente
  t <- table(factor(x, levels = niv))
  col_v <- unname(COL_FAMILIA[FAMILIA[v]])
  relleno <- if (v %in% ORDINALES) degradado(col_v, length(niv))
             else if (v == "Rotación") ifelse(niv == "Si", COL_NAR, COL_GRIS)
             else rep(col_v, length(niv))
  data.frame(var = v, nivel = niv, n = as.integer(t), relleno = relleno,
             clave = paste(v, etq_ord(v, niv), sep = "___"))
})) %>%
  mutate(var = factor(var, levels = CUALI_EDA),
         clave = factor(clave, levels = unique(clave)))
ggplot(bar_eda, aes(n, clave, fill = relleno)) +
  geom_col(width = 0.7) +
  geom_text(aes(label = pct(n / N, 0)), hjust = -0.1, size = 2.3, colour = COL_TXT) +
  scale_fill_identity() +
  scale_y_discrete(labels = function(k) sub("^.*___", "", k)) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.35))) +
  facet_wrap(~ var, scales = "free", ncol = 3) +
  labs(title = "Varias cualitativas tienen una categoría dominante",
       subtitle = "El color identifica la familia temática; las ordinales van en degradado",
       x = "Empleados", y = NULL) +
  theme(strip.text = element_text(size = 8), axis.text.y = element_text(size = 7),
        axis.text.x = element_text(size = 6), panel.grid.major.y = element_blank())

Figura 3.3: Distribución de las 13 variables cualitativas. Las nominales se ordenan por frecuencia; las ordinales, en su orden natural con un degradado del matiz de su familia. En Rotación, el naranja marca a quienes rotaron.

La Figura 3.3 muestra la estructura de la plantilla: dos tercios de los empleados están en Investigación y Desarrollo, los cargos más comunes son los de ventas y los técnicos o investigadores, y la mayoría viaja raramente. Entre las ordinales, las dos escalas de satisfacción se reparten de forma parecida entre sus cuatro niveles; Equilibrio_Trabajo_Vida se concentra en el nivel 3 (medio) y Rendimiento_Laboral solo toma los valores 3 (alto) y 4 (muy alto). Las frecuencias exactas están en las Tablas 5.4 y 5.5.

3.5 Asociaciones entre variables

Se calculan dos matrices de asociación: la correlación de Spearman entre las variables cuantitativas y ordinales, y la V de Cramér entre las nominales. Spearman usa solo el orden de los valores, así que es válida para las ordinales codificadas con números y no la afectan los atípicos ni la asimetría. Viaje de Negocios entra con los puntajes 1, 2, 3 de su orden.

NUM_EDA <- c(CUANTITATIVAS, ORDINALES)
M_num <- crudo %>% mutate(`Viaje de Negocios` = match(`Viaje de Negocios`,
                                   c("No_Viaja", "Raramente", "Frecuentemente"))) %>%
  select(all_of(NUM_EDA))
R_all <- cor(M_num, method = "spearman")
# Orden por conglomerado jerárquico de 1 - |r| para agrupar las variables afines
orden <- hclust(as.dist(1 - abs(R_all)), method = "average")$order
niv_h <- NUM_EDA[orden]
heat <- as.data.frame(as.table(R_all)) %>%
  setNames(c("v1", "v2", "r")) %>%
  mutate(v1 = factor(v1, levels = niv_h), v2 = factor(v2, levels = rev(niv_h)))
marco <- heat %>% filter(v1 %in% c("Edad", "Ingreso_Mensual", "Antigüedad_Cargo"),
                         v2 %in% c("Edad", "Ingreso_Mensual", "Antigüedad_Cargo"))
ggplot(heat, aes(v1, v2, fill = r)) +
  geom_tile(colour = "white", linewidth = 0.4) +
  geom_tile(data = marco, fill = NA, colour = COL_TXT, linewidth = 0.6) +
  geom_text(aes(label = num(r, 2), colour = abs(r) > 0.5), size = 2.1) +
  scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = COL_TXT), guide = "none") +
  scale_fill_gradient2(low = COL_NAR, mid = "white", high = COL_JAV, midpoint = 0,
                       limits = c(-1, 1)) +
  coord_equal() +
  labs(title = "Las variables de tiempo forman un bloque muy correlacionado",
       x = NULL, y = NULL) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 6.5),
        axis.text.y = element_text(size = 6.5), panel.grid = element_blank())

Figura 3.4: Correlación de Spearman entre las 11 variables cuantitativas y las 6 ordinales. Azul: correlación positiva; naranja: negativa. Los recuadros marcan las tres cuantitativas seleccionadas.

ut <- which(upper.tri(R_all), arr.ind = TRUE)
pares <- data.frame(a = NUM_EDA[ut[, 1]], b = NUM_EDA[ut[, 2]], r = R_all[ut]) %>%
  arrange(desc(abs(r)))
n_fuertes <- sum(abs(pares$r) > 0.7)
r_max_ord <- max(abs(R_all[ORDINALES, ORDINALES][upper.tri(R_all[ORDINALES, ORDINALES])]))
R_o <- abs(R_all[ORDINALES, setdiff(NUM_EDA, ORDINALES)])
R_o["Rendimiento_Laboral", "Porcentaje_aumento_salarial"] <- NA
r_ord_resto <- max(R_o, na.rm = TRUE)

La Figura 3.4 ordena las variables por afinidad y deja ver un bloque de variables de tiempo: Antigüedad, Antigüedad_Cargo, Años_acargo_con_mismo_jefe, Años_ultima_promoción y Años_Experiencia, junto con la edad y el ingreso. Hay 4 pares con \(|r_s|>0{,}7\); el mayor es Antigüedad–Antigüedad_Cargo (\(r_s=0,854\)). Este bloque es el que justifica elegir una sola variable de antigüedad para el modelo (Sección 4). Las ordinales son casi independientes entre sí (la mayor correlación entre ellas es \(|r_s|=0,030\)) y del resto, con una excepción: Rendimiento_Laboral y Porcentaje_aumento_salarial tienen \(r_s=0,629\), lo que es esperable si el aumento se asigna según la calificación de rendimiento. Sin ese par, ninguna ordinal supera \(|r_s|=0,205\) con otra variable.

cramer <- function(a, b) {
  t <- table(a, b)
  ch <- suppressWarnings(chisq.test(t, correct = FALSE))
  sqrt(unname(ch$statistic) / (sum(t) * (min(dim(t)) - 1)))
}
V <- outer(NOMINALES, NOMINALES, Vectorize(function(i, j)
  if (i == j) 1 else cramer(crudo[[i]], crudo[[j]])))
dimnames(V) <- list(NOMINALES, NOMINALES)
heat_v <- as.data.frame(as.table(V)) %>% setNames(c("v1", "v2", "V")) %>%
  mutate(v1 = factor(v1, levels = NOMINALES), v2 = factor(v2, levels = rev(NOMINALES)))
ggplot(heat_v, aes(v1, v2, fill = V)) +
  geom_tile(colour = "white", linewidth = 0.4) +
  geom_text(aes(label = num(V, 2), colour = V > 0.5), size = 2.5) +
  scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = COL_TXT), guide = "none") +
  scale_fill_gradient(low = "white", high = COL_TEAL, limits = c(0, 1)) +
  coord_equal() +
  labs(title = "Departamento y cargo son casi redundantes",
       x = NULL, y = NULL) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7),
        axis.text.y = element_text(size = 7), panel.grid = element_blank())

Figura 3.5: V de Cramér entre las 7 variables nominales, incluida la respuesta.

V_rot <- sort(V["Rotación", setdiff(NOMINALES, "Rotación")], decreasing = TRUE)

En la Figura 3.5, Departamento y Cargo tienen \(V=0,939\): cada cargo pertenece casi siempre a un solo departamento, así que incluir las dos en un modelo sería redundante. Campo_Educación tiene una asociación moderada con el departamento (\(V=0,590\)) y más débil con el cargo (\(V=0,343\)). Con la respuesta, la asociación nominal más fuerte es la de Horas_Extra (\(V=0,246\)), seguida de Cargo (\(V=0,242\)) y Estado_Civil (\(V=0,177\)). Cargo no se eligió como covariable a pesar de su asociación con la rotación: sus 9 niveles exigen 8 parámetros, lo mismo que las seis covariables juntas, y duplica la información de Departamento. Es una lectura descriptiva, sin contraste formal; las pruebas de asociación con la rotación se hacen en la Sección 6 para las covariables seleccionadas.

4 Selección de variables e hipótesis

Se seleccionan tres covariables categóricas y tres cuantitativas con un mecanismo plausible de relación con la rotación (Tabla 4.1). La referencia de cada factor es el nivel frente al cual se interpretan los odds ratio (OR).

HIP <- data.frame(
  Variable = c("Horas_Extra", "Estado_Civil", "Viaje_Negocios",
               "Edad", "Ingreso_Mensual", "Antigüedad_Cargo"),
  Tipo = c("Nominal (2 niveles)", "Nominal (3 niveles)", "Ordinal (3 niveles)",
           "Continua", "Continua", "Discreta"),
  Referencia = c("No", "Casado", "No_Viaja", "—", "—", "—"),
  Mecanismo = c(
    "La sobrecarga produce desgaste y conflicto entre trabajo y vida personal.",
    "Menos cargas familiares reducen el costo de cambiar de empleo.",
    "Los viajes aumentan el desgaste y el tiempo fuera de casa.",
    "Los trabajadores jóvenes tienen más movilidad y menos arraigo.",
    "Un salario mayor eleva el costo de oportunidad de irse.",
    "La permanencia refleja ajuste al puesto y capital específico."),
  `Dirección esperada` = c("Sí > No", "Soltero > Casado",
                           "Creciente con la frecuencia", "Negativa", "Negativa",
                           "Negativa"),
  check.names = FALSE)
HIP$Signo <- c(+1, +1, +1, -1, -1, -1)   # signo esperado del efecto principal
tabla(HIP[, 1:5], "Covariables seleccionadas, mecanismo y dirección esperada.",
      font = 8) %>%
  column_spec(4, width = "5.4cm") %>% column_spec(5, width = "2.6cm")

Tabla 4.1: Covariables seleccionadas, mecanismo y dirección esperada.

Variable	Tipo	Referencia	Mecanismo	Dirección esperada
Horas_Extra	Nominal (2 niveles)	No	La sobrecarga produce desgaste y conflicto entre trabajo y vida personal.	Sí > No
Estado_Civil	Nominal (3 niveles)	Casado	Menos cargas familiares reducen el costo de cambiar de empleo.	Soltero > Casado
Viaje_Negocios	Ordinal (3 niveles)	No_Viaja	Los viajes aumentan el desgaste y el tiempo fuera de casa.	Creciente con la frecuencia
Edad	Continua	—	Los trabajadores jóvenes tienen más movilidad y menos arraigo.	Negativa
Ingreso_Mensual	Continua	—	Un salario mayor eleva el costo de oportunidad de irse.	Negativa
Antigüedad_Cargo	Discreta	—	La permanencia refleja ajuste al puesto y capital específico.	Negativa

En Estado_Civil la hipótesis direccional se plantea para solteros frente a casados. Para los divorciados no se formula una dirección, porque dos mecanismos opuestos compiten: menos cargas que un casado, pero en muchos casos con responsabilidades económicas por separado.

Hipótesis y regla de decisión

Todas las pruebas son bilaterales con \(\alpha=0{,}05\).

• Factores (bivariado): \(H_0\): Rotación y \(X_j\) son independientes; \(H_1\): no lo son.
• Viaje_Negocios (bivariado): \(H_0\): la proporción de rotación no tiene tendencia lineal en los puntajes \(1,\,2,\,3\) de la frecuencia; \(H_1\): sí la tiene.
• Cuantitativas (bivariado): \(H_0\): \(P(X_1>X_0)+\tfrac12P(X_1=X_0)=\tfrac12\), donde \(X_1\) y \(X_0\) son observaciones al azar de quienes rotan y de quienes no; \(H_1\): \(\neq\tfrac12\).
• Modelo: \(H_0\): \(\beta_j=0\) (o todos los coeficientes del término son nulos, en los factores); \(H_1\): al menos uno es distinto de cero.

La dirección se lee en el signo de \(\hat\beta_j\) o en el OR. Veredicto: confirmada si \(p<\alpha\) y el signo observado coincide con el esperado; significativa en sentido contrario si \(p<\alpha\) y el signo es el opuesto; no confirmada si \(p\ge\alpha\).

Mapa de las pruebas de hipótesis. La Tabla 4.2 reúne todas las pruebas del informe: qué hipótesis nula contrasta cada una, con qué estadístico y distribución de referencia, por qué es la adecuada y qué supuestos exige. En cada sección, antes de la tabla de resultados, un recuadro detalla la prueba en el contexto de las variables, y cada tabla de resultados incluye la decisión con \(\alpha=0{,}05\).

PRUEBAS <- data.frame(
  Prueba = c("Shapiro-Wilk y Anderson-Darling", "$\\chi^2$ de independencia de Pearson",
             "Tendencia de Cochran-Armitage", "Residuos estandarizados ajustados",
             "Mann-Whitney (suma de rangos)", "Razón de verosimilitudes (logística simple)",
             "Box-Tidwell", "Wald", "Razón de verosimilitudes tipo II",
             "Razón de verosimilitudes global", "Hosmer-Lemeshow",
             "Razón de verosimilitudes: lineal vs. splines", "AUC = 0,5 (Mann-Whitney sobre $\\hat\\pi$)",
             "Binomial exacta: exactitud vs. NIR", "McNemar"),
  Sección = c("5.2", "6.1", "6.1", "6.1", "6.2", "6.3", "6.4", "7.1", "7.1", "7.1", "7.5",
              "7.6", "8.4", "8.4", "8.4"),
  `Hipótesis nula` = c(
    "La variable sigue una distribución normal",
    "Rotación y la covariable son independientes",
    "Sin tendencia lineal de la proporción de rotación en los puntajes 1, 2, 3",
    "La celda se ajusta a la independencia",
    "$A=P(X_1>X_0)+\\tfrac12P(X_1=X_0)=\\tfrac12$",
    "$\\beta_1=0$ en el modelo simple",
    "$\\gamma=0$: el logit es lineal en la covariable",
    "$\\beta_j=0$ (un coeficiente)",
    "Todos los coeficientes del término son nulos",
    "$\\beta_1=\\dots=\\beta_8=0$",
    "El modelo está bien calibrado",
    "Los términos no lineales son nulos",
    "El modelo no discrimina mejor que el azar",
    "Exactitud $\\le$ NIR",
    "$P(FP)=P(FN)$"),
  `Estadístico y distribución bajo $H_0$` = c(
    "$W$ y $A^2$ (distribuciones propias)",
    "$X^2=\\sum(O-E)^2/E\\sim\\chi^2_{(r-1)(c-1)}$",
    "$Z^2\\sim\\chi^2_1$",
    "$r_{ij}\\sim N(0,\\,1)$",
    "$U$; $z\\sim N(0,\\,1)$ con corrección por empates",
    "$G^2=D_0-D_1\\sim\\chi^2_{\\text{gl}}$",
    "$G^2\\sim\\chi^2_1$",
    "$z=\\hat\\beta_j/\\widehat{\\mathrm{EE}}\\sim N(0,\\,1)$",
    "$G^2\\sim\\chi^2_{\\text{gl}}$",
    "$G^2\\sim\\chi^2_8$",
    "$\\hat C\\sim\\chi^2_{g-2}$",
    "$G^2\\sim\\chi^2_6$",
    "$U$; $z\\sim N(0,\\,1)$",
    "Aciertos $\\sim\\text{Bin}(n,\\text{NIR})$",
    "$(|FP-FN|-1)^2/(FP+FN)\\sim\\chi^2_1$"),
  `Por qué esta prueba` = c(
    "Describe la forma; Anderson-Darling pondra más las colas",
    "Dos variables categóricas",
    "Covariable ordinal: usa el orden y concentra la potencia en 1 gl",
    "Localiza las celdas que generan la asociación",
    "Covariable cuantitativa no normal, con atípicos y empates",
    "Más fiable que Wald; da la dirección del efecto",
    "Contrasta el supuesto de linealidad en el logit",
    "Contrasta cada nivel por separado",
    "Factores de varios niveles; no depende del orden de entrada",
    "Significancia conjunta del modelo",
    "Bondad de ajuste con covariables continuas",
    "Compara modelos anidados",
    "Discriminación en datos nuevos",
    "Compara con la regla trivial de la clase mayoritaria",
    "Detecta si el corte sesga los errores hacia un lado"),
  Supuestos = c(
    "Observaciones independientes",
    "Independencia; condición de Cochran (si falla, Fisher)",
    "Independencia; puntajes fijados a priori",
    "Muestra grande",
    "Dos muestras independientes",
    "Observaciones independientes; MV regular",
    "Covariable positiva",
    "Normalidad asintótica del EMV",
    "Modelo sin interacciones",
    "Observaciones independientes",
    "Grupos con frecuencias esperadas suficientes",
    "Modelos anidados en los mismos datos",
    "Muestra de prueba independiente",
    "Muestra de prueba independiente",
    "Datos pareados (real y predicho del mismo empleado)"),
  check.names = FALSE)
PRUEBAS$`Por qué esta prueba`[1] <- "Describe la forma; Anderson-Darling pondera más las colas"
tabla(PRUEBAS, "Pruebas de hipótesis del informe: hipótesis nula, estadístico, justificación y supuestos.",
      align = c("l", "c", "l", "l", "l", "l")) %>%
  pack_rows("Análisis univariado", 1, 1) %>% pack_rows("Análisis bivariado", 2, 7) %>%
  pack_rows("Modelo de regresión logística", 8, 12) %>% pack_rows("Evaluación predictiva", 13, 15)

Tabla 4.2: Pruebas de hipótesis del informe: hipótesis nula, estadístico, justificación y supuestos.

Prueba	Sección	Hipótesis nula	Estadístico y distribución bajo \(H_0\)	Por qué esta prueba	Supuestos
Análisis univariado
Shapiro-Wilk y Anderson-Darling	5.2	La variable sigue una distribución normal	\(W\) y \(A^2\) (distribuciones propias)	Describe la forma; Anderson-Darling pondera más las colas	Observaciones independientes
Análisis bivariado
\(\chi^2\) de independencia de Pearson	6.1	Rotación y la covariable son independientes	\(X^2=\sum(O-E)^2/E\sim\chi^2_{(r-1)(c-1)}\)	Dos variables categóricas	Independencia; condición de Cochran (si falla, Fisher)
Tendencia de Cochran-Armitage	6.1	Sin tendencia lineal de la proporción de rotación en los puntajes 1, 2, 3	\(Z^2\sim\chi^2_1\)	Covariable ordinal: usa el orden y concentra la potencia en 1 gl	Independencia; puntajes fijados a priori
Residuos estandarizados ajustados	6.1	La celda se ajusta a la independencia	\(r_{ij}\sim N(0,\,1)\)	Localiza las celdas que generan la asociación	Muestra grande
Mann-Whitney (suma de rangos)	6.2	\(A=P(X_1\gt X_0)+\tfrac12P(X_1=X_0)=\tfrac12\)	\(U\); \(z\sim N(0,\,1)\) con corrección por empates	Covariable cuantitativa no normal, con atípicos y empates	Dos muestras independientes
Razón de verosimilitudes (logística simple)	6.3	\(\beta_1=0\) en el modelo simple	\(G^2=D_0-D_1\sim\chi^2_{\text{gl}}\)	Más fiable que Wald; da la dirección del efecto	Observaciones independientes; MV regular
Box-Tidwell	6.4	\(\gamma=0\): el logit es lineal en la covariable	\(G^2\sim\chi^2_1\)	Contrasta el supuesto de linealidad en el logit	Covariable positiva
Modelo de regresión logística
Wald	7.1	\(\beta_j=0\) (un coeficiente)	\(z=\hat\beta_j/\widehat{\mathrm{EE}}\sim N(0,\,1)\)	Contrasta cada nivel por separado	Normalidad asintótica del EMV
Razón de verosimilitudes tipo II	7.1	Todos los coeficientes del término son nulos	\(G^2\sim\chi^2_{\text{gl}}\)	Factores de varios niveles; no depende del orden de entrada	Modelo sin interacciones
Razón de verosimilitudes global	7.1	\(\beta_1=\dots=\beta_8=0\)	\(G^2\sim\chi^2_8\)	Significancia conjunta del modelo	Observaciones independientes
Hosmer-Lemeshow	7.5	El modelo está bien calibrado	\(\hat C\sim\chi^2_{g-2}\)	Bondad de ajuste con covariables continuas	Grupos con frecuencias esperadas suficientes
Razón de verosimilitudes: lineal vs. splines	7.6	Los términos no lineales son nulos	\(G^2\sim\chi^2_6\)	Compara modelos anidados	Modelos anidados en los mismos datos
Evaluación predictiva
AUC = 0,5 (Mann-Whitney sobre \(\hat\pi\))	8.4	El modelo no discrimina mejor que el azar	\(U\); \(z\sim N(0,\,1)\)	Discriminación en datos nuevos	Muestra de prueba independiente
Binomial exacta: exactitud vs. NIR	8.4	Exactitud \(\le\) NIR	Aciertos \(\sim\text{Bin}(n,\text{NIR})\)	Compara con la regla trivial de la clase mayoritaria	Muestra de prueba independiente
McNemar	8.4	\(P(FP)=P(FN)\)	\((|FP-FN|-1)^2/(FP+FN)\sim\chi^2_1\)	Detecta si el corte sesga los errores hacia un lado	Datos pareados (real y predicho del mismo empleado)

Decisión metodológica

Satisfacción laboral no entra entre las seis covariables. Es categórica ordinal y podría entrar como factor, igual que Viaje_Negocios. Se excluye por su papel causal plausible: la satisfacción puede ser mediadora entre las condiciones de trabajo (horas extra, viajes) y la rotación, y ajustar por ella absorbería parte del efecto de esas variables. Además, es una medida autorreportada, a diferencia de las seleccionadas, que son registros administrativos. Se describe en el análisis univariado.

Colinealidad entre candidatas. Años_Experiencia y Antigüedad se descartan como covariables por su correlación con las seleccionadas (Tabla 4.3). La experiencia tiene \(r_s=0,710\) con el ingreso y \(r_s=0,657\) con la edad; la antigüedad en la empresa tiene \(r_s=0,854\) con la antigüedad en el cargo. Entre las tres cuantitativas seleccionadas la correlación máxima es moderada, y la colinealidad en el modelo se verifica con el GVIF (Sección 7.5).

vars_sp <- c("Edad", "Ingreso_Mensual", "Antigüedad_Cargo", "Antigüedad", "Años_Experiencia")
R_sp <- cor(crudo[vars_sp], method = "spearman")
tabla(data.frame(Variable = vars_sp, R_sp, check.names = FALSE),
      "Correlaciones de Spearman entre las cuantitativas seleccionadas y las descartadas.",
      digits = 3, font = 8.5)

Tabla 4.3: Correlaciones de Spearman entre las cuantitativas seleccionadas y las descartadas.

Variable	Edad	Ingreso_Mensual	Antigüedad_Cargo	Antigüedad	Años_Experiencia
Edad	1,000	0,472	0,198	0,252	0,657
Ingreso_Mensual	0,472	1,000	0,395	0,464	0,710
Antigüedad_Cargo	0,198	0,395	1,000	0,854	0,493
Antigüedad	0,252	0,464	0,854	1,000	0,594
Años_Experiencia	0,657	0,710	0,493	0,594	1,000

max_sp_sel <- max(abs(R_sp[1:3, 1:3][upper.tri(R_sp[1:3, 1:3])]))

La correlación máxima entre las tres seleccionadas es \(|r_s|=0,472\).

5 Análisis univariado

5.1 Variable respuesta: Rotación

w <- prop.test(N1, N, correct = FALSE)$conf.int    # intervalo de Wilson
tabla(data.frame(Rotación = c("No (Y = 0)", "Sí (Y = 1)", "Total"),
                 n = c(N0, N1, N), `Proporción` = c(N0, N1, N) / N,
                 check.names = FALSE),
      "Distribución de la variable respuesta.", digits = c(0, 0, 4))

Tabla 5.1: Distribución de la variable respuesta.

Rotación	n	Proporción
No (Y = 0)	1.233	0,8388
Sí (Y = 1)	237	0,1612
Total	1.470	1,0000

data.frame(g = factor(c("No rotó", "Rotó"), levels = c("Rotó", "No rotó")),
           n = c(N0, N1)) %>%
  ggplot(aes(n, g, fill = g)) + geom_col(width = 0.6) +
  geom_text(aes(label = paste0(ent(n), " (", pct(n / N), ")")), hjust = -0.08,
            size = 3.3, colour = COL_TXT) +
  scale_fill_manual(values = c("Rotó" = COL_NAR, "No rotó" = COL_GRIS)) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.25))) +
  labs(title = "Uno de cada seis empleados rotó", x = "Empleados", y = NULL) +
  theme(panel.grid.major.y = element_blank())

Figura 5.1: Distribución de la rotación.

Rotaron 237 de los 1.470 empleados, una proporción de \(\hat p=0,1612\) con IC 95 % de Wilson \([0,1433;\ 0,1809]\) (Tabla 5.1, Figura 5.1). El intervalo de Wilson se prefiere al de Wald porque mantiene la cobertura nominal con proporciones alejadas de 0,5 (Agresti, 2013, §1.4).

Nota

Consecuencia del desbalance. Con una prevalencia cercana al 16 %, un modelo que asignara “No” a todos acertaría en el 83,9 % de los casos. Por eso la exactitud no sirve como medida principal, y un corte de 0,5 tenderá a clasificar a casi todos como “No”, con baja sensibilidad. El corte se decide en la Sección 9 con criterios explícitos.

5.2 Variables cuantitativas

desc_cuanti <- function(v) {
  x <- crudo[[v]]
  data.frame(Variable = v, n = length(x), Media = mean(x), DE = sd(x),
             Mediana = median(x), RIC = IQR(x), Mín = min(x), Máx = max(x),
             Asimetría = skewness(x, type = 2), `Curtosis (exceso)` = kurtosis(x, type = 2),
             check.names = FALSE)
}
tabla(bind_rows(lapply(CUANTITATIVAS, desc_cuanti)),
      "Estadísticos descriptivos de las variables cuantitativas.", font = 8,
      digits = c(0, 0, 2, 2, 1, 1, 0, 0, 3, 3))

Tabla 5.2: Estadísticos descriptivos de las variables cuantitativas.

Variable	n	Media	DE	Mediana	RIC	Mín	Máx	Asimetría	Curtosis (exceso)
Edad	1.470	36,92	9,14	36	13	18	60	0,413	-0,405
Ingreso_Mensual	1.470	6.502,93	4.707,96	4.919	5.468	1.009	19.999	1,370	1,005
Distancia_Casa	1.470	9,19	8,11	7	12	1	29	0,958	-0,225
Trabajos_Anteriores	1.470	2,69	2,50	2	3	0	9	1,026	0,010
Porcentaje_aumento_salarial	1.470	15,21	3,66	14	6	11	25	0,821	-0,301
Años_Experiencia	1.470	11,28	7,78	10	9	0	40	1,117	0,918
Capacitaciones	1.470	2,80	1,29	3	1	0	6	0,553	0,495
Antigüedad	1.470	7,01	6,13	5	6	0	40	1,765	3,936
Antigüedad_Cargo	1.470	4,23	3,62	3	5	0	18	0,917	0,477
Años_ultima_promoción	1.470	2,19	3,22	1	3	0	15	1,984	3,613
Años_acargo_con_mismo_jefe	1.470	4,12	3,57	3	5	0	17	0,833	0,171

La asimetría y la curtosis se calculan con los estimadores de tipo 2 (Joanes y Gill, 1998), que son los que reportan los programas estadísticos habituales; la curtosis se reporta en exceso, de modo que 0 corresponde a la normal. Todas las cuantitativas son asimétricas a la derecha salvo la edad, que es casi simétrica. El ingreso mensual es la más asimétrica de las seleccionadas: su media supera a su mediana en 1.584 unidades.

SEL_CUANTI <- c("Edad", "Ingreso_Mensual", "Antigüedad_Cargo")
largo <- crudo %>% select(all_of(SEL_CUANTI)) %>%
  pivot_longer(everything(), names_to = "var", values_to = "x") %>%
  mutate(var = factor(var, levels = SEL_CUANTI))
largo <- largo %>% mutate(col = unname(COL_FAMILIA[FAMILIA[as.character(var)]]))
g1 <- ggplot(largo, aes(x)) +
  geom_histogram(aes(y = after_stat(density), fill = col), bins = 30, colour = "white") +
  geom_density(colour = COL_TXT, linewidth = 0.6) +
  scale_fill_identity() +
  facet_wrap(~ var, scales = "free", nrow = 1) +
  labs(title = "Ingreso y antigüedad en el cargo, asimétricos a la derecha; edad, casi simétrica",
       x = NULL, y = "Densidad")
g2 <- ggplot(largo, aes(x, y = 0)) +
  geom_boxplot(aes(fill = col), width = 0.5, colour = COL_TXT, alpha = 0.85,
               outlier.colour = COL_TXT, outlier.size = 0.9, outlier.alpha = 0.5) +
  scale_fill_identity() +
  facet_wrap(~ var, scales = "free_x", nrow = 1) +
  labs(x = NULL, y = NULL) +
  theme(axis.text.y = element_blank(), strip.text = element_blank())
gridExtra::grid.arrange(g1, g2, heights = c(2.3, 1))

Figura 5.2: Distribución de las tres covariables cuantitativas seleccionadas. Arriba, histograma con densidad; abajo, diagrama de caja (los puntos son atípicos según Tukey).

ggplot(largo, aes(sample = x)) +
  stat_qq(aes(colour = col), size = 0.5, alpha = 0.6) + scale_colour_identity() +
  stat_qq_line(colour = COL_TXT, linewidth = 0.5) +
  facet_wrap(~ var, scales = "free", nrow = 1) +
  labs(title = "Las tres se apartan de la normalidad en las colas",
       x = "Cuantiles teóricos", y = "Cuantiles observados")

Figura 5.3: Gráficos cuantil-cuantil normales de las tres covariables cuantitativas seleccionadas.

Pruebas de normalidad de Shapiro-Wilk y Anderson-Darling

\(H_0\): la covariable sigue una distribución normal; \(H_1\): no la sigue. Shapiro-Wilk usa \(W=\left(\sum a_i x_{(i)}\right)^2/\sum(x_i-\bar x)^2\), que compara los estadísticos de orden con sus valores esperados bajo normalidad; Anderson-Darling usa \(A^2=n\int(F_n-F)^2\,[F(1-F)]^{-1}dF\), que pondera más las discrepancias en las colas. Se rechaza \(H_0\) si \(p<0{,}05\).

Por qué estas pruebas. Shapiro-Wilk (Shapiro y Wilk, 1965) es la prueba general de mayor potencia, y Anderson-Darling (Anderson y Darling, 1954) complementa la lectura de las colas, donde están los atípicos de estas variables. Se usan solo para describir la forma: la regresión logística no supone normalidad de las covariables, y la conclusión orienta la elección de pruebas no paramétricas en el análisis bivariado.

norm_tab <- bind_rows(lapply(SEL_CUANTI, function(v) {
  x <- crudo[[v]]; sw <- shapiro.test(x); ad <- ad.test(x)
  data.frame(Variable = v,
             `W (Shapiro-Wilk)` = sw$statistic, `p (SW)` = pval(sw$p.value),
             `A (Anderson-Darling)` = ad$statistic, `p (AD)` = pval(ad$p.value),
             Decisión = ifelse(sw$p.value < 0.05 & ad$p.value < 0.05,
                               "Se rechaza $H_0$ con ambas", "Ver valores p"),
             check.names = FALSE)
}))
tabla(norm_tab, "Pruebas de normalidad de las tres covariables cuantitativas seleccionadas.",
      digits = c(0, 4, 0, 3, 0, 0))

Tabla 5.3: Pruebas de normalidad de las tres covariables cuantitativas seleccionadas.

Variable	W (Shapiro-Wilk)	p (SW)	A (Anderson-Darling)	p (AD)	Decisión
Edad	0,9774	\(2,03\times10^{-14}\)	9,975	\(\lt 2\times10^{-16}\)	Se rechaza \(H_0\) con ambas
Ingreso_Mensual	0,8279	\(\lt 2\times10^{-16}\)	85,391	\(\lt 2\times10^{-16}\)	Se rechaza \(H_0\) con ambas
Antigüedad_Cargo	0,8962	\(\lt 2\times10^{-16}\)	50,263	\(\lt 2\times10^{-16}\)	Se rechaza \(H_0\) con ambas

Las pruebas de Shapiro-Wilk y Anderson-Darling rechazan la normalidad en las tres variables (Tabla 5.3; Figuras 5.2 y 5.3). Se reportan solo como descripción de la forma: con \(n=1.470\) estas pruebas detectan desviaciones de magnitud trivial, y la normalidad de las covariables no es un supuesto de la regresión logística. La forma sí orienta dos decisiones: el uso de pruebas no paramétricas en el bivariado y la evaluación de la escala logarítmica del ingreso.

5.3 Variables cualitativas

frec_de <- function(vars) bind_rows(lapply(vars, function(v) {
  x <- crudo[[v]]
  t <- if (v == "Viaje de Negocios")
         table(factor(x, levels = c("No_Viaja", "Raramente", "Frecuentemente")))
       else if (v %in% ORDINALES) table(factor(x, levels = sort(unique(x))))
       else sort(table(x), decreasing = TRUE)
  data.frame(Variable = c(v, rep("", length(t) - 1)), Nivel = etq_ord(v, names(t)),
             n = as.integer(t), `%` = 100 * as.numeric(t) / length(x), check.names = FALSE)
}))
tabla(frec_de(setdiff(NOMINALES, "Rotación")),
      "Frecuencias de las variables cualitativas nominales.",
      digits = c(0, 0, 0, 1), font = 8)

Tabla 5.4: Frecuencias de las variables cualitativas nominales.

Variable	Nivel	n	%
Departamento	IyD	961	65,4
	Ventas	446	30,3
	RH	63	4,3
Campo_Educación	Ciencias	606	41,2
	Salud	464	31,6
	Mercadeo	159	10,8
	Tecnicos	132	9,0
	Otra	82	5,6
	Humanidades	27	1,8
Genero	M	882	60,0
	F	588	40,0
Cargo	Ejecutivo_Ventas	326	22,2
	Investigador_Cientifico	292	19,9
	Tecnico_Laboratorio	259	17,6
	Director_Manofactura	145	9,9
	Representante_Salud	131	8,9
	Gerente	102	6,9
	Representante_Ventas	83	5,6
	Director_Investigación	80	5,4
	Recursos_Humanos	52	3,5
Estado_Civil	Casado	673	45,8
	Soltero	470	32,0
	Divorciado	327	22,2
Horas_Extra	No	1.054	71,7
	Si	416	28,3

tabla(frec_de(ORDINALES), "Frecuencias de las variables cualitativas ordinales.",
      digits = c(0, 0, 0, 1), font = 8)

Tabla 5.5: Frecuencias de las variables cualitativas ordinales.

Variable	Nivel	n	%
Viaje de Negocios	No_Viaja	150	10,2
	Raramente	1.043	71,0
	Frecuentemente	277	18,8
Educación	1 · Primaria	170	11,6
	2 · Secundaria	282	19,2
	3 · Técnico/tecnólogo	572	38,9
	4 · Pregrado	398	27,1
	5 · Posgrado	48	3,3
Satisfacción_Ambiental	1 · Muy insatisfecho	284	19,3
	2 · Insatisfecho	287	19,5
	3 · Satisfecho	453	30,8
	4 · Muy satisfecho	446	30,3
Satisfación_Laboral	1 · Muy insatisfecho	289	19,7
	2 · Insatisfecho	280	19,0
	3 · Satisfecho	442	30,1
	4 · Muy satisfecho	459	31,2
Rendimiento_Laboral	3 · Alto	1.244	84,6
	4 · Muy alto	226	15,4
Equilibrio_Trabajo_Vida	1 · Muy bajo	80	5,4
	2 · Bajo	344	23,4
	3 · Medio	893	60,7
	4 · Alto	153	10,4

ord_num <- setdiff(ORDINALES, "Viaje de Negocios")
moda <- function(x) { t <- table(x); names(t)[which.max(t)] }
tabla(data.frame(Variable = c(ord_num, "Viaje de Negocios"),
                 Mediana = c(sapply(ord_num, function(v) etq_ord(v, median(crudo[[v]]))),
                             levels(datos$Viaje_Negocios)[median(as.integer(datos$Viaje_Negocios))]),
                 Moda = c(sapply(ord_num, function(v) etq_ord(v, moda(crudo[[v]]))),
                          moda(crudo$`Viaje de Negocios`)),
                 row.names = NULL),
      "Mediana y moda de las variables ordinales.")

Tabla 5.6: Mediana y moda de las variables ordinales.

Variable	Mediana	Moda
Educación	3 · Técnico/tecnólogo	3 · Técnico/tecnólogo
Satisfacción_Ambiental	3 · Satisfecho	3 · Satisfecho
Satisfación_Laboral	3 · Satisfecho	4 · Muy satisfecho
Rendimiento_Laboral	3 · Alto	3 · Alto
Equilibrio_Trabajo_Vida	3 · Medio	3 · Medio
Viaje de Negocios	Raramente	Raramente

En las ordinales se reportan la mediana y la moda, que solo usan el orden de las categorías (Tabla 5.6). Satisfación_Laboral tiene mediana 3 · Satisfecho: al menos la mitad de los empleados está satisfecha o muy satisfecha, aunque el 38,7 % está insatisfecho o muy insatisfecho. La moda de la educación es el pregrado y la del equilibrio trabajo-vida, el nivel medio.

SEL_CUALI <- c("Horas_Extra", "Estado_Civil", "Viaje_Negocios")
bar <- bind_rows(lapply(SEL_CUALI, function(v) {
  t <- table(datos[[v]])
  data.frame(var = v, nivel = names(t), n = as.integer(t))
})) %>%
  mutate(var = factor(var, levels = SEL_CUALI),
         nivel = factor(nivel, levels = unique(c(levels(datos$Horas_Extra),
                                                  levels(datos$Estado_Civil),
                                                  levels(datos$Viaje_Negocios)))),
         relleno = case_when(nivel == "No_Viaja" ~ degradado(COL_TERR, 3)[1],
                             nivel == "Raramente" ~ degradado(COL_TERR, 3)[2],
                             nivel == "Frecuentemente" ~ COL_TERR,
                             nivel %in% c("No", "Si") ~ COL_TERR,
                             TRUE ~ COL_VIO))
ggplot(bar, aes(nivel, n, fill = relleno)) + geom_col(width = 0.65) +
  geom_text(aes(label = paste0(ent(n), "\n", pct(n / N))), vjust = -0.15, size = 2.7,
            colour = COL_TXT, lineheight = 0.85) +
  scale_fill_identity() +
  scale_x_discrete(labels = ETQ_NIVEL) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.3))) +
  facet_wrap(~ var, scales = "free_x", nrow = 1) +
  labs(title = "La mayoría no trabaja horas extra, está casada y viaja raramente",
       x = NULL, y = "Empleados")

Figura 5.4: Distribución de las tres covariables categóricas seleccionadas. Los niveles de viaje de negocios se muestran en su orden natural, con un degradado del matiz de su familia.

La mayoría de los empleados no trabaja horas extra (71,7 %), está casada (45,8 %) y viaja raramente (71,0 %). El nivel menos frecuente entre las seleccionadas es No_Viaja, con 150 casos: suficientes para estimar su efecto, aunque es la referencia con menos información.

6 Análisis bivariado

La respuesta es \(Y\) (\(1\) = rotó, \(0\) = no rotó). Para cada covariable se usa la prueba que corresponde a su tipo, se reporta un tamaño de efecto y se estima una regresión logística simple, cuyo signo responde a la pregunta del enunciado sobre la dirección del efecto.

6.1 Covariables categóricas

prop_nivel <- function(v) datos %>% group_by(nivel = .data[[v]]) %>%
  summarise(n = n(), Rotan = sum(y), .groups = "drop") %>%
  mutate(var = v, prop = Rotan / n)
pn <- bind_rows(lapply(SEL_CUALI, prop_nivel))
tabla(pn %>% transmute(Variable = ifelse(duplicated(var), "", var), Nivel = as.character(nivel),
                       n, `Rotan (n)` = Rotan, `Proporción de rotación` = prop),
      "Proporción de rotación por nivel de las covariables categóricas.",
      digits = c(0, 0, 0, 0, 4))

Tabla 6.1: Proporción de rotación por nivel de las covariables categóricas.

Variable	Nivel	n	Rotan (n)	Proporción de rotación
Horas_Extra	No	1.054	110	0,1044
	Si	416	127	0,3053
Estado_Civil	Casado	673	84	0,1248
	Divorciado	327	33	0,1009
	Soltero	470	120	0,2553
Viaje_Negocios	No_Viaja	150	12	0,0800
	Raramente	1.043	156	0,1496
	Frecuentemente	277	69	0,2491

ggplot(pn %>% mutate(var = factor(var, levels = SEL_CUALI)),
       aes(nivel, prop)) +
  geom_col(width = 0.6, fill = COL_NAR) +
  geom_hline(yintercept = N1 / N, linetype = "dashed", colour = "grey50") +
  geom_text(aes(label = pct(prop)), vjust = -0.3, size = 2.9, colour = COL_TXT) +
  scale_y_continuous(labels = function(x) pct(x, 0), expand = expansion(mult = c(0, 0.18))) +
  scale_x_discrete(labels = ETQ_NIVEL) +
  facet_wrap(~ var, scales = "free_x", nrow = 1) +
  labs(title = "Rotan más: con horas extra, solteros y viajeros frecuentes",
       x = NULL, y = "Proporción que rota")

Figura 6.1: Proporción de rotación por nivel de las covariables categóricas. La línea punteada es la proporción global.

Prueba \(\chi^2\) de independencia de Pearson

\(H_0\): Rotación y la covariable \(X\) son independientes, \(\pi_{ij}=\pi_{i+}\pi_{+j}\) en todas las celdas; \(H_1\): al menos una celda cumple \(\pi_{ij}\neq\pi_{i+}\pi_{+j}\). El estadístico es \(X^2=\sum_{i,j}(O_{ij}-E_{ij})^2/E_{ij}\), con \(E_{ij}=n_{i+}n_{+j}/n\); bajo \(H_0\), \(X^2\overset{a}{\sim}\chi^2_{(r-1)(c-1)}\). Se rechaza \(H_0\) si \(p<0{,}05\).

Por qué esta prueba. Las dos variables son categóricas y cada empleado aporta una sola observación, así que las observaciones son independientes. La aproximación \(\chi^2\) requiere frecuencias esperadas suficientes; se verifica la condición de Cochran (ninguna \(E_{ij}<1\) y al menos el 80 % de las \(E_{ij}\ge5\)) y, si fallara, se usaría la prueba exacta de Fisher. No se aplica la corrección de Yates, que es conservadora. Como con \(n=1.470\) el valor \(p\) no mide la magnitud de la asociación, se reporta la \(V\) de Cramér, \(V=\sqrt{X^2/[n(\min(r,c)-1)]}\in[0,\,1]\).

prueba_chi <- function(v) {
  t  <- table(datos[[v]], datos$y)
  ch <- suppressWarnings(chisq.test(t, correct = FALSE))
  cochran <- all(ch$expected >= 1) && mean(ch$expected >= 5) >= 0.8
  p <- if (cochran) ch$p.value else fisher.test(t)$p.value
  data.frame(Variable = v, Prueba = if (cochran) "Chi-cuadrado" else "Fisher",
             `Chi-cuadrado` = unname(ch$statistic), gl = unname(ch$parameter),
             `Valor p` = p, `Mín. esperado` = min(ch$expected),
             `V de Cramér` = sqrt(unname(ch$statistic) / (sum(t) * (min(dim(t)) - 1))),
             Decisión = decide(p), check.names = FALSE)
}
chi_tab <- bind_rows(lapply(SEL_CUALI, prueba_chi))
tabla(chi_tab %>% mutate(`Valor p` = pval(`Valor p`)),
      "Pruebas de independencia entre cada covariable categórica y la rotación.",
      digits = c(0, 0, 3, 0, 0, 1, 3))

Tabla 6.2: Pruebas de independencia entre cada covariable categórica y la rotación.

Variable	Prueba	Chi-cuadrado	gl	Valor p	Mín. esperado	V de Cramér	Decisión
Horas_Extra	Chi-cuadrado	89,044	1	\(\lt 2\times10^{-16}\)	67,1	0,246	Se rechaza \(H_0\)
Estado_Civil	Chi-cuadrado	46,164	2	\(9,46\times10^{-11}\)	52,7	0,177	Se rechaza \(H_0\)
Viaje_Negocios	Chi-cuadrado	24,182	2	\(5,61\times10^{-6}\)	24,2	0,128	Se rechaza \(H_0\)

# Cochran-Armitage para la tendencia en Viaje_Negocios (puntajes 1, 2, 3)
tv <- table(datos$Viaje_Negocios, datos$y)
ca <- prop.trend.test(x = tv[, "1"], n = rowSums(tv), score = 1:3)

Las tres covariables categóricas cumplen la condición de Cochran (todas las frecuencias esperadas son mayores que 5; la mínima es 24,2), así que se usa la prueba \(\chi^2\) de independencia sin corrección y no hace falta la de Fisher. Las tres rechazan la independencia (Tabla 6.2). Por la V de Cramér, la asociación más fuerte es la de las horas extra (\(V=0,246\)).

Prueba de tendencia de Cochran-Armitage (Armitage, 1955)

\(H_0\): la proporción de rotación no cambia linealmente con la frecuencia de viaje, \(\beta_1=0\) en \(\pi_k=\beta_0+\beta_1 s_k\) con puntajes \(s_k=1,\,2,\,3\); \(H_1\): \(\beta_1\neq0\). El estadístico es \[ Z^2=\frac{\left[\sum_k s_k\,(y_k-n_k\bar p)\right]^2}{\bar p(1-\bar p)\sum_k n_k(s_k-\bar s)^2}\ \overset{a}{\sim}\ \chi^2_1 , \] donde \(y_k\) es el número de rotaciones en el nivel \(k\), \(\bar p\) la proporción global y \(\bar s=\sum_k n_ks_k/n\).

Por qué esta prueba. Viaje_Negocios es ordinal. La \(\chi^2\) general ignora el orden y reparte la evidencia en 2 gl; esta prueba la concentra en la dirección monótona que plantea la hipótesis, y por eso es más potente cuando la asociación es monótona. Los puntajes se fijan a priori por el orden natural de los niveles.

Para Viaje_Negocios, que es ordinal, la prueba principal es la de tendencia de Cochran-Armitage con puntajes \(1,\,2,\,3\): \(\chi^2_1=23,712\), \(p=1,12\times10^{-6}\). La proporción de rotación crece con la frecuencia de los viajes, de 8,0 % entre quienes no viajan a 24,9 % entre quienes viajan con frecuencia. La prueba de tendencia usa un solo grado de libertad y es más potente que la \(\chi^2\) general cuando la asociación es monótona (Agresti, 2013, §3.4.6).

res_aj <- bind_rows(lapply(SEL_CUALI, function(v) {
  ch <- suppressWarnings(chisq.test(table(datos[[v]], datos$y), correct = FALSE))
  data.frame(Variable = v, Nivel = rownames(ch$stdres),
             `Residuo (Y = 0)` = ch$stdres[, "0"], `Residuo (Y = 1)` = ch$stdres[, "1"],
             check.names = FALSE)
})) %>% mutate(Variable = ifelse(duplicated(Variable), "", Variable))
tabla(res_aj, "Residuos estandarizados ajustados de las tablas de contingencia.",
      digits = c(0, 0, 2, 2))

Tabla 6.3: Residuos estandarizados ajustados de las tablas de contingencia.

Variable	Nivel	Residuo (Y = 0)	Residuo (Y = 1)
Horas_Extra	No	9,44	-9,44
	Si	-9,44	9,44
Estado_Civil	Casado	3,49	-3,49
	Divorciado	3,36	-3,36
	Soltero	-6,73	6,73
Viaje_Negocios	No_Viaja	2,85	-2,85
	Raramente	1,90	-1,90
	Frecuentemente	-4,41	4,41

Los residuos estandarizados ajustados siguen aproximadamente una \(N(0,\,1)\) bajo independencia, así que \(|r|>1{,}96\) señala las celdas que generan la asociación (Agresti, 2013, §2.4.5). En la Tabla 6.3, las celdas con exceso de rotación son las horas extra, los solteros y los viajes frecuentes; los divorciados y los casados tienen rotación por debajo de la esperada, y quienes no viajan también.

6.2 Covariables cuantitativas

Prueba de Mann-Whitney (Mann y Whitney, 1947)

\(H_0\): \(A=P(X_1>X_0)+\tfrac12P(X_1=X_0)=\tfrac12\), es decir, el valor de la covariable de un empleado que rota no tiende a ser mayor ni menor que el de uno que no rota; \(H_1\): \(A\neq\tfrac12\). El estadístico es \(U=\sum_{i\in G_1}R_i-n_1(n_1+1)/2\), con \(R_i\) los rangos en la muestra combinada; bajo \(H_0\), \(z=(U-n_1n_0/2)/\sigma_U\overset{a}{\sim}N(0,\,1)\), con \(\sigma_U\) corregida por empates. Se rechaza \(H_0\) si \(p<0{,}05\).

Por qué esta prueba. Las tres covariables no son normales (Tabla 5.3), tienen atípicos y, en Antigüedad_Cargo, muchos empates; la prueba \(t\) compararía medias, que son sensibles a las colas. Mann-Whitney solo usa el orden de los datos y exige dos muestras independientes, lo que se cumple porque cada empleado está en un solo grupo. Con \(n_1=237\) y \(n_0=1.233\) la aproximación normal es precisa.

mw_tab <- bind_rows(lapply(SEL_CUANTI, function(v) {
  x1 <- datos[[v]][datos$y == 1]; x0 <- datos[[v]][datos$y == 0]
  wt <- wilcox.test(x1, x0, alternative = "two.sided", exact = FALSE)
  A  <- unname(wt$statistic) / (length(x1) * length(x0))
  data.frame(Variable = v, `Mediana (Y = 1)` = median(x1), `Mediana (Y = 0)` = median(x0),
             U = unname(wt$statistic), `Valor p` = wt$p.value, `A estimado` = A,
             `r biserial` = 2 * A - 1, Decisión = decide(wt$p.value), check.names = FALSE)
}))
tabla(mw_tab %>% mutate(`Valor p` = pval(`Valor p`)),
      "Prueba de Mann-Whitney para las covariables cuantitativas según la rotación.",
      digits = c(0, 1, 1, 0, 0, 3, 3))

Tabla 6.4: Prueba de Mann-Whitney para las covariables cuantitativas según la rotación.

Variable	Mediana (Y = 1)	Mediana (Y = 0)	U	Valor p	A estimado	r biserial	Decisión
Edad	32	36	106.855	\(5,28\times10^{-11}\)	0,366	-0,269	Se rechaza \(H_0\)
Ingreso_Mensual	3.202	5.204	100.620	\(2,95\times10^{-14}\)	0,344	-0,311	Se rechaza \(H_0\)
Antigüedad_Cargo	2	3	105.214	\(4,43\times10^{-12}\)	0,360	-0,280	Se rechaza \(H_0\)

Se usa la prueba de Mann-Whitney bilateral, por la asimetría y los atípicos de las tres variables. El estadístico \(U\) se calcula para el grupo \(Y=1\), de modo que \(\hat A=U/(n_1n_0)\) estima \(P(X_1>X_0)+\tfrac12P(X_1=X_0)\); la prueba contrasta la dominancia estocástica (\(A=\tfrac12\)), no la igualdad de medianas, porque las formas de las dos distribuciones pueden diferir. El tamaño de efecto es la correlación biserial de rangos \(r_{rb}=2\hat A-1\in[-1;1]\): un valor negativo indica valores menores entre quienes rotan. Las tres pruebas rechazan \(H_0\) con efectos negativos y de magnitud similar (Tabla 6.4): quienes rotan son más jóvenes, ganan menos y llevan menos tiempo en el cargo. La mediana del ingreso de quienes rotan es 3.202 frente a 5.204 de quienes no rotan.

datos %>% select(y, all_of(SEL_CUANTI)) %>%
  pivot_longer(-y, names_to = "var", values_to = "x") %>%
  mutate(var = factor(var, levels = SEL_CUANTI),
         g = factor(ifelse(y == 1, "Rotó", "No rotó"), levels = c("No rotó", "Rotó"))) %>%
  ggplot(aes(g, x, fill = g)) +
  geom_boxplot(width = 0.55, colour = COL_TXT, outlier.size = 0.6, outlier.colour = "grey60") +
  scale_fill_manual(values = c("No rotó" = COL_GRIS, "Rotó" = COL_NAR)) +
  facet_wrap(~ var, scales = "free_y", nrow = 1) +
  labs(title = "Quienes rotan son más jóvenes, ganan menos y llevan menos tiempo en el cargo",
       x = NULL, y = NULL)

Figura 6.2: Distribución de las covariables cuantitativas según la rotación.

6.3 Regresión logística simple

Prueba de razón de verosimilitudes en cada regresión logística simple

\(H_0\): \(\beta_1=0\) (o \(\beta_1=\beta_2=0\) en los factores de tres niveles) en \(\operatorname{logit}\pi_i=\beta_0+\beta_1x_i\); \(H_1\): algún coeficiente es distinto de cero. El estadístico es la diferencia de devianzas entre el modelo nulo y el simple, \(G^2=D_0-D_1\overset{a}{\sim}\chi^2_{\text{gl}}\). Se prefiere a la de Wald porque es invariante a la parametrización y más fiable en muestras finitas (Hauck y Donner, 1977). Además de la significancia, el modelo simple da el signo del efecto, que el enunciado pide interpretar.

simple <- function(term, v) {
  f  <- reformulate(term, "y")
  m1 <- glm(f, binomial, datos); m0 <- glm(y ~ 1, binomial, datos)
  ic <- suppressMessages(confint(m1))
  lr <- anova(m0, m1, test = "LRT")
  cf <- coef(m1)[-1]
  data.frame(Variable = v, Coeficiente = etiqueta(names(cf)), beta = cf,
             OR = exp(cf), `IC 2,5 %` = exp(ic[-1, 1]), `IC 97,5 %` = exp(ic[-1, 2]),
             `p (LR del término)` = c(lr$`Pr(>Chi)`[2], rep(NA, length(cf) - 1)),
             check.names = FALSE, row.names = NULL)
}
TERMINOS <- c(Horas_Extra = "Horas_Extra", Estado_Civil = "Estado_Civil",
              Viaje_Negocios = "Viaje_Negocios", Edad = "Edad",
              Ingreso_Mensual = "log(Ingreso_Mensual)", Antigüedad_Cargo = "Antigüedad_Cargo")
sl <- bind_rows(Map(simple, TERMINOS, names(TERMINOS)))
tabla(sl %>% mutate(Variable = ifelse(duplicated(Variable), "", Variable),
                    `p (LR del término)` = ifelse(is.na(`p (LR del término)`), "",
                                                   pval(`p (LR del término)`))) %>%
        rename(`$\\hat\\beta$` = beta),
      "Regresiones logísticas simples: coeficiente, OR con IC 95 % por verosimilitud perfilada y prueba de razón de verosimilitudes.",
      digits = c(0, 0, 4, 3, 3, 3, 0), font = 8.5)

Tabla 6.5: Regresiones logísticas simples: coeficiente, OR con IC 95 % por verosimilitud perfilada y prueba de razón de verosimilitudes.

Variable	Coeficiente	\(\hat\beta\)	OR	IC 2,5 %	IC 97,5 %	p (LR del término)
Horas_Extra	Horas extra: Sí vs. No	1,3274	3,771	2,832	5,032	\(\lt 2\times10^{-16}\)
Estado_Civil	Estado civil: Divorciado vs. Casado	-0,2395	0,787	0,508	1,194	\(2,79\times10^{-10}\)
	Estado civil: Soltero vs. Casado	0,8772	2,404	1,769	3,281
Viaje_Negocios	Viaje: Raramente vs. No viaja	0,7044	2,023	1,139	3,931	\(6,93\times10^{-6}\)
	Viaje: Frecuentemente vs. No viaja	1,3389	3,815	2,061	7,641
Edad	Edad (por año)	-0,0523	0,949	0,933	0,965	\(3,25\times10^{-10}\)
Ingreso_Mensual	ln(Ingreso mensual)	-0,9053	0,404	0,317	0,512	\(4,44\times10^{-15}\)
Antigüedad_Cargo	Antigüedad en el cargo (por año)	-0,1463	0,864	0,823	0,905	\(6,38\times10^{-11}\)

b_ing_s <- sl$beta[sl$Coeficiente == "ln(Ingreso mensual)"]

Cada modelo simple es \(\operatorname{logit}\pi_i=\beta_0+\beta_1x_i\) (o con indicadoras para los factores), y la prueba de razón de verosimilitudes compara la devianza con la del modelo nulo. Los signos coinciden con los de las pruebas anteriores (Tabla 6.5):

• Horas extra: \(\hat\beta>0\). Los odds de rotación de quienes hacen horas extra son 3,77 veces los de quienes no las hacen.
• Estado civil: los solteros tienen odds 2,40 veces los de los casados. Los divorciados no difieren significativamente de los casados (su IC contiene el 1).
• Viaje de negocios: los OR frente a quienes no viajan crecen con la frecuencia, de 2,02 (raramente) a 3,81 (frecuentemente).
• Edad y antigüedad en el cargo: \(\hat\beta<0\). Cada año adicional multiplica los odds por 0,949 y 0,864, respectivamente.
• Ingreso: \(\hat\beta<0\) en la escala logarítmica. Duplicar el ingreso multiplica los odds por \(2^{\hat\beta}=0,534\).

6.4 Linealidad en el logit

La regresión logística supone que el logit es lineal en cada covariable cuantitativa. Se evalúa de dos formas: con el logit empírico por deciles, \(\ln\{\bar y_g/(1-\bar y_g)\}\), y con la prueba de Box-Tidwell, que añade el término \(x\ln x\) al modelo de las seis covariables (con \(x\) en la escala en que entra la variable) y lo contrasta con una prueba de razón de verosimilitudes (Hosmer et al., 2013, §4.2). En Antigüedad_Cargo, que tiene ceros, se usa \(x+1\).

le <- bind_rows(lapply(SEL_CUANTI, function(v) {
  g <- cut(datos[[v]], unique(quantile(datos[[v]], 0:10 / 10)), include.lowest = TRUE)
  datos %>% mutate(g = g) %>% group_by(g) %>%
    summarise(x = mean(.data[[v]]), p = mean(y), .groups = "drop") %>%
    mutate(var = v, logit = qlogis(p))
})) %>% mutate(var = factor(var, levels = SEL_CUANTI),
               col = unname(COL_FAMILIA[FAMILIA[as.character(var)]]))
ggplot(le, aes(x, logit)) +
  geom_smooth(method = "lm", se = FALSE, colour = COL_GRIS, linewidth = 0.8, formula = y ~ x) +
  geom_line(aes(colour = col), linewidth = 0.5) + geom_point(aes(colour = col), size = 1.6) +
  scale_colour_identity() +
  facet_wrap(~ var, scales = "free_x", nrow = 1) +
  labs(title = "El riesgo cae rápido en los valores bajos y se aplana después",
       x = NULL, y = "Logit empírico")

Figura 6.3: Logit empírico de la rotación por deciles de cada covariable cuantitativa. Cada punto es un grupo; la línea es un ajuste lineal de referencia.

Prueba de Box-Tidwell

\(H_0\): \(\gamma=0\) en \(\operatorname{logit}\pi=\mathbf x^\top\boldsymbol\beta+\gamma\,z\ln z\), es decir, el logit es lineal en \(z\); \(H_1\): \(\gamma\neq0\) (el logit es curvo en \(z\)). El estadístico es \(G^2=D_{\text{sin }z\ln z}-D_{\text{con }z\ln z}\overset{a}{\sim}\chi^2_1\).

Por qué esta prueba. La linealidad en el logit es el supuesto de forma de la regresión logística para las covariables cuantitativas. El término \(z\ln z\) es la derivada de \(z^\lambda\) respecto a \(\lambda\) en \(\lambda=1\), así que la prueba detecta si una potencia distinta de 1 ajusta mejor (Box y Tidwell, 1962). Requiere \(z>0\); por eso en la antigüedad en el cargo se usa \(x+1\).

# Comparación de la escala del ingreso en el modelo completo
m_ing_lin <- update(modelo, . ~ . - log(Ingreso_Mensual) + Ingreso_Mensual)
aic_lin <- AIC(m_ing_lin); aic_log <- AIC(modelo)

# Box-Tidwell: se añade z ln z al modelo base que contiene a z en esa escala
bt <- function(v, desplaza = 0, escala = identity, nombre = v, base = modelo) {
  z <- escala(datos[[v]] + desplaza)
  d2 <- datos %>% mutate(bt_ = z * log(z))
  m2 <- update(base, . ~ . + bt_, data = d2)
  G  <- deviance(base) - deviance(m2)
  data.frame(Variable = nombre,
             `Modelo base` = if (identical(base, modelo)) "ln(ingreso)" else "ingreso lineal",
             `Chi-cuadrado LR` = G, `Valor p` = pchisq(G, 1, lower.tail = FALSE),
             Decisión = decide(pchisq(G, 1, lower.tail = FALSE)), check.names = FALSE)
}
bt_tab <- bind_rows(bt("Edad"),
                    bt("Ingreso_Mensual", nombre = "Ingreso_Mensual (escala original)", base = m_ing_lin),
                    bt("Ingreso_Mensual", escala = log, nombre = "ln(Ingreso_Mensual)"),
                    bt("Antigüedad_Cargo", desplaza = 1, nombre = "Antigüedad_Cargo + 1"))
tabla(bt_tab %>% mutate(`Valor p` = pval(`Valor p`)),
      "Prueba de Box-Tidwell en el modelo de las seis covariables.",
      digits = c(0, 0, 3, 0))

Tabla 6.6: Prueba de Box-Tidwell en el modelo de las seis covariables.

Variable	Modelo base	Chi-cuadrado LR	Valor p	Decisión
Edad	ln(ingreso)	8,566	0,003	Se rechaza \(H_0\)
Ingreso_Mensual (escala original)	ingreso lineal	7,269	0,007	Se rechaza \(H_0\)
ln(Ingreso_Mensual)	ln(ingreso)	6,080	0,014	Se rechaza \(H_0\)
Antigüedad_Cargo + 1	ln(ingreso)	8,050	0,005	Se rechaza \(H_0\)

Decisión metodológica

El ingreso entra en escala logarítmica. Con el ingreso en su escala original el modelo tiene AIC \(=1.100,4\); con \(\ln(\text{Ingreso})\), AIC \(=1.092,6\), con el mismo número de parámetros. La escala logarítmica también corresponde a la asimetría de la variable y permite una lectura en términos relativos (efecto de duplicar el ingreso).

La linealidad es una aproximación. Box-Tidwell rechaza la linealidad de la edad, de la antigüedad en el cargo y del ingreso en las dos escalas (Tabla 6.6); la escala logarítmica reduce la curvatura del ingreso pero no la elimina. El logit empírico muestra la forma de cada variable (Figura 6.3). En el ingreso y la antigüedad en el cargo la tendencia global es decreciente, con irregularidades locales; en la edad, el logit empírico baja hasta cerca de los 40 años y luego se estabiliza o sube, así que la relación no es monótona y el coeficiente lineal resume solo la tendencia promedio. Se conserva la forma lineal para que los coeficientes sean interpretables, como pide el enunciado, y en la Sección 7.6 se ajusta un modelo con splines para verificar que las conclusiones no dependen de esta decisión.

6.5 Contraste con las hipótesis

dir_obs <- c(
  Horas_Extra      = sign(sl$beta[sl$Coeficiente == "Horas extra: Sí vs. No"]),
  Estado_Civil     = sign(sl$beta[sl$Coeficiente == "Estado civil: Soltero vs. Casado"]),
  Viaje_Negocios   = sign(glm(y ~ as.integer(Viaje_Negocios), binomial, datos)$coef[2]),
  Edad             = sign(mw_tab$`r biserial`[1]),
  Ingreso_Mensual  = sign(mw_tab$`r biserial`[2]),
  Antigüedad_Cargo = sign(mw_tab$`r biserial`[3]))
p_princ <- c(chi_tab$`Valor p`[1:2], ca$p.value, mw_tab$`Valor p`)
b_tend  <- coef(glm(y ~ as.integer(Viaje_Negocios), binomial, datos))[2]
efecto  <- c(paste0("V = ", num(chi_tab$`V de Cramér`[1:2], 3)),
             paste0("OR por nivel = ", num(exp(b_tend), 2)),
             paste0("r = ", num(mw_tab$`r biserial`, 3)))
veredicto <- function(p, s_obs, s_esp) ifelse(p >= 0.05, "No confirmada",
  ifelse(s_obs == s_esp, "Confirmada", "Significativa en sentido contrario"))
contraste <- data.frame(
  Variable = HIP$Variable,
  Prueba = c("Chi-cuadrado", "Chi-cuadrado", "Cochran-Armitage", rep("Mann-Whitney", 3)),
  `Dirección esperada` = HIP$`Dirección esperada`,
  `Dirección observada` = c(
    paste0("OR = ", num(sl$OR[sl$Coeficiente == "Horas extra: Sí vs. No"], 2)),
    paste0("OR (Soltero) = ", num(sl$OR[sl$Coeficiente == "Estado civil: Soltero vs. Casado"], 2)),
    ifelse(dir_obs[3] > 0, "Creciente", "Decreciente"),
    ifelse(dir_obs[4:6] < 0, "Negativa", "Positiva")),
  `Valor p` = pval(p_princ), `p (Holm)` = pval(p.adjust(p_princ, "holm")),
  Efecto = efecto,
  Veredicto = veredicto(p_princ, dir_obs, HIP$Signo), check.names = FALSE)
tabla(contraste, "Contraste de las hipótesis con los resultados del análisis bivariado.",
      font = 7.5) %>% column_spec(2, width = "2cm") %>% column_spec(3, width = "2cm") %>%
  column_spec(4, width = "2.2cm") %>% column_spec(8, width = "1.7cm")

Tabla 6.7: Contraste de las hipótesis con los resultados del análisis bivariado.

Variable	Prueba	Dirección esperada	Dirección observada	Valor p	p (Holm)	Efecto	Veredicto
Horas_Extra	Chi-cuadrado	Sí > No	OR = 3,77	\(\lt 2\times10^{-16}\)	\(\lt 2\times10^{-16}\)	V = 0,246	Confirmada
Estado_Civil	Chi-cuadrado	Soltero > Casado	OR (Soltero) = 2,40	\(9,46\times10^{-11}\)	\(1,89\times10^{-10}\)	V = 0,177	Confirmada
Viaje_Negocios	Cochran-Armitage	Creciente con la frecuencia	Creciente	\(1,12\times10^{-6}\)	\(1,12\times10^{-6}\)	OR por nivel = 1,92	Confirmada
Edad	Mann-Whitney	Negativa	Negativa	\(5,28\times10^{-11}\)	\(1,58\times10^{-10}\)	r = -0,269	Confirmada
Ingreso_Mensual	Mann-Whitney	Negativa	Negativa	\(2,95\times10^{-14}\)	\(1,48\times10^{-13}\)	r = -0,311	Confirmada
Antigüedad_Cargo	Mann-Whitney	Negativa	Negativa	\(4,43\times10^{-12}\)	\(1,77\times10^{-11}\)	r = -0,280	Confirmada

n_conf_biv <- sum(contraste$Veredicto == "Confirmada")
p_solt  <- summary(glm(y ~ Estado_Civil, binomial, datos))$coefficients["Estado_CivilSoltero", 4]
r_solt  <- suppressWarnings(chisq.test(table(datos$Estado_Civil, datos$y), correct = FALSE))$stdres["Soltero", "1"]

Las 6 hipótesis se confirman en el análisis bivariado (Tabla 6.7). La decisión principal usa \(\alpha=0{,}05\) sin ajuste, porque cada hipótesis tiene un mecanismo y una dirección definidos a priori; ahora bien, la selección de las covariables se hizo después del análisis exploratorio (Sección 3.5), así que las seis hipótesis no son estrictamente preespecificadas y su confirmación no es independiente de esa selección; los valores \(p\) ajustados por Holm (1979), que controlan la tasa de error por familia de los seis contrastes, llevan a la misma decisión. En Estado_Civil la prueba \(\chi^2\) es ómnibus (2 gl); la hipótesis direccional se verifica con el contraste específico de solteros frente a casados, que en la regresión simple tiene \(p=2,54\times10^{-8}\) (Wald), y con el residuo ajustado de la celda Soltero–rotó, \(r=6,73\). Para los divorciados no se planteó dirección.

7 Modelo de regresión logística

7.1 Especificación y estimación

Sea \(Y_i\in\{0,\,1\}\) la rotación del empleado \(i\) y \(\pi_i=P(Y_i=1\mid\mathbf x_i)\). Las \(Y_i\) se suponen Bernoulli independientes con \[ \operatorname{logit}\pi_i=\ln\frac{\pi_i}{1-\pi_i}=\mathbf x_i^\top\boldsymbol\beta , \qquad \pi_i=\frac{1}{1+e^{-\mathbf x_i^\top\boldsymbol\beta}}, \] donde \(\mathbf x_i\) contiene el intercepto, las indicadoras de Horas_Extra (Sí), Estado_Civil (Divorciado, Soltero) y Viaje_Negocios (Raramente, Frecuentemente), y las variables Edad, \(\ln(\texttt{Ingreso\_Mensual})\) y Antigüedad_Cargo. El estimador de máxima verosimilitud maximiza \[ \ell(\boldsymbol\beta)=\sum_{i=1}^{n}\left[y_i\ln\pi_i+(1-y_i)\ln(1-\pi_i)\right], \] cuyo gradiente \(\nabla\ell=\mathbf X^\top(\mathbf y-\boldsymbol\pi)\) no tiene raíz en forma cerrada; se resuelve por mínimos cuadrados iterativamente reponderados (IRLS), equivalente al método de Newton-Raphson con el enlace canónico (Agresti, 2013, §4.6). Con enlace logit, \(e^{\beta_j}\) es el OR asociado a un aumento unitario de \(x_j\) con las demás covariables fijas.

Decisión metodológica

El modelo inferencial se estima con la base completa (\(n=1.470\)), porque es el que se interpreta y el que usa toda la información disponible. La partición en entrenamiento y prueba se reserva para evaluar el poder predictivo (Sección 8). La relación entre eventos y parámetros es \(\text{EPV}=237/8=29,6\), por encima del umbral de 10 de Peduzzi et al. (1996).

Pruebas sobre los coeficientes del modelo

• Wald, para cada coeficiente. \(H_0\): \(\beta_j=0\); \(H_1\): \(\beta_j\neq0\). El estadístico \(z_j=\hat\beta_j/\widehat{\mathrm{EE}}(\hat\beta_j)\) sigue aproximadamente una \(N(0,\,1)\) bajo \(H_0\), por la normalidad asintótica del estimador de máxima verosimilitud.
• Razón de verosimilitudes tipo II, para cada término. \(H_0\): todos los coeficientes del término son nulos, dados los demás términos; \(H_1\): alguno no lo es. \(G^2=D_{\text{sin el término}}-D_{\text{completo}}\overset{a}{\sim}\chi^2_{\text{gl}}\), con gl igual al número de coeficientes del término.
• Razón de verosimilitudes global. \(H_0\): \(\beta_1=\dots=\beta_8=0\); \(H_1\): algún \(\beta_j\neq0\). \(G^2=D_0-D\overset{a}{\sim}\chi^2_8\).

Por qué estas pruebas. La de Wald contrasta cada nivel por separado, pero un factor de tres niveles necesita una prueba conjunta de sus dos coeficientes. Para ella se usa la razón de verosimilitudes, que además es más fiable que la de Wald cuando un coeficiente es grande (Hauck y Donner, 1977). Se usa el tipo II porque el modelo no tiene interacciones: cada término se contrasta después de los demás y el resultado no depende del orden en que entran. Los IC de los OR se obtienen por verosimilitud perfilada, coherentes con la prueba de razón de verosimilitudes. Se rechaza \(H_0\) si \(p<0{,}05\).

cf  <- summary(modelo)$coefficients
icp <- suppressMessages(confint(modelo))       # IC por verosimilitud perfilada
coef_tab <- data.frame(
  Término = etiqueta(rownames(cf)), beta = cf[, 1], EE = cf[, 2], z = cf[, 3],
  `p (Wald)` = pval(cf[, 4]), OR = exp(cf[, 1]),
  `IC 2,5 %` = exp(icp[rownames(cf), 1]), `IC 97,5 %` = exp(icp[rownames(cf), 2]),
  check.names = FALSE, row.names = NULL)
tabla(coef_tab %>% rename(`$\\hat\\beta$` = beta),
      "Coeficientes del modelo logístico: log-odds, error estándar, prueba de Wald, OR e IC 95 % por verosimilitud perfilada.",
      digits = c(0, 4, 4, 3, 0, 3, 3, 3), font = 8.5)

Tabla 7.1: Coeficientes del modelo logístico: log-odds, error estándar, prueba de Wald, OR e IC 95 % por verosimilitud perfilada.

Término	\(\hat\beta\)	EE	z	p (Wald)	OR	IC 2,5 %	IC 97,5 %
Intercepto	3,0273	1,2014	2,520	0,012	20,641	1,990	222,147
Horas extra: Sí vs. No	1,4685	0,1593	9,216	\(\lt 2\times10^{-16}\)	4,343	3,183	5,948
Estado civil: Divorciado vs. Casado	-0,2968	0,2305	-1,288	0,198	0,743	0,468	1,158
Estado civil: Soltero vs. Casado	0,7920	0,1721	4,602	\(4,18\times10^{-6}\)	2,208	1,578	3,101
Viaje: Raramente vs. No viaja	0,6800	0,3317	2,050	0,040	1,974	1,069	3,962
Viaje: Frecuentemente vs. No viaja	1,3333	0,3542	3,764	\(1,67\times10^{-4}\)	3,793	1,955	7,914
Edad (por año)	-0,0239	0,0100	-2,384	0,017	0,976	0,957	0,995
ln(Ingreso mensual)	-0,6036	0,1547	-3,902	\(9,55\times10^{-5}\)	0,547	0,402	0,739
Antigüedad en el cargo (por año)	-0,0884	0,0277	-3,187	0,001	0,915	0,866	0,966

b <- coef(modelo); OR <- exp(b)
ORic <- exp(icp)

av <- Anova(modelo, type = "II", test.statistic = "LR")
av_tab <- data.frame(Término = c("Horas_Extra", "Estado_Civil", "Viaje_Negocios", "Edad",
                                 "ln(Ingreso_Mensual)", "Antigüedad_Cargo"),
                     `Chi-cuadrado LR` = av$`LR Chisq`, gl = av$Df,
                     `Valor p` = av$`Pr(>Chisq)`, check.names = FALSE, row.names = NULL) %>%
  mutate(`Chi-cuadrado - gl` = `Chi-cuadrado LR` - gl) %>%
  arrange(desc(`Chi-cuadrado - gl`)) %>% mutate(Rango = row_number(), Decisión = decide(`Valor p`))
tabla(av_tab %>% mutate(`Valor p` = pval(`Valor p`)),
      "Pruebas de razón de verosimilitudes tipo II por término, ordenadas por importancia.",
      digits = c(0, 3, 0, 0, 3, 0))

Tabla 7.2: Pruebas de razón de verosimilitudes tipo II por término, ordenadas por importancia.

Término	Chi-cuadrado LR	gl	Valor p	Chi-cuadrado - gl	Rango	Decisión
Horas_Extra	86,470	1	\(\lt 2\times10^{-16}\)	85,470	1	Se rechaza \(H_0\)
Estado_Civil	32,573	2	\(8,45\times10^{-8}\)	30,573	2	Se rechaza \(H_0\)
Viaje_Negocios	20,132	2	\(4,25\times10^{-5}\)	18,132	3	Se rechaza \(H_0\)
ln(Ingreso_Mensual)	15,696	1	\(7,44\times10^{-5}\)	14,696	4	Se rechaza \(H_0\)
Antigüedad_Cargo	10,737	1	0,001	9,737	5	Se rechaza \(H_0\)
Edad	5,881	1	0,015	4,881	6	Se rechaza \(H_0\)

lr_glob <- modelo$null.deviance - modelo$deviance
gl_glob <- modelo$df.null - modelo$df.residual
p_glob  <- pchisq(lr_glob, gl_glob, lower.tail = FALSE)
mcf     <- 1 - as.numeric(logLik(modelo)) / as.numeric(logLik(update(modelo, . ~ 1)))

Significancia global. La prueba de razón de verosimilitudes frente al modelo nulo da \(G^2=224,00\) con 8 grados de libertad (\(p<2\times10^{-16}\)): al menos una covariable está asociada con la rotación. El pseudo-\(R^2\) de McFadden es \(0,1725\); en la escala de la log-verosimilitud, valores entre 0,2 y 0,4 se consideran ya un ajuste muy bueno (McFadden, 1979), así que el valor obtenido es moderado.

Significancia por término. Las pruebas de Wald (Tabla 7.1) evalúan cada coeficiente por separado; para los factores de tres niveles se requiere la prueba conjunta de razón de verosimilitudes tipo II, que compara el modelo con y sin el término (Tabla 7.2). Los seis términos son significativos al 5 %. Dentro de Estado_Civil, solo el contraste Soltero vs. Casado es significativo (\(p=4,18\times10^{-6}\)); el de Divorciado vs. Casado no lo es (\(p=0,198\)).

7.2 Interpretación de los coeficientes

Todas las interpretaciones se hacen con las demás covariables constantes y describen asociaciones, no efectos causales.

• Horas extra (\(\hat\beta=1,468\)): los odds de rotar de quien hace horas extra son 4,34 veces los de quien no las hace (IC 95 %: 3,18–5,95).
• Estado civil: los odds de un soltero son 2,21 veces los de un casado (IC: 1,58–3,10). El OR de los divorciados, 0,74, tiene un IC que contiene el 1.
• Viaje de negocios: frente a quien no viaja, los odds se multiplican por 1,97 si viaja raramente y por 3,79 si viaja con frecuencia. El gradiente respalda la hipótesis de tendencia creciente.
• Edad (\(\hat\beta=-0,0239\)): cada año adicional multiplica los odds por 0,9764, una reducción de 2,4 %; cinco años más los multiplican por \(e^{5\hat\beta}=0,887\).
• Ingreso (\(\hat\beta=-0,6036\) en escala logarítmica): como \(\ln(k\,x)=\ln k+\ln x\), multiplicar el ingreso por \(k\) multiplica los odds por \(k^{\hat\beta}\). Duplicarlo los multiplica por \(2^{\hat\beta}=0,658\), y aumentarlo un 10 %, por \(1{,}1^{\hat\beta}=0,944\).
• Antigüedad en el cargo (\(\hat\beta=-0,0884\)): cada año adicional multiplica los odds por 0,915, una reducción de 8,5 %.
• Intercepto (\(\hat\beta_0=3,027\)): es el log-odds de un casado sin horas extra que no viaja, con edad 0, ingreso 1 y antigüedad 0. Ese perfil está fuera del rango de los datos, así que el intercepto no tiene lectura sustantiva; solo fija el nivel de la curva.

7.3 Jerarquía de los factores

El enunciado pide identificar qué factores inciden “en mayor proporción”. Comparar coeficientes crudos no sirve, porque dependen de las unidades de cada variable y los factores de tres niveles tienen dos coeficientes. Se usan dos criterios complementarios. El primero es \(\chi^2_{LR}-\text{gl}\) por término (Harrell, 2015, §5.4), que descuenta el valor esperado del estadístico bajo \(H_0\) y así corrige la ventaja de los términos con más grados de libertad (Tabla 7.2). El segundo expresa las cuantitativas como el OR entre el tercer y el primer cuartil, para ponerlas en una escala comparable con los contrastes entre niveles de los factores.

or_ric <- bind_rows(lapply(SEL_CUANTI, function(v) {
  q <- quantile(datos[[v]], c(.25, .75))
  cb <- if (v == "Ingreso_Mensual") "log(Ingreso_Mensual)" else v
  delta <- if (v == "Ingreso_Mensual") diff(log(q)) else diff(q)
  data.frame(Variable = v, Q1 = q[1], Q3 = q[2],
             `OR (Q3 vs. Q1)` = exp(b[cb] * delta),
             `IC 2,5 %` = exp(icp[cb, 1] * delta), `IC 97,5 %` = exp(icp[cb, 2] * delta),
             check.names = FALSE, row.names = NULL)
}))
tabla(or_ric, "OR entre el tercer y el primer cuartil de las covariables cuantitativas.",
      digits = c(0, 0, 0, 3, 3, 3))

Tabla 7.3: OR entre el tercer y el primer cuartil de las covariables cuantitativas.

Variable	Q1	Q3	OR (Q3 vs. Q1)	IC 2,5 %	IC 97,5 %
Edad	30	43	0,733	0,565	0,943
Ingreso_Mensual	2.911	8.379	0,528	0,382	0,726
Antigüedad_Cargo	2	7	0,643	0,487	0,840

Por \(\chi^2_{LR}-\text{gl}\) el orden es: Horas_Extra, Estado_Civil, Viaje_Negocios, ln(Ingreso_Mensual), Antigüedad_Cargo, Edad. Las horas extra son, con amplia diferencia, el factor de mayor peso. Entre las cuantitativas, pasar del primer al tercer cuartil del ingreso multiplica los odds por 0,528, y de la antigüedad en el cargo, por 0,643 (Tabla 7.3). El IC del OR por RIC se obtiene escalando el IC perfilado de \(\beta\), lo cual es exacto porque la transformación es monótona.

7.4 Comparación entre el análisis bivariado y el modelo ajustado

comp <- data.frame(Término = etiqueta(names(b)[-1]),
                   `OR simple` = sl$OR[match(etiqueta(names(b)[-1]), sl$Coeficiente)],
                   `OR ajustado` = OR[-1],
                   `Cambio (%)` = 100 * (OR[-1] / sl$OR[match(etiqueta(names(b)[-1]), sl$Coeficiente)] - 1),
                   `p (Wald, ajustado)` = pval(cf[-1, 4]),
                   check.names = FALSE, row.names = NULL)
tabla(comp, "OR de las regresiones simples frente a los OR del modelo ajustado.",
      digits = c(0, 3, 3, 1, 0))

Tabla 7.4: OR de las regresiones simples frente a los OR del modelo ajustado.

Término	OR simple	OR ajustado	Cambio (%)	p (Wald, ajustado)
Horas extra: Sí vs. No	3,771	4,343	15,1	\(\lt 2\times10^{-16}\)
Estado civil: Divorciado vs. Casado	0,787	0,743	-5,6	0,198
Estado civil: Soltero vs. Casado	2,404	2,208	-8,2	\(4,18\times10^{-6}\)
Viaje: Raramente vs. No viaja	2,023	1,974	-2,4	0,040
Viaje: Frecuentemente vs. No viaja	3,815	3,793	-0,6	\(1,67\times10^{-4}\)
Edad (por año)	0,949	0,976	2,9	0,017
ln(Ingreso mensual)	0,404	0,547	35,2	\(9,55\times10^{-5}\)
Antigüedad en el cargo (por año)	0,864	0,915	6,0	0,001

cambia_signo <- sum(sign(log(comp$`OR simple`)) != sign(log(comp$`OR ajustado`)))

Ningún coeficiente cambia de signo al ajustar (0 cambios; Tabla 7.4). En la escala de \(\beta\), los mayores cambios están en la edad, la antigüedad en el cargo y el ingreso, cuyos coeficientes se acercan a 0: es compatible con confusión parcial, porque están correlacionadas entre sí (\(r_s=0,472\) entre la edad y el ingreso). El cambio no se explica solo por confusión: el OR no es colapsable, así que el OR ajustado puede alejarse de 1 sin que haya confusión, como ocurre con las horas extra (Greenland, Robins y Pearl, 1999). La edad sigue siendo significativa en el modelo ajustado, aunque con el menor peso de los seis términos.

7.5 Diagnóstico del modelo

gv <- car::vif(modelo)   # rms también exporta vif(); se usa la de car (GVIF)
tabla(data.frame(Término = rownames(gv), GVIF = gv[, 1], gl = gv[, 2],
                 `$\\mathrm{GVIF}^{1/(2\\,\\mathrm{gl})}$` = gv[, 3], check.names = FALSE, row.names = NULL) %>%
        mutate(Término = sub("log\\(Ingreso_Mensual\\)", "ln(Ingreso_Mensual)", Término)),
      "Factores de inflación de varianza generalizados.", digits = c(0, 3, 0, 3))

Tabla 7.5: Factores de inflación de varianza generalizados.

Término	GVIF	gl	\(\mathrm{GVIF}^{1/(2\,\mathrm{gl})}\)
Horas_Extra	1,036	1	1,018
Estado_Civil	1,030	2	1,007
Viaje_Negocios	1,013	2	1,003
Edad	1,269	1	1,126
ln(Ingreso_Mensual)	1,403	1	1,185
Antigüedad_Cargo	1,154	1	1,074

Colinealidad. Para términos con varios grados de libertad se usa el GVIF de Fox y Monette (1992), y \(\text{GVIF}^{1/(2\,\text{gl})}\) se compara con \(\sqrt{5}\approx2{,}24\), el equivalente del umbral \(\text{VIF}=5\) para un solo coeficiente. El máximo es 1,185 (Tabla 7.5): no hay colinealidad preocupante.

celdas_min <- min(sapply(SEL_CUALI, function(v) min(table(datos[[v]], datos$y))))
max_ee <- max(cf[-1, 2])

Separación. Hay separación completa o cuasicompleta cuando una combinación lineal de covariables predice perfectamente la respuesta; entonces el estimador de máxima verosimilitud no existe y los errores estándar crecen sin límite. El algoritmo convergió en 5 iteraciones, la celda menos poblada de los cruces entre factores y respuesta tiene 12 casos y el mayor error estándar de las pendientes es 0,354. No hay indicios de separación.

Prueba de Hosmer-Lemeshow

\(H_0\): el modelo está bien calibrado, \(E(Y\mid\mathbf x)=\pi(\mathbf x;\boldsymbol\beta)\); \(H_1\): no lo está. Las observaciones se ordenan por \(\hat\pi\) y se agrupan en \(g\) grupos de tamaño similar; el estadístico es \[ \hat C=\sum_{k=1}^{g}\frac{(O_k-E_k)^2}{E_k\,(1-E_k/n_k)}\ \overset{a}{\sim}\ \chi^2_{g-2}, \] donde \(O_k\) y \(E_k\) son las rotaciones observadas y esperadas en el grupo \(k\), y \(n_k\) su tamaño.

Por qué esta prueba. Con covariables continuas casi cada empleado tiene un patrón de covariables único, y entonces la devianza y la \(\chi^2\) de Pearson no tienen una distribución \(\chi^2\) de referencia; agrupar por \(\hat\pi\) resuelve ese problema. Como el resultado depende de \(g\), se reporta con tres valores.

hosmer <- function(y, p, g) {
  cortes <- unique(quantile(p, probs = seq(0, 1, length.out = g + 1)))
  grp <- cut(p, cortes, include.lowest = TRUE)
  o1 <- tapply(y, grp, sum); e1 <- tapply(p, grp, sum); ng <- tapply(y, grp, length)
  C <- sum((o1 - e1)^2 / (e1 * (1 - e1 / ng)))
  gg <- nlevels(grp)
  data.frame(Grupos = gg, `C de Hosmer-Lemeshow` = C, gl = gg - 2,
             `Valor p` = pchisq(C, gg - 2, lower.tail = FALSE),
             Decisión = decide(pchisq(C, gg - 2, lower.tail = FALSE)), check.names = FALSE)
}
hl_tab <- bind_rows(lapply(c(8, 10, 12), function(g) hosmer(datos$y, fitted(modelo), g)))
tabla(hl_tab %>% mutate(`Valor p` = pval(`Valor p`)),
      "Prueba de Hosmer-Lemeshow con distinto número de grupos.", digits = c(0, 3, 0, 0))

Tabla 7.6: Prueba de Hosmer-Lemeshow con distinto número de grupos.

Grupos	C de Hosmer-Lemeshow	gl	Valor p	Decisión
8	21,776	6	0,001	Se rechaza \(H_0\)
10	20,951	8	0,007	Se rechaza \(H_0\)
12	24,480	10	0,006	Se rechaza \(H_0\)

Bondad de ajuste. El estadístico de Hosmer-Lemeshow agrupa las observaciones por cuantiles de \(\hat\pi\) y compara eventos observados y esperados; bajo \(H_0\) sigue aproximadamente una \(\chi^2_{g-2}\) (Hosmer et al., 2013, §5.2.2). Su resultado depende del número de grupos, por eso se reporta con \(g\in\{8,\,10,\,12\}\) (Tabla 7.6). Con \(g=10\), \(p=0,007\). La prueba rechaza el ajuste con los tres números de grupos.

cb10 <- unique(quantile(fitted(modelo), seq(0, 1, 0.1)))
g10  <- cut(fitted(modelo), cb10, include.lowest = TRUE)
hl_grp <- data.frame(Decil = seq_len(nlevels(g10)), n = as.integer(table(g10)),
                     `Observados` = as.numeric(tapply(datos$y, g10, sum)),
                     `Esperados` = as.numeric(tapply(fitted(modelo), g10, sum)),
                     check.names = FALSE) %>%
  mutate(`Observados - esperados` = Observados - Esperados)

Las rotaciones observadas y esperadas por decil (Anexo B) muestran de dónde viene la falta de ajuste: el modelo subestima la rotación en los deciles extremos y la sobreestima en los intermedios. Es el patrón que produce un logit que no es lineal en las covariables, y coincide con lo que mostró Box-Tidwell. La prueba tiene además mucha potencia con \(n=1.470\), así que detecta desvíos de magnitud moderada (Hosmer et al., 2013, §5.2.2). La magnitud práctica de la descalibración se examina con la curva de calibración de la Sección 8, y el modelo con splines (Sección 7.6) indica cuánto se corrige al relajar la linealidad.

cook <- cooks.distance(modelo); hat <- hatvalues(modelo); dfb <- dfbetas(modelo)
p_c  <- length(b)
infl <- data.frame(
  Medida = c("Distancia de Cook $> 4/n$", "Distancia de Cook $> 1$",
             "Apalancamiento $> 2p/n$", "Algún $|\\mathrm{DFBETAS}| > 2/\\sqrt{n}$"),
  Umbral = c(4 / N, 1, 2 * p_c / N, 2 / sqrt(N)),
  Observaciones = c(sum(cook > 4 / N), sum(cook > 1), sum(hat > 2 * p_c / N),
                    sum(apply(abs(dfb) > 2 / sqrt(N), 1, any))))
dfb_tab <- data.frame(Término = etiqueta(colnames(dfb)),
                      `Máximo |DFBETAS|` = apply(abs(dfb), 2, max),
                      `Observaciones con |DFBETAS| > 2/raíz(n)` = colSums(abs(dfb) > 2 / sqrt(N)),
                      check.names = FALSE, row.names = NULL)
rot_cook <- mean(datos$y[cook > 4 / N] == 1)
pi_cook  <- median(fitted(modelo)[cook > 4 / N & datos$y == 1])
alto_h   <- hat > 2 * p_c / N

Influencia. Ninguna observación tiene distancia de Cook mayor que 1 (la máxima es 0,0158), el umbral usual de influencia grave. Con el umbral más exigente, \(4/n\), se marcan 115 observaciones; el 95,7 % son rotaciones de bajo riesgo estimado, cuyo residuo es grande por construcción en un evento poco frecuente. El máximo \(|\mathrm{DFBETAS}|\) es 0,268: ninguna observación desplaza un coeficiente ni un tercio de su error estándar, y se conservan todas. El detalle está en el Anexo C.

7.6 Sensibilidad a la forma funcional

modelo_ns <- glm(y ~ Horas_Extra + Estado_Civil + Viaje_Negocios +
                   splines::ns(Edad, 3) + splines::ns(log(Ingreso_Mensual), 3) +
                   splines::ns(Antigüedad_Cargo, 3), binomial, datos)
lr_nl <- deviance(modelo) - deviance(modelo_ns)
gl_nl <- modelo$df.residual - modelo_ns$df.residual
p_nl  <- pchisq(lr_nl, gl_nl, lower.tail = FALSE)
auc_lin <- as.numeric(auc(roc(datos$y, fitted(modelo), levels = c(0, 1), direction = "<", quiet = TRUE)))
auc_ns  <- as.numeric(auc(roc(datos$y, fitted(modelo_ns), levels = c(0, 1), direction = "<", quiet = TRUE)))
hl_ns <- bind_rows(lapply(c(8, 10, 12), function(g) hosmer(datos$y, fitted(modelo_ns), g)))
base_pred <- data.frame(Horas_Extra = factor("No", levels(datos$Horas_Extra)),
                        Estado_Civil = factor("Casado", levels(datos$Estado_Civil)),
                        Viaje_Negocios = factor("Raramente", levels(datos$Viaje_Negocios)),
                        Edad = median(datos$Edad), Ingreso_Mensual = median(datos$Ingreso_Mensual),
                        Antigüedad_Cargo = median(datos$Antigüedad_Cargo))
curvas <- bind_rows(lapply(SEL_CUANTI, function(v) {
  r <- quantile(datos[[v]], c(0.02, 0.98))
  g <- base_pred[rep(1, 80), ]; g[[v]] <- seq(r[1], r[2], length.out = 80)
  bind_rows(data.frame(var = v, x = g[[v]], eta = predict(modelo, g), Modelo = "Lineal"),
            data.frame(var = v, x = g[[v]], eta = predict(modelo_ns, g), Modelo = "Splines"))
})) %>% mutate(var = factor(var, levels = SEL_CUANTI))
edad_min <- curvas %>% filter(var == "Edad", Modelo == "Splines") %>%
  slice_min(eta, n = 1) %>% pull(x)

Como Box-Tidwell rechazó la linealidad, se reajustó el modelo reemplazando cada cuantitativa por un spline cúbico natural con 3 grados de libertad (Anexo D). Los splines mejoran el ajuste (\(G^2=24,91\), 6 gl, \(p=3,55\times10^{-4}\)), pero la ganancia en discriminación es pequeña: el AUC aparente pasa de 0,778 a 0,785. El ingreso y la antigüedad en el cargo conservan un efecto decreciente; la edad tiene forma de U, con mínimo cerca de los 43 años. El modelo lineal se mantiene como principal por su interpretabilidad; su limitación es que promedia un efecto que no es uniforme a lo largo del rango.

8 Evaluación del poder predictivo

8.1 Partición y esquema de validación

La base se divide de forma estratificada por \(Y\) en entrenamiento (70 %) y prueba (30 %), con semilla 2026, para que las dos partes conserven la prevalencia de rotación. El modelo de la Sección 7 se reajusta solo con entrenamiento y se evalúa en prueba, que el modelo no vio.

tabla(data.frame(Conjunto = c("Entrenamiento", "Prueba", "Total"),
                 n = c(nrow(entren), nrow(prueba), N),
                 Rotaciones = c(sum(entren$y), sum(prueba$y), N1),
                 Prevalencia = c(mean(entren$y), mean(prueba$y), N1 / N)),
      "Partición estratificada de la base.", digits = c(0, 0, 0, 4))

Tabla 8.1: Partición estratificada de la base.

Conjunto	n	Rotaciones	Prevalencia
Entrenamiento	1.029	166	0,1613
Prueba	441	71	0,1610
Total	1.470	237	0,1612

En entrenamiento hay 166 rotaciones, así que \(\text{EPV}=166/8=20,8\), también por encima de 10.

8.2 Curva ROC y AUC

Para un corte \(c\), la sensibilidad es \(P(\hat\pi\ge c\mid Y=1)\) y la especificidad es \(P(\hat\pi<c\mid Y=0)\). La curva ROC grafica la sensibilidad frente a \(1-\text{especificidad}\) al recorrer todos los cortes. El área bajo ella tiene una lectura probabilística directa: \[ \text{AUC}=P\left(\hat\pi_{(1)}>\hat\pi_{(0)}\right)+\tfrac12\,P\left(\hat\pi_{(1)}=\hat\pi_{(0)}\right), \] donde \(\hat\pi_{(1)}\) y \(\hat\pi_{(0)}\) son las probabilidades estimadas de un rotante y de un no rotante elegidos al azar. Es el mismo \(\hat A\) de Mann-Whitney aplicado a \(\hat\pi\), y no depende del corte (Fawcett, 2006).

ETQ_TR <- paste0("Entrenamiento (AUC = ", num(auc(roc_tr), 3), ")")
ETQ_TE <- paste0("Prueba (AUC = ", num(auc(roc_te), 3), ")")
roc_df <- bind_rows(
  data.frame(fpr = 1 - roc_tr$specificities, tpr = roc_tr$sensitivities, Conjunto = ETQ_TR),
  data.frame(fpr = 1 - roc_te$specificities, tpr = roc_te$sensitivities, Conjunto = ETQ_TE)) %>%
  arrange(Conjunto, fpr, tpr)
pt_y <- data.frame(fpr = 1 - m_youden_te$Esp, tpr = m_youden_te$Sens)
ggplot(roc_df, aes(fpr, tpr, colour = Conjunto)) +
  geom_abline(linetype = "dashed", colour = "grey75") +
  geom_step(linewidth = 0.8) +
  geom_point(data = pt_y, aes(fpr, tpr), inherit.aes = FALSE, colour = COL_NAR, size = 2.6) +
  scale_colour_manual(values = setNames(c(COL_TEAL, COL_AZUL), c(ETQ_TR, ETQ_TE)), name = NULL) +
  coord_equal() +
  theme(legend.position = "bottom", legend.text = element_text(size = 8),
        legend.key.size = unit(0.4, "cm")) +
  theme(plot.title.position = "plot") +
  labs(title = "Discriminación aceptable en datos nuevos",
       x = "1 - especificidad", y = "Sensibilidad")

Figura 8.1: Curvas ROC del modelo en entrenamiento y en prueba. El punto marca el corte de Youden elegido en entrenamiento.

ic_tr <- ci.auc(roc_tr, method = "delong")
# Bootstrap de optimismo (Harrell, 2015, §5.3) sobre la base completa
fit_lrm <- lrm(FORMULA, data = datos, x = TRUE, y = TRUE)
set.seed(2026)
val <- validate(fit_lrm, method = "boot", B = 500)
auc_orig   <- 0.5 + val["Dxy", "index.orig"] / 2
slope_corr <- val["Slope", "index.corrected"]
auc_corr <- 0.5 + val["Dxy", "index.corrected"] / 2
optim_auc <- val["Dxy", "optimism"] / 2
tabla(data.frame(
  Estimación = c("AUC en entrenamiento", "AUC en prueba",
                 "AUC aparente (base completa)", "AUC corregido por optimismo (bootstrap, B = 500)"),
  AUC = c(auc(roc_tr), auc(roc_te), auc_orig, auc_corr),
  `IC 95 % (DeLong)` = c(paste0(num(ic_tr[1], 3), " – ", num(ic_tr[3], 3)),
                         paste0(num(ic_te[1], 3), " – ", num(ic_te[3], 3)), "—", "—"),
  check.names = FALSE),
  "Poder discriminante del modelo.", digits = c(0, 4, 0))

Tabla 8.2: Poder discriminante del modelo.

Estimación	AUC	IC 95 % (DeLong)
AUC en entrenamiento	0,7918	0,751 – 0,833
AUC en prueba	0,7402	0,671 – 0,810
AUC aparente (base completa)	0,7776	—
AUC corregido por optimismo (bootstrap, B = 500)	0,7700	—

El AUC en prueba es 0,740, con IC 95 % de DeLong 0,671–0,810 (Tabla 8.2). El IC de DeLong et al. (1988) se basa en la teoría de los estadísticos \(U\) y no supone ninguna distribución para \(\hat\pi\). Según la escala de Hosmer et al. (2013, §5.2.4), la discriminación es aceptable (entre 0,7 y 0,8): si se elige al azar un empleado que rotó y uno que no, el modelo asigna mayor probabilidad al primero en el 74,0 % de los casos.

El AUC de prueba (0,740) es 0,052 menor que el de entrenamiento (0,792). Esa brecha mezcla el sobreajuste con la variabilidad de una sola partición: la muestra de prueba tiene solo 71 rotaciones, y el IC de su AUC (0,671–0,810) contiene al AUC de entrenamiento. Para separar los dos efectos se usa el bootstrap de optimismo, que no depende de una partición particular: se reajusta el modelo en 500 remuestras, se mide cuánto se degrada su AUC en la base original y se descuenta ese optimismo del AUC aparente. El optimismo estimado es 0,0077 y el AUC corregido, 0,770: con 8 parámetros y 237 eventos, el modelo casi no está sobreajustado, y la brecha observada en la partición se debe sobre todo al azar de la división. El AUC corregido es la mejor estimación del desempeño esperado en datos nuevos de la misma población.

8.3 Calibración

La discriminación no garantiza que las probabilidades sean correctas en magnitud, y el corte por costos de la Sección 9 necesita que lo sean. La calibración se evalúa en prueba con el puntaje de Brier, la curva por deciles y la recta de calibración \(\operatorname{logit}P(Y=1)=a+b\,\operatorname{logit}\hat\pi\), cuyo ideal es \(a=0\), \(b=1\).

brier <- function(y, p) mean((p - y)^2)
br_te <- brier(prueba$y, p_te); br_ref <- mean(prueba$y) * (1 - mean(prueba$y))
cal <- glm(prueba$y ~ qlogis(p_te), family = binomial)
cal_int <- glm(prueba$y ~ 1, offset = qlogis(p_te), family = binomial)   # a con b = 1
tabla(data.frame(
  Medida = c("Brier del modelo (prueba)", "Brier de referencia: $\\bar y(1-\\bar y)$",
             "Brier escalado: $1-B/B_{\\mathrm{ref}}$",
             "Pendiente de calibración $b$ (prueba)",
             "Pendiente corregida por optimismo (bootstrap, base completa)",
             "Intercepto de calibración $a$ (con $b=1$)"),
  Valor = c(br_te, br_ref, 1 - br_te / br_ref, coef(cal)[2], slope_corr, coef(cal_int)[1])),
  "Medidas de calibración (muestra de prueba y bootstrap en la base completa).", digits = 4)

Tabla 8.3: Medidas de calibración (muestra de prueba y bootstrap en la base completa).

Medida	Valor
Brier del modelo (prueba)	0,1142
Brier de referencia: \(\bar y(1-\bar y)\)	0,1351
Brier escalado: \(1-B/B_{\mathrm{ref}}\)	0,1542
Pendiente de calibración \(b\) (prueba)	0,7979
Pendiente corregida por optimismo (bootstrap, base completa)	0,9639
Intercepto de calibración \(a\) (con \(b=1\))	0,0147

gcal <- cut(p_te, unique(quantile(p_te, 0:10 / 10)), include.lowest = TRUE)
cal_df <- data.frame(y = prueba$y, p = p_te, g = gcal) %>% group_by(g) %>%
  summarise(pm = mean(p), obs = mean(y), n = n(), k = sum(y), .groups = "drop") %>%
  rowwise() %>%
  mutate(li = prop.test(k, n, correct = FALSE)$conf.int[1],
         ls = prop.test(k, n, correct = FALSE)$conf.int[2]) %>% ungroup()
lim <- max(c(cal_df$ls, cal_df$pm)) * 1.05
ggplot(cal_df, aes(pm, obs)) +
  geom_abline(linetype = "dashed", colour = "grey70") +
  geom_errorbar(aes(ymin = li, ymax = ls), width = 0, colour = COL_GRIS, linewidth = 0.9) +
  geom_line(colour = COL_TEAL, linewidth = 0.7) + geom_point(colour = COL_TEAL, size = 2) +
  coord_equal(xlim = c(0, lim), ylim = c(0, lim)) +
  labs(title = "Las probabilidades estimadas siguen a las observadas",
       x = "Probabilidad estimada (media del decil)", y = "Proporción observada")

Figura 8.2: Curva de calibración en la muestra de prueba: proporción observada de rotación frente a la probabilidad media estimada, por decil de probabilidad, con IC 95 % de Wilson.

El puntaje de Brier en prueba es 0,1142, frente a 0,1351 de un modelo que asignara a todos la prevalencia: una reducción de 15,4 % (Tabla 8.3). La pendiente de calibración es 0,798, menor que 1: en esta partición las predicciones son algo más extremas de lo que deberían, lo que concuerda con la brecha de AUC de la misma partición; el intercepto de calibración, 0,015, indica que el nivel medio de riesgo está bien estimado. La pendiente corregida por optimismo sobre la base completa es 0,964, más cercana a 1 que la de esta partición, lo que es compatible con que el valor de 0,798 refleje sobre todo el azar de la división. Las dos no son estrictamente comparables: la pendiente corregida corresponde al modelo de la base completa (\(n=1.470\)) y la de prueba, al ajustado solo con entrenamiento (\(n=1.029\)), que se espera algo más optimista. Los intervalos de 8 de los 10 deciles contienen la probabilidad estimada (Figura 8.2). La falta de ajuste que detectó Hosmer-Lemeshow en la base completa es, en la escala de la decisión, moderada: las probabilidades son utilizables para fijar cortes por costos.

8.4 Matriz de confusión y métricas de clasificación

El AUC, el puntaje de Brier y la calibración evalúan las probabilidades sin fijar un corte. Para clasificar a un empleado como “rota” o “no rota” hay que fijar uno, y entonces el desempeño se resume en la matriz de confusión: verdaderos positivos (VP, rotantes señalados), falsos positivos (FP, no rotantes señalados), verdaderos negativos (VN) y falsos negativos (FN, rotantes que el modelo deja escapar). La clase positiva es la rotación (\(Y=1\)). Se comparan dos cortes en la muestra de prueba: el convencional de 0,5 y el de Youden (\(c_J=0,25\)), elegido en entrenamiento; la justificación del corte está en la Sección 9.

cm_tab <- function(m, etiqueta) data.frame(
  Corte = c(etiqueta, "", ""), Predicción = c("Rota (1)", "No rota (0)", "Total"),
  `Real: rota (1)` = c(m$cm["VP"], m$cm["FN"], m$cm["VP"] + m$cm["FN"]),
  `Real: no rota (0)` = c(m$cm["FP"], m$cm["VN"], m$cm["FP"] + m$cm["VN"]),
  Total = c(m$cm["VP"] + m$cm["FP"], m$cm["FN"] + m$cm["VN"], sum(m$cm)),
  check.names = FALSE, row.names = NULL)
tabla(bind_rows(cm_tab(met_05, "c = 0,50 (convencional)"),
                cm_tab(met_yt, paste0("c = ", num(c_youden, 2), " (Youden)"))),
      "Matrices de confusión en la muestra de prueba con los dos cortes.") %>%
  row_spec(c(3, 6), italic = TRUE)

Tabla 8.4: Matrices de confusión en la muestra de prueba con los dos cortes.

Corte	Predicción	Real: rota (1)	Real: no rota (0)	Total
c = 0,50 (convencional)	Rota (1)	15	6	21
	No rota (0)	56	364	420
	Total	71	370	441
c = 0,25 (Youden)	Rota (1)	36	55	91
	No rota (0)	35	315	350
	Total	71	370	441

cm_largo <- bind_rows(lapply(list(met_05, met_yt), function(m) {
  et <- if (m$corte == 0.5) "c = 0,50 (convencional)" else paste0("c = ", num(m$corte, 2), " (Youden)")
  data.frame(Corte = et,
             Real = factor(c("Rota", "No rota", "Rota", "No rota"), levels = c("Rota", "No rota")),
             Pred = factor(c("Rota", "Rota", "No rota", "No rota"), levels = c("No rota", "Rota")),
             tipo = c("VP", "FP", "FN", "VN"),
             n = unname(m$cm[c("VP", "FP", "FN", "VN")]))
})) %>% group_by(Corte, Real) %>% mutate(pc = n / sum(n)) %>% ungroup() %>%
  mutate(acierto = tipo %in% c("VP", "VN"), Corte = factor(Corte, levels = unique(Corte)))
ggplot(cm_largo, aes(Real, Pred, fill = acierto, alpha = pc)) +
  geom_tile(colour = "white", linewidth = 1) +
  geom_text(aes(label = paste0(tipo, "\n", ent(n), " (", pct(pc, 0), ")")), alpha = 1,
            size = 3.4, colour = COL_TXT, lineheight = 0.9) +
  scale_fill_manual(values = c(`TRUE` = COL_AZUL, `FALSE` = COL_NAR)) +
  scale_alpha(range = c(0.25, 0.85)) +
  facet_wrap(~ Corte) +
  labs(title = "El corte de Youden recupera más rotaciones a cambio de más falsas alarmas",
       x = "Clase real", y = "Clase predicha") +
  theme(panel.grid = element_blank())

Figura 8.3: Matrices de confusión en la muestra de prueba. Cada celda muestra el conteo y el porcentaje de su columna (clase real): en la columna «Rota», la diagonal es la sensibilidad; en la columna «No rota», la especificidad.

Con el corte de 0,5 el modelo señala solo 21 de los 441 empleados de prueba y deja escapar 56 de las 71 rotaciones (Tabla 8.4, Figura 8.3). Con el corte de Youden detecta 36 rotaciones, 21 más, a cambio de 49 falsas alarmas adicionales.

Pruebas sobre la clasificación

• Exactitud frente a la tasa de no información. \(H_0\): \(\text{Exactitud}\le\text{NIR}\); \(H_1\): \(\text{Exactitud}>\text{NIR}\), donde \(\text{NIR}=\max(\bar y,1-\bar y)\) es la exactitud de asignar a todos la clase mayoritaria. Bajo \(H_0\) el número de aciertos sigue una \(\text{Bin}(n,\text{NIR})\), y se usa la prueba binomial exacta unilateral. Responde a si el clasificador supera a la regla trivial, que con clases desbalanceadas ya tiene una exactitud alta.
• McNemar (1947). \(H_0\): \(P(\text{FP})=P(\text{FN})\), es decir, el clasificador no se equivoca más hacia un lado que hacia el otro; \(H_1\): \(P(\text{FP})\neq P(\text{FN})\). Estadístico \(\chi^2=(|FP-FN|-1)^2/(FP+FN)\overset{a}{\sim}\chi^2_1\). Usa solo las celdas discordantes, porque las concordantes no informan sobre el sesgo; detecta si el corte produce un desequilibrio sistemático entre los dos errores.
• AUC frente a 0,5. \(H_0\): \(\text{AUC}=0{,}5\) (el modelo no discrimina mejor que el azar); \(H_1\): \(\text{AUC}>0{,}5\). Como el AUC es el \(\hat A\) de Mann-Whitney aplicado a \(\hat\pi\), se contrasta con la prueba de Mann-Whitney unilateral sobre las probabilidades estimadas de los dos grupos.

Los intervalos de la exactitud son de Clopper y Pearson (1934), exactos, y los de la sensibilidad, la especificidad y los valores predictivos son de Wilson (1927), adecuados para proporciones cercanas a 0 o a 1. El kappa de Cohen (1960) corrige la exactitud por el acuerdo esperado por azar, \(p_e\), y el coeficiente de Matthews (1975) es la correlación de Pearson entre la clase real y la predicha.

f_ic <- function(campo) function(m) paste0(num(m[[campo]], 3), " [", num(m[[paste0(campo, "_ic")]][1], 3),
                                           "; ", num(m[[paste0(campo, "_ic")]][2], 3), "]")
f_v  <- function(campo, d = 3) function(m) num(m[[campo]], d)
FILAS <- list(
  list("Exactitud (*accuracy*)", "$(VP+VN)/n$", f_ic("acc")),
  list("Tasa de no información (NIR)", "$\\max(\\bar y,\\,1-\\bar y)$", f_v("nir")),
  list("Valor p: exactitud > NIR", "Binomial exacta unilateral", function(m) pval(m$p_nir)),
  list("Kappa de Cohen", "$(\\text{Exact.}-p_e)/(1-p_e)$", f_v("kappa")),
  list("Sensibilidad (*recall*, TVP)", "$VP/(VP+FN)$", f_ic("sens")),
  list("Especificidad (TVN)", "$VN/(VN+FP)$", f_ic("esp")),
  list("Precisión (VPP)", "$VP/(VP+FP)$", f_ic("vpp")),
  list("Valor predictivo negativo (VPN)", "$VN/(VN+FN)$", f_ic("vpn")),
  list("$F_1$", "$2\\,\\text{VPP}\\cdot\\text{Sens}/(\\text{VPP}+\\text{Sens})$", f_v("f1")),
  list("Exactitud balanceada", "$(\\text{Sens}+\\text{Esp})/2$", f_v("bal")),
  list("Índice de Youden $J$", "$\\text{Sens}+\\text{Esp}-1$", f_v("J")),
  list("Correlación de Matthews (MCC)",
       "$\\dfrac{VP\\cdot VN-FP\\cdot FN}{\\sqrt{(VP+FP)(VP+FN)(VN+FP)(VN+FN)}}$", f_v("mcc")),
  list("Tasa de falsos positivos", "$FP/(FP+VN)$", f_v("fpr")),
  list("Tasa de falsos negativos", "$FN/(FN+VP)$", f_v("fnr")),
  list("Razón de verosimilitud positiva $LR^+$", "$\\text{Sens}/(1-\\text{Esp})$", f_v("lrp", 2)),
  list("Razón de verosimilitud negativa $LR^-$", "$(1-\\text{Sens})/\\text{Esp}$", f_v("lrn")),
  list("Prevalencia", "$(VP+FN)/n$", f_v("prev")),
  list("Tasa de detección", "$VP/n$", f_v("detec")),
  list("Proporción señalada", "$(VP+FP)/n$", f_v("senal")),
  list("Valor p de McNemar", "$H_0:\\ P(FP)=P(FN)$", function(m) pval(m$p_mcnemar)))
met_tab <- bind_rows(lapply(FILAS, function(f)
  data.frame(Métrica = f[[1]], Definición = f[[2]],
             a = f[[3]](met_05), b = f[[3]](met_yt), c = f[[3]](met_ytr))))
names(met_tab)[3:5] <- c("Prueba, c = 0,50", paste0("Prueba, c = ", num(c_youden, 2)),
                         paste0("Entrenamiento, c = ", num(c_youden, 2)))
met_tab$Métrica <- gsub("\\*", "", met_tab$Métrica)
tabla(met_tab, "Métricas de clasificación con IC 95 % (Clopper-Pearson para la exactitud; Wilson para las proporciones condicionales).",
      align = c("l", "l", "r", "r", "r")) %>%
  column_spec(3:5, extra_css = "white-space: nowrap;") %>%
  pack_rows("Exactitud global", 1, 4) %>%
  pack_rows("Desempeño por clase", 5, 8) %>%
  pack_rows("Medidas resumen para clases desbalanceadas", 9, 12) %>%
  pack_rows("Errores y razones de verosimilitud", 13, 16) %>%
  pack_rows("Prevalencia y sesgo del clasificador", 17, 20)

Tabla 8.5: Métricas de clasificación con IC 95 % (Clopper-Pearson para la exactitud; Wilson para las proporciones condicionales).

Métrica	Definición	Prueba, c = 0,50	Prueba, c = 0,25	Entrenamiento, c = 0,25
Exactitud global
Exactitud (accuracy)	\((VP+VN)/n\)	0,859 [0,823; 0,890]	0,796 [0,755; 0,833]	0,825 [0,800; 0,848]
Tasa de no información (NIR)	\(\max(\bar y,\,1-\bar y)\)	0,839	0,839	0,839
Valor p: exactitud > NIR	Binomial exacta unilateral	0,135	0,993	0,890
Kappa de Cohen	\((\text{Exact.}-p_e)/(1-p_e)\)	0,273	0,322	0,426
Desempeño por clase
Sensibilidad (recall, TVP)	\(VP/(VP+FN)\)	0,211 [0,132; 0,320]	0,507 [0,393; 0,620]	0,614 [0,539; 0,685]
Especificidad (TVN)	\(VN/(VN+FP)\)	0,984 [0,965; 0,993]	0,851 [0,811; 0,884]	0,866 [0,841; 0,887]
Precisión (VPP)	\(VP/(VP+FP)\)	0,714 [0,500; 0,862]	0,396 [0,301; 0,498]	0,468 [0,403; 0,534]
Valor predictivo negativo (VPN)	\(VN/(VN+FN)\)	0,867 [0,831; 0,896]	0,900 [0,864; 0,927]	0,921 [0,900; 0,938]
Medidas resumen para clases desbalanceadas
\(F_1\)	\(2\,\text{VPP}\cdot\text{Sens}/(\text{VPP}+\text{Sens})\)	0,326	0,444	0,531
Exactitud balanceada	\((\text{Sens}+\text{Esp})/2\)	0,598	0,679	0,740
Índice de Youden \(J\)	\(\text{Sens}+\text{Esp}-1\)	0,195	0,358	0,480
Correlación de Matthews (MCC)	\(\dfrac{VP\cdot VN-FP\cdot FN}{\sqrt{(VP+FP)(VP+FN)(VN+FP)(VN+FN)}}\)	0,337	0,325	0,432
Errores y razones de verosimilitud
Tasa de falsos positivos	\(FP/(FP+VN)\)	0,016	0,149	0,134
Tasa de falsos negativos	\(FN/(FN+VP)\)	0,789	0,493	0,386
Razón de verosimilitud positiva \(LR^+\)	\(\text{Sens}/(1-\text{Esp})\)	13,03	3,41	4,57
Razón de verosimilitud negativa \(LR^-\)	\((1-\text{Sens})/\text{Esp}\)	0,802	0,579	0,445
Prevalencia y sesgo del clasificador
Prevalencia	\((VP+FN)/n\)	0,161	0,161	0,161
Tasa de detección	\(VP/n\)	0,034	0,082	0,099
Proporción señalada	\((VP+FP)/n\)	0,048	0,206	0,212
Valor p de McNemar	\(H_0:\ P(FP)=P(FN)\)	\(4,88\times10^{-10}\)	0,045	\(1,44\times10^{-4}\)

La exactitud no discrimina entre modelos con este desbalance. Con el corte de 0,5 la exactitud es 85,9 %, pero la tasa de no información es 83,9 %: la diferencia no es significativa (\(p=0,135\)). Con el corte de Youden la exactitud baja a 79,6 %, por debajo de la NIR, aunque el clasificador es más útil: detecta más del doble de rotaciones. La exactitud premia acertar en la clase mayoritaria, que es la que menos interesa a la empresa, y por eso no se usa para elegir el corte.

Métricas por clase. El corte de 0,5 tiene alta precisión (71,4 % de los señalados rotan) pero una sensibilidad de solo 21,1 %. El de Youden invierte el balance: sensibilidad de 50,7 % (IC 95 %: 39,3 %–62,0 %) y precisión de 39,6 %, que sigue siendo 2,5 veces la prevalencia. El VPN es alto con ambos cortes (86,7 % y 90,0 %): un empleado no señalado tiene baja probabilidad de rotar.

Medidas resumen. El \(F_1\) sube de 0,326 a 0,444, la exactitud balanceada de 0,598 a 0,679 y el kappa de 0,273 a 0,322. El MCC es casi igual con los dos cortes (0,337 y 0,325): resume la correlación entre la clase real y la predicha, y en este caso los dos cortes intercambian errores sin cambiar mucho esa correlación. La elección entre ellos depende, entonces, del costo relativo de cada error, que se analiza en la Sección 9.4.

Sesgo del clasificador. La prueba de McNemar rechaza la simetría de errores con el corte de 0,5 (\(p=4,88\times10^{-10}\)): hay 56 falsos negativos frente a 6 falsos positivos, así que ese corte subestima sistemáticamente la rotación. Con el corte de Youden el desequilibrio se invierte (55 FP frente a 35 FN, \(p=0,045\)): el modelo señala a más empleados de los que rotan, que es el sesgo buscado cuando perder a un empleado cuesta más que intervenirlo.

Entrenamiento frente a prueba. Con el mismo corte, las métricas de entrenamiento son algo mejores (sensibilidad 61,4 % frente a 50,7 %; \(F_1\) 0,531 frente a 0,444), en línea con la brecha de AUC de la Sección 8, que el bootstrap atribuye sobre todo al azar de la partición.

ll    <- function(y, p) -mean(y * log(p) + (1 - y) * log(1 - p))
ks    <- function(r) max(r$sensitivities + r$specificities - 1)
p_auc <- function(y, p) wilcox.test(p[y == 1], p[y == 0], alternative = "greater")$p.value
prob_tab <- data.frame(
  Métrica = c("AUC [IC 95 % de DeLong]", "Coeficiente de Gini", "Estadístico KS",
              "Puntaje de Brier", "Pérdida logarítmica (*log-loss*)", "Valor p: $H_0$ AUC = 0,5"),
  Definición = c("$P(\\hat\\pi_{(1)}>\\hat\\pi_{(0)})$", "$2\\,\\text{AUC}-1$",
                 "$\\max_c\\,[\\text{Sens}(c)+\\text{Esp}(c)-1]$",
                 "$\\frac1n\\sum(\\hat\\pi_i-y_i)^2$",
                 "$-\\frac1n\\sum[y_i\\ln\\hat\\pi_i+(1-y_i)\\ln(1-\\hat\\pi_i)]$",
                 "Mann-Whitney unilateral sobre $\\hat\\pi$"),
  Entrenamiento = c(paste0(num(auc(roc_tr), 3), " [", num(ic_tr[1], 3), "; ", num(ic_tr[3], 3), "]"),
                    num(2 * auc(roc_tr) - 1, 3), num(ks(roc_tr), 3), num(brier(entren$y, p_tr), 4),
                    num(ll(entren$y, p_tr), 4), pval(p_auc(entren$y, p_tr))),
  Prueba = c(paste0(num(auc(roc_te), 3), " [", num(ic_te[1], 3), "; ", num(ic_te[3], 3), "]"),
             num(2 * auc(roc_te) - 1, 3), num(ks(roc_te), 3), num(br_te, 4),
             num(ll(prueba$y, p_te), 4), pval(p_auc(prueba$y, p_te))),
  check.names = FALSE)
prob_tab$Métrica <- gsub("\\*", "", prob_tab$Métrica)
tabla(prob_tab, "Métricas que no dependen del corte, en entrenamiento y en prueba.",
      align = c("l", "l", "r", "r"))

Tabla 8.6: Métricas que no dependen del corte, en entrenamiento y en prueba.

Métrica	Definición	Entrenamiento	Prueba
AUC [IC 95 % de DeLong]	\(P(\hat\pi_{(1)}\gt \hat\pi_{(0)})\)	0,792 [0,751; 0,833]	0,740 [0,671; 0,810]
Coeficiente de Gini	\(2\,\text{AUC}-1\)	0,584	0,480
Estadístico KS	\(\max_c\,[\text{Sens}(c)+\text{Esp}(c)-1]\)	0,485	0,396
Puntaje de Brier	\(\frac1n\sum(\hat\pi_i-y_i)^2\)	0,1072	0,1142
Pérdida logarítmica (log-loss)	\(-\frac1n\sum[y_i\ln\hat\pi_i+(1-y_i)\ln(1-\hat\pi_i)]\)	0,3580	0,3854
Valor p: \(H_0\) AUC = 0,5	Mann-Whitney unilateral sobre \(\hat\pi\)	\(\lt 2\times10^{-16}\)	\(7,03\times10^{-11}\)

Las métricas que no dependen del corte completan la evaluación (Tabla 8.6). El AUC rechaza con holgura la hipótesis de discriminación aleatoria (\(p=7,03\times10^{-11}\) en prueba). El coeficiente de Gini, \(2\,\text{AUC}-1=0,480\), es la misma información en escala \([0,\,1]\), y el estadístico de Kolmogorov-Smirnov (0,396) es la máxima separación entre las distribuciones acumuladas de \(\hat\pi\) en los dos grupos, que coincide con el máximo del índice de Youden sobre la curva ROC. El puntaje de Brier y la pérdida logarítmica miden el error de las probabilidades; ambas penalizan más las predicciones seguras y equivocadas, y la pérdida logarítmica es la devianza media del modelo.

9 Predicción y sensibilidad del punto de corte

9.1 Empleado hipotético

tabla(data.frame(
  Covariable = c("Edad", "Ingreso_Mensual", "Antigüedad_Cargo", "Horas_Extra", "Estado_Civil",
                 "Viaje_Negocios", "Probabilidad estimada de rotar", "IC 95 %"),
  Valor = c(ent(PERFIL$Edad), ent(PERFIL$Ingreso_Mensual), ent(PERFIL$Antigüedad_Cargo),
            as.character(PERFIL$Horas_Extra), as.character(PERFIL$Estado_Civil),
            as.character(PERFIL$Viaje_Negocios), num(pi_hat, 4),
            paste0(num(pi_ic[1], 4), " – ", num(pi_ic[2], 4))),
  `Posición en los datos` = c(
    paste0("Percentil ", num(100 * mean(datos$Edad <= PERFIL$Edad), 0)),
    paste0("Percentil ", num(100 * mean(datos$Ingreso_Mensual <= PERFIL$Ingreso_Mensual), 0)),
    paste0("Percentil ", num(100 * mean(datos$Antigüedad_Cargo <= PERFIL$Antigüedad_Cargo), 0)),
    paste0(pct(mean(datos$Horas_Extra == "No")), " de la base"),
    paste0(pct(mean(datos$Estado_Civil == "Soltero")), " de la base"),
    paste0(pct(mean(datos$Viaje_Negocios == "Raramente")), " de la base"), "", ""),
  check.names = FALSE),
  "Perfil del empleado hipotético y probabilidad estimada de rotación.")

Tabla 9.1: Perfil del empleado hipotético y probabilidad estimada de rotación.

Covariable	Valor	Posición en los datos
Edad	32	Percentil 35
Ingreso_Mensual	2.500	Percentil 15
Antigüedad_Cargo	1	Percentil 20
Horas_Extra	No	71,7 % de la base
Estado_Civil	Soltero	32,0 % de la base
Viaje_Negocios	Raramente	71,0 % de la base
Probabilidad estimada de rotar	0,2541
IC 95 %	0,2003 – 0,3165

El perfil corresponde a un empleado joven que empieza en su cargo, con un ingreso en la parte baja de la distribución, soltero, sin horas extra y que viaja raramente (Tabla 9.1). Todos los valores están dentro del rango observado, así que la predicción es una interpolación. Se calcula con el modelo de la base completa. El IC se construye en la escala del predictor lineal, \(\hat\eta\pm z_{0,975}\,\widehat{\mathrm{EE}}(\hat\eta)\) con \(\widehat{\mathrm{EE}}(\hat\eta)=\sqrt{\mathbf x_0^\top\widehat{\mathrm{Var}}(\hat{\boldsymbol\beta})\mathbf x_0}\), y se transforma con \(\pi=1/(1+e^{-\eta})\); así el intervalo queda dentro de \([0,\,1]\), lo que no garantiza el método delta en la escala de \(\pi\). Resulta \(\hat\pi_0=0,2541\), con IC 95 % \([0,2003;\ 0,3165]\): 1,58 veces la prevalencia de la base.

9.2 Métricas sobre la rejilla de cortes

Se evalúan los cortes \(c\in\{0{,}05;0{,}06;\dots;0{,}95\}\). Para cada uno se calculan la sensibilidad, la especificidad, el valor predictivo positivo (VPP: proporción de rotantes entre los señalados), el negativo (VPN), la exactitud y \(F_1=2\,\mathrm{VPP}\cdot\mathrm{Sens}/(\mathrm{VPP}+\mathrm{Sens})\), en entrenamiento y en prueba.

R_COSTOS <- c(2, 5, 10)
c_cost <- 1 / R_COSTOS          # c* = C_I / C_L (Elkan, 2001)
rej_te %>% select(corte, Sens, Esp, VPP, F1) %>%
  pivot_longer(-corte, names_to = "Métrica", values_to = "v") %>%
  ggplot(aes(corte, v, colour = Métrica)) +
  geom_vline(xintercept = c_cost, colour = "grey80", linetype = "dashed") +
  geom_vline(xintercept = c_youden, colour = COL_NAR, linewidth = 0.8) +
  geom_line(linewidth = 0.8) +
  geom_text(data = . %>% group_by(Métrica) %>% filter(abs(corte - 0.60) < 1e-9),
            aes(label = Métrica), vjust = -0.6, size = 3.2, show.legend = FALSE) +
  scale_colour_manual(values = c(Sens = COL_AZUL, Esp = COL_JAV, VPP = COL_TEAL, F1 = COL_NAR)) +
  scale_x_continuous(breaks = seq(0, 1, 0.1)) +
  labs(title = "Subir el corte gana especificidad a costa de dejar escapar rotaciones",
       x = "Punto de corte c", y = NULL)

Figura 9.1: Métricas de clasificación en la muestra de prueba según el punto de corte. Las líneas verticales marcan el corte de Youden (naranja) y los cortes por costos (gris).

CORTES_SEL <- sort(unique(round(c(0.10, 0.15, 0.20, 0.30, 0.40, 0.50, c_youden, c_cost), 3)))
sel <- bind_rows(lapply(CORTES_SEL, function(c)
  bind_cols(metricas(entren$y, p_tr, c) %>% select(corte, Sens, Esp) %>%
              rename(`Sens (ent.)` = Sens, `Esp (ent.)` = Esp),
            metricas(prueba$y, p_te, c) %>% select(-corte, -Youden))))
tabla(sel %>% rename(Corte = corte), "Métricas de clasificación en cortes seleccionados, en entrenamiento (ent.) y en prueba.",
      digits = c(3, 3, 3, 3, 3, 3, 3, 3, 3), font = 8)

Tabla 9.2: Métricas de clasificación en cortes seleccionados, en entrenamiento (ent.) y en prueba.

Corte	Sens (ent.)	Esp (ent.)	Sens	Esp	VPP	VPN	Exactitud	F1
0,100	0,831	0,562	0,761	0,532	0,238	0,921	0,569	0,362
0,150	0,753	0,710	0,662	0,695	0,294	0,915	0,689	0,407
0,200	0,669	0,793	0,549	0,784	0,328	0,901	0,746	0,411
0,250	0,614	0,866	0,507	0,851	0,396	0,900	0,796	0,444
0,300	0,506	0,908	0,451	0,903	0,471	0,895	0,830	0,460
0,400	0,325	0,951	0,310	0,962	0,611	0,879	0,857	0,411
0,500	0,217	0,978	0,211	0,984	0,714	0,867	0,859	0,326

m05 <- metricas(prueba$y, p_te, 0.5)

La Tabla 9.2 y la Figura 9.1 muestran el intercambio entre sensibilidad y especificidad. Con el corte convencional de 0,5 el modelo detecta solo el 21,1 % de las rotaciones de prueba, aunque su exactitud (85,9 %) apenas supera el 83,9 % que se obtendría clasificando a todos como “No”. Es la consecuencia del desbalance anticipada en la Sección 5. El corte de 0,5 no es arbitrario —equivale a una razón de costos \(R=2\), como se ve en la Sección 9.4—, pero con esa exactitud alta y esa sensibilidad baja resulta poco útil si perder a un empleado cuesta bastante más que retenerlo.

9.3 Criterio de Youden

El índice \(J(c)=\text{Sens}(c)+\text{Esp}(c)-1\) es la distancia vertical entre la curva ROC y la diagonal; su máximo da el corte que mejor separa los grupos si los dos errores pesan igual (Youden, 1950). El corte se elige en entrenamiento y se evalúa en prueba, para que las métricas reportadas no estén sesgadas a favor del corte.

El corte de Youden en entrenamiento es \(c_J=0,25\) (\(J=0,480\)). En prueba detecta el 50,7 % de las rotaciones con una especificidad del 85,1 % (Tabla 8.4; métricas completas en la Tabla 8.5). El VPP es 39,6 %: de cada 100 empleados señalados, unos 40 efectivamente rotan, frente a 16 si se eligieran al azar.

9.4 Criterio de costos

Sean \(C_I\) el costo de una acción de retención y \(C_L\) el costo de perder a un empleado. Se supone que la intervención evita la rotación (\(e=1\)) y que su costo se paga siempre, rote o no el empleado. Si \(\hat\pi\) está calibrada, intervenir tiene costo esperado \(C_I\) y no intervenir, \(\hat\pi\,C_L\). Conviene intervenir cuando \[ \hat\pi\,C_L\ \ge\ C_I \iff \hat\pi\ \ge\ c^*=\frac{C_I}{C_L}=\frac{1}{R}, \qquad R=\frac{C_L}{C_I} . \] Con eficacia parcial \(e<1\), el corte es \(c^*=1/(eR)\) (Elkan, 2001). El supuesto es que perder a un empleado cuesta más que intervenirlo (\(R>1\)). Como \(R\) no se conoce, se evalúa \(R\in\{2,\,5,\,10\}\), que da \(c^*=0,500;\ 0,200;\ 0,100\). Con \(R=2\) el corte coincide con el convencional de 0,5: usar 0,5 equivale a suponer que perder a un empleado cuesta solo el doble que intervenirlo. La calibración verificada en la Sección 8 respalda el uso de este criterio.

costo_tab <- bind_rows(lapply(seq_along(R_COSTOS), function(i) {
  m <- metricas(prueba$y, p_te, c_cost[i])
  data.frame(R = R_COSTOS[i], Corte = c_cost[i], Sens = m$Sens, Esp = m$Esp, VPP = m$VPP,
             `% de empleados intervenidos` = 100 * mean(p_te >= c_cost[i]), check.names = FALSE)
}))
tabla(costo_tab, "Cortes por costos y su desempeño en la muestra de prueba.",
      digits = c(0, 3, 3, 3, 3, 1))

Tabla 9.3: Cortes por costos y su desempeño en la muestra de prueba.

R	Corte	Sens	Esp	VPP	% de empleados intervenidos
2	0,500	0,211	0,984	0,714	4,8
5	0,200	0,549	0,784	0,328	27,0
10	0,100	0,761	0,532	0,238	51,5

A mayor \(R\), más bajo el corte: se interviene a más empleados para no perder rotaciones. Con \(R=10\) se detecta el 76,1 % de las rotaciones, pero se interviene al 51,5 % de la plantilla (Tabla 9.3).

9.5 Decisión sobre el empleado hipotético

CRIT <- data.frame(Criterio = c("Convencional", "Youden", paste0("Costos, R = ", R_COSTOS)),
                   Corte = c(0.5, c_youden, c_cost))
pi_hat_tr <- unname(predict(modelo_tr, PERFIL, type = "response"))
dec <- CRIT %>% mutate(
  `pi estimada` = pi_hat,
  Decisión = ifelse(pi_hat >= Corte, "Intervenir", "No intervenir"),
  `pi entrenamiento` = pi_hat_tr,
  `Decisión (entren.)` = ifelse(pi_hat_tr >= Corte, "Intervenir", "No intervenir"),
  `Estable en c ± 0,05` = ifelse(abs(pi_hat - Corte) > 0.05 & abs(pi_hat_tr - Corte) > 0.05,
                                 "Sí", "No"),
  `IC contiene c` = ifelse(pi_ic[1] <= Corte & Corte <= pi_ic[2], "Sí (incierta)", "No"))
names(dec)[c(3, 5)] <- c("$\\hat\\pi_0$ (completa)", "$\\hat\\pi_0$ (entren.)")
tabla(dec, "Decisión de intervenir al empleado hipotético según cada criterio de corte, con los dos modelos.",
      digits = c(0, 3, 4, 0, 4, 0, 0, 0), font = 8)

Tabla 9.4: Decisión de intervenir al empleado hipotético según cada criterio de corte, con los dos modelos.

Criterio	Corte	\(\hat\pi_0\) (completa)	Decisión	\(\hat\pi_0\) (entren.)	Decisión (entren.)	Estable en c ± 0,05	IC contiene c
Convencional	0,500	0,2541	No intervenir	0,2425	No intervenir	Sí	No
Youden	0,250	0,2541	Intervenir	0,2425	No intervenir	No	Sí (incierta)
Costos, R = 2	0,500	0,2541	No intervenir	0,2425	No intervenir	Sí	No
Costos, R = 5	0,200	0,2541	Intervenir	0,2425	Intervenir	No	No
Costos, R = 10	0,100	0,2541	Intervenir	0,2425	Intervenir	Sí	No

n_interv <- sum(pi_hat >= CRIT$Corte)

La decisión depende del criterio (Tabla 9.4). Una decisión es estable si \(|\hat\pi_0-c|>0{,}05\) con los dos modelos (base completa y entrenamiento), es decir, si no cambia al mover el corte en \(\pm0{,}05\) con ninguno de ellos; es incierta si el IC de \(\hat\pi_0\) contiene al corte, porque entonces los datos no permiten afirmar de qué lado está la probabilidad verdadera.

El corte de Youden se elige con el modelo de entrenamiento, mientras que \(\hat\pi_0\) se calcula con el de la base completa. Con el modelo de entrenamiento, \(\hat\pi_0=0,2425\): las decisiones coinciden en los cortes por costos, pero no en el de Youden, que el empleado supera con un modelo y no con el otro. La decisión con ese criterio depende, entonces, de qué modelo se use, y no solo del corte; la recomendación se apoya en el criterio de costos.

Con \(\hat\pi_0=0,254\), el empleado supera el corte de costos con \(R=5\), aunque el límite inferior de su IC (0,2003) queda casi sobre el corte (0,20); no alcanza el de \(R=2\), que es el corte convencional de 0,5. Intervenir conviene cuando \(R\ge 1/\hat\pi_0\), es decir, cuando perder al empleado cuesta al menos 3,9 veces lo que cuesta la intervención (4,1 con el modelo de entrenamiento). Se recomienda intervenirlo si la empresa acepta esa razón de costos: su riesgo es 1,58 veces el promedio, y lo explican variables modificables (un ingreso bajo y poco tiempo en el cargo) y una no modificable (ser soltero). Si la empresa estima que la rotación cuesta menos de esa razón —por ejemplo, el doble que la intervención—, la recomendación se invierte. La acción sugerida es la que corresponde a sus factores de riesgo: revisar su trayectoria salarial y darle un plan de desarrollo en el cargo.

10 Estrategia de retención y conclusiones

10.1 Balance de las hipótesis

p_mod <- av$`Pr(>Chisq)`
s_mod <- sign(c(b["Horas_ExtraSi"], b["Estado_CivilSoltero"], b["Viaje_NegociosFrecuentemente"],
                b["Edad"], b["log(Ingreso_Mensual)"], b["Antigüedad_Cargo"]))
balance <- data.frame(Variable = HIP$Variable, `Dirección esperada` = HIP$`Dirección esperada`,
                      `Bivariado` = contraste$Veredicto,
                      `Modelo ajustado` = veredicto(p_mod, s_mod, HIP$Signo),
                      `p (LR, modelo)` = pval(p_mod),
                      `Rango de importancia` = match(c("Horas_Extra", "Estado_Civil", "Viaje_Negocios",
                                                      "Edad", "ln(Ingreso_Mensual)", "Antigüedad_Cargo"),
                                                    av_tab$Término),
                      check.names = FALSE)
tabla(balance, "Balance de las hipótesis en el análisis bivariado y en el modelo ajustado.", font = 8.5)

Tabla 10.1: Balance de las hipótesis en el análisis bivariado y en el modelo ajustado.

Variable	Dirección esperada	Bivariado	Modelo ajustado	p (LR, modelo)	Rango de importancia
Horas_Extra	Sí > No	Confirmada	Confirmada	\(\lt 2\times10^{-16}\)	1
Estado_Civil	Soltero > Casado	Confirmada	Confirmada	\(8,45\times10^{-8}\)	2
Viaje_Negocios	Creciente con la frecuencia	Confirmada	Confirmada	\(4,25\times10^{-5}\)	3
Edad	Negativa	Confirmada	Confirmada	0,015	6
Ingreso_Mensual	Negativa	Confirmada	Confirmada	\(7,44\times10^{-5}\)	4
Antigüedad_Cargo	Negativa	Confirmada	Confirmada	0,001	5

n_ambos <- sum(balance$Bivariado == "Confirmada" & balance$`Modelo ajustado` == "Confirmada")

Las 6 hipótesis se confirman en ambos análisis (Tabla 10.1): ninguna resultó significativa en sentido contrario ni dejó de ser significativa al ajustar. La única precisión es que el efecto del estado civil corresponde a los solteros; los divorciados no difieren de los casados.

10.2 Estrategia

Las acciones siguen el orden de importancia del modelo. Solo se proponen sobre variables significativas, y se distingue entre las que la empresa puede modificar y las que solo sirven para focalizar.

1. Horas extra (modificable; factor de mayor peso). Redistribuir la carga: contratar personal o crear turnos de refuerzo en las áreas con horas extra recurrentes, poner un tope mensual por empleado y compensarlas con tiempo libre. Es el factor con mayor asociación: los odds de quienes hacen horas extra son 4,34 veces los de quienes no las hacen. Que reducirlas disminuya la rotación es una hipótesis de intervención, no un resultado del modelo.
2. Estado civil (no modificable; focaliza). Los solteros rotan más. No es una variable sobre la que se actúe, pero orienta el diseño: planes de carrera visibles y beneficios valorados por empleados sin cargas familiares (formación, movilidad interna), que compitan con las ofertas externas.
3. Viajes de negocios (modificable). Reducir la frecuencia de los viajes con reuniones virtuales cuando sea posible, rotar los viajes entre más personas y compensar los días fuera de casa.
4. Ingreso (modificable). Revisar la política salarial en la parte baja de la escala. En el modelo principal el efecto es relativo: un mismo aumento porcentual multiplica los odds por el mismo factor en cualquier nivel de ingreso. El modelo con splines matiza esa lectura, porque la curva del ingreso es más empinada en los valores bajos (Anexo D). Como además un mismo aumento porcentual cuesta menos en términos absolutos en los salarios bajos, ahí se concentra la mejor relación entre costo y efecto esperado.
5. Antigüedad en el cargo (modificable de forma indirecta). Reforzar el acompañamiento en los primeros años en el cargo: inducción, mentoría y metas de desarrollo, para acelerar el ajuste al puesto.
6. Edad (no modificable; focaliza). Priorizar a los empleados jóvenes en los programas anteriores. El modelo con splines indica que los empleados de mayor edad también tienen riesgo elevado; en ellos la rotación puede corresponder a retiro.

El modelo permite, además, priorizar la intervención: calcular \(\hat\pi\) para cada empleado y aplicar el corte que corresponda a la razón de costos que estime la empresa (Sección 9).

10.3 Conclusiones

1. La rotación afecta al 16,1 % de los empleados. Con este desbalance, la exactitud no sirve para evaluar el modelo, y el corte de 0,5 solo es adecuado si perder a un empleado cuesta apenas el doble de intervenirlo (\(R=2\)).
2. Las seis covariables están asociadas con la rotación, en la dirección esperada, tanto en el análisis bivariado como en el modelo ajustado. Ningún coeficiente cambia de signo al ajustar.
3. Las horas extra son el factor de mayor peso, seguidas del estado civil, los viajes de negocios y el ingreso; la antigüedad en el cargo y la edad tienen efectos protectores menores.
4. El modelo discrimina de forma aceptable en datos no usados para estimarlo (AUC = 0,740; Gini = 0,480), con sobreajuste escaso (AUC corregido por optimismo = 0,770), y sus probabilidades están razonablemente calibradas. Con el corte de Youden, en prueba alcanza una sensibilidad de 50,7 %, una especificidad de 85,1 %, una precisión de 39,6 % y un \(F_1\) de 0,444; la exactitud no es informativa con este desbalance, porque no supera la de la regla trivial.
5. La decisión de intervenir depende de la razón de costos \(R\). Para el empleado hipotético (\(\hat\pi_0=0,254\)) conviene intervenir cuando \(R\ge 1/\hat\pi_0=3,9\), es decir, cuando perderlo cuesta al menos 3,9 veces lo que cuesta la intervención (4,1 con el modelo de entrenamiento).

11 Limitaciones

• Asociación, no causalidad. Los datos son observacionales: los OR describen asociaciones condicionales. Por ejemplo, las horas extra podrían ser consecuencia de otros problemas del puesto y no su causa.
• Forma funcional. Box-Tidwell y Hosmer-Lemeshow indican que la linealidad en el logit es una aproximación. En el ingreso y la antigüedad el efecto se concentra en los valores bajos; en la edad la relación tiene forma de U, y el OR lineal no representa a los empleados de mayor edad.
• Qué mide Rotación. La variable no distingue entre un cambio interno de cargo y un retiro de la empresa, que podrían tener determinantes distintos.
• Carácter de la base. Los datos son una adaptación didáctica de Weiers (2006); por su tamaño (1.470 registros), su número de rotaciones (237) y sus primeros registros, coinciden con el conjunto IBM HR Analytics Employee Attrition, que es ficticio. Las conclusiones ilustran el método y no deben trasladarse sin más a una organización concreta.
• Costos. Los cortes por costos dependen de una razón \(R\) que la empresa debe estimar; aquí solo se evaluó su sensibilidad.
• Validación. La partición y el bootstrap validan internamente el modelo; su desempeño en otra organización o en otro periodo no está garantizado.

Referencias

Agresti, A. (2013). Categorical Data Analysis (3.ª ed.). Wiley.

Anderson, T. W. y Darling, D. A. (1954). A test of goodness of fit. Journal of the American Statistical Association, 49(268), 765–769.

Armitage, P. (1955). Tests for linear trends in proportions and frequencies. Biometrics, 11(3), 375–386.

Box, G. E. P. y Tidwell, P. W. (1962). Transformation of the independent variables. Technometrics, 4(4), 531–550.

Clopper, C. J. y Pearson, E. S. (1934). The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika, 26(4), 404–413.

Cohen, J. (1960). A coefficient of agreement for nominal scales. Educational and Psychological Measurement, 20(1), 37–46.

DeLong, E. R., DeLong, D. M. y Clarke-Pearson, D. L. (1988). Comparing the areas under two or more correlated receiver operating characteristic curves: A nonparametric approach. Biometrics, 44(3), 837–845.

Elkan, C. (2001). The foundations of cost-sensitive learning. Proceedings of the Seventeenth International Joint Conference on Artificial Intelligence (IJCAI), 973–978.

Fawcett, T. (2006). An introduction to ROC analysis. Pattern Recognition Letters, 27(8), 861–874.

Fox, J. y Monette, G. (1992). Generalized collinearity diagnostics. Journal of the American Statistical Association, 87(417), 178–183.

Greenland, S., Robins, J. M. y Pearl, J. (1999). Confounding and collapsibility in causal inference. Statistical Science, 14(1), 29–46.

Harrell, F. E. (2015). Regression Modeling Strategies (2.ª ed.). Springer.

Hauck, W. W. y Donner, A. (1977). Wald’s test as applied to hypotheses in logit analysis. Journal of the American Statistical Association, 72(360), 851–853.

Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6(2), 65–70.

Hosmer, D. W., Lemeshow, S. y Sturdivant, R. X. (2013). Applied Logistic Regression (3.ª ed.). Wiley.

Joanes, D. N. y Gill, C. A. (1998). Comparing measures of sample skewness and kurtosis. Journal of the Royal Statistical Society, Series D, 47(1), 183–189.

Mann, H. B. y Whitney, D. R. (1947). On a test of whether one of two random variables is stochastically larger than the other. The Annals of Mathematical Statistics, 18(1), 50–60.

Matthews, B. W. (1975). Comparison of the predicted and observed secondary structure of T4 phage lysozyme. Biochimica et Biophysica Acta, 405(2), 442–451.

McFadden, D. (1979). Quantitative methods for analysing travel behaviour of individuals: Some recent developments. En D. A. Hensher y P. R. Stopher (Eds.), Behavioural Travel Modelling (pp. 279–318). Croom Helm.

McNemar, Q. (1947). Note on the sampling error of the difference between correlated proportions or percentages. Psychometrika, 12(2), 153–157.

Peduzzi, P., Concato, J., Kemper, E., Holford, T. R. y Feinstein, A. R. (1996). A simulation study of the number of events per variable in logistic regression analysis. Journal of Clinical Epidemiology, 49(12), 1373–1379.

Shapiro, S. S. y Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611.

Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.

Weiers, R. M. (2006). Introducción a la estadística para negocios (5.ª ed.). Thomson.

Wilson, E. B. (1927). Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association, 22(158), 209–212.

Youden, W. J. (1950). Index for rating diagnostic tests. Cancer, 3(1), 32–35.

Anexos

Anexo A. Instalación del entorno

install.packages("devtools")                       # solo la primera vez
devtools::install_github("centromagis/paqueteMODELOS", force = TRUE)
install.packages(c("dplyr", "tidyr", "ggplot2", "kableExtra", "car", "pROC",
                   "rms", "e1071", "nortest", "gridExtra", "bookdown", "rmarkdown"))

El código de cada tabla y de cada figura se despliega con el botón Code que aparece a su derecha; el menú Code del encabezado los muestra u oculta todos a la vez. La compilación del archivo .Rmd reproduce todas las cifras con la semilla 2026.

Anexo B. Bondad de ajuste por deciles

tabla(hl_grp, "Rotaciones observadas y esperadas por decil de probabilidad estimada.",
      digits = c(0, 0, 0, 1, 1))

Tabla 11.1: Rotaciones observadas y esperadas por decil de probabilidad estimada.

Decil	n	Observados	Esperados	Observados - esperados
1	147	7	2,8	4,2
2	147	9	5,4	3,6
3	147	11	8,0	3,0
4	147	9	10,7	-1,7
5	147	7	13,9	-6,9
6	147	17	18,0	-1,0
7	147	19	24,0	-5,0
8	147	24	32,7	-8,7
9	147	51	45,2	5,8
10	147	83	76,3	6,7

La Tabla 11.1 desagrega la prueba de Hosmer-Lemeshow con \(g=10\). En los tres primeros deciles y en los dos últimos se observan más rotaciones de las esperadas; en los deciles intermedios, menos.

Anexo C. Observaciones influyentes

tabla(infl, "Diagnóstico de observaciones influyentes.", digits = c(0, 4, 0))

Tabla 11.2: Diagnóstico de observaciones influyentes.

Medida	Umbral	Observaciones
Distancia de Cook \(\gt 4/n\)	0,0027	115
Distancia de Cook \(\gt 1\)	1,0000	0
Apalancamiento \(\gt 2p/n\)	0,0122	125
Algún \(|\mathrm{DFBETAS}| \gt 2/\sqrt{n}\)	0,0522	364

tabla(dfb_tab, "Máximo valor absoluto de los DFBETAS por coeficiente.", digits = c(0, 3, 0))

Tabla 11.3: Máximo valor absoluto de los DFBETAS por coeficiente.

Término	Máximo |DFBETAS|	Observaciones con |DFBETAS| > 2/raíz(n)
Intercepto	0,147	121
Horas extra: Sí vs. No	0,097	202
Estado civil: Divorciado vs. Casado	0,140	98
Estado civil: Soltero vs. Casado	0,090	197
Viaje: Raramente vs. No viaja	0,268	35
Viaje: Frecuentemente vs. No viaja	0,247	74
Edad (por año)	0,235	132
ln(Ingreso mensual)	0,150	123
Antigüedad en el cargo (por año)	0,222	104

Con el umbral \(4/n\) se marcan 115 observaciones (Tabla 11.2); el 95,7 % son rotaciones, con una mediana de \(\hat\pi\) de 0,149: empleados que rotaron con un perfil de bajo riesgo, no errores de registro. Las 125 observaciones con apalancamiento mayor que \(2p/n\) corresponden a combinaciones poco frecuentes de las covariables; su apalancamiento máximo es 0,0341, y ninguna de ellas tiene distancia de Cook cercana a 1. El DFBETAS mide el cambio en \(\hat\beta_j\) al excluir la observación, en unidades de su error estándar; el máximo por coeficiente está en la Tabla 11.3.

Nota

Excluir en bloque las observaciones con Cook \(>4/n\) no es un análisis de sensibilidad válido aquí: casi todas son rotaciones, así que eliminarlas vaciaría categorías con pocos eventos (en No_Viaja solo rotaron 12 empleados) y produciría separación cuasicompleta, con OR artificialmente grandes.

Anexo D. Sensibilidad a la forma funcional

tabla(data.frame(Modelo = c("Lineal (principal)", "Splines naturales"),
                 Parámetros = c(length(coef(modelo)), length(coef(modelo_ns))),
                 AIC = c(AIC(modelo), AIC(modelo_ns)),
                 `AUC aparente` = c(auc_lin, auc_ns), check.names = FALSE),
      "Comparación del modelo lineal con el modelo con splines.", digits = c(0, 0, 1, 4))

Tabla 11.4: Comparación del modelo lineal con el modelo con splines.

Modelo	Parámetros	AIC	AUC aparente
Lineal (principal)	9	1.092,6	0,7776
Splines naturales	15	1.079,7	0,7847

ggplot(curvas, aes(x, eta, colour = Modelo)) + geom_line(linewidth = 0.8) +
  scale_colour_manual(values = c(Lineal = COL_VIO, Splines = COL_TEAL)) +
  facet_wrap(~ var, scales = "free_x", nrow = 1) +
  labs(title = "Misma dirección; con splines, el riesgo se concentra en valores bajos",
       x = NULL, y = "Logit estimado")

Figura 11.1: Efecto parcial de cada covariable cuantitativa sobre el logit, en el modelo lineal (violeta) y en el modelo con splines (verde azulado). Las demás covariables se fijan en su mediana o en su nivel de referencia.

El modelo alternativo reemplaza cada cuantitativa por un spline cúbico natural con 3 grados de libertad (splines::ns). Mejora el ajuste (\(G^2=24,91\), 6 gl, \(p=3,55\times10^{-4}\)), pero la ganancia en discriminación es pequeña (Tabla 11.4). Las curvas del ingreso y de la antigüedad en el cargo son decrecientes; la de la edad tiene forma de U, con mínimo cerca de los 43 años y un aumento del riesgo en las edades mayores, compatible con retiros (Figura 11.1). La dirección lineal negativa de la edad describe la tendencia promedio, no un efecto monótono. Con los splines, la prueba de Hosmer-Lemeshow da \(p=0,101\), \(p=0,005\) y \(p=0,148\) para \(g=8,\,10,\,12\): la no linealidad explica parte de la falta de ajuste del modelo lineal, aunque con \(g=10\) la prueba sigue rechazando.

Anexo E. Información de la sesión

si <- sessionInfo()
paq <- c("paqueteMODELOS", "dplyr", "tidyr", "ggplot2", "kableExtra", "car", "pROC",
         "rms", "e1071", "nortest")
tabla(data.frame(
  Elemento = c("Versión de R", "Plataforma", "Semilla global",
               paste("Paquete", paq)),
  Valor = c(si$R.version$version.string, si$R.version$platform, "2026",
            sapply(paq, function(p) as.character(packageVersion(p)))),
  row.names = NULL),
  "Entorno de cómputo utilizado.")

Tabla 11.5: Entorno de cómputo utilizado.

Elemento	Valor
Versión de R	R version 4.5.2 (2025-10-31 ucrt)
Plataforma	x86_64-w64-mingw32
Semilla global	2026
Paquete paqueteMODELOS	0.1.0
Paquete dplyr	1.2.1
Paquete tidyr	1.3.2
Paquete ggplot2	4.0.3
Paquete kableExtra	1.4.0
Paquete car	3.1.3
Paquete pROC	1.19.0.1
Paquete rms	8.1.1
Paquete e1071	1.7.17
Paquete nortest	1.0.4

Oriana Giraldo Arcia · Luis Javier Rubio Hernández
Pontificia Universidad Javeriana Cali