Monthly Archives: February 2014

Estimación no paramétrica de la función de supervivencia

Teniendo un estudio de estudio de supervivencia entre manos, puede ser complicado ajustar un modelo paramétrico a nuestros datos. Existen alternativas no paramétricas, como el estimador de Kaplan-Meier de la función de supervivencia \(\hat{S(t)}\), que se define de la siguiente manera:

\(\hat{S(t)}=\prod_{t_{i}\leq t}(1-\frac{d_{i}}{Y_{i}})\)

Suficiente teoría por hoy. Pasemos a un ejemplo práctico con R. Supongamos que tenemos una serie de individuos con tiempos de supervivencia y la típica variable de estado indicando si es un dato censurado o no.

time <- c(1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)
status <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) 

Cargamos el paquete survival.

library(survival)

A continuación usamos el comando survfit, encargado de realizar la estimación de Kaplan-Meier.

fit2 <- survfit(Surv(time, status) ~ 1) 
plot(fit2,xlab="Tiempo",ylab="Survival p.",conf.int=FALSE,col='blue')
title("Kaplan-Meier estimation")

Que produce el siguiente gráfico:

km

Ahora planteamos lo siguiente. Dadas dos estimaciones de probabilidades de supervivencia para dos grupos, ¿podemos afirmar que las curvas reales de supervivencia son iguales?. Es decir, nos gustaría plantear la hipótesis:

\(
H_{0}=\textrm{S(t) es igual para los dos grupos}
\)

Primero vamos a limitarnos a introducir los datos para nuestros grupos de la siguiente manera:

time1 <- c(6,6,6,7,10,13,16,22,23,6,9,10,11,17,19,20,25,32,32,34,35)
status1 <- c(1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0)

time2 <- c(1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)
status2 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) 

Ajustamos un modelo equivalentemente a lo anterior para cada uno de los grupos:

fit1 <- survfit(Surv(time1, status1) ~ 1)
fit2 <- survfit(Surv(time2, status2) ~ 1)

plot(fit1,xlab="Tiempo",ylab="Survival p.",conf.int=FALSE,col='blue')
lines(fit2,col='red',conf.int=FALSE)
legend(21,1,c('Group 1', 'Group 2'), col = c('blue','red'), lty = 1) 

El resultado es el que cabe esperar:
km2

Pero claro, aunque visualmente queda claro, no es suficiente evidencia como para afirmar que son distintas. Recurrimos a lo que se denomina Log-Rank test, que no es más que un test chi-cuadrado para contrastar nuestra hipótesis. Una expresión más detallada de este estadístico puede encontrarse al final del documento. En R este test puede hacerse de la siguiente manera (consideramos una variable auxiliar que distinga entre ambos grupos):

time <- c(6,6,6,7,10,13,16,22,23,6,9,10,11,17,19,20,25,32,32,34,35,1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,
12,12,15,17,22,23)
status <- c(1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) 
treatment <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2)

fit <- survdiff(Surv(time, status) ~ treatment) 

La salida de ‘fit’ nos da lo que nos interesa saber:

Call:
survdiff(formula = Surv(time, status) ~ treatment)

             N Observed Expected (O-E)^2/E (O-E)^2/V
treatment=1 21        9     19.3      5.46      16.8
treatment=2 21       21     10.7      9.77      16.8

 Chisq= 16.8  on 1 degrees of freedom, p= 4.17e-05 

Para un \(\alpha=0.05\) rechazaríamos nuestra hipótesis nula. Este test puede generalizarse para más grupos, incluso considerando estratos.

Tiempo residual medio

Supongamos que tenemos un estudio de análisis de supervivencia entre manos, y dada una variable \(X\) que define el tiempo de vida de algo, estamos interesados en saber en cuál es el tiempo de vida esperado de un individuo dado que ya ha vivido una cantidad de tiempo \(t\). Para responder a eso recurrimos a lo que se denomina tiempo residual medio, que se define como:

\(e_{x}(t)=E[X-t|X>t]=\int_{t}^{\infty}\frac{1-F_{x}(x)}{1-F_{x}(t)}dx\)

Donde \(F\) representa la función de distribución de la v.a. en cuestión. ¿Y cómo es esta función para algunos modelos conocidos? (Como somos un poco vagos para integrar a mano, usamos Wolfram Mathematica)

    • \(X\sim Exp(\lambda)\)


Integrate[(1 - CDF[ExponentialDistribution[\[Lambda]], x])/(1 -
CDF[ExponentialDistribution[\[Lambda]], t]), {x, t, Infinity},
Assumptions -> \[Lambda] > 0 && t > 0 && x > 0]

Resultado: \(\frac{1}{\lambda}\). Es decir, para el modelo exponencial el tiempo residual medio es constante. Esto es consecuencia directa de la propiedad de pérdida de memoria.

 

    • \(X\sim Weib(\alpha,\beta)\)


Integrate[(1 - CWeibull[x, \[Alpha], \[Beta]])/(1 -
CWeibull[t, \[Alpha], \[Beta]]), {x, t, Infinity},
Assumptions -> t > 0 && x > 0 && \[Alpha] > 0 && \[Beta] > 0 ]

Resultado: \(\frac{e^{(t/\alpha)^{\beta}} \alpha \gamma(1/\beta,(t/\alpha)^{\beta})}{\beta}\), donde \(\gamma\) es la función gamma incompleta.

Para el modelo Weibull dependiendo de los parámetros \((\alpha,\beta)\), la función en cuestión sería creciente (cuánto más tiempo vive un elemento, más tiempo esperado de supervivencia tiene) o decreciente (recíproco). Para comprobar este hecho, podríamos estudiar el signo de la primera derivada de la función igualándola a 0 y despejando t, pero hacer esto no es fácil. Para comprobar el comportamiento de dicha función empíricamente podemos recurrir de nuevo a Mathematica:


Manipulate[
Plot[(E^(t/\[Alpha])^\[Beta] \[Alpha] Gamma[
1/\[Beta], (t/\[Alpha])^\[Beta]])/\[Beta], {t, 0,
8}], {\[Alpha], -8, 8}, {\[Beta], -2, 2}]

manipulate