Capítulo 1 ¿Qué es el análisis de sobrevivencia?

El análisis de sobrevivencia es un conjunto de técnicas estadísticas en las que la variable respuesta es el tiempo que transcurre entre el comienzo de seguimiento del individuo en el estudio y la aparición del evento de interés (tolerancia a la leche, muerte, etc.).

En el análisis de sobrevivencia, la variable respuesta, es generalmente, el tiempo hasta la ocurrencia de un evento de interés.

1.1 Tiempos de fallas.

Tiempo hasta la muerte de un paciente.
Tiempo hasta la recaída de una enfermedad.
Tiempo hasta la cura del paciente.

Ejemplo: En estudios de cáncer es usual el registro de fechas correspondientes al diágnostico de la enfermedad, a la remisión (después del tratamiento, el paciente queda libre de los síntomas de la enfermedad), a la recurrencia de la enfermedad (recaída) y hasta la muerte del paciente.

1.1.1 Aplicaciones.

Asumamos que estamos interesados en el tiempo hasta que ocurra un cierto evento de interés.

Medicina: Tiempo hasta que un tratamiento surga efecto.
Ingeniería: Tiempo hasta que un componente electrónico se deteriore.
Educación: Tiempo hasta que un alumno adquiera cierto conocimiento.
Ecología: Tiempo hasta que un cierto animal alcance un determinado tamaño, etc.

1.1.2 Conceptos.

Es necesario definir los siguientes elementos: El tiempo de inicio del estudio debe ser definido de forma precisa, de forma que las observaciones puedan ser comparables en el origen del estudio.

La escala de medida: Usualmente el tiempo (cronológico).
El evento de interés: Nos enfocaremos en el caso en que sólo puede ocurrir debido a una única causa.

1.1.3 Datos completos.

Indican que todos los individuos en el estudio fallaron o el tiempo de falla de cada individuo es conocido. Por lo tanto asumimos informaciones completas.

1.2 Censura.

Los estudios que involucran una respuesta temporal usualmente son de una “larga” duraci ́on. Sin embargo, a pesar de ser largos en duración, muchos estudios clínicos terminan antes de que todos los individuos dentro del estudio “fallen”.

De esta forma, tenemos observaciones “incompletas” dentro de nuestro estudio. Estas observaciones llamadas censuras, pueden ocurrir por una variedad de razones.

Siendo algunas de ellas:

Muertes por una causa diferente a la relacionada con el estudio.
Término del estudio.
Pérdida de contacto con el paciente.
Rechazar continuar participando del estudio y consecuentemente hace abandono por parte del paciente.
Abandono del estudio debido a efectos colaterales del tratamiento.

Nota: Los tiempos de censura, deben ser incorporados en el análisis porque traen información sobre un individuo o componente.

1.3 Tipos de censura.

Censura a la derecha: El tiempo de ocurrencia del evento de interés está a la derecha del tiempo registrado.
Censura a la izquierda: Cuando el tiempo resgistrado es mayor al tiempo de falla. Es decir, el evento ya aconteció cuando el individuo fue observado.
Censura a la derecha e izquierda: Datos doblemente censurados.
Censura Intervalar: El evento ocurrió en un cierto intervalo de tiempo.

Para efecto del desarrollo de este documento se trabajará sólo con censura a la derecha.

1.3.1 Censura a la derecha.

Existen tres mecanismos de censura a la derecha.

Censura Tipo I: El estudio será determinado después de un período pre-establecido de tiempo. Las observaciones cuyo evento de interés no fueron observadas hasta este tiempo son dichas censuradas. Es más común en estudios médicos.
Censura Tipo II: El término del estudios será determinado una vez ocurrido el evento de interés en un número pre-establecido de individuos, es decir, en este esquema de censura, se fija en número r de fallas y el experimento es conducido con n unidades en el estudio, n > r, hasta que r de esas unidades hayan fallado.
Censura Aleátoria: Ocurre si una observación fue retirada del estudio sin haber ocurrido el evento de interés, por una razón diferente a la estudiada.

1.4 Truncación.

La truncación ocurre cuando parte de la población no puede ser observada.

Es caracterizado por una condición que excluye ciertos individuos del estudio. En estos estudios, los pacientes no son acompañados a partir del tiempo inicial, pero solo después de experimentar un cierto evento.

1.5 Tiempo de Sobrevivencia.

La v.a. no negativa T, usualmente continua, que representa el tiempo de falla, es generalmente representada por su función de sobrevivencia o por la función de tasa de falla (o de riesgo).

1.5.1 Función de Sobrevivencia.

La función de sobrevivencia es definada como la probabilidad de una observación no fallar hasta un cierto tiempo t, es decir, la probabilidad de una observación sobrevivir hasta un cierto tiempo t. En terminos probabilísticos, esto es descrito como:

\[\begin{equation} S(t) = P(T > t) = 1 − P(T < t) \end{equation}\]

Propiedades

$S(t) = 0$, cuando $t \to \infty$.
$S(t) = 1$, cuando $t \to 0$.

1.5.2 Función de riesgo.

Función de riesgo o tasa de falla: Es la tasa de falla instantánea en el tiempo t condicional a que la observación sobrevivió hasta el tiempo t.

\[\begin{equation} h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t} \end{equation}\]

Propiedades

$h(t) \to 0$, cuando $t \to \infty$.
$h(t)$ es una tasa y no una probabilidad.

1.5.3 Función de riesgo acumulada.

Función de riesgo acumulada:

\[\begin{equation} H(t) = \int_{0}^{t} h(u) \, du \end{equation}\]

1.5.4 Tiempo medio y Vida media residual.

El área bajo la curva de sobrevivencia es definida como el tiempo medio, es decir,

\[\begin{equation} tm = \int_{0}^{\infty} S(t) \, dt \end{equation}\]

donde $S(t)$ es la función de sobrevivencia. Además, se define el $vmr(t)$ como:

\[\begin{equation} vmr(t) = \int_{t}^{\infty} (u - t) f(u) \, du \end{equation}\]

y finalmente, otra expresión para $S(t)$ en términos de $S(u)$:

\[\begin{equation} S(t) = \int_{t}^{\infty} S(u) \, du \end{equation}\]

1.5.5 Algunas Relaciones Importantes.

Tabla 1

\[ \begin{array}{c|c} \hline \text{Función requerida} & \text{Función dada} \\ \hline f(t) & S(t) \\ - \int_{t}^{\infty} f(u) \, du & - \frac{dS(t)}{dt} \\ h(t) \exp \left( - \int_{0}^{t} h(u) \, du \right) & H(t) \\ \frac{dH(t)}{dt} & \exp \{ -H(t) \} \\ \hline \end{array} \]

Capítulo 2 Instalación de paquetes.

Este manual fue realizado utilizando la versión 4.2.1 (2022-06-23 ucrt) de R.

Las librerías se instalan por una única vez y se cargan antes de iniciar cualquier análisis. La instalación y las librerías utilizadas se realizan en el siguiente script:

#install.packages('tinytex')
#install.packages("survival")
#install.packages("KMsurv")
#install.packages("survMisc")
#install.packages("survminer")
#install.packages("flexsurv")
#install.packages("actuar")
#install.packages("dplyr")

La carga de las librerías se realiza con la función library():

suppressPackageStartupMessages({
  library(survival)
  library(KMsurv)
  library(survMisc)
  library(survminer)
  library(flexsurv)
  library(actuar)
  library(dplyr)
  library(tinytex)
})

2.1 Bases de datos

Al menos que se indique lo contrarío las bases de datos utilizadas son extraídas de Klein and Moeschberger (1997), Survival Analysis, Techniques for Censored and Truncated Data, Springer obtenidas a través de la librería KMSurv.

Las bases de da tos de Klein and Moeschberger se pueden accesar utilizando la función data(name), en donde name toma el nombre de la base de datos que se desea extraer, esta función crea una variable del mismo nombre que contiene la base de datos.

Ejemplo de la base de datos.

library(survival)
data("ovarian")
head(ovarian)

##   futime fustat     age resid.ds rx ecog.ps
## 1     59      1 72.3315        2  1       1
## 2    115      1 74.4932        2  1       1
## 3    156      1 66.4658        2  1       2
## 4    421      0 53.3644        2  2       1
## 5    431      1 50.3397        2  1       1
## 6    448      0 56.4301        1  1       2

2.2 Preparación de los datos.

R utiliza el directorio de trabajo para leer y escribir archivos. Para saber cuál es este directorio puede utilizar el comando getwd(). Se recomienda cambiar el directorio de trabajo utilizando la función setwd() por ejemplo, setwd(“C:/data”) o setwd(“/home/R”).

Es necesario proporcionar la dirección completa del archivo si este no se encuentra en el directorio de trabajo, para efectos del desarrollo del presente documento se trabajaran sólo con bases de datos que están presentes en el directorio.

Antes de iniciar con los análisis se requiere un pre-procesamiento de los datos, estos deberán tener el siguiente formato:

La primera fila, conocida como encabezado, corresponde al nombre de cada variable.
Las filas posteriores al encabezado representan las observaciones o individuos.
Los valores faltantes se deberán capturar como espacios vacíos.
Las variables categóricas se deberán codificar como números enteros.
El archivo debe estar en formato de “valores separados por comas” (.csv).
El archivo debe estar guardado en el directorio de trabajo.

Este proceso se puede realizar en un software de hojas de cálculo y posteriormente guardar el archivo en formato csv.

2.3 El Objeto Surv.

La mayoría de las funciones utilizan un objeto tipo Surv como argumento, un objeto Surv no es más que la combinación de información entre los tiempos y su censura. Existen varios tipos de censura, sin embargo, en esta oportunidad solo se ejemplificarán con datos con censura por la derecha.

La creación de un objetivo Surv se realiza con la función Surv(Time,Event), que toma como argumentos: Time como el tiempo hasta la observación y Event un indicador en donde es 1 si se observó el evento y 0 si es censura.

Por ejemplo:

Surv(5, 1) #Observación del evento

## [1] 5

Surv(5, 0) #Observación Censurada

## [1] 5+

Surv(c(1, 2, 3, 4), c(1, 0, 1, 0))

## [1] 1  2+ 3  4+

Los tiempos censurados se caracterizan por tener el símbolo + acompañado a su derecha.

Como se verá posteriormente, Surv(Time,Event) en combinación de otras funciones tomará como argumento el Nombre de la columna que hace referencia tanto a los tiempos de observación, así como su censura.

Capítulo 3 Estimación No - Paramétrica de la Función de Sobrevivencia.

3.1 Estimador de Kaplan-Meier.

Los estimadores de Kaplan-Meier para la función de sobrevivencia es obtenido a través paquete estadístico survival mediante la función survfit(). Esta función en su forma más sencilla, solo requiere un objeto de sobrevivencia creado por la función Surv(). Los argumentos de la función survfit() son los siguientes:

formula: Un objeto fórmula y ~ x, que debe tener un objeto Surv como variable respuesta a la izquierda del ~ y, si se desea, el nombre de las covariables por la derecha. Uno de los términos puede ser un objeto strata. Para una sola curva de sobrevivencia del lado derecho se coloca ~ 1.
data: Objeto data.frame donde están los datos.
type: Tipo de estimador: "kaplan-meier".

Ejemplo:

La base de datos tongue con tiene un total de 80 observaciones de 3 variables:

time. Tiempo hasta la muerte o en estudio en semanas.
delta. Indicador de muerte(0=Vivo,1=Muerte).
type. Perfil de ADN tumoral.(1=Tumor Aneuploide, 2=Tumor Diploide).

library(KMsurv)
data("tongue")
head(tongue)

##   type time delta
## 1    1    1     1
## 2    1    3     1
## 3    1    3     1
## 4    1    4     1
## 5    1   10     1
## 6    1   13     1

#Guardando el Objeto Surv
tongue.surv <- Surv(tongue$time, tongue$delta) #Objeto tipo Surv
tongue.km <- survfit(tongue.surv ~ 1, data = tongue,
                     type = "kaplan-meier") #Estimación KM

#Sin Guardar el Objeto Surv
tongue.km <- survfit(Surv(time, delta) ~ 1, data = tongue, type = "kaplan-meier")

La estimación de la función de sobrevivencia se lleva a cabo con la función summary().

summary(tongue.km)

La estimación devuelve los siguientes valores:

time : Tiempo de la observación
n.risk : Número de sujetos en riesgo de experimentar el evento en cada tiempo.
n.evento : Número de eventos (como muertes) que ocurrieron en cada punto de tiempo.
survival : Estimación de la probabilidad de sobrevivencia hasta cada punto de tiempo.
std.err : Error estándar de la estimación de la sobrevivencia.
lower y upper CI : Los intervalos de confianza para la estimación.

La función survfit() devuelve un resumen de la estimación, la información se puede acceder agregando el símbolo “$” seguido del nombre del elemento de la lista.

Una mejor manera de extraer la información es utilizando la función fortify() sobre el objeto survfit, esta función devuelve un data.frame con la información. Al tener presencia de covariables se anexa la columna llamada strata.

Ejemplo, utilizando la variable delta como covariable.

tongue.km <- survfit(Surv(time, delta) ~ type, data = tongue, type = "kaplan-meier") #Estimación KM por grupos
summary(tongue.km) #Resumen estadístico

tongue.km[2]$surv #Elementos columna surv delta=1

library(broom)
tongue_ovarian <- tidy(tongue.km)
#fortify(tongue.km)

3.2 Desviación estándar e Intervalos de confianza de la estimación de la función de sobrevivencia.

Tanto como la desviación estándar y los Intervalos de confianza de la curva de sobrevivencia es estimada mediante la función survfit() con los siguientes argumentos:

error: Tipo de estimación para las desviaciones, los posibles valores son “greenwood”(defecto) para la fórmula de Greenwodd o “tsiatis” para la fórmula de Tsiatis/Aalen.
conf.type: Tipo de transformación para calcular los intervalos de confianza, “plain”, “log”(defecto), *“log-log”.
conf.int: El nivel de confianza para el intervalo de confianza (.95 por defecto).

Ejemplo

tongue.km <- survfit(Surv(time, delta) ~ 1, data = tongue, type = "kaplan-meier",
error = "tsiatis", conf.type = "log-log", conf.int = 0.99)
summary(tongue.km)

3.3 Graficos de la curva de sobrevivencia.

La curva de sobrevivencia estimada se gráfica con la paquetería survminer, está gráfica esta hecha utilizando la librería ggplot2 y contiene un numero grande de parámetros, por lo que solamente ilustraremos los más importantes.

fit: Objeto tipo survfit.
data: Un conjunto de datos utilizado para ajustar curvas de sobrevivencia.
fun: Transformación de la curva de sobrevivencia(Opcional), las posibles opciones son: “event” para los eventos acumulados, “cumhaz” para el riesgo acumulado y “pct” para la curva de sobrevivencia en porcentaje.
conf.int: Indicador para graficar los intervalos de confianza.
title: Título
xlab: Eje x
ylab: Eje Y
legends.lab: Vector de nombres para identificar las curvas.
legend.title: Título de la leyenda.

Ejemplo

library(survminer)
ggsurvplot(fit = tongue.km, data = tongue, conf.int = TRUE,

title = "Curva de Sobrevivencia",
xlab = "Tiempo",
ylab = "Probabilidad de sobrevivencia",
legend.title = "Estimación",)

3.4 Estimación de la función de riesgo acumulado.

La estimación de la función de riesgo acumulado se puede estimar de dos formas, la primera es por Máxima Verosimilitud: tomando el menos logaritmo de la función de sobrevivencia y el segundo es utilizando el estimador Nelson-Aalen.

3.5 Estimación Máxima Verosimilitud.

La información del Tiempo y Probabilidad de sobrevivencia es extraída del objeto survfit, se calcula el riesgo acumulado utilizando mutate y finalmente se guarda en la variable R.

Las funciones %>%, group_by, mutate, pertenecen a la librería dplyr, una librería para la gramática de manipulación de datos.

3.5.1 Ejemplo Sin Covariables.

tongue.km <- survfit(Surv(time, delta) ~ type, data = tongue)

#R <- tongue.km %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))

#Convertir el objeto survfit a un dataframe
tongue_ovarian <- tidy(tongue.km)

#Variable CumHaz
R <- tongue_ovarian %>%
  mutate(CumHaz = cumsum(estimate))
R

## # A tibble: 69 × 10
##     time n.risk n.event n.censor estimate std.error conf.high conf.low strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>     <dbl>    <dbl> <chr> 
##  1     1     52       1        0    0.981    0.0194     1        0.944 type=1
##  2     3     51       2        0    0.942    0.0343     1        0.881 type=1
##  3     4     49       1        0    0.923    0.0400     0.998    0.853 type=1
##  4    10     48       1        0    0.904    0.0452     0.988    0.827 type=1
##  5    13     47       2        0    0.865    0.0547     0.963    0.777 type=1
##  6    16     45       2        0    0.827    0.0634     0.936    0.730 type=1
##  7    24     43       1        0    0.808    0.0677     0.922    0.707 type=1
##  8    26     42       1        0    0.788    0.0718     0.908    0.685 type=1
##  9    27     41       1        0    0.769    0.0760     0.893    0.663 type=1
## 10    28     40       1        0    0.75     0.0801     0.877    0.641 type=1
## # ℹ 59 more rows
## # ℹ 1 more variable: CumHaz <dbl>

3.5.2 Ejemplo Con Covariables

tongue_ovarian <- tidy(tongue.km)

#Variable CumHaz agrupado por estratos
R <- tongue_ovarian %>%
group_by(strata) %>%
mutate(CumHaz = cumsum(n.event / n.risk))

print(R, 60)

## # A tibble: 69 × 10
## # Groups:   strata [2]
##     time n.risk n.event n.censor estimate std.error
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>
##  1     1     52       1        0    0.981    0.0194
##  2     3     51       2        0    0.942    0.0343
##  3     4     49       1        0    0.923    0.0400
##  4    10     48       1        0    0.904    0.0452
##  5    13     47       2        0    0.865    0.0547
##  6    16     45       2        0    0.827    0.0634
##  7    24     43       1        0    0.808    0.0677
##  8    26     42       1        0    0.788    0.0718
##  9    27     41       1        0    0.769    0.0760
## 10    28     40       1        0    0.75     0.0801
## # ℹ 59 more rows
## # ℹ 4 more variables: conf.high <dbl>, conf.low <dbl>,
## #   strata <chr>, CumHaz <dbl>

La estimación queda guardada en la variable R, y se puede graficar directamente usando el objeto survfit utilizando plot() o ggsurvplot.

plot(tongue.km, fun = "cumhaz", conf.int = F, main = "Riesgo Acumulado", col = 1:2,
xlab = "Tiempo", ylab = "Riesgo Acumulado")

ggsurvplot(tongue.km, fun = "cumhaz", xlab = "Tiempo", censor = T,
           ylab = "Riesgo Acumulado", title = "Riesgo Acumulado",
           legend.title = "Perfil de ADN tumoral",
legend.labs = c("Aneuploide", "Diploide"))

3.6 Estimación Nelson-Aalen.

La información es extraída de la misma manera que la Estimación Máxima Verosimilitud, con la diferencia que se calcula la suma acumulada del número de eventos entre las personas en riesgo.

La gráfica de riesgo acumulado se realiza manualmente utilizando la función qplot.

3.6.1 Ejemplo sin covariables.

tongue_ovarian <- tidy(tongue.km)
#Crear variable CumHaz (
R <- tongue_ovarian %>%
mutate(CumHaz = cumsum(n.event / n.risk))
R

## # A tibble: 69 × 10
##     time n.risk n.event n.censor estimate std.error conf.high conf.low strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>     <dbl>    <dbl> <chr> 
##  1     1     52       1        0    0.981    0.0194     1        0.944 type=1
##  2     3     51       2        0    0.942    0.0343     1        0.881 type=1
##  3     4     49       1        0    0.923    0.0400     0.998    0.853 type=1
##  4    10     48       1        0    0.904    0.0452     0.988    0.827 type=1
##  5    13     47       2        0    0.865    0.0547     0.963    0.777 type=1
##  6    16     45       2        0    0.827    0.0634     0.936    0.730 type=1
##  7    24     43       1        0    0.808    0.0677     0.922    0.707 type=1
##  8    26     42       1        0    0.788    0.0718     0.908    0.685 type=1
##  9    27     41       1        0    0.769    0.0760     0.893    0.663 type=1
## 10    28     40       1        0    0.75     0.0801     0.877    0.641 type=1
## # ℹ 59 more rows
## # ℹ 1 more variable: CumHaz <dbl>

qplot(time, CumHaz, data = R, geom = "step", xlab = "Tiempo (Semanas)", ylab = "Riesgo Acumulado",
main = "Riesgo Acumulado")

## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

3.6.2. Ejemplo con covariables.

tongue_ovarian <- tidy(tongue.km)

# Crear la variable CumHaz (Hazard acumulado) agrupando por estratos
R <- tongue_ovarian %>%
group_by(strata) %>%
mutate(CumHaz = cumsum(n.event / n.risk))
R

## # A tibble: 69 × 10
## # Groups:   strata [2]
##     time n.risk n.event n.censor estimate std.error conf.high conf.low strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>     <dbl>    <dbl> <chr> 
##  1     1     52       1        0    0.981    0.0194     1        0.944 type=1
##  2     3     51       2        0    0.942    0.0343     1        0.881 type=1
##  3     4     49       1        0    0.923    0.0400     0.998    0.853 type=1
##  4    10     48       1        0    0.904    0.0452     0.988    0.827 type=1
##  5    13     47       2        0    0.865    0.0547     0.963    0.777 type=1
##  6    16     45       2        0    0.827    0.0634     0.936    0.730 type=1
##  7    24     43       1        0    0.808    0.0677     0.922    0.707 type=1
##  8    26     42       1        0    0.788    0.0718     0.908    0.685 type=1
##  9    27     41       1        0    0.769    0.0760     0.893    0.663 type=1
## 10    28     40       1        0    0.75     0.0801     0.877    0.641 type=1
## # ℹ 59 more rows
## # ℹ 1 more variable: CumHaz <dbl>

qplot(time, CumHaz, col = strata, data = R, geom = "step", xlab = "Tiempo (Semanas)",
ylab = "Riesgo Acumulado", main = "Riesgo Acumulado")

Si solo se desea ver la estimación del riesgo acumulado en un gráfico, la función plot() o ggsurvplot bajo el argumento fun=“cumhaz” hace esto posible, esta estimación utiliza el método Máxima Verosimilitud.

3.7 Estimación de la media, mediana y percentiles de los tiempos de Sobrevivencia.

La estimación de la media se realiza integrando la función sobrevivencia estimada hasta un valor t. En R la estimación se realiza con la función print() usando como argumento una estimación survfit, bajo los siguientes argumentos:

print.rmean: Indicador para mostrar la estimación de la media (TRUE/FALSE).
rmean: valor de t (Opcional).

Ejemplo:

print(tongue.km, print.rmean = TRUE)

## Call: survfit(formula = Surv(time, delta) ~ type, data = tongue)
## 
##         n events rmean* se(rmean) median 0.95LCL 0.95UCL
## type=1 52     31  146.6      27.7     93      67      NA
## type=2 28     22   86.8      24.0     42      23     112
##     * restricted mean with upper limit =  400

Tanto como para la estimación de la mediana, tanto como los percentiles se utiliza la función quantile(), los argumentos de la función son los siguientes:

x : Objeto tipo Survfit.
probs : Vector de cuantiles que se desea estimar.
conft.int* : Indicador(TRUE/FALSE) para estimar el intervalo de confianza para los cuantiles.

NOTA: La estimación del intervalo de confianza se utiliza en el mismo nivel de confianza que se usó en el survfit.

Ejemplo:

quantile(tongue.km, c(0.05, 0.5, 0.95))

## $quantile
##         5 50 95
## type=1  3 93 NA
## type=2  3 42 NA
## 
## $lower
##         5 50  95
## type=1  1 67  NA
## type=2  1 23 181
## 
## $upper
##         5  50 95
## type=1 16  NA NA
## type=2 12 112 NA

3.8 Prueba de hipótesis para igualdad de dos o más funciones de sobrevivencia.

Estas pruebas de hipótesis se realizan en R utilizando la función survdiff() bajo los siguientes argumentos:

formula: Un objeto fórmula y ~ x, que debe tener un objeto Surv como la respuesta a la izquierda del ~ y, si se desea, las covariables por la derecha. Uno de los términos puede ser un objeto strata. Para una sola curva de sobrevivencia del lado derecho se coloca ~1.
data: Objeto data.frame donde están los datos.
rho: Un valor numérico que asigna los pesos de acuerdo a $S(t)^p$. Estas pruebas están basadas en la familia G-rho de Harrington and Fleming. De acuerdo a los valores particulares que se le asigna a rho, se tienen las siguientes pruebas de hipótesis:
- rho = 0, Prueba log-rank (Defecto).
- rho = 1, Prueba Peto & Peto de la modificación de la prueba Gehan-Wilcoxon.

Ejemplo:

Se utilizará la base de datos “Ovarian” que contiene 26 observaciones de 6 variables, trata la sobrevivencia en un ensayo aleatorizado comparando dos tratamientos para el cáncer de ovario:

futime: Tiempo de sobrevivencia o de censura.
delta: Indicador (0=Vivo, 1=Muerte).
state: Estado del deceso (1, 2, 3 y 4).
age: Edad al diagnóstico de cáncer de laringe.
diagyr: Año del diagnóstico del cáncer.

library(survival)
data("ovarian")
head(ovarian)

##   futime fustat     age resid.ds rx ecog.ps
## 1     59      1 72.3315        2  1       1
## 2    115      1 74.4932        2  1       1
## 3    156      1 66.4658        2  1       2
## 4    421      0 53.3644        2  2       1
## 5    431      1 50.3397        2  1       1
## 6    448      0 56.4301        1  1       2

Utilizando esta información se compara si existe alguna diferencia de las curvas de sobrevivencia entre los estados del deceso.

survdiff(Surv(futime, fustat) ~ resid.ds, data = ovarian, rho = 0) #Prueba log-rank

## Call:
## survdiff(formula = Surv(futime, fustat) ~ resid.ds, data = ovarian, 
##     rho = 0)
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## resid.ds=1 11        3     6.26      1.70      3.62
## resid.ds=2 15        9     5.74      1.85      3.62
## 
##  Chisq= 3.6  on 1 degrees of freedom, p= 0.06

survdiff(Surv(futime, fustat) ~ resid.ds, data = ovarian, rho = 1) #Prueba Peto & Peto

## Call:
## survdiff(formula = Surv(futime, fustat) ~ resid.ds, data = ovarian, 
##     rho = 1)
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## resid.ds=1 11     1.95     4.79      1.68      4.31
## resid.ds=2 15     7.45     4.61      1.75      4.31
## 
##  Chisq= 4.3  on 1 degrees of freedom, p= 0.04

3.9 Comparaciones Múltiples por parejas.

Las combinaciones o pruebas posibles se realizan con la función pairwise_survdiff(), esta devuelve una matriz con los p-valores de las comparaciones a pares. Es importante asegurar que la covariable a evaluar sea de tipo factor.

Los argumentos son los siguientes:

formula: Un objeto fórmula y ~ x, que debe tener un objeto Surv como la respuesta a la izquierda del ~ y, la covariable por la derecha.
data: Objeto data.frame donde están los datos.
p.adjust.method: Método para ajustar los p-valores. Valores incluidos son "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none". Si no se desea ajustar el p-valor (Práctica no recomendada), use p.adjust.method = "none".
rho: Un valor numérico que asigna los pesos de acuerdo a $S(t)^p$.

Ejemplo:

En la base de datos ovarian, la covariable resid.ds no es de tipo factor por lo que se tiene que convertir en factor() antes de utilizar la prueba.

class(ovarian$resid.ds) #Tipo Numérico

## [1] "numeric"

ovarian$resid.ds <- factor(ovarian$resid.ds)
class(ovarian$resid.ds) #Tipo factor

## [1] "factor"

pairwise_survdiff(Surv(futime, fustat) ~ resid.ds, data = ovarian, p.adjust.method = "bonferroni",
rho = 1)

## 
##  Pairwise comparisons using Peto & Peto test 
## 
## data:  ovarian and resid.ds 
## 
##   1    
## 2 0.038
## 
## P value adjustment method: bonferroni

3.9.1 Prueba Estratificada.

La prueba estratificada se hace similarmente a la anterior con la distinción que se le agrega +strata(estrato) en la formula.

Ejemplo:

Para nuestro ejemplo utilizaremos la base de datos ovarian y crearemos una nueva variable llamada GrupoEdad en donde se agrupe a las personas en 3 tipos de edades: 39-54, 54-60, 60+ utilizando la función cut(), esta variable será utilizada como estrato.

ovarian$GrupoEdad <- cut(ovarian$age, breaks = c(39, 54, 60, 90)) #Nueva variable
survdiff(Surv(futime, fustat) ~ resid.ds + strata(GrupoEdad), data = ovarian, rho = 1)

## Call:
## survdiff(formula = Surv(futime, fustat) ~ resid.ds + strata(GrupoEdad), 
##     data = ovarian, rho = 1)
## 
## n=25, 1 observation deleted due to missingness.
## 
##             N Observed Expected (O-E)^2/E (O-E)^2/V
## resid.ds=1 11     1.71     3.48     0.893      2.14
## resid.ds=2 14     7.17     5.41     0.574      2.14
## 
##  Chisq= 2.1  on 1 degrees of freedom, p= 0.1

3.9.2 Prueba de Tendencia.

La función comp() de la librería survMisc realiza comparaciones entre curvas de sobrevivencia con diferentes métodos. Cuando el número de grupos es mayor a 2, esta función realiza la prueba de tendencia. Los argumentos de la función son los siguientes:

x: la función ten() aplicado a una formula Surv.
p: Valor de p para la prueba Fleming-Harrington.
q: Valor de q para la prueba Fleming-Harrington.
scores: Un vector de pesos para la prueba de tendencia.

Los argumentos p y q pueden ser ignorados, el argumento scores si no se le es indicado toma como valores un vector 1 hasta n, donde n es el numero de grupos.

La función devuelve como valores dos tablas, la primera una serie de pruebas de hipótesis de igualdad de curvas mientras que la segunda una serie de pruebas de hipótesis de tendencia de curvas.

Ejemplo:

library(survMisc)
model <- ten(Surv(futime, fustat) ~ resid.ds, data = ovarian)
comp(model, scores = c(1, 2, 3, 4))

##                     Q        Var      Z pNorm
## 1          3.2607e+00 2.9343e+00 1.9035     5
## n          7.2000e+01 1.1610e+03 2.1131     1
## sqrtN      1.5321e+01 5.6551e+01 2.0374     4
## S1         2.7303e+00 1.7227e+00 2.0802     3
## S2         2.6122e+00 1.5615e+00 2.0904     2
## FH_p=1_q=1 3.5717e-01 8.3477e-02 1.2362     6
##               maxAbsZ        Var      Q pSupBr
## 1          3.9577e+00 2.9343e+00 2.3104      5
## n          8.0000e+01 1.1610e+03 2.3479      4
## sqrtN      1.7682e+01 5.6551e+01 2.3513      3
## S1         3.1070e+00 1.7227e+00 2.3672      2
## S2         2.9587e+00 1.5615e+00 2.3677      1
## FH_p=1_q=1 5.2755e-01 8.3477e-02 1.8259      6

Capítulo 4 Estimación Paramétrica de la Función de Sobrevivencia.

La función de sobrevivencia puede ser estimada por medio de distribuciones paramétricas con la función flexsurvreg() de la paquetería flexsurv. Los parámetros son estimados por máxima verosimilitud utilizando los algoritmos disponibles en la función de optim() de R.

Los parámetros definidos como positivos son estimados en la escala logarítmica. Los intervalos de confianza se calculan a partir de la matriz Hessiana, y se transforman de nuevo a la escala original de los parámetros.

formula: Un objeto fórmula y ~ 1, que debe tener un objeto Surv como variable respuesta y a la izquierda del ” ~” y del lado derecho se coloca ~1.

data: Objeto data frame donde están los datos.
cl: Ancho de los intervalos de confianza simétricos para estimaciones de máxima verosimilitud, por defecto 0.95.
dist: Nombre de la distribución a ajustar, en la siguiente tabla se muestran las posibles distribuciones que por defecto puede ajustar.

Tabla 2

Nombre de la distribución	Función	Número de parámetros
Gamma generalizada (estable)	`gengamma`	3
Gamma generalizada (original)	`gengamma.orig`	3
F generalizada (estable)	`genf`	4
F generalizada (original)	`genf.orig`	4
Weibull	`weibull`	2
Gamma	`gamma`	2
Exponencial	`exp`	1
Log-logística	`llogis`	2
Log-normal	`lnorm`	2
Gompertz	`gompertz`	2

Las distribuciones tal como la Weibull, Gamma, Exponencial, etc. en su gran mayoría distribuciones pro- gramadas en R cuentan con sus respectivas funciones para calcular su densidad, su función de distribución, su función de cuantiles y la generación aleatoria de acuerdo a su distribución.

Estas funciones tienen como nombre la distribución y se le agrega al inicio la letra según la función deseada.

p. Función de Distribución.
d. Función de Densidad.
q. Función Cuantiles.
r. Generación Aleatoria.

Ejemplo:

pgamma(10, shape = 5, scale = 4) #Función de Distribucion Gamma

## [1] 0.108822

dweibull(1, shape = 0.5, scale = 2) #Función de Densidad Weibull

## [1] 0.1743261

qlogis(0.5, location = 1, scale = 3) #Función Cuantil de log-Logística

## [1] 1

rgompertz(n = 10, shape = 1, rate = 1/50) #Función de Gompertz

##  [1] 1.7467062 4.4651631 4.0142776 4.2894442 4.9889105 0.8915778 3.0524591
##  [8] 4.8136235 3.6007778 2.8415390

4.1.1 Ejemplo.

Retomando la base de datos anterior “ovarian”.

library(survival)
data("ovarian")

## Warning in data("ovarian"): data set 'ovarian' not found

head(ovarian)

##   futime fustat     age resid.ds rx ecog.ps GrupoEdad
## 1     59      1 72.3315        2  1       1   (60,90]
## 2    115      1 74.4932        2  1       1   (60,90]
## 3    156      1 66.4658        2  1       2   (60,90]
## 4    421      0 53.3644        2  2       1   (39,54]
## 5    431      1 50.3397        2  1       1   (39,54]
## 6    448      0 56.4301        1  1       2   (54,60]

flex <- flexsurvreg(Surv(futime, fustat) ~ 1, data = ovarian, dist = "exp") #Ajuste exponencial
flex

## Call:
## flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, 
##     dist = "exp")
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  0.000770  0.000437  0.001356  0.000222
## 
## N = 26,  Events: 12,  Censored: 14
## Total time at risk: 15588
## Log-likelihood = -98.0322, df = 1
## AIC = 198.0644

La salida de la función son las siguientes:

est: La estimación máxima verosimilitud.
L y U: Intervalos de confianza inferior y superior.
se: La desviación estándar de la estimación.
N: Numero de datos.
Events: Numero de datos observados.
Censored: Numero de datos censurados.
Total time at risk: Tiempo total en riesgo.
Log-likelihood: Log-verosimilitud.
ovarian: Grados de libertad.
AIC: Índice AIC.

4.2 Graficación del ajuste.

4.2.1 R Base.

Se puede comparar el ajuste paramétrico con el ajuste no-paramétrico utilizando la función plot(), con los siguientes argumentos:

flexsurv: Objeto flexsurvreg.
type: tipo de curva a comprar, “survival” para Kaplan-Meier, “cumhaz” para riesgo acumulado y “hazard” para la función de riesgo.
ci: Indicador (TRUE/FALSE) para mostrar las bandas de confianza para distribución ajustada.
conf.int: Indicador(TRUE/FALSE) para mostrar las bandas de confianza para la estimación no paramétrica para el caso de la función de sobrevivencia o riesgo acumulado.
cl: Nivel de confianza para el intervalo.
col: Color de la curva paramétrica.
col.obs: Color de la curva no paramétrica.
main: Titulo del gráfico.
xlab: Nombre del eje x.
ylab: Nombre del eje y.

Ejemplo:

plot(flex, type = "cumhaz", ci = F, main = "Riesgo Acumulado", conf.int = F,
col = 3, col.obs = 4, xlab = "Tiempo", ylab = "Riesgo Acumulado")

plot(flex, type = "survival", ci = F, main = "Probabilidad de sobrevivencia",
conf.int = F, col = 1, col.obs = 2,
xlab = "Tiempo", ylab = "Probabilidad de sobrevivencia")

4.2.2 Ggplot2.

Por el momento no se tiene implementando ningún código para visualizar este gráfico en ggplot2, sin embargo, con el siguiente script para realizar la comparación entre el ajuste paramétrico y ajuste no-paramétrico en la Curva de sobrevivencia y Riesgo Acumulado.

El script calcula las estimaciones de la probabilidad de sobrevivencia (O Riesgo Acumulado) mediante la función summary() (Paramétrico) y fortify() o tidy() (No Paramétrico) en caso de no funcionar fortify según la versión de R utilizada y los gráfica mediante geom_line() y geom_step() respectivamente.

Los cambios que se realizan en el script son los argumentos de la función flexsurvreg, survfit.

4.2.3 Curva de Sobrevivencia.

library(survival) # Para survfit y Surv
library(flexsurv) # Para flexsurvreg
library(broom) # Para tidy
library(ggplot2) # Para graficar
library(dplyr) # Para manipulación de datos

#Ajuste modelo paramétrico
flexgg <- flexsurvreg(Surv(futime, fustat) ~ 1, data = ovarian, dist = "exp") %>%
summary(type = "survival") %>%
as.data.frame()

#Renombrar columnas para asegurar compatibilidad con ggplot
names(flexgg) <- c("time", "est", "lower", "upper")

#Ajuste del modelo no paramétrico y conversión con tidy
kmgg <- survfit(Surv(futime, fustat) ~ 1, data = ovarian) %>%
tidy()

#Comprobar nombres de columnas
print(names(flexgg))

## [1] "time"  "est"   "lower" "upper"

print(names(kmgg))

## [1] "time"      "n.risk"    "n.event"   "n.censor"  "estimate"  "std.error"
## [7] "conf.high" "conf.low"

#
kmgg <- survfit(Surv(futime, fustat) ~ 1, data = ovarian) %>% tidy()
#
if (!"surv" %in% names(kmgg)) {
kmgg$surv <- kmgg$estimate # Asumiendo que 'estimate' es el equivalente a 'surv'
}
#
ggplot() +
geom_line(aes(x = time, y = est, col = "Paramétrico"), data = flexgg) +
geom_step(aes(x = time, y = surv, col = "No Paramétrico"), data = kmgg) +
labs(x = "Time",
y = "Probabilidad de Sobrevivencia",
colour = "Ajustes",
title = "Comparación entre curvas de Sobrevivencia") +
scale_color_manual(values = c("Paramétrico" = "blue", "No Paramétrico" = "red")) +
theme_minimal()

4.2.4 Riesgo Acumulado.

La gráfica de riesgo acumulado se realiza utilizando el mismo script utilizando y=−log(Sest(t))

library(survival) # Para survfit y Surv
library(flexsurv) # Para flexsurvreg
library(broom) # Para tidy
library(ggplot2) # Para graficar
library(dplyr) # Para manipulación de datos

#Ajuste del modelo paramétrico
flexgg <- flexsurvreg(Surv(futime, fustat) ~ 1,
                      data = ovarian, dist = "exp") %>%
summary(type = "survival") %>%
as.data.frame()

#Renombrar columnas para asegurar compatibilidad con ggplot
names(flexgg) <- c("time", "est", "lower", "upper")

#Ajuste del modelo no paramétrico y conversión con tidy
kmgg <- survfit(Surv(futime, fustat) ~ 1, data = ovarian) %>%
tidy()

#Nombres de columnas
print(names(flexgg))

## [1] "time"  "est"   "lower" "upper"

print(names(kmgg))

## [1] "time"      "n.risk"    "n.event"   "n.censor"  "estimate"  "std.error"
## [7] "conf.high" "conf.low"

#
kmgg <- survfit(Surv(futime, fustat) ~ 1, data = ovarian) %>% tidy()
#
if (!"surv" %in% names(kmgg)) {
kmgg$surv <- kmgg$estimate
}
#
ggplot() +
geom_line(aes(x = time, y = -log(est), col = "Paramétrico"), data = flexgg) +
geom_step(aes(x = time, y = -log(surv), col = "No Paramétrico"), data = kmgg) +
labs(x = "Time",
y = "Probabilidad de Sobrevivencia",
colour = "Ajustes",
title = "Comparación entre curvas de Sobrevivencia") +
scale_color_manual(values = c("Paramétrico" = "blue", "No Paramétrico" = "red")) +
theme_minimal()

4.2.4 Comparación de Curvas de Sobrevivencia por Grupos.

Dados dos o más grupos, la comparación entre las curvas paramétricas se realiza mediante la función flex- survreg indicando el nombre de la columna de los grupos igualada a un nivel del mismo en el argumento subset.

Ejemplo:

flex <- flexsurvreg(Surv(futime, fustat) ~ 1, data = ovarian, dist = "exp", subset = ecog.ps ==
1)
flex

## Call:
## flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, 
##     subset = ecog.ps == 1, dist = "exp")
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  0.000600  0.000250  0.001442  0.000268
## 
## N = 14,  Events: 5,  Censored: 9
## Total time at risk: 8332
## Log-likelihood = -42.0921, df = 1
## AIC = 86.18421

Bajo este método, el análisis se realiza individualmente para cada nivel de la variable. Otro método para analizar es ajustar una misma distribución a cada grupo. Se deja el siguiente script con los siguientes parámetros modificables:

ovarian: data.frame que contiene los datos.
ecog.ps: Nombre de la covariable de los grupos.
exp: Distribución a ajustar.

El script funciona agrupando la base de datos (ovarian) por la covariable (ecog.ps) y aplicando a cada grupo un ajuste paramétrico (gompertz), guardando los resultados en la variable model.

head(ovarian)

##   futime fustat     age resid.ds rx ecog.ps GrupoEdad
## 1     59      1 72.3315        2  1       1   (60,90]
## 2    115      1 74.4932        2  1       1   (60,90]
## 3    156      1 66.4658        2  1       2   (60,90]
## 4    421      0 53.3644        2  2       1   (39,54]
## 5    431      1 50.3397        2  1       1   (39,54]
## 6    448      0 56.4301        1  1       2   (54,60]

model <- ovarian %>% group_by(resid.ds) %>% do(flex = flexsurvreg(Surv(futime, fustat) ~
1, data = ., dist = "gompertz"))
model

## # A tibble: 2 × 2
## # Rowwise: 
##   resid.ds flex      
##   <fct>    <list>    
## 1 1        <flxsrvrg>
## 2 2        <flxsrvrg>

model$flex

## [[1]]
## Call:
## flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ., dist = "gompertz")
## 
## Estimates: 
##        est        L95%       U95%       se       
## shape   9.00e-04  -2.14e-03   3.94e-03   1.55e-03
## rate    2.38e-04   2.72e-05   2.07e-03   2.62e-04
## 
## N = 11,  Events: 3,  Censored: 8
## Total time at risk: 8243
## Log-likelihood = -26.63719, df = 2
## AIC = 57.27437
## 
## 
## [[2]]
## Call:
## flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ., dist = "gompertz")
## 
## Estimates: 
##        est        L95%       U95%       se       
## shape  -0.000743  -0.003359   0.001874   0.001335
## rate    0.001555   0.000565   0.004281   0.000803
## 
## N = 15,  Events: 9,  Censored: 6
## Total time at risk: 7345
## Log-likelihood = -69.17407, df = 2
## AIC = 142.3481

Se utilizar model$flex[[]] para acceder a los ajustes individuales.

Ejemplo:

model$flex[[1]]

## Call:
## flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ., dist = "gompertz")
## 
## Estimates: 
##        est        L95%       U95%       se       
## shape   9.00e-04  -2.14e-03   3.94e-03   1.55e-03
## rate    2.38e-04   2.72e-05   2.07e-03   2.62e-04
## 
## N = 11,  Events: 3,  Censored: 8
## Total time at risk: 8243
## Log-likelihood = -26.63719, df = 2
## AIC = 57.27437

4.3 Gráficas Curvas No-Paramétricas con las Curvas Paramétricas.

El siguiente script gráfica las curvas no paramétricas con las curvas paramétricas. Únicamente se tiene que realizar primero el ajuste no paramétrico en la variable KM.

Este script se modifica la variable KM de acuerdo a los parámetros deseados.

#Ajuste Curva de Sobrevivencia No-Paramétrica con KM
KM <- survfit(Surv(futime, fustat) ~ rx, data = ovarian)

#Ajuste modelos paramétricos para cada grupo en 'rx'
models_flex <- lapply(unique(ovarian$rx), function(x) {
data_subset <- ovarian[ovarian$rx == x, ]
flexsurvreg(Surv(futime, fustat) ~ 1, data = data_subset, dist = "exp")
})

#
plot(KM, col = 1:length(models_flex), main = "Comparación Paramétrica",
     xlab = "Tiempo",ylab = "Probabilidad de Sobrevivencia")

#Líneas para cada modelo paramétrico
for (i in 1:length(models_flex)) {
curve_data <- summary(models_flex[[i]], type = "survival")
print(curve_data)
}

4.3.1 Riesgo Acumulado.

#Ajuste del modelo Kaplan-Meier con estratificación por 'rx'
KM <- survfit(Surv(futime, fustat) ~ rx, data = ovarian)

#
plot(KM, col = 1:length(KM$strata), main = "Comparación Paramétrica", xlab = "Tiempo",
ylab = "Riesgo Acumulado", fun = "cumhaz")

#Verificar estructura KM
print(KM)

## Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian)
## 
##       n events median 0.95LCL 0.95UCL
## rx=1 13      7    638     268      NA
## rx=2 13      5     NA     475      NA

4.3.2 ggplot2.

Para realizar la gráfica en ggplot2 se deja el siguiente script, los parámetros a modificar son survfit de la variable KM y ejecutar previamente el script de comparación entre curvas.

4.3.3 Curva de Sobrevivencia.

El script extrae de la lista los valores estimados de la curva de sobrevivencia.

KM <- survfit(Surv(futime, fustat) ~ resid.ds, data = ovarian)
plot(KM, col = 1:length(model$flex), main = "Comparación Paramétrica", xlab = "Tiempo",
ylab = "Probabilidad de Sobrevivencia")
for (i in 1:length(model$flex)) lines(model$flex[[i]], ci = F, col = i)

4.3.4 Riesgo Acumulado

KM <- survfit(Surv(futime, fustat) ~ resid.ds, data = ovarian)

plot(KM, col = 1:length(model$flex), main = "Comparación Paramétrica", xlab = "Tiempo",
ylab = "Riesgo Acumulado", fun = "cumhaz")
for (i in 1:length(model$flex)) lines(model$flex[[i]], ci = F, col = i, type = "cumhaz")

4.4 Comparación de Modelos.

La comparación entre los modelos ajustados se puede realizar visualmente entre la curva de sobrevivencia Kaplan-Meier y las curvas paramétricas ajustadas o mediante el índice AIC.

4.4.1 Comparación grafica.

La comparación gráfica se realiza con el siguiente script, únicamente se debe definir en la variable Dist las distribuciones que se desean comparar, el objeto tipo Surv en la variable data.Surv y guardar la base de datos como datosflex, si se desea utilizar por grupos se debe utilizar el argumento subset en la función flexsurvreg o filtrar los datos antes de ejecutar.

datosflex <- ovarian
Dist <- c("exp", "weibull", "llogis", "gompertz")
data.Surv <- Surv(datosflex$futime, datosflex$fustat)
model <- sapply(Dist, function(x) flexsurvreg(data.Surv ~ 1, dist = x), USE.NAMES = T,
simplify = F)
model

## $exp
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  0.000770  0.000437  0.001356  0.000222
## 
## N = 26,  Events: 12,  Censored: 14
## Total time at risk: 15588
## Log-likelihood = -98.0322, df = 1
## AIC = 198.0644
## 
## 
## $weibull
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape     1.108     0.674     1.822     0.281
## scale  1225.419   690.421  2174.979   358.714
## 
## N = 26,  Events: 12,  Censored: 14
## Total time at risk: 15588
## Log-likelihood = -97.9539, df = 2
## AIC = 199.9078
## 
## 
## $llogis
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape     1.376     0.842     2.249     0.345
## scale   837.753   470.311  1492.268   246.769
## 
## N = 26,  Events: 12,  Censored: 14
## Total time at risk: 15588
## Log-likelihood = -97.3547, df = 2
## AIC = 198.7094
## 
## 
## $gompertz
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est        L95%       U95%       se       
## shape  -0.000509  -0.002568   0.001550   0.001050
## rate    0.000930   0.000371   0.002330   0.000436
## 
## N = 26,  Events: 12,  Censored: 14
## Total time at risk: 15588
## Log-likelihood = -97.90959, df = 2
## AIC = 199.8192

Cada modelo se puede acceder de manera individual utilizado $ sobre la variable model.

model$gompertz

## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est        L95%       U95%       se       
## shape  -0.000509  -0.002568   0.001550   0.001050
## rate    0.000930   0.000371   0.002330   0.000436
## 
## N = 26,  Events: 12,  Censored: 14
## Total time at risk: 15588
## Log-likelihood = -97.90959, df = 2
## AIC = 199.8192

4.4.2 R Base

El siguiente script realiza la comparación gráfica utilizando los gráficos bases de R.

#Curva de sobrevivencia

plot(model[[1]], ci = F, conf.int = F, lty = 2, main = "Ajuste Paramétrico",
xlab = "Time", ylab = "Probabilidad de sobrevivencia")
for (i in 2:length(Dist)) plot(model[[i]], ci = F, conf.int = F, add = T, col = i +
1, lty = i)
legend("topright", c("KM", Dist), lty = 1:(length(Dist) + 1), col = 1:(length(Dist) +
1))

#Riesgo Acumulado

plot(model[[1]], ci = F, conf.int = F, lty = 2, main = "Ajuste Paramétrico",
xlab = "Time", ylab = "Riesgo Acumulado", type = "cumhaz")
for (i in 2:length(Dist)) plot(model[[i]], ci = F, conf.int = F, add = T, col = i +
1, lty = i, type = "cumhaz")

legend("topright", c("KM", Dist), lty = 1:(length(Dist) + 1), col = 1:(length(Dist) +
1))

4.4.3 Índice AIC, BIC y Log-verosimilitud.

Los índices AIC, BIC y Log-verosimilitud, se extrae utilizando las funciones: AIC(), BIC() y logLik() tomando como único argumento un objeto flexsurvreg.

model <- flexsurvreg(Surv(futime, fustat) ~ 1, data = ovarian, dist = "lnorm")
AIC(model)

## [1] 198.2435

BIC(model)

## [1] 200.7597

logLik(model)

## 'log Lik.' -97.12174 (df=2)

Cuando se realiza una comparación entre modelos se puede utilizar el siguiente script para generar una tabla con los respectivos índices para cada distribución. Los parámetros a modificar son los datos datosflex así como data.Surv y las distribuciones a comparar en Dist.

datosflex <- ovarian
data.Surv <- Surv(datosflex$futime, datosflex$fustat)
Dist <- c("exp", "weibull", "llogis", "gompertz")
#Inicio Script
model <- sapply(Dist, function(x) flexsurvreg(data.Surv ~ 1, data = datosflex,
dist = x), USE.NAMES = T, simplify = F)
AIC.model <- sapply(model, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = logLik(x)),
simplify = T)
AIC.model

##             exp  weibull   llogis  gompertz
## AIC    198.0644 199.9078 198.7094 199.81919
## BIC    199.3225 202.4240 201.2256 202.33538
## LogLik -98.0322 -97.9539 -97.3547 -97.90959

#Ordenados por menor AIC
AIC.model[, order(AIC.model["AIC", ])]

##             exp   llogis  gompertz  weibull
## AIC    198.0644 198.7094 199.81919 199.9078
## BIC    199.3225 201.2256 202.33538 202.4240
## LogLik -98.0322 -97.3547 -97.90959 -97.9539

4.4.4 Pruebas Gráficas para Modelos Específicos.

Se compara mediante una serie de transformaciones si la estimación no-paramétrica puede ser razonablemente estimada por algún modelo paramétrico. Los argumentos de la función son los siguientes:

dist: Nombre de la distribución (Ver tabla).
KM: Estimación no paramétrica de los datos utilizando survfit().
lm: Indicador (TRUE/FALSE) si se desea ajustar una recta de mínimos cuadrados.

Distribución y sus códigos:

Tabla 3

Tabla de Distribuciones

Distribución	Función
Exponencial	`exp`
Weibull	`weibull`
Gumbel	`gumbel`
Log-Normal	`lnorm`
Gamma	`gamma`
Normal	`norm`
Logística	`logis`
Log-logística	`llogis`
Pareto	`pareto`

La función es un objeto ggplot2 y se le puede agregar títulos y valor a los ejes con labs().

Ejemplo:

#Ajuste no paramétrico con Kaplan-Meier
KM <- survfit(Surv(futime, fustat) ~ 1, data = ovarian)
plot(KM)

#Ajuste paramétrico con una distribución de Weibull
parametric_model <- flexsurvreg(Surv(futime, fustat) ~ 1,
                                data = ovarian, dist = "weibull")
plot(parametric_model)

4.4.5 Gráfico de Bondad de Ajuste.

El gráfico de bondad de ajuste se realiza mediante un gráfico de dispersión entre $1 - S_{KM}(t_i)$ y $F_{EMV}(t_i; \hat{\theta})$, donde $S_{KM}(t_i)$ es la estimación no paramétrica de la curva de sobrevivencia y $F_{EMV}(t_i; \hat{\theta})$ es la estimación de Máxima Verosimilitud suponiendo una distribución específica.

El script extrae estos valores, guardándolos en las variables F_est y KM_est y los gráfica utilizando la función plot() o qplot().

Los parámetros modificables de este script son los argumentos de las funciones flexsurvreg() y survfit().

flex <- flexsurvreg(Surv(futime, fustat) ~ 1, 
                    data = ovarian, dist = "gompertz") # Ajuste Paramétrico
F_est <- 1 - as.data.frame(summary(flex))[, "est"]
KM_est <- survfit(Surv(futime, fustat) ~ 1, ovarian) 

#Ajuste no-paramétrico
KM_est <- 1 - KM_est$surv

plot(F_est, KM_est)

qplot(F_est, KM_est)

4.4.6 Gráfico de Probabilidad Estabilizado.

El gráfico de probabilidad estabilizado se realiza utilizando el anterior script, anexando las transformaciones en ambos ejes: $2\pi \sin^{-1}(\sqrt{F(t_i)})$ en las funciones plot() o qplot().

plot((2/pi) * asin(sqrt(KM_est)), (2/pi) * asin(sqrt(F_est)), xlab = "F-km-stab",
ylab = "F-emv-stab", main = "Gráfica de Probabilidad Estabilizada")

qplot((2/pi) * asin(sqrt(KM_est)), (2/pi) * asin(sqrt(F_est))) + labs(x = "F-km-stab",
y = "F-emv-stab", title = "Gráfica de Probabilidad Estabilizada")

Capítulo 5 Modelo de Riesgos Proporcionales de Cox.

Supongamos que en un estudio que consiste en la comparación de los tiempos de falla de dos grupos. Los pacientes fueron seleccionados aleatoriamente para recibir el tratamiento estándar (grupo 0) y un nuevo tratamiento (grupo 1). La función de tasa de falla del primer grupo será representada por $h_0(t)$ y la del segundo grupo por $h_1(t)$. Asumiendo proporcionalidad entre estas funciones, se tiene que:

\[ \frac{h_1(t)}{h_0(t)} = K \]

Dónde $K$ es la razón de tasa de falla o riesgo relativo, constante para todo tiempo $t$ de acompañamiento del estudio. Si $X$ es la variable indicadora del grupo, dónde

\[ X = \begin{cases} 0 & \text{si es del grupo 0} \\ 1 & \text{si es del grupo 1} \end{cases} \]

Si $K = \exp(\beta x)$, entonces

\[ h(t) = h_0(t) \exp(\beta x) \]

esto es,

\[ h(t) = \begin{cases} h_1(t) = h_0(t) \exp(\beta) & \text{si } x = 1 \\ h_0(t) & \text{si } x = 0 \end{cases} \]

De una forma genérica consideremos $p$ covariables de modo que $x$ es un vector con los componentes $x = (x_1, \ldots, x_p)$. Así que la expresión general del modelo de regresión de Cox será:

\[ h(t) = h_0(t) g(x^\top \beta) \]

Dónde $g$ es una función que debe ser especificada, tal que $g(0) = 1$. El componente no paramétrico, $h_0(t)$, no es especificado y es una función no negativa del tiempo. El componente paramétrico es frecuentemente usado de la siguiente forma multiplicativa:

\[ g(x^\top \beta) = \exp(x^\top \beta) = \exp(\beta_1 x_1 + \cdots + \beta_p x_p) \]

En que $\beta$ es el vector de parámetros asociado a las covariables. Esta forma garantiza que $h(t)$ será siempre positiva. Este modelo es también llamado modelo de riesgos proporcionales, esto debido a que la razón de las tasas de falla de dos individuos diferentes es constante en el tiempo. Esto es, la razón de las tasas de falla para dos individuos diferentes, $i$ y $j$, es:

\[ \frac{h_i(t)}{h_j(t)} = \frac{h_0(t) \exp(x_i \beta)}{h_0(t) \exp(x_j \beta)} \]

que no depende del tiempo. La suposición básica para el uso del modelo de regresión de Cox, por tanto, es que las tasas de falla sean proporcionales.

5.1 Ajustando el modelo de Cox

El modelo de regresión de Cox utiliza coeficientes β para medir los efectos de las covariables en la tasa de falla. Para determinar el modelo, estos coeficientes deben estimarse a partir de datos muestrales.

A menudo se usa el método de máxima verosimilitud; sin embargo, la inclusión de un componente no paramétrico h0(t) en la función de verosimilitud hace que este método no sea el más adecuado.

Cox propuso una solución en su trabajo original, y la formalizó más adelante en 1975, el método de máxima verosimilitud parcial, que condiciona la verosimilitud para excluir el componente perturbador.

5.2 Método de Máxima Verosimilitud Parcial.

Consideremos una muestra de $n$ individuos donde existen $k \leq n$ fallas distintas en los tiempos $t_1 \leq t_2 \leq \ldots \leq t_k$. Además, considerando el siguiente argumento condicional: la probabilidad condicional de que la i-ésima observación venga a fallar en el tiempo $t_i$ dado las observaciones que están en riesgo en $t_i$ es

\[ \frac{h_i(t)}{\sum_{j \in R(t_i)} h_j(t)} = \frac{h_0(t) \exp\{x_i \beta\}}{\sum_{j \in R(t_i)} h_0(t) \exp\{x_j \beta\}} = \frac{\exp\{x_i \beta\}}{\sum_{j \in R(t_i)} \exp\{x_j \beta\}} \]

donde $R(t_i)$ es el conjunto de índices de las observaciones sobre el riesgo en el tiempo $t_i$.

5.3 La Función de Verosimilitud.

La función de verosimilitud a ser utilizada para hacer inferencias en el modelo, es entonces, formada por el producto de todos los términos representados de la ecuación anterior asociados a los tiempos diferentes de falla, esto es,

\[ L(\beta) = \prod_{i=1}^{k} \left( \frac{\exp\{x_i \beta\}}{\sum_{j \in R(t_i)} \exp\{x_j \beta\}} \right)^{\delta_i} \]

con $\delta_i$ el indicador de falla.

Los valores de $\beta$ que maximizan $L(\beta)$ son obtenidos resolviendo el sistema de ecuaciones $U(\beta) = 0$, donde $U(\beta)$ es el vector de escore de las primeras derivadas de la función $l(\beta) = \log(L(\beta))$.

5.4 Interpretación de los Coeficientes.

A partir del modelo general de regresión de Cox, se puede observar que el efecto de las covariables es de acelerar o desacelerar la función de riesgo. Sin embargo, la propiedad de riesgos proporcionales del modelo debe ser usada para interpretar los coeficientes estimados. Tomando la razón de tasas de falla de dos individuos (i y j) que tienen los mismos valores para las covariables con excepción de la l-ésima, se tiene:

\[ \frac{h_i(t)}{h_j(t)} = \exp\{\beta_l (x_{il} - x_{jl})\} \]

que puede ser interpretado como la razón de riesgos o de riesgo relativo instantáneo en el tiempo t. Sin embargo, como esta razón es constante para todo el acompañamiento, se puede suprimir la palabra “instantánea” de la interpretación.

Por ejemplo, suponga que $x_l$ sea una covariable dicotómica indicando pacientes hipertensos. El riesgo de muerte entre los hipertensos es $\exp\{\beta_l\}$ veces el riesgo de pacientes con presión normal, manteniendo fijas las otras covariables.

5.5 Estimando Funciones Relacionadas a h0(t).

Función de Tasa de Falla Acumulada

\[ H_0(t) = \int_{0}^{t} h_0(u) \, du \]

y la correspondiente función de sobrevivencia

\[ S_0(t) = \exp(-H_0(t)) \]

Una estimación para $H_0(t)$ es la propuesta por Breslow (1972), es una función escalonada con saltos en tiempos diferentes de falla y expresada por:

\[ H_0(t_i) = \sum_{j : t_j \leq t} \frac{d_j}{\sum_{l \in R_j} \exp(x_l \beta)} \]

con $d_j$ el número de fallas en $t_j$.

Consecuentemente, las funciones de sobrevivencia $S_0(t)$ y $S(t)$ pueden ser estimadas por:

\[ S(t) = \exp[-H_0(t)] \]

\[ S(t) = [S_0(t)]^{\exp(x \beta)} \]

Tanto $S_0(t)$ como $S(t)$ son funciones escalonadas decrecientes en el tiempo. Notemos que en la ausencia de covariables, la expresión de $H_0(t_i)$ se reduce a

\[ H_0(t_i) = \sum_{j : t_j \leq t} \frac{d_j}{n_j} \]

que es el estimador de Nelson-Aalen. Notemos que por ese hecho, el estimador $H_0(t_i)$ es también referenciado en la literatura como estimador de Nelson-Aalen-Breslow.

5.6 Elección del Modelo de Cox.

El modelo de regresión de Cox es bastante flexıble debido a la presencia del componente no paramétrico. Incluso así, el no se ajusta a cualquier situación clínica y como cualquier otro modelo estadístico, requiere de técnicas para validar su elección.

En particular, como mencionado anteriormente, el tiene una suposición básica que es la de riesgos propor- cionales. La violación de esta suposición puede acarriar serios vicios en la estimación de los coeficientes del modelo (Struthers e Kalbfleisch, 1986).

5.7 Verificando la suposición de riesgos proporcionales.

5.7.1 Métodos Gráficos.

Para verificar la suposición de riesgos proporcionales en el modelo de Cox, un gráfico simple y bastante usado es propuesto para esta finalidad. La obtención de este gráfico consiste, inicialmente, en dividir los datos en m estratos, usualmente de acuerdo con alguna covariable.

Por ejemplo, dividir los datos en dos estratos de acuerdo a la covariable sexo. En seguida, se debe estimar $H_{0j}(t)$ para cada estrato. Si la suposición fuese válida, las curvas del logaritmo de $H_{0j}(t)$ versus $t$, o $\log(t)$, deben presentar diferencias aproximadamente constantes en el tiempo. Curvas no paralelas significan una contradicción en la suposición de riesgos proporcionales. Es razonable construir este gráfico para cada covariable incluida en el estudio. Si la covariable fuese de naturaleza continua, una sugerencia es agruparla en un pequeño número de categorías.

Una ventaja de esta técnica gráfica es la de indicar la covariable que está generando una violación de la suposición, en caso de que eso ocurra. Una propuesta adicional que viene siendo usada para verificar la suposición de los riesgos proporcionales en el modelo de Cox es analizar los residuos de Schoenfeld (1982).

Para definir los residuos, consideremos que existen $k \leq n$ tiempos diferentes de falla $t_1, \ldots, t_k$. Si el individuo $i$ con vector de covariables $x_i = (x_{1i}, \ldots, x_{2i}, \ldots, x_{pi})$ falló, entonces se tiene para ese individuo un vector de residuos Schoenfeld $r_i = (r_{1i}, \ldots, r_{2i}, \ldots, r_{pi})$ en que cada componente $r_{qi} (q = 1, \ldots, p)$ es definido por:

\[ r_{qi} = x_{qi} - \frac{\sum_{j \in R(t_i)} x_{qj} \exp\{x_j \beta\}}{\sum_{j \in R(t_i)} \exp\{x_j \beta\}} \]

Note que para cada una de las p covariables consideradas en el modelo, se tiene, para el individuo $i$, un correspondiente residuo Schoenfeld. Como los residuos son definidos en cada falla, el conjunto de residuos Schoenfeld es una matriz con k líneas y p columnas. Cada línea corresponde a un tiempo distinto de falla y cada columna a una de las p covariables consideradas en el modelo. La k-ésima línea de esta matriz es obtenida por la ecuación anterior. Los residuos estandarizados están definidos por:

\[ r_i^* = I(\beta)^{-1} r_i \]

donde $I(\beta)$ es la matriz de información de Fisher observada.

5.8 Ejemplo en R

5.8.1 Preprocesamiento de los datos.

Para estimar correctamente el modelo de riesgos proporcionales se requiere que las variables a utilizar cumplan con ciertos requisitos dependiendo su tipo.

Se trabajará con la base de datos lung (package=“survival”). Esta trata sobre la sobrevivencia en pacientes con cáncer de pulmón avanzado del Grupo de Tratamiento del Cáncer Centro Norte. Las puntuaciones de rendimiento califican qué tan bien el paciente puede realizar las actividades diarias habituales. La data contiene 228 observaciones y 10 variables.

datos <- lung
str(datos)

## 'data.frame':    228 obs. of  10 variables:
##  $ inst     : num  3 3 3 5 1 12 7 11 1 7 ...
##  $ time     : num  306 455 1010 210 883 ...
##  $ status   : num  2 2 1 2 2 1 2 2 2 2 ...
##  $ age      : num  74 68 56 57 60 74 68 71 53 61 ...
##  $ sex      : num  1 1 1 1 1 1 2 2 1 1 ...
##  $ ph.ecog  : num  1 0 0 1 0 1 2 2 1 2 ...
##  $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ...
##  $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ...
##  $ meal.cal : num  1175 1225 NA 1150 NA ...
##  $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

head(datos)

##   inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1    3  306      2  74   1       1       90       100     1175      NA
## 2    3  455      2  68   1       0       90        90     1225      15
## 3    3 1010      1  56   1       0       90        90       NA      15
## 4    5  210      2  57   1       1       90        60     1150      11
## 5    1  883      2  60   1       0      100        90       NA       0
## 6   12 1022      1  74   1       1       50        80      513       0

Variables Numéricas.

Para el uso de variables numéricas, no se requiere ningún formato en específico, su uso simplemente consiste en agregarla en la fórmula.

coxph(Surv(time, status) ~ age, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ age, data = datos)
## 
##         coef exp(coef) se(coef)     z      p
## age 0.018720  1.018897 0.009199 2.035 0.0419
## 
## Likelihood ratio test=4.24  on 1 df, p=0.03946
## n= 228, number of events= 165

Las variables numéricas se pueden aplicar transformaciones en la expresión de la formula tales por ejemplo; log, exp, de manera general se usa la función para realizar estas transformaciones.

coxph(Surv(time, status) ~ log(age) + I(age^2) + I(age > 60), data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ log(age) + I(age^2) + I(age > 
##     60), data = datos)
## 
##                       coef  exp(coef)   se(coef)      z     p
## log(age)        -0.3847904  0.6805933  3.1873225 -0.121 0.904
## I(age^2)         0.0002760  1.0002760  0.0004346  0.635 0.525
## I(age > 60)TRUE -0.1844504  0.8315611  0.2863351 -0.644 0.519
## 
## Likelihood ratio test=4.82  on 3 df, p=0.1855
## n= 228, number of events= 165

5.8.3 Variables Nominales u Ordinales.

Para el uso de variables nominales u ordinales se debe codificar cada factor con un numero entero, posteri- ormente de leer los datos se debe indicar las variables que son Nominales u Ordinales.

Este procedimiento se realiza utilizando la función factor() sobre la columna del factor con los siguientes argumentos:

x: Vector de datos.
levels: Un vector opcional de los valores que x pudo haber tomado. El valor predeterminado es el conjunto único de valores tomados por x, ordenados en orden creciente de x.
labels: Un vector de caracteres de etiquetas para los niveles (en el mismo orden de levels).
ordered: Indicador (TRUE/FALSE) para determinar si los niveles deben considerarse con ordenados en el orden indicado en labels(Opcional, por defecto FALSE).

En caso de factores ordenados, las estimaciones en el modelo de Cox se usarán utilizando polinomios ortog- onales, mientras que en factores no ordenado se utilizará variables dummies.

Ejemplo:

En nuestra data lung no tenemos ninguna variable categórica por ello para fines demostrativos transformare- mos la variable ‘sex’ en un factor con etiquetas específicas.

Descripción del Conjunto de Datos

Tabla de Variables

Variable	Descripción
inst	Código de la institución
time	Tiempo de sobrevivencia en días
status	Estado de censura: 1=censurado, 2=fallecido
age	Edad en años
sex	Sexo: Masculino=1, Femenino=2
ph.ecog	Puntuación del estado funcional ECOG por el médico: 0=asintomático, 1=sintomático pero completamente ambulatorio
ph.karno	Puntuación de desempeño Karnofsky (malo=0-bueno=100) evaluada por el médico
pat.karno	Puntuación de desempeño Karnofsky evaluada por el paciente
meal.cal	Calorías consumidas en las comidas
wt.loss	Pérdida de peso en los últimos seis meses (libras)

Sexo: 1- Male, 2- Female. ph.ecog: 0=asintomático, 1=sintomático.

datos$sex <- factor(datos$sex, levels = c(1, 2), labels = c("Male", "Female"))
datos$ph.ecog <- factor(datos$ph.ecog, levels = c(0, 1),

labels = c("asintomático", "sintomático."))


head(datos) #Datos tranformados

##   inst time status age  sex      ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1    3  306      2  74 Male sintomático.       90       100     1175      NA
## 2    3  455      2  68 Male asintomático       90        90     1225      15
## 3    3 1010      1  56 Male asintomático       90        90       NA      15
## 4    5  210      2  57 Male sintomático.       90        60     1150      11
## 5    1  883      2  60 Male asintomático      100        90       NA       0
## 6   12 1022      1  74 Male sintomático.       50        80      513       0

coxph(Surv(time, status) ~ sex, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ sex, data = datos)
## 
##              coef exp(coef) se(coef)      z       p
## sexFemale -0.5310    0.5880   0.1672 -3.176 0.00149
## 
## Likelihood ratio test=10.63  on 1 df, p=0.001111
## n= 228, number of events= 165

coxph(Surv(time, status) ~ ph.ecog, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ ph.ecog, data = datos)
## 
##                       coef exp(coef) se(coef)     z      p
## ph.ecogsintomático. 0.3673    1.4439   0.1987 1.849 0.0645
## 
## Likelihood ratio test=3.56  on 1 df, p=0.05912
## n= 176, number of events= 119 
##    (52 observations deleted due to missingness)

Por defecto se toma como etiqueta base a la etiqueta con código más pequeño. De manera general se puede utilizar la función I( ) para crear las variables dummies, su uso consiste en crear indicadoras.

Ejemplo:

coxph(Surv(time, status) ~ I(sex == "Male") + I(sex == "Female"), data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ I(sex == "Male") + I(sex == 
##     "Female"), data = datos)
## 
##                          coef exp(coef) se(coef)     z       p
## I(sex == "Male")TRUE   0.5310    1.7007   0.1672 3.176 0.00149
## I(sex == "Female")TRUE     NA        NA   0.0000    NA      NA
## 
## Likelihood ratio test=10.63  on 1 df, p=0.001111
## n= 228, number of events= 165

5.8.4 Interacción.

Para realizar la interacción con variables categóricas se requiere que estas variables sean codificadas con la función factor, tal como en el procedimiento anterior. Existen 3 formas principales para manejar la interacción en los modelos de regresión, y son mediante el uso de los siguientes operadores:

El operador * agrega las variables principales y todas las posibles interacciones entre ellas en un modelo.

Por ejemplo, para un modelo de tasa de riesgo donde

\[ \lambda(t \mid X_i) = \lambda_0(t) \exp(\beta_0 \text{Edad} + \beta_1 \text{Sexo} + \beta_3 \text{Edad} \times \text{Sexo}) \]

la fórmula en R de la función podría ser representada de la siguiente manera:

coxph(Surv(time, status) ~ age * sex, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ age * sex, data = datos)
## 
##                    coef exp(coef)  se(coef)      z      p
## age            0.020343  1.020551  0.011395  1.785 0.0742
## sexFemale      0.095322  1.100013  1.228718  0.078 0.9382
## age:sexFemale -0.009689  0.990358  0.019414 -0.499 0.6177
## 
## Likelihood ratio test=14.37  on 3 df, p=0.002441
## n= 228, number of events= 165

Todas las posibles interacciones entre age, sex y ph.ecog:

coxph(Surv(time, status) ~ age * sex * ph.ecog, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ age * sex * ph.ecog, data = datos)
## 
##                                        coef exp(coef)  se(coef)      z      p
## age                                0.055106  1.056653  0.026798  2.056 0.0397
## sexFemale                         -0.176416  0.838269  3.252224 -0.054 0.9567
## ph.ecogsintomático.                3.773425 43.528902  1.976268  1.909 0.0562
## age:sexFemale                     -0.004511  0.995499  0.051763 -0.087 0.9306
## age:ph.ecogsintomático.           -0.052495  0.948859  0.030368 -1.729 0.0839
## sexFemale:ph.ecogsintomático.      0.151380  1.163439  3.634085  0.042 0.9668
## age:sexFemale:ph.ecogsintomático. -0.005490  0.994525  0.058147 -0.094 0.9248
## 
## Likelihood ratio test=20.38  on 7 df, p=0.004799
## n= 176, number of events= 119 
##    (52 observations deleted due to missingness)

El operador : se utiliza para especificar únicamente la interacción entre variables en un modelo es- tadístico, sin incluir los efectos principales de las variables involucradas entre sex, age y las variables principales.

coxph(Surv(time, status) ~ age + sex + age:sex, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ age + sex + age:sex, data = datos)
## 
##                    coef exp(coef)  se(coef)      z      p
## age            0.020343  1.020551  0.011395  1.785 0.0742
## sexFemale      0.095322  1.100013  1.228718  0.078 0.9382
## age:sexFemale -0.009689  0.990358  0.019414 -0.499 0.6177
## 
## Likelihood ratio test=14.37  on 3 df, p=0.002441
## n= 228, number of events= 165

Modelo utilizando solo la interacción entre sex y age:

coxph(Surv(time, status) ~ age:sex, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ age:sex, data = datos)
## 
##                   coef exp(coef) se(coef)     z     p
## age:sexMale   0.019821  1.020019 0.009191 2.157 0.031
## age:sexFemale 0.011624  1.011692 0.009537 1.219 0.223
## 
## Likelihood ratio test=14.37  on 2 df, p=0.0007596
## n= 228, number of events= 165

El operador ^ en un modelo estadístico se utiliza para incluir interacciones entre variables hasta el nivel indicado.

Modelo utilizando las variables principales y sus interacciones dobles y triples de las variables age, meal.cal y sex.

coxph(Surv(time, status) ~ (age + meal.cal + sex)^3, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ (age + meal.cal + sex)^3, 
##     data = datos)
## 
##                              coef  exp(coef)   se(coef)      z     p
## age                     1.607e-02  1.016e+00  3.333e-02  0.482 0.630
## meal.cal               -2.692e-04  9.997e-01  2.010e-03 -0.134 0.893
## sexFemale              -8.997e-01  4.067e-01  3.190e+00 -0.282 0.778
## age:meal.cal            1.171e-06  1.000e+00  3.079e-05  0.038 0.970
## age:sexFemale           3.311e-03  1.003e+00  4.959e-02  0.067 0.947
## meal.cal:sexFemale      8.775e-04  1.001e+00  3.113e-03  0.282 0.778
## age:meal.cal:sexFemale -1.057e-05  1.000e+00  5.043e-05 -0.210 0.834
## 
## Likelihood ratio test=11.04  on 7 df, p=0.1369
## n= 181, number of events= 134 
##    (47 observations deleted due to missingness)

Modelo utilizando las variables principales y las interacciones dobles de las variables age, meal.cal y ph.ecog.

coxph(Surv(time, status) ~ (age + meal.cal + ph.ecog)^2, data = datos)

## Call:
## coxph(formula = Surv(time, status) ~ (age + meal.cal + ph.ecog)^2, 
##     data = datos)
## 
##                                    coef  exp(coef)   se(coef)      z     p
## age                           5.312e-02  1.055e+00  3.637e-02  1.461 0.144
## meal.cal                     -3.464e-04  9.997e-01  1.983e-03 -0.175 0.861
## ph.ecogsintomático.           1.763e+00  5.831e+00  1.982e+00  0.889 0.374
## age:meal.cal                 -5.646e-06  1.000e+00  2.955e-05 -0.191 0.848
## age:ph.ecogsintomático.      -3.879e-02  9.620e-01  2.853e-02 -1.359 0.174
## meal.cal:ph.ecogsintomático.  9.362e-04  1.001e+00  6.787e-04  1.379 0.168
## 
## Likelihood ratio test=8.17  on 6 df, p=0.2261
## n= 136, number of events= 94 
##    (92 observations deleted due to missingness)

De manera general se puede utilizar la función I( ) para crear nuevas variables y/o interacciones deseadas.

5.8.5 Intervalos de Confianza y Pruebas de Significancia.

Los intervalos de confianza y las Pruebas de Significancia se realizan con la función summary() a un modelo previamente ajustado, esta función tiene los siguientes argumentos:

object: Modelo coxph.
conf.int: Nivel de confianza (95% Defecto)

modelo <- coxph(Surv(time, status) ~ age:sex, data = datos)
summary(modelo)

## Call:
## coxph(formula = Surv(time, status) ~ age:sex, data = datos)
## 
##   n= 228, number of events= 165 
## 
##                   coef exp(coef) se(coef)     z Pr(>|z|)  
## age:sexMale   0.019821  1.020019 0.009191 2.157    0.031 *
## age:sexFemale 0.011624  1.011692 0.009537 1.219    0.223  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## age:sexMale       1.020     0.9804     1.002     1.039
## age:sexFemale     1.012     0.9884     0.993     1.031
## 
## Concordance= 0.603  (se = 0.026 )
## Likelihood ratio test= 14.37  on 2 df,   p=8e-04
## Wald test            = 14.05  on 2 df,   p=9e-04
## Score (logrank) test = 14.34  on 2 df,   p=8e-04

Los resultados de la función coxph en R proporcionan una tabla detallada que incluye los coeficientes estimados de $e^{\beta_i}$, junto con sus respectivos intervalos de confianza. Este enfoque ayuda a interpretar más directamente los factores de riesgo en términos de riesgos relativos.

Adicionalmente, se proporcionan estadísticos de Concordancia y $R^2$, que son medidas de la calidad del ajuste del modelo. Estas métricas son útiles para evaluar qué tan bien el modelo predice los datos observados.

Los resultados también incluyen estadísticas de pruebas importantes para la validación del modelo:

Prueba de Razón de Verosimilitud: Evalúa si el modelo con más parámetros (más complejo) proporciona un mejor ajuste a los datos que un modelo más simple.
Prueba de Wald: Se utiliza para probar la significancia de los coeficientes individuales del modelo.
Prueba de Puntajes (Score Test): Examina la derivada de la función de verosimilitud y es útil para evaluar la contribución de los coeficientes al modelo.

Cada una de estas pruebas proporciona un valor de p-valor y los grados de libertad asociados, permitiendo determinar la significancia estadística de los hallazgos.

5.8.6 Comparación entre Modelos.

La comparación entre modelos de regresión se puede realizar utilizando la función anova(), específicamente diseñada para efectuar la prueba de Razón de Verosimilitud entre dos modelos ajustados. Esta prueba es útil para determinar si la inclusión de variables adicionales en un modelo mejora significativamente el ajuste.

Argumentos de la función anova()

Los modelos ajustados que se desean comparar se pasan como argumentos a la función.

Salida de la función

La función anova() produce una tabla que incluye los siguientes elementos para cada uno de los modelos comparados:

Log-verosimilitud (Log-Likelihood): Muestra el logaritmo de la función de verosimilitud para cada modelo, que indica qué tan bien el modelo se ajusta a los datos.
Estadístico Chi-cuadrado (Chisq): Muestra el valor del estadístico Chi-cuadrado para la comparación, que es la medida de cuánto mejora el modelo más complejo en comparación con el más simple.
Grados de Libertad (Df): Indica el número de parámetros adicionales en el modelo más complejo.
P-valor (P(>|Chi|)): Proporciona el p-valor asociado a la prueba, que indica la significancia estadística de las diferencias entre los modelos.

Ejemplo

modeloCompleto <- coxph(Surv(time, status) ~ age + sex, data = datos)
modeloReducido <- coxph(Surv(time, status) ~ age, data = datos)

anova(modeloCompleto, modeloReducido)

## Analysis of Deviance Table
##  Cox model: response is  Surv(time, status)
##  Model 1: ~ age + sex
##  Model 2: ~ age
##    loglik  Chisq Df Pr(>|Chi|)   
## 1 -742.85                        
## 2 -747.79 9.8822  1   0.001669 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.9 Validación de Supuestos.

5.9.1 Riesgos Proporcionales.

La validación de los supuestos de proporcionalidad de riesgos en un modelo de regresión de Cox se lleva a cabo mediante el análisis de los residuos de Schoenfeld. Estos residuos son utilizados para examinar si las covariables interactúan con el tiempo de una manera que podría violar el supuesto de proporcionalidad. La función cox.zph() en R facilita esta validación realizando una prueba de correlación entre los residuos de Schoenfeld y alguna transformación del tiempo.

Argumentos de cox.zph(): - fit: Un objeto coxph resultante de ajustar un modelo de regresión de Cox. - transform: Define el tipo de transformación del tiempo a aplicar en el análisis. Los valores posibles incluyen “km” para la transformación Kaplan-Meier de los tiempos de falla, “rank” para el uso de rangos de los tiempos, “identity” que mantiene los tiempos sin cambios, o puede ser una función definida por el usuario que transforme los tiempos.

Salida de la función cox.zph():

La función produce una tabla que incluye los siguientes elementos para cada covariable en el modelo:

Rho: El coeficiente de correlación de Schoenfeld, que mide la dependencia entre los residuos y el tiempo transformado.
Chisq: El valor del estadístico Chi-cuadrado para la prueba de correlación.
P-valor: El p-valor asociado con el estadístico Chi-cuadrado, que ayuda a determinar si existe una correlación significativa. Además, la función ofrece una prueba global para evaluar la proporcionalidad de los riesgos en todo el modelo.

Ejemplo

modelo <- coxph(Surv(time, status) ~ log(age) * sex, data = datos)
riesgo <- cox.zph(modelo)
riesgo

##              chisq df    p
## log(age)     0.138  1 0.71
## sex          2.631  1 0.10
## log(age):sex 2.696  1 0.10
## GLOBAL       3.346  3 0.34

Para realizar una validación gráfica del supuesto de proporcionalidad de riesgos en un modelo de regresión de Cox, la función ggcoxzph() de la librería survminer es una herramienta efectiva. Esta función permite visualizar los residuos de Schoenfeld en contra del tiempo, facilitando la identificación de posibles violaciones del supuesto de proporcionalidad.

Uso de ggcoxzph():

cox.zph object: Se debe proporcionar un objeto resultante de la función cox.zph().
var: Un vector que indica las covariables cuyos residuos de Schoenfeld se desean graficar. Si se omite, se graficarán todos.

La librería survminer está diseñada específicamente para facilitar la visualización de resultados de análisis de sobrevivencia, incluyendo las pruebas de proporcionalidad de riesgos.

Ejemplo

riesgo

##              chisq df    p
## log(age)     0.138  1 0.71
## sex          2.631  1 0.10
## log(age):sex 2.696  1 0.10
## GLOBAL       3.346  3 0.34

ggcoxzph(riesgo)

ggcoxzph(riesgo, var = c("log(age)", "sex"))

5.9.2 Análisis Residual.

El análisis residual para modelos de regresión de Cox se puede facilitar mediante la función ggcoxdiagnostics de la librería survminer. Esta función proporciona gráficos de diagnóstico útiles para evaluar la adecuación del modelo ajustado con coxph. A continuación se describen los argumentos más relevantes de la función ggcoxdiagnostics:

Argumentos de ggcoxdiagnostics:

fit: Un objeto modelo calculado mediante la función coxph de la librería survival.
type: Define el tipo de residual que se mostrará en el eje Y. Los valores posibles incluyen:
- "martingale": Residuos de martingala.
- "deviance": Residuos de desviación.
- "score": Residuos de score.
- "schoenfeld": Residuos de Schoenfeld.
- "dfbeta": Residuos DFBETA.
- "dfbetas": Residuos DFBETAS.
- "scaledsch": Schoenfeld escalado.
x.scale: Especifica cómo se representará el eje X. Las opciones son:
- "linear.predictions": Predicciones lineales.
- "observation.id": Identificador de la observación.
- "time": Tiempo.
title: Título del gráfico.
subtitle: Subtítulo del gráfico.
caption: Leyenda del gráfico.

Ejemplo

modelo <- coxph(Surv(time, status) ~ log(age) , data = lung)
ggcoxdiagnostics(modelo, type = "schoenfeld", ox.scale = "time")

## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## ℹ Please use `gather()` instead.
## ℹ The deprecated feature was likely used in the survminer package.
##   Please report the issue at <https://github.com/kassambara/survminer/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

ggcoxdiagnostics(modelo, type = "dfbetas", ox.scale = "observation.id",
                 title = "Residuales dfbetas",
subtitle = "dfbetas vs id", caption = "Análisis residual")

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

Al utilizar la función ggcoxdiagnostics para el análisis de residuos en modelos de regresión de Cox, se puede encontrar que por defecto, esta gráfica los residuos para todas las covariables del modelo. Para obtener gráficos individuales por cada covariable o un subconjunto específico de covariables, se puede emplear una función adicional llamada ggcoxdiagnostics_var. Esta función permite filtrar y visualizar los residuos para covariables seleccionadas.

Argumentos de ggcoxdiagnostics_var:

*: Objeto ggplot generado con ggcoxdiagnostics.
var: Vector de nombres de las covariables que se desean mostrar en el gráfico.

Ejemplo

p <- ggcoxdiagnostics(modelo, type = "dfbetas", ox.scale = "observation.id",
title = "Residuales dfbetas", subtitle = "dfbetas vs id", caption = "Análisis residual")

p

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

6 Predicción de Curva de Sobrevivencia.

La función survfit se utiliza en el contexto de la sobrevivencia para calcular la función de sobrevivencia pronosticada basándose en un modelo de riesgos proporcionales de Cox ajustado previamente con la función coxph. A continuación se describen los argumentos principales de survfit:

Argumentos de survfit:

formula: Un objeto modelo de Cox, el resultado de coxph.
newdata: Un dataframe que contiene las mismas variables utilizadas en el modelo ajustado. La función de sobrevivencia que se calcula representa la cohorte definida por los valores específicos en newdata.
conf.int: Nivel del intervalo de confianza, con un valor predeterminado de 0.95.
se.fit: Indicador booleano para determinar si se deben calcular los errores estándares de las estimaciones de la sobrevivencia.
conf.type: Define el método para calcular los intervalos de confianza. Las opciones incluyen "none", "plain", "log" (predeterminado), y "log-log".

Salida de survfit:*

La salida incluye estimaciones de la función de sobrevivencia que son comparables a las obtenidas mediante el estimador de Kaplan-Meier y los métodos de Fleming-Harrington. Esto incluye:

Estimaciones de la probabilidad de sobrevivencia en cada tiempo observado.
Intervalos de confianza para estas estimaciones, dependiendo de conf.type.
Los tiempos y los eventos asociados a cada estimación.

Ejemplo

(modelo <- coxph(Surv(time, status) ~ age + meal.cal, data = datos))

## Call:
## coxph(formula = Surv(time, status) ~ age + meal.cal, data = datos)
## 
##                coef  exp(coef)   se(coef)      z      p
## age       1.892e-02  1.019e+00  1.044e-02  1.813 0.0699
## meal.cal -3.217e-05  1.000e+00  2.349e-04 -0.137 0.8911
## 
## Likelihood ratio test=3.67  on 2 df, p=0.1592
## n= 181, number of events= 134 
##    (47 observations deleted due to missingness)

datosnuevos <- data.frame(age = c(50, 60), meal.cal = c(1500, 2000)) #Datos a Estimar
datosnuevos

##   age meal.cal
## 1  50     1500
## 2  60     2000

curva <- survfit(modelo, newdata = datosnuevos)
curva

## Call: survfit(formula = modelo, newdata = datosnuevos)
## 
##     n events median 0.95LCL 0.95UCL
## 1 181    134    361     291     558
## 2 181    134    329     239     574

Las curvas de sobrevivencia estimadas generadas por la función survfit en R son objetos del tipo Survfit. Estos objetos permiten realizar análisis adicionales, como la estimación de medias y cuantiles de las funciones de sobrevivencia. Se pueden utilizar las funciones print() y quantile() para extraer esta información.

print(curva, print.rmean = T)

## Call: survfit(formula = modelo, newdata = datosnuevos)
## 
##     n events rmean* se(rmean) median 0.95LCL 0.95UCL
## 1 181    134    436      31.8    361     291     558
## 2 181    134    389      25.2    329     239     574
##     * restricted mean with upper limit =  1022

quantile(curva)

## $quantile
##    25  50  75
## 1 199 361 643
## 2 175 329 567
## 
## $lower
##    25  50  75
## 1 156 291 477
## 2 118 239 390
## 
## $upper
##    25  50 75
## 1 288 558 NA
## 2 301 574 NA

6.1 Graficar la Curva de Sobrevivencia.

La visualización de las curvas de sobrevivencia en R puede realizarse eficazmente utilizando la función ggsurvplot de la librería survminer. Esta función es una interfaz avanzada que utiliza ggplot2 para la creación de gráficos de sobrevivencia, ofreciendo una amplia gama de opciones para personalizar la visualización. A continuación, se describen algunos de los parámetros más importantes:

Parámetros de ggsurvplot:

fit: Objeto de tipo survfit. Este objeto contiene las estimaciones de la función de sobrevivencia que se desean graficar.
data: Conjunto de datos utilizado para ajustar las curvas de sobrevivencia. Es opcional pero útil para etiquetar eventos o censuras directamente desde los datos.
fun: Transformación aplicada a la curva de sobrevivencia. Las opciones incluyen:
- "event": para mostrar los eventos acumulados.
- "cumhaz": para el riesgo acumulado.
- "pct": para representar la curva de sobrevivencia en porcentaje.
conf.int: Un indicador booleano que, si se establece en TRUE, permite la visualización de los intervalos de confianza alrededor de las estimaciones de sobrevivencia.
title: Título del gráfico.
xlab y ylab: Etiquetas para los ejes X e Y, respectivamente.
legend.labs: Vector de nombres para identificar las curvas, útil cuando se comparan varias curvas de sobrevivencia.
legend.title: Título para la leyenda del gráfico.

(modelo <- coxph(Surv(time, status) ~ age + meal.cal, data = datos))

## Call:
## coxph(formula = Surv(time, status) ~ age + meal.cal, data = datos)
## 
##                coef  exp(coef)   se(coef)      z      p
## age       1.892e-02  1.019e+00  1.044e-02  1.813 0.0699
## meal.cal -3.217e-05  1.000e+00  2.349e-04 -0.137 0.8911
## 
## Likelihood ratio test=3.67  on 2 df, p=0.1592
## n= 181, number of events= 134 
##    (47 observations deleted due to missingness)

datosnuevos <- data.frame(age = c(50, 60), meal.cal = c(1500, 2000)) #Datos a Estimar

curva <- survfit(modelo, newdata = datosnuevos, conf.type = "log-log")

ggsurvplot(fit = curva, data = datos, conf.int = T, title = "Curva de sobrevivencia",
           xlab = "time", ylab = "Probabilidad de sobrevivencia",
           legend.labs = c("age=50 y meal.cal=1500",
"age=60 y meal.cal=2000"), legend.title = "Grupos: ")

6.2 Estimación de Riesgo Base.

La estimación de riesgo base representa la particular predicción de la curva de sobrevivencia utilizando las covariables igual a 0.

(modelo <- coxph(Surv(time, status) ~ age + meal.cal, data = datos))

datosnuevo <- data.frame(age = 0, meal.cal = 0)
RiesgoBase <- survfit(modelo, datosnuevo)
summary(RiesgoBase)

cbind(Tiempo = RiesgoBase$time, RiesgoAcumulado = RiesgoBase$cumhaz)

Análisis de Sobrevivencia

Yolanda Gómez & Camila Baeza

2025-01-23

Capítulo 1 ¿Qué es el análisis de sobrevivencia?

1.1 Tiempos de fallas.

1.1.1 Aplicaciones.

1.1.2 Conceptos.

1.1.3 Datos completos.

1.2 Censura.

1.3 Tipos de censura.

1.3.1 Censura a la derecha.

1.4 Truncación.

1.5 Tiempo de Sobrevivencia.

1.5.1 Función de Sobrevivencia.

1.5.2 Función de riesgo.

1.5.3 Función de riesgo acumulada.

1.5.4 Tiempo medio y Vida media residual.

1.5.5 Algunas Relaciones Importantes.

Capítulo 2 Instalación de paquetes.

2.1 Bases de datos

2.2 Preparación de los datos.

2.3 El Objeto Surv.

Capítulo 3 Estimación No - Paramétrica de la Función de Sobrevivencia.

3.1 Estimador de Kaplan-Meier.

3.2 Desviación estándar e Intervalos de confianza de la estimación de la función de sobrevivencia.

3.3 Graficos de la curva de sobrevivencia.

3.4 Estimación de la función de riesgo acumulado.

3.5 Estimación Máxima Verosimilitud.

3.5.1 Ejemplo Sin Covariables.

3.5.2 Ejemplo Con Covariables

3.6 Estimación Nelson-Aalen.

3.6.1 Ejemplo sin covariables.

3.6.2. Ejemplo con covariables.

3.7 Estimación de la media, mediana y percentiles de los tiempos de Sobrevivencia.

3.8 Prueba de hipótesis para igualdad de dos o más funciones de sobrevivencia.

3.9 Comparaciones Múltiples por parejas.

3.9.1 Prueba Estratificada.

3.9.2 Prueba de Tendencia.

Capítulo 4 Estimación Paramétrica de la Función de Sobrevivencia.

4.1.1 Ejemplo.

4.2 Graficación del ajuste.

4.2.1 R Base.

4.2.2 Ggplot2.

4.2.3 Curva de Sobrevivencia.

4.2.4 Riesgo Acumulado.

4.2.4 Comparación de Curvas de Sobrevivencia por Grupos.

4.3 Gráficas Curvas No-Paramétricas con las Curvas Paramétricas.

4.3.1 Riesgo Acumulado.

4.3.2 ggplot2.

4.3.3 Curva de Sobrevivencia.

4.3.4 Riesgo Acumulado

4.4 Comparación de Modelos.

4.4.1 Comparación grafica.

4.4.2 R Base

4.4.3 Índice AIC, BIC y Log-verosimilitud.

4.4.4 Pruebas Gráficas para Modelos Específicos.

4.4.5 Gráfico de Bondad de Ajuste.

4.4.6 Gráfico de Probabilidad Estabilizado.

Capítulo 5 Modelo de Riesgos Proporcionales de Cox.

5.1 Ajustando el modelo de Cox

5.2 Método de Máxima Verosimilitud Parcial.

5.3 La Función de Verosimilitud.

5.4 Interpretación de los Coeficientes.

5.5 Estimando Funciones Relacionadas a h0(t).

Función de Tasa de Falla Acumulada

5.6 Elección del Modelo de Cox.

5.7 Verificando la suposición de riesgos proporcionales.

5.7.1 Métodos Gráficos.

5.8 Ejemplo en R

5.8.1 Preprocesamiento de los datos.

Variables Numéricas.

5.8.3 Variables Nominales u Ordinales.

5.8.4 Interacción.

5.8.5 Intervalos de Confianza y Pruebas de Significancia.

5.8.6 Comparación entre Modelos.

5.9 Validación de Supuestos.

5.9.1 Riesgos Proporcionales.

5.9.2 Análisis Residual.

6 Predicción de Curva de Sobrevivencia.

6.2 Estimación de Riesgo Base.