Modelo de Cox con covariantes dependientes del tiempo

FACULTAD DE MATEMÁTICAS ESTADÍSTICA E INVESTIGACIÓN OPERATIVA Trabajo Fin de Máster Modelo de Cox con covariantes dependientes del tiempo Elena de las Heras Jiménez Máster en Matemáticas Dirigido por: D. Juan Manuel Muñoz Pichardo Sevilla, noviembre de 2023 Índice general Resumen ........................................ iii Abstract ........................................ iv Índice de Figuras ................................... v Índice de Tablas ....................................vii 1. Introducción 1 1.1. Breve referencia histórica ............................ 1 1.2. Estructura del trabajo ............................. 3 2. Modelo de Cox 4 2.1. Formulación del modelo ............................ 4 2.2. Características del modelo ........................... 5 2.3. Estimación de parámetros ........................... 6 2.4. Contrastes de hipótesis ............................ 9 2.5. Diagnosis del modelo .............................. 10 2.5.1. Procedimiento gráfico ......................... 11 2.5.2. Contraste de bondad de ajuste (GOF) ................ 12 2.5.3. Uso de variables dependientes del tiempo .............. 13 2.6. Evaluación del modelo ............................. 14 2.7. Variantes del modelo .............................. 15 2.8. Técnicas de selección de variables ....................... 17 3. Modelo de Cox con covariantes dependientes del tiempo 18 3.1. Tipos de covariantes dependientes del tiempo ................ 18 3.2. Aplicaciones de las covariantes dependientes del tiempo ........... 19 3.3. Formulación del modelo ............................ 20 3.4. Características del modelo ........................... 20 3.5. Estimación de parámetros ........................... 21 3.5.1. Propiedades asintóticas ......................... 25 i 3.6. Contrastes de hipótesis ............................. 26 3.7. Evaluación del modelo ............................. 27 3.8. Técnicas de selección de variables ....................... 28 4. Modelo de Cox en R 29 4.1. Introducción ................................... 29 4.2. Modelo de Cox ................................. 30 4.2.1. Diagnosis del modelo de Cox ...................... 34 4.2.2. Evaluación del modelo de Cox ..................... 36 4.2.3. Selección de variables del modelo de Cox ............... 37 4.2.4. Otras consideraciones del modelo de Cox ............... 39 4.3. Aplicaciones ................................... 40 4.3.1. Modelo de Cox con covariantes independientes del tiempo ..... 41 4.3.1.1. Diagnosis del modelo .................... 44 4.3.1.2. Estimación de parámetros .................. 51 4.3.1.3. Contrastes de hipótesis ................... 55 4.3.1.4. Evaluación del modelo .................... 58 4.3.1.5. Selección de variables del modelo .............. 59 4.3.1.6. Otras consideraciones del modelo .............. 64 4.3.2. Modelo de Cox con covariantes dependientes del tiempo ....... 70 4.3.2.1. Estimación de parámetros .................. 74 4.3.2.2. Contrastes de hipótesis ................... 76 4.3.2.3. Evaluación del modelo .................... 78 4.3.2.4. Selección de variables del modelo .............. 82 4.3.2.5. Otras consideraciones del modelo .............. 87 Conclusión 90 Bibliografía 95 ii Resumen El análisis de supervivencia es una metodología de análisis estadístico muy útil en cualquier ámbito de investigación (industria, medicina, economía, biología, demografía...) y tiene como objetivo analizar o modelizar el tiempo que tardan en ocurrir uno o más eventos de interés. El modelo más usado en este campo es el modelo semiparamétrico de regresión de Cox, conocido también como modelo de riesgos proporcionales. Una extensión del mismo con grandes aplicaciones prácticas en la ciencia, surge ante la presencia de covariantes dependientes del tiempo. De esta forma, el propósito de este trabajo es el estudio de este último modelo extendido de Cox que incluirá la formulación general para permitir variables dependientes del tiempo, una discusión de las características del modelo y la inferencia estadística sobre el mismo. Además, con objeto de ilustrar la utilidad de este modelo, se usará R y se estudiará dicho modelo sobre datos reales. Palabras claves: análisis de supervivencia, modelo de regresión de Cox, función de riesgo, hipótesis de riesgos proporcionales, covariante dependiente del tiempo. iii Abstract Survival analysis is a very useful statistical analysis methodology in any field of research (industry, medicine, economics, biology, demography. ..) and aims to analyze or model the time it takes for one or more events of interest to occur. The most widely used model in this field is the semi-parametric Cox regression model, also known as the proportional hazards model. An extension of this model with great practical applications in science arises in the presence of time-dependent covariates. Thus, the purpose of this paper is the study of the latter extended Cox model which will include the general formulation to allow for time-dependent variables, a discussion of the characteristics of the model and statistical inference on it. In addition, in order to illustrate the usefulness of this model, R will be used and the model will be studied on real data. Keywords: survival analysis, Cox regression model, hazard function, proportional hazards hypothesis, time-dependent covariates. iv Índice de figuras 4.1. Gráficas de los residuos de Schoenfeld con la orden plot(cox.zph()) ..... 47 4.2. Gráficas de los residuos de Schoenfeld con la orden ggcoxdiagnostics() de la librería survminer .............................. 48 4.3. Comparación de las curvas log-log de supervivencia ............. 49 4.4. Gráficas obtenidas con la orden profLik() ................... 55 4.5. Captura de la salida summary() para el estudio del contraste de hipótesis conjunto del modelo de Cox .......................... 56 4.6. Resumen gráfico del análisis del modelo de Cox con la función ggforest() . 58 4.7. Evolución del valor de log p λ q en función de la verosimilitud parcial ..... 61 4.8. Evolución del valor de los coeficientes frente a la norma L1 en el método Lasso ....................................... 64 4.9. Estimación de la función de supervivencia para el modelo creado ...... 68 4.10. Estimación de la función de supervivencia para el modelo creado ...... 69 4.11. Estimación de la función de supervivencia para la variable site ....... 69 4.12. Resumen gráfico del análisis del modelo de Cox con covariantes dependientes del tiempo realizado con la función ggforest() ............... 78 4.13. Gráficas de los residuos escalados de Schoenfeld con la orden plot(cox.zph()) para el modelo de Cox creado con covariantes dependientes del tiempo . . 81 4.14. Estudio gráfico y analítico de los residuos de Schoenfeld con la orden ggcoxzph() de la librería survminer ....................... 81 4.15. Evolución del valor de log p λ q en función de la verosimilitud parcial del modelo de Cox con covariantes dependientes del tiempo ........... 84 4.16. Evolución del valor de los coeficientes frente a la norma L1 en el método Lasso para un modelo de Cox con covariantes dependientes del tiempo . . 87 4.17. Estimación de la función de supervivencia para el modelo creado con covariantes dependientes del tiempo ........................ 89 v Índice de tablas 4.1. Muestra de datos del fichero uissurv...................... 42 4.2. Estudio de la hipótesis de riesgos proporcionales del modelo de Cox . . . . 45 4.3. Residuos escalados de Schoenfeld asociados al modelo de Cox ........ 46 4.4. Estudio de la hipótesis de riesgos proporcionales del nuevo modelo ..... 51 4.5. Estimación de los coeficientes asociados a cada una de las variables explicativas del modelo ............................... 52 4.6. Estimación de la exponencial de los coeficientes asociados a cada una de las variables explicativas del modelo ...................... 53 4.7. Estimación de los coeficientes asociados a cada una de las variables explicativas del modelo junto a su intervalo de confianza al 95 % ......... 54 4.8. Contrastes de hipótesis individuales asociados al modelo de Cox ...... 56 4.9. Test de razón de verosimilitud del modelo de Cox .............. 57 4.10. Test de Wald del modelo de Cox ........................ 57 4.11. Test score del modelo de Cox ......................... 57 4.12. Estimación de los coeficientes usando λm´ın en el método Lasso ....... 62 4.13. Estimación de los coeficientes usando λ1se en el método Lasso ....... 63 4.14. Predicción sobre un nuevo conjunto de datos ................. 65 4.15. Estimación de la función de supervivencia a partir del modelo creado . . . 66 4.16. Estimación de la función de supervivencia para el programa ’A’ ...... 67 4.17. Estimación de la función de supervivencia para el programa ’B’ ...... 68 4.18. Muestra de datos del conjunto creado pbc_unif ............... 72 4.19. Datos heart del paquete survival ........................ 73 4.20. Estimación de los coeficientes asociados a cada una de las variables explicativas del modelo con covariantes dependientes del tiempo ......... 74 4.21. Estimación de la exponencial de los coeficientes asociados a cada una de las variables explicativas del modelo ...................... 75 4.22. Estimación de los coeficientes asociados a cada una de las variables explicativas del modelo con covariantes dependientes del tiempo junto a su intervalo de confianza al 95 % ......................... 75 vi 4.23. Contrastes de hipótesis individuales asociados al modelo de Cox con covariantes dependientes del tiempo ........................ 76 4.24. Contraste de hipótesis conjunto asociado al modelo de Cox con covariantes dependientes del tiempo ............................ 77 4.25. Estudio de la hipótesis de riesgos proporcionales con covariantes dependientes del tiempo para el modelo de Cox ................... 79 4.26. Residuos escalados de Schoenfeld asociados al modelo de Cox con covariantes dependientes del tiempo .......................... 80 4.27. Ajuste de los parámetros en la selección de variables en el modelo de Cox con covariantes dependientes del tiempo ................... 84 4.28. Estimación de los coeficientes usando λm´ın en el método Lasso para un modelo de Cox con covariantes dependientes del tiempo ........... 85 4.29. Estimación de los coeficientes usando λ1se en el método Lasso para un modelo de Cox con covariantes dependientes del tiempo ........... 86 4.30. Estimación de la función de supervivencia a partir del modelo creado con covariantes dependientes del tiempo ...................... 88 vii Capítulo 1 Introducción En estadística, el análisis de supervivencia se plantea como principal objetivo analizar y modelizar el tiempo transcurrido hasta que se produce un evento o suceso de interés. Dicho análisis tiene utilidad en muchos ámbitos de investigación como pueden ser la medicina, la ingeniería o la sociología. En la medicina, por ejemplo, es habitual el estudio de la variable tiempo hasta que se produce la recaída de un paciente tras suministrarle un tratamiento médico específico para una enfermedad determinada. En el análisis de supervivencia, la variable objetivo es el tiempo de vida y es imprescindible que se conozca el momento en el que se inicia el estudio de dicha variable en cada caso para así poder realizar un estudio detallado y fiable. Nótese que en algunas ocasiones, puede suceder que los individuos abandonen el estudio antes de que se produzca el evento en cuestión o que se incorporen después del inicio del estudio, de modo que en dichos casos, la información conocida del individuo no es completa, sino parcial. Esta peculiaridad caracteriza al análisis de supervivencia y los datos que la presentan, se conocen como datos censurados. Existen diferentes técnicas estadísticas para modelar la variable objetivo o de interés, sin embargo, este trabajo se centrará en el estudio del modelo de regresión de Cox o de riesgos proporcionales para predecir el tiempo de fallo o probabilidad de supervivencia a través de diferentes variables que pueden ser influyentes en la variable objetivo y que incluso pueden llegar a depender del tiempo. Además, permite identificar y evaluar la relación que existe entre el conjunto de variables explicativas y la función de riesgo en estudio. Esta técnica multivariante del análisis de supervivencia puede ser útil, por ejemplo, para observar o detectar las variables que influyen a la hora de delinquir de nuevo un preso que ha sido puesto en libertad recientemente (véase Berk et al. (1980)). Otra técnica usual como la comparación de curvas de supervivencia y la presentación de conceptos básicos del análisis de supervivencia que podrían ser de ayuda para el presente trabajo, se pueden consultar en el Trabajo Fin de Grado Comparación de curvas de supervivencia con datos censurados de de las Heras (2021). 1.1. Breve referencia histórica El origen del análisis de supervivencia se sitúa en el siglo xvii con la creación de tablas de mortalidad que se usaban para realizar estudios demográficos. No obstante, el 1 2.3. ESTIMACIÓN DE PARÁMETROS instante ti, de modo que se verifica que card t U p i q u  m p i q ;r p i q  card t R p ti qu siendo R p ti q el conjunto de individuos en riesgo en el instante ti(incluye a los casos censurados con censura mayor que ti) y z p i q  ¸ j P U p i q x p j q . En base a la notación anterior, los dos primeros modelos citados se usan cuando el tiempo es continuo y su complejidad computacional es menor que la del método exacto. A continuación, se enuncian cada uno de los métodos considerados: Método de Breslow: L  B p β q  D ¹ i  1 exp p z 1 p i q β q ¸ j P R p ti q exp p x 1 p j q β q  m p i q , donde el numerador tiene en cuenta las covariantes que fallan en el instante ti, mientras que el denominador, considera la suma de la función de riesgo parcial del conjunto de riesgo R p ti q elevada a la multiplicidad de cada uno de los instantes ti considerados. Método de Efron: L  E p β q  D ¹ i  1 exp p z 1 p i q β q m p i q ¹ j  1 ¸ k P R p ti q exp p x 1 p k q β q  j  1 m p i q ¸ k P U p i q exp p x 1 p k q β q  , aproximación más precisa que los métodos exacto y de Breslow y se usa preferiblemente cuando la muestra en estudio es pequeña. Por último, cuando los tiempos de supervivencia se observan en tiempo discreto, Cox propone el método Exact Partial Likelihood que se define de la siguiente manera: L  C p β q  D ¹ i  1 exp p z 1 p i q β q ¸ u P U p i q exp #¸ j P u x p j q  1 β + . Tal y como se puede consultar en este apartado, a la hora de estimar los parámetros, solo se ha tenido en cuenta el orden de los tiempos de fallo y censura. Como consecuencia de ello, Harrell (2001) concluye que el modelo de Cox es menos sensible a los outliers del estudio que cualquier modelo parámetrico de análisis de supervivencia. Esta característica favorece el uso del modelo de Cox en problemas de regresión en el ámbito del análisis de supervivencia. Para estudiar en más detalle estas técnicas, consultar Ramírez et al. (2017). No obstante, la estimación de los parámetros del modelo, también se pueden realizar mediante métodos de remuestreo como las técnicas Bootstrap7o Jackknife8. Dichas 7Técnica de remuestreo que calcula el estimador de un parámetro de interés usando muestras con reemplazamiento a partir de la muestra original y con el mismo tamaño que esta. Se utiliza cuando la distribución de los datos en estudio es difícil de obtener analíticamente o cuando las hipótesis clásicas del problema no se cumplen. 8Técnica de remuestreo cuyo objetivo es el cálculo del estimador de un parámetro de interés usando la media de las estimaciones obtenidas al eliminar, para cada una de ellas, una observación de la muestra original en estudio. Se utiliza para la estimación del sesgo y la varianza de diferentes estimadores. CAPÍTULO 2. MODELO DE COX 8 2.4. CONTRASTES DE HIPÓTESIS técnicas poseen la ventaja de que reducen la necesidad de asumir determinados modelos probabilísticos para las observaciones. Sin embargo, tiene la desventaja de que, en general, requieren un número elevado de cálculos. 2.4. Contrastes de hipótesis Como es habitual en los modelos de regresión cuando se estiman los parámetros, se realizan test sobre los parámetros estimados para estudiar si dichos parámetros pueden considerarse significativos o no. Para ello, considérese que, asintóticamente, los parámetros estimados p βsiguen una distribución normal de media βy matriz de varianzas y covarianzas Σ  Φ  1 p β q de dimensión p  p, i.e., p βa  Np p β, Φ  1 p β qq . A su vez, la media y la matriz de varianzas y covarianzas de la distribución normal en cuestión, pueden ser estimadas de la siguiente manera: p β  β y p Σ  Φ  1 p p β q , donde: Φ p β qB 2l  p β q B β B β 1 . El primer contraste de hipótesis que se estudiará es el contraste de hipótesis individual: " H0:βj  0 H1:βj  0, para cada j  1, . . . , p. Para estudiar la significación de la j-ésima covariante del modelo de Cox considerado, se puede utilizar el estadístico de Wald que se presenta a continuación: Z  p βj by V ar p p βj q , donde y V ar p p βj q es el j-ésimo elemento de la diagonal de p Σ. Bajo hipótesis nula, Zse distribuye asintóticamente según una distribución normal estándar. En consecuencia, el intervalo de confianza para cada uno de los coeficientes p βjse obtiene de la siguiente forma: p βj  zα { 2 by V ar p p βj q , dado que bajo la hipótesis nula, el estadístico en estudio Zsigue una distribución normal estándar. Asimismo, es posible estudiar también el contraste de hipótesis conjunto: " H0:β  p β10, . . . , βp0 q 1 H1:βj  βj0para algún j  1, . . . , p. Por su parte, este último contraste, se puede resolver usando tres métodos o test diferentes: el test de razón de verosimilitud, el test de Wald o el test score. Cada uno de ellos tiene unas propiedades distintas y se basan en estadísticos diferentes. A continuación, se enuncian cada uno de los test citados y recogidos en Cox y Oakes (1984): CAPÍTULO 2. MODELO DE COX 9 2.5. DIAGNOSIS DEL MODELO El test de Wald: se basa en la media estimada ( p β) que sigue asintóticamente una distribución normal y el estadístico considerado es: WW ald  p p β  β0 q 1 Φ p p β qp p β  β0 q . Dicho estadístico, bajo la hipótesis nula del contraste, sigue una distribución χ2con pgrados de libertad (WW ald  χ2 p). El test de razón de verosimilud: compara la función de verosimilitud parcial evaluada en la estimación del parámetro p β, i.e., L  p p β q , con esa misma función evaluada en el caso en el que la hipótesis nula del contraste es cierta, L  p β0 q . El estadístico asociado a este test es el siguiente: WLR  2 p log L  p p β q  log L  p β0 qq . Bajo hipótesis nula, el estadístico enunciado sigue una distribución χ2con pgrados de libertad (WLR  χ2 p). Test score: usa las derivadas del logaritmo de la función de verosimilitud parcial evaluada en la hipótesis nula del contraste y el estadístico asociado se enuncia de la siguiente forma: Wscore   B l  p β0 q B β  1  B l  2 p β0 q B β B β 1   1 B l  p β0 q B β. Bajo hipótesis nula, el estadístico enunciado sigue una distribución χ2con pgrados de libertad (Wscore  χ2 p). Tanto el test de razón de verosimilitud como el test score son contrastes invariantes ante diferentes parametrizaciones, sin embargo, el test de Wald tiene una interpretación más directa que el resto de los test en estudio. Además, el test de razón de verosimilitud converge de manera más rápida a una distribución normal mientras que el test score es más rápido computacionalmente cuando el test a contrastar tiene varios parámetros. No obstante, el test que usualmente se utiliza más es el test de razón de verosimilitud. 2.5. Diagnosis del modelo La hipótesis fundamental que se debe verificar en el modelo de regresión de Cox es la hipótesis de riesgos proporcionales. Esta hipótesis debe verificarse tanto para cada una de las variables como para el conjunto de variables involucradas en el modelo. En el caso de que no se verifique el contraste para alguna de las variables, se deberá considerar un nuevo modelo de Cox sin dicha variable o realizar alguna transformación para que la variable en cuestión cumpla la hipótesis en estudio. Asimismo, en el caso de que todas las variables cumplan la hipótesis de riesgos proporcionales pero no se cumpla de manera global, será necesario replantear el modelo o realizar alguna transformación para que así se verifique la hipótesis fundamental del modelo. Tal y como se puede consultar en Kleinbaum y Klein (2012), existen diferentes procedimientos para estudiar el cumplimiento de la hipótesis de riesgos proporcionales: un CAPÍTULO 2. MODELO DE COX 10 2.5. DIAGNOSIS DEL MODELO procedimiento gráfico, un contraste de bondad de ajuste y un procedimiento que involucra el uso de variables dependientes del tiempo. Se recomienda el uso de al menos dos de los procedimientos que se detallan a continuación a la hora de evaluar si la hipótesis de riesgos proporcionales es válida para el modelo en estudio, pues cada una de las técnicas propuestas tienen sus ventajas e inconvenientes. 2.5.1. Procedimiento gráfico Como procedimientos gráficos, se encuentran dos posibles técnicas que permiten evaluar si la hipótesis de riesgos proporcionales se cumple o no en el modelo. La primera de ellas consiste en la comparación de las curvas log-log de supervivencia y la segunda, en la comparación de las curvas de supervivencia observadas frente a las esperadas. Curvas log-log de supervivencia: se trata de una transformación logarítmica de la curva de supervivencia estimada, en concreto, consiste en tomar dos veces el logaritmo natural de la curva en cuestión. Matemáticamente, la curva log-log de supervivencia se escribe como  ln p ln p S q y su signo es desconocido. Normalmente, este tipo de gráficos utilizan como estimación de la curva de supervivencia la estimación de Kaplan-Meier. Dado que esta estimación de la curva de supervivencia es escalonada, la curva log-log también lo es. Véase de manera genérica, la construcción de este tipo de curvas: sea la función de supervivencia S p t | X q tal y como se define a continuación: S p t | X qr S0 p t qs exp p X 1 β q  r S0 p t qs exp  ° p j  1βjXj  . Al tomar un logaritmo natural en dicha función, se obtiene que: ln S p t | X q  exp  p ¸ j  1 βjXj   ln S0 p t q . Como es bien sabido, se cumple que 0 ¤ S p t | X q ¤ 1, por lo que al tomar un segundo logaritmo, se tendrá que: ln r ln S p t | X qs  ln   exp  p ¸ j  1 βjXj   ln S0 p t q   ln  exp  p ¸ j  1 βjXj   ln r ln S0 p t qs  p ¸ j  1 βjXj  ln r ln S0 p t qs , (2.3) de modo que el rango de la curva  ln p ln q está comprendido entre 8 y 8 . Tal y como se puede consultar en la expresión (2.3), la curva log-log de supervivencia se puede expresar como la suma de dos términos, por un lado, un predictor lineal y por otro, la transformación log-log de la función de supervivencia base. En el caso de que un modelo de Cox sea adecuado para un conjunto de predictores en estudio, las gráficas empíricas esperadas de las curvas log-log de supervivencia para diferentes individuos de la muestra, deben ser aproximadamente paralelas. CAPÍTULO 2. MODELO DE COX 11 2.5. DIAGNOSIS DEL MODELO Curvas de supervivencia observadas frente a las esperadas: se basa en la estratificación de los datos por categorías del predictor a evaluar y posteriormente, la estimación de la curva de supervivencia asociada a cada una de las categorías consideradas. Si para cada categoría del predictor evaluado se tiene que los gráficos observados y esperados están próximos entre sí y no discrepan mucho, entonces se puede concluir que la hipótesis de riesgos proporcionales se cumple. En contraposición, si los gráficos de una o más categorías discrepan entre sí, se llegará a la conclusión de que se está incumpliendo la hipótesis en estudio. 2.5.2. Contraste de bondad de ajuste (GOF) Se trata de una opción más objetiva que el procedimiento gráfico, dado que con el contraste que se realiza para evaluar si la hipótesis de riesgos proporcionales es cierta o no, se obtiene un estadístico de prueba y un p-valor. En esta sección, se presenta el test de Harrell y Lee (1986) que se trata de una variación de la prueba propuesta originalmente por Schoenfeld (1982) y la cual se basa en los residuos denominados de Schoenfeld asociados a cada uno de los predictores involucrados en el modelo de Cox en estudio y a cada uno de los sujetos que tiene un evento. Dichos residuos se definen como los valores observados de las covariantes o variables explicativas menos los valores esperados en cada instante de fallo y para cada uno de los casos no censurados. Además, se verifica que este tipo de residuos para cada una de las variables predictoras suma cero y que un valor positivo de dicho residuo significa que el valor de la variable asociada es más alto de lo esperado en ese instante de fallo. A continuación, se enuncia el cálculo de los residuos de Schoenfeld para cada variable predictora o explicativa y caso no censurado: r p s q ij  xij  ° k P R p ti q exp p x 1 k p β q xkj ° k P R p ti q exp p x 1 k p β q xij  p s x p wi q j, siendo R p ti q el conjunto de todos los individuos o sujetos muestrales que se encuentran en riesgo justo antes del instante ti, tal y como se había especificado al comienzo del Apartado 2.3. Se obtiene así una serie de residuos para cada una de las pvariables predictoras (j  1, . . . , p). Si estos residuos mantienen un patrón aleatorio o no sistemático, se puede concluir que existen evidencias de que el efecto de la covariable no cambia respecto del tiempo, hecho que presupone el modelo de Cox. En contraposición, si existiese algún tipo de patrón sistemático en los residuos, se concluirá que el efecto de la covariable cambia a lo largo del tiempo. De esta forma, se considerará que la hipótesis de riesgos proporcionales es cierta cuando los residuos no muestran tendencias temporales y en la gráfica de los residuos frente al tiempo la pendiente sea nula. El test en cuestión puede realizarse en tres pasos y estudia si el efecto de cada una de las variables es constante en el tiempo: 1. Creación o ejecución del modelo de Cox y obtención de los residuos de Schoenfeld asociados a cada uno de los predictores, i.e., para cada j  1, . . . , p, se consideran t r p s q ij u icon i  1, . . . , n. CAPÍTULO 2. MODELO DE COX 12 2.5. DIAGNOSIS DEL MODELO 2. Creación de una variable que es capaz de clasificar en orden los distintos sujetos en función de los tiempos de fallo: U p i q  rango de tien la colección de tiempos de fallo. 3. Estudio de la correlación entre las variables creadas en los dos primeros pasos y se considera como hipótesis nula que la correlación entre los residuos de Schoenfeld y el tiempo de fallo censurado es nulo (H0:ρk  0). En el caso de que se rechace la hipótesis nula, se concluirá también que la hipótesis de riesgos proporcionales no se cumple. Dicho test se realiza para cada una de las variables por separado, así como para el conjunto de variables involucradas en el modelo de Cox creado. Además, el estadístico de prueba asociado a este contraste de hipótesis se define de la siguiente forma: Z  ρk ? nu  2 a 1  ρ2 k , donde ρkes la correlación entre los residuos de Schoenfeld y el orden de los tiempos de fallo; y nues el número total de observaciones no censuradas en la muestra estudiada. Dicho estadístico, bajo la hipótesis nula del contraste, sigue una distribución χ2con un grado de libertad cuando se estudia la hipótesis de riesgos proporcionales de manera individual para cada una de las variables involucradas en el modelo y sgrados de libertad cuando el estudio de la hipótesis de riesgos proporcionales se hace de manera conjunta sobre todas las variables consideradas en el modelo de Cox creado, siendo sel número de variables involucradas en el modelo bajo estudio. Nótese que el p-valor obtenido tras realizar este test de bondad de ajuste puede verse influenciado por el tamaño de la muestra con la que se esté trabajando. De esta forma, si la muestra es pequeña, una violación grave de la hipótesis nula puede no ser estadísticamente significativa. Por el contrario, una violación leve de la hipótesis nula puede ser muy significativa. A pesar de este inconveniente, el uso de una prueba estadística ofrece un enfoque más objetivo para evaluar la hipótesis de riesgos proporcionales que el uso de un enfoque gráfico. No obstante, en ocasiones puede ser de utilidad usar un enfoque gráfico dado que permite detectar tipos específicos de desviaciones de la hipótesis de riesgos proporcionales. Se recomienda así el uso tanto de una prueba estadística como de un procedimiento gráfico para concluir de manera definitiva una decisión sobre la violación o no de la hipótesis de riesgos proporcionales. 2.5.3. Uso de variables dependientes del tiempo En ocasiones, puede ser útil el uso de variables dependientes del tiempo para evaluar la hipótesis de riesgos proporcionales de un modelo de Cox en el que las variables predictoras son independientes del tiempo. Cuando se hace uso de este procedimiento, es necesario considerar el modelo de Cox ampliado en el que se tienen en cuenta las interacciones entre las variables predictoras en estudio y el tiempo, es decir, un modelo de Cox en el que haya involucrados productos de las variables explicativas con alguna función del tiempo (por ejemplo, t,log t...). CAPÍTULO 2. MODELO DE COX 13 2.6. EVALUACIÓN DEL MODELO Una vez considerado el modelo de Cox ampliado, el estudio de si la hipótesis de riesgos proporcionales es cierta o no para una variable predictora jconcreta, se basa en la evaluación de la significación del coeficiente asociado al término producto del modelo de la variable en estudio. En este caso, el modelo se definirá de la siguiente forma: h p t | X q  h0 p t q exp p βX  γ p Xj  g p t qqq , donde g p t q es la función del tiempo asociada a la variable predictora Xjen estudio. Para concluir si la hipótesis de riesgos proporcionales es cierta para la variable predictora en cuestión, será necesario estudiar el contraste cuya hipótesis nula es H0:γ  0. En el caso de que el coeficiente sea significativo (γ  0), se concluirá que la hipótesis de riesgos proporcionales para la variable en estudio no se cumple. Nótese que para el estudio de la significación del coeficiente γse puede usar el estadístico de Wald o uno de razón de verosimilitud. Ambos, bajo la hipótesis nula siguen una distribución χ2con un grado de libertad. No obstante, también es posible estudiar la validez de la hipótesis de riesgos proporcionales de manera conjunta para varios predictores. Para ello, es necesario estudiar de manera simultánea si todos los coeficientes γjdel modelo que se presenta a continuación, son todos nulos: h p t | X q  h0 p t q exp  p ¸ j  1 r βjXj  γj p Xj  gj p t qqs  , donde gj p t q es la función del tiempo asociada a cada uno de ppredictores (j  1, . . . , p). Debe tenerse en cuenta que diferentes variables predictoras pueden llevar asociadas distintas funciones del tiempo y que la hipótesis nula del contraste a estudiar es H0:γ1    γp  0. Este contraste requiere el uso de un estadístico de razón de verosimilitud que bajo hipótesis nula sigue una distribución χ2con pgrados de libertad. Además, en caso de que no haya evidencias significativas en contra de la hipótesis nula, se concluirá que la hipótesis de riesgos proporcionales conjunta del modelo de Cox no se cumple. Para más detalles de esta aproximación al problema de la diagnosis del modelo, véase Kleinbaum y Klein (2012). 2.6. Evaluación del modelo Una vez comprobadas y verificadas las hipótesis de riesgos proporcionales de los predictores involucrados en el modelo, tanto individualmente como conjuntamente, es recomendable evaluar la calidad del modelo creado. Para ello, al igual que sucede en cualquier modelo de regresión, se puede utilizar una medida análoga al coeficiente de determinación R2que mida el rendimiento del modelo como la proporción de varianza explicada por dicho modelo. Según Schemper y Stare (1996), en un modelo de riesgos proporcionales no existe una medida única, sencilla, útil y fácil de interpretar y calcular como coeficiente de determinación, dado que tal y como señalan Schemper (1990), Korn y Simon (1990), el coeficiente de determinación R2usual es muy sensible a la distribución de los tiempos censurados. Por ello, se presentan a continuación diferentes posibles medidas que actúan como el R2de un modelo de regresión lineal. CAPÍTULO 2. MODELO DE COX 14 2.7. VARIANTES DEL MODELO En primer lugar, se adjunta el coeficiente sugerido por Nagelkerke (1991): R2 p  1  exp  2 n p l  0  l  p q  ,(2.4) donde l  pes la log-verosimilitud parcial para el modelo ajustado con todas las covariantes (p covariantes), l  0es la log-verosimilitud parcial para el modelo sin covariantes y nel número de sujetos o individuos existentes en el estudio. Dicho coeficiente tiene las siguientes propiedades: Es consistente con la definición del R2clásico de una regresión lineal. Es consistente con el método de estimación de máxima verosimilitud. Es asintóticamente independiente del tamaño de la muestra, n. Puede ser interpretado como la proporción de variación explicada por el modelo. Es adimensional, es decir, no depende de las unidades utilizadas. Otras medidas posibles, son las propuestas por O’Quigley et al. (2005)yRoyston (2006). La sugerida por O’Quigley et al. (2005) reemplaza en (2.4) el número de individuos o sujetos npor el número de eventos observados m, de modo que, esta nueva medida depende menos del porcentaje de observaciones censuradas y se conoce como medida de aleatoriedad explicada: R2 p,e  1  exp  2 m p l  0  l  p q  . Por último, la medida propuesta por Royston (2006) dependiente de la anterior, tiene como objetivo parecerse más al coeficiente de determinación R2para la regresión lineal y sugiere la siguiente formulación: R2 p,v  R2 p,e R2 p,e  π2 6 p 1  R2 p,e q . 2.7. Variantes del modelo Además del modelo de Cox estudiado, existen variaciones o modificaciones del mismo para la adaptación de dicho modelo a diferentes situaciones. Entre las modificaciones, destacan el modelo de Cox estratificado, los modelos de fragilidad y el modelo de Cox con covariantes dependientes del tiempo. A continuación, se incluye una breve introducción a cada una de las variantes citadas y posteriormente, en el capítulo siguiente, se profundizará en el modelo de Cox con covariantes dependientes del tiempo, tema principal del presente trabajo: CAPÍTULO 2. MODELO DE COX 15 2.7. VARIANTES DEL MODELO a) Modelo de Cox estratificado: tal y como se ha presentado, se trata de una modificación del modelo de Cox en el que se controla con el uso de una estratificación de un predictor o predictores, el hecho de que dicho predictor o dichos predictores no verifiquen la hipótesis de riesgos proporcionales. Este tipo de modelo puede ser aplicado para la estratificación de varias variables explicativas en varios estratos o para la estratificación de un único predictor. Partiendo del hecho de que en un modelo de Cox existan rvariables (Z1, . . . , Zr) que no satisfacen la hipótesis de riesgos proporcionales y svariables (X1, . . . , Xs) que sí las satisfacen de las pvariables en estudio, se considerará una nueva variable Z  a partir de las Z’s que se utilizarán para la estratificación. Para cada Zise crearán categorías y las combinaciones de dichas categorías conformarán los diferentes estratos correspondientes a las categorías de la variable creada Z  . De esta forma, para cada estrato, existirá un modelo general de la forma: hg p t | X q  h0g p t q exp r β 1 X s , donde g  1, . . . , k  , con k  el número de estratos de la variable Z  ,β  p β1, . . . , βs q el vector de parámetros del modelo, X  p X1, . . . , Xs q el vector de variables que cumplen la hipótesis de riesgos proporcionales y h0gla función de riesgo base para cada uno de los estratos considerados. En este caso, se obtienen tantas curvas de supervivencia como estratos tiene la variable Z  y la forma de estimar los parámetros β’s del modelo será a través de la maximización de una función de verosimilitud parcial resultado del producto de todas las k  funciones de verosimilitud parcial de cada uno de los estratos. Adicionalmente, en este tipo de modelo, se debe estudiar si se cumple o no lo que se conoce como la hipótesis de la no interacción, supuesto que se refiere al hecho de que los sparámetros βdel modelo (β  p β1, . . . , βs q ) no varían en los diferentes estratos. b) Modelos de fragilidad: el modelo de Cox asume que la población en estudio es homogénea, o en otras palabras, que todos los individuos poseen el mismo riesgo de sufrir el evento de interés o independencia de los tiempos de supervivencia de los diferentes sujetos. No obstante, dicha hipótesis no es válida en muchos estudios de forma que es necesario introducir la heterogeneidad no observable (efecto aleatorio) con el uso de una componente de fragilidad. De esta forma, esta componente es capaz de dividir la variabilidad total del estudio en una parte que depende de las variables predictoras, y por tanto, que se puede predecir y otra teóricamente impredecible que es explicada por la componente de fragilidad. En concreto, la fragilidad es una medida de riesgo relativo, ya que cuanto mayor sea la fragilidad de un individuo con respecto al evento o suceso de interés, mayor probabilidad habrá de que se produzca ese evento en dicho sujeto. Este tipo de método es un modelo multiplicativo de riesgo compuesto por tres elementos: la componente de fragilidad, la función de riesgo base y la componente que tiene en cuenta la influencia de las covariables usadas. La componente de fragilidad, representada por Z, asume una distribución de probabilidad (normalmente es considerada la distribución Gamma, dado que ayuda a la hora de maximizar la función de verosimilitud parcial necesaria para obtener los estimadores de los parámetros del modelo) pues se trata de una variable aleatoria y la formulación general del modelo queda de la siguiente manera: hi p t | Xi, Z q  Zh0 p t q exp p β 1 Xi q . CAPÍTULO 2. MODELO DE COX 16 2.8. TÉCNICAS DE SELECCIÓN DE VARIABLES c) Modelo de Cox con covariantes dependientes del tiempo: ampliación del modelo de Cox en el que algunas o todas las variables explicativas o predictoras dependen del tiempo. En este caso, la formulación del modelo es: hi p t | Xi p t qq  h0 p t q exp  p1 ¸ j  1 βjXij  p2 ¸ k  1 γkXik p t q  , que difiere con la formulación del modelo de Cox sin covariantes dependientes del tiempo en el hecho de que en la parte exponencial, se incluyen tanto los predictores independientes del tiempo, Xi  p Xi1, . . . , Xip1 q , como los predictores dependientes del tiempo, Xi p t qp Xi1 p t q , . . . , Xip2 p t qq . Al igual que en el modelo de Cox, los parámetros del modelo se estiman mediante un procedimiento de máxima verosimilitud aunque en este tipo de modelo, los cálculos son más complicados dado que los conjuntos de riesgo no son triviales. Además, en este tipo de modelos, no suele cumplirse la hipótesis de riesgos proporcionales y el riesgo en el instante tde que un individuo sufra el evento de interés, depende del valor de la variable Xi p t q en ese mismo momento. En el siguiente capítulo, se profundizará en este modelo, dado que es el tema principal del presente Trabajo Fin de Máster. 2.8. Técnicas de selección de variables Tal y como indican Guyon y Elisseeff (2003), la selección de variables en un modelo de regresión posee diferentes ventajas en el modelo como la mejora del rendimiento de la predicción de las variables explicativas, la obtención de predictores más rápidos, la reducción de la dimensionalidad de un problema o el evitar la multicolinealidad. En general, el aumento del número de variables predictoras consideradas en un modelo de regresión, lleva consigo un mayor número de parámetros a estimar y una disminución de la precisión individual de cada una de las estimaciones que se traduce en una mayor varianza y en un sobreajuste de la función de regresión estimada. En contraposición, si se consideran menos variables de las que debiera en el modelo, los sesgos aumentan mientras que las varianzas disminuyen, obteniéndose así una mala descripción de los datos. Otro aspecto a considerar, es la correlación existente entre las variables explicativas: si se consideran predictores muy correlacionados entre ellos, la confianza de la predicción del modelo se ve perjudicada. De este modo, los métodos de selección de variables ayudan a encontrar un modelo apropiado que se ajuste bien a los datos y que a su vez, posea un equilibrio entre su bondad de ajuste y simplicidad. En el ámbito del análisis de supervivencia, se pueden aplicar las mismas técnicas de selección de variables que se usan en cualquier otro modelo de regresión: selección paso a paso, selección de los mejores subconjuntos de covariantes, método backward. .. En el capítulo 5 del libro de Lemeshow et al. (2008) se detallan tres métodos útiles para la selección de variables entre los que se encuentran la selección paso a paso o la selección de los mejores subconjuntos de predictores. CAPÍTULO 2. MODELO DE COX 17 3.5. ESTIMACIÓN DE PARÁMETROS Existen numerosos algoritmos para la resolución de este tipo de problemas de optimización con restricciones como el método de Newton. Al igual que se expuso en el Apartado 2.3 del Capítulo 2, dichos métodos son iterativos, ya que a partir de una estimación inicial de cada parámetro, el algoritmo va iterando hasta su convergencia. Sin embargo, cuando mes grande, existe un algoritmo de fácil implementación y muy eficiente para la resolución de las ecuaciones (3.5), (3.6) y (3.7). Dicho algoritmo, se presenta a continuación y comienza estimando β p k q ,γ p k q yθ p k q . Una vez obtenidas dichas estimaciones, la iteración k  1se calcula siguiendo los siguientes pasos: (1) Calcular β p k  1 q tal que l p β p k  1 q , γ p k q , θ p k q q ¥ l p β p k q , γ p k q , θ p k q q . (2) Calcular γ p k  1 q tal que l p β p k  1 q , γ p k  1 q , θ p k q q ¥ l p β p k  1 q , γ p k q , θ p k q q . (3) Calcular θ p k  1 q ¥ 0tal que l p β p k  1 q , γ p k  1 q , θ p k  1 q q ¥ l p β p k  1 q , γ p k  1 q , θ p k q q . Estas condiciones aseguran que l p β p k  1 q , γ p k  1 q , θ p k  1 q q ¥ l p β p k q , γ p k q , θ p k q q al final de la iteración k  1, lo cual es un requisito clave para la convergencia del algoritmo en estudio. Los pasos 1 y 2 se resuelven usando el algoritmo de Newton incorporando pasos de búsqueda lineal12, y el paso 3, se resuelve mediante un algoritmo multiplicativo-iterativo (MI) diseñado para respetar las restricciones no negativas de θ. Para más información sobre este último tipo de algoritmo, consultar Chan y Ma (2012) y Ma (2010). El algoritmo presentado se conoce con el nombre de algoritmo Newton-MI. Para actualizar el valor de β, se emplea una iteración del algoritmo de Newton con búsqueda lineal. Comenzando con β p k q y usando una búsqueda lineal con tamaño de paso ω p k q 1 P p 0,1 s , se tiene: β p k  1 q  β p k q  ω p k q 1 p X 1 A p k q X q  1X 1 p A p k q 1n  δ q , donde Xes la matriz n  p1compuesta por todas las covariantes no dependientes del tiempo, Aes una matriz diagonal definida como A  diag p H  01 p y1 q exp p x 1 1β q , . . . , H  0n p yn q exp p x 1 nβ qq , 1nes un vector de 1’s de tamaño nyδes un vector de ncomponentes para los δi’s. Además, se tiene que A p k q es Acon β  β p k q ,γ  γ p k q yθ  θ p k q . La matriz X 1 AX es el hessiano negativo de l p β, γ, θ q con respecto a β, y el parámetro de búsqueda lineal ω p k q 1 ayuda a conseguir que l p β p k  1 q , γ p k q , θ p k q q ¥ l p β p k q , γ p k q , θ p k q q . Análogamente, el método de Newton con búsqueda lineal se aplica para la actualización de γ. Primero, se considera la matriz Z p t qr Z 1 1 p t q , . . . , Z 1 n p t qs 1 donde Zi p t q  r zi1, . . . , zi,ni  1 s 1 yzia está definido con la Ecuación 3.3. La matriz Z p t q es la matriz del modelo asociada a las covariantes dependientes del tiempo y su dimensión es N  p2siendo N  ° ini. Sea también B  diag p exp p x 1 1β q B1,...,exp p x 1 nβ q Bn q , 12Enfoque iterativo y básico en los problemas de optimización que se basa en la búsqueda de una dirección de descenso a lo largo de la cual la función objetivo del problema se reducirá. CAPÍTULO 3. MODELO DE COX CON COVARIANTES DEPENDIENTES DEL TIEMPO 24 3.5. ESTIMACIÓN DE PARÁMETROS donde Bies una matriz diagonal de tamaño ni  nicuyos elementos de la diagonal son r H0 p ri,a  1 q  H0 p ri,a qs exp p z 1 iaγ q para a  0, . . . , ni  1. Sea ζun N-vector dado por ζ  r ζ11, . . . , ζ1,n1  1, . . . , ζn1, . . . , ζn,nn  1 s 1 donde ζia  1solo si a  ni  1yδi  1, en cualquier otro caso, ζia  0. De esta forma, el algoritmo de Newton con búsqueda lineal actualiza γde acuerdo a la ecuación siguiente: γ p k  1 q  γ p k q  ω p k q 2 p Z 1 p t q B p k q Z p t qq  1Z 1 p B p k q 1N  ζ q , donde ω p k q 2 P p 0,1 s se corresponde con el tamaño del paso de búsqueda lineal y B p k q denota Bpero con β  β p k  1 q ,γ  γ p k q yθ  θ p k q . La matriz Z 1 p t q BZ p t q es el hessiano negativo de l p β, γ, θ q con respecto a γ. Nótese que ω p k q 2se usa para conseguir que se verifique que l p β p k  1 q , γ p k  1 q , θ p k q q ¥ l p β p k  1 q , γ p k q , θ p k q q . Finalmente, para actualizar el parámetro θ, se usa un algoritmo de tipo multiplicativoiterativo (MI) que respeta la restricción no negativa de θy, además, es fácil de implementar. Sean CyC  matrices n  mcuyos elementos p i, u q son ψu p yi q yΨ  ui p yi q , respectivamente. Sea δel vector de dimensión npara δiyfel n-vector para exp p xiβ q . El algoritmo MI actualiza θde acuerdo con la siguiente ecuación: θ p k  1 q  θ p k q  ω p k q 3S p k q p C 1 r D p k q s  1δ  r C p k q s 1 f p k q q , donde DySson matrices diagonales cuyos elementos son h0 p yi q yθu {p ° iΨ  ui p yi q exp p x 1 iβ q ε q , respectivamente. Además, εes un pequeño umbral utilizado para evitar que el denominador correspondiente sea nulo y ω p k q 3 P p 0,1 s para garantizar que se verifique que l p β p k  1 q , γ p k  1 q , θ p k  1 q q ¥ l p β p k  1 q , γ p k  1 q , θ p k q q . Todas las búsquedas lineales utilizadas en el algoritmo, pueden calcularse eficientemente con el uso de la regla de Armijo13 (Luenberger y Ye (1984)). Asimismo, cuando A1 { 2Xy B1 { 2Z p t q tienen rangos máximos, las matrices X 1 AX yZ 1 p t q BZ p t q son definidas positivas de modo que las actualizaciones de βyγestán bien definidas. Por último, destacar que se puede demostrar que si θ p k q es no negativa, entonces θ p k  1 q es también no negativa y que bajo ciertas condiciones de regularidad, el algoritmo expuesto converge a la solución satisfaciendo las condiciones de KKT. Para profundizar en estos últimos resultados, consultar Chan y Ma (2012). 3.5.1. Propiedades asintóticas Los resultados asintóticos son fundamentales para realizar inferencias sin la necesidad de usar métodos computacionales complejos. Sean p β0, γ0, h00 p t qq los verdaderos parámetros y p βy p γlas estimaciones por máxima verosimilitud de βyγrespectivamente. Además, se toma la notación siguiente: hn p t q  m ¸ u  1 θuψu p t q y p hn p t q  m ¸ u  1 p θuψu p t q , 13Método esencial en la optimización numérica y búsqueda de línea en problemas de optimización no lineal, especialmente en algoritmos iterativos como el descenso de gradiente, ya que determina el tamaño adecuado del paso a lo largo de una dirección de búsqueda, evitando pasos excesivamente grandes y garantizando la convergencia de los algoritmos de optimización hacia un mínimo local en una función objetivo. CAPÍTULO 3. MODELO DE COX CON COVARIANTES DEPENDIENTES DEL TIEMPO 25 3.6. CONTRASTES DE HIPÓTESIS donde p θues la estimación de máxima verosimilitud del parámetro θuymaumenta con n. Sea también θel vector de todos los θuy p θel vector de los p θu. Según Xu et al. (2018), se pueden obtener resultados con consistencia fuerte para p β, p γy p hn p t q cuando m Ñ 8 yn Ñ 8 pero m { n Ñ 0. Además, existen resultados que aseguran la normalidad asintótica de las estimaciones de p β, p γy p θpara un mfijo. Asimismo, se puede demostrar que la varianza asintótica es exacta cuando se compara con la varianza obtenida a partir del método de Monte Carlo. Estos últimos resultados se pueden consultar y profundizar en Thackham y Ma (2020). 3.6. Contrastes de hipótesis Al igual que ocurre en el modelo de Cox con covariantes independientes del tiempo, una vez estimados los parámetros βyγdel modelo en estudio, es necesario realizar contrastes para estudiar si dichos parámetros pueden considerarse significativos o no. Dadas las propiedades asintóticas anteriores, se puede considerar que tanto p βcomo p γsiguen una distribución normal de media βyγ, y matriz de varianzas y covarianzas Σβ  Φ  1 p β q de dimensión p1  p1yΣγ  Φ  1 p γ q de dimensión p2  p2, respectivamente. Por lo tanto, se concluye que p βa  Np1 p β, Φ  1 p β qq y p γa  Np2 p γ, Φ  1 p γ qq . Asimismo, la media y la matriz de varianzas y covarianzas de las distribuciones normales anteriores, se pueden estimar de la siguiente forma: p i q p β  βyΣβ  Φ  1 p p β q con Φ p β qB 2l  p β q B β B β 1 , p ii q p γ  γy p Σγ  Φ  1 p p γ q con Φ p γ q  B 2l  p γ q B γ B γ 1 . En este modelo, habrá que considerar contrastes de hipótesis individuales para cada uno de los parámetros βjcon j  1, . . . , p1yγkcon k  1, . . . , p2. Estos contrastes se expresan en los siguientes términos: p i q " H0:βj  0 H1:βj  0 p ii q " H0:γk  0 H1:γk  0 para cada j  1, . . . , p1yk  1, . . . , p2, respectivamente. Dichos contrastes se resuelven utilizando el estadístico de Wald tal y como se indica en el Apartado 2.4 del Capítulo 2. Además, también es posible considerar el contraste de hipótesis conjunto para cada uno de los parámetros βyγ. Estos contrastes por separado, se pueden resolver usando los tres métodos y las mismas consideraciones citadas en el Apartado 2.4 anteriormente indicado. A continuación, se presentan los dos contrastes de interés: p i q " H0:β  p β10, . . . , βp10 q 1 H1:βj  βj0para algún j  1, . . . , p1, p ii q " H0:γ  p γ10, . . . , γp20 q 1 H1:γ  γj0para algún j  1, . . . , p2. CAPÍTULO 3. MODELO DE COX CON COVARIANTES DEPENDIENTES DEL TIEMPO 26 3.7. EVALUACIÓN DEL MODELO Asimismo, es posible considerar el siguiente contraste de hipótesis conjunto de βyγ: " H0: p β, γ q  p β10, . . . , βp10, γ10, . . . , γp20 q 1 H1:βj  βj0oγk  γk0para algún j  1, . . . , p1ok  1, . . . , p2. Este último contraste se suele utilizar en los diferentes software de programación para mostrar la significación del modelo en estudio. Puede resolverse usando por ejemplo, el test de razón de verosimilitud de forma que el estadístico asociado a dicho contraste sería: XLR  2 p log L  p ˆ β, ˆγ q  log L  p β0, γ0 qq . Bajo hipótesis nula, el estadístico enunciado sigue una distribución χ2con p1  p2grados de libertad (XLR  χ2 p1  p2). 3.7. Evaluación del modelo Tal y como ocurre en cualquier modelo de regresión, una vez contrastadas la significación o no significación de los parámetros involucrados en el modelo, es conveniente evaluar la capacidad predictiva del modelo creado. Análogamente a lo que sucedía con el modelo de Cox con covariantes independientes del tiempo, una posible medida para evaluar la calidad del modelo es usando un coeficiente de determinación R2adaptado al tipo de datos con el que se trabaja en el análisis de supervivencia. En el modelo de Cox con covariantes dependientes del tiempo, se pueden considerar las mismas medidas de evaluación sugeridas en el Apartado 2.6 adaptadas al modelo objeto de estudio de este capítulo. Entre ellos, se encuentran: El coeficiente sugerido por Nagelkerke (1991) definido como: R2 p  1  exp  2 n p l  0  l  p q  , donde l  pes la log-verosimilitud parcial para el modelo ajustado con todas las covariantes (p  p1  p2covariantes), l  0es la log-verosimilitud parcial para el modelo sin covariantes y nel número de sujetos o individuos existentes en el estudio. El coeficiente propuesto por O’Quigley et al. (2005) y conocido como medida de aleatoriedad explicada: R2 p,e  1  exp  2 m p l  0  l  p q  , siendo mel número de eventos observados, l  pla log-verosimilitud parcial para el modelo ajustado con todas las covariantes y l  0la log-verosimilitud parcial para el modelo sin covariantes y nel número de sujetos o individuos existentes en el estudio. El coeficiente sugerido por Royston (2006) depende del anterior y se define como: R2 p,v  R2 p,e R2 p,e  π2 6 p 1  R2 p,e q . CAPÍTULO 3. MODELO DE COX CON COVARIANTES DEPENDIENTES DEL TIEMPO 27 3.8. TÉCNICAS DE SELECCIÓN DE VARIABLES 3.8. Técnicas de selección de variables Por último, en ocasiones, puede ser interesante considerar el modelo de Cox con un menor número de variables dado que el aumento del número de variables predictoras consideradas en un modelo de regresión, lleva consigo un mayor número de parámetros a estimar y una disminución de la precisión individual de cada una de las estimaciones, lo que se traduce en una mayor varianza y en un sobreajuste de la función de regresión estimada. De esta forma, el considerar un menor número de variables en el modelo podría evitar la multicolinealidad entre las variables explicativas o la obtención de un modelo con una mayor capacidad predictiva. La solución a estas casuísticas se encuentra con el uso de técnicas dedicadas a la selección de variables como ocurre en cualquier tipo de regresión. Entre estas técnicas, tal y como se comentó en el Apartado 2.8, se encuentran por ejemplo, la selección paso a paso o el método forward. Tal y como se especificó en el Capítulo 2, en el libro de Lemeshow et al. (2008) se detallan tres métodos útiles para la selección de variables entre los que se encuentran la selección paso a paso o la selección de los mejores subconjuntos de predictores. Asimismo, se puede hacer uso del enfoque del conjunto activo primal-dual para resolver el problema de la selección de subconjuntos de variables y propuesto por Ito y Kunisch (2013). Para profundizar en esta técnica, consultar Jareño (2022). En general, los métodos de selección de variables ayudan a encontrar un modelo apropiado que se ajuste bien a los datos y que a su vez, posea un equilibrio entre su bondad de ajuste y simplicidad. CAPÍTULO 3. MODELO DE COX CON COVARIANTES DEPENDIENTES DEL TIEMPO 28 Capítulo 4 Modelo de Cox en R En este capítulo, se detallará cómo se puede usar el modelo de Cox tanto con covariantes independientes como dependientes del tiempo en el software libre R. Para ello, se irán exponiendo las librerías o paquetes necesarios para realizar dicho estudio y las funciones imprescindibles para obtener las conclusiones oportunas. Para mostrar cada uno de los modelos de Cox considerados, se detallarán los argumentos y la salida de cada una de las funciones más relevantes. Asimismo, se adjuntarán ejemplos con datos disponibles en los propios paquetes o con datos extraídos de la asignatura Modelado y Predicción Estadística. Cabe destacar que se usará un nivel de significación del 5 % en todos los contrastes a realizar. 4.1. Introducción La librería más habitual a la hora de tratar datos de análisis de supervivencia es la librería survival (Therneau (2023)). Dicha librería posee funciones esenciales para el análisis de supervivencia y trabaja con objetos de tipo Surv. Este tipo de objetos se caracteriza por combinar información temporal y de censura para cada uno de los individuos u observaciones en estudio; y se puede crear con la función Surv(time, time2, event, type, origin) cuyos argumentos se enuncian a continuación: time: en el caso de que los datos estén censurados a la derecha, en este argumento se anotará el tiempo de seguimiento. En contraposición, cuando la censura sea por intervalos, se introducirá el instante inicial del intervalo. event: indicador del estado de la observación o individuo en estudio, 0  vivo y1  muerto. No obstante, cuando se trabaja con datos censurados por intervalos, el estado puede ser 0  censura a la derecha,1  evento observado,2  censura a la izquierda o3  censura por intervalos. time2: argumento necesario cuando se trabaja con datos censurados por intervalos (solo admite intervalos semiabiertos de la forma p inicio, fin s ) y debe anotar el instante final del intervalo. Además, cuando se trabaja con datos de procesos de conteo el argumento event indica si se produjo el evento al final del intervalo. 29 4.2. MODELO DE COX type: indicador del tipo de censura y admite los siguientes valores “right”,“left”, “counting”,“interval”,“interval2” o“mstate”. Para una mayor explicación de los posibles valores del argumento en estudio, consúltese Therneau (2023). origin: útil cuando se trabaja con datos de procesos de conteo e indica la función de riesgo al comienzo del estudio. No obstante, dicha librería no va a ser la única que se va a utilizar a lo largo del presente capítulo. Otras librerías útiles para el estudio del modelo de Cox en R son la librería survMisc (Dardis (2022)), la librería survminer (Biecek et al. (2021)) o la librería glmnet (Friedman et al. (2023)). Por ejemplo, la primera de ellas permite consultar, entre otros conceptos del modelo de Cox, los diferentes coeficientes R2expuestos a lo largo del trabajo; y, la segunda, ayuda a presentar gráficos detallados. 4.2. Modelo de Cox Para ajustar un modelo de Cox en R a partir de ciertas variables explicativas o covariantes, se usa la función coxph() de la librería survival (Therneau (2023)). Esta función coxph() permite crear el modelo de Cox con las covariantes que se le especifique en dicha función a partir de una fórmula. Los argumentos de dicha función se detallan a continuación: library(survival) args(coxph) ## function (formula, data, weights, subset, na.action, init, control, ## ties = c("efron", "breslow", "exact"), singular.ok = TRUE, ## robust, model = FALSE, x = FALSE, y = TRUE, tt, method = ties, ## id, cluster, istate, statedata, nocenter = c(-1, 0, 1), ...) ## NULL donde: formula: objeto de tipo fórmula en el que la parte de la derecha contiene un objeto de tipo Surv() seguido por ~ y las variables explicativas o predictoras que se quieran utilizar. data: conjunto de datos que contiene las variables objeto de estudio para el modelo de Cox. weights: vector de ponderaciones que se utilizará en el caso que corresponda. subset: subconjunto de datos que debe ser utilizado en el ajuste a realizar. Es muy útil cuando se trabaja con conjunto de datos de entrenamiento y test. na.action: función que filtra los datos de tipo NA. CAPÍTULO 4. MODELO DE COX EN R 30 4.2. MODELO DE COX init: vector de valores iniciales de la iteración del método que calcula los coeficientes del modelo. control: objeto de clase coxph.control que especifica el límite de iteraciones y otro tipo de opciones de control del ajuste del modelo de Cox. ties: cadena de caracteres que especifica el método a utilizar en caso de empates en los tiempos de fallo. En el caso de que no existan tiempos de fallo coincidentes, cualquier método es equivalente. Por defecto, se utiliza el método de Efron “efron” (más preciso y computacionalmente más eficiente) pero otras opciones posibles son el método de Breslow “breslow” o el método Exact Partial Likelihood “exact”, aconsejable cuando se trabaja con un conjunto pequeño de valores discretos. singular.ok: valor lógico que indica qué hacer en caso de colinealidad en la matriz del modelo. Por defecto, toma el valor TRUE de manera que el programa no considera las columnas de la matriz Xque son combinaciones lineales de las columnas restantes. En dicho caso, los coeficientes para dichas columnas serán NA y la matriz de varianzas contendrá ceros. Además, para los cálculos auxiliares como el predictor lineal, los valores perdidos se tratan como valores nulos. robust: valor lógico que garantiza el cálculo de una varianza robusta en el caso de tomar el valor TRUE (valor por defecto). id: variable opcional que se utiliza para identificar a los sujetos o individuos de la muestra cuando existen varias filas que hacen referencia al mismo individuo. cluster: variable opcional que agrupa las observaciones en clústeres y se utiliza cuando el argumento robust toma el valor TRUE. Normalmente, esta variable se encuentra en los datos de estudio. Para un estudio más detallado sobre el resto de argumentos que admite la función en estudio, consultar la función coxph() y el apartado coxph.control del manual de Therneau (2023). En la salida que proporciona la función se podrá consultar una tabla con la siguiente información para cada una de las variables explicativas involucradas en el modelo: coef : coeficientes estimados asociados a cada uno de las variables predictoras. exp(coef): exponencial de los coeficientes estimados y se interpretan como efectos multiplicativos sobre la función de riesgo h p t q . se(coef): errores de estimación asociados a los coeficientes estimados. zyp: muestran la significación de cada una de las variables consideradas en el modelo creado. El parámetro zmuestra el valor del estadístico de Wald asociado a cada variable junto con el pque se corresponde con el p-valor asociado. En función del valor de pse podrán considerar nulos o no los coeficientes de las variables del modelo en estudio. CAPÍTULO 4. MODELO DE COX EN R 31 4.2. MODELO DE COX Asimismo, en la salida, se puede consultar también el estadístico de razón de verosimilitud del contraste de hipótesis conjunto del modelo junto con su grado de libertad correspondiente y p-valor. De esta forma, de manera sencilla se puede concluir si el modelo de Cox considerado es óptimo o se debe ajustar con menos variables explicativas u otras variables diferentes. La salida también adjunta el tamaño de la muesta en estudio, n, el número de eventos que se han producido a lo largo del tiempo de estudio y el número de observaciones que se han eliminado por ser valores perdidos. No obstante, un estudio más completo de la salida se puede obtener ejecutando la orden summary() del modelo creado. En este caso, se muestra una tabla en la que además de la información expuesta anteriormente, se recoge para cada una de las variables explicativas consideradas en el modelo, el valor de exp(-coef), el extremo inferior (lower .95) y superior (upper .95) del intervalo de confianza asociado a cada una de las exponenciales de los coeficientes estimados cambiados de signo e información relevante de cada uno de los test (valor del estadístico, grados de libertad asociado y p-valor) de Wald, de razón de verosimilitud y test score correspondientes al contraste de hipótesis conjunto del modelo. Por defecto, el intervalo de confianza se calcula al 95%. En el caso de que se quiera conocer el valor del intervalo de confianza asociado directamente a los coeficientes estimados, se debe usar la función confint(). La estructura de dicha función es: confint(x, parm, level), donde xes el objeto coxph que contiene el modelo de Cox creado, parm es el vector que contiene el nombre de las variables para las cuales se quiere estudiar el intervalo de confianza de su coeficiente estimado y level indica el nivel de confianza para el cual se quiere consultar el intervalo de confianza (por defecto, toma el valor 0.95). La salida de la función es una matriz (o vector) con dos columnas: en la primera de ellas se muestra el límite inferior del intervalo de confianza y en la segunda, el límite superior de este. Ambas columnas se etiquetan como (1 - nivel)/2 y 1 - (1 - nivel)/2 en porcentaje. En ocasiones, puede ser de interés observar el comportamiento de la estimación de los coeficientes del modelo de Cox creado, para ello, existe la orden profLik() de la librería survMisc. Dicha orden presenta como salida una gráfica por cada una de las variables explicativas contempladas en el modelo de Cox bajo estudio en la que se puede consultar la evolución del valor de la función de verosimilitud parcial frente al coeficiente estimado para cada una de las iteraciones realizadas en el ajuste del modelo. Además, en cada gráfica se puede consultar el intervalo de confianza asociado al coeficiente estimado (máximo de la función dibujada) a través de dos circunferencias que denotan el extremo inferior y superior de este. El modo en el que se usa la función en estudio es: profLik(x, CI, interval, devNew, ...), siendo xel objeto coxph que contiene el modelo de Cox bajo estudio, CI el nivel a usar para mostrar el intervalo de confianza (por defecto, toma el valor 0.95), interval el número de puntos sobre los que evaluar el coeficiente (por defecto, toma el valor 50), devNew el valor lógico para indicar si las gráficas se quieren mostrar en otra pestaña o no y parámetros gráficos adicionales como col para indicar el color de la gráfica a dibujar. Por último, destacar también que la clase de objeto coxph contiene más información que se puede consultar para la evaluación del modelo como la matriz de covarianzas de los estimadores, el valor de la log-verosimilitud del modelo con los valores iniciales y finales de los coeficientes, las medias de las variables predictoras, el número de iteraciones CAPÍTULO 4. MODELO DE COX EN R 32 4.2. MODELO DE COX realizadas para la estimación de los coeficientes14 o el vector de predictores lineales de las observaciones del conjunto de estudio. Para conocer con más detalle la información que contiene el objeto de tipo coxph, consultar el apartado coxph.object del manual de Therneau (2023). Además, es posible realizar un estudio de los residuos de Schoenfeld asociados a cada observación del conjunto de datos bajo estudio con la orden residuals(object, type = c("schoenfeld")), siendo object el modelo de Cox creado. En el caso de que se quiera realizar un estudio gráfico del modelo de Cox creado, existen funciones como ggforest() del paquete survminer que pueden ayudar a realizarlo. Esta función tiene la peculiaridad de mostrar gráficamente la exponencial de los coeficientes obtenidos en el modelo ajustado junto con el intervalo de confianza obtenido para este, así como presentar un pequeño resumen que contiene el p-valor asociado a cada uno de los contrastes individuales de los coeficientes del modelo y al contraste conjunto, el AIC del modelo o el número de observaciones asociadas a cada una de las variables predictoras del modelo. La estructura de la función es: ggforest(model, data, main, cpositions, fontsize, refLabel, noDigits) y cuyos argumentos se definen a continuación: model: objeto de tipo coxph. data: conjunto de datos usado para ajustar el modelo en estudio. main: título del gráfico a realizar. cpositions: posición relativa en el eje OX de las tres primeras columnas que se presentan en la gráfica. fontsize: tamaño de la fuente utilizada para realizar las anotaciones que aparecen en la gráfica. refLabel: etiqueta de referencia para las variables de tipo factor. noDigits: número de dígitos para las estimaciones y los p-valores presentados en el estudio realizado. Una vez creado un modelo de Cox, también se puede hacer un estudio de las curvas de supervivencia ajustadas para el modelo de regresión creado. Para dicho objetivo, existen dos órdenes en el paquete survminer:ggadjustedcurves() ysurv_adjustedcurves(). La primera de ellas, realiza una gráfica de la estimación de las curvas de supervivencia del modelo de Cox creado con la función coxph(), mientras que la segunda de ellas, la estima y la muestra de manera tabular en cada instante de tiempo t. Ambas órdenes comparten argumentos y estos, se detallan en la siguiente lista: fit: objeto de tipo coxph que contiene el modelo de Cox creado. 14El ajuste del modelo de Cox se puede controlar con diferentes parámetros numéricos de la orden coxph.control() como el número máximo de iteraciones para la convergencia del modelo. CAPÍTULO 4. MODELO DE COX EN R 33 4.3. APLICACIONES stype: tipo de cálculo de la curva de supervivencia. Puede tomar los valores 1(cálculo de la curva de supervivencia) o 2(cálculo de la exponencial del riesgo acumulado). Por defecto, toma el valor 2. ctype: indica si el cálculo de la tasa de riesgo acumulativo debe calcularse con una corrección. Los valores posibles son 1(no) o 2(sí). Por defecto, toma el valor 2. conf.type: indicador del tipo de estimación del intervalo a realizar. Los posibles valores son “log” (valor por defecto y genera intervalos basados en la tasa de riesgo acumulada o el logaritmo de la curva de supervivencia) “log-log” (genera el logaritmo de la opción “log”), “plain” (genera intervalos de tipo estándar), “none” (no genera intervalos de confianza) o “logit” (genera el logaritmo del cociente de la curva de supervivencia entre la unidad menos la curva de supervivencia). censor: valor lógico que indica si se incluyen o no en el resultado las observaciones sobre las que no se ha producido el evento en estudio. Por defecto, toma el valor TRUE. La estimación de la función de supervivencia se puede consultar con la orden summary() y la información que devuelve es la siguiente: time: instante de la observación o tiempos observados. n.risk: número de individuos en riesgo en cada instante de tiempo. n.event: número de individuos que presentan el evento en cada instante de tiempo. survival: estimación de la función de supervivencia. std.err: desviación estándar de la estimación realizada. lower 95 % CI yupper 95 % CI: extremo inferior y superior del intervalo de confianza de la estimación al 95%. Además, se puede consultar también el número de observaciones censuradas en cada instante de tiempo en estudio (n.censor). Para ampliar la información de manera más detallada, consultar los argumentos de esta orden en el apartado survfit.coxph del manual de Therneau (2023). De manera gráfica, se puede estudiar la estimación de la función de supervivencia calculada junto con el intervalo de confianza obtenido para dicha estimación. La mencionada gráfica se puede realizar con la ayuda del paquete survminer (Biecek et al. (2021)) y la función ggsurvplot(). 4.3. Aplicaciones En los próximos dos subapartados se pondrán en práctica las funciones y órdenes descritas en el apartado anterior. En el primero de ellos, las funciones se aplicarán a un conjunto de datos llamado uissurv.csv extraído de la asignatura de Modelado y Predicción Estadística y que está formado por variables independientes del tiempo. Dicho conjunto de datos CAPÍTULO 4. MODELO DE COX EN R 40 4.3. APLICACIONES permitirá mostrar el ajuste de un modelo de Cox con ese tipo de variables, las cuales se presentan y describen en el siguiente subapartado. Por otro lado, en el segundo de los subapartados se mostrará la forma con la que se debe trabajar cuando el conjunto de datos posee covariantes dependientes del tiempo. Para tal caso y con la ayuda de Therneau et al. (2022), se usarán los datos pbc yheart disponibles en el paquete survival. A lo largo del segundo subapartado de esta sección, se mostrarán ambos conjuntos de datos, explicando sus objetivos y detallando las variables de estudio de cada uno de ellos. Además de poder mostrar todo el potencial de las funciones que proporciona R, se mostrarán las conclusiones lógicas que derivan de las salidas obtenidas a lo largo del estudio de los modelos de Cox que aportarán información relevante para la realización de un análisis completo de dichos modelos. 4.3.1. Modelo de Cox con covariantes independientes del tiempo Tal y como se ha presentado en la introducción anterior, para la construcción y análisis de un modelo de Cox con covariantes independientes del tiempo se usará el conjunto de datos procedente del archivo uissurv.csv. Dicho fichero posee una serie de datos procedentes de un estudio titulado UMARU Impact Study (UIS), resultado de un proyecto colaborativo de investigación llevado a cabo entre los años 1989 y 1994 por la Unidad de Investigación sobre el SIDA de la Universidad de Massachusetts (UMARU) y dirigido por los doctores Jane McCusker, Carol Bigelow y Anne Stoddard (Lemeshow et al. (2008)). El objetivo del estudio fue comparar dos tratamientos diferentes diseñados para reducir el abuso de drogas y prevenir conductas de alto riesgo de contagio de VIH. Para ello, se diseñaron dos programas diferentes y el UIS buscaba determinar si la eficacia de los tratamientos en estudio dependía de la duración del programa diseñado. A continuación, se enuncian los dos programas contemplados en el estudio: a) Programa A: con una duración comprendida entre 3 y 6 meses, los 444 individuos sometidos a este programa debían asistir a un curso de educación sanitaria y prevención de recaídas. Durante dicho periodo de tiempo, los individuos aprenden a reconocer situaciones de alto riesgo que podían provocar una recaída y adquirían las habilidades necesarias para no usar las drogas en dichos momentos. b) Programa B: 184 individuos eran sometidos a un programa que duraba de 6 a 12 meses y en el que aprendían un estilo de vida muy estructurado en un entorno social de comunidad. Las variables involucradas en el estudio son: id: identificador de cada uno de los individuos de la muestra considerada en el estudio. age: edad de cada individuo expresada en años. beck: puntuación que el sujeto ha obtenido en el test de depresión de Beck (BDI)18. Los valores están comprendidos entre 0y54. 18El test de depresión de Beck es un cuestionario que consta de 21 preguntas de respuesta múltiple utilizado para medir la severidad de una depresión y fue desarrollado por el psiquiatra, investigador y fundador de la terapia cognitiva, Aaron T. Beck. CAPÍTULO 4. MODELO DE COX EN R 41 4.3. APLICACIONES hercoc: variable categórica que especifica si el sujeto ha usado heroína/cocaína durante los 3 meses previos a la participación en el estudio. Los posibles valores son: 1  Heroína & Cocaína (Heroin & Cocain), 2  Solo Heroína (Heroin only), 3  Solo Cocaína (Cocain only) o 4  Ni Heroína ni Cocaína (Neither). ivhx: indicador de consumo de drogas de uso intravenoso (drogas IV) antes de la admisión al programa. Los posibles valores son: 1  Nunca (Never), 2  Previamente (Previous) o 3  Recientemente (Recent). ndrugtx: número de tratamientos previos a la incorporación al programa. El rango de valores está comprendido entre 0y40. race: raza del individuo. Los posibles valores son: 0  Blanca (White) o 1  No blanca (Non-White). treat: variable que denota el tratamiento aleatorizado asignado a cada uno de los individuos. Los posibles valores son 0  Corto (Short q o1  Largo (Long). site: tipo de programa asociado a cada uno de los individuos. Los valores contemplados son 0  A o 1  B. los: tiempo de estancia en el tratamiento expresado en días y medido desde el día de admisión hasta el día de salida del programa. time: tiempo hasta la recaída en las drogas expresado en días y medido desde el día de admisión del individuo en el programa. censor: estado de la censura (1  Recaída en las drogas o 0  censurado). Una vez presentadas cada una de las variables, se leerá el fichero que contiene dichos datos: uissurv <- read.csv("uissurv.csv",header = TRUE, sep = ";") Además, se mostrará una muestra aleatoria del objeto uissurv que será un data.frame que contendrá la información necesaria para construir posteriormente un modelo de Cox: id age beck hercoc ivhx ndrugtx race treat site los time censor 284 33 17.00 Cocain only Never 3 Non-White Short A 55 115 1 623 36 9.00 Neither Previous NA White Short B 37 51 1 98 26 11.00 Heroin & Cocain Recent 1 White Long A 40 73 1 104 25 6.00 Heroin only Recent 5 White Long A 90 655 0 5 24 5.00 Heroin only Never 5 Non-White Long A 173 551 0 33 37 9.45 Neither Recent 1 White Short A 90 259 1 454 43 9.00 Heroin & Cocain Recent 0 Non-White Long B 6 6 1 499 30 8.40 Cocain only Previous 40 White Short B 36 36 1 Tabla 4.1: Muestra de datos del fichero uissurv. CAPÍTULO 4. MODELO DE COX EN R 42 4.3. APLICACIONES Una vez comprendida la estructura de los datos, se procede a construir un modelo de Cox para la variable tiempo de supervivencia frente al resto de variables que serán las variables explicativas del modelo: library(survival) mod_cox_indp_t <- coxph( Surv(time, censor) ~ age + beck + hercoc + ivhx + ndrugtx + race + treat + site + los, data = uissurv ) A continuación, se adjunta la información que recoge el objeto coxph anterior: mod_cox_indp_t ## Call: ## coxph(formula = Surv(time, censor) ~ age + beck + hercoc + ivhx + ## ndrugtx + race + treat + site + los, data = uissurv) ## ## coef exp(coef) se(coef) z p ## age -0.0209245 0.9792929 0.0082448 -2.538 0.01115 ## beck 0.0053053 1.0053194 0.0049233 1.078 0.28122 ## hercocHeroin & Cocain 0.0143785 1.0144823 0.1687127 0.085 0.93208 ## hercocHeroin only 0.1660367 1.1806165 0.1646092 1.009 0.31313 ## hercocNeither 0.0125511 1.0126302 0.1228395 0.102 0.91862 ## ivhxPrevious 0.1882694 1.2071587 0.1387700 1.357 0.17488 ## ivhxRecent 0.3986325 1.4897860 0.1482168 2.690 0.00716 ## ndrugtx 0.0262922 1.0266409 0.0085909 3.060 0.00221 ## raceWhite 0.3063290 1.3584291 0.1160987 2.639 0.00833 ## treatShort -0.1381856 0.8709370 0.0969330 -1.426 0.15399 ## siteB 0.4437656 1.5585651 0.1112148 3.990 6.6e-05 ## los -0.0094204 0.9906238 0.0008215 -11.467 < 2e-16 ## ## Likelihood ratio test=193.3 on 12 df, p=< 2.2e-16 ## n= 575, number of events= 464 ## (53 observations deleted due to missingness) Tal y como se especifica en la presentación de la función coxph() del Apartado 4.3.1, en la salida anterior se pueden consultar los coeficientes estimados (columna coef ) de cada una de las variables explicativas consideradas en el modelo y evaluar su significación con la ayuda de la columna p. Es destacable que, dado que en el conjunto de datos considerado existen variables categóricas con kmodalidades, se muestra un coeficiente estimado para cada una de las k  1modalidades de la variable en cuestión. A la vista del resultado obtenido, se sabe que el modelo se ha creado a partir de 575 observaciones de las cuales, en 464 de ellas, se ha producido el evento de interés (recaída en las drogas del individuo) y53 observaciones no se han usado por poseer valores perdidos. En el caso de que se quiera tener una salida más detallada del modelo creado, se usará la orden summary(): CAPÍTULO 4. MODELO DE COX EN R 43 4.3. APLICACIONES summary(mod_cox_indp_t) ## Call: ## coxph(formula = Surv(time, censor) ~ age + beck + hercoc + ivhx + ## ndrugtx + race + treat + site + los, data = uissurv) ## ## n= 575, number of events= 464 ## (53 observations deleted due to missingness) ## ## coef exp(coef) se(coef) z Pr(>|z|) ## age -0.0209245 0.9792929 0.0082448 -2.538 0.01115 * ## beck 0.0053053 1.0053194 0.0049233 1.078 0.28122 ## hercocHeroin & Cocain 0.0143785 1.0144823 0.1687127 0.085 0.93208 ## hercocHeroin only 0.1660367 1.1806165 0.1646092 1.009 0.31313 ## hercocNeither 0.0125511 1.0126302 0.1228395 0.102 0.91862 ## ivhxPrevious 0.1882694 1.2071587 0.1387700 1.357 0.17488 ## ivhxRecent 0.3986325 1.4897860 0.1482168 2.690 0.00716 ** ## ndrugtx 0.0262922 1.0266409 0.0085909 3.060 0.00221 ** ## raceWhite 0.3063290 1.3584291 0.1160987 2.639 0.00833 ** ## treatShort -0.1381856 0.8709370 0.0969330 -1.426 0.15399 ## siteB 0.4437656 1.5585651 0.1112148 3.990 6.6e-05 *** ## los -0.0094204 0.9906238 0.0008215 -11.467 < 2e-16 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## exp(coef) exp(-coef) lower .95 upper .95 ## age 0.9793 1.0211 0.9636 0.9952 ## beck 1.0053 0.9947 0.9957 1.0151 ## hercocHeroin & Cocain 1.0145 0.9857 0.7288 1.4121 ## hercocHeroin only 1.1806 0.8470 0.8551 1.6301 ## hercocNeither 1.0126 0.9875 0.7960 1.2883 ## ivhxPrevious 1.2072 0.8284 0.9197 1.5845 ## ivhxRecent 1.4898 0.6712 1.1142 1.9920 ## ndrugtx 1.0266 0.9741 1.0095 1.0441 ## raceWhite 1.3584 0.7361 1.0820 1.7055 ## treatShort 0.8709 1.1482 0.7202 1.0532 ## siteB 1.5586 0.6416 1.2533 1.9382 ## los 0.9906 1.0095 0.9890 0.9922 ## ## Concordance= 0.741 (se = 0.011 ) ## Likelihood ratio test= 193.3 on 12 df, p=<2e-16 ## Wald test = 170.5 on 12 df, p=<2e-16 ## Score (logrank) test = 176.2 on 12 df, p=<2e-16 A lo largo de las siguientes secciones, se irá analizando el modelo de Cox creado con las covariantes independientes del tiempo contenidas en el conjunto de datos bajo estudio. En primer lugar, se realizará la diagnosis del modelo de forma que se pueda ajustar un modelo que no viole las hipótesis de riesgos tanto global como individualmente. Una vez obtenido dicho modelo, se analizarán los parámetros estimados y los contrastes de hipótesis asociados a estos; y se terminará con la evaluación del mismo. 4.3.1.1. Diagnosis del modelo Para obtener un modelo de Cox válido y con capacidad predictiva es fundamental que las variables involucradas en él cumplan las hipótesis de riesgos proporcionales tanto de CAPÍTULO 4. MODELO DE COX EN R 44 4.3. APLICACIONES manera conjunta como individual. El estudio de si la hipótesis de riesgos proporcionales se verifica o no en el modelo se puede realizar de manera gráfica y analítica mediante el test de bondad de ajuste (GOF). Tras esta breve introducción, se presentan ambos enfoques y sus correspondientes conclusiones. A continuación, se procede a analizar los diferentes test de bondad de ajuste asociados a cada una de las variables predictoras del modelo de Cox creado. Para extraer dichos resultados de manera analítica, se usará la orden cox.zph() del paquete survival con una transformación de los tiempos de supervivencia de rango (transform = “rank”). En la siguiente tabla, se puede consultar el resultado obtenido a través de la orden en cuestión. Dicha tabla contiene para cada una de las variables predictoras consideradas en el modelo, el valor del estadístico χ2asociado al contraste realizado junto a sus grados de libertad y su p-valor: cox.zph(mod_cox_indp_t, transform = "rank") Variable chisq df p age 0.397 1 0.5288 beck 2.979 1 0.0843 hercoc 1.179 3 0.7581 ivhx 0.326 2 0.8496 ndrugtx 0.346 1 0.5564 race 0.503 1 0.4781 treat 8.526 1 0.0035 site 2.414 1 0.1202 los 207.403 1 <2e-16 GLOBAL 213.953 12 <2e-16 Tabla 4.2: Estudio de la hipótesis de riesgos proporcionales tanto de manera individual como global para el modelo de Cox creado. A la vista de los resultados obtenidos, dado que los p-valores asociados a las variables predictoras treat ylos son menores que 0.05, se concluye que existen evidencias significativas para determinar que dichas variables violan la hipótesis de riesgos proporcionales. Es más, el modelo de Cox creado no es válido dado que globalmente se viola también dicha hipótesis. Para intentar solucionar este problema, se deben eliminar dichas variables, replantear el modelo y analizarlo de nuevo. Otra forma de obtener la tabla anterior es mediante la orden que se adjunta a continuación: cox.zph(mod_cox_indp_t, transform = "rank")$table El estudio analítico de la hipótesis de riegos proporcionales se puede hacer de manera gráfica a través del estudio de los residuos de Schoenfeld tal y como se especificó en el CAPÍTULO 4. MODELO DE COX EN R 45 4.3. APLICACIONES Apartado 2.5.2. Antes de proceder a hacer las gráficas, se especificará la manera en la que se pueden extraer los residuos de Schoenfeld asociados a cada uno de los individuos de la muestra y a cada variable predictora considerada. Para ello, se usará la función residuals() si se quieren obtener los residuos sin escalar y cox.zph()$y en caso contrario: residuals(mod_cox_indp_t, type = c("schoenfeld")) cox.zph(mod_cox_indp_t, transform = "rank")$y time age beck hercoc ivhx ndrugtx race treat site los 4 0.0823 -0.0916 29.973 1.7194 -0.0126 1.259 -2.8373 0.4223 -0.0341 4 0.2277 0.0220 -0.2571 -10.1006 -0.069 3.1586 -2.5264 -1.2056 -0.0264 6 -0.1071 -0.1726 -2.7891 3.3803 -0.0205 -4.7547 2.0453 -0.9294 -0.0194 6 0.2863 -0.1350 -10.6887 14.3498 -0.3027 -4.4125 -3.4116 6.223 -0.0372 7 0.1487 0.057 33.2189 -6.2191 0.4165 0.9228 1.7559 -0.3232 -0.0264 7 -0.15683 0.0809 49.0868 -13.681 0.0496 1.8397 -1.986 -0.9337 -0.0289 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 502 0.0273 0.008 -16.256 1.9374 0.12689 1.2539 -1.259 1.6857 0.0643 516 0.2697 0.0981 -11.126 -1.7983 -0.0675 2.8468 -1.8483 -1.6535 0.0162 519 -0.1844 -0.0224 2.1318 2.330 -0.0213 -4.2965 -1.7249 3.7895 0.0016 559 -0.2312 -0.1011 0.0223 -0.989 0.0421 2.216 -2.338 -0.574 -0.027 568 0.1275 -0.15836 29.0773 -1.8809 0.4715 0.4334 1.7649 -0.0166 -0.0148 659 -0.1485 -0.0033 -2.5585 3.9765 0.0151 -3.5266 -2.348 0.8089 -0.0059 Tabla 4.3: Residuos escalados de Schoenfeld asociados al modelo de Cox creado. Una vez mostrada la forma en la que se pueden extraer los residuos de Schoenfeld asociados a cada una de las variables explicativas de un modelo de Cox, se procede a realizar unas gráficas de estos para su posterior análisis. Para ello, se puede usar la función plot(cox.zph()) o directamente la función ggcoxdiagnostics() aplicada al objeto coxph que contiene el modelo de Cox. La primera de las opciones tiene la peculiaridad de que las gráficas que hace para las variables de tipo factor no se realizan para cada una de las categorías consideradas en el modelo como variables predictoras, sino que hace una única gráfica que resume todas las categorías. En cambio, en el caso de la función ggcoxdiagnostics() se dibujan tantas gráficas como coeficientes se han estimado a la hora de ajustar el modelo de Cox creado. A continuación, se presentan los resultados obtenidos con cada una de las opciones estudiadas y las conclusiones obtenidas de ellas: par(mfrow = c(1,2)) plot(cox.zph(mod_cox_indp_t, transform = "rank"), col = "red") par(mfrow = c(1,1)) CAPÍTULO 4. MODELO DE COX EN R 46 4.3. APLICACIONES Figura 4.1: Gráficas de los residuos de Schoenfeld con la orden plot(cox.zph()). CAPÍTULO 4. MODELO DE COX EN R 47 4.3. APLICACIONES Tal y como se especificó en el Apartado 2.5.2, se dirá que una variable explicativa del modelo cumple la hipótesis de riesgos proporcionales cuando su gráfica de los residuos frente al tiempo no muestra ningún patrón o, la curva que se dibuja en rojo es de pendiente nula. Teniendo en cuenta dicha apreciación, se podría considerar que efectivamente las variables treat ylos no cumplen la hipótesis de riesgos proporcionales, aunque no de manera muy clara. Dado que las conclusiones de manera gráfica no son muy certeras, es mejor hacer uso de un enfoque más objetivo como son los contrastes de hipótesis. Tras el análisis gráfico anterior, se procede a mostrar la salida de la segunda de las opciones consideradas para las gráficas de los residuos de Schoenfeld: library(survminer) ggcoxdiagnostics(mod_cox_indp_t, type = "schoenfeld" ) ndrugtx raceWhite siteB treatShort hercocNeither ivhxPrevious ivhxRecent los age beck hercocHeroin & Cocain hercocHeroin only 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0 100200300400 0.0 0.3 0.6 0.9 −100 0 100 200 −0.50 −0.25 0.00 0.25 0.50 0.0 0.3 0.6 −0.50 −0.25 0.00 0.25 0.50 0.75 −0.2 0.0 0.2 0.4 0.6 −20 −10 0 10 20 0.0 0.3 0.6 −0.6 −0.3 0.0 0.3 −10 0 10 20 0.0 0.5 0 10 20 30 Observation Id Residuals (type = schoenfeld) Fuente propia Figura 4.2: Gráficas de los residuos de Schoenfeld con la orden ggcoxdiagnostics() de la librería survminer. Al igual que se ha comentado anteriormente, el análisis de las gráficas de los residuos anteriores para cada una de las variables explicativas consideradas en el modelo no es objetivo por lo que no es aconsejable sacar muchas conclusiones de manera gráfica. Otra manera de conseguir gráficas de los residuos de Schoenfeld es a través de la función ggcoxzph(). Sin embargo, cuando el modelo tiene muchas variables, las gráficas obtenidas no se ven correctamente por lo que solo se adjunta la orden que permite conseguir dichas gráficas: CAPÍTULO 4. MODELO DE COX EN R 48 4.3. APLICACIONES ggcoxzph(cox.zph(mod_cox_indp_t)) Asimismo, en el Apartado 2.5.1 se expuso que otra manera gráfica de diagnosticar si una variable cumple o no la hipótesis de riesgos proporcionales es mediante la comparación de las curvas log-log de supervivencia. Este método es útil para las variables categóricas o de tipo factor y se concluirá que dicha variable cumple la hipótesis de riesgos proporcionales cuando las correspondientes curvas son aproximadamente paralelas. A continuación, para mostrar el uso de esta técnica se adjunta un ejemplo realizado sobre la variable race y el uso de la función ggsurvplot() del paquete survminer sobre un objeto de tipo survfit: ggsurvplot(survfit(Surv(time, censor) ~ race, data = uissurv), fun = "cloglog", xlab = "t (días)", ylab = "Curva log-log de supervivencia", lab = c(10,10,7), size = 1,legend.title = "Race", palette = c("chartreuse","cyan4"), conf.int = FALSE, legend.labs = c("White","Non-White"), ggtheme = theme_grey(), font.caption = c(10,"italic","black"), font.x = c(12,"bold.italic","black"), font.y = c(12,"bold.italic","black") ) + ++++++++++++++++++++++++++ +++ + + ++++++++++++++++++++++ +++++ +++++++++++++++++++++ + −6 −4 −2 0 1 10 100 1000 10000 t (días) Curva log−log de supervivencia Race ++ White Non−White Fuente propia Figura 4.3: Comparación de las curvas log-log de supervivencia para la variable predictora race. CAPÍTULO 4. MODELO DE COX EN R 49 4.3. APLICACIONES summary(mod_cox_indp_t_ajus)$coef[, c(4,5)] Variable z Pr(>|z|) age -3.4946910 0.0004746109 *** beck 1.6932172 0.0904141405 . hercocHeroin & Cocain 0.6544510 0.5128212840 hercocHeroin only 0.9058276 0.3650271098 hercocNeither 0.8032623 0.4218231408 ivhxPrevious 1.2877376 0.1978373357 ivhxRecent 1.9978879 0.0457288101 * ndrugtx 3.3094180 0.0009349015 *** raceWhite 1.8524678 0.0639586787 . siteB -0.7105864 0.4773405668 Tabla 4.8: Contrastes de hipótesis individuales asociados al modelo de Cox creado. A la vista de los resultados obtenidos, se concluye que las variables significativas del modelo son age,ivhxRecent yndrugtx dado que el p-valor asociado a ellas es menor que 0.05; y las variables no significativas son beck,hercocHeroin & Cocain,hercocHeroin only, hercocNeither,ivhxPrevious,raceWhite ysiteB. Estas últimas variables no significativas podrían eliminarse del modelo y este no sufriría ningún cambio. Además, en el caso de que se considerase un nivel de significación del 10%, también pasarían a ser significativas las variables beck yraceWhite. Una vez realizados y analizados los contrastes de hipótesis individuales, es momento de analizar lo que ocurre con el contraste de hipótesis conjunto de todas las covariantes usadas en el modelo: " H0:β  p β10, . . . , β100 q 1  0 H1:βj  βj0para algún j  1,...,10. Para dicho cometido, existen tres formas de resolverlo: mediante el test de razón de verosimilitud, el test de Wald o el test Score, los cuales son asintóticamente equivalentes. Tanto el estadístico asociado a cada uno de dichos test como su p-valor, se pueden consultar en el summary() del modelo creado y se adjunta a continuación: Figura 4.5: Captura de la salida summary(mod_cox_indp_t_ajus) para el estudio del contraste de hipótesis conjunto del modelo de Cox. Dado que los p-valores obtenidos son todos significativos, i.e., menores que 0.05, se concluye que el modelo creado permite explicar la variable tiempo hasta la recaída en CAPÍTULO 4. MODELO DE COX EN R 56 4.3. APLICACIONES las drogas. Además, el valor de los estadísticos χ2asociados a los test en estudio junto a sus grados de libertad y p-valores se pueden consultar directamente con las siguientes órdenes: Test de razón de verosimilitud: summary(mod_cox_indp_t_ajus)$logtest test df pvalue 4.160053e+01 1.000000e+01 8.830632e-06 Tabla 4.9: Test de razón de verosimilitud del modelo de Cox creado. Test de Wald: summary(mod_cox_indp_t_ajus)$waldtest test df pvalue 4.26600e+01 1.00000e+01 5.72741e-06 Tabla 4.10: Test de Wald del modelo de Cox creado. Test score: summary(mod_cox_indp_t_ajus)$sctest test df pvalue 4.309274e+01 1.000000e+01 4.786700e-06 Tabla 4.11: Test score del modelo de Cox creado. Por último, es destacable que todo el estudio analítico del modelo que se ha realizado a través de la estimación de los parámetros y de los contrastes de hipótesis individuales y conjunto, se puede resumir gráficamente con la orden ggforest() de la librería survminer. En dicha salida y, tal y como se especificó en el Apartado 4.3.1, se puede consultar el coeficiente estimado de cada una de las variables explicativas junto con su intervalo de confianza, el número de observaciones asociadas a cada una de ellas o los p-valores asociados a cada uno de los parámetros, entre otra información relevante. La orden en cuestión, se muestra en la siguiente línea de código: CAPÍTULO 4. MODELO DE COX EN R 57 4.3. APLICACIONES ggforest(mod_cox_indp_t_ajus, main = "Resumen del modelo") site race ndrugtx ivhx hercoc beck age B A White Non−White Recent Previous Never Neither Heroin only Heroin & Cocain Cocain only (N=184) (N=444) (N=467) (N=155) (N=628) (N=262) (N=115) (N=233) (N=207) (N=114) (N=111) (N=178) (N=628) (N=628) 0.93 reference 1.24 reference 1.03 1.34 1.20 reference 1.10 1.16 1.11 reference 1.01 0.97 (0.75 − 1.14) (0.99 − 1.56) (1.01 − 1.04) (1.01 − 1.78) (0.91 − 1.57) (0.87 − 1.40) (0.84 − 1.59) (0.81 − 1.54) (1.00 − 1.02) (0.96 − 0.99) 0.477 0.064 <0.001 *** 0.046 * 0.198 0.422 0.365 0.513 0.09 <0.001 *** # Events: 464; Global p−value (Log−Rank): 8.8306e−06 AIC: 5304.7; Concordance Index: 0.59 0.8 1 1.2 1.4 1.6 1.8 Resumen del modelo de Cox Figura 4.6: Resumen gráfico del análisis del modelo de Cox creado con la función ggforest() del paquete survminer. En el resumen anterior, se observan cinco columnas bien diferenciadas: la primera de ellas contiene el nombre de las variables explicativas consideradas en el modelo de Cox creado; la segunda, especifica los diferentes valores o categorías que toman las variables de tipo factor junto con el número de individuos u observaciones que presentan cada una de ellas; en la tercera, se muestra la exponencial del coeficiente estimado de cada variable junto con su intervalo de confianza al 95% (valor por defecto); en la cuarta, se dibuja dicha estimación con su intervalo de confianza; y por último, en la quinta, se presenta el p-valor asociado a las estimaciones de los parámetros realizadas. Además, en el extremo inferior izquierdo también se puede consultar un resumen global del modelo creado que incluye el p-valor del contraste conjunto o el número de sujetos sobre los que se ha producido el evento de interés. Dado que con esta gráfica se está analizando el mismo modelo de Cox construido al principio de la sección, las conclusiones que se obtienen consultando esta última gráfica son las mismas que se han ido destacando a lo largo de los Apartados 4.3.1.2 y4.3.1.3. 4.3.1.4. Evaluación del modelo Para evaluar la calidad de un modelo de Cox, se usarán las medidas expuestas en el Apartado 2.6 que incluyen el coeficiente sugerido por Nagelkerke (1991), el propuesto CAPÍTULO 4. MODELO DE COX EN R 58 4.3. APLICACIONES por O’Quigley et al. (2005) y el sugerido por Royston (2006). Todos estos coeficientes se pueden obtener con la orden rsq() de la librería survMisc. A continuación, se adjuntan dichos coeficientes para el nuevo modelo de Cox creado: library(survMisc) rsq(mod_cox_indp_t_ajus, sigD = NULL) ## $cod ## [1] 0.06979357 ## ## $mer ## [1] 0.08575466 ## ## $mev ## [1] 0.05394637 En vista a los resultados obtenidos, no parece que la calidad del modelo creado sea excesivamente buena para la predicción del tiempo de supervivencia pues, los coeficientes obtenidos son muy cercanos al 0e interesaría que fueran próximos a 1. No obstante, el modelo obtenido es válido para determinar aquellas variables que influyen en dicho tiempo y el sentido de dicha influencia. 4.3.1.5. Selección de variables del modelo En esta sección, se realizará un estudio de selección de variables en un modelo de Cox con la librería glmnet. Para llevarlo a cabo, se considerará el modelo de Cox que incluye todas las posibles variables explicativas (excepto la variable los) del conjunto de datos uissurv con el que se ha estado trabajando en las secciones anteriores. En primer lugar, se hará un estudio de los λ’s óptimos (λm´ın yλ1se) con la ayuda de la función cv.glmnet() y una vez encontrados dichos valores, se crearán los modelos correspondientes con la orden glmnet(). Antes de proceder a la programación del modelo, es preciso resaltar que las funciones anteriores operan con estructuras de datos de tipo matrix en las que todas sus variables son numéricas y no existen valores perdidos. Debido a esta peculiaridad, es necesario convertir las variables de tipo factor cuyas categorías son cadenas de caracteres en otras variables que toman solo valores numéricas. Como consecuencia de ello y dado que el conjunto de datos uissurv a utilizar posee variables de tipo factor sin codificación numérica y con datos perdidos, se procederá a modificarlo para que sea posible su posterior aplicación con las funciones citadas anteriormente. En primer lugar, se procede a eliminar aquellas observaciones o individuos que poseen valores perdidos en cualquiera de sus variables. El nuevo conjunto de datos se renombrará como uissurv_glmnet: uissurv_glmnet <- na.omit(uissurv) Tras eliminar los valores perdidos, se procede a recodificar cada una de las categorías de las variables de tipo factor (hercoc,ivhx,race,treat ysite) para que cumplan el requisito CAPÍTULO 4. MODELO DE COX EN R 59 4.3. APLICACIONES de ser variables numéricas. Para ello, es necesario el uso de la función recode() de la librería dplyr (Wickham et al. (2023b)): library(dplyr) uissurv_glmnet$hercoc <- recode(uissurv_glmnet$hercoc, "Heroin & Cocain" =1, "Heroin only" =2,"Cocain only" =3, "Neither" =4 ) uissurv_glmnet$ivhx <- recode(uissurv_glmnet$ivhx, "Never" =1, "Previous" =2, "Recent" =3 ) uissurv_glmnet$race <- recode(uissurv_glmnet$race, "White" =0, "Non-White" =1 ) uissurv_glmnet$treat <- recode(uissurv_glmnet$treat, "Short" =0, "Long" =1 ) uissurv_glmnet$site <- recode(uissurv_glmnet$site, "A" =0, "B" =1 ) Una vez preparadas las variables categóricas, se procede a crear dos nuevos conjuntos de datos: el primero de ellos, x, contiene las variables explicativas del modelo y el segundo, y, incluye las variables tiempo y censura. Nótese también que para trabajar un modelo de Cox con la librería glmnet, es necesario que la columna que contiene los datos de censura se llame status. Por lo tanto, será necesario aplicar dicho cambio también al conjunto de datos utilizado. A continuación, se presenta el código necesario para terminar la preparación de los datos: x<- uissurv_glmnet[, c( "age","beck","hercoc","ivhx","ndrugtx", "race","treat","site" )] y<- uissurv_glmnet[, c("time","censor")] colnames(y)[2]<- "status" En el momento en el que se tienen los datos preparados, se procede a la construcción del ajuste con la librería glmnet y la búsqueda de los λ’s óptimos. Recuérdese que para ello, se usará la función cv.glmnet() con el argumento family = “cox”,type.measure = “deviance” para que el ajuste del modelo de Cox sea a través de la verosimilitud parcial, mientras CAPÍTULO 4. MODELO DE COX EN R 60 4.3. APLICACIONES que el resto de argumentos que toma por defecto (alpha = 119 ynlambda = 100). De esta forma, el estudio que se realizará será una selección de variables por el método Lasso para así conseguir que algunos coeficientes estimados sean nulos. Se empleará también una semilla para que los resultados obtenidos mediante validación cruzada sean reproducibles posteriormente. Asimismo, se añadirá una gráfica en la que se puede observar la evolución del logaritmo del λutilizado en función de la verosimilitud parcial. A continuación, se procede a realizar el estudio: library(glmnet) set.seed(1) cvfit <- cv.glmnet(as.matrix(x), as.matrix(y), family = "cox", type.measure = "deviance" ) Se adjunta la gráfica en cuestión en la que también pueden consultarse en el extremo superior el número de variables consideradas en los ajustes realizados con los diferentes λ’s y con líneas verticales discontinuas se indican los valores óptimos de λ: plot(cvfit) −6 −5 −4 −3 −2 10.90 10.95 11.00 11.05 Log(λ) Partial Likelihood Deviance 888777777766652 Fuente propia Figura 4.7: Evolución del valor de log p λ q en función de la verosimilitud parcial obtenida con la función cv.glmnet(). Se observa que para ambos parámetros óptimos se ha producido una selección de variables, consecuencia directa con el hecho de que se está realizando un ajuste usando el método Lasso. Una vez consultada la gráfica, se extraen los valores óptimos de λ: 19El valor del argumento alpha determina el método de regularización deseado y toma valores entre 0 y1. En el caso del valor 0, aplica el método Ridge; y el 1, realiza el método Lasso. CAPÍTULO 4. MODELO DE COX EN R 61 4.3. APLICACIONES Valor obtenido de λm´ın: cvfit$lambda.min ## [1] 0.01223744 Valor obtenido de λ1se: cvfit$lambda.1se ## [1] 0.09475009 Tras la extracción de los parámetros λóptimos, se procede a construir los modelos correspondientes con la orden glmnet(). El primer modelo se creará con el valor del λm´ın; y el segundo, con el valor del λ1se: 1. Modelo usando λm´ın: fit <- glmnet(as.matrix(x), as.matrix(y), family = "cox",lambda = cvfit$lambda.min ) Una vez definido el modelo correspondiente con su λm´ın, se estudia el valor de sus coeficientes estimados: fit$beta Variable s0 age -0.02421305 beck 0.00664753 hercoc . ivhx 0.15320376 ndrugtx 0.02529730 race -0.17576828 treat -0.19918352 site -0.05699427 Tabla 4.12: Estimación de los coeficientes asociados a cada una de las variables explicativas usando λm´ın en el método Lasso. Tal y como se puede observar, con este valor de λse obtiene un modelo de Cox en el que el coeficiente de la variable hercoc es nula. De esta forma, se ha conseguido un nuevo modelo de Cox más simple. Una vez obtenido este modelo, se debería realizar un estudio de diagnosis del modelo para de manera objetiva, saber si el modelo creado es válido o no. CAPÍTULO 4. MODELO DE COX EN R 62 4.3. APLICACIONES 2. Modelo usando λ1se: fit1 <- glmnet(as.matrix(x), as.matrix(y), family = "cox",lambda = cvfit$lambda.1se ) Al igual que se ha realizado en el modelo anterior, se estudian cada uno de los coeficientes estimados: fit1$beta Variable s0 age . beck . hercoc . ivhx 0.048330501 ndrugtx 0.006150463 race 0.006891106 treat . site . Tabla 4.13: Estimación de los coeficientes asociados a cada una de las variables explicativas usando λ1se en el método Lasso. A la vista de los coeficientes, se ha obtenido un nuevo modelo en el que solo se consideran tres variables explicativas no nulas (ivhx,ndrugtx yrace) para modelizar el tiempo hasta la recaída en las drogas. Al igual que se especificó en el caso anterior, para garantizar que el modelo obtenido es válido, será necesario hacer un estudio de diagnosis de dicho modelo. En este caso, el modelo obtenido mediante el método Lasso con λ1se es mucho más simple que el creado con λm´ın. Por último, se realiza una gráfica en la que se puede observar la evolución de los valores de los coeficientes de las variables explicativas consideradas frente a la norma L1 a medida que el valor de λvaría. Además, en el eje superior de esta se puede consultar el número de coeficientes distintos de 0: plot(glmnet(as.matrix(x), as.matrix(y), family = "cox")) CAPÍTULO 4. MODELO DE COX EN R 63 4.3. APLICACIONES 0.0 0.2 0.4 0.6 −0.2 −0.1 0.0 0.1 L1 Norm Coefficients 0 6 6 7 Fuente propia Figura 4.8: Evolución del valor de los coeficientes de las variables predictoras frente a la norma L1 en el método Lasso. 4.3.1.6. Otras consideraciones del modelo Para terminar con el análisis completo de un modelo de Cox, se añaden órdenes que pueden ser de utilidad a la hora de programar un modelo de regresión de este tipo. A lo largo de esta sección, el modelo que se usará es el mod_cox_indp_t_ajus que se ajustó en apartados anteriores. En primer lugar, se muestra la estimación de la matriz de varianzas y covarianzas de las variables predictoras usando la orden coxph()$var, la cual permite evaluar la relación existente entre las variables explicativas del modelo: mod_cox_indp_t_ajus$var ## [,1] [,2] [,3] [,4] [,5] ## [1,] 6.618911e-05 2.423426e-06 -6.345488e-05 -8.855666e-05 -3.738056e-05 ## [2,] 2.423426e-06 2.451767e-05 2.693840e-05 2.899191e-05 4.293907e-05 ## [3,] -6.345488e-05 2.693840e-05 2.708693e-02 1.556921e-02 8.178327e-03 ## [4,] -8.855666e-05 2.899191e-05 1.556921e-02 2.647426e-02 8.290126e-03 ## [5,] -3.738056e-05 4.293907e-05 8.178327e-03 8.290126e-03 1.473712e-02 ## [6,] -2.031930e-04 -5.588681e-05 -1.889368e-03 -1.904810e-03 -2.545384e-03 ## [7,] -2.520241e-04 -7.524165e-05 -1.148694e-02 -1.151066e-02 -1.443276e-03 ## [8,] -1.173992e-05 1.579762e-06 1.749931e-05 6.985647e-05 4.541880e-05 ## [9,] 7.498269e-05 2.613260e-05 9.140565e-04 -8.209515e-04 -1.366308e-04 ## [10,] 1.381050e-05 2.908462e-05 7.206944e-04 5.188565e-04 -4.136942e-04 ## [,6] [,7] [,8] [,9] [,10] ## [1,] -2.031930e-04 -2.520241e-04 -1.173992e-05 7.498269e-05 1.381050e-05 ## [2,] -5.588681e-05 -7.524165e-05 1.579762e-06 2.613260e-05 2.908462e-05 ## [3,] -1.889368e-03 -1.148694e-02 1.749931e-05 9.140565e-04 7.206944e-04 ## [4,] -1.904810e-03 -1.151066e-02 6.985647e-05 -8.209515e-04 5.188565e-04 ## [5,] -2.545384e-03 -1.443276e-03 4.541880e-05 -1.366308e-04 -4.136942e-04 ## [6,] 1.926279e-02 8.321383e-03 -1.205006e-04 -1.520755e-03 -1.095199e-03 ## [7,] 8.321383e-03 2.141026e-02 -2.283340e-04 -3.168170e-03 1.992492e-03 ## [8,] -1.205006e-04 -2.283340e-04 6.819688e-05 -6.367003e-05 1.111069e-05 CAPÍTULO 4. MODELO DE COX EN R 64 4.3. APLICACIONES ## [9,] -1.520755e-03 -3.168170e-03 -6.367003e-05 1.351003e-02 -2.506650e-03 ## [10,] -1.095199e-03 1.992492e-03 1.111069e-05 -2.506650e-03 1.172941e-02 Cuando un modelo está ajustado es posible predecir el tiempo hasta la recaída en las drogas de nuevas observaciones y para ello, se usará la orden predict(). Nótese que para predecir, ajustar y evaluar el rendimiento del modelo, así como evitar el sobreajuste, hubiera sido interesante haber considerado dos conjuntos de datos, el de entrenamiento y el test. A pesar de ello, se procede a hacer una predicción sobre un conjunto de datos generado aleatoriamente: set.seed(12) age <- sample(18:60,10,replace = TRUE) beck <- sample(0:54,10,replace = TRUE) hercoc <- sample(c( "Heroin & Cocain","Heroin only", "Cocain only","Neither" ), 10,replace = T) ivhx <- sample(c("Never","Recent","Previous"), 10,replace = T) ndrugtx <- sample(0:40,10,replace = TRUE) race <- sample(c("Non-White","White"), 10,replace = T) site <- sample(c("A","B"), 10,replace = T) datos_pred <- data.frame( age, beck, hercoc, ivhx, ndrugtx, race, site ) Una vez generados los datos, se procede a la predicción de estos: predict(mod_cox_indp_t_ajus, newdata = datos_pred, type = "lp") Observación Predicción 1 0.6017720 2 0.3122243 3 1.3195960 4 0.5770842 5 0.9976317 6 1.2379449 7 0.5594878 8 0.6526373 9 0.8530893 10 0.7766950 Tabla 4.14: Predicción sobre un nuevo conjunto de datos generado de manera aleatoria. CAPÍTULO 4. MODELO DE COX EN R 65 4.3. APLICACIONES Al aplicar los argumentos tdc para crear las nuevas covariantes dependientes del tiempo, se han generado además dos variables nuevas, tstart ytstop, a partir de las variables time yday. En dichas variables se encuentra la información de los nuevos intervalos de estudio donde varían las covariantes dependientes del tiempo y cuya estructura es de la forma p tstart, tstop s . De esta forma, se ha obtenido un único conjunto de datos con la información necesaria para realizar un estudio de análisis de supervivencia con algunas variables dependientes del tiempo. A continuación, se adjunta una muestra del conjunto de datos creado donde la variable id indica el paciente al que corresponden las mediciones realizadas y donde se han seleccionado las variables más relevantes a lo largo de la construcción del nuevo conjunto de datos: id status ascites bili albumin protime alk.phos tstart tstop death 1 2 1 14.5 2.60 12.2 1718 0 192 0 1 2 1 21.3 2.94 11.2 1612 192 400 2 84 0 0 0.4 3.76 11.2 1345 0 184 0 84 0 0 0.6 3.6 10.5 784 184 371 0 84 0 0 0.4 3.62 10.4 771 371 4032 0 Tabla 4.18: Muestra de datos del conjunto creado pbc_unif con la función tmerge(). Dado que el conjunto de datos creado con la función tmerge() es complejo y posee muchas variables, se ha estimado más conveniente usar el conjunto de datos heart (Crowley y Hu (1977)) de la librería survival para la construcción de un modelo de Cox con covariantes dependientes del tiempo. Dicho conjunto de datos estudia la supervivencia de los pacientes en lista de espera del programa de trasplantes cardíacos de la ciudad de Stanford y posee las siguientes variables: start ystop: tiempo de inicio y fin del intervalo de tiempo en estudio. event: estado en el intervalo asociado. Puede tomar los valores 0  vivo y 1  muerto. age: edad del paciente menos 48 años. year: año de aceptación en el programa. Medida expresada en años. surgery: variable que indica si el paciente ha tenido una cirugía de bypass previa o no. Los posibles valores son: 1si el paciente ha tenido cirugía previa o 0en caso contrario. Se trata de una variable dependiente del tiempo. transplant: variable binaria que especifica si el paciente ha recibido un transplante o no. Esta variable tomará el valor 1cuando el paciente ha recibido un transplante y0en caso contrario. Se trata de una variable dependiente del tiempo. id: identificador numérico de cada uno de los pacientes. Una vez presentadas cada una de las variables del conjunto de datos, se extraen los datos directamente desde R con la siguiente orden: CAPÍTULO 4. MODELO DE COX EN R 72 4.3. APLICACIONES data("heart",package = "survival") En la siguiente tabla se adjunta una muestra del data.frame en estudio: start stop event age year surgery transplant id 0 46 0 0.9253936 2.507871 0 0 36 46 100 1 0.9253936 2.507871 0 1 36 0 78 0 -6.617385 3.994524 1 0 59 78 916 0 -6.617385 3.994524 1 1 59 0 6 0 4.892539 5.284052 0 0 81 6 445 0 4.892539 5.284052 0 1 81 0 21 0 -0.9007529 3.340178 0 0 47 21 72 1 -0.9007529 3.340178 0 1 47 0 60 0 -1.746749 5.470226 0 0 87 60 110 1 -1.746749 5.470226 0 1 87 Tabla 4.19: Muestra de datos del conjunto heart contenido en el paquete survival. Tras la comprensión de la estructura de los datos, se procede a la construcción de un modelo de Cox para la variable tiempo de supervivencia expresada mediante intervalos de tiempo de la forma (tstart, tstop] frente al resto de variables disponibles en el conjunto de datos bajo estudio. Entre estas variables disponibles, se encuentran dos variables no dependientes del tiempo (age yyear) y dos variables dependientes del tiempo o variantes en el tiempo (surgery ytransplant). Todas estas variables actuarán como variables explicativas o predictoras en el modelo de regresión a construir: mod_cox_dep_t <- coxph(Surv(start, stop, event) ~ transplant + age + year + surgery, data = heart) A continuación, se muestra la información que posee el objeto coxph creado con la orden anterior: summary(mod_cox_dep_t) ## Call: ## coxph(formula = Surv(start, stop, event) ~ transplant + age + ## year + surgery, data = heart) ## ## n= 172, number of events= 75 ## ## coef exp(coef) se(coef) z Pr(>|z|) ## transplant1 -0.01025 0.98980 0.31375 -0.033 0.9739 ## age 0.02717 1.02754 0.01371 1.981 0.0476 * ## year -0.14635 0.86386 0.07047 -2.077 0.0378 * ## surgery -0.63721 0.52877 0.36723 -1.735 0.0827 . CAPÍTULO 4. MODELO DE COX EN R 73 4.3. APLICACIONES ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## exp(coef) exp(-coef) lower .95 upper .95 ## transplant1 0.9898 1.0103 0.5352 1.8307 ## age 1.0275 0.9732 1.0003 1.0555 ## year 0.8639 1.1576 0.7524 0.9918 ## surgery 0.5288 1.8912 0.2574 1.0860 ## ## Concordance= 0.636 (se = 0.033 ) ## Likelihood ratio test= 15.11 on 4 df, p=0.004 ## Wald test = 14.49 on 4 df, p=0.006 ## Score (logrank) test = 15.03 on 4 df, p=0.005 La información que se puede consultar en la salida obtenida es la misma que se resaltó en el Apartado 4.3.1 del presente trabajo. Entre ella, es destacable el número de observaciones utilizadas a la hora de construir el modelo, en este caso, 172; o el número de observaciones para las que se ha producido el evento de interés (muerte del paciente que se encuentra a la espera de recibir un transplante de corazón), 75. El resto de información relevante se analizará en las próximas secciones. A continuación, se procede al estudio detallado del modelo con la estimación de parámetros, la significación de estos y la evaluación de su calidad predictiva. Recuérdese que en el modelo de Cox con covariantes dependientes del tiempo no hace falta estudiar si se verifica la hipótesis de riesgos proporcionales pues, en general, el modelo de Cox extendido para este tipo de variables no satisface dicha hipótesis (véase el Apartado 3.4). 4.3.2.1. Estimación de parámetros El estudio de los parámetros o coeficientes estimados para cada una de las variables involucradas en el modelo se realizará con la siguiente orden: mod_cox_dep_t$coefficients Variable Estimación transplant1 -0.01025077 age 0.02716664 year -0.14634635 surgery -0.63720989 Tabla 4.20: Estimación de los coeficientes asociados a cada una de las variables explicativas del modelo creado con covariantes dependientes del tiempo. Para una interpretación de los coeficientes anteriores y el efecto que tienen sobre el tiempo de supervivencia modelizado en intervalos de tiempo, es necesario estudiar la exponencial de dichas estimaciones. A continuación, se adjunta una tabla en la que se pueden consultar dichas exponenciales para cada uno de los parámetros estimados junto a su intervalo de confianza al 95 %: CAPÍTULO 4. MODELO DE COX EN R 74 4.3. APLICACIONES summary(mod_cox_dep_t)$conf.int Variable exp(coef) lower .95 upper .95 transplant1 0.9898016 0.5351550 1.8306980 age 1.0275390 0.7524197 0.9918021 year 0.8638585 0.8066522 1.5376930 surgery 0.5287657 0.2574423 1.0860419 Tabla 4.21: Estimación de la exponencial de los coeficientes asociados a cada una de las variables explicativas del modelo creado con covariantes dependientes del tiempo junto a su intervalo de confianza al 95 %. Para evidenciar el efecto que tienen las variables involucradas en el modelo sobre la función de riesgo de muerte de los pacientes que se encuentran a la espera de un transplante cardíaco, se analizarán dos variables: Variable age: se ha obtenido que ˆ β2  0.02716664 yeˆ β2  1.0275390, de modo que, manteniendo constantes el resto de variables involucradas en el modelo de Cox creado, el aumento en un año de la edad de un paciente supone que la función riesgo de muerte de este aumente de manera aproximada en un 2.75 %, ya que p 1.0275390  1 q  100%  2.75%. Variable year: se tiene que ˆ β3   0.14634635 yeˆ β3  0.8638585, por lo que, si se mantienen constantes el resto de variables incluidas en el modelo de Cox construido, el entrar en el programa un año más tarde le supone a un paciente que la función riesgo de muerte disminuya en un factor de 0.8638585, que en promedio se corresponde con un 13.6 % aproximadamente, dado que p 1  0.8638585 q 100%  13.61%. Nótese que esta interpretación detallada de los coeficientes sobre el efecto que tienen en la función de riesgo de muerte asociada a los pacientes debe realizarse por un experto en la materia. Además, se puede estudiar el intervalo de confianza asociado a cada uno de los parámetros estimados de la siguiente manera: confint(mod_cox_dep_t) Variable coef 2.5 % 97.5 % transplant1 -0.01025077 -0.6251988771 0.604697332 age 0.02716664 0.0002874691 0.054045813 year -0.14634635 -0.2844610476 -0.008231644 surgery -0.63720989 -1.3569596167 0.082539837 Tabla 4.22: Estimación de los coeficientes asociados a cada una de las variables explica- CAPÍTULO 4. MODELO DE COX EN R 75 4.3. APLICACIONES tivas del modelo creado con covariantes dependientes del tiempo junto a su intervalo de confianza al 95%. En el estudio del modelo, es posible consultar también el número de iteraciones, en este caso 4, que ha tenido que realizar el algoritmo para que el método converja: mod_cox_dep_t$iter ## [1] 4 Asimismo, puede consultarse el valor de la log-verosimilitud del modelo junto al mismo valor bajo la hipótesis nula H0:β  0: mod_cox_dep_t$loglik ## [1] -298.1214 -290.5656 4.3.2.2. Contrastes de hipótesis Tras la estimación de los parámetros asociados a las variables predictoras consideradas, es necesario estudiar si estos son significativos o no para el modelo de regresión construido. Es decir, se deben resolver los siguientes contrastes de hipótesis para cada una de las variables explicativas consideradas en el modelo creado; i.e.: " H0:βj  0 H1:βj  0, para todo j  1,...,4donde, por simplicidad, los índices j  1yj  4hacen referencia a las variables dependientes del tiempo. Para alcanzar dicho objetivo, se analiza la columna Pr(>|z|) de la siguiente tabla que contiene los p-valores asociados a cada una de las estimaciones obtenidas mediante el estadístico de Wald que se adjunta también en la columna z: summary(mod_cox_dep_t)$coef[, c(4,5)] Variable z Pr(>|z|) transplant1 -0.03267128 0.97393672 age 1.98092553 0.04759963 * year -2.07677794 0.03782206 * surgery -1.73519821 0.08270570 . Tabla 4.23: Contrastes de hipótesis individuales asociados al modelo de Cox creado con covariantes dependientes del tiempo. CAPÍTULO 4. MODELO DE COX EN R 76 4.3. APLICACIONES Consultando los p-valores obtenidos, se concluye que las variables significativas del modelo creado a un nivel de significación del 5 % son age yyear, ambas variables no dependientes del tiempo. En contraposición, las variables dependientes del tiempo consideradas en el modelo (transplant ysurgery) son no significativas. Sin embargo, en el caso de que se considerase el estudio a un nivel de significación del 10%, la variable surgery sería también significativa en el modelo. Una vez analizados los contrastes de hipótesis individuales, debe estudiarse también lo que ocurre con el contraste de hipótesis conjunto de todas las variables explicativas usadas en el modelo de Cox. Este se enuncia a continuación: " H0:β  p β10, . . . , β40 q 1  0 H1:βj  βj0para algún j  1,...,4, donde al igual que se enunció anteriormente, los índices j  1yj  4hacen referencia a las variables dependientes del tiempo. Dicho contraste de hipótesis se puede resolver a través de tres test diferentes (test de razón de verosimilitud, test de Wald y test score), los cuales son asintóticamente equivalentes. En la tabla siguiente se pueden examinar los estadísticos χ2con sus grados de libertad y p-valores asociados a cada uno de los test expuestos: t(data.frame( "Test de razón de verosimilitud" =summary(mod_cox_dep_t)$logtest, "Test de Wald" =summary(mod_cox_dep_t)$waldtest, "Test score" =summary(mod_cox_dep_t)$sctest )) test df pvalue Test de razón de verosimilitud 15.11148 4 0.004475500 Test de Wald 14.49000 4 0.005876870 Test score 15.03420 4 0.004630811 Tabla 4.24: Contraste de hipótesis conjunto asociado al modelo de Cox creado con covariantes dependientes del tiempo. Como todos los p-valores obtenidos son significativos, se concluye que el modelo creado permite explicar la variable tiempo hasta la muerte del paciente. Para finalizar el estudio de los parámetros y los contrastes de hipótesis asociados a ellos, se adjunta una gráfica que reúne toda la información relevante para el modelo. En dicha gráfica se observarán las cuatro variables predictoras consideradas junto a la exponencial de sus coeficientes estimados. Asimismo, para cada uno de los coeficientes se mostrará un gráfico de la estimación realizada junto a su intervalo de confianza y p-valor asociado. Una explicación más detallada del contenido de la gráfica se puede consultar en la aplicación expuesta en la página 58 para la Figura 4.6. Finalmente, la orden necesaria para obtener la gráfica en cuestión se enuncia a continuación: CAPÍTULO 4. MODELO DE COX EN R 77 4.3. APLICACIONES ggforest(mod_cox_dep_t, main = "Resumen del modelo") surgery year age transplant 1 0 (N=172) (N=172) (N=172) (N=69) (N=103) 0.53 0.86 1.03 0.99 reference (0.26 − 1.09) (0.75 − 0.99) (1.00 − 1.06) (0.54 − 1.83) 0.083 0.038 * 0.048 * 0.974 # Events: 75; Global p−value (Log−Rank): 0.0044755 AIC: 589.13; Concordance Index: 0.64 0.4 0.6 0.8 1 1.2 1.41.61.8 2 Resumen del modelo de Cox Figura 4.12: Resumen gráfico del análisis del modelo de Cox con covariantes dependientes del tiempo realizado con la función ggforest() del paquete survminer. 4.3.2.3. Evaluación del modelo Tal y como se indicó en el Apartado 3.7 de la parte teórica de la evaluación del modelo de Cox con covariantes dependientes del tiempo, la calidad predictiva del modelo se puede estudiar a través de medidas análogas al coeficiente de determinación R2. Estas medidas incluyen el coeficiente sugerido por Nagelkerke (1991), el propuesto por O’Quigley et al. (2005) y el sugerido por Royston (2006). Al igual que se especificó en Apartado 4.3.1.4 de la aplicación anterior del modelo de Cox que incluía solo covariantes independientes del tiempo, dichos coeficientes se pueden calcular y mostrar de la siguiente forma: rsq(mod_cox_dep_t, sigD = NULL) ## $cod ## [1] 0.08410856 ## ## $mer ## [1] 0.1824853 ## ## $mev ## [1] 0.1194867 CAPÍTULO 4. MODELO DE COX EN R 78 4.3. APLICACIONES Dado que los coeficientes obtenidos son muy cercanos al 0, se deduce que la calidad del modelo creado con covariantes dependientes del tiempo no es excesivamente buena para la predicción del tiempo de supervivencia. Sin embargo, debe tenerse en cuenta que el modelo obtenido es válido para determinar aquellas variables que influyen en dicho tiempo y el sentido de dicha influencia. A pesar de que en un modelo de Cox con covariantes dependientes del tiempo no hace falta estudiar si las variables involucradas en el modelo cumplen la hipótesis de riesgos proporcionales o no, a continuación, se presenta un pequeño estudio para analizar qué es lo que ocurre con dicha hipótesis en el modelo bajo estudio. Al igual que se mostró en la diagnosis realizada sobre el modelo de Cox con covariantes independientes del tiempo del apartado anterior, el estudio de si la hipótesis de riesgos proporcionales se viola o no, se puede hacer de manera gráfica y analítica mediante el test de bondad de ajuste (GOF). En la siguiente tabla, se pueden analizar los resultados obtenidos tras resolver los diferentes contrastes de bondad de ajuste asociados a cada una de las variables predictoras del modelo de Cox creado: cox.zph(mod_cox_dep_t, transform = "rank") Variable chisq df p transplant 0.25244 1 0.62 age 0.83072 1 0.36 year 1.05187 1 0.31 surgery 0.00385 1 0.95 GLOBAL 2.68433 4 0.61 Tabla 4.25: Estudio de la hipótesis de riesgos proporcionales de manera individual y global para el modelo de Cox creado con covariantes dependientes del tiempo. A la vista de los p-valores obtenidos para cada contraste realizado, se concluye que no existen evidencias significativas en contra de la hipótesis nula de dichos contrastes. Como consecuencia de ello, todas las variables predictoras del modelo creado verifican la hipótesis de riesgos proporcionales a nivel individual. Asimismo, dicha hipótesis se cumple de manera conjunta para todas las variables involucradas en el modelo. Nótese que en un modelo de Cox con covariantes dependientes del tiempo, el hecho de que se cumpla la hipótesis de riesgos proporcionales para cada una de las variables explicativas no es usual. Sin embargo, en el caso bajo estudio, esta condición se cumple ya que, tal y como se ha comprobado en el Apartado 4.3.2.2, las variables dependientes del tiempo son no significativas y pueden considerarse nulos sus coeficientes estimados. Para el estudio gráfico de las hipótesis de riesgos proporcionales, es necesario calcular los residuos de Schoenfeld asociados a cada uno de los individuos del conjunto de datos bajo estudio y a cada variable explicativa considerada en el modelo. A continuación, se muestra el cálculo de estos residuos escalados: CAPÍTULO 4. MODELO DE COX EN R 79 4.3. APLICACIONES cox.zph(mod_cox_dep_t, transform = "rank")$y time transplant age year surgery 1 -0.8012696 0.09340004 -1.1193608 9.347928 2 -0.3913630 0.11085396 0.7490845 -2.239041 2 -0.4647442 0.10960172 0.4462945 -1.915566 2 -0.1436648 -0.03385701 -0.4402613 8.687208 3 -0.3362103 0.04168269 0.2541290 -1.707941 3 0.4353129 -0.12542568 -0.1590904 -1.309205 . . . . . . . . . . . . . . . 733 -0.5020444 0.204263442 -0.19367658 -2.618850 852 0.4910750 0.035798724 -0.13991515 -3.136782 980 0.3940582 0.051646341 0.09952175 6.115883 996 0.2148273 0.090662022 0.09884397 6.524236 1032 0.7810621 0.009864047 -0.25340874 -2.322435 1387 0.2522059 0.149259359 -0.30396670 -1.402432 Tabla 4.26: Residuos escalados de Schoenfeld asociados al modelo de Cox creado con covariantes dependientes del tiempo. Una vez obtenidos los residuos de Schoenfeld, se procede a realizar unas gráficas de estos para su posterior análisis: par(mfrow = c(1,2)) plot(cox.zph(mod_cox_dep_t, transform = "rank"), col = "red") par(mfrow = c(1,1)) CAPÍTULO 4. MODELO DE COX EN R 80 4.3. APLICACIONES Figura 4.13: Gráficas de los residuos escalados de Schoenfeld con la orden plot(cox.zph()) para el modelo de Cox creado con covariantes dependientes del tiempo. Observando las gráficas obtenidas, no se tiene muy claro si las variables usadas en el ajuste del modelo cumplen o no la hipótesis de riesgos proporcionales pues no se ve claramente si existen patrones en los residuos. Ante esta dificultad, se presenta otra forma de obtener las gráficas de los residuos escalados de Schoenfeld mucho más interpretable. Esta alternativa hace uso de la función ggcoxzph() y se aplica de la siguiente forma: ggcoxzph(cox.zph(mod_cox_dep_t), font.main = 10,font.submain = 8,font.x = 10, font.y = 10,font.tickslab = 10,font.legend = 8 ) −20 −10 0 10 20 3.8 17 37 62 83 190340 990 Time Beta(t) for transplant Schoenfeld Individual Test p: 0.6524 −1.0 −0.5 0.0 0.5 1.0 3.8 17 37 62 83 190340 990 Time Beta(t) for age Schoenfeld Individual Test p: 0.3518 −3 0 3 3.8 17 37 62 83 190340 990 Time Beta(t) for year Schoenfeld Individual Test p: 0.1902 −20 −10 0 10 20 30 3.8 17 37 62 83 190340 990 Time Beta(t) for surgery Schoenfeld Individual Test p: 0.755 Global Schoenfeld Test p: 0.482 Fuente propia Figura 4.14: Estudio gráfico y analítico de los residuos de Schoenfeld con la orden ggcoxzph() de la librería survminer. CAPÍTULO 4. MODELO DE COX EN R 81 4.3. APLICACIONES Otra forma equivalente de obtener dicha estimación es mediante la siguiente orden: survfit(mod_cox_dep_t, stype = 1) Tras el cálculo de su estimación, se presenta un resumen en el que se puede consultar para cada tiempo observado el número de individuos en riesgo (n.risk), el número de individuos que presentan el evento en dicho instante (n.event), la estimación de S p t q (survival) junto con su intervalo de confianza al 95 % (lower 95 % CI yupper 95 % CI ) y la desviación estándar de la estimación realizada (std.err): summary(func_superv_dep) time n.risk n.event survival std.err lower 95 % CI upper 95 % CI 1 103 1 0.990 0.00956 0.9718 1.000 2 102 3 0.962 0.01895 0.9253 1.000 3 99 3 0.933 0.02479 0.8856 0.983 5 96 2 0.914 0.02797 0.8605 0.970 6 94 2 0.894 0.03074 0.8361 0.957 8 92 1 0.885 0.03205 0.8239 0.950 . . . . . . . . . . . . . . . . . . . . . 733 16 1 0.236 0.08937 0.1126 0.496 852 14 1 0.217 0.08878 0.0974 0.484 980 11 1 0.196 0.08791 0.0812 0.472 996 10 1 0.176 0.08608 0.0671 0.459 1032 9 1 0.155 0.08296 0.0545 0.442 1387 6 1 0.130 0.07748 0.0405 0.418 Tabla 4.30: Estimación de la función de supervivencia a partir del modelo de Cox creado con covariantes dependientes del tiempo. Por último, para completar el estudio de la función de supervivencia, se hace una gráfica de su estimación para todo el conjunto de datos involucrado en el modelo de Cox creado con covariantes dependientes del tiempo: ggsurvplot(survfit(mod_cox_dep_t), data = uissurv, conf.int = T, xlab = "t (días)",ylab = "S(t) estimada", legend.title = "Estimación", legend.labs = "Kaplan-Meier", ggtheme = theme_grey() ) CAPÍTULO 4. MODELO DE COX EN R 88 4.3. APLICACIONES + + +++ + + + ++ ++ + + +++ +++++ + + ++ + ++ + + + + + ++ + + + + + + ++ + +++++++++++++++ +++ + + ++ ++ + 0.00 0.25 0.50 0.75 1.00 0 500 1000 1500 2000 t (días) S(t) estimada Estimación +Kaplan−Meier Fuente propia Figura 4.17: Estimación de la función de supervivencia para el modelo creado y ajustado con covariantes dependientes del tiempo. CAPÍTULO 4. MODELO DE COX EN R 89 Conclusión El análisis de supervivencia es una técnica muy útil cuando el objetivo del problema o análisis a realizar es el estudio y modelización de la variable tiempo que transcurre hasta que se produce un determinado suceso de interés. Se trata también de una rama de la estadística en continua evolución y estudio, a la que todavía le queda mucho por avanzar tal y como se puede comprobar en los artículos recientes que se han utilizado como bibliografía a la hora de la realización de este trabajo. Asimismo, el análisis de supervivencia destaca por las aplicaciones que tiene hoy en día sobre todo en disciplinas como la medicina o la ingeniería. A lo largo del trabajo, se ha realizado un estudio detallado sobre el modelo de Cox con covariantes dependientes e independientes del tiempo. Sin embargo, dada la magnitud de la materia, no se ha podido profundizar en algunos aspectos citados en el presente Trabajo Fin de Máster. A continuación, se enumeran diferentes conceptos y temas que en un posterior estudio se podrían profundizar: Propiedades asintóticas tanto del modelo de Cox con covariantes dependientes del tiempo como dicho modelo con covariantes independientes del tiempo. Estudio de aplicaciones del modelo de Cox con covariantes dependientes del tiempo en el ámbito científico, médico.. . Estimación de parámetros para el modelo de Cox con covariantes dependientes del tiempo y tipos de censuras diferentes a la censura a la derecha. Nótese que para la censura a la izquierda bastaría con hacer una pequeña transformación. Tal y como determinan Goel y Klein (1992), la obtención de datos censurados a la derecha partiendo de datos censurados a la izquierda se consigue con el producto de cada uno de los datos por  1. De esta forma, cualquier técnica o estudio realizado para los datos censurados a la derecha se pueden extrapolar a los datos censurados a la izquierda. Estimación de parámetros del modelo de Cox con covariantes dependientes del tiempo cuando existen tiempos de fallos coincidentes o empates. A la hora de resolver la estimación de los parámetros de cada uno de los modelos de Cox expuestos a lo largo del trabajo, se necesita resolver un sistema de ecuaciones no lineales cuya expresión no se detalla en el presente trabajo. Se podría ampliar con el análisis de los procedimientos iterados de aproximación a la solución de los problemas de optimización 90 Bibliografía O. Aalen, O. Borgan y H. Gjessing. Survival and event history analysis: a process point of view. Springer Science & Business Media, 2008. M.P. Allen. Understanding regression analysis. Springer Science & Business Media, 1997. doi: https://doi.org/dnkkn9. A. Allignol y A. Latouche. CRAN Task View: Survival Analysis. R Foundation for Statistical Computing, Version 2022-03-07. Disponible en: https://cran.r-project.org/ web/views/Survival.html. E.E. Alvarez y J. Ferrario. Revisión de estimación robusta en modelos semiparamétricos de supervivencia. Estadistica, 64, 2012. C.A. Bellera, G. MacGrogan, M. Debled, C.T. De Lara, V. Brouste y S. Mathoulin- Pélissier. Variables with time-varying effects and the Cox model: some statistical concepts illustrated with a prognostic factor study in breast cancer. BMC medical research methodology, 10(1):1–12, 2010. R.A. Berk, P.H. Rossi y K.J. Lenihan. Money, Work, and Crime: Experimental Evidence. New York: Academic Press, 1980. P. Biecek, S. Fabian, A. Kassambara y M. Kosinski. survminer: Drawing Survival Curves using ’ggplot2’. R package version 0.4.9, 2021. Disponible en: https://CRAN.R-project. org/package=survminer. N. Breslow. Covariance analysis of censored survival data. Biometrics, 30:89–99, 1974. R.H. Chan y J. Ma. A multiplicative iterative algorithm for box-constrained penalized likelihood image restoration. IEEE transactions on image processing, 21(7):3168–3181, 2012. D. Cox y D. Oakes. Analysis of survival data, volume 21. CRC press, 1984. D.R. Cox. Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological), 34(2):187–220, 1972. doi: https://doi.org/gf2fbn. D.R. Cox. Partial likelihood. Biometrika, 62(2):269–276, 1975. J. Crowley y M. Hu. Covariance analysis of heart transplant survival data. Journal of the American Statistical Association, 72(357):27–36, 1977. C. Dardis. survMisc: Miscellaneous Functions for Survival Data. R package version 0.5.6, 2022. Disponible en: https://CRAN.R-project.org/package=survMisc. 91 Bibliografía E. de las Heras. Comparación de curvas de supervivencia con datos censurados. Trabajo Fin de Grado, Universidad de Sevilla. Disponible en https://hdl.handle.net/11441/ 134450, 2021. P. Deuflhard. Newton methods for nonlinear problems: affine invariance and adaptive algorithms, volume 35. Springer Science & Business Media, 2011. doi: https://doi.org/ dv8qcx. B. Efron. The efficiency of Cox’s likelihood function for censored data. Journal of the American statistical Association, 72(359):557–565, 1977. L. Fahrmeir, G. Tutz, W. Hennevogl y E. Salem. Multivariate statistical modelling based on generalized linear models, volume 425. Springer, 1994. L.D. Fisher y D.Y. Lin. Time-dependent covariates in the Cox proportional-hazards regression model. Annual review of public health, 20(1):145–157, 1999. J. Fox y S. Weisberg. An R companion to applied regression. Sage publications, 2018. J. Fox y S. Weisberg. Cox proportional-hazards regression for survival data. An R companion to applied regression, 2018, 2023. J. Friedman, T. Hastie, R. Tibshirani, B. Narasimhan, K. Tay, N. Simon, J. Qian y J. Yang. glmnet: Lasso and Elastic-Net Regularized Generalized Linear Models. R package version 4.1.8., 2023. Disponible en: https://CRAN.R-project.org/package= glmnet. P. Goel y J.P. Klein. Survival Analysis: State of the Art. Springer Science & Business Media, First edition, 1992. A. González. Selección de variables: Una revisión de métodos existentes. Trabajo Fin de Máster, Universidad de los Andes. Disponible en http://eio.usc.es/pub/mte/index. php/es/trabajos-fin-de-master/finalizados, 2015. I. Guyon y A. Elisseeff. An introduction to variable and feature selection. Journal of machine learning research, 3(Mar):1157–1182, 2003. G. Gómez y C. Cadarso-Suárez. El modelo de riesgos proporcionales de Cox y sus extensiones. Impacto en Estadística y Biomedicina. La Gaceta de la RSME, 20(3):513–538, 2017. M.A. Gómez. Karl Pearson, el creador de la estadística matemática. Historia de la Probabilidad y la estadística (IV), páginas 351–356, 2009. F.E. Harrell. Regression modeling strategies: with applications to linear models, logistic regression, and survival analysis, volume 608. Springer, 2001. F.E. Harrell y K.L. Lee. Verifying assumptions of the Cox proportional hazards model. In Proceedings of the eleventh annual SAS Users group international conference, páginas 823–828. SAS Institute Inc, Cary, NC, 1986. F.E. Harrell, K.L. Lee y D.B. Mark. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in medicine, 15(4):361–387, 1996. 92 Bibliografía E. Harrison, T. Drake y R. Ots. finalfit: Quickly Create Elegant Regression Results Tables and Plots when Modelling. R package version 1.0.6, 2023. Disponible en: https:// CRAN.R-project.org/package=finalfit. P. Hougaard. Analysis of multivariate survival data, volume 564. Springer, 2000. K. Ito y K. Kunisch. A variational approach to sparsity optimization based on Lagrange multiplier theory. Inverse problems, 30(1):015001, 2013. C.M. Jareño. Selección de subconjuntos de variables en modelos estadísticos. Trabajo Fin de Grado, Universidad de Sevilla. Disponible en https://hdl.handle.net/11441/142918, 2022. J.D. Kalbfleisch y R.L. Prentice. The statistical analysis of failure time data. John Wiley & Sons, 2011. W. Karush. Minima of functions of several variables with inequalities as side constraints. M. Sc. Dissertation. Dept. of Mathematics, Univ. of Chicago, 1939. D.G. Kleinbaum y M. Klein. Survival Analysis: A Self-Learning Text. Springer, Third edition, 2012. doi: https://doi.org/cmnq7s. E.L. Korn y R. Simon. Measures of explained variation for survival data. Statistics in medicine, 9(5):487–503, 1990. H. Kuhn y A. Tucker. Nonlinear programming. In Proceedings of 2nd Berkeley symposium, páginas 481–492. Berkeley: University of California Press, 1951. S. Lemeshow, S. May y D.W. Hosmer Jr. Applied survival analysis: regression modeling of time-to-event data. John Wiley & Sons, 2008. doi: https://doi.org/ftds85. J. Li y S. Ma. Survival analysis in medicine and genetics. CRC Press, 2013. X. Liu. Survival analysis: models and applications. John Wiley & Sons, 2012. doi: https://doi.org/kcd7. J.S. López. Selección de modelos de supervivencia para pacientes con cáncer gastrointestinal-un enfoque bajo modelos de fragilidad. Trabajo Fin de Máster, Universidad de los Andes. Disponible en http://hdl.handle.net/1992/12348, 2014. D.G. Luenberger y Y. Ye. Linear and nonlinear programming, volume 2. Springer, 1984. J. Ma. Positively constrained multiplicative iterative algorithm for maximum penalized likelihood tomographic reconstruction. IEEE Transactions on Nuclear Science, 57(1): 181–192, 2010. D.F. Moore. Applied survival analysis using R, volume 473. Springer, 2016. N.J.D. Nagelkerke. A note on a general definition of the coefficient of determination. biometrika, 78(3):691–692, 1991. J. O’Quigley, R. Xu y J. Stare. Explained randomness in proportional hazards models. Statistics in medicine, 24(3):479–489, 2005. 93 Bibliografía A.L. Palmer. Modelo de regresión de Cox: ejemplo numérico del proceso de estimación de parámetros. Psicothema, páginas 387–402, 1993. A. Pettitt y I.B Daud. Investigating time dependence in Cox’s proportional hazards model. Journal of the Royal Statistical Society Series C: Applied Statistics, 39(3):313– 329, 1990. M.D. Pérez, J. Basulto y J.A. Camúñez. Un antecedente histórico de regresión lineal: la estimación mediana propuesta por Boscovich. Gaceta de la Real Sociedad Matematica Española, 22(2):351–364, 2019. 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, E. Regino y S.Y. Guerrero. Comparación de métodos de estimación en regresión de Cox. Comunicaciones en Estadística, 10(1):101–112, 2017. James O Ramsay. Monotone regression splines in action. Statistical science, páginas 425–441, 1988. P. Royston. Explained variation for survival models. The Stata Journal, 6(1):83–96, 2006. RStudio Team. RStudio: Integrated Development Environment for R. RStudio, PBC., Boston, MA, 2020. Disponible en: http://www.rstudio.com/. G. Ruiz et al. Los orígenes del método de mínimos cuadrados. Suma, 2003. N.M. Russo. Aplicaciones del análisis de supervivencia en la biotecnología. Trabajo Fin de Grado, Universidad de Almería. Disponible en http://hdl.handle.net/10835/13488, 2021. B. San José, E. Pérez y R. Madero. Métodos estadísticos en estudios de supervivencia. Anales de Pediatría Continuada, 7(1):55–59, 2009. J.P. Sánchez. Modelos semiparamétricos en el análisis de supervivencia aplicado a la mejora genética. Dep. de Ciencia Animal, UPV, 2005. M. Schemper. The explained variation in proportional hazards regression. Biometrika, 77 (1):216–218, 1990. M. Schemper y J. Stare. Explained variation in survival analysis. Statistics in medicine, 15(19):1999–2012, 1996. D. Schoenfeld. Partial residuals for the proportional hazards model. Biometrika, 69:51–55, 1982. M. Spreafico, F. Ieva y M. Fiocco. Modelling time-varying covariates effect on survival via functional data analysis: application to the MRC BO06 trial in osteosarcoma. Statistical Methods & Applications, 32(1):271–298, 2023. J. M. Stanton. Galton, Pearson, and the peas: A brief history of linear regression for statistics instructors. Journal of Statistics Education, 9(3), 2001. doi: https://doi.org/ gd82dx. 94 Bibliografía M. Tableman y J.S. Kim. Survival analysis using S: analysis of time-to-event data. CRC press, 2003. K. Tay, N. Simon, J. Friedman, T. Hastie, R. Tibshirani y B. Narasimhan. Regularized Cox Regression, 2023. M. Thackham y J. Ma. On maximum likelihood estimation of the semi-parametric Cox model with time-varying covariates. Journal of Applied Statistics, 47(9):1511–1528, 2020. T. Therneau, C. Crowson y E. Atkinson. Using time dependent covariates and time dependent coefficients in the Cox model. Survival Vignettes, 2(3):1–31, 2022. T.M. Therneau. A Package for Survival Analysis in R. R package version 3.5-7, 2023. Disponible en: https://CRAN.R-project.org/package=survival. T.M. Therneau y P.M. Grambsch. Modeling Survival Data: Extending the Cox Model. Springer, 2000. doi: https://doi.org/j8kf. T.M. Therneau, P.M. Grambsch y T.R. Fleming. Martingale-based residuals for survival models. Biometrika, 77(1):147–160, 1990. A.A. Tsiatis. Semiparametric Theory and Missing Data. Springer, 2006. doi: https: //doi.org/dx8pkp. P. Velasco. Modelo de regresión de Cox y sus aplicaciones biosanitarias. Trabajo Fin de Grado, Universidad de Sevilla. Disponible en http://hdl.handle.net/11441/43493, 2016. A. Webb y J. Ma. Cox models with time-varying covariates and partly-interval censoring– A maximum penalised likelihood approach. Statistics in medicine, 42(6):815–833, 2023. 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.4.3, 2023a. Disponible en: https://CRAN. R-project.org/package=ggplot2. H. Wickham, R. François, L. Henry, K. Müller y D. Vaughan. dplyr: A Grammar of Data Manipulation. R package version 4.1.8., 2023b. Disponible en: https://CRAN. R-project.org/package=glmnet. W.H. Wong y T.A. Severini. On maximum likelihood estimation in infinite dimensional parameter spaces. The Annals of Statistics, páginas 603–632, 1991. J. Xu, J. Ma, M.H. Connors y H. Brodaty. Proportional hazard model estimation under dependent censoring using copulas and penalized likelihood. Statistics in medicine, 37 (14):2238–2251, 2018. Z. Zhang, J. Reinikainen, K.A. Adeleke, M.E. Pieterse y C.G.M. Groothuis-Oudshoorn. Time-varying covariates and coefficients in Cox regression models. Annals of translational medicine, 6(7), 2018. R. Álvarez. Sir Francis Galton, padre de la eugenesia, volume 4. Editorial CSIC-CSIC Press, 1985. 95