Apunte 24 — Modelo de Cox y hazard ratios

Analítica de Personas · Semestre otoño 2026 · Semana 9 · Prof. René Gempp

1. Del log-rank al modelo multivariable

En el Apunte 23 aprendiste a comparar curvas de Kaplan-Meier entre grupos con el log-rank test. Pero hay un problema: el log-rank compara curvas sin controlar por otras variables.

Si la curva de Soporte Técnico cae más rápido que la de Desarrollo, ¿es porque el departamento en sí causa más rotación, o porque Soporte Técnico tiene empleados más jóvenes, menos satisfechos, con menores salarios? El log-rank no puede distinguir — es como comparar salarios brutos entre hombres y mujeres sin controlar por cargo ni antigüedad (Clase 5).

El modelo de riesgos proporcionales de Cox (Cox, 1972) resuelve esto: es el análogo de la regresión múltiple para datos de supervivencia. Permite estimar el efecto de cada covariable sobre el riesgo de salida, controlando por las demás.

2. coxph(): sintaxis y ajuste

La sintaxis de coxph() es idéntica a la de lm() y glm(), con la diferencia de que el lado izquierdo de la fórmula es un objeto Surv():

# Modelo de Cox
modelo_cox <- coxph(
  Surv(tiempo_meses, evento) ~
    satisfaccion_laboral +
    evaluacion_desempeno +
    departamento +
    genero +
    ingreso_mensual +
    horas_capacitacion +
    nivel_jerarquico,
  data = datos_surv
)

# Resumen del modelo
summary(modelo_cox)

Call:
coxph(formula = Surv(tiempo_meses, evento) ~ satisfaccion_laboral + ...

  n= 1200, number of events= 326

                          coef  exp(coef)  se(coef)    z     Pr(>|z|)
satisfaccion_laboral   -0.1523   0.8588    0.0412  -3.696   0.000219 ***
evaluacion_desempeno   -0.0847   0.9188    0.0398  -2.128   0.033345 *
departamentoFinanzas    0.0912   1.0955    0.2148   0.425   0.671120
departamentoMarketing   0.1574   1.1705    0.2534   0.621   0.534520
departamentoRRHH        0.3821   1.4654    0.1945   1.965   0.049420 *
departamentoSoporte     0.4265   1.5320    0.1621   2.631   0.008512 **
departamentoVentas      0.2041   1.2264    0.1589   1.284   0.199020
generoMujer            -0.0312   0.9693    0.1145  -0.273   0.785190
ingreso_mensual        -0.0001   0.9999    0.0001  -1.412   0.157900
...

Concordance= 0.621  (se = 0.016)
Likelihood ratio test= 42.3  on 11 df,   p=1.5e-05

Las columnas clave del output son:

3. Hazard ratio: interpretación

El hazard ratio (HR) es la medida de efecto del modelo de Cox. Es el exponencial del coeficiente β: HR = exp(β). Para obtener los HR directamente con broom::tidy():

library(broom)

hr_tabla <- tidy(modelo_cox,
                 exponentiate = TRUE,   # exp(β) = HR
                 conf.int     = TRUE) |>
  select(term, HR = estimate, IC_bajo = conf.low,
         IC_alto = conf.high, p_valor = p.value)

print(hr_tabla)

¿Cómo se interpreta un hazard ratio?

Variables continuas vs. categóricas

Para una variable continua (como satisfaccion_laboral, escala 1–5), el HR se interpreta por cada unidad de aumento. Si HR = 0,86, cada punto adicional en la escala reduce el riesgo en un 14 %. Dos puntos de diferencia: 0,86² ≈ 0,74 → reduce el riesgo en un 26 %.

Para una variable categórica (como departamento), R automáticamente crea variables dummy. El HR compara cada categoría con la categoría de referencia (por defecto, la primera en orden alfabético). Verifica cuál es la referencia en tu output.

4. OR vs. HR: cuándo divergen

En la Clase 3 usaste odds ratios (OR) de la regresión logística. Ahora tienes hazard ratios (HR) del modelo de Cox. Son métricas conceptualmente distintas:

En InnovaCo, la tasa de rotación es 27,2 %. Eso no es un evento raro. Por tanto, el OR sobreestimará el efecto relativo comparado con el HR. Para ver esto empíricamente, puedes estimar ambos modelos con las mismas covariables y comparar lado a lado:

# Modelo logístico (mismo que Clase 3)
modelo_log <- glm(evento ~ satisfaccion_laboral + evaluacion_desempeno +
                          departamento + genero + ingreso_mensual,
                  data = datos_surv, family = binomial)

# Tabla comparativa
or_tabla <- tidy(modelo_log, exponentiate = TRUE, conf.int = TRUE) |>
  select(term, OR = estimate)

hr_tabla2 <- tidy(modelo_cox, exponentiate = TRUE, conf.int = TRUE) |>
  select(term, HR = estimate)

left_join(or_tabla, hr_tabla2, by = "term")

5. El forest plot: visualizar los hazard ratios

El forest plot es la visualización estándar de los hazard ratios: una fila por variable, con un punto (HR) y una barra horizontal (IC 95 %). La línea vertical en HR = 1 es la referencia: si la barra cruza el 1, el efecto no es significativo al 5 %.

Forest plot manual con ggplot2 (recomendado)

Partimos de la tabla de hazard ratios generada con broom::tidy() y la visualizamos directamente. Este enfoque es más robusto y nos da control total sobre la estética:

# 1. Limpiar los nombres de coeficientes para legibilidad
hr_plot <- hr_tabla |>
  mutate(
    etiqueta = term |>
      str_replace("^departamento", "Depto: ") |>
      str_replace("^nivel_jerarquico", "Nivel: ") |>
      str_replace("^genero", "Género: ") |>
      str_replace("satisfaccion_laboral", "Satisfacción laboral") |>
      str_replace("evaluacion_desempeno", "Evaluación desempeño") |>
      str_replace("ingreso_mensual", "Ingreso mensual") |>
      str_replace("horas_capacitacion", "Horas capacitación"),
    direccion = if_else(HR > 1, "Riesgo", "Protector")
  ) |>
  mutate(etiqueta = fct_reorder(etiqueta, HR))

# 2. Forest plot
ggplot(hr_plot, aes(x = HR, y = etiqueta, color = direccion)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = IC_bajo, xmax = IC_alto), height = 0.25) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
  scale_color_manual(
    values = c("Riesgo" = "#B85042", "Protector" = "#0D9488"),
    guide  = "none"
  ) +
  labs(
    title = "¿Qué factores aceleran o frenan la salida?",
    subtitle = "Hazard ratios con IC 95 % — Modelo de Cox",
    x = "Hazard Ratio", y = NULL,
    caption = "HR > 1 = mayor riesgo · HR < 1 = menor riesgo"
  ) +
  theme_minimal() +
  theme(panel.grid.major.y = element_blank())

6. cox.zph(): verificar proporcionalidad

El supuesto fundamental del modelo de Cox es que los hazard ratios son constantes en el tiempo. Si el HR de Soporte Técnico es 1,53, debe serlo a los 6 meses y a los 60 meses. Si el efecto cambia con el tiempo, el modelo viola su supuesto principal.

El test de Grambsch y Therneau (1994) evalúa esto formalmente:

# Test de proporcionalidad
test_ph <- cox.zph(modelo_cox)
print(test_ph)

                          chisq  df       p
satisfaccion_laboral      5.23   1    0.022
evaluacion_desempeno      0.87   1    0.351
departamentoFinanzas      0.12   1    0.725
...
GLOBAL                   14.5  11    0.205

Un p-valor significativo (p < 0,05) para una variable indica que su HR no es constante en el tiempo. El test GLOBAL evalúa el supuesto para todo el modelo conjuntamente.

# Visualizar los residuos de Schoenfeld
ggcoxzph(test_ph)

Los gráficos de residuos de Schoenfeld muestran el efecto estimado de cada variable a lo largo del tiempo. Si la línea ajustada es horizontal, el supuesto se cumple. Si tiene pendiente, el efecto cambia con el tiempo.

¿Qué hacer si hay violación?

# Ejemplo de modelo estratificado
# Si "departamento" viola proporcionalidad:
modelo_strat <- coxph(
  Surv(tiempo_meses, evento) ~
    satisfaccion_laboral + evaluacion_desempeno +
    genero + ingreso_mensual +
    strata(departamento),  # ← no estima HR para depto, pero lo controla
  data = datos_surv
)
summary(modelo_strat)

7. Concordance index (C-index)

El C-index (o concordancia) es la medida de capacidad discriminativa del modelo de Cox. Es conceptualmente análogo al AUC de la regresión logística: mide qué tan bien el modelo distingue entre personas que se van rápido y personas que se quedan.

# Extraer el C-index
summary(modelo_cox)$concordance

En datos de rotación, C-indexes de 0,60–0,65 son habituales: los factores que medimos (satisfacción, salario, departamento) explican parte de la rotación, pero hay factores no observables (ofertas externas, situación personal, relación con el jefe directo) que el modelo no captura.

8. Curvas ajustadas por Cox

Las curvas de Kaplan-Meier del Apunte 23 son "crudas": muestran la supervivencia de cada grupo sin controlar por nada. Pero a veces quieres ver la curva de un grupo después de controlar por las demás variables. survminer::ggadjustedcurves() hace esto:

# Curvas de supervivencia ajustadas por departamento
ggadjustedcurves(
  modelo_cox,
  data     = datos_surv,
  variable = "departamento",
  method   = "average"
) +
  labs(
    title = "Supervivencia ajustada por departamento",
    subtitle = "Controlando por satisfacción, desempeño, salario y género",
    x = "Tiempo (meses)",
    y = "Probabilidad de permanencia"
  )

La diferencia entre la curva cruda (Kaplan-Meier) y la ajustada (Cox) te dice cuánto de la diferencia entre grupos se explica por las covariables. Si las curvas ajustadas convergen (se acercan), las diferencias crudas se debían a confounders. Si siguen separadas, hay algo propio del grupo que no capturan las covariables.

9. Extensiones que quedan fuera

El modelo de Cox es la herramienta estándar, pero la estadística de supervivencia moderna tiene extensiones poderosas que este curso no cubre. Las nombramos para que sepas dónde buscar si las necesitas:

10. Errores frecuentes