# devtools::install_github("centromagis/paqueteMODELOS", force = TRUE)
library(paqueteMODELOS)
library(tidyverse)
library(broom)
library(corrplot)
library(car)
library(pROC)
library(knitr)
library(kableExtra)

theme_set(theme_minimal(base_size = 12))
# Formatos para el texto y las tablas
f2  <- function(x, d = 2) formatC(x, format = "f", digits = d)
pct <- function(x, d = 1) paste0(formatC(100 * x, format = "f", digits = d), "%")
fp  <- function(p) ifelse(p < 0.001, "< 0.001", formatC(p, format = "f", digits = 3))
fpi <- function(p) ifelse(p < 0.001, "p < 0.001", paste0("p = ", formatC(p, format = "f", digits = 3)))
red <- function(or) pct(1 - or, 0)   # reducción porcentual de los odds

tabla <- function(x, full_width = FALSE, ...) {
  kable(x, ...) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = full_width)
}

# Métricas de clasificación para un punto de corte (clase positiva: y = 1, rota)
metricas_corte <- function(prob, y, corte) {
  pred <- as.integer(prob >= corte)
  VP <- sum(pred == 1 & y == 1); FP <- sum(pred == 1 & y == 0)
  FN <- sum(pred == 0 & y == 1); VN <- sum(pred == 0 & y == 0)
  tibble(Corte = corte, VP = VP, FP = FP, FN = FN, VN = VN,
         Sensibilidad  = VP / (VP + FN),
         Especificidad = VN / (VN + FP),
         VPP           = ifelse(VP + FP > 0, VP / (VP + FP), NA_real_),
         Exactitud     = (VP + VN) / length(y),
         Intervenidos  = (VP + FP) / length(y))
}

# Lectura convencional del AUC (Hosmer y Lemeshow)
leer_auc <- function(a) {
  if (a < 0.6) "pobre" else if (a < 0.7) "débil" else if (a < 0.8) "aceptable"
  else if (a < 0.9) "muy buena" else "sobresaliente"
}

1 El caso

Una empresa quiere anticipar qué empleados van a rotar de cargo y saber qué factores pesan más en esa decisión. Tiene datos históricos de 1470 trabajadores. Con ellos estimamos un modelo de regresión logística que da la probabilidad de rotación de cada empleado, medimos qué tan bien predice, fijamos un punto de corte para decidir a quién intervenir y, con las variables que resultan significativas, proponemos cómo reducir la rotación.

2 Datos y preparación

data("rotacion")

n_na  <- sum(is.na(rotacion))
n_dup <- sum(duplicated(rotacion))
c(filas = nrow(rotacion), columnas = ncol(rotacion),
  valores_faltantes = n_na, filas_duplicadas = n_dup)
##             filas          columnas valores_faltantes  filas_duplicadas 
##              1470                24                 0                 0

La base rotacion tiene 1470 registros y 24 variables, sin valores faltantes ni filas duplicadas. La preparamos así: la respuesta queda como \(y=1\) si el empleado rota y \(y=0\) si no; las variables ordinales que vienen como números pasan a factores con sus etiquetas; las categorías de referencia son “no trabaja horas extra”, “casado” y “muy insatisfecho”; y el ingreso se expresa en miles para que su coeficiente se lea por cada 1000.

df_full <- rotacion %>% mutate(
  y = ifelse(Rotación == "Si", 1, 0),

  Rendimiento_Laboral = factor(Rendimiento_Laboral, levels = 1:4,
                               labels = c("Bajo", "Medio", "Alto", "Muy alto")),
  Educación = factor(Educación, levels = 1:5,
                     labels = c("Primaria", "Secundaria", "Técnico/Tecnólogo", "Pregrado", "Posgrado")),
  Satisfacción_Ambiental = factor(Satisfacción_Ambiental, levels = 1:4,
                                  labels = c("Muy insatisfecho", "Insatisfecho", "Satisfecho", "Muy satisfecho")),
  Satisfación_Laboral = factor(Satisfación_Laboral, levels = 1:4,
                               labels = c("Muy insatisfecho", "Insatisfecho", "Satisfecho", "Muy satisfecho")),
  Equilibrio_Trabajo_Vida = factor(Equilibrio_Trabajo_Vida, levels = 1:4,
                                   labels = c("Muy bajo", "Bajo", "Medio", "Alto")),
  across(where(is.character), as.factor),

  Horas_Extra         = relevel(Horas_Extra, ref = "No"),
  Estado_Civil        = relevel(Estado_Civil, ref = "Casado"),
  Satisfación_Laboral = relevel(Satisfación_Laboral, ref = "Muy insatisfecho"),

  Ingreso_Miles = Ingreso_Mensual / 1000
)

stopifnot(sum(is.na(df_full)) == 0)   # la recodificación no genera faltantes

3 Selección de variables e hipótesis

num_vars <- c("Edad", "Distancia_Casa", "Ingreso_Mensual", "Trabajos_Anteriores",
              "Años_Experiencia", "Antigüedad", "Antigüedad_Cargo",
              "Años_ultima_promoción", "Años_acargo_con_mismo_jefe",
              "Porcentaje_aumento_salarial", "Capacitaciones")

mat_cor <- cor(df_full[num_vars])

corrplot(mat_cor, method = "color", type = "upper",
         col = colorRampPalette(c(col_alerta, "white", col_primario))(200),
         tl.cex = 0.7, tl.col = "black", addCoef.col = "black", number.cex = 0.55)

Figura 1. Correlaciones entre las variables cuantitativas de la base.

Varias variables cuantitativas miden casi lo mismo. La experiencia total se mueve con el ingreso (r = 0.77) y la antigüedad en la empresa con la antigüedad en el cargo (r = 0.76). Ponerlas juntas en el modelo inflaría los errores estándar, como vimos con la multicolinealidad en la Unidad 2, así que elegimos cuantitativas poco correlacionadas entre sí.

Variable Tipo Por qué se relacionaría con la rotación Hipótesis
Horas_Extra Cualitativa nominal Las horas extra desgastan y quitan tiempo personal. Quien trabaja horas extra rota más.
Satisfación_Laboral Cualitativa ordinal Quien está a gusto tiene menos motivos para cambiar. A mayor satisfacción, menor rotación.
Estado_Civil Cualitativa nominal Con compromisos familiares se valora más la estabilidad. Los solteros rotan más que los casados.
Edad Cuantitativa Con los años la carrera se asienta. A mayor edad, menor rotación.
Ingreso_Mensual Cuantitativa Un mejor salario encarece la decisión de irse. A mayor ingreso, menor rotación.
Antigüedad_Cargo Cuantitativa El tiempo en el cargo genera arraigo. A mayor antigüedad en el cargo, menor rotación. 
cat_vars   <- c("Horas_Extra", "Satisfación_Laboral", "Estado_Civil")
cuant_vars <- c("Edad", "Ingreso_Miles", "Antigüedad_Cargo")
vars       <- c(cat_vars, cuant_vars)

df    <- df_full %>% select(y, Rotación, all_of(vars))
p_rot <- mean(df$y)

cor_sel <- cor(df[cuant_vars])

La mayor correlación entre las tres cuantitativas elegidas es 0.50, entre edad e ingreso.

4 Análisis univariado

Las variables cualitativas se describen con frecuencias y barras; las cuantitativas, con medidas de centro, dispersión y forma, e histogramas.

Tabla 1. Frecuencias de la respuesta y de las variables cualitativas.

uni_cat <- map_dfr(c("Rotación", cat_vars), function(v) {
  df %>% count(Categoría = as.character(.data[[v]])) %>%
    mutate(Variable = v, Porcentaje = pct(n / sum(n))) %>%
    select(Variable, Categoría, Frecuencia = n, Porcentaje)
})

tabla(uni_cat) %>% collapse_rows(columns = 1, valign = "top")
Variable Categoría Frecuencia Porcentaje
Rotación No 1233 83.9%
Si 237 16.1%
Horas_Extra No 1054 71.7%
Si 416 28.3%
Satisfación_Laboral Insatisfecho 280 19.0%
Muy insatisfecho 289 19.7%
Muy satisfecho 459 31.2%
Satisfecho 442 30.1%
Estado_Civil Casado 673 45.8%
Divorciado 327 22.2%
Soltero 470 32.0%

Rotaron 237 de los 1470 empleados, el 16.1%. La respuesta está desbalanceada y eso condiciona dos decisiones posteriores. Un modelo que dijera que nadie rota acertaría el 83.9% de las veces, de modo que la exactitud dice poco. Y como las probabilidades estimadas van a quedar cerca de 0.16, el corte de 0.5 casi nunca se alcanzaría.

El 28.3% de la planta trabaja horas extra, el 19.7% se declara muy insatisfecho con su trabajo y el 32.0% es soltero.

Tabla 2. Medidas descriptivas de las variables cuantitativas (ingreso en miles).

asimetria <- function(x) mean((x - mean(x))^3) / sd(x)^3

uni_cuant <- df %>%
  select(all_of(cuant_vars)) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  group_by(Variable) %>%
  summarise(Media = mean(valor), Mediana = median(valor), DE = sd(valor),
            `CV (%)` = 100 * sd(valor) / mean(valor),
            Mín = min(valor), Q1 = quantile(valor, 0.25),
            Q3 = quantile(valor, 0.75), Máx = max(valor),
            Asimetría = asimetria(valor))

tabla(uni_cuant, digits = 2)
Variable Media Mediana DE CV (%) Mín Q1 Q3 Máx Asimetría
Antigüedad_Cargo 4.23 3.00 3.62 85.67 0.00 2.00 7.00 18 0.92
Edad 36.92 36.00 9.14 24.74 18.00 30.00 43.00 60 0.41
Ingreso_Miles 6.50 4.92 4.71 72.40 1.01 2.91 8.38 20 1.37 
largo_cuant <- df %>%
  select(all_of(cuant_vars)) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "valor")

ggplot(largo_cuant, aes(valor)) +
  geom_histogram(bins = 20, fill = col_secundario, color = "white") +
  geom_vline(data = largo_cuant %>% group_by(variable) %>% summarise(m = median(valor)),
             aes(xintercept = m), linetype = "dashed", color = col_primario) +
  facet_wrap(~variable, scales = "free") +
  labs(x = NULL, y = "Empleados")

Figura 2. Histogramas de las variables cuantitativas. La línea punteada es la mediana.

La edad es casi simétrica, con mediana de 36 años. El ingreso y la antigüedad en el cargo tienen cola a la derecha (asimetría de 1.37 y 0.92): la mitad de los empleados gana menos de 4.92 miles y lleva 3 años o menos en su cargo, mientras unos pocos llegan a 20 miles y 18 años. Por esa asimetría comparamos los grupos con una prueba no paramétrica en la sección siguiente. Los valores extremos se revisan en el Anexo 10.1.

5 Análisis bivariado

Cada variable se cruza con la rotación usando la prueba que corresponde a su tipo, todas al 5%.

Para una variable cualitativa con \(r\) categorías usamos la prueba de independencia de Pearson de la Unidad 1. Si \(O_{ij}\) son las frecuencias observadas y \(E_{ij}=n_{i\cdot}\,n_{\cdot j}/n\) las esperadas bajo independencia,

\[ \chi^2=\sum_{i=1}^{r}\sum_{j=1}^{2}\frac{(O_{ij}-E_{ij})^2}{E_{ij}}\;\sim\;\chi^2_{(r-1)}, \qquad V=\sqrt{\frac{\chi^2}{n}} \]

donde \(V\) es la V de Cramér, que mide la fuerza de la asociación entre 0 y 1. Para una variable cuantitativa usamos Wilcoxon-Mann-Whitney, que compara los rangos de los dos grupos y no supone normalidad.

tasas_cat <- map_dfr(cat_vars, function(v) {
  df %>% group_by(categoria = as.character(.data[[v]])) %>%
    summarise(n = n(), tasa = mean(y), .groups = "drop") %>%
    mutate(variable = v)
}) %>%
  mutate(categoria = factor(categoria,
                            levels = c("No", "Si", levels(df$Satisfación_Laboral), levels(df$Estado_Civil))))

ggplot(tasas_cat, aes(categoria, tasa)) +
  geom_col(fill = col_alerta) +
  geom_hline(yintercept = p_rot, linetype = "dashed", color = col_primario) +
  geom_text(aes(label = pct(tasa)), vjust = -0.4, size = 3.2) +
  facet_wrap(~variable, scales = "free_x") +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.12))) +
  labs(x = NULL, y = "Empleados que rotan") +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

Figura 3. Porcentaje de empleados que rota en cada categoría. La línea punteada es la rotación global.

Tabla 3. Prueba chi-cuadrado de independencia frente a la rotación.

pruebas_cat <- map_dfr(cat_vars, function(v) {
  tb <- table(df[[v]], df$Rotación)
  ch <- chisq.test(tb, correct = FALSE)
  tibble(Variable = v, `Chi-cuadrado` = unname(ch$statistic), gl = unname(ch$parameter),
         p = ch$p.value,
         `V de Cramér` = sqrt(unname(ch$statistic) / sum(tb)),
         `Esperada mínima` = min(ch$expected))
})

pruebas_cat %>%
  mutate(`Valor p` = fp(p)) %>%
  select(Variable, `Chi-cuadrado`, gl, `Valor p`, `V de Cramér`, `Esperada mínima`) %>%
  tabla(digits = 2)
Variable Chi-cuadrado gl Valor p V de Cramér Esperada mínima
Horas_Extra 89.04 1 < 0.001 0.25 67.07
Satisfación_Laboral 17.51 3 < 0.001 0.11 45.14
Estado_Civil 46.16 2 < 0.001 0.18 52.72 
df %>%
  select(Rotación, all_of(cuant_vars)) %>%
  pivot_longer(-Rotación, names_to = "variable", values_to = "valor") %>%
  ggplot(aes(Rotación, valor, fill = Rotación)) +
  geom_boxplot(show.legend = FALSE, outlier.size = 0.8, outlier.alpha = 0.5) +
  scale_fill_manual(values = pal_rot) +
  facet_wrap(~variable, scales = "free") +
  labs(x = "¿El empleado rotó?", y = NULL)

Figura 4. Variables cuantitativas según la rotación.

Tabla 4. Comparación entre quienes rotan y quienes no (Wilcoxon-Mann-Whitney).

pruebas_cuant <- map_dfr(cuant_vars, function(v) {
  x <- df[[v]]
  tibble(Variable = v,
         `Mediana (no rota)` = median(x[df$y == 0]),
         `Mediana (rota)`    = median(x[df$y == 1]),
         `Shapiro-Wilk (p)`  = fp(shapiro.test(x)$p.value),
         p = wilcox.test(x ~ df$Rotación)$p.value)
})

pruebas_cuant %>%
  mutate(`Wilcoxon (p)` = fp(p)) %>% select(-p) %>%
  tabla(digits = 2)
Variable Mediana (no rota) Mediana (rota) Shapiro-Wilk (p) Wilcoxon (p)
Edad 36.0 32.0 < 0.001 < 0.001
Ingreso_Miles 5.2 3.2 < 0.001 < 0.001
Antigüedad_Cargo 3.0 2.0 < 0.001 < 0.001

Las seis pruebas rechazan la hipótesis nula. Las frecuencias esperadas superan 5 en todas las celdas, así que la aproximación \(\chi^2\) es válida, y Shapiro-Wilk rechaza la normalidad de las tres cuantitativas, lo que justifica Wilcoxon. La asociación más fuerte es la de las horas extra (V = 0.25): rota el 30.5% de quienes las trabajan frente al 10.4% del resto. Los solteros rotan el doble que los casados (25.5% frente a 12.5%). En satisfacción la diferencia está en los extremos, 22.8% en los muy insatisfechos y 11.3% en los muy satisfechos, con los dos niveles intermedios casi iguales. Quienes rotan son más jóvenes, ganan menos (mediana de 3.20 frente a 5.20 miles) y llevan menos tiempo en el cargo.

5.1 Signo del coeficiente y contraste con las hipótesis

Para leer el signo estimamos una regresión logística con cada variable por separado. En el caso de una variable binaria, el resultado se puede anticipar desde la tabla de contingencia. Con horas extra el modelo es \(\ln\frac{P}{1-P}=\beta_0+\beta_1x\), con \(x=1\) si el empleado las trabaja, y entonces

\[ \beta_1=\ln\frac{P_1}{1-P_1}-\ln\frac{P_0}{1-P_0} \quad\Longrightarrow\quad e^{\beta_1}=\frac{P_1/(1-P_1)}{P_0/(1-P_0)}=\frac{a\,d}{b\,c} \]

donde \(a\) y \(b\) son los que rotan y no rotan entre quienes trabajan horas extra, y \(c\) y \(d\) los mismos conteos entre quienes no. El odds ratio de la tabla debe coincidir con \(e^{\hat\beta_1}\).

tb_he <- table(df$Horas_Extra, df$y)
a <- tb_he["Si", "1"]; b <- tb_he["Si", "0"]; c0 <- tb_he["No", "1"]; d <- tb_he["No", "0"]

or_tabla <- (a * d) / (b * c0)
or_logit <- exp(coef(glm(y ~ Horas_Extra, family = binomial, data = df))[2])

round(c(a = a, b = b, c = c0, d = d, OR_tabla = or_tabla, OR_logit = unname(or_logit)), 3)
##        a        b        c        d OR_tabla OR_logit 
##  127.000  289.000  110.000  944.000    3.771    3.771

Se cumple: \(\frac{127\times944}{289\times110}=\) 3.771, igual al \(e^{\hat\beta_1}\) del modelo.

Tabla 5. Regresiones logísticas simples y contraste con el signo esperado.

esperado <- c(Horas_Extra = 1, Satisfación_Laboral = -1, Estado_Civil = 1,
              Edad = -1, Ingreso_Miles = -1, Antigüedad_Cargo = -1)

biv <- map_dfr(vars, function(v) {
  m <- glm(as.formula(paste("y ~", v)), family = binomial, data = df)
  tidy(m) %>%
    filter(term != "(Intercept)") %>%
    transmute(Variable = v, term, Término = str_remove(term, v),
              Coeficiente = estimate, OR = exp(estimate), p = p.value)
}) %>%
  mutate(Término = ifelse(Término == "", "(por unidad)", Término),
         `Signo esperado` = ifelse(esperado[Variable] > 0, "+", "−"),
         `Signo estimado` = ifelse(Coeficiente > 0, "+", "−"),
         Significativo    = ifelse(p < 0.05, "Sí", "No"))

biv %>%
  mutate(`Valor p` = fp(p)) %>%
  select(Variable, Término, Coeficiente, OR, `Valor p`, `Signo esperado`, `Signo estimado`, Significativo) %>%
  tabla(digits = 3) %>%
  collapse_rows(columns = 1, valign = "top")
Variable Término Coeficiente OR Valor p Signo esperado Signo estimado Significativo
Horas_Extra Si 1.327 3.771 < 0.001
Sí
Satisfación_Laboral Insatisfecho -0.409 0.664 0.055 − − No
Satisfecho -0.403 0.668 0.034 − − Sí
Muy satisfecho -0.840 0.432 < 0.001 − − Sí
Estado_Civil Divorciado -0.239 0.787 0.271
− No
Soltero 0.877 2.404 < 0.001
Sí
Edad (por unidad) -0.052 0.949 < 0.001 − − Sí
Ingreso_Miles (por unidad) -0.127 0.881 < 0.001 − − Sí
Antigüedad_Cargo (por unidad) -0.146 0.864 < 0.001 − − Sí

Los signos coinciden con las seis hipótesis. Trabajar horas extra y ser soltero aumentan la probabilidad de rotar; la satisfacción, la edad, el ingreso y la antigüedad en el cargo la bajan. El único término no significativo es “divorciado” (p = 0.271), que no se distingue de “casado”. Como varias de estas variables están relacionadas entre sí, falta ver si los efectos se sostienen cuando entran juntas.

6 Estimación del modelo

6.1 Planteamiento y estimación

La rotación de cada empleado es \(Y_i\sim\text{Bernoulli}(P_i)\). La regresión lineal de la Unidad 2 no aplica porque \(x_i'\beta\) puede salirse de \([0,1]\), y por eso se modela el logaritmo de los odds:

\[ \ln\frac{P_i}{1-P_i}=x_i'\beta \quad\Longleftrightarrow\quad P_i=\frac{1}{1+e^{-x_i'\beta}} \]

con \(x_i'\beta=\beta_0+\beta_1HE_{Si}+\beta_2Sat_{Insat}+\beta_3Sat_{Satisf}+\beta_4Sat_{Muy\,satisf}+\beta_5EC_{Divorc}+\beta_6EC_{Soltero}+\beta_7Edad+\beta_8Ingreso+\beta_9AntCargo\). Las cualitativas entran como variables indicadoras frente a su categoría de referencia.

Los parámetros se estiman por máxima verosimilitud. La log-verosimilitud y su gradiente son

\[ \ell(\beta)=\sum_{i=1}^{n}\left[y_i\,x_i'\beta-\ln\!\left(1+e^{x_i'\beta}\right)\right], \qquad \frac{\partial\ell}{\partial\beta}=X'(y-p)=0 \]

Estas ecuaciones se parecen a las ecuaciones normales de MCO, \(X'(y-X\beta)=0\), pero \(p\) depende de \(\beta\) de forma no lineal y no tienen solución cerrada. Se resuelven con Newton-Raphson, que en cada paso equivale a un mínimo cuadrado ponderado con \(W=\text{diag}\{p_i(1-p_i)\}\):

\[ \beta^{(t+1)}=\beta^{(t)}+\left(X'WX\right)^{-1}X'(y-p), \qquad \widehat{V}(\hat\beta)=\left(X'WX\right)^{-1} \]

Igual que en la Unidad 2 comparamos el MCO matricial con lm(), aquí programamos el algoritmo y lo comparamos con glm().

modelo <- glm(y ~ Horas_Extra + Satisfación_Laboral + Estado_Civil +
                Edad + Ingreso_Miles + Antigüedad_Cargo,
              family = binomial(link = "logit"), data = df)
X  <- model.matrix(modelo)
yv <- df$y
beta <- rep(0, ncol(X))

for (iter in 1:25) {
  p_i  <- as.vector(1 / (1 + exp(-X %*% beta)))
  W    <- p_i * (1 - p_i)
  paso <- solve(t(X) %*% (W * X), t(X) %*% (yv - p_i))
  beta <- beta + as.vector(paso)
  if (max(abs(paso)) < 1e-10) break
}

p_i   <- as.vector(1 / (1 + exp(-X %*% beta)))
V_man <- solve(t(X) %*% (p_i * (1 - p_i) * X))
logv  <- sum(yv * (X %*% beta) - log(1 + exp(X %*% beta)))

dif_coef <- max(abs(beta - coef(modelo)))
dif_ee   <- max(abs(sqrt(diag(V_man)) - sqrt(diag(vcov(modelo)))))

data.frame(iteraciones = iter,
           log_verosimilitud = round(logv, 4),
           logLik_glm = round(as.numeric(logLik(modelo)), 4),
           dif_max_coef = formatC(dif_coef, format = "e", digits = 1),
           dif_max_ee   = formatC(dif_ee, format = "e", digits = 1))

El algoritmo converge en 7 iteraciones y reproduce los coeficientes y los errores estándar de glm(); las diferencias aparecen después de la quinta cifra decimal y vienen de que glm() se detiene con una tolerancia menos estricta.

Tabla 6. Coeficientes del modelo, prueba de Wald y odds ratios con intervalo de confianza al 95%.

coefs <- tidy(modelo) %>%
  bind_cols(as_tibble(confint.default(modelo)) %>% set_names(c("li", "ls"))) %>%
  mutate(manual = beta, OR = exp(estimate), or_li = exp(li), or_ls = exp(ls))

coefs %>%
  transmute(Término = term, `Coef. (manual)` = manual, `Coef. (glm)` = estimate,
            `Error est.` = std.error, z = statistic, `Valor p` = fp(p.value),
            OR, `IC 2.5%` = or_li, `IC 97.5%` = or_ls) %>%
  tabla(digits = 3)
Término Coef. (manual) Coef. (glm) Error est. z Valor p OR IC 2.5% IC 97.5%
(Intercept) -0.165 -0.165 0.375 -0.441 0.659 0.848 0.406 1.768
Horas_ExtraSi 1.516 1.516 0.160 9.497 < 0.001 4.554 3.330 6.227
Satisfación_LaboralInsatisfecho -0.465 -0.465 0.236 -1.972 0.049 0.628 0.396 0.997
Satisfación_LaboralSatisfecho -0.480 -0.480 0.208 -2.303 0.021 0.619 0.412 0.931
Satisfación_LaboralMuy satisfecho -1.079 -1.079 0.223 -4.843 < 0.001 0.340 0.220 0.526
Estado_CivilDivorciado -0.346 -0.346 0.230 -1.503 0.133 0.707 0.450 1.111
Estado_CivilSoltero 0.828 0.828 0.172 4.821 < 0.001 2.288 1.634 3.204
Edad -0.028 -0.028 0.010 -2.834 0.005 0.972 0.953 0.991
Ingreso_Miles -0.072 -0.072 0.025 -2.866 0.004 0.930 0.885 0.977
Antigüedad_Cargo -0.099 -0.099 0.027 -3.673 < 0.001 0.905 0.859 0.955 
coefs %>%
  filter(term != "(Intercept)") %>%
  mutate(term = fct_reorder(term, OR),
         efecto = case_when(or_li > 1 ~ "Aumenta la rotación",
                            or_ls < 1 ~ "Reduce la rotación",
                            TRUE ~ "No significativo")) %>%
  ggplot(aes(OR, term, color = efecto)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey40") +
  geom_errorbarh(aes(xmin = or_li, xmax = or_ls), height = 0.2, linewidth = 0.8) +
  geom_point(size = 3) +
  scale_x_log10(breaks = c(0.25, 0.5, 1, 2, 4, 6)) +
  scale_color_manual(values = c("Aumenta la rotación" = col_alerta,
                                "Reduce la rotación" = col_primario,
                                "No significativo" = "grey55")) +
  labs(x = "Odds ratio (IC 95%)", y = NULL, color = NULL) +
  theme(legend.position = "bottom")

Figura 5. Odds ratios del modelo con su intervalo al 95%, en escala logarítmica.

6.2 Interpretación y significancia

Si \(x_j\) aumenta en una unidad y lo demás queda fijo, el logit sube \(\beta_j\) y los odds se multiplican por \(e^{\beta_j}\):

\[ \frac{odds(x_j+1)}{odds(x_j)}=\frac{e^{\beta_0+\dots+\beta_j(x_j+1)+\dots}}{e^{\beta_0+\dots+\beta_jx_j+\dots}}=e^{\beta_j} \]

Cada coeficiente se contrasta con la prueba de Wald, \(z_j=\hat\beta_j/ee(\hat\beta_j)\sim N(0,1)\) bajo \(H_0:\beta_j=0\), y el intervalo del OR es \(\exp\{\hat\beta_j\pm1.96\,ee(\hat\beta_j)\}\). Para decidir si una variable completa aporta, o el modelo entero, se usa la razón de verosimilitud, que hace aquí el papel de la prueba F parcial de la Unidad 2:

\[ G=-2\left[\ell(\hat\beta_{reducido})-\ell(\hat\beta_{completo})\right]\sim\chi^2_{(q)} \]

con \(q\) igual al número de parámetros que se anulan.

Con las demás variables fijas, trabajar horas extra multiplica los odds de rotar por 4.55 (IC 95%: 3.33 a 6.23). Es el efecto más grande del modelo, y es mayor que el 3.77 del análisis bivariado, así que las otras variables no lo explican. Los solteros tienen odds 2.29 veces los de los casados (IC 95%: 1.63 a 3.20), y los divorciados no se diferencian de los casados (p = 0.133).

Frente a un empleado muy insatisfecho, los odds son 37% menores en uno insatisfecho, 38% menores en uno satisfecho y 66% menores en uno muy satisfecho. El nivel “insatisfecho” queda en el borde (p = 0.049) y lo leemos con cautela.

Cada año de edad reduce los odds en 3%, cada 1000 de ingreso en 7% y cada año en el cargo en 9%. Son reducciones menores que las de las regresiones simples (5%, 12% y 14%), lo esperable cuando las tres están correlacionadas y antes cada una cargaba con parte del efecto de las otras.

Tabla 7. Razón de verosimilitud por covariable y para el modelo completo.

anova_lr <- car::Anova(modelo, type = "II", test.statistic = "LR")

# Prueba global calculada con la log-verosimilitud del modelo nulo:
# l0 = n [ ybar ln(ybar) + (1 - ybar) ln(1 - ybar) ]
l0  <- nrow(df) * (p_rot * log(p_rot) + (1 - p_rot) * log(1 - p_rot))
G   <- -2 * (l0 - logv)
gl  <- ncol(X) - 1
r2_mcfadden <- 1 - logv / l0

an <- as.data.frame(anova_lr) %>%
  rownames_to_column("Contraste") %>%
  set_names(c("Contraste", "G", "gl", "p")) %>%
  arrange(desc(G)) %>%
  bind_rows(tibble(Contraste = "Modelo completo (todas las β = 0)", G = G, gl = gl,
                   p = pchisq(G, gl, lower.tail = FALSE)))

an %>% mutate(`Valor p` = fp(p)) %>% select(Contraste, G, gl, `Valor p`) %>% tabla(digits = 2)
Contraste G gl Valor p
Horas_Extra 92.22 1 < 0.001
Estado_Civil 37.01 2 < 0.001
Satisfación_Laboral 24.52 3 < 0.001
Antigüedad_Cargo 14.40 1 < 0.001
Ingreso_Miles 8.93 1 0.003
Edad 8.37 1 0.004
Modelo completo (todas las β = 0) 220.55 9 < 0.001

Las seis covariables son significativas al 5%, y Estado_Civil lo es como variable aunque “divorciado” no lo sea como categoría. Para el modelo completo, \(\ell(\hat\beta_0)\) = -649.29 y \(\ell(\hat\beta)\) = -539.02, de donde \(G\) = 220.55 con 9 grados de libertad (p < 0.001). Es la misma diferencia de desviaciones que reporta glm() (220.55). El pseudo-\(R^2\) de McFadden, \(1-\ell(\hat\beta)/\ell(\hat\beta_0)\), es 0.170: el modelo explica una parte de la rotación y deja espacio para otras variables de la base. El VIF del Anexo 10.2 descarta multicolinealidad.

7 Evaluación del poder predictivo

Evaluar el modelo con los mismos datos que lo estimaron lo favorece. Como en la validación cruzada de la Unidad 2, separamos al azar el 70% de la base para estimar y el 30% para evaluar, con la misma proporción de rotación en las dos partes.

set.seed(123)
idx <- unlist(lapply(split(seq_len(nrow(df)), df$y),
                     function(i) sample(i, round(0.7 * length(i)))))
train <- df[idx, ]
test  <- df[-idx, ]

modelo_tr <- glm(formula(modelo), family = binomial, data = train)
p_train   <- predict(modelo_tr, type = "response")
p_test    <- predict(modelo_tr, newdata = test, type = "response")

c(n_entrenamiento = nrow(train), rotan_entrenamiento = sum(train$y),
  n_prueba = nrow(test), rotan_prueba = sum(test$y))
##     n_entrenamiento rotan_entrenamiento            n_prueba        rotan_prueba 
##                1029                 166                 441                  71

Para un corte \(c\) se interviene al empleado si \(\hat P\ge c\), y de ahí salen

\[ \text{Sensibilidad}(c)=P(\hat P\ge c\mid Y=1),\qquad \text{Especificidad}(c)=P(\hat P<c\mid Y=0) \]

La curva ROC dibuja la sensibilidad contra \(1-\)especificidad para todos los cortes, y el AUC es el área bajo esa curva. Tiene una lectura directa: es la probabilidad de que un empleado que rota reciba una probabilidad estimada mayor que uno que no rota. Por eso se puede calcular sin la curva, con el estadístico \(U\) de Mann-Whitney:

\[ AUC=P\!\left(\hat P_{rota}>\hat P_{no\ rota}\right)=\frac{U}{n_1\,n_0} \]

roc_train <- roc(train$y, p_train, quiet = TRUE)
roc_test  <- roc(test$y,  p_test,  quiet = TRUE)

auc_train <- as.numeric(auc(roc_train))
auc_test  <- as.numeric(auc(roc_test))
ic_auc    <- as.numeric(ci.auc(roc_test))   # IC 95% de DeLong

# Comprobación: AUC = U / (n1 * n0)
n1 <- sum(test$y == 1); n0 <- sum(test$y == 0)
U  <- unname(wilcox.test(p_test[test$y == 1], p_test[test$y == 0])$statistic)
auc_manual <- U / (n1 * n0)

ggroc(list(Entrenamiento = roc_train, Prueba = roc_test), linewidth = 1) +
  geom_abline(intercept = 1, slope = 1, linetype = "dashed", color = "grey50") +
  scale_color_manual(values = c(Entrenamiento = col_secundario, Prueba = col_primario),
                     labels = c(paste0("Entrenamiento (AUC = ", f2(auc_train, 3), ")"),
                                paste0("Prueba (AUC = ", f2(auc_test, 3), ")"))) +
  labs(x = "Especificidad", y = "Sensibilidad", color = NULL) +
  coord_equal() +
  theme(legend.position = "bottom")

Figura 6. Curva ROC en entrenamiento y prueba. La diagonal es un modelo que clasifica al azar.

En la muestra de prueba, \(U\) = 20404 con \(n_1\) = 71 y \(n_0\) = 370, lo que da un AUC de 0.777, el mismo valor que entrega pROC. Su intervalo de DeLong al 95% va de 0.715 a 0.838 y no toca el 0.5 del azar. En la escala de Hosmer y Lemeshow es una discriminación aceptable. El AUC de entrenamiento es 0.779, casi igual, así que no hay sobreajuste; la validación cruzada del Anexo 10.3 da lo mismo. El modelo sirve para ordenar a los empleados por riesgo, y con seis variables no alcanza para dar por seguro lo que hará una persona en particular.

8 Predicción y punto de corte

8.1 Sensibilidad del punto de corte

El corte se puede deducir antes de mirar los datos. Sea \(C_{FN}\) el costo de no intervenir a alguien que rota y \(C_{FP}\) el de intervenir a alguien que no iba a rotar. Para un empleado con probabilidad \(\hat P\), no intervenir tiene un costo esperado de \(\hat P\,C_{FN}\) e intervenir uno de \((1-\hat P)\,C_{FP}\). Conviene intervenir cuando

\[ \hat P\,C_{FN}\ \ge\ (1-\hat P)\,C_{FP} \quad\Longleftrightarrow\quad \hat P\ \ge\ c^{*}=\frac{C_{FP}}{C_{FP}+C_{FN}}=\frac{1}{1+k},\qquad k=\frac{C_{FN}}{C_{FP}} \]

El corte de 0.5 solo es óptimo si los dos errores cuestan lo mismo (\(k=1\)). Perder a un empleado cuesta más que una conversación de seguimiento, así que \(k>1\) y el corte debe quedar por debajo de 0.5. Cuando no se conocen los costos, el criterio usual es el índice de Youden, que da igual peso a los dos errores en términos de tasas:

\[ J(c)=\text{Sensibilidad}(c)+\text{Especificidad}(c)-1 \]

Elegimos el corte que maximiza \(J\) en entrenamiento y lo evaluamos en prueba, para no calificarlo con los datos que lo definieron.

youden <- coords(roc_train, "best", best.method = "youden",
                 ret = c("threshold", "sensitivity", "specificity"))
corte  <- round(youden$threshold[1], 2)
k_impl <- (1 - corte) / corte

c(corte_youden = corte, k_implicito = k_impl)
## corte_youden  k_implicito 
##     0.210000     3.761905
rejilla <- map_dfr(seq(0.02, 0.70, by = 0.01), ~ metricas_corte(p_test, test$y, .x))

rejilla %>%
  select(Corte, Sensibilidad, Especificidad, Exactitud) %>%
  pivot_longer(-Corte, names_to = "Medida", values_to = "valor") %>%
  ggplot(aes(Corte, valor, color = Medida)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = corte, linetype = "dashed", color = "grey30") +
  geom_vline(xintercept = 0.5, linetype = "dotted", color = "grey30") +
  annotate("text", x = corte, y = 1.04, label = paste0("Corte elegido (", corte, ")"), size = 3.2, hjust = -0.05) +
  annotate("text", x = 0.5, y = 1.04, label = "Corte 0.5", size = 3.2, hjust = -0.1) +
  scale_color_manual(values = c(Sensibilidad = col_alerta, Especificidad = col_primario, Exactitud = "#F2B134")) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1.06), breaks = seq(0, 1, 0.25)) +
  labs(x = "Punto de corte", y = NULL, color = NULL) +
  theme(legend.position = "bottom")

Figura 7. Sensibilidad, especificidad y exactitud en la muestra de prueba según el corte.

Tabla 8. Desempeño en la muestra de prueba para distintos cortes. La columna \(k\) es la relación de costos \(C_{FN}/C_{FP}\) para la que ese corte sería el óptimo. VPP: proporción de intervenidos que sí rotan.

cortes_tabla <- sort(unique(c(0.10, 0.15, 0.20, corte, 0.25, 0.30, 0.40, 0.50)))
sens_corte   <- map_dfr(cortes_tabla, ~ metricas_corte(p_test, test$y, .x)) %>%
  mutate(k = (1 - Corte) / Corte, J = Sensibilidad + Especificidad - 1)

sens_corte %>%
  mutate(across(c(Sensibilidad, Especificidad, VPP, Exactitud, Intervenidos), ~ pct(.x)),
         k = f2(k, 1), J = f2(J, 3)) %>%
  select(Corte, k, VP, FP, FN, VN, Sensibilidad, Especificidad, J, VPP, Exactitud, Intervenidos) %>%
  tabla() %>%
  row_spec(which(sens_corte$Corte == corte), bold = TRUE, background = "#E8F1FA")
Corte k VP FP FN VN Sensibilidad Especificidad J VPP Exactitud Intervenidos
0.10 9.0 62 178 9 192 87.3% 51.9% 0.392 25.8% 57.6% 54.4%
0.15 5.7 54 113 17 257 76.1% 69.5% 0.455 32.3% 70.5% 37.9%
0.20 4.0 45 86 26 284 63.4% 76.8% 0.401 34.4% 74.6% 29.7%
0.21 3.8 44 81 27 289 62.0% 78.1% 0.401 35.2% 75.5% 28.3%
0.25 3.0 40 57 31 313 56.3% 84.6% 0.409 41.2% 80.0% 22.0%
0.30 2.3 33 43 38 327 46.5% 88.4% 0.349 43.4% 81.6% 17.2%
0.40 1.5 17 18 54 352 23.9% 95.1% 0.191 48.6% 83.7% 7.9%
0.50 1.0 11 8 60 362 15.5% 97.8% 0.133 57.9% 84.6% 4.3%

El corte de Youden es 0.21, cerca de la proporción de rotación de la base, y equivale a suponer que un falso negativo cuesta 3.8 veces lo que un falso positivo. Con él, el modelo detecta en prueba a 44 de los 71 empleados que rotaron (sensibilidad de 62.0%) y habría que intervenir al 28.3% de la planta. De cada 100 intervenidos, 35 iban a rotar, 2.2 veces lo que se lograría eligiendo al azar.

Con 0.5 la exactitud sube a 84.6%, pero solo se detectan 11 de los 71. Esa exactitud es prácticamente la de decir que nadie rota (83.9%), tal como se anticipó en el análisis univariado. Bajar el corte a 0.10 detecta a casi todos y obliga a intervenir a más de la mitad de la planta. La zona útil está entre 0.15 y 0.30, y dentro de ella la empresa puede mover el corte según su propia relación de costos \(k\).

8.2 Predicción para un empleado hipotético

Tomamos un empleado de 28 años, soltero, que trabaja horas extra, se declara insatisfecho, gana 3000 al mes y lleva un año en el cargo. La predicción sale de reemplazar su perfil \(x_0\) en el modelo estimado con toda la base, y el intervalo se construye en la escala logit, donde \(\hat\eta_0\) es aproximadamente normal:

\[ \hat\eta_0=x_0'\hat\beta,\qquad \hat P_0=\frac{1}{1+e^{-\hat\eta_0}},\qquad ee(\hat\eta_0)=\sqrt{x_0'\,\widehat V(\hat\beta)\,x_0},\qquad IC_{95\%}(P_0)=\frac{1}{1+e^{-(\hat\eta_0\pm1.96\,ee)}} \]

nuevo <- data.frame(
  Horas_Extra         = factor("Si", levels = levels(df$Horas_Extra)),
  Satisfación_Laboral = factor("Insatisfecho", levels = levels(df$Satisfación_Laboral)),
  Estado_Civil        = factor("Soltero", levels = levels(df$Estado_Civil)),
  Edad                = 28,
  Ingreso_Miles       = 3,
  Antigüedad_Cargo    = 1
)

# Cálculo matricial con los coeficientes y la varianza obtenidos a mano
x0      <- as.vector(model.matrix(delete.response(terms(modelo)), nuevo))
eta0    <- sum(x0 * beta)
ee0     <- sqrt(as.numeric(t(x0) %*% V_man %*% x0))
p_nuevo <- 1 / (1 + exp(-eta0))
ic_nuevo <- 1 / (1 + exp(-(eta0 + c(-1, 1) * qnorm(0.975) * ee0)))

decidir <- function(p) ifelse(p >= corte, "Intervenir", "No intervenir")

round(c(eta0 = eta0, ee = ee0, P = p_nuevo, IC_inf = ic_nuevo[1], IC_sup = ic_nuevo[2],
        P_predict_glm = unname(predict(modelo, newdata = nuevo, type = "response"))), 4)
##          eta0            ee             P        IC_inf        IC_sup 
##        0.6000        0.2362        0.6457        0.5342        0.7432 
## P_predict_glm 
##        0.6457

Al reemplazar,

\[ \hat\eta_0=-0.165+1.516-0.465+0.828-0.0285(28)-0.0722(3)-0.0993(1)=0.600 \quad\Longrightarrow\quad \hat P_0=\frac{1}{1+e^{-0.600}}=0.646 \]

La probabilidad de que este empleado rote es 64.6% (IC 95%: 53.4% a 74.3%), 4.0 veces la rotación promedio. Todo el intervalo queda por encima del corte de 0.21, así que la decisión es intervenir.

¿Con qué medida? El mismo modelo deja comparar. Edad y estado civil se quedan fijos porque la empresa no los gestiona.

Tabla 9. Probabilidad de rotación del empleado hipotético según la medida que se le aplique.

escenario <- function(nombre, ...) {
  d <- nuevo
  cambios <- list(...)
  for (v in names(cambios)) {
    d[[v]] <- if (is.factor(d[[v]])) factor(cambios[[v]], levels = levels(d[[v]])) else cambios[[v]]
  }
  tibble(Escenario = nombre, p = as.numeric(predict(modelo, newdata = d, type = "response")))
}

escenarios <- bind_rows(
  escenario("Situación actual"),
  escenario("Aumento salarial de 1000", Ingreso_Miles = 4),
  escenario("Pasa a muy satisfecho", Satisfación_Laboral = "Muy satisfecho"),
  escenario("Deja las horas extra", Horas_Extra = "No"),
  escenario("Deja las horas extra y pasa a muy satisfecho",
            Horas_Extra = "No", Satisfación_Laboral = "Muy satisfecho")
) %>%
  mutate(Decisión = decidir(p))

escenarios %>%
  transmute(Escenario, `Probabilidad de rotar` = pct(p), Decisión) %>%
  tabla()
Escenario Probabilidad de rotar Decisión
Situación actual 64.6% Intervenir
Aumento salarial de 1000 62.9% Intervenir
Pasa a muy satisfecho 49.7% Intervenir
Deja las horas extra 28.6% Intervenir
Deja las horas extra y pasa a muy satisfecho 17.8% No intervenir

Quitarle las horas extra es lo que más baja su riesgo, de 64.6% a 28.6%. El aumento de 1000 casi no lo mueve. Ninguna medida sola lo deja por debajo del corte; hace falta combinar la carga de trabajo con la satisfacción para lograrlo. Estas cifras describen a empleados con ese perfil en datos observados. No prueban que la medida cause la reducción, y eso habría que comprobarlo con un piloto.

9 Conclusiones y estrategia

Las seis hipótesis se confirmaron, tanto variable por variable como en el modelo conjunto. Las horas extra son el factor que más pesa (OR = 4.55), por encima del salario, y le sigue ser soltero (OR = 2.29). El modelo discrimina de forma aceptable (AUC = 0.777 en datos que no vio) y con el corte de 0.21 detecta al 62% de quienes rotan revisando al 28% de la planta. Tiene límites claros: cerca de 65 de cada 100 empleados señalados no iban a rotar, usa seis de las 23 variables disponibles y lo que mide son asociaciones.

Con base en las variables significativas, la estrategia que proponemos es esta:

  1. Empezar por las horas extra. El 28.3% de la planta las trabaja y en ese grupo rota el 30.5%. Hay que ubicar los cargos donde son habituales, reforzarlos con más personal o más turnos y poner un tope mensual por persona.
  2. Atender a los muy insatisfechos, donde rota el 22.8%. Medir la satisfacción con regularidad y conversar con ese grupo para saber si el problema es el jefe, las tareas o la falta de crecimiento. Salir del nivel más bajo ya reduce los odds en 37%.
  3. Revisar los salarios más bajos y no subirlos a todos. Cada 1000 adicionales reducen los odds apenas 7%, de modo que un aumento general sería caro para lo que consigue; la rotación está en la franja de menores ingresos.
  4. Acompañar los dos primeros años en el cargo, con inducción y seguimiento. Quienes rotan llevan una mediana de 2 años en el puesto.
  5. Usar la edad y el estado civil solo para entender al grupo que más rota. Los empleados jóvenes y solteros valoran el plan de carrera, la formación y la movilidad interna. Estas dos variables no deben pesar en contrataciones, ascensos ni salarios.
  6. Calcular cada período la probabilidad de rotación de cada empleado y hacer seguimiento a quienes superen 0.21. Cuando la empresa conozca sus costos de reemplazo y de intervención, el corte se ajusta con \(c^{*}=1/(1+k)\).

10 Anexos

10.1 Valores atípicos y observaciones influyentes

Tabla 10. Valores atípicos por variable según la regla \([Q_1-1.5\,IQR,\ Q_3+1.5\,IQR]\).

atipicos <- df_full %>%
  select(where(is.numeric), -y, -Ingreso_Miles) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  group_by(Variable) %>%
  summarise(Q1 = quantile(valor, 0.25), Q3 = quantile(valor, 0.75),
            `Límite inferior` = Q1 - 1.5 * (Q3 - Q1),
            `Límite superior` = Q3 + 1.5 * (Q3 - Q1),
            Atípicos = sum(valor < `Límite inferior` | valor > `Límite superior`),
            `% atípicos` = 100 * Atípicos / n()) %>%
  arrange(desc(`% atípicos`))

tabla(atipicos, digits = 1)
Variable Q1 Q3 Límite inferior Límite superior Atípicos % atípicos
Capacitaciones 2 3 0.5 4.5 238 16.2
Ingreso_Mensual 2911 8379 -5291.0 16581.0 114 7.8
Años_ultima_promoción 0 3 -4.5 7.5 107 7.3
Antigüedad 3 9 -6.0 18.0 104 7.1
Años_Experiencia 6 15 -7.5 28.5 63 4.3
Trabajos_Anteriores 1 4 -3.5 8.5 52 3.5
Antigüedad_Cargo 2 7 -5.5 14.5 21 1.4
Años_acargo_con_mismo_jefe 2 7 -5.5 14.5 14 1.0
Distancia_Casa 2 14 -16.0 32.0 0 0.0
Edad 30 43 10.5 62.5 0 0.0
Porcentaje_aumento_salarial 12 18 3.0 27.0 0 0.0 
cook <- cooks.distance(modelo)

tibble(obs = seq_along(cook), cook = cook) %>%
  ggplot(aes(obs, cook)) +
  geom_segment(aes(xend = obs, yend = 0), color = col_primario, alpha = 0.6) +
  labs(x = "Observación", y = "Distancia de Cook")

Figura 8. Distancia de Cook de cada observación en el modelo logístico.

8 de las 11 variables cuantitativas tienen valores atípicos. Corresponden a salarios altos y trayectorias largas, que son datos plausibles, y no los eliminamos. Tampoco afectan el modelo: la mayor distancia de Cook es 0.0196 (observación 46), muy lejos del umbral de 1.

10.2 Multicolinealidad

Tabla 11. Factor de inflación de varianza generalizado.

vif_mod <- car::vif(modelo)

as.data.frame(vif_mod) %>%
  rownames_to_column("Variable") %>%
  set_names(c("Variable", "GVIF", "gl", "GVIF^(1/(2·gl))")) %>%
  tabla(digits = 3)
Variable GVIF gl GVIF^(1/(2·gl))
Horas_Extra 1.045 1 1.022
Satisfación_Laboral 1.028 3 1.005
Estado_Civil 1.036 2 1.009
Edad 1.272 1 1.128
Ingreso_Miles 1.333 1 1.154
Antigüedad_Cargo 1.100 1 1.049

El mayor GVIF es 1.33, muy por debajo del umbral de 10 de la Unidad 2. No hay multicolinealidad.

10.3 Validación cruzada

Tabla 12. Validación cruzada con 10 pliegues estratificados.

set.seed(2026)
k_pl <- 10
pliegue <- integer(nrow(df))
for (clase in unique(df$y)) {
  i <- which(df$y == clase)
  pliegue[i] <- sample(rep(1:k_pl, length.out = length(i)))
}

cv <- map_dfr(1:k_pl, function(j) {
  m  <- glm(formula(modelo), family = binomial, data = df[pliegue != j, ])
  pj <- predict(m, newdata = df[pliegue == j, ], type = "response")
  yj <- df$y[pliegue == j]
  met <- metricas_corte(pj, yj, corte)
  tibble(AUC = as.numeric(auc(roc(yj, pj, quiet = TRUE))),
         Sensibilidad = met$Sensibilidad, Especificidad = met$Especificidad)
})

cv %>%
  pivot_longer(everything(), names_to = "Medida", values_to = "v") %>%
  group_by(Medida) %>%
  summarise(Media = mean(v), `Desviación estándar` = sd(v), Mínimo = min(v), Máximo = max(v)) %>%
  tabla(digits = 3)
Medida Media Desviación estándar Mínimo Máximo
AUC 0.773 0.048 0.702 0.846
Especificidad 0.821 0.032 0.764 0.862
Sensibilidad 0.641 0.072 0.542 0.792

El AUC promedio de los 10 pliegues es 0.773, prácticamente el de la muestra de prueba (0.777). El resultado no depende de la partición. La sensibilidad y la especificidad se calculan con el corte de 0.21.

10.4 Información de sesión

sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.4.1     pROC_1.19.1          car_3.1-5           
##  [4] carData_3.0-6        corrplot_0.95        lubridate_1.9.5     
##  [7] forcats_1.0.1        stringr_1.6.0        dplyr_1.2.1         
## [10] purrr_1.2.2          readr_2.2.0          tidyr_1.3.2         
## [13] tibble_3.3.1         tidyverse_2.0.0      paqueteMODELOS_0.1.0
## [16] summarytools_1.1.5   knitr_1.52           gridExtra_2.3.1     
## [19] GGally_2.4.0         ggplot2_4.0.3        broom_1.0.13        
## [22] boot_1.3-32         
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.60          bslib_0.12.0       tzdb_0.5.0        
##  [5] vctrs_0.7.3        tools_4.6.1        generics_0.1.4     pkgconfig_2.0.3   
##  [9] checkmate_2.3.4    RColorBrewer_1.1-3 S7_0.2.2           lifecycle_1.0.5   
## [13] compiler_4.6.1     farver_2.1.2       textshaping_1.0.5  rapportools_1.2   
## [17] htmltools_0.5.9    sass_0.4.10        yaml_2.3.12        Formula_1.2-6     
## [21] pillar_1.11.1      jquerylib_0.1.4    MASS_7.3-65        cachem_1.1.0      
## [25] magick_2.9.1       abind_1.4-8        ggstats_0.14.0     tidyselect_1.2.1  
## [29] digest_0.6.39      stringi_1.8.9      reshape2_1.4.5     pander_0.6.6      
## [33] bookdown_0.48      labeling_0.4.3     fastmap_1.2.0      grid_4.6.1        
## [37] cli_3.6.6          magrittr_2.0.5     base64enc_0.1-6    withr_3.0.3       
## [41] scales_1.4.0       backports_1.5.1    timechange_0.4.0   rmarkdown_2.32    
## [45] matrixStats_1.5.0  otel_0.2.0         hms_1.1.4          evaluate_1.0.5    
## [49] tcltk_4.6.1        viridisLite_0.4.3  rlang_1.3.0        Rcpp_1.1.2        
## [53] glue_1.8.1         xml2_1.6.0         svglite_2.2.2      rstudioapi_0.19.0 
## [57] jsonlite_2.0.0     R6_2.6.1           plyr_1.8.9         systemfonts_1.3.2