Análisis de Supervivencia

El objetivo de este RMarkDown es mostrar como se puede realizar análisis con el ejemplo del capitulo 27 con The Epidemiologist R Handbook

https://www.epirhandbook.com/en/

27.1 Descripción general

El análisis de supervivencia se centra en describir, para un individuo o grupo de individuos, un punto definido de evento denominado fracaso (aparición de una enfermedad, curación, fallecimiento, recaída tras la respuesta al tratamiento, etc.), que ocurre tras un período denominado tiempo de fracaso (o tiempo de seguimiento en estudios de cohortes/poblacionales) durante el cual se observa a los individuos. Para determinar el tiempo de fracaso, es necesario definir un momento de origen (que puede ser la fecha de inclusión, la fecha de diagnóstico, etc.).

El objetivo de la inferencia para el análisis de supervivencia es, por lo tanto, el tiempo transcurrido entre un origen y un evento. En la investigación médica actual, se utiliza ampliamente en estudios clínicos para evaluar el efecto de un tratamiento, por ejemplo, o en epidemiología del cáncer para evaluar una amplia variedad de medidas de supervivencia.

Generalmente se expresa a través de la probabilidad de supervivencia , que es la probabilidad de que el evento de interés no haya ocurrido en una duración t.

Censura : La censura ocurre cuando, al final del seguimiento, algunos individuos no han presentado el evento de interés y, por lo tanto, se desconoce el tiempo real hasta el evento. Nos centraremos principalmente en la censura correcta, pero para más detalles sobre la censura y el análisis de supervivencia en general, pueden consultar las referencias.

27.2 Preparación Cargar paquetes Para ejecutar análisis de supervivencia en R, uno de los paquetes más utilizados es survival . Primero lo instalamos y luego lo cargamos, junto con los demás paquetes que se usarán en esta sección:

En este manual, destacamos p_load()el uso de pacman , que instala el paquete si es necesario y lo carga para su uso. También puede cargar paquetes instalados library()desde R base. Consulte la página sobre conceptos básicos de R para obtener más información sobre los paquetes de R. https://www.epirhandbook.com/en/new_pages/basics.html

Paquetes y librerias necesarias

Paquetes usados:

  • pacman: gestor que instala y carga paquetes en un solo paso.
  • survival: funciones básicas de supervivencia (Surv, survfit, coxph…).
  • rio: importación simple de datos (import de .rds, .csv…).
  • dplyr: manipulación de datos con tuberías (%>%).
  • forcats: manejo y reordenamiento de factores.
  • janitor: limpieza de tablas y tabyl() para frecuencias.
  • tibble: creación y manipulación de tibbles; rowid_to_column().
  • survminer: gráficos de supervivencia y diagnóstico de modelos.
  • ggplot2: esquema de graficación base para survminer.
  • RColorBrewer: paletas de color para gráficas.
  • SemiCompRisks: dataset BMT de trasplante de médula ósea.
# Instalar pacman si no está instalado
if (!require("pacman")) install.packages("pacman")
## Cargando paquete requerido: pacman
# Cargar (y si no están, instalar) todos los paquetes necesarios de una sola llamada
pacman::p_load(
  survival,
  rio,
  dplyr,
  forcats,
  janitor,
  tibble,
  survminer,
  ggplot2,
  RColorBrewer,
  SemiCompRisks,
  tidyverse,
  modelsummary,
  readxl,
  openxlsx,
  plotly,
  car,
  knitr,
  DT
)
## Installing package into 'C:/Users/crist/AppData/Local/R/win-library/4.5'
## (as 'lib' is unspecified)
## also installing the dependency 'snakecase'
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5:
##   no fue posible abrir la URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5/PACKAGES'
## package 'snakecase' successfully unpacked and MD5 sums checked
## package 'janitor' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\crist\AppData\Local\Temp\RtmpSIRblS\downloaded_packages
## 
## janitor installed
## Installing package into 'C:/Users/crist/AppData/Local/R/win-library/4.5'
## (as 'lib' is unspecified)
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5:
##   no fue posible abrir la URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5/PACKAGES'
## package 'SemiCompRisks' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\crist\AppData\Local\Temp\RtmpSIRblS\downloaded_packages
## 
## SemiCompRisks installed
## Installing package into 'C:/Users/crist/AppData/Local/R/win-library/4.5'
## (as 'lib' is unspecified)
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5:
##   no fue posible abrir la URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5/PACKAGES'
## package 'plotly' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\crist\AppData\Local\Temp\RtmpSIRblS\downloaded_packages
## 
## plotly installed
## Installing package into 'C:/Users/crist/AppData/Local/R/win-library/4.5'
## (as 'lib' is unspecified)
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5:
##   no fue posible abrir la URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.5/PACKAGES'
## package 'DT' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\crist\AppData\Local\Temp\RtmpSIRblS\downloaded_packages
## 
## DT installed

Forma controlada de cargar librerias y dependencias

Importar conjunto de datos

Importamos el conjunto de datos de casos de una epidemia simulada de ébola. Si desea seguir los pasos, haga clic para descargar la lista de líneas «limpia» (como archivo .rds). Importe los datos con la función import()del paquete rio (admite muchos tipos de archivos, como .xlsx, .csv, .rds; consulte la página Importar y exportar para obtener más información).

https://github.com/appliedepi/epirhandbook/blob/master/inst/extdata/case_linelists/linelist_cleaned.rds

https://www.epirhandbook.com/en/new_pages/importing.html

Variable Descripción
case_id Identificador único de cada caso en el linelist; permite enlazar registros de un mismo individuo. (epirhandbook.com)
generation Generación de transmisión: número de saltos desde el caso índice hasta el caso actual. (epirhandbook.com)
date_infection Fecha (estimada o conocida) en que ocurrió la infección. (epirhandbook.com)
date_onset Fecha de inicio de síntomas. (epirhandbook.com)
date_hospitalisation Fecha de ingreso hospitalario. (epirhandbook.com)
date_outcome Fecha en que se registró el desenlace final (recuperación, muerte, etc.). (epirhandbook.com)
outcome Resultado observado: categoría del desenlace (p. ej., “Death”, “Recover”). (epirhandbook.com)
gender Sexo del caso (“m” para masculino, “f” para femenino). (epirhandbook.com)
age Edad del caso en la unidad original (número). (epirhandbook.com)
age_unit Unidad en la que se midió la edad (p. ej., “years”, “months”). (epirhandbook.com)
age_years Edad convertida a años (numérico). (epirhandbook.com)
age_cat Categoría de edad original (rangos definidos en el linelist). (epirhandbook.com)
age_cat5 Categorías de edad en intervalos de 5 años. (epirhandbook.com)
hospital Nombre del hospital o centro donde fue ingresado el caso. (epirhandbook.com)
lon Longitud geográfica de residencia o ubicación del caso. (epirhandbook.com)
lat Latitud geográfica de residencia o ubicación del caso. (epirhandbook.com)
infector Identificador de la persona que infectó a este caso (vinculado a case_id). (epirhandbook.com)
source Fuente o tipo de exposición que causó la infección (p. ej., “community”, “nosocomial”). (epirhandbook.com)
wt_kg Peso del caso en kilogramos. (epirhandbook.com)
ht_cm Talla (altura) del caso en centímetros. (epirhandbook.com)
ct_blood Valor de ciclo umbral (Ct) en prueba PCR de sangre. (epirhandbook.com)
fever Indica si el caso presentó fiebre (“yes”/“no”). (epirhandbook.com)
chills Indica si el caso presentó escalofríos. (epirhandbook.com)
cough Indica si el caso presentó tos. (epirhandbook.com)
aches Indica si el caso presentó dolores o molestias. (epirhandbook.com)
vomit Indica si el caso presentó vómito. (epirhandbook.com)
temp Temperatura corporal registrada (°C). (epirhandbook.com)
time_admission Hora o marca temporal del ingreso al hospital. (epirhandbook.com)
bmi Índice de masa corporal: peso (kg) dividido por talla al cuadrado (m²). (epirhandbook.com)
days_onset_hosp Días transcurridos entre inicio de síntomas y hospitalización. (epirhandbook.com)
event Indicador de evento para análisis de supervivencia: 1 = muerte (evento observado), 0 = censurado (recuperación o desenlace desconocido). (epirhandbook.com)
futime Tiempo de seguimiento en días: diferencia entre date_outcome y date_onset. (epirhandbook.com)
age_cat_small Categorías simplificadas de edad (“0–4”, “5–19”, “20+”) para análisis estratificado. (epirhandbook.com)

El objeto linelist_surv proviene de un conjunto de datos de un brote simulado de Ébola, utilizado en el capítulo de análisis de supervivencia del Epidemiologist R Handbook.

En este análisis, el evento de interés es la muerte por Ébola, y los pacientes que se recuperan o cuyo desenlace no se conoce se consideran censurados (event = 0). Además, la variable de tiempo de seguimiento (futime) se calcula como los días transcurridos entre el inicio de los síntomas y la fecha del desenlace (muerte o recuperación), permitiendo así estimar la probabilidad de supervivencia tras la infección por Ébola.

# Toma como ruta setwd la misma ruta en la que se guarde el archivo.RMD
#path <- dirname(rstudioapi::getActiveDocumentContext()$path)
# importar linelist
linelist_case_data <- rio::import("linelist_cleaned.rds")

Gestión y transformación de datos

En análisis de supervivencia, los datos deben estructurarse en torno a tres características:

  1. Tiempo hasta evento: la variable respuesta es el número de días hasta que ocurre un evento definido.
  2. Censura: algunas observaciones no han registrado el evento al final del seguimiento.
  3. Predictores: variables explicativas cuyo efecto sobre el tiempo de espera queremos evaluar.

Por ello, a partir de linelist_case_data crearemos:

  • Un nuevo data frame linelist_surv.
  • La variable de evento (event), codificando muerte como 1 y todo lo demás (recuperados o desconocidos) como 0.
  • El tiempo de seguimiento (futime), en días, como la diferencia entre fecha de inicio y fecha de resultado.
  • Eliminaremos filas con fechas faltantes o inconsistentes (inicio posterior a resultado).

PRECAUCIÓN: Dado que en un estudio de cohorte real se conoce la información sobre el momento de origen y la finalización del seguimiento, dado que se observa a los individuos, eliminaremos las observaciones cuya fecha de inicio o de resultado se desconozca. Asimismo, se eliminarán los casos en que la fecha de inicio sea posterior a la fecha de resultado, ya que se consideran incorrectos.

SUGERENCIA: Dado que filtrar por valores mayores que (>) o menores que (<) una fecha puede eliminar filas con valores faltantes, aplicar el filtro en las fechas incorrectas también eliminará las filas con fechas faltantes.

  • Finalmente, definiremos age_cat_small con tres categorías de edad usando case_when().
#crea linelist_surv a partir de linelist_case_data

linelist_surv <-  linelist_case_data %>% 
     
  dplyr::filter(
       # eliminar observaciones con fechas de inicio o de resultado erróneas o ausentes
       date_outcome > date_onset) %>% 
  
  dplyr::mutate(
       # crear la var evento que es 1 si el paciente falleció y 0 si fue correctamente censurado
       event = ifelse(is.na(outcome) | outcome == "Recover", 0, 1), 
    
       # crear la var de tiempo de seguimiento en días
       futime = as.double(date_outcome - date_onset), 
    
       # crear una nueva variable de categoría de edad con sólo 3 estratos
       age_cat_small = dplyr::case_when( 
            age_years < 5  ~ "0-4",
            age_years >= 5 & age_years < 20 ~ "5-19",
            age_years >= 20   ~ "20+"),
       
       # el paso anterior creó la variable age_cat_small como carácter.
       # ahora se convierte en factor y se especifican los niveles.
       # Nótese que los valores NA siguen siendo NA y no se ponen en un nivel "desconocido" por ejemplo,
       # ya que en los siguientes análisis tienen que ser eliminados.
       age_cat_small = fct_relevel(age_cat_small, "0-4", "5-19", "20+")
       )

CONSEJO: Podemos verificar las nuevas columnas que hemos creado haciendo un resumen en futimey una tabulación cruzada entre eventy outcomea partir de la cual se creó. Además de esta verificación, es una buena costumbre comunicar el tiempo medio de seguimiento al interpretar los resultados del análisis de supervivencia.

summary(linelist_surv$futime)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    6.00   10.00   11.98   16.00   64.00

Tabulación cruzada entre outcome (variable original) y event (nueva variable) para verificar que event se creó correctamente según lo esperado.

# cruzar la nueva variable de evento y la de resultado a partir de la que se creó
# para asegurarse de que el código ha hecho lo que debía
linelist_surv %>% 
  tabyl(outcome, event)
##  outcome    0    1
##    Death    0 1952
##  Recover 1547    0
##     <NA> 1040    0

Ahora cruzamos la nueva variable age_cat_small y la antigua columna age_cat para asegurarnos de que las asignaciones son correctas.

linelist_surv %>% 
  tabyl(age_cat_small, age_cat)
##  age_cat_small 0-4 5-9 10-14 15-19 20-29 30-49 50-69 70+ NA_
##            0-4 834   0     0     0     0     0     0   0   0
##           5-19   0 852   717   575     0     0     0   0   0
##            20+   0   0     0     0   862   554    69   5   0
##           <NA>   0   0     0     0     0     0     0   0  71

Ahora revisamos las 10 primeras observaciones de los datos linelist_surv, analizando variables específicas (incluidas las de nueva creación).

# Carga explícita de dplyr (y magrittr para %>% si no lo trae)
library(dplyr)

linelist_surv %>%
  # asegurar que usamos select de dplyr
  dplyr::select(
    case_id,
    age_cat_small,
    date_onset,
    date_outcome,
    outcome,
    event,
    futime
  ) %>%
  head(10)
##    case_id age_cat_small date_onset date_outcome outcome event futime
## 1   8689b7           0-4 2014-05-13   2014-05-18 Recover     0      5
## 2   11f8ea           20+ 2014-05-16   2014-05-30 Recover     0     14
## 3   893f25           0-4 2014-05-21   2014-05-29 Recover     0      8
## 4   be99c8          5-19 2014-05-22   2014-05-24 Recover     0      2
## 5   07e3e8          5-19 2014-05-27   2014-06-01 Recover     0      5
## 6   369449           0-4 2014-06-02   2014-06-07   Death     1      5
## 7   f393b4           20+ 2014-06-05   2014-06-18 Recover     0     13
## 8   1389ca           20+ 2014-06-05   2014-06-09   Death     1      4
## 9   2978ac          5-19 2014-06-06   2014-06-15   Death     1      9
## 10  fc15ef          5-19 2014-06-16   2014-07-09 Recover     0     23

También podemos tabular las columnas age_cat_smally genderpara obtener más detalles sobre la distribución de esta nueva columna por género. Utilizamos tabyl()y las funciones adorn de janitor, tal y como se describe en la página Tablas descriptivas.

linelist_surv %>% 
  tabyl(gender, age_cat_small, show_na = F) %>% 
  adorn_totals(where = "both") %>% 
  adorn_percentages() %>% 
  adorn_pct_formatting() %>% 
  adorn_ns(position = "front")
##  gender         0-4          5-19           20+          Total
##       f 482 (22.4%) 1,184 (54.9%)   490 (22.7%) 2,156 (100.0%)
##       m 325 (15.0%)   880 (40.6%)   960 (44.3%) 2,165 (100.0%)
##   Total 807 (18.7%) 2,064 (47.8%) 1,450 (33.6%) 4,321 (100.0%)
names(linelist_surv)
##  [1] "case_id"              "generation"           "date_infection"      
##  [4] "date_onset"           "date_hospitalisation" "date_outcome"        
##  [7] "outcome"              "gender"               "age"                 
## [10] "age_unit"             "age_years"            "age_cat"             
## [13] "age_cat5"             "hospital"             "lon"                 
## [16] "lat"                  "infector"             "source"              
## [19] "wt_kg"                "ht_cm"                "ct_blood"            
## [22] "fever"                "chills"               "cough"               
## [25] "aches"                "vomit"                "temp"                
## [28] "time_admission"       "bmi"                  "days_onset_hosp"     
## [31] "event"                "futime"               "age_cat_small"
table(linelist_surv$hospital)
## 
##                     Central Hospital                    Military Hospital 
##                                  352                                  689 
##                              Missing                                Other 
##                                 1130                                  658 
##                        Port Hospital St. Mark's Maternity Hospital (SMMH) 
##                                 1396                                  314
summary(linelist_surv$ct_blood)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.00   20.00   22.00   21.18   22.00   26.00
linelist_surv %>%
  tbl_summary(statistic = list(
      all_continuous() ~ "{mean} ± {sd}"), include = c(gender,outcome,wt_kg, wt_kg,age_years,ht_cm,ct_blood,fever,cough, chills, aches,temp,bmi,event, futime),
              by= gender,
              label = list(outcome ~ "Desenlace",
                           wt_kg ~ "Peso",
                           ht_cm ~ "Talla",
                           ct_blood ~ "PCR dE sangre",
                           fever ~ "Fiebre",
                           cough ~ "Tos",
                           chills ~ "Escalofríos",
                           aches ~ "Dolores o Molestias",
                           temp ~ "Temperatura Comporal",
                           bmi ~ "IMC",
                           event ~ "Muerto",
                           age_years  ~ "Edad",
                           futime ~ "Tiempo de Supervivencia"),
              missing_text = "Sin dato"
              )  %>%
  modify_header(label ~ "**Variable**")
## 218 missing rows in the "gender" column have been removed.
Variable f
N = 2,1561		 m
N = 2,1651
Desenlace

    Death 924 (55%) 929 (56%)
    Recover 743 (45%) 738 (44%)
    Sin dato 489 498
Peso 46 ± 17 59 ± 18
Edad 12 ± 9 19 ± 14
Talla 109 ± 45 141 ± 48
PCR dE sangre 21.16 ± 1.70 21.15 ± 1.74
Fiebre 1,666 (81%) 1,685 (81%)
    Sin dato 91 84
Tos 1,775 (86%) 1,794 (86%)
    Sin dato 91 84
Escalofríos 406 (20%) 410 (20%)
    Sin dato 91 84
Dolores o Molestias 184 (8.9%) 203 (9.8%)
    Sin dato 91 84
Temperatura Comporal 38.58 ± 0.98 38.55 ± 0.97
    Sin dato 59 58
IMC 54 ± 69 38 ± 28
Muerto 924 (43%) 929 (43%)
Tiempo de Supervivencia 12 ± 9 12 ± 9
1 n (%); Mean ± SD

27.3 Conceptos básicos del análisis de supervivencia

Creación de un objeto de surv-type

Primero utilizaremos Surv()de survival para construir un objeto de supervivencia a partir de las columnas de tiempo de seguimiento y eventos. El resultado de este paso es un objeto de tipo Surv que condensa la información temporal y si se ha observado el evento de interés (muerte). Este objeto se utilizará finalmente en el lado derecho de las fórmulas del modelo posteriores (véase la documentación).

https://cran.r-project.org/web/packages/survival/vignettes/survival.pdf

# Utilizar la sintaxis Suv() para datos correctamente censurados
survobj <- Surv(time = linelist_surv$futime,
                event = linelist_surv$event)

A modo de resumen, aquí están las primeras 10 filas de los datos linelist_surv, mostrando solo algunas columnas importantes.

linelist_surv %>%
  # asegurar que usamos select de dplyr
  dplyr::select(
   case_id, 
   date_onset, 
   date_outcome, 
   futime, 
   outcome, 
   event
  ) %>%
  head(10)
##    case_id date_onset date_outcome futime outcome event
## 1   8689b7 2014-05-13   2014-05-18      5 Recover     0
## 2   11f8ea 2014-05-16   2014-05-30     14 Recover     0
## 3   893f25 2014-05-21   2014-05-29      8 Recover     0
## 4   be99c8 2014-05-22   2014-05-24      2 Recover     0
## 5   07e3e8 2014-05-27   2014-06-01      5 Recover     0
## 6   369449 2014-06-02   2014-06-07      5   Death     1
## 7   f393b4 2014-06-05   2014-06-18     13 Recover     0
## 8   1389ca 2014-06-05   2014-06-09      4   Death     1
## 9   2978ac 2014-06-06   2014-06-15      9   Death     1
## 10  fc15ef 2014-06-16   2014-07-09     23 Recover     0

Y aquí están los primeros 10 elementos de survobj. Se imprime esencialmente como un vector de tiempo de seguimiento, con «+» para representar si una observación fue censurada por la derecha. Observe cómo se alinean los números arriba y abajo.

#print the 50 first elements of the vector to see how it presents
head(survobj, 10)
##  [1]  5+ 14+  8+  2+  5+  5  13+  4   9  23+

Ejecución de análisis iniciales

A continuación, iniciamos nuestro análisis utilizando la survfit()función para generar un objeto de ajuste de supervivencia (survfit) , que se ajusta a los cálculos predeterminados de las estimaciones de Kaplan-Meier (KM) de la curva de supervivencia global (marginal), que, de hecho, son una función escalonada con saltos en los tiempos de los eventos observados. El objeto de ajuste de supervivencia final contiene una o más curvas de supervivencia y se crea utilizando el objeto Surv como variable de respuesta en la fórmula del modelo.

NOTA: La estimación de Kaplan-Meier es una estimación de máxima verosimilitud (MLE) no paramétrica de la función de supervivencia. (consulte los recursos para obtener más información).

El resumen de este objeto de ajuste de supervivencia generará lo que se denomina una tabla de vida . Para cada intervalo de tiempo del seguimiento ( time) donde ocurrió un evento (en orden ascendente):

Ajustamos las estimaciones de KM utilizando la fórmula donde el objeto de supervivencia anterior “survobj” es la variable de respuesta. “~ 1” precisa que ejecutamos el modelo para la supervivencia general.

# ajustar las estimaciones de KM utilizando una fórmula donde el objeto Surv "survobj" es la variable de respuesta.
# "~ 1" significa que ejecutamos el modelo para la supervivencia global  
linelistsurv_fit <-  survival::survfit(survobj ~ 1)

#imprimir su resumen para más detalles
summary(linelistsurv_fit)
## Call: survfit(formula = survobj ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1   4539      30    0.993 0.00120        0.991        0.996
##     2   4500      69    0.978 0.00217        0.974        0.982
##     3   4394     149    0.945 0.00340        0.938        0.952
##     4   4176     194    0.901 0.00447        0.892        0.910
##     5   3899     214    0.852 0.00535        0.841        0.862
##     6   3592     210    0.802 0.00604        0.790        0.814
##     7   3223     179    0.757 0.00656        0.745        0.770
##     8   2899     167    0.714 0.00700        0.700        0.728
##     9   2593     145    0.674 0.00735        0.660        0.688
##    10   2311     109    0.642 0.00761        0.627        0.657
##    11   2081     119    0.605 0.00788        0.590        0.621
##    12   1843      89    0.576 0.00809        0.560        0.592
##    13   1608      55    0.556 0.00823        0.540        0.573
##    14   1448      43    0.540 0.00837        0.524        0.556
##    15   1296      31    0.527 0.00848        0.511        0.544
##    16   1152      48    0.505 0.00870        0.488        0.522
##    17   1002      29    0.490 0.00886        0.473        0.508
##    18    898      21    0.479 0.00900        0.462        0.497
##    19    798       7    0.475 0.00906        0.457        0.493
##    20    705       4    0.472 0.00911        0.454        0.490
##    21    626      13    0.462 0.00932        0.444        0.481
##    22    546       8    0.455 0.00948        0.437        0.474
##    23    481       5    0.451 0.00962        0.432        0.470
##    24    436       4    0.447 0.00975        0.428        0.466
##    25    378       4    0.442 0.00993        0.423        0.462
##    26    336       3    0.438 0.01010        0.419        0.458
##    27    297       1    0.436 0.01017        0.417        0.457
##    29    235       1    0.435 0.01030        0.415        0.455
##    38     73       1    0.429 0.01175        0.406        0.452
#imprimir los percentiles 25, 50 y 75
qs <- quantile(linelistsurv_fit, probs = c(0.0, 0.05, 0.25, 0.50, 0.75, 0.95, 1.0))
qs
## $quantile
##   0   5  25  50  75  95 100 
##   0   3   8  17  NA  NA  NA 
## 
## $lower
##   0   5  25  50  75  95 100 
##   0   3   7  16  NA  NA  NA 
## 
## $upper
##   0   5  25  50  75  95 100 
##   0   4   8  18  NA  NA  NA

Al utilizar summary(), podemos añadir la opción timesy especificar determinados momentos en los que queremos ver la información de supervivencia.

#imprimir su resumen en momentos concretos
summary(linelistsurv_fit, times = c(5,10,20,30,60))
## Call: survfit(formula = survobj ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5   3899     656    0.852 0.00535        0.841        0.862
##    10   2311     810    0.642 0.00761        0.627        0.657
##    20    705     446    0.472 0.00911        0.454        0.490
##    30    210      39    0.435 0.01030        0.415        0.455
##    60      2       1    0.429 0.01175        0.406        0.452

Se puede tavla más bonita la tabla Flextable.

# 1. Obtener el resumen en los tiempos de interés
surv_sum <- summary(linelistsurv_fit, times = c(5, 10, 15, 20, 25, 30, 35, 40, 45,50,55,60))

# 2. Extraer la tabla de resultados y convertir a data.frame
tab <- with(
  surv_sum,
  data.frame(
    time      = time,
    n.risk    = n.risk,
    n.event   = n.event,
    surv      = round(surv, 3),
    std.error = round(std.err, 3),
    lower     = round(lower, 3),
    upper     = round(upper, 3)
  )
)

# 3. Crear la flextable
library(flextable)

ft <- flextable(tab) %>%
  set_header_labels(
    time      = "Tiempo (días)",
    n.risk    = "En riesgo",
    n.event   = "Eventos",
    surv      = "Supervivencia",
    std.error = "Error estándar",
    lower     = "IC 95% (inf.)",
    upper     = "IC 95% (sup.)"
  ) %>%
  colformat_num(
    j = c("surv","std.error","lower","upper"),
    digits = 3
  ) %>%
  theme_zebra(
    odd_header = "transparent", 
    odd_body   = "#f9f9f9"
  ) %>%
  bold(part = "header") %>%
  autofit()

# 4. Mostrar la tabla
ft

Tiempo (días)

En riesgo

Eventos

Supervivencia

Error estándar

IC 95% (inf.)

IC 95% (sup.)

5

3,899

656

0.852

0.005

0.841

0.862

10

2,311

810

0.642

0.008

0.627

0.657

15

1,296

337

0.527

0.008

0.511

0.544

20

705

109

0.472

0.009

0.454

0.490

25

378

34

0.442

0.010

0.423

0.462

30

210

5

0.435

0.010

0.415

0.455

35

114

0

0.435

0.010

0.415

0.455

40

61

1

0.429

0.012

0.406

0.452

45

36

0

0.429

0.012

0.406

0.452

50

18

0

0.429

0.012

0.406

0.452

55

9

0

0.429

0.012

0.406

0.452

60

2

0

0.429

0.012

0.406

0.452

También podemos usar la print()función. El print.rmean = TRUEargumento se utiliza para obtener el tiempo medio de supervivencia y su error estándar (EE).

NOTA: El tiempo medio de supervivencia restringido (RMST) es una medida de supervivencia específica cada vez más utilizada en el análisis de supervivencia del cáncer y que a menudo se define como el área bajo la curva de supervivencia, dado que observamos a los pacientes hasta el tiempo restringido T (más detalles en la sección Recursos).

# imprimir el objeto linelistsurv_fit con el tiempo medio de supervivencia y su se. 
print(linelistsurv_fit, print.rmean = TRUE)
## Call: survfit(formula = survobj ~ 1)
## 
##         n events rmean* se(rmean) median 0.95LCL 0.95UCL
## [1,] 4539   1952   33.1     0.539     17      16      18
##     * restricted mean with upper limit =  64

SUGERENCIA: Podemos crear el objeto surv directamente en lasurvfit()función y guardar una línea de código. Se verá asílinelistsurv_quick <- survfit(Surv(futime, event) ~ 1, data=linelist_surv):

Riesgo acumulado

Además de la función summary(), también podemos utilizar la función str() que da más detalles sobre la estructura del objeto survfit(). Se trata de una lista de 16 elementos.

Entre estos elementos hay uno importante: el cumhaz, que es un vector numérico. Se puede trazar para mostrar el riesgo acumulado, siendo el riesgo la tasa instantánea de ocurrencia del evento (ver referencias).

str(linelistsurv_fit)
## List of 17
##  $ n        : int 4539
##  $ time     : num [1:59] 1 2 3 4 5 6 7 8 9 10 ...
##  $ n.risk   : num [1:59] 4539 4500 4394 4176 3899 ...
##  $ n.event  : num [1:59] 30 69 149 194 214 210 179 167 145 109 ...
##  $ n.censor : num [1:59] 9 37 69 83 93 159 145 139 137 121 ...
##  $ surv     : num [1:59] 0.993 0.978 0.945 0.901 0.852 ...
##  $ std.err  : num [1:59] 0.00121 0.00222 0.00359 0.00496 0.00628 ...
##  $ cumhaz   : num [1:59] 0.00661 0.02194 0.05585 0.10231 0.15719 ...
##  $ std.chaz : num [1:59] 0.00121 0.00221 0.00355 0.00487 0.00615 ...
##  $ type     : chr "right"
##  $ logse    : logi TRUE
##  $ conf.int : num 0.95
##  $ conf.type: chr "log"
##  $ lower    : num [1:59] 0.991 0.974 0.938 0.892 0.841 ...
##  $ upper    : num [1:59] 0.996 0.982 0.952 0.91 0.862 ...
##  $ t0       : num 0
##  $ call     : language survfit(formula = survobj ~ 1)
##  - attr(*, "class")= chr "survfit"

Representar curvas de Kaplan-Meir

Una vez ajustadas las estimaciones de KM, podemos visualizar la probabilidad de estar vivo a lo largo de un tiempo determinado utilizando la función básica plot() que dibuja la “curva de Kaplan-Meier”. En otras palabras, la curva de abajo es una ilustración convencional de la experiencia de supervivencia en todo el grupo de pacientes.

Podemos verificar rápidamente el tiempo de seguimiento mínimo y máximo en la curva.

Una forma fácil de interpretarlo es decir que en el momento cero, todos los participantes están vivos y la probabilidad de supervivencia es entonces del 100%. Esta probabilidad disminuye con el tiempo a medida que los pacientes mueren. La proporción de participantes que sobreviven más allá de los 60 días de seguimiento se sitúa en torno al 40%.

plot(linelistsurv_fit, 
     xlab = "Days of follow-up",    # etiqueta del eje-x
     ylab="Survival Probability",   # etiqueta del eje-y
     main= "Overall survival curve" # título de la figura
     )

El intervalo de confianza de las estimaciones de supervivencia de KM también se representa por defecto y puede descartarse añadiendo la opción conf.int = FALSE al comando plot().

Dado que el evento de interés es “death”, dibujar una curva que describa los complementos de las proporciones de supervivencia llevará a dibujar las proporciones de mortalidad acumulada. Esto puede hacerse con lines(), que añade información a un gráfico existente.

# gráfico original
plot(
  linelistsurv_fit,
  xlab = "Days of follow-up",       
  ylab = "Survival Probability",       
  mark.time = TRUE,              # marca los eventos en la curva: se imprime un "+" en cada evento
  conf.int = FALSE,              # no representa el intervalo de confianza
  main = "Overall survival curve and cumulative mortality"
  )

# dibujar una curva adicional a la anterior
lines(
  linelistsurv_fit,
  lty = 3,             # use different line type for clarity
  fun = "event",       # dibuja los eventos acumulativos en lugar de la supervivencia 
  mark.time = FALSE,
  conf.int = FALSE
  )

# añade una leyenda al gráfico
legend(
  "topright",                               # posición de la leyenda
  legend = c("Survival", "Cum. Mortality"), # texto de la leyenda  
  lty = c(1, 3),                            # tipos de línea a utilizar en la leyenda
  cex = .85,                                # parámetro que define el tamaño del texto de la leyenda
  bty = "n"                                 # no se dibujará ningún tipo de recuadro para la leyenda
  )

27.4 Comparación de las curvas de supervivencia

Para comparar la supervivencia dentro de los diferentes grupos de nuestros participantes o pacientes observados, es posible que tengamos que observar primero sus respectivas curvas de supervivencia y luego realizar pruebas para evaluar la diferencia entre grupos independientes. Esta comparación puede referirse a grupos basados en el género, la edad, el tratamiento, la comorbilidad…

Test Log rank

El test Log rank (de rango logarítmico) es una prueba popular que compara toda la experiencia de supervivencia entre dos o más grupos independientes y puede considerarse como una prueba de si las curvas de supervivencia son idénticas (se superponen) o no (hipótesis nula de no diferencia de supervivencia entre los grupos). La función survdiff() del paquete survival permite ejecutar el test Log rank cuando especificamos rho = 0 (que es el valor predeterminado). Los resultados de la prueba dan un estadístico chi-cuadrado junto con un valor-p, ya que el estadístico log rank se distribuye aproximadamente como un test estadístico de chi-cuadrado.

En primer lugar, tratamos de comparar las curvas de supervivencia por grupos de género. Para ello, primero intentamos visualizarlo (comprobar si las dos curvas de supervivencia se superponen). Se creará un nuevo objeto survfit con una fórmula ligeramente diferente. Luego se creará el objeto survdiff.

Al suministrar ~ gender como lado derecho de la fórmula, ya no trazamos la supervivencia global sino por género.

# crear el nuevo objeto survfit basado en el género
linelistsurv_fit_sex <-  survfit(Surv(futime, event) ~ gender, data = linelist_surv)

Ahora podemos trazar las curvas de supervivencia por género. Observa el orden de los niveles de los estratos en la columna de género antes de definir los colores y la leyenda.

# establecer colores
col_sex <- c("lightgreen", "darkgreen")

# crear gráfico
plot(
  linelistsurv_fit_sex,
  col = col_sex,
  xlab = "Days of follow-up",
  ylab = "Survival Probability")

# añadir leyenda
legend(
  "topright",
  legend = c("Female","Male"),
  col = col_sex,
  lty = 1,
  cex = .9,
  bty = "n")

Y ahora podemos calcular la prueba de la diferencia entre las curvas de supervivencia utilizando survdiff()

#calcular el test de la diferencia entre las curvas de supervivencia
survival::survdiff(
  Surv(futime, event) ~ gender, 
  data = linelist_surv
  )
## Call:
## survival::survdiff(formula = Surv(futime, event) ~ gender, data = linelist_surv)
## 
## n=4321, 218 observations deleted due to missingness.
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## gender=f 2156      924      909     0.255     0.524
## gender=m 2165      929      944     0.245     0.524
## 
##  Chisq= 0.5  on 1 degrees of freedom, p= 0.5

Vemos que la curva de supervivencia de las mujeres y la de los hombres se superponen y la prueba de rango logarítmico no da pruebas de una diferencia de supervivencia entre mujeres y hombres.

Algunos otros paquetes de R permiten ilustrar curvas de supervivencia para diferentes grupos y probar la diferencia de una sola vez. Utilizando la función ggsurvplot() del paquete survminer, también podemos incluir en nuestra curva las tablas de riesgo impresas para cada grupo, así como el valor p del test log-rank.

PRECAUCIÓN: las funciones de survminer requieren que especifiques el objeto de supervivencia y que vuelvas a especificar los datos utilizados para ajustar el objeto de supervivencia. Recuerda hacer esto para evitar mensajes de error no específicos.

survminer::ggsurvplot(
    linelistsurv_fit_sex, 
    data = linelist_surv,          # especifica de nuevo los datos usados para ajustar linelistsurv_fit_sex 
    conf.int = FALSE,              # no muestra el intervalo de confianza de las estimaciones de KM
    surv.scale = "percent",        # presenta las probabilidades en el eje y en %
    break.time.by = 10,            # presenta el eje temporal con un incremento de 10 días
    xlab = "Follow-up days",
    ylab = "Survival Probability",
    pval = T,                      # imprime el valor p de la prueba Log-rank 
    pval.coord = c(40,.91),        # imprime el valor p en estas coordenadas del gráfico
    risk.table = T,                # imprime la tabla de riesgos en la parte inferior 
    legend.title = "Gender",       # características de la leyenda
    legend.labs = c("Female","Male"),
    font.legend = 10, 
    palette = "Dark2",             # especifica la paleta de colores 
    surv.median.line = "hv",       # dibuja líneas horizontales y verticales a las medianas de supervivencia
    ggtheme = theme_light()        # simplifica el fondo del gráfico
)

Ajustar para que se vea mejor la leyenda

library(grid)  # para unit()

ggsurvplot(
  linelistsurv_fit_sex,
  data            = linelist_surv,
  conf.int        = FALSE,
  surv.scale      = "percent",
  break.time.by   = 10,
  xlab            = "Follow-up days",
  ylab            = "Survival Probability",
  pval            = TRUE,
  pval.coord      = c(40, .91),
  risk.table      = TRUE,
  legend.title    = "Gender",
  legend.labs     = c("Female", "Male"),
  surv.median.line= "hv",
  palette         = "Dark2",
  ggtheme         = theme_light() +
    theme(
      # Texto de leyenda
      legend.title    = element_text(size = 14),
      legend.text     = element_text(size = 12),
      # Tamaño de los recuadros de leyenda (key) para que el número no se corte
      legend.key.width  = unit(1.5, "cm"),
      legend.key.height = unit(1.2, "cm"),
      # Espacio interno entre ítems de leyenda
      legend.spacing.y  = unit(0.5, "lines")
    )
)

También podríamos querer comprobar si hay diferencias en la supervivencia según la fuente de infección (fuente de contaminación).

En este caso, la prueba de rangos logarítmicos proporciona evidencia suficiente de una diferencia en las probabilidades de supervivencia en [número] alpha= 0.005. Las probabilidades de supervivencia de los pacientes infectados en funerales son mayores que las de los pacientes infectados en otros lugares, lo que sugiere una ventaja en la supervivencia.

linelistsurv_fit_source <-  survfit(
  Surv(futime, event) ~ source,
  data = linelist_surv
  )

# gráfico
ggsurvplot( 
  linelistsurv_fit_source,
  data = linelist_surv,
  size = 1, linetype = "strata",   # line types
  conf.int = T,
  surv.scale = "percent",  
  break.time.by = 10, 
  xlab = "Follow-up days",
  ylab= "Survival Probability",
  pval = T,
  pval.coord = c(40,.91),
  risk.table = T,
  legend.title = "Source of \ninfection",
  legend.labs = c("Funeral", "Other"),
  font.legend = 10,
  palette = c("#E7B800","#3E606F"),
  surv.median.line = "hv", 
  ggtheme = theme_light()
)

### Pruebas complementarios al log-rank

Log-RANK

Gehan-Breslow

Tarone-Wane

Peto-Peto

Fleming-Harrington

logrank_test <- survdiff(Surv(futime, event) ~ source, data = linelist_surv, rho = 0)
logrank_test
## Call:
## survdiff(formula = Surv(futime, event) ~ source, data = linelist_surv, 
##     rho = 0)
## 
## n=2906, 1633 observations deleted due to missingness.
## 
##                   N Observed Expected (O-E)^2/E (O-E)^2/V
## source=funeral  438      172      201     4.129      5.17
## source=other   2468     1072     1043     0.795      5.17
## 
##  Chisq= 5.2  on 1 degrees of freedom, p= 0.02
ggsurvplot(linelistsurv_fit_source, data = linelist_surv, pval = TRUE, pval.method = TRUE,break.time.by = 10, legend.title = "Source of \ninfection",palette = c("#E7B800","#3E606F"),risk.table = T,surv.scale = "percent", conf.int = T,
           legend.labs = c("Funeral", "Other"),ggtheme = theme_minimal())

ggsurvplot(linelistsurv_fit_source, data = linelist_surv, pval = TRUE, pval.method = TRUE,  log.rank.weights = "n", pval.method.coord = c(5, 0.1),   pval.method.size = 3,
           break.time.by = 10, legend.title = "Source of \ninfection",palette = c("#E7B800","#3E606F"),risk.table = T,surv.scale = "percent", conf.int = T,
           legend.labs = c("Funeral", "Other"),ggtheme = theme_minimal())

ggsurvplot(linelistsurv_fit_source, data = linelist_surv, pval = TRUE, pval.method = TRUE, log.rank.weights = "sqrtN", pval.method.coord = c(3, 0.1), pval.method.size = 4,
           break.time.by = 10, legend.title = "Source of \ninfection",palette = c("#E7B800","#3E606F"),risk.table = T,surv.scale = "percent", conf.int = T,
           legend.labs = c("Funeral", "Other"),ggtheme = theme_minimal())

ggsurvplot(linelistsurv_fit_source, data = linelist_surv, pval = TRUE, pval.method = TRUE, log.rank.weights = "S1", pval.method.coord = c(5, 0.1), pval.method.size = 3,
           break.time.by = 10, legend.title = "Source of \ninfection",palette = c("#E7B800","#3E606F"),risk.table = T,surv.scale = "percent", conf.int = T,
           legend.labs = c("Funeral", "Other"),ggtheme = theme_minimal())

ggsurvplot(linelistsurv_fit_source, data = linelist_surv, pval = TRUE, pval.method = TRUE, log.rank.weights = "FH_p=1_q=1", pval.method.coord = c(5, 0.1), pval.method.size = 4,
           break.time.by = 10, legend.title = "Source of \ninfection",palette = c("#E7B800","#3E606F"),risk.table = T,surv.scale = "percent", conf.int = T,
           legend.labs = c("Funeral", "Other"),ggtheme = theme_minimal())

# 3. Definir los tiempos que quieres mostrar
times_to_show <- c(10, 20, 30, 40, 50, 60)

# 4. Obtener el summary en esos tiempos
sum_km <- summary(linelistsurv_fit_source, times = times_to_show)

# 5. Convertir a data.frame y limpiar
df_km <- data.frame(
  Source  = sub("source=", "", sum_km$strata),
  time    = sum_km$time,
  n.risk  = sum_km$n.risk,
  n.event = sum_km$n.event,
  survival= sum_km$surv,
  std.err = sum_km$std.err,
  lower   = sum_km$lower,
  upper   = sum_km$upper,
  stringsAsFactors = FALSE
)

# 6. Dar formato numérico consistente
df_km <- df_km %>%
  mutate_at(vars(survival, std.err, lower, upper), ~ round(., 3))

# 7. Crear la flextable
ft <- flextable(df_km) %>%
  # Renombrar cabeceras
  set_header_labels(
    cancer   = "Source",
    time     = "Tiempo",
    n.risk   = "n.riesgo",
    n.event  = "n.eventos",
    survival = "Supervivencia",
    std.err  = "Std.Err",
    lower    = "Lower 95% CI",
    upper    = "Upper 95% CI"
  ) %>%
  # Agrupar filas por tipo de cáncer
  merge_v(j = "cancer") %>%
  # Formato numérico (columnas 4:8)
  colformat_num(j = c("survival","std.err","lower","upper"), digits = 3) %>%
  # Estilo limpio
  theme_booktabs() %>%
  bold(part = "header") %>%
  autofit()
## Warning: Some column(s) can not be found in the original data.frame: "cancer".
# 8. Imprimir la tabla
ft

Source

Tiempo

n.riesgo

n.eventos

Supervivencia

Std.Err

Lower 95% CI

Upper 95% CI

funeral

10

248

132

0.673

0.024

0.629

0.721

funeral

20

89

35

0.549

0.027

0.498

0.605

funeral

30

26

5

0.490

0.035

0.426

0.564

funeral

40

9

0

0.490

0.035

0.426

0.564

funeral

50

3

0

0.490

0.035

0.426

0.564

other

10

1,240

792

0.642

0.010

0.622

0.663

other

20

362

259

0.458

0.013

0.434

0.483

other

30

100

21

0.419

0.014

0.392

0.448

other

40

28

0

0.419

0.014

0.392

0.448

other

50

7

0

0.419

0.014

0.392

0.448

other

60

1

0

0.419

0.014

0.392

0.448

# 1. Instala e importa officer
#install.packages("officer")      # si no lo tienes instalado

# 2. Crea tu flextable (por ejemplo, ft de tu código anterior)
#    Supongamos que ya tienes: ft

# 3. Crea el documento Word y añade la tabla
doc <- read_docx() %>%
  body_add_par("Tabla de supervivencia Kaplan–Meier", style = "heading 1") %>%
  body_add_flextable(ft) %>%
  body_add_par("", style = "Normal")  # salto de línea

# 4. Guarda el documento
print(doc, target = "KM_Survival_Table.docx")

Tasa de riesgo instantánea (hazard)

La tasa de riesgo instantánea, o hazard \(h(t)\), mide la probabilidad condicional de que ocurra el evento de interés (por ejemplo, la muerte) en un instante inmediatamente posterior a \(t\), dado que el sujeto ha sobrevivido hasta ese momento. Formalmente,
\[ h(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}, \]
y se interpreta como el número esperado de eventos por unidad de tiempo entre los sujetos todavía en riesgo. A diferencia de una probabilidad pura, \(h(t)\) es una tasa que puede exceder 1 si se considera un intervalo infinitesimal muy pequeño :contentReferenceoaicite:0.

Función de riesgo acumulado (cumulative hazard)

La función de riesgo acumulado \(H(t)\) cuantifica el “acúmulo” de peligros desde el inicio del estudio hasta el tiempo \(t\). Se define como la integral de la tasa de riesgo:
\[ H(t) = \int_{0}^{t} h(u)\,du. \]
Esta función crece monotonamente y está relacionada directamente con la función de supervivencia \(S(t)\) mediante
\[ S(t) = \exp\bigl(-H(t)\bigr), \]
de modo que a mayor riesgo acumulado, menor probabilidad de que el sujeto siga con vida en \(t\) :contentReferenceoaicite:1.

haz2 <- bshazard(Surv(futime, event) ~ source,  data=linelist_surv)
## NOTE: entry.status has been set to 0 for all.
## Iterations: relative error in phi-hat = 1e-04 
## phi= 0.6932073   sv2= 0.1613212   df= 8.625101   lambda= 4.297062 
## phi= 0.6363475   sv2= 0.2272495   df= 9.986668   lambda= 2.800215 
## phi= 0.6127985   sv2= 0.2749406   df= 10.65197   lambda= 2.22884 
## phi= 0.6017291   sv2= 0.3052854   df= 10.99368   lambda= 1.971038 
## phi= 0.5961656   sv2= 0.3232501   df= 11.17348   lambda= 1.844286 
## phi= 0.5932665   sv2= 0.3334759   df= 11.26941   lambda= 1.779039 
## phi= 0.5917262   sv2= 0.3391741   df= 11.321   lambda= 1.744609 
## phi= 0.5908993   sv2= 0.3423132   df= 11.34888   lambda= 1.726195 
## phi= 0.590453   sv2= 0.3440316   df= 11.36399   lambda= 1.716275 
## phi= 0.5902112   sv2= 0.3449692   df= 11.37218   lambda= 1.71091 
## phi= 0.5900802   sv2= 0.3454799   df= 11.37663   lambda= 1.708002 
## phi= 0.590009   sv2= 0.3457577   df= 11.37904   lambda= 1.706423 
## phi= 0.5899703   sv2= 0.3459087   df= 11.38036   lambda= 1.705567 
## phi= 0.5899494   sv2= 0.3459908   df= 11.38107   lambda= 1.705101 
## phi= 0.5899379   sv2= 0.3460355   df= 11.38146   lambda= 1.704848 
## phi= 0.5899317   sv2= 0.3460597   df= 11.38167   lambda= 1.704711
plot(haz2, conf.int=F,xlab = "Tiempo", ylab = "Hazard instantaneo", main = "Estimacion del Hazard Instantaneo (bshazard)")

## A medida que pasa el tiempo incrementa el peligro de fallecer

ggsurvplot(linelistsurv_fit_source, fun="cumhaz", xlab="Tiempo", ylab="Hazard", title="Grafico de Hazard Acumulado", conf.int = T,
            break.time.by = 10,palette = c("#E7B800","#3E606F"),
           risk.table = T)

27.5 Análisis de regresión de Cox

La regresión de riesgos proporcionales de Cox es una de las técnicas de regresión más populares para el análisis de supervivencia. También se pueden utilizar otros modelos, ya que el modelo de Cox requiere supuestos importantes que deben verificarse para un uso adecuado, como el supuesto de riesgos proporcionales: véanse las referencias.

En un modelo de regresión de riesgos proporcionales de Cox, la medida del efecto es la tasa de riesgo (HR), que es el riesgo de fracaso (o el riesgo de muerte en nuestro ejemplo), dado que el participante ha sobrevivido hasta un momento específico. Normalmente, nos interesa comparar grupos independientes con respecto a sus riesgos, y utilizamos una razón de riesgo, que es análoga a una razón de probabilidades en el entorno del análisis de regresión logística múltiple. La función cox.ph() del paquete de supervivencia se utiliza para ajustar el modelo. La función cox.zph() del paquete survival puede utilizarse para probar la suposición de riesgos proporcionales para un ajuste del modelo de regresión de Cox.

NOTA: Una probabilidad debe estar en el rango de 0 a 1. Sin embargo, el peligro representa el número esperado de eventos por una unidad de tiempo.

\[ h(t \mid \mathbf{X}) \;=\; h_0(t)\,\exp\bigl(\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p\bigr) \]
o, equivalentemente, en escala logarítmica:

\[ \log h(t \mid \mathbf{X}) \;=\; \log h_0(t)\;+\;\sum_{i=1}^{p}\beta_i\,X_i \]

Primero podemos ajustar un modelo para evaluar el efecto de la edad y el sexo en la supervivencia. Con sólo imprimir el modelo, tenemos la información sobre:

La función summary() aplicada al objeto del modelo de Cox ofrece más información, como el intervalo de confianza de la RR estimada y las diferentes puntuaciones de la prueba.

El efecto de la primera covariable, gender, se presenta en la primera fila. Se imprime genderm (masculino), lo que implica que el primer nivel de estrato (“f”), es decir, el grupo femenino, es el grupo de referencia para el género. Por lo tanto, la interpretación del parámetro de la prueba es la de los hombres en comparación con las mujeres. El valor p indica que no hay pruebas suficientes de un efecto del género sobre el peligro esperado o de una asociación entre el género y la mortalidad por todas las causas.

La misma falta de pruebas se observa en relación con el grupo de edad.

#ajustar el modelo cox
linelistsurv_cox_sexage <-  survival::coxph(
              Surv(futime, event) ~ gender + age_cat_small, 
              data = linelist_surv
              )


#impresión del modelo ajustado
linelistsurv_cox_sexage
## Call:
## survival::coxph(formula = Surv(futime, event) ~ gender + age_cat_small, 
##     data = linelist_surv)
## 
##                       coef exp(coef) se(coef)      z     p
## genderm           -0.03149   0.96900  0.04767 -0.661 0.509
## age_cat_small5-19  0.09400   1.09856  0.06454  1.456 0.145
## age_cat_small20+   0.05032   1.05161  0.06953  0.724 0.469
## 
## Likelihood ratio test=2.8  on 3 df, p=0.4243
## n= 4321, number of events= 1853 
##    (218 observations deleted due to missingness)
#resumen del modelo
summary(linelistsurv_cox_sexage)
## Call:
## survival::coxph(formula = Surv(futime, event) ~ gender + age_cat_small, 
##     data = linelist_surv)
## 
##   n= 4321, number of events= 1853 
##    (218 observations deleted due to missingness)
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)
## genderm           -0.03149   0.96900  0.04767 -0.661    0.509
## age_cat_small5-19  0.09400   1.09856  0.06454  1.456    0.145
## age_cat_small20+   0.05032   1.05161  0.06953  0.724    0.469
## 
##                   exp(coef) exp(-coef) lower .95 upper .95
## genderm               0.969     1.0320    0.8826     1.064
## age_cat_small5-19     1.099     0.9103    0.9680     1.247
## age_cat_small20+      1.052     0.9509    0.9176     1.205
## 
## Concordance= 0.514  (se = 0.007 )
## Likelihood ratio test= 2.8  on 3 df,   p=0.4
## Wald test            = 2.78  on 3 df,   p=0.4
## Score (logrank) test = 2.78  on 3 df,   p=0.4

Fue interesante ejecutar el modelo y observar los resultados, pero un primer vistazo para verificar si se respetan los supuestos de riesgos proporcionales podría ayudar a ahorrar tiempo.

test_ph_sexage <- survival::cox.zph(linelistsurv_cox_sexage)
test_ph_sexage
##               chisq df    p
## gender        0.454  1 0.50
## age_cat_small 0.838  2 0.66
## GLOBAL        1.399  3 0.71

NOTA: Se puede especificar un segundo argumento llamado método cuando se calcula el modelo de Cox, que determina cómo se manejan los empates. El valor por defecto es “efron”, y las otras opciones son “breslow” y “exact”.

En otro modelo añadimos más factores de riesgo, como el origen de la infección y el número de días entre la fecha de inicio y el ingreso. Esta vez, primero verificamos la hipótesis de riesgos proporcionales antes de seguir adelante.

En este modelo, hemos incluido un predictor continuo (days_onset_hosp). En este caso, interpretamos las estimaciones de los parámetros como el aumento del logaritmo esperado del riesgo relativo por cada aumento de una unidad en el predictor, manteniendo los demás predictores constantes. Primero verificamos el supuesto de riesgos proporcionales.

#ajustar el modelo
linelistsurv_cox <-  coxph(
                        Surv(futime, event) ~ gender + age_years+ source + days_onset_hosp,
                        data = linelist_surv
                        )


#comprobar el modelo de riesgo proporcional
linelistsurv_ph_test <- cox.zph(linelistsurv_cox)
linelistsurv_ph_test
##                    chisq df       p
## gender           0.45062  1    0.50
## age_years        0.00199  1    0.96
## source           1.79622  1    0.18
## days_onset_hosp 31.66167  1 1.8e-08
## GLOBAL          34.08502  4 7.2e-07

La verificación gráfica de esta suposición puede realizarse con la función ggcoxzph() del paquete survminer mediante los residuos de Schoenfield.

survminer::ggcoxzph(linelistsurv_ph_test)

#Prueba de bondad de Ajuste-proporcionalidad
linelistsurv_ph_test <- cox.zph(linelistsurv_cox)
linelistsurv_ph_test
##                    chisq df       p
## gender           0.45062  1    0.50
## age_years        0.00199  1    0.96
## source           1.79622  1    0.18
## days_onset_hosp 31.66167  1 1.8e-08
## GLOBAL          34.08502  4 7.2e-07
# grafico pruebas
plot(linelistsurv_ph_test)

Los resultados del modelo indican que existe una asociación negativa entre la duración del inicio del ingreso y la mortalidad por todas las causas. El riesgo esperado es 0,9 veces menor en una persona que ingresa un día más tarde que otra, manteniendo el género constante. O, en una explicación más directa, un aumento de una unidad en la duración del inicio al ingreso se asocia con una disminución del 10,7% (coef *100) en el riesgo de muerte.

Los resultados muestran también una asociación positiva entre la fuente de infección y la mortalidad por todas las causas. Es decir, hay un mayor riesgo de muerte (1,21 veces) para los pacientes que tuvieron una fuente de infección distinta de los funerales.

#imprimir el resumen del modelo
summary(linelistsurv_cox)
## Call:
## coxph(formula = Surv(futime, event) ~ gender + age_years + source + 
##     days_onset_hosp, data = linelist_surv)
## 
##   n= 2772, number of events= 1180 
##    (1767 observations deleted due to missingness)
## 
##                      coef exp(coef)  se(coef)      z Pr(>|z|)    
## genderm          0.004710  1.004721  0.060827  0.077   0.9383    
## age_years       -0.002249  0.997753  0.002421 -0.929   0.3528    
## sourceother      0.178393  1.195295  0.084291  2.116   0.0343 *  
## days_onset_hosp -0.104063  0.901169  0.014245 -7.305 2.77e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                 exp(coef) exp(-coef) lower .95 upper .95
## genderm            1.0047     0.9953    0.8918    1.1319
## age_years          0.9978     1.0023    0.9930    1.0025
## sourceother        1.1953     0.8366    1.0133    1.4100
## days_onset_hosp    0.9012     1.1097    0.8764    0.9267
## 
## Concordance= 0.566  (se = 0.009 )
## Likelihood ratio test= 71.31  on 4 df,   p=1e-14
## Wald test            = 59.22  on 4 df,   p=4e-12
## Score (logrank) test = 59.54  on 4 df,   p=4e-12

Podemos comprobar esta relación con una tabla:

linelist_case_data %>% 
  tabyl(days_onset_hosp, outcome) %>% 
  adorn_percentages() %>%  
  adorn_pct_formatting()
##  days_onset_hosp Death Recover   NA_
##                0 44.3%   31.4% 24.3%
##                1 46.6%   32.2% 21.2%
##                2 43.0%   32.8% 24.2%
##                3 45.0%   32.3% 22.7%
##                4 41.5%   38.3% 20.2%
##                5 40.0%   36.2% 23.8%
##                6 32.2%   48.7% 19.1%
##                7 31.8%   38.6% 29.5%
##                8 29.8%   38.6% 31.6%
##                9 30.3%   51.5% 18.2%
##               10 16.7%   58.3% 25.0%
##               11 36.4%   45.5% 18.2%
##               12 18.8%   62.5% 18.8%
##               13 10.0%   60.0% 30.0%
##               14 10.0%   50.0% 40.0%
##               15 28.6%   42.9% 28.6%
##               16 20.0%   80.0%  0.0%
##               17  0.0%  100.0%  0.0%
##               18  0.0%  100.0%  0.0%
##               22  0.0%  100.0%  0.0%
##               NA 52.7%   31.2% 16.0%

Habría que considerar e investigar por qué existe esta asociación en los datos. Una posible explicación podría ser que los pacientes que viven lo suficiente como para ser ingresados más tarde tenían una enfermedad menos grave para empezar. Otra explicación, quizá más probable, es que, dado que utilizamos unos datos falsos simulados, este patrón no refleja la realidad.

Forest plots

A continuación, podemos visualizar los resultados del modelo cox utilizando los prácticos gráficos de bosque con la función ggforest() del paquete survminer.

ggforest(linelistsurv_cox, data = linelist_surv)

# Models Summary

Con modelsummary se puede organizar mejor

https://modelsummary.com/vignettes/modelsummary.html

# Ajustar los tres modelos
cox_model1 <- coxph(Surv(futime, event) ~ gender + age_years + source + days_onset_hosp, data = linelist_surv)

cox_model2 <- coxph(Surv(futime, event) ~ ct_blood + bmi + temp + source + days_onset_hosp, data = linelist_surv)

cox_model3 <- coxph(Surv(futime, event) ~ fever + chills + cough + aches + vomit + source + days_onset_hosp, data = linelist_surv)

cox_model4 <- coxph(Surv(futime, event) ~ bmi + age_years + cough + aches + temp + source + days_onset_hosp, data = linelist_surv)

cox_model5 <- coxph(Surv(futime, event) ~ gender + chills + cough + aches + age_years + source + days_onset_hosp, data = linelist_surv)

cox_model6 <- coxph(Surv(futime, event) ~ ct_blood + bmi + temp + age_years + gender + source + days_onset_hosp, data = linelist_surv)
# Crear la tabla comparativa
modelsummary(
  list(
    "Model 1" = cox_model1,
    "Model 2" = cox_model2,
    "Model 3" = cox_model3,
    "Model 4" = cox_model4,
    "Model 5" = cox_model5,
    "Model 6" = cox_model6
    
  ),
  exponentiate = TRUE, # Para mostrar HR en lugar de coeficientes
  statistic = "conf.int", # Puedes usar también "{std.error} ({p.value})"
  stars = TRUE, # Para mostrar significancia
  output = "html" # Puedes cambiar a "latex", "html", etc. según lo necesites
)
Model 1 Model 2 Model 3 Model 4 Model 5 Model 6
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
genderm 1.005 1.001 1.013
[0.892, 1.132] [0.886, 1.130] [0.898, 1.143]
age_years 0.998 0.997 0.998 0.996
[0.993, 1.002] [0.992, 1.002] [0.993, 1.003] [0.991, 1.001]
sourceother 1.195* 1.189* 1.164+ 1.162+ 1.176+ 1.206*
[1.013, 1.410] [1.010, 1.400] [0.988, 1.372] [0.983, 1.373] [0.994, 1.392] [1.020, 1.426]
days_onset_hosp 0.901*** 0.889*** 0.900*** 0.900*** 0.902*** 0.888***
[0.876, 0.927] [0.858, 0.920] [0.876, 0.926] [0.875, 0.926] [0.877, 0.927] [0.857, 0.920]
ct_blood 0.974 0.969
[0.934, 1.017] [0.928, 1.013]
bmi 0.999+ 0.999+ 0.999*
[0.998, 1.000] [0.998, 1.000] [0.998, 1.000]
temp 0.992 0.983 1.010
[0.936, 1.050] [0.924, 1.046] [0.952, 1.072]
feveryes 0.927
[0.804, 1.068]
chillsyes 1.022 1.012
[0.886, 1.179] [0.874, 1.172]
coughyes 1.007 0.988 1.027
[0.853, 1.189] [0.836, 1.169] [0.866, 1.219]
achesyes 0.973 1.002 0.962
[0.799, 1.185] [0.820, 1.223] [0.784, 1.179]
vomityes 1.113+
[0.993, 1.247]
Num.Obs. 2772 2832 2791 2672 2663 2699
AIC 17409.5 17990.3 17640.8 16856.5 16676.2 16930.8
BIC 17433.3 18020.0 17682.3 16897.7 16717.4 16972.1
RMSE 0.66 0.66 0.66 0.66 0.66 0.66 
# Versión para guardar en word
modelsummary(
  list(
    "Model 1" = cox_model1,
    "Model 2" = cox_model2,
    "Model 3" = cox_model3
  ),
  exponentiate = TRUE, # Para mostrar HR en lugar de coeficientes
  statistic = "conf.int", # Puedes usar también "{std.error} ({p.value})"
  stars = TRUE, # Para mostrar significancia
  output = "Modelos.docx" # Puedes cambiar a "latex", "html", etc. según lo necesites
)

Models Summary saca gráficas para comparar modelos de regresión

# Graficar comparacion de modelos
modelplot(
  list(
    "Modelo 1" = cox_model1,
    "Modelo 2" = cox_model2,
    "Model 3" = cox_model3,
    "Model 4" = cox_model4,
    "Model 5" = cox_model5,
    "Model 6" = cox_model6
  ),
  exponentiate = TRUE,  # Mostrar HR en lugar de coeficientes
  coef_omit = "(Intercept)",  # Omitir intercepto si aparece
  draw = TRUE
) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray") +
  labs(
    title = "Comparacion de Modelos de Cox",
    x = "Hazard Ratio (HR)",
    y = "Covariable",
    caption = "Batra, Neale, et al. Manual de R para Epidemiologia. 2021" ) +
  theme_minimal(base_size = 12)

cox_model1 <- coxph(Surv(futime, event) ~ gender + age_years + source + days_onset_hosp, data = linelist_surv)
cox_model4 <- coxph(Surv(futime, event) ~ gender + age_years + source , data = linelist_surv)
# Comparación de Modelos
anova(cox_model1,cox_model4)
## Analysis of Deviance Table
##  Cox model: response is  Surv(futime, event)
##  Model 1: ~ gender + age_years + source + days_onset_hosp
##  Model 2: ~ gender + age_years + source
##    loglik  Chisq Df Pr(>|Chi|)    
## 1 -8700.8                         
## 2 -8733.2 64.863  1  8.029e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

27.6 Covariables tiempo-dependientes en modelos de supervivencia

Algunas de las siguientes secciones han sido adaptadas con permiso de la excelente introducción al análisis de supervivencia en R por la Dra. Emily Zabor

https://www.emilyzabor.com/survival-analysis-in-r.html

En la última sección hemos tratado el uso de la regresión de Cox para examinar las asociaciones entre las covariables de interés y los resultados de supervivencia, pero estos análisis dependen de que la covariable se mida en la línea de base, es decir, antes de que comience el tiempo de seguimiento del evento.

¿Qué ocurre si tienes interés en una covariable que se mide después del tiempo de seguimiento? O, ¿qué pasa si tienes una covariable que puede cambiar con el tiempo?

Por ejemplo, tal vez estés trabajando con datos clínicos en los que se repiten medidas de valores de laboratorio del hospital que pueden cambiar con el tiempo. Este es un ejemplo de una covariable dependiente del tiempo. Para abordar esto se necesita una configuración especial, pero afortunadamente el modelo de Cox es muy flexible y este tipo de datos también puede ser modelado con herramientas del paquete survival.

Configuración de covariables dependientes del tiempo El análisis de covariables dependientes del tiempo en R requiere la configuración de unos datos especial. Si tienes interés, mira el documento del autor del paquete survival Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model.

Para ello, utilizaremos un nuevo conjunto de datos del paquete SemiCompRisks denominado BMT, que incluye datos de 137 pacientes de trasplante de médula ósea. Las variables en las que nos centraremos son:

  • T1 - tiempo (en días) hasta la muerte o el último seguimiento
  • delta1 - indicador de muerte; 1-Muerto, 0-Vivo
  • TA - tiempo (en días) hasta la enfermedad aguda de injerto contra huésped
  • deltaA - indicador de la enfermedad aguda de injerto contra huésped;
    • 1 - Desarrolló la enfermedad aguda de injerto contra huésped
    • 0 - Nunca desarrolló la enfermedad aguda de injerto contra huésped

Cargaremos este conjunto de datos del paquete survival utilizando el comando DE R base data(), que puede utilizarse para cargar datos que ya están incluidos en un paquete de R que se ha cargado. El dataframe BMT aparecerá en tu entorno de R.

data(BMT, package = "SemiCompRisks")

Añadir un identificador único de paciente

No hay una columna de identificación única en los datos de BMT, que es necesaria para crear el tipo de conjunto de datos que queremos. Así que utilizamos la función rowid_to_column() del paquete tibble de tidyverse para crear una nueva columna de identificación llamada my_id (añade una columna al principio del dataframe con identificadores de fila secuenciales, empezando por el 1). Llamamos al dataframe bmt.

bmt <- rowid_to_column(BMT, "my_id")
names(bmt)
##  [1] "my_id"  "g"      "T1"     "T2"     "delta1" "delta2" "delta3" "TA"    
##  [9] "deltaA" "TC"     "deltaC" "TP"     "deltaP" "Z1"     "Z2"     "Z3"    
## [17] "Z4"     "Z5"     "Z6"     "Z7"     "Z8"     "Z9"     "Z10"

Consulta de Tabla de Base de Datos

bmt <- rowid_to_column(BMT, "my_id")
datatable(bmt, options = list(pageLength = 5))  # tabla navegable

Ampliar las filas de pacientes

A continuación, utilizaremos la función tmerge() con las funciones de ayuda event() y tdc() para crear el conjunto de datos reestructurado. Nuestro objetivo es reestructurar el conjunto de datos para crear una fila separada para cada paciente por cada intervalo de tiempo en el que tengan un valor diferente de deltaA. En este caso, cada paciente puede tener como máximo dos filas dependiendo de si desarrollaron la enfermedad aguda de injerto contra huésped durante el periodo de recogida de datos. Llamaremos a nuestro nuevo indicador para el desarrollo de la enfermedad aguda de injerto contra huésped agvhd.

  • tmerge() crea unos datos largos con múltiples intervalos de tiempo para los diferentes valores de las covariables de cada paciente
  • event() crea el nuevo indicador de eventos para que vaya con los intervalos de tiempo recién creados
  • tdc() crea la columna de covarianza dependiente del tiempo, agvhd, para que vaya con los intervalos de tiempo recién creados
td_dat <- 
  tmerge(
    data1 = bmt %>% select(my_id, T1, delta1), 
    data2 = bmt %>% select(my_id, T1, delta1, TA, deltaA), 
    id = my_id, 
    death = event(T1, delta1),
    agvhd = tdc(TA)
    )

Para ver qué hace esto, veamos los datos de los 5 primeros pacientes individuales.

Las variables de interés en los datos originales tenían este aspecto:

bmt %>% 
  select(my_id, T1, delta1, TA, deltaA) %>% 
  filter(my_id %in% seq(1, 5))
##   my_id   T1 delta1   TA deltaA
## 1     1 2081      0   67      1
## 2     2 1602      0 1602      0
## 3     3 1496      0 1496      0
## 4     4 1462      0   70      1
## 5     5 1433      0 1433      0

Ahora algunos de nuestros pacientes tienen dos filas en el conjunto de datos correspondientes a intervalos en los que tienen un valor diferente de nuestra nueva variable, agvhd. Por ejemplo, el paciente 1 tiene ahora dos filas con un valor de agvhd de cero desde el tiempo 0 hasta el tiempo 67, y un valor de 1 desde el tiempo 67 hasta el tiempo 2081.

Regresión de Cox con covariables dependientes del tiempo

Ahora que hemos remodelado nuestros datos y añadido la nueva variable aghvd dependiente del tiempo, vamos a ajustar un modelo de regresión cox simple de una sola variable. Podemos utilizar la misma función coxph() que antes, sólo tenemos que cambiar nuestra función Surv() para especificar tanto el tiempo de inicio como el de finalización de cada intervalo utilizando los argumentos time1 = y time2 =.

bmt_td_model = coxph(
  Surv(time = tstart, time2 = tstop, event = death) ~ agvhd, 
  data = td_dat
  )

summary(bmt_td_model)
## Call:
## coxph(formula = Surv(time = tstart, time2 = tstop, event = death) ~ 
##     agvhd, data = td_dat)
## 
##   n= 163, number of events= 80 
## 
##         coef exp(coef) se(coef)    z Pr(>|z|)
## agvhd 0.3351    1.3980   0.2815 1.19    0.234
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## agvhd     1.398     0.7153    0.8052     2.427
## 
## Concordance= 0.535  (se = 0.024 )
## Likelihood ratio test= 1.33  on 1 df,   p=0.2
## Wald test            = 1.42  on 1 df,   p=0.2
## Score (logrank) test = 1.43  on 1 df,   p=0.2

De nuevo, visualizaremos los resultados de nuestro modelo de Cox utilizando la función ggforest() del paquete urvminer.:

ggforest(bmt_td_model, data = td_dat)

Como se puede ver en el gráfico de forest, el intervalo de confianza y el valor-p, no parece haber una fuerte asociación entre la muerte y la enfermedad aguda de injerto contra huésped en el contexto de nuestro modelo simple.

27.7 Recursos