L’analisi de supervivència es una branca de l’estadística que s’utilitza per estudiar el temps que passa fins que ocorre un esdeveniment. També se’n diu anàlisi de temps a l’esdeveniment.
Les aplicacions clíniques habituals son l’anàlisi del temps fins a la mort, recidiva, ingres hospitalàri o perdua de seguiment d’un pacient.
En estudis reals, sovint no coneixem el temps exacte de l’esdeveniment per a tothom, ni tothom comença i acaba un estudi al mateix temps. Per fer calculs amb els seguiments dels que se disposa, existeix el concepte de censura.
📌 Censura per la dreta vol dir que sabem que l’esdeveniment no ha passat fins a cert moment, però no sabem què passa després.
Un subjecte pot estar censurat per diversos motius:
Aquests són casos de censura per la dreta, que és la més habitual. Existeixen censura per l’esquerra o per intervals que no tractarem en aquest curs.
Per il·lustrar la censura i la supervivència, observem la següent figura:
Pregunta: quina proporció de pacients es mantenen lliures de l’esdeveniment al cap de 10 anys?
📌 Tot i que no sabem què ha passat amb els censurats, sí sabem que han estat seguits durant un temps determinat sense l’esdeveniment.
Per fer una anàlisi de supervivència en R, necessitem un conjunt de dades que inclogui com a mínim:
| Variable | Descripció |
|---|---|
temps |
Temps de seguiment per a cada subjecte (pot estar en dies, mesos, anys…) |
event |
Indicador de si ha passat l’esdeveniment (1) o si és censurat (0) |
📌 Aquestes dues variables s’utilitzen per construir un objecte especial de tipus Surv(), que encapsula tota la informació necessària per als mètodes de supervivència: el temps, l’estat de censura i (més endavant) grups comparatius.
Surv()?Un objecte Surv() és un format especial que R utilitza per representar dades de supervivència.
Conté, per a cada cas:
Aquest objecte és imprescindible per poder fer servir funcions com survfit(), coxph() o survdiff() del paquet survival.
Abans de calcular el temps, cal assegurar-se que les dates estiguin en format correcte. A R això pot ser complicat, però hi ha un paquet molt útil per treballar amb dates:
lubridate és un paquet de la família tidyverse que facilita molt la conversió de textos en dates, així com càlculs, extracció d’anys, mesos, etc.
Per exemple, ymd("2020-01-01") crea una data en format “any-mes-dia”.
library(lubridate)
library(survival)
ingres <- ymd(c("2020-01-01", "2020-02-15", "2020-03-01", "2020-01-20", "2020-01-10"))
alta <- ymd(c("2020-01-12", "2020-03-05", "2020-03-20", "2020-01-28", "2020-01-25"))
event <- c(1, 0, 1, 1, 0)
temps <- as.numeric(alta - ingres)
dades <- data.frame(subjecte = 1:5, ingres, alta, temps, event)
dadesSurv()surv_obj <- Surv(time = temps, event = event)
surv_obj[1] 11 19+ 19 8 15+plot(surv_obj)
Ara ja tenim un objecte Surv preparat per utilitzar-lo en l’estimació de corbes de supervivència amb survfit().
Treballarem amb el conjunt de dades lung, que forma part del paquet survival de R.
Aquest dataset prové d’un assaig clínic realitzat per la Veterans Administration i recull informació de pacients amb càncer de pulmó, incloent:
📦 El dataset està disponible automàticament quan carreguem el paquet survival.
library(survival)
head(lung) # Mostra les primeres filesstatus indica si el pacient ha mort (2) o està censurat (1)sex està codificat com 1 = home, 2 = donalung <- lung %>%
mutate(
status = ifelse(status == 2, 1, 0), # 1 = mort, 0 = censurat
sex = factor(sex, levels = c(1, 2), labels = c("Home", "Dona"))
)
table(lung$status)
0 1
63 165 table(lung$sex)
Home Dona
138 90 Surv()surv_lung <- Surv(time = lung$time, event = lung$status)
head(surv_lung)[1] 306 455 1010+ 210 883 1022+Ara tenim preparades les dades per a estimar la funció de supervivència.
Continuarem amb la corba de Kaplan-Meier i la seva representació visual.
Ara que tenim el nostre objecte Surv() creat, podem estimar la funció de supervivència amb la funció survfit().
km_total <- survfit(surv_lung ~ 1, data = lung)
summary(km_total)$table records n.max n.start events rmean se(rmean) median 0.95LCL
228.00000 228.00000 228.00000 165.00000 376.27475 19.70779 310.00000 285.00000
0.95UCL
363.00000 La taula ens mostra que la mediana de supervivència és de 310 dies (IC 95%: 285–363 dies), que correspon al moment en què la probabilitat de supervivència arriba al 50%.
A la pràctica clínica, però, sovint necessitem respondre a una altra pregunta clau: quina és la probabilitat de supervivència a un temps determinat (p. ex., als 6 mesos o a l’any)? Per extreure aquestes xifres sense haver d’inspeccionar tota la taula d’esdeveniments, utilitzem l’argument times dins de summary():
# Supervivència acumulada estimada als 6 mesos (180 dies) i a 1 any (365 dies)
summary(km_total, times = c(180, 365))Call: survfit(formula = surv_lung ~ 1, data = lung)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
180 160 63 0.722 0.0298 0.666 0.783
365 65 58 0.409 0.0358 0.345 0.486Per representar la corba utilitzarem el paquet survminer, que facilita gràfics clars i publicables.
library(survminer)
ggsurvplot(km_total,
conf.int = TRUE,
xlab = "Dies de seguiment",
ylab = "Probabilitat de supervivència",
title = "Corba de supervivència (Kaplan-Meier)",
surv.median.line = "hv")
La línia de punts mostra la mediana de supervivència estimada.
Podem estimar i visualitzar la supervivència segons el sexe:
km_sexe <- survfit(surv_lung ~ sex, data = lung)
ggsurvplot(km_sexe,
data = lung,
conf.int = TRUE,
pval = TRUE,
risk.table = TRUE,
xlab = "Dies de seguiment",
ylab = "Supervivència estimada",
title = "Supervivència segons el sexe",
legend.labs = c("Home", "Dona"))
Aquest gràfic ens permet observar si hi ha diferències clíniques rellevants i si són estadísticament significatives (valor p del test de log-rank).
Un cop hem representat les corbes de supervivència per grups (per exemple, per sexe), volem saber si la diferència observada és estadísticament significativa.
El test de log-rank (o test de Mantel-Haenszel) és la prova estàndard per comparar dues o més corbes de Kaplan-Meier.
Comprova la hipòtesi nul·la que les funcions de supervivència són iguals entre grups. És a dir:
\[ H_0: S_1(t) = S_2(t) \quad \text{per a tot } t \]
Si el valor p és petit, rebutgem ( H_0 ) → hi ha diferències estadísticament significatives entre grups
survdiff()logrank <- survdiff(surv_lung ~ sex, data = lung)
logrankCall:
survdiff(formula = surv_lung ~ sex, data = lung)
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=Home 138 112 91.6 4.55 10.3
sex=Dona 90 53 73.4 5.68 10.3
Chisq= 10.3 on 1 degrees of freedom, p= 0.001 La sortida ens dona un estadístic de tipus \(\chi^2\) amb els graus de llibertat i el p-valor associat.
El model de Cox ens permet estimar l’efecte de diverses variables (edat, sexe, tractament…) sobre el risc instantani de tenir l’esdeveniment (mort, recidiva, etc.).
És l’equivalent de la regressió logística per a dades de supervivència, però tenint en compte el temps de seguiment.
Modelitza el hazard (risc instantani de morir) en funció de predictors:
\(h(t) = h_0(t) \cdot \exp(\beta_1 X_1 + \beta_2 X_2 + \dots)\)
coxph()cox_model <- coxph(surv_lung ~ sex + age, data = lung)
summary(cox_model)Call:
coxph(formula = surv_lung ~ sex + age, data = lung)
n= 228, number of events= 165
coef exp(coef) se(coef) z Pr(>|z|)
sexDona -0.513219 0.598566 0.167458 -3.065 0.00218 **
age 0.017045 1.017191 0.009223 1.848 0.06459 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
sexDona 0.5986 1.6707 0.4311 0.8311
age 1.0172 0.9831 0.9990 1.0357
Concordance= 0.603 (se = 0.025 )
Likelihood ratio test= 14.12 on 2 df, p=9e-04
Wald test = 13.47 on 2 df, p=0.001
Score (logrank) test = 13.72 on 2 df, p=0.001Els coeficients s’interpreten com a hazard ratios:
Per veure els HR i els seus intervals de confiança:
exp(coef(cox_model)) # HR sexDona age
0.598566 1.017191 exp(confint(cox_model)) # IC 95% dels HR 2.5 % 97.5 %
sexDona 0.4310936 0.8310985
age 0.9989686 1.0357467A partir de la sortida del nostre model (summary(cox_model)), podem interpretar clínicament els coeficients exponenciats:
sexDona, \(\text{HR} = 0.60\), IC 95%: 0.43–0.83, \(p = 0.002\)):Home), un HR de 0.60 indica que les dones tenen un 40% menys de risc instantani de mort que els homes durant el seguiment, mantenint l’edat constant.
exp(-coef) de la sortida d’R), els homes presenten un 67% més de risc de morir que les dones.age, \(\text{HR} = 1.017\), IC 95%: 0.999–1.036, \(p = 0.065\)):📌 El model de Cox és útil per:
El model de Cox assumeix que l’efecte de cada predictor és constant en el temps. És a dir:
El hazard ratio (HR) per a cada variable no canvia al llarg del seguiment.
Això s’anomena l’assumpció de riscos proporcionals.
Si l’HR canvia amb el temps (per exemple, un tractament només és efectiu al principi), el model pot ser inadequat o enganyós.
cox.zph()test_proporcionalitat <- cox.zph(cox_model)
test_proporcionalitat chisq df p
sex 2.608 1 0.11
age 0.209 1 0.65
GLOBAL 2.771 2 0.25plot(test_proporcionalitat)
| Resultat | Què vol dir |
|---|---|
| p > 0.05 | L’assumpció de riscos proporcionals es compleix |
| p < 0.05 | Podria haver-hi un efecte no constant en el temps |
📌 Si es detecta violació, es poden considerar models amb efectes dependents del temps (nivell més avançat) o estratificar per aquella variable.
En aquest capítol hem après a treballar amb dades de temps fins a esdeveniment, i hem vist com analitzar-les i representar-les correctament.
| Eina / concepte | Què permet fer |
|---|---|
Surv() |
Codificar temps i censura per a cada individu |
survfit() + ggsurvplot() |
Estimar i visualitzar corbes de supervivència |
survdiff() |
Comparar grups amb el test de log-rank |
coxph() |
Ajustar un model amb múltiples predictors |
cox.zph() |
Comprovar l’assumpció de riscos proporcionals |
Posa a prova els conceptes apresos sobre censura, estimació no paramètrica de Kaplan-Meier, contrast de log-rank i regressió de Cox amb el teu propi script o document Quarto[cite: 1]:
Surv manual: Un assaig clínic segueix 6 pacients després d’una intervenció coronària percussió (temps en mesos):
Surv() i comprova quins casos queden marcats amb censura a la dreta (+)[cite: 1].lung del paquet survival[cite: 1]:
status de manera que l’esdeveniment (mort) sigui 1 i la censura sigui 0 (recorda que originàriament ve codificat com 1=censura, 2=èxitus)[cite: 1].sex en un factor amb les etiquetes "Home" i "Dona"[cite: 1].lung[cite: 1].ph.ecog mesura l’escala funcional de l’ECOG (0 = asimptomàtic fins a 3 = enllitat la major part del dia)[cite: 1].
ph.ecog[cite: 1].ph.ecog incorporant la taula de pacients a risc amb survminer::ggsurvplot()[cite: 1].survdiff) per contrastar si les diferències de supervivència entre graus d’ECOG són estadísticament significatives[cite: 1].coxph) per avaluar la mortalitat en funció de l’edat (age), el sexe (sex) i l’estat funcional (ph.ecog)[cite: 1].ph.ecog mantenint constants l’edat i el sexe[cite: 1].cox.zph()[cite: 1].library(tidyverse)
library(survival)
library(survminer)
# -------------------------------------------------------------
# Exercici 1: Codificació de censura i objecte Surv manual
# -------------------------------------------------------------
cohort_manual <- tibble(
pacient = 1:6,
temps_mesos = c(6, 18, 24, 12, 3, 24),
esdeveniment = c(1, 0, 0, 1, 1, 0)
)
surv_manual <- Surv(time = cohort_manual$temps_mesos, event = cohort_manual$esdeveniment)
surv_manual[1] 6 18+ 24+ 12 3 24+# -------------------------------------------------------------
# Exercici 2: Carrega neta de dades 'lung'
# -------------------------------------------------------------
# Carreguem directament les dades originals del paquet survival per evitar col·lisions
df_lung_raw <- as_tibble(survival::lung)
lung_exercicis <- df_lung_raw %>%
mutate(
status_net = ifelse(status == 2, 1, 0), # 1 = mort, 0 = censurat
sexe = factor(sex, levels = c(1, 2), labels = c("Home", "Dona")),
ecog = factor(ph.ecog)
)
# Recompte d'esdeveniments i censures
table(Estat = lung_exercicis$status_net)Estat
0 1
63 165 table(Sexe = lung_exercicis$sexe, Estat = lung_exercicis$status_net) Estat
Sexe 0 1
Home 26 112
Dona 37 53# -------------------------------------------------------------
# Exercici 3: Kaplan-Meier global i probabilitat a 1 any (365 dies)
# -------------------------------------------------------------
surv_obj_ex <- Surv(time = lung_exercicis$time, event = lung_exercicis$status_net)
km_global <- survfit(surv_obj_ex ~ 1, data = lung_exercicis)
summary(km_global)$table # Mediana i IC 95% records n.max n.start events rmean se(rmean) median 0.95LCL
228.00000 228.00000 228.00000 165.00000 376.27475 19.70779 310.00000 285.00000
0.95UCL
363.00000 # Estimació puntual de supervivència als 365 dies
summary(km_global, times = 365)Call: survfit(formula = surv_obj_ex ~ 1, data = lung_exercicis)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
365 65 121 0.409 0.0358 0.345 0.486# -------------------------------------------------------------
# Exercici 4: Corbes per ECOG i contrast de Log-Rank
# -------------------------------------------------------------
lung_ecog_valids <- lung_exercicis %>%
filter(!is.na(ecog))
km_ecog <- survfit(Surv(time, status_net) ~ ecog, data = lung_ecog_valids)
ggsurvplot(
km_ecog,
data = lung_ecog_valids,
pval = TRUE,
risk.table = TRUE,
conf.int = FALSE,
title = "Supervivència segons l'estatus ECOG (ph.ecog)",
xlab = "Dies de seguiment",
ylab = "Supervivència acumulada",
risk.table.height = 0.28
)
# Test de log-rank
survdiff(Surv(time, status_net) ~ ecog, data = lung_ecog_valids)Call:
survdiff(formula = Surv(time, status_net) ~ ecog, data = lung_ecog_valids)
N Observed Expected (O-E)^2/E (O-E)^2/V
ecog=0 63 37 54.153 5.4331 8.2119
ecog=1 113 82 83.528 0.0279 0.0573
ecog=2 50 44 26.147 12.1893 14.6491
ecog=3 1 1 0.172 3.9733 4.0040
Chisq= 22 on 3 degrees of freedom, p= 7e-05 # -------------------------------------------------------------
# Exercici 5: Model multivariable de Cox
# -------------------------------------------------------------
mod_cox <- coxph(Surv(time, status_net) ~ age + sexe + ecog, data = lung_exercicis)
summary(mod_cox)Call:
coxph(formula = Surv(time, status_net) ~ age + sexe + ecog, data = lung_exercicis)
n= 227, number of events= 164
(1 observation deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.010795 1.010853 0.009312 1.159 0.24637
sexeDona -0.545831 0.579360 0.168228 -3.245 0.00118 **
ecog1 0.410048 1.506890 0.199604 2.054 0.03995 *
ecog2 0.903303 2.467741 0.228078 3.960 7.48e-05 ***
ecog3 1.954543 7.060694 1.029701 1.898 0.05767 .
---
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.9926 1.0295
sexeDona 0.5794 1.7260 0.4166 0.8056
ecog1 1.5069 0.6636 1.0190 2.2284
ecog2 2.4677 0.4052 1.5782 3.8587
ecog3 7.0607 0.1416 0.9383 53.1289
Concordance= 0.637 (se = 0.025 )
Likelihood ratio test= 30.87 on 5 df, p=1e-05
Wald test = 31.43 on 5 df, p=8e-06
Score (logrank) test = 33.56 on 5 df, p=3e-06# Hazard Ratios (HR) i intervals de confiança del 95%
exp(cbind(HR = coef(mod_cox), confint(mod_cox))) HR 2.5 % 97.5 %
age 1.0108531 0.9925710 1.0294721
sexeDona 0.5793604 0.4166324 0.8056467
ecog1 1.5068901 1.0190093 2.2283583
ecog2 2.4677406 1.5781881 3.8586933
ecog3 7.0606942 0.9383484 53.1288860# -------------------------------------------------------------
# Exercici 6: Diagnòstic de riscos proporcionals (Schoenfeld)
# -------------------------------------------------------------
test_prop <- cox.zph(mod_cox)
print(test_prop) chisq df p
age 0.163 1 0.69
sexe 2.406 1 0.12
ecog 5.650 3 0.13
GLOBAL 7.777 5 0.17# Diagnòstic visual
par(mfrow = c(2, 2))
plot(test_prop)
par(mfrow = c(1, 1))