Primer ejemplo

Retorno al uso de drogas

Información inicial

El objetivo de este estudio es analizar el tiempo de retorno al uso de drogas de pacientes enrolados en dos programas residenciales de tratamiento, que difieren en su periodo de duracion treat:

  • treat = 0 es el programa corto
  • treat = 1 es el programa largo

Los pacientes fueron asignados aleatoriamente a dos sitios diferentes:

  • site = 0 es el sitio A

  • site = 1 es el sitio B

  • age es la edad a la que se enrolaron.

  • hercoc indica el consumo de heroína o cocaína en los tres meses anteriores con codigos:

    • hercoc=1: uso de heroina y cocaina
    • hercoc=2: solo heroina
    • hercoc=3: Solo cocaina
    • hercoc=4: No uso de heroina ni de cocaina

  • ndrugtx indica el número de tratamientos previos de drogas

  • race es la raza de los pacientes con:

    • race=0: Blanca
    • race=1: No blanca

  • ivhx indica la historia de uso de drogas intravenosas:

    • ivhx=1. Nunca
    • ivhx=2: Previo
    • ivhx=3: Reciente

  • time contiene el tiempo de retorno al uso de drogas

  • censor indica si el sujeto ha retornado o no al uso con codigo:

    • censor=1: retorno al uso de drogas
    • censor=0: no retorno

Análisis de los datos

Comenzamos por cargar y ver los datos

uis <- read.csv("uis.csv")
head(uis)
##   id age  beck hercoc ivhx ndrugtx race treat site time censor
## 1  1  39  9.00      4    3       1    0     1    0  188      1
## 2  2  33 34.00      4    2       8    0     1    0   26      1
## 3  3  33 10.00      2    3       3    0     1    0  207      1
## 4  4  32 20.00      4    3       1    0     0    0  144      1
## 5  5  24  5.00      2    1       5    1     1    0  551      0
## 6  6  30 32.55      3    3       1    0     1    0   32      1
dim(uis)
## [1] 575  11
attach(uis)

Vemos un gráfico de barras respecto a la variable hercoc:

T1 <- table(hercoc)
names(T1) <- c("H&C","H","C","NINGUNA")
b <- barplot(T1, col=rainbow(15:18), main="Gráfica de barras: Tipo de drogas", col.main="red", ylab="Frecuencia absoluta", col.lab="green")
text(x=b,y=c(T1), labels=c(T1), pos=3, col="darkgreen", cex=1.5, xpd=TRUE)

de donde se observa que hay un mayor número de personas que previamente no han consumido ni heroína ni cocaína. Podemos ver el número de fallas y censuras en nuestros datos:

T2 <- table(censor)
names(T2) <- c("censura","falla")
b <- barplot(T2, col=rainbow(2), main="Censura", col.main="red", ylab="Frecuencia absoluta", col.lab="green")
text(x=b, y=c(T2), labels=c(T2), pos=3, col="darkgreen", cex=1.5, xpd=TRUE)

Por otro lado tenemos:

T3 <- table(ndrugtx)
b <- barplot(T3, col=rainbow(30), main="Intentos previos de dejar las drogas", col.main="red", ylab="Frecuencia absoluta", col.lab="green")
text(x=b, y=c(T3), labels=c(T3), pos=3, col="darkgreen", cex=1.5, xpd=TRUE)

vemos que ha habido al menos una persona que ha intentado dejar las drogas más de 31 veces; hay 80 personas que están por primera vez en tratamiento para dejar las drogas.

T5 <- table(race)
names(T5) <- c("White","Non White")
b <- barplot(T5,col=rainbow(2), main="Raza", col.main="red", ylab="Frecuencia absoluta", col.lab="green")
text(x=b,y=c(T5), labels=c(T5), pos=3, col="darkgreen", cex=1.5, xpd=TRUE)

vemos que hay más personas de raza blanca dentro del estudio.

T6 <- table(site)
names(T6) <- c("Site A","Site B")
b <- barplot(T6,col=rainbow(2), main="Sitio", col.main="red", ylab="Frecuencia absoluta", col.lab="green")
text(x=b,y=c(T6), labels=c(T6), pos=3, col="darkgreen", cex=1.5, xpd=TRUE)

Se han asignado aleatoriamente más persona al sitio a respecto del B.

Finalmente vemos

hist(time, probability=T, col="green", border="white", main="Histograma del tiempo de supervencia") 
lines(density(time), col="red",lwd=2)

sin embargo, no es correcto realizar un histograma del tiempo de supervivencia puesto que el tiempo registrado en los datos, no siempre es el tiempo de supervivencia, pues puede ser también (por ejemplo) el tiempo en el cual se observó a in individuo hasta la censura. En otras palabras, no sabemos cuáles tiempos corresponden a fallas y cuáles a censuras.

Tabla actuarial de vida: Intervalos de 6 meses (180 dias)

Comenzamos creandando la tabla actuarial de vida como sigue:

attach(uis)
uis$interval <- floor(uis$time/180)*180

# Tabla para registrar fallas (1) y censuras (0) 
tab <- table(uis$interval, uis$censor)

# Tiempo observado a_0, a_1,... con longitud de 180 días
intEndpts <- (0:6)*180

# Número de sujetos que inicia en el estudio
ntotal <- dim(uis)[1]

# Vector de números referentes a las censuras
cens <- tab[,1]

# Vector de números referentes a las fallas
fallas <- tab[,2]

# Definición de la tabla de vida
fitLifeTable <- lifetab(tis = intEndpts, ninit = ntotal, nlost = cens, nevent = fallas)

round(fitLifeTable,8)
##          nsubs nlost nrisk nevent      surv        pdf     hazard    se.surv
## 0-180      575     0 575.0    296 1.0000000 0.00285990 0.00385116 0.00000000
## 180-360    279     0 279.0    128 0.4852174 0.00123671 0.00330749 0.02084233
## 360-540    151    16 143.0     37 0.2626087 0.00037749 0.00165105 0.01835142
## 540-720     98    90  53.0      3 0.1946610 0.00006121 0.00032362 0.01665955
## 720-900      5     4   3.0      0 0.1836424 0.00000000 0.00000000 0.01688753
## 900-1080     1     1   0.5      0 0.1836424         NA         NA 0.01688753
##              se.pdf  se.hazard
## 0-180    0.00011579 0.00020997
## 180-360  0.00009638 0.00027909
## 360-540  0.00005959 0.00026842
## 540-720  0.00003472 0.00018677
## 720-900         NaN        NaN
## 900-1080         NA         NA
names(fitLifeTable)
##  [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"      
##  [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

donde vemos en la columna surv el cálculo estimado de la función de supervivencia. Podemos graficar la función de supervivencia asociada a esta tabla:

x = rep(intEndpts,rep(2,7))[2:13]
y = rep(fitLifeTable$surv, rep(2,6))
plot(x, y, type="l", col="darkred", xlab="Time", ylab="Supervivencia (Life Table)", col.lab="darkblue", lwd=2)
title(main="Duración de recaída", cex=0.6, col.main="darkblue")

y la función de riesgo:

y = rep(fitLifeTable$hazard, rep(2,6))
plot(x, y, type="l", col="darkgreen", xlab="Time", ylab="Riesgo (Life Table)",col.lab="darkblue",lwd=2)
title(main="Duración de recaída", cex=.6,col.main="orange")

Ahora, obtenemos la gráfica de la función de riesgo y de supervivencia del tratamiento corto:

# Obtenemos la tabla de vida filtrada al tratamiento A
uis0 <- subset(uis,treat==0)
attach(uis0)

uis0$interval <- floor(uis0$time/180)*180
tab <- table(uis0$interval,uis0$censor)
intEndpts <- (0:5)*180
ntotal <- dim(uis0)[1]
cens <- tab[,1]
fallas <- tab[,2]
fitLifeTable0 <- lifetab(tis = intEndpts, ninit=ntotal, nlost=cens, nevent=fallas)
round(fitLifeTable0,4)
##         nsubs nlost nrisk nevent   surv    pdf hazard se.surv se.pdf se.hazard
## 0-180     289     0   289    169 1.0000 0.0032 0.0046  0.0000  2e-04     3e-04
## 180-360   120     0   120     57 0.4152 0.0011 0.0035  0.0290  1e-04     4e-04
## 360-540    63    10    58     12 0.2180 0.0003 0.0013  0.0243  1e-04     4e-04
## 540-720    41    36    23      1 0.1729 0.0000 0.0002  0.0225  0e+00     2e-04
## 720-900     4     4     2      0 0.1654     NA     NA  0.0227     NA        NA
# Función de riesgo
par(mfrow=c(1,2))
x = rep(intEndpts,rep(2,6))[2:11]
y = rep(fitLifeTable0$surv,rep(2,5))
plot(x, y, type="l", col="darkred", xlab="Time", ylab="Supervivencia (Life Table)", col.lab="darkblue", lwd=2)
title(main = "Duración de recaída tratamiento corto",  cex=.6, col.main="darkblue")

# Función de riesgo
y=rep(fitLifeTable0$hazard, rep(2,5))
plot(x, y, type="l", col="darkgreen", xlab="Time", ylab="Riesgo (Life Table)", col.lab="darkblue", lwd=2)
title(main="Duración de recaída tratamiento corto",  cex=.6, col.main="orange")

y análogamente respecto al tratamiento largo:

# Obtenemos la tabla de vida filtrada al tratamiento B
uis1 <- subset(uis,treat==1)
attach(uis1)

uis1$interval <- floor(uis1$time/180)*180
tab <- table(uis1$interval,uis1$censor)
intEndpts <- (0:5)*180
ntotal <- dim(uis1)[1]
cens <- tab[,1]
fallas <- tab[,2]
fitLifeTable1 <- lifetab(tis = intEndpts, ninit=ntotal, nlost=cens, nevent=fallas)
round(fitLifeTable1,4)
##         nsubs nlost nrisk nevent   surv    pdf hazard se.surv se.pdf se.hazard
## 0-180     286     0 286.0    127 1.0000 0.0025 0.0032  0.0000  2e-04     3e-04
## 180-360   159     0 159.0     71 0.5559 0.0014 0.0032  0.0294  1e-04     4e-04
## 360-540    88     6  85.0     25 0.3077 0.0005 0.0019  0.0273  1e-04     4e-04
## 540-720    57    54  30.0      2 0.2172 0.0001 0.0004  0.0245  1e-04     3e-04
## 720-900     1     1   0.5      0 0.2027     NA     NA  0.0250     NA        NA
# Función de supervivencia
par(mfrow=c(1,2))
x=rep(intEndpts,rep(2,6))[2:11]
y=rep(fitLifeTable1$surv,rep(2,5))
plot(x, y, type="l", col="darkred", xlab="Time", ylab="Supervivencia (Life Table)", col.lab="darkblue", lwd=2)
title(main = "Duración de recaída tratamiento largo",  cex=.6, col.main="darkblue")

# Función de riesgo
y=rep(fitLifeTable1$hazard, rep(2,5))
plot(x, y, type="l", col="darkgreen", xlab="Time", ylab="Riesgo (Life Table)", col.lab="darkblue", lwd=2)
title(main="Duración de recaída tratamiento largo",  cex=.6, col.main="orange")

Comparacion de supervivencias entre tratamiento largo y corto

x=rep(intEndpts,rep(2,6))[2:11]
y=rep(fitLifeTable0$surv,rep(2,5))
plot(x, y, type="l", col="darkred", xlab="Time", ylab="", col.lab="darkblue", lwd=2,yaxt="n")
title(main = "Comparación de supervivencias: tratamiento corto vs. largo", cex=.6,col.main="darkblue")

par(new=T)
x=rep(intEndpts,rep(2,6))[2:11]
y=rep(fitLifeTable1$surv,rep(2,5))
plot(x, y, type="l", lty=2, col="darkgreen", xlab="Time", ylab="", col.lab="darkblue", lwd=2, yaxt="n")
legend("topright",c("Tratamiento corto","Tratamiento largo"), col=c("darkred","darkgreen"), pch=15, bty="n")

vemos que el tratamiento largo tiene una mayor supervivencia respecto al tratamiento corto.

Comparacion de riesgos entre tratamiento largo y corto:

x=rep(intEndpts,rep(2,6))[2:11]
y=rep(fitLifeTable0$hazard,rep(2,5))
plot(x, y, type="l", col="darkred", xlab="Time", ylab="",col.lab="darkblue", lwd=2, yaxt="n")
title(main = "Comparación de riesgos: tratamiento corto vs. largo", cex=.6, col.main="darkblue")

par(new=T)
x=rep(intEndpts,rep(2,6))[2:11]
y=rep(fitLifeTable1$hazard,rep(2,5))
plot(x, y, type="l",lty=2, col="darkgreen", xlab="Time", ylab="",col.lab="darkblue",lwd=2,yaxt="n")
legend("topright",c("Tratamiento 0","Tratamiento 1"), col=c("darkred","darkgreen"), pch=15, bty="n")

vemos que el tratamiento largo presenta mayor riesgo respecto al tratamiento corto.