## ------------------------------------------------------------------------ ##
##  R scripts examples pour sharing                                         ##
## ------------------------------------------------------------------------ ##
##  Scripts       Estimation of a competing risks                           ## 
## ------------------------------------------------------------------------ ##
##  Author        Quentin Guibert                                           ##
## ------------------------------------------------------------------------ ##
##  Version       01 - 03/01/2013                                           ##
## ------------------------------------------------------------------------ ##
#----------------------------------------------------------

#----------------------------------------------------------
# Initialization
#----------------------------------------------------------
# Clean memory
rm(list=ls())

# Define the path
dossier <- getwd()
# Load personal functions
source(paste(dossier, "/R/ToolsMstate.R", sep=""))
source(paste(dossier, "/R/GraphicTools.R", sep=""))
source(paste(dossier, "/R/ToolsRetraitementMstate.R",sep=""))

#----------------------------------------------------------
# Load packages
library(mstate)
library(dplyr)

#----------------------------------------------------------
# Load data (toy dataset)
load( file = paste(dossier, "/data/data_example.RData", sep=""))

#----------------------------------------------------------
# Define transition matrix for a competing risks model
tmat <- transMat(x = list(c(2,3,4,5,6,7), c(), c(),  c(), c(),c(),c()), names = c("Cotisant", "Dépendant_1", "Dépendant_2",
                                                                                 "Dépendant_3","Dépendant_4","Décés","Résilié"))
attr(t,"trans") <- tmat
attr(t,"class") <- c("msdata","data.frame")

#----------------------------------------------------------
# Parameters
conf <- 95 # Level of the confidence interval
coef_conf <- qnorm((1-conf/100)/2)
AgeMin <- 65 # Age min
AgeMax <- 90 # Age max
DureeMax <- 60

#----------------------------------------------------------
#1. Basic mstate model - Transition intensities from the state "Cotisant"
#----------------------------------------------------------
# Estimate cox model without covariate
c_a <- coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data=t, method = "breslow")

# Estimate transition intensities
msf_a <- msfit(object = c_a, vartype = "aalen", trans = tmat, variance=TRUE)
#----------------------------------------------------------
# Sex effect
#Cox
c_a_s <- estimcox(t, "Sexe", TRUE, "clock.forward")

# Check coefficients
summary(c_a_s)

# Test the PH assumption
cox.zph(c_a_s)
# ==> the PH assumption is not satisfyed ==> we consider two subpopulation: Male and Female

# Estimation transition intensities for males and females
msf_a_m <- estimintens(t, tmat , c("Sexe"), c("M"))
msf_a_f <- estimintens(t, tmat, c("Sexe"), c("F"))

# Save results
save(c_a_s, file = paste(dossier,"/data/c_a_s.RData",sep=""))
save(msf_a_m, file = paste(dossier,"/data/msf_a_m.RData",sep=""))
save(msf_a_f, file = paste(dossier,"/data/msf_a_f.RData",sep=""))
#----------------------------------------------------------
# Plot transition intensities with confidence interval
# Create a data.frame with the results
msf_a_graph <- msf_a$Haz
msf_a_graph$Haz_m <- msf_a_m$Haz$Haz
msf_a_graph$Haz_f <- msf_a_f$Haz$Haz
msf_a_graph$se <- subset(msf_a$varHaz,trans1==trans2)$varHaz^0.5
msf_a_graph$se_m <- subset(msf_a_m$varHaz,trans1==trans2)$varHaz^0.5
msf_a_graph$se_f <- subset(msf_a_f$varHaz,trans1==trans2)$varHaz^0.5
msf_a_graph$borneup <- msf_a_graph$Haz*exp(-coef_conf*msf_a_graph$se/msf_a_graph$Haz)
msf_a_graph$bornedown <- msf_a_graph$Haz*exp(coef_conf*msf_a_graph$se/msf_a_graph$Haz)
msf_a_graph$borneup_m <- msf_a_graph$Haz_m*exp(-coef_conf*msf_a_graph$se_m/msf_a_graph$Haz_m)
msf_a_graph$bornedown_m <- msf_a_graph$Haz_m*exp(coef_conf*msf_a_graph$se_m/msf_a_graph$Haz_m)
msf_a_graph$borneup_f <- msf_a_graph$Haz_f*exp(-coef_conf*msf_a_graph$se_f/msf_a_graph$Haz_f)
msf_a_graph$bornedown_f <- msf_a_graph$Haz_f*exp(coef_conf*msf_a_graph$se_f/msf_a_graph$Haz_f)
msf_a_graph$trans <- as.factor(msf_a_graph$trans)
msf_a_graph$trans <- factor(msf_a_graph$trans,labels=c("Pathologie 1","Pathologie 2","Pathologie 3","Pathologie 4","Deces","Resiliation"))

# Age interval 65-90
msf_a_graph <- subset(msf_a_graph,msf_a_graph$time>=AgeMin & msf_a_graph$time<=AgeMax) 
# Create a ggplot objet
plot_a <- ggplot(msf_a_graph,aes(x=time,y=Haz),colour="black")+
  geom_line()+
  geom_ribbon(aes(x=time,ymin=bornedown,ymax=borneup),colour="darkgrey",alpha=0.2,linetype=0)
plot_a <- plot_a+scale_x_continuous("Age")+
  scale_y_continuous("Intensité")+
  facet_wrap(~trans,3,2,scales="free_y")
# Print the results
png(file="figure/intensity.png", width  = 8000, height = 4800,res=750,antialias = "cleartype" )
plot_a
dev.off()

# Plot by sex
plot_a_s <- ggplot(msf_a_graph,aes(x=time,y=Haz_m))+
  geom_line(aes(color="Homme"))+
  geom_ribbon(aes(x=time,ymin=bornedown_m,ymax=borneup_m),colour="darkgrey",alpha=0.2,linetype=0)+
  geom_line(aes(x=time,y=Haz_f,color="Femme"))+
  geom_ribbon(aes(x=time,ymin=bornedown_f,ymax=borneup_f),colour="darkgrey",alpha=0.2,linetype=0)
plot_a_s <- plot_a_s+scale_x_continuous("Age")+
  scale_y_continuous("Intensité")+
  scale_colour_manual("", 
                      breaks = c("Homme", "Femme"),
                      values = c("blue","red"))+
  facet_wrap(~trans,3,2,scales="free_y")
plot_a_s <- plot_a_s+scale_colour_hue("Sexe")+theme(legend.position="bottom",
                                                 legend.background = element_rect(colour = "black"))+
  theme(legend.key = element_rect(colour = "white"))
png(file="figure/intensity_sex.png", width  = 8000, height = 4800,res=750,antialias = "cleartype" )
plot_a_s
dev.off()


#----------------------------------------------------------
# Calculate annual transition probabilities

pt_a <- probtrans(msf_a,predt=AgeMin,method="aalen",variance=T,direction="forward") # Transition probabilities P(65,t)
qx_a <- probtrans2(msf_a,predt=AgeMin,pas=1,method="aalen",variance=T,direction="ratecurve") # incidence rates P(x, x + 1)

#Sauvegarde des résultats
save(pt_a,file = paste(dossier,"/data/pt_a.RData",sep=""))
save(qx_a,file = paste(dossier,"/data/pt_qx_a.RData",sep=""))

#----------------------------------------------------------
# plot (crude) transition probabilities with confidence intervals

time=rep(pt_a[[1]]$time,8)
prob=c(pt_a[[1]]$pstate1,
       pt_a[[1]]$pstate2,
       pt_a[[1]]$pstate3,
       pt_a[[1]]$pstate4,
       pt_a[[1]]$pstate5,
       pt_a[[1]]$pstate6,
       pt_a[[1]]$pstate7,
       pt_a[[1]]$pstate2+pt_a[[1]]$pstate3+pt_a[[1]]$pstate4+pt_a[[1]]$pstate5) # For all pathologies

se=c(pt_a[[1]]$se1,
     pt_a[[1]]$se2,
     pt_a[[1]]$se3,
     pt_a[[1]]$se4,
     pt_a[[1]]$se5,
     pt_a[[1]]$se6,
     pt_a[[1]]$se7,
     rep(0,nrow(pt_a[[1]])))

borneup=prob*exp(-coef_conf*se/prob)
bornedown=prob*exp(coef_conf*se/prob)
trans=c(rep(0,nrow(pt_a[[1]])),
        rep(1,nrow(pt_a[[1]])),
        rep(2,nrow(pt_a[[1]])),
        rep(3,nrow(pt_a[[1]])),
        rep(4,nrow(pt_a[[1]])),
        rep(5,nrow(pt_a[[1]])),
        rep(6,nrow(pt_a[[1]])),
        rep(7,nrow(pt_a[[1]]))
)

trans=as.factor(trans)
trans=factor(trans,labels=c("Survie","Pathologie 1","Pathologie 2","Pathologie 3","Pathologie 4","Deces","Resiliation","Dependance"))

pt_a_graph <-  data.frame(time,prob,se,borneup,bornedown,trans)
# Plot on the age interval 65-90
pt_a_graph <- subset(pt_a_graph,pt_a_graph$time>=AgeMin & pt_a_graph$time<=AgeMax)
#Traçage du graphique
plot_pt_a <- ggplot(pt_a_graph,aes(x=time,y=prob),colour="black")+
  geom_line()+
  geom_ribbon(aes(x=time,ymin=bornedown,ymax=borneup),colour="darkgrey",alpha=0.2,linetype=0)
plot_pt_a <- plot_pt_a+scale_x_continuous("Age")+
  scale_y_continuous("Probabilité")+
  facet_wrap(~trans,ncol=2,scales="free_y")

# Print
png(file="figure/proba.png", width  = 8000, height = 4800,res=750,antialias = "cleartype" )
plot_pt_a
dev.off()

rm(time);rm(prob);rm(se);rm(borneup);rm(bornedown);rm(trans)
#----------------------------------------------------------
# plot (crude) transition probabilities and incidence rates with confidence intervals

# The initial data corresponds to the transitions probabilities between t-1 and t (not t and t+1). The layout is updated to have usual definition of a incidence rate 

time <- rep(qx_a[[1]]$time,6)
prob <- c(qx_a[[1]]$pstate2,
       qx_a[[1]]$pstate3,
       qx_a[[1]]$pstate4,
       qx_a[[1]]$pstate5,
       qx_a[[1]]$pstate6,
       qx_a[[1]]$pstate7)

se <- c(qx_a[[1]]$se2,
     qx_a[[1]]$se3,
     qx_a[[1]]$se4,
     qx_a[[1]]$se5,
     qx_a[[1]]$se6,
     qx_a[[1]]$se7)

borneup <- prob*exp(-coef_conf*se/prob)
bornedown <- prob*exp(coef_conf*se/prob)
trans <- c(rep(1,nrow(qx_a[[1]])),
        rep(2,nrow(qx_a[[1]])),
        rep(3,nrow(qx_a[[1]])),
        rep(4,nrow(qx_a[[1]])),
        rep(5,nrow(qx_a[[1]])),
        rep(6,nrow(qx_a[[1]])))

trans <- as.factor(trans)
trans <- factor(trans,labels=c("Pathologie neurologique","Pathologie diverse","Cancer","Demence","Deces","Resiliation"))
qx_a_brut <- data.frame(time,prob,se,borneup,bornedown,trans)

#Age interval: 65-90

qx_a_brut$time <- qx_a_brut$time-1
qx_a_brut <- qx_a_brut[qx_a_brut$time!=AgeMin-1,]
qx_a_brut <- subset(qx_a_brut,qx_a_brut$time<=AgeMax)
qx_a_brut$time <- rep(AgeMin:AgeMax,length(levels(trans)))

# Print
plot_qx_a <- ggplot(qx_a_brut,aes(x=time,y=prob),colour="black")+
  geom_line()+
  geom_ribbon(aes(x=time,ymin=bornedown,ymax=borneup),colour="darkgrey",alpha=0.2,linetype=0)
plot_qx_a <- plot_qx_a+scale_x_continuous("Age")+
  scale_y_continuous("Taux d'incidence")+
  facet_wrap(~trans,ncol=2,scales="free_y")

png(file="figure/incidence_rates.png", width  = 10000, height = 7500,res=1200,antialias = "cleartype" )
plot_qx_a
dev.off()
rm(time);rm(prob);rm(se);rm(borneup);rm(bornedown);rm(trans)

# Save results
save(qx_a_brut,file = paste(dossier,"/data/qx_a_brut.RData",sep=""))

#----------------------------------------------------------
# Incidence rates using KM estimator
surv_a <- survfit(Surv(Tstart,Tstop,status)~trans,data=t,type = "kaplan-meier")
qx_a_brut_marg <- getqxmarg(surv_a,AgeMin,AgeMax,strata=1:6,pas=1)
qx_a_brut_marg$trans  <- as.factor(qx_a_brut_marg$trans)
qx_a_brut_marg$trans <- factor(qx_a_brut_marg$trans,labels=c("Pathologie 1","Pathologie 2","Pathologie 3","Pathologie 4","Deces","Resiliation"))

# Save
save(qx_a_brut_marg,file = paste(dossier,"//qx_a_brut_marg.RData",sep=""))

# Merge all pathologies (i.e. 3 states = Décès, Résilié, Dépendant)

t2=t
levels(t2$to)= c("Décès","Dependant","Dependant","Dependant","Dependant","Résilié")
t2$trans[t2$trans==2]=1
t2$trans[t2$trans==3]=1
t2$trans[t2$trans==4]=1
t2$trans[t2$trans==5]=2
t2$trans[t2$trans==6]=3
t2$duplic=duplicated(t2)
t2=subset(t2,t2$duplic==FALSE)
t2$duplic=duplicated(t2[,1:7])
t2$status[t2$id %in% t2$id[which(t2$duplic)] & t2$trans==1]=1
t2$duplic=duplicated(t2)
t2=subset(t2,t2$duplic==FALSE)
t2$duplic=duplicated(t2)
t2=subset(t2,t2$duplic==FALSE)

# Estimation
surv_a2=survfit(Surv(Tstart,Tstop,status)~trans,data=t2,type = "kaplan-meier")
qx_a_brut_marg2=getqxmarg(surv_a2,AgeMin,AgeMax,strata=1:3,pas=1)
qx_a_brut_marg2$trans=as.factor(qx_a_brut_marg2$trans)
qx_a_brut_marg2$trans=factor(qx_a_brut_marg2$trans,labels=c("Dependance","Deces","Resiliation"))
save(qx_a_brut_marg2,file = paste(dossier,"/02 - Data/4.qx_a_brut_marg2.RData",sep=""))

