# Paquetes necesarios
library(survival) # Análisis de supervivencia
library(survminer) # Visualización
library(flexsurv) # Modelos paramétricos flexibles
library(KMsurv) # Datasets de supervivencia
library(ggplot2) # Gráficos
library(dplyr) # Manipulación de datos
library(knitr) # Tablas
library(rms) # Regresión y modelización
library(cmprsk) # Riesgos competitivos
# Configuración de gráficos
theme_set(theme_minimal())Tutorial de Análisis de Supervivencia en R
Métodos estadísticos para datos de tiempo hasta el evento
Resumen
Este tutorial presenta una introducción exhaustiva al análisis de supervivencia, cubriendo desde conceptos fundamentales hasta técnicas avanzadas de modelización. El análisis de supervivencia es una rama de la estadística que estudia el tiempo hasta que ocurre un evento de interés, considerando la presencia de datos censurados. Se presentan métodos no paramétricos (Kaplan-Meier), semiparamétricos (Cox) y paramétricos (Weibull, exponencial), con ejemplos prácticos implementados en R.
Palabras clave
análisis de supervivencia, survival analysis, Kaplan-Meier, modelo de Cox, censura, riesgo proporcional, modelos paramétricos, bioestadística, epidemiología, tiempo hasta el evento, riesgos competitivos, R
1 Introducción
1.1 ¿Qué es el análisis de supervivencia?
El análisis de supervivencia es un conjunto de métodos estadísticos diseñados para analizar datos donde la variable de interés es el tiempo hasta que ocurre un evento. Aunque el término proviene de estudios médicos donde el evento de interés es la muerte, estos métodos se aplican ampliamente en:
- Medicina: tiempo hasta recaída, muerte, recuperación
- Ingeniería: tiempo hasta fallo de un componente
- Economía: duración del desempleo, permanencia en un trabajo
- Ciencias sociales: duración de matrimonios, tiempo hasta reincidencia criminal
1.2 Características especiales de los datos de supervivencia
Los datos de supervivencia presentan dos características distintivas:
- Positividad: Los tiempos de supervivencia son siempre positivos (\(T > 0\))
- Censura: No siempre observamos el evento de interés para todos los sujetos
1.3 Censura
La censura es la característica más importante que distingue al análisis de supervivencia de otros métodos estadísticos. Un dato está censurado cuando no observamos el tiempo exacto del evento.
1.3.1 Tipos de censura
1.3.1.1 Censura por la derecha
El tipo más común. Ocurre cuando:
- El estudio termina antes de que ocurra el evento
- El sujeto se pierde durante el seguimiento
- El sujeto experimenta un evento competitivo
Notación: Si \(T_i\) es el verdadero tiempo de supervivencia y \(C_i\) el tiempo de censura, observamos:
\[ Y_i = \min(T_i, C_i) \]
y el indicador de censura:
\[ \delta_i = \begin{cases} 1 & \text{si } T_i \leq C_i \text{ (evento observado)} \\ 0 & \text{si } T_i > C_i \text{ (censurado)} \end{cases} \]
1.3.1.2 Censura por la izquierda
Ocurre cuando el evento ya ha ocurrido antes del inicio del estudio. Menos común en práctica.
1.3.1.3 Censura por intervalo
El evento ocurre en un intervalo conocido \((L_i, R_i]\) pero no se conoce el tiempo exacto.
1.4 Configuración del entorno
2 Conceptos fundamentales
2.1 Función de supervivencia
La función de supervivencia \(S(t)\) define la probabilidad de que un individuo sobreviva más allá del tiempo \(t\):
\[ S(t) = P(T > t) = 1 - F(t) \]
donde \(F(t) = P(T \leq t)\) es la función de distribución acumulada.
Propiedades:
- \(S(0) = 1\) (todos están vivos al inicio)
- \(S(\infty) = 0\) (eventualmente todos experimentan el evento)
- \(S(t)\) es una función decreciente y continua por la derecha
2.2 Función de densidad
La función de densidad de probabilidad se define como:
\[ f(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t)}{\Delta t} = \frac{dF(t)}{dt} = -\frac{dS(t)}{dt} \]
2.3 Función de riesgo (hazard)
La función de riesgo \(h(t)\) (o hazard function) representa la tasa instantánea de ocurrencia del evento en el tiempo \(t\), dado que el individuo ha sobrevivido hasta ese momento:
\[ h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t} = \frac{f(t)}{S(t)} = -\frac{d \log S(t)}{dt} \]
Interpretación: \(h(t)\) es la probabilidad condicional instantánea del evento.
2.4 Función de riesgo acumulado
La función de riesgo acumulado \(H(t)\) se define como:
\[ H(t) = \int_0^t h(u) \, du = -\log S(t) \]
De donde se deduce:
\[ S(t) = \exp\left(-H(t)\right) = \exp\left(-\int_0^t h(u) \, du\right) \]
2.5 Relaciones entre funciones
Las cuatro funciones (\(S(t)\), \(f(t)\), \(h(t)\), \(H(t)\)) están interrelacionadas:
\[ \begin{aligned} S(t) &= \exp\left(-\int_0^t h(u) \, du\right) \\ f(t) &= h(t) S(t) \\ h(t) &= \frac{f(t)}{S(t)} = -\frac{d \log S(t)}{dt} \end{aligned} \]
3 Métodos no paramétricos
3.1 Estimador de Kaplan-Meier
El estimador de Kaplan-Meier (1) es el método no paramétrico más utilizado para estimar la función de supervivencia en presencia de censura.
3.1.1 Construcción del estimador
Sean \(t_1 < t_2 < \cdots < t_k\) los tiempos únicos de eventos observados. En cada tiempo \(t_j\):
- \(d_j\) = número de eventos
- \(n_j\) = número de individuos en riesgo (no han experimentado el evento ni han sido censurados antes de \(t_j\))
El estimador de Kaplan-Meier es:
\[ \hat{S}(t) = \prod_{j: t_j \leq t} \left(1 - \frac{d_j}{n_j}\right) \]
3.1.2 Varianza del estimador
La varianza se estima mediante la fórmula de Greenwood (2):
\[ \widehat{\text{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{j: t_j \leq t} \frac{d_j}{n_j(n_j - d_j)} \]
El error estándar es:
\[ \widehat{\text{SE}}[\hat{S}(t)] = \sqrt{\widehat{\text{Var}}[\hat{S}(t)]} \]
3.1.3 Intervalos de confianza
El intervalo de confianza al \((1-\alpha) \times 100\%\) para \(S(t)\) se puede construir usando:
Método lineal (puede producir límites fuera de \([0,1]\)):
\[ \hat{S}(t) \pm z_{\alpha/2} \times \widehat{\text{SE}}[\hat{S}(t)] \]
Método log-log (recomendado, garantiza \(\hat{S}(t) \in [0,1]\)):
\[ \hat{S}(t)^{\exp\left(\pm \frac{z_{\alpha/2}}{\log \hat{S}(t) \times \widehat{\text{SE}}[\hat{S}(t)]}\right)} \]
3.1.4 Ejemplo: datos de leucemia
Utilizamos los datos clásicos de Freireich et al. (1963) sobre tiempos de remisión en leucemia aguda.
# Cargar datos
data(aml, package = "survival")
# Crear objeto de supervivencia
km_fit <- survfit(Surv(time, status) ~ x, data = aml)
# Resumen
summary(km_fit)Call: survfit(formula = Surv(time, status) ~ x, data = aml)
x=Maintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
9 11 1 0.909 0.0867 0.7541 1.000
13 10 1 0.818 0.1163 0.6192 1.000
18 8 1 0.716 0.1397 0.4884 1.000
23 7 1 0.614 0.1526 0.3769 0.999
31 5 1 0.491 0.1642 0.2549 0.946
34 4 1 0.368 0.1627 0.1549 0.875
48 2 1 0.184 0.1535 0.0359 0.944
x=Nonmaintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
5 12 2 0.8333 0.1076 0.6470 1.000
8 10 2 0.6667 0.1361 0.4468 0.995
12 8 1 0.5833 0.1423 0.3616 0.941
23 6 1 0.4861 0.1481 0.2675 0.883
27 5 1 0.3889 0.1470 0.1854 0.816
30 4 1 0.2917 0.1387 0.1148 0.741
33 3 1 0.1944 0.1219 0.0569 0.664
43 2 1 0.0972 0.0919 0.0153 0.620
45 1 1 0.0000 NaN NA NA# Gráfico
ggsurvplot(
km_fit,
data = aml,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
xlab = "Tiempo (semanas)",
ylab = "Probabilidad de supervivencia",
title = "Curvas de Kaplan-Meier por grupo de tratamiento",
legend.title = "Grupo",
legend.labs = c("Mantenido", "No mantenido"),
ggtheme = theme_minimal()
)
3.1.5 Tiempo de supervivencia mediano
El tiempo de supervivencia mediano es el tiempo \(t_{0.5}\) tal que:
\[ \hat{S}(t_{0.5}) = 0.5 \]
# Calcular medianas
print(km_fit, print.rmean = TRUE)Call: survfit(formula = Surv(time, status) ~ x, data = aml)
n events rmean* se(rmean) median 0.95LCL 0.95UCL
x=Maintained 11 7 52.6 19.83 31 18 NA
x=Nonmaintained 12 11 22.7 4.18 23 8 NA
* restricted mean with upper limit = 161 3.2 Comparación de curvas de supervivencia
3.2.1 Test de log-rank
El test de log-rank (3) compara dos o más curvas de supervivencia. Es el más utilizado y da peso igual a todos los tiempos.
Hipótesis:
- \(H_0\): Las funciones de riesgo son iguales en todos los grupos (\(h_1(t) = h_2(t) = \cdots = h_K(t)\))
- \(H_1\): Al menos una función de riesgo es diferente
Estadístico de prueba: Sea \(O_{kj}\) el número observado de eventos en el grupo \(k\) en el tiempo \(t_j\), y \(E_{kj}\) el número esperado bajo \(H_0\):
\[ E_{kj} = n_{kj} \times \frac{d_j}{n_j} \]
El estadístico de log-rank es:
\[ \chi^2 = \sum_{k=1}^K \frac{(O_k - E_k)^2}{E_k} \sim \chi^2_{K-1} \]
donde \(O_k = \sum_j O_{kj}\) y \(E_k = \sum_j E_{kj}\).
# Test de log-rank
survdiff(Surv(time, status) ~ x, data = aml)Call:
survdiff(formula = Surv(time, status) ~ x, data = aml)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 7 10.69 1.27 3.4
x=Nonmaintained 12 11 7.31 1.86 3.4
Chisq= 3.4 on 1 degrees of freedom, p= 0.07 3.2.2 Test de Wilcoxon-Breslow-Gehan
Variante que da más peso a tiempos tempranos:
\[ W = \sum_{k=1}^K \sum_j n_j (O_{kj} - E_{kj}) \]
# Test de Wilcoxon
survdiff(Surv(time, status) ~ x, data = aml, rho = 1)Call:
survdiff(formula = Surv(time, status) ~ x, data = aml, rho = 1)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 3.85 6.14 0.859 2.78
x=Nonmaintained 12 7.18 4.88 1.081 2.78
Chisq= 2.8 on 1 degrees of freedom, p= 0.1 3.2.3 Test de Tarone-Ware
Compromiso entre log-rank y Wilcoxon, usando pesos \(w_j = \sqrt{n_j}\).
3.3 Estimador de Nelson-Aalen
Alternativa al estimador de Kaplan-Meier, estima directamente el riesgo acumulado:
\[ \hat{H}(t) = \sum_{j: t_j \leq t} \frac{d_j}{n_j} \]
La función de supervivencia se estima como:
\[ \hat{S}(t) = \exp(-\hat{H}(t)) \]
# Estimador de Nelson-Aalen
na_fit <- survfit(Surv(time, status) ~ 1, data = aml, type = "fleming-harrington")
# Comparación
par(mfrow = c(1, 2))
plot(km_fit, conf.int = FALSE, main = "Kaplan-Meier", xlab = "Tiempo", ylab = "S(t)")
plot(na_fit, conf.int = FALSE, main = "Nelson-Aalen", xlab = "Tiempo", ylab = "S(t)")
4 Modelos paramétricos
Los modelos paramétricos asumen que los tiempos de supervivencia siguen una distribución de probabilidad específica.
4.1 Distribución exponencial
La distribución más simple. Asume riesgo constante \(h(t) = \lambda\).
4.1.1 Funciones
\[ \begin{aligned} f(t) &= \lambda \exp(-\lambda t) \\ S(t) &= \exp(-\lambda t) \\ h(t) &= \lambda \\ H(t) &= \lambda t \end{aligned} \]
Propiedad clave: No tiene memoria. \(P(T > s + t \mid T > s) = P(T > t)\).
4.1.2 Media y varianza
\[ E[T] = \frac{1}{\lambda}, \quad \text{Var}(T) = \frac{1}{\lambda^2} \]
4.1.3 Estimación
El estimador de máxima verosimilitud para \(\lambda\) con censura es:
\[ \hat{\lambda} = \frac{\sum_{i=1}^n \delta_i}{\sum_{i=1}^n y_i} \]
donde \(\delta_i\) es el indicador de evento e \(y_i\) el tiempo observado.
# Ajuste exponencial
exp_fit <- survreg(Surv(time, status) ~ x, data = aml, dist = "exponential")
summary(exp_fit)
Call:
survreg(formula = Surv(time, status) ~ x, data = aml, dist = "exponential")
Value Std. Error z p
(Intercept) 4.101 0.378 10.85 <2e-16
xNonmaintained -0.958 0.483 -1.98 0.048
Scale fixed at 1
Exponential distribution
Loglik(model)= -81.3 Loglik(intercept only)= -83.3
Chisq= 4.06 on 1 degrees of freedom, p= 0.044
Number of Newton-Raphson Iterations: 4
n= 23 # Extraer lambda
lambda <- 1 / exp(coef(exp_fit)[1])
cat("Lambda estimado:", lambda, "\n")Lambda estimado: 0.01654846 cat("Tiempo medio:", 1/lambda, "\n")Tiempo medio: 60.42857 4.2 Distribución de Weibull
Generalización de la exponencial que permite riesgo creciente, decreciente o constante.
4.2.1 Funciones
\[ \begin{aligned} f(t) &= \lambda \gamma t^{\gamma-1} \exp(-\lambda t^\gamma) \\ S(t) &= \exp(-\lambda t^\gamma) \\ h(t) &= \lambda \gamma t^{\gamma-1} \\ H(t) &= \lambda t^\gamma \end{aligned} \]
Parámetros:
- \(\lambda > 0\): parámetro de escala
- \(\gamma > 0\): parámetro de forma
- \(\gamma = 1\): exponencial (riesgo constante)
- \(\gamma < 1\): riesgo decreciente
- \(\gamma > 1\): riesgo creciente
4.2.2 Parametrización alternativa
En R, survreg usa la parametrización:
\[ \log T = \mu + \sigma W \]
donde \(W\) sigue una distribución de valor extremo. La relación es:
\[ \gamma = \frac{1}{\sigma}, \quad \lambda = \exp(-\mu/\sigma) \]
# Ajuste Weibull
weib_fit <- survreg(Surv(time, status) ~ x, data = aml, dist = "weibull")
summary(weib_fit)
Call:
survreg(formula = Surv(time, status) ~ x, data = aml, dist = "weibull")
Value Std. Error z p
(Intercept) 4.109 0.300 13.70 <2e-16
xNonmaintained -0.929 0.383 -2.43 0.015
Log(scale) -0.235 0.178 -1.32 0.188
Scale= 0.791
Weibull distribution
Loglik(model)= -80.5 Loglik(intercept only)= -83.2
Chisq= 5.31 on 1 degrees of freedom, p= 0.021
Number of Newton-Raphson Iterations: 5
n= 23 # Extraer parámetros
gamma <- 1 / weib_fit$scale
lambda <- exp(-coef(weib_fit)[1] / weib_fit$scale)
cat("Gamma (forma):", gamma, "\n")Gamma (forma): 1.264295 cat("Lambda (escala):", lambda, "\n")Lambda (escala): 0.005543889 # Verificación gráfica: log-log plot
# Si Weibull es apropiada, debe ser lineal
km_all <- survfit(Surv(time, status) ~ 1, data = aml)
plot(log(km_all$time), log(-log(km_all$surv)),
xlab = "log(tiempo)", ylab = "log(-log(S(t)))",
main = "Diagnóstico Weibull: gráfico log-log")
abline(lm(log(-log(km_all$surv)) ~ log(km_all$time)), col = "red")
4.3 Distribución log-normal
Los logaritmos de los tiempos siguen una distribución normal.
4.3.1 Funciones
Si \(\log T \sim N(\mu, \sigma^2)\):
\[ \begin{aligned} f(t) &= \frac{1}{t \sigma \sqrt{2\pi}} \exp\left(-\frac{(\log t - \mu)^2}{2\sigma^2}\right) \\ S(t) &= 1 - \Phi\left(\frac{\log t - \mu}{\sigma}\right) \end{aligned} \]
donde \(\Phi\) es la función de distribución acumulada normal estándar.
Característica: La función de riesgo \(h(t)\) primero aumenta y luego disminuye (forma de campana).
# Ajuste log-normal
lognorm_fit <- survreg(Surv(time, status) ~ x, data = aml, dist = "lognormal")
summary(lognorm_fit)
Call:
survreg(formula = Surv(time, status) ~ x, data = aml, dist = "lognormal")
Value Std. Error z p
(Intercept) 3.579 0.285 12.57 <2e-16
xNonmaintained -0.724 0.380 -1.90 0.057
Log(scale) -0.145 0.170 -0.86 0.391
Scale= 0.865
Log Normal distribution
Loglik(model)= -78.9 Loglik(intercept only)= -80.7
Chisq= 3.49 on 1 degrees of freedom, p= 0.062
Number of Newton-Raphson Iterations: 4
n= 23 4.4 Distribución log-logística
Similar a log-normal pero con función de supervivencia de forma cerrada.
4.4.1 Funciones
\[ \begin{aligned} S(t) &= \frac{1}{1 + \lambda t^\gamma} \\ h(t) &= \frac{\lambda \gamma t^{\gamma-1}}{1 + \lambda t^\gamma} \end{aligned} \]
Característica: También tiene función de riesgo que primero aumenta y luego disminuye cuando \(\gamma > 1\).
# Ajuste log-logístico
loglog_fit <- survreg(Surv(time, status) ~ x, data = aml, dist = "loglogistic")
summary(loglog_fit)
Call:
survreg(formula = Surv(time, status) ~ x, data = aml, dist = "loglogistic")
Value Std. Error z p
(Intercept) 3.503 0.288 12.18 <2e-16
xNonmaintained -0.604 0.393 -1.54 0.1243
Log(scale) -0.667 0.192 -3.48 0.0005
Scale= 0.513
Log logistic distribution
Loglik(model)= -79.4 Loglik(intercept only)= -80.6
Chisq= 2.41 on 1 degrees of freedom, p= 0.12
Number of Newton-Raphson Iterations: 3
n= 23 4.5 Comparación de modelos
4.5.1 Criterios de información
AIC (Akaike Information Criterion):
\[ \text{AIC} = -2 \log L + 2p \]
BIC (Bayesian Information Criterion):
\[ \text{BIC} = -2 \log L + p \log n \]
donde \(L\) es la verosimilitud, \(p\) el número de parámetros y \(n\) el tamaño muestral.
Regla: Menor valor indica mejor ajuste.
# Comparar modelos
modelos <- list(
Exponencial = exp_fit,
Weibull = weib_fit,
LogNormal = lognorm_fit,
LogLogistic = loglog_fit
)
comparacion <- data.frame(
Modelo = names(modelos),
AIC = sapply(modelos, AIC),
BIC = sapply(modelos, BIC),
LogLik = sapply(modelos, logLik)
)
kable(comparacion, digits = 2, caption = "Comparación de modelos paramétricos")| Modelo | AIC | BIC | LogLik | |
|---|---|---|---|---|
| Exponencial | Exponencial | 166.57 | 168.85 | -81.29 |
| Weibull | Weibull | 167.04 | 170.45 | -80.52 |
| LogNormal | LogNormal | 163.86 | 167.26 | -78.93 |
| LogLogistic | LogLogistic | 164.71 | 168.11 | -79.35 |
4.5.2 Gráfico comparativo
# Función para predecir supervivencia
pred_surv <- function(fit, newdata, times) {
predict(fit, newdata = newdata, type = "quantile", p = 1 - times)
}
# Tiempos
tiempos <- seq(0.1, 0.99, by = 0.01)
# KM para comparación
km_ref <- survfit(Surv(time, status) ~ 1, data = aml)
# Gráfico
plot(km_ref, conf.int = FALSE, xlab = "Tiempo", ylab = "S(t)",
main = "Comparación de modelos")
legend("topright",
legend = c("Kaplan-Meier", "Exponencial", "Weibull", "Log-normal", "Log-logística"),
col = c("black", "red", "blue", "green", "purple"),
lty = 1)
5 Modelo de riesgos proporcionales de Cox
El modelo de Cox (4) es el modelo semiparamétrico más utilizado en análisis de supervivencia.
5.1 Especificación del modelo
El modelo asume que la función de riesgo para el individuo \(i\) es:
\[ h_i(t \mid \mathbf{X}_i) = h_0(t) \exp(\boldsymbol{\beta}' \mathbf{X}_i) \]
donde:
- \(h_0(t)\) es la función de riesgo basal (no especificada)
- \(\mathbf{X}_i = (X_{i1}, \ldots, X_{ip})'\) es el vector de covariables
- \(\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)'\) es el vector de coeficientes
5.1.1 Razón de riesgos (Hazard Ratio)
Comparando dos individuos con covariables \(\mathbf{X}_i\) y \(\mathbf{X}_j\):
\[ \text{HR}_{ij}(t) = \frac{h_i(t \mid \mathbf{X}_i)}{h_j(t \mid \mathbf{X}_j)} = \frac{h_0(t) \exp(\boldsymbol{\beta}' \mathbf{X}_i)}{h_0(t) \exp(\boldsymbol{\beta}' \mathbf{X}_j)} = \exp[\boldsymbol{\beta}'(\mathbf{X}_i - \mathbf{X}_j)] \]
Propiedad clave: La razón de riesgos es constante en el tiempo (asunción de riesgos proporcionales).
5.1.2 Interpretación de coeficientes
Para una covariable continua \(X_k\):
\[ \text{HR} = \exp(\beta_k) \]
- \(\beta_k > 0\) (\(\text{HR} > 1\)): Aumento del riesgo
- \(\beta_k < 0\) (\(\text{HR} < 1\)): Disminución del riesgo
- \(\beta_k = 0\) (\(\text{HR} = 1\)): Sin efecto
Ejemplo: Si \(\beta_1 = 0.5\), entonces \(\text{HR} = \exp(0.5) = 1.65\), indicando un aumento del 65% en el riesgo por cada unidad de aumento en \(X_1\).
5.2 Verosimilitud parcial
La estimación de \(\boldsymbol{\beta}\) se realiza mediante verosimilitud parcial (5):
\[ L(\boldsymbol{\beta}) = \prod_{i: \delta_i = 1} \frac{\exp(\boldsymbol{\beta}' \mathbf{X}_i)}{\sum_{j \in R(t_i)} \exp(\boldsymbol{\beta}' \mathbf{X}_j)} \]
donde \(R(t_i)\) es el conjunto de riesgo en el tiempo \(t_i\) (individuos que no han experimentado el evento ni han sido censurados antes de \(t_i\)).
5.3 Ejemplo: datos de cáncer de pulmón
# Datos
data(lung, package = "survival")
lung <- na.omit(lung)
# Ajuste del modelo
cox_fit <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno +
pat.karno + meal.cal + wt.loss,
data = lung)
summary(cox_fit)Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + ph.karno +
pat.karno + meal.cal + wt.loss, data = lung)
n= 167, number of events= 120
coef exp(coef) se(coef) z Pr(>|z|)
age 1.080e-02 1.011e+00 1.160e-02 0.931 0.35168
sex -5.536e-01 5.749e-01 2.016e-01 -2.746 0.00603 **
ph.ecog 7.395e-01 2.095e+00 2.250e-01 3.287 0.00101 **
ph.karno 2.244e-02 1.023e+00 1.123e-02 1.998 0.04575 *
pat.karno -1.207e-02 9.880e-01 8.116e-03 -1.488 0.13685
meal.cal 2.835e-05 1.000e+00 2.594e-04 0.109 0.91298
wt.loss -1.420e-02 9.859e-01 7.766e-03 -1.828 0.06748 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0109 0.9893 0.9881 1.0341
sex 0.5749 1.7395 0.3872 0.8534
ph.ecog 2.0950 0.4773 1.3479 3.2560
ph.karno 1.0227 0.9778 1.0004 1.0455
pat.karno 0.9880 1.0121 0.9724 1.0038
meal.cal 1.0000 1.0000 0.9995 1.0005
wt.loss 0.9859 1.0143 0.9710 1.0010
Concordance= 0.653 (se = 0.029 )
Likelihood ratio test= 28.16 on 7 df, p=2e-04
Wald test = 27.5 on 7 df, p=3e-04
Score (logrank) test = 28.31 on 7 df, p=2e-045.3.1 Intervalos de confianza
Intervalos de confianza al 95% para HR:
# IC para HR
exp(confint(cox_fit)) 2.5 % 97.5 %
age 0.9881389 1.0341075
sex 0.3872369 0.8534079
ph.ecog 1.3479247 3.2559956
ph.karno 1.0004241 1.0454545
pat.karno 0.9724065 1.0038410
meal.cal 0.9995201 1.0005369
wt.loss 0.9710067 1.00102195.4 Estratificación
Cuando la asunción de riesgos proporcionales no se cumple para una covariable, se puede estratificar por ella:
\[ h_{ij}(t \mid \mathbf{X}_{ij}) = h_{0j}(t) \exp(\boldsymbol{\beta}' \mathbf{X}_{ij}) \]
donde \(j\) indexa el estrato.
# Modelo estratificado por sexo
cox_strat <- coxph(Surv(time, status) ~ age + strata(sex) + ph.ecog,
data = lung)
summary(cox_strat)Call:
coxph(formula = Surv(time, status) ~ age + strata(sex) + ph.ecog,
data = lung)
n= 167, number of events= 120
coef exp(coef) se(coef) z Pr(>|z|)
age 0.007446 1.007474 0.011088 0.671 0.501902
ph.ecog 0.458227 1.581268 0.138465 3.309 0.000935 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.007 0.9926 0.9858 1.030
ph.ecog 1.581 0.6324 1.2054 2.074
Concordance= 0.619 (se = 0.031 )
Likelihood ratio test= 13.45 on 2 df, p=0.001
Wald test = 13.51 on 2 df, p=0.001
Score (logrank) test = 13.73 on 2 df, p=0.0015.5 Diagnóstico del modelo
5.5.1 Verificación de riesgos proporcionales
5.5.1.1 Test de Schoenfeld
Los residuos de Schoenfeld permiten verificar la asunción de riesgos proporcionales:
\[ r_{ik} = \delta_i \left[X_{ik} - \frac{\sum_{j \in R(t_i)} X_{jk} \exp(\boldsymbol{\beta}' \mathbf{X}_j)}{\sum_{j \in R(t_i)} \exp(\boldsymbol{\beta}' \mathbf{X}_j)}\right] \]
Si la asunción se cumple, los residuos escalados deben tener pendiente cero al graficarlos contra el tiempo.
# Test de riesgos proporcionales
test_ph <- cox.zph(cox_fit)
print(test_ph) chisq df p
age 0.0581 1 0.810
sex 1.3567 1 0.244
ph.ecog 3.1105 1 0.078
ph.karno 5.9656 1 0.015
pat.karno 3.2117 1 0.073
meal.cal 6.5599 1 0.010
wt.loss 0.0833 1 0.773
GLOBAL 13.8427 7 0.054# Gráfico
par(mfrow = c(2, 2))
plot(test_ph)
Interpretación: Un p-valor < 0.05 sugiere violación de la asunción de riesgos proporcionales.
5.5.1.2 Gráficos log-log
Para covariables categóricas:
km_sex <- survfit(Surv(time, status) ~ sex, data = lung)
plot(km_sex, fun = "cloglog",
xlab = "log(Tiempo)", ylab = "log(-log(S(t)))",
main = "Verificación de riesgos proporcionales")
legend("topleft", legend = c("Hombre", "Mujer"), col = 1:2, lty = 1)
Si las líneas son paralelas, la asunción se cumple.
5.5.2 Residuos de martingala
Los residuos de martingala detectan la forma funcional de covariables continuas:
\[ M_i = \delta_i - \hat{H}_0(t_i) \exp(\boldsymbol{\beta}' \mathbf{X}_i) \]
mart_res <- residuals(cox_fit, type = "martingale")
plot(lung$age, mart_res,
xlab = "Edad", ylab = "Residuos de martingala",
main = "Verificación de forma funcional")
lines(lowess(lung$age, mart_res), col = "red", lwd = 2)
abline(h = 0, lty = 2)
5.5.3 Residuos de desviación
Los residuos de desviación son transformaciones simétricas de los residuos de martingala:
\[ D_i = \text{sign}(M_i) \sqrt{-2[M_i + \delta_i \log(\delta_i - M_i)]} \]
Útiles para detectar observaciones atípicas.
dev_res <- residuals(cox_fit, type = "deviance")
plot(dev_res,
xlab = "Índice", ylab = "Residuos de desviación",
main = "Detección de observaciones atípicas")
abline(h = c(-2, 0, 2), lty = c(2, 1, 2))
5.5.4 Estadísticos de influencia
DFBETAs: Miden el cambio en \(\hat{\boldsymbol{\beta}}\) al eliminar la observación \(i\).
dfb <- residuals(cox_fit, type = "dfbeta")
par(mfrow = c(2, 2))
for (j in 1:min(4, ncol(dfb))) {
plot(dfb[, j], ylab = paste("DFBETA", j),
main = colnames(dfb)[j])
abline(h = 0, lty = 2)
}
5.6 Extensiones del modelo de Cox
5.6.1 Covariables dependientes del tiempo
Permite que las covariables cambien con el tiempo:
\[ h_i(t \mid \mathbf{X}_i(t)) = h_0(t) \exp(\boldsymbol{\beta}' \mathbf{X}_i(t)) \]
# Ejemplo: efecto de tratamiento que cambia después de 6 meses
lung$trt_after6 <- lung$time > 180
# Formato tstart-tstop para covariables dependientes del tiempo
lung_td <- survSplit(Surv(time, status) ~ .,
data = lung,
cut = 180,
episode = "interval")
# Modelo
cox_td <- coxph(Surv(tstart, time, status) ~ age + sex + interval,
data = lung_td)
summary(cox_td)Call:
coxph(formula = Surv(tstart, time, status) ~ age + sex + interval,
data = lung_td)
n= 287, number of events= 120
coef exp(coef) se(coef) z Pr(>|z|)
age 0.01737 1.01752 0.01084 1.603 0.1089
sex -0.44669 0.63974 0.19754 -2.261 0.0237 *
interval NA NA 0.00000 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0175 0.9828 0.9961 1.0394
sex 0.6397 1.5631 0.4344 0.9422
interval NA NA NA NA
Concordance= 0.601 (se = 0.032 )
Likelihood ratio test= 8.88 on 2 df, p=0.01
Wald test = 8.47 on 2 df, p=0.01
Score (logrank) test = 8.6 on 2 df, p=0.015.6.2 Términos de interacción
Permiten que el efecto de una covariable dependa de otra:
# Interacción entre edad y sexo
cox_int <- coxph(Surv(time, status) ~ age * sex + ph.ecog, data = lung)
summary(cox_int)Call:
coxph(formula = Surv(time, status) ~ age * sex + ph.ecog, data = lung)
n= 167, number of events= 120
coef exp(coef) se(coef) z Pr(>|z|)
age 0.04141 1.04228 0.03245 1.276 0.201924
sex 1.03506 2.81528 1.39970 0.739 0.459611
ph.ecog 0.48187 1.61909 0.13868 3.475 0.000512 ***
age:sex -0.02444 0.97586 0.02211 -1.105 0.268992
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0423 0.9594 0.9781 1.111
sex 2.8153 0.3552 0.1812 43.747
ph.ecog 1.6191 0.6176 1.2337 2.125
age:sex 0.9759 1.0247 0.9345 1.019
Concordance= 0.648 (se = 0.03 )
Likelihood ratio test= 21.23 on 4 df, p=3e-04
Wald test = 21.6 on 4 df, p=2e-04
Score (logrank) test = 22 on 4 df, p=2e-045.6.3 Términos no lineales (splines)
Para capturar efectos no lineales:
# Efecto no lineal de edad usando splines
library(splines)
cox_spline <- coxph(Surv(time, status) ~ ns(age, df = 3) + sex + ph.ecog,
data = lung)
summary(cox_spline)Call:
coxph(formula = Surv(time, status) ~ ns(age, df = 3) + sex +
ph.ecog, data = lung)
n= 167, number of events= 120
coef exp(coef) se(coef) z Pr(>|z|)
ns(age, df = 3)1 -0.1268 0.8809 0.4256 -0.298 0.76581
ns(age, df = 3)2 1.6708 5.3164 1.3340 1.252 0.21040
ns(age, df = 3)3 0.8402 2.3169 0.5354 1.569 0.11660
sex -0.5187 0.5953 0.1985 -2.613 0.00899 **
ph.ecog 0.4691 1.5986 0.1386 3.386 0.00071 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
ns(age, df = 3)1 0.8809 1.1352 0.3825 2.0287
ns(age, df = 3)2 5.3164 0.1881 0.3891 72.6311
ns(age, df = 3)3 2.3169 0.4316 0.8112 6.6173
sex 0.5953 1.6798 0.4034 0.8785
ph.ecog 1.5986 0.6256 1.2184 2.0973
Concordance= 0.649 (se = 0.032 )
Likelihood ratio test= 22.12 on 5 df, p=5e-04
Wald test = 21.76 on 5 df, p=6e-04
Score (logrank) test = 22.32 on 5 df, p=5e-046 Tópicos avanzados
6.1 Riesgos competitivos
Los riesgos competitivos ocurren cuando un sujeto puede experimentar uno de varios eventos mutuamente excluyentes.
6.1.1 Función de incidencia acumulada
La función de incidencia acumulada (CIF) para el evento tipo \(k\) es:
\[ F_k(t) = P(T \leq t, \text{tipo } = k) = \int_0^t S(u) h_k(u) \, du \]
donde \(S(t)\) es la supervivencia global y \(h_k(t)\) es la causa-específica hazard.
Propiedad: \(\sum_{k=1}^K F_k(\infty) = 1\) si todos eventualmente experimentan algún evento.
6.1.2 Test de Gray
El test de Gray (6) compara CIFs entre grupos considerando los eventos competitivos.
# Simular datos con riesgos competitivos
set.seed(123)
n <- 200
tiempo <- rexp(n, 0.01)
evento <- sample(0:2, n, replace = TRUE, prob = c(0.3, 0.4, 0.3))
# 0 = censura, 1 = evento tipo 1, 2 = evento tipo 2
grupo <- rbinom(n, 1, 0.5)
datos_cr <- data.frame(tiempo, evento, grupo)
# Estimar CIF
cif_fit <- cuminc(ftime = datos_cr$tiempo,
fstatus = datos_cr$evento,
group = datos_cr$grupo)
# Gráfico
plot(cif_fit, lty = 1:2, col = 1:2,
xlab = "Tiempo", ylab = "Incidencia acumulada",
main = "Funciones de incidencia acumulada por grupo")
legend("topleft",
legend = c("Grupo 0, Evento 1", "Grupo 0, Evento 2",
"Grupo 1, Evento 1", "Grupo 1, Evento 2"),
lty = rep(1:2, 2), col = rep(1:2, each = 2))
# Test de Gray
print(cif_fit)Tests:
stat pv df
1 0.7391073 0.38994699 1
2 3.3594887 0.06681881 1
Estimates and Variances:
$est
200 400 600
0 1 0.5611684 0.6742783 0.6742783
1 1 0.4600541 0.5299475 NA
0 2 0.2228946 0.2228946 0.2743081
1 2 0.3709500 0.4116342 NA
$var
200 400 600
0 1 0.004210873 0.005942215 0.005942215
1 1 0.002719739 0.003628017 NA
0 2 0.002681503 0.002681503 0.005927133
1 2 0.002613671 0.003064468 NA
6.1.3 Modelo de subdistribución de Fine-Gray
El modelo de Fine-Gray (7) modela directamente la CIF:
\[ h_k^*(t \mid \mathbf{X}) = h_{0k}^*(t) \exp(\boldsymbol{\beta}_k' \mathbf{X}) \]
donde \(h_k^*(t)\) es el subdistribution hazard.
# Modelo de Fine-Gray
fg_fit <- crr(ftime = datos_cr$tiempo,
fstatus = datos_cr$evento,
cov1 = datos_cr$grupo,
failcode = 1)
summary(fg_fit)Competing Risks Regression
Call:
crr(ftime = datos_cr$tiempo, fstatus = datos_cr$evento, cov1 = datos_cr$grupo,
failcode = 1)
coef exp(coef) se(coef) z p-value
datos_cr$grupo1 -0.166 0.847 0.207 -0.799 0.42
exp(coef) exp(-coef) 2.5% 97.5%
datos_cr$grupo1 0.847 1.18 0.565 1.27
Num. cases = 200
Pseudo Log-likelihood = -426
Pseudo likelihood ratio test = 0.61 on 1 df,6.2 Modelos de fragilidad
Los modelos de fragilidad introducen heterogeneidad no observada mediante un efecto aleatorio:
\[ h_i(t \mid \mathbf{X}_i, Z_i) = Z_i h_0(t) \exp(\boldsymbol{\beta}' \mathbf{X}_i) \]
donde \(Z_i\) es la fragilidad (usualmente \(Z_i \sim \text{Gamma}(\theta, \theta)\) con \(E[Z_i] = 1\)).
Aplicaciones: Datos agrupados (pacientes dentro de hospitales) o eventos recurrentes.
# Modelo de fragilidad gamma
data(rats, package = "survival")
frailty_fit <- coxph(Surv(time, status) ~ rx + frailty(litter),
data = rats)
summary(frailty_fit)Call:
coxph(formula = Surv(time, status) ~ rx + frailty(litter), data = rats)
n= 300, number of events= 42
coef se(coef) se2 Chisq DF p
rx 0.7271 0.3183 0.3126 5.22 1.00 0.022
frailty(litter) 63.05 47.67 0.067
exp(coef) exp(-coef) lower .95 upper .95
rx 2.069 0.4833 1.109 3.861
Iterations: 6 outer, 23 Newton-Raphson
Variance of random effect= 2.021201 I-likelihood = -217.5
Degrees of freedom for terms= 1.0 47.7
Concordance= 0.923 (se = 0.014 )
Likelihood ratio test= 88.54 on 48.63 df, p=4e-046.3 Eventos recurrentes
Para sujetos que pueden experimentar el evento múltiples veces.
6.3.1 Procesos de conteo
Modelamos el número de eventos \(N_i(t)\) usando modelos de intensidad:
\[ \lambda_i(t \mid \mathbf{X}_i) = Y_i(t) h_0(t) \exp(\boldsymbol{\beta}' \mathbf{X}_i) \]
donde \(Y_i(t)\) indica si el sujeto \(i\) está en riesgo en el tiempo \(t\).
6.3.2 Modelos para eventos recurrentes
- Andersen-Gill (AG): Todos los eventos se tratan igualmente
- Prentice-Williams-Peterson (PWP): Estratifica por número de evento
- Wei-Lin-Weissfeld (WLW): Modelo marginal para cada evento
# Ejemplo con datos simulados de eventos recurrentes
datos_rec <- data.frame(
id = rep(1:50, each = 3),
evento = sample(0:1, 150, replace = TRUE),
tiempo = abs(rnorm(150, 10, 3)),
x1 = rep(rnorm(50), each = 3)
)
# Modelo Andersen-Gill
ag_fit <- coxph(Surv(tiempo, evento) ~ x1 + cluster(id),
data = datos_rec)
summary(ag_fit)Call:
coxph(formula = Surv(tiempo, evento) ~ x1, data = datos_rec,
cluster = id)
n= 150, number of events= 71
coef exp(coef) se(coef) robust se z Pr(>|z|)
x1 0.08985 1.09401 0.12400 0.13615 0.66 0.509
exp(coef) exp(-coef) lower .95 upper .95
x1 1.094 0.9141 0.8378 1.429
Concordance= 0.525 (se = 0.045 )
Likelihood ratio test= 0.53 on 1 df, p=0.5
Wald test = 0.44 on 1 df, p=0.5
Score (logrank) test = 0.52 on 1 df, p=0.5, Robust = 0.4 p=0.5
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).6.4 Análisis de concordancia
El índice de concordancia (C-index) mide la capacidad discriminativa del modelo:
\[ C = P(\hat{T}_i > \hat{T}_j \mid T_i > T_j) \]
- \(C = 0.5\): Sin capacidad discriminativa (aleatorio)
- \(C = 1\): Discriminación perfecta
# Calcular C-index
c_index <- concordance(cox_fit)
print(c_index)Call:
concordance.coxph(object = cox_fit)
n= 167
Concordance= 0.6532 se= 0.0289
concordant discordant tied.x tied.y tied.xy
6900 3664 0 11 0 # Usando paquete survival
survConcordance(Surv(time, status) ~ predict(cox_fit), data = lung)$concordance
concordant
0.6531617
$stats
concordant discordant tied.risk tied.time std(c-d)
6900.000 3664.000 0.000 11.000 648.288
$n
[1] 167
$std.err
std(c-d)
0.03068383
$call
survConcordance(formula = Surv(time, status) ~ predict(cox_fit),
data = lung)
attr(,"class")
[1] "survConcordance"6.5 Validación del modelo
6.5.1 Validación cruzada
# Validación cruzada 10-fold
set.seed(456)
k <- 10
folds <- sample(1:k, nrow(lung), replace = TRUE)
c_indices <- numeric(k)
for (i in 1:k) {
train <- lung[folds != i, ]
test <- lung[folds == i, ]
modelo <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = train)
pred <- predict(modelo, newdata = test)
c_indices[i] <- concordance(Surv(time, status) ~ pred, data = test)$concordance
}
cat("C-index medio (CV):", mean(c_indices), "\n")C-index medio (CV): 0.3761658 cat("Desviación estándar:", sd(c_indices), "\n")Desviación estándar: 0.09477215 6.5.2 Calibración
La calibración evalúa si las probabilidades predichas coinciden con las frecuencias observadas.
# Gráfico de calibración
library(rms)
dd <- datadist(lung)
options(datadist = "dd")
cox_rms <- cph(Surv(time, status) ~ age + sex + ph.ecog,
data = lung, x = TRUE, y = TRUE, surv = TRUE, time.inc = 365)
cal <- calibrate(cox_rms, u = 365, B = 50)Using Cox survival estimates at 365 Daysplot(cal, main = "Gráfico de calibración")
7 Predicción de supervivencia
7.1 Predicción individual
Para un nuevo paciente con covariables \(\mathbf{X}_{\text{new}}\), la función de supervivencia predicha es:
\[ \hat{S}(t \mid \mathbf{X}_{\text{new}}) = \hat{S}_0(t)^{\exp(\hat{\boldsymbol{\beta}}' \mathbf{X}_{\text{new}})} \]
donde \(\hat{S}_0(t)\) es la función de supervivencia basal estimada.
# Definir perfiles de pacientes
perfiles <- data.frame(
age = c(50, 50, 70, 70),
sex = c(1, 2, 1, 2),
ph.ecog = c(0, 0, 2, 2),
ph.karno = c(100, 100, 100, 100),
pat.karno = c(100, 100, 100, 100),
meal.cal = c(1000, 1000, 1000, 1000),
wt.loss = c(0, 0, 0, 0)
)
# Predicciones
pred_surv <- survfit(cox_fit, newdata = perfiles)
# Gráfico
plot(pred_surv, col = 1:4, lwd = 2,
xlab = "Tiempo (días)", ylab = "Probabilidad de supervivencia",
main = "Curvas de supervivencia predichas")
legend("topright",
legend = c("50 años, Hombre, ECOG 0",
"50 años, Mujer, ECOG 0",
"70 años, Hombre, ECOG 2",
"70 años, Mujer, ECOG 2"),
col = 1:4, lwd = 2)
7.2 Probabilidades de supervivencia
# Probabilidad de sobrevivir 1 año (365 días)
summary(pred_surv, times = 365)Call: survfit(formula = cox_fit, newdata = perfiles)
time n.risk n.event survival1 survival2 survival3 survival4
365 49 87 0.535 0.698 0.0332 0.141# Probabilidad de sobrevivir 2 años
summary(pred_surv, times = 730)Call: survfit(formula = cox_fit, newdata = perfiles)
time n.risk n.event survival1 survival2 survival3 survival4
730 10 116 0.198 0.395 0.000149 0.006327.3 Nomogramas
Los nomogramas permiten cálculos gráficos de probabilidades de supervivencia.
library(rms)
dd <- datadist(lung)
options(datadist = "dd")
# Modelo usando rms
cox_nom <- cph(Surv(time, status) ~ age + sex + ph.ecog,
data = lung, x = TRUE, y = TRUE, surv = TRUE, time.inc = 365)
# Crear perfiles de pacientes para predicción
perfiles_nom <- data.frame(
age = c(40, 50, 60, 70, 80),
sex = rep(1, 5),
ph.ecog = rep(1, 5)
)
# Calcular predicciones de supervivencia a 1 año y 2 años
survest_1yr <- survest(cox_nom, newdata = perfiles_nom, times = 365)
survest_2yr <- survest(cox_nom, newdata = perfiles_nom, times = 730)
# Crear tabla de predicciones
pred_table <- data.frame(
Edad = perfiles_nom$age,
Sexo = "Masculino",
ECOG = perfiles_nom$ph.ecog,
Supervivencia_1_año = round(survest_1yr$surv, 3),
Supervivencia_2_años = round(survest_2yr$surv, 3)
)
kable(pred_table,
caption = "Probabilidades de supervivencia predichas según edad del paciente",
col.names = c("Edad", "Sexo", "ECOG", "1 año", "2 años"))| Edad | Sexo | ECOG | 1 año | 2 años |
|---|---|---|---|---|
| 40 | Masculino | 1 | 0.406 | 0.109 |
| 50 | Masculino | 1 | 0.377 | 0.091 |
| 60 | Masculino | 1 | 0.347 | 0.074 |
| 70 | Masculino | 1 | 0.318 | 0.060 |
| 80 | Masculino | 1 | 0.288 | 0.047 |
8 Consideraciones prácticas
8.1 Tamaño muestral
El tamaño muestral requerido depende del número de eventos, no del número de sujetos.
Regla general: Se necesitan al menos 10 eventos por covariable para evitar sobreajuste.
\[ n_{\text{eventos}} \geq 10 \times p \]
donde \(p\) es el número de covariables.
8.2 Manejo de datos faltantes
8.2.1 Análisis de casos completos
Elimina sujetos con valores faltantes. Problema: Pérdida de información y posible sesgo.
8.2.2 Imputación múltiple
Genera varios conjuntos de datos imputados y combina los resultados.
# Ejemplo de imputación múltiple con mice
library(mice)
# Crear valores faltantes artificiales
lung_na <- lung
lung_na$age[sample(1:nrow(lung_na), 20)] <- NA
# Imputación múltiple
imp <- mice(lung_na, m = 5, maxit = 10, seed = 123, printFlag = FALSE)
# Ajustar modelo en cada conjunto imputado
fit_imp <- with(imp, coxph(Surv(time, status) ~ age + sex + ph.ecog))
# Combinar resultados
pool_fit <- pool(fit_imp)
summary(pool_fit)8.3 Selección de variables
8.3.1 Método backward
# Modelo completo
full_model <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno +
pat.karno + meal.cal + wt.loss,
data = lung)
# Selección backward usando AIC
step_model <- step(full_model, direction = "backward", trace = 0)
summary(step_model)Call:
coxph(formula = Surv(time, status) ~ sex + ph.ecog + ph.karno +
pat.karno + wt.loss, data = lung)
n= 167, number of events= 120
coef exp(coef) se(coef) z Pr(>|z|)
sex -0.558190 0.572244 0.199202 -2.802 0.00508 **
ph.ecog 0.742983 2.102197 0.227604 3.264 0.00110 **
ph.karno 0.020366 1.020575 0.011080 1.838 0.06604 .
pat.karno -0.012401 0.987675 0.007978 -1.554 0.12008
wt.loss -0.014494 0.985611 0.007693 -1.884 0.05957 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
sex 0.5722 1.7475 0.3873 0.8456
ph.ecog 2.1022 0.4757 1.3457 3.2841
ph.karno 1.0206 0.9798 0.9987 1.0430
pat.karno 0.9877 1.0125 0.9724 1.0032
wt.loss 0.9856 1.0146 0.9709 1.0006
Concordance= 0.658 (se = 0.029 )
Likelihood ratio test= 27.28 on 5 df, p=5e-05
Wald test = 26.89 on 5 df, p=6e-05
Score (logrank) test = 27.64 on 5 df, p=4e-058.3.2 Regularización (LASSO)
library(glmnet)
# Preparar datos
x <- model.matrix(~ age + sex + ph.ecog + ph.karno +
pat.karno + meal.cal + wt.loss - 1, data = lung)
y <- Surv(lung$time, lung$status)
# LASSO
lasso_fit <- cv.glmnet(x, y, family = "cox", alpha = 1)
# Coeficientes
coef(lasso_fit, s = "lambda.min")7 x 1 sparse Matrix of class "dgCMatrix"
1
age 0.001091615
sex -0.264579494
ph.ecog 0.305668826
ph.karno .
pat.karno -0.001755917
meal.cal .
wt.loss . # Gráfico
plot(lasso_fit)
8.4 Reporte de resultados
8.4.1 Tabla de resultados
# Función para crear tabla de resultados
crear_tabla <- function(modelo) {
coef_df <- data.frame(
Variable = names(coef(modelo)),
Coeficiente = coef(modelo),
HR = exp(coef(modelo)),
IC_inf = exp(confint(modelo)[, 1]),
IC_sup = exp(confint(modelo)[, 2]),
p_valor = summary(modelo)$coefficients[, "Pr(>|z|)"]
)
rownames(coef_df) <- NULL
return(coef_df)
}
tabla <- crear_tabla(cox_fit)
kable(tabla, digits = 3,
caption = "Resultados del modelo de Cox: Hazard Ratios e intervalos de confianza al 95%")| Variable | Coeficiente | HR | IC_inf | IC_sup | p_valor |
|---|---|---|---|---|---|
| age | 0.011 | 1.011 | 0.988 | 1.034 | 0.352 |
| sex | -0.554 | 0.575 | 0.387 | 0.853 | 0.006 |
| ph.ecog | 0.740 | 2.095 | 1.348 | 3.256 | 0.001 |
| ph.karno | 0.022 | 1.023 | 1.000 | 1.045 | 0.046 |
| pat.karno | -0.012 | 0.988 | 0.972 | 1.004 | 0.137 |
| meal.cal | 0.000 | 1.000 | 1.000 | 1.001 | 0.913 |
| wt.loss | -0.014 | 0.986 | 0.971 | 1.001 | 0.067 |
8.4.2 Gráfico Forest Plot
library(forestplot)
# Extraer datos
coefs <- coef(cox_fit)
hrs <- exp(coefs)
ci <- exp(confint(cox_fit))
# Preparar datos
tabletext <- cbind(
c("Variable", names(coefs)),
c("HR", sprintf("%.2f", hrs)),
c("95% IC", sprintf("(%.2f-%.2f)", ci[, 1], ci[, 2]))
)
# Forest plot
forestplot(tabletext,
mean = c(NA, hrs),
lower = c(NA, ci[, 1]),
upper = c(NA, ci[, 2]),
xlog = TRUE,
xlab = "Hazard Ratio",
col = fpColors(box = "royalblue", line = "darkblue"))
9 Casos de estudio
9.1 Caso 1: Trasplante de médula ósea
# Datos simulados de trasplante
set.seed(789)
n_bmt <- 137
bmt <- data.frame(
time = abs(rnorm(n_bmt, 500, 300)),
status = sample(0:1, n_bmt, replace = TRUE, prob = c(0.4, 0.6)),
age = rnorm(n_bmt, 35, 10),
donor = factor(sample(c("Relacionado", "No relacionado"), n_bmt, replace = TRUE)),
disease = factor(sample(c("AML", "ALL", "CML"), n_bmt, replace = TRUE))
)
# KM por tipo de donante
km_bmt <- survfit(Surv(time, status) ~ donor, data = bmt)
ggsurvplot(km_bmt, data = bmt,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
xlab = "Días desde trasplante",
ylab = "Probabilidad de supervivencia",
title = "Supervivencia por tipo de donante",
ggtheme = theme_minimal())
# Modelo multivariado
cox_bmt <- coxph(Surv(time, status) ~ age + donor + disease, data = bmt)
summary(cox_bmt)Call:
coxph(formula = Surv(time, status) ~ age + donor + disease, data = bmt)
n= 137, number of events= 76
coef exp(coef) se(coef) z Pr(>|z|)
age -0.01675 0.98339 0.01247 -1.342 0.1795
donorRelacionado -0.48890 0.61330 0.24781 -1.973 0.0485 *
diseaseAML 0.09609 1.10086 0.29521 0.325 0.7448
diseaseCML 0.48582 1.62551 0.28686 1.694 0.0903 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 0.9834 1.0169 0.9596 1.0077
donorRelacionado 0.6133 1.6305 0.3773 0.9968
diseaseAML 1.1009 0.9084 0.6172 1.9634
diseaseCML 1.6255 0.6152 0.9264 2.8521
Concordance= 0.598 (se = 0.038 )
Likelihood ratio test= 8.03 on 4 df, p=0.09
Wald test = 8.07 on 4 df, p=0.09
Score (logrank) test = 8.17 on 4 df, p=0.09
9.2 Caso 2: Cáncer de mama
# Datos simulados de cáncer de mama con recaída
set.seed(321)
n_mama <- 200
mama <- data.frame(
id = 1:n_mama,
time = abs(rnorm(n_mama, 1000, 400)),
status = sample(0:1, n_mama, replace = TRUE, prob = c(0.3, 0.7)),
age = rnorm(n_mama, 55, 12),
size = exp(rnorm(n_mama, log(2.5), 0.5)),
nodes = rpois(n_mama, 3),
grade = factor(sample(1:3, n_mama, replace = TRUE)),
er = factor(sample(c("Positivo", "Negativo"), n_mama, replace = TRUE, prob = c(0.7, 0.3)))
)
# Modelo de Cox
cox_mama <- coxph(Surv(time, status) ~ age + size + nodes + grade + er,
data = mama)
summary(cox_mama)Call:
coxph(formula = Surv(time, status) ~ age + size + nodes + grade +
er, data = mama)
n= 200, number of events= 134
coef exp(coef) se(coef) z Pr(>|z|)
age -0.006931 0.993093 0.007690 -0.901 0.367
size -0.007479 0.992549 0.068401 -0.109 0.913
nodes -0.033845 0.966721 0.051670 -0.655 0.512
grade2 -0.076670 0.926195 0.209897 -0.365 0.715
grade3 -0.128404 0.879498 0.225124 -0.570 0.568
erPositivo -0.354447 0.701562 0.187682 -1.889 0.059 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 0.9931 1.007 0.9782 1.008
size 0.9925 1.008 0.8680 1.135
nodes 0.9667 1.034 0.8736 1.070
grade2 0.9262 1.080 0.6138 1.398
grade3 0.8795 1.137 0.5657 1.367
erPositivo 0.7016 1.425 0.4856 1.013
Concordance= 0.548 (se = 0.03 )
Likelihood ratio test= 4.68 on 6 df, p=0.6
Wald test = 4.79 on 6 df, p=0.6
Score (logrank) test = 4.82 on 6 df, p=0.6# Validación de riesgos proporcionales
test_ph_mama <- cox.zph(cox_mama)
print(test_ph_mama) chisq df p
age 0.0401 1 0.84
size 0.0133 1 0.91
nodes 0.1069 1 0.74
grade 1.6304 2 0.44
er 0.6103 1 0.43
GLOBAL 2.3862 6 0.8810 Software y recursos
10.1 Paquetes de R
10.2 Recursos adicionales
10.2.1 Libros
- Hosmer, Lemeshow & May (2008). Applied Survival Analysis
- Kleinbaum & Klein (2012). Survival Analysis: A Self-Learning Text
- Therneau & Grambsch (2000). Modeling Survival Data
- Klein & Moeschberger (2003). Survival Analysis: Techniques for Censored and Truncated Data
10.2.2 Cursos online
10.2.3 Documentación
11 Conclusiones
El análisis de supervivencia proporciona herramientas estadísticas poderosas para analizar datos de tiempo hasta el evento en presencia de censura. Los puntos clave son:
- Censura: Requiere métodos especializados que no ignoran la información parcial
- Kaplan-Meier: Método no paramétrico para estimar funciones de supervivencia
- Test de log-rank: Prueba estándar para comparar curvas de supervivencia
- Modelos paramétricos: Útiles cuando se conoce la distribución subyacente
- Modelo de Cox: Modelo semiparamétrico más utilizado, flexible y robusto
- Diagnóstico: Verificación de asunciones es crucial (riesgos proporcionales, forma funcional)
- Extensiones: Riesgos competitivos, fragilidad y eventos recurrentes para situaciones más complejas
La elección del método depende de:
- Objetivo del análisis (descripción vs predicción)
- Estructura de los datos (censura, covariables, eventos competitivos)
- Asunciones que estemos dispuestos a hacer
- Tamaño muestral disponible
12 Referencias
1.
Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association. 1958;53(282):457-81.
2.
Greenwood M. The errors of sampling of the survivorship tables. Reports on Public Health and Statistical Subjects. 1926;33:1-26.
3.
Mantel N. Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemotherapy Reports. 1966;50(3):163-70.
4.
Cox DR. Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological). 1972;34(2):187-202.
5.
Cox DR. Partial likelihood. Biometrika. 1975;62(2):269-76.
6.
Gray RJ. A class of K-sample tests for comparing the cumulative incidence of a competing risk. The Annals of Statistics. 1988;16(3):1141-54.
7.
Fine JP, Gray RJ. A proportional hazards model for the subdistribution of a competing risk. Journal of the American Statistical Association. 1999;94(446):496-509.
8.
Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. Springer Science & Business Media; 2000.
9.
Kassambara A, Kosinski M, Biecek P. survminer: Drawing Survival Curves using ’ggplot2’ [Internet]. 2021. Disponible en: https://CRAN.R-project.org/package=survminer
10.
Jackson C. flexsurv: A Platform for Parametric Survival Modeling in R. Journal of Statistical Software. 2016.
11.
Harrell Jr FE. Regression modeling strategies: with applications to linear models, logistic and ordinal regression, and survival analysis. 2.ª ed. Springer; 2015.
12.
Gray B. cmprsk: Subdistribution Analysis of Competing Risks [Internet]. 2014. Disponible en: https://CRAN.R-project.org/package=cmprsk
13.
Therneau TM. coxme: Mixed Effects Cox Models [Internet]. 2022. Disponible en: https://CRAN.R-project.org/package=coxme
13 Apéndice A: Código R completo
# ============================================================
# ============================================================
# 1. INSTALACIÓN DE PAQUETES --------------------------------
install.packages(c("survival", "survminer", "flexsurv",
"KMsurv", "cmprsk", "rms", "ggplot2"))
# 2. CARGA DE LIBRERÍAS -------------------------------------
library(survival)
library(survminer)
library(flexsurv)
library(KMsurv)
library(cmprsk)
library(rms)
library(ggplot2)
# 3. ANÁLISIS NO PARAMÉTRICO --------------------------------
# Datos de ejemplo
data(aml, package = "survival")
# Kaplan-Meier
km_fit <- survfit(Surv(time, status) ~ x, data = aml)
summary(km_fit)
# Gráfico
ggsurvplot(km_fit, data = aml, pval = TRUE, conf.int = TRUE)
# Test de log-rank
survdiff(Surv(time, status) ~ x, data = aml)
# 4. MODELOS PARAMÉTRICOS -----------------------------------
# Exponencial
exp_fit <- survreg(Surv(time, status) ~ x, data = aml,
dist = "exponential")
# Weibull
weib_fit <- survreg(Surv(time, status) ~ x, data = aml,
dist = "weibull")
# Log-normal
lognorm_fit <- survreg(Surv(time, status) ~ x, data = aml,
dist = "lognormal")
# Comparar modelos
AIC(exp_fit, weib_fit, lognorm_fit)
# 5. MODELO DE COX ------------------------------------------
data(lung, package = "survival")
lung <- na.omit(lung)
# Ajuste
cox_fit <- coxph(Surv(time, status) ~ age + sex + ph.ecog,
data = lung)
summary(cox_fit)
# Diagnóstico
cox.zph(cox_fit)
plot(cox.zph(cox_fit))
# 6. PREDICCIÓN ---------------------------------------------
# Definir perfil
nuevo_paciente <- data.frame(age = 60, sex = 1, ph.ecog = 1)
# Predecir
pred <- survfit(cox_fit, newdata = nuevo_paciente)
plot(pred)
# Probabilidad a 1 año
summary(pred, times = 365)14 Apéndice B: Fórmulas matemáticas
14.1 Distribuciones comunes
| Distribución | \(h(t)\) | \(S(t)\) | Parámetros |
|---|---|---|---|
| Exponencial | \(\lambda\) | \(\exp(-\lambda t)\) | \(\lambda > 0\) |
| Weibull | \(\lambda \gamma t^{\gamma-1}\) | \(\exp(-\lambda t^\gamma)\) | \(\lambda, \gamma > 0\) |
| Log-normal | \(\frac{\phi\left(\frac{\log t - \mu}{\sigma}\right)}{\sigma t \left[1 - \Phi\left(\frac{\log t - \mu}{\sigma}\right)\right]}\) | \(1 - \Phi\left(\frac{\log t - \mu}{\sigma}\right)\) | \(\mu \in \mathbb{R}, \sigma > 0\) |
| Log-logística | \(\frac{\lambda \gamma t^{\gamma-1}}{1 + \lambda t^\gamma}\) | \(\frac{1}{1 + \lambda t^\gamma}\) | \(\lambda, \gamma > 0\) |
14.2 Tests estadísticos
14.2.1 Test de log-rank
\[ \chi^2 = \frac{(O_1 - E_1)^2}{E_1} + \frac{(O_2 - E_2)^2}{E_2} \sim \chi^2_1 \]
14.2.2 Test de Wald
Para \(H_0: \beta_j = 0\):
\[ z = \frac{\hat{\beta}_j}{\widehat{\text{SE}}(\hat{\beta}_j)} \sim N(0, 1) \]
14.2.3 Test de razón de verosimilitud
\[ -2 \log \Lambda = -2[\log L(\text{modelo restringido}) - \log L(\text{modelo completo})] \sim \chi^2_q \]
donde \(q\) es la diferencia en número de parámetros.
Información de la sesión:
sessionInfo()R version 4.6.0 (2026-04-24)
Platform: x86_64-apple-darwin20
Running under: macOS Monterey 12.7.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Madrid
tzcode source: internal
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] forestplot_3.2.0 abind_1.4-8 checkmate_2.3.4 glmnet_5.1
[5] Matrix_1.7-6 mice_3.19.0 cmprsk_2.2-12 rms_8.1-1
[9] Hmisc_5.3-0 knitr_1.52 dplyr_1.2.1 KMsurv_0.1-6
[13] flexsurv_2.3.2 survminer_0.5.2 ggpubr_1.0.0 ggplot2_4.0.3
[17] survival_3.8-12
loaded via a namespace (and not attached):
[1] Rdpack_2.6.6 gridExtra_2.3.1 sandwich_3.1-3
[4] rlang_1.3.0 magrittr_2.0.5 multcomp_1.4-32
[7] otel_0.2.0 polspline_1.1.25 compiler_4.6.0
[10] vctrs_0.7.3 quantreg_6.1 quadprog_1.5-8
[13] stringr_1.6.0 pkgconfig_2.0.3 shape_1.4.6.1
[16] fastmap_1.2.0 backports_1.5.1 labeling_0.4.3
[19] deSolve_1.42 rmarkdown_2.32 nloptr_2.2.1
[22] MatrixModels_0.5-4 purrr_1.2.2 xfun_0.60
[25] jomo_2.7-6 jsonlite_2.0.0 mstate_0.3.3
[28] pan_2.0 broom_1.0.13 cluster_2.1.8.3
[31] R6_2.6.1 stringi_1.8.9 RColorBrewer_1.1-3
[34] car_3.1-5 boot_1.3-32 rpart_4.1.27
[37] numDeriv_2016.8-1.1 Rcpp_1.1.2 iterators_1.0.14
[40] zoo_1.9-0 base64enc_0.1-6 splines_4.6.0
[43] nnet_7.3-21 tidyselect_1.2.1 rstudioapi_0.19.0
[46] yaml_2.3.12 ggtext_0.2.0 codetools_0.2-20
[49] lattice_0.23-1 tibble_3.3.1 withr_3.0.3
[52] S7_0.2.2 evaluate_1.0.5 foreign_0.8-91
[55] xml2_1.6.0 muhaz_1.2.6.5 pillar_1.11.1
[58] carData_3.0-6 stats4_4.6.0 foreach_1.5.2
[61] reformulas_0.4.4 generics_0.1.4 scales_1.4.0
[64] minqa_1.2.8 glue_1.8.1 tools_4.6.0
[67] data.table_1.18.6.1 lme4_2.0-6 SparseM_1.84-2
[70] ggsignif_0.6.4 mvtnorm_1.4-2 tidyr_1.3.2
[73] rbibutils_2.4.1 colorspace_2.1-3 nlme_3.1-171
[76] htmlTable_2.5.0 Formula_1.2-6 cli_3.6.6
[79] gtable_0.3.6 rstatix_1.1.0 digest_0.6.39
[82] TH.data_1.1-5 htmlwidgets_1.6.4 farver_2.1.2
[85] htmltools_0.5.9 lifecycle_1.0.5 mitml_0.4-5
[88] statmod_1.5.2 gridtext_0.1.6 MASS_7.3-66