Tabla de superviviencia
Edson Rolando Tmayo Castañeda
13 de mayo de 2017
Los siguientes datos(Descargar) fueron tomados de la tabla 3.1 de la tercera edición del libro Statistical Methods for Survival Data Analyse de Elisa T. Lee y John Wenyu Wang. Las columnas “Censura Rd” y “Censura ST” fueron agregadas para indicar sensura en los datos de Remission Duration y Survival Time respectivamente, donde “s” indica que el dato esta censurado y “m” lo contrario.
Son datos de 30 pacientes con melanoma en 4 estados, donde se realizo un estudio para comparar dos tratamientos BCG(1) y Corynebacterium parvun(2). SE midio el tiempo supervivencia y la duracion de la remisión.
#leemos los datos
library(dplyr)
library(xtable)
library(openxlsx)
sup<-read.xlsx("superviv.xlsx")
#Tomamos los datos necesarios para estimar la funcion de superviencia s(t)
sup<-select(sup,Patient,Tratament.Received,Survival.Time,Censura.ST)
#Ordenamos los datos por el tiempo de supervivencia
sup<-sup[order(sup$Survival.Time),]Ya que tenemos los datos ordenados por tiempo de supervivencia estimaresmo una funcion para la misma, por separado pra cada uno de los dos tratamientos.
Tratamiento 1 BCG
#Filtramos los datos para el tratamiento BCG
BCG<-filter(sup,Tratament.Received==1)
#Calculamos r "posicion de los datos sin censura"
r<-vector()
for(i in 1:nrow(BCG)){
if(BCG$Censura.ST[i]=="m"){
r[i]<-i
}else{
r[i]<-NA
}
}
BCG<-cbind(BCG,r)Calcularemos la proporcion de pacientes que sobrevie a cada intervalo como
\[sup\_int=\dfrac{n-r}{n-r+1}\]
donde \(n\) es el tamaño de la muestra y \(r\) es el rango no censurado.
#Calcular el tamaño de la muestra
n=nrow(BCG)
#Calcular la proporcion de supervivienca por intervalo.
BCG<-mutate(BCG,sup_int=(n-r)/(n-r+1))Ahora calcularemos \[\hat{s}(t)=\hat{s}(t-1)*sup\_int_t\]
denotaremos en la tabla a \(\hat{s}(t)\) simplemente como “S”.
BCG$S<-NA
for(i in 1:n){
if((!is.na(BCG$r[i]))&(BCG$r[i]==1)){
BCG$S[i]<-BCG$sup_int[i]
a<-i
}
if((!is.na(BCG$r[i]))&(!BCG$r[i]==1)){
BCG$S[i]<-BCG$sup_int[i]*BCG$S[a]
a<-i
}
}Calcularemos la desviación estandar de \(\hat{S}(t)\) con la formula de Greenwood.
\[s.e.\{\hat{S}(t)\}\approx [\hat{S}(t)]\left\{\sum_{j=1}^{k}\dfrac{d_j}{(n-r+1)(n-r)}\right\}^{1/2}\] Tomando en cuenta que para todo j \(d_j=1\)
g<-BCG[BCG$Censura.ST=="m",]
aux<-vector()
se<-vector()
for(i in 1:nrow(g)){
for(j in 1:i){
aux[j]<-1/((n-g$r[j])*(n-g$r[j]+1))
}
se[i]<-(g$S[i])*(sum(aux)^(.5))
}
g<-cbind(g,se)Calculemos los intervalos de confianza para S al 95%. Dados por
\[\hat{S}(t)\pm1.96s.e.\{\hat{S}(t)\}\]
#Calcula los limites del intervalo de confianza
g<-mutate(g, li=S-1.96*se,ls=S+1.96*se)
#Une los datos de g con los datos completos
BCG<-full_join(BCG,select(g,r,se,li,ls),by="r")
#imprime la tabla de resultados
print(xtable(BCG),type="html")| Patient | Tratament.Received | Survival.Time | Censura.ST | r | sup_int | S | se | li | ls | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2.00 | 1.00 | 3.90 | m | 1 | 0.91 | 0.91 | 0.09 | 0.74 | 1.08 |
| 2 | 4.00 | 1.00 | 5.40 | m | 2 | 0.90 | 0.82 | 0.12 | 0.59 | 1.05 |
| 3 | 7.00 | 1.00 | 7.90 | m | 3 | 0.89 | 0.73 | 0.13 | 0.46 | 0.99 |
| 4 | 3.00 | 1.00 | 10.50 | m | 4 | 0.88 | 0.64 | 0.15 | 0.35 | 0.92 |
| 5 | 9.00 | 1.00 | 16.60 | s | ||||||
| 6 | 8.00 | 1.00 | 16.90 | s | ||||||
| 7 | 11.00 | 1.00 | 17.10 | s | ||||||
| 8 | 5.00 | 1.00 | 19.50 | m | 8 | 0.75 | 0.48 | 0.18 | 0.13 | 0.82 |
| 9 | 6.00 | 1.00 | 23.80 | s | ||||||
| 10 | 1.00 | 1.00 | 33.70 | s | ||||||
| 11 | 10.00 | 1.00 | 33.70 | s |
Grafia de superviencia para el tratamiento BCG
plot(g$S~g$Survival.Time,data=g,type="s"
,xlim=c(0,max(g$Survival.Time)),ylim=c(0,max(g$ls))
,xlab="Survival Time",ylab="S(t)")
lines(g$li~g$Survival.Time,data=g,type="s")
lines(g$ls~g$Survival.Time,data=g,type="s")
Tratamiento 2 Corynebacterium
#Filtramos los datos para el tratamiento BCG
BCG<-filter(sup,Tratament.Received==2)
#Calculamos r "posicion de los datos sin censura"
r<-vector()
for(i in 1:nrow(BCG)){
if(BCG$Censura.ST[i]=="m"){
r[i]<-i
}else{
r[i]<-NA
}
}
BCG<-cbind(BCG,r)Calcularemos la proporcion de pacientes que sobrevie a cada intervalo como
\[sup\_int=\dfrac{n-r}{n-r+1}\]
donde \(n\) es el tamaño de la muestra y \(r\) es el rango no censurado.
#Calcular el tamaño de la muestra
n=nrow(BCG)
#Calcular la proporcion de supervivienca por intervalo.
BCG<-mutate(BCG,sup_int=(n-r)/(n-r+1))Ahora calcularemos \[\hat{s}(t)=\hat{s}(t-1)*sup\_int_t\]
denotaremos en la tabla a \(\hat{s}(t)\) simplemente como “S”.
BCG$S<-NA
for(i in 1:n){
if((!is.na(BCG$r[i]))&(BCG$r[i]==1)){
BCG$S[i]<-BCG$sup_int[i]
a<-i
}
if((!is.na(BCG$r[i]))&(!BCG$r[i]==1)){
BCG$S[i]<-BCG$sup_int[i]*BCG$S[a]
a<-i
}
}Calcularemos la desviación estandar de \(\hat{S}(t)\) con la formula de Greenwood.
\[s.e.\{\hat{S}(t)\}\approx [\hat{S}(t)]\left\{\sum_{j=1}^{k}\dfrac{d_j}{(n-r+1)(n-r)}\right\}^{1/2}\] Tomando en cuenta que para todo j \(d_j=1\)
g<-BCG[BCG$Censura.ST=="m",]
aux<-vector()
se<-vector()
for(i in 1:nrow(g)){
for(j in 1:i){
aux[j]<-1/((n-g$r[j])*(n-g$r[j]+1))
}
se[i]<-(g$S[i])*(sum(aux)^(.5))
}
g<-cbind(g,se)Calculemos los intervalos de confianza para S al 95%. Dados por
\[\hat{S}(t)\pm1.96s.e.\{\hat{S}(t)\}\]
#Calcula los limites del intervalo de confianza
g<-mutate(g, li=S-1.96*se,ls=S+1.96*se)
#Une los datos de g con los datos completos
BCG<-full_join(BCG,select(g,r,se,li,ls),by="r")
#imprime la tabla de resultados
print(xtable(BCG),type="html")| Patient | Tratament.Received | Survival.Time | Censura.ST | r | sup_int | S | se | li | ls | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 17.00 | 2.00 | 6.90 | m | 1 | 0.95 | 0.95 | 0.05 | 0.85 | 1.05 |
| 2 | 26.00 | 2.00 | 7.70 | m | 2 | 0.94 | 0.89 | 0.07 | 0.76 | 1.03 |
| 3 | 30.00 | 2.00 | 7.80 | s | ||||||
| 4 | 12.00 | 2.00 | 8.00 | m | 4 | 0.94 | 0.84 | 0.09 | 0.67 | 1.01 |
| 5 | 28.00 | 2.00 | 8.20 | s | ||||||
| 6 | 29.00 | 2.00 | 8.20 | s | ||||||
| 7 | 21.00 | 2.00 | 8.30 | m | 7 | 0.92 | 0.77 | 0.10 | 0.58 | 0.97 |
| 8 | 22.00 | 2.00 | 10.80 | s | ||||||
| 9 | 18.00 | 2.00 | 11.00 | s | ||||||
| 10 | 23.00 | 2.00 | 12.20 | s | ||||||
| 11 | 24.00 | 2.00 | 12.50 | s | ||||||
| 12 | 27.00 | 2.00 | 14.80 | s | ||||||
| 13 | 16.00 | 2.00 | 16.00 | s | ||||||
| 14 | 15.00 | 2.00 | 18.10 | s | ||||||
| 15 | 14.00 | 2.00 | 21.40 | s | ||||||
| 16 | 20.00 | 2.00 | 23.00 | s | ||||||
| 17 | 25.00 | 2.00 | 24.40 | m | 17 | 0.67 | 0.52 | 0.22 | 0.08 | 0.95 |
| 18 | 19.00 | 2.00 | 24.80 | s | ||||||
| 19 | 13.00 | 2.00 | 26.90 | s |
Grafia de superviencia para el tratamiento BCG
plot(g$S~g$Survival.Time,data=g,type="s"
,xlim=c(0,max(g$Survival.Time)),ylim=c(0,max(g$ls))
,xlab="Survival Time",ylab="S(t)")
lines(g$li~g$Survival.Time,data=g,type="s")
lines(g$ls~g$Survival.Time,data=g,type="s")
Sin distinción de tratamientos
#Filtramos los datos para el tratamiento BCG
BCG<-sup
#Calculamos r "posicion de los datos sin censura"
r<-vector()
for(i in 1:nrow(BCG)){
if(BCG$Censura.ST[i]=="m"){
r[i]<-i
}else{
r[i]<-NA
}
}
BCG<-cbind(BCG,r)Calcularemos la proporcion de pacientes que sobrevie a cada intervalo como
\[sup\_int=\dfrac{n-r}{n-r+1}\]
donde \(n\) es el tamaño de la muestra y \(r\) es el rango no censurado.
#Calcular el tamaño de la muestra
n=nrow(BCG)
#Calcular la proporcion de supervivienca por intervalo.
BCG<-mutate(BCG,sup_int=(n-r)/(n-r+1))Ahora calcularemos \[\hat{s}(t)=\hat{s}(t-1)*sup\_int_t\]
denotaremos en la tabla a \(\hat{s}(t)\) simplemente como “S”.
BCG$S<-NA
for(i in 1:n){
if((!is.na(BCG$r[i]))&(BCG$r[i]==1)){
BCG$S[i]<-BCG$sup_int[i]
a<-i
}
if((!is.na(BCG$r[i]))&(!BCG$r[i]==1)){
BCG$S[i]<-BCG$sup_int[i]*BCG$S[a]
a<-i
}
}Calcularemos la desviación estandar de \(\hat{S}(t)\) con la formula de Greenwood.
\[s.e.\{\hat{S}(t)\}\approx [\hat{S}(t)]\left\{\sum_{j=1}^{k}\dfrac{d_j}{(n-r+1)(n-r)}\right\}^{1/2}\] Tomando en cuenta que para todo j \(d_j=1\)
g<-BCG[BCG$Censura.ST=="m",]
aux<-vector()
se<-vector()
for(i in 1:nrow(g)){
for(j in 1:i){
aux[j]<-1/((n-g$r[j])*(n-g$r[j]+1))
}
se[i]<-(g$S[i])*(sum(aux)^(.5))
}
g<-cbind(g,se)Calculemos los intervalos de confianza para S al 95%. Dados por
\[\hat{S}(t)\pm1.96s.e.\{\hat{S}(t)\}\]
#Calcula los limites del intervalo de confianza
g<-mutate(g, li=S-1.96*se,ls=S+1.96*se)
#Une los datos de g con los datos completos
BCG<-full_join(BCG,select(g,r,se,li,ls),by="r")
#imprime la tabla de resultados
print(xtable(BCG),type="html")| Patient | Tratament.Received | Survival.Time | Censura.ST | r | sup_int | S | se | li | ls | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2.00 | 1.00 | 3.90 | m | 1 | 0.97 | 0.97 | 0.03 | 0.90 | 1.03 |
| 2 | 4.00 | 1.00 | 5.40 | m | 2 | 0.97 | 0.93 | 0.05 | 0.84 | 1.02 |
| 3 | 17.00 | 2.00 | 6.90 | m | 3 | 0.96 | 0.90 | 0.05 | 0.79 | 1.01 |
| 4 | 26.00 | 2.00 | 7.70 | m | 4 | 0.96 | 0.87 | 0.06 | 0.75 | 0.99 |
| 5 | 30.00 | 2.00 | 7.80 | s | ||||||
| 6 | 7.00 | 1.00 | 7.90 | m | 6 | 0.96 | 0.83 | 0.07 | 0.70 | 0.97 |
| 7 | 12.00 | 2.00 | 8.00 | m | 7 | 0.96 | 0.80 | 0.07 | 0.65 | 0.94 |
| 8 | 28.00 | 2.00 | 8.20 | s | ||||||
| 9 | 29.00 | 2.00 | 8.20 | s | ||||||
| 10 | 21.00 | 2.00 | 8.30 | m | 10 | 0.95 | 0.76 | 0.08 | 0.60 | 0.92 |
| 11 | 3.00 | 1.00 | 10.50 | m | 11 | 0.95 | 0.72 | 0.08 | 0.56 | 0.89 |
| 12 | 22.00 | 2.00 | 10.80 | s | ||||||
| 13 | 18.00 | 2.00 | 11.00 | s | ||||||
| 14 | 23.00 | 2.00 | 12.20 | s | ||||||
| 15 | 24.00 | 2.00 | 12.50 | s | ||||||
| 16 | 27.00 | 2.00 | 14.80 | s | ||||||
| 17 | 16.00 | 2.00 | 16.00 | s | ||||||
| 18 | 9.00 | 1.00 | 16.60 | s | ||||||
| 19 | 8.00 | 1.00 | 16.90 | s | ||||||
| 20 | 11.00 | 1.00 | 17.10 | s | ||||||
| 21 | 15.00 | 2.00 | 18.10 | s | ||||||
| 22 | 5.00 | 1.00 | 19.50 | m | 22 | 0.89 | 0.64 | 0.11 | 0.43 | 0.85 |
| 23 | 14.00 | 2.00 | 21.40 | s | ||||||
| 24 | 20.00 | 2.00 | 23.00 | s | ||||||
| 25 | 6.00 | 1.00 | 23.80 | s | ||||||
| 26 | 25.00 | 2.00 | 24.40 | m | 26 | 0.80 | 0.51 | 0.14 | 0.23 | 0.79 |
| 27 | 19.00 | 2.00 | 24.80 | s | ||||||
| 28 | 13.00 | 2.00 | 26.90 | s | ||||||
| 29 | 1.00 | 1.00 | 33.70 | s | ||||||
| 30 | 10.00 | 1.00 | 33.70 | s |
Grafia de superviencia para el tratamiento BCG
plot(g$S~g$Survival.Time,data=g,type="s"
,xlim=c(0,max(g$Survival.Time)),ylim=c(0,max(g$ls))
,xlab="Survival Time",ylab="S(t)")
lines(g$li~g$Survival.Time,data=g,type="s")
lines(g$ls~g$Survival.Time,data=g,type="s")