Taller 2: Cox PH.

1. Paquetes necesarios.

library(naniar)

## Warning: package 'naniar' was built under R version 4.4.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.4.3

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(survival)

## Warning: package 'survival' was built under R version 4.4.3

library(survminer)

## Warning: package 'survminer' was built under R version 4.4.3

## Cargando paquete requerido: ggplot2

## Cargando paquete requerido: ggpubr

## Warning: package 'ggpubr' was built under R version 4.4.2

## 
## Adjuntando el paquete: 'survminer'

## The following object is masked from 'package:survival':
## 
##     myeloma

library(MASS)

## 
## Adjuntando el paquete: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(rms)

## Warning: package 'rms' was built under R version 4.4.3

## Cargando paquete requerido: Hmisc

## Warning: package 'Hmisc' was built under R version 4.4.3

## 
## Adjuntando el paquete: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

2. Cargar el conjunto de datos lung (viene precargado en survival).

lung <- lung

3. Eliminar:

a.- Datos faltantes de variables que me interesa trabajar.

vis_miss()

→ Es una función que genera un gráfico de barras horizontales que muestra, para cada observación y cada variable, si hay o no datos faltantes.

Negro: datos faltantes (NA).
Gris claro: datos presentes.

vis_miss(lung[,c(1:5)])

vis_miss(lung[,c(6:10)])

lung <- lung %>% dplyr::filter(!is.na(ph.ecog))
lung <- lung %>% dplyr::filter(!is.na(ph.karno))
lung <- lung %>% dplyr::filter(!is.na(pat.karno))

b.- Columnas de las variables meal.cal y wt.loss porque tienen demasiados datos faltantes.

vis_miss(lung[,c(1:5)])

lung <- lung |> dplyr::select(-c(meal.cal,wt.loss,inst))
vis_miss(lung)

4. Convertimos a factor con etiquetas:

sex

lung$sex <- factor(lung$sex, levels = c(1, 2), labels = c("Hombre", "Mujer"))
levels(lung$sex)

## [1] "Hombre" "Mujer"

ph.ecog

lung$ph.ecog <- factor(lung$ph.ecog, levels = c(0, 1, 2, 3), labels = c("Normal", "Limitado", "Restriccion Moderada", "Restriccion Total"))
levels(lung$ph.ecog)

## [1] "Normal"               "Limitado"             "Restriccion Moderada"
## [4] "Restriccion Total"

5. Crear objeto con Surv: variable de respuesta “Tiempo al Evento”: tiempo_evento (time, status).

tiempo_evento <- Surv(lung$time,lung$status)

6. Selección de variables.

Vamos a usar un método automático asumiendo que se trata de un modelo PREDICTIVO donde buscamos los mejores predictores sin necesidad de que sean variables explicativas o un control por confundidores clásico.

Ajustar modelo full.

cox_model_full <- coxph(tiempo_evento ~ age+sex+ph.ecog+ph.karno+pat.karno, data = lung)
options(scipen = 999)
summary(cox_model_full)

## Call:
## coxph(formula = tiempo_evento ~ age + sex + ph.ecog + ph.karno + 
##     pat.karno, data = lung)
## 
##   n= 223, number of events= 160 
## 
##                                  coef exp(coef)  se(coef)      z Pr(>|z|)   
## age                          0.011550  1.011617  0.009565  1.208  0.22722   
## sexMujer                    -0.557191  0.572816  0.171146 -3.256  0.00113 **
## ph.ecogLimitado              0.598343  1.819103  0.239681  2.496  0.01255 * 
## ph.ecogRestriccion Moderada  1.076914  2.935606  0.380776  2.828  0.00468 **
## ph.ecogRestriccion Total     2.367312 10.668674  1.085446  2.181  0.02919 * 
## ph.karno                     0.015599  1.015721  0.009842  1.585  0.11297   
## pat.karno                   -0.011234  0.988828  0.007227 -1.555  0.12005   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## age                            1.0116    0.98852    0.9928    1.0308
## sexMujer                       0.5728    1.74576    0.4096    0.8011
## ph.ecogLimitado                1.8191    0.54972    1.1372    2.9099
## ph.ecogRestriccion Moderada    2.9356    0.34065    1.3918    6.1918
## ph.ecogRestriccion Total      10.6687    0.09373    1.2711   89.5454
## ph.karno                       1.0157    0.98452    0.9963    1.0355
## pat.karno                      0.9888    1.01130    0.9749    1.0029
## 
## Concordance= 0.648  (se = 0.025 )
## Likelihood ratio test= 33.41  on 7 df,   p=0.00002
## Wald test            = 33.97  on 7 df,   p=0.00002
## Score (logrank) test = 36.23  on 7 df,   p=0.000007

Aplicamos seleccion automatica mediante stepAIC. Esta función está en el paquete MASS. Admite en sentido forward (desde modelo vacío), backward (desde modelo completo) y both (ambas direcciones)

Selección de variables para modelo final.

👉 El algoritmo elige eliminar la variable que baje más el AIC, en este caso age, porque pasa de 1426.5 a 1425.9.

cox_model <- stepAIC(cox_model_full, direction = "both")

## Start:  AIC=1426.45
## tiempo_evento ~ age + sex + ph.ecog + ph.karno + pat.karno
## 
##             Df    AIC
## - age        1 1425.9
## <none>         1426.5
## - pat.karno  1 1426.8
## - ph.karno   1 1427.1
## - ph.ecog    3 1430.3
## - sex        1 1435.6
## 
## Step:  AIC=1425.93
## tiempo_evento ~ sex + ph.ecog + ph.karno + pat.karno
## 
##             Df    AIC
## <none>         1425.9
## - ph.karno   1 1426.0
## - pat.karno  1 1426.3
## + age        1 1426.5
## - ph.ecog    3 1429.7
## - sex        1 1434.9

summary(cox_model)

## Call:
## coxph(formula = tiempo_evento ~ sex + ph.ecog + ph.karno + pat.karno, 
##     data = lung)
## 
##   n= 223, number of events= 160 
## 
##                                  coef exp(coef)  se(coef)      z Pr(>|z|)   
## sexMujer                    -0.552095  0.575743  0.170772 -3.233  0.00123 **
## ph.ecogLimitado              0.583002  1.791408  0.239930  2.430  0.01510 * 
## ph.ecogRestriccion Moderada  1.086390  2.963555  0.384499  2.825  0.00472 **
## ph.ecogRestriccion Total     2.414081 11.179488  1.086390  2.222  0.02628 * 
## ph.karno                     0.014072  1.014171  0.009837  1.430  0.15258   
## pat.karno                   -0.011207  0.988855  0.007242 -1.548  0.12174   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## sexMujer                       0.5757    1.73689    0.4120    0.8046
## ph.ecogLimitado                1.7914    0.55822    1.1193    2.8670
## ph.ecogRestriccion Moderada    2.9636    0.33743    1.3948    6.2965
## ph.ecogRestriccion Total      11.1795    0.08945    1.3295   94.0067
## ph.karno                       1.0142    0.98603    0.9948    1.0339
## pat.karno                      0.9889    1.01127    0.9749    1.0030
## 
## Concordance= 0.651  (se = 0.025 )
## Likelihood ratio test= 31.93  on 6 df,   p=0.00002
## Wald test            = 33.06  on 6 df,   p=0.00001
## Score (logrank) test = 35.16  on 6 df,   p=0.000004

Vemos que la selección automática dejó 2 variables que no tienen Wald test significativo. Pero podemos probar que si las retiráramos índice C aumentaría. Por eso la función stepAIC las dejó. Porque están aportando información.

🔹 Relación entre Concordance e índice C:

El índice C (o índice de concordancia) mide la capacidad discriminativa del modelo: qué tan bien predice el orden de ocurrencia de los eventos.
En la salida de summary(cox_model) de R, el valor reportado como “Concordance” es justamente ese índice C.
El SE (error estándar) indica la precisión de esa estimación.

🔹 Interpretación rápida:

0.5 = sin poder discriminativo (equivalente a azar).
1.0 = predicción perfecta.
En tu caso: 0.651 → el modelo discrimina mejor que el azar, pero todavía lejos de perfecto.

🔹 Wald test en la salida de summary(cox_model):

Es una prueba global para el modelo ajustado.
Contrasta la H0 (ninguna covariable tiene efecto sobre el riesgo) con H1 (al menos una covariable tiene un coeficiente distinto de cero).

📊 Interpretación en tu caso:

Valor del estadístico: 33.06
Grados de libertad (df): 6 (todas las variables que quedaron en el modelo).
p-valor: 0.00001 → Rechaza la hipótesis nula global, es decir, el modelo contiene al menos una variable significativa.

💡 Importante:

Este Wald test global no es lo mismo que los p-valores individuales que aparecen al lado de cada variable.
Los p-valores individuales son tests de Wald para cada coeficiente por separado.
Podés tener un Wald test global muy significativo (el modelo sirve) aunque algunas variables individuales no lo sean (como ocurre aquí).

7. Chequeo de asunciones del modelo final.

Residuos de Schoenfeld

Los residuos de Schoenfeld son una herramienta central para evaluar el cumplimiento del supuesto de riesgos proporcionales en el modelo de Cox. Se calculan solo en los tiempos en que ocurre un evento (no en censuras), y hay uno por evento y por covariable.

La evaluación se realiza mediante un test de hipótesis + visualización de los residuos.

El test de Schoenfeld evalúa si el efecto de una covariable (β) cambia o no con el tiempo. Es decir si la relación entre esa covariable y el riesgo es constante en el tiempo (lo que nos exige el modelo de Cox), o si varía a lo largo del seguimiento.

Este test admite evaluar variables numéricas (continuas y discretas) y categóricas (binomiales, multinomiales y ordinales).

Como es un test de hipótesis tiene H₀ y H₁:

H₀ (hipótesis nula):

El efecto de la covariable no cambia con el tiempo → el riesgo es proporcional

H₁ (hipótesis alternativa):

El efecto de la covariable cambia con el tiempo → el riesgo no es proporcional

Cómo interpretamos el test entonces?

Si el p-valor es > o = 0.05, no se rechaza H₀
Si el p-valor es < 0.05, se rechaza H₀

El test nos da 2 informaciones:

Evaluación por covariable, evalua una a una en contexto multiple.

Evaluación global, evalua si el conjunto de covariables viola el supuesto.

Visualización:

Una por cada variable en contexto múltiple.

Test de proporcionalidad de riesgos de Schoenfeld.

test_ph <- cox.zph(cox_model)
test_ph

##           chisq df     p
## sex        1.77  1 0.184
## ph.ecog    6.03  3 0.110
## ph.karno   5.34  1 0.021
## pat.karno  4.58  1 0.032
## GLOBAL    11.74  6 0.068

📌 Evaluación individual:

sex → p = 0.184 ✅ (no viola el supuesto)
ph.ecog → p = 0.110 ✅ (no viola el supuesto)
ph.karno → p = 0.021 ❌ (posible violación)
pat.karno → p = 0.032 ❌ (posible violación)

📌 Evaluación global:

GLOBAL: p = 0.068 → ✅ No se rechaza H₀ → el modelo en conjunto no muestra violación significativa del supuesto de riesgos proporcionales.

Gráficos por variables de residuos de Schoenfeld en el tiempo.

Si el supuesto de Proporcionalidad del Hazard se cumple, estos residuos deberían fluctuar aleatoriamente alrededor de cero a lo largo del tiempo.

ggcoxzph(test_ph)

En el test individual de Schoenfeld, para que el supuesto de proporcionalidad del hazard se cumpla, la distribución visual de los residuos para cada covariable debería verse así:

Los puntos rojos (residuos estandarizados β(t) a lo largo del tiempo) deberían dispersarse de forma aleatoria alrededor de la línea horizontal en cero.
Sin tendencias claras: no debería apreciarse una pendiente marcada hacia arriba o hacia abajo.
La línea negra ajustada (tendencia suavizada) debería ser prácticamente horizontal y cercana a cero.
Las bandas punteadas (intervalos de confianza) deberían incluir el cero durante la mayor parte del rango de tiempo.
No deberían existir patrones sistemáticos ni oscilaciones amplias que sugieran que el efecto de la covariable cambia con el tiempo.

En resumen:

🔹 Visualmente, buscamos ruido aleatorio alrededor de cero.

🔹 Si la curva suavizada se aleja mucho de cero o muestra tendencia marcada, eso sugiere violación del supuesto.

En tu gráfico, aunque ph.karno y pat.karno dieron p < 0.05, la tendencia visual es leve, lo que coincide con la nota del PEC: “Si bien hay dos variables que individualmente no cumplen supuesto en el test de Schoenfeld, el gráfico no da señales de violación grosera.”

Gráficos Log-log.

Estos gráficos se utilizan cada vez menos ya que solo se puede evaluar variables categóricas.
Si se dibujan rectas más o menos paralelas sin entrecruzamientos groseros se asume que el supuesto se cumple.

Gráficos log(-log) para la variable sex.

ggsurvplot(survfit(tiempo_evento ~ sex, data = lung), fun = "cloglog",
           xlab = "Tiempo (dias)",
           ylab = "Log(-Log(S(t)))",
           palette = "simpsons")

Gráficos log(-log) para la variable ph.ecog.

ggsurvplot(survfit(tiempo_evento ~ ph.ecog, data = lung), fun = "cloglog",
           xlab = "Tiempo (dias)",
           ylab = "Log(-Log(S(t)))",
           palette = "simpsons")

Respecto a los gráficos log-log:

Se espera que las curvas sean más o menos paralelas y que no haya entrecruzamientos grandes.
El primer gráfico (sex) cumple: las líneas son casi paralelas y no se cruzan de forma grosera.
El segundo (ph.ecog) muestra un cruce claro entre categorías → eso indica que para esta variable el supuesto podría estar comprometido, lo que coincide con tu resultado del test de Schoenfeld (p global ≈ 0,11 para esa variable).

Residuos de Martingale y Deviance.

Básicamente lo que miden estos residuos es la diferencia entre lo observado (si ocurrió el evento) y lo esperado por el modelo. Evaluan variables continuas y categóricas como los de Schoenfeld.

Evaluan si la relación entre covariable y log hazard es lineal. Pueden identificar necesidad de hacer transformaciones o términos no lineales para mejorar el ajuste:

Si los puntos se dispersan de forma aleatoria alrededor de la línea cero → buena especificación funcional.
Si hay patrones curvos o sistemáticos → posible comportamiento no lineal del modelo para la variable que se está evaluado.

Dado que los residuos de Martingale son asimétricos, se los puede transformar en residuos de Deviance, que son más comparables a los residuos studentizados que usamos en otros modelos de regresión. Son simétricos y con media cero.

Al tener un aspecto más “estandarizado” también son más fáciles de usar para identificar observaciones influyentes.

Patrones curvos o con tendencia señalan mala especificación (no lineal) como los de Martingale mientras que puntos extremos alejados del resto pueden indicar individuos con influencia desproporcionada sobre el modelo.

Gráficos de residuos.

a.- Martingale

ggcoxdiagnostics(cox_model, type = "martingale", linear.predictions = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

Martingale residuals:

Se usan principalmente para evaluar la forma funcional de covariables continuas.
Si el modelo está bien especificado, la línea suavizada (azul) debería oscilar alrededor de la línea roja en cero sin mostrar patrones sistemáticos.
En tu gráfico, la curva azul se mantiene relativamente cerca de cero, aunque hay ligera curvatura al inicio: podría sugerir que alguna variable continua no está perfectamente modelada de forma lineal.

b.- Deviance.

ggcoxdiagnostics(cox_model, type = "deviance", linear.predictions = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

Deviance residuals:

Son más simétricos y útiles para detectar outliers o valores atípicos que no se ajustan bien.
Valores grandes en valor absoluto (muy alejados de cero, arriba o abajo) indican observaciones con mal ajuste.
En tu gráfico, hay algunos puntos que se alejan bastante, pero no parece haber un patrón fuerte sistemático.
Lo importante es revisar si hay observaciones con residuales extremos (ej. > ±3) para evaluar si son influyentes.

8. Evaluación de Discriminación y Calibración del modelo.

Discriminación.

Indice C de Harrell.

summary(cox_model)$concordance

##          C      se(C) 
## 0.65099677 0.02516653

✅ ¿Qué significa?

El índice C = 0.65 indica que el modelo tiene una capacidad de discriminación del 65%, es decir:

En el 65% de las veces, el modelo predice correctamente que un paciente con peor pronóstico (mayor riesgo) tendrá el evento antes que otro con mejor pronóstico.

🟡 Valores típicos y su interpretación:

0.5 = Sin discriminación (como tirar una moneda)

0.6–0.7 = Discriminación pobre/modesta

0.7–0.8 = Buena discriminación

mayor a 0.8 = Excelente discriminación

Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete rms.

dd <- datadist(lung); options(datadist = "dd")
cox_metrics <- cph(tiempo_evento ~ age + sex + ph.ecog + ph.karno + pat.karno, 
                   data = lung, x = TRUE, y = TRUE, surv = TRUE)

Validamos métricas con bootstrap.

invisible(capture.output({ #esta linea es para evitar que se imprima el output del bootstrap)
  cox_rms <- validate(cox_metrics, method = "boot", B = 200)
}))

## Warning in fitter(..., strata = strata, rownames = rownames, offset = offset, :
## Loglik converged before variable 1 ; coefficient may be infinite.

⚙️ Se están generando 200 muestras bootstrap (con reemplazo) para validar internamente el modelo de Cox y calcular métricas corregidas por optimismo como:

Dxy (Índice de Somers)
R²
Slope
D, U, Q
g (índice de discriminación)

🧠 ¿Por qué corregir por optimismo?

Cuando uno ajusta un modelo y lo evalúa sobre la misma base de datos, es común que las métricas de desempeño estén sobreestimadas.

🔁 El bootstrap:

Toma una muestra aleatoria con reemplazo (bootstrap sample).
Ajusta el modelo sobre esa muestra.
Lo evalúa tanto en esa muestra (entrenamiento) como en la base original (test).
Calcula la diferencia (optimismo).
Repite este proceso B = 200 veces.
Calcula un promedio del optimismo y lo resta al desempeño original → obteniendo el index.corrected.

📊 En resumen:

Acá el bootstrap se usa para validar internamente el modelo de Cox y corregir las métricas por sobreajuste.
Es el mismo concepto que en otros modelos, pero aplicado a métricas propias de supervivencia.

Métricas de discriminación y calibración (validadas con bootstrap).

Revisamos las metricas validadas (en la columna index.corrected).

cox_rms

##       index.orig training   test optimism index.corrected   n
## Dxy       0.2959   0.3221 0.2824   0.0397          0.2562 142
## R2        0.1393   0.1666 0.1205   0.0461          0.0932 142
## Slope     1.0000   1.0000 0.8391   0.1609          0.8391 142
## D         0.0224   0.0278 0.0191   0.0087          0.0137 142
## U        -0.0014  -0.0014 0.0021  -0.0035          0.0021 142
## Q         0.0238   0.0292 0.0170   0.0122          0.0116 142
## g         0.5667   0.6296 0.5118   0.1178          0.4489 142

Para visualizar más facil podemos pedir el output solo la columna que nos interesa y redondear a 3 decimales.

round(cox_rms[,5],3)

##   Dxy    R2 Slope     D     U     Q     g 
## 0.256 0.093 0.839 0.014 0.002 0.012 0.449

¿Como lo interpretamos?

Dxy es el Indice de Somers. El valor es positivo pero relativamente bajo como el índice de Harrel.
R2 es el pseudoR2. El valor es muy bajo.
Slope es la pendiente de calibración. Menor a 1 sobrestima cuando el riesgo es bajo y subestima cuando el riesgo es alto.
D U Q g es la descomposicion de información. El valor de D cercano a cero señala que el modelo aporta pero no discrimina bien, si bien el g es aceptable como discriminacion promedio.

En el contexto de la validación con bootstrap (como en el output de validate.cph del paquete rms):

La columna optimism mide la diferencia entre lo que el modelo logra en el conjunto de entrenamiento y lo que se espera en validación repetida con bootstrap.

Es, en otras palabras, una estimación del sobreajuste.

📌 Interpretación práctica:

Optimism ≈ 0 → el modelo casi no está sobreajustado, generaliza bien.
Optimism pequeño (ej. <0.02–0.03) → indica un buen comportamiento (poco sobreajuste).
Optimism más alto (>0.05) → alerta 🚨, significa que el modelo rinde bastante mejor en entrenamiento que en validación, y puede no generalizar bien.

Grafiquemos la curva de calibración.

Calibración. Evaluación Gráfica.

invisible(capture.output({
cal <- calibrate(cox_metrics, method = "boot", B = 200, u = 120)
}))
plot(cal, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Gráfico de Calibración",
     cex.main = 1.2,
     cex.axis = 0.7,
     cex.lab = 0.9,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")

📈 ¿Qué muestra el gráfico?

Este gráfico de calibración compara la probabilidad de sobrevida predicha por el modelo (eje X) con la sobrevida observada (eje Y) a los 120 días. Sirve para evaluar cuán bien calibrado está el modelo de Cox: es decir, si las probabilidades predichas coinciden con lo que realmente ocurrió.

📏 Líneas importantes:

Línea negra punteada: es la línea de identidad (perfecta calibración). Cada punto sobre esta línea implica que lo observado coincide exactamente con lo predicho.
Línea azul: es la curva de calibración generada por el modelo a partir de validación con bootstrap.

✅ Interpretación práctica:

En general, la curva azul sigue bastante bien la línea negra punteada, lo que indica que el modelo tiene buena calibración.
Hay una leve subestimación de la sobrevida observada en algunos tramos (la línea azul está por encima de la diagonal), lo que sugiere que el modelo predice una probabilidad algo más baja de lo que realmente se observó.
Como está basado en validación bootstrap (B = 200), es más confiable porque reduce el optimismo del modelo entrenado sobre la misma base.

📌 Conclusión:

El modelo de Cox calibrado muestra un ajuste adecuado entre la sobrevida predicha y la observada a los 120 días, con ligera subestimación en los valores altos de sobrevida, pero en general es un modelo bien calibrado.

Resumen explicativo:

En esta primera parte, trabajamos con un modelo de Cox para analizar datos de supervivencia.

El flujo general fue:

Partimos de un modelo inicial (cox_model_full) con un conjunto de variables seleccionadas según criterio clínico o de interés.
Aplicamos selección automática de variables usando stepAIC() para quedarnos con las que aportan más al modelo según el criterio de información AIC.

Este paso puede eliminar variables no relevantes (por ejemplo, se quitó age en este caso).
El modelo reducido resultante se llama cox_model.

Evaluamos la capacidad de discriminación del modelo:

Índice C de Harrell: nos dice con qué frecuencia el modelo predice correctamente quién tendrá antes el evento.
- 0.5 es azar,
- 0.65–0.7 es discriminación modesta,
- mayor a 0.8 es excelente.
Índice D de Somers (Dxy) y g: miden también la discriminación (qué tan bien el modelo separa alto de bajo riesgo).

Evaluamos la calibración (qué tan bien las predicciones coinciden con la realidad):

Slope: pendiente ideal ≈ 1. Valores menores indican sobreestimación/subestimación de riesgo.
R²: proporción de variabilidad explicada (pseudo R² en supervivencia).
D, U, Q: medidas de descomposición de la información (valores cercanos a 0 indican poco aporte a la discriminación).

Validamos el modelo con bootstrap:

Repetimos el ajuste muchas veces en muestras simuladas para estimar el optimismo (sobreajuste).
Restamos ese optimismo del desempeño original → obtenemos métricas corrigidas (columna index.corrected).

🎰 BONUS!

Explorando las variables problema del dataset, esas que no son significativas en el Wald Test pero quedaron por mejorar el AIC.

cox_ph.karno <- coxph(tiempo_evento ~ ph.karno, data = lung)
lung$deviance.ph <- resid(cox_ph.karno, type = "deviance")
ggplot(lung, aes(x = ph.karno, y = deviance.ph)) +
  geom_point() +
  geom_smooth(method = "loess", se = TRUE, col = "blue", linetype = "dashed") +
  geom_hline(yintercept = 0, col = "red") +
  labs(title = "Residuos de Deviance vs ph.karno",
       y = "Residuo Deviance", x = "ph.karno")

## `geom_smooth()` using formula = 'y ~ x'

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 100.25

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 20.25

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 8.0122e-16

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 100

## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 100.25

## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 20.25

## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 8.0122e-16

## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 100

Contexto del análisis

Se está evaluando la variable ph.karno en un modelo de Cox.
Aunque esta variable no es significativa en el Wald test, el algoritmo de selección automática la mantuvo porque ayudaba a mejorar el AIC.

Qué muestra el gráfico

Eje x: valores de ph.karno (escala de estado funcional del paciente).
Eje y: residuo deviance (mide la discrepancia entre lo observado y lo esperado por el modelo para cada observación).
Línea roja: nivel de referencia en cero (residuo esperado si el modelo ajusta bien).
Línea azul discontinua: tendencia suavizada de los residuos (método loess).
Banda gris: intervalo de confianza de la tendencia.

Hallazgos clave:

La línea azul se mantiene cerca de cero, pero con cierta curvatura, sobre todo en valores bajos de ph.karno (50–70), lo que podría indicar una relación no estrictamente lineal con el riesgo.
Se observa que ph.karno no es verdaderamente continua en esta base: los valores aparecen “escalonados” en múltiplos de 10, lo que se denomina discretización.
Esto puede hacer que el modelo trate mal la variable como continua, afectando su ajuste.

Posibles soluciones planteadas:

Opción 1: Categorizar ph.karno en tramos (p. ej. ≤60, 70, 80, etc.) para que el modelo no asuma linealidad donde no la hay.
Opción 2: Verificar colinealidad con pat.karno (escala reportada por el médico), ya que podrían estar muy correlacionadas.

En pocas palabras:

ph.karno quedó en el modelo por mejorar el AIC, pero los residuos sugieren que su relación con la supervivencia podría no ser lineal y que la medición en saltos de 10 puntos la hace más apta para análisis como variable categórica o para evaluar redundancia con otra escala similar.

cox_pat.karno <- coxph(Surv(time, status) ~ pat.karno, data = lung)
lung$deviance.pat <- resid(cox_ph.karno, type = "deviance")
ggplot(lung, aes(x = pat.karno, y = deviance.pat)) +
  geom_point() +
  geom_smooth(method = "loess", se = TRUE, col = "blue", linetype = "dashed") +
  geom_hline(yintercept = 0, col = "red") +
  labs(title = "Residuos de Deviance vs pat.karno",
       y = "Residuo Deviance", x = "pat.karno")

## `geom_smooth()` using formula = 'y ~ x'

Interpretación visual:

La línea azul se mantiene bastante cerca de cero, pero presenta una leve curvatura, sobre todo en valores bajos (<60), sugiriendo que la relación con el riesgo podría no ser perfectamente lineal.
Igual que con ph.karno, los valores están “escalonados” en múltiplos de 10, lo que indica que la variable está discretizada aunque se almacene como numérica continua.
Esto podría estar introduciendo pequeños problemas de ajuste si el modelo la trata como verdaderamente continua.

Conclusión práctica:

Tanto ph.karno como pat.karno presentan patrones residuales similares y escalas discretizadas.
Es probable que exista colinealidad entre ambas, ya que ambas miden el estado funcional (una autorreportada y la otra por el médico).
Antes de decidir mantener ambas, conviene evaluar la redundancia y, si es necesario, quedarse con una o bien categorizarlas.

Vemos que si bien son escalas continuas, ambas están “discretizadas” ya que no hay valores intermedios entre decenas.

1️⃣ Una opción es categorizarlas en forma ordinal y juntas las categorías más bajas (menor o igual a 60) para tener todas con una cantidad suficiente de observaciones.

2️⃣ Otra opción es evaluar su colinealidad ya que pareciera ser bastante coincidente lo que reportan medicos y pacientes en esta escala.

Categorizamos ph.karno y pat.karno

lung <- lung |> mutate(ph.karno.cat=as.factor(case_when(ph.karno>90~1,
                                           ph.karno<=90 & ph.karno>80~2,
                                           ph.karno<=80 & ph.karno>70~3,
                                           ph.karno<=70 & ph.karno>60~4,
                                           ph.karno<=60~5
                                           )))

lung <- lung |> mutate(pat.karno.cat=as.factor(case_when(pat.karno>90~1,
                                           pat.karno<=90 & pat.karno>80~2,
                                           pat.karno<=80 & pat.karno>70~3,
                                           pat.karno<=70 & pat.karno>60~4,
                                           pat.karno<=60~5
                                           )))

Reajustamos el modelo.

cox_model_cat <- coxph(tiempo_evento ~ sex+ph.ecog+ph.karno.cat+pat.karno.cat, data = lung)
summary(cox_model_cat)

## Call:
## coxph(formula = tiempo_evento ~ sex + ph.ecog + ph.karno.cat + 
##     pat.karno.cat, data = lung)
## 
##   n= 223, number of events= 160 
## 
##                                 coef exp(coef) se(coef)      z Pr(>|z|)   
## sexMujer                    -0.55055   0.57663  0.17260 -3.190  0.00142 **
## ph.ecogLimitado              0.46241   1.58790  0.28566  1.619  0.10551   
## ph.ecogRestriccion Moderada  0.89987   2.45928  0.42717  2.107  0.03515 * 
## ph.ecogRestriccion Total     2.40453  11.07326  1.10894  2.168  0.03013 * 
## ph.karno.cat2               -0.01466   0.98545  0.31893 -0.046  0.96334   
## ph.karno.cat3               -0.01839   0.98178  0.39165 -0.047  0.96254   
## ph.karno.cat4               -0.18070   0.83468  0.45440 -0.398  0.69088   
## ph.karno.cat5               -0.33976   0.71194  0.48586 -0.699  0.48437   
## pat.karno.cat2               0.08940   1.09352  0.27849  0.321  0.74820   
## pat.karno.cat3               0.12092   1.12854  0.27955  0.433  0.66533   
## pat.karno.cat4               0.14473   1.15573  0.31723  0.456  0.64822   
## pat.karno.cat5               0.52789   1.69535  0.34448  1.532  0.12542   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## sexMujer                       0.5766    1.73420    0.4111    0.8088
## ph.ecogLimitado                1.5879    0.62976    0.9071    2.7796
## ph.ecogRestriccion Moderada    2.4593    0.40662    1.0646    5.6808
## ph.ecogRestriccion Total      11.0733    0.09031    1.2599   97.3209
## ph.karno.cat2                  0.9855    1.01476    0.5274    1.8412
## ph.karno.cat3                  0.9818    1.01856    0.4557    2.1154
## ph.karno.cat4                  0.8347    1.19806    0.3426    2.0338
## ph.karno.cat5                  0.7119    1.40461    0.2747    1.8450
## pat.karno.cat2                 1.0935    0.91448    0.6335    1.8875
## pat.karno.cat3                 1.1285    0.88610    0.6525    1.9519
## pat.karno.cat4                 1.1557    0.86526    0.6206    2.1522
## pat.karno.cat5                 1.6954    0.58985    0.8630    3.3303
## 
## Concordance= 0.662  (se = 0.025 )
## Likelihood ratio test= 32.73  on 12 df,   p=0.001
## Wald test            = 35.11  on 12 df,   p=0.0004
## Score (logrank) test = 37.63  on 12 df,   p=0.0002

test_cat <- cox.zph(cox_model_cat)
test_cat

##               chisq df     p
## sex            2.08  1 0.149
## ph.ecog        6.34  3 0.096
## ph.karno.cat   6.41  4 0.171
## pat.karno.cat  6.16  4 0.188
## GLOBAL        15.09 12 0.237

ggcoxzph(test_cat)

Como vemos, ahora se cumplen supuestos para todas las variables

¿Pero el rendimiento es mejor?

summary(cox_model)$concordance

##          C      se(C) 
## 0.65099677 0.02516653

summary(cox_model_cat)$concordance

##          C      se(C) 
## 0.66228650 0.02490197

La diferencia en C de Harrell es mínima pero mayor en el modelo nuevo.

Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete rms

cox_metrics_cat <- cph(tiempo_evento ~ sex + ph.ecog + ph.karno.cat + pat.karno.cat, data = lung, x = TRUE, y = TRUE, surv = TRUE)

Validacion con bootstrap de las metricas calculadas.

invisible(capture.output({
cox_rms_cat <-validate(cox_metrics_cat, method = "boot", B = 200) 
}))

## Warning in fitter(..., strata = strata, rownames = rownames, offset = offset, :
## Loglik converged before variable 5,11 ; coefficient may be infinite.

Revisamos las metricas validadas de ambos modelos.

round(cox_rms[,5],3)

##   Dxy    R2 Slope     D     U     Q     g 
## 0.256 0.093 0.839 0.014 0.002 0.012 0.449

round(cox_rms_cat[,5],3)

##    Dxy     R2  Slope      D      U      Q      g 
##  0.245  0.046  0.681  0.005  0.006 -0.001  0.293

🔥 La métricas de discriminación empeoraron un poco excepto la Dxy y la C de Harrel… Evaluemos gráficamente la calibración aunque ya vimos la pendiente.

Calibración. Evaluación Gráfica.

invisible(capture.output({
cal_cat <- calibrate(cox_metrics_cat, method = "boot", B = 200, u = 120)
}))

## Warning in fitter(X, Y, strata = Strata, offset = offset, weights = weights, :
## Loglik converged before variable 7,8 ; coefficient may be infinite.

plot(cal_cat, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Gráfico de Calibración",
     cex.main = 1.2,
     cex.axis = 0.7,
     cex.lab = 0.9,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")

Comparemos gráficos.

par(mfrow=c(1,2))
    
plot(cal, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Modelo Cont",
     cex.main = 1,
     cex.axis = 0.5,
     cex.lab = 0.5,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")
plot(cal_cat, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Modelo Cat",
     cex.main = 1,
     cex.axis = 0.5,
     cex.lab = 0.5,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")

🔥 La calibración del segundo modelo es peor respecto al primero.

Veamos residuos.

Gráficos de residuos

ggcoxdiagnostics(cox_model_cat, type = "martingale", linear.predictions = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

ggcoxdiagnostics(cox_model_cat, type = "deviance", linear.predictions = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

✅ Los residuos no muestran patrones importantes

Entonces… ¿A qué se deben la pérdida de calibración y capacidad discriminativa con este modelo?

Básicamente a la pérdida de granularidad producto de la categorización de las variables ph.karno y pat.karno.

Tambien es posible que la categorización haya afectado la colinealidad de variables. Veamos…

vif(cox_model_cat)

##                    sexMujer             ph.ecogLimitado 
##                    1.037053                    3.241072 
## ph.ecogRestriccion Moderada    ph.ecogRestriccion Total 
##                    5.534007                    1.193442 
##               ph.karno.cat2               ph.karno.cat3 
##                    3.406899                    4.969823 
##               ph.karno.cat4               ph.karno.cat5 
##                    4.709075                    3.913406 
##              pat.karno.cat2              pat.karno.cat3 
##                    2.215041                    2.283554 
##              pat.karno.cat4              pat.karno.cat5 
##                    2.366891                    3.058447

vif(cox_model)

##                    sexMujer             ph.ecogLimitado 
##                    1.014802                    2.286385 
## ph.ecogRestriccion Moderada    ph.ecogRestriccion Total 
##                    4.482448                    1.145403 
##                    ph.karno                   pat.karno 
##                    2.543566                    1.678548

🚩 Tenemos un problema:

El modelo con variables categorizadas cumple todos los supuestos y el modelo continuo no.
El modelo categorizado pierde calibración y capacidad discriminativa y suma problemas de moderada multicolinearidad (VIF >5).

Podríamos probar retirar una de las variables karno y comparar modelos…

Preservamos la escala de pacientes.

cox_model_cat2 <- coxph(tiempo_evento ~ sex+ph.ecog+pat.karno.cat, data = lung)
summary(cox_model_cat2)

## Call:
## coxph(formula = tiempo_evento ~ sex + ph.ecog + pat.karno.cat, 
##     data = lung)
## 
##   n= 223, number of events= 160 
## 
##                                 coef exp(coef) se(coef)      z Pr(>|z|)   
## sexMujer                    -0.54094   0.58220  0.17221 -3.141  0.00168 **
## ph.ecogLimitado              0.42175   1.52462  0.21268  1.983  0.04736 * 
## ph.ecogRestriccion Moderada  0.67449   1.96303  0.29584  2.280  0.02261 * 
## ph.ecogRestriccion Total     2.09360   8.11407  1.04978  1.994  0.04612 * 
## pat.karno.cat2               0.07634   1.07933  0.27693  0.276  0.78281   
## pat.karno.cat3               0.07341   1.07617  0.27513  0.267  0.78962   
## pat.karno.cat4               0.10038   1.10560  0.30598  0.328  0.74286   
## pat.karno.cat5               0.55324   1.73888  0.33195  1.667  0.09559 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## sexMujer                       0.5822     1.7176    0.4154    0.8159
## ph.ecogLimitado                1.5246     0.6559    1.0049    2.3131
## ph.ecogRestriccion Moderada    1.9630     0.5094    1.0993    3.5055
## ph.ecogRestriccion Total       8.1141     0.1232    1.0367   63.5057
## pat.karno.cat2                 1.0793     0.9265    0.6272    1.8573
## pat.karno.cat3                 1.0762     0.9292    0.6276    1.8453
## pat.karno.cat4                 1.1056     0.9045    0.6069    2.0140
## pat.karno.cat5                 1.7389     0.5751    0.9072    3.3329
## 
## Concordance= 0.663  (se = 0.025 )
## Likelihood ratio test= 31.8  on 8 df,   p=0.0001
## Wald test            = 33.74  on 8 df,   p=0.00005
## Score (logrank) test = 36.15  on 8 df,   p=0.00002

Preservamos la escala de médicos (ph=phisicians)

cox_model_cat3 <- coxph(tiempo_evento ~ sex+ph.ecog+ph.karno.cat, data = lung)
summary(cox_model_cat3)

## Call:
## coxph(formula = tiempo_evento ~ sex + ph.ecog + ph.karno.cat, 
##     data = lung)
## 
##   n= 223, number of events= 160 
## 
##                                  coef exp(coef)  se(coef)      z Pr(>|z|)    
## sexMujer                    -0.564114  0.568864  0.171365 -3.292 0.000995 ***
## ph.ecogLimitado              0.460496  1.584859  0.277192  1.661 0.096656 .  
## ph.ecogRestriccion Moderada  1.161863  3.195883  0.394961  2.942 0.003264 ** 
## ph.ecogRestriccion Total     2.518990 12.416056  1.101542  2.287 0.022208 *  
## ph.karno.cat2                0.006539  1.006560  0.314437  0.021 0.983410    
## ph.karno.cat3                0.061227  1.063140  0.376871  0.162 0.870943    
## ph.karno.cat4               -0.133046  0.875425  0.450984 -0.295 0.767984    
## ph.karno.cat5               -0.406575  0.665927  0.475499 -0.855 0.392524    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## sexMujer                       0.5689    1.75789    0.4066    0.7959
## ph.ecogLimitado                1.5849    0.63097    0.9205    2.7286
## ph.ecogRestriccion Moderada    3.1959    0.31290    1.4737    6.9308
## ph.ecogRestriccion Total      12.4161    0.08054    1.4333  107.5519
## ph.karno.cat2                  1.0066    0.99348    0.5435    1.8642
## ph.karno.cat3                  1.0631    0.94061    0.5079    2.2253
## ph.karno.cat4                  0.8754    1.14230    0.3617    2.1188
## ph.karno.cat5                  0.6659    1.50167    0.2622    1.6911
## 
## Concordance= 0.644  (se = 0.026 )
## Likelihood ratio test= 29.95  on 8 df,   p=0.0002
## Wald test            = 31.34  on 8 df,   p=0.0001
## Score (logrank) test = 33.46  on 8 df,   p=0.00005

Comparamos con modelo categórico base.

anova(cox_model_cat,cox_model_cat2)

## Analysis of Deviance Table
##  Cox model: response is  tiempo_evento
##  Model 1: ~ sex + ph.ecog + ph.karno.cat + pat.karno.cat
##  Model 2: ~ sex + ph.ecog + pat.karno.cat
##    loglik  Chisq Df Pr(>|Chi|)
## 1 -706.56                     
## 2 -707.03 0.9315  4       0.92

AIC(cox_model_cat)

## [1] 1437.125

AIC(cox_model_cat2)

## [1] 1430.057

AIC(cox_model_cat3)

## [1] 1431.91

Si bien las diferencias no son significativas el modelo cox_model_cat2 (que preserva pat.karno_cat) es superior en AIC al que mantiene ambas variables categoricas y al que preserva solo ph.karno_cat.

👀 Veamos supuestos y rendimiento de este modelo:

test_cat2 <- cox.zph(cox_model_cat2)
test_cat2

##               chisq df    p
## sex            2.36  1 0.12
## ph.ecog        6.15  3 0.10
## pat.karno.cat  5.87  4 0.21
## GLOBAL        10.97  8 0.20

🎊 El modelo cumple supuestos.

¿Multicolinearidad?

vif(cox_model_cat2)

##                    sexMujer             ph.ecogLimitado 
##                    1.033582                    1.798675 
## ph.ecogRestriccion Moderada    ph.ecogRestriccion Total 
##                    2.653884                    1.069492 
##              pat.karno.cat2              pat.karno.cat3 
##                    2.193361                    2.210704 
##              pat.karno.cat4              pat.karno.cat5 
##                    2.202081                    2.834860

🙌 No hay evidencias de multicolinearidad moderada o severa.

Métricas de rendimiento.

summary(cox_model_cat2)$concordance

##          C      se(C) 
## 0.66263022 0.02481438

Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete rms.

cox_metrics_cat2 <- cph(tiempo_evento ~ sex + ph.ecog +pat.karno.cat, data = lung, x = TRUE, y = TRUE, surv = TRUE)

Validación bootstrap.

invisible(capture.output({
cox_rms_cat2 <-validate(cox_metrics_cat2, method = "boot", B = 200) # Validacion con bootstrap de las metricas calculadas
}))

Revisamos las metricas validadas y C de Harrel.

summary(cox_model_cat2)$concordance

##          C      se(C) 
## 0.66263022 0.02481438

round(cox_rms_cat2[,5],3)

##   Dxy    R2 Slope     D     U     Q     g 
## 0.269 0.072 0.802 0.010 0.003 0.008 0.382

🎉 Mejores métricas respecto al modelo categórico base e incluso una C de Harrel y una D de Somers algo mejores que en el modelo continuo

Vemos el grafico de calibracion.

Calibración. Evaluación Gráfica.

invisible(capture.output({
cal_cat2 <- calibrate(cox_metrics_cat2, method = "boot", B = 200, u = 120)
}))

plot(cal_cat2, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Gráfico de Calibración",
     cex.main = 1.2,
     cex.axis = 0.7,
     cex.lab = 0.9,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")

Comparemos con la calibración del modelo que preserva las dos variables categóricas.

par(mfrow=c(1,2))
    
  plot(cal_cat, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Modelo Cat 2 variables",
     cex.main = 1,
     cex.axis = 0.5,
     cex.lab = 0.5,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")
plot(cal_cat2, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Modelo Cat 1 variable",
     cex.main = 1,
     cex.axis = 0.5,
     cex.lab = 0.5,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")

✅ La calibración es algo mejor que la del modelo categorico original.

¿Y respecto al modelo original?

par(mfrow=c(1,2))
    
  plot(cal, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Modelo Cont",
     cex.main = 1,
     cex.axis = 0.5,
     cex.lab = 0.5,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")
plot(cal_cat2, 
     subtitles = FALSE,
     lwd = 2, 
     lty = 2, 
     main = "Modelo Cat 1 variable",
     cex.main = 1,
     cex.axis = 0.5,
     cex.lab = 0.5,
     xlab = "Probabilidad Predicha de Sobrevida a 120 días",
     ylab = "Sobrevida Observada")

😎 Definitivamente este es nuestro mejor modelo, que sufre en calibración pero menos que el modelo con 2 variables, cumple supuestos que es el punto de partida y tiene mejor rendimiento que el modelo categórico base y que el modelo que preserva la escala ph en lugar de la escala pat.

Taller 2: Cox PH.

Ceci Cocucci + ChatGPT

2025-08-06