## ------------------------------------------------------------------------ ##
##  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)
library(gss)

#----------------------------------------------------------
# Load data and results
load( file = paste(dossier, "/data/data_example.RData", sep=""))
load(file = paste(dossier,"/data/qx_a_brut.RData",sep=""))
load(file = paste(dossier,"/data/qx_a_brut_marg.RData",sep=""))
load(file = paste(dossier,"/data/qx_a_brut_marg2.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_trans2=c("Dependance","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) 

#----------------------------------------------------------
# Whittaker Henderson: Find the best parameters
crit_WH<- data.frame(alpha=rep(0,length(names_trans)),ordre=rep(0,length(names_trans)),trans=names_trans)

# Pathologie 1
y <- qx_a_brut$prob[qx_a_brut$trans==names_trans[1]]
# Best parameters W-H 
preselect_W_H1(y, poids)
compar_W_H1(AgeMin:AgeMax,y,poids,c(2,3,4),c(2,3,4,5,6))
# Select manually the best parameters
crit_WH$alpha[crit_WH$trans==names_trans[1]]=4
crit_WH$ordre[crit_WH$trans==names_trans[1]]=2

# Pathologie 2 ...
# IDEM for each transition

#----------------------------------------------------------
# Analyze the results graphically 
# A dataframe with the results
# !!! the 'crit_WH' list should be completed

qx_a_res=data.frame(Age=qx_a_brut$time,qx.brut=qx_a_brut$prob,trans=qx_a_brut$trans)
qx_a_res=qx_a_res[qx_a_res$trans %in% names_trans,]
qx_a_res$qx.fit=0
qx_a_res$infl=0
qx_a_res$resid=0
qx_a_res$resid.std=0
qx_a_res$resid.stud=0

# Optimal parameter and the compute the effect of each point
for (i in names_trans){
  y=qx_a_res$qx.brut[qx_a_res$trans==i]
  ordre=crit_WH$ordre[crit_WH$trans==i]
  alpha=crit_WH$alpha[crit_WH$trans==i]
  temp.WH=W_H1(y,poids,ordre,alpha)
  qx_a_res$qx.fit[qx_a_res$trans==i]=temp.WH$fit
  qx_a_res$infl[qx_a_res$trans==i]=temp.WH$infl
  temp.resid=resid(y,temp.WH$fit,temp.WH$infl)
  qx_a_res$resid[qx_a_res$trans==i]=temp.resid$resid
  qx_a_res$resid.std[qx_a_res$trans==i]=temp.resid$resid.std
  qx_a_res$resid.stud[qx_a_res$trans==i]=temp.resid$resid.stud
  rm(y,ordre,alpha,temp.WH,temp.resid)
}
test.qx.fit=test_process.fit(qx_a_res,names_trans,t_stat_a)

