## ------------------------------------------------------------------------ ##
##  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=""))

#----------------------------------------------------------
# Parameters
conf=95
coef_conf=qnorm((1-conf/100)/2)
AgeMin=65
AgeMax=90
DureeMax=60

names_trans=c("Pathologie 1","Pathologie 2","Pathologie 3","Pathologie 4","Deces","Resiliation")
names_trans_bis=c("Pathologie neurologique","Pathologie diverse","Cancer","Demence","Deces","Resiliation")

#----------------------------------------------------------
# Calculate weight
# As an example, we consider weights equal to 1. In the paper, they are calculed using the exposure at risk. 
poids <- rep(1, AgeMax - AgeMin + 1) 

#----------------------------------------------------------
#Bootstrap process
qx_res_boot=boot.process.qx(process.qx, data=t, B = 1000, id = "id", AgeMin=AgeMin, AgeMax=AgeMax, crit_WH=crit_WH, crit_WH_marg= crit_WH_marg, poids=poids) 

save(qx_res_boot,file = paste(dossier,"/Data/qx_res_boot.RData",sep=""))
#----------------------------------------------------------
# Calculate mean, variance ... 
stat_boot=data.frame(Age=qx_a_res$Age,trans=qx_a_res$trans,qx.fit=qx_a_res$qx.fit)
stat_boot$trans=factor(stat_boot$trans,labels=names_trans_bis)
n=length(qx_res_boot)
temp=0
temp2=0
temp3=matrix(0,n,dim(stat_boot)[1])
for (i in 1:n){
  temp=temp+qx_res_boot[[i]]$qx.fit
  temp2=temp2+(qx_res_boot[[i]]$qx.fit-stat_boot$qx.fit)^2
  temp3[i,]=qx_res_boot[[i]]$qx.fit
}
temp4=apply(temp3,2,quantile,probs=c(1-conf/100,conf/100))
stat_boot$mean=temp/n
stat_boot$var=temp2/(n-1)
stat_boot$disp=sqrt(stat_boot$var*(n-1)/n)/stat_boot$qx.fit
stat_boot$bornedown=temp4[1,]
stat_boot$borneup=temp4[2,]
rm(n,temp,temp2,temp3,temp4)

# Plot the results
plot_stat_boot=ggplot(stat_boot,aes(x=Age,y=qx.fit),colour="black")+
  geom_point(aes(x=Age,y=qx.fit,shape="Lois estimées"))+
  geom_point(aes(x=Age,y=mean,shape="Lois simulées"))+
  geom_line(aes(x=Age,y=qx.fit,linetype="Lois estimées"))+
  geom_line(aes(x=Age,y=mean,linetype="Lois simulées"))+
  geom_ribbon(aes(x=Age,ymin=bornedown,ymax=borneup),colour="darkgrey",alpha=0.2,linetype=0)
plot_stat_boot=plot_stat_boot+scale_x_continuous("Age")+
  scale_y_continuous("Taux d'incidence")+
  scale_shape_manual("", breaks = c("Lois estimées","Lois simulées"),values = c(1,32))+
  scale_linetype_manual("", breaks = c("Lois estimées","Lois simulées"),values = c(0,1))+
  theme(legend.position="bottom",
        legend.background = element_rect(colour = "black"))+
  theme(legend.key = element_rect(colour = "white"))+
  facet_wrap(~trans,ncol=2,scales="free_y")

