###################################################
##  APLICACIÓN PRÁCTICA ANÁLISIS DE SUPERVIVENCIA##
###################################################

install.packages('survminer')
install.packages('prodlim')
install.packages('pec')
install.packages("ggplot2")
install.packages("ggpubr")
install.packages("survMisc")
install.packages("flexsurv")
install.packages("actuar")
install.packages("dplyr")
install.packages("ggfortify")
library(survival)
library(survminer)
library(prodlim)
library(pec)
library(KMsurv)
library(ggplot2)
library(ggpubr)
library(survMisc)
library(flexsurv)
library(actuar)
library(dplyr)
library(ggfortify)

##IMPORTAT DATOS:
#################

data('psych')
head(psych)
oldNames<-colnames(psych)
colnames(psych)<-c('Sexo','Edad','Tiempo','Estado')
dim(psych)
str(psych) #Comprobar que están ordenados
head(psych)
summary(psych)

psych$Sexo<-factor(psych$Sexo)
levels(psych$Sexo)<-c('M','F')
summary(psych)
str(psych)
classes<-sapply(psych,class)
for (i in 1:ncol(psych)) {
  if(classes[i]=='integer') psych[[i]]= as.numeric(psych[[i]])
  
}
classes=sapply(psych,class)
classes
str(psych)


## HISTOGRAMA DE LA VARIABLE TIEMPO:

hist(psych$Tiempo, freq = F,
     ylab = 'Densidad',
     xlab = 'Tiempo', main = 'Histograma variable tiempo',
     col = 'pink')
lines(density(psych$Tiempo), col ='black')
curve(dnorm(x,mean(psych$Tiempo),sd(psych$Tiempo)),
      col = 'orange',lwd=3, add = TRUE) #Normal
shapiro.test(psych$Tiempo)

##KAPLAN-MEIER:
###############
# PRIMERO SIN COVARIABLE
Surv(psych$Tiempo, psych$Estado) # muestran los tiempos en los que  
                                 # se dieron los eventos y las censuras.
                                 # Los tiempos censurado se representan con el signo +

psych.KM<-survfit(Surv(Tiempo,Estado)~1, data=psych) #~1 para no distinguir entre niveles
print(psych.KM)
summary(psych.KM)
median(psych.KM)

 # CURVA DE SUPERVIVENCIA:

ggsurvplot(fit = psych.KM, data = psych, 
           conf.int = T, 
           title = 'Curva de supervivencia',
           xlab = 'Tiempo en semanas',
           ylab = 'Probabilidad de supervivencia',
           legend.title = 'Estimación',
           legend.labs = 'Kaplan-Meier',
           surv.median.line = 'hv',
           ggtheme = theme_light())

 # ESTIMACIÓN DE LA FUNCIÓN DE RIESGO ACUMULADO:

psych.km<-survfit(Surv(Tiempo, Estado)~1, psych)
psych.Racum<-psych.km %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
psych.Racum
ggsurvplot(psych.km, fun = 'cumhaz',
           xlab = 'Tiempo en semanas',
           censor = T,
           ylab = 'Riesgo Acumulado',
           title = 'Riesgo Acumulado',
           legend.title = 'Estimación',
           legend.labs = 'Kaplan-Meier',
           ggtheme = theme_light())

# SEGUNDO CON LA COVARIABLE SEXO:

psych.KM<-survfit(Surv(Tiempo,Estado)~Sexo, data=psych) 
print(psych.KM)
summary(psych.KM)
median(psych.KM)

# TEST LOG-RANk (COMPARAR DOS FUNCIONES DE SUPERVIVENCIA)

survival::survdiff(Surv(Tiempo, Estado)~Sexo, data = psych, rho = 0) # Test log-rank

 # CURVA DE SUPERVIVENCIA

ggsurvplot(fit = psych.KM, data = psych, 
           conf.int = T, 
           title = 'Curva de supervivencia por Sexo',
           xlab = 'Tiempo en semanas',
           ylab = 'Probabilidad de supervivencia',
           legend.title = 'Género',
           legend.labs = c('Masculino', 'Femenino'),
           surv.median.line = 'hv',
           ggtheme = theme_light())
  
 # ESTIMACIÓN DE LA FUNCIÓN DE RIESGO ACUMULADO:

psych.km<-survfit(Surv(Tiempo, Estado==1)~Sexo, psych)
psych.Racum<-psych.km %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
psych.Racum
ggsurvplot(psych.km, fun = 'cumhaz',
           xlab = 'Tiempo en semanas',
           censor = T,
           ylab = 'Riesgo Acumulado',
           title = 'Riesgo Acumulado',
           legend.title = 'Sexo',
           legend.labs = c('Masculino','Femenino'),
           ggtheme = theme_light())



### NELSON AALEN###
###################

# PRIMERO SIN COVARIABLES:

psych.na<-survfit(Surv(Tiempo, Estado)~1, psych)
psych.Racum.na<-psych.na %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
psych.Racum.na

qplot(time, CumHaz, data = psych.Racum.na,
      geom = 'step',
      xlab = 'Tiempo en semanas',
      ylab = 'Riesgo Acumulado',
      main = 'Riesgo Acumulado con Nelson Aalen')

# SEGUNDO CON COVARIABLES:

psych.na<-survfit(Surv(Tiempo, Estado==1)~Sexo, psych)
psych.Racum.na<-psych.na %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
psych.Racum.na

qplot(time, CumHaz, col = strata, data = psych.Racum.na,
      geom = 'step',
      xlab = 'Tiempo en semanas',
      ylab = 'Riesgo Acumulado',
      main = 'Riesgo Acumulado con Nelson Aalen')