Al final del capítulo 7 hay un ejercicio que pide construir la función de supervivencia de 10 pacientes que han sido aleatoriamente asignados a los tratamientos A y B. Con esosmismos datos en el capítulo 9.3 se pide hacer la prueba de comparar dichas funciones.

CAPITULO 7

Suponga que disponemos de los datos de supervivencia de 10 pacientes que han sido aleatoriamente asignados a los tratamientos A y B.

A: 3, 5, 7, 9+, 18

B: 12, 19, 20, 20+, 33+

Construya la función de supervivencia para cada tratamiento y grafíquelas. ¿Qué se puede decir de los tratamientos a partir de las gráficas?

librerias

library(ggplot2)
library(latex2exp)
library(knitr)
library(kableExtra)

Datos

# tratamiento A
tratamientoA = data.frame(Tiempo = c(3,5,7,9,18),
                          Censura= c(0,0,0,1,0),
                          Status = c(1,1,1,0,1))
# tratamiento B
tratamientoB = data.frame(Tiempo = c(12,19,20,20,33),
                          Censura= c(0,0,0,1,1),
                          Status = c(1,1,1,0,0))

Para calcular Kaplan-Meier usamos la funcion

El método K-M se basa en el uso de la Log-Verosimilitud para estimar los parametros que permitan calcular la función de supervivencia.

Recordemos que para una muestra aleatoria \(\{T_i\}_{i=1}^n\) con soporte discreto en \(\{u_1,u_2,...\}\) la función de supervivencia usando K-M se estima como sigue:

\[\hat S(t) = \prod_{k:u_k\leq t}\left(1-\frac{d_k}{n_k}\right)\] Donde:

\[d_k= \text{número de tiempos de fallo iguales a } u_k \\ n_k= \text{número de individuos en riesgo al tiempo } u_k\] Usamos la siguiente funcion

#Funcion que calcula la supervivencia mediante el metodo Kaplan-Meier

K_M = function(d){#d  datos(dataframe)
  #soporte de los tiempos(uk)
  t<-sort(unique(d$Tiempo))
  l<-length(t)
  #calculamos
  dk = rep(0,l)#dk:numero de tiempos de fallo iguales a uk
  nk = rep(0,l)#nk:numero de individuos en riesgo a tiempo uk
  ck = rep(0,l)#ck:numero de censuras a tiempo uk
  Sk = rep(1,l)#Sk:supervivencia
  pk = rep(0,l)
  for(i in 1:(l)){
    dk[i] = nrow(subset(d, Tiempo==t[i] & Status==1))# si no hay censura
    ck[i] = nrow(subset(d, Tiempo==t[i] & Status==0))# si hay censura
    nk[i] = nrow(subset(d, Tiempo>=t[i]))
    pk[i] = 1-dk[i]/nk[i]
    Sk[i] = prod(pk[1:i])# el productos
  }
  resultados<-data.frame(t,dk,nk,ck,Sk)  
  return(resultados)
}

Aplicamos la funcion los tratammientos A y B

tratamientoA_KM = K_M(tratamientoA)
tratamientoB_KM = K_M(tratamientoB)

Imprimimos tablas

tratamientoA_KM

tratamientoB_KM

graficamos

ggplot()+ 
  geom_step(aes(x=tratamientoA_KM$t,y=tratamientoA_KM$Sk, color ="Tratamiento A"),size=2)+
  geom_step(aes(x=tratamientoB_KM$t,y=tratamientoB_KM$Sk, color ="Tratamiento B"),size=2)+
  theme_classic()+
  labs(title="K-M", x = "Tiempo", y = "Supervivencia estimada") +
  scale_color_manual(breaks=c("Tratamiento A","Tratamiento B") ,values = c("Tratamiento A"="green","Tratamiento B"="blue")) +
  theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
  labs(color = "")

CAPITULO 9

Pruebas de bondad de ajuste con Wilcoxon

Utilizaremos la prueba Wilcoxon para comparar dos funciones de supervivencia

\[H_0: S_1(t)=S_2(t)~~~\forall t>0\\ vs \\ H_a: S_1(t)\neq S_2(t)~~p.a.~t>0\] El estadistico de esta prueba es:

\[W^2=\frac{ (\sum_{i=1}^kn_{i}(d_{1i}-e_{1i}))^2}{\sum_{i=1}^kn_{i}^2V_{1i}}\] con:

\[e_{1i}=n_{1i}\frac{d_i}{n_i}\\ V_{1i}=\frac{n_{1i}n_{2i}d_i(n_i-d_i)}{n_i^2(n_i-1)}\] Se rechaza \(H_0\) cuando: \[W^2>J_{1-\alpha}\] Donde \(J_{1-\alpha}\) es el cuantil \(1-\alpha\) de una distribución \(\chi^2_{(1)}\).

Los resultados son los siguientes:

Creamos una funcion que haga la prueba Wilcoxon

#recibe las 2 funciones de supervivencia que serán testeadas

Wilcoxon = function(S1,S2){# funciones
  n1i = S1$nk
  n2i = S2$nk
  d1i = S1$dk
  di = S1$dk + S2$dk 
  ni = S1$nk + S2$nk
  
  e1i = n1i*di/ni
  v1i = n1i*n2i*di*(ni-di)/((ni^2)*(ni-1))
  
  # estadístico
  W <- (sum(ni*(d1i-e1i)))^2/(sum(ni^2*v1i))
  
  # Regla de decisión
  if (W > qchisq(0.95, df = 1)){conclusion <- "Se rechaza Ho"
  }else{  conclusion <- "No se rechaza Ho"}
  
  resultado = c("Estadístico: ",(sum(ni*(d1i-e1i)))^2/(sum(ni^2*v1i)),"\ncuantil al 95%:", qchisq(0.95, df = 1),conclusion)
  return(resultado)
}

aplicamos la funcion

result = Wilcoxon(tratamientoA_KM,tratamientoB_KM)

## Warning in S1$dk + S2$dk: longitud de objeto mayor no es múltiplo de la longitud
## de uno menor

## Warning in S1$nk + S2$nk: longitud de objeto mayor no es múltiplo de la longitud
## de uno menor

## Warning in n1i * n2i: longitud de objeto mayor no es múltiplo de la longitud de
## uno menor

result

## [1] "Estadístico: "     "0.169719827586207" "\ncuantil al 95%:"
## [4] "3.84145882069412"  "No se rechaza Ho"

EJERCICIOS CAP 7 Y CAP 9.3

Cruz Mateo David