Análisis de supervivencia a través de árboles de decisión

Trabajo Fin de Grado Doble Grado en Matem´ aticas y Estad ´ ıstica Facultad de Matem´aticas An´alisis de supervivencia a trav´es de ´arboles de decisi´on Estudiante: Jos´e Antonio Tirado Carballo Tutor: Juan M. Mu˜noz Pichardo Departamento: Estad´ıstica e Investigaci´on Operativa Fecha: 5 de junio de 2024 A mi madre y a mi padre, por siempre dejarme escoger mi camino y ense˜narme que el talento no es nada sin esfuerzo. A mi familia y amigos, por estar en los momentos dif´ıciles. Espero que todos est´en orgullosos de lo que estoy construyendo. Resumen El an´alisis de datos es una herramienta esencial en diversas disciplinas, permitiendo extraer informaci´on valiosa y realizar predicciones basadas en datos hist´oricos. Entre las m´ultiples t´ecnicas disponibles, los ´arboles de decisi´on se destacan por su simplicidad, interpretabilidad y versatilidad. Este trabajo de fin de grado se centra en los ´arboles de supervivencia, abordando su base en los ´arboles de decisi´on. En el cap´ıtulo 1. Introducci´on, se realiza una breve introducci´on a los ´arboles de decisi´on y su aplicabilidad en el estudio del an´alisis de supervivencia. En el cap´ıtulo 2. ´ Arboles de decisi´on, comenzamos con una introducci´on general que da paso a la discusi´on sobre ´arboles de regresi´on y ´arboles de clasificaci´on, detallando los procesos de construcci´on y poda de estos modelos. Adem´as, exploramos las ventajas y desventajas inherentes a los ´arboles de decisi´on, ofreciendo una visi´on equilibrada sobre cu´ando y c´omo utilizarlos de manera efectiva. Esta secci´on tambi´en incluye una introducci´on al an´alisis con datos censurados, preparando el terreno para la discusi´on m´as especializada en ´arboles de supervivencia. En el cap´ıtulo 3. Bosques aleatorios de supervivencia, se ampl´ıa el horizonte de los m´etodos basados en ´arboles introduciendo los m´etodos de ensemble, como bagging yrandom forests. Aqu´ı, se analiza c´omo la combinaci´on de m´ultiples ´arboles puede llevar a una mejora significativa en la precisi´on de las predicciones. Adem´as, se presentan conceptos cruciales como el error out-of-bag y la importancia de los predictores, con subsecciones espec´ıficas que abordan la funci´on de riesgo acumulada y el error de predicci´on en bosques de supervivencia. A pesar de que los ´arboles de supervivencia constituyen el n´ucleo de este trabajo, considero oportuno incluir una discusi´on sobre los bosques. Los bosques, al combinar m´ultiples ´arboles y agregar aleatoriedad en su construcci´on, no solo mejoran la robustez y la precisi´on de los modelos, sino que tambi´en permiten manejar mejor la variabilidad en los datos y las interacciones complejas entre variables. Esta secci´on es fundamental para entender c´omo se pueden superar algunas limitaciones de los ´arboles de decisi´on individuales. Finalmente, en el cap´ıtulo 4. Ilustraciones con R Program, proporcionamos una serie de ejemplos pr´acticos que ilustran la implementaci´on de ´arboles y bosques de supervivencia en R, aplicando los conceptos discutidos. Abstract Data analysis is an essential tool in various disciplines, enabling the extraction of valuable information and making predictions based on historical data. Among the many available techniques, decision trees stand out for their simplicity, interpretability, and versatility. This thesis focuses on survival trees, addressing their foundation in decision trees. In chapter 1. Introduction, a brief introduction to decision trees and their applicability in the study of survival analysis is provided. In chapter 2. Decision trees, we begin with a general introduction that leads to a discussion on regression trees and classification trees, detailing the processes of constructing and pruning these models. Additionally, we explore the inherent advantages and disadvantages of decision trees, providing a balanced view of when and how to use them effectively. This section also includes an introduction to censored data analysis, setting the stage for the more specialized discussion on survival trees. In chapter 3. Random Survival Forests, we broaden the scope of treebased methods by introducing ensemble methods such as bagging and random forests. Here, we analyze how the combination of multiple trees can lead to a significant improvement in prediction accuracy. Additionally, crucial concepts such as out-of-bag error and predictor importance are presented, with specific subsections addressing the cumulative hazard function and prediction error in survival forests. Although survival trees constitute the core of this work, I find it appropriate to include a discussion on forests. By combining multiple trees and introducing randomness in their construction, forests not only enhance the robustness and accuracy of the models but also allow better handling of data variability and complex interactions between variables. This section is fundamental to understanding how some limitations of individual decision trees can be overcome. Finally, in chapter 4. Illustrations with R Program, we provide a series of practical examples illustrating the implementation of survival trees and forests in R, applying the discussed concepts. ´ Indice general 1. Introducci´on 1 1.1. Unpocodehistoria........................ 2 2. ´ Arboles de decisi´on 5 2.1. Introducci´on............................ 5 2.1.1. ´ Arboles de regresi´on . . . . . . . . . . . . . . . . . . . 6 2.1.1.1. Proceso de construcci´on de un ´arbol de regre- si´on ....................... 7 2.1.1.2. Algoritmo b´asico de poda de un ´arbol de decisi´on ...................... 12 2.1.2. ´ Arboles de clasificaci´on . . . . . . . . . . . . . . . . . . 15 2.1.2.1. Proceso de construcci´on de un ´arbol de clasificaci´on ..................... 16 2.1.3. Ventajas y desventajas de los ´arboles . . . . . . . . . . 21 2.2. Introducci´on al an´alisis con datos censurados . . . . . . . . . . 22 2.3. ´ Arboles de supervivencia . . . . . . . . . . . . . . . . . . . . . 24 2.3.1. Definiciones previas . . . . . . . . . . . . . . . . . . . . 25 2.3.2. Un conjunto de cuestiones “binarias” para hacer las divisiones ......................... 27 2.3.3. Un criterio para evaluar la bondad de las divisiones . . 27 2.3.3.1. Estad´ıstico log-rank . . . . . . . . . . . . . . 27 2.3.3.2. Reducci´on de impureza . . . . . . . . . . . . . 29 2.3.4. Una regla de parada . . . . . . . . . . . . . . . . . . . 31 2.3.5. Una regla para predecir en cada nodo terminal . . . . . 32 2.3.5.1. Estimador de Kaplan-Meier . . . . . . . . . . 32 2.3.5.2. Estimador de Nelson-Aalen . . . . . . . . . . 33 2.4. ElmodelodeCox......................... 35 2.4.1. Formulaci´on del modelo . . . . . . . . . . . . . . . . . 35 3 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON las m´as cl´asicas el enfoque de los ´arboles de clasificaci´on y regresi´on, ideado principalmente por Breiman et al. (1984) 2.1.1. ´ Arboles de regresi´on Un ´arbol de regresi´on consiste en una serie de reglas de divisi´on del espacio definido por las variables predictoras. ´ Estos nos ofrecen una alternativa al an´alisis tradicional de discriminaci´on o a la predicci´on tradicional de regresi´on. Se trata de una t´ecnica especialmente relevante en el an´alisis de datos de supervivencia, es decir, datos relativos al tiempo hasta ocurrir un determinado suceso. Supongamos que queremos predecir una variable Yempleando dos variables predictoras cuantitativas X1yX2. La variable X1puede tomar cualquier valor positivo, mientras que X2se encuentra en el intervalo [0,100]. Tenemos entonces que el espacio predictor χest´a definido como R+×[0,100]. La divi- si´on de este espacio se lleva a cabo a trav´es de un ´arbol de decisi´on, similar al siguiente: Figura 2.1: Ejemplo de ´arbol. La predicci´on de la variable Yen cada regi´on Rjse realiza utilizando ¯y(Rj), que es la media muestral de Ypara los puntos observados en Rj. En estas t´ecnicas, las regiones Rjse corresponden con los nodos terminales o hojas del ´arbol. Los puntos del ´arbol donde se efect´ua la divisi´on del espacio predictor se denominan nodos, y las conexiones entre estos nodos se conocen como ramas del ´arbol. Como se puede observar, los ´arboles de regresi´on presentan una simplificaci´on en la relaci´on entre variables, lo que permite una f´acil interpretabilidad. De esta forma, un ´arbol de regresi´on con Jnodos se puede expresar como: 6 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON f(x) = J X j=1 cj·IRj(x) donde IRj(x) = (1 si x∈Rj, 0 c.c. La f´ormula presenta la funci´on de predicci´on fde un ´arbol de regresi´on que tiene Jnodos. Cada nodo o regi´on se representa por Rj, y tiene asociado un valor constante cj, cuyo valor depender´a del enfoque de optimalidad que empleemos. Gracias a la funci´on indicadora, tenemos que, dado un valor x, f(x) ser´a igual al valor cjde la regi´on Rjen la que xcae. Si lo que deseamos es minimizar la suma de cuadrados sobre la muestra o conjunto de entrenamiento, tenemos que es inmediato m´ın c1...cJ n X i=1 [yi−f(xi)]2⇒ˆcj= ¯y(Rj) = media muestral de Yen Rj. El error sobre el conjunto de entrenamiento (error de entrenamiento o p´erdida de entrenamiento) es una estimaci´on del error esperado: E[Loss(Y, f(X)], donde Loss() es la funci´on de error seleccionada. En el caso de m´ınimos cuadrados tenemos que E[Loss(Y, f(X) es el error cuadr´atico esperado. 2.1.1.1. Proceso de construcci´on de un ´arbol de regresi´on Esencialmente el proceso consta de dos pasos: Paso 1: Se construye una partici´on del espacio predictor en Jregiones Rj (j= 1, . . . , J) Paso 2: Asignamos a cada observaci´on que cae en la regi´on Rjla predicci´on dada por la media de los valores de la respuesta para las observaciones de la muestra de Rj El problema radica en desarrollar el Paso 1. Por simplicidad, se divide el espacio predictor en rect´angulos o cajas, y se busca el conjunto de cajas que minimice la suma total de los cuadrados de los residuos RSS, que se define como RSS = J X j=1 X i∈Rj (yi−ˆyRj)2, 7 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON donde ˆyRjcorresponde a la media muestral de las respuestas observadas en la muestra de entrenamiento para la j-´esima caja. Debido a la inviabilidad computacional de evaluar todas las particiones posibles, se adopta un m´etodo de divisi´on binaria recursiva de arriba hacia abajo (top-down). As´ı, una vez establecida una partici´on, esta permanece fija para las divisiones siguientes. Inicialmente se considera todo el espacio predictor en una ´unica regi´on llamada nodo ra´ız y se procede a dividir en nodos hijos. Cada nodo se divide en dos nuevos nodos hijos. Cada divisi´on se realiza de manera ´optima, exclusivamente a trav´es de las observaciones incluidas en cada subgrupo y sin considerar el futuro. La divisi´on planteada anteriormente es v´alida para predictores Xknum´ericos. En caso de predictores cualitativos se procede con cuestiones binarias, como se indicar´a posteriormente. En cada nodo, para llevar a cabo la divisi´on binaria recursiva, tomamos un predictor Xj, y un punto de corte sde manera que la divisi´on del espacio predictor en las regiones {X:Xk< s}y{X:Xk≥s}conduzca a la mayor reducci´on posible de RSS. Este procedimiento se realiza considerando todas las variables predictoras X1, . . . , Xpy todos los valores posibles de s. As´ı, para cada kysse define el par de semiplanos: R1(k, s) = {X:Xk< s}y R2(k, s) = {X:Xk≥s}y se seleccionan: k, s : minimizan X i∈R1(k,s) (yi−ˆyR1(k,s))2+X i∈R2(k,s) (yi−ˆyR2(k,s))2. Una vez encontrados estos valores de kys, se repite el proceso, buscando el mejor predictor y el mejor punto de corte con el fin de dividir los datos adicionales a fin de minimizar el RSS dentro de cada una de las regiones resultantes. Eso s´ı, en las divisiones que siguen se divide una de las regiones previamente identificadas, en lugar de todo el espacio predictor. Tengamos en cuenta que, previo a la primera divisi´on, la suma residual de cuadrados o estad´ıstico desviaci´on de m´etodo de predicci´on basado en la media es D= n X i=1 (yi−¯y)2. 8 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Tras la divisi´on mediante la variable predictora Xky el punto de corte s se tiene que D= n X i=1 (yi−¯y)2=SCW+SCB, SCW=X i∈R1(k,s) (yi−ˆyR1(k,s))2+X i∈R2(k,s) (yi−ˆyR2(k,s))2, SCB=n1(ˆyR1(k,s)−¯y)2+n2(ˆyR2(k,s)−¯y)2. Estas ecuaciones definen la descomposici´on de la suma total de cuadrados Den la suma de cuadrados dentro de los grupos SCWy la suma de cuadrados entre los grupos SCB, respectivamente. Las n1yn2representan el n´umero de observaciones en las regiones R1(k, s) y R2(k, s), mientras que ˆyR1(k,s)y ˆyR2(k,s)son las medias de las respuestas en esas regiones. De esta forma, m´ın k,s SCW(k, s)⇐⇒ m´ax k,s SCB(k, s). Esto concuerda con el prop´osito de que, cuando se construye un ´arbol de regresi´on, buscamos dividir el conjunto de datos de tal manera que los datos dentro de cada grupo (o “hoja” del ´arbol) sean lo m´as similares posible (minimizando SCW), y que haya una clara distinci´on entre los grupos (maximizando SCB) . Por tanto, SCBmide cu´anta variabilidad hay entre los diferentes grupos creados por las divisiones en el ´arbol. Si SCBes grande en comparaci´on con D, significa que la divisi´on que hemos hecho captura una gran parte de la variabilidad de los datos y por lo tanto es probable que sea una divisi´on ´util. En otras palabras, queremos encontrar las divisiones que maximicen la diferencia en las medias de las respuestas entre los grupos, lo cual es equivalente a maximizar SCB. Tras el primer proceso de divisi´on, se tendr´an dos regiones. De nuevo, se vuelve a seleccionar kyspara dividir una de estas dos regiones de forma que se minimice RSS. El proceso contin´ua hasta que se alcanza un criterio de parada. Llegados a este punto, debemos mencionar que el proceso descrito, aunque puede ser buen predictor sobre la muestra de entrenamiento, es posible que sobreajuste los datos, lo que se traduce en un rendimiento pobre por la posible complejidad elevada del ´arbol resultante. 9 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Un ´arbol m´as peque˜no con menor n´umero de divisiones podr´ıa conducir a menor varianza y mejor interpretaci´on a costa de un leve aumento de sesgo. De manera general, una divisi´on aleatoria de los datos en distintos subconjuntos proporcionan distintos ´arboles de regresi´on cuyas primeras particiones tienden a ser similares pero difieren sustancialmente m´as cerca de los nodos terminales. Estos nodos m´as profundos tienden a sobreajustarse a atributos espec´ıficos de los datos de la muestra; en consecuencia, muestras ligeramente diferentes dar´an como resultado predicciones variables en los nodos terminales. Introducimos aqu´ı el t´ermino de poda, que consiste en podar estos nodos de niveles inferiores con el objetivo de estabilizar las predicciones y generalizar mejor para datos nuevos, a costa de introducir sesgo. Seg´un lo descrito anteriormente, la construcci´on de un ´arbol de regresi´on requiere de los siguientes elementos: Un conjunto de cuestiones “binarias” para hacer las divisiones. Un criterio para evaluar la bondad de las divisiones. Una regla de parada. Una regla para predecir en cada nodo terminal. A continuaci´on se describen brevemente estos elementos. Un conjunto de cuestiones “binarias” para hacer las divisiones Cada divisi´on del ´arbol se basa en una ´unica variable predictora. En el caso de las variables binarias (0/1, Si/No) solo existe una divisi´on posible. Sin embargo, para variables cuantitativas u ordinales, existen te´oricamente infinitas posibles divisiones, como por ejemplo, del tipo {x≤c}. Dado un conjunto de nobservaciones de la muestra, si se ordenan de menor a mayor como x(1) ≤x(2) ≤. . . ≤x(n), las divisiones posibles se reducen a las que est´an entre estas observaciones. Por lo tanto, las posibles divisiones ser´ıan {x≤cj}donde cj=1 2[x(j)+x(j+1)] para cada par de observaciones consecutivas. En el caso de las variables cualitativas nominales que tienen Bmodalidades, las posibles divisiones pueden expresarse como variaciones con repetici´on, 10 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON lo que lleva a un n´umero total de divisiones posibles de VB R 2−2 donde VB R representa las variaciones con repetici´on de Belementos. Esto se simplifica a2B−2 2= 2(B−1) −1, proporcionando un m´etodo para calcular el n´umero de divisiones ´unicas para variables categ´oricas nominales. Un criterio para evaluar la bondad de las divisiones Por otra parte, para evaluar la bondad de las divisiones, un criterio est´andar consiste en buscar la m´ınima suma de los cuadrados de los errores dentro de los grupos o nodos. En t´erminos relativos, de forma global, se puede considerar el coeficiente de determinaci´on R2. Una regla de parada Esto se tratar´a m´as adelante cuando se presenten los algoritmos. Aun as´ı, se puede adelantar que, aunque en principio se puede lograr un ´arbol con un nodo final por cada observaci´on, esto resultar´ıa en un modelo muy complejo, con capacidad de generalizaci´on muy baja. Por ello, estas reglas se basan en la profundidad del ´arbol y en el n´umero de casos en cada nodo del mismo. Una regla para predecir en cada nodo terminal Predecir en la regi´on Rtasociada al nodo tcon ˆyRt=1 #(Rt)X i∈Rt yi. Es frecuente que un ´arbol de regresi´on demasiado detallado se ajuste en exceso a los datos, resultando un ajuste insatisfactorio al ser aplicado a nuevos conjuntos de datos; lo que se conoce como sobreajuste. Un ´arbol m´as simple, con menos divisiones, puede ofrecer un mejor equilibrio entre la varianza y sesgo como mencionamos anteriormente, favoreciendo a una mayor generalizaci´on del modelo y a una mayor dependencia continua de los datos. La t´ecnica com´unmente aceptada para controlar el sobreajuste implica la construcci´on inicial de un ´arbol grande T0seguido de un proceso de poda para obtener una versi´on m´as reducida. Un sub´arbol T⊂T0se define como cualquier ´arbol que se pueda derivar de T0eliminando algunas de sus ramas, es decir, suprimiendo nodos terminales seleccionados para formar un ´arbol con menos nodos terminales que el original. 11 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON La cuesti´on que subyace es determinar el m´etodo m´as adecuado para podar el ´arbol, equilibrando la necesidad de simplificaci´on con la de retener la suficiente estructura para capturar las relaciones entre los datos. La respuesta a esta pregunta implica tanto la teor´ıa como la pr´actica de t´ecnicas estad´ısticas y computacionales. 2.1.1.2. Algoritmo b´asico de poda de un ´arbol de decisi´on Con el objetivo de aplicar el algoritmo de poda posteriormente para un ´arbol de supervivencia, presentamos a continuaci´on el algoritmo base en el que se apoyan todas las posibles variaciones algor´ıtmicas a aplicar en ´arboles. Entrada: Conjunto de datos de entrenamiento con la variable objetivo y las predictoras. 1. Construir un ´arbol inicial T0. 2. Hacer T∗=T0yk= 0. 3. Mientras |T∗|>1 (n´umero de nodos) hacer: a)k=k+ 1. b) Seleccionar un nodo ven T∗y podar la rama con ra´ız en v:T∗ v. c) Hacer T(k)=T∗−T∗ v(resultado de la poda del paso anterior). d) Hacer T∗=T(k). 4. Salida: T0, T(1), . . . , T(k):T(1) ⊃T(2) ⊃. . . ⊃T(K). 5. Seleccionar el mejor ´arbol entre todos los incluidos en la colecci´on {T0, T(1), . . . , T(K)}. Para realizar los dos procesos de selecci´on, se emplea un criterio basado en “costes” de error y la “complejidad” del ´arbol. Criterio Coste-Complejidad: “poda por el eslab´on m´as d´ebil” Se busca el sub´arbol ´optimo usando un par´ametro Coste-Complejidad (α) que penaliza la funci´on objetivo con el n´umer de nodos terminales del ´arbol. En lugar de considerar cada posible sub´arbol, se considera una secuencia de ´arboles indexada por dicho par´ametro αno negativo. A cada αle corresponde un ´arbol Tα⊂T0tal que 12 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON C(Tα) = m´ın T⊆T0{Cα(T)} con Cα(T) = |T| X j=1 X xi∈Rj (yi−ˆyRj)2+α|T|, donde |T|es el tama˜no de T(n´umero de nodos terminales) y Rjla regi´on del espacio predictor asociada al j-´esimo nodo terminal. El par´ametro αse introduce para equilibrar el tama˜no o complejidad del ´arbol y el ajuste a los datos. Valores altos de αconducen a ´arboles peque˜nos, y viceversa. A continuaci´on se detallan algunas consideraciones sobre el algoritmo: Cuando α= 0, entonces el sub´arbol Tαser´a simplemente igual T0. Cuando αaumenta, hay una penalizaci´on por la complejidad del ´arbol, por lo que C(T, α) tiende a ser m´as minimizado para un sub´arbol m´as peque˜no. Con el aumento de α, se podan las ramas del ´arbol en una forma anidada y predecible, por lo que la obtenci´on de la secuencia de sub´arboles viables como una funci´on de αes factible. Para cada αexiste un ´unico sub´arbol m´ınimo Tαque minimiza C(T, α). Se puede obtener una familia anidada de sub´arboles donde cada ´arbol sea sub´arbol del anterior y una serie de valores α1< ... < αmtales que T(h)=Tαpara todo α∈[αh, αh+1). Una vez generados los sub´arboles, debe elegirse uno de ellos mediante un estimador insesgado del error (riesgo) esperado de cada sub´arbol, lo que puede realizarse a trav´es de un conjunto test o mediante validaci´on cruzada. Selecci´on del ´arbol podado ´optimo El objetivo es seleccionar uno de los ´arboles de la secuencia de ´arboles podados como el ´arbol ´optimo, que constituir´a el ´arbol soluci´on. Por tanto, se trata de asociar una medida de error o coste C(·) a cada ´arbol Tde la secuencia y escoger aquel que tenga asociado el menor error, o bien, aquel sub´arbol que tenga un error “muy similar” al ´optimo pero con un tama˜no (o complejidad) menor. 13 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Elegir T(k0)tal que: C(T(k0)) = m´ın kC(T(k)) con C(T) = |T| X j=1 X xi∈Rj (yi−ˆyRj)2=C(T, 0), Para estimar “honestamente” C(T(k)), se plantean dos posibles soluciones: Cts(T), estimador por conjunto de prueba (test set), que se calcula dividiendo el conjunto de aprendizaje Den dos subconjuntos disjuntos, aprendizaje DLy test DT. Su aplicaci´on consiste en: 1. Se obtiene un ´arbol T0inicial con DL, y se obtiene la secuencia de ´arboles a trav´es del procedimiento de poda: T0⊇T(1) ⊇T(2) ⊇ . . . ⊇T(k). 2. Con el conjunto DTse obtiene la secuencia Cts(T(k)). El ´arbol ´optimo ser´a T(k0):Cts(T(k0)) = m´ınkCts(T(k)). Ccv(T), estimador por validaci´on cruzada. Suele utilizarse validaci´on cruzada con K= 10 pliegues o iteraciones. En la validaci´on cruzada de Kpliegues (K-fold cross-validation), el conjunto de datos Dse dividen en Ksubconjuntos. Uno de los subconjuntos se utiliza como datos de prueba y el resto (K−1) como datos de entrenamiento. Este proceso de validaci´on cruzada es repetido durante Kiteraciones, con cada uno de los posibles subconjuntos de datos de prueba. En este caso Ccv(T) se obtiene como la media de las Mmedidas obtenidas a trav´es de las iteraciones. Para evitar la inestabilidad de un estimador de este tipo resulta recomendable calibrar la incertidumbre de Cts(T) o Ccv(T) calculando su error est´andar (SE). Se define el error est´andar asociado al estimador Cts(T), o Ccv(T), como: SE(Cts(T)) = sCts(T)1−Cts(T)) |DT|,con |DT|= card{DT}. Los estimadores de C(T) en funci´on del n´umero de nodos terminales, |T(k)|, se comportan de manera similar a la figura: r´apido decrecimiento inicial 14 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Figura 2.2: Estimadores de C(T) en funci´on del n´umero de nodos terminales, |T(k)|. seguido de un largo valle, casi plano, y un posible crecimiento leve para valores altos de |T(k)|. As´ı, el m´ınimo es muy inestable, por lo que, se propone la regla 1 −SE como alternativa a la selecci´on propuesta del m´ınimo con el objetivo de reducir esta inestabilidad a la vez que nos aseguramos que el ´arbol seleccionado sea el m´as simple con valor criterio cercano a dicho m´ınimo. La regla 1−SE consiste en seleccionar el ´arbol m´as simple (menos nodos) que no se desv´ıa m´as de una SE del m´ınimo: Si T(k0):Cts(T(k0)) = m´ın kCts(T(k)) entonces seleccionar T(k1)tal que |T(k1)|= m´ın{|T(k)|:Cts(T(k))≤Cts(T(k0)) + SE(Cts(T(k0)))}. Ver imagen en la p´agina siguiente. 2.1.2. ´ Arboles de clasificaci´on Mientras que en los ´arboles de regresi´on la variable objetivo es cuantitativa, en los ´arboles de clasificaci´on la variable objetivo es cualitativa. B´asi- 15 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Mas a´un, los ´arboles tienden al sobreajuste, lo que puede reducir la capacidad de generalizaci´on del modelo. Tambi´en son conocidos por su inestabilidad ante peque˜nas variaciones en los datos de entrada. No obstante, estas desventajas pueden reducirse o eliminarse mediante t´ecnicas como el bagging, el boosting o la combinaci´on de m´ultiples ´arboles en un modelo de bosque aleatorio (Random Forest). Estos m´etodos tienen como objetivo aumentar la precisi´on y robustez del modelo, mejorando el rendimiento predictivo. 2.2. Introducci´on al an´alisis con datos censurados Introducimos esta secci´on con el objetivo de presentar la idea de censura en los datos, ya que es una caracter´ıstica frecuente durante an´alisis de supervivencia, que se ver´a posteriormente. Como se discute en el trabajo H¨ardle et al. (2006), la caracter´ıstica principal inherente a estos datos es la censura. El t´ermino censura se refiere a datos parcialmente conocidos, a los que llamaremos datos censurados u observaciones censuradas, los cuales pueden surgir por m´ultiples razones: finalizaci´on del estudio antes de que el evento de inter´es ocurra para todos los sujetos, la p´erdida de seguimiento de participantes, o el inicio de la observaci´on despu´es de que el evento haya ocurrido para algunos casos. Distinguimos entre censura por la derecha y censura por la izquierda, seg´un el tiempo de supervivencia exceda o no el tiempo observado. Adem´as, podemos hacer una distinci´on para los datos censurados por la derecha: Censura de Tipo I: las unidades bajo estudio comienzan el estudio en el mismo punto temporal, en una fecha fijada o plazo de tiempo fijado el estudio finaliza y en algunas de ellas no se ha seguido el estudio o no ha ocurrido el evento. No obstante, podemos distinguir entre: •Fija: tiempo inicial y final del estudio es el mismo para todos los sujetos bajo estudio. •Progresiva: se estratifica a la poblaci´on bajo estudio y se le asigna a cada grupo un tiempo de estudio concreto. •Generalizada: se le asigna a cada sujeto un tiempo inicial y final espec´ıfico. 22 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Censura de Tipo II: las unidades bajo estudio comienzan el estudio en el mismo punto temporal, el final del estudio no es fijado inicialmente y es continuado hasta que el evento ocurre para una determinada proporci´on de unidades bajo estudio. Censura de Tipo III: el tiempo de censura lo determina alg´un fen´omeno aleatorio no controlado o esperado, que ocurre durante el transcurso del estudio. Entre las posibles causas se incluyen: p´erdida del seguimiento del sujeto sin mayor informaci´on o abandono voluntario del sujeto. En todos estos casos, el tiempo de censura es aleatorio. Figura 2.4: Tipos de censura. Fuente: H¨ardle et al. (2006). Notemos que la censura de tipo II es una censura de tipo no aleatoria, mientras que la tipo I y tipo III son de tipo aleatoria. Por otro lado, las observaciones censuradas por la izquierda son aquellas en las que el evento de inter´es ha ocurrido antes del inicio del estudio, de modo que el tiempo de censura para estas observaciones ser´a el momento del inicio del estudio, puesto que se conoce que el suceso ha tenido lugar previamente. Esta circunstancia de que el tiempo de supervivencia no puede ser completamente observado para todas los unidades bajo estudio puede ser expresada 23 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON como sigue. Sean YyTel tiempo observado y el tiempo de supervivencia, respectivamente. Sin censura, Y=T, es decir, el tiempo observado es el verdadero tiempo de supervivencia. Con censura, Y < T, es decir, el tiempo observado es el tiempo censurado denotado por U. Podemos emplear un indicador δque nos permitir´a saber si un dato est´a o no censurado, esto es, δ= 1 si Y=Tyδ= 0 en caso contrario. Para el ´ultimo, tendremos que Y=min(T, U). Ejemplo: Si un paciente muere a los 3 a˜nos (T= 3), aunque el estudio termina a los 5 a˜nos, entonces Y=T= 3 (no hay censura). Pero si el paciente sigue vivo al final del estudio, entonces T es desconocido y mayor que 5, y Y=U= 5 (censurado por la derecha). Se tiene que para el an´alisis de este tipo de datos se puede recurrir a m´ultiples t´ecnicas, bien param´etricas, bien semiparam´etricas o no param´etricas. Cada una de ellas presentar´a ventajas e inconvenientes. Nosotros nos centraremos en los ´arboles de supervivencia, t´ecnica que se incluye en el grupo de las no param´etricas, pues no requiere realizar suposiciones sobre la distribuci´on de las funciones de supervivencia. 2.3. ´ Arboles de supervivencia Siguiendo con lo que se discute en H¨ardle et al. (2006), cualquier m´etodo basado en ´arboles incluye dos pasos cruciales: 1. Crecimiento del ´arbol, es decir, particionar los datos (nodos internos) de acuerdo a un criterio de divisi´on oportuno. Este criterio de divisi´on nos permite seleccionar la mejor covariable para dividir cualquier nodo, junto a su respectivo punto de corte. 2. Poda del ´arbol, es decir, la eliminaci´on de algunas ramas con el objetivo de simplificar el ´arbol para obtener conclusiones m´as precisas y claras. Puesto que el tiempo de supervivencia no se observa completamente para todas las unidades cuando estamos trabajando con datos censurados, surgen de manera natural dos variables de respuesta: el tiempo observado Yy el indicador de censura δ. Podemos entonces representar cada dato por una tripleta (yi, δi, xi), donde yies el tiempo observado para la i-´esima unidad experimental, δiindica si el dato yiest´a censurado y xi= (xi1, . . . , xip) es el 24 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON vector de las pcovariables asociadas con la i-´esima unidad experimental. La igualdad yi=tise denomina tiempo de ocurrencia, tiempo de fallo o muerte. Se tiene que, en este enfoque, se asume que la censura es aleatoria (tipos I y III), lo que implica que las distribuciones de tiempo de supervivencia y el tiempo de censura son independientes. 2.3.1. Definiciones previas A continuaci´on se definen las funciones que se emplean en el contexto de an´alisis de supervivencia, y que, al aplicar t´ecnicas de ´arboles se estiman con el objetivo de predecir, como veremos posteriormente. La variable aleatoria bajo estudio ser´a el tiempo Tde supervivencia (tambi´en conocida como tiempo de fallo). Definici´on 1. Se definen la funci´on de densidad de T,f(t) y la de distribuci´on, F(t), como F(t) = P(T≤t) = Zt 0 f(x)dx, f(t) = F′(t). Propiedades de F(t) 1. R∞ 0f(x)dx = 1. 2. F(t) es continua. 3. F(t) es mon´otona creciente para todo ten el soporte de T. Definici´on 2. Se define la funci´on de supervivencia, denotada como S(t), como la probabilidad de que el tiempo de supervivencia Tsea mayor que un tiempo tespec´ıfico. Esto es, S(t) = P(T > t). N´otese que la funci´on de supervivencia es la complementaria de la funci´on de distribuci´on F(t) de T, puesto que F(t) mide la probabilidad de que el sujeto “muera” o “falle” antes de t. Es decir, S(t) = 1 −F(t). Propiedades de S(t) 25 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON 1. S(0) = 1. 2. S(+∞) = 0. 3. S(t) es continua y mon´otona decreciente para todo ten el soporte de T. Definici´on 3. Se define la funci´on de riesgo instant´aneo h(t) como la tasa de fallo en un tiempo tdado que el sujeto ha sobrevivido hasta ese momento sin experimentar el evento de inter´es, h(t) = l´ım ∆t→0 P(t≤T < t + ∆t|T≥t) ∆t. Esta funci´on describe la probabilidad instant´anea de que ocurra el evento en el tiempo t, condicionada a la supervivencia hasta ese momento. Se tiene como propiedad que h(t) = f(t) S(t)=f(t) 1−F(t). Formalmente, la funci´on de riesgo h(t) puede presentarse de manera creciente, decreciente o constante. No obstante, en el contexto del an´alisis de supervivencia esta funci´on suele ser creciente, lo que indicar´ıa que es m´as probable que los sujetos fallezcan con el tiempo. Definici´on 4. Se define la funci´on de riesgo acumulado H(t) como el acumulado de la funci´on de riesgo h(t) desde el comienzo del estudio hasta el tiempo t. Esto es, H(t) = Zt 0 h(u)du. Proporciona una medida acumulativa del riesgo de fallo hasta el tiempo t y es ´util para describir la carga total de riesgo experimentada por un sujeto hasta un tiempo t. En el contexto del an´alisis de supervivencia, al representar H(t), nos aparecer´a una funci´on que crece m´as r´apido que una recta por ser h(t) creciente. Relaciones entre las funciones de riesgo y de supervivencia 1. S′(t) = −f(t). 2. h(t) = −d dt logS(t). 26 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON 3. S(t) = e−H(t)=⇒H(t) = −lnS(t). La demostraci´on de las propiedades aportadas en esta secci´on pueden consultarse en Kalbfleisch y Prentice (2002). Siguiendo con el esquema empleado para presentar los ´arboles de regresi´on y clasificaci´on, se muestran a continuaci´on los elementos necesarios en un ´arbol de supervivencia. Esta vez lo haremos por secciones. 2.3.2. Un conjunto de cuestiones “binarias” para hacer las divisiones Similar a los ´arboles de clasificaci´on y regresi´on, en los ´arboles de supervivencia se utilizan preguntas binarias para determinar las divisiones de los datos. Cada divisi´on se basa en una covariable y un punto de corte, formulando preguntas del tipo “¿Es el valor de esta covariable mayor/menor que un umbral?”. Estas divisiones se emplean para estratificar el conjunto de datos en subgrupos que son m´as homog´eneos en cuanto a la supervivencia de los sujetos. El proceso de divisi´on busca nodos donde las diferencias en la supervivencia, medidas mediante el tiempo hasta el evento de inter´es o el estado de censura, sean estad´ısticamente significativas. 2.3.3. Un criterio para evaluar la bondad de las divisiones El criterio para evaluar la bondad de las divisiones en ´arboles de supervivencia debe considerar la naturaleza de los datos de supervivencia, los cuales incluyen censuras y tiempos de eventos. Los criterios comunes son el estad´ıstico log-rank y la reducci´on de impureza. 2.3.3.1. Estad´ıstico log-rank El estad´ıstico basado en el test log-rank se usa normalmente para analizar datos censurados con el objetivo de comparar la distribuciones de supervivencia de dos grupos diferentes. Este m´etodo se basa en la hip´otesis nula de que no hay diferencia entre las poblaciones en t´erminos de supervivencia, seg´un fue originalmente propuesto por Mantel y Haenszel (1959), que formalizaron este test en un contexto biom´edico. 27 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON Se trata de un test muy potente para calcular las diferencias mencionadas cuando los logaritmos de las curvas de supervivencias se consideran proporcionales. No obstante, suele tener dificultades para detectar diferencias significativas cuando las curvas se cruzan. El estad´ıstico log-rank es fundamental en la construcci´on de ´arboles de supervivencia para comparar las funciones de supervivencia entre dos grupos de datos. Este m´etodo se basa en evaluar si existen diferencias estad´ısticamente significativas entre los grupos en t´erminos de supervivencia. Supongamos que deseamos comparar la supervivencia de dos grupos A1 yA2en un nodo del ´arbol. Definimos la siguiente tabla de contingencia para cada tiempo tjde fallo observado en los datos: Tabla 2.1: Tabla de contingencia para el estad´ıstico log-rank Grupo Muertes No muertes En riesgo A1d1jn1j−d1jn1j A2d2jn2j−d2jn2j Total dj=d1j+d2jnj−dj= (n1j+n2j)−(d1j+d2j)nj=n1j+n2j Definimos el contraste de hip´otesis del test log-rank de la siguiente manera: H0:SA1(t) = SA2(t) contra H1:SA1(t)=SA2(t), donde SA1(t) y SA2(t) son las funciones de supervivencia para los grupos A1 yA2, respectivamente. Bajo la hip´otesis nula, es decir, cuando las funciones de supervivencia coinciden, la muerte es independiente del grupo. Dados nj(n´umero de casos en riesgo en el instante tj), nj1(n´umero de casos del grupo A1) y dj(n´umero de fallos o muertes en el instante tj), la variable Vj= n´umero de muertes en A1en el instante tjverifica Vj|nj,nj1,dj∼H[nj, nj1, dj]. Por ello, su valor esperado es e1j=n1j dj nj . Y su varianza 28 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON V ar[Vj] = nj1nj2dj(nj−dj) n2 j(nj−1) . Adem´as las distribuciones Vj|nj,nj1,djson asint´oticamente independientes ∀j= 1, . . . , K. Podemos definir los fallos observados menos los esperados: u1= k X j=1 (d1j−e1j). En el caso de que el valor de ksea lo suficientemente grande se tiene por el teorema central del l´ımite que el siguiente cociente se aproxima a una N(0,1): u √ν=Pk j=1(d1j−eij) qPk j=1 vj . donde νj=nj1nj2dj(nj−dj) n2 j(nj−1) . Tambi´en podemos construir el test del estad´ıstico log-rank como: u2 √ν H0 ∼χ2 1. Ejemplo: Imaginemos que estamos construyendo un ´arbol de supervivencia para analizar datos de pacientes con una enfermedad espec´ıfica. Usar el criterio de log-rank para dividir los nodos significa que buscaremos dividir los pacientes en grupos que muestren diferencias significativas en sus patrones de supervivencia, por ejemplo, bas´andonos en su edad, tipo de tratamiento, etc. 2.3.3.2. Reducci´on de impureza Se utiliza una medida de impureza introducida por Breiman et al. (1984) adaptada que combina el tiempo observado y el estado de censura de los datos. Se selecciona la divisi´on que maximiza la disminuci´on de la impureza, lo que refleja un aumento en la homogeneidad de los tiempos de supervivencia dentro de los nodos hijos. Con el fin de introducir el criterio de divisi´on basado en la impureza, es conveniente recordar algunos conceptos b´asicos y notaci´on en la clasificaci´on 29 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON de una respuesta con varias posibles categor´ıas. Consideremos un criterio de divisi´on sde un nodo vcandidato a ser dividido en dos “hijos” (o “subnodos”) vlyvr, basado en una regla que involucra una o m´as variables, y sean p(vl) y p(vr) las proporciones de las observaciones de vque han ca´ıdo en los nuevos nodos hijos vlyvr, respectivamente. La impureza en el nodo v, denotada i(v), mide la impureza basada en las probabilidades entre clases. Entonces, parece natural evaluar el rendimiento de la divisi´on como el cambio en la impureza dado por: △i(s, v) = i(v)−(p(vl)i(vl) + p(vr)i(vr)). Esta cantidad se usa como criterio de partici´on. No obstante, esta noci´on de impureza en el caso de datos censurados no conviene usarla tal y como la hemos presentado, pues aunque estamos interesados en el tiempo de supervivencia, recordemos que est´an involucradas dos variables respuesta: el tiempo observado yi(continua) y el indicador de censura δi(binaria). Entonces, un nodo inicial contendr´a unidades cuyos tiempos observados sean similares y que est´an censurados o no censurados. B´asicamente, una medida de impureza razonable para datos censurados debe tener en cuenta tanto el tiempo observado como la censura. De esta forma, la impureza de un nodo la podemos expresar como una suma ponderada de medidas de impureza basadas en el tiempo observado y en la censura, es decir, i(v) = w1(v)iy(v) + w2(v)iδ(v), donde w1(v) y w2(v) son pesos pre-fijados y iy(v) y iδ(v) denotan la impureza del nodo vpara el tiempo observado y la censura, respectivamente. En particular, puede considerar la impureza para el tiempo como una medida del error cuadr´atico iy(v) = n(v) X i=1 (yi−¯y(v))2 Py2 i , donde n(v) es el n´umero de observaciones en el nodo vy ¯y(v) es la media del tiempo en las observaciones del nodo. El denominador se usa para normalizar respecto a la otra componente de la impureza. Para medir la impureza del indicador de censura, usamos la medida de entrop´ıa: iδ(v) = −pvlog(pv)−(1 −pv)log(1 −pv), 30 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON donde pvdenota la proporci´on de datos censurados en el nodo v. De entre todos los criterios de divisi´on posible de un nodo, se escoge uno para maximizar la reducci´on de impureza medida por la f´ormula inicial i(v). Un buen criterio de divisi´on reducir´a la impureza, lo que significa que las observaciones en cada nodo hijo son m´as “puras” (es decir, m´as propensas a pertenecer a una sola clase). As´ı, si el nodo original vtiene una alta mezcla de clases (alta impureza) y despu´es de la divisi´on, los nodos hijos vlyvrson cada uno muy homog´eneos (baja impureza), entonces la reducci´on de impureza △i(s, v) ser´a grande, indicando que hemos usado un buen criterio de divisi´on. En resumen, esta f´ormula nos ayuda a cuantificar cu´anto mejora la homogeneidad (o pureza) de los nodos despu´es de aplicar un criterio de divisi´on. En la construcci´on de ´arboles de decisi´on, buscamos continuamente divisiones que maximicen este cambio en impureza, ya que eso mejora la capacidad del ´arbol para distinguir claramente entre las diferentes clases. Tanto el criterio del estad´ıstico log-rank como el de la reducci´on de la impureza buscan mejorar la homogeneidad interna de los nodos con respecto a la supervivencia. Sin embargo, es fundamental que la medida de impureza refleje adecuadamente la heterogeneidad en los tiempos de supervivencia y el estado de censura, lo que puede requerir la ponderaci´on o la combinaci´on de diferentes componentes de impureza para lograr este objetivo. 2.3.4. Una regla de parada En los ´arboles de supervivencia, el criterio de parada es fundamental para asegurar que el modelo final sea tanto representativo de los datos como predictivo para nuevas observaciones. Estos son los aspectos clave que se consideran para detener el crecimiento del ´arbol: 1. N´umero m´ınimo de observaciones: Se puede establecer un l´ımite para el n´umero m´ınimo de observaciones que cada nodo debe tener, especialmente en los nodos terminales. Esto asegura que cada grupo tenga suficientes datos para que la estimaci´on de la supervivencia sea confiable. 2. Complejidad del ´arbol: Un ´arbol demasiado complejo puede llevar al sobreajuste. Una m´etrica de complejidad, como el coste-complejidad de Breiman et al. (1984), puede utilizarse para equilibrar la precisi´on y la simplicidad. 31 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON R(ti) es el conjunto de riesgo en el tiempo ti, es decir, el conjunto de todos los sujetos que a´un est´an bajo estudio justo antes del tiempo ti. ˆ βson los coeficientes estimados del modelo. Xjson las covariables para el sujeto j. Este m´etodo asume que los tiempos de fallo son discretos y que los fallos ocurren en los tiempos observados. La estimaci´on de h0(t) se puede obtener diferenciando ˆ H0(t), es decir, calculando la diferencia entre estimaciones consecutivas de ˆ H0(t). Una vez estimada la parte param´etrica exp(ˆ β′X) y la no param´etrica ˆ h0(t), tendremos la estimaci´on del modelo de Cox completo: ˆ h(t, X) = ˆ h0(t)exp(ˆ β′X). 2.4.3. Integraci´on del modelo de Cox en ´arboles de supervivencia A continuaci´on, describimos c´omo se puede integrar este modelo en la metodolog´ıa de ´arboles de supervivencia. En ´arboles de supervivencia, cada nodo terminal puede estar asociado con un modelo de Cox espec´ıfico, donde las covariables son aquellas que definen la ruta desde la ra´ız hasta el nodo. La funci´on de riesgo en cada nodo terminal, por lo tanto, se modela como: h(t, X) = h0,nodo(t) exp(β1X1+β2X2+. . . +βpXp), donde h0,nodo(t) es la funci´on de riesgo base para el nodo espec´ıfico, que podr´ıa ser estimada utilizando los datos que caen dentro de ese nodo, reflejando una estrategia m´as localizada en comparaci´on con el enfoque global habitual del modelo de Cox. La funci´on de riesgo base en cada nodo del ´arbol, h0,nodo(t), puede ser estimada utilizando el estimador de Breslow, adaptado a las observaciones en ese nodo. Esto implica calcular: ˆ H0,nodo(t) = X ti≤t di Pj∈Rnodo(ti)exp(ˆ β′Xj), 38 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON donde Rnodo(ti) es el conjunto de riesgo en el tiempo ti, pero limitado a las observaciones que terminan en el nodo en cuesti´on. Este enfoque permite que cada nodo del ´arbol adapte su propia estimaci´on de la funci´on de riesgo base bas´andose en su subconjunto de datos. Integrar el modelo de Cox en los ´arboles de supervivencia permite obtener una interpretaci´on de los efectos de las covariables en la supervivencia dentro de subgrupos espec´ıficos definidos por las reglas del ´arbol. Esto es particularmente ´util en estudios m´edicos donde las interacciones entre variables pueden tener impactos significativos en los resultados de supervivencia. Adem´as, esta metodolog´ıa facilita la visualizaci´on de la igualdad de los efectos del tratamiento o de otros factores de riesgo a trav´es de diferentes subgrupos en la poblaci´on de estudio. Definimos ahora el coeficiente de Cox (Cox Score) para un modelo de ´arbol de Cox como la media ponderada de las log-verosimilitudes parciales m´aximas de los modelos de Cox en cada uno de los nodos terminales del ´arbol: CS(T) = 1 MX j∈nodos terminales mjmaxβjl(βj|T), donde mjes el n´umero de elementos del nodo terminal jyMes el n´umero total de elementos en el conjunto de datos. Para facilitar la interpretaci´on, tambi´en definimos el cociente del coeficiente Cox (CSR) como la reducci´on porcentual en el coeficiente para el ´arbol Ten relaci´on a un modelo nulo, CSR(T) = 1 −CS(T) maxβl(β|T0). El CSR compara el ajuste de un ´arbol de decisi´on espec´ıfico, T, con el de un modelo nulo, T0, que no incluye ninguna de las covariables predictivas; es decir, es el modelo m´as simple posible. Un CSR cercano a 0 indicar´ıa que la estructura del ´arbol no mejora el ajuste del modelo m´as all´a de lo que se podr´ıa obtener con el modelo nulo. Por otra parte, un CSR cercano a 1 indicar´ıa que el modelo de ´arbol Testar´ıa capturando mejor las relaciones entre las covariables y la variable objetivo, comparado con el modelo nulo. 39 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON 2.4.4. Residuos Llevar a cabo an´alisis de residuos ha sido una de las ventajas que han surgido a trav´es de an´alisis de supervivencia, ya que ´estos se pueden emplear para: 1. Determinar la forma funcional apropiada para un predictor continuo, lo que implica evaluar c´omo se relaciona cada predictor con la probabilidad de ocurrencia del evento de inter´es. 2. Identificar los individuos que pueden estar mal pronosticados por el modelo, lo cual es crucial para mejorar la precisi´on y la generalidad del modelo. 3. Identificar puntos o individuos de influencia, cuya presencia puede afectar de manera significativa los estimados del modelo, sugiriendo la necesidad de investigar m´as a fondo o posiblemente de ajustar el modelo. 4. Verificar la suposici´on de riesgo proporcional, asegurando que la relaci´on entre las covariables y la funci´on de riesgo se mantenga constante a lo largo del tiempo. Therneau et al. (1990) sugieren utilizar los residuos martingala como la variable respuesta en un algoritmo de ´arbol de regresi´on. Los residuos martingala, en el contexto del modelo de Cox, se definen para un individuo ien el tiempo tcomo la diferencia entre el evento observado y la integral del riesgo acumulado hasta ese momento, ajustado por las covariables. Para el individuo i, el residuo martingala Mise puede expresar matem´aticamente como: Mi=δi−Zti 0 Yi(u)h(u, Xi)du, donde δies el indicador de evento, siendo 1 si el evento de inter´es ocurri´o y 0 en caso contrario. ties el tiempo de observaci´on o censura para el individuo i. Yi(t) es una variable indicadora que vale 1 si el individuo iest´a en riesgo en el tiempo ty 0 en caso contrario. 40 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON h(u, Xi) es la funci´on de riesgo en el tiempo upara un individuo con covariables Xi, como se defini´o anteriormente. Esto permite que los m´etodos est´andar de ´arboles de regresi´on se apliquen directamente a los datos de supervivencia, utilizando el software existente de ´arboles de regresi´on con los residuos martingala como resultados modificados. Esta aproximaci´on facilita la identificaci´on de relaciones entre variables y de riesgos en el an´alisis de supervivencia utilizando la metodolog´ıa de ´arboles de decisi´on. Existen otro tipo de residuos que pueden ser empleados, y en los que no entraremos en detalle en este trabajo: Residuos de Cox-Snell. Residuos basados en el estad´ıstico Deviance. Residuos de Schoenfeld. Residuos escalados de Schoenfeld. Residuos dfbeta. 41 CAP´ ITULO 2. ´ ARBOLES DE DECISI´ ON 42 Cap´ıtulo 3 Bosques aleatorios de supervivencia Los bosques de decisi´on son una extensi´on natural de los ´arboles de decisi´on y proporcionan una mejora significativa en t´erminos de precisi´on y robustez al combinar m´ultiples ´arboles. De manera similar, los bosques de supervivencia, que aplican el enfoque de los bosques de decisi´on al an´alisis de supervivencia, permiten una mejor modelizaci´on de los datos censurados y proporcionan estimaciones m´as estables y precisas de las funciones de riesgo y supervivencia. Esta inclusi´on proporciona una visi´on m´as completa del campo y destaca la evoluci´on y el potencial de las t´ecnicas basadas en ´arboles para el an´alisis de datos de supervivencia. Antes de adentrarnos en el estudio de los bosques de supervivencia, es esencial entender los fundamentos y mecanismos de los bosques de decisi´on, tambi´en conocidos como Random Forests. Este modelo es una evoluci´on natural de los ´arboles de decisi´on y constituye la base sobre la cual se construyen los bosques de supervivencia. Los bosques de decisi´on son una t´ecnica de aprendizaje ensemble que mejora la capacidad predictiva y la estabilidad de los modelos mediante la combinaci´on de m´ultiples ´arboles de decisi´on durante el proceso de entrenamiento. Cada ´arbol en el bosque es construido a partir de una muestra aleatoria del conjunto de datos de entrenamiento, utilizando un subconjun- to de las variables en cada divisi´on, lo que asegura la diversidad dentro del modelo y reduce el riesgo de sobreajuste. Este m´etodo no solo es potente por su robustez y precisi´on, sino tambi´en por su capacidad para manejar grandes vol´umenes de datos y su eficacia 43 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA frente a conjuntos de datos con muchas variables explicativas, posiblemente correlacionadas. Adem´as, los bosques de decisi´on proporcionan una medida intr´ınseca de la importancia de las variables, lo que permite una evaluaci´on efectiva de qu´e factores son m´as significativos al predecir la variable respuesta. La introducci´on y comprensi´on de los bosques de decisi´on son fundamentales para abordar los bosques de supervivencia, ya que comparten muchos principios t´ecnicos pero se aplican en contextos diferentes, especialmente en an´alisis donde el tiempo hasta un evento es crucial, como es el caso en los estudios de supervivencia. En la siguiente secci´on, exploraremos c´omo los principios de los bosques de decisi´on se adaptan y extienden en los bosques de supervivencia para abordar de manera eficiente y efectiva los an´alisis de tiempo hasta el fallo o muerte. Como idea intuitiva, un modelo Random Forest est´a formado por un conjunto de ´arboles de decisi´on individuales, cada uno entrenado con una muestra aleatoria que se extrae de los datos de entrenamiento originales mediante la t´ecnica bootstrap. Por tanto, cada ´arbol se entrena con unos datos ligeramente distintos. En cada ´arbol individual, las observaciones se van distribuyendo por los nodos hasta alcanzar un nodo terminal. Como cabe esperar, la predicci´on de una nueva observaci´on se obtiene agregando las predicciones de todos los ´arboles individuales que forman el modelo. Para entender c´omo funcionan los modelos Random Forest es necesario conocer primero los conceptos de ensemble ybagging. Los contenidos que se siguen en esta secci´on se tratan en el art´ıculo de Rodrigo (2020) 3.1. M´etodos de ensemble Todos los modelos referentes al aprendizaje autom´atico y la estad´ıstica, enfrentan el desaf´ıo com´un del manejo del equilibrio entre el sesgo y la varianza, caracter´ısticas elementales para medir la eficacia de cualquier modelo predictivo. El sesgo representa el error introducido por aproximar un problema real y su soluci´on mediante un modelo simplificado. Indica la diferencia entre las predicciones promedio de un modelo y los valores verdaderos. Adem´as, refleja c´omo de capaz es el modelo de aprender la relaci´on real que existe entre la variable objetivo y las predictoras. Por ejemplo, si la relaci´on sigue un patr´on no lineal, al aplicar un modelo de regresi´on lineal se tendr´a un 44 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA sesgo alto puesto que el modelo no podr´a capturar correctamente la relaci´on. Por otro lado, la varianza se refiere a cu´anto cambia el modelo las predicciones cuando se entrena con diferentes conjuntos de datos. Un modelo ideal debe ser estable y no variar demasiado con peque˜nos cambios en el conjun- to de entrenamiento. Sin embargo, modelos demasiado complejos, como un ´arbol de decisi´on profundo, pueden ajustarse excesivamente a los datos de entrenamiento, memorizando en lugar de aprender, lo que se traduce en una alta varianza. Incrementar la complejidad de un modelo suele disminuir el sesgo, ya que el modelo puede ajustarse mejor a los datos. No obstante, esto puede llevar a un sobreajuste, donde el modelo se adapta tan estrechamente a los datos de entrenamiento que falla al predecir nuevos datos. Es por ello que se habla de equilibrio entre sesgo y varianza. Los m´etodos de ensemble combinan m´ultiples modelos en uno nuevo con el objetivo de lograr un equilibrio entre sesgo y varianza. Dos de los tipos de ensemble m´as utilizados son: Bagging: se ajustan m´ultiples modelos, cada uno con un subconjunto distinto de los datos de entrenamiento. Finalmente, para predecir se toman todos los modelos que forman el agregado y se toma la media de todas las predicciones en el caso de las variables continuas, o la clase m´as frecuente en el caso de variables categ´oricas. Los modelos Random Forest est´an dentro de esta categor´ıa. Boosting: modelos sencillos se ajustan secuencialmente, de manera que cada modelo “aprende” de los errores del anterior. Como valor final se toma la media de todas las predicciones (variables continuas) o la clase m´as frecuente (variables cualitativas). Tres de los m´etodos de boosting m´as empleados son AdaBoost,Gradient Boosting yStochastic Gradient Boosting, aunque no entraremos en detalle pues no forman parte de Random Forest. En bagging, se emplean modelos con muy poco sesgo pero mucha varianza, de forma que agreg´andolos se consigue reducir la varianza sin apenas incrementar el sesgo. Adem´as, cada modelo es distinto del resto porque cada uno se entrena con una muetra distinta obtenida mediante bootstrapping. A continuaci´on, se describe con m´as detalle la t´ecnica de bagging. 45 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA 3.2. Bagging El t´ermino bagging es el diminutivo de bootstrap aggregation, y hace referencia al empleo del muestreo repetido con reposici´on bootstrapping, con el objetivo de reducir la varianza de algunos modelos de aprendizaje estad´ıstico, entre ellos los basados en ´arboles. Un resultado cl´asico de inferencia estad´ıstica es que dadas nmuestras de observaciones independientes X1, . . . , Xn, con var(Xi) = σ2,∀i= 1, . . . , n, la varianza de la media de las observaciones ¯ Xes σ2 n. Es decir, se consigue reducir la varianza mediante el promedio de un conjunto de observaciones. En esta idea se basa el bagging: se simulan pseudomuestras mediante bootstrapping, con los que se ajustan diferentes modelos y despu´es se agregan, reduciendo la varianza. La forma de aplicarlo es la siguiente: 1. Generar Bconjuntos de prueba mediante bootstrapping a partir de la muestra de entrenamiento original. 2. Entrenar un ´arbol con cada una de las Bmuestras del paso 1. Cada ´arbol se crea sin apenas restricciones y no se somete a poda, por lo que tiene varianza alta pero poco sesgo. En la mayor´ıa de casos, la ´unica regla de parada es el n´umero m´ınimo de observaciones que deben tener los nodos terminales. El valor ´optimo de este hiperpar´ametro puede obtenerse comparando el out-of-bag error (se ver´a posteriormente) o por validaci´on cruzada. 3. Para cada nueva observaci´on, obtener la predicci´on de cada uno de los B´arboles. El valor final de la predicci´on se obtiene como la media de las Bpredicciones en el caso de variables cuantitativas y como la clase predicha m´as frecuente (moda) para variables cualitativas. En an´alisis de supervivencia, podr´ıa ser la media de las Bestimaciones de la funci´on de supervivencia. 46 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA 3.3. Entrenamiento de Random Forest Conviene llegados a este punto aclarar que la t´ecnica de bagging para reducir varianza mediante el promedio de modelos se consigue siempre cuando los modelos agregados no tengan una alta correlaci´on. Por ejemplo, si en un conjunto de datos tenemos una variable predictora con una alta influencia, todos o casi todos los ´arboles creados en el proceso de bagging tendr´an como dominador esta variable. En consecuencia, los ´arboles ser´an parecidos entre ellos, no consiguiendo reducir la varianza. La t´ecnica de Random Forest evita este problema al seleccionar aleatoriamente mvariables predictoras antes de evaluar cada divisi´on. De esta forma, un promedio de p−m pdivisiones no contendr´an la variable influyente. A˜nadiendo este paso extra se consigue reducir la correlaci´on entre ´arboles a´un m´as. El algoritmo de Random Forest es id´entico al de bagging con la diferencia de que antes de cada divisi´on se seleccionan aleatoriamente mvariables predictoras. Si m=pentonces las dos t´ecnicas son equivalentes. Se recomiendan: La ra´ız cuadrada del n´umero total de predictores para problemas de clasificaci´on: m≈√p. Un tercio del n´umero de predictores para regresi´on: m≈p 3. Aun as´ı, la mejor forma para encontrar el valor ´optimo de mes evaluar el out-of-bag-error o recurrir a validaci´on cruzada. 3.4. Out-of-Bag Error El proceso de bagging ofrece la posibilidad de estimar el error de test sin necesidad de emplear m´etodos de validaci´on cruzada. Esto se debe a que, al ajustar los ´arboles utilizando muestras generadas por bootstrapping, en promedio, cada ajuste utiliza aproximadamente dos tercios de las observaciones originales, dejando un tercio restante denominado out-of-bag (OOB). Registrando las observaciones empleadas para cada ´arbol ajustado, es posible predecir la respuesta de una observaci´on excluida haciendo uso de los ´arboles en los que no particip´o su ajuste y promedi´andolos (o calculando la moda en el caso de ´arboles de clasificaci´on). Este proceso permite obtener predicciones para todas las observaciones y calcular el OOB-mean square error 47 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA 3.6.3.2. Error de predicci´on OOB Calcular Crequiere calcular las predicciones. Utilizamos la estimaci´on de la funci´on de riesgo acumulada del conjunto OOB para definir las predicciones. Debido a que este valor se deriva de los datos OOB, puede usarse para obtener una estimaci´on OOB para C, y, en consecuencia, una tasa de error OOB. Sean to 1, . . . , to mlos puntos de tiempo ´unicos preseleccionados (usamos los tiempos de eventos ´unicos, t1, . . . , tN). Para clasificar dos casos iyj, decimos que itiene un peor resultado predicho que jsi m X l=1 H∗∗ e(to l|xi)> m X l=1 H∗∗ e(to l|xj). Usando esta regla, se calcula Ccomo se describi´o anteriormente. Denotamos la estimaci´on OOB por C∗∗. El error de predicci´on OOB, PE∗∗, se define como 1 −C∗∗. N´otese que 0 ≤PE∗∗ ≤1 y que un valor de PE∗∗ = 0,5 indica que la predicci´on no mejora el “adivinar del azar”. En la pr´oxima secci´on, presentaremos ejemplos pr´acticos sobre c´omo implementar ´arboles y bosques de supervivencia en el software R, destacando tanto su aplicabilidad como su implementaci´on paso a paso. 3.6.4. Funci´on de supervivencia S(t) Se sigue ahora un an´alisis similar al de la secci´on anterior para H(t). La estimaci´on de S(t) para un nodo terminal h∈ T se puede obtener usando la relaci´on entre la funci´on de supervivencia y la funci´on de riesgo acumulada: ˆ Sh(t) = exp −ˆ Hh(t). Todos los casos dentro de htienen la misma estimaci´on de Sh(t). Cada caso itiene una covariable xide dimensi´on d. Sea S(t|xi) la funci´on de supervivencia para i. Para predecir este valor, hacemos descender xipor el ´arbol. Debido a la naturaleza binaria de un ´arbol de supervivencia, xicaer´a en un nodo terminal ´unico h∈ T . La funci´on de supervivencia del sujeto i es el estimador de supervivencia para el nodo terminal de xi: ˆ S(t|xi) = ˆ Sh(t),si xi∈h. 54 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA Esta identidad define la funci´on de supervivencia para todos los casos, y para el ´arbol en s´ı. Para obtener una estimaci´on del conjunto, promediamos sobre B´arboles de supervivencia. Describimos tanto una estimaci´on OOB como bootstrap. Recordemos que cada ´arbol en el bosque se construye usando una muestra bootstrap independiente. Definimos Ii,b = 1 si ies un caso OOB para b; de lo contrario, definimos Ii,b = 0. Sea ˆ S∗ b(t|x) la estimaci´on para un ´arbol construido a partir de la muestra bootstrap b. La estimaci´on de la funci´on de supervivencia para el conjunto OOB para ies ˆ S∗∗ e(t|xi) = PB b=1 Ii,b ˆ S∗ b(t|xi) PB b=1 Ii,b . Observemos que la identidad anterior es un promedio sobre muestras bootstrap en las que ies OOB. Equivalente, S∗∗ e(t|xi) se puede calcular de la siguiente manera. Hacemos descender los datos OOB por un ´arbol de supervivencia construido a partir de datos in-bag (bootstrap). Llevamos a cabo un registro del nodo terminal de iy su predicci´on de la funci´on de supervivencia. El siguiente paso ser´ıa promediar estas predicciones. Esto producir´ıa la identidad anterior. En contraste con la identidad anterior, la predicci´on del conjunto bootstrap para ies S∗ e(t|xi) = 1 B B X b=1 S∗ b(t|xi). Observemos que usa todos los ´arboles de supervivencia y no solo aquellos donde ies OOB. Cabe destacar que en lugar de la relaci´on ˆ Sh(t) = exp −ˆ Hh(t), en el caso de no tener calculada ˆ Hh(t), se puede emplear en su lugar el estimador de Kaplan-Meier (de la misma forma que se emple´o el estimador de Nelson-Aalen en la secci´on anterior para Hh(t)): ˆ Sh(t) = Y tl,h≤t1−dl,h Yl,h . 55 CAP´ ITULO 3. BOSQUES ALEATORIOS DE SUPERVIVENCIA 3.6.5. Error de predicci´on para S(t) La estimaci´on del error de predicci´on para S(t) se realiza de manera an´aloga al de H(t). Se basa en el ´ındice C. De esta forma, de manera an´aloga, sean to 1, . . . , to mlos puntos de tiempo ´unicos preseleccionados (usamos los tiempos de eventos ´unicos, t1, . . . , tN). Para clasificar dos casos iyj, decimos que itiene un peor resultado predicho que jsi m X l=1 S∗∗ e(to l|xi)< m X l=1 S∗∗ e(to l|xj). Usando esta regla, se calcula Ccomo se describi´o anteriormente. Denotamos la estimaci´on OOB por C∗∗. El error de predicci´on OOB, PE∗∗, se define como 1 −C∗∗. N´otese que 0 ≤PE∗∗ ≤1 y que un valor de PE∗∗ = 0,5 indica que la predicci´on no mejora el “adivinar del azar”. 56 Cap´ıtulo 4 Ilustraciones con R Program 4.1. ´ Arbol de supervivencia En esta secci´on veremos un caso pr´actico de ´arbol de supervivencia con el paquete LongCART desarrollado y actualizado por ´ultima vez por Kundu (2022). Este paquete construye ´arboles para datos longitudinales continuos y datos de supervivencia utilizando covariables basales como variables de partici´on seg´un los algoritmos ‘LongCART’ y ‘SurvCART’, respectivamente. Posteriormente, tambi´en se incluyen funciones para calcular la potencia condicional y la potencia predictiva de ´exito basadas en resultados intermedios y la probabilidad de ´exito para un ensayo prospectivo. Vamos a usar el conjunto de datos GBSCG2 presentado por Sauerbrei y Royston (1999) y disponible en el paquete TH.data desarrollado por Hothorn (2023). Se trata de un estudio alem´an sobre c´ancer de mama. Tenemos un conjunto de 686 observaciones con 10 variables. Las variables medidas son: horTh: Terapia hormonal, un factor con niveles no ys´ı. age: Edad en a˜nos. menostat: Estado menop´ausico, un factor con niveles Pre yPost. tsize: Tama˜no del tumor (en mm). tgrade: Un factor ordenado con niveles I <II <III. pnodes: N´umero de ganglios positivos. 57 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM progrec: Receptor de progesterona (en fmol). estrec: Receptor de estr´ogeno (en fmol). time: Tiempo de supervivencia libre de recurrencia (en d´ıas). cens: Indicador de censura (0 - censurado, 1 - evento). A continuaci´on se detalla paso a paso el proceso seguido en R para el an´alisis de supervivencia mediante ´arbol de decisi´on para este conjunto de datos. # Paso 0. Importar las librer´ıas necesarias library(LongCART) library(TH.data) library(survival) library(rpart) library(survminer) # Paso 1. Cargar el conjunto de datos. data(GBSG2) str(GBSG2) 'data.frame': 686 obs. of 10 variables: $ horTh : Factor w/ 2 levels "no","yes": 1 2 2 2 1 1 2 1 1 1 ... $ age : int 70 56 58 59 73 32 59 65 80 66 ... $ menostat: Factor w/ 2 levels "Pre","Post": 2 2 2 2 2 1 2 2 2 2 ... $ tsize : int 21 12 35 17 35 57 8 16 39 18 ... $ tgrade : Ord.factor w/ 3 levels "I"<"II"<"III": 2 2 2 2 2 3 2 2 2 2 ... $ pnodes : int 3 7 9 4 1 24 2 1 30 7 ... $ progrec : int 48 61 52 60 26 0 181 192 0 0 ... $ estrec : int 66 77 271 29 65 13 0 25 59 3 ... $ time : int 1814 2018 712 1807 772 448 2172 2161 471 2014 ... $cens :int 1111110010... # Paso 2. Convertir variables categ´oricas a num´ericas y establecer una columna 'ID'. GBSG2$horTh1 <- as.numeric(GBSG2$horTh) GBSG2$tgrade1 <- as.numeric(GBSG2$tgrade) GBSG2$menostat1 <- as.numeric(GBSG2$menostat) GBSG2$subjid <- 1:nrow(GBSG2) # Paso 3. Crear el modelo SurvCART. 58 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM DESCRIPCI´ ON DE LOS PAR´ AMETROS # Especifica el conjunto de datos que se va a utilizar. En este caso, GBSG2, que contiene informaci´on sobre pacientes con c´ancer de mama. data=GBSG2 # Especifica la columna que contiene los identificadores ´unicos de los pacientes. Aqu´ı, la columna se llama subjid. patid="subjid" # Especifica la columna que indica el estado de censura. Esta columna contiene valores binarios (0 para censurados y 1 para evento). censorvar="cens" # Especifica la columna que contiene los tiempos de seguimiento hasta el evento o censura (en d´ıas en este conjunto de datos). timevar="time" # Lista de variables de partici´on que se utilizar´an para dividir los nodos del ´arbol. Estas variables no cambian con el tiempo y pueden ser categ´oricas o continuas. gvars=c('horTh1','age','menostat1','tsize','tgrade1','pnodes', 'progrec','estrec') # Tipos de las variables de partici´on especificadas en gvars. Un valor de 0 indica una variable categ´orica y un valor de 1 indica una variable continua. La longitud de tgvars debe coincidir con la de gvars. tgvars=c(0,1,0,1,0,1,1,1) # Indica el valor en censorvar que representa la ocurrencia de un evento (en lugar de la censura). Aqu´ı, 1 indica que ocurri´o un evento. event.ind=1 # Nivel de significancia nominal para la prueba de inestabilidad de par´ametros. Este es el umbral para determinar si una divisi´on en el ´arbol es estad´ısticamente significativa. alpha=0.05 # N´umero m´ınimo de observaciones que deben existir en un nodo para que se intente una divisi´on. Esto ayuda a prevenir la sobreajuste al asegurarse de que los nodos no sean demasiado peque~nos. 59 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM minsplit=80 # N´umero m´ınimo de observaciones en cualquier nodo terminal. Similar a minsplit, esto ayuda a prevenir nodos terminales muy peque~nos. minbucket=40 # Si se establece en TRUE, imprime un resumen del modelo ajustado, incluyendo el n´umero de sujetos en riesgo, el n´umero de eventos, el tiempo mediano de eventos y el tiempo mediano de censura para cada nodo del ´arbol. print=TRUE # MODELO. Suponiendo que el tiempo hasta el evento, es decir, el tiempo hasta que se detecta c´ancer, sigue una exponencial. out <- SurvCART(data=GBSG2, patid="subjid", censorvar="cens", timevar="time", gvars=c('horTh1','age','menostat1','tsize', 'tgrade1','pnodes','progrec','estrec'), tgvars=c(0,1,0,1,0,1,1,1), event.ind=1, alpha=0.05, minsplit=80, minbucket=40, print=TRUE) ID n D Q1.T median.T Q3.T Q1.C median.C Q3.C 1 686 299 727 1807 NA 1100 1645 1976 2 583 226 857 2030 NA 1113 1645 1979 3 376 119 1164 NA NA 1152 1632 1984 4 239 91 918 1990 NA 1091 1604 1976 5 137 28 2039 NA NA 1169 1692 2009 6 86 24 1279 NA NA 1119 1624 1984 7 51 4 NA NA NA 1329 1743 2009 8 207 107 577 1337 2456 1090 1666 1979 9 117 72 417 893 2456 974 1637 1977 10 90 35 982 2015 NA 1095 1702 2010 11 103 73 420 747 1493 1062 1722 1846 12 60 51 359 515 855 973 1826 1856 60 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM ID n D Q1.T median.T Q3.T Q1.C median.C Q3.C 13 43 22 730 1170 1975 1062 1505 1834 yval loglik AIC var index p Instability Terminal 1807 -2647.8 5297.6 pnodes 10 0.000 23.8 FALSE 2030 -2039.8 4081.6 pnodes 4 0.000 13.8 FALSE NA -1106.8 2215.7 progrec 90 0.033 6.9 FALSE 1990 -824.9 1651.8 tgrade1 NA 0.936 NA TRUE NA -275.0 552.1 horTh1 2 0.023 4.9 FALSE NA -226.6 455.1 tsize NA 0.928 NA TRUE NA -43.6 89.1 NA NA NA NA TRUE 1337 -919.2 1840.4 progrec 56 0.008 7.8 FALSE 893 -594.8 1191.6 horTh1 NA 0.984 NA TRUE 2015 -316.6 635.2 horTh1 NA 0.473 NA TRUE 747 -584.2 1170.4 progrec 25 0.003 8.8 FALSE 515 -386.3 774.7 NA NA NA NA TRUE 1170 -189.1 380.1 NA NA NA NA TRUE logLikelihood (root)=-2647.8 logLikelihood (tree)=-2581.9 AIC (root)=5297.6 AIC (tree)=5177.6 Descripci´on del proceso llevado a cabo: 1. El modelo SurvCART realiza una partici´on recursiva de los datos bas´andose en las variables especificadas en gvars. Para cada posible divisi´on en los datos, el modelo eval´ua si la divisi´on mejora la homogeneidad dentro de los nodos resultantes. 2. Se utilizan las pruebas de inestabilidad de par´ametros para decidir si una variable de partici´on (y su punto de corte) proporciona una mejora estad´ısticamente significativa en la separaci´on de los datos. 3. Los par´ametros minsplit y minbucket aseguran que las divisiones resultantes no conduzcan a nodos muy peque˜nos, lo que podr´ıa llevar a un sobreajuste. 4. El par´ametro alpha establece el umbral de significancia para aceptar una divisi´on. 61 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM 5. Al final, el ´arbol de supervivencia resultante tiene nodos terminales que representan subgrupos de pacientes con diferentes distribuciones de tiempo hasta el evento, lo que permite una interpretaci´on detallada de los factores que influyen en la supervivencia. # Paso 4. Visualizaci´on del ´arbol de supervivencia. par(xpd =TRUE) plot(out, compress =TRUE) text(out, use.n =TRUE) Figura 4.1: ´ Arbol de Supervivencia generado por SurvCART. El ´arbol de supervivencia generado mediante el modelo SurvCART en el conjunto de datos GBSG2 muestra c´omo diferentes variables afectan el tiempo de supervivencia sin recurrencia en pacientes con c´ancer de mama. A continuaci´on, se realiza un an´alisis detallado de cada nodo del ´arbol y de las decisiones de partici´on: Nodo ra´ız (Nodo 1) Variable de partici´on: pnodes (n´umero de nodos positivos). Punto de corte: 10. N´umero de sujetos: 686 (N). N´umero de eventos: 150 (D). Este nodo ra´ız divide la muestra en dos grupos bas´andose en si el n´umero de nodos positivos es menor que 10. 62 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Sub´arbol izquierdo del nodo ra´ız (Nodo 2) Condici´on: pnodes <10. Variable de partici´on: pnodes. Punto de corte: 4. Este nodo se divide nuevamente en funci´on de si el n´umero de nodos positivos es menor que 4. Sub´arbol izquierdo del nodo 2 (Nodo 4) Condici´on: pnodes <4. Variable de partici´on: progrec (receptor de progesterona). Punto de corte: 90. Este nodo se divide en funci´on de los niveles del receptor de progesterona. Sub´arbol izquierdo del nodo 4 (Nodo 8) Condici´on: progrec <90. Estad´ısticas: Mediana de tiempo de supervivencia (q0.5): 1990 d´ıas. N´umero de sujetos (N): 239. N´umero de eventos (D): 91. Este nodo es terminal, lo que indica que los pacientes con menos de 4 nodos positivos y menos de 90 unidades del receptor de progesterona tienen una mediana de supervivencia de aproximadamente 1990 d´ıas. Sub´arbol derecho del nodo 4 (Nodo 9) Condici´on: progrec >= 90. Variable de partici´on: horTh1 (terapia hormonal). Punto de corte: 2. Este nodo se divide en funci´on de si los pacientes recibieron terapia hormonal. 63 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM # Paso 6. Visualizaci´on de la funci´on de riesgo H(t). KMPlot(out, xscale=365.25, type=2, ylab='Probabilidad de muerte', xlab='Tiempo (a~nos)') Figura 4.3: Funci´on de riesgo H(t) por nodos finales. No se aprecia bien cada grupo por separado al tener curvas de riesgo similares. Vamos a visualizar las gr´aficas individualmente, por separado. KMPlot(out, xscale=365.25, type=2, overlay=FALSE, mfrow=c(2,2), xlab="Tiempo (a~nos)", ylab="Probabilidad de censura") Figura 4.4: Funci´on H(t) para los nodos finales 8, 18, 19 y 10. 70 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Figura 4.5: Funci´on H(t) para los nodos finales 11, 6 y 7. Las gr´aficas de probabilidad de censura muestran la probabilidad de que los datos sean censurados (es decir, que el evento no haya ocurrido durante el per´ıodo de observaci´on) a lo largo del tiempo para cada uno de los nodos terminales del ´arbol de supervivencia. A continuaci´on se proporciona un an´alisis detallado de cada una de las 7 gr´aficas: Nodo 8: pnodes <10 & pnodes <4 & progrec <90 Condiciones: Pacientes con menos de 10 nodos positivos, menos de 4 nodos positivos y niveles del receptor de progesterona menores de 90. Curva de censura: La probabilidad de censura disminuye gradualmente a lo largo del tiempo. Al final del per´ıodo de observaci´on (7 a˜nos), la probabilidad de censura es aproximadamente entre 1-5 Interpretaci´on: Una baja probabilidad de censura sugiere que la mayor´ıa de los pacientes en este subgrupo experimentaron el evento durante el per´ıodo de seguimiento. Nodo 18: pnodes <10 & pnodes <4 & progrec >=90 & horTh1 <2 Condiciones: Pacientes con menos de 10 nodos positivos, menos de 4 nodos positivos, niveles del receptor de progesterona mayores o iguales a 90 y no recibieron terapia hormonal. 71 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Curva de censura: La probabilidad de censura disminuye de manera similar a la del Nodo 8, con una probabilidad de censura al final del per´ıodo de observaci´on de aproximadamente entre el 1-5 Interpretaci´on: Este subgrupo tiene una probabilidad de censura comparable a la del Nodo 8, aunque cae de manera m´as lenta, lo que sugiere que la ausencia de terapia hormonal puede afecta ligeramente la probabilidad de censura en este contexto. Nodo 19: pnodes <10 & pnodes <4 & progrec >=90 & horTh1 >=2 Condiciones: Pacientes con menos de 10 nodos positivos, menos de 4 nodos positivos, niveles del receptor de progesterona mayores o iguales a 90 y que recibieron terapia hormonal. Curva de censura: La probabilidad de censura disminuye gradualmente a lo largo del tiempo, con una probabilidad de censura de aproximadamente entre 1-5 Interpretaci´on: Este subgrupo tiene una probabilidad similar de censura en comparaci´on con los Nodos 8 y 18, lo que sugiere que la terapia hormonal podr´ıa estar no estar asociada a la probabilidad de ocurrencia de c´ancer. Nodo 10: pnodes <10 & pnodes >=4 & progrec <56 Condiciones: Pacientes con menos de 10 nodos positivos, 4 o m´as nodos positivos y niveles del receptor de progesterona menores de 56. Curva de censura: La probabilidad de censura disminuye gradualmente, con una probabilidad de censura de aproximadamente el 10 Interpretaci´on: Este subgrupo tiene una probabilidad de censura m´as alta que los anteriores nodos, lo que sugiere que una parte significativa de los pacientes experiment´o el evento (c´ancer) durante el per´ıodo de seguimiento. Nodo 11: pnodes <10 & pnodes >=4 & progrec >=56 Condiciones: Pacientes con menos de 10 nodos positivos, 4 o m´as nodos positivos y niveles del receptor de progesterona mayores o iguales a 56. Curva de censura: La probabilidad de censura disminuye gradualmente, con una probabilidad de censura de aproximadamente el 10 72 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Interpretaci´on: Este subgrupo tiene una probabilidad de censura comparable a la del Nodo 10, lo que sugiere que niveles m´as altos del receptor de progesterona no afectan significativamente la probabilidad de censura en pacientes con 4 o m´as nodos positivos. Nodo 6: pnodes >=10 & progrec <25 Condiciones: Pacientes con 10 o m´as nodos positivos y niveles del receptor de progesterona menores de 25. Curva de censura: La probabilidad de censura disminuye gradualmente, con una probabilidad de censura de aproximadamente entre 20-60 Interpretaci´on: Este subgrupo tiene una probabilidad de censura moderada, lo que sugiere que una parte significativa de los pacientes experiment´o el evento (c´ancer) durante el per´ıodo de seguimiento. Nodo 7: pnodes >=10 & progrec >=25 Condiciones: Pacientes con 10 o m´as nodospositivos y niveles del receptor de progesterona mayores o iguales a 25. Curva de censura: La probabilidad de censura disminuye gradualmente, con una probabilidad de censura de aproximadamente el 10 Interpretaci´on: Este subgrupo tiene la mayor probabilidad de censura entre todos los nodos analizados, lo que sugiere que una mayor proporci´on de pacientes no experiment´o el evento durante el per´ıodo de seguimiento. Este subgrupo tiene una probabilidad de censura similar a los nodos anteriores que tambi´en ten´ıan el 10observaci´on es menor, lo que sugiere que los pacientes sufrieron c´ancer en un periodo m´as corto de tiempo. Las gr´aficas de probabilidad de censura proporcionan una visualizaci´on clara de c´omo diferentes combinaciones de variables afectan la probabilidad de censura en pacientes con c´ancer de mama. Estas gr´aficas ayudan a identificar subgrupos de pacientes con diferentes perfiles de riesgo de censura, lo que puede guiar la interpretaci´on de los resultados de supervivencia y la planificaci´on de futuros estudios. # Paso 7. Validaci´on del modelo: validaci´on cruzada. #Validaci´on cruzada usando crossval 73 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM library(crossval) library(survivalROC) # Definir la f´ormula formula <- Surv(time, cens) ~horTh1 +age +menostat1 +tsize + tgrade1 +pnodes +progrec +estrec # Funci´on de validaci´on cruzada sin visualizaci´on de curvas ROC crossval_cox <- function(data, formula, k = 10) { set.seed(49084092)# Para reproducibilidad folds <- sample(rep(1:k, length =nrow(data))) results <- lapply(1:k, function(i) { train_data <- data[folds != i, ] test_data <- data[folds == i, ] fit <- coxph(formula, data =train_data) pred <- predict(fit, newdata =test_data, type ="risk") surv_roc <- survivalROC(Stime =test_data$time, status = test_data$cens, marker =pred, predict.time = median(test_data$time), method ="KM") return(surv_roc$AUC) }) list(mean_auc =mean(unlist(results))) } # Ejecutar validaci´on cruzada auc_mean <- crossval_cox(data =GBSG2, formula =formula, k = 10) print(auc_mean$mean_auc) [1] 0.7247845 El AUC promedio obtenido es aproximadamente 0.725. Un AUC de 0.7247845 indica que el modelo tiene una buena capacidad para distinguir entre los sujetos que experimentan el evento (p.ej., recurrencia del c´ancer) y aquellos que no lo experimentan. En general, un AUC entre 0.7 y 0.8 se considera aceptable, mientras que un AUC superior a 0.8 se considera excelente. Por lo tanto, este resultado sugiere que el modelo de Cox ajustado proporciona predicciones razonablemente precisas sobre la supervivencia de los pacientes con c´ancer de mama en el conjunto de datos GBSG2. # Paso 8. Predicci´on. pred <- predict.SurvCART(object=out, newdata=GBSG2) head(pred) 74 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM horTh age menostat tsize tgrade pnodes progrec estrec time cens horTh1 no 70 Post 21 II 3 48 66 1814 1 1 yes 56 Post 12 II 7 61 77 2018 1 2 yes 58 Post 35 II 9 52 271 712 1 2 yes 59 Post 17 II 4 60 29 1807 1 2 no 73 Post 35 II 1 26 65 772 1 1 no 32 Pre 57 III 24 0 13 448 1 1 tgrade1 menostat1 subjid node median.T median.C Q1.T Q3.T Q1.C Q3.C 2 2 1 8 1990 1604 918 NA 1091 1976 2 2 2 11 2015 1702 982 NA 1095 2010 2 2 3 10 893 1637 417 2456 974 1977 2 2 4 11 2015 1702 982 NA 1095 2010 2 2 5 8 1990 1604 918 NA 1091 1976 3 1 6 6 515 1826 359 855 973 1856 A modo de ejemplo, supongamos que nos llegan tres nuevos pacientes, vamos a predecir tras medirles las variables bajo estudio. nuevos_pacientes <- data.frame( horTh1 =c(1,2,2), # No, S´ı, S´ı age =c(45,60,50), # Edad: 45, 60, 50 menostat1 =c(1,2,2), # Premenop´ausica, Postmenop´ausica, Postmenop´ausica tsize =c(20,30,15), # Tama~no del tumor: 20, 30, 15 mm tgrade1 =c(1,3,2), # Grado del tumor: I, III, II pnodes =c(0,4,10), # Ganglios positivos: 0, 4, 10 progrec =c(50,150,80), # Receptor de progesterona: 50, 150, 80 estrec =c(200,50,75)# Receptor de estr´ogeno: 200, 50, 75 ) # Usar el modelo SurvCART para predecir en qu´e nodo se encasillar´ıan estos pacientes pred_nuevos <- predict.SurvCART(object =out, newdata = nuevos_pacientes) # Mostrar las predicciones id_pred = 1:nrow(pred_nuevos) pred_nuevos =cbind(id_pred,pred_nuevos) 75 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM id pred horTh1 age menostat1 tsize tgrade1 pnodes progrec 1 1 45 1 20 1 0 50 2 2 60 2 30 3 4 150 3 2 50 2 15 2 10 80 estrec node median.T median.C Q1.T Q3.T Q1.C Q3.C 200 8 1990 1604 918 NA 1091 1976 50 11 2015 1702 982 NA 1095 2010 75 7 1170 1505 730 1975 1062 1834 Interpretaci´on de los resultados Paciente 1: Este paciente, que no recibe terapia hormonal, tiene un tumor peque˜no (20 mm) y no tiene nodos positivos, se clasifica en el nodo 8. La mediana de supervivencia para pacientes en este nodo es alta (1990 d´ıas), lo que sugiere un pron´ostico favorable. Sin embargo, hay censura significativa (1604 d´ıas), lo que indica que muchos pacientes no tuvieron el evento durante el seguimiento. Paciente 2: Este paciente, que recibe terapia hormonal, es postmenop´ausica, tiene un tumor m´as grande (30 mm) y varios ganglios positivos (4), se clasifica en el nodo 11. La mediana de supervivencia para este nodo tambi´en es alta (2015 d´ıas), lo que sugiere un buen pron´ostico, pero nuevamente con una mediana de censura alta (1702 d´ıas). Paciente 3: Este paciente, que tambi´en recibe terapia hormonal, es postmenop´ausica, tiene un tumor m´as peque˜no (15 mm) pero muchos ganglioss positivos (10), se clasifica en el nodo 7. La mediana de supervivencia es menor (1170 d´ıas) comparada con los otros nodos, indicando un pron´ostico menos favorable. La censura es alta (1505 d´ıas), pero a´un muestra una mayor ocurrencia de eventos comparado con los otros dos nodos. Conclusi´on El paciente 1 tiene un pron´ostico m´as favorable con una mediana de supervivencia de aproximadamente 5.5 a˜nos, y una alta proporci´on de censura. El paciente 2 tambi´en tiene un pron´ostico favorable con una mediana de supervivencia de aproximadamente 5.5 a˜nos, pero con caracter´ısticas de mayor riesgo (edad avanzada, tumor m´as grande y mayor grado). El paciente 3 tiene un pron´ostico menos favorable con una mediana de supervivencia de aproximadamente 3.2 a˜nos, principalmente debido a la alta cantidad de nodos positivos. 76 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM # Paso 9. Prueba de estabilidad de par´ametros para la variable de partici´on categ´orica. stab_cat <- StabCat.surv(data=GBSG2, timevar="time", censorvar="cens", splitvar="horTh1", event.ind=1) print(stab_cat$pval) [1] 0.003431809 Puesto que el p valor es menor que 0.05, hay evidencia de heterogeneidad en los par´ametros del tiempo hasta el evento entre los pacientes que reciben y no reciben terapia hormonal. Esto sugiere que la terapia hormonal tiene un impacto significativo en la distribuci´on del tiempo hasta el evento y puede ser un factor importante a considerar en el an´alisis de supervivencia. # Paso 10. Prueba de estabilidad de par´ametros para las variables de partici´on continua. # Prueba de estabilidad de par´ametros para la variable continua 'progrec' stab_cont_progrec <- StabCont.surv(data =GBSG2, timevar ="time", censorvar ="cens", splitvar ="progrec", event.ind = 1) print(stab_cont_progrec$pval) [1] 2.688536e-09 Hemos obtenido un p-valor muy bajo: Heterogeneidad: Hay evidencia de heterogeneidad en los par´ametros del tiempo hasta el evento entre los diferentes niveles de progrec. Impacto Cl´ınico: Esto sugiere que los niveles de receptor de progesterona (progrec) tienen un impacto significativo en la distribuci´on del tiempo hasta el evento y podr´ıan justificar una partici´on en el modelo de ´arbol. # Prueba de estabilidad de par´ametros para la variable continua 'pnodes' stab_cont_pnodes <- StabCont.surv(data =GBSG2, timevar ="time", 77 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM censorvar ="cens", splitvar ="pnodes", event.ind = 1) print(stab_cont_pnodes$pval) [1] 2.172706e-13 Hemos obtenido un p-valor muy bajo: Heterogeneidad: Hay evidencia de heterogeneidad en los par´ametros del tiempo hasta el evento entre los diferentes niveles de pnodes. Impacto Cl´ınico: Esto sugiere que el n´umero de nodos positivos (pnodes) tiene un impacto significativo en la distribuci´on del tiempo hasta el evento y podr´ıan justificar una partici´on en el modelo de ´arbol. # Paso 11. Evaluar la probabilidad de ´exito cl´ınico. Supongamos que queremos evaluar la probabilidad de ´exito cl´ınico de un ensayo sobre la supervivencia sin recurrencia en pacientes con c´ancer de mama. n_treatment <- sum(GBSG2$horTh1 == 2)# N´umero de sujetos en el grupo de tratamiento n_control <- sum(GBSG2$horTh1 == 1)# N´umero de sujetos en el grupo de control a<- n_treatment /n_control # hr.prior y D.prior se obtiene de estudios previos o meta-an´alisis. Por este motivo, asumiremos valores hipot´eticos para estos par´ametros. pos_result <- PoS(type ="surv", nsamples = 2,# Dos grupos de pacientes: tratamiento / sin control null.value = 1,# H_0: No hay diferencia en el hazard ratio entre los grups (HR = 1) alternative ="less",# Estamos interesados en saber si el hazard ratio es menor que 1, # lo que sugerir´ıa que el tratamiento reduce el riesgo c´ancer en comparaci´on con el contrl N= 686,# Tama~no de la muestra D=sum(GBSG2$cens == 1), # N´umero total de eventos (personas con c´ancer) al final del an´alisis. a, succ.crit ="clinical", 78 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM clin.succ.threshold = 0.8,# Umbral cl´ınicamente significatio para el HR. Un HR menor que 0.8 # indicar´ıa una reducci´on del riesgo cl´ınicamente significativa hr.prior = 0.75,#Valor medio de la distribuci´on previa del hazard ratio. D.prior = 133)#N´umero de eventos en la distribuci´on previa. Based on the prior distribution the probability of success is calculated as 0.62 Figura 4.6: Distribuci´on predictiva del cociente de riesgo dada la distribuci´on inicial. Interpretaci´on del Resultado La probabilidad de ´exito cl´ınico obtenida es 0.62 o 62 %. Esto significa que, basado en la distribuci´on previa del hazard ratio y los par´ametros especificados, hay un 62 % de probabilidad de que el ensayo cl´ınico muestre un HR menor que 0.8, indicando una reducci´on del riesgo cl´ınicamente significativa. Un HR menor que 0.8 sugiere una reducci´on significativa en el riesgo de recurrencia del c´ancer debido al tratamiento. Por lo tanto, un 62 % de probabilidad de ´exito cl´ınico indica que el tratamiento tiene una buena probabilidad de ser eficaz. 79 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM plot.survival(v.all) Figura 4.9: Supervivencia OOB, Puntuaci´on de Brier OOB, Puntuaci´on CRPS OOB, y Mortalidad OOB vs. Tiempo . La figura 4.9 muestra cuatro gr´aficos que eval´uan el rendimiento del modelo de bosque aleatorio de supervivencia utilizando el conjunto de datos veteran. A continuaci´on se presenta un an´alisis conciso de cada gr´afico: Supervivencia OOB (Out-of-Bag): •Descripci´on: Este gr´afico muestra las curvas de supervivencia estimadas para cada individuo en el conjunto de datos OOB. •Interpretaci´on: La tendencia descendente de las curvas indica una disminuci´on de la probabilidad de supervivencia con el tiempo, lo cual es coherente con la naturaleza de los datos de supervivencia. La curva roja promedio proporciona una visi´on general del comportamiento de la supervivencia en la poblaci´on estudiada. 86 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Puntuaci´on de Brier OOB: La puntuaci´on Brier (1950) es una funci´on que mide la precisi´on de las predicciones probabil´ıstica, equivalente al error cuadr´atico medio aplicado a las probabilidades ajustadas. La puntuaci´on de de Brier es aplicable a tareas en las que las predicciones deben asignar probabilidades a un conjunto de clases o resultados discretos mutuamente excluyentes. •Descripci´on: Representa la puntuaci´on de Brier, una medida de la precisi´on de las predicciones probabil´ısticas, a lo largo del tiempo. •Interpretaci´on: Valores m´as bajos de la puntuaci´on de Brier indican mejores predicciones. Las curvas muestran una mayor preci- si´on en las predicciones durante los primeros d´ıas, con una disminuci´on gradual en la precisi´on a medida que avanza el tiempo. Las etiquetas azules indican los cuartiles de la puntuaci´on de Brier. Puntuaci´on CRPS OOB: •Descripci´on: Muestra la puntuaci´on CRPS (Continuous Ranked Probability Score) para diferentes grupos de pacientes a lo largo del tiempo. •Interpretaci´on: Al igual que la puntuaci´on de Brier, valores m´as bajos de CRPS indican una mejor calidad de las predicciones. Este gr´afico proporciona una visi´on adicional de la precisi´on del modelo, mostrando c´omo el CRPS var´ıa con el tiempo. Mortalidad OOB vs. Tiempo: •Descripci´on: Este gr´afico presenta la mortalidad estimada para cada individuo en funci´on del tiempo. •Interpretaci´on: La dispersi´on de los puntos azules indica la variabilidad en las predicciones de mortalidad del modelo. La l´ınea punteada muestra una tendencia general de disminuci´on de la mortalidad a medida que aumenta el tiempo. Podr´ıamos preguntarnos qu´e n´umero de ´arboles es adecuado tomar para construir el bosque. Para ello tenemos una gr´afica que nos proporciona el error OOB en funci´on del n´umero de ´arboles, desde 1 hasta ntree. 87 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Vamos a dibujar el gr´afico con la siguiente orden: plot(rfsrc(Surv(time, status) ~., data =veteran, ntree = 1000, block.size = 1) Figura 4.10: Evoluci´on del error OOB en funci´on del n´umero de ´arboles. Podemos visualizar las gr´aficas parciales para las variables predictoras, tanto para la mortalidad como para la probabilidad de supervivencia: plot.variable(v.all,partial=T) plot.variable(v.all, surv.type ="surv", time = 200,partial =TRUE, smooth.lines =TRUE) 88 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Figura 4.11: Gr´aficas parciales para la mortalidad seg´un las variables predictoras. Figura 4.12: Gr´aficas parciales para la probabilidad de supervivencia seg´un las variables predictoras. 89 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Podemos ahora estudiar la importancia de las variables mediante permutaci´on. # Establecemos una semilla set.seed(49084092) print(100*vimp(v.all, importance ="permute")$importance) Variable Importancia (*100) trt -0.0897 celltype 5.1099 karno 12.6002 diagtime 0.0386 age -0.0903 prior -0.0340 Tabla 4.7: Importancia de las variables en veteran. Vamos a hacer un ajuste de cada valor para que la suma total sea 100, manteniendo las proporciones y poder visualizarlo en un gr´afico de barras de manera m´as intuitiva. # Datos de importancia originales importance <- c(-0.07234318,5.01158028,13.03408674,0.07073978, -0.14105759,-0.02370235) # Sumar los valores de importancia total_importance <- sum(importance) # Ajustar los valores para que sumen 100, manteniendo la proporci´on scaled_importance <- 100 * importance /total_importance # Verificar que los valores ajustados sumen 100 sum(scaled_importance) scaled_importance variables <- c("trt","celltype","karno","diagtime","age", "prior") # Crear el gr´afico de barras barplot(scaled_importance, names.arg =variables, col ='skyblue', xlab ='Variables', ylab ='Importancia (%)', main ='Importancia de las Variables en el Modelo de Bosque de Supervivencia') 90 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM # Agregar l´ınea horizontal en y = 0 para referencia abline(h = 0, col ="red", lwd = 2, lty = 2) Figura 4.13: Gr´afico de barras de la importancia de las variables predictoras en veteran. Se observa claramente que las variables m´as relevantes son karno ycelltype. Ya que la puntuaci´on de Karnofsky es significativamente la m´as relevante, podemos observar las curvas de supervivencia por grupos para esta variable: veteran$karno_group <- cut(veteran$karno, breaks = quantile(veteran$karno, probs =seq(0,1,0.25)), include.lowest = TRUE) surv_fit <- survfit(Surv(time, status) ~karno_group, data = veteran) ggsurvplot(surv_fit, data =veteran, p.val =F, conf.int =F, title ="Curvas de Supervivencia por Grupos de Karnofsky", legend.title ="Grupo de Karnofsky", xlab ="Tiempo (d´ıas)", ylab ="Probabilidad de Supervivencia") 91 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM Figura 4.14: Curvas de supervivencia para diferentes subgrupos de la variable karno. Ahora consideramos el efecto incremental de cada predictor utilizando un an´alisis anidado. Ordenamos los predictores por sus valores de importancia y consideramos la secuencia anidada de modelos comenzando con la variable m´as importante, seguido del modelo con las dos variables principales, luego el modelo con las tres variables principales, y as´ı sucesivamente: imp <- vimp(v.all)$importance pnames <- vimp(v.all)$xvar.names pnames.order <- pnames[rev(order(imp))] n.pred <- length(pnames) vet.err <- rep(0, n.pred) for (k in 1:n.pred){ rsf.f <- "Surv(time,status)~" rsf.f <- as.formula(paste(rsf.f, paste(pnames.order[1:k],collapse="+" ⌋ ))) vet.err[k] <- rfsrc(rsf.f, veteran, ntree=ntree, splitrule="logrank")$err.rate[ntree] } vet.imp.out <- as.data.frame( cbind(round(rev(sort(imp)),4), round(vet.err,4), 92 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM round(-diff(c(0.5,vet.err)),4)), row.names=pnames.order) colnames(vet.imp.out) <- c("Imp","Err","Drop Err") print(vet.imp.out) Variable Importancia (*100) Error Disminuci´on del Error karno 12.6055 0.3580 0.1420 celltype 5.0805 0.3093 0.0486 age 0.0893 0.3038 0.0055 diagtime -0.0079 0.2952 0.0086 trt -0.0224 0.2920 0.0032 prior -0.0345 0.2934 -0.0014 Tabla 4.8: Importancia de las Variables, Error y Disminuci´on del Error. La primera columna es el valor de importancia de un predictor en el modelo completo. El valor k-´esimo en la segunda columna es la tasa de error para el k-´esimo modelo anidado, mientras que el valor k-´esimo en la tercera columna es la diferencia entre la tasa de error para el k-´esimo y el (k−1)- ´esimo modelo anidado, donde la tasa de error para el modelo nulo, k= 0, es 0.5. Se puede observar que no se gana mucho al usar m´as de 3 predictores y que los 2 predictores principales explican gran parte del poder predictivo. Por ´ultimo, utilizamos el modelo entrenado para realizar predicciones en el conjunto de datos y evaluamos la adecuaci´on de las predicciones utilizando el ´ındice C, para evaluar qu´e tan bueno es el modelo. # Realizar predicciones en el conjunto de datos completo (modo restauraci´on) v.pred <- predict(v.all, veteran) print(v.pred) Sample size of test (predict) data: 137 Number of grow trees: 1000 Average no. of grow terminal nodes: 6.136 Total no. of grow variables: 6 Resampling used to grow trees: swor Resample size used to grow trees: 87 Analysis: RSF Family: surv CRPS: 0.05372582 Requested performance error: 0.23809793 93 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM El ´ındice C es de 0.2381. Este valor indica que el modelo tiene una tasa de error del 23.81 % en sus predicciones. A continuaci´on, se presentan algunas observaciones y conclusiones sobre este valor: Un ´ındice C de 0.2381 sugiere que el modelo tiene un rendimiento aceptablemente bueno en la predicci´on de los tiempos de supervivencia. Este valor indica que el modelo ordena correctamente los tiempos de supervivencia en aproximadamente el 76.19 % de los pares comparables. Comparado con un modelo que no mejorar´ıa m´as all´a del azar (´ındice C = 0.5), este modelo muestra una mejora significativa en la capacidad predictiva. La tasa de error del 23.81 % es relativamente baja, lo que sugiere que el modelo es robusto y realiza predicciones precisas sobre la supervivencia de los pacientes en el conjunto de datos veteran.La m´etrica CRPS complementa esta interpretaci´on, indicando un buen ajuste del modelo con un valor de 0.0537, que tambi´en sugiere alta precisi´on en las predicciones. Tambi´en se podr´ıa realizar validaci´on cruzada para estudiar la calibraci´on del modelo. La calibraci´on eval´ua la relaci´on entre las probabilidades predichas por el modelo y las proporciones observadas de eventos en los datos. Si un modelo predice que la probabilidad de un evento es del 80 % para un conjunto de casos, y efectivamente el 80 % de esos casos experimentan el evento, el modelo est´a bien calibrado. Para evaluar la calibraci´on del modelo, utilizamos la funci´on pec de la librer´ıa pec desarrollada por Gerds (2023) que genera la curva de calibraci´on: # Calibraci´on del Modelo usando pec library(pec) cal_model <- pec(object =list("RSF" =v.all), formula =Surv(time, status) ~karno +trt + celltype +diagtime +age +prior, data =veteran, times =seq(0,max(veteran$time), length.out = 100), exact =TRUE, splitMethod ="BootCv", B= 1000, cens.model ="none") plot(cal_model, xlim =c(0,max(veteran$time)), ylim =c(0,1)) La funci´on pec nos permite evaluar la precisi´on del modelo a trav´es de la puntuaci´on de Brier y generar una curva de calibraci´on que muestra c´omo se ajustan las predicciones a los datos observados. La gr´afica resultante nos 94 CAP´ ITULO 4. ILUSTRACIONES CON R PROGRAM da una idea visual de la calibraci´on del modelo, complementando el an´alisis realizado con el ´ındice C. Figura 4.15: Gr´afica de calibraci´on. Ejes: •Eje X (Tiempo): Representa el tiempo en d´ıas. •Eje Y (Error de Predicci´on): Representa el error de predicci´on, con valores que var´ıan de 0 a 1. Curvas: •L´ınea Negra (Kaplan-Meier): Representa el error de predicci´on del modelo de referencia Kaplan-Meier. •L´ınea Roja (RSF): Representa el error de predicci´on del modelo de bosque aleatorio de supervivencia (RSF). Error Inicial: Ambos modelos comienzan con un error de predicci´on bajo, que r´apidamente aumenta. Esto es com´un en modelos de supervivencia debido a que al principio del seguimiento hay menos eventos y la predicci´on es m´as f´acil. 95