Comparación de curvas de supervivencia con datos censurados

FACULTAD DE MATEMÁTICAS ESTADÍSTICA E INVESTIGACIÓN OPERATIVA Trabajo Fin de Grado Comparación de curvas de supervivencia con datos censurados Elena de las Heras Jiménez Doble Grado en Matemáticas y Estadística Dirigido por: D. Juan Manuel Muñoz Pichardo Sevilla, Julio de 2021 Índice general Resumen ....................................... iv Abstract ....................................... v Índice de Figuras .................................. vi Índice de Tablas ................................... viii 1. Introducción 1 1.1. Breve referencia histórica ........................... 2 1.2. Estructura del trabajo ............................ 2 2. Conceptos básicos 4 2.1. Función de supervivencia ........................... 4 2.1.1. Estimador de Kaplan-Meier ..................... 5 2.1.2. Estimador actuarial ......................... 6 2.1.3. Estimador de Turnbull ........................ 9 2.1.4. Estimador de Fleming-Harrington ................. 11 2.2. Tasas de fallo ................................. 12 2.2.1. Tasa de fallo instantánea o función de riesgo ............ 12 2.2.2. Tasa de fallo acumulado ....................... 12 2.3. Tipo de datos censurados ........................... 13 2.3.1. Censura de tipo I ........................... 14 2.3.2. Censura de tipo II .......................... 14 2.3.3. Censura aleatoria ........................... 14 2.3.3.1. Censura a la derecha .................... 14 2.3.3.2. Censura a la izquierda ................... 15 2.3.3.3. Datos doblemente censurados ............... 16 2.3.4. Datos censurados por intervalo ................... 16 2.4. Teoría de proceso de conteo (counting process theory)........... 16 2.5. Test de permutación ............................. 17 i 2.6. Teoría de Martingala ............................. 18 2.7. Algoritmo EM ................................ 19 2.8. Variables aleatorias intercambiables ..................... 19 2.9. Método de Montecarlo ............................ 20 3. Comparación de curvas de supervivencia 22 3.1. Técnicas generales ............................... 22 3.1.1. El test log-rank ponderado ..................... 22 3.1.2. Prueba de tendencia (trend-testing)................ 24 3.2. Técnicas específicas para datos censurados ................. 26 3.2.1. Datos censurados a la derecha .................... 27 3.2.1.1. Enfoque mediante el proceso de conteo ......... 27 3.2.1.2. Enfoque de permutación ................. 28 3.2.2. Datos censurados por intervalos ................... 29 3.2.2.1. Enfoque de permutación ................. 29 3.2.2.2. Enfoque basado en un vector de puntuación ....... 31 4. Comparación de curvas de supervivencia en R 33 4.1. Librerías para el análisis de supervivencia .................. 33 4.1.1. Librería survival ........................... 33 4.1.1.1. Estimador de Kaplan-Meier ................ 34 4.1.1.2. Estimador de Fleming-Harrington ............. 38 4.1.1.3. Estimador de la tasa de fallo acumulado ......... 42 4.1.1.4. Comparación de curvas de supervivencia ........ 46 4.1.1.4.1. Comparaciones múltiples ............. 49 4.1.2. Librería survMisc ........................... 51 4.1.2.1. Comparación de curvas de supervivencia ......... 51 4.1.3. Librería KMsurv ........................... 58 4.1.3.1. Estimador actuarial .................... 58 4.1.4. Librería gte .............................. 63 4.1.4.1. Estimador de Turnbull ................... 63 4.1.5. Librería coin ............................. 66 4.1.5.1. Comparación de curvas de supervivencia ......... 66 4.2. Librería FHtest ................................ 71 4.2.1. Comparación de curvas de supervivencia .............. 72 ii 4.2.1.1. Datos censurados a la derecha ............... 72 4.2.1.1.1. Enfoque mediante el proceso de conteo ..... 72 4.2.1.1.2. Enfoque de permutación ............. 80 4.2.1.2. Datos censurados por intervalos .............. 83 4.2.1.2.1. Enfoque de permutación ............. 84 4.2.1.2.2. Enfoque basado en un vector de puntuación . . 92 Conclusión 99 Bibliografía 105 iii Resumen El análisis de supervivencia es una rama de la Estadística que tiene como objetivo analizar o modelizar el tiempo que tarda en ocurrir uno o más eventos de interés como la muerte de un paciente con cáncer o los fallos en los motores de unas máquinas de una industria textil. Para conseguir dicho objetivo, existe una extensa gama de herramientas entre las que se encuentra la comparación de curvas de supervivencia. El propósito de este trabajo es describir la teoría y la metodología de diversas técnicas de comparación de curvas de supervivencia con datos censurados. Además, con objeto de ilustrar la utilidad de estas técnicas se usarán librerías de R que contienen funciones que llevan a cabo dichas técnicas y se aplicarán sobre datos reales. Palabras clave: análisis de supervivencia, datos censurados por intervalos, datos censurados a la derecha, función de supervivencia, prueba de tendencia, test log-rank ponderado, test de permutación, procesos de conteo. iv Abstract Survival analysis is a branch of Statistics that aims to analyse or model the time it takes for one or more events of interest to occur, such as the death of a cancer patient or the failure of machine motors in a textile industry. In order to achieve this goal, a wide range of tools are available, including the comparison of survival curves. The purpose of this dissertation is to describe the theory and methodology of tests for comparing survival curves with censored data. Furthermore, in order to illustrate the usefulness of these techniques, R libraries containing functions that perform these techniques will be used and applied on real data. Keywords: survival analysis, interval-censored data, right-censored data, survival function, trend-testing, weighted log-rank test, permutation tests, counting processes. v Índice de figuras 2.1. Ejemplo de funciones de supervivencia ................... 4 2.2. Estimador de Kaplan-Meier de la función de supervivencia ........ 6 2.3. Estimador actuarial de la función de supervivencia ............ 7 2.4. Estimador de Turnbull de la función de supervivencia ........... 11 2.5. Estimador de Fleming-Harrington de la función de supervivencia ..... 12 2.6. Ejemplo de funciones de riesgo ........................ 13 4.1. Estimador de Kaplan-Meier para los datos AML .............. 38 4.2. Estimador de Fleming-Harrington para los datos oropharynx ....... 42 4.3. Estimador de la tasa de fallo acumulado para los datos rats ........ 45 4.4. Comparación de curvas de supervivencia con ggsurvplot .......... 49 4.5. Prueba de tendencia con la función comp .................. 53 4.6. Prueba de tendencia de manera gráfica ................... 57 4.7. Estimador actuarial para los datos Nurshome ............... 62 4.8. Estimador de Turnbull para los datos simul ................ 66 4.9. Curvas de supervivencia de los pacientes de lung .............. 69 4.10. Curvas de supervivencia de los pacientes de lung por género ....... 71 4.11. Curvas de supervivencia de los pacientes de ovarian ............ 76 4.12. Curvas de supervivencia de los pacientes de bmt .............. 79 4.13. Curvas de supervivencia de los pacientes de bcos .............. 88 4.14. Curvas de supervivencia de los pacientes de illust3 ............ 91 4.15. Curvas de supervivencia de los usuarios de duser por género ....... 95 4.16. Curvas de supervivencia asociadas a cada período de duser ........ 97 vi Índice de tablas 2.1. Datos usados para el estimador actuarial .................. 7 2.2. Tabla de vida o estimador actuarial ..................... 8 3.1. Ejemplos de funciones pesos ω(ti)...................... 24 4.1. Datos AML .................................. 35 4.2. Datos AML implementados en R ...................... 36 4.3. Estimador de Kaplan-Meier para el grupo Maintained ........... 37 4.4. Estimador de Kaplan-Meier para el grupo Nonmaintained ........ 37 4.5. Datos del fichero de oropharynx ....................... 39 4.6. Estimador de Fleming-Harrington para el grupo de pacientes con radioterapia sola (estándar) ............................. 41 4.7. Estimador de Fleming-Harrington para el grupo de pacientes con radioterapia junto a un agente quimioterapéutico (test) .............. 41 4.8. Datos rats del paquete survival ....................... 43 4.9. Estimación de la función de riesgo acumulado para el grupo tratado con fármaco .................................... 44 4.10. Estimación de la función de riesgo acumulado para el grupo de control . 44 4.11. Datos colon del paquete survival ....................... 47 4.12. Datos larynx del paquete KMsurv ...................... 53 4.13. Primera tabla de la salida comp ....................... 54 4.14. Segunda tabla de la salida comp ....................... 55 4.15. Datos Nurshome ............................... 59 4.16. Estimador actuarial del conjunto de datos Nurshome ........... 61 4.17. Datos simul del paquete gte ......................... 64 4.18. Estimador de Turnbull del conjunto de datos simul ............ 65 4.19. Datos lung del paquete survival ....................... 68 4.20. Datos ovarian del paquete survival ..................... 74 vii 4.21. Datos bmt del paquete KMsurv ....................... 77 4.22. Datos bcos del paquete interval ....................... 86 4.23. Datos illust3 del paquete FHtest ...................... 89 4.24. Datos duser del paquete FHtest ....................... 93 viii 2.1. FUNCIÓN DE SUPERVIVENCIA ˆ S(ti+1) = ˆ S(ti)ˆ S(ti+1 |ti) = i Y j=0 rj−dj rj : probabilidad de estar vivo en el instante ti+1, i.e., al inicio del siguiente intervalo. Al igual que ocurría con el estimador de Kaplan-Meier, debe hallarse un intervalo de confianza a partir del error estándar que se calculará tal y como se ha expuesto en la Sección 2.1.1 mediante la aproximación de Greenwood. A continuación se adjuntará una representación de la función de supervivencia estimada mediante el método actuarial de los siguientes datos: nsubs nlost nevent 0-2 927 2 77 2-3 848 3 71 3-5 774 6 119 5-7 649 9 75 7-11 565 7 109 11-17 449 5 148 17-25 296 3 107 25-37 186 0 74 37-53 112 0 85 53-NA 27 0 27 Tabla 2.1: Datos extraídos del ejemplo que aparece en la ayuda de R de la orden lifetab() del paquete KMsurv (Klein y Moeschberger (2012)) La estimación gráfica siguiente se ha llevado a cabo a partir de la Tabla 2.2 que muestra la información necesaria para la construcción de una tabla de vida. 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 t S(t) estimada Fuente propia Figura 2.3: Estimación de la función de supervivencia mediante el método actuarial realizada con R CAPÍTULO 2. CONCEPTOS BÁSICOS 7 2.1. FUNCIÓN DE SUPERVIVENCIA Tabla 2.2: Tabla de vida o estimador actuarial nsubs nlost nrisk nevent surv pdf hazard se.surv se.pdf se.hazard 0-2 927 2 926.0 77 1.0000000 0.0415767 0.0433803 0.0000000 0.0045368 0.0049390 2-3 848 3 846.5 71 0.9168467 0.0769003 0.0875462 0.0090737 0.0087684 0.0103799 3-5 774 6 771.0 119 0.8399463 0.0648208 0.0836261 0.0120584 0.0055430 0.0076391 5-7 649 9 644.5 75 0.7103048 0.0413288 0.0617792 0.0149472 0.0045695 0.0071200 7-11 565 7 561.5 109 0.6276472 0.0304602 0.0537475 0.0159669 0.0027313 0.0051183 11-17 449 5 446.5 148 0.5058065 0.0279430 0.0662192 0.0165929 0.0020898 0.0053347 17-25 296 3 294.5 107 0.3381483 0.0153573 0.0554979 0.0158122 0.0013853 0.0052313 25-37 186 0 186.0 74 0.2152897 0.0071377 0.0413870 0.0138261 0.0007904 0.0046604 37-53 112 0 112.0 85 0.1296368 0.0061491 0.0764388 0.0113583 0.0006305 0.0065601 53-NA 27 0 27.0 27 0.0312517 NA NA 0.0059119 NA NA CAPÍTULO 2. CONCEPTOS BÁSICOS 8 2.1. FUNCIÓN DE SUPERVIVENCIA Lo que representan cada una de las columnas de la tabla anterior es (véase Klein y Moeschberger (2012)): nsubs número de sujetos que entran en los intervalos y que no han experimentado el suceso en estudio nlost número de individuos censurados nrisk número de individuos en riesgo de sufrir el suceso en estudio nevent número de individuos que experimentaron el suceso surv función de supervivencia estimada al principio de los intervalos pdf función de densidad estimada en el punto medio de los intervalos hazard función de riesgo estimada en el punto medio de los intervalos se.surv desviación estándar de la función de supervivencia estimada al principio de los intervalos se.pdf desviación estándar de la función de densidad estimada en el punto medio de los intervalos se.hazard desviación estándar de la función de riesgo estimada en el punto medio de los intervalos 2.1.3. Estimador de Turnbull Cuando se trabaja con datos censurados por intervalos, el método de Kaplan-Meier no se puede aplicar. Ante esta situación, Turnbull (1974,1976) (véase Jaramillo y Lopera (2017)) desarrolló un algoritmo iterativo que realiza una estimación no paramétrica de máxima verosimilitud de la función de supervivencia de los datos. Sea una muestra de intervalos observados independientes [Li, Ri], i = 1, . . . , n en los cuales se producen los tiempos de fallos T1, . . . , Tn(variables aleatorias independientes). En estas condiciones, la función de verosimilitud que debe ser maximizada es: L(S) = n Y i=1 [S(R+ i)−S(L− i)] (2.1) donde S(R+ i)es el valor de la función de supervivencia justo después del extremo derecho RiyS(L− i)es el valor de la función de supervivencia justo antes del extremo izquierdo Lide cada uno de los iintervalos observados. Para la resolución del problema (2.1) es necesario definir un conjunto Aque contenga los extremos izquierdos y otro conjunto Bcon los extremos derechos de la muestra inicial de intervalos observados: A={Li, i = 1, . . . , n} B={Ri, i = 1, . . . , n} Una vez definidos los conjuntos anteriores, se pueden considerar nuevos intervalos de tiempo [aj, bj], j = 1, . . . , m con m≤n,aj∈Aybj∈Bde tal forma que una función que maximice (2.1) es constante entre los intervalos [aj, bj]e indefinida dentro de ellos, lo CAPÍTULO 2. CONCEPTOS BÁSICOS 9 2.1. FUNCIÓN DE SUPERVIVENCIA que implica que ˆ P[T∈(bj−1, aj)] = 0 para cualquier j. Denótense los incrementos de S dentro de los intervalos [aj, bj]por sj, j = 1, . . . , m. Por tanto, (2.1) debe ser maximizada como función de s1, . . . , smsujeto a que:        sj≥0 sm= 1 − m−1 X j=1 sj Para resolver el problema de optimización, Turnbull propone el uso del algoritmo de auto-consistencia. Sean αij =I{[aj,bj]∈[Li,Ri]}, i = 1, . . . , n, j = 1, . . . , m las variables indicadoras cuyo objetivo es informar de si el intervalo [aj, bj]se encuentra dentro del intervalo [Li, Ri]o no. Entonces, la probabilidad de que Tiesté dentro del intervalo [aj, bj] dado un vector s= (s1, s2, . . . , sm)′viene dado por: µij(s) = αijsj m X k=1 αiksk (2.2) debido a que la estimación de Ses constante fuera de los intervalos [aj, bj]. Además, la proporción de observaciones en [aj, bj]es igual a: πj(s) = 1 n n X i=1 µij(s)(2.3) y se dice que un vector s= (s1, . . . , sm)′es auto-consistente, si sj=πj(s), j = 1, . . . , m. Una vez definidas las herramientas necesarias, el algoritmo de auto-consistencia de Turnbull para calcular el estimador no paramétrico de S(t)se puede implementar siguiendo los siguientes pasos: 1. Obtención de las estimaciones iniciales de s; por ejemplo, s(0) j=1 m, j = 1, . . . , m. 2. Para i= 1, . . . , n,j= 1, . . . , m, hallar µij(s(0))conforme a (2.2), y luego πj(s(0)) según (2.3). 3. Obtención de estimaciones mejoradas para scalculando s(1) j=πj(s(0)). 4. Volver al paso 2, reemplazando s(0) por s(1) y continuar hasta lograr la convergencia. CAPÍTULO 2. CONCEPTOS BÁSICOS 10 2.1. FUNCIÓN DE SUPERVIVENCIA 0 10 20 30 40 50 0.0 0.2 0.4 0.6 0.8 1.0 t S(t) estimada Fuente propia Figura 2.4: Estimación de la función de supervivencia mediante el estimador de Turnbull de los datos simul del paquete gte (Dehghan et al. (2015)) realizada con R 2.1.4. Estimador de Fleming-Harrington El estimador de Fleming-Harrington (véase Fleming y Harrington (1984)) se trata de un estimador no paramétrico de la función de supervivencia y puede obtenerse a partir del estimador de la función de riesgo acumulado o tasa de fallo acumulado de Nelson y Aalen (véanse Ramírez-Montoya (2015), Nelson (1969)yAalen y Johansen (1978)), el cual se puede calcular de la siguiente manera: ˆ H(t) = n X i=1 di ri representando dial número de muertes (o de ocurrencia del evento de interés) en el momento tiyrial número de individuos en riesgo en el instante ti. El estimador también se puede calcular bajo el enfoque de procesos de conteo (véase la Sección 2.4) y toma la siguiente forma: ˆ H(t) = n X i=1 Zt 0 dNi(s) r(s)ds Este estimador suele ser útil cuando los tamaños de la muestra en estudio son pequeños. Se deduce así que el estimador de Fleming-Harrington es: ˆ S(t) = e−ˆ H(t) Además, una forma de calcular la varianza de la estimación de la función de supervivencia ˆ S(t)sería: d V ar(ˆ S(t)) = ˆ S2(t) n X i=1 di r2 i Asimismo, se podría calcular el error cuadrático medio del estimador en estudio: ECM(ˆ S(ti)) = E[ˆ S(ti)−S(ti)]2= [ˆ S(ti)−S(ti)]2+V ar[ˆ S(ti)] CAPÍTULO 2. CONCEPTOS BÁSICOS 11 2.2. TASAS DE FALLO 0 50 100 150 200 0.0 0.2 0.4 0.6 0.8 1.0 t S(t) estimada Fuente propia Figura 2.5: Estimación de la función de supervivencia mediante el estimador de Fleming-Harrington de los datos NursingHome (Pešta (2021)) realizada con R 2.2. Tasas de fallo En análisis de supervivencia se conoce como tasa de fallo o función de riesgo a aquella función que tiene la capacidad de medir la probabilidad de que a un individuo le ocurra el suceso o evento en estudio a lo largo de un tiempo predeterminado. 2.2.1. Tasa de fallo instantánea o función de riesgo La función de riesgo h(t)(también conocida como tasa de fallo instantánea o hazard function) se define como la tasa de riesgo infinitesimal en T=tcondicionada a que el individuo ha sobrevivido al instante t: h(t) = l´ım ∆t→0 P[t < T ≤t+ ∆t|T >t] ∆t= l´ım ∆t→0 F(t+ ∆t)−F(t) ∆t[1 −F(t)] =f(t) S(t) Equivalentemente, la expresión de dicha función es: h(t) = −1 S(t) d dtS(t) = −d dt ln[S(t)] Se trata de una medida del potencial de fallo en el instante t, depende de la distribución de la variable Ty puede tomar cualquier valor no negativo. 2.2.2. Tasa de fallo acumulado La tasa de fallo acumulado H(t)se define como la acumulación de la tasa de fallo instantánea hasta el instante t. Por tanto, para su cálculo basta considerar la integral de CAPÍTULO 2. CONCEPTOS BÁSICOS 12 2.3. TIPO DE DATOS CENSURADOS la función de riesgo tal y como se indica a continuación: H(t) = Zt 0h(x)dx Además, la tasa de fallo acumulado se relaciona directamente con la función de supervivencia mediante las siguientes expresiones: H(t) = −ln[S(t)] ⇒S(t) = e−H(t) Para obtener un estimador de la tasa de fallo acumulado se puede trabajar de dos maneras distintas: 1. Estimación usando el método de máxima verosimilitud: este método consiste en tomar el menos logaritmo de la función de supervivencia. ˆ H(t) = −log( ˆ S(t)) 2. Estimación de Nelson-Aalen: consultar la Sección 2.1.4. 0 2 4 6 8 10 12 14 0.0 0.5 1.0 1.5 2.0 2.5 t h(t) Fuente propia Figura 2.6: Ejemplo de funciones de riesgo 2.3. Tipo de datos censurados Cuando se realiza un análisis de supervivencia a un conjunto de individuos, puede ocurrir que el suceso de interés no se observe durante el tiempo de estudio o de supervivencia. Esto podría suceder por varios motivos entre los que se encuentra que el suceso se haya producido en un instante anterior al momento en el que comienza el estudio, que el suceso ocurra después de que termine el estudio o simplemente porque un individuo decide o se ve obligado a abandonar el estudio. En todos esos casos, la información que aporta el individuo al estudio es parcial y dichos datos se conocen como datos censurados (censored data). En la práctica existen varios modelos de censura y a continuación, se presentan algunos de ellos (véanse Gijbels (2010), Fahrmeir y Tutz (2013)yCox y Oakes (1984)). CAPÍTULO 2. CONCEPTOS BÁSICOS 13 2.3. TIPO DE DATOS CENSURADOS 2.3.1. Censura de tipo I Este tipo de censura se produce cuando el experimento comienza en un momento determinado t= 0 para un conjunto de individuos concreto y se registra el tiempo de fallo de cada uno de ellos. El estudio termina en un tiempo prefijado t=tc. Si Tfuera el tiempo de fallo y nel número de individuos en estudio entonces las realizaciones de la variable observada serían: Yi=   Tisi Ti≤tc tcsi Ti> tc Normalmente, se considera una variable cuya función es indicar si el tiempo es observado o censurado: δi=   1si Ti≤tc(observado) 0si Ti> tc(censurado) En este caso, el tiempo de censura tces fijo y el número de tiempos de fallo observados es una variable aleatoria que toma valores en el conjunto {0,1, . . . , n}. 2.3.2. Censura de tipo II La censura de tipo II se produce cuando el experimento se interrumpe cuando ha fallado una fracción prefijada r/n de individuos en estudio, siendo run número entero tal que 1≤r≤n. Esta situación se suele dar por ejemplo en pruebas de fiabilidad cuando, para un determinado lote de productos o piezas, el experimento se detiene en cuanto ha fallado una fracción determinada de las nunidades. Si los tiempos de fallos son T1, . . . , Tn y se considera la muestra ordenada T(1) ≤ ··· ≤ T(n), entonces las nrealizaciones de la variable observada serán: Yi=   T(i)si i≤r T(r)si i>r Existen así n−robservaciones de las que solo se conoce que su tiempo de fallo es mayor que T(r)y que se consideran observaciones censuradas. De esta forma, el número de tiempos observados res fijo, mientras que el punto final Tres aleatorio. 2.3.3. Censura aleatoria 2.3.3.1. Censura a la derecha Este tipo de censura suele darse en estudios de medicina y se caracteriza por la existencia de individuos que no experimentan el suceso de interés antes de que acabe el tiempo de seguimiento prefijado. Considérese un estudio clínico con npacientes (i= 1, . . . , n), cuyo objetivo es el estudio de la eficacia de un nuevo tratamiento médico y supóngase que el seguimiento de los pacientes será durante un período de tiempo fijado de antemano. Cada paciente recibió el tratamiento y sea la cura del paciente el evento de interés del estudio. Denótese: CAPÍTULO 2. CONCEPTOS BÁSICOS 14 2.3. TIPO DE DATOS CENSURADOS Ei: momento en el que el i-ésimo paciente entra en el estudio Di: momento en el que el i-ésimo paciente se cura (evento de interés) Siguiendo dicha notación, se tendrá entonces que el tiempo de supervivencia del i-ésimo individuo será: Ti=Di−Ei Sin embargo, existen situaciones en las que los tiempos de supervivencia de todos los individuos no pueden ser observados y la información que se puede recabar de dichos pacientes es parcial. Esta censura se puede deber a alguna de las siguientes causas: Fin del estudio: termina el ensayo y no se ha observado el evento de interés. Abandono del individuo: el individuo deja de participar en el estudio. Salida forzosa: los efectos del tratamiento obligan al paciente a abandonar el estudio. En los individuos en los que se produce alguna de las situaciones anteriores se observa que Ci=Ti−Ei, siendo Ti> CiyCiel tiempo de censura. Se dice así que el tiempo de supervivencia Tiestá censurado por el tiempo de censura Ci. Una modelización estadística del problema contará con la participación de dos variables aleatorias: T≡tiempo de supervivencia con función de distribución Fy función de supervivencia SfyC≡tiempo de censura con función de distribución Gy función de supervivencia Sg. De esta forma, en cada individuo se observa la dupla (Yi, δi)siendo Yi= m´ın{Ti, Ci};δi=   1si Ti≤Ci(observado) 0si Ti> Ci(censurado) Normalmente se suele suponer independencia entre las variables TiyCi. No obstante, esta suposición es violada cuando el individuo debe salir forzosamente del estudio. 2.3.3.2. Censura a la izquierda La censura a la izquierda se produce cuando solo se conoce que el tiempo de supervivencia de algún individuo es anterior a un momento determinado. Por ejemplo, este tipo de censura puede darse cuando el tiempo final del estudio ocurre antes del momento que se considera como origen del tiempo de estudio. Considérese un estudio pediátrico en el que se desea estudiar la edad en la que los bebés son capaces de dar sus primeros pasos. En este caso, es claro que el suceso de interés es empezar a caminar de manera independiente y la muestra de los bebés puede ser tomada de una guardería cualquiera. Un dato será censurado a la izquierda si al inicio del estudio el niño ya es capaz de caminar solo. Estas observaciones se caracterizan por la siguiente propiedad Ti≤Cidonde Ciserá la edad del bebé al inicio del estudio. En un problema con censura a la izquierda interviene un tiempo de censura aleatorio Ciasociado a cada tiempo de supervivencia Tiy cada una de las observaciones viene determinada por el par (Yi, δi)donde: Yi= m´ax{Ti, Ci};δi=   1si Ti≥Ci(observado) 0si Ti< Ci(censurado) CAPÍTULO 2. CONCEPTOS BÁSICOS 15 2.4. TEORÍA DE PROCESO DE CONTEO (COUNTING PROCESS THEORY ) Al igual que ocurría con la censura aleatoria a la derecha, una suposición que se suele hacer en este tipo de censura es considerar que las variables TiyCison independientes. 2.3.3.3. Datos doblemente censurados Los datos doblemente censurados se producen cuando existen a la vez tiempos de censura a la derecha y a la izquierda asociados a cada uno de los tiempos de supervivencia Ti. Cabe señalar que frecuentemente en estudios donde existe censura por la izquierda, también existe censura por la derecha. En estos casos, las observaciones vienen determinadas por la dupla (Yi, δi)donde: Yi= m´ın{m´ax(Ti, Cl,i), Cr,i};δi=       1si Cl,i ≤Ti≤Cr,i (observado) 0si Cr,i < Ti(censurado por la derecha) −1si Ti< Cl,i (censurado por la izquierda) En otras palabras, el tiempo de supervivencia real solo se observa cuando supera el tiempo de censura izquierdo Cl,i (tiempo antes del cual algunos individuos presentan el evento de interés) y a la misma vez no supera el tiempo de censura derecho Cr,i (tiempo después del cual algunos individuos presentan el evento de interés). Obsérvese que fijando Cl=−∞ yCr= +∞, el caso de la censura doblemente aleatoria incluye como casos especiales la censura aleatoria a la derecha y la censura aleatoria a la izquierda, respectivamente. 2.3.4. Datos censurados por intervalo La censura por intervalo o censura intermitente se produce cuando el individuo que está participando en el ensayo experimenta el punto final del estudio dentro de un intervalo determinado. Este tipo de censura suele darse cuando los individuos en estudio no son monitorizados de manera continua sino que solo se observan en momentos discretos de tiempo. Además, la censura a la izquierda es un caso especial de este tipo de censura, debido a que el tiempo final se producirá en el intervalo posterior a la entrada del individuo en el estudio y antes de que la investigación llegue a su fin. Considérese el individuo i-ésimo del estudio que no ha presentado el evento de interés antes del tiempo Lia pesar de que en el siguiente tiempo de observación Riel individuo ya presentó el evento en estudio. En ese caso, se dice que Ties censurado en el intervalo (Li, Ri], i.e., Li< Ti≤Ri. 2.4. Teoría de proceso de conteo (counting process theory) Tal y como se enuncia en Andersen et al. (1985), un proceso de conteo multivariante N={N1(t), N2(t), . . . , Nk(t)}, t ≥0es un proceso estocástico con kcomponentes que puede considerarse como un recuento de sucesos a lo largo del tiempo tsiendo klos sucesos posibles que pueden ocurrir. De esta forma, Nh(t)representará el número de sucesos de tipo hque se producen en el intervalo de tiempo [0, t]. Se supondrá que cada componente del proceso Nh(t)tiene saltos de tamaño 1. Además, no se permite que se produzan CAPÍTULO 2. CONCEPTOS BÁSICOS 16 3.1. TÉCNICAS GENERALES sea la naturaleza de los datos) en el instante tiyri, i = 1, . . . , n el número de individuos en riesgo en el instante ti, se tiene que: ZLR g= n X i=1 ω(ti)rig dig rig −di ri!, g = 1, . . . , k, donde dig es el número de fallos o muertes en el grupo o población ge instante ti,rig el número de individuos en riesgo en el grupo o población ge instante tiyω(ti)la función peso. Asimismo, se cumplen las siguientes relaciones: di= k X g=1 dig ri= k X g=1 rig Además, si se considera un vector al que se llamará ZLR que contenga los kestadísticos dados por cada uno de los ZLR gy sea ΣLR la matriz de varianza-covarianza del vector ZLR, entonces el estadístico asociado al test vendrá dado por: QLR = (ZLR)′(ΣLR)−ZLR (3.3) donde (ΣLR)−es la inversa generalizada4de ΣLR y los elementos de la matriz de varianzacovarianza del vector ZLR son: σLR gj = n X i=1 ω(ti)2rig riδgj −rij ri ri−di ri−1!di siendo δgj =   1si g=j 0si g=j yg, j = 1, . . . , k. El estadístico QLR sigue aproxidamente una distribución χ2con k−1grados de libertad. A continuación, se enumeran algunas de las funciones peso más habituales. Sin embargo, el test que se suele utilizar más a menudo es aquel que usa pesos iguales y se conoce como el test log-rank. En la siguiente tabla se pueden consultar algunos de los pesos ω(ti)más utilizados a la hora de aplicar el test log-rank ponderado. Nótese que en dicha tabla se usará la siguiente notación: ˆ S(t) = m Y i=1 1−di ri ˜ S(t) = m Y i=1 1−di ri+ 1 donde mrepresenta el m-ésimo instante de ocurrencia del evento en estudio. 4Sea Xuna matriz n×m, entonces se dice que la matriz X−de dimensión m×nes una inversa generalizada de Xsi se cumple que XX−X=X. CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 23 3.1. TÉCNICAS GENERALES Test Peso Log-rank 1 Gehan-Breslow ri Tarone-Ware √ri Peto-Peto ˜ S(ti) Peto-Peto modificado ˜ S(ti)ri/(ri+ 1) Fleming-Harrington (0,1) 1 −ˆ S(ti−1) Fleming-Harrington (1,0) ˆ S(ti−1) Fleming-Harrington (1,1) (1 −ˆ S(ti−1))ˆ S(ti−1) Fleming-Harrington (0.5,2) (1 −ˆ S(ti−1))2qˆ S(ti−1) Tabla 3.1: Ejemplos de funciones pesos ω(ti) 3.1.2. Prueba de tendencia (trend-testing) En muchos experimentos y estudios cuando se desea comparar la eficacia de varios tratamientos o modalidades, es necesaria la existencia de pruebas sensibles. Un enfoque clásico de análisis es el análisis de la varianza (ANOVA) y las pruebas chi-cuadrado de grupos múltiples. Sin embargo, un enfoque más sensible puede ser el análisis de tendencias, entendiéndose por tendencia cualquier asociación entre el orden del tratamiento y la magnitud de la respuesta. Se suele usar cuando se trabaja con modelos de regresión y cuando se comparan tres o más curvas de supervivencia en las que las columnas están en un orden natural (edades o estadios del cáncer5). Además, en esencia, el test comprueba si existe una tendencia lineal entre el orden de las columnas y la mediana de supervivencia6. Tal y como especifican Kirkwood y Sterne (2010), una prueba de tendencia es una aproximación (basada en una prueba de puntuación7) a un test de razón de verosimilitud en el que la hipótesis nula es el coeficiente de regresión para un efecto lineal igual a cero. Este tipo de análisis es adecuado para datos binarios mediante la prueba χ2, aunque se pueden obtener resultados similares a partir de un modelo de regresión logística. Sin embargo, en general el test de tendencia suele ser más sensible que las pruebas estándar. Cuando se trabaja con dos tratamientos, la prueba de tendencia y la prueba estándar obtienen resultados idénticos ya que en cada caso, el estadístico correspondiente tiene los mismos grados de libertad. Sin embargo, con tres o más tratamientos, los grados de 5La determinación del estadio de un cáncer ayuda a definir la ubicación del cáncer, si se ha diseminado o hacia dónde, y si está afectando a otras zonas del cuerpo. 6Concepto que se suele usar en medicina y que denota el tiempo que pasa desde la fecha del diagnóstico o el comienzo del tratamiento de una enfermedad, como el cáncer, durante el cual la mitad de los pacientes con la enfermedad siguen vivos. En un ensayo clínico, se mide la mediana de supervivencia como una forma de determinar la eficacia de un tratamiento nuevo. 7La prueba de puntuación, score test o prueba del multiplicador de Lagrange es un test utilizado para comprobar si una restricción impuesta a un modelo estimado por máxima verosimilitud es violada por los datos. Dicha prueba evalúa las restricciones de los parámetros estadísticos basándose en el gradiente de la función de verosimilitud evaluada en el valor del parámetro hipotético bajo la hipótesis nula. Además, la distribución asintótica bajo la hipótesis nula es una distribución χ2. CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 24 3.1. TÉCNICAS GENERALES libertad del estadístico asociado al análisis estándar aumentan rápidamente, mientras que con el análisis de tendencias no lo hacen, dando lugar a p-valores más pequeños. No obstante, se deben tener en cuenta las limitaciones que presenta el análisis de tendencias. En primer lugar, si no existe tendencia en los datos, el método de análisis estándar puede ser más sensible que la prueba de tendencia. Como consecuencia, se deben realizar pruebas estándar además de las pruebas de tendencia. En segundo lugar, las pruebas de tendencia suponen una tendencia lineal en la respuesta de los datos con los tratamientos posteriores. Esto significa que, con datos continuos, las medias de los grupos de tratamiento aumentan linealmente, y con datos binarios, las probabilidades de respuesta aumentan exponencialmente, i.e., los logaritmos de las probabilidades aumentan linealmente. El efecto lineal es una suposición simplificadora que debe comprobarse. Dicha comprobación es fácil cuando se tienen como máximo tres categorías. De esta forma, se puede concluir que las pruebas de tendencia pueden ser más sensibles que los métodos estándar de comparación de tratamientos. Las pruebas de tendencia están diseñadas para detectar diferencias ordenadas en las curvas de supervivencia, i.e., para al menos un grupo: S1(t)≥S2(t)≥S3(t)≥ ··· ≥ Sk(t)t≤τ con al menos una desigualdad estricta y donde τes la mayor t∈Ten la que todos los grupos involucrados en el estudio tienen al menos un sujeto en riesgo y k= 1, . . . , M con Mel número de poblaciones en estudio. Además, la hipótesis nula del test es: S1(t) = S2(t) = ··· =Sk(t)t≤τ Según Dardis (2018) las puntuaciones utilizadas para construir la prueba suelen ser s= 1,2, . . . , k, pero pueden darse como un vector sque represente una característica numérica del grupo correspondiente. Además, la prueba se basa en los siguientes cálculos: ZTT g(ti) = X ti≤τ ω(ti)heig −rig ei rii, g = 1,2, . . . , k, donde ω(ti)es la función peso (véase la tabla de pesos 3.1), ries el número de individuos en riesgo en el instante ti,rig es el número de individuos en riesgo en el grupo o población ge instante ti,eies el número de eventos (por tiempo) y eig es el número de eventos (por tiempo) en el grupo o población g. De esta forma, el estadístico asociado al test viene definido por: ZTT = k X g=1 sgZTT g(τ) v u u t k X g=1 k X j=1 sgsjσTT gj (3.4) donde σTT gj es el elemento apropiado de la matriz de varianza-covarianza asociada a dicho estadístico y sgysjson los elementos g-ésimo y j-ésimo del vector srespectivamente. Además, cuando la hipótesis nula del contraste se rechaza, se cumple que el estadístico ZTT es mayor que el α-ésimo percentil de una distribución normal. CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 25 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS A continuación, se especifican los elementos de la matriz de varianza-covarianza de ZT T para al menos k≥2poblaciones: σTT gj =               X ti≤τ rig (ri−rig)ei(ri−ei) r2 i(ri−1) g, j = 1, . . . , k yg=j, −X ti≤τ rig rij ei(ri−ei) r2 i(ri−1) g, j = 1, . . . , k con g=j. 3.2. Técnicas específicas para datos censurados Tal y como comentan Calle et al. (2009), el comportamiento asintótico de los estadísticos asociados a los test basados en datos censurados por intervalos es más difícil de estudiar que en el caso de la censura a la derecha porque por ejemplo, la teoría del proceso de recuento no puede ser aplicada. Esta dificultad se afronta con métodos de remuestreo o con métodos de inferencia estándar que requieren supuestos adicionales. Entre los métodos de remuestreo se encuentran los procedimientos de permutación, imputación múltiple y procedimientos bootstrap. A continuación, se enumerarán diferentes técnicas para la comparación de curvas de supervivencia de dos o más poblaciones diferentes con la presencia de datos censurados a la derecha (RC ≡right censored data) o por intervalo (IC ≡interval censoring o interval-censored data) (véase Oller y Langohr (2017)). Para cada caso, se describirán los estadísticos de prueba de la clase de Fleming-Harrington Gρ,λ, concepto que introdujeron Fleming y Harrington (1982) y que se trata de un subconjunto de lo que se conoce como la prueba log-rank ponderada (Sección 3.1.1). Observación: No se tratarán técnicas para datos censurados a la izquierda debido a que dichos datos pueden transformarse fácilmente en datos censurados a la derecha. En el libro de Goel y Klein (1992) se determina que para obtener datos censurados a la derecha partiendo de una muestra de datos con censura a la izquierda basta con multiplicar cada uno de los datos por −1. De esta forma, pueden extrapolarse las técnicas que se detallan a continuación a los datos censurados a la izquierda. ▲ Cuando se trata con muestras de dos o más poblaciones diferentes, puede ser de interés comparar las curvas de supervivencia asociadas a cada una de las poblaciones a través de los siguientes contrastes: Test con 2 poblaciones (bilateral): H0:S1(t) = S2(t)∀t∈T H1:S1(t)=S2(t)para algún t∈T(3.5) También puede contemplarse el contraste en el que la hipótesis alternativa es unilateral. Test con kpoblaciones: H0:S1(t) = ··· =Sk(t)∀t∈T H1:∃i, j, t tal que Si(t)=Sj(t)(3.6) CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 26 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS Prueba de tendencia con hipótesis alternativa creciente: H0:S1(t) = ··· =Sk(t)∀t∈T H1:S1(t)≤ ··· ≤ Sk(t)∀t∈T(3.7) con al menos una desigualdad estricta. También puede contemplarse la hipótesis alternativa decreciente H1:S1(t)≥ ··· ≥ Sk(t)∀t∈Tcon al menos una desigualdad estricta. Bajo condiciones generales, el estadístico básico de cada uno de los contrastes anteriores es: Z= (Z1, . . . , Zk)′ el cual es asintóticamente normal de media cero y matriz de varianza-covarianza Σ. La forma explícita de ZyΣse indica a continuación y depende del tipo de censura y de la aproximación asintótica. El estadístico de prueba asociado al contraste (3.5) es: V1=Z2 √σ22 (3.8) donde σ22 es el segundo elemento de la diagonal de Σy la distribución asintótica de V1es la distribución normal estándar. Por otro lado, el estadístico asociado a la hipótesis del test con kpoblaciones (3.6) es: V2=Z′Σ−1Z(3.9) donde Σ−1es la inversa generalizada de Σy la distribución asintótica de V2es una χ2con k−1grados de libertad. Por último, el test de tendencia (3.7) con la hipótesis alternativa creciente puede resolverse con el estadístico: V3=a′Z √a′Σa(3.10) donde a= (a1, . . . , ak)′,a1<··· < ak, son los valores de las covariables de cada una de las poblaciones en estudio. La distribución asintótica de V3es una distribución normal estándar. Un test unilateral por la izquierda puede mostrar evidencias de una tendencia creciente y uno por la derecha, de una tendencia decreciente. 3.2.1. Datos censurados a la derecha 3.2.1.1. Enfoque mediante el proceso de conteo Considérese una muestra aleatoria de datos censurados a la derecha [Li, Ri], i = 1, . . . , n con Li=Ri(datos no censurados) o Ri= +∞(datos censurados a la derecha). Sean t1<··· < tmlos mtiempos no censurados y supóngase que existen n−m tiempos censurados a la derecha. En estas condiciones, las componentes del estadístico del test log-rank ponderado ZRC son de la forma: ZRC g= m X q=1 ω(tq)(dqg −rqg rq dq), g = 1, . . . , k, (3.11) CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 27 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS donde dqrepresenta al número de muertes en el instante tq,rqel número de individuos en riesgo en el instante tq, es decir, el número de individuos vivos justo antes del instante tqyω(tq)es el peso sugerido por Fleming-Harrington. Tanto dqg como rqg corresponden a las mismas definiciones pero referidas a las observaciones de la población g-ésima. Las ponderaciones o pesos sugeridos por Fleming-Harrington son: ω(tq) = ˆ S(t− q)ρ(1 −ˆ S(t− q))λ,(3.12) donde ˆ S(t− q)es el valor del estimador de Kaplan-Meier de la función de supervivencia justo antes del instante tq. Cuando (ρ, λ) = (0,0) y(ρ, λ) = (1,0), el estadístico ZRC corresponde a los estadísticos de prueba log-rank y Prentice-Wilcoxon, respectivamente. La elección de los pesos debe hacerse antes del examen de los datos teniendo en cuenta que deben proporcionar la mayor potencia estadística, que a su vez depende de cómo se cree que se viola la nulidad. No obstante, en los ensayos clínicos, el poder estadístico puede no ser la única consideración para determinar el esquema de ponderación, ya que también deben contemplarse los aspectos clínicos. Asintóticamente, la distribución de (3.11) se obtiene mediante el teorema central del límite de Martingala. Además, bajo la hipótesis nula, ZRC sigue asintóticamente una distribución normal con media nula y matriz de varianza-covarianza ΣRCCount.. Los elementos de dicha matriz vienen dados por: σRCCount. gj =                    m X q=1 ω(tq)2rqg (rq−rqg) r2 q (rq−dq)dq rq−1g, j = 1, . . . , k yg=j, − m X q=1 ω(tq)2rqg rqj r2 q (rq−dq)dq rq−1g, j = 1, . . . , k con g=j. 3.2.1.2. Enfoque de permutación Según Oller y Langohr (2017), el estadístico (3.11) puede escribirse de forma lineal tal y como se indica a continuación: ZRCP erm. = n X i=1 zici(3.13) donde zi= (Z1i, . . . , ZKi)′es un vector de covariables indicadoras del grupo o población al que pertenece cada observación, i.e., Zgi =   1si la i-ésima observación está en la g-ésima población 0en caso contrario (3.14) ycies un valor de puntuación asociado a cada una de las observaciones y toma la siguiente forma: ci=           ω(tq)− q X s=1 ω(ts)ds rs si Li=Ri=tq, − q X s=1 ω(ts)ds rs si Ri= +∞ytq≤Li< tq+1. (3.15) CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 28 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS El enfoque de permutación puede aplicarse fácilmente a esta forma lineal. Si la hipótesis nula es cierta y el proceso de censura subyacente es idéntico en todas las poblaciones, entonces las etiquetas de las puntuaciones cison intercambiables. La distribución de permutación de ZRCP erm. se obtiene entonces permutando las etiquetas y volviendo a calcular el estadístico de prueba para todas las posibles permutaciones. Además, los p-valores de la distribución de permutación pueden calcularse utilizando la enumeración completa, un algoritmo de red, un método de aproximación como el de Montecarlo o asintóticamente utilizando una versión del teorema del límite central para variables aleatorias intercambiables. Nótese que este último método produce una aproximación a una distribución normal de media 0 y matriz de varianza-covarianza definida por: ΣRCP erm. =1 n−1 n X i=1 c2 i! n X i=1 ziz′ i−¯zi¯z′ i! donde ¯z =1 n n X i=1 zies la media muestral. 3.2.2. Datos censurados por intervalos 3.2.2.1. Enfoque de permutación Considérese una muestra aleatoria de datos censurados por intervalos (Li, Ri], dados t1<··· < tmdenótense los elementos ordenados {Li, Ri;i= 1, . . . , n}y sea ˆ S(t)la estimación no paramétrica de la función de supervivencia a través del estimador de Turnbull. En dichas condiciones, las componentes del estadístico ZIC son de la forma ZIC g= m X q=1 ωIC(tq)dIC qg −rIC qg rIC q dIC q, g = 1, . . . , k, (3.16) donde dIC q=n(ˆ S(t− q)−ˆ S(tq)) representa el número estimado de muertes en el instante tqyrIC q=nˆ S(t− q)expresa el número estimado de individuos en riesgo en el instante tq. Además, dIC qg yrIC qg se definen como: dIC qg =n 1 Ng N X j=1 Zgj ˆ Sj(t− q)−1 Ng N X j=1 Zgj ˆ Sj(tq)! rIC qg =n1 Ng N X j=1 Zgj ˆ Sj(t− q) siendo Ngel tamaño de la muestra de la g-ésima población, Zgj es la variable indicadora definida en (3.14), ˆ Sjla estimación de la función de supervivencia del j-ésimo individuo usando el algoritmo EM y definida como la estimación de la función de supervivencia global de Turnbull ˆ S(t)truncada en el j-ésimo intervalo observado (véase Gómez y Oller (2008,2012), i.e., ˆ Sj(t) = Pˆ S[T∈(t, +∞)|(Li, Ri]] = ˆ S(m´ax{Li, t})−ˆ S(m´ax{Ri, t}) ˆ S(Li)−ˆ S(Ri) CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 29 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS siendo Pˆ Sla probabilidad dada por la función de supervivencia ˆ S(t)y N= k X g=1 Ng. Observación: La estimación de ˆ Sj(t)obtenida en el paso de maximización del algoritmo EM posee las siguientes propiedades: 1. ˆ S(t) = 1 n n X j=1 ˆ Sj(t) 2. ˆ Sj(t) = I{t<tj}cuando el j-ésimo individuo no está censurado y T=tj. ▲ Por último, los pesos o ponderaciones propuestos por Oller y Langohr son: ωIC(tq) = ˆ S(t− q)B(1 −ˆ S(tq); λ+ 1, ρ)−B(1 −ˆ S(t− q); λ+ 1, ρ) ˆ S(t− q)−ˆ S(tq),(3.17) donde B(x;p, r) = Zx 0tp−1(1 −t)r−1dt;p, r ∈Ctales que Re(p)>0yRe(r)>0 denota a la función beta incompleta8. Cabe destacar que estas ponderaciones derivan del enfoque del vector de puntuación que se verá en la Sección 3.2.2.2. Cuando la dupla (ρ, λ) toma los valores (0,0), el estadístico ZIC se reduce al estadístico log-rank mientras que si toma los valores (1,0), el estadístico se reduce al estadístico de Wilcoxon-Peto. Además, como ˆ Ses asintóticamente continua en el instante tq, la función peso ωIC(tq)se parece a[ˆ S(t− q)]ρ[1 −ˆ S(t− q)]λy como consecuencia, la interpretación de los parámetros λyρ reproducen la interpretación que posee la familia de datos censurados a la derecha. Esta clase de estadísticos ZIC puede escribirse en forma lineal tal y como se especifica a continuación: ZICPerm. = n X i=1 zIC icIC i(3.18) donde zIC i= (ZIC 1i, . . . , ZIC Ki )′es un vector de covariables indicadoras del grupo o población al que pertenece cada observación, i.e., ZIC gi =   1si la i-ésima observación está en la g-ésima población 0en caso contrario (3.19) ycIC ies un valor de puntuación asociado a cada una de las observaciones: cIC i=       ˆ S(Ri)B(1−ˆ S(Ri); λ+1, ρ)−ˆ S(Li)B(1−ˆ S(Li); λ+1, ρ) ˆ S(Li)−ˆ S(Ri)si ˆ S(Ri)= 0, −B(1 −ˆ S(Li); λ+ 1, ρ)si ˆ S(Ri)=0. 8La función beta incompleta es una generalización de la función beta. Además, coincide con la función beta cuando x= 1. Se define como B(x;a, b) = Zx 0 ta−1(1 −t)b−1dt con a, b ∈Ctales que Re(a)>0y Re(b)>0. CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 30 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS Nótese que la equivalencia entre la forma log-rank ponderada (3.16) y la forma lineal (3.18) se basa en las ponderaciones (3.17). Además, en el hipotético caso de que el tipo de censura sea el mismo en cada una de las poblaciones en estudio, la forma lineal ZICP erm. dada en (3.18) permite el uso de una distribución de permutación con media 0 y matriz de varianza-covarianza definida de la siguiente forma: ΣICPerm. =1 n−1 n X i=1 (cIC i)2! n X i=1 zIC i(zIC i)′−¯zIC i(¯zIC i)′! donde ¯zIC =1 n n X i=1 zIC ies la media muestral. Observación: Para la construcción del estadístico ZICP erm. se ha considerado que los intervalos observados son semiabiertos aunque con pequeñas modificaciones, se puede hacer extensible a intervalos cerrados o abiertos. Además, dicha metodología puede aplicarse a los datos censurados a la derecha como un caso particular de censura por intervalos. En ese caso, los resultados obtenidos por las formas lineales ZRCP erm. dada por (3.13) y ZICPerm. dada por (3.18) son parecidas pero no idénticas. ▲ 3.2.2.2. Enfoque basado en un vector de puntuación Cuando se restringe la clase de estadísticos ZIC de Fleming-Harrington a la subclase λ= 0 (suele denominarse como familia Gρ), cada estadístico del test dado por el valor de ρ puede considerarse como una prueba de puntuación eficiente asociado a un modelo lineal. Para tratar este tipo de enfoque es necesario que los datos censurados se muestren como datos discretos o agrupados en intervalos de tiempo continuo, i.e., L, R ∈ {t0, t1, . . . , tm} donde 0 = t0< t1<··· < tm< tm+1 = +∞y además, se asume un modelo lineal de la forma: g(Ti) = −(zIC i)′β+εi siendo guna función creciente desconocida, βun vector de constantes desconocidas y εi una variable que tiene asociada como función de supervivencia lo siguiente: Sε(t) = 1 + ρ et−1 ρ En tales condiciones, se tiene que el estadístico asociado al contraste es: ZICV ec. ="∂log(L(β,θ)) ∂β#β=0,θ=ˆ θ0 (3.20) donde L(β,θ)es la función de verosimilitud y θ= (θq)m q=1 son parámetros molestos o de ruido9tales que θq=g(tq). Considerando datos discretos censurados por intervalos y suponiendo cierta la hipótesis nula del contraste de igualdad de curvas de supervivencia, el comportamiento asintótico del vector de puntuación ZICV ec. procede de la teoría de máxima verosimilitud. Usando 9Los parámetros molestos o de ruido (nuisance parameters) son parámetros de un modelo estadístico que no tienen interés científico alguno en el estudio realizado, pero cuyos valores son necesarios para hacer inferencia sobre otros parámetros de interés. CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 31 3.2. TÉCNICAS ESPECÍFICAS PARA DATOS CENSURADOS este enfoque es necesario que ninguno de los valores estimados del vector θse aproxime a la frontera del espacio que contiene a dichos parámetros, en otras palabras, al espacio definido por 1>ˆ S(t1)>··· >ˆ S(tm)>0. En dicho caso, la distribución del estadístico es asintóticamente normal de media 0 y matriz de varianza-covarianza ΣICV ec. especificada por Gómez y Oller (2008) y que se muestra a continuación: ΣICV ec. =−"∂2log(L(S)) ∂β∂β′− ∂2log(L(S)) ∂β∂θ′! ∂2log(L(S)) ∂θ∂θ′!−1 ∂2log(L(S)) ∂β∂θ′!#β=0,θ=ˆ θ0 donde L(S) = n Y i=1{S(Li)−S(Ri)} con S(Li)el valor de la función de supervivencia en el extremo izquierdo Li,S(Ri)el valor de la función de supervivencia en el extremo derecho Riy para cualesquiera parámetros ϕxeϕyarbitrarios se cumple: ∂2log(L(S)) ∂ϕx∂ϕy = n X i=1 1 ∆i(S)" ∂2∆i(S) ∂ϕx∂ϕy!−1 ∆i(S) ∂∆i(S) ∂ϕx! ∂∆i(S) ∂ϕy!# con ∆i(S) = S(Ri|zIC i)−S(Li|zIC i) donde zIC ies el vector de covariables indicadoras de la población a la que pertenece cada observación y cuyos elementos están definidos en (3.19). Además, se cumple que: "∂S(tj|zIC i) ∂βx#β=0,θ=ˆ θ =−1 ρ{ˆ S(tj)}{1−(ˆ S(tj))ρ}ZIC ix "∂S(tj|zIC i) ∂θx#β=0,θ=ˆ θ =−1 ρ{ˆ S(tj)}{1−(ˆ S(tj))ρ}I{j=x} "∂2S(tj|zIC i) ∂βx∂βy#β=0,θ=ˆ θ =1 ρ2{ˆ S(tj)}{1−(ˆ S(tj))ρ}{1−(ρ+ 1)(ˆ S(tj))ρ}ZIC ix ZIC iy "∂2S(tj|zIC i) ∂βx∂θy#β=0,θ=ˆ θ =1 ρ2{ˆ S(tj)}{1−(ˆ S(tj))ρ}{1−(ρ+ 1)(ˆ S(tj))ρ}ZIC ix I{j=y} "∂2S(tj|zIC i) ∂θx∂θy#β=0,θ=ˆ θ =1 ρ2{ˆ S(tj)}{1−(ˆ S(tj))ρ}{1−(ρ+ 1)(ˆ S(tj))ρ}I{j=x=y} Observación: el estadístico ZRC y su correspondiente matriz de varianza-covarianza ΣRCCount. pueden obtenerse a partir del enfoque del vector de puntuación cuando λ= 0. En dicho caso, en vez de utilizar una función de verosimilitud completa, se usa una función de verosimilitud parcial. Debido a esta apreciación, la combinación ZRC yΣRCCount. o la combinación ZICV ec. yΣICV ec. pueden utilizarse para el análisis de datos censurados a la derecha a pesar de que los resultados obtenidos difieran un poco entre ellos. ▲ CAPÍTULO 3. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA 32 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA pacientes abandonaron el seguimiento debido a un traslado o a que sobrevivieron hasta la finalización del estudio. Es destacable el hecho de que en los datos se incluyeron mediciones de covariables muy relacionadas con el análisis de supervivencia. Entre ellas se encuentra la variable Condition que proporciona una medida de la capacidad funcional del paciente en el momento del diagnóstico y la variable Grade que ofrece el grado en que la célula tumoral se asemeja a la célula huésped. En primer lugar, se leerá el fichero que contiene dichos datos y se mostrará la cabecera del objeto oropharynx que será un data.frame que contendrá la información necesaria para realizar posteriormente una estimación de la función de supervivencia: oropharynx <- read.table("oropharynx.txt",header = T) Inst Sex TrtGrp Grade Age Cond Site T N entrydate Status Time 2 2 1 1 51 1 2 3 1 2468 1 631 2 1 2 1 65 1 4 2 3 2968 1 270 2 1 1 2 64 2 1 3 3 3368 1 327 2 1 1 1 73 1 1 4 0 5768 1 243 5 1 2 2 64 1 1 4 3 9568 1 916 4 1 2 1 61 1 2 3 0 10668 0 1823 4 1 1 2 65 1 2 4 3 10768 1 637 4 1 2 3 84 1 4 1 3 12068 1 235 6 1 1 2 54 2 1 3 3 13368 1 255 Tabla 4.5: Primeros datos del fichero oropharynx El significado de las variables en estudio son: Inst: indica la institución a la que pertenece cada uno de los pacientes. Los posibles valores que puede tomar dicha variable son 1,2,3,4,5o6. Sex: variable sexo que puede tomar los valores 1 = hombre o 2 = mujer. TrtGrp: grupo de tratamiento que puede tomar los valores 1 = estándar o 2 = test. Grade: grado en que la célula tumoral se asemeja a la célula huésped y puede tomar los valores 1 = bien diferenciado, 2 = moderadamente diferenciado o 3 = pobremente diferenciado. Age: edad del paciente en el momento del diagnóstico. Condition: medida de la capacidad funcional del paciente en el momento del diagnóstico y puede tomar los valores 1 = ninguna incapacidad, 2 = trabajo restringido, 3 = requiere ayuda con los cuidados personales o 4 = confinado en cama. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 39 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA Site: localización del carcinoma 1 = arco facial, 2 = fosa amigdalina, 3 = pilar posterior, 4 = lengua faríngea o 5 = pared posterior. T:1 = primaria, 2 cm o menos, 2 = primarios, 2 a 4 cm, 3 = primaria, más de 4 cm o 4 = tumor invasivo masivo. N:0 = no metástasis, 1 = un solo nodo de 3 cm o menos, 2 = un solo nodo de 3 cm o más o 3 = múltiples nodos. Entrydate: día del año y año, dddyy. Status: estado de la censura (0 = censurado o 1 = muerte). Time: tiempo en días desde el primer día del diagnóstico. Una vez presentadas cada una de las variables, se procede a estimar la función de supervivencia con respecto a la variable TrtGrp: FH.oropharynx <- survfit(Surv(Time, Status) ~ TrtGrp, type = "fleming-harrington", data = oropharynx, conf.type = "plain") Al igual que en la sección anterior, se adjuntan dos tablas (Tabla 4.6 yTabla 4.7) en las que se pueden consultar las estimaciones obtenidas para cada uno de los tratamientos que se contemplan en el conjunto de datos oropharynx: time n.risk n.event survival std.err lower 95 % CI upper 95 % CI 10 100 0 1.000 0.0000 1.000 1.000 50 100 0 1.000 0.0000 1.000 1.000 100 94 5 0.950 0.0218 0.907 0.993 150 88 6 0.890 0.0313 0.828 0.951 200 82 6 0.830 0.0377 0.756 0.903 250 76 6 0.769 0.0422 0.686 0.852 300 70 6 0.709 0.0455 0.620 0.798 400 58 12 0.588 0.0493 0.492 0.685 500 50 7 0.517 0.0502 0.418 0.615 600 43 7 0.445 0.0500 0.347 0.543 700 38 4 0.404 0.0494 0.307 0.501 800 28 7 0.326 0.0479 0.232 0.420 900 25 0 0.326 0.0479 0.232 0.420 1000 16 3 0.286 0.0473 0.193 0.379 CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 40 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA (continúa) time n.risk n.event survival std.err lower 95 % CI upper 95 % CI 1250 10 2 0.248 0.0481 0.153 0.342 1575 1 2 0.117 0.0692 0.000 0.253 Tabla 4.6: Estimador de Fleming-Harrington para el grupo de pacientes con radioterapia sola (estándar) time n.risk n.event survival std.err lower 95 % CI upper 95 % CI 10 95 0 1.000 0.0000 1.000 1.000 50 91 4 0.958 0.0204 0.918 0.998 100 90 1 0.948 0.0227 0.903 0.992 150 84 6 0.885 0.0326 0.821 0.949 200 76 7 0.812 0.0400 0.733 0.890 250 68 8 0.727 0.0457 0.637 0.816 300 58 9 0.631 0.0496 0.534 0.728 400 46 12 0.502 0.0516 0.401 0.603 500 38 7 0.425 0.0512 0.325 0.526 600 26 9 0.323 0.0490 0.227 0.419 700 21 3 0.285 0.0479 0.191 0.379 800 18 1 0.270 0.0476 0.177 0.363 900 16 1 0.255 0.0472 0.163 0.348 1000 13 1 0.239 0.0470 0.147 0.331 1250 10 0 0.239 0.0470 0.147 0.331 1575 3 0 0.239 0.0470 0.147 0.331 Tabla 4.7: Estimador de Fleming-Harrington para el grupo de pacientes con radioterapia junto a un agente quimioterapéutico (test) Se adjunta también una gráfica en la que se pueden observar las estimaciones de Fleming-Harrington de las funciones de supervivencia asociadas a las poblaciones de pacientes tratadas con radioterapia sola o con radioterapia junto a un agente quimioterapéutico: ggsurvplot( FH.oropharynx, data = oropharynx, xlab = "t (días)", ylab = "S(t) estimada", lab = c(10,10,7), size = 1,legend.title = "Tratamientos", palette = c("darkslategray4","darkorange"), CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 41 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA conf.int = TRUE, legend.labs = c("Estándar","Test"), ggtheme = theme_grey(), caption = "Fuente propia", font.caption = c(10,"italic","black"), font.x = c(12,"bold.italic","black"), font.y = c(12,"bold.italic","black") ) + + +++ ++++ + + ++++ ++ + +++ ++++ ++ + + ++ +++ +++ + + + ++ + ++ + ++ + ++ 0.00 0.25 0.50 0.75 1.00 0 500 1000 1500 2000 t (días) S(t) estimada Tratamientos ++ Estándar Test Fuente propia Figura 4.2: Estimador de Fleming-Harrington para cada uno de los tratamientos estudiados en el conjunto de pacientes disponibles en oropharynx 4.1.1.3. Estimador de la tasa de fallo acumulado Un estimador de la función de riesgo acumulado puede obtenerse a partir de un estimador de la función de supervivencia. Se mostrará un ejemplo con el uso del conjunto de datos rats del paquete survival. Dicho conjunto procede del estudio realizado por Mantel et al. (1977) y en él se recoge información sobre la aplicación de dos tratamientos a 300 ratas. Se partió de 100 camadas de ratas y de manera aleatoria se eligieron 3 ratas de cada una de las camadas, una de las cuales fue tratada con un fármaco. Posteriormente, se observó la incidencia de tumores en cada uno de los individuos de la muestra. Las variables involucradas en la investigación fueron: litter: camada a la que pertenece la rata y toma los valores del 1 al 100. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 42 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA rx: tratamiento administrado (1 = fármaco o 0 = control). time: tiempo hasta que aparece un tumor o el último seguimiento realizado. status: estado del evento (1 = tumor o 0 = censura). sex: sexo del animal (m=macho o f=hembra). En la Tabla 4.8 se pueden consultar algunos de los datos del conjunto rats con el que se trabajará: litter rx time status sex 1 1 101 0 f 2 1 91 0 m 3 1 104 0 f 4 0 104 0 m 5 1 104 0 f 6 1 98 0 m 7 1 77 0 f 8 0 76 0 m 9 0 104 0 f 10 1 104 0 m 11 1 88 1 f Tabla 4.8: Cabecera de los datos rats del paquete survival Una vez expuestos los datos, se procede a la obtención del estimador deseado. En primer lugar, se debe estimar la función de supervivencia asociada a cada una de las poblaciones y luego, basta considerar la suma acumulada del cociente entre las variables n.event (individuos que presentaron el evento) y n.risk (individuos en riesgo): ajuste.rats <- survfit(Surv(time, status) ~ rx, data = rats) Se adjuntan dos tablas en las que se pueden observar el resultado de la estimación para cada uno de los tratamientos en estudio: time n.risk n.event survival hazard 34 99 1 0.990 0.0101 39 98 1 0.980 0.0203 45 97 1 0.970 0.0306 67 89 1 0.959 0.0419 70 86 1 0.948 0.0535 72 85 1 0.937 0.0652 CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 43 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA (continúa) time n.risk n.event survival hazard 73 84 1 0.925 0.0771 80 78 2 0.902 0.1028 86 72 1 0.889 0.1167 88 67 1 0.876 0.1316 89 64 2 0.848 0.1629 92 54 1 0.833 0.1814 94 50 1 0.816 0.2014 96 47 1 0.799 0.2226 102 43 1 0.780 0.2459 103 41 3 0.723 0.3191 104 38 1 0.704 0.3454 Tabla 4.9: Estimación de la función de riesgo acumulado para el grupo tratado con fármaco time n.risk n.event survival hazard 40 198 1 0.995 0.0051 49 196 1 0.990 0.0102 50 195 1 0.985 0.0153 54 191 1 0.980 0.0205 55 188 1 0.974 0.0258 64 184 1 0.969 0.0313 66 182 1 0.964 0.0368 68 181 1 0.958 0.0423 71 176 1 0.953 0.0480 73 173 1 0.948 0.0538 75 168 1 0.942 0.0597 77 164 1 0.936 0.0658 78 158 1 0.930 0.0721 79 156 1 0.924 0.0785 81 149 2 0.912 0.0920 84 142 2 0.899 0.1060 96 111 1 0.891 0.1151 101 98 1 0.882 0.1253 102 96 1 0.873 0.1357 Tabla 4.10: Estimación de la función de riesgo acumulado para el grupo de control CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 44 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA Por último, se añade una gráfica en la que se muestra la estimación. Para ello, basta añadir el argumento fun = “cumhaz” en la función ggsurvplot() con la que se ha trabajado en secciones anteriores: ggsurvplot( ajuste.rats, fun = "cumhaz", xlab = "t", ylab = "H(t) estimada", lab = c(10,10,7), size = 1,legend.title = "Tratamientos", palette = c("burlywood3","lightslateblue"), conf.int = TRUE, legend.labs = c("Fármaco","Control"), ggtheme = theme_grey(), caption = "Fuente propia", font.caption = c(10,"italic","black"), font.x = c(12,"bold.italic","black"), font.y = c(12,"bold.italic","black") ) + + +++ + + ++ +++ ++ ++ ++ + + ++++ +++++++++++ +++++ + +++ +++ + ++++++++++ + ++ + +++ ++ ++ ++ + + + 0.0 0.1 0.2 0.3 0.4 0.5 0 25 50 75 100 t H(t) estimada Tratamientos ++ Fármaco Control Fuente propia Figura 4.3: Estimador de la tasa de fallo acumulado para cada uno de los tratamientos estudiados en el conjunto de ratas disponibles en rats del paquete survival CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 45 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA 4.1.1.4. Comparación de curvas de supervivencia Uno de los objetivos del análisis de supervivencia es el estudio de la igualdad de dos o más curvas de supervivencia asociadas a distintas poblaciones en estudio (véase el contraste (3.1)). En el paquete survival, se encuentra la función survdiff(formula, data, subset, na.action, rho), la cual es muy útil para dicho objetivo, ya que comprueba si existen diferencias significativas entre dos o más curvas de supervivencia con el uso de la familia Gρdefinida por Fleming y Harrington (1982). Además, la función en cuestión puede comparar una sola curva de supervivencia con una curva alternativa y conocida. Tal y como se ha especificado, esta función survdiff() implementa la familia de test Gρde Harrington y Fleming con pesos en cada instante de fallo de {SKM (t)}ρsiendo SKM (t)el estimador de Kaplan-Meier de la función de supervivencia. Los argumentos de la función en estudio son: formula: objeto fórmula del tipo Surv y de la forma Surv(time, status) ∼predictor. data: objeto data.frame en el cual se encuentran los datos. subset: argumento que indica qué subconjunto de las filas de datos debe utilizarse en el ajuste (vector lógico, vector numérico o vector de caracteres). Si no se indica nada, se incluyen todas las observaciones. na.action: argumento para filtrar los datos perdidos. rho: valor numérico que controla el tipo de test a aplicar. Si toma el valor 0 (por defecto) se usará el test log-rank y cuando toma el valor 1 se usará la modificación de Peto-Peto de la prueba de Gehan-Breslow, i.e., test log-rank ponderado con los pesos de Peto-Peto modificados (véase la tabla de pesos 3.1). A continuación, se aplicará la función al conjunto de datos colon (Laurie et al. (1989)) del paquete survival. Los datos pertenecen a uno de los primeros ensayos que se realizaron sobre la aplicación a pacientes con cáncer de colon de distintos tratamientos. A una parte de la muestra de datos se le aplicó levamisol (compuesto de baja toxicidad que se usaba para tratar infecciones en animales y representado por Lev), a otra parte se le aplicó una combinación de levamisol y 5-FU (agente quimioterapéutico moderadamente tóxico y representado por Lev+5FU) y a la última parte de la muestra no se le aplicó ningún tratamiento (representado por Obs). Además, para cada paciente existen dos registros, uno para la recurrencia y otro para la muerte. Solo se trabajará con los datos de recurrencia y con las variables rx (tratamiento), time (días hasta que se produce el evento o la censura) y status (estado de la censura). colon.rec <- colon[colon$etype == 1,] time status rx 968 1 Lev+5FU 3087 0 Lev+5FU CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 46 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA (continúa) time status rx 542 1 Obs 245 1 Lev+5FU 523 1 Obs 904 1 Lev+5FU 229 1 Lev 3192 0 Obs 3173 0 Lev 3308 0 Lev+5FU 2908 0 Lev 3309 0 Lev+5FU 1130 1 Obs 2231 1 Lev Tabla 4.11: Muestra de datos del conjunto colon contenido en el paquete survival El objetivo que se persigue es estudiar el contraste siguiente: H0:SObs(t) = SLev(t) = SLev+5F U (t)para todo t∈T H1:Si(t)=Sj(t)para algún i=j, t ∈T(4.1) donde i, j =Obs, Lev, Lev + 5FU. A continuación, se resuelve dicho contraste en R para ρ= 0 (test log-rank ponderado con pesos iguales): survdiff(Surv(time, status) ~ rx, data = colon.rec, rho = 0) ## Call: ## survdiff(formula = Surv(time, status) ~ rx, data = colon.rec, ## rho = 0) ## ## N Observed Expected (O-E)^2/E (O-E)^2/V ## rx=Obs 315 177 151 4.62 6.82 ## rx=Lev 310 172 149 3.69 5.41 ## rx=Lev+5FU 304 119 169 14.69 23.04 ## ## Chisq= 23.1 on 2 degrees of freedom, p= 1e-05 Antes de concluir si el contraste es significativo o no, se detalla el significado de cada una de las componentes de la salida obtenida: N: número de individuos u observaciones existentes en cada grupo o población en estudio. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 47 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA Observed: número ponderado de eventos observados en cada uno de los grupos en estudio. Expected: número ponderado y esperado de eventos en cada uno de los grupos en estudio. (O−E)2/E: representa el cálculo del estadístico χ2para cada uno de los grupos en estudio. (O−E)2/V : valor numérico que se calcula elevando al cuadrado la diferencia de las componentes Observed yExpected y dividiendo entre la varianza. chisq: valor del estadístico χ2del test junto con los grados de libertad asociados. var: matriz de varianza del test. p: p-valor asociado al contraste. Se concluye así que el contraste (4.1) realizado tiene asociado un p-valor de 0.00001, por lo que, existen evidencias significativas en contra de la hipótesis nula o igualdad de las curvas de supervivencia de los tres tratamientos (Obs,Lev yLev+5FU) en estudio. En el caso de que se quisiera obtener un gráfico en el que se pudiera consultar la estimación de la curva de supervivencia asociada a cada tratamiento junto con el p-valor vinculado al contraste de hipótesis, se usará la función ggsurvplot() del paquete survminer utilizada en apartados anteriores. En este caso, a dicha función es necesario añadirle el argumento pval = TRUE ylog.rank.weights para indicarle el tipo de peso que se quiere usar para el test log-rank ponderado. Además, con el argumento risk.table = TRUE se adjunta una tabla en la que se puede consultar el número de individuos en riesgo en cada uno de los instantes considerados. Para una consulta más detallada sobre los posibles argumentos que admite la función ggsurvplot() véase el apartado ggsurvplot_arguments del manual de Biecek et al. (2021). A continuación, se muestra la gráfica asociada al contraste (4.1) en el que se puede visualizar el p-valor obtenido y la tabla de individuos en riesgo: fit <- survfit(Surv(time, status) ~ rx, data = colon.rec, type = "kaplan-meier",conf.type = "plain") ggsurvplot( fit, data = colon.rec, size = 1, conf.int = TRUE, pval = TRUE, legend.labs = c("Obs","Lev","Lev+5FU"), ylim = c(0.1,1), xlim = c(0,3250), lab = c(10,10,7), xlab = "t (días)", ylab = "S(t) estimada", CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 48 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA Por otro lado, los resultados obtenidos por las pruebas de tendencia son: Q Var Z pNorm 1 -25.8061 48.1505 -3.7190 0.00020005 *** n -1939.0000 210644.6590 -4.2248 2.3919e-05 *** sqrtN -221.9185 2990.6169 -4.0580 4.9493e-05 *** S1 -21.3895 26.8311 -4.1293 3.6380e-05 *** S2 -21.0942 26.0075 -4.1363 3.5292e-05 *** FH_p=1_q=1 -2.9412 1.5056 -2.3970 0.01653027 * Tabla 4.14: Segunda tabla de la salida comp en la que se contemplan diferentes pruebas de tendencia En cada una de las filas se realiza una prueba de tendencia concreta con distintos pesos. Dichos pesos son los mismos que aparecen en la Tabla 4.13 y que se han especificado antes. Por otro lado, en la columna Zse muestra el estadístico (3.4) de la prueba de tendencia y en la columna pNorm se consulta el p-valor asociado al contraste. Nótese que el contraste que se está realizando es: H0:S1(t) = S2(t) = S3(t) = S4(t)t≤τ H1:S1(t)≥S2(t)≥S3(t)≥S4(t)t≤τ siendo τla mayor t∈Ten la que todos los estados involucrados en el estudio tienen al menos un sujeto en riesgo. Si se observan los p-valores, todos ellos son menores que 0.05, así que claramente se rechaza la hipótesis nula. Otra forma de obtener el p-valor asociado a la prueba de tendencia es de manera gráfica con la función ggsurvplot() de la librería survminer con la que ya se ha trabajado. En este caso, para visualizar el p-valor es necesario añadirle a la función los argumentos pval = TRUE ytest.for.trend = TRUE. Además, con el argumento log.rank.weights se puede modificar el tipo de peso que se quiere usar en el test y con la orden arrange_ggsurvplots() se mostrarán varias gráficas en la misma figura. En las siguientes gráficas se muestran los mismos p-valores que se habían obtenido en la Tabla 4.14: # Ajuste de la función de supervivencia: ajuste1 <- survfit(Surv(time, delta) ~ stage, data = larynx, type = "kaplan-meier") splots <- list() # Prueba de tendencia con pesos iguales: splots[[1]] <- ggsurvplot( ajuste1, data = larynx, pval = TRUE, test.for.trend = TRUE, pval.method = TRUE, log.rank.weights = "1" ) CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 55 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA # Test de Gehan-Breslow: splots[[2]] <- ggsurvplot( ajuste1, data = larynx, pval = TRUE, test.for.trend = TRUE, pval.method = TRUE, log.rank.weights = "n" ) # Test de Tarone-Ware: splots[[3]] <- ggsurvplot( ajuste1, data = larynx, pval = TRUE, test.for.trend = TRUE, pval.method = TRUE, log.rank.weights = "sqrtN" ) # Test de Peto-Peto: splots[[4]] <- ggsurvplot( ajuste1, data = larynx, pval = TRUE, test.for.trend = TRUE, pval.method = TRUE, log.rank.weights = "S1" ) # Test de Peto-Peto modificado: splots[[5]] <- ggsurvplot( ajuste1, data = larynx, pval = TRUE, test.for.trend = TRUE, pval.method = TRUE, log.rank.weights = "S2" ) # Test de Fleming-Harrington con p = 1 y q = 1: splots[[6]] <- ggsurvplot( ajuste1, data = larynx, pval = TRUE, test.for.trend = TRUE, pval.method = TRUE, log.rank.weights = "FH_p=1_q=1" ) arrange_ggsurvplots(splots, print = TRUE, ncol = 2,nrow = 3) CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 56 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA ++ + +++ ++ ++++ + + + ++ + ++ + ++ + + ++ ++ ++ + + + +p = 2e−04 Log−rank, tft 0.00 0.25 0.50 0.75 1.00 0 2.5 5 7.5 10 t (mes) S(t) estimada Estado ++++ S1 S2 S3 S4 ++ + +++ ++ ++++ + + + ++ + ++ + ++ + + ++ ++ ++ + + + +p < 0.0001 Gehan−Breslow, tft 0.00 0.25 0.50 0.75 1.00 0 2.5 5 7.5 10 t (mes) S(t) estimada Estado ++++ S1 S2 S3 S4 ++ + +++ ++ ++++ + + + ++ + ++ + ++ + + ++ ++ ++ + + + +p < 0.0001 Tarone−Ware, tft 0.00 0.25 0.50 0.75 1.00 0 2.5 5 7.5 10 t (mes) S(t) estimada Estado ++++ S1 S2 S3 S4 ++ + +++ ++ ++++ + + + ++ + ++ + ++ + + ++ ++ ++ + + + +p < 0.0001 Peto−Peto, tft 0.00 0.25 0.50 0.75 1.00 0 2.5 5 7.5 10 t (mes) S(t) estimada Estado ++++ S1 S2 S3 S4 ++ + +++ ++ ++++ + + + ++ + ++ + ++ + + ++ ++ ++ + + + +p < 0.0001 modified Peto−Peto, tft 0.00 0.25 0.50 0.75 1.00 0 2.5 5 7.5 10 t (mes) S(t) estimada Estado ++++ S1 S2 S3 S4 ++ + +++ ++ ++++ + + + ++ + ++ + ++ + + ++ ++ ++ + + + +p = 0.016 Fleming−Harrington (p=1, q=1), tft 0.00 0.25 0.50 0.75 1.00 0 2.5 5 7.5 10 t (mes) S(t) estimada Estado ++++ S1 S2 S3 S4 Figura 4.6: Obtención de los p-valores asociados a distintos pesos en pruebas de tendencia de manera gráfica con la librería survminer CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 57 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA 4.1.3. Librería KMsurv La librería KMsurv (Klein y Moeschberger (2012)) se caracteriza por poseer una gran selección de conjuntos de datos procedentes de estudios realizados y ejercicios propuestos en Klein y Moeschberger (1997). Además, incluye una función que es capaz de generar una tabla de vida o estimador actuarial a partir de un par de variables, el número de censuras y el número de eventos que se han producido en cada instante de tiempo. library(KMsurv) 4.1.3.1. Estimador actuarial Tal y como se especificó en la Sección 2.1.2 del Capítulo 2, una forma de calcular un estimador de la función de supervivencia cuando los tiempos vienen agrupados en intervalos es mediante el estimador actuarial o tabla de vida. En el paquete en estudio, existe la función lifetab(tis, ninit, nlost, nevent) capaz de calcular dicho estimador. A continuación, se especifica el significado de cada uno de los argumentos de la función anterior: tis: vector que contiene los instantes finales de los intervalos de tiempo involucrados en el estudio. Dicho vector debe contener un elemento más de los que poseen los vectores nlost ynevent. ninit: número de individuos que intervienen en el inicio del estudio. nlost: vector que contiene el número de individuos que abandonan el estudio o cuyos datos se pierden a lo largo del estudio. nevent: vector que posee el número de individuos que experimentaron el evento en estudio. Además, como salida se obtiene una tabla con diez columnas en las que se puede consultar la estimación de la función de densidad, de la función de supervivencia o de la función de riesgo, entre otra información disponible. El significado de cada una de las columnas se puede consultar en la página 9 del Capítulo 2. Se adjunta un ejemplo de aplicación de la función en cuestión. Para ello, se usarán los datos nurshome que se pueden consultar en el siguiente enlace https://www2.karlin.mff. cuni.cz/~pesta/prednasky/NMFM404/Data/NursingHome.dat.The National Center for Health Services Research estudió 36 residencias de ancianos con el objetivo de evaluar los efectos que tenían ciertos incentivos económicos en la duración de la estancia de cada uno de los residentes. Las residencias “incentivadas” recibían más dinero por los pacientes de Medicaid y bonificaciones por mejorar la salud de un paciente y por darles el alta. En el estudio participaron 1601 pacientes ingresados en dichas residencias entre el 1 de mayo de 1981 y el 30 de abril de 1982. Las variables que se consideraron para el estudio fueron: los: duración de la estancia en días. age: edad del residente. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 58 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA rx: indica el tipo de residencia. Si está incentivada toma el valor 1 y en caso contrario toma el valor 0. gender: especifica el género del residente. El valor 1 indica que el residente es hombre y el 0 que es mujer. married: indicador de si el residente estaba casado o no. Toma el valor 1 si está casado y el 0 en caso contrario. health: indicador del estado de salud del residente. Se indica a través de una escala numérica del 2 al 5 donde el 2 indica que el estado de salud es el mejor y 5 indica el peor estado de salud. censor: indicador de la censura (1 = censurado o 0 = no censurado). Se procede a la lectura de los datos y a la muestra de algunos de ellos: nurshome <- read.table("http://www.karlin.mff.cuni.cz/~pesta/prednasky/ NMFM404/Data/NursingHome.dat", header = FALSE, skip = 14, col.names = c("los","age","rx", "gender","married", "health","censor")) los age rx gender married health censor 864 94 1 1 0 2 1 18 86 1 1 0 4 0 1084 75 1 0 0 3 1 1092 77 1 0 1 2 1 23 86 1 0 0 4 0 971 82 1 1 1 4 1 176 93 0 0 0 5 0 396 86 0 0 1 3 1 533 89 0 0 0 5 1 25 65 0 0 0 4 0 2 84 0 1 0 5 0 38 95 0 1 0 3 0 Tabla 4.15: Muestra de datos del conjunto Nurshome localizado en https://www2.karlin. mff.cuni.cz/~pesta/prednasky/NMFM404/Data/NursingHome.dat Para calcular un estimador actuarial, los datos se restringirán a aquellos en los que la variable rx toma el valor 1, es decir, aquellas residencias de ancianos que han sido incentivadas: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 59 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA nurshome_inc <- nurshome[nurshome$rx == 1, ] Además, como ya se ha indicado anteriormente, se necesita que los datos estén agrupados por intervalos, así que lo que se hará será agrupar la duración de la estancia de cada uno de los residentes en intervalos de 50 días. Para ello, se hará uso de la función floor() y se creará una nueva variable int que contendrá el extremo izquierdo del intervalo al que pertenece cada uno de los residentes o individuos de la muestra: nurshome_inc$int <- floor(nurshome_inc$los/50)*50 A partir de esos extremos, se deben calcular los argumentos necesarios para introducir en la función lifetab(): Cálculo del vector tis, el cual debe contener el extremo derecho de cada uno de los intervalos considerados. Como el tiempo máximo que ha estado un paciente ingresado ha sido 1092, basta considerar el vector que comienza en el 0 y acaba en 1050 siendo los saltos de amplitud 50. El último intervalo a considerar será (1050,+∞): extr_derecho <- c(seq(0,1050,by = 50), Inf) Cálculo de ninit: nTotal <- nrow(nurshome_inc) Cálculo de los vectores nlost ynevent. Se necesita una tabla de frecuencias en la que se contabilice para cada intervalo de tiempo considerado, el número de individuos censurados y completos. Para ello, se hará uso de la función table(): tabla <- table(nurshome_inc$int, nurshome_inc$censor) Una vez calculada la tabla, el vector que contabiliza el número de individuos censurados (censura) y no censurados (evento) en cada uno de los intervalos será: censura <- tabla[, 2] evento <- tabla[, 1] Una vez que se tienen los datos preparados, se procede a calcular la tabla de vida o estimador actuarial: tabla_vida <- lifetab(tis = extr_derecho, ninit = nTotal, nlost = censura, nevent = evento) En la siguiente página se puede consultar el resultado obtenido del estimador calculado: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 60 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA Tabla 4.16: Estimador actuarial del conjunto de datos Nurshome nsubs nlost nrisk nevent surv pdf hazard se.surv se.pdf se.hazard 0-50 712 0 712.0 241 1.00000 0.00677 0.00815 0.00000 0.00035 0.00051 50-100 471 0 471.0 89 0.66152 0.00250 0.00417 0.01773 0.00025 0.00044 100-150 382 0 382.0 48 0.53652 0.00135 0.00268 0.01869 0.00019 0.00039 150-200 334 0 334.0 38 0.46910 0.00107 0.00241 0.01870 0.00017 0.00039 200-250 296 0 296.0 35 0.41573 0.00098 0.00251 0.01847 0.00016 0.00042 250-300 261 0 261.0 30 0.36657 0.00084 0.00244 0.01806 0.00015 0.00044 300-350 231 0 231.0 15 0.32444 0.00042 0.00134 0.01755 0.00011 0.00035 350-400 216 0 216.0 23 0.30337 0.00065 0.00225 0.01723 0.00013 0.00047 400-450 193 1 192.5 17 0.27107 0.00048 0.00185 0.01666 0.00011 0.00045 450-500 175 0 175.0 15 0.24713 0.00042 0.00179 0.01617 0.00011 0.00046 500-550 160 0 160.0 5 0.22595 0.00014 0.00063 0.01568 0.00006 0.00028 550-600 155 0 155.0 8 0.21889 0.00023 0.00106 0.01550 0.00008 0.00037 600-650 147 0 147.0 3 0.20759 0.00008 0.00041 0.01521 0.00005 0.00024 650-700 144 0 144.0 10 0.20335 0.00028 0.00144 0.01509 0.00009 0.00045 700-750 134 10 129.0 7 0.18923 0.00021 0.00112 0.01469 0.00008 0.00042 750-800 117 20 107.0 3 0.17896 0.00010 0.00057 0.01440 0.00006 0.00033 800-850 94 13 87.5 1 0.17394 0.00004 0.00023 0.01428 0.00004 0.00023 850-900 80 16 72.0 3 0.17196 0.00014 0.00085 0.01426 0.00008 0.00049 900-950 61 15 53.5 4 0.16479 0.00025 0.00155 0.01425 0.00012 0.00078 950-1000 42 15 34.5 0 0.15247 0.00000 0.00000 0.01446 NaN NaN 1000-1050 27 14 20.0 0 0.15247 0.00000 0.00000 0.01446 NaN NaN 1050-Inf 13 13 6.5 0 0.15247 NA NA 0.01446 NA NA CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 61 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA A partir de la Tabla 4.16, se puede graficar la estimación de la función de supervivencia asociada a aquellas residencias incentivadas: plot(tabla_vida$surv, type = "s",xlab = "t (días)", ylab = "S(t) estimada",col = "darkslateblue") 5 10 15 20 0.2 0.4 0.6 0.8 1.0 t (días) S(t) estimada Fuente propia Figura 4.7: Estimador actuarial para el grupo de residencias de ancianos incentivadas del conjunto de datos Nurshome Cabe destacar también, que en el paquete biostat3 (Tillander et al. (2020)) existe otra función capaz de calcular la tabla de vida. Esta función se llama lifetab2 (formula, data, breaks) y la principal diferencia que presenta con respecto a la función lifetab() del paquete KMsurv es que no es necesario preparar los datos en intervalos. Basta introducir como primer argumento un objeto Surv como formula; como segundo argumento, el data.frame de donde proceden los datos; y como último argumento, un vector con los extremos izquierdos de los intervalos a considerar. Con el siguiente código se obtiene la misma tabla de vida que la Tabla 4.16. Nótese que se ha tenido que recodificar la variable censor debido a que la función Surv considera como dato censurado aquellos que toman el 0 en la variable censor en vez del valor 1 como considera el conjunto de datos Nurshome. library(biostat3) nurshome_inc$censor <- car::Recode(nurshome_inc$censor, "1 ='0'; 0 ='1'") tabla_vida1 <- lifetab2(Surv(los, censor) ~ 1,data = nurshome_inc, breaks = seq(0,1050,by = 50)) CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 62 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA En este caso, para mostrar la gráfica del estimador de la función de supervivencia, basta realizar lo siguiente: plot(tabla_vida1, type = "s") 4.1.4. Librería gte La librería gte (Dehghan et al. (2015)) fue creada para obtener el estimador de Turnbull propuesto por Dehghan y Duchesne (2011) con una función capaz de realizar el cálculo de dicho estimador denominada gte. A continuación, se procede a cargar la librería: library(gte) 4.1.4.1. Estimador de Turnbull Tal y como se ha especificado, la función que calcula el estimador de Turnbull de la función de supervivencia es gte(formula, data, z, itermax, tole). Los argumentos de dicha función son: formula: objeto de la forma y∼xdonde yserá un objeto de tipo Surv que actúa como variable respuesta y xserá el nombre de las covariables. data: objeto data.frame que contiene el conjunto de datos en estudio. z: debe ser una matriz en la que cada fila contenga los valores de un vector de covariables con los que se calculan la función de supervivencia condicionada. En el caso de que exista una sola covariable, zserá un vector. itermax: indica el número máximo de iteraciones que debe realizar el algoritmo (por defecto es 100000). tole: distancia máxima permitida entre iteraciones sucesivas antes de considerar la convergencia del método. Por defecto toma el valor 0.0005. El estimador en cuestión es una generalización del estimador de Turnbull y consiste en una estimación no paramétrica de la función de supervivencia de la variable aleatoria T condicionada a un valor z0del vector de covariables Z, es decir, se pretende construir un estimador no paramétrico de S(t|z0) = P[T > t |Z=z0]. Para ver el funcionamiento de la función, se aplicará al conjunto de datos simul del paquete en estudio. Dicho conjunto de datos contiene 100 observaciones censuradas por intervalos y tres variables: L: variable que contiene los extremos izquierdos de cada uno de los intervalos de censura. R: variable que posee cada uno de los extremos derechos de los intervalos de censura. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 63 4.1. LIBRERÍAS PARA EL ANÁLISIS DE SUPERVIVENCIA Z: covariable continua. Nótese que cuando R=NA, la observación en cuestión está censurada a la derecha y cuando L=NA, la observación está censurada por la izquierda. Además, si se cumple que L=R, entonces el instante en el que se ha producido el evento se conoce con exactitud. A continuación, se leerá el conjunto de datos y se adjuntará una muestra de ellos: data(simul) L R Z 41.8099110 50.510700 18.122569 NA 36.926040 13.889860 2.8240208 37.949810 13.859484 NA 10.669690 5.353621 15.8330992 21.621300 17.423702 8.0705407 19.397390 5.345955 NA 22.475340 13.843842 46.2031487 NA 21.199203 51.6190973 NA 23.548849 12.7220895 28.013940 21.302348 17.4136333 25.578350 14.021004 Tabla 4.17: Muestra de datos del conjunto simul del paquete gte Una vez mostrados los datos, se procede a generar el código necesario para obtener el estimador de Turnbull de los datos. En el objeto Surv que se debe introducir como primer argumento en la función gte() hay que especificar que el type es interval2 debido a que los datos están censurados por intervalos y son indicados por el valor NA: turnbull <- gte(Surv(L, R, type = "interval2") ~ Z, data = simul, z = c(15,30,45)) Como resultado se obtiene una tabla en la que se pueden consultar cada una de las estimaciones realizadas. Se adjunta una tabla del resultado obtenido: survival|Z=z[1] survival|Z=z[2] survival|Z=z[3] [8.6042993, 8.6042993] 8.260578e-01 1.0000000 1.0000000 (10.0582938, 10.26988] 8.260488e-01 0.9999999 1.0000000 (10.5500817, 10.66969] 8.260488e-01 0.9999998 1.0000000 [13.5215962, 13.5215962] 5.895147e-01 0.9999678 1.0000000 (14.5131108, 14.53361] 5.895147e-01 0.9999367 1.0000000 CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 64 4.2. LIBRERÍA FHTEST + ++ ++ ++ + + +++ + ++ + + ++ ++ ++++ + +++++ ++ ++ + ++ ++ ++++ ++++ + +++ +++ +++ ++ + 0.00 0.25 0.50 0.75 1.00 0 250 500 750 1000 t (días) S(t) estimada Género ++ Masculino Femenino Fuente propia Figura 4.10: Curvas de supervivencia del conjunto de datos lung del paquete survival por género Según la gráfica, el tiempo de supervivencia asociado a las mujeres es mayor durante todo el período en estudio aunque cabe destacar que al final del estudio, la diferencia entre ambas curvas de supervivencia es pequeña. 4.2. Librería FHtest La librería FHtest (Oller y Langohr (2020)) fue implementada en el año 2017 por Ramon Oller y Klaus Langohr debido a la falta de funciones en R para test log-rank ponderados con datos censurados por intervalos basados en la clase de Fleming-Harrington. Además, incorpora funciones para trabajar con datos censurados por intervalos y datos censurados a la derecha. Por otro lado, proporciona un enfoque integrado para realizar pruebas de dos muestras, kmuestras y de tendencia basadas en la teoría de proceso de conteo, teoría de verosimilitud o distribuciones de permutación. Una peculiaridad que posee el paquete es que necesita un paquete externo llamado Icens (Bioconductor (2021b)) y que para ser instalado es conveniente consultar el siguiente enlace http://www.bioconductor.org/packages/release/bioc/html/Icens.html. Una forma de descargar dicho paquete es ejecutando las siguientes órdenes: source("http://bioconductor.org/biocLite.R") if (!requireNamespace("BiocManager",quietly = TRUE)) install.packages("BiocManager") BiocManager::install(version = "3.12") BiocManager::install(c("GenomicFeatures","AnnotationDbi")) biocLite("Icens") CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 71 4.2. LIBRERÍA FHTEST Una vez descargado el paquete necesario, se procede a cargar el paquete FHtest: library(FHtest) 4.2.1. Comparación de curvas de supervivencia La librería en estudio posee cuatro funciones destinadas a llevar a cabo las técnicas específicas de comparación de curvas de supervivencia con datos censurados por intervalos y a la derecha estudiadas en la Sección 3.2 del Capítulo 3. Las cuatro funciones permiten realizar pruebas o test de dos poblaciones o grupos, kpoblaciones y de tendencia. Además, en los contrastes a realizar también se pueden contemplar el uso de hipótesis alternativas con desigualdades. Debido a que dichas técnicas se basan en los pesos de Fleming-Harrington, la selección de los parámetros λyρdebe ser adecuada. Por ejemplo, si se está trabajando en ensayos clínicos y se quiere evaluar si el efecto de un tratamiento o terapia sobre la supervivencia es mayor en las primeras fases de la terapia, se deberá elegir λ= 0 yρde tal forma que cuánto mayor sea su valor, más se enfatizará en las diferencias tempranas. Sin embargo, en el caso de que existan evidencias médicas para creer que el efecto de la terapia está más pronunciado en la mitad o el final del período de seguimiento, se deberá tomar ρ=λ > 0 oρ= 0, respectivamente. Además, cuanto más elevado sea el valor de λmás se enfatizará en las diferencias medias o tardías. 4.2.1.1. Datos censurados a la derecha A la hora de abordar la comparación de curvas de supervivencia con datos censurados a la derecha, se puede utilizar la función FHtestrcc() o la función FHtestrcp(). Según el enfoque con el que se quiera resolver el contraste se utilizará una u otra función. La primera de ellas aplica un enfoque basado en el proceso de conteo mientras que la segunda aplica un enfoque de permutación. 4.2.1.1.1. Enfoque mediante el proceso de conteo La función asociada a este enfoque se llama FHtestrcc y tal y como se especificó en la Sección 3.2.1.1, el método usa la familia de estadísticos Gρ,λ de Fleming-Harrington para detectar las diferencias entre dos o más curvas de supervivencia y se puede aplicar de dos formas distintas. La primera de ellas sigue el siguiente esquema FHtestrcc(L, R, group, rho, lambda, alternative, method) mientras que el segundo esquema es FHtestrcc(formula, data, subset, na.action, rho, lambda, alternative, method). A continuación, se especifica el significado de cada uno de los argumentos: L: vector numérico que contiene los extremos izquierdos de los intervalos con los que se va a trabajar. R: vector numérico que contiene los extremos derechos de los intervalos de la muestra en estudio. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 72 4.2. LIBRERÍA FHTEST group: vector que denota la variable para la que se desea efectuar el contraste. Además, determina el tipo de test a realizar. En el caso de que group sea de tipo factor con k≥2posibles valores, entonces se llevará a cabo un test de kpoblaciones. Sin embargo, si group es un vector numérico, entonces se llevará a cabo una prueba de tendencia. rho: valor de la constante ρdel peso de Fleming-Harrington (véase (3.12)). Por defecto, toma el valor 0. lambda: valor de la constante λdel peso de Fleming-Harrington. Por defecto, toma el valor 0. alternative: carácter que proporciona el tipo de hipótesis alternativa (H1) que se quiere utilizar para los test con dos poblaciones o pruebas de tendencia. Los posibles valores son "different" (test bilateral), "increasing" (test unilateral con ≤) o "decreasing" (test unilateral con ≥). method: indicador del tipo de método usado para obtener el p-valor del test. Por defecto toma el valor "pclt" y usa el teorema central del límite con permutaciones. Otros posibles valores pueden ser "exact.mc" que usa un método exacto vía Montecarlo, "exact.network" que utiliza un método exacto a través de un algoritmo de red o "exact.ce" que a su vez usa un método exacto pero en este caso usando enumeración completa, lo que provoca que solo sea apropiada para muestras de tamaño pequeño. Nótese que los dos últimos posibles valores solo son compatibles con test de 2 poblaciones. formula: objeto de la forma y∼xdonde yserá un objeto de tipo Surv yxserá la variable que contiene los grupos a los que pertenece cada individuo u observación en estudio. data:data.frame donde se encuentran los datos con los que se va a trabajar. subset: vector opcional que indica el subconjunto de individuos u observaciones que debe usar la función. na.action: función cuyo fin es indicar cómo se van a tratar aquellos valores perdidos NA. Observación: cuando un tiempo iestá censurado a la derecha, el elemento i-ésimo del vector Rtomará el valor +∞, de esta forma, el intervalo asociado a dicho tiempo será de la forma (a, +∞). Sin embargo, cuando el elemento i-ésimo del vector Ly el correspondiente al vector Rcoinciden, se dice que se conoce de manera exacta el instante en el que se ha producido el evento en estudio y el intervalo es de la forma (a, a). ▲ En la salida que proporciona la función se podrá consultar la siguiente información: N: número de observaciones existentes en cada una de las poblaciones o grupos en estudio. Observed: número ponderado de eventos observados en cada uno de los grupos en estudio. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 73 4.2. LIBRERÍA FHTEST Expected: número ponderado de eventos esperados en cada uno de los grupos en estudio. statistic: estadístico asociado al contraste. var: matriz de varianza del test. pvalue: p-valor asociado al contraste que se está realizando. Además, se especifican los valores de los parámetros ρyλusados para los pesos de Fleming-Harrington, el tipo de test que se está aplicando y la descripción de la hipótesis alternativa del contraste. También se pueden observar las cantidades (O−E)2/E y(O− E)2/V necesarias para el cálculo del estadístico. Para ilustrar su funcionamiento se mostrarán tres ejemplos, el primero de ellos consistirá en realizar un test de 2 poblaciones, el segundo constará de un contraste con k > 2 poblaciones y el tercero consistirá en una prueba de tendencia. Test con 2 poblaciones Se trabajará con el conjunto de datos ovarian (Edmonson et al. (1979)) disponible en el paquete survival y que contiene los tiempos de supervivencia de 26 pacientes con cáncer de ovario tratados cada uno de ellos con uno de los dos tratamientos disponibles con el fin de comparar la eficacia de ambos tratamientos. Entre las variables en estudio se encuentran: futime: tiempo de supervivencia. fustat: estado de la censura (0 = censurado o 1 = no censurado). rx: tratamiento administrado al paciente. Si toma el valor 1, se hablará de Tratamiento 1 y si toma el valor 2, se hablará entonces de Tratamiento 2. A continuación, se muestran algunos de los datos del conjunto ovarian: futime fustat rx 59 1 1 115 1 1 464 1 2 563 1 2 770 0 2 1106 0 1 769 0 2 855 0 1 Tabla 4.20: Muestra de datos del conjunto ovarian del paquete survival CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 74 4.2. LIBRERÍA FHTEST Para poder comparar la eficacia de los tratamientos se debe resolver el siguiente contraste: H0:STrat.1(t) = STrat.2(t)∀t > 0 H1:STrat.1(t)=STrat.2(t)para algún t > 0(4.5) Una vez modelado el problema, se procede a resolverlo con la función en estudio mediante el test log-rank ponderado con pesos iguales (test log-rank): FHtestrcc(Surv(futime, fustat) ~ rx, data = ovarian) ## ## Two-sample test for right-censored data ## ## Parameters: rho=0, lambda=0 ## Distribution: counting process approach ## ## Data: Surv(futime, fustat) by rx ## ## N Observed Expected O-E (O-E)^2/E (O-E)^2/V ## rx=1 13 7 5.23 1.77 0.596 1.06 ## rx=2 13 5 6.77 -1.77 0.461 1.06 ## ## Statistic Z= -1, p-value= 0.303 ## Alternative hypothesis: survival functions not equal Analizando los resultados obtenidos del contraste, se concluye que debido a que el pvalor es 0.303, no existen evidencias en contra de la igualdad de curvas de supervivencia. Para resolver el contraste sin necesidad de incluir una formula en la función es necesario que en el vector Lse encuentren los tiempos observados mientras que Rdebe contener los tiempos de supervivencia no censurados o el +∞en caso de que el tiempo esté censurado a la derecha. Siguiendo las indicaciones expuestas para datos no censurados y censurados a la derecha la orden sería: with(ovarian, FHtestrcc(L = futime, R = futime / fustat, group = rx)) ## ## Two-sample test for right-censored data ## ## Parameters: rho=0, lambda=0 ## Distribution: counting process approach ## ## Data: {futime,futime/fustat} by rx ## ## N Observed Expected O-E (O-E)^2/E (O-E)^2/V ## rx=1 13 7 5.23 1.77 0.596 1.06 ## rx=2 13 5 6.77 -1.77 0.461 1.06 ## ## Statistic Z= -1, p-value= 0.303 ## Alternative hypothesis: survival functions not equal CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 75 4.2. LIBRERÍA FHTEST Tal y como era de esperar, el resultado obtenido es el mismo. La prueba expuesta utiliza pesos iguales en todos los sumandos de (3.11), pero en algunas ocasiones puede ser de interés detectar diferencias entre las curvas de supervivencia a corto o largo plazo. En el caso de que se esté interesado en las diferencias a largo plazo, es conveniente dar más peso a las últimas diferencias entre el número de eventos observados y esperados. Para ello, puede tomarse por ejemplo λ= 1 yρ= 0: FHtestrcc(Surv(futime, fustat) ~ rx, data = ovarian, rho = 0,lambda = 1) ## ## Two-sample test for right-censored data ## ## Parameters: rho=0, lambda=1 ## Distribution: counting process approach ## ## Data: Surv(futime, fustat) by rx ## ## N Observed Expected O-E (O-E)^2/E (O-E)^2/V ## rx=1 13 1.11 1.11 -0.00447 1.80e-05 0.000102 ## rx=2 13 1.50 1.49 0.00447 1.34e-05 0.000102 ## ## Statistic Z= 0, p-value= 0.992 ## Alternative hypothesis: survival functions not equal En este caso, el estadístico de prueba es el recogido en la ecuación (3.8) y el p-valor es cercano a 1 por lo que se puede concluir que las diferencias a largo plazo entre las curvas de supervivencia asociadas a cada uno de los tratamientos son bastantes pequeñas. Véase gráficamente dicha diferencia: ++ + + + + + + +++ + 0.00 0.25 0.50 0.75 1.00 0 1 2 3 t (años) S(t) estimada Tratamientos ++ Trat.1 Trat.2 Fuente propia Figura 4.11: Curvas de supervivencia asociadas a cada uno de los tratamientos destinados a los pacientes de cáncer de ovario del conjunto de datos ovarian del paquete survival CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 76 4.2. LIBRERÍA FHTEST Se observa así que existió una mejor supervivencia a corto plazo entre los pacientes con el tratamiento 2. Sin embargo, después de un año y medio, las diferencias entre las curvas de supervivencia disminuyeron. Test con kpoblaciones Para ilustrar la prueba con k > 2poblaciones se utilizará el conjunto de datos bmt (Klein y Moeschberger (1997)) del paquete KMsurv. El conjunto de datos procede de un estudio sobre la supervivencia de 137 pacientes con leucemia tras realizarle un trasplante de médula ósea a cada uno de ellos y está compuesto por 22 variables, entre las que se encuentran: group: grupo de enfermedad al que pertenece cada paciente. Puede tomar el valor 1 para referirse a la leucemia linfoblástica aguda (ALL), el valor 2 para hacer referencia a la leucemia mieloide aguda de bajo riesgo (AML Low Risk) o el valor 3 para referirse a la leucemia mieloide aguda de alto riesgo (AML High Risk). t2: tiempo de supervivencia de cada individuo libre de la enfermedad, es decir, el tiempo hasta la recaída, la muerte o el final del estudio. d3: indicador del tipo de censura. Si la variable toma el valor 1, se produce la muerte o la recaída del individuo en cuestión y si toma el valor 0, el tiempo de supervivencia estará censurado. Se procede a la carga de los datos y se muestran algunos de ellos: data("bmt",package = "KMsurv") group t2 d3 1 2081 0 1 1496 0 1 609 1 1 172 1 2 860 0 2 1324 0 2 390 1 2 79 1 3 2430 0 3 1345 0 3 100 1 3 2 1 Tabla 4.21: Muestra de datos del conjunto bmt del paquete KMsurv CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 77 4.2. LIBRERÍA FHTEST Con el fin de que la interpretación de la salida sea más clara, se modificarán las etiquetas de la variable group: bmt$group <- factor(bmt$group, labels = c("ALL","AML low risk", "AML high risk")) Una vez que se tienen los datos preparados, se pretende comparar si el tiempo de supervivencia hasta la recaída de los pacientes es igual para todos los grupos de enfermedades o no. En ese caso, el contraste que se debe resolver es: H0:SALL(t) = SAML low risk(t) = SAML high risk ∀t > 0 H1:Si(t)=Sj(t)para algún t, i, j (4.6) donde i, j toman como posibles valores ALL,AML low risk yAML high risk. Para resolver el contraste anterior, se elegirá ρ=λ= 1 de tal forma que se dará más importancia a las diferencias existentes entre las curvas en torno a la mediana global: FHtestrcc(Surv(t2, d3) ~ group, data = bmt, rho = 1,lambda = 1) ## ## K-sample test for right-censored data ## ## Parameters: rho=1, lambda=1 ## Distribution: counting process approach ## ## Data: Surv(t2, d3) by group ## ## N Observed Expected O-E (O-E)^2/E (O-E)^2/V ## group=ALL 38 4.55 3.79 0.769 0.156 1.02 ## group=AML high risk 45 5.41 3.54 1.864 0.981 6.28 ## group=AML low risk 54 4.87 7.50 -2.633 0.924 8.99 ## ## Chisq= 9.9 on 2 degrees of freedom, p-value= 0.00697 ## Alternative hypothesis: survival functions not equal En la salida obtenida se observa que el estadístico de prueba (3.9) sigue asintóticamente una distribución χ2con k−1=2grados de libertad. Además, se muestra un p-valor muy pequeño, por lo que se debe rechazar la hipótesis nula. Así, se puede suponer que hay diferencias entre los tres grupos de enfermos de leucemia con respecto al tiempo de supervivencia libre hasta la recaída después de recibir un trasplante. Véanse dichas diferencias en la siguiente gráfica: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 78 4.2. LIBRERÍA FHTEST + ++ ++++ +++++ + + ++++++ ++++++++++ ++ ++++++ ++ + ++ + + + + + + ++ + + + 0.00 0.25 0.50 0.75 1.00 02468 t (años) S(t) estimada Enfermedad +++ ALL AML de bajo riesgo AML de alto riesgo Fuente propia Figura 4.12: Curvas de supervivencia asociadas a cada uno de los pacientes de leucemia del conjunto de datos bmt del paquete KMsurv En promedio, los tiempos de supervivencia libre de enfermedad fueron mayores entre los pacientes con AML de bajo riesgo y menores entre los pacientes con AML de alto riesgo. Prueba de tendencia Para poder aplicar una prueba de tendencia es imprescindible que la variable de agrupación group sea numérica. Este test puede ser más útil que la prueba con kpoblaciones porque la hipótesis alternativa del contraste puede contener desigualdades. Si se sigue trabajando con el conjunto de datos bmt, puede ser de interés resolver el siguiente contraste: H0:SALL(t) = SAML low risk(t) = SAML high risk ∀t > 0 H1:SAML low risk(t)≥SALL(t)≥SAML high risk ∀t > 0 con al menos una desigualdad estricta debido a que en la Figura 4.12 se observa que la curva de supervivencia asociada al grupo AML de bajo riesgo está por encima de la de ALL y esta a su vez por encima de la curva del grupo AML de alto riesgo. Como se necesita que la variable group sea numérica, se reconvierte a tipo numeric y se recodifica para que el valor 1 esté asociado al grupo AML de bajo riesgo; el 2 al grupo ALL; y el 3, al grupo AML de alto riesgo: bmt$group <- as.numeric(relevel(bmt$group, ref = "AML low risk")) Una vez preparados los datos, se procede a resolver el contraste. Para ello, se debe indicar que la hipótesis alternativa es decreciente (alternative = “decreasing”): CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 79 4.2. LIBRERÍA FHTEST FHtestrcc(Surv(t2, d3) ~ group, data = bmt, rho = 1,lambda = 1, alternative = "decreasing") ## ## Trend FH test for right-censored data ## ## Parameters: rho=1, lambda=1 ## Distribution: counting process approach ## ## Data: Surv(t2, d3) by group ## ## N Observed Expected O-E ## group=1 54 4.87 7.50 -2.633 ## group=2 38 4.55 3.79 0.769 ## group=3 45 5.41 3.54 1.864 ## ## Statistic Z= 3.1, p-value= 0.000883 ## Alternative hypothesis: decreasing survival functions (higher group implies earlier event times) Nótese que en la salida de este test se suprimen las columnas (O−E)2/E y(O−E)2/V ya que estas no son necesarias para el cálculo del p-valor. La función ha utilizado el estadístico de prueba (3.10) y el p-valor obtenido es 0.000883 por lo que existen evidencias en contra de H0. Obsérvese que el signo negativo en la columna O−Eindica que existen menos eventos observados que esperados bajo la hipótesis nula, lo que quiere decir que la supervivencia es mejor de la esperada. Por lo tanto, el aumento de las diferencias O−E del primer al último grupo avalan la hipótesis alternativa de funciones de supervivencia decrecientes. 4.2.1.1.2. Enfoque de permutación La función asociada a este enfoque se llama FHtestrcp y utiliza la familia de estadísticos Gρ,λ para detectar las diferencias entre dos o más curvas de supervivencia. Al igual que la función estudiada en la sección anterior, admite dos tipos de esquemas. El primero de ellos es FHtestrcp(L, R, group, rho, lambda, alternative, method, exact) y el segundo FHtestrcp(formula, data, subset, na.action, rho, lambda, alternative, method, exact). El significado de cada uno de los argumentos se puede consultar en la descripción detallada de la función FHtestrcc() a falta del argumento exact, el cual puede tomar el valor TRUE oFALSE según se desee resolver el contraste de manera exacta o no. En el caso de que method sea no nulo, este argumento es ignorado. Otro posible argumento es permcontrol que se encarga de controlar aspectos relacionados con los test de permutación. Al ejecutar la función se obtiene la siguiente información: N: número de observaciones en cada grupo de estudio. O-E: número ponderado de eventos observados menos el esperado en cada uno de los grupos en estudio. var: matriz de varianza del test. statistic: estadístico asociado al contraste. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 80 4.2. LIBRERÍA FHTEST ## Data: Surv(left, right, type = "interval2") by treatment ## ## N O-E ## treatment=Rad 46 -0.598 ## treatment=RadChem 48 0.598 ## ## Statistic Z= 3.9, p-value= 0.000113 ## Alternative hypothesis: survival functions not equal El p-valor asociado al contraste usando λ=ρ= 2 es 0.000113, por lo que se concluye que la hipótesis de igualdad de curvas debe ser rechazada. Si se quisiera estudiar la razón por la que se está rechazando la H0, se podría resolver un contraste con hipótesis alternativa de tipo decreasing: H0:SRad(t) = SRadChem(t)∀t > 0 H1:SRad(t)≥SRadChem(t)∀t > 0 FHtesticp(Surv(left, right, type = "interval2") ~ treatment, data = bcos, Rin = TRUE, Lin = TRUE, alternative = "decreasing", lambda = 2,rho = 2) ## ## Two-sample test for interval-censored data ## ## Parameters: rho=2, lambda=2 ## Distribution: permutation approach (asymptotic approximation using central limit theorem) ## ## Data: Surv(left, right, type = "interval2") by treatment ## ## N O-E ## treatment=Rad 46 -0.598 ## treatment=RadChem 48 0.598 ## ## Statistic Z= 3.9, p-value= 5.67e-05 ## Alternative hypothesis: decreasing survival functions (group=RadChem has earlier event times) En este contraste, el p-valor obtenido es mucho más pequeño que el nivel de significación prefijado, así que existen suficientes evidencias en contra de la hipóteis nula y por tanto, se puede concluir que en general los tiempos de supervivencia asociados al tratamiento de radioterapia son mayores que los asociados al tratamiento de radioterapia más quimioterapia, sobre todo en la mitad del período de seguimiento. Por consiguiente, da la impresión de que la efectividad del tratamiento en radioterapia es mayor. Véanse gráficamente las estimaciones de las curvas de supervivencia: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 87 4.2. LIBRERÍA FHTEST + ++ + + +++ +++ +++ ++ + + ++ ++ + 0.00 0.25 0.50 0.75 1.00 0 10 20 30 40 50 t (meses) S(t) estimada Tratamientos ++ Radioterapia Radioterapia+Quimio Fuente propia Figura 4.13: Curvas de supervivencia asociadas a cada uno de los tratamientos destinados a los pacientes de cáncer de mama del conjunto de datos bcos del paquete interval Se observa así que aproximadamente en los primeros siete meses de estudio, las dos curvas de supervivencia son prácticamente iguales. Sin embargo, en los siguientes diez meses la curva de supervivencia asociada al tratamiento de radioterapia más quimioterapia se sitúa por encima de la asociada al tratamiento de radioterapia. Finalmente, a largo plazo se observa que la curva de supervivencia asociada al tratamiento de radioterapia se sitúa por encima de la otra curva, por lo que, en general, la efectividad del tratamiento de radioterapia es mayor en dicho período de tiempo. Test con kpoblaciones El conjunto de datos illust3 del paquete FHtest que contiene datos procedentes de un ensayo clínico sobre el sida se utilizará para la ilustración del test con kpoblaciones y datos censurados por intervalos. El marco de datos anterior posee 1607 observaciones divididas en tres variables y como se ha especificado, procede de un ensayo clínico sobre el sida que comenzó a realizarse hacia el año 1987 por Volberding et al. (1995). Fue diseñado para estudiar los beneficios de la terapia con zidovudina en pacientes en la fase inicial de la infección por VIH y las tres variables en estudio son: left: extremo izquierdo de cada uno de los intervalos de tiempo en estudio. right: extremo derecho de cada uno de los intervalos de tiempo en estudio. group: grupo de tratamiento en el que el valor 1 hace referencia a terapia diferida (deferred therapy), el 2 hace referencia a dosis de 500 mg/día (500 mg/day dosage) y el 3 hace referencia a dosis de 1500 mg/día (1500 mg/day dosage). CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 88 4.2. LIBRERÍA FHTEST La variable censurada por intervalos es el tiempo (en meses) desde que el paciente comienza el tratamiento hasta que muestra síntomas del sida o muere. A continuación, se cargan los datos y se muestran algunos de ellos: data("illust3",package = "FHtest") left right group 7 999 1 9 35 1 2 2 1 0 999 2 2 8 2 16 45 2 8 9 3 46 55 3 70 999 3 Tabla 4.23: Muestra de datos del conjunto illust3 del paquete FHtest En este conjunto de datos, existen tiempos censurados por intervalos, tiempos censurados a la derecha, tiempos perdidos y tiempos sin censura. Los tiempos sin censura son aquellos en los que la variable left yright coinciden; los tiempos perdidos se dan cuando la variable left toma el valor 0 y right toma el valor 999; los censurados a la derecha toman el valor 999 en la variable right y un número distinto de 0 en la variable left y por último, los censurados por intervalo son aquellos que no verifican ninguna de las condiciones anteriores. Se convertirán los valores 999 a Inf para distinguir de manera más sencilla los tiempos censurados: library(dplyr) illust3 <- illust3 %> % mutate(right = ifelse(right == 999, Inf, right)) Una vez presentados y arreglados los datos, se procede a resolver el siguiente contraste: H0:S1(t) = S2(t) = S3(t)∀t > 0 H1:S1(t)=S2(t)=S3(t)para algún t > 0 donde S1hace referencia a la función de supervivencia asociada a la terapia diferida; S2 a la función de supervivencia asociada a la dosis de 500 mg/día; y S3, a la dosis de 1500 mg/día. Dado que la variable group es numérica, se necesita recodificar la variable: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 89 4.2. LIBRERÍA FHTEST illust3$group <- factor(illust3$group, labels = c("Deferred therapy", "500 mg/day dosage", "1500 mg/day dosage")) El contraste anterior se resolverá para los valores de ρ= 1 yλ= 1 de forma que se dé más importancia a las diferencias intermedias que a las diferencias a largo o corto plazo entre las curvas de supervivencia: FHtesticp(Surv(left, right, type = "interval2") ~ group, data = illust3, rho = 1,lambda = 1,Lin = TRUE, Rin = TRUE) ## ## K-sample test for interval-censored data ## ## Parameters: rho=1, lambda=1 ## Distribution: permutation approach (asymptotic approximation using central limit theorem) ## ## Data: Surv(left, right, type = "interval2") by group ## ## N O-E ## group=1500 mg/day dosage 528 -7.63 ## group=500 mg/day dosage 538 -1.72 ## group=Deferred therapy 541 9.35 ## ## Chisq= 15.7 on 2 degrees of freedom, p-value= 4e-04 ## Alternative hypothesis: survival functions not equal El estadístico asociado a la prueba sigue asintóticamente una distribución χ2con 2 grados de libertad y el p-valor obtenido es 0.0004, por lo que existen evidencias suficientes en contra de la igualdad de curvas de supervivencia. Se adjunta una gráfica de la estimación de las curvas de supervivencia asociadas a cada uno de los tratamientos en estudio. En ella se puede observar que aproximadamente durante el primer año de tratamiento, las curvas de supervivencia son prácticamente iguales. Sin embargo, posteriormente la curva asociada a la terapia diferida se sitúa por debajo de las otras dos, aunque la del grupo de las dosis de 1500 mg/día es algo mayor que la de 500 mg/día. Además, con la ayuda del contraste realizado se ha podido comprobar que realmente existen diferencias entre los tres tratamientos con respecto a la supervivencia de cada uno de los pacientes. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 90 4.2. LIBRERÍA FHTEST ++ + +++ + ++++ ++++ + + ++ ++++ + ++++ +++++++++++++++++ ++++++++++++++ +++ ++ + + + ++ +++ +++ + + + ++ + ++ +++ + ++ + ++ ++++ + + + +++++ ++ ++ + + ++++++ ++++ ++++++++++ ++ + +++ + + ++ +++ ++++ + +++++ + + + + +++ ++ ++ ++ + + ++++ + +++ ++++ +++++ + ++ +++++++++ + 0.00 0.25 0.50 0.75 1.00 0246 t (años) S(t) estimada Tratamientos +++ Terapia dif. dosis de 500 mg/día dosis de 1500 mg/día Fuente propia Figura 4.14: Curvas de supervivencia asociadas a cada uno de los tratamientos destinados a los pacientes de sida del conjunto de datos illust3 del paquete FHtest Prueba de tendencia Con la ayuda de una prueba de tendencia se estudiará si puede considerarse que el tiempo de supervivencia asociado al tratamiento de la dosis de 1500 mg/día es superior al resto. Se considerará el siguiente contraste: H0:S1(t) = S2(t) = S3(t)∀t > 0 H1:S1(t)≤S2(t)≤S3(t)∀t > 0 con al menos una desigualdad estricta. A la hora de resolver dicho contraste, se introducirá el argumento alternative = “increasing” y se convertirá a una variable numérica la variable de agrupación group: FHtesticp(Surv(left, right, type = "interval2") ~ as.numeric(group), data = illust3, Lin = TRUE, Rin = TRUE) ## ## Trend FH test for interval-censored data ## ## Parameters: rho=0, lambda=0 ## Distribution: permutation approach (asymptotic approximation using central limit theorem) ## ## Data: Surv(left, right, type = "interval2") by as.numeric(group) ## ## N O-E ## as.numeric(group)=1 541 51.4 ## as.numeric(group)=2 538 -8.0 CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 91 4.2. LIBRERÍA FHTEST ## as.numeric(group)=3 528 -43.4 ## ## Statistic Z= -4.2, p-value= 3.15e-05 ## Alternative hypothesis: survival functions not equal Se ha obtenido un p-valor cercano a 0, así que se puede suponer que efectivamente el tiempo de supervivencia de los pacientes con una dosis de 1500 mg/día es mayor que el tiempo de supervivencia asociado a los pacientes con una dosis de 500 mg/día y este tiempo es a su vez mayor que el asociado a la terapia diferida. 4.2.1.2.2. Enfoque basado en un vector de puntuación La función asociada a este enfoque se llama FHtestics, utiliza la familia de estadísticos Gρ,λ con λ= 0 para encontrar las diferencias entre dos o más curvas de supervivencia y al igual que todas las funciones del paquete en estudio, su aplicación puede llevarse a cabo de dos formas distintas. El primer esquema de aplicación es FHtestics(L, R, group, rho, alternative, tol, icFIT, icontrol, Lin, Rin) y el segundo FHtestics(formula, data, subset, na.action, rho, alternative, tol, icFIT, icontrol, Lin, Rin). El significado de cada uno de los argumentos se puede consultar en la descripción detallada de la función FHtesticp() a falta del argumento tol, el cual denota la tolerancia permitida para el cálculo de la inversa de g. En la salida que proporciona la función se podrá consultar la siguiente información: N: número de observaciones en cada grupo de estudio. O-E: número ponderado de eventos observados menos el esperado en cada uno de los grupos en estudio. statistic: estadístico asociado al contraste. pvalue: p-valor asociado al contraste que se está realizando. var: matriz de varianza del test. d2L.dB2: segunda derivada de la log-verosimilitud con respecto a β. d2L.dgam2: segunda derivada de la log-verosimilitud con respecto a θ. d2L.dBdgam: segunda derivada de la log-verosimilitud con respecto a βyθ. Además, se especifica el valor de ρusado en los pesos de Fleming-Harrington, el tipo de test aplicado y la descripción de la hipótesis alternativa del contraste. Por último, en el caso de que se utilice el método exact.mc, se puede consultar un intervalo de confianza del p-valor. En las siguientes páginas se mostrarán ejemplos en los que se usará la función en estudio, entre ellos, habrá test con k≥2poblaciones y una prueba de tendencia. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 92 4.2. LIBRERÍA FHTEST Test con 2 poblaciones Se trabajará con el conjunto de datos duser (Gómez et al. (2000)) que contiene 940 observaciones y está disponible en el paquete FHtest. Dicho conjunto de datos contiene los tiempos desde que un usuario usa por primera vez una serie de drogas inyectables hasta la seroconversión del VIH. Los datos proceden de la unidad de desintoxicación del hospital universitario Germans Trias i Pujol de Badalona y las variables en estudio son: left: extremo izquierdo de cada uno de los intervalos en estudio expresado en meses. right: extremo derecho de cada uno de los intervalos en estudio expresado en meses. zper: período del calendario que se denota por un valor del 1 al 4. El primer período va desde el año 1972 al 1980, el segundo desde 1981 a 1985, el tercero desde 1986 hasta 1991 y el último desde 1992 hasta 1997. zgen: género del usuario (0 = Hombre y 1 = Mujer) age: edad a la que comienza el usuario a consumir drogas inyectables. La variable censurada por intervalos es el tiempo (en meses) desde el inicio del consumo de drogas inyectables hasta la seroconversión, síntoma que puede utilizarse como indicador de la infección del sida. A continuación, se cargan los datos y se muestran algunos de ellos: data("duser",package = "FHtest") left right zper zgen age 201 9999 1 0 20 0 211 1 1 12 83 9999 2 0 35 0 106 2 1 17 47 61 2 0 17 19 28 3 0 25 35 9999 3 1 14 18 23 4 0 49 24 9999 4 1 31 Tabla 4.24: Muestra de datos del conjunto duser del paquete FHtest Los tiempos sin censura se caracterizan por la coincidencia de las variables right yleft; los tiempos censurados a la izquierda son aquellos en los que la variable left toma el valor 0 y los tiempos censurados a la derecha son aquellos en los que la variable right toma el valor 9999. Por último, aquellos tiempos que no verifiquen ninguna de las condiciones expuestas se considerarán que están censurados por intervalos. A continuación, se convertirán los valores 9999 a Inf para facilitar el estudio: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 93 4.2. LIBRERÍA FHTEST duser <- duser %> % mutate(right = ifelse(right == 9999, Inf, right)) Una vez preparados los datos, se procede a comparar los tiempos de supervivencia según el género del usuario, es decir, se debe resolver el siguiente contraste: H0:SHombre(t) = SMujer(t)∀t > 0 H1:SHombre(t)=SMujer(t)para algún t > 0 Se tomará como ρel valor 1 para dar más importancia a las diferencias tempranas y se añadirán los parámetros Lin yRin para considerar intervalos de tiempo cerrados. Además, se aumentará el número máximo de iteraciones que utiliza el programa internamente para estimar las funciones de supervivencia por el método EM debido a que surgen problemas de convergencia cuando se mantiene el valor por defecto de 10000 iteraciones. Para ello, se añadirá el argumento icontrol: FHtestics(Surv(left, right, type = "interval2") ~ zgen, data = duser, rho = 1,Rin = TRUE, Lin = TRUE, icontrol = icfitControl(maxit = 50000)) ## ## Two-sample test for interval-censored data ## ## Parameters: rho=1, lambda=0 ## Distribution: score vector approach ## ## Data: Surv(left, right, type = "interval2") by zgen ## ## N O-E ## zgen=0 759 -12.4 ## zgen=1 181 12.4 ## ## Statistic Z= 2.3, p-value= 0.0214 ## Alternative hypothesis: survival functions not equal El p-valor obtenido es 0.0214, por lo que existen evidencias en contra de la hipótesis nula y por tanto, no se puede considerar que los tiempos de supervivencia hasta la seroconversión sean iguales para los hombres y las mujeres. Ante el resultado obtenido, podría plantearse resolver un contraste unilateral como el siguiente: H0:SHombre(t) = SMujer(t)∀t > 0 H1:SHombre(t)≥SMujer(t)∀t > 0 El contraste anterior se resuelve de la siguiente forma: FHtestics(Surv(left, right, type = "interval2") ~ zgen, data = duser, Rin = TRUE, Lin = TRUE, rho = 1,alternative = "decreasing") CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 94 4.2. LIBRERÍA FHTEST ## ## Two-sample test for interval-censored data ## ## Parameters: rho=1, lambda=0 ## Distribution: score vector approach ## ## Data: Surv(left, right, type = "interval2") by zgen ## ## N O-E ## zgen=0 759 -12.4 ## zgen=1 181 12.4 ## ## Statistic Z= 2.3, p-value= 0.0107 ## Alternative hypothesis: decreasing survival functions (group=1 has earlier event times) Dado que el p-valor es significativo, se concluye que se rechaza la hipóteis nula y que por tanto puede considerarse que la curva de supervivencia de los hombres es mayor que la de las mujeres. Se adjunta una gráfica de las estimaciones de las correspondientes funciones de supervivencia de hombres y mujeres consumidores de drogas: ++++++ +++ + ++ + ++ + + + ++ + + +++++++ +++++ + ++ +++++ +++++++++++ + ++++ ++ ++ + + + ++ + ++ +++ + + ++ +++++++ +++++ + ++ ++ + ++ ++ 0.00 0.25 0.50 0.75 1.00 0 5 10 15 20 t (años) S(t) estimada Sexo ++ Hombre Mujer Fuente propia Figura 4.15: Curvas de supervivencia del conjunto de datos duser del paquete FHtest por género Según la gráfica, la probabilidad de seroconversión durante los tres o cuatro primeros meses es mayor entre las mujeres. Después, la función de supervivencia asociada a los hombres se sitúa por encima de la asociada a las mujeres a excepción de que, al comienzo del quinto año, la curva que se sitúa por encima es la asociada a las mujeres. Es destacable también que después de 13 años, no se observa más seroconversión entre los consumidores femeninos. CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 95 4.2. LIBRERÍA FHTEST Test con kpoblaciones Para ilustrar la prueba con k > 2poblaciones o grupos se usará el conjunto de datos duser que se ha utilizado en el apartado anterior. En este caso, la variable de agrupación que se tomará es zper y el estudio se centrará en los consumidores de drogas masculinos. En concreto, se resolverá el siguiente contraste: H0:SP1,H(t) = SP2,H (t) = SP3,H =SP4,H ∀t > 0 H1:SP1,H(t)=SP2,H (t)=SP3,H =SP4,H para algún t > 0 A la hora de aplicar dicho contraste en R, debe recordarse que la variable de agrupación no puede ser numérica y para trabajar solo con los datos de usuarios masculinos se añade el argumento subset junto con la condición zgen == 0: FHtestics(Surv(left, right, type = "interval2") ~ as.character(zper), data = duser, subset = (zgen == 0), Lin = TRUE, Rin = TRUE, rho = 5) ## ## K-sample test for interval-censored data ## ## Parameters: rho=5, lambda=0 ## Distribution: score vector approach ## ## Data: Surv(left, right, type = "interval2") by as.character(zper) ## ## N O-E ## as.character(zper)=1 146 -0.524 ## as.character(zper)=2 300 2.632 ## as.character(zper)=3 240 -0.997 ## as.character(zper)=4 73 -1.111 ## ## Chisq= 1.8 on 3 degrees of freedom, p-value= 0.605 ## Alternative hypothesis: survival functions not equal Tomando como ρel valor 5, se obtiene un p-valor de 0.605, por lo que no existen evidencias suficientes en contra de H0. De esta forma, puede admitirse que las curvas de supervivencia asociadas a cada uno de los períodos son iguales, sobre todo en las fases tempranas del estudio. A continuación, se adjunta una gráfica de las estimaciones de las funciones de supervivencia asociadas a cada uno de los períodos en estudio para la población masculina: CAPÍTULO 4. COMPARACIÓN DE CURVAS DE SUPERVIVENCIA EN R 96 Bibliografía M.C. Jaramillo y C.M. Lopera. Estudio del efecto de la imputación de fallas en la estimación de la curva de supervivencia bajo censura a intervalo. Ciencia en Desarrollo, 8(1):21–28, 2017. H.W. Jones. John Graunt and his Bills of Mortality. Bulletin of the Medical Library Association, 33(1):3–4, 1945. J.D. Kalbfleisch y R.L. Prentice. The Statistical Analysis of Failure Time Data. John Wiley & Sons, Second edition, 2002. E.L. Kaplan y P. Meier. Nonparametric Estimation from Incomplete Observations. Journal of the American Statistical Association, 53(282):457–481, 1958. doi: https: //doi.org/gdz3gq. C. Kartsonaki. Survival Analysis. Diagnostic Histopathology, 22(7):263–270, 2016. ISSN 1756-2317. doi: https://doi.org/gf7fph. Mini-Symposium: Medical Statistics. B.R. Kirkwood y J.A.C. Sterne. Essential Medical Statistics. John Wiley & Sons, 2010. M. Klass, H. Teicher, et al. The Central Limit Theorem for Exchangeable Random Variables without Moments. The Annals of Probability, 15(1):138–153, 1987. J.P. Klein y M.L. Moeschberger. Survival Analysis Techniques for Censored and Truncated Data. Springer, 1997. J.P. Klein y M.L. Moeschberger. KMsurv: Data sets from Klein and Moeschberger (1997), Survival Analysis. R package version 0.1-5, 2012. Disponible en: https: //CRAN.R-project.org/package=KMsurv. D.G. Kleinbaum y M. Klein. Survival Analysis: A Self-Learning Text. Springer, Third edition, 2012. A. Knezevic y S. Patil. "Combination Weighted Log-rank Tests for Survival Analysis with Non-Proportional Hazards". Disponible en https://www.sas.com/content/dam/ SAS/support/en/sas-global-forum-proceedings/2020/5062-2020.pdf, 2020. T.L. Lai y H. Xing. Statistical Models and Methods for Financial Markets. Springer, 2008. J.A. Laurie, C.G. Moertel, T.R. Fleming, H.S. Wieand, J.E. Leigh, J. Rubin, G.W. Mc- Cormack, J.B. Gerstner, J.E. Krook, y J. Malliard. Surgical Adjuvant Therapy of Large-Bowel Carcinoma: an Evaluation of Levamisole and the Combination of Levamisole and Fluorouracil. The North Central Cancer Treatment Group and the Mayo Clinic. Journal of Clinical Oncology, 7(10):1447–1456, 1989. E.T. Lee y J. Wang. Statistical Methods for Survival Data Analysis, volume 476. John Wiley & Sons, 2003. C.L. Loprinzi, J.A. Laurie, H.S. Wieand, J.E. Krook, P.J. Novotny, J.W. Kugler, J. Bartel, M. Law, M. Bateman, y N.E. Klatt. Prospective Evaluation of Prognostic Variables from Patient-Completed Questionnaires. North Central Cancer Treatment Group. Journal of Clinical Oncology, 12(3):601–607, 1994. 103 Bibliografía P.L. Luque. Escribir un Trabajo Fin de Estudios con R Markdown. Disponible en http://destio.us.es/calvo, 2017. P.L. Luque. Cómo crear Tablas de información en R Markdown. Disponible en http: //destio.us.es/calvo, 2019. M. Maathuis. MLEcens: Computation of the MLE for bivariate (interval) censored data. R package version 0.1-4, 2013. Disponible en: https://CRAN.R-project.org/package= MLEcens. N. Mantel, N.R. Bohidar, y J.L. Ciminera. Mantel-Haenszel Analyses of Litter-Matched Time-to-Response Data, with Modifications for Recovery of Interlitter Information. Cancer Research, 37(11):3863–3868, 1977. C.M. Martínez. Bonferroni’s Method to Compare Survival Curves with Recurrent Events. Ingeniería Industrial. Actualidad y Nuevas Tendencias, 3(10):105–114, 2013. D.F. Moore. Applied Survival Analysis Using R. Springer, 2016. W. Nelson. Hazard Plotting for Incomplete Failure Data. Journal of Quality Technology, 1(1):27–52, 1969. doi: https://doi.org/gfgm. R. Oller y K. Langohr. FHtest: An R Package for the Comparison of Survival Curves with Censored Data. Journal of Statistical Software, 81(15):1–25, 2017. doi: https: //doi.org/fz35. R. Oller y K. Langohr. FHtest: Tests for Right and Interval-Censored Survival Data Based on the Fleming-Harrington Class. R package version 1.5., Barcelona, España, 2020. Disponible en: https://CRAN.R-project.org/package=FHtest. M. Pešta. "Survival Data Analysis". Disponible en https://www2.karlin.mff.cuni.cz/ ~pesta/NMFM404/survival.html, 2021. R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2020. Disponible en: https://www. R-project.org/. J. Ramírez-Montoya. Comparación de estimadores no paramétricos frente a los paramétricos para la función de confiabilidad. Ingeniería Investigación y Desarrollo, 15(2): 15–24, 2015. P. Rebasa. Conceptos básicos del análisis de supervivencia. Cirugía española, 78(4): 222–230, 2005. ISSN 0009-739X. doi: https://doi.org/c87m3p. B. Ripley, B. Venables, D.M. Bates, K. Hornik, A. Gebhardt, y D. Firth. MASS: Support Functions and Datasets for Venables and Ripley’s MASS. R package version 7.3-54, 2021. Disponible en: https://CRAN.R-project.org/package=MASS. RStudio Team. RStudio: Integrated Development Environment for R. RStudio, Inc., Boston, MA, 2015. Disponible en: http://www.rstudio.com/. P. Sedgwick. Kaplan-Meier Survival Analysis: Types of Censored Observations. Bmj, 347, 2013. doi: https://doi.org/gh7rp5. 104 Bibliografía A.N. Shiryaev. "Martingale. Encyclopedia of Mathematics.". Disponible en http: //encyclopediaofmath.org/index.php?title=Martingale&oldid=49256. NCSS Statistical Software. "Kaplan-Meier Curves (Logrank Tests)". Disponible en https: //www.ncss.com/software/ncss/ncss-documentation/. M. Tableman y J.S. Kim. Survival Analysis Using S: Analysis of Time-to-Event Data. CRC press, 2003. T.M. Therneau. A Package for Survival Analysis in R. R package version 3.2-11, 2021. Disponible en: https://CRAN.R-project.org/package=survival. A. Tillander, A. Karlsson, J. Zetterqvist, P. Strom, B. Delcoigne, y M. Clements. biostat3: Utility Functions, Datasets and Extended Examples for Survival Analysis. R package version 0.1.5, 2020. Disponible en: https://CRAN.R-project.org/package=biostat3. B.W. Turnbull. Nonparametric Estimation of a Survivorship Function with Doubly Censored Data. Journal of the American Statistical Association, 69(345):169–173, 1974. doi: https://doi.org/gg2nc6. B.W. Turnbull. The Empirical Distribution Function with Arbitrarily Grouped, Censored and Truncated Data. Journal of the Royal Statistical Society: Series B (Methodological), 38(3):290–295, 1976. doi: https://doi.org/ggqzts. P.A. Volberding, S.W. Lagakos, J.M. Grimes, D.S. Stein, J. Rooney, T. Meng, M.A. Fischl, A.C. Collier, J.P. Phair, M.S. Hirsch, et al. A Comparison of Immediate with Deferred Zidovudine Therapy for Asymptomatic HIV-Infected Adults with CD4 Cell Counts of 500 or More per Cubic Millimeter. New England Journal of Medicine, 333 (7):401–407, 1995. doi: https://doi.org/dk6b89. H. Wickham y Rstudio. tidyverse: Easily Install and Load the ’Tidyverse’. R package version 1.3.1, 2021. Disponible en: https://CRAN.R-project.org/package=tidyverse. H. Wickham, W. Chang, L. Henry, T.L. Pedersen, K. Takahashi, C. Wilke, K. Woo, H. Yutani, y D. Dunnington. ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. R package version 3.3.3, 2020. Disponible en: https://CRAN. R-project.org/package=ggplot2. 105