Técnicas no paramétricas y modelos de regresión para datos de tiempo de vida

UNIVERSIDAD DE SEVILLA -CURSO 2017/2018 -GRADO EN MATEM´ ATICAS TRABAJO FIN DE GRADO T´ecnicas no param´etricas y modelos de regresi´on para datos de tiempo de vida Departamento : Estad´ıstica e Investigaci´on Operativa Realizado por : Lorena Barrenechea L´opez Tutora : Inmaculada Barranco Chamorro 2 Resumen La caracter´ıstica principal en un estudio de supervivencia es que los sujetos bajo estudio se observan durante un tiempo estipulado, denominado tiempo de seguimiento, y el objetivo es conocer el tiempo en el que tiene lugar el evento o suceso de inter´es (tiempo de ocurrencia). Dado que este an´alisis puede aplicarse a distintos campos de investigaci´on, el evento de estudio ser´a el propio de cada disciplina: en medicina puede ser la muerte de un paciente que sufre una patolog´ıa concreta o bien un episodio relacionado con la enfermedad; en ingenier´ıa puede referirse al fallo de una m´aquina o de una de sus piezas y en econom´ıa podr´ıa ser el inicio de un nuevo empleo despu´es de un periodo de desempleo. A lo largo de este trabajo, explicaremos algunas de las t´ecnicas utilizadas para el estudio de supervivencia. En el Cap´ıtulo 1, estableceremos algunas definiciones para poder entender las funciones usadas en cada m´etodo, y que desarrollaremos en los siguientes cap´ıtulos. Hablaremos sobre la censura y el truncamiento, que es una manera de poder continuar con el estudio, cuando un individuo (o una pieza) muere (o falla). Y por ´ultimo, las t´ecnicas param´etricas que se utilizan. En el Cap´ıtulo 2, presentaremos las t´ecnicas no param´etricas que se utilizan cuando los datos no se ajustan a ninguna distribuci´on conocida. Estas t´ecnicas son Kaplan-Meier, Nelson-Aalen y el test Log-Rank. En el Cap´ıtulo 3, nos centramos en un modelo en el que adem´as de relacionar la tasa de supervivencia con el tiempo, se a˜naden diferentes covariables explicativas. Es el denominado modelo de regresi´on de Cox. Vemos los residuos que se utilizan para comprobar que ese modelo es v´alido. Por ´ultimo, introducimos lo que se llama “modelo de Cox estratificado” que surge cuando la hip´otesis de riesgos proporcionales no se cumple. Finalmente, en el Cap´ıtulo 4, utilizamos el software R para aplicar lo expuesto en este trabajo a un conjunto de datos reales. 3 4 Abstract The main feature in a survival study is that, the subjects under study are observed during a fixed period of time, called follow-up time, and the aim is to determine the time in which the event of interest takes place (time of occurrence). Since this analysis can be applied to different fields of research, we will adapt the study to each discipline: in the case of medicine, it can be the death of a patient suffering from a specific pathology, or an episode related to the disease; in engineering it can refer to the failure of a machine or of one of its parts, and in economy it could be the start of a new job after a period of unemployment. Throughout this essay, we will explain some of the techniques used for the survival study. In Chapter 1, we will establish some definitions to understand the functions used in each method, and we will develop them in the next chapters. We will talk about censoring and truncation, which are two methods that allow us to continue with the study. Both happen when an individual (or an item) dies or has a failure. And finally, we will focus on the parametric techniques that are used. In Chapter 2, we will show the non-parametric techniques required when the data doesn’t fit to any of the known distributions.These techniques are KaplanMeier, Nelson-Aalen and the Log-Rank test. In Chapter 3, we will focus on a model which, apart from associating the survival rate and time, adds different explanatory covariates. This is called Cox regression model. This chapter also analyses the residues that are used to validate the model. Later, we will introduce what is called “stratified Cox model”, which arises when the proportional risks hypothesis does not meet the conditions required. Finally, in Chapter 4, we will use R software to apply what has been previously mentioned in this essay to a set of real data. 5 6 ´ Indice general 1. Generalidades y Modelos 11 1.1. Funcionesb´asicas.......................... 11 1.2. Censura y Truncamiento . . . . . . . . . . . . . . . . . . . . . . 17 1.2.1. Censura ........................... 17 1.2.2. Truncamiento . . . . . . . . . . . . . . . . . . . . . . . . 19 1.3. T´ecnicas param´etricas . . . . . . . . . . . . . . . . . . . . . . . 19 1.3.1. Distribuci´on Exponencial . . . . . . . . . . . . . . . . . . 20 1.3.2. Distribuci´on Weibull . . . . . . . . . . . . . . . . . . . . 21 1.3.3. Distribuci´on Log-normal . . . . . . . . . . . . . . . . . . 22 2. T´ecnicas no param´etricas 25 2.1. Introducci´on............................. 25 2.2. M´etodo Kaplan-Meier . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.1. Varianza del estimador de Kaplan-Meier ˆ S(t) ...... 34 2.3. M´etodo Nelson-Aalen . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.1. Varianza del estimador de Nelson-Aalen H........ 37 2.3.2. Estimador de Nelson-Aalen, ˆ H(t), a partir del estimador de Kaplan-Meier, ˆ S(t) ................... 37 2.4. TestLog-Rank ........................... 39 3. Modelo de Regresi´on de Cox 45 3.1. Formulaci´on del modelo . . . . . . . . . . . . . . . . . . . . . . 46 3.2. Hip´otesis de riesgos proporcionales, PH . . . . . . . . . . . . . . 48 3.2.1. Diagrama de diagn´ostico para PH . . . . . . . . . . . . . 51 3.3. Estimaci´on de los coeficientes . . . . . . . . . . . . . . . . . . . 52 3.3.1. Contrastes de hip´otesis . . . . . . . . . . . . . . . . . . . 56 3.4. Residuos en el an´alisis de supervivencia . . . . . . . . . . . . . . 57 3.4.1. Residuos de Cox-Snell . . . . . . . . . . . . . . . . . . . 59 3.4.2. Residuos de martingala . . . . . . . . . . . . . . . . . . . 59 3.4.3. Residuos basados en el estad´ıstico Deviance . . . . . . . 60 3.4.4. Residuos de Schoenfeld . . . . . . . . . . . . . . . . . . . 60 3.4.5. Residuos escalados de de Schoenfeld . . . . . . . . . . . . 61 3.4.6. Residuos dfbeta ....................... 61 7 8´ INDICE GENERAL 3.5. Modelo estratificado . . . . . . . . . . . . . . . . . . . . . . . . 61 4. Aplicaci´on pr´actica con el software R 63 4.1. Conclusiones al estudio de los datos de metadona ........ 82 Referencias 86 ´ Indice de figuras 1.1. Funci´on de riesgo creciente. . . . . . . . . . . . . . . . . . . . . 14 1.2. Funci´on de riesgo decreciente. . . . . . . . . . . . . . . . . . . . 15 1.3. Funci´on de riesgo constante. . . . . . . . . . . . . . . . . . . . . 15 1.4. Funci´on de riesgo con forma de ba˜nera. . . . . . . . . . . . . . . 16 1.5. Representaci´on de una Exponencial. . . . . . . . . . . . . . . . . 21 1.6. Representaci´on de una Weibull. . . . . . . . . . . . . . . . . . . 22 1.7. Representaci´on de una Log-normal. . . . . . . . . . . . . . . . . 24 2.1. Representaci´on de la funci´on emp´ırica de una muestra de valores. 27 2.2. Estimaci´on utilizando el m´etodo Kaplan-Meier con datos de metadona................................. 28 2.3. Estudio de 6 unidades. . . . . . . . . . . . . . . . . . . . . . . . 30 2.4. M´etodo de Kaplan-Meier estimando S(t). . . . . . . . . . . . . . 33 2.5. Informaci´on obtenida con el programa R, con los datos de n´umero de millones de revoluciones de bolas de cer´amica. . . . . . . . 35 2.6. M´etodo de Kaplan-Meier estimando S(t). . . . . . . . . . . . . . 36 2.7. Estimaci´on utilizando el m´etodo Nelson-Aalen. . . . . . . . . . . 39 2.8. Estimaci´on utilizando el m´etodo Kaplan-Meier para dos grupos. 43 3.1. Gr´afico en el que se representan las curvas de riesgo para el grupo sin cirug´ıa y el grupo con cirug´ıa. . . . . . . . . . . . . . 50 3.2. Riesgos proporcionales del rodamiento de bolas cer´amicas. . . . 52 4.1. Estimador de Kaplan-Meier. . . . . . . . . . . . . . . . . . . . . 68 4.2. Estimador de Nelson-Aalen. . . . . . . . . . . . . . . . . . . . . 70 4.3. Modelo con clinica como estrato. . . . . . . . . . . . . . . . . . 74 4.4. Salida de los residuos de Cox-Snell. . . . . . . . . . . . . . . . . 78 4.5. Salida de los residuos de martingala para prision. . . . . . . . . 79 4.6. Salida de los residuos del estad´ıstico Deviance. . . . . . . . . . . 80 4.7. Salida de los residuos escalados de Schoenfeld. . . . . . . . . . . 81 4.8. Salida de los residuos dfbeta. . . . . . . . . . . . . . . . . . . . . 82 9 16 CAP´ ITULO 1. GENERALIDADES Y MODELOS Figura 1.4: Funci´on de riesgo con forma de ba˜nera. Definici´on 1.1.4 Se define la funci´on de riesgo acumulada (cumulative hazard rate), como H(t) = Zt 0 h(x)dx, t > 0 Nota 1.1.1 Al representar H(t), nos aparecer´a una l´ınea recta si h(t) es constante, una funci´on que crece m´as r´apido que una recta si h(t) es creciente, y una que crece m´as despacio si h(t) es decreciente. Proposici´on 1.1.3 Se tienen las siguientes relaciones con respecto a la funci´on de riesgo y la de supervivencia: 1. S0(t) = −f(t) 2. h(t) = −d dtlogS(t) 3. S(t) = e−H(t) Demostraci´on 1.1.3 1. S(t)=1−F(t) =⇒S0(t) = −F0(t)(1.1) =−f(t). 2. d dt{logS(t)}=S0(t) S(t)=−f(t) S(t) (1.6) =−h(t) =⇒h(t) = −d dtlogS(t). 3. Del anterior punto, si integramos de 0 a xrespecto de t, tendremos Zx 0 d dtlogS(t)dt =logS(t)it=x t=00 = logS(x)−logS(0) = logS(x)−log1 = =logS(x) = −Zx 0 h(t)dt =⇒logS(x) = −Zx 0 h(t)dt =−H(x) =⇒ =⇒S(x) = e−H(x).  1.2. CENSURA Y TRUNCAMIENTO 17 Definici´on 1.1.5 Se denomina cuantil de orden p,tp, al valor que cumple P(T≤tp) = p, 0<p<1. Es decir, tp=F−1(p). Nota 1.1.2 Si hablamos en t´erminos de percentiles,tpes el percentil 100p. Nota 1.1.3 El percentil 50, o cuantil de orden 0.50, t0.50 es la mediana. 1.2. Censura y Truncamiento 1.2.1. Censura Se dice que una observaci´on en un estudio,(Lawless 2003, [15]), est´a censurada cuando el individuo no presenta el evento de inter´es durante el tiempo de seguimiento (duraci´on del estudio). Adem´as, como es l´ogico, el estudio debe tener un punto de inicio y final adaptado a los recursos disponibles. Seg´un si el evento ha tenido lugar antes del punto de inicio del estudio o bien, no se ha ocasionado una vez dado por finalizado, el tiempo de censura ser´a la fecha de inicio o final del estudio. Dependiendo del ´ambito de estudio, el dise˜no experimental y los resultados obtenidos para la variable, las observaciones censuradas se pueden agrupar en 3 grandes grupos: censuradas por la derecha, censuradas por la izquierda y censuradas en un intervalo. En los siguientes subapartados se explica cada uno de estos grupos de forma detallada. 1.Variables censuradas por la derecha Una observaci´on est´a censurada por la derecha cuando el evento de inter´es no tiene lugar durante el per´ıodo de observaci´on, de modo que no es posible determinar el tiempo de ocurrencia. En funci´on de la causa que ha originado la censura por la derecha, podemos encontrar distintos tipos: Censura de tipo I La finalizaci´on del per´ıodo de observaci´on de los individuos tiene lugar a un tiempo predeterminado, el cual se fija durante el dise˜no del estudio; as´ı el tiempo de censura es conocido y fijo. La censura de tipo I puede ser: 18 CAP´ ITULO 1. GENERALIDADES Y MODELOS - FIJA: el tiempo de inicio y final del estudio (censura) es el mismo para todos los individuos. - PROGRESIVA: los individuos se dividen en grupos y cada uno de los grupos tiene un tiempo de censura concreto. - GENERALIZADA: cada uno de los individuos tiene un tiempo de entrada al estudio y un tiempo de censura espec´ıfico. Censura de tipo II La finalizaci´on del per´ıodo de observaci´on de los sujetos no ocurre en un tiempo prefijado, sino que ´este contin´ua hasta que ocurre el suceso de estudio en una proporci´on establecida de individuos respecto al total. Por tanto, el tiempo de censura no se conoce a priori, sino que se trata de una variable aleatoria, dado que la proporci´on que se estipula durante el dise˜no del estudio es la proporci´on de fallos. Censura de tipo III Este grupo est´a formado por la censura aleatoria, tambi´en llamada no informativa. Se denomina as´ı, porque el tiempo de censura lo determina un fen´omeno aleatorio no esperado, que tiene lugar durante la consecuci´on del estudio, e impide seguir con la observaci´on del individuo hasta el tiempo final. Estos sucesos no esperados pueden ser: la p´erdida del sujeto sin m´as informaci´on, el abandono voluntario del estudio o la experimentaci´on de un evento de competencia con el evento de inter´es, que obliga a eliminar al individuo del estudio. Por tanto, en todos estos casos, el tiempo de censura es aleatorio y tiene lugar antes del tiempo de finalizaci´on del estudio que se ha estipulado. 2.Variables censuradas por la izquierda Las observaciones censuradas por la izquierda son aquellas en las que el evento de inter´es ha tenido lugar antes del punto de inicio del estudio. As´ı, el tiempo de censura en este caso, ser´a el tiempo de inicio del per´ıodo de seguimiento ya que se conoce que el suceso ha ocurrido previamente, pero no puede saberse con exactitud cu´ando (no es cuantificable). Tal y como se ha indicado en la censura por la derecha, la censura por la izquierda tambi´en la encontramos en mediciones anal´ıticas de la variable. En este caso, el valor obtenido en la medici´on es inferior a un umbral determinado y, por este motivo no es cuantificable. Un ´ambito de estudio donde este tipo de censura es recurrente, es el orientado a la investigaci´on medioambiental, dado que los instrumentos de medida tienen un l´ımite de detecci´on espec´ıfico y a menudo se detectan observaciones que no lo alcanzan. 1.3. T´ ECNICAS PARAM ´ ETRICAS 19 3.Variables censuradas en un intervalo Las observaciones que presentan censura en un intervalo, suelen ocasionarse en aquellos estudios donde las mediciones de la variable de inter´es se realizan de forma peri´odica, de modo que es posible que el suceso de estudio haya tenido lugar en un tiempo entre dos de las mediciones. En este caso, se sabe que el tiempo de ocurrencia se sit´ua entre dos tiempos de censura (el m´aximo y el m´ınimo valor que conforman el intervalo), pero se desconoce el tiempo exacto y esto no permite su cuantificaci´on. Est´a documentado, que uno de los tipos de estudios donde este fen´omeno es m´as frecuente, son los estudios de vida ´util de alimentos, ya que el producto puede deteriorarse entre dos tiempos consecutivos de evaluaci´on. Tambi´en en problemas dentales, ocurren en un instante indeterminado, entre dos revisiones consecutivas. Aunque es cierto, que en ellos tambi´en podr´ıamos encontrar observaciones censuradas por la derecha (por ejemplo, si un individuo acepta el producto a´un superando el tiempo m´aximo de almacenaje) o por la izquierda (por ejemplo, si se testan productos de la competencia y se desconoce la fecha de producci´on, pero el individuo no lo acepta desde el primer momento). 1.2.2. Truncamiento El truncamiento tiene lugar cuando s´olo aquellos sujetos que manifiestan el evento dentro de una ventana observacional se observan, del resto no se realiza ning´un seguimiento y, por tanto, no se obtiene informaci´on sobre ellos. El ejemplo m´as claro de truncamiento lo encontramos en el campo de la astronom´ıa: en una parte del espacio, s´olo los elementos suficientemente brillantes pueden observarse desde la Tierra; aquellos cuya intensidad lum´ınica es inferior a un cierto nivel, no es posible saber de su existencia. 1.3. T´ecnicas param´etricas Existen numerosos modelos param´etricos,(Lawless 2003, [15]), que se usan en el an´alisis de tiempos de vida, y en problemas relacionados con la modelizaci´on del envejecimiento y el proceso de fallo. Entre los modelos univariantes, son unas pocas distribuciones las que toman un papel fundamental, dada su demostrada utilidad en casos pr´acticos. As´ı se tiene, entre otras: la exponencial, la Weibull y la log-normal. El m´etodo consiste en estimar, por m´etodos num´ericos (m´axima verosimilitud o m´ınimos cuadrados), los par´ametros caracter´ısticos de la distribuci´on, y usar su normalidad asint´otica para realizar la estimaci´on por intervalos y resolver contrastes de hip´otesis. 20 CAP´ ITULO 1. GENERALIDADES Y MODELOS 1.3.1. Distribuci´on Exponencial La distribuci´on exponencial tiene un papel fundamental en el An´alisis de Fiabilidad, ya que se trata de la distribuci´on m´as b´asica en el an´alisis de datos de tiempo de fallo. Se utiliza para modelar el tiempo transcurrido entre dos sucesos aleatorios, no muy frecuentes, cuando la tasa de ocurrencia, λ, se supone constante. En fiabilidad, se usa para describir los tiempos de fallo de un dispositivo durante su etapa de vida ´util, en la cu´al la tasa de fallo es (aproximadamente) constante, es decir, h(t) = λ. Una tasa de fallo constante, significa que, para un dispositivo que no haya fallado con anterioridad, la probabilidad de fallar en el siguiente intervalo infinitesimal es independiente de la edad del dispositivo. La expresi´on de la funci´on de densidad que sigue una distribuci´on exponencial es: f(t) = λ exp{−λt},0<t<∞, λ > 0. La funci´on de distribuci´on es: F(t) = Zt 0 f(u)du = 1 −exp{−λt},0<t<∞, λ > 0. La funci´on de supervivencia queda como: S(t)=1−F(t) = exp{−λt},0<t<∞, λ > 0. La funci´on de riesgo es: h(t) = f(t) S(t)=λ, 0<t<∞, λ > 0. A continuaci´on, en la Figura 1.5, hemos representado la funci´on de densidad, funci´on de distribuci´on, supervivencia y tasa de riesgo para una distribuci´on exponencial para varios valores de λ. 1.3. T´ ECNICAS PARAM ´ ETRICAS 21 Figura 1.5: Representaci´on de una Exponencial. 1.3.2. Distribuci´on Weibull Un inconveniente de la distribuci´on exponencial es, que no sirve como modelo para tiempos de vida en los que la raz´on de fallo no es una funci´on constante, sino que la probabilidad condicional de fallo instant´aneo var´ıa con el tiempo. Es decir, mientras que la distribuci´on exponencial supone una raz´on de fallo constante, la familia de distribuciones Weibull incluyen razones de fallo crecientes y decrecientes. Como muchos fallos que encontramos en la pr´actica presentan una tendencia creciente, debido al envejecimiento o desgaste, esta distribuci´on es ´util para describir los patrones de este tipo de fallo. La expresi´on de la funci´on de densidad que sigue una distribuci´on Weibull es: f(t) = β αt αβ−1 exp (−t αβ),0<t<∞, α > 0, β > 0. La funci´on de distribuci´on es: F(t) = Zt 0 f(u)du = 1 −exp (−t αβ),0< t < ∞, α > 0, β > 0. La funci´on de supervivencia queda como: S(t) = 1 −F(t) = exp (−t αβ),0<t<∞, α > 0, β > 0. La funci´on de riesgo es: h(t) = f(t) S(t)=β αt αβ−1 ,0< t < ∞, α > 0, β > 0. 22 CAP´ ITULO 1. GENERALIDADES Y MODELOS Esta distribuci´on viene caracterizada por dos par´ametros: α(escala) y β (forma). Figura 1.6: Representaci´on de una Weibull. En la Figura 1.6 se han representado f(t), F(t), S(t) y h(t) para una escala fija (α= 10), y distintos valores del par´ametro de forma. 1.3.3. Distribuci´on Log-normal La distribuci´on normal, es sin duda la m´as importante de las distribuciones estad´ısticas, sin embargo no resulta de mucho inter´es a la hora de modelar tiempos de fallo. Esto es debido al hecho de que la distribuci´on normal admite valores negativos, lo cual contrasta con el hecho de que los tiempos transcurridos hasta el fallo sean siempre valores positivos. Una forma de solventar esta dificultad, es recurrir a la distribuci´on log-normal, relacionada con la normal, y s´olo considera valores positivos. Se dice que una variable aleatoria T > 0 tiene un comportamiento log-normal, de par´ametros µyσ, si su logaritmo es una variable aleatoria con distribuci´on normal, es decir, si ln(T)∼ N(µ, σ). Si T sigue una distribuci´on log-normal se representa por T∼ LN(µ, σ), donde µyσson los par´ametros de localizaci´on ydispersi´on respectivamente, de la distribuci´on de ln(T). 1.3. T´ ECNICAS PARAM ´ ETRICAS 23 La expresi´on de la funci´on de densidad que sigue una distribuci´on log-normal es: f(t) = 1 σ t√2πexp −1 2σ2(ln(t)−µ)2, t > 0. La funci´on de distribuci´on es: F(t) = Zt 0 f(u)du = Φ lnt −µ σ, t > 0. donde Φ(z) representa la funci´on de distribuci´on de una normal est´andar, cuyo c´alculo se obtiene con la integral Φ(z) = Zz −∞ φ(u)du donde φ(u) es la funci´on de densidad de una Normal(0,1). La funci´on de supervivencia depende tambi´en de Φ(z). Se tiene que S(t)=1−F(t) = P(T > t) = PlnT −µ σ>lnt −µ σ= 1 −Φlnt −µ σ, para t > 0. La funci´on de riesgo h(t) = f(t) S(t)tiene valor cero en t= 0, es creciente hasta un m´aximo y despu´es decrece muy lentamente sin llegar a 0. A continuaci´on, en la Figura 1.7 se han representado f(t), F(t), S(t) y h(t) para una distribuci´on log-normal con par´ametro de localizaci´on cero, µ= 0, y distintos valores del par´ametro σ. Estos gr´aficos se han realizado con el paquete ‘eha’ de R [9]. 24 CAP´ ITULO 1. GENERALIDADES Y MODELOS Figura 1.7: Representaci´on de una Log-normal. Cap´ıtulo 2 T´ecnicas no param´etricas En este cap´ıtulo se introducen las t´ecnicas no param´etricas que se utilizan para solventar el problema de que nuestros datos no se ajusten a ninguna distribuci´on conocida (Exponencial, Weibull, Log-normal, entre otras.), y por tanto, no dependen de ning´un par´ametro. Si recordamos lo visto anteriormente en asignaturas como Inferencia Estad´ıstica, se suele usar la funci´on de densidad f(t) y la funci´on de distribuci´on F(t) para analizar y resolver un problema no param´etrico. Sin embargo, en el estudio del an´alisis de tiempos de vida, se usar´an las siguientes funciones: S(t) que es la funci´on de supervivencia para el m´etodo Kaplan-Meier; h(t) que es la funci´on de hazard (o de riesgo) . Esta es muy irregular y dif´ıcil de entender, y por ´ultimo, H(t) que es la funci´on cumulative hazard (o riesgo acumulado) para el m´etodo Nelson-Aalen. 2.1. Introducci´on En general, las t´ecnicas no param´etricas de datos de tiempo de vida son ´utiles para hacer un an´alisis preliminar, que nos puede ayudar a elegir un modelo param´etrico, haciendo uso de las caracter´ısticas de dicho modelo.De esta manera, es mucho m´as f´acil estudiarlo, porque conocemos a qu´e distribuci´on se ajustan los datos, o bien, en el caso de no poder elegir un modelo, las t´ecnicas no param´etricas nos proporcionan herramientas para estudiar datos de supervivencia que no se ajusten a ning´un modelo. Un ejemplo de ello, es lo que ocurre con datos relacionados con el comportamiento humano. Estos no suelen seguir las mismas pautas, y ser´a m´as dif´ıcil ajustar un modelo, que en el caso de estudios relacionados con caracter´ısticas de m´aquinas. 25 32 CAP´ ITULO 2. T´ ECNICAS NO PARAM ´ ETRICAS En general, para i: ˆ P(T > t(i)/T > t(i−1)) = ˆ S(t(i)) = 1 −pi donde pi=di ni es la estimaci´on de la proporci´on de fallos en t(i), sabiendo que niunidades est´an funcionando y fallan di. Por tanto, ˆ P(T > t(i)/T > t(i−1)) = 1 −di ni =ni−di ni . La relaci´on anterior motiva la propuesta del m´etodo Kaplan-Meier para estimar S(t) de manera: Para t < t(1) =⇒ˆ S(t) = 1 ya que antes de t(1) no ha fallado ninguno, siempre suponiendo que ning´un dato antes de t(1) haya sido censurado. Para t≥t(1) =⇒tomamos ital que t(i)≤t≤t(i+1) y analizamos cada caso: •Si i= 1 =⇒t(1) ≤t < t(2) ˆ S(t(1)) = ˆ P(T > t(1)) = 1 −p1=n1−d1 n1 •Si i= 2 =⇒t(2) ≤t<t(3) ˆ S(t(2)) = ˆ P(T > t(1))ˆ P(T > t(2)/T > t(1)) = n1−d1 n1·n2−d2 n2 •En general, sea i=⇒t(i)≤t ˆ S(t(i)) = P(T > t(1))P(T > t(2)/T > t(1)). . . P(T > t(i)/T > t(i−1)) = =n1−d1 n1·n2−d2 n2····· ni−di ni En resumen: ˆ S(t) =     Y j:t(j)≤t nj−dj nj ,para t≥t(1). 1,para t<t(1). 2.2. M´ ETODO KAPLAN-MEIER 33 Calculamos el estimador de Kaplan-Meier de la funci´on de supervivencia, ˆ S(t), en el siguiente ejemplo. Ejemplo 2.2.2 Consideramos los datos dados en el Ejemplo 2.2.1 donde aparecen los tiempos y fallos observados junto con las unidades que fallan. t(j)njdj nj−dj nj ˆ S t(1) 6 1 n1−d1 n1 = 5/6 5/6 (c1) t(2) 4 2 n2−d2 n2 = 1/2 5/6*1/2=5/12 t(3) 2 1 n3−d3 n3 = 1/2 5/6*1/2*1/2=5/24 (c2) Los datos que han sido censurados se representan por ci. En este caso, hay un dato censurado c1entre t(1) yt(2) y otro dato c2despu´es de t(3). Veamos gr´aficamente el estimador de Kaplan-Meier que se ajusta a estos datos: Figura 2.4: M´etodo de Kaplan-Meier estimando S(t). Para dibujarlo hemos supuesto que el valor de los tiempos de fallo observados son: t(1) = 2, t(2) = 4 y t(3) = 6. 34 CAP´ ITULO 2. T´ ECNICAS NO PARAM ´ ETRICAS Propiedad 2.2.1 El estimador Kaplan – Meier es el estimador no param´etrico, m´aximo verosimil de la funci´on de supervivencia, lo cu´al hace que presente algunas de las propiedades de un buen estimador, como ser insesgado, consistente, eficiente y suficiente. Este estimador cumple con varias de estas propiedades ya que es consistente (Efron, 1967) y eficiente (Wellner, 1982). Adem´as, es un estimador de m´axima verosimilitud para datos censurados (Peterson, 1997) y es asint´oticamente normal (Breslow & Crowley, 1974). Estas propiedades facilitan el c´alculo de este estimador y la utilizaci´on en problemas con datos censurados por la derecha. Cuando las estimaciones de S(t) se realizan con datos no censurados, coincide con el estimador no param´etrico de la funci´on de supervivencia. Lo anterior hace que el estimador Kaplan–Meier sea un buen estimador para muestras grandes. Para muestras muy peque˜nas, hemos de tener en cuenta, que las propiedades asint´oticas ya no se cumplen. 2.2.1. Varianza del estimador de Kaplan-Meier ˆ S(t) El estimador de Kaplan-Meier da una estimaci´on puntual, o un ´unico valor de la funci´on de supervivencia en cualquier instante t. Por lo tanto, si se desea tener una precisi´on de este estimador, en diferentes instantes de tiempo, o sobre diferentes muestras, es necesario contar con un buen estimador de la varianza. Para ello usaremos la f´ormula de Greenwood que nos da una expresi´on aproximada para la varianza: d V ar(ˆ S(t)) = ˆ S2(t)  X t(j)≤t dj nj(nj−dj)  (2.4) Adem´as, la desviaci´on est´andar o error est´andar, se calcula como la ra´ız cuadrada de la varianza dada en (2.4). s.e.(ˆ S(t)) = ˆ S(t)  X t(j)≤t dj nj(nj−dj)   1/2 Si ahora suponemos que ˆ S(t), en un momento determinado t, sigue una distribuci´on normal, con media µigual al valor verdadero de S(t), y varianza σ2, entonces un intervalo de confianza aproximado, que es asint´otico en el sentido de que ksea grande, ser´ıa I.C (S(t); 1 −α) = ˆ S(t)−Z1−α/2s.e.(ˆ S(t)) ; ˆ S(t) + Z1−α/2s.e.(ˆ S(t)) donde Z1−α/2representa el cuantil al nivel 1 −α/2 de N(0,1). 2.2. M´ ETODO KAPLAN-MEIER 35 Si lo queremos al 95 % : 1−α= 0.95 →α= 0.05 // Z1−α/2=Z0.975 = 1.96 I.C (S(t); 0.95) = ˆ S(t)−1.96 s.e.(ˆ S(t)) ; ˆ S(t)+1.96 s.e.(ˆ S(t)) En el siguiente ejemplo realizamos una aplicaci´on donde se calcula ˆ S(t),s.e.(ˆ S(t)) eI.C al 95 %, entre otras. Ejemplo 2.2.3 Consideramos como tiempos de fallo, (obtenidos de Caroni [6]), T= “n´umero de millones de revoluciones”. Se estudian 25 rodamientos de bolas de cer´amica antes de que se produzca un fallo (n=25). En este caso, tenemos 6 observaciones censuradas por la derecha (representadas con un *): 17.88 28.92 33.00 41.52 42.12 45.60 48.48 51.84 51.96 54.12 55.56 67.80 67.80* 67.80* 68.64 68.64* 68.88* 84.12 93.12 98.64 105.12 105.84* 127.92 128.04 173.40* Usando el m´etodo de Kaplan-Meier, sabemos que ˆ S(t) es una funci´on escalonada, con saltos s´olo en los tiempos de fallo observados t(i). En este caso, hay 19 saltos (ya que 6 de las observaciones han sido censuradas, por lo tanto, no entran en el an´alisis), y ˆ Sno cambia en los 6 tiempos en los que las observaciones han sido censuradas. Figura 2.5: Informaci´on obtenida con el programa R, con los datos de n´umero de millones de revoluciones de bolas de cer´amica. 36 CAP´ ITULO 2. T´ ECNICAS NO PARAM ´ ETRICAS Veamos gr´aficamente el estimador de Kaplan-Meier de la curva de supervivencia que se ajusta a estos datos: Figura 2.6: M´etodo de Kaplan-Meier estimando S(t). Para poder dibujarlo, adem´as de la variable tiempo (dado en la tabla), hemos introducido censura (aquellos valores que est´an censurados y aquellos que no). 2.3. M´etodo Nelson-Aalen Este estimador fue propuesto, por primera vez por Nelson W. Aalen (1969), y luego por Altschuler (1970), qui´en lo descubri´o utilizando t´ecnicas de conteo con animales. Recordemos que, por definici´on, la funci´on de riesgo acumulada (o cumulative hazard rate) es de la forma: H(t) = Zt 0 h(u)du (2.5) con h(t) la funci´on de riesgo (o fallo) en el instante t. La expresi´on (2.5) tambi´en se considera como la suma de las probabilidades de fallo en el intervalo (0, t] por lo que se sugiere el siguiente estimador: ˆ H(t) = X j:t(j)≤t dj nj donde t(j)representa los tiempos observados, djel n´umero de fallos ocurridos en el instante t(j), y njel n´umero de individuos en riesgo antes de t(j). 2.3. M´ ETODO NELSON-AALEN 37 Propiedad 2.3.1 El estimador de Nelson-Aalen es otra funci´on de escalonada, luego ˆ Hpuede ser ´util para ayudar a entender lo que est´a sucediendo, y as´ı poder elegir un modelo. Ejemplo 2.3.1 Si ˆ Hes aproximadamente lineal entonces hes aproximadamente constante, como corresponder´ıa al modelo exponencial. 2.3.1. Varianza del estimador de Nelson-Aalen H Se propone: d V ar(ˆ H(t)) = X j:t(j)≤t dj n2 j A continuaci´on vamos a construir el estimador de Nelson-Aalen, a partir de ˆ S(t), ya que existe un relaci´on entre ellos. 2.3.2. Estimador de Nelson-Aalen, ˆ H(t), a partir del estimador de Kaplan-Meier, ˆ S(t) Usando la Propiedad 3 de la Proposici´on 1.1.3, vista en el Cap´ıtulo 1, se tiene: S(t) = e−H(t)=⇒H(t) = −ln S(t) (2.6) Para el caso de una variable continua, este estimador ha tenido mucha discusi´on por parte de Nelson (1972), Breslow y Crowley (1974), Efron (1977) y Altschuler (1979). Los cu´ales llegaron a la conclusi´on que en este caso ˆ H(t) y ˆ S(t) son asint´oticamente equivalentes, con la excepci´on de valores altos de t, en los que las estimaciones son menos estables. La diferencia entre ellos ser´a, por lo general peque˜na, luego no existe ninguna raz´on suficiente para escoger alguna de estas. Una de las utilidades de las estimaciones ˆ H(t) y ˆ S(t) es, en la construcci´on de gr´aficas para evaluar la selecci´on de una determinada familia param´etrica de distribuciones. Comenzamos con la construcci´on del estimador de Nelson-Aalen, paso a paso. Por (2.6) sabemos que H(t) = −ln S(t) luego si: t < t(1) =⇒ˆ H(t) = −ln ˆ S(t) = −ln 1 = 0 38 CAP´ ITULO 2. T´ ECNICAS NO PARAM ´ ETRICAS t≥t(1) =⇒ˆ H(t) = −ln ˆ S(t) = −ln Y j:tj≤t nj−dj nj =−X j:tj≤t ln nj−dj nj= =−X j:tj≤t ln 1−dj nj Lema 2.3.1 (Aproximaci´on por Taylor) Dado uj<1, se tiene que: −ln(1 −uj) = ∞ X l=1 ul j l. Por tanto, se puede aproximar −ln(1 −u)≈upara u < 1, y en nuestro caso se aplicar´a uj=dj nj <1, porque dj< nj. Finalmente podemos aproximar −ln 1−dj nj≈dj nj . En resumen, definimos el estimador de Nelson-Aalen como: ˆ H(t) =     X j:t(j)≤t dj nj , para t≥t(1). 0,para t < t(1). (2.7) luego este estimador, ˆ H, nos ofrece una alternativa al estimador ˆ S. Por otra parte, sabiendo por (2.7) que, ˆ H(t) = X j:t(j)≤t dj nj y por (2.6) que S(t) = e−H(t), se tiene: ˆ S(t) = exp n−ˆ H(t)o=exp   −X j:t(j)≤t dj nj  =Y j:t(j)≤t exp −dj nj que es el llamado estimador de Altshuler de S. Propiedad 2.3.2 Una de las propiedades m´as importantes de este estimador es: Altshuler ≥Kaplan-Meier,para cada t. La diferencia entre los dos es muy peque˜na, y es equivalente, para tpeque˜nos. En el siguiente ejemplo demostraremos de qu´e manera se utiliza el m´etodo de Nelson-Aalen, aplicado a unos datos. 2.4. TEST LOG-RANK 39 Ejemplo 2.3.2 Consideramos los datos utilizados en el Ejemplo 2.2.3 teniendo en cuenta todas las variables que intervienen. Figura 2.7: Estimaci´on utilizando el m´etodo Nelson-Aalen. 2.4. Test Log-Rank Para comparar la supervivencia de dos o m´as grupos de observaciones, necesitamos un test estad´ıstico apropiado. El fin es determinar si todos los grupos presentan la misma supervivencia, o qu´e grupos son distintos. Las representaciones gr´aficas de las curvas de supervivencia dan una alerta sobre posibles diferencias entre estas, pero con las pruebas estad´ısticas, se obtiene si existen diferencias m´as significativas entre las curvas, indic´andonos que el factor considerado influye de forma importante en el riesgo de que falle. Ejemplo 2.4.1 Un ejemplo de esta situaci´on ser´ıa analizar, si los pacientes sobreviven m´as tiempo con el nuevo tratamiento que con el establecido o si las unidades fabricadas en A duran m´as que las fabricadas en B. Para comparar la igualdad de dos o m´as funciones de supervivencia (fiabilidad) con datos censurados, se presentan los siguientes contrastes no param´etricos: Test de Log-Rank (riesgos proporcionales): es muy potente para calcular diferencias cuando los logaritmos de las curvas de supervivencias son proporcionales, pero tiene muchos problemas para detectar las diferencias cuando las curvas de supervivencias se cruzan. 40 CAP´ ITULO 2. T´ ECNICAS NO PARAM ´ ETRICAS Test de Breslow (test de Wilcoxon generalizado): detecta las diferencias cuando las curvas de supervivencia se cruzan o se cortan, pero solamente al principio, por lo cual no es recomendable para un estudio a largo plazo. Test de Tarone-Ware : es un test intermedio a los otros dos. En este caso, nos centraremos en el Test de Log-Rank. Es el test m´as potente cuando el cociente de las funciones de riesgo es aproximadamente constante. Supongamos que se va a comparar la supervivencia de dos grupos A1yA2 donde la funci´on de supervivencia es, respectivamente, S1yS2. Adem´as, se tiene una muestra de cada poblaci´on, con su respectivo tama˜no n1yn2, donde n=n1+n2es el n´umero total de datos en la muestra combinada. Los tiempos de fallo se definen como t(1) < t(2) < . . . < t(k). La comparaci´on de las dos curvas de supervivencias se efect´ua a trav´es de contrastes basados en tablas de contingencia, como la siguiente: GRUPO EVENTO A1A2TOTAL Muerte d1jd2jdj No muerte n1j−d1jn2j−d2jnj−dj EN RIESGO n1jn2jnj Definimos el contraste de hip´otesis de la siguiente manera: H0:S1(t) = S2(t) H1:S1(t)6=S2(t) donde tes el tiempo total de observaci´on de la muestra. En este test de Log-Rank se compara el n´umero de fallos observados dentro de cada grupo A1yA2, adem´as del n´umero de fallos esperados bajo la hip´otesis nula. Cuando la hip´otesis nula es cierta, es decir la funci´on de supervivencia es igual en ambas poblaciones, la probabilidad condicional de fallo en t(j)es igual para los dos grupos λi, por lo tanto, la distribuci´on de probabilidad de (d1j, d2j) est´a dada de la siguiente manera: 2.4. TEST LOG-RANK 41 2 Y i=1 nij dij λdij j(1 −λi)nj−dij = 2 Y i=1 nij dij λdij j(1 −λi)nj−dij donde: dj= n´umero total de fallos ocurridos en el tiempo t(j). nj= n´umero total de unidades en riesgos antes de t(j). dij = el n´umero de fallos ocurridos en el tiempo t(j)entre los individuos del grupo i(= 1,2). nij = el n´umero de unidades en riesgo al principio de t(j)entre los individuos del grupo i(= 1,2). Como las funciones de supervivencia coinciden (debido a la hip´otesis nula), la funci´on de riesgo es la misma en ambas poblaciones, por lo que el fallo es independiente del grupo, lo que implica que los fallos esperados en el grupo 1 viene dado por: ei1=n1jdj nj Por tanto, se define el estad´ıstico de Log-Rank: ui= k X j=1 (dij −eij) es decir, los fallos observados menos los esperados. En el caso de que el valor de ksea lo suficientemente grande se tiene por el teorema central del l´ımite lo siguiente: u √v= k X j=1 (d1j−e1j) v u u t k X j=1 vj ∼ N(0,1) donde vjes la varianza de d1jusando la distribuci´on hipergeom´etrica de la forma: vj=n1jn2jdj(nj−dj) n2 j(nj−1) Por lo tanto, el test estad´ıstico Log-Rank se define: u2 √v∼ X2 1 48 CAP´ ITULO 3. MODELO DE REGRESI ´ ON DE COX Ejemplo 3.1.1 El sexo, la raza o el grupo de tratamiento son variables fijas, s´olo toman un valor, el inicial. Tambi´en podr´ıamos considerar variables como el hecho de ser o no fumador (estado de fumador) como variable independiente del tiempo, ya que aunque el estado de fumador puede variar en el tiempo, para el estudio no var´ıa, porque se parte de un estado inicial, y se supone que no cambia hasta el final, y por lo tanto s´olo toma un valor por individuo. El modelo de Cox, definido en (3.1), se considera un modelo semiparam´etrico, debido a que incluye una parte param´etrica y otra parte no param´etrica: 1. La parte param´etrica se corresponde con exp(β0X), es decir, con la exponencial del predictor lineal η=β0X. En esta parte del modelo se estiman los par´ametros o coeficientes βde la regresi´on mediante la maximizaci´on de la denominada funci´on de verosimilitud parcial que estudiaremos con detalle m´as adelante. 2. La parte no param´etrica es la funci´on de riesgo basal h0(t). Esta es una funci´on arbitraria condicionada a la estimaci´on de los par´ametros β de la regresi´on. Es por esta componente no param´etrica de la f´ormula que el modelo de Cox se considera semiparam´etrico. Una vez estimada la parte param´etrica, exp(ˆ β0X), y posteriormente la no param´etrica ˆ h0(t), tendremos la estimaci´on del modelo semiparam´etrico completo: ˆ h(t;X) = ˆ h0(t)exp(ˆ β0X) 3.2. Hip´otesis de riesgos proporcionales, PH En el modelo de Cox se busca como primer paso, la relaci´on entre los riesgos de muerte de dos individuos expuestos a factores de riesgo diferentes. Para ello, el modelo parte de una hip´otesis fundamental, la de que los riesgos son proporcionales. Definici´on 3.2.1 Se define la raz´on de riesgo(o hazard ratio) entre dos sujetos con diferente vector de covariables X= (X1, . . . , Xp)0yX∗= (X∗ 1, . . . , X∗ p) como h(t;X∗) h(t;X).(3.3) 3.2. HIP ´ OTESIS DE RIESGOS PROPORCIONALES, PH 49 Se eval´ua en el numerador el grupo de mayor riesgo, definido por X∗, y en el denominador el grupo de menor riesgo, definido por X. Se espera que la raz´on de riesgo sea mayor que 1, ya que h(t;X∗)> h(t;X) cuantifica cu´antas veces es mayor el riesgo de morir con perfil X∗que con X, luego el tiempo de supervivencia disminuye. Si sustituimos en la expresi´on (3.3), el modelo dado en (3.1) obtenemos h(t;X∗) h(t;X)=h0(t)exp(β0X∗) h0(t)exp(β0X)= exp p X j=1 βjX∗ j! exp p X j=1 βjXj!= =exp p X j=1 βj(X∗ j−Xj)!=exp((X∗−X)0β) Por tanto exp((X∗−X)0β) es igual al cociente (3.3) h(t, X∗) h(t, X).(3.4) Observamos que el resultado de la raz´on de riesgo dada en (3.4) no depende de la funci´on de riesgo basal h0(t), tan solo del valor de los predictores y de las betas estimadas, es decir, no depende del tiempo. En resumen, en el modelo de Cox se supone la hip´otesis de que los riesgos son proporcionales, ya que las covariables se suponen no dependientes del tiempo. Veamos un ejemplo en el que no se cumple la hip´otesis de riesgos proporcionales. Ejemplo 3.2.1 Consideremos un estudio en el que los pacientes con c´ancer se reparten aleatoriamente a radioterapia con o sin cirug´ıa. Definimos una variable que toma los valores 0 y 1 que indica si ha habido cirug´ıa o no, respectivamente, y supongamos que ´esta es la ´unica variable predictora de inter´es en el modelo de Cox. La raz´on de riesgo se calcular´ıa como ˆ h(t, X = 1) ˆ h(t, X = 0) Cuando un paciente recibe cirug´ıa para eliminar un tumor canceroso hay, por lo general, un alto riesgo de complicaciones por la propia cirug´ıa, o quiz´as de muerte temprana, en el proceso de recuperaci´on. 50 CAP´ ITULO 3. MODELO DE REGRESI ´ ON DE COX Una vez que el paciente pasa este periodo cr´ıtico, es cuando se pueden observar las ventajas de la cirug´ıa. Supongamos que se dibujan las curvas de riesgo para los dos grupos y que se observa que se cruzan aproximadamente en el d´ıa 3, que antes del d´ıa 3 el riesgo para el grupo de cirug´ıa es m´as alto que el riesgo para el grupo sin cirug´ıa, mientras que despu´es de 3 d´ıas, el riesgo para el grupo de cirug´ıa es inferior que el riesgo para el grupo sin cirug´ıa. Figura 3.1: Gr´afico en el que se representan las curvas de riesgo para el grupo sin cirug´ıa y el grupo con cirug´ıa. Escrito en t´erminos de la raz´on de riesgo se tendr´ıa lo siguiente: t=2: ˆ h(2; X= 1) ˆ h(2, X = 0) <1 t=3: ˆ h(3; X= 1) ˆ h(3, X = 0) = 1 t=5: ˆ h(5; X= 1) ˆ h(5, X = 0) >1 Se observa que la raz´on de riesgo que obtenemos no es constante a lo largo del tiempo t. Se encontrar´a alguna soluci´on a este problema cuando se estudie el modelo de Cox estratificado, el cual permite la utilizaci´on de predictores dependientes del tiempo, y que veremos m´as adelante. 3.2. HIP ´ OTESIS DE RIESGOS PROPORCIONALES, PH 51 En resumen, si los riesgos del numerador y denominador de la raz´on de riesgo se cruzan, la hip´otesis de riesgos proporcionales no se cumple, y por lo tanto el modelo (3.1) de riesgos proporcionales no es adecuado. 3.2.1. Diagrama de diagn´ostico para PH En esta secci´on introducimos una forma de decidir si el modelo es adecuado o no. Para ello utilizaremos el diagrama de diagn´ostico para riesgos proporcionales (PH), el cual explicamos a continuaci´on. Lema 3.2.1 Bajo la hip´otesis (PH), la funci´on de supervivencia para cualquier t es: S(t;X) = S0(t)ˆ{exp(β0X)}(3.5) Demostraci´on 3.2.1 Por la Propiedad 3 de la Proposici´on 1.1.3 del Cap´ıtulo 1 se tiene S(t;X) = exp{−H(t;X)} Sustituyendo el modelo dado en (3.1) en la expresi´on anterior S(t;X) = exp{−H0(t)exp(β0X)}={exp(−H0(t))}ˆ{exp(β0X)} Si llamamos S0(t) = exp(−H0(t)) (3.6) nos queda lo siguiente S(t;X) = S0(t)ˆ{exp(β0X)}  Una vez demostrada la expresi´on (3.5), continuamos aplicando logaritmo neperiano a ambos lados de la igualdad: lnS(t, X) = ln[S0(t)ˆ{exp(β0X)}] = {exp(β0X)}lnS0(t) Como S(t, X)<1 entonces lnS(t, X)<0. Para solucionar este problema, multiplicamos todo por (-1), para convertirlo en positivo, (−lnS(t, X)<0), y obtenemos as´ı: −lnS(t, X) = exp(β0X){(−lnS0(t)} Volvemos a aplicar logaritmo neperiano: ln{−lnS(t, X)}=ln[exp(β0X){(−lnS0(t))}] = =ln{exp(β0X)}+ln{−lnS0(t)}=β0X+ln{−lnS0(t)} 52 CAP´ ITULO 3. MODELO DE REGRESI ´ ON DE COX Por la definici´on dada en (2.6) sabemos que H0(t) = −lnS0(t) luego podemos escribir: ln{−lnS(t, X)}=β0X+ln{H0(t)} lo que significa que la curva ln{−lnS(t, X)}para cada Xdebe ser paralela aln{H0(t)}en el tiempo, luego si se tiene la hip´otesis (PH) se mantiene las curvas de supervivencia para diferentes X0sdeben ser paralelas. Veamos un ejemplo en el que se cumple la hip´otesis de riesgos proporcionales. Ejemplo 3.2.2 Consideramos los tiempos de fallo (en millones de revoluciones) de rodamientos de bolas de cer´amica dividi´endolas, seg´un hayan salido de la f´abrica 1 o de la f´abrica 2 (datos extra´ıdos del Ejemplo 2.2.3): Figura 3.2: Riesgos proporcionales del rodamiento de bolas cer´amicas. En este caso, el gr´afico sugiere, que la hip´otesis PH es satisfactoria, ya que ambas curvas son paralelas. 3.3. Estimaci´on de los coeficientes En el modelo de regresi´on de Cox los par´ametros β= (β1, . . . , βp)0se estiman maximizando el logaritmo de la denominada funci´on de verosimilitud parcial. La maximizaci´on de dicha funci´on se realiza mediante m´etodos num´ericos, obteniendo de esta forma la estimaci´on ˆ β= (ˆ β1,...,ˆ βp)0. Con la estimaci´on de estos par´ametros ya tendremos la componente param´etrica totalmente especificada en el modelo h(t, X) = h0(t)exp(ˆ β0X) 3.3. ESTIMACI ´ ON DE LOS COEFICIENTES 53 y consecuentemente, podremos hacer inferencia sobre dicho vector de par´ametros, para calcular la raz´on de riesgos, de inter´es en el estudio. La funci´on de verosimilitud parcial, que a continuaci´on vamos a definir, se denomina parcial, debido a que, se tiene en cuenta ´unicamente, en la funci´on de verosimilitud, las probabilidades de los tiempos de muerte/fallo, y no incluye las probabilidades de los tiempos de datos censurados. Sin embargo, en el c´alculo de las probabilidades de los tiempos, se tiene en cuenta a todos los individuos (censurados o no a posteriori) en riesgo, al inicio de los diferentes tiempos de muerte/fallo. Denominamos L≡L(β1, . . . , βp) a la funci´on de verosimilitud parcial. Supongamos que tenemos ktiempos de muerte, y que no hay empates. As´ı, tendremos n−ktiempos censurados. Notaci´on 3.3.1 Sea t(1), . . . , t(k), los tiempos de muerte ordenados. Sea Ripara i= 1, . . . , k, el conjunto de los individuos que est´an en riesgo antes del tiempo t(i). Definimos por Li≡ Lt(i)(β1, . . . , βp) para i= 1, . . . , k a cada porci´on de la verosimilitud parcial anterior perteneciente a los diferentes tiempos de muerte t(1), . . . , t(k). Construiremos la funci´on de verosimilitud parcial como el producto de cada una de las aportaciones de los ktiempos de muerte: i= 1 −→ L1≡ Lt(1) (β1, . . . , βp) i= 2 −→ L2≡ Lt(2) (β1, . . . , βp) · · · i=k−→ Lk≡ Lt(k)(β1, . . . , βp)                =⇒ L = k Y i=1 Li Veamos cuanto vale exactamente cada una de las Li≡ Lt(i)(β1, . . . , βp) para i= 1, . . . , k. Para cada unidad l∈Ri h(t(i);Xl)δt =P(fallo en t(i), t(i)+δt) dado que la funci´on de riesgo en t(i)da la probabilidad instant´anea de fallo condicionada a haber alcanzado t(i)cuyas unidades en Rilo han hecho 54 CAP´ ITULO 3. MODELO DE REGRESI ´ ON DE COX P(de que un individuo ifalle/a que haya un fallo en t(i)) = =h(t(i), Xi) X l∈Ri h(t(i), Xl)=h0(t(i))exp(β0Xi) X l∈Ri h0(t(i))exp(β0Xl)=h0(t(i))exp(β0Xi) h0(t(i))X l∈Ri exp(β0Xl) Esta ´ultima igualdad se tiene porque h0(t(i)) no depende de l, lo que implica que se anule tanto en el numerador, como en el denominador. Xies el vector de covariables para el individuo con tiempo de muerte t(i)yXl, para l∈Ri, el vector de covariables de cada uno de los individuos de Ri. Por tanto nos queda la expresi´on: Lt(i)(β1, . . . , βp) = exp(β0Xi) X l∈Ri exp(β0Xl) Se observa que la funci´on de verosimilitud parcial, as´ı calculada, no depende de las cuant´ıas de los tiempos, sino tan solo de su ordenaci´on, y de si el dato estaba o no censurado. Como consecuencia se podr´ıa obtener las mismas estimaciones de βpara distintos datos, siempre que estos tengan el mismo patr´on de orden y censura en los tiempos de supervivencia. En resumen, la funci´on de verosimilitud parcial queda de la siguiente manera: L(β) = k Y i=1        exp(β0Xi) X l∈Ri exp(β0Xl)       Una vez que tenemos la funci´on de verosimilitud parcial construida, procedemos a calcular ˆ β, las estimaciones de los coeficientes del modelo β: Sea ˆ β= m´ax ln L(β),calculamos: (3.7) ln L(β) = k X i=1 ln        exp(β0Xi) X l∈Ri exp(β0Xl)       = k X i=1 (ln exp(β0Xi)−ln X l∈Ri exp(β0Xl))= = k X i=1 β0Xi− k X i=1 ln (X l∈Ri exp(β0Xl)) 3.3. ESTIMACI ´ ON DE LOS COEFICIENTES 55 A continuaci´on derivamos: ∂ln L ∂βj =∂ ∂βj k X i=1 β0Xi− k X i=1 ln (X l∈Ri exp(β0Xl))!= =∂ ∂βj k X i=1 β0Xi! | {z } (1) −∂ ∂βj k X i=1 ln (X l∈Ri exp(β0Xl))! | {z } (2) Calculamos (1): ∂ ∂βj k X i=1 β0Xi!=∂ ∂βj k X i=1 βjXij!= k X i=1 ∂ ∂βj (βjXij) = k X i=1 Xij Calculamos (2): ∂ ∂βj k X i=1 ln (X l∈Ri exp(β0Xl))!= k X i=1 ∂ ∂βj ln (X l∈Ri exp(β0Xl))!= = k X i=1     X l∈Ri ∂ ∂βj{exp (βjXjl)} X l∈Ri exp(β0Xl)     = k X i=1    X l∈Ri Xjl exp(β0Xl) X l∈Ri exp(β0Xl)    Por tanto: ∂ln L ∂βj = k X i=1 Xij − k X i=1    X l∈Ri Xjl exp(β0Xl) X l∈Ri exp(β0Xl)   (3.8) Derivando por segunda vez ∂2ln L ∂βi∂βj ,(3.9) Si igualamos (3.8) a 0, para j= 1, . . . , p, obtendremos las ecuaciones que nos permitir´an obtener las estimaciones de ˆ β= ( ˆ β1,..., ˆ βp) mediante la utilizaci´on de alg´un m´etodo num´erico. De (3.9), comprobamos que realmente es un m´aximo, y podemos obtener, como ocurre cuando se trabaja en general con una funci´on de verosimilitud, la matriz de informaci´on (observada),I(β), donde cada elemento se iguala a: 56 CAP´ ITULO 3. MODELO DE REGRESI ´ ON DE COX Iij(β) = −∂2ln L ∂βi∂βj . As´ı, la matriz de varianzas y covarianzas estimada (pxp) es ˆ Σ = I−1(ˆ β). Cabe notar que este estimador, obtenido a partir de la maximizaci´on de la funci´on de verosimilitud parcial, es asint´oticamente no sesgado, eficiente y normal. Adem´as, aunque el estimador ˆ βestime consistentemente el vector de par´ametros βno es completamente eficiente. Finalmente, la distribuci´on de ˆ β= ( ˆ β1,..., ˆ βp) es aproximadamente normal, de media (β1, . . . , βp) y matriz de varianzas y covarianzas Σ. 3.3.1. Contrastes de hip´otesis Tras el ajuste del modelo de Cox, se ha de comprobar si las variables del modelo son significativas. Para ello, existen pruebas que se encargan de validar las correspondientes hip´otesis. En ellas, se considera el vector de par´ametros estimados ˆ β= ( ˆ β1,..., ˆ βp)0, y la matriz de varianzas y covarianzas estimada, ˆ Σ, lo que nos permite utilizar tests an´alogos a los utilizados en un modelo lineal, o lineal generalizado. Para contrastar la hip´otesis H0:βj= 0 vs. H1:βj6= 0, podemos utilizar el estad´ıstico de Wald dado por: z=ˆ βj qˆ V ar(ˆ βj) =ˆ βj s.e.(ˆ βj)∼N(0,1) asint´oticamente. La f´ormula para calcular un intervalo de confianza aproximado al nivel (1−α) para el coeficiente βjes la siguiente: ˆ βj±Z1−α/2qˆ V ar(ˆ βj) = ˆ βj±Z1−α/2s.e.(ˆ βj). En cambio, si lo que queremos hacer es el test H0:β=β0vs. H1:β6=β0se suelen utilizar tres contrastes: 1. El contraste de Wald. Este contraste se basa en que los coeficientes ˆ β= ( ˆ β1,..., ˆ βp) siguen una distribuci´on aproximadamente normal de media (β1, . . . , βp) y matriz de varianzas y covarianzas ˆ Σ = I−1(ˆ β). El estad´ıstico se define como: X2 W= (ˆ β−β0)0I(ˆ β)(ˆ β−β0), que bajo la hip´otesis nula sigue una distribuci´on X2con pgrados de libertad. 3.4. RESIDUOS EN EL AN ´ ALISIS DE SUPERVIVENCIA 57 2. El contraste de la raz´on de verosimilitud. En este contraste se utiliza el valor de la funci´on de verosimilitud parcial evaluada en ˆ β(L(ˆ β)) y evaluada en β0(L(ˆ β0)): X2 LR =−2{lnL(β0)−lnL(ˆ β)}, que bajo la hip´otesis nula sigue una distribuci´on X2. 3. El contraste del “score” (Log Rank). En este contraste se utiliza el gradiente (derivadas) del logaritmo de la funci´on de verosimilitud parcial evaluada en la hip´otesis nula y supone que bajo la hip´otesis nula el vector scores: X2 SC =∂L(β0) ∂β 0 −∂L2(β0) ∂β∂β0!−1∂L(β0) ∂β , sigue una distribuci´on aproximada X2. Nota 3.3.1 El contraste de Wald tiene una interpretaci´on m´as directa, que el contraste de verosimilitud y el del score, sin embargo no es invariante ante diferentes parametrizaciones y los otros dos, s´ı. Con el contraste del score s´olo hace falta maximizar bajo la hip´otesis nula, con lo que si hay que realizar el test para varios par´ametros, es m´as r´apido computacionalmente. Sin embargo, el test de la m´axima verosimilitud converge m´as r´apido hacia la distribuci´on normal. Ante la duda de cu´al utilizar, es recomendable decantarse por el test de la raz´on de verosimilitudes. 3.4. Residuos en el an´alisis de supervivencia Una de las ventajas que han surgido del enfoque del an´alisis de supervivencia es la posibilidad de efectuar an´alisis de residuos. Los residuos se pueden utilizar para: 1. Descubrir la forma funcional correcta de un predictor continuo. 2. Identificar los sujetos que est´an pobremente pronosticados por el modelo. 3. Identificar los puntos o individuos de influencia. 4. Verificar el supuesto de riesgo proporcional. 64 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R Preparaci´on de los datos En primer lugar, para realizar el estudio sobre dicho conjunto de datos, los cargamos con la funci´on read.table(): cdat<-read.table("heroina.txt",header=TRUE,row.names =1) head(cdat) ## cliente clinica censura tiempo prision dosis ## 1 244 0 0 2 1 60 ## 2 202 0 1 7 1 40 ## 3 190 0 1 17 1 40 ## 4 247 0 1 19 1 40 ## 5 220 0 0 28 0 50 ## 6 230 0 0 28 0 50 Veamos, con la funci´on str(), el aspecto de estos datos: attach(cdat) str(cdat) ## 'data.frame': 238 obs. of 6 variables: ## $ cliente: num 244 202 190 247 220 230 203 212 261 248 ... ## $ clinica: num 0 0 0 0 0 0 0 0 0 0 ... ## $ censura: num 0 1 1 1 0 0 1 1 1 1 ... ## $ tiempo : num 2 7 17 19 28 28 29 30 33 35 ... ## $ prision: num 1 1 1 1 0 0 1 0 1 0 ... ## $ dosis : num 60 40 40 40 50 50 60 60 60 60 ... La informaci´on que nos devuelve es que es una muestra de 238 datos con 6 variables. El objeto Surv Hemos de preparar los datos para realizar un estudio de supervivencia con el paquete estad´ıstico survival mediante la funci´on survfit() install.packages("survival") library(survival) Un objeto Surv no es m´as que la combinaci´on de informaci´on entre los tiempos y su censura. 65 Es necesario trabajar con los datos en este formato para m´as tarde aplicar algunas t´ecnicas. sup<-Surv(tiempo,censura) Vemos como obtener si la observaci´on del evento no est´a censurada Surv(5,1) ## [1] 5 Si la observaci´on del evento est´a censurada Surv(5,0) ## [1] 5+ Si queremos los primeros 6 t´erminos head(sup) ## [1] 2+ 7 17 19 28+ 28+ como vemos, se representa con un “+” a la derecha del dato, aquel que est´a censurado. Estimaci´on no param´etrica de la funci´on de supervivencia Estimador de Kaplan-Meier Recordamos por teor´ıa, visto en el Cap´ıtulo 2, la finalidad de este m´etodo es que la proporci´on acumulada que se mantiene en el tratamiento se calcula para el tiempo de supervivencia individual de cada paciente y no se agrupan los tiempos en intervalos. Adem´as, matem´aticamente su f´ormula viene dada por: ˆ S(t) = Y j:t(j)≤t nj−dj nj El estimador de Kaplan-Meier para la funci´on de supervivencia se obtiene a trav´es del paquete estad´ıstico survival (cargado anteriormente) mediante la funci´on survfit(). Para ello, consideramos los estratos (clinica A y B) para poder compararlos. 66 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R outp<-survfit(Surv(tiempo,censura)~strata(clinica), type="kaplan-meier",data=cdat) attach(outp) Veamos un resumen estad´ıstico, mostrando los 10 primeros t´erminos de cada estrato summary(outp) ## Call: survfit(formula = Surv(tiempo, censura) ~ strata(clinica), ## data = cdat, type = "kaplan-meier") ## ## ## strata(clinica)=clinica=0 ## time n.risk n.event survival std.err lower 95% CI upper 95% CI ## 7 162 1 0.9938 0.00615 0.98184 1.000 ## 17 161 1 0.9877 0.00868 0.97080 1.000 ## 19 160 1 0.9815 0.01059 0.96094 1.000 ## 29 157 1 0.9752 0.01223 0.95155 0.999 ## 30 156 1 0.9690 0.01366 0.94258 0.996 ## 33 155 1 0.9627 0.01493 0.93390 0.992 ## 35 154 1 0.9565 0.01609 0.92545 0.989 ## 37 153 1 0.9502 0.01716 0.91719 0.984 ## 41 152 1 0.9440 0.01815 0.90907 0.980 ## 47 151 1 0.9377 0.01907 0.90107 0.976 ## ## strata(clinica)=clinica=1 ## time n.risk n.event survival std.err lower 95% CI upper 95% CI ## 13 74 1 0.986 0.0134 0.961 1.000 ## 26 73 1 0.973 0.0189 0.937 1.000 ## 35 72 1 0.959 0.0229 0.916 1.000 ## 41 71 1 0.946 0.0263 0.896 0.999 ## 79 68 1 0.932 0.0294 0.876 0.991 ## 109 66 1 0.918 0.0321 0.857 0.983 ## 122 65 1 0.904 0.0346 0.838 0.974 ## 143 64 1 0.890 0.0368 0.820 0.965 ## 149 62 1 0.875 0.0389 0.802 0.955 ## 161 61 1 0.861 0.0408 0.785 0.945 Esta salida devuelve los siguientes valores: time: tiempo de supervivencia de cada cliente dado en d´ıas. 67 n.risk: n´umero de elementos en riesgo en ese instante. n.event: n´umero de elementos que fallan en ese momento. survival: es la estimaci´on de Kaplan-Meier de la funci´on de supervivencia en el instante correspondiente, ˆ S(t). std.err: es el error est´andar asociado a cada ˆ S(t). lower 95 % CI: es el extremo inferior del intervalo de confianza para S(t) al nivel 95 %. upper 95 % CI: es el extremo superior del intervalo de confianza para S(t) al nivel 95 %. Representamos la estimaci´on calculada de S(t). Para ello, usaremos la funci´on ggsurvplot() contenida en el paquete survminer, (Kasambara y Kosinski [14]). install.packages("survminer") library(survminer) ggsurvplot(fit = outp, data = cdat, conf.int =T,title ="Curva de Supervivencia",xlab ="Tiempo",ylab ="Probabilidad de supervivencia") En esta figura presentamos las gr´aficas de las funciones de supervivencia estimadas, con las bandas de confianza asociadas por el m´etodo de Kaplan-Meier, para los dos tipos de cl´ınicas. Notamos que las curvas estimadas son de tipo escalonada, ambas con datos censurados (prueba de ello es la aparici´on de s´ımbolos “+” en cada una de ellas). Observamos diferencias en las curvas de supervivencia de las dos cl´ınicas. En general, ˆ S(t) es mayor para la cl´ınica B (codificada con el valor 1). 68 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R ++ + ++++++ + +++ + +++ + ++ ++ ++++ ++++++++ + ++ +++ +++++ + + + ++++ +++++++++++++ ++ ++ ++++ ++ +++ + + 0.00 0.25 0.50 0.75 1.00 0 250 500 750 1000 Tiempo Probabilidad de supervivencia Strata ++ strata(clinica)=clinica=0 strata(clinica)=clinica=1 Curva de Supervivencia Figura 4.1: Estimador de Kaplan-Meier. Estimador de Nelson-Aalen Recordar que, H(t) es la funci´on de hazard acumulada (o cumulative hazard) definida como: H(t) = Zt 0 h(u)du que representa la suma de las probabilidades de fallo (en nuestro caso, de abandonar el tratamiento) en el intervalo (0, t]. Matem´aticamente, su f´ormula viene dada por: ˆ H(t) = X j:t(j)≤t dj nj Instalamos los siguientes paquetes que nos har´an falta para calcular este estimador. install.packages("ggfortify") library(ggfortify) 69 install.packages("dplyr") library(dplyr) En este caso, la informaci´on se extrae calculando la suma acumulada del n´umero de eventos entre las personas en riesgo como se vi´o en el Cap´ıtulo 2. mod<- survfit(Surv(tiempo,censura)~clinica, cdat) R<- mod %>% fortify %>% group_by(strata) %>% mutate(CumHaz = cumsum(n.event/n.risk)) R ## # A tibble: 215 x 10 ## # Groups: strata [2] ## time n.risk n.event n.censor surv std.err upper lower strata ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> ## 1 2 163 0 1 1 0 1 1 1 ## 2 7 162 1 0 0.994 0.00619 1 0.982 1 ## 3 17 161 1 0 0.988 0.00878 1 0.971 1 ## 4 19 160 1 0 0.981 0.0108 1 0.961 1 ## 5 28 159 0 2 0.981 0.0108 1 0.961 1 ## 6 29 157 1 0 0.975 0.0125 0.999 0.952 1 ## 7 30 156 1 0 0.969 0.0141 0.996 0.943 1 ## 8 33 155 1 0 0.963 0.0155 0.992 0.934 1 ## 9 35 154 1 0 0.956 0.0168 0.989 0.925 1 ## 10 37 153 1 0 0.950 0.0181 0.984 0.917 1 ## # ... with 205 more rows ## ## ## # A tibble: 215 x 10 ## # Groups: strata [2] ## CumHaz ## <dbl> ## 0 ## 0.00617 ## 0.0124 ## 0.0186 ## 0.0186 ## 0.0250 ## 0.0314 ## 0.0379 ## 0.0444 ## 0.0509 ## # ... with 205 more rows 70 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R Cuya gr´afica asociada es la siguiente: 0 1 2 3 0 300 600 900 Tiempo (Dias) Riesgo Acumulado strata 0 1 Riesgo Acumulado Figura 4.2: Estimador de Nelson-Aalen. De nuevo observamos bastantes diferencias entre las cl´ınicas. Se aprecia que hay un menor ˆ H(t) para la cl´ınica B (codificada con el valor 1) lo que es coherente con los resultados obtenidos al aplicar Kaplan-Meier puesto que ˆ H(t) = −ln ˆ S(t). Para formalizar estas apreciaciones realizamos, a continuaci´on, un test de hip´otesis. Test de Log-Rank Para aplicar este test, definimos el siguiente contraste de hip´otesis: H0:S1(t) = S2(t) H1:S1(t)6=S2(t) Para comparar ambas funciones de supervivencia, aplicamos el test de LogRank: 71 out1<-survdiff(Surv(tiempo, censura) ~clinica,data=cdat) out1 ## Call: ## survdiff(formula = Surv(tiempo, censura) ~ clinica, data = cdat) ## ## N Observed Expected (O-E)^2/E (O-E)^2/V ## clinica=0 163 122 90.9 10.6 27.9 ## clinica=1 75 28 59.1 16.4 27.9 ## ## Chisq= 27.9 on 1 degrees of freedom, p= 1.28e-07 En base a los valores obtenidos en la tabla anterior, el p-valor p= 1.28·10−07 < 0.05 luego se rechaza la hip´otesis nula de igualdad de funciones de supervivencia (para un nivel de significaci´on del 5 %). En consecuencia, podemos concluir que existe una clara evidencia de desigualdad entre las curvas de supervivencia. Adem´as, X2 1= 27.9. Como son distintas, los datos generados permiten a su vez realizar una estimaci´on del riesgo hr =O0/E0 O1/E1 donde: O0 representa el valor “Oberved” en la cl´ınica A (valor 0) E0 representa el valor “Expected” en la cl´ınica A (valor 0) O1 representa el valor “Oberved” en la cl´ınica B (valor 1) E1 representa el valor “Oberved” en la cl´ınica B (valor 1) luego hr<-(122/90.9)/(28/59.1) hr ## [1] 2.832862 Por tanto, los clientes que est´an en la cl´ınica A se mantienen en el tratamiento 2.832862 veces m´as que los de la cl´ınica B. 72 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R Ajuste del modelo Cox Consideramos un modelo de Cox con las covariables prision ydosis, usando la covariable clinica como estrato para diferenciar el comportamiento entre cl´ınicas. El modelo es de la forma: hm(t;X) = h0m(t)exp(β0X), m = 1,2. donde h1(t;X) referencia a la cl´ınica A y h2(t;X) referencia a la cl´ınica B. Hay que tener en cuenta que observando dicha expresi´on, los coeficientes correspondientes a cada covariable β= (β1, β2) son iguales en ambos modelos. Lo que var´ıa respecto a las cl´ınicas, es la funci´on baseline hazard h0m(t). Creamos el modelo de Cox con prision ydosis como covariables y clinica como estrato, d´andole el nombre mod(). mod<-coxph(Surv(tiempo,censura)~prision+dosis+strata(clinica), data=cdat, method="breslow") Para ver si el modelo es correcto, veamos si se cumple la hip´otesis de riesgos proporcionales. Al hacerlo, podemos obtener uno de estos dos casos: 1. Si las curvas se cruzan =⇒se rechaza la hip´otesis de riesgos proporcionales y el modelo resultante es con la covariable cl´ınica como estrato. 2. Si las curvas son paralelas =⇒se acepta la hip´otesis de riesgos proporcionales y el modelo resultante es introduciendo cl´ınica como covariable. Para ello, calculamos las funciones baseline hazard para cada tiempo, mostrando los 10 primeros t´erminos para cada estrato: bh <- basehaz(mod, centered=TRUE) bh ## hazard time strata ## 1 0.000000000 2 clinica=0 ## 2 0.005302894 7 clinica=0 ## 3 0.010677622 17 clinica=0 ## 4 0.016126156 19 clinica=0 ## 5 0.016126156 28 clinica=0 ## 6 0.021724919 29 clinica=0 ## 7 0.027363076 30 clinica=0 ## 8 0.033028253 33 clinica=0 ## 9 0.038733766 35 clinica=0 73 ## 10 0.044466952 37 clinica=0 ## ## ## 148 0.000000000 2 clinica=1 ## 149 0.012484863 13 clinica=1 ## 150 0.025167324 26 clinica=1 ## 151 0.038130669 35 clinica=1 ## 152 0.051531857 41 clinica=1 ## 153 0.051531857 53 clinica=1 ## 154 0.051531857 72 clinica=1 ## 155 0.065994531 79 clinica=1 ## 156 0.065994531 86 clinica=1 ## 157 0.081399362 109 clinica=1 Calculamos los logaritmos: lnhazard<-log(bh[,1]) lntime<-log(bh[,2]) Veamos si son paralelos o no, para ello instalamos el paquete, que se necesita para dibujarlo: install.packages("lattice") library(lattice) xyplot(lnhazard~lntime, group=strata,auto.key =TRUE,data=bh, xlab ="ln(t)",ylab ="ln(-ln(S(t))",main ="Clinica como estrato",type="l") 80 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R devresi <- resid(mod2, type="deviance") plot(mod2$linear.predictor, devresi, ylab="Residuos de Deviance", main='Residuos de deviance') abline(h=0,lty=2,col='black') −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 −2 −1 0 1 2 Residuos de deviance mod2$linear.predictor Residuos de Deviance Figura 4.6: Salida de los residuos del estad´ıstico Deviance. Como se comprueba, no apreciamos patrones definidos pero si valores alejados del origen. Residuos escalados de Schoenfeld Ahora estamos interesados en evaluar la hip´otesis de riesgos proporcionales del modelo de CPH, examinando si el impacto de una o m´as covariables sobre el riesgo de los adictos a la hero´ına puede variar con el tiempo. Calculamos los residuos escalados de Schoenfeld para nuestro caso de la forma: par(mfrow=c(1,2)) plot(cox.zph(mod2)) 81 Time Beta(t) for prision 45 220 470 740 −2 −1 0 1 2 3 4 Time Beta(t) for dosis 45 220 470 740 −0.2 −0.1 0.0 0.1 0.2 Figura 4.7: Salida de los residuos escalados de Schoenfeld. Las tendencias en los diagramas de dispersi´on de los residuos escalados de Schoenfeld a menudo son dif´ıciles de determinar, especialmente con las covariables binarias (como prision en nuestro caso) donde s´olo hay dos bandas horizontales de residuos presentes. Residuos dfbeta dfbeta <- residuals(mod2, type="dfbetas") par(mfrow=c(1,2)) for (j in 1:2){ plot(dfbeta[,j], ylab=names(coef(mod2))[j]) abline(h=0,lty=2,col='black') lines(c(0,0),c(0,0)) } 82 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R 0 50 100 150 200 −0.2 −0.1 0.0 0.1 0.2 Index prision 0 50 100 150 200 −0.2 −0.1 0.0 0.1 0.2 0.3 Index dosis Figura 4.8: Salida de los residuos dfbeta. En estas figuras se nos muestran los residuos dfbeta del modelo. Como vemos estos residuos est´an centrados con respecto al origen, y no presentan patrones definidos. Se nos presentan algunos datos demasiados alejados del origen en ambas figuras. 4.1. Conclusiones al estudio de los datos de metadona Los datos analizados, corresponden a programas de rehabilitaci´on realizados en distintas cl´ınicas australianas durante los a˜nos 1970. Individuos adictos a la hero´ına fueron sometidos a tratamientos con metadona. Trabajamos con los datos del estudio realizado por Caplehorn y Bell (1991), cuyo objetivo era verificar si exist´ıan diferencias significativas entre el tipo de tratamiento aplicado a los pacientes en las etiquetadas como cl´ınica A y B (la pol´ıtica de estas cl´ınicas fue diferente) . Se consideran como covariables de inter´es en el estudio: dosis (dosis diaria de metadona administrada durante el tratamiento), y prisi´on (si el individuo ha estado antes en prisi´on o no). 4.1. CONCLUSIONES AL ESTUDIO DE LOS DATOS DE METADONA83 Aplicando la metodolog´ıa propuesta en el trabajo hemos obtenido los siguientes resultados: 1. Se han estimado las funciones de supervivencia, S(t), y de hazard acumulada, H(t), en cada cl´ınica. 2. Hemos enconttrado que existen diferencias significativas entre las funciones de supervivencia en las dos cl´ınicas, aplicando el test de log-rank. 3. Se ha discutido por qu´e es adecuado utilizar un modelo de Cox estratificado en este caso. 4. El modelo estratificado de Cox se ha aplicado considerando como estratos las cl´ınicas A y B. Esto nos ha permitio estimar las funciones basales para cada cl´ınica, y estimar el efecto de las covariables (dosis y prisi´on) en el tiempo de permanencia en tratamiento. Se incluye el an´alisis de los residuos, para ilustrar c´omo se calculan estos en el modelo de Cox. 84 CAP´ ITULO 4. APLICACI ´ ON PR ´ ACTICA CON EL SOFTWARE R Bibliograf´ıa [1] Boj del Val,Eva. “El modelo de regresi´on de Cox”. Departamento de Matem´atica Econ´omica, Financiera y Actuarial de Mayo de 2017 (Universidad de Barcelona). [2] Borges P.,Rafael Eduardo. “An´alisis de supervivencia de pacientes con di´alisis peritoneal”. Revista Colombiana de Estad´ıstica de Diciembre de 2005 (Volumen 28 No 2. pp. 243 a 259). [3] Caplehorn, J.R.M. and Bell. “Methadone dosage and retention of clients in maintenance treatment”. Med.J.Aust. 154, pp. 195-199. (1991). [4] Cardona Hurtado,Diego Alejandro and Trujillo Bonilla,Jenny Carolina. “Aspectos b´asicos de estimaci´on no param´etrica en an´alisis de sobrevivencia. Una aplicaci´on a un estudio de deserci´on estudiantil”. Trabajo para optar al t´ıtulo de Profesional en Matem´aticas con ´ Enfasis en Estad´ıstica de 2013 (Universidad Del Tolima) Ibagu´e, Colombia. [5] Caroni,Chrys. “Lifetime data analysis. Reliability and survival analysis”. (2002). [6] Caroni,Chrys. “The Correct “Ball Bearings” Data”. Lifetime Data Analysis, 8, 395-399. Kluwer Academic Publishers. (2002). [7] Garc´ ıa-Hinojosa,Cristina Pruenza. “Estudio de an´alisis de supervivencia”. Trabajo Fin de Grado de Mayo 2014 (Universidad Aut´onoma de Madrid). [8] Grambsch, P. and Therneau, T.M.. “Proportional hazards tests and diagnostics based on weighted residuals”. Biometrika. 81,515-26. (1994). [9] Brostr¨ om, G¨ oran. “eha: Event History Analysis”. R package version 2.5.1.(2017). [10] Hern´ andez Dom´ ınguez,Ana Mar´ ıa. “An´alisis estad´ıstico de datos de tiempos de fallo en R”. M´aster Universitario en Estad´ıstica Aplicada de 2010 (Universidad de Granada). 85 86 BIBLIOGRAF´ IA [11] Hess, K.R.. “Graphical Methods for assessing violations of the proportional hazards assumption in Cox regression”. Statistics in Medicine, 14, 1707-1723. (1995). [12] Jim´ enez,Pablo Moreno. “Inferencia estad´ıstica para datos censurados. M´etodos y aplicaciones”. Trabajo Fin de Grado de Junio de 2014 (Universidad de Sevilla). [13] Kalbfleisch, John D. and Prentice, Ross L. “The Statistical Analysis of Failure Time Data”. Second Edition. Wiley & Sons. (2002). [14] Kassambara, Alboukadel and Kosinski, Marcin. “survminer: Drawing Survival Curves using ’ggplot2’ ”. (2018). [15] Lawless, J.F. ”Statistical Models and Methods for Lifetime Data”. Second Edition. Wiley & Sons. (2003). [16] L´ opez Montoya,Antonio Jes´ us. “Comparaci´on de dos modelos de Regresi´on en fiabilidad”. M´aster Universitario en Estad´ıstica Aplicada de 2011 (Universidad de Granada). [17] Martinez,Javier. “An´alisis de Supervivencia en R”. 22 de mayo de 2017. [18] Quintanilla Casas,Beatriz. “Estad´ıstica en variables con censura: Aplicaci´on a datos medioambientales”. M´aster en Bioinform´atica y Bioestad´ıstica de Junio de 2017 (Universitat Oberta de Catalunya). [19] R Core Team, “R: A Language and Environment for Statistical Computing”. R Foundation for Statistical Computing. Vienna, Austria. (2018). [20] Roque Roque,Daniel Octavio, “Forma funcional de covariables en el modelo de Cox”. Tesis de 2009 (Universidad Nacional Mayor de San Marcos) Lima-Per´u. [21] Sarkar, Deepayan. “lattice: Multivariate Data Visualization with R”. Springer. New York. (2008). [22] Tang, Yuan; Horikoshi, Masaaki and Li, Wenxuan “ggfortify: Unified Interface to Visualize Statistical Result of Popular R Packages”. (2016). [23] Therneau, Terry M. and Grambsch, Patricia M.. “survival: Modeling Survival Data: Extending the Cox Model”. Springer, New York.(2000). [24] Wickham, Hadley ; Franc¸ois, Romain ; Henry, Lionel and M¨ uller, Kirill. “dplyr: A Grammar of Data Manipulation”. (2018).