Análisis de Supervivencia en R.
Javier Martinez
22 de mayo de 2017
1 Introducción.
Los métodos estadísticos utilizado en el análisis de supervivencia frecuentemente presentan cierta complejidad en los cálculos, esto debido a la naturaleza de las observaciones las cuales suelen presentar censura y/o truncamiento. La utilización de software estadístico para realizar dicho análisis se ha vuelto indispensable en la práctica.
El presente manual tiene como finalidad dar un acercamiento el uso del software R en técnicas comúnmente usadas en un análisis de supervivencia.
El manual está dirigido para personas que han tomado un curso de Análisis de Supervivencia y tengan conocimientos básicos con el software R.
El manual está dividido por temas en los cuales se presentarán los comandos, argumentos y su forma de emplearlos para ejecutar sus correspondientes técnicas.
La sintaxis utilizada en el manual es el siguiente:
Los script’s, entradas, salidas y comentarios están escritas en forma de bloques de la siguiente manera:
# Esto es un Comentario
x <- 2 * 4 #esto es un Entrada
x #Salida## [1] 8Las salidas del programa son diferenciadas por el símbolo “#” al inicio de la línea, mientras que los comentarios se mostraran de un color naranja oscuro. Todos los comandos o funciones utilizadas en el texto se verán de esta manera.
1.1 Instalación de paquetes.
El software R, así como sus librerías están en constante actualización por lo que se recomienda revisar si las librerías utilizadas soportan la versión actualmente usada de R. Este manual fue realizado utilizando la versión 3.3.2 (2016-10-31) “Sincere Pumpkin Patch” de R.
Las librerías utilizadas, así como sus versiones fueron las en este manual son las siguientes:
## [1] "survival (2.41.3)"
## [1] "KMsurv (0.1.5)"
## [1] "survMisc (0.5.4)"
## [1] "survminer (0.4.0)"
## [1] "flexsurv (1.1)"
## [1] "actuar (2.1.1)"Las librerías se instalan por una única vez y se cargan antes de iniciar cualquier análisis.
La instalación se realiza con el siguiente script
install.packages("survival")
install.packages("KMsurv")
install.packages("survMisc")
install.packages("survminer")
install.packages("flexsurv")
install.packages("actuar")
install.packages("dplyr")Mientras que la carga de las librerías se realiza con la función library()
library(survival)
library(KMsurv)
library(survMisc)
library(survminer)
library(ggfortify)
library(flexsurv)
library(actuar)
library(dplyr)1.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/Estadistica/R"). Es necesario proporcionar la dirección completa del archivo si este no se encuentra en el directorio de trabajo.
Antes de iniciar con los análisis comando 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 formado 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.
Así un estudio con 40 pacientes y 8 variables medidas corresponderá a una tabla de 40 + 1 = 41 Filas y 8 columnas.
Es importante que el nombre de cada variable en el encabezado sea corto y sencillo, ya que con ese mismo nombre se hará referencia en los análisis.
La lectura de los datos se realiza con la función read.csv(file,header) utilizando los siguientes argumentos:
- file: El nombre del archivo del que se deben leer los datos.
- header: Valor lógico (TRUE/FALSE) que indica si el archivo contiene el nombre de las variables en la primera fila.
Ejemplo:
datos <- read.csv(file = "datos/cox.csv", header = T) #Leer Datos
head(datos)## Tiempo Censura Sexo Edad Ciudad IMC
## 1 12.04 1 0 42 1 19
## 2 0.72 1 1 29 4 20
## 3 44.65 1 0 37 3 23
## 4 56.80 0 1 50 4 26
## 5 29.29 1 1 38 1 17
## 6 4.51 1 0 48 3 251.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. Bases de datos complementaría se podrá encontrar en http://www.dropbox.com/sh/40lwu3hcvlhfcwa/AAARdsePFmpOoOnmOLL1yPl6a
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 tongue.
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 11.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 este manual 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.
Ejemplo:
Surv(5, 1) #Observación del evento## [1] 5Surv(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. El objeto, se puede guardar en una variable para utilizar posteriormente, sin embargo, en este manual no se hará.
2 Estimación no-paramétrica de la función de supervivencia
2.1 Estimador de Kaplan-Meier y Fleming-Harrington.
Los estimadores de Kaplan-Meier y Fleming-Harrington para la función de supervivencia 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 supervivencia 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
Survcomo 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 estrato. Para una sola curva de supervivencia del lado derecho se coloca~1. - data. objeto data frame donde están los datos.
- type. Tipo de estimador: “kaplan-meier” o “fleming-harrington”.
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).
# tongue<-read.csv(file = 'datos/tongue.csv',header = TRUE)
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) #Creando objeto tipo Surv
tongue.km <- survfit(tongue.surv ~ 1, data = tongue, type = "kaplan-meier") #Estimación Kaplan Meier
# 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 supervivencia se lleva a cabo con la función summary().
summary(tongue.km)## Call: survfit(formula = Surv(time, delta) ~ 1, data = tongue, type = "kaplan-meier")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 80 2 0.975 0.0175 0.9414 1.000
## 3 78 3 0.938 0.0271 0.8859 0.992
## 4 75 2 0.913 0.0316 0.8526 0.977
## 5 73 2 0.888 0.0353 0.8209 0.960
## 8 71 1 0.875 0.0370 0.8054 0.951
## 10 69 1 0.862 0.0386 0.7900 0.941
## 12 68 1 0.850 0.0400 0.7747 0.932
## 13 67 3 0.812 0.0438 0.7301 0.902
## 16 64 2 0.786 0.0460 0.7011 0.882
## 18 62 1 0.774 0.0470 0.6868 0.871
## 23 61 1 0.761 0.0479 0.6726 0.861
## 24 60 1 0.748 0.0487 0.6585 0.850
## 26 59 2 0.723 0.0503 0.6307 0.828
## 27 57 2 0.697 0.0516 0.6033 0.806
## 28 55 1 0.685 0.0522 0.5897 0.795
## 30 54 3 0.647 0.0537 0.5496 0.761
## 32 51 1 0.634 0.0541 0.5363 0.750
## 41 50 1 0.621 0.0545 0.5232 0.738
## 42 49 1 0.609 0.0549 0.5101 0.726
## 51 48 1 0.596 0.0552 0.4971 0.715
## 56 47 1 0.583 0.0554 0.4842 0.703
## 62 45 1 0.570 0.0557 0.4710 0.691
## 65 44 1 0.557 0.0559 0.4579 0.679
## 67 43 1 0.544 0.0561 0.4449 0.666
## 69 41 1 0.531 0.0563 0.4315 0.654
## 70 40 1 0.518 0.0564 0.4183 0.641
## 72 39 1 0.505 0.0565 0.4051 0.628
## 73 38 1 0.491 0.0566 0.3921 0.616
## 77 35 1 0.477 0.0567 0.3782 0.602
## 91 27 1 0.460 0.0573 0.3600 0.587
## 93 26 1 0.442 0.0577 0.3421 0.571
## 96 24 1 0.424 0.0582 0.3236 0.554
## 100 22 1 0.404 0.0586 0.3042 0.537
## 104 20 3 0.344 0.0594 0.2449 0.482
## 112 13 1 0.317 0.0604 0.2184 0.461
## 129 11 1 0.288 0.0614 0.1900 0.438
## 157 8 1 0.252 0.0634 0.1541 0.413
## 167 7 1 0.216 0.0638 0.1213 0.386
## 181 5 1 0.173 0.0640 0.0838 0.357La estimación devuelve los siguientes valores:
- time : Tiempo de la observación
- n.risk : El número de sujetos en riesgo.
- n.evento : El número de sujetos que presentaron el evento.
- survival : La estimación de la función de supervivencia.
- std.err : La desviación estándar de la estimación.
- 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 Kaplan Meier por grupos
summary(tongue.km) #Resumen estadístico ## Call: survfit(formula = Surv(time, delta) ~ type, data = tongue, type = "kaplan-meier")
##
## type=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 52 1 0.981 0.0190 0.944 1.000
## 3 51 2 0.942 0.0323 0.881 1.000
## 4 49 1 0.923 0.0370 0.853 0.998
## 10 48 1 0.904 0.0409 0.827 0.988
## 13 47 2 0.865 0.0473 0.777 0.963
## 16 45 2 0.827 0.0525 0.730 0.936
## 24 43 1 0.808 0.0547 0.707 0.922
## 26 42 1 0.788 0.0566 0.685 0.908
## 27 41 1 0.769 0.0584 0.663 0.893
## 28 40 1 0.750 0.0600 0.641 0.877
## 30 39 2 0.712 0.0628 0.598 0.846
## 32 37 1 0.692 0.0640 0.578 0.830
## 41 36 1 0.673 0.0651 0.557 0.813
## 51 35 1 0.654 0.0660 0.537 0.797
## 65 33 1 0.634 0.0669 0.516 0.780
## 67 32 1 0.614 0.0677 0.495 0.762
## 70 31 1 0.594 0.0683 0.475 0.745
## 72 30 1 0.575 0.0689 0.454 0.727
## 73 29 1 0.555 0.0693 0.434 0.709
## 77 27 1 0.534 0.0697 0.414 0.690
## 91 19 1 0.506 0.0715 0.384 0.667
## 93 18 1 0.478 0.0728 0.355 0.644
## 96 16 1 0.448 0.0741 0.324 0.620
## 100 14 1 0.416 0.0754 0.292 0.594
## 104 12 1 0.381 0.0767 0.257 0.566
## 157 5 1 0.305 0.0918 0.169 0.550
## 167 4 1 0.229 0.0954 0.101 0.518
##
## type=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 28 1 0.9643 0.0351 0.8979 1.000
## 3 27 1 0.9286 0.0487 0.8379 1.000
## 4 26 1 0.8929 0.0585 0.7853 1.000
## 5 25 2 0.8214 0.0724 0.6911 0.976
## 8 23 1 0.7857 0.0775 0.6475 0.953
## 12 21 1 0.7483 0.0824 0.6031 0.929
## 13 20 1 0.7109 0.0863 0.5603 0.902
## 18 19 1 0.6735 0.0895 0.5190 0.874
## 23 18 1 0.6361 0.0921 0.4790 0.845
## 26 17 1 0.5986 0.0939 0.4402 0.814
## 27 16 1 0.5612 0.0952 0.4025 0.783
## 30 15 1 0.5238 0.0959 0.3658 0.750
## 42 14 1 0.4864 0.0961 0.3302 0.716
## 56 13 1 0.4490 0.0957 0.2956 0.682
## 62 12 1 0.4116 0.0948 0.2621 0.646
## 69 10 1 0.3704 0.0938 0.2255 0.608
## 104 8 2 0.2778 0.0904 0.1468 0.526
## 112 5 1 0.2222 0.0877 0.1025 0.482
## 129 4 1 0.1667 0.0815 0.0639 0.435
## 181 2 1 0.0833 0.0717 0.0155 0.449tongue.km[2]$surv #Elementos de la columna surv de delta=1## [1] 0.964 0.929 0.893 0.821 0.786 0.748 0.711 0.673 0.636 0.599 0.561
## [12] 0.524 0.486 0.449 0.412 0.412 0.370 0.370 0.278 0.222 0.167 0.167
## [23] 0.083 0.083fortify(tongue.km)## time n.risk n.event n.censor surv std.err upper lower strata
## 1 1 52 1 0 0.981 0.019 1.00 0.944 1
## 2 3 51 2 0 0.942 0.034 1.00 0.881 1
## 3 4 49 1 0 0.923 0.040 1.00 0.853 1
## 4 10 48 1 0 0.904 0.045 0.99 0.827 1
## 5 13 47 2 0 0.865 0.055 0.96 0.777 1
## 6 16 45 2 0 0.827 0.063 0.94 0.730 1
## 7 24 43 1 0 0.808 0.068 0.92 0.707 1
## 8 26 42 1 0 0.788 0.072 0.91 0.685 1
## 9 27 41 1 0 0.769 0.076 0.89 0.663 1
## 10 28 40 1 0 0.750 0.080 0.88 0.641 1
## 11 30 39 2 0 0.712 0.088 0.85 0.598 1
## 12 32 37 1 0 0.692 0.092 0.83 0.578 1
## 13 41 36 1 0 0.673 0.097 0.81 0.557 1
## 14 51 35 1 0 0.654 0.101 0.80 0.537 1
## 15 61 34 0 1 0.654 0.101 0.80 0.537 1
## 16 65 33 1 0 0.634 0.105 0.78 0.516 1
## 17 67 32 1 0 0.614 0.110 0.76 0.495 1
## 18 70 31 1 0 0.594 0.115 0.74 0.475 1
## 19 72 30 1 0 0.575 0.120 0.73 0.454 1
## 20 73 29 1 0 0.555 0.125 0.71 0.434 1
## 21 74 28 0 1 0.555 0.125 0.71 0.434 1
## 22 77 27 1 0 0.534 0.130 0.69 0.414 1
## 23 79 26 0 1 0.534 0.130 0.69 0.414 1
## 24 80 25 0 1 0.534 0.130 0.69 0.414 1
## 25 81 24 0 1 0.534 0.130 0.69 0.414 1
## 26 87 23 0 2 0.534 0.130 0.69 0.414 1
## 27 88 21 0 1 0.534 0.130 0.69 0.414 1
## 28 89 20 0 1 0.534 0.130 0.69 0.414 1
## 29 91 19 1 0 0.506 0.141 0.67 0.384 1
## 30 93 18 1 1 0.478 0.152 0.64 0.355 1
## 31 96 16 1 0 0.448 0.165 0.62 0.324 1
## 32 97 15 0 1 0.448 0.165 0.62 0.324 1
## 33 100 14 1 0 0.416 0.181 0.59 0.292 1
## 34 101 13 0 1 0.416 0.181 0.59 0.292 1
## 35 104 12 1 1 0.381 0.201 0.57 0.257 1
## 36 108 10 0 1 0.381 0.201 0.57 0.257 1
## 37 109 9 0 1 0.381 0.201 0.57 0.257 1
## 38 120 8 0 1 0.381 0.201 0.57 0.257 1
## 39 131 7 0 1 0.381 0.201 0.57 0.257 1
## 40 150 6 0 1 0.381 0.201 0.57 0.257 1
## 41 157 5 1 0 0.305 0.301 0.55 0.169 1
## 42 167 4 1 0 0.229 0.417 0.52 0.101 1
## 43 231 3 0 1 0.229 0.417 0.52 0.101 1
## 44 240 2 0 1 0.229 0.417 0.52 0.101 1
## 45 400 1 0 1 0.229 0.417 0.52 0.101 1
## 46 1 28 1 0 0.964 0.036 1.00 0.898 2
## 47 3 27 1 0 0.929 0.052 1.00 0.838 2
## 48 4 26 1 0 0.893 0.065 1.00 0.785 2
## 49 5 25 2 0 0.821 0.088 0.98 0.691 2
## 50 8 23 1 1 0.786 0.099 0.95 0.648 2
## 51 12 21 1 0 0.748 0.110 0.93 0.603 2
## 52 13 20 1 0 0.711 0.121 0.90 0.560 2
## 53 18 19 1 0 0.673 0.133 0.87 0.519 2
## 54 23 18 1 0 0.636 0.145 0.84 0.479 2
## 55 26 17 1 0 0.599 0.157 0.81 0.440 2
## 56 27 16 1 0 0.561 0.170 0.78 0.402 2
## 57 30 15 1 0 0.524 0.183 0.75 0.366 2
## 58 42 14 1 0 0.486 0.198 0.72 0.330 2
## 59 56 13 1 0 0.449 0.213 0.68 0.296 2
## 60 62 12 1 0 0.412 0.230 0.65 0.262 2
## 61 67 11 0 1 0.412 0.230 0.65 0.262 2
## 62 69 10 1 0 0.370 0.253 0.61 0.225 2
## 63 76 9 0 1 0.370 0.253 0.61 0.225 2
## 64 104 8 2 1 0.278 0.325 0.53 0.147 2
## 65 112 5 1 0 0.222 0.395 0.48 0.103 2
## 66 129 4 1 0 0.167 0.489 0.43 0.064 2
## 67 176 3 0 1 0.167 0.489 0.43 0.064 2
## 68 181 2 1 0 0.083 0.860 0.45 0.015 2
## 69 231 1 0 1 0.083 0.860 0.45 0.015 22.2 Desviación estándar e Intervalos de confianza de la estimación de la función de supervivencia.
Tanto como la desviación estándar y los Intervalos de confianza de la curva de supervivencia 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)## Call: survfit(formula = Surv(time, delta) ~ 1, data = tongue, type = "kaplan-meier",
## error = "tsiatis", conf.type = "log-log", conf.int = 0.99)
##
## time n.risk n.event survival std.err lower 99% CI upper 99% CI
## 1 80 2 0.975 0.0172 0.8582 0.996
## 3 78 3 0.938 0.0266 0.8184 0.979
## 4 75 2 0.913 0.0311 0.7876 0.965
## 5 73 2 0.888 0.0348 0.7572 0.950
## 8 71 1 0.875 0.0364 0.7421 0.942
## 10 69 1 0.862 0.0380 0.7269 0.934
## 12 68 1 0.850 0.0395 0.7119 0.925
## 13 67 3 0.812 0.0432 0.6687 0.897
## 16 64 2 0.786 0.0453 0.6404 0.878
## 18 62 1 0.774 0.0463 0.6263 0.869
## 23 61 1 0.761 0.0472 0.6124 0.859
## 24 60 1 0.748 0.0480 0.5986 0.849
## 26 59 2 0.723 0.0495 0.5717 0.828
## 27 57 2 0.697 0.0508 0.5451 0.807
## 28 55 1 0.685 0.0515 0.5319 0.797
## 30 54 3 0.647 0.0528 0.4935 0.764
## 32 51 1 0.634 0.0533 0.4806 0.753
## 41 50 1 0.621 0.0537 0.4679 0.742
## 42 49 1 0.609 0.0540 0.4553 0.731
## 51 48 1 0.596 0.0543 0.4428 0.720
## 56 47 1 0.583 0.0546 0.4304 0.709
## 62 45 1 0.570 0.0549 0.4177 0.697
## 65 44 1 0.557 0.0551 0.4051 0.685
## 67 43 1 0.544 0.0553 0.3926 0.673
## 69 41 1 0.531 0.0555 0.3799 0.661
## 70 40 1 0.518 0.0556 0.3672 0.649
## 72 39 1 0.505 0.0557 0.3547 0.637
## 73 38 1 0.491 0.0558 0.3422 0.624
## 77 35 1 0.477 0.0559 0.3290 0.611
## 91 27 1 0.460 0.0564 0.3111 0.596
## 93 26 1 0.442 0.0569 0.2937 0.580
## 96 24 1 0.424 0.0573 0.2756 0.564
## 100 22 1 0.404 0.0577 0.2569 0.547
## 104 20 3 0.344 0.0573 0.2024 0.490
## 112 13 1 0.317 0.0583 0.1766 0.468
## 129 11 1 0.288 0.0591 0.1493 0.443
## 157 8 1 0.252 0.0606 0.1156 0.415
## 167 7 1 0.216 0.0604 0.0863 0.384
## 181 5 1 0.173 0.0594 0.0547 0.3472.3 Graficación de la curva de supervivencia.
La curva de supervivencia estimada se gráfica con la función ggsurvplot() de 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 y se recomienda revisar los demás utilizando el comando help("ggsurvplot").
- fit: Objeto tipo survfit.
- data: Un conjunto de datos utilizado para ajustar curvas de supervivencia.
- fun: Transformación de la curva de supervivencia(Opcional), las posibles opciones son: “event” para los eventos acumulados, “cumhaz” para el riesgo acumulado y “pct” para la curva de supervivencia en porcentaje.
- conf.int: Indicador para graficar los intervalos de confianza.
- title: Titulo
- xlab: Eje x
- ylab: Eje Y
- legends.lab: Vector de nombres para identificar las curvas.
- legend.title : Titulo de la leyenda.
2.3.1 Ejemplo
ggsurvplot(fit = tongue.km, data = tongue, conf.int = T, title = "Curva de Supervivencia",
xlab = "Tiempo", ylab = "Probabilidad de supervivencia", legend.title = "Estimación",
legend.labs = "Kaplan-Meier")
2.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 supervivencia y el segundo es utilizando el estimador Nelson-Aalen.
2.4.1 Estimación Máxima Verosimilitud.
La información del Tiempo y Probabilidad de supervivencia es extraída del objeto survfit (tongue.km) utilizando fortify, posteriormente si existe covariable se agrupa utilizando group_by, 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, se puede encontrar más información en http://dplyr.tidyverse.org/.
2.4.1.1 Ejemplo Sin Covariables
tongue.km <- survfit(Surv(time, delta) ~ 1, tongue)
R <- tongue.km %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))## Warning: package 'bindrcpp' was built under R version 3.3.3R## time n.risk n.event n.censor surv std.err upper lower CumHaz
## 1 1 80 2 0 0.97 0.018 1.00 0.941 0.025
## 2 3 78 3 0 0.94 0.029 0.99 0.886 0.063
## 3 4 75 2 0 0.91 0.035 0.98 0.853 0.090
## 4 5 73 2 0 0.89 0.040 0.96 0.821 0.118
## 5 8 71 1 1 0.88 0.042 0.95 0.805 0.132
## 6 10 69 1 0 0.86 0.045 0.94 0.790 0.146
## 7 12 68 1 0 0.85 0.047 0.93 0.775 0.161
## 8 13 67 3 0 0.81 0.054 0.90 0.730 0.206
## 9 16 64 2 0 0.79 0.058 0.88 0.701 0.237
## 10 18 62 1 0 0.77 0.061 0.87 0.687 0.253
## 11 23 61 1 0 0.76 0.063 0.86 0.673 0.269
## 12 24 60 1 0 0.75 0.065 0.85 0.659 0.286
## 13 26 59 2 0 0.72 0.070 0.83 0.631 0.320
## 14 27 57 2 0 0.70 0.074 0.81 0.603 0.355
## 15 28 55 1 0 0.68 0.076 0.80 0.590 0.373
## 16 30 54 3 0 0.65 0.083 0.76 0.550 0.429
## 17 32 51 1 0 0.63 0.085 0.75 0.536 0.448
## 18 41 50 1 0 0.62 0.088 0.74 0.523 0.468
## 19 42 49 1 0 0.61 0.090 0.73 0.510 0.489
## 20 51 48 1 0 0.60 0.093 0.71 0.497 0.510
## 21 56 47 1 0 0.58 0.095 0.70 0.484 0.531
## 22 61 46 0 1 0.58 0.095 0.70 0.484 0.531
## 23 62 45 1 0 0.57 0.098 0.69 0.471 0.553
## 24 65 44 1 0 0.56 0.100 0.68 0.458 0.576
## 25 67 43 1 1 0.54 0.103 0.67 0.445 0.599
## 26 69 41 1 0 0.53 0.106 0.65 0.432 0.623
## 27 70 40 1 0 0.52 0.109 0.64 0.418 0.648
## 28 72 39 1 0 0.50 0.112 0.63 0.405 0.674
## 29 73 38 1 0 0.49 0.115 0.62 0.392 0.700
## 30 74 37 0 1 0.49 0.115 0.62 0.392 0.700
## 31 76 36 0 1 0.49 0.115 0.62 0.392 0.700
## 32 77 35 1 0 0.48 0.119 0.60 0.378 0.729
## 33 79 34 0 1 0.48 0.119 0.60 0.378 0.729
## 34 80 33 0 1 0.48 0.119 0.60 0.378 0.729
## 35 81 32 0 1 0.48 0.119 0.60 0.378 0.729
## 36 87 31 0 2 0.48 0.119 0.60 0.378 0.729
## 37 88 29 0 1 0.48 0.119 0.60 0.378 0.729
## 38 89 28 0 1 0.48 0.119 0.60 0.378 0.729
## 39 91 27 1 0 0.46 0.125 0.59 0.360 0.766
## 40 93 26 1 1 0.44 0.131 0.57 0.342 0.804
## 41 96 24 1 0 0.42 0.137 0.55 0.324 0.846
## 42 97 23 0 1 0.42 0.137 0.55 0.324 0.846
## 43 100 22 1 0 0.40 0.145 0.54 0.304 0.892
## 44 101 21 0 1 0.40 0.145 0.54 0.304 0.892
## 45 104 20 3 2 0.34 0.173 0.48 0.245 1.042
## 46 108 15 0 1 0.34 0.173 0.48 0.245 1.042
## 47 109 14 0 1 0.34 0.173 0.48 0.245 1.042
## 48 112 13 1 0 0.32 0.190 0.46 0.218 1.119
## 49 120 12 0 1 0.32 0.190 0.46 0.218 1.119
## 50 129 11 1 0 0.29 0.213 0.44 0.190 1.209
## 51 131 10 0 1 0.29 0.213 0.44 0.190 1.209
## 52 150 9 0 1 0.29 0.213 0.44 0.190 1.209
## 53 157 8 1 0 0.25 0.251 0.41 0.154 1.334
## 54 167 7 1 0 0.22 0.295 0.39 0.121 1.477
## 55 176 6 0 1 0.22 0.295 0.39 0.121 1.477
## 56 181 5 1 0 0.17 0.370 0.36 0.084 1.677
## 57 231 4 0 2 0.17 0.370 0.36 0.084 1.677
## 58 240 2 0 1 0.17 0.370 0.36 0.084 1.677
## 59 400 1 0 1 0.17 0.370 0.36 0.084 1.6772.4.1.2 Ejemplo Con Covariables
tongue.km <- survfit(Surv(time, delta) ~ type, tongue)
R <- tongue.km %>% fortify %>% group_by(strata) %>% mutate(CumHaz = cumsum(n.event/n.risk))
print(R, 60)## # A tibble: 69 x 10
## # Groups: strata [2]
## time n.risk n.event n.censor surv std.err upper lower strata CumHaz
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fctr> <dbl>
## 1 1 52 1 0 0.98 0.019 1.00 0.94 1 0.019
## 2 3 51 2 0 0.94 0.034 1.00 0.88 1 0.058
## 3 4 49 1 0 0.92 0.040 1.00 0.85 1 0.079
## 4 10 48 1 0 0.90 0.045 0.99 0.83 1 0.100
## 5 13 47 2 0 0.87 0.055 0.96 0.78 1 0.142
## 6 16 45 2 0 0.83 0.063 0.94 0.73 1 0.187
## 7 24 43 1 0 0.81 0.068 0.92 0.71 1 0.210
## 8 26 42 1 0 0.79 0.072 0.91 0.68 1 0.234
## 9 27 41 1 0 0.77 0.076 0.89 0.66 1 0.258
## 10 28 40 1 0 0.75 0.080 0.88 0.64 1 0.283
## # ... with 59 more rowsLa 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 (Semanas)", ylab = "Riesgo Acumulado")
ggsurvplot(tongue.km, fun = "cumhaz", xlab = "Tiempo (Semanas)", censor = T,
ylab = "Riesgo Acumulado", title = "Riesgo Acumulado", legend.title = "Perfil de ADN tumoral",
legend.labs = c("Aneuploide", "Diploide"))
2.4.2 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.
2.4.2.1 Ejemplo sin covariables
tongue.km <- survfit(Surv(time, delta) ~ 1, tongue)
R <- tongue.km %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
R## time n.risk n.event n.censor surv std.err upper lower CumHaz
## 1 1 80 2 0 0.97 0.018 1.00 0.941 0.025
## 2 3 78 3 0 0.94 0.029 0.99 0.886 0.063
## 3 4 75 2 0 0.91 0.035 0.98 0.853 0.090
## 4 5 73 2 0 0.89 0.040 0.96 0.821 0.118
## 5 8 71 1 1 0.88 0.042 0.95 0.805 0.132
## 6 10 69 1 0 0.86 0.045 0.94 0.790 0.146
## 7 12 68 1 0 0.85 0.047 0.93 0.775 0.161
## 8 13 67 3 0 0.81 0.054 0.90 0.730 0.206
## 9 16 64 2 0 0.79 0.058 0.88 0.701 0.237
## 10 18 62 1 0 0.77 0.061 0.87 0.687 0.253
## 11 23 61 1 0 0.76 0.063 0.86 0.673 0.269
## 12 24 60 1 0 0.75 0.065 0.85 0.659 0.286
## 13 26 59 2 0 0.72 0.070 0.83 0.631 0.320
## 14 27 57 2 0 0.70 0.074 0.81 0.603 0.355
## 15 28 55 1 0 0.68 0.076 0.80 0.590 0.373
## 16 30 54 3 0 0.65 0.083 0.76 0.550 0.429
## 17 32 51 1 0 0.63 0.085 0.75 0.536 0.448
## 18 41 50 1 0 0.62 0.088 0.74 0.523 0.468
## 19 42 49 1 0 0.61 0.090 0.73 0.510 0.489
## 20 51 48 1 0 0.60 0.093 0.71 0.497 0.510
## 21 56 47 1 0 0.58 0.095 0.70 0.484 0.531
## 22 61 46 0 1 0.58 0.095 0.70 0.484 0.531
## 23 62 45 1 0 0.57 0.098 0.69 0.471 0.553
## 24 65 44 1 0 0.56 0.100 0.68 0.458 0.576
## 25 67 43 1 1 0.54 0.103 0.67 0.445 0.599
## 26 69 41 1 0 0.53 0.106 0.65 0.432 0.623
## 27 70 40 1 0 0.52 0.109 0.64 0.418 0.648
## 28 72 39 1 0 0.50 0.112 0.63 0.405 0.674
## 29 73 38 1 0 0.49 0.115 0.62 0.392 0.700
## 30 74 37 0 1 0.49 0.115 0.62 0.392 0.700
## 31 76 36 0 1 0.49 0.115 0.62 0.392 0.700
## 32 77 35 1 0 0.48 0.119 0.60 0.378 0.729
## 33 79 34 0 1 0.48 0.119 0.60 0.378 0.729
## 34 80 33 0 1 0.48 0.119 0.60 0.378 0.729
## 35 81 32 0 1 0.48 0.119 0.60 0.378 0.729
## 36 87 31 0 2 0.48 0.119 0.60 0.378 0.729
## 37 88 29 0 1 0.48 0.119 0.60 0.378 0.729
## 38 89 28 0 1 0.48 0.119 0.60 0.378 0.729
## 39 91 27 1 0 0.46 0.125 0.59 0.360 0.766
## 40 93 26 1 1 0.44 0.131 0.57 0.342 0.804
## 41 96 24 1 0 0.42 0.137 0.55 0.324 0.846
## 42 97 23 0 1 0.42 0.137 0.55 0.324 0.846
## 43 100 22 1 0 0.40 0.145 0.54 0.304 0.892
## 44 101 21 0 1 0.40 0.145 0.54 0.304 0.892
## 45 104 20 3 2 0.34 0.173 0.48 0.245 1.042
## 46 108 15 0 1 0.34 0.173 0.48 0.245 1.042
## 47 109 14 0 1 0.34 0.173 0.48 0.245 1.042
## 48 112 13 1 0 0.32 0.190 0.46 0.218 1.119
## 49 120 12 0 1 0.32 0.190 0.46 0.218 1.119
## 50 129 11 1 0 0.29 0.213 0.44 0.190 1.209
## 51 131 10 0 1 0.29 0.213 0.44 0.190 1.209
## 52 150 9 0 1 0.29 0.213 0.44 0.190 1.209
## 53 157 8 1 0 0.25 0.251 0.41 0.154 1.334
## 54 167 7 1 0 0.22 0.295 0.39 0.121 1.477
## 55 176 6 0 1 0.22 0.295 0.39 0.121 1.477
## 56 181 5 1 0 0.17 0.370 0.36 0.084 1.677
## 57 231 4 0 2 0.17 0.370 0.36 0.084 1.677
## 58 240 2 0 1 0.17 0.370 0.36 0.084 1.677
## 59 400 1 0 1 0.17 0.370 0.36 0.084 1.677qplot(time, CumHaz, data = R, geom = "step", xlab = "Tiempo (Semanas)", ylab = "Riesgo Acumulado",
main = "Riesgo Acumulado")
2.4.2.2 Ejemplo con covariables
tongue.km <- survfit(Surv(time, delta) ~ type, tongue)
R <- tongue.km %>% fortify %>% group_by(strata) %>% mutate(CumHaz = cumsum(n.event/n.risk))
R## # A tibble: 69 x 10
## # Groups: strata [2]
## time n.risk n.event n.censor surv std.err upper lower strata CumHaz
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fctr> <dbl>
## 1 1 52 1 0 0.98 0.019 1.00 0.94 1 0.019
## 2 3 51 2 0 0.94 0.034 1.00 0.88 1 0.058
## 3 4 49 1 0 0.92 0.040 1.00 0.85 1 0.079
## 4 10 48 1 0 0.90 0.045 0.99 0.83 1 0.100
## 5 13 47 2 0 0.87 0.055 0.96 0.78 1 0.142
## 6 16 45 2 0 0.83 0.063 0.94 0.73 1 0.187
## 7 24 43 1 0 0.81 0.068 0.92 0.71 1 0.210
## 8 26 42 1 0 0.79 0.072 0.91 0.68 1 0.234
## 9 27 41 1 0 0.77 0.076 0.89 0.66 1 0.258
## 10 28 40 1 0 0.75 0.080 0.88 0.64 1 0.283
## # ... with 59 more rowsqplot(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, y coincide con el método de Nelson-Aalen si se estimó la función de supervivencia vía Fleming-Harrington.
2.5 Estimación de la media, mediana y percentiles de los tiempos de supervivencia.
La estimación de la media se realiza integrando la función supervivencia estimada hasta un valor \(\tau\).
\[ \hat\mu = \int_{0}^{\tau}\hat{S}(t)dt \]
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 \(\tau\) (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 127.3 20.1 93 67 NA
## type=2 28 22 79.7 18.9 42 23 112
## * restricted mean with upper limit = 316Tanto 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 NA2.6 Prueba de hipótesis para igualdad de dos o más funciones de supervivencia.
Dado dos o más grupos se tener la siguiente la siguiente prueba de hipótesis:
\[H_0 : h_1(t)=h_2(t)=...=h_n(t)\text{ para todo }t \leq \tau \] \[H_1: h_i(t_0) \neq h_j(t_0)\text{ para al menos un par i,j y tiempo } t_0 \]
O en su versión equivalente
\[H_0 : S_1(t)=S_2(t)=...=S_k(t)\text{ para todo }t \leq \tau \] \[H_1: S_i(t_0) \neq S_j(t_0)\text{ para almenos un par i,j y tiempo }t_0 \]
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 estratos. Para una sola curva de supervivencia 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 \(\hat{S}(t)^\rho\)
Estas pruebas están basadas en la familia G-rho de Harrington and Fleming, acuerdo a los valores particulares que se le asigna a rho se tiene 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:
La base de datos larynx contiene 90 observaciones de 5 variables:
- time. Tiempo hasta la muerte o en estudio en meses
- delta. Indicador(0=Vivo,1=Muerte)
- state. Estado del deceso(1,2,3 y 4)
- age. Edad al diagnóstico de cáncer de laringue.
- diagyr. Año del diagnóstico del cáncer
data("larynx")
head(larynx)## stage time age diagyr delta
## 1 1 0.6 77 76 1
## 2 1 1.3 53 71 1
## 3 1 2.4 45 71 1
## 4 1 2.5 57 78 0
## 5 1 3.2 58 74 1
## 6 1 3.2 51 77 0Utilizando esta información se compara si existe alguna diferencia de las curvas de supervivencia entre los estados del deceso.
survdiff(Surv(time, delta) ~ stage, data = larynx, rho = 0) #Prueba log-rank## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx,
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=1 33 15 22.57 2.537 4.741
## stage=2 17 7 10.01 0.906 1.152
## stage=3 27 17 14.08 0.603 0.856
## stage=4 13 11 3.34 17.590 19.827
##
## Chisq= 22.8 on 3 degrees of freedom, p= 4.53e-05survdiff(Surv(time, delta) ~ stage, data = larynx, rho = 1) #Prueba Peto & Peto## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx,
## rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=1 33 9.26 15.78 2.694 6.48
## stage=2 17 4.66 7.18 0.884 1.47
## stage=3 27 12.77 10.06 0.730 1.35
## stage=4 13 9.15 2.82 14.219 18.78
##
## Chisq= 23.1 on 3 degrees of freedom, p= 3.85e-052.6.1 Comparaciones Múltiples por parejas.
Las \((^n_2)\) combinaciones 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 (Practica no recomendada), use
p.adjust.method = "none". - rho. Un valor numérico que asigna los pesos de acuerdo a \(\hat{S}(t)^\rho\)
Ejemplo:
En la base de datos larunx, la covariable stage no es de tipo factor por lo que se procederá a convertirlo con factor() antes de utilizar la prueba.
class(larynx$stage) #Tipo Numérico## [1] "integer"larynx$stage <- factor(larynx$stage)
class(larynx$stage) #Tipo factor## [1] "factor"pairwise_survdiff(Surv(time, delta) ~ stage, data = larynx, p.adjust.method = "bonferroni",
rho = 1)##
## Pairwise comparisons using Peto & Peto test
##
## data: larynx and stage
##
## 1 2 3
## 2 1.000 - -
## 3 0.201 0.945 -
## 4 8e-06 0.007 0.265
##
## P value adjustment method: bonferroni2.6.2 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 larynx y crearemos una nueva variable llamada GrupoEdad en donde se agrupe a las personas en 3 tipos de edades: 40-55, 55-75, 75+ utilizando la función cut(), esta variable será utilizada como estrato.
larynx$GrupoEdad <- cut(larynx$age, breaks = c(40, 55, 75, 90)) #Creando nueva variable
survdiff(Surv(time, delta) ~ stage + strata(GrupoEdad), data = larynx, rho = 1)## Call:
## survdiff(formula = Surv(time, delta) ~ stage + strata(GrupoEdad),
## data = larynx, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=1 33 8.90 16.42 3.44 8.68
## stage=2 17 4.78 6.55 0.48 0.77
## stage=3 27 12.98 9.78 1.05 1.89
## stage=4 13 9.17 3.08 12.03 16.29
##
## Chisq= 21.7 on 3 degrees of freedom, p= 7.65e-052.6.3 Prueba de Tendencia.
La función comp() de la librería survMisc realiza comparaciones entre curvas de supervivencia 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(Ver Ejemplo). - 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
model <- ten(Surv(time, delta) ~ stage, data = larynx)
comp(model, scores = c(1, 2, 3, 4))## chiSq df pChisq
## 1 22.8 3 5e-05 ***
## n 23.2 3 4e-05 ***
## sqrtN 23.1 3 4e-05 ***
## S1 23.2 3 4e-05 ***
## S2 23.2 3 4e-05 ***
## FH_p=1_q=1 16.7 3 8e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## $tft
## Q Var Z pNorm
## 1 -4.46e+00 7.28e+01 -0.52 0.6
## n -3.77e+02 2.82e+05 -0.71 0.5
## sqrtN -4.01e+01 4.19e+03 -0.62 0.5
## S1 -4.09e+00 3.70e+01 -0.67 0.5
## S2 -4.03e+00 3.58e+01 -0.67 0.5
## FH_p=1_q=1 -2.22e-01 2.59e+00 -0.14 0.9
##
## $scores
## [1] 1 2 3 43 Estimación Paramétrica de la función de supervivencia
La función de supervivencia 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 respuestaya 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. Puede consultar sus respectivas parametrizaciones en https://cran.r-project.org/web/packages/flexsurv/flexsurv.pdf.
| Nombre | dist | Numero de parámetros |
|---|---|---|
| Gamma generalizada (stable) | “gengamma” | 3 |
| Gamma generalizada (original) | “gengamma.orig” | 3 |
| F generalizada (stable) | “genf” | 4 |
| F generalizada (original) | “genf.orig” | 4 |
| Weibull | “weibull” | 2 |
| Gamma | “gamma” | 2 |
| Exponencial | “exp” | 1 |
| Log-logistica | “llogis” | 2 |
| Log-normal | “lnorm” | 2 |
| Gompertz | “gompertz” | 2 |
Todas las distribuciones anteriores y la gran mayoría programadas en R cuenta 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.11dweibull(1, shape = 0.5, scale = 2) #Funcion de Densidad Weibull## [1] 0.17qllogis(0.5, shape = 1, scale = 3) #Funcion Cuantil de log-Logistica## [1] 3rgompertz(10, shape = 1, rate = 1/50) # Simulación Gompertz## [1] 2.3 4.6 4.5 1.9 3.7 4.9 2.5 2.5 4.1 3.43.1 Ejemplo
La base de datos bfeed contiene 927 observaciones de 10 variables, alguna de ellas son los siguientes:
- duration. Duración de la lactancia materna, en Semanas.
- delta. Indicador de la lactancia materna completada(1=Si,0=No)
- race. Raza de la madre(1=Blanco,2=Negro,3=Otro)
- smoke. Indicador si la madre fumo al nacer el bebé (1=Si,0=No)
Se ajustará la distribución exponencial.
data("bfeed") #Carga Los Datos
head(bfeed)## duration delta race poverty smoke alcohol agemth ybirth yschool pc3mth
## 1 16 1 1 0 0 1 24 82 14 0
## 2 1 1 1 0 1 0 26 85 12 0
## 3 4 0 1 0 0 0 25 85 12 0
## 4 3 1 1 0 1 1 21 85 9 0
## 5 36 1 1 0 1 0 22 82 12 0
## 6 36 1 1 0 0 0 18 82 11 0flex <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = "exp") #Ajuste exponencial
flex## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = bfeed,
## dist = "exp")
##
## Estimates:
## est L95% U95% se
## rate 0.05948 0.05570 0.06352 0.00199
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3409, df = 1
## AIC = 6821La 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.
- df: Grados de libertad.
- AIC: Índice AIC.
3.2 Graficación del ajuste
3.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 supervivencia 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 Supervivencia",
conf.int = F, col = 1, col.obs = 2, xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
3.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 Supervivencia y Riesgo Acumulado.
El script calcula las estimaciones de la probabilidad de supervivencia (O Riesgo Acumulado) mediante la función summary() (Paramétrico) y fortify()(No Paramétrico) 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.
3.2.2.1 Curva de Supervivencia
flexgg <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = "exp") %>%
summary(type = "survival") %>% data.frame
kmgg <- survfit(Surv(duration, delta) ~ 1, data = bfeed) %>% fortify
ggplot() + geom_line(aes(time, est, col = "Paramétrico"), data = flexgg) + geom_step(aes(time,
surv, col = "No Paramétrico"), data = kmgg) + labs(x = "Tiempo (Semanas)",
y = "Probabilidad de Supervivencia", col = "Ajustes", title = "Comparación entre curvas de supervivencia")
3.2.2.2 Riesgo Acumulado
La gráfica de riesgo acumulado se realiza utilizando el mismo script utilizando \(y=-log(S_{est}(t))\)
flexch <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = "exp") %>%
summary() %>% data.frame
kmch <- survfit(Surv(duration, delta) ~ 1, data = bfeed) %>% fortify
ggplot() + geom_line(aes(time, -log(est), col = "Paramétrico"), data = flexch) +
geom_step(aes(time, -log(surv), col = "No Paramétrico"), data = kmch) +
labs(x = "Tiempo (Semanas)", y = "Riesgo Acumulado", col = "Ajustes", title = "Comparación entre Riesgo acumulado")
3.3 Comparación de curvas de Supervivencia por grupos.
Dados dos o más grupos, la comparación entre las curvas paramétricas se realiza mediante la función flexsurvreg indicando el nombre de la columna de los grupos igualada a un nivel del mismo en el argumento subset.
Ejemplo.
flex <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = "exp", subset = smoke ==
1)
flex## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = bfeed,
## subset = smoke == 1, dist = "exp")
##
## Estimates:
## est L95% U95% se
## rate 0.0713 0.0632 0.0805 0.0044
##
## N = 270, Events: 263, Censored: 7
## Total time at risk: 3687
## Log-likelihood = -957, df = 1
## AIC = 1917Bajo este método el análisis se realiza individualmente para cada nivel de la variable, y la graficación se realiza utilizando el script de la sección 3.2
Otro método para analizar es ajustar una misma distribución a cada grupos, se deja el siguiente script con los siguientes parámetros modificables:
- bfeed. data.frame que contiene los datos.
- race. Nombre de la covariable de los grupos
- exp. Distribución a Ajustar.
El script funciona agrupando la base de datos(bfeed) por la covariable(race) y aplicando a cada grupo un ajusto paramétrico(gompertz) guardando los resultados en la variable model.
head(bfeed)## duration delta race poverty smoke alcohol agemth ybirth yschool pc3mth
## 1 16 1 1 0 0 1 24 82 14 0
## 2 1 1 1 0 1 0 26 85 12 0
## 3 4 0 1 0 0 0 25 85 12 0
## 4 3 1 1 0 1 1 21 85 9 0
## 5 36 1 1 0 1 0 22 82 12 0
## 6 36 1 1 0 0 0 18 82 11 0model <- bfeed %>% group_by(race) %>% do(flex = flexsurvreg(Surv(duration, delta) ~
1, data = ., dist = "gompertz"))
model## Source: local data frame [3 x 2]
## Groups: <by row>
##
## # A tibble: 3 x 2
## race flex
## * <int> <list>
## 1 1 <S3: flexsurvreg>
## 2 2 <S3: flexsurvreg>
## 3 3 <S3: flexsurvreg>model$flex## [[1]]
## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = ., dist = "gompertz")
##
## Estimates:
## est L95% U95% se
## shape -0.00295 -0.00709 0.00119 0.00211
## rate 0.05928 0.05330 0.06592 0.00321
##
## N = 662, Events: 634, Censored: 28
## Total time at risk: 11278
## Log-likelihood = -2458, df = 2
## AIC = 4920
##
##
## [[2]]
## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = ., dist = "gompertz")
##
## Estimates:
## est L95% U95% se
## shape 0.00425 -0.00987 0.01837 0.00720
## rate 0.05964 0.04483 0.07935 0.00869
##
## N = 117, Events: 113, Censored: 4
## Total time at risk: 1777
## Log-likelihood = -424, df = 2
## AIC = 852
##
##
## [[3]]
## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = ., dist = "gompertz")
##
## Estimates:
## est L95% U95% se
## shape -0.01305 -0.02327 -0.00283 0.00521
## rate 0.09114 0.07375 0.11263 0.00985
##
## N = 148, Events: 145, Censored: 3
## Total time at risk: 1941
## Log-likelihood = -517, df = 2
## AIC = 1039Se utilizar model$flex[[]] para acceder a los ajustes individuales.
Ejemplo
model$flex[[1]]## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = ., dist = "gompertz")
##
## Estimates:
## est L95% U95% se
## shape -0.00295 -0.00709 0.00119 0.00211
## rate 0.05928 0.05330 0.06592 0.00321
##
## N = 662, Events: 634, Censored: 28
## Total time at risk: 11278
## Log-likelihood = -2458, df = 2
## AIC = 49203.3.1 Graficación
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.
3.3.1.1 R-Base
Este script se modifica la variable KM de acuerdo a los parámetros deseados.
3.3.1.1.1 Curva de Supervivencia
KM <- survfit(Surv(duration, delta) ~ race, data = bfeed)
plot(KM, col = 1:length(model$flex), main = "Comparación Paramétrica", xlab = "Tiempo",
ylab = "Probabilidad de Supervivencia")
for (i in 1:length(model$flex)) lines(model$flex[[i]], ci = F, col = i)
3.3.1.1.2 Riesgo Acumulado
KM <- survfit(Surv(duration, delta) ~ race, data = bfeed)
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")
3.3.1.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.
3.3.1.2.1 Curva de Supervivencia
El script extrae de la lista los valores estimados de la curva de supervivencia.
KM <- survfit(Surv(duration, delta) ~ race, data = bfeed)
KMgg <- KM %>% fortify()
Flexgg <- model[[1]] %>% sapply(function(x) cbind(Grupo = levels(KMgg$strata)[x],
data.frame(summary(model[[x, 2]]))), simplify = F) %>% do.call("rbind",
.)
ggplot() + geom_step(aes(time, surv, col = strata, linetype = "No Parametricos"),
KMgg) + geom_line(aes(time, est, col = factor(Grupo), linetype = "Paramétrico"),
Flexgg, size = 0.9) + scale_linetype_manual(values = c("dashed", "solid")) +
labs(col = "Grupos", linetype = "Ajustes", x = "Tiempo(Semanas)", y = "Probabilidad de Supervivencia",
title = "Ajustes Paramétrico vs No Paramétrico")
3.3.1.2.2 Riesgo Acumulado
KM <- survfit(Surv(duration, delta) ~ race, data = bfeed)
KMgg <- KM %>% fortify()
Flexgg <- model[[1]] %>% sapply(function(x) cbind(Grupo = levels(KMgg$strata)[x],
data.frame(summary(model[[x, 2]]))), simplify = F) %>% do.call("rbind",
.)
ggplot() + geom_step(aes(time, -log(surv), col = strata, linetype = "No Parametricos"),
KMgg) + geom_line(aes(time, -log(est), col = factor(Grupo), linetype = "Paramétrico"),
Flexgg, size = 0.9) + scale_linetype_manual(values = c("dashed", "solid")) +
labs(col = "Grupos", linetype = "Ajustes", x = "Tiempo(Semanas)", y = "Riesgo Acumulado",
title = "Ajustes Paramétrico vs No Paramétrico")
3.4 Comparación de modelos.
La comparación entre los modelos ajustados se puede realizar visualmente entre la curva de supervivencia Kaplan-Meier y las curvas paramétricas ajustadas o mediante el índice AIC.
3.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 <- bfeed
Dist <- c("exp", "weibull", "llogis", "gompertz", "genf")
data.Surv <- Surv(datosflex$duration, datosflex$delta)
# Inicia SCRIPT#
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.05948 0.05570 0.06352 0.00199
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3409, df = 1
## AIC = 6821
##
##
## $weibull
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## shape 0.9701 0.9229 1.0197 0.0247
## scale 16.5888 15.4514 17.8100 0.6012
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3409, df = 2
## AIC = 6821
##
##
## $llogis
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## shape 1.4385 1.3634 1.5176 0.0393
## scale 9.8064 9.0547 10.6205 0.3990
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3429, df = 2
## AIC = 6863
##
##
## $gompertz
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## shape -0.004255 -0.007957 -0.000552 0.001889
## rate 0.064008 0.058516 0.070014 0.002929
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3407, df = 2
## AIC = 6817
##
##
## $genf
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## mu 2.49e+00 2.35e+00 2.63e+00 7.12e-02
## sigma 1.13e+00 1.07e+00 1.19e+00 3.12e-02
## Q 4.33e-01 2.26e-01 6.39e-01 1.05e-01
## P 8.36e-04 1.54e-09 4.54e+02 5.63e-03
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3395, df = 4
## AIC = 6797Cada 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.004255 -0.007957 -0.000552 0.001889
## rate 0.064008 0.058516 0.070014 0.002929
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3407, df = 2
## AIC = 68173.4.1.1 R Base
El siguiente script realiza la comparación gráfica utilizando los gráficos bases de R.
3.4.1.1.1 Curva de Supervivencia
plot(model[[1]], ci = F, conf.int = F, lty = 2, main = "Ajuste Paramétrico",
xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
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))
3.4.1.1.2 Riesgo Acumulado
plot(model[[1]], ci = F, conf.int = F, lty = 2, main = "Ajuste Paramétrico",
xlab = "Tiempo", 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))
3.4.1.2 ggplot2
En ggplot2 es un poco más complejo graficarlo, se deja el siguiente script en donde solo se debe haber ejecutado el primer script de comparación gráfica.
3.4.1.2.1 Curva de Supervivencia
KM <- survfit(data.Surv ~ 1, data = datosflex)
KMdist <- KM %>% fortify()
Flexdist <- names(model) %>% sapply(function(x) cbind(Dist = x, data.frame(summary(model[[x]]))),
simplify = F) %>% do.call("rbind", .)
x <- (1:length(Dist)) + 1 #Colores
names(x) <- Dist #Cambiar si se Desea
ggplot() + geom_line(aes(time, est, col = factor(Dist)), Flexdist, size = 0.9,
alpha = 0.5) + geom_step(aes(time, surv, col = "No Parametrico", group = 1),
KMdist, size = 1, alpha = 0.9) + scale_color_manual(values = c(`No Parametrico` = "black",
x)) + labs(col = "Distribuciones", linetype = "Ajustes", x = "Tiempo(Semanas)",
y = "Probabilidad de Supervivencia", title = "Ajustes Paramétrico vs No Paramétrico")
3.4.1.2.2 Riesgo Acumulado
KM <- survfit(data.Surv ~ 1, data = datosflex)
KMdist <- KM %>% fortify()
Flexdist <- names(model) %>% sapply(function(x) cbind(Dist = x, data.frame(summary(model[[x]]))),
simplify = F) %>% do.call("rbind", .)
x <- (1:length(Dist)) + 1 #Colores
names(x) <- Dist #Cambiar si se Desea
ggplot() + geom_line(aes(time, -log(est), col = factor(Dist)), Flexdist, size = 0.9,
alpha = 0.5) + geom_step(aes(time, -log(surv), col = "No Parametrico", group = 1),
KMdist, size = 1, alpha = 0.9) + scale_color_manual(values = c(`No Parametrico` = "black",
x)) + labs(col = "Distribuciones", linetype = "Ajustes", x = "Tiempo(Semanas)",
y = "Riesgo Acumulado", title = "Ajustes Paramétrico vs No Paramétrico")
3.4.2 Í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(duration, delta) ~ 1, data = bfeed, dist = "lnorm")
AIC(model)## [1] 6810BIC(model)## [1] 6819logLik(model)## 'log Lik.' -3403 (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 <- bfeed
data.Surv <- Surv(datosflex$duration, datosflex$delta)
Dist <- c("exp", "weibull", "llogis", "gompertz", "genf")
# Inicia 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 genf
## AIC 6821 6821 6863 6817 6797
## BIC 6825 6831 6872 6827 6816
## LogLik -3409 -3409 -3429 -3407 -3395# Ordenados por menor AIC
AIC.model[, order(AIC.model["AIC", ])]## genf gompertz exp weibull llogis
## AIC 6797 6817 6821 6821 6863
## BIC 6816 6827 6825 6831 6872
## LogLik -3395 -3407 -3409 -3409 -34293.5 Pruebas gráficas para modelos especificos.
La función PPlot() (Ver anexos) 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 | dist |
|---|---|
| Exponencial | exp |
| Weibull | weibull |
| Gumbel | gumbel |
| Log-Normal | lnorm |
| Gamma | gamma |
| Normal | norm |
| Logistica | logis |
| Log-logistica | 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
KM <- survfit(Surv(duration, delta) ~ 1, data = bfeed)
PPlot("pareto", KM)## Warning: Removed 1 rows containing non-finite values (stat_smooth).
PPlot("exp", KM, lm = T) + labs(title = "Papel de Probabilidad Exponencial",
y = "ln(S(t))", x = "t")## Warning: Removed 1 rows containing non-finite values (stat_smooth).
3.6 Gráfico de Bondad de ajuste.
El gráfico de bondad de ajuste se realiza mediante un gráfico de dispersión entre
\[1-\hat{S}_{KM}(t_i) \text{ vs } F_{EMV}(t_i;\hat{\theta})\]
En donde \(\hat{S}_{KM}(t_i)\) es la estimación no paramétrica de la curva de supervivencia y \(F_{EMV}(t_i;\hat{\theta})\) es la estimación 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 argumento de la función flexsurvreg() y survfit()
flex <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = "gompertz") # Ajuste Paramétrico
F_est <- 1 - as.data.frame(summary(flex))[, "est"]
KM_est <- survfit(Surv(duration, delta) ~ 1, bfeed) #Ajuste no-paramétrico
KM_est <- 1 - KM_est$surv
plot(F_est, KM_est)
qplot(F_est, KM_est)
3.6.1 Gráfico de Probabilidad Estabilizado.
El gráfico de probabilidad estabilizado se realiza utilizando el anterior script, anexando las transformaciones en ambos ejes: \(( \frac{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")
3.7 Otros Modelos
El paquete flexsurv no es limitado en cuanto a distribuciones, el usuario puede utilizar sus distribuciones personalizadas si este le puede proveer la función de densidad o la función de riesgo. En R la paquetería actuar contiene muchas distribuciones usadas en la ciencia actuarial, se mostrara de manera de ejemplo como ajustar una distribución Pareto Tipo II o Lomax, cuya función de probabilidad acumulada es:
\[F(x)= 1- \left[1 + \frac{x}{\sigma}\right]^{-\alpha} \qquad x,\alpha,\sigma > 0, \]
Para utilizar una distribución personalizada el usuario debe indicar en el argumento dist en la función flexsurvreg() un objeto de R tipo Lista, el cual debe tener los siguientes elementos:
- name: Nombre de la distribución (pareto) En el que se asume por lo menos tenga una función de densidad (dpareto) o una función de riesgo.
Igual deberán existir una función de distribución acumulada (ppareto), o en su caso H para la función de riego acumulado, en caso de que las formas analíticas no estén disponibles flexsurvreg calculará por medio de integración numérica estas funciones.
Estas funciones deberán ser vectorizadas, y la función de densidad debe aceptar el argumento log, en donde TRUE, devuelve la log-densidad.
pars: Vector de nombres de los parámetros de la variable, \((\alpha,\sigma)\), Deben coincidir con los argumentos de su distribución en R, en el mismo orden
c("shape","scale").location: Nombre del parámetro de localización
c("shape").transforms: una lista de funciones que transforman los parámetros de su rango natural al rango de los reales, por ejemplo
c(log,identity)si el primero solo toma valores positivos y el segundo no tiene restricción.inv.transforms: una lista de sus correspondientes funciones inversas.
inits: Función que proporcionará valores iniciales para el algoritmo de maximización de la estimación máxima verosimilitud(Opcional), Si esta no es proporcionada se deberá indicar cuando se realiza la estimación en flexsurvreg .
Para nuestro ejemplo el objeto tipo lista quedaría de la siguiente manera:
custom.pareto <- list(name = "pareto", pars = c("shape", "scale"), location = c("scale"),
transforms = c(log, log), inv.transforms = c(exp, exp))Para realizar la estimación se utiliza en el argumento dist de la función flexsurvreg nuestra lista definida, seguido del parámetro inits(En caso de no haberse definido), el cual será un vector que proporcionará valores iniciales para el algoritmo de maximización de la estimación máxima verosimilitud.
library(actuar)
flex <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = custom.pareto,
inits = c(3, 3))## Forming integrated rmst function...## Forming integrated mean function...flex## Call:
## flexsurvreg(formula = Surv(duration, delta) ~ 1, data = bfeed,
## dist = custom.pareto, inits = c(3, 3))
##
## Estimates:
## est L95% U95% se
## shape 12.10 4.83 30.31 5.67
## scale 187.01 69.49 503.23 94.45
##
## N = 927, Events: 892, Censored: 35
## Total time at risk: 14996
## Log-likelihood = -3406, df = 2
## AIC = 6817plot(flex, main = "Estimación Pareto Tipo 2 (Lomax)")
plot(flex, main = "Estimación Pareto Tipo 2 (Lomax)", type = "cumhaz")
4 Modelo de Riesgo proporcionales de Cox.
El modelo semiparamétrico de riesgo proporcionales de Cox es realizado con la función coxph() de la paquetería survival, los argumentos de esta función son los siguientes:
formula: Un objeto fórmula
y ~ x, que debe tener un objeto Surv como variable respuesta y a la izquierda del “~”, y del lado derecho las covariables del modelo. Las diferentes formas de incluir covariables se explican a detalle más adelante.data: Objeto data frame que contiene los datos y covariables nombradas en la formula.
method: Nombre del métodos usado para los datos repetidos(“efron”,“breslow”,“exact”), el método “efron” por defecto.
singular.ok: Valor lógico(TRUE/FALSE) que indica el tratamiento de la colinealidad en la matriz del modelo. si es verdadero(TRUE) el programa automáticamente omite las columnas de la matriz que son combinación lineal de otras, en este caso las columnas mostraran NA, y la matriz de covarianzas contienen ceros(Opcional).
La salida de la función es un resumen de la siguiente forma:
Proporciona una tabla que para cada variable del modelo muestra:
- coef: Coeficiente estimado de la beta.
- exp(coef): Exponencial elevado al coeficiente beta estimado.
- se(coef): La desviación estándar de la estimación.
- z: Valor del estadístico para prueba Wald de beta igual a cero.
- p: P-valor de la prueba Wald, beta igual a cero.
Debajo de esta tabla se encuentra el resume el estadístico de la prueba de razón de verosimilitud con sus respectivos grados de libertad y o valor, el numero de datos y observaciones del evento.
4.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.
En los ejemplos utilizaremos una base de datos simulada que se puede descargar en https://www.dropbox.com/s/26p7qj2786qy6va/cox.csv?dl=0
datos <- read.csv("datos/cox.csv") #Leer Datos
str(datos)## 'data.frame': 100 obs. of 6 variables:
## $ Tiempo : num 12.04 0.72 44.65 56.8 29.29 ...
## $ Censura: int 1 1 1 0 1 1 1 1 1 1 ...
## $ Sexo : int 0 1 0 1 1 0 1 1 0 0 ...
## $ Edad : int 42 29 37 50 38 48 28 42 42 41 ...
## $ Ciudad : int 1 4 3 4 1 3 3 3 3 4 ...
## $ IMC : num 19.4 19.9 22.6 26.2 16.5 ...head(datos)## Tiempo Censura Sexo Edad Ciudad IMC
## 1 12.04 1 0 42 1 19
## 2 0.72 1 1 29 4 20
## 3 44.65 1 0 37 3 23
## 4 56.80 0 1 50 4 26
## 5 29.29 1 1 38 1 17
## 6 4.51 1 0 48 3 254.1.1 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 formula.
coxph(Surv(Tiempo, Censura) ~ Edad + IMC, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad + IMC, data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.0377 0.9630 0.0188 -2.00 0.045
## IMC -0.0514 0.9499 0.0243 -2.11 0.035
##
## Likelihood ratio test=9.29 on 2 df, p=0.00963
## n= 100, number of events= 77Las 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 I( ) para realizar estas transformaciones.
Ejemplo
coxph(Surv(Tiempo, Censura) ~ log(Edad) + I(IMC^2) + I(IMC > 20), data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ log(Edad) + I(IMC^2) +
## I(IMC > 20), data = datos)
##
## coef exp(coef) se(coef) z p
## log(Edad) -1.337005 0.262631 0.704341 -1.90 0.058
## I(IMC^2) -0.001245 0.998755 0.000835 -1.49 0.136
## I(IMC > 20)TRUE 0.042861 1.043793 0.411025 0.10 0.917
##
## Likelihood ratio test=9.1 on 3 df, p=0.028
## n= 100, number of events= 774.1.2 Variables Nominales u Ordinales.
Para el uso de variables nominales u ordinales se debe codificar cada factor con un numero entero, posteriormente 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 ortogonales, mientras que en factores no ordenado se utilizará variables dummies.
Para seleccionar la columna del factor se utiliza el operador $, su uso es el siguiente: Nombre_de_la_tabla$Columna.
Ejemplo:
En el siguiente ejemplo tenemos la siguiente 2 variables categóricas codificadas de la siguiente manera:
Sexo: 0-Hombre, 1-Mujer Ciudad: 1-MID, 2-MTY, 3-CDMX, 4-VER
datos <- read.csv("datos/cox.csv") #Leer Datos
head(datos) #Datos Originales## Tiempo Censura Sexo Edad Ciudad IMC
## 1 12.04 1 0 42 1 19
## 2 0.72 1 1 29 4 20
## 3 44.65 1 0 37 3 23
## 4 56.80 0 1 50 4 26
## 5 29.29 1 1 38 1 17
## 6 4.51 1 0 48 3 25datos$Sexo <- factor(x = datos$Sexo, levels = c(0, 1), labels = c("H", "M"))
datos$Ciudad <- factor(x = datos$Ciudad, levels = c(1, 2, 3, 4), labels = c("MID",
"MTY", "CDMX", "VER"))
head(datos) #Datos tranformados## Tiempo Censura Sexo Edad Ciudad IMC
## 1 12.04 1 H 42 MID 19
## 2 0.72 1 M 29 VER 20
## 3 44.65 1 H 37 CDMX 23
## 4 56.80 0 M 50 VER 26
## 5 29.29 1 M 38 MID 17
## 6 4.51 1 H 48 CDMX 25coxph(Surv(Tiempo, Censura) ~ Sexo, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Sexo, data = datos)
##
## coef exp(coef) se(coef) z p
## SexoM -0.723 0.485 0.234 -3.09 0.002
##
## Likelihood ratio test=9.52 on 1 df, p=0.00203
## n= 100, number of events= 77coxph(Surv(Tiempo, Censura) ~ Ciudad, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Ciudad, data = datos)
##
## coef exp(coef) se(coef) z p
## CiudadMTY -0.129 0.879 0.319 -0.41 0.69
## CiudadCDMX 0.183 1.201 0.351 0.52 0.60
## CiudadVER -0.117 0.890 0.325 -0.36 0.72
##
## Likelihood ratio test=1.01 on 3 df, p=0.8
## n= 100, number of events= 77Por 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 de la forma: I(Variable=="Etiqueta").
Ejemplo:
coxph(Surv(Tiempo, Censura) ~ I(Ciudad == "MTY") + I(Ciudad == "MID"), data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ I(Ciudad == "MTY") +
## I(Ciudad == "MID"), data = datos)
##
## coef exp(coef) se(coef) z p
## I(Ciudad == "MTY")TRUE -0.127587 0.880217 0.271562 -0.47 0.64
## I(Ciudad == "MID")TRUE -0.000284 0.999716 0.289904 0.00 1.00
##
## Likelihood ratio test=0.26 on 2 df, p=0.879
## n= 100, number of events= 774.1.3 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:
- ** * **: Este operador agregar las variables principales y todas las posibles interacciones posibles en las variables.
Ejemplo para el modelo \(\lambda(t|X_i)=\lambda_0(t)*exp(\beta_0 Edad+\beta_1Sexo+\beta_3Edad*Sexo)\) la formula en la función sería la siguiente:
coxph(Surv(Tiempo, Censura) ~ Edad * Sexo, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad * Sexo, data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.0651 0.9370 0.0294 -2.21 0.027
## SexoM -1.8001 0.1653 1.5089 -1.19 0.233
## Edad:SexoM 0.0252 1.0255 0.0382 0.66 0.509
##
## Likelihood ratio test=17 on 3 df, p=0.000708
## n= 100, number of events= 77Todas las posibles interacciones entre Edad, IMC y Ciudad:
coxph(Surv(Tiempo, Censura) ~ Edad * IMC * Ciudad, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad * IMC * Ciudad,
## data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -3.82e-01 6.82e-01 3.80e-01 -1.01 0.31
## IMC -4.39e-01 6.45e-01 5.67e-01 -0.77 0.44
## CiudadMTY -6.56e+00 1.42e-03 1.87e+01 -0.35 0.73
## CiudadCDMX 6.82e+00 9.18e+02 2.32e+01 0.29 0.77
## CiudadVER -4.82e-01 6.18e-01 1.72e+01 -0.03 0.98
## Edad:IMC 1.11e-02 1.01e+00 1.45e-02 0.77 0.44
## Edad:CiudadMTY 1.60e-01 1.17e+00 4.76e-01 0.34 0.74
## Edad:CiudadCDMX -7.15e-02 9.31e-01 5.85e-01 -0.12 0.90
## Edad:CiudadVER 6.88e-02 1.07e+00 4.55e-01 0.15 0.88
## IMC:CiudadMTY 1.60e-01 1.17e+00 7.09e-01 0.23 0.82
## IMC:CiudadCDMX -3.22e-01 7.25e-01 9.47e-01 -0.34 0.73
## IMC:CiudadVER -1.22e-01 8.85e-01 6.70e-01 -0.18 0.86
## Edad:IMC:CiudadMTY -3.87e-03 9.96e-01 1.79e-02 -0.22 0.83
## Edad:IMC:CiudadCDMX 4.33e-03 1.00e+00 2.37e-02 0.18 0.86
## Edad:IMC:CiudadVER 5.76e-04 1.00e+00 1.75e-02 0.03 0.97
##
## Likelihood ratio test=23.9 on 15 df, p=0.066
## n= 100, number of events= 77- “:” : Este operador realizar únicamente la interacción señalada.
Ejemplo para el modelo en donde se considera la interacción entre Sexo y Edad y las variables principales:
coxph(Surv(Tiempo, Censura) ~ Edad + Sexo + Edad:Sexo, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad + Sexo + Edad:Sexo,
## data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.0651 0.9370 0.0294 -2.21 0.027
## SexoM -1.8001 0.1653 1.5089 -1.19 0.233
## Edad:SexoM 0.0252 1.0255 0.0382 0.66 0.509
##
## Likelihood ratio test=17 on 3 df, p=0.000708
## n= 100, number of events= 77Modelo utilizando solo la interacción entre Sexo y Edad:
coxph(Surv(Tiempo, Censura) ~ Edad:Sexo, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad:Sexo, data = datos)
##
## coef exp(coef) se(coef) z p
## Edad:SexoH -0.0379 0.9628 0.0186 -2.04 0.0412
## Edad:SexoM -0.0576 0.9440 0.0198 -2.91 0.0036
##
## Likelihood ratio test=15.6 on 2 df, p=0.000413
## n= 100, number of events= 77- “^” : Este operador realiza las posibles interacciones hasta el nivel indicado.
Modelo utilizando las variables principales y sus interacciones dobles y triples de las variables Edad, IMC y Sexo.
coxph(Surv(Tiempo, Censura) ~ (Edad + IMC + Sexo)^3, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ (Edad + IMC + Sexo)^3,
## data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.30321 0.73845 0.18869 -1.61 0.11
## IMC -0.47459 0.62214 0.30423 -1.56 0.12
## SexoM 2.35107 10.49676 9.88937 0.24 0.81
## Edad:IMC 0.01004 1.01009 0.00762 1.32 0.19
## Edad:SexoM -0.10877 0.89693 0.25710 -0.42 0.67
## IMC:SexoM -0.12145 0.88564 0.39541 -0.31 0.76
## Edad:IMC:SexoM 0.00420 1.00421 0.01011 0.42 0.68
##
## Likelihood ratio test=28.4 on 7 df, p=0.00019
## n= 100, number of events= 77Modelo utilizando las variables principales y las interacciones dobles de las variables Edad, IMC y Ciudad.
coxph(Surv(Tiempo, Censura) ~ (Edad + IMC + Ciudad)^2, data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ (Edad + IMC + Ciudad)^2,
## data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.36499 0.69420 0.16864 -2.16 0.030
## IMC -0.41248 0.66201 0.24319 -1.70 0.090
## CiudadMTY -2.40837 0.08996 2.55753 -0.94 0.346
## CiudadCDMX 2.76508 15.88031 3.18734 0.87 0.386
## CiudadVER -0.93221 0.39368 2.59589 -0.36 0.720
## Edad:IMC 0.01048 1.01053 0.00608 1.72 0.085
## Edad:CiudadMTY 0.05575 1.05733 0.06103 0.91 0.361
## Edad:CiudadCDMX 0.02916 1.02959 0.06871 0.42 0.671
## Edad:CiudadVER 0.08350 1.08708 0.06552 1.27 0.203
## IMC:CiudadMTY 0.00492 1.00493 0.07034 0.07 0.944
## IMC:CiudadCDMX -0.14752 0.86285 0.09663 -1.53 0.127
## IMC:CiudadVER -0.10391 0.90130 0.07322 -1.42 0.156
##
## Likelihood ratio test=23.8 on 12 df, p=0.0219
## n= 100, number of events= 77De manera general se puede utilizar la función I( ) para crear nuevas variables y/o interacciones deseadas.
Ejemplo utilizando como covariables el logaritmo de Edad, el IMC-Centrado y la variable dummie de Mérida-Hombre
coxph(Surv(Tiempo, Censura) ~ log(Edad) + I(IMC - mean(IMC)) + I((Ciudad ==
"MID") * (Sexo == "H")), data = datos)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ log(Edad) + I(IMC - mean(IMC)) +
## I((Ciudad == "MID") * (Sexo == "H")), data = datos)
##
## coef exp(coef) se(coef) z
## log(Edad) -1.5389 0.2146 0.7170 -2.15
## I(IMC - mean(IMC)) -0.0541 0.9473 0.0246 -2.20
## I((Ciudad == "MID") * (Sexo == "H")) 0.5348 1.7072 0.3259 1.64
## p
## log(Edad) 0.032
## I(IMC - mean(IMC)) 0.028
## I((Ciudad == "MID") * (Sexo == "H")) 0.101
##
## Likelihood ratio test=11.3 on 3 df, p=0.0104
## n= 100, number of events= 774.2 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(Tiempo, Censura) ~ Edad:Sexo, data = datos)
summary(modelo)## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad:Sexo, data = datos)
##
## n= 100, number of events= 77
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Edad:SexoH -0.0379 0.9628 0.0186 -2.04 0.0412 *
## Edad:SexoM -0.0576 0.9440 0.0198 -2.91 0.0036 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Edad:SexoH 0.963 1.04 0.928 0.998
## Edad:SexoM 0.944 1.06 0.908 0.981
##
## Concordance= 0.656 (se = 0.038 )
## Rsquare= 0.144 (max possible= 0.997 )
## Likelihood ratio test= 15.6 on 2 df, p=0.000413
## Wald test = 14.1 on 2 df, p=0.000853
## Score (logrank) test = 14.6 on 2 df, p=0.000678Los resultados son una versión más detallada de la función coxph, en esta se encuentran agregada una tabla con los coeficientes estimados de \(e^{\beta_i}\) con sus respectivos intervalos de confianza.
Además se agregan estadísticos de Concordancia y \(R^2\); así como los estadísticos de prueba, grados de libertad y p-valores de las pruebas de hipótesis de Razón de Verosimilitud, Prueba de Wald y Prueba de Puntajes (score test).
4.3 Comparación entre modelos.
La comparación entre los modelos se realiza con la función anova(), esta función realiza la prueba de Razón de Verosimilitud entre dos modelos. Los argumentos de la función son los dos modelos ajustados a comparar.
La salida de la función es una tabla en donde muestra para los dos modelos la log-verosimilitud, el estadísticos Chisq, los grados de libertad Df y el p-valor P(>|Chi|) de la prueba.
Ejemplo
modeloCompleto <- coxph(Surv(Tiempo, Censura) ~ Edad + IMC, data = datos)
modeloReducido <- coxph(Surv(Tiempo, Censura) ~ Edad, data = datos)
anova(modeloCompleto, modeloReducido)## Analysis of Deviance Table
## Cox model: response is Surv(Tiempo, Censura)
## Model 1: ~ Edad + IMC
## Model 2: ~ Edad
## loglik Chisq Df P(>|Chi|)
## 1 -284
## 2 -286 4.38 1 0.036 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.4 Validación de Supuestos.
4.4.1 Riesgos Proporcionales
La validación de riesgos proporcionales se realiza analizando los residuos Schoenfeld en contra de alguna transformación del tiempo, la función cox.zph() realiza la prueba de hipótesis de correlación para los residuales Schoenfeld de cada covariable con alguna transformación del tiempo, los argumentos de la función son:
fit: Un objeto coxph resultado de un ajuste del modelo de regresión de Cox.
tranform: Tipo de transformación del tiempo, posibles valores son “km”, “rank”, “identity” o una función.
La salida de la función es una tabla en donde para cada covariable calcula el coeficiente rho, el estadístico de prueba chisq y su respectivo p-valor, además incluye una prueba global.
Ejemplo
modelo <- coxph(Surv(Tiempo, Censura) ~ log(Edad) * Sexo, data = datos)
riesgo <- cox.zph(modelo)
riesgo## rho chisq p
## log(Edad) 0.0736 0.483 0.487
## SexoM -0.0588 0.257 0.612
## log(Edad):SexoM 0.0628 0.292 0.589
## GLOBAL NA 3.707 0.295La validación del supuesto se puede hacer de forma gráfica utilizando la función ggcoxzph() al objeto cox.zph, la selección de la covariable a mostrar se le indica en forma de vector al parámetro var. Esta función pertenece a la librería survminer.
Ejemplo
riesgo## rho chisq p
## log(Edad) 0.0736 0.483 0.487
## SexoM -0.0588 0.257 0.612
## log(Edad):SexoM 0.0628 0.292 0.589
## GLOBAL NA 3.707 0.295ggcoxzph(riesgo)
ggcoxzph(riesgo, var = c("log(Edad)", "SexoM"))
4.4.2 Análisis Residual
El análisis residual se lleva a cabo gráficamente con la función ggcoxdiagnostics de la librería survminer, esta función muestra los gráficos de diagnóstico de un modelo calculado con la función coxph. Los argumentos de la función son:
- fit: Un modelo calculado con la función
cox.ph. - type: Tipo de residual representado en el eje Y. Los posibles valores son: “martingale”, “deviance”,“score”,“schoenfeld”,“dfbeta”,“dfbetas” y “scaledsch”.
- ox.scale : Forma de representar el eje X. los posibles valores son : “linear.predictions” para Y Predichos, “observation.id” para id de la observación y “time” para Tiempo.
- title: Titulo del gráfico.
- subtitle: Subtitulo del gráfico.
- caption: Leyenda del gráfico.
Ejemplo.
modelo <- coxph(Surv(Tiempo, Censura) ~ log(Edad) * Sexo, data = datos)
ggcoxdiagnostics(modelo, type = "schoenfeld", ox.scale = "time")
ggcoxdiagnostics(modelo, type = "dfbetas", ox.scale = "observation.id", title = "Residuales dfbetas",
subtitle = "dfbetas vs id", caption = "Análisis residual")
Un problema de esta función es que gráfica todas las variables, si se desea obtener una gráfica individual para cada covariable o ciertas covariables se deberá utilizar la función en los anexos llamada ggcoxdiagnostics_var esta función toma los siguientes argumentos:
- ** *: Objeto ggplot generado con
ggcoxdiagnostics() - var : Un vector de las variables a mostrar
Ejemplo.
p <- ggcoxdiagnostics(modelo, type = "dfbetas", ox.scale = "observation.id",
title = "Residuales dfbetas", subtitle = "dfbetas vs id", caption = "Análisis residual")
ggcoxdiagnostics_var(p, c("SexoM", "log(Edad)"))
ggcoxdiagnostics_var(p, "log(Edad):SexoM")
4.5 Predicción de curva de supervivencia.
La función survfit calcula la función de supervivencia pronosticada para un modelo de riesgos proporcionales de Cox. Los argumentos de esta función son:
- formula: Un modelo calculado con la función coxph.
- newdata: Un data frame con los mismo nombres que aparecen en el modelo ajustado. La curva producida representara la cohorte con los valores proporcionados.
- conf.int: Nivel del intervalo de confianza(.95 por defecto).
- se.fit: Indicador si se deben calcular el error estándar.
- conf.type: Tipo de calculado para los intervalos de confianza (“none”,“plain”,“log”(defecto),“log-log”).
La salida de esta función son las misma que se ven en la sección Estimador de Kaplan-Meier y Fleming-Harrington.
Ejemplo:
Utilizando el modelo con variables Edad y IMC, se calcularan las curva de supervivencia para una persona de Edad 50 con IMC de 15 y de otra persona con Edad 60 y con IMC de 20.
(modelo <- coxph(Surv(Tiempo, Censura) ~ Edad + IMC, data = datos))## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad + IMC, data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.0377 0.9630 0.0188 -2.00 0.045
## IMC -0.0514 0.9499 0.0243 -2.11 0.035
##
## Likelihood ratio test=9.29 on 2 df, p=0.00963
## n= 100, number of events= 77datosnuevos <- data.frame(Edad = c(50, 60), IMC = c(15, 20)) #Datos a Estimar
datosnuevos## Edad IMC
## 1 50 15
## 2 60 20curva <- survfit(modelo, newdata = datosnuevos)
curva## Call: survfit(formula = modelo, newdata = datosnuevos)
##
## n events median 0.95LCL 0.95UCL
## 1 100 77 16.8 10.4 76.7
## 2 100 77 31.1 16.8 NALas curvas estimadas son de tipo Survfit por lo que se puede estimar su media y cuantiles utilizando las funciones print() con el argumento print.rmean=True y quantile().
print(curva, print.rmean = T)## Call: survfit(formula = modelo, newdata = datosnuevos)
##
## n events *rmean *se(rmean) median 0.95LCL 0.95UCL
## 1 100 77 28.3 3.15 16.8 10.4 76.7
## 2 100 77 49.1 7.55 31.1 16.8 NA
## * restricted mean with upper limit = 110quantile(curva)## $quantile
## 25 50 75
## 1 7.2 17 36
## 2 12.7 31 90
##
## $lower
## 25 50 75
## 1 4.6 10 20
## 2 7.2 17 36
##
## $upper
## 25 50 75
## 1 20 77 NA
## 2 NA NA NA4.5.1 Graficar la curva de supervivencia.
La curva de supervivencia se gráfica con la función ggsurvplot de 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 y se recomienda revisar los demás utilizando el comando help("ggsurvplot").
- fit: Objeto tipo survfit.
- data: Un conjunto de datos utilizado para ajustar curvas de supervivencia.
- fun: Transformación de la curva de supervivencia(Opcional), las posibles opciones son: “event” para los eventos acumulados, “cumhaz” para el riesgo acumulado y “pct” para la curva de supervivencia en porcentaje.
- conf.int: Indicador para graficar los intervalos de confianza.
- title: Titulo
- xlab: Eje x
- ylab: Eje Y
- legends.lab: Vector de nombres para identificar las curvas.
- legend.title : Titulo de la leyenda.
(modelo <- coxph(Surv(Tiempo, Censura) ~ Edad + IMC, data = datos))## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad + IMC, data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.0377 0.9630 0.0188 -2.00 0.045
## IMC -0.0514 0.9499 0.0243 -2.11 0.035
##
## Likelihood ratio test=9.29 on 2 df, p=0.00963
## n= 100, number of events= 77datosnuevos <- data.frame(Edad = c(50, 60), IMC = c(15, 20)) #Datos a Estimar
curva <- survfit(modelo, newdata = datosnuevos, conf.type = "log-log")
ggsurvplot(fit = curva, data = datos, conf.int = T, title = "Curva de Supervivencia",
xlab = "Tiempo", ylab = "Probabilidad de supervivencia", legend.labs = c("EDAD=50 y IMC=15",
"EDAD=60 y IMC=20"), legend.title = "Grupos: ")
4.5.2 Estimación de Riesgo Base
La estimación de riesgo base representa la particular predicción de la curva de supervivencia utilizando las covariables igual a 0.
(modelo <- coxph(Surv(Tiempo, Censura) ~ Edad + IMC, data = datos))## Call:
## coxph(formula = Surv(Tiempo, Censura) ~ Edad + IMC, data = datos)
##
## coef exp(coef) se(coef) z p
## Edad -0.0377 0.9630 0.0188 -2.00 0.045
## IMC -0.0514 0.9499 0.0243 -2.11 0.035
##
## Likelihood ratio test=9.29 on 2 df, p=0.00963
## n= 100, number of events= 77datosnuevo <- data.frame(Edad = 0, IMC = 0)
RiesgoBase <- survfit(modelo, datosnuevo)
summary(RiesgoBase)## Call: survfit(formula = modelo, newdata = datosnuevo)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.38 100 1 8.71e-01 1.59e-01 6.08e-01 1
## 0.72 99 1 7.57e-01 2.36e-01 4.11e-01 1
## 0.78 97 1 6.56e-01 2.89e-01 2.77e-01 1
## 0.96 96 1 5.67e-01 3.23e-01 1.86e-01 1
## 1.81 95 1 4.89e-01 3.43e-01 1.24e-01 1
## 1.97 94 1 4.21e-01 3.51e-01 8.22e-02 1
## 2.47 93 1 3.62e-01 3.50e-01 5.44e-02 1
## 3.00 92 1 3.11e-01 3.42e-01 3.60e-02 1
## 3.02 91 1 2.67e-01 3.30e-01 2.37e-02 1
## 3.15 90 1 2.28e-01 3.13e-01 1.54e-02 1
## 3.83 89 1 1.94e-01 2.95e-01 1.00e-02 1
## 4.25 88 1 1.65e-01 2.74e-01 6.42e-03 1
## 4.30 87 1 1.41e-01 2.53e-01 4.11e-03 1
## 4.41 86 1 1.19e-01 2.32e-01 2.62e-03 1
## 4.51 85 1 1.01e-01 2.11e-01 1.67e-03 1
## 4.60 84 1 8.54e-02 1.91e-01 1.06e-03 1
## 4.76 83 1 7.20e-02 1.72e-01 6.64e-04 1
## 5.36 82 2 5.08e-02 1.37e-01 2.56e-04 1
## 5.37 80 1 4.25e-02 1.21e-01 1.57e-04 1
## 5.42 79 2 2.94e-02 9.37e-02 5.75e-05 1
## 6.17 76 1 2.42e-02 8.14e-02 3.36e-05 1
## 7.18 74 1 1.98e-02 7.02e-02 1.92e-05 1
## 7.19 73 1 1.62e-02 6.02e-02 1.09e-05 1
## 7.65 71 1 1.31e-02 5.12e-02 6.08e-06 1
## 7.82 70 1 1.06e-02 4.34e-02 3.37e-06 1
## 8.42 69 1 8.51e-03 3.66e-02 1.85e-06 1
## 8.55 68 1 6.82e-03 3.07e-02 1.00e-06 1
## 8.97 67 1 5.45e-03 2.57e-02 5.39e-07 1
## 9.00 66 1 4.35e-03 2.13e-02 2.87e-07 1
## 9.63 64 1 3.44e-03 1.76e-02 1.50e-07 1
## 10.37 62 1 2.69e-03 1.44e-02 7.62e-08 1
## 10.54 61 1 2.10e-03 1.17e-02 3.82e-08 1
## 10.55 60 1 1.63e-03 9.43e-03 1.87e-08 1
## 12.03 59 1 1.25e-03 7.57e-03 9.09e-09 1
## 12.04 58 1 9.61e-04 6.04e-03 4.34e-09 1
## 12.18 57 1 7.33e-04 4.79e-03 2.04e-09 1
## 12.34 56 1 5.55e-04 3.76e-03 9.37e-10 1
## 12.67 55 1 4.16e-04 2.93e-03 4.19e-10 1
## 13.02 54 1 3.10e-04 2.27e-03 1.84e-10 1
## 13.68 53 1 2.29e-04 1.74e-03 7.78e-11 1
## 15.12 50 1 1.65e-04 1.31e-03 3.11e-11 1
## 15.58 49 1 1.18e-04 9.71e-04 1.20e-11 1
## 16.10 48 1 8.39e-05 7.16e-04 4.53e-12 1
## 16.27 47 1 5.90e-05 5.23e-04 1.67e-12 1
## 16.81 46 1 4.12e-05 3.79e-04 6.01e-13 1
## 17.12 45 1 2.86e-05 2.73e-04 2.14e-13 1
## 17.38 44 1 1.97e-05 1.95e-04 7.45e-14 1
## 17.60 43 1 1.34e-05 1.37e-04 2.49e-14 1
## 18.02 42 1 9.02e-06 9.59e-05 8.10e-15 1
## 19.83 39 1 5.85e-06 6.45e-05 2.35e-15 1
## 20.44 38 1 3.75e-06 4.30e-05 6.60e-16 1
## 21.50 37 1 2.38e-06 2.83e-05 1.81e-16 1
## 21.54 36 1 1.50e-06 1.84e-05 4.80e-17 1
## 22.66 34 1 9.15e-07 1.17e-05 1.19e-17 1
## 24.15 33 1 5.52e-07 7.32e-06 2.82e-18 1
## 24.99 31 1 3.21e-07 4.43e-06 6.04e-19 1
## 25.43 30 1 1.83e-07 2.62e-06 1.22e-19 1
## 25.77 29 1 1.02e-07 1.52e-06 2.33e-20 1
## 29.20 28 1 5.63e-08 8.68e-07 4.25e-21 1
## 29.29 26 1 2.97e-08 4.76e-07 6.92e-22 1
## 30.55 25 1 1.50e-08 2.50e-07 9.80e-23 1
## 31.10 24 1 7.42e-09 1.29e-07 1.33e-23 1
## 33.22 22 1 3.46e-09 6.25e-08 1.52e-24 1
## 36.33 20 1 1.52e-09 2.86e-08 1.48e-25 1
## 44.65 16 1 5.79e-10 1.14e-08 9.85e-27 1
## 45.04 15 1 2.03e-10 4.20e-09 5.22e-28 1
## 49.65 14 1 6.78e-11 1.47e-09 2.42e-29 1
## 50.34 13 1 2.09e-11 4.76e-10 9.15e-31 1
## 53.84 11 1 4.69e-12 1.13e-10 1.31e-32 1
## 56.31 10 1 8.33e-13 2.15e-11 8.93e-35 1
## 65.52 7 1 7.03e-14 1.98e-12 6.74e-38 1
## 76.70 6 1 4.73e-15 1.46e-13 2.92e-41 1
## 81.30 5 1 2.20e-16 7.40e-15 4.91e-45 1
## 90.21 4 1 6.89e-18 2.54e-16 3.20e-49 1
## 110.19 1 1 5.87e-23 2.79e-21 2.06e-63 1cbind(Tiempo = RiesgoBase$time, RiesgoAcumulado = RiesgoBase$cumhaz)## Tiempo RiesgoAcumulado
## [1,] 0.38 0.14
## [2,] 0.72 0.28
## [3,] 0.78 0.42
## [4,] 0.96 0.57
## [5,] 1.81 0.72
## [6,] 1.97 0.86
## [7,] 2.47 1.02
## [8,] 3.00 1.17
## [9,] 3.02 1.32
## [10,] 3.15 1.48
## [11,] 3.83 1.64
## [12,] 4.25 1.80
## [13,] 4.30 1.96
## [14,] 4.41 2.13
## [15,] 4.51 2.29
## [16,] 4.60 2.46
## [17,] 4.76 2.63
## [18,] 5.36 2.98
## [19,] 5.37 3.16
## [20,] 5.42 3.53
## [21,] 6.16 3.53
## [22,] 6.17 3.72
## [23,] 6.73 3.72
## [24,] 7.18 3.92
## [25,] 7.19 4.13
## [26,] 7.43 4.13
## [27,] 7.65 4.34
## [28,] 7.82 4.55
## [29,] 8.42 4.77
## [30,] 8.55 4.99
## [31,] 8.97 5.21
## [32,] 9.00 5.44
## [33,] 9.58 5.44
## [34,] 9.63 5.67
## [35,] 9.72 5.67
## [36,] 10.37 5.92
## [37,] 10.54 6.17
## [38,] 10.55 6.42
## [39,] 12.03 6.68
## [40,] 12.04 6.95
## [41,] 12.18 7.22
## [42,] 12.34 7.50
## [43,] 12.67 7.78
## [44,] 13.02 8.08
## [45,] 13.68 8.38
## [46,] 13.78 8.38
## [47,] 14.16 8.38
## [48,] 15.12 8.71
## [49,] 15.58 9.04
## [50,] 16.10 9.39
## [51,] 16.27 9.74
## [52,] 16.81 10.10
## [53,] 17.12 10.46
## [54,] 17.38 10.83
## [55,] 17.60 11.22
## [56,] 18.02 11.62
## [57,] 18.07 11.62
## [58,] 18.53 11.62
## [59,] 19.83 12.05
## [60,] 20.44 12.49
## [61,] 21.50 12.95
## [62,] 21.54 13.41
## [63,] 22.66 13.90
## [64,] 24.15 14.41
## [65,] 24.75 14.41
## [66,] 24.99 14.95
## [67,] 25.43 15.51
## [68,] 25.77 16.09
## [69,] 29.20 16.69
## [70,] 29.27 16.69
## [71,] 29.29 17.33
## [72,] 30.55 18.02
## [73,] 31.10 18.72
## [74,] 33.13 18.72
## [75,] 33.22 19.48
## [76,] 34.19 19.48
## [77,] 36.33 20.30
## [78,] 37.97 20.30
## [79,] 42.22 20.30
## [80,] 43.93 20.30
## [81,] 44.65 21.27
## [82,] 45.04 22.32
## [83,] 49.65 23.42
## [84,] 50.34 24.59
## [85,] 50.81 24.59
## [86,] 53.84 26.08
## [87,] 56.31 27.81
## [88,] 56.80 27.81
## [89,] 57.40 27.81
## [90,] 65.52 30.29
## [91,] 76.70 32.99
## [92,] 81.30 36.05
## [93,] 90.21 39.52
## [94,] 98.35 39.52
## [95,] 100.51 39.52
## [96,] 110.19 51.195 Ejemplos - Tiempos de muerte de pacientes de trasplante de riñón
Analizaremos la base de datos kidtran. Los datos sobre la hora de la muerte de 863 pacientes con trasplante de riñón. A todos los pacientes se les realizó su trasplante en el Centro de trasplantes de la Universidad Estatal de Ohio durante el período 1982-1992. El tiempo máximo de seguimiento para este estudio fue de 9.47 años. Los pacientes fueron censurados si se mudaron de Colón (perdidos durante el seguimiento) o si estaban vivos el 30 de junio de 1992. En la muestra, había 432 hombres blancos, 92 hombres negros, 280 mujeres blancas y 59 mujeres negras. Las edades de los pacientes en el trasplante variaron de 9,5 meses a 74,5 años con una edad promedio de 42,8 años. Setenta y tres (16.9%) de los hombres blancos, 14 (15.2%) de los hombres negros, 39 (13.9%) de las mujeres blancas y 14 (23.7%) de las mujeres negras murieron antes del final del estudio.
Las variables contiene las siguientes variables:
- obs. Numero de observación.
- time. Tiempo hasta la muerte o en estudio.
- delta. Indicador de muerte(1=Vivo,0=Muerte)
- gender. 1=Hombre, 2=Mujer
- race. Raza (1=Blanco,2=Negro)
- age. Edad en Años.
data("kidtran")
str(kidtran)## 'data.frame': 863 obs. of 6 variables:
## $ obs : int 1 2 3 4 5 6 7 8 9 10 ...
## $ time : int 1 5 7 9 13 13 17 20 26 26 ...
## $ delta : int 0 0 1 0 0 0 1 0 1 1 ...
## $ gender: int 1 1 1 1 1 1 1 1 1 1 ...
## $ race : int 1 1 1 1 1 1 1 1 1 1 ...
## $ age : int 46 51 55 57 45 43 47 65 55 44 ...head(kidtran)## obs time delta gender race age
## 1 1 1 0 1 1 46
## 2 2 5 0 1 1 51
## 3 3 7 1 1 1 55
## 4 4 9 0 1 1 57
## 5 5 13 0 1 1 45
## 6 6 13 0 1 1 435.1 Análisis Exploratorio
Se calculará medidas descriptivas sobre el número de observaciones censuradas de acuerdo a las variables. Estas estadísticas se utilizaran utilizando las funciones table y prob.table en conjunto a with para utilizar las columnas de la base de datos como variables.
with(kidtran, table(delta)) #Censura## delta
## 0 1
## 723 140with(kidtran, prop.table(table(delta)))## delta
## 0 1
## 0.84 0.16El numero de observaciones por muerte por trasplante de riñón es bajo a comparación con el número de observaciones perdidas o vivos.
Las covariables que representan grupos las convertiremos a tipo factor para un mejor uso y visualización.
kidtran$gender <- factor(kidtran$gender, labels = c("Hombre", "Mujer"))
kidtran$race <- factor(kidtran$race, labels = c("Blanco", "Negro"))with(kidtran, table(gender, delta)) #genero## delta
## gender 0 1
## Hombre 437 87
## Mujer 286 53with(kidtran, table(race, delta)) #Raza## delta
## race 0 1
## Blanco 600 112
## Negro 123 28with(kidtran, hist(age))
La variable edad se analizará creando una nueva variable edadg en donde se dividirá en 3 grupos: Menores hasta 25, Mayores a 25 hasta 50 y mayores a 50.
kidtran$edadg <- cut(kidtran$age, c(0, 25, 50, 75))
with(kidtran, table(edadg, delta))## delta
## edadg 0 1
## (0,25] 100 7
## (25,50] 406 69
## (50,75] 217 645.2 Estimación No-Parametríca
5.2.1 Análisis Exploratorio.
Se estimará la curva de supervivencia para diferentes variables con la finalidad de encontrar alguna diferencia notable entre estas.
survfit(Surv(time, delta) ~ gender, kidtran, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Género",
conf.int = T, legend.title = "Género", legend.labs = c("Hombres", "Mujeres"))
survfit(Surv(time, delta) ~ race, kidtran, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Raza",
conf.int = T, legend.title = "Raza", legend.labs = c("Blanco", "Negro"))
survfit(Surv(time, delta) ~ edadg, kidtran, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Edad",
conf.int = T, legend.title = "Edad", legend.labs = c("0-25", "25-50", "50+"))
A primera vista los grupos de edad presentan una clara diferencia entre sus curvas de supervivencia, en donde las personas con mayor edad presentan un mayor riesgo, la curva por Raza presenta al final una ligera diferencia.
Ahora se analizarán las curvas de supervivencia por covariables utilizando subgrupos de otra covariable utilizando el argumento facet.by de ggsurvplot.
survfit(Surv(time, delta) ~ gender + edadg, kidtran, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia entre Género por Edad", conf.int = T,
facet.by = "edadg", legend.title = "Género", panel.labs = list(edadg = c("0-25",
"25-50", "75+")), short.panel.labs = T)
survfit(Surv(time, delta) ~ race + edadg, kidtran, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia entre Raza por Edad", conf.int = T, facet.by = "edadg",
legend.title = "Raza", panel.labs = list(edadg = c("0-25", "25-50",
"75+")), short.panel.labs = T)
survfit(Surv(time, delta) ~ race + gender, kidtran, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia entre Género por Raza", conf.int = T,
facet.by = "gender", legend.title = "Raza", short.panel.labs = T)
survfit(Surv(time, delta) ~ race + gender + edadg, kidtran, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia entre Género por Raza", conf.int = F,
facet.by = "edadg")
Las gráficas sugieren que el género y la raza no son factores de riesgo para la muerte, sin embargo, la edad gráficamente parecer un factor de riesgo, esta hipótesis se pondrá a prueba a continuación.
5.2.2 Pruebas de Hipótesis.
A continuación, realizaremos las respectivas pruebas de hipótesis:
- Género y Raza
- Género y Raza estratificando los grupos de Edad.
- Grupos de Edad.
- Tendencia entre los Grupos de Edad.
Todas las pruebas de hipótesis se harán con un nivel de confianza del 95%.
survdiff(Surv(time, delta) ~ gender + race, kidtran, rho = 1)## Call:
## survdiff(formula = Surv(time, delta) ~ gender + race, data = kidtran,
## rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## gender=Hombre, race=Blanco 432 66.3 62.71 0.2114 0.462
## gender=Hombre, race=Negro 92 12.5 13.16 0.0342 0.042
## gender=Mujer, race=Blanco 280 35.2 42.62 1.2777 2.131
## gender=Mujer, race=Negro 59 12.5 8.04 2.4179 2.827
##
## Chisq= 4.3 on 3 degrees of freedom, p= 0.227survdiff(Surv(time, delta) ~ gender + race + strata(edadg), kidtran, rho = 1)## Call:
## survdiff(formula = Surv(time, delta) ~ gender + race + strata(edadg),
## data = kidtran, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## gender=Hombre, race=Blanco 432 64.9 63.91 0.0149 0.034
## gender=Hombre, race=Negro 92 12.3 14.01 0.1972 0.248
## gender=Mujer, race=Blanco 280 35.3 38.74 0.3009 0.490
## gender=Mujer, race=Negro 59 12.4 8.25 2.0389 2.451
##
## Chisq= 2.9 on 3 degrees of freedom, p= 0.413survdiff(Surv(time, delta) ~ edadg, kidtran, rho = 1)## Call:
## survdiff(formula = Surv(time, delta) ~ edadg, data = kidtran,
## rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## edadg=(0,25] 107 6.43 18.7 8.08 10.60
## edadg=(25,50] 475 61.87 72.8 1.63 4.22
## edadg=(50,75] 281 58.24 35.0 15.34 23.40
##
## Chisq= 27.8 on 2 degrees of freedom, p= 9.19e-07ten(Surv(time, delta) ~ edadg, data = kidtran) %>% comp()## chiSq df pChisq
## 1 28.5 2 7e-07 ***
## n 22.9 2 1e-05 ***
## sqrtN 25.6 2 3e-06 ***
## S1 27.8 2 9e-07 ***
## S2 27.8 2 9e-07 ***
## FH_p=1_q=1 23.6 2 8e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## $tft
## Q Var Z pNorm
## 1 -2.35e+00 7.61e+01 -0.27 0.8
## n 1.13e+03 2.89e+07 0.21 0.8
## sqrtN 5.80e+00 4.36e+04 0.03 1.0
## S1 -1.40e+00 6.14e+01 -0.18 0.9
## S2 -1.37e+00 6.11e+01 -0.18 0.9
## FH_p=1_q=1 -7.15e-01 7.87e-01 -0.81 0.4
##
## $scores
## [1] 1 2 3Las pruebas de hipótesis no rechazan la hipótesis de diferencia de curvas por lo que podemos suponer que las variables género y raza no son factores de riesgo para la muerte en el trasplante de riñón. Mientras que se rechaza la hipótesis de igualdad de curvas de supervivencia entre los grupos de edad, sin embargo, no se rechazó la prueba de tendencia entre estas.
El análisis paramétrico se llevará a cabo de dos maneras, la primera haciendo un ajuste de distribución para para grupo de edades y la segunda realizando un modelo semiparamétrico de Cox.
5.3 Modelo Paramétrico
El análisis Paramétrico se llevará a cabo separando la base de datos en los 3 grupos y aplicando para cada uno el script de comparación de modelos.
Para cada grupo se estimará vía máxima verosimilitud los parámetros de las siguientes distribuciones: Weibull, Loglogistica, Gompertz, LogNormal. Eligiéremos al modelo que presente un mejor AIC, la validación se hará de manera gráfica comparando las curvas de supervivencia y el gráfico estabilizado de Bondad de Ajuste.
5.3.1 Grupo (0-25]
datosflex <- kidtran %>% filter(edadg == "(0,25]") #Filtrar Base de datos
Dist <- c("weibull", "llogis", "gompertz", "lnorm")
data.Surv <- Surv(datosflex$time, datosflex$delta)
# Inicia SCRIPT#
model <- sapply(Dist, function(x) flexsurvreg(data.Surv ~ 1, dist = x), USE.NAMES = T,
simplify = F)
AIC.model <- sapply(model, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = x$loglik),
simplify = T)
AIC.model[, order(AIC.model["AIC", ])]## gompertz lnorm llogis weibull
## AIC 155 155 155 156
## BIC 160 161 161 161
## LogLik -75 -76 -76 -76El Modelo con menor AIC es el Gompertz.
Ahora validaremos con el Gráfico de Probabilidad Estabilizada.
flex <- flexsurvreg(Surv(time, delta) ~ 1, data = datosflex, dist = "gompertz") # Ajuste Paramétrico
F_est <- 1 - as.data.frame(summary(flex))[, "est"]
KM_est <- survfit(Surv(time, delta) ~ 1, datosflex) #Ajuste no-paramétrico
KM_est <- 1 - KM_est$surv
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")
Se tiene tendencia de línea recta, sin embargo al tener un gran porcentaje de observaciones censuradas se presentan líneas verticales.
Ahora visualizaremos la curva de supervivencia y riesgo acumulado.
flexgg <- flexsurvreg(Surv(time, delta) ~ 1, data = datosflex, dist = "gompertz") %>%
summary(type = "survival") %>% data.frame
kmgg <- survfit(Surv(time, delta) ~ 1, data = datosflex) %>% fortify
ggplot() + geom_line(aes(time, est, col = "Paramétrico"), data = flexgg) + geom_step(aes(time,
surv, col = "No Paramétrico"), data = kmgg) + labs(x = "Tiempo (Semanas)",
y = "Probabilidad de Supervivencia", col = "Ajustes", title = "Comparación entre curvas de supervivencia")
ggplot() + geom_line(aes(time, -log(est), col = "Paramétrico"), data = flexgg) +
geom_step(aes(time, -log(surv), col = "No Paramétrico"), data = kmgg) +
labs(x = "Tiempo (Días)", y = "Riesgo Acumulado", col = "Ajustes", title = "Comparación entre curvas de Riesgo Acumulado")
Gráficamente el modelo parece adecuado por lo cual optaremos por no desecharlo. Las estimaciones de los parámetros y sus intervalos de confianza son:
round(model$gompertz$res, 8)## est L95% U95% se
## shape -0.00157 -2.8e-03 -0.00038 6.1e-04
## rate 0.00013 4.7e-05 0.00035 6.6e-055.3.2 Grupo (25-50]
datosflex <- kidtran %>% filter(edadg == "(25,50]") #Filtrar Base de datos
Dist <- c("weibull", "llogis", "gompertz", "lnorm")
data.Surv <- Surv(datosflex$time, datosflex$delta)
# Inicia SCRIPT#
model <- sapply(Dist, function(x) flexsurvreg(data.Surv ~ 1, dist = x), USE.NAMES = T,
simplify = F)
AIC.model <- sapply(model, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = x$loglik),
simplify = T)
AIC.model[, order(AIC.model["AIC", ])]## weibull lnorm llogis gompertz
## AIC 1400 1400 1400 1408
## BIC 1408 1408 1408 1417
## LogLik -698 -698 -698 -702El Modelo con menor AIC es Weibull, aunque no con gran diferencia.
Ahora validaremos con el Gráfico de Probabilidad Estabilizada.
flex <- flexsurvreg(Surv(time, delta) ~ 1, data = datosflex, dist = "weibull") # Ajuste Paramétrico
F_est <- 1 - as.data.frame(summary(flex))[, "est"]
KM_est <- survfit(Surv(time, delta) ~ 1, datosflex) #Ajuste no-paramétrico
KM_est <- 1 - KM_est$surv
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 Estabiizada")
El ajuste tiene tendencia de línea aunque se puede notar ligeras desviaciones en la cola.
flexgg <- flexsurvreg(Surv(time, delta) ~ 1, data = datosflex, dist = "weibull") %>%
summary(type = "survival") %>% data.frame
kmgg <- survfit(Surv(time, delta) ~ 1, data = datosflex) %>% fortify
ggplot() + geom_line(aes(time, est, col = "Paramétrico"), data = flexgg) + geom_step(aes(time,
surv, col = "No Paramétrico"), data = kmgg) + labs(x = "Tiempo (Semanas)",
y = "Probabilidad de Supervivencia", col = "Ajustes", title = "Comparación entre curvas de Supervivencia")
ggplot() + geom_line(aes(time, -log(est), col = "Paramétrico"), data = flexgg) +
geom_step(aes(time, -log(surv), col = "No Paramétrico"), data = kmgg) +
labs(x = "Tiempo (Días)", y = "Riesgo Acumulado", col = "Ajustes", title = "Comparación entre curvas de Riesgo Acumulado")
El ajuste es bueno a pesar de que parece subestimar la probabilidad de supervivencia. Los parámetros finales del modelo son:
round(model$weibull$res, 5)## est L95% U95% se
## shape 6.8e-01 5.5e-01 8.4e-01 7.4e-02
## scale 2.3e+04 1.2e+04 4.5e+04 7.7e+035.3.3 Grupo (50-75]
datosflex <- kidtran %>% filter(edadg == "(50,75]") #Filtrar Base de datos
Dist <- c("weibull", "llogis", "gompertz", "lnorm")
data.Surv <- Surv(datosflex$time, datosflex$delta)
# Inicia SCRIPT#
model <- sapply(Dist, function(x) flexsurvreg(data.Surv ~ 1, dist = x), USE.NAMES = T,
simplify = F)
AIC.model <- sapply(model, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = x$loglik),
simplify = T)
AIC.model[, order(AIC.model["AIC", ])]## lnorm weibull llogis gompertz
## AIC 1197 1197 1197 1206
## BIC 1204 1205 1205 1214
## LogLik -597 -597 -597 -601El Modelo con menor AIC es Lognormal, aunque no con gran diferencia.
Ahora validaremos con el Gráfico de Probabilidad Estabilizada.
flex <- flexsurvreg(Surv(time, delta) ~ 1, data = datosflex, dist = "lnorm") # Ajuste Paramétrico
F_est <- 1 - as.data.frame(summary(flex))[, "est"]
KM_est <- survfit(Surv(time, delta) ~ 1, datosflex) #Ajuste no-paramétrico
KM_est <- 1 - KM_est$surv
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 Estabiizada")
El ajuste tiene tendencia de línea aunque se puede notar ligeras desviaciones en la cola.
flexgg <- flexsurvreg(Surv(time, delta) ~ 1, data = datosflex, dist = "lnorm") %>%
summary(type = "survival") %>% data.frame
kmgg <- survfit(Surv(time, delta) ~ 1, data = datosflex) %>% fortify
ggplot() + geom_line(aes(time, est, col = "Paramétrico"), data = flexgg) + geom_step(aes(time,
surv, col = "No Paramétrico"), data = kmgg) + labs(x = "Tiempo (Semanas)",
y = "Probabilidad de Supervivencia", col = "Ajustes", title = "Comparación entre curvas de Supervivencia")
ggplot() + geom_line(aes(time, -log(est), col = "Paramétrico"), data = flexgg) +
geom_step(aes(time, -log(surv), col = "No Paramétrico"), data = kmgg) +
labs(x = "Tiempo (Días)", y = "Riesgo Acumulado", col = "Ajustes", title = "Comparación entre curvas de Riesgo Acumulado")
El ajuste se subestima ligeramente en los primeros 2200 días, y a partir de ahí el estimador no paramétrico presenta un disminución notable aunque no demasiado.
Los parámetros del modelo final son:
round(model$lnorm$res, 3)## est L95% U95% se
## meanlog 9.0 8.3 9.6 0.33
## sdlog 2.7 2.2 3.3 0.265.4 Modelo Semiparamétrico
El modelo semiparamétrico de Cox permite utilizar la variable age sin categorizar por grupos, utilizaremos la variable original para modelar y validaremos utilizando las gráficas de residuales.
5.4.1 Selección de Variables
Comenzaremos con un modelo completo contemplando las variables age, race y gender, e iremos eliminando las variables no significativas.
coxph(Surv(time, delta) ~ age + race + gender, kidtran)## Call:
## coxph(formula = Surv(time, delta) ~ age + race + gender, data = kidtran)
##
## coef exp(coef) se(coef) z p
## age 0.05104 1.05236 0.00718 7.11 1.2e-12
## raceNegro 0.11636 1.12341 0.21151 0.55 0.58
## genderMujer 0.02654 1.02689 0.17490 0.15 0.88
##
## Likelihood ratio test=57.1 on 3 df, p=2.5e-12
## n= 863, number of events= 140coxph(Surv(time, delta) ~ age + race, kidtran)## Call:
## coxph(formula = Surv(time, delta) ~ age + race, data = kidtran)
##
## coef exp(coef) se(coef) z p
## age 0.05095 1.05227 0.00716 7.12 1.1e-12
## raceNegro 0.11671 1.12379 0.21150 0.55 0.58
##
## Likelihood ratio test=57 on 2 df, p=4.12e-13
## n= 863, number of events= 140coxph(Surv(time, delta) ~ age, kidtran)## Call:
## coxph(formula = Surv(time, delta) ~ age, data = kidtran)
##
## coef exp(coef) se(coef) z p
## age 0.05107 1.05239 0.00714 7.16 8.3e-13
##
## Likelihood ratio test=56.7 on 1 df, p=4.97e-14
## n= 863, number of events= 140cox.final <- coxph(Surv(time, delta) ~ age, kidtran)
summary(cox.final)## Call:
## coxph(formula = Surv(time, delta) ~ age, data = kidtran)
##
## n= 863, number of events= 140
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.05107 1.05239 0.00714 7.16 8.3e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.05 0.95 1.04 1.07
##
## Concordance= 0.671 (se = 0.026 )
## Rsquare= 0.064 (max possible= 0.87 )
## Likelihood ratio test= 56.7 on 1 df, p=4.97e-14
## Wald test = 51.2 on 1 df, p=8.3e-13
## Score (logrank) test = 53.4 on 1 df, p=2.67e-13La variable age fue la única significativa en el modelo.
5.4.2 Validación de Supuestos.
5.4.2.1 Riesgos Proporcionales.
cox.zph(cox.final)## rho chisq p
## age 0.0941 1.14 0.285ggcoxzph(cox.zph(cox.final))
El valor-p no es menor a .05 por lo que no podemos rechazar la hipótesis de riesgos proporcionales.
5.4.2.2 Linealidad.
ggcoxdiagnostics(cox.final, type = "martingale")
Al gráficas los residuales martingala y graficar la recta a través de ellos vemos como se asemeja a una línea recta muy cercana a cero, la robustez de las pruebas consideramos que la desviaciones que se tienen en realidad no son significativas y que se cumple el supuesto de linealidad de la variable.
5.4.2.3 Residuales Dfbeta
ggcoxdiagnostics(cox.final, type = "dfbeta")
La gráfica se observa una serie de puntos que se alejan de los demás puntos, sin embargo el cambio en el coeficiente es muy pequeño.
5.4.3 Predicción
Por último, se presenta la estimación de la curva de supervivencia para las edades de 25, 45, 55 y 75.
datosnuevo <- data.frame(age = c(25, 45, 55, 75))
Pred <- survfit(cox.final, newdata = datosnuevo)
ggsurvplot(Pred, title = "Curva de Supervivencia por Edad")
6 Anexos
6.1 Función PPlot
PPlot <- function(dist, KM, lm = F, ...) {
time <- KM$time
st <- KM$surv
x <- switch(dist, exp = time, weibull = log(time), gumbel = time, lnorm = log(time),
gamma = sqrt(time), norm = time, logis = time, llogis = log(time), pareto = log(time))
y <- switch(dist, exp = log(st), weibull = log(-log(st)), gumbel = log(-log(st)),
lnorm = qnorm(1 - st), gamma = qnorm(1 - st), norm = qnorm(1 - st),
logis = log((1 - st)/st), llogis = log((1 - st)/st), pareto = log(st))
qplot(x, y, main = paste("Papel de Probabilidad", dist)) + geom_line(stat = "smooth",
method = "lm", alpha = lm)
}6.2 Función coxdiagnostics_var
ggcoxdiagnostics_var <- function(p, var) {
g <- p
g$data <- g$data[g$data$covariate == var, ]
g
}6.3 Leyendas en R-Base
Para añadir leyendas a los gráficos de R-Base, puede usarse la función R legend() con los siguientes argumentos:
- x : La posición de la leyenda, “bottomright”, “bottom”, “bottomleft”, “left”, “topleft”, “top”, “topright”, “right” and “center”.
- legend : Texto de la leyenda, puede ser un vector.
- title: Titulo de la leyenda.
- fill : Color del borde de la leyenda.
- col : Color de las lineas o puntos, Puede ser vector.
- bg : Color del fondo del cuadro.
- lty: Tipo de linea, puede ser un vector.
Ejemplo
flex <- flexsurvreg(Surv(duration, delta) ~ 1, data = bfeed, dist = "exp") #Ajuste exponencial
plot(flex, ci = F, type = "cumhaz", conf.int = F, main = "Riesgo Acumulado")
legend("topright", legend = c("Kaplan Meier", "Exponencial"), col = c(1, 2),
lty = 1)
6.4 Base de datos (datos)
| Tiempo | Censura | Sexo | Edad | Ciudad | IMC |
|---|---|---|---|---|---|
| 12.04 | 1 | H | 42 | MID | 19 |
| 0.72 | 1 | M | 29 | VER | 20 |
| 44.65 | 1 | H | 37 | CDMX | 23 |
| 56.80 | 0 | M | 50 | VER | 26 |
| 29.29 | 1 | M | 38 | MID | 17 |
| 4.51 | 1 | H | 48 | CDMX | 25 |
| 10.54 | 1 | M | 28 | CDMX | 23 |
| 20.44 | 1 | M | 42 | CDMX | 30 |
| 1.97 | 1 | H | 42 | CDMX | 24 |
| 24.99 | 1 | H | 41 | VER | 22 |
| 0.72 | 0 | H | 50 | VER | 19 |
| 7.18 | 1 | H | 46 | MID | 22 |
| 100.51 | 0 | M | 37 | VER | 31 |
| 6.16 | 0 | M | 29 | VER | 17 |
| 45.04 | 1 | H | 50 | CDMX | 23 |
| 8.55 | 1 | H | 42 | MID | 24 |
| 22.66 | 1 | M | 39 | CDMX | 25 |
| 25.43 | 1 | M | 40 | MTY | 25 |
| 33.13 | 0 | H | 33 | MTY | 29 |
| 18.02 | 1 | H | 46 | MTY | 26 |
| 53.84 | 1 | M | 28 | MTY | 26 |
| 0.38 | 1 | M | 40 | CDMX | 17 |
| 8.42 | 1 | H | 38 | VER | 21 |
| 6.17 | 1 | H | 32 | MID | 22 |
| 3.02 | 1 | H | 33 | MID | 17 |
| 6.73 | 0 | H | 41 | CDMX | 17 |
| 1.81 | 1 | H | 36 | MTY | 23 |
| 12.03 | 1 | M | 31 | VER | 25 |
| 24.75 | 0 | M | 44 | CDMX | 19 |
| 34.19 | 0 | M | 41 | CDMX | 26 |
| 76.70 | 1 | M | 46 | MTY | 27 |
| 5.36 | 1 | H | 27 | MID | 26 |
| 16.10 | 1 | M | 38 | MTY | 21 |
| 18.07 | 0 | H | 42 | MID | 18 |
| 12.34 | 1 | H | 35 | VER | 16 |
| 7.19 | 1 | H | 34 | VER | 27 |
| 0.78 | 1 | H | 30 | CDMX | 20 |
| 36.33 | 1 | M | 39 | MTY | 30 |
| 57.40 | 0 | M | 55 | CDMX | 26 |
| 12.67 | 1 | H | 42 | MID | 20 |
| 19.83 | 1 | H | 34 | VER | 30 |
| 8.97 | 1 | H | 44 | MTY | 20 |
| 50.34 | 1 | M | 37 | MTY | 17 |
| 29.20 | 1 | H | 46 | VER | 28 |
| 5.42 | 1 | H | 29 | MTY | 21 |
| 0.96 | 1 | M | 31 | MTY | 23 |
| 42.22 | 0 | M | 56 | MTY | 30 |
| 98.35 | 0 | H | 39 | VER | 24 |
| 16.27 | 1 | H | 42 | CDMX | 23 |
| 9.58 | 0 | M | 42 | MTY | 19 |
| 49.65 | 1 | M | 39 | MID | 26 |
| 4.41 | 1 | H | 47 | MTY | 31 |
| 21.54 | 0 | H | 43 | MID | 22 |
| 10.37 | 1 | H | 33 | VER | 24 |
| 90.21 | 1 | M | 43 | MTY | 16 |
| 25.77 | 1 | M | 41 | MTY | 30 |
| 17.38 | 1 | H | 32 | VER | 22 |
| 4.25 | 1 | M | 43 | MTY | 27 |
| 9.63 | 1 | M | 35 | MID | 18 |
| 3.00 | 1 | H | 40 | MTY | 21 |
| 4.76 | 1 | M | 41 | VER | 17 |
| 4.30 | 1 | H | 40 | MTY | 26 |
| 33.22 | 1 | M | 46 | VER | 29 |
| 12.18 | 1 | M | 40 | VER | 16 |
| 21.50 | 1 | H | 40 | MID | 25 |
| 5.37 | 1 | H | 32 | MID | 21 |
| 7.65 | 1 | M | 44 | VER | 29 |
| 43.93 | 0 | H | 45 | VER | 29 |
| 24.15 | 1 | H | 41 | MTY | 27 |
| 17.60 | 1 | M | 35 | MID | 30 |
| 7.82 | 1 | H | 37 | MID | 23 |
| 50.81 | 0 | M | 34 | MID | 24 |
| 13.78 | 0 | M | 36 | CDMX | 28 |
| 13.02 | 1 | M | 33 | VER | 16 |
| 17.12 | 1 | M | 38 | VER | 29 |
| 31.10 | 1 | M | 47 | VER | 25 |
| 7.43 | 0 | M | 38 | MTY | 17 |
| 5.36 | 1 | H | 44 | CDMX | 19 |
| 15.58 | 1 | M | 40 | MID | 30 |
| 5.42 | 1 | H | 37 | CDMX | 17 |
| 18.53 | 0 | M | 34 | MTY | 21 |
| 37.97 | 0 | M | 41 | MTY | 22 |
| 16.81 | 1 | H | 46 | VER | 29 |
| 81.30 | 1 | H | 47 | MID | 29 |
| 15.12 | 1 | M | 36 | VER | 17 |
| 3.83 | 1 | H | 31 | VER | 19 |
| 13.68 | 1 | M | 34 | CDMX | 30 |
| 29.27 | 0 | M | 41 | MID | 22 |
| 3.15 | 1 | M | 38 | CDMX | 20 |
| 110.19 | 1 | M | 27 | VER | 28 |
| 56.31 | 1 | M | 38 | MID | 17 |
| 9.00 | 1 | M | 41 | CDMX | 20 |
| 14.16 | 0 | M | 30 | MTY | 27 |
| 4.60 | 1 | M | 30 | VER | 19 |
| 10.55 | 1 | H | 37 | MTY | 29 |
| 65.52 | 1 | M | 49 | MTY | 30 |
| 30.55 | 1 | H | 47 | MID | 30 |
| 9.72 | 0 | M | 47 | MTY | 25 |
| 21.54 | 1 | H | 41 | MTY | 31 |
| 2.47 | 1 | H | 45 | MTY | 26 |