PRACTICA DE SUPERVIVENCIA
Cruz Mateo David
Cargamos las siguientes librerias
library(survival)
library(survminer)## Loading required package: ggplot2## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomalibrary(flexsurv)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.4 v purrr 0.3.4
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()Evento Ser rearestado despues de salir de prision
El seguimiento de los casos duro tres años
Importamos datos
henning <- read.csv("C:/Users/doria/OneDrive/Escritorio/Cuestionario_2/PRACTICA DE SUPERVIVENCIA/data/Henning.txt", sep="")Caracteristicas
summary(henning)## id months censor personal
## Min. : 1.00 Min. : 0.06571 Min. :0.0000 Min. :0.0000
## 1st Qu.: 49.25 1st Qu.: 4.27926 1st Qu.:0.0000 1st Qu.:0.0000
## Median : 97.50 Median :10.23409 Median :0.0000 Median :0.0000
## Mean : 97.50 Mean :13.80647 Mean :0.4536 Mean :0.3144
## 3rd Qu.:145.75 3rd Qu.:19.92608 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :194.00 Max. :36.00000 Max. :1.0000 Max. :1.0000
## property cage
## Min. :0.0000 Min. :-12.487
## 1st Qu.:1.0000 1st Qu.: -5.935
## Median :1.0000 Median : -1.916
## Mean :0.8144 Mean : 0.000
## 3rd Qu.:1.0000 3rd Qu.: 3.998
## Max. :1.0000 Max. : 32.868id numero de ientificacion del preso del 1 al 194
months cuantos meses tarda la persona en ser rearrestada de 0 a 36
las variables censor, personal, property son indicadoras
censor indica si se censuro o no la persona ,1 obserbaste y no ocurrio el evento, 0 si ocurrio el evento
naturaleza del registro criminal
personal, 1 ,antecednetes criminales con respecto a personas
property, 1, antecednetes criminales con respecto hacia propiedades
cage, edad centrada, diferencia entre edad y edad promedio
Evento de interes ser arrestado
Creamos la variable event (evento)
1 si fue rearrestado
0 si no fue rearrestado
henning <- henning %>%
mutate(event = case_when(
censor == 0 ~ 1,
TRUE ~ 0
))Visualizamos de nuevo encabezados
head(henning)## id months censor personal property cage event
## 1 1 0.06570842 0 1 1 -1.675198 1
## 2 2 0.13141684 0 0 1 -10.482864 1
## 3 3 0.22997947 0 1 1 -4.426738 1
## 4 4 0.29568789 0 0 1 -11.328860 1
## 5 5 0.29568789 0 1 1 -7.164589 1
## 6 6 0.32854209 0 1 0 -2.868901 1Agrupamos
henning %>%
group_by(event) %>% #agrupamos por evento
summarise(avg_months = mean(months), # meses promedio de rearresto
avg_personal = mean(personal), #promedio registro de incidentes con personas
avg_property = mean(property),#promedio de registro de incidentes con propiedades
avg_cage = mean(cage))# promedio de edad centrada en la que se ometen los incidentes## # A tibble: 2 x 5
## event avg_months avg_personal avg_property avg_cage
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 19.6 0.25 0.693 3.25
## 2 1 9.02 0.368 0.915 -2.70Observamos para que una persona vuele a ser rearrestada en promedio son 9 meses, el numero promedio de registros es mayor , tanto en incidentes con personas, como en incidentes con propiedades y la edad media delictiva es mucho menor.
Funcion intro_surv
Contruimos una funcion que nos indica dada una muestra de tamaño n si la observacion fue censurada, o el evento ocurrio ## intro_surv
intro_surv<- function(data=NULL,n=NULL){
muestra=sample_n(data,n) # seleccionamos una muestra de n
muestra$id=seq(1,n)# renombramos los id
time=muestra$months# meses de la muestra
event=muestra$event# si fue rearrestado
ggplot(muestra)+
geom_point(aes(x =id, y = time, shape= as.factor(event)))+
geom_linerange(aes(x =id, ymin = 0, ymax = time), linetype = "dotdash")+
theme_bw()+
labs(title ='Survival times plot',
shape = 'Evento',
x = 'Muestra',
y='Tiempo')+
coord_flip() #giramos los ejes
}llamamos la funcion, n=10
intro_surv(data=henning,n=10)
En su mayoria se observo que habia mas sujetos rearrestados , estos ocurrian entre los 5 y 15 meses , en el estudio, entonces la censura es menor
Objeto de sobrevivencia
Regresion entre el logaritmo de la supervivencia y los tiempos
# creamos el objeto de tipo sobrevencia
# sfit ajuste de supervivencia
# Tiempo = months y evento = event
sfit <- survfit(Surv(months, event) ~ 1, data = henning)Graficamos
Probabilidad sobrevivencia S(t)
ggsurvplot(sfit, data = henning,
size=0.4,
palette = c("darkblue"),
title = "Reincidencia criminal",
xlab = "Tiempo",
ylab = "Probabilidad de supervivencia, S(t)",
legend.title = "Persona y Propiedad")
Obsevamos una tasa baja debido a que la probabilidad de ser rearrestados durante los primeros 10 meses es mayor, Apartir de mes 20 la probabilidad de ser rearrestado parece permanecer constante en torno al 30 porciento ### Riesgo acumulada de hazard H(t)
ggsurvplot(sfit, data = henning,
fun = "cumhaz", # riesgo acumulada
size=0.4,
palette = c("darkblue"),
title = "Reincidencia criminal",
xlab = "Tiempo",
ylab = "Riesgo Acumulada H(t)",
legend.title = "Persona y Propiedad")
Mide la frecuencia con que ocurren los fallos en el tiempo y en el análisis de residuos para el ajuste de algunos modelos.
Estimacion por medio de Modelos parametricos
Regresion entre el logaritmo de la supervivencia y los tiempos
Usamos la siguiente funcion de Diagnostico paraver si los datos podrian ajustarse a una distribucion conocida exponencial, gompertz, Weibull, Log-logistica, Lognormal.
diagnos_surv <- function(times= NULL, surv=NULL) {# parmetros
cuantil <- qnorm(1-exp(-log(surv)))# para la lognormal
# creamos dataframe para el ggplot2
datos=data.frame(times = times, surv = surv,cuantil= cuantil)
g1<-ggplot(datos)+# exponecial
geom_point(aes(x = log(surv), y = times))+
geom_smooth(aes(x = log(surv), y = times), method = "lm")+
labs(title="Diagnostico exponencial",subtitle = "Se espera una linea negativa", x="S(t)", y="t")+
theme_bw()
g2 <- ggplot(datos)+# Gompertz
geom_point(aes(x = log(-log(surv)+1), y = times))+
geom_smooth(aes(x = log(-log(surv)+1), y = times), method = "lm")+
labs(title="Diagnostico Gompertz",subtitle = "Se espera una linea positiva", x="S(t)", y="t")+
theme_bw()
g3 <- ggplot(datos)+ # Weibull
geom_point(aes(x = log(-log(surv)), y = times))+
geom_smooth(aes(x = log(-log(surv)), y = times), method = "lm")+
labs(title="Diagnostico weibull",subtitle = "Se espera una linea positiva", x="S(t)", y="t")+
theme_bw()
g4 <- ggplot(datos)+# Log-logistic
geom_point(aes(x = log(exp(-log(surv)-1)), y = times))+
geom_smooth(aes(x = log(exp(-log(surv)-1)), y = times), method = "lm")+
labs(title="Diagnostico Log-logistic",subtitle = "Se espera una linea positiva", x="S(t)", y="t")+
theme_bw()
g5 <- ggplot(datos)+# Lognormal
geom_point(aes(x =cuantil , y = times))+
geom_smooth(aes(x = cuantil, y = times), method = "lm")+
theme_bw()
return(ggarrange(g1,g2,g3,g4,ncol=2,nrow=2))
}llamamos la funcion
diagnos_surv(times=sfit$time,surv=sfit$surv)## Warning in qnorm(1 - exp(-log(surv))): NaNs produced## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Por lo tanto vemos que los modelos que podriamos usar serian una Exponencial, Gommpertz o Log-logistico
Pero no podriamos usar una Weibull o Log normal
Seguimiento, relacionado con personas
0 Aquellos que no cometen crimenes en contra de personas 1 Aquellos que cometen crimenes en contra de personas
# objeto de ajuste
# dividimos las poblaciones
sfit_personal <- survfit(Surv(months, event) ~ personal, data = henning)
ggsurvplot(sfit_personal,
data = henning,
palette = c("darkgreen","red"),
size=0.4,
title = "Reincidencia criminal",
xlab = "Tiempo",
ylab = "Supervivencia S(t)",
legend.title = "Personal")
Aquellos que cometen crimenes en contra de personas son arrestados mas pronto.
Seguimiento, relacionado con propiedad
0 Aquellos que no cometen crimenes en contra de propiedades 1 Aquellos que cometen crimenes en contra de propiedades
# dividimos la poblacion
sfit_property <- survfit(Surv(months, event) ~ property, data = henning)
ggsurvplot(sfit_property,
data = henning,
size=0.4,
palette = c("darkgreen","red"),
title = "Reincidencia criminal",
xlab = "Tiempo",
ylab = "Supervivencia S(t)",
legend.title = "Propiedad")
Aquellos que cometen crimenes en contra de Propiedad son rearrestados mas pronto incluso que aquellos que cometen crimenes encontra de personas,vemos un decrecimiento mas rapido, tal vez debido, a que el daño a propiedad la condena es menor, y por eso es las personas cometen mas crimenes de este tipo, el indide rearrestos es menor
Distancia entre las curvas
Las curvas con discretas, utilizamos la funcion aproxfun para poder Usamos la funcion max_distancia para poder medir la maxima distancia entre las dos graficas de supervivencia
max_distancia <- function(sfitp = NULL,p=NULL) {
f1 <- sfitp$strata[1] # Strata 0
f2 <- length(sfitp$surv) # Strata 1
t1<-sfitp$time[1:f1]# tiempo 0
s1<-sfitp$surv[1:f1]# sobrevivencia 0
t2<-sfitp$time[(f1+1):f2]# tiempo 1
s2<-sfitp$surv[(f1+1):f2]# sobrevivencia 1
#suavizamos funciones de supervivencia y las evaluamos
s1aprox <- approxfun(s1,t1)
s2aprox <- approxfun(s2,t2)
# creamos dataframe para ggplot2
df1 <- data.frame(s=s2,ta1=s1aprox(s2), ta2=s2aprox(s2))
# diferencia entre los tiempos para las funciones
df1$distancia<-abs(df1$ta1-df1$ta2)
# Eliminamos NAs para graficar
df1 <- df1[!is.na(df1$distancia),]
# Conservamos las S(t)>p
df<-df1[df1$s>p,]
# Obtenemos la mayor de las distancias
maxima <- max(df1$distancia)
# guardamos la distancia, en un nuevo df
d<-df1[df1$distancia==maxima,]
max<-as.character(round(maxima,6))
grafica<-ggplot()+
geom_step(aes(t1,s1, color="darkred"),size=0.6)+# personal=0
geom_step(aes(t2,s2, color="darkgreen"),size=0.6)+# persoaml=1
labs(title=paste("Máxima Distancia entre","personal=1","y","personal=0", sep = " "),
x="Tiempo", y="S(t)")+
theme_bw()+ # graficamos puntos maximo
geom_hline(yintercept = d$s, linetype="dashed", color = "gray50", size=0.7)+
labs(color = "red") +
geom_vline(xintercept = d$ta1, linetype="dashed", color = "gray50", size=0.7)+
geom_vline(xintercept = d$ta2, linetype="dashed", color = "gray50", size=0.7)+
annotate("text", x =min(d$ta1,d$ta2), y = d$s, label = paste("Distancia=",max,sep =""))
return(grafica)
}Con respecto a personas
max_distancia(sfitp =sfit_personal,p=0)## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
Entre las personas que cometen crimenes relacionados con personas y las que no estan relacionadas con personas, la distancia maxima es de 11 meses, debido a que el delito es mas castigado
Con respecto a propiedad
max_distancia(sfitp =sfit_property,p=0)## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
Entre las personas que cometen crimenes relacionados con propiedades y las que no estan relacionadas con propiedades, la distancia maxima es de 22 meses, debido a que el delito es menos penado