Supervivencia

MODELOS DE SUPERVIVENCIA

El análisis de supervivencia estudia el tiempo hasta que ocurre un evento y los factores que influyen en los mismos. Podríamos decir que es una mezcla de regresión lineal y logística.

Algunos ejemplos de tiempos de supervivencia (o tiempo hasta el evento): + tiempo desde el nacimiento hasta el fallecimiento + tiempo desde que se entra en un ensayo clínico hasta el abandono o fallecimiento + tiempo desde que se fabrica una pieza hasta que se estropea

En la regresión logística, yo estoy fijando el tiempo. Pero el analisis de supervivencia lo que nos permite es modelizar el tiempo que pasa hasta que ocurre un evento. Entonces, yo estoy teniendo en cuenta dos cosas. Por un lado el hecho de que ocurra el evento, y por otro lado, el tiempo que ha pasado hasta que ese evento ocurre.

La variable respuesta es el tiempo, con lo cual es una variable continua.(parecido a lo que nos pasaba en poisson, solo que en poisson lavariable respuesta en una variable discreta),y clarisimamente no negativa.

Una característica clave del análisis de supervivencia es que la variable respuesta es una variable aleatoria no-negativa con distribución continua o discreta. Esta variable define el tiempo que transcurre desde el origen (bien definido) hasta que ocurre el evento (bien definido).

La segunda característica importante es la censura, que se da cuando los momentos en los que ocurren el evento origen o evento final no son observados.

La censura más común es la censura por la derecha, lo que implica que sólo se sabe que el evento de darse, se ha dado a partir de un cierto momento. Matemáticamente sea $T^∗$el tiempo hasta el evento, y U una variable aleatoria que representa el tiempo hasta un evento censura, lo que observamos es $T = min(T^∗,U)$ y un indicador de censura $δ = I[T^∗ < U$]. Es decir, δ toma los valores 0 y 1, en función de si el tiempo observado está censurado o si es un evento observado, respectivamente.

La censura se puede clasificar en tres tipos: tipo I, tipo II y aleatoria.

Tipo I:

los tiempos de censura están preestablecidos. Por ejemplo, en un experimento animal, una cohorte de animales puede comenzar en un momento específico, y todos son seguidos hasta una hora final preestablecida. Animales que no han experimentado el evento de interés antes del final del estudio se censuran en ese momento.

Tipo 2:

ocurre cuando los objetos experimentales son seguidos hasta que una fracción preespecificada ha fallado. Este tipo de diseño es poco frecuente en los estudios biomédicos, pero puede ser utilizado en entornos industriales, donde el tiempo hasta el fallo de un dispositivo es de interés primordial. Un ejemplo sería aquel en el que el estudio se detiene después de que, por ejemplo, se observa que 25 de 100 dispositivos han fallado. Los 75 dispositivos restantes serían censurados. En este ejemplo, se observa el 25 % más pequeño de los tiempos de fallo ordenados, y el resto son censurados.

Aleatorio:

se debe tener especial cuidado con este tipo de censura para evitar resultados sesgados. Una de las causas de la censura aleatoria es el abandono de los pacientes del estudio. Si el abandono ocurre verdaderamente al azar, y no está relacionado con el proceso de la enfermedad, tal censura podría no causar ningún problema de sesgo en el análisis. Pero si los pacientes que están cerca de la muerte son más propensos a abandonar el estudio que otros pacientes, pueden surgir sesgos graves. Nosotros trabajaremos asumiendo que las censuras son independientes de los tiempos de supervivencia, lo cual tendremos que comprobar previamente.

FUNCIÓN DE SUPERVIVENCIA.

No podemos entender un modelo de supervivencia si no podemos entender las funciones de probabilidad que tenemos por detrás del modelo de supervivencia.

Los métodos de análisis de supervivencia dependen de la distribución de la supervivencia, y las dos maneras clave de especificarlo son mediante la función de supervivencia y la función de riesgo (o hazard function).

La función de supervivencia define la probabilidad de sobrevivir hasta el punto t. Formalmente, \[S(t) = P(T > t) \hspace{1em} para \hspace{1em} 0 < t < ∞.\] Es decir, que probabilidad tiene un individuo de haber sobrevivido hasta ese tiempo.Es decir, que después de ese tiempo este vivo, aunque fallezca justo en ese instante después.

Esta función toma el valor 1 en el tiempo 0, disminuye (o permanece constante) en el tiempo, y por supuesto nunca cae por debajo de 0.(Puede disminuir, o puede permanecer constante). La integral me tiene que dar 1.

La probabilidad de que una variable continua tome un valro concreto es cero. imposible.

FUNCIÓN DE RIESGO

La función de supervivencia se define a menudo en términos de la función de riesgo o función hazard, la cual es la tasa de eventos instantáneos. Es decir, cual es la probabilidad de que en un tiempo phi haya un número de eventos dividido entre lo que dura ese tiempo, en una fracción infinitesimal de ese tiempo.

Es la probabilidad de que, dado que un sujeto ha sobrevivido hasta el tiempo t, tiene el evento (fallece) en el siguiente pequeño intervalo de tiempo, dividido entre la duración de ese intervalo. Formalmente, esto puede expresarse como \[h(t) = \lim_{φ\to 0} \frac{P(t < T < t + φ|T > t)}{φ}\] Las funciones de riesgo y supervivencia son dos maneras de especificar la distribución de la supervivencia.

FUNCIÓN DE RIESGO ACUMULADA

Además de las funciones de supervivencia y de riesgo, existen otras maneras de definir la supervivencia. La función de distribución acumulada (FDA), que es comúnmente utilizada en el área de teoría de la probabilidad y sus aplicaciones, viene dada por, \[F(t) = P(T ≤ t), \hspace{1em} 0 < t < ∞\]

Es decir, el complementario de la función de supervivencia. Y, al igual que esta, es continua por la derecha. En el análisis de supervivencia, esta función se conoce como la función de riesgo acumulada (no confundir con función hazard acumulada que se define más adelante).

RELACIÓN ENTRE LA FUNCIÓN DE DENSIDAD Y LA DE SUPERVIVENCIA

Una vez tenemos la función acumulada, podemos definir la función de densidad.( la derivada de la función acumulada) La función de densidad se define como: \[f(t) = −\frac{d}{dt} S(t) = \frac{d}{dt} F(t)\] La función hazard h(t) es precisamente el ratio entre la función de densidad y la función de supervivencia: \[h(t) = \frac{f (t)}{S(t)}\] Es decir, la hazard en el tiempo t es la probabilidad de que un evento ocurra en el vecindario del tiempo t dividido por la probabilidad de que el sujeto esté vivo en el tiempo t.

FUNCIÓN HAZARD ACUMULADA

La función hazard acumulada se define como el área bajo la función hazard h(t) descrita anteriormente, \[H(t) = \int_{0}^{t}h(u)du\] Y la función de supervivencia S(t) se puede escribir en términos de esta: \[S(t) = exp(−\int_{0}^{t}h(u)du)= exp(−H(t))\]

MEDIA DE LOS TIEMPOS DE SUPERVIVENCIA

La supervivencia media es la esperanza matemática de la variable aleatoria tiempo de supervivencia (T), \[µ = E(T) = \int_{0}^{∞}tf (t)dt\] lo que también se puede escribir como \[µ =\int_{0}^{∞}S(t)dt\] NOTA: El tiempo medio de supervivencia sólo está definido si S(∞) = 0, es decir, si todos los sujetos eventualmente tienen el evento de interés.

Este podría no ser el caso si, por ejemplo, el interés está en el tiempo hasta la recurrencia de un cáncer, y alguna fracción c de los sujetos se curan y por lo tanto no tienen recurrencia.

En ese caso, S(∞) = c y el área bajo la curva de supervivencia es infinita. Si se requiere un tiempo medio de supervivencia en esta situación, una solución es especificar un tiempo de supervivencia máximo posible, de modo que la integral se vuelva finita.

DISTRIBUCIÓN DE SUPERVIVENCIA PARAMÉTRICA: EXPONENCIAL

Existen varias distribuciones de supervivencia disponibles para modelar los datos de supervivencia.

La distribución exponencial, la distribución de supervivencia más sencilla, tiene una función hazard constante, $h(t) = λ$. Y por lo tanto:

La función hazard acumulada viene dada por: $H(t) = λt$ (se obtiene facilmente integrando la función hazard $h(t)$).
La función de supervivencia viene dada por: $S(t) = e^{−H(t)} = e^{−λt}$
La función de densidad es: $f (t) = h(t)S(t) = λe^{−λt}$

Es fácil trabajar con la distribución exponencial, pero la asunción de riesgo constante a menudo no es apropiado para describir la vida de los seres humanos o animales.

DISTRIBUCIÓN DE SUPERVIVENCIA PARAMÉTRICA: WEIBULL

La distribución Weibull por su parte ofrece mayor flexibilidad para modelizar datos de supervivencia. + función hazard: $h(t) = αλ(λt)^{α−1} = αλ^{α}t^{α−1}$. + función hazard acumulada: $H(t) = (λt)α$. + función de supervivencia: $S(t) = e^{−(λt)^{α}}$.

ESTIMACIÓN NO PARAMÉTRICA

ESTIMACIÓN NO PARAMÉTRICA DE LA FUNCIÓN DE SUPERVIVENCIA

Hemos visto que hay una gran variedad de formas entre las que elegir para la función hazard si uno quiere modelar unos datos utilizando un modelo paramétrico. Pero, ¿cuáles son los parámetros que debería utilizarse para una aplicación en particular?¿y cuál es esa distribucion en concreto a la que yo me tendría que ajustar?

Cuando se modelan datos reales sobre seres vivos por ejemplo, es difícil saber qué familia paramétrica elegir, y a menudo ninguna de las familias disponibles tiene la flexibilidad suficiente para modelar la forma real de la distribución. Así, en aplicaciones médicas y sanitarias, utilizar métodos no paramétricos tienen ventajas considerables.

ESTIMADOR KAPLAN MEIER

El estimador Kaplan Meier es el estimador no-paramétrico más común, sencillo e intuitivo para estimar la función de supervivencia, y viene dado por la siguiente expresión: \[\hat{S}(t) = \prod_{t_{k}≤t}(1 − \hat{q}_{k} ) = \prod_{t_{k}≤t}\left(1 −\frac{dk}{nk}\right)\] donde $n_k$ es el número de individuos en riesgo en el tiempo $t_k$ y $d_k$ es el número de individios para los que se ha observado el evento en ese tiempo.

EJEMPLO

library(survival)# cargamos la librería survival

## Warning: package 'survival' was built under R version 4.0.3

tt <- c(7,6,6,5,2,4)# indicamos para cada uno de los individuos su tiempo de supervivencia
cens <- c(0,1,0,0,1,1)# y le indicamos si hemos observado el evento o no lo hemos observado. Es decir, si sabemos que es un evento o si sabemos que es una censura. 0= censura, 1=evento.
Surv(tt,cens)# tenemos que crear un objeto con la función surv, al cual le tengo que pasar el tiempo de supervivencia, y el vector de censura asociado, el orden es importante.

## [1] 7+ 6  6+ 5+ 2  4

Para hacer el etimador kaplan meier lo que utilizamos es la función survfit que basicamente nos ajusta una función de supervivencia kaplan meier.

result.km <- survfit(Surv(tt,cens)~1, conf.type="log-log") ## como ajusto sin covariables, gusanito 1.
summary(result.km)

## Call: survfit(formula = Surv(tt, cens) ~ 1, conf.type = "log-log")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     2      6       1    0.833   0.152       0.2731        0.975
##     4      5       1    0.667   0.192       0.1946        0.904
##     6      3       1    0.444   0.222       0.0662        0.785

Podemos calcular el intervalo de confianza para esa curva, para ese estimador. Tenemos diferentes opciones para ese intervalo de confianza pero el que más se utiliza, es el log-log, y es por que ajustamos el logaritmo del logaritmo.

Intervalo de confianza para el etimador K.M(I de II)

En primer lugar estimamos la varianza de log(Sˆ(t)) mediante el método Delta, \[var(log\hat{S}(t_k)) = \sum_{t_i≤t}var(log(1 − \hat{q}_i)) ≈ \sum_{t_i≤t}\frac{di}{n_i(n_i − di)}\] Y para obtener la varianza de $var(Sˆ(t))$, utilizamos de nuevo el método Delta \[var(\hat{S}(t)) ≈ [\hat{S}(t)]^2\sum_{t_i≤t}\frac{di}{n_i(n_i − di)}\]

Intervalo de confianza para el etimador K.M(II de II)

Sin embargo, los intervalos de confianza basados en esta varianza pueden exceder de 1 o estar por debajo de 0. Por lo tanto, una aproximación más satisfactoria es obtener intervalos de confianza para la transformación complementaria log-log de $\hat{S}(t)$, \[var(log[−log\hat{S}(t)] ≈\frac{1}{[log\hat{S}(t)]^2}\sum_{t_i≤t}\frac{d_i}{n_i(n_i − d_i)}\] Volvémos al ejemplo:

result.km <- survfit(Surv(tt,cens)~1, conf.type="log-log") ## como ajusto sin covariables, gusanito 1.
summary(result.km)

## Call: survfit(formula = Surv(tt, cens) ~ 1, conf.type = "log-log")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     2      6       1    0.833   0.152       0.2731        0.975
##     4      5       1    0.667   0.192       0.1946        0.904
##     6      3       1    0.444   0.222       0.0662        0.785

Calculamos la función de supervivencia, y lo que tenemos aqui es el error estándar asociado a cada uno de ellos. El error estándar aumenta, por que mi tamaño muestral disminuye. A medida que me voy quedando en el final, tengo menos individuos, aumeta mi error estándar y aumenta mi intervalo de confianza.

Y de ahi obtenemos la función kaplan-meier.

plot(result.km)

En la práctica, la función Kaplan-Meier se representa como una función escalonada. El tiempo medio de supervivencia está en t = 6, que es es el tiempo t más pequeño tal que S(t) ≤ 0.5.

MODELO DE SUPERVIVENCIA DE COX

comparación de dos distribuciones de supervivencia

Supongamos que queremos comparar dos funciones de supervivencia $S_1(t)$ y $S_0(t)$. ¿Qué podemos hacer? \[\left\{ \begin{array}{ll} H_0: & \mbox{} S_1(t) = S_0(t)\\ H_1: & \mbox{} ??? \end{array} \right.\] Una alternativa es introducir un parámetro $ψ$ de forma que $S_1(t) = [S_0(t)]^ψ$, o equivalentemente, $h_1(t) = ψh_0(t)$ (recordad que $S(t) = exp(H(t))$, y al mismo $h(t) =\frac{dH(t)}{dt} =\frac{S^\prime(t)}{S(t)}$$. De esta forma, el contraste se convertiría en:

\[\left\{ \begin{array}{ll} H_0: & \mbox{} ψ = 1 \\ H_1: & \mbox{} ψ \neq 1 \end{array} \right.\]

Si toma un valor 1, esas funciones van a ser iguales; y si es distinto de 1, esas funciones van a ser distintas.

Vamos a utilizar este concepto como base para introducir los modelos de regresión de riesgos proporcionales. Yo asumo que la relacion entre dos funciones de riesgo es siempre esa constante ψ, no depende del tiempo. Para más información sobre comparación de distribuciones de supervivencia ver el Capítulo 4 de Moore, 2006.

La ecuación $h_1(t) = ψh_0(t)$ es la clave para cuantificar la diferencia entre dos funciones de riesgo (dos funciones hazard), y el modelo de riesgos proporcionales es ampliamente utilizado en este caso. (Más adelante veremos cómo evaluar la validez de esta suposición). Además tambien nos interesa modelizar esas varaible de interés en función de otros parámetros. Osea, yo queiro saber si el sexo cambia mi supervivencia, si la edad cambia mi supervivencia…

Además, podemos ampliar el modelo para incluir un vector de q covariables Z de la siguiente manera: \[ψ = e^{Zβ}\].

Eso quiere decir que yo puedo plantearme que ese parametro phi(yo tengo dos funciones de supervivencia que dependen del tiempo, que las estoy relacionando con un parámetro, pero es que ahora lo que voy a decir es que ese parámetro está relacionado con una serie de covariables que yo le voy a meter al modelo).

A priori, y en los modelos que nosotros vamosa ver, esas covariables no dependen del tiempo.

Es posible relacionar $ψ$ y el vector de covariables $Z$ con otras relaciones funcionales, sin embargo, esta es la relación más común en la práctica.

El modelo de riesgos proporcionales nos permitirá ajustar los modelos de regresión a datos de supervivencia con censura, disponiendo así de modelos muy similares a los que obtenemos con la regresión lineal y la regresión logística. Sin embargo, no asumir una forma paramétrica particular para $h_0(t)$, junto con la presencia de censura, hace que la modelización de datos de supervivencia sea particularmente complicada.

Para ello vamos a introducir la verosimilitud parcial, desarrollada inicialmente por D.R. Cox, y, por lo tanto, a menudo se le conoce como el modelo de riesgos proporcionales de Cox.

Estimación delModelo de Cox de R.P

Consideramos un modelo donde la covariable Z representa el grupo de un ensayo clínico (casos=1, controles=0).

Sea k el indicador del tiempo hasta el evento (ordenado).

La función hazard para el individuo $i$ en el tiempo de evento $t_k$ es $h_i(t_k)$. Bajo el modelo de riesgos proporcionales, podemos escribir esta función hazard como: \[h_i(t_k) = h_0(t_k)ψ_i = h_0(t_k )e^{βz_i}\] $h_0$ se le llamariesgo basal y no depende del individuo, lo que quiere decir que es fgenerico. es decir, que va a ser el mismo para todos. Es cómo e riesgo base, cómo la media. Lo mismo que cuando tenemos un modelo de regresion lo que representa la $\beta_0$, el intercept.

Entonces, para un indivuduo i,su función de riesgo en el tiempo $t_k$, la puedo representar en función del riesgo base y un parametro phi que si depende de ese individuo.Pero es que ademas depende de ese individuo en fucion de los valores que ese individuo tome en una serie de variables que van a ser las $Z_i$

Esperamos que los pacientes en el grupo experimental tengan menor riesgo de experimentar el evento que los pacientes en el grupo control, por lo tanto esperamos que $β < 0$, y consecuentemente $ψ < 1$. De forma que la hazard para ese paciente al que se le esta dando la vacuna tiene menos riesgo, porque su funcion de riesgo, como phi es menor que 1, es menor que la funcion de reisgo media para los que están en el control.

Sea $X = (X_1, . . . , X_q)$ un vector de q covariables, que pueden ser continuas, ordinarias, dicotómicas, etc.

El modelo de Cox de riesgos proporcionales se escribe como: \[h(t) = h_0(t)ψ = h_0(t)e^{\sum_{j=1}^{q}β_jX_j}\] Los parámetros $β = (β_1, . . . , β_q)$ se estiman maximizando el logaritmo de la denominada función de verosimilitud parcial. La maximización de dicha función se realiza mediante métodos numéricos (por lo que nos tenemos que asegurar de que obtenemos la convergencia), obteniendo de esta forma la estimación $\hat{β} = (\hat{β_1}, . . . , \hat{β_q})$.

Signifiación estadística de las variables del modelo

De forma similar a lo que ocurría en los modelos GLM, existen diferentes alternativas para estudiar la significación estadísticas de los parámetros $ , para $j = 1, . . . , q$ y/o las variables predictoras $X_j$. + Test de Wald + Test de la Razón de verosimilitud + Score Test

Test de Wall

Los coeficientes del modelo estimados $(\hat{β}_j)$ tienen una distribución normal aproximada cuando hay un número adecuado de eventos en la muestra. La aproximación normal se verifica mejor para las estimaciones de los coeficientes que los hazard ratio, por lo que los contrastes de hipótesis y los intervalos de confianza se construyen para los coeficientes y sus errores estándar. El test de Wald es quizás el test más utilizado, y se obtiene directamente en la salida de R.

El estadístico de contraste se construye de la misma manera que en un modelo GLM: \[z_p =\frac{\hat{β}_j}{s.e(\hat{β}_j)}, \hspace{2em} j = 1, . . . , q\] donde s.e representa el error estándar de la estimación del parámetro $_j. Rechazaremos $H_0 : β_j = 0$, si $2 ∗ P(Z > |z_p|) < α$ (equivalentemente, si $|z_p| > z_{α/2}$).

También podemos construir el intervalo de confianza al (1 − α)% como \[\hat{β}_j ± z_{α/2} · s.e(\hat{β}_j)\] ##### Teste de la zon de verosimilitud (univariante)

El test de la razón de verosimilitud utiliza el resultado de la teoría estadística de que $2(l(β = \hat{β}) − l(β = 0))$ sigue aproximadamente una distribución de chi-cuadrado con un grado de libertad. La ventaja clave de este test sobre el test de Wald es que es invariable a las transformaciones monótonas de $β$. Por ejemplo, el hecho de que el test se aplique en términos de $β$ o en términos de D $ψ = e^{β}$ no tiene ningún efecto en el p-valor para contrastar $H_0 : β = 0$.

Ejemplo: Tabaquismo

El propósito de este estudio fue evaluar la duración de una combinación de triple medicación versus terapia sólo con el parche de nicotina en fumadores con enfermedades médicas.
Los pacientes con un historial de tabaquismo fueron aleatoriamente asignados a la terapia de triple combinación o de parche y se les hizo un seguimiento a lo largo de seis meses.
La variable de resultado primaria fue el tiempo transcurrido desde la asignación al azar hasta la recaída (volver a fumar);
Las personas que permanecieron como no fumadores durante seis meses fueron censuradas.

Las variables de las que se dispone en la base de datos son las siguientes:

ttr: número de días sin fumar (tiempo hasta la recaida)
relapse: 1= recae y vuelve a fumar 0= censura
grp: tipo de tratamiento (patchonly vs combination)
employment: ft = jornada completa pt = jornada parcial

Objetivos principales:

comparar los dos tipos de tratamiento con respecto al tiempo hasta la recaida
estudiar otros factores que puedan influir en la variable respuesta

# Lectura de los datos
setwd("C:/Users/Alex/Desktop/MÁSTER MATEMÁTICAS/MODELIZACIÓN ESTADÍSTICA/IRANZU BARRIO")
tabaquismo <- read.table("tabaquismo.txt",header=TRUE)
summary(tabaquismo[,1:9])

##        id              ttr            relapse          grp           
##  Min.   :  1.00   Min.   :  0.00   Min.   :0.000   Length:125        
##  1st Qu.: 33.00   1st Qu.:  8.00   1st Qu.:0.000   Class :character  
##  Median : 67.00   Median : 49.00   Median :1.000   Mode  :character  
##  Mean   : 66.15   Mean   : 77.44   Mean   :0.712                     
##  3rd Qu.: 99.00   3rd Qu.:182.00   3rd Qu.:1.000                     
##  Max.   :130.00   Max.   :182.00   Max.   :1.000                     
##       age           gender              race            employment       
##  Min.   :22.00   Length:125         Length:125         Length:125        
##  1st Qu.:41.00   Class :character   Class :character   Class :character  
##  Median :49.00   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :48.84                                                           
##  3rd Qu.:56.00                                                           
##  Max.   :86.00                                                           
##   yearsSmoking  
##  Min.   : 9.00  
##  1st Qu.:22.00  
##  Median :30.00  
##  Mean   :30.88  
##  3rd Qu.:39.00  
##  Max.   :56.00

attach(tabaquismo)

Ajustamos el primer modelo univariante para la variable tipo de tratamiento

modelo.cox1<-coxph(Surv(ttr,relapse)~grp)## le tengo que decir el tiempo y la variable que me representa si ehe observado el tiempo o si es una censura. Muy importante: el 1 siempre es el evento, el 0 siemrpe es una censura. 
summary(modelo.cox1)

## Call:
## coxph(formula = Surv(ttr, relapse) ~ grp)
## 
##   n= 125, number of events= 89 
## 
##                coef exp(coef) se(coef)   z Pr(>|z|)   
## grppatchOnly 0.6050    1.8313   0.2161 2.8  0.00511 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## grppatchOnly     1.831     0.5461     1.199     2.797
## 
## Concordance= 0.581  (se = 0.027 )
## Likelihood ratio test= 7.99  on 1 df,   p=0.005
## Wald test            = 7.84  on 1 df,   p=0.005
## Score (logrank) test = 8.07  on 1 df,   p=0.004

El tratamiento es significativo.

el b1 es positivo, lo que quiere decir que noseque el mayor que uno. entonces, el riesgo de recaer que tienen los que tienen parche es mayor. Es decir, es mejor no usar el parche.

Y el riesgo es 1.8 veces mayor. (1:33H video)

Hoy vamos a empezar con la funcion de verosimilitud parcial, donde lo dejamos el otro día. Vamos a ver como ajustamos un modelo multiple y como seleccionar las variables que entran en el modelo.

Tenemos diferentes variables que podemos introducir, tenemos que ver cuales son significativas de forma univariante, en el modelo multiple tenemos la opcion de comparar modelos anidados mediante el test de la razon de verosimilitud, pero esta vez en vez de usar la función de verosimilitud, utilizará la funci´n de verosimilitud parcial, y vamos a acabar viendo cuales son las asunciones que se tienen que verificar en un modelo de cox, que las asumimos a priori antes de ajustar un modelo pero una vez que tenemos el modelo ajustado tenemos que comprobar que efectivamente se cumple. Vamos a ver cómo comprobamos que se cumplen, y en aquellos casos que no se cumplan vamos a discutir lo que podríamos hacer. y con esto terminamos los modelos de cox y supervivencia.

Función de verosimilitud parcial

La función de verosimilitud parcial, se denomina parcial debido a que tiene en cuenta únicamente en la función de verosimilitud las probabilidades de los tiempos de evento (muerte) y no incluye las probabilidades de los tiempos de datos censurados. Sin embargo, en el cálculo de las probabilidades de los tiempos de evento sí tiene en cuenta a todos los sujetos (censurados o no a posteriori) objeto de riesgo al inicio de los diferentes tiempos de evento.

Denominamos $L ≡ L(\hat{β})$, a la función de verosimilitud parcial. Supongamos que tenemos K tiempos de evento y que no hay empates. Así, tendremos n −K tiempos censurados. Los tiempos de evento ordenados los denotamos por t1, . . . ,tK .

Denominamos por $L_k ≡ L_{t_k}(\hat{β})$ para $k = 1, . . . ,K$ a las porciones de la verosimilitud total anterior debidas a la aportación de los diferentes tiempos de evento $t_1, . . . ,t_K$. Construiremos la función de verosimilitud total como el producto de cada una de las aportaciones de los K tiempos de evento: \[L(β) = \prod_{k=1}^{K}L_k(β)\] Pero, ¿cómo se calcula cada una de las aportaciones $L_k$ ?

Denotamos por $R(t_k)$ para $k = 1, . . . ,K$, al conjunto de los sujetos en riesgo en el tiempo $t_k$. \[L_k =\frac{h_0(t_k)e^{\sum_{j=1}^{q}β_jX_j^{(k)}}}{\sum_{l∈R(t_k)}h_0(t_k)e^{\sum_{j=1}^{q}β_jX^l_j}}=\frac{e^{\sum_{j=1}^{q}β_jX^{(k)}_j}}{\sum_{l∈R(tk)}e^{\sum_{j=1}^{q}β_jX^l_j}}\]

Si aplicamos el logaritmo sobre la verosimilitud parcial L(β), obtenemos la log-verosimilitud parcial: \[l(β) = \sum^{K}_{k=1}βX^{(k)} − \sum^{K}_{k=1}log \left(\sum_{l∈R(tk)} e^{βX^{(l)}} \right)\]

Derivando la log-verosimilitud parcial con respecto a los $β_j, j = 1, . . . , q$, obtenemos los estimadores del modelo.

Estimación Hazard basal

Una estimación de la función hazard basal viene dada por: \[h_0(t_k ) =\frac{dk}{\sum_{l∈R(t_k)}e^{\hat{β}X^{(l)}}}\] donde $d_k$ es el número de eventos en el tiempo $t_k$ . El modelo de Cox de riesgos proporcionales se llama modelo semi-paramétrico porque no asume una distribución conocida para la función hazard basal.

una vez que conocemos la VSparcial, conocemos, hemos estimado el hazard basal, lo siguiente es obtener los estimadores de los parametros para mi modelo, que seran aquellos que maximicen mi función de VS parcial, y luego tendre que ver si esos parámetros son significativos o no.(Test de Wald, Test de la Razón de verosimilitud, Score Test)

#####Test de la razón de verosimilitud (univariante)

En el test de la razon de VS recordad que lo que hacemos es ver la significación de mi variable.Independientemente de que sea una variable continua, sea una variable categórica con 2 categorías o con mas categorías.en la nula asumo que mi beta es igual a cero, y en la alternativa que mi beta es distinto de cero.

entonces lo que hacemos es comparar es la logVs del modelo con todos mis parametros con respecto a la hipótesis nula.

El test de la razón de verosimilitud utiliza el resultado de la teoría estadística de que $2 \left(l(β = \hat{β}) − l(β = 0)\right)$ sigue aproximadamente una distribución de chi-cuadrado con un grado de libertad. La ventaja clave de este test sobre el test de Wald es que es invariable a las transformaciones monótonas de $β$. Por ejemplo, el hecho de que el test se aplique en términos de $β$ o en términos de D $ψ = e^β$ no tiene ningún efecto en el p-valor para contrastar $H_0 : β = 0$.

modelo.cox1 <-coxph(Surv(ttr,relapse)~ grp)
summary(modelo.cox1)

## Call:
## coxph(formula = Surv(ttr, relapse) ~ grp)
## 
##   n= 125, number of events= 89 
## 
##                coef exp(coef) se(coef)   z Pr(>|z|)   
## grppatchOnly 0.6050    1.8313   0.2161 2.8  0.00511 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## grppatchOnly     1.831     0.5461     1.199     2.797
## 
## Concordance= 0.581  (se = 0.027 )
## Likelihood ratio test= 7.99  on 1 df,   p=0.005
## Wald test            = 7.84  on 1 df,   p=0.005
## Score (logrank) test = 8.07  on 1 df,   p=0.004

Ajustábamos el modelo para la variable tipo de tratamiento. Es muy importante en el caso de los modelos de supervivencia en general, mi varaible respuesta no es ni el tiempo solo ni el evento, es una combinación de ambas. porque es el tiempo hasta que yo observo el enveto o el tiempo hasta que tengo la censura.con lo cual necesito la informacion combinada de ambas variables. para decirle eso necesito un objeto tipo survival como el de la orden, y luego gusanito y mi variable.

que vemos en la salida:

En este caso tenemos solo dos categorias, el parche o el medicamento. En este caso el medicamento sería la categoria de referencia, y obtenemos el beta asociado al parche (0.60), su riesgo relativo (1.8, recordad que el riesgo relativo lo calculabamos con la exponencial del beta) y el p-valor esta asociado al test de wall. luego abajo lo tengo igual.

luego tenemmos el test de lar razon de VS. volvemos a tener un pvalor muy similar al test de wall lo cual quiere decir que mi variable es significativa, es decir, que el riesgo de volver a fumar o de tener una recaida cambia en función delt tipo de medicamento que este utilizando.Pues esto seria elmodelo univariante

####INTEPRETACIÖN:

Los hazard ratio se pueden interpretar en términos de cambios porcentuales en el riesgo.
En el ejemplo anterior, el estimador $β = 0.6050$ para la tratamiento “parche”, y por lo tanto, su hzard ratio correspondiente será $e^β = e^{0.6050} = 1.8313$, lo que significa que el riesgo de recaer y volver a fumar incrementa un 83% en los individuos con parche con respecto a los individos con el tratamiento triple.
La hipótesis de riesgos proporcionales asume que el riesgo no varía con el tiempo, y por tanto, la relación de riesgo entre los dos grupos en este caso, consideramos que se mantiene constante el tiempo, independientemente de cuando haga el estudio.
Por otro lado, se observa la estimación del parámetro $β$ obtenido es significativa (test de Wald y razón de verosimilitud).

###MODELO MULTIVARIANTE

Una vez que hemos establecido un conjunto de covariables q, algunas de las cuales son variables dummies y algunas variables continuas, podemos escribir el modelo de la siguiente manera: \[h_1(t) = h_0(t)ψ = h_0(t)e^{\sum_{j=1}^{q}β_jX_j}\] \[log(ψ) = β_1X_1 + . . . + β_qX_q\] Para cada covariable, el parámetro $β_j$ es el log hazard ratio para el efecto de la variable sobre la supervivencia, ajustado por las otras covariables. Para covariables continuas, representa el efecto de un cambio de unidad en la covariable; para variables dummies por contra, representa el efecto del nivel correspondiente en comparación con la categoría de referencia.

Test de la razón de verosimilitud (multivariable)

Si queremos comparar dos modelos anidados podemos utilizar el test de la razón de verosimilitud.

El estadístico se basa en las diferencias de la log verosimilitud parcial de los dos modelos: \[2\left(l(\hat{β}reduced ) − l(\hat{β}full)\right)\] Este estadístico sigue una distribución Chi-cuadrado, siendo los grados de libertad la diferencia en parámetros entre $\hat{β}reduced$ y $\hat{β}full$.

Para obtener la log verosimilitud parcial de un modelo de Cox en R podemos utilizar la función logLik.

entonces voy a rechazar la hipotesis nula, es decir, me va a merecer la pena meter dos variables mas en el modelo si la diferencia entre mis vSpar es lo suficientmente grande y obtengo un p-valor pequeño.

Comparación de modelos no- anidados

A menudo tenemos un gran número de covariables potenciales y necesitamos podar el modelo para que sólo se incluyan las covariables necesarias.
Hay una variedad de herramientas disponibles para ayudar en este proceso.
Un método muy utilizado consiste en realizar la selección del modelo mediante un proceso “stepwise”.
En la versión “forward” de este método, primero ajustamos los modelos univariantes para cada covariable. Se selecciona la covariable con el p-valor más pequeño y se añade al modelo base.
Luego, con esa covariable incluida, se ajusta un modelo adicional con cada una de las covariables restantes por separado, y se selecciona a la mejor (segunda covariable, la que tiene el p-valor más pequeño), de modo que tenemos un modelo con dos covariables.
Continuamos hasta que ninguna covariable adicional tenga un p-valor menor que un cierto valor crítico α;

+Otra versión, conocida como el procedimiento “backward”, comenzamos con todas las covariables en el modelo, las eliminamos una a una, en función del mayor p-valor (no significativo), hasta que todas las variables incluidas en el modelo son significativas.

####AIC para comparar modelos no-anidados

Sin embargo, los procedimientos stepwise, aunque es cierto que nos pueden dar una idea de las variables a incorporar en un modelo, tienen algunas limitaciones:

múltiples comparaciones
el p-valor obtenido corresponde a la comparación de modelos anidados

Una alternativa es comparar modelos no anidados mediante el criterio de información de Akaike (AIC). \[AIC = −2 · l(\hat{β}) + 2q\], siendo q el número de parámetros del modelo.

Al comparar modelos de AICs, el tamaño muestral tiene que ser el mismo.

# la variable edad que la teníamos antes de manera continua la ehmos categorizado en 4 niveles
tabaquismo$ageGroup4=factor(tabaquismo$ageGroup4)
levels(tabaquismo$ageGroup4)

## [1] "21-34" "35-49" "50-64" "65+"

# y luego tenemos la varibale tipo de trabajo,jornada parcial, jornada completa u otro 
tabaquismo$employment=factor(tabaquismo$employment)
levels(tabaquismo$employment)

## [1] "ft"    "other" "pt"

## entonces nos planeamos la necesidad de ajustar un modelo con dos variables.
modelo.cox2 <-coxph(Surv(ttr,relapse)~ ageGroup4 + employment)
modelo.cox3 <-coxph(Surv(ttr,relapse)~ ageGroup4)
modelo.cox4 <-coxph(Surv(ttr,relapse)~ employment)

¿Cómo saber si en este modelo multivariante, cualqueira de esas dos variables es significativa?- comparar con los univariantes.

si yo comparo cox2 con cox3 lo que estoy comparando es el modelo pequeño que solo tiene edad, con el modelo grande que tiene edad y tipo de rtabajo. Si obtengo un pvalor significativo quiere decir que la variable employment es significativa, que necesito meterla en el modelo. Porque aqui obtengo une stadistico G, una diferencia entre las logverosimilitudes parciales lo suficientemente grande para que me merezca meter esas variables.

summary(modelo.cox2)$coefficients

##                       coef exp(coef)  se(coef)          z    Pr(>|z|)
## ageGroup435-49  -0.1299275 0.8781591 0.3213101 -0.4043679 0.685942230
## ageGroup450-64  -1.0238705 0.3592019 0.3585151 -2.8558649 0.004291978
## ageGroup465+    -0.7824610 0.4572793 0.5046472 -1.5505109 0.121018946
## employmentother  0.5257333 1.6916990 0.2748459  1.9128297 0.055769852
## employmentpt     0.5000966 1.6488805 0.3315223  1.5084854 0.131430353

veos el modelo de las dos covariables y vemos lo pvalores asociados al test de wall. cada uno de los pvalores lo que me esta diciendo es si esa categoria es estadísticamente diferente con respecto a la categoría de referencia.

en este caso tenemos diferencias estadisticamente significativas en el grupo de 50-64 años con repspecto al de menos de 35 años. Ese, es diferente al de menos de 35 años, pero no lo es el de 35-49, ni el de 65+.

Con esta información y estos 3 pvalores yo no se si la variable edad es significativa o no. yo lo único que se es que hay dos categorías entre las que hay diferencias, pero no si si son lo suficiente para que la variable en su conjunto sea significativa.

y con los empleos para lo mismo. employmentother esta ahi muy al limite con un alpha justo de 0.06 casi, y el otro no es significativo. ademas sus betas son muy parecidos, lo que me da para pensar que entre otro tipo de trabajo y el trabajo aparcial no hay evidencias ES SIG. porque los betas asociados son muy parecidos. con lo cual con esto me da por pensar que el tipo de trabajo no va a ser significativo, pero no lo se.

summary(modelo.cox2)$logtest

##         test           df       pvalue 
## 16.787204170  5.000000000  0.004921562

Este p-valor de aqui lo que me dice es que este modelo, cuando lo comparo con el modelo nulo, es decir, con no hacer nada, si es significativo. ahora, necesito las dos variables, o me vale con una? Para eso comparo las verosimilitudes de los 2 modelos. Puedo utilizar las funciones loglik. Calculo la logVS de uno y de otro. importante, GL, la diferencia va en funcion de cual compare (5vs3). cuando comparo esos dos, la diferecia en grados de libertad es dos.

logLik(modelo.cox2)

## 'log Lik.' -377.7597 (df=5)

logLik(modelo.cox3)

## 'log Lik.' -380.043 (df=3)

logLik(modelo.cox4)

## 'log Lik.' -385.1232 (df=2)

\[2 \left(l(\hat{β}reduced) − l(\hat{β}full)\right) = 2(-377.7597 + 380.043) = 4.567\]

pchisq(4.567, df=2,lower.tail=FALSE)

## [1] 0.1019268

pero no es necesario estar haciendo todos los loglik de todos los modelos. Podemos utilizar la funcion anova

anova(modelo.cox2,modelo.cox3)

## Analysis of Deviance Table
##  Cox model: response is  Surv(ttr, relapse)
##  Model 1: ~ ageGroup4 + employment
##  Model 2: ~ ageGroup4
##    loglik  Chisq Df P(>|Chi|)
## 1 -377.76                    
## 2 -380.04 4.5666  2    0.1019

la funcion anova al modelo.cox2 sabe que es una objeto tipo coxph entonces ya saber que lo que tiene que compara son las logverosimilitudes parciales. Y entonces compara el 2 con el 3. lo que estoy viendo aqui es si el tipo de trabajo es significativo o no. obtenemos un pvlor de 0.1, con lo que comluimos que en ese modelo multiarante la variable tipod e trabajo no es significativa.

Sin embargo si comparamos el modelo multivariante y miro la significacion de la variable edad categorizada, (fijaros que aqui tengo los 3 grados de libertad) y en este caso la diferencia si es estadisticamente significativa.

anova(modelo.cox2, modelo.cox4)

## Analysis of Deviance Table
##  Cox model: response is  Surv(ttr, relapse)
##  Model 1: ~ ageGroup4 + employment
##  Model 2: ~ employment
##    loglik  Chisq Df P(>|Chi|)   
## 1 -377.76                       
## 2 -385.12 14.727  3  0.002065 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Es decir, que en este modelo multivariantela variable edad es significativa pero la variable tipo de trabajo no lo es. Con lo cual, teniendo esas dos variables, nos tendriamos que quedar solo con el modelo que tuviese la varaible edad categorizada.

AIC(modelo.cox2)

## [1] 765.5194

AIC(modelo.cox3)

## [1] 766.086

AIC(modelo.cox4)

## [1] 774.2464

Fijaros que si comparamos los modelos, curiosamente, en base al aic no obtengo el mismo resultado que en base al p-valor. porque en base al AIC, el modelo 2, que es el modelo multivariante, es el que me da menor AIC, a pesar de que la variable no es significativa.

MIconsejo, o en la forma en la que yo entiendo que hay q ajustar los modelos y ver que son validos. A mi no me interesa ajustar el aic de un modelo done todas mis variables no son significativas. El primer criterio es que la variable sea significativa. luego ya mirare los demás paramtros que le puedan afectar como le afectan.

A nivel estadistico no tiene sentido quedarnos con una variable que no sea significativa. Y una vez que tengamos modelos distintos donde todo sea signicativo, usaré el aic para comparalos.

modelo.cox.all<-coxph(Surv(ttr,relapse) ~ grp + gender + race +  employment + yearsSmoking + levelSmoking + ageGroup4 + priorAttempts + longestNoSmoke )

result.step<- step(modelo.cox.all, scope = list(upper = ~ grp + gender + race +  employment + yearsSmoking + levelSmoking + ageGroup4 + priorAttempts + longestNoSmoke, lower = ~ grp ))

## Start:  AIC=770.2
## Surv(ttr, relapse) ~ grp + gender + race + employment + yearsSmoking + 
##     levelSmoking + ageGroup4 + priorAttempts + longestNoSmoke
## 
##                  Df    AIC
## - race            3 766.98
## - yearsSmoking    1 768.20
## - gender          1 768.20
## - priorAttempts   1 768.24
## - levelSmoking    1 768.47
## - longestNoSmoke  1 769.04
## <none>              770.20
## - employment      2 772.45
## - ageGroup4       3 774.11
## 
## Step:  AIC=766.98
## Surv(ttr, relapse) ~ grp + gender + employment + yearsSmoking + 
##     levelSmoking + ageGroup4 + priorAttempts + longestNoSmoke
## 
##                  Df    AIC
## - levelSmoking    1 764.98
## - gender          1 765.00
## - priorAttempts   1 765.01
## - yearsSmoking    1 765.04
## - longestNoSmoke  1 766.29
## <none>              766.98
## - employment      2 768.37
## - ageGroup4       3 770.16
## + race            3 770.20
## 
## Step:  AIC=764.98
## Surv(ttr, relapse) ~ grp + gender + employment + yearsSmoking + 
##     ageGroup4 + priorAttempts + longestNoSmoke
## 
##                  Df    AIC
## - gender          1 763.00
## - priorAttempts   1 763.01
## - yearsSmoking    1 763.06
## - longestNoSmoke  1 764.29
## <none>              764.98
## - employment      2 766.37
## + levelSmoking    1 766.98
## - ageGroup4       3 768.18
## + race            3 768.47
## 
## Step:  AIC=763
## Surv(ttr, relapse) ~ grp + employment + yearsSmoking + ageGroup4 + 
##     priorAttempts + longestNoSmoke
## 
##                  Df    AIC
## - priorAttempts   1 761.02
## - yearsSmoking    1 761.08
## - longestNoSmoke  1 762.31
## <none>              763.00
## - employment      2 764.42
## + gender          1 764.98
## + levelSmoking    1 765.00
## - ageGroup4       3 766.32
## + race            3 766.48
## 
## Step:  AIC=761.02
## Surv(ttr, relapse) ~ grp + employment + yearsSmoking + ageGroup4 + 
##     longestNoSmoke
## 
##                  Df    AIC
## - yearsSmoking    1 759.10
## - longestNoSmoke  1 760.34
## <none>              761.02
## - employment      2 762.42
## + priorAttempts   1 763.00
## + gender          1 763.01
## + levelSmoking    1 763.02
## - ageGroup4       3 764.50
## + race            3 764.52
## 
## Step:  AIC=759.1
## Surv(ttr, relapse) ~ grp + employment + ageGroup4 + longestNoSmoke
## 
##                  Df    AIC
## - longestNoSmoke  1 758.42
## <none>              759.10
## - employment      2 760.42
## + yearsSmoking    1 761.02
## + gender          1 761.08
## + levelSmoking    1 761.08
## + priorAttempts   1 761.08
## + race            3 762.52
## - ageGroup4       3 766.90
## 
## Step:  AIC=758.42
## Surv(ttr, relapse) ~ grp + employment + ageGroup4
## 
##                  Df    AIC
## <none>              758.42
## + longestNoSmoke  1 759.10
## - employment      2 760.31
## + yearsSmoking    1 760.34
## + gender          1 760.39
## + priorAttempts   1 760.40
## + levelSmoking    1 760.41
## + race            3 761.53
## - ageGroup4       3 767.24

Y el tema de la funcion etep de la que hablamos. Ajustamos un modelo con todas las variables. y caundo hago un fucnion step lo que le digo es cual es mi modelo base, de dode parto, mi lista desde donde empiezo(cual sería mi modelo posible mas grande, que seria mi modelo upper), y el mas pequeño ( aqui le vamos a decir: y ya se que el tipo de parche que vamos a utilizar es significativo, y ademas es lo que me interesa analizar en este estudio, asi que esa variable dejamela. )

Partiendo de la base que en el modelo mas pequeño esa va a estar, de todo eso, ¿con qué me quedo?

Entones vamos a ir viendo todo el proceso y vamos a ir viendo en cada una de las iteraciones que decisiones toma y por que.

Partimos del modelo con otdas las variables y me va diciendo cual es el AIC.recordad que antes hemos dicho que la funcion step la podemos calcular en funcion del p-valor, pero que eso tenia problemas asociados, con lo cual lo tenemos encuenta en función del aic. Esto es una herramienta para que os hagais una idea cuando tengais muchas muchas variables.

Se reduce el problema a las mas importantes y despues incluis de una en una las que realmente importan.

Entonces esto que quiere decir. Yo parto de un Aic, y me va sacando los aic que deja el modelo tras quitar la variable . hasta conseguir el ultimo modelo en el que ya me empieza a subir el Aic.

Selección de variables del modelo

Una vez que nos hemos quedado con un modelo final, tenemos que comprobar que se cumplen las hipotesisi que hemos asumido incialmente.

Del mismo modo que ocurre en otros modelos de regresión, debemos de comprobar que se cumplen las siguientes condiciones con las variables continuas que introducimos en el modelo:

Relación lineal entre las covariables continuas y el log hazard.
Incorporar en caso necesario términos de interacción.

¿Cómo podemos verificar la relacion lineal?

Puedo ajustar un modelo con una spline, que no tiene por que tener una forma polinomica por defecto. en este caso, cuando le digo df=4, lo que le estoy diciendo es que me ajuste una curva donde sus grados de libertad vayan hasta el orden 4. una vez que yo le meto esto lo que le estoy diciendo es que esta variable edad de continuo me la ajuste con una funcion suave.

modelo.cox.spline <- coxph(Surv(ttr,relapse) ~ grp + employment + pspline(age, df=4))
modelo.cox.spline

## Call:
## coxph(formula = Surv(ttr, relapse) ~ grp + employment + pspline(age, 
##     df = 4))
## 
##                              coef se(coef)     se2   Chisq   DF       p
## grppatchOnly               0.6507   0.2210  0.2194  8.6738 1.00 0.00323
## employmentother            0.6330   0.2774  0.2750  5.2068 1.00 0.02250
## employmentpt               0.5700   0.3403  0.3332  2.8051 1.00 0.09396
## pspline(age, df = 4), lin -0.0339   0.0102  0.0102 11.0668 1.00 0.00088
## pspline(age, df = 4), non                           4.2026 3.08 0.25164
## 
## Iterations: 3 outer, 9 Newton-Raphson
##      Theta= 0.709 
## Degrees of freedom for terms= 1.0 2.0 4.1 
## Likelihood ratio test=27.3  on 7.02 df, p=3e-04
## n= 125, number of events= 89

En la salidda obtengo 2 lineas. una la linela, y otra la correspondiente a la no lineal. Entonces lo que me viene a decir aqui es que la lineal es significativa, y la no lineal no es significativa.

lo cual quiere decirq ue asumiendo una relacion lineal entre mi edad y el logaritmo de pi, lo estoy haciendo bien, no necesito meter nada mas.

Pero imaginaros que me hubiese dado que la funcion suave es significativa, que hago?

Introducir la variable al cuadrado en el modelo.(la facil) ¿qué se hace muchas veces? categorizar la variable. y una vez la he categorizado me puedo olvidar de si la variable es lineal o no lo es.

Importante en el modelo de cox de riesgo proporcionales cuadno yo estoy asumiendo que los riesgos son proporcionales, si me equivoco en la relacion funcional de las variables, es decir, si yo asumo que es lineal cuando no lo es, puedo influir en los riesgos proporcionales. es decir, puedo estar obteniendo que mis riesgos son proporcionales cuando mi probema de fondo es que la relacion que he asumido no es la correcta. entonces la opcion de categorizar es muy sencilla por que me permite comparar grupos de riesgo, y ese riesgo va a depender de los grupos.

en cualquier caso lo importante es que os asegureis que se esta dando esa linealidad.

termplot(modelo.cox.spline, se=T, terms = 3, ylabs = "Log hazard")

En este caso podemos dibujar esa relacion. esa seria esta funcion suave que estamos estimando, le digo que me dibuje el intervalo de confianza, y cuando le digo terms = 3 es por que quiero ver el efecto de mi tercera variable, es decir, la edad.

Y lo que vemos es que mas o menos (no es una linea recta rectisima pero siempre baja), salvo aqui al final. pero que pasa aqui al final, que como veis, tengo un intervalo de confianza super ancho. Eso quiere decir que tengo muy pocos individuos y si que sube un poco pero realmente podria subir hasta arriba o bajar hasta abajo. no tengo la suficiente exactitud en ese cambi de ahi para decir o como para rechazar que sea lineal. por lo tanto asumo que es lineal, por que no tengo datos suficientes como para decir que no lo es.

Asunción de riesgos proporcionales

La suposición de riesgos proporcionales es clave para la construcción de la verosimilitud parcial, ya que es esta propiedad la que permite cancelar el riesgo basal de los factores de verosimilitud parciales.
Si uno tiene una variable predictiva binaria, como el tratamiento experimental vs. el tratamiento estándar, lo que esta suposición significa es que las funciones de riesgo son proporcionales y, por lo tanto, los log hazards están separados por una constante en todos los puntos de tiempo.
Del mismo modo, una variable categórica con muchos niveles resultan en funciones logaritmo de las hazard paralelas.

Esta suposición es, en el mejor de los casos, una aproximación en la práctica, y es poco probable que un pequeño desvio en esta suposición tenga efectos importantes en la inferencia de los parámetros del modelo.

Log Cumulative Hazard Plots

¿Cómo vamos a ver si se cumple la proporcionalidad de los riesgos o no?

Si comparamos los tiempos de supervivencia entre dos grupos, hay un gráfico simple que puede ayudarnos a evaluar el supuesto de riesgos proporcionales.

Bajo la asunción de riesgos proporcionales, tenemos \[S_1(t) = [S_0(t)]^{exp(β)}\] donde exp(β) es la constante de los riesgos proporcionales (proportional hazards). Tomando logaritmos, \[log [S_1(t)] = exp(β) · log [S_0(t)]\]

Por esta razón, los test estadísticos de proporcionalidad tienen a menudo un valor limitado. Aún así, es útil evaluar, en los datos disponible, si esta suposición es razonable.

Dado que las funciones de supervivencia son inferiores a 1, sus logaritmos son negativos. Por lo tanto, debemos negarlas antes de tomar un segundo logaritmo, \[log {−log [S_1(t)]} = β + log {−log [S_0(t)]}\] La función $g(u) = log {−log(u)}$ es conocida como la transformación complementaria log-log, y tiene el efecto de transformar valores de u en el intervalo $(0, 1)$ al intervalo $(−∞, ∞)$ para $g(u)$.

Un gráfico de $g [S_1(t)]$ y $g [S_0(t)]$ frente a $t$ o $log(t)$ producirá dos curvas paralelas separados por β si la suposición de riesgos proporcionales es correcta.

result.surv.patch<-survfit(Surv(ttr) ~ grp, subset = {grp == "patchOnly"})
time.patch <-result.surv.patch$time
surv.patch <-result.surv.patch$surv
cloglog.patch <-log(-log(surv.patch))
logtime.patch <-log(time.patch)

result.surv.comb <-survfit(Surv(ttr) ~ grp, subset = {grp == "combination"})
time.comb <-result.surv.comb$time
surv.comb <-result.surv.comb$surv
cloglog.comb<-log(-log(surv.comb))
logtime.comb<-log(time.comb)

plot(cloglog.patch  ~ logtime.patch, type ="s", col = "blue", lwd = 2, ylab = "supervivencia log-log",xlab="Log time")
lines(cloglog.comb ~ logtime.comb,  col = "red", lwd = 2,  type ="s")
legend("bottomright", legend =c("patchonly","combination"), col=c("blue","red"),lwd=2)

lo que estoy asjutando es esa log log cuando tengo el grupo del parche y cuando tengo el grupo del medicamento combinado. Y estas curvas son mas o menos paralelas. es decir, el paso de uno a otro a lo largo del tiempo se me mantiene constante. Si se separasen querria decir que varia el riesgo.

Y otra forma de ver esa proporcionalidad es mediante los residuos.

Residuos de Schoenfeld

Los gráficos de residuos de Schoenfeld proporcionan una forma útil de evaluar la asunción de riesgos proporcionales.

Los residuos se definen como \[\hat{r}_i = x_i − \sum_{l∈R_i}x_l· p(\hat{β}, x_l) = xi − \bar{x}(t_i).\] donde $x¯(t_i)$ es el valor esperado de la variable $X$ para el individuo i en el tiempo $t_i$.

Un gráfico de estos residuos frente a la covariable $X$ producirá un patrón de puntos que se centran en cero, si la suposición de riesgos proporcionales es correcta.

####Residuos transformados de Gramsch and Therneau

Gramsch y Therneau propusieron reescalar cada residuo con una estimación de su varianza. Este residuo reescalado puede ser convenientemente aproximado de la siguiente manera: \[r^∗_i = r_i· d · var(\hat{β})\] donde d es el número total de eventos, y $var(\hat{β})$ es la varianza del parámetro estimado. En R podemos evaluar estos residuos con la función cox.zph

modelo.cox.final <- coxph(Surv(ttr,relapse) ~ grp + employment
+ageGroup4)
result.res <- cox.zph(modelo.cox.final, transform="rank")
result.res

##            chisq df    p
## grp        0.114  1 0.74
## employment 0.624  2 0.73
## ageGroup4  3.203  3 0.36
## GLOBAL     4.008  6 0.68

lo que me interesa aqui es ver si se cumplen los riesgos proporcionales para diferntes variables que yo meto en el modelo. Entonces, hemos dicho antes que nos quedabamos con un modelo que teniamos la varaible que me repsentaba el tipo de medicamento, el trabajo y al variable edad categorizada en 4 grupos. entonces si a este modelo final el aplica esta funcion coxzph lo que me devuelve es esta salida de aqui. y esta columna final me represetna el p-valor. Enotnces tengo un pvalor asociado a cada una de als covariables, y un pvalor asociado de forma global al modelo.

Entonces, lo que yo stoy asumiendo, la hipotesis nula en este caso es que la diferencia

En este caso, como todos los pvalores que obtengo son mayores que alfa, no rechazo la hipotesis nula. es decir, puedo asumir que se cumplen la proporcionalidad y los riesgos para todas las variables.

R- cuadrado

El estadístico “r-squared” es una adaptación del R2 de la regresión lineal al análisis de supervivencia.Se define como: \[R^2 = 1 −\left({\frac{l(0)}{l(\hat{β})}}\right)^{2/n}\]$ y refleja la mejora en el ajuste del modelo con la covariable en comparación con el modelo nulo.