##########################################################################################
##========================================================================================
## R PREPARATION
##========================================================================================
##########################################################################################

### SET UP A WORKING DIRECTORY (NOT NECESSARY BUT RECOMMENDED)

setwd("C:/CEA with R tutorial")
setwd("~/C/CEA with R tutorial")

### THE FOLLOWING PACKAGES NEED TO BE INSTALLED INTO R IF THEY AREN'T ALREADY:
### mstate eha caTools mc2d msm cmprsk flexsurv BBmisc

install.packages("mstate")
install.packages("eha")
install.packages("caTools")
install.packages("mc2d")
install.packages("msm")
install.packages("cmprsk")
install.packages("flexsurv")
install.packages("BBmisc")

### READ IN REQUIRED PACKAGES INTO R

library(mstate)
library(eha)
library(caTools)
library(mc2d)
library(msm)
library(cmprsk)
library(flexsurv)
library(BBmisc)

### READ FUNCTIONS INTO R

source("function modelparam.R")
source("function Markov.R")
source("function semiMarkov.R")
source("function visualMarkov.R")
source("function visualsemiMarkov.R")
#source("function proportion.R") not run - requires actual data
source("function meanLY.R")
source("function PSAprob.R")
source("function PSAmeanLY.R")
source("function PSAQALY.R")
source("function CEplane.R")
source("function CEAC.R")



### READ IN DATA
#==============================================================================
# NOTE - THE DATASET USED IN THIS ILLUSTRATION IS NOT ACTUAL TRIAL DATA, BUT 
# HAS BEEN GENERATED BY DIGITISING PUBLISHED KAPLAN-MEIER CURVES. THEREFORE 
# PATIENT IDENTIFIERS (PATID) ARE NOT AVAILABLE. ALSO VARIABLES ON THE SAME ROW,
# SUCH AS TIME AND STATUS FOR DIFFERENT EVENTS, DO NOT NECESSARILY RELATE TO 
# THE SAME PATIENT. IT IS NOT POSSIBLE TO ACCURATELY REFLECT A
# "ONE PATIENT - ONE ROW" DATASET FOR THE MODEL IN THIS ILLUSTRATION WITHOUT THE
# ACTUAL TRIAL DATA.

# THE  VARIABLE HAS BEEN ADDED FROM THE REAL DATA TO ALLOW STATE-ARRIVAL
# EXTENDED MODELLING
#==============================================================================


data<-read.table("data.txt",sep="\t",skip=1) 
names(data)=c("patid","treat","prog","death", "prog_t","death_t",
              "progdeath", "progdeath_t")

### variables in the dataset are:
### patid          patient identification number
### treat          treatment    0=FC, 1=RFC
### prog           progression  0=no, 1=yes
### death          death        0=no, 1=yes
### prog_t         time (in months) to progression or last known not to have a progression
### death_t        time (in months) to death or last known to be alive
### progdeath      progression or death 0=no, 1=yes
### progdeath_t    time (in months) to progression or death or last known to be alive  
###                without progression

### convert times from months to years

data$prog_ty<-data$prog_t/12
data$death_ty<-data$death_t/12
data$progdeath_ty<-data$progdeath_t/12

### create data subsets based on treatment arm
RFCdata<-subset(data, treat==1)
FCdata<-subset(data, treat==0)

### create competing risk indicatior
### 0=did not die or have a progression, 1=had a progression, 2=died without progression
RFCdata$crstatus<-0
RFCdata$crstatus[RFCdata$prog==1]<-1
RFCdata$crstatus[RFCdata$prog==0 & RFCdata$death==1]<-2

FCdata$crstatus<-0
FCdata$crstatus[FCdata$prog==1]<-1
FCdata$crstatus[FCdata$prog==0 & FCdata$death==1]<-2

### create data subset for those who had progression
RFCdataprog<-subset(RFCdata, prog==1)
RFCdataprog$time<-RFCdataprog$death_ty-RFCdataprog$prog_ty


#######  ILLNESS-DEATH MODEL DATA PREPARATION WITH MSTATE

### create transition matrix
### sets up a matrix that shows the states that can be reached from each state. In this
### example the states progression and death can be reached from the progression-free 
### state, the death state can be reached from the progression state and no further 
### movement is possible from the death state
tmat<- transMat(list(c(2,3), 3, c()), 
                names = c("progression-free", "progression", "death"))
tmat
### create transition matrix with death split into two states
tmat2 <- transMat(x = list(c(2, 4), c(3), c(), c()), 
                  names=c("PFS", "prog", 
                          "death after prog", "death without prog"))
tmat2

### create dataset that can be used for multi-state modelling
covs<-c("treat","prog_ty")
msmcancer<-msprep(time=c(NA, "prog_ty", "death_ty"),
                  status=c(NA, "prog", "death"), data=data, trans=tmat,keep=covs)

### create datasets specific to each transition
msmcancer1<-subset(msmcancer,trans==1)
msmcancer2<-subset(msmcancer,trans==2)
msmcancer3<-subset(msmcancer,trans==3)

### create data subsets based on treatment arm

msmcancer1FC<-subset(msmcancer1,treat==0)
msmcancer2FC<-subset(msmcancer2,treat==0)
msmcancer3FC<-subset(msmcancer3,treat==0)
msmcancer1RFC<-subset(msmcancer1,treat==1)
msmcancer2RFC<-subset(msmcancer2,treat==1)
msmcancer3RFC<-subset(msmcancer3,treat==1)

##########################################################################################
##========================================================================================
## CONSIDERATION OF THE APPROPRIATENESS OF THE PROPORTIONAL HAZARDS ASSUMPTION
##========================================================================================
##########################################################################################

###================================================================
### LOG CUMULATIVE HAZARD VS LOG TIME PLOT progression-free -> progression
###================================================================

pdf('log-log plots.pdf', width=12, height=12)
par(mfrow=c(2,2))
plot(log(survfit(Surv(time,status)~ 1, data=msmcancer1RFC)$time), 
     log(-log(survfit(Surv(time,status)~ 1, data=msmcancer1RFC)$surv)), type="l",
     xlab="log(time)", ylab="log(Cumulative hazard)", xaxs="i",yaxs="i",
     las=1,cex.lab=1.5, cex.axis=1.5)
lines(log(survfit(Surv(time,status)~ 1, data=msmcancer1FC)$time), 
      log(-log(survfit(Surv(time,status)~ 1, data=msmcancer1FC)$surv)), lty=2)
legend(-4,-2, c("FC", "RFC"), bty="n", lty=c(2,1), cex=1.5)
title(main="(a) progression-free -> progression", cex.main=1.5)


###================================================================
### LOG CUMULATIVE HAZARD PLOT VS LOG TIME progression-free -> death
###================================================================

plot(log(survfit(Surv(time,status)~ 1, data=msmcancer2RFC)$time), 
     log(-log(survfit(Surv(time,status)~ 1, data=msmcancer2RFC)$surv)), type="l",
     xlab="log(time)", ylab="log(Cumulative hazard)", xaxs="i",yaxs="i", ylim=c(-6,-2), 
     las=1,cex.lab=1.5, cex.axis=1.5)
lines(log(survfit(Surv(time,status)~ 1, data=msmcancer2FC)$time), 
      log(-log(survfit(Surv(time,status)~ 1, data=msmcancer2FC)$surv)), lty=2)
legend(-4,-3, c("FC", "RFC"), bty="n", lty=c(2,1), cex=1.5)
title(main="(b) progression-free -> death", cex.main=1.5)

###================================================================
### LOG CUMULATIVE HAZARD PLOT VS LOG TIME progression -> death
###================================================================

plot(log(survfit(Surv(time,status)~ 1, data=msmcancer3RFC)$time), 
     log(-log(survfit(Surv(time,status)~ 1, data=msmcancer3RFC)$surv)), type="l",
     xlab="log(time)", ylab="log(Cumulative hazard)", xaxs="i",yaxs="i", 
     las=1,cex.lab=1.5, cex.axis=1.5)
lines(log(survfit(Surv(time,status)~ 1, data=msmcancer3FC)$time), 
      log(-log(survfit(Surv(time,status)~ 1, data=msmcancer3FC)$surv)), lty=2)
legend(-4,-2, c("FC", "RFC"), bty="n", lty=c(2,1), cex=1.5)
title(main="(c) progression -> death", cex.main=1.5)

dev.off()

###================================================================
### CUMULATIVE HAZARD PLOT:  I.E. -LOG(KM SURVIVAL) progression-free -> progression
###================================================================

pdf('cum haz plots.pdf', width=12, height=12)
par(mfrow=c(2,2))
par(mar= c(5, 7, 4, 2)+0.1)
plot(survfit(Surv(time,status)~ 1, data=msmcancer1RFC)$time, 
     -log(survfit(Surv(time,status)~ 1, data=msmcancer1RFC)$surv), type="l",
     xlab="time", ylab="Cumulative hazard", xaxs="i",yaxs="i", 
     las=1,cex.lab=1.5, cex.axis=1.5, col=1)
lines(survfit(Surv(time,status)~ 1, data=msmcancer1FC)$time, 
      -log(survfit(Surv(time,status)~ 1, data=msmcancer1FC)$surv), lty=2)
legend(1, 0.6, c("FC", "RFC"), bty="n", lty=c(2,1), cex=1.5)
title(main="(a) progression-free -> progression", cex.main=1.5)

###================================================================
### CUMULATIVE HAZARD PLOT:  I.E. -LOG(KM SURVIVAL) progression-free -> death
###================================================================

plot(survfit(Surv(time,status)~ 1, data=msmcancer2RFC)$time, 
     -log(survfit(Surv(time,status)~ 1, data=msmcancer2RFC)$surv), type="l",
     xlab="time", ylab="", xaxs="i",yaxs="i", 
     las=1,cex.lab=1.5, cex.axis=1.5, col=1)
lines(survfit(Surv(time,status)~ 1, data=msmcancer2FC)$time, 
      -log(survfit(Surv(time,status)~ 1, data=msmcancer2FC)$surv), lty=2)
title(main="(b) progression-free -> death", cex.main=1.5)
legend(2,0.02, c("FC", "RFC"), bty="n", lty=c(2,1), cex=1.5)
title(ylab="Cumulative hazard",cex.lab = 1.5,line = 4.5)

###================================================================
### CUMULATIVE HAZARD PLOT:  I.E. -LOG(KM SURVIVAL) progression-> death
###================================================================

plot(survfit(Surv(time,status)~ 1, data=msmcancer3RFC)$time, 
     -log(survfit(Surv(time,status)~ 1, data=msmcancer3RFC)$surv), type="l",
     xlab="time", ylab="Cumulative hazard", xaxs="i",yaxs="i", 
     las=1,cex.lab=1.5, cex.axis=1.5, col=1)
lines(survfit(Surv(time,status)~ 1, data=msmcancer3FC)$time, 
      -log(survfit(Surv(time,status)~ 1, data=msmcancer3FC)$surv), lty=2)
legend(2,0.2, c("FC", "RFC"), bty="n", lty=c(2,1), cex=1.5)
title(main="(c) progression -> death", cex.main=1.5)

dev.off()

##########################################################################################
##========================================================================================
## DECIDING WHETHER THE MARKOV PROPERTY IS REASONABLE
##========================================================================================
##########################################################################################

### build a Cox Markov model including time in previous state as a covariate

# code shown for illustration only. Markov modelling not appropriate with digitised data.
# results in tutorial paper were based on actual trial data
CoxMarkov<-coxph(Surv(Tstart,Tstop,status)~treat+prog_ty,
                 data=msmcancer3,method = "breslow") 
summary(CoxMarkov)

##########################################################################################
##========================================================================================
## BUILDING PARAMETRIC MULTI-STATE MODELS:
## DISPLAYING THE MODEL PARAMETERS
##========================================================================================
##########################################################################################

### semi-Markov Weibull model for transition 3 with treatment as only covariate 
modelparam()
### Markov Weibull model for transition 3 with treatment as only covariate 
modelparam(Markov=TRUE) # code shown for illustration only. Markov modelling not 
                        # appropriate with digitised data
### expontential model for transition 3 with treatment as only covariate 
modelparam(dist="exp")
### semi-Markov gompertz  model for transition 3 with treatment as only covariate 
modelparam(dist="gom")
### semi-Markov loglogistic  model for transition 3 with treatment as only covariate 
modelparam(dist="logl")
### semi-Markov lognormal  model for transition 3 with treatment as only covariate
modelparam(dist="logn")
### semi-Markov generalised gamma  model for transition 3 with treatment as only covariate
modelparam(dist="gam")
### semi-Markov Weibull model for transition 2 with treatment as only covariate 
modelparam(transnum=2)


##########################################################################################
##========================================================================================
## ASSESSMENT OF THE FIT FOR PROGRESSION-> DEATH
##========================================================================================
##########################################################################################

### THE FIT FOR PROGRESSION->DEATH CAN BE ASSESSED VISUALLY AND USING AICs BY RUNNING THE 
### FOLLOWING LINE OF CODE. THIS WILL RUN THE INDICATED R PROGRAM.

source("assessment of fit for prog to death.R", echo=TRUE)
### or alternatively open the program above and run it manually

##########################################################################################
##========================================================================================
## CALCULATING STATE OCCUPANCY PROBABILITIES
##========================================================================================
##########################################################################################

### RUNNING THE FOLLOWING LINES OF CODE WILL CREATE THE STATE OCCUPANCY PROBABILITIES.
### MODELS FOR ALL COMBINATIONS OF USING EACH OF THE SIX DISTRIBUTIONS FOR EACH 
### OF THE TRANSITIONS ARE CONSIDERED.

# ***** WARNING : THIS FILE BUILDS 432 MODELS AND TAKES MANY HOURS AND HENCE IS  
# ***** BEST RUN OVERNIGHT. THE RESULTANT DATA FILE ("allcombos.RData) IS INCLUDED ON THE 
# ***** WEB PAGE FOR CONVENIENCE. THEREFORE THE CODE DOES NOT NEED TO BE RUN. 


source("building all combinations of msms.R",echo=TRUE)
#source("building all combinations of msms PPS.R") # not run

### or alternatively open the programs above and run them manually

save.image("allcombos.RData") 
load("allcombos.RData") 

##########################################################################################
##========================================================================================
## VISUAL ASSESSMENT OF FITS FOR PROGRESSION-FREE TO PROGRESSION 
##========================================================================================
##########################################################################################

### RUNNING THE FOLLOWING LINE OF CODE WILL CREATE PLOTS ALLOWING ASSESSMENT OF THE 
### DISTRIBUTION TO USE FOR PROGRESSION-FREE -> PROGRESSION (PLOTS SHOWN IN APPENDIX 1)

source("visual assessment of fit for progression-free to progression.R", echo=TRUE)
### or alternatively open the program above and run it manually

##########################################################################################
##========================================================================================
## VISUAL ASSESSMENT OF FITS FOR PROGRESSION-FREE TO DEATH 
##========================================================================================
##########################################################################################

### RUNNING THE FOLLOWING LINE OF CODE WILL CREATE PLOTS ALLOWING ASSESSMENT OF THE 
### DISTRIBUTION TO USE FOR PROGRESSION-FREE -> DEATH (PLOTS SHOWN IN APPENDIX 1)

source("visual assessment of fit for progression-free to death.R", echo=TRUE)
### or alternatively open the program above and run it manually

##########################################################################################
##========================================================================================
## MEASURING CLINICAL EFFECTIVENESS -BASE CASE
##========================================================================================
##########################################################################################

#############################################################
####### building the BASE CASE  multi-state model
#############################################################

### ALL THE COMPUTATIONAL TIMES INCLUDED AS COMMENTS IN THIS R SCRIPT AND THE OTHERS ON 
### THE WEBSITE RELATE TO A PC RUNNING WINDOWS 7 64-BIT, WITH 4 GB OF RAM, AN INTEL
### CORE 2 QUAD CPU 2.83 GHz PROCESSOR.

#### THIS SECTION OF CODE DOES NOT NEED TO BE RUN IF allcombos.RData (SEE ABOVE) HAS
#### BEEN LOADED INTO R


### RFC 

ptm <- proc.time() #start counting the time it takes
gom1gam2gom3simRFC<-semiMarkov(ncovs=c(1,1,1),covs=rbind("treat", "treat", "treat"),
                                                 coveval=rbind(1,1,1),  
                                                 dist=cbind("gom", "gam", "gom"),
                                                 timeseq=seq(0,4,1/12), 
                                                 timeseq_ext=c(seq(49/12,100/12,1/12), seq(100/12+1/144, 12, 1/144), 
                                                               seq(12+1/600,15, 1/600)),
                                                 trans=tmat2, M=5000)
(proc.time() - ptm) #stop counting the time it took (seconds)
(proc.time() - ptm)/60 #stop counting the time it took (minutes)
(proc.time() - ptm)/3600 #stop counting the time it took (hours)
#### TOOK 2.2 MINS

### FC 

ptm <- proc.time() #start counting the time it takes
gom1gam2gom3simFC<-semiMarkov(ncovs=c(1,1,1),covs=rbind("treat", "treat", "treat"),
                                                coveval=rbind(0,0,0),  
                                                dist=cbind("gom", "gam", "gom"),
                                                timeseq=seq(0,4,1/12), 
                                                timeseq_ext=c(seq(49/12,100/12,1/12),
                                                              seq(100/12+1/144, 12, 1/144), 
                                                              seq(12+1/600,15, 1/600)),
                                                trans=tmat2, M=5000)
(proc.time() - ptm) #stop counting the time it took (seconds)
(proc.time() - ptm)/60 #stop counting the time it took (minutes)
(proc.time() - ptm)/3600 #stop counting the time it took (hours)
#### TOOK 2.2 MINS

#### SAVE DATA (IF DESIRED)

save.image("basecase15gom1gam2gom3.RData")


#############################################################
####### calculating base case mean Life Years
#############################################################

load("basecase15gom1gam2gom3.RData")
######## remove times that don't equate to months before calculating area under curves 

is.wholenumber <-
  function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

gom1gam2gom3simRFCmon<-gom1gam2gom3simRFC
gom1gam2gom3simRFCmon<-subset(gom1gam2gom3simRFCmon,is.wholenumber(gom1gam2gom3simRFCmon[,1]*12)==1)
gom1gam2gom3simFCmon<-gom1gam2gom3simFC
gom1gam2gom3simFCmon<-subset(gom1gam2gom3simFCmon,is.wholenumber(gom1gam2gom3simFCmon[,1]*12)==1)

### RFC treatment arm 
### mean life years in initial progression-free state 
(RFCmeanLYPFSdis<-meanLY(object=gom1gam2gom3simRFCmon, state=1,instate=TRUE,discounted=TRUE, rate=0.035))

### mean life years in progression
(RFCmeanLYprogdis<-meanLY(object=gom1gam2gom3simRFCmon, state=2,instate=TRUE,discounted=TRUE, rate=0.035))


### FC treatment arm 
### mean life years in initial progression-free state 
(FCmeanLYPFSdis<-meanLY(object=gom1gam2gom3simFCmon,state=1,instate=TRUE,discounted=TRUE, rate=0.035))

### mean life years in progression
(FCmeanLYprogdis<-meanLY(object=gom1gam2gom3simFCmon,state=2,instate=TRUE,discounted=TRUE, rate=0.035))

################################################################################
##========================================================================================
## INCORPORATING COSTS AND MEASURING COST-EFFECTIVENESS -BASE CASE
##========================================================================================
##########################################################################################

### costs

cost_support_PFS_RFC<- 28*12*RFCmeanLYPFSdis
cost_support_PFS_FC<- 28*12*FCmeanLYPFSdis

cost_support_prog_RFC<- 84*12*RFCmeanLYprogdis
cost_support_prog_FC<- 84*12*FCmeanLYprogdis

cost_therapy_RFC<-(5179/((RFCmeanLYprogdis + FCmeanLYprogdis)/2*12))*12*RFCmeanLYprogdis
cost_therapy_FC<-(5179/((RFCmeanLYprogdis + FCmeanLYprogdis)/2*12))*12*FCmeanLYprogdis

mean_cost_PFS_RFC<-10113+1224+2776+1109+21+1109+cost_support_PFS_RFC+592+640
mean_cost_PFS_FC<- 2790+1115+22+1115+cost_support_PFS_FC+360+507

mean_cost_prog_RFC<-cost_support_prog_RFC+cost_therapy_RFC
mean_cost_prog_FC<-cost_support_prog_FC+cost_therapy_FC

(mean_cost_RFC<-mean_cost_PFS_RFC+mean_cost_prog_RFC)
(mean_cost_FC<-mean_cost_PFS_FC+mean_cost_prog_FC)


##########################################################################################
##========================================================================================
## SELECTED ONE-WAY SENSITIVITY ANALYSES
##========================================================================================
##########################################################################################

### RUNNING THE FOLLOWING LINES OF CODE WILL CALCULATE THE COST PER QALYS GAINED FOR EACH
### OF THE ONE-WAY SENSITIVITY ANALYSES SHOWN IN THE TUTORIAL ARTICLE. IN PARTICULAR, 
### ASSUMING NO TREATMENT EFFECT IN THE UNOBSERVED PERIOD IS EVALUATED.

source("one way sens.R", echo=TRUE)
### or alternatively open the program above and run it manually

##########################################################################################
##========================================================================================
## ONE-WAY SENSITIVITY ANALYSIS: ALTERNATIVE FITS FOR EACH OF THE TRANSITIONS
##========================================================================================
##########################################################################################

### RUNNING THE FOLLOWING LINE OF CODE WILL CALCULATE THE COST PER QALYS GAINED FOR 432
### OF THE MODELS CREATED IN THE 'CALCULATING STATE OCCUPANCY PROBABILITIES' STEP ABOVE

source("sensitivity analysis - alternative fits ICERs table.R", echo=TRUE)
### or alternatively open the program above and run it manually

##########################################################################################
##========================================================================================
## PSA: UNCERTAINTY IN THE SURVIVAL REGRESSION MODEL PARAMETERS
##========================================================================================
##########################################################################################

# ***** WARNING : THIS STEP CAN BE TIME-CONSUMING THEREFORE THE RESULTANT DATA FILE 
# ("PSAdata15gom1gam2gom3.RData") IS INCLUDED ON THE WEB PAGE FOR CONVENIENCE. THEREFORE
# THE CODE IN THIS SECTION DOES NOT NEED TO BE RUN. 


source("function PSAprob.R")

##### RFC 
ptm <- proc.time() #start counting the time it takes
PSASMSARFCwitherrid<-PSAprob(ntrans=3, ncovs=c(1,1,1), 
                                     covs=rbind("treat", "treat", "treat"),
                                     coveval=rbind(1,1,1), 
                                     dist=cbind("gom", "gam", "gom"),
                                     timeseq=seq(0,4,1/12), 
                                     timeseq_ext=c(seq(49/12,97/12,1/12), seq(97/12+1/144, 12, 1/144), 
                                     seq(12+1/600,14, 1/600),  seq(14+1/1200,15, 1/1200)),
                                     data=msmcancer, trans=tmat, Markov=FALSE,
                                     nruns=1000,seedno=12345, M=100)
(proc.time() - ptm) #stop counting the time it took (seconds)
(proc.time() - ptm)/60 #stop counting the time it took (minutes)
(proc.time() - ptm)/3600 #stop counting the time it took (hours)
### 1000 runs took 20 mins

##### FC 
ptm <- proc.time() #start counting the time it takes
PSASMSAFCwitherrid<-PSAprob(ntrans=3, ncovs=c(1,1,1), 
                                    covs=rbind("treat", "treat", "treat"),
                                    coveval=rbind(0,0,0), 
                                    dist=cbind("gom", "gam", "gom"),
                                    timeseq=seq(0,4,1/12),  
                                    timeseq_ext=c(seq(49/12,97/12,1/12), seq(97/12+1/144, 12, 1/144), 
                                                  seq(12+1/600,14, 1/600),  seq(14+1/1200,15, 1/1200)),
                                    data=msmcancer, trans=tmat, Markov=FALSE,
                                    nruns=1000,seedno=12345, M=100)
(proc.time() - ptm) #stop counting the time it took (seconds)
(proc.time() - ptm)/60 #stop counting the time it took (minutes)
(proc.time() - ptm)/3600 #stop counting the time it took (hours)
### 1000 runs  took  20 mins


##### SAVE DATA - RECOMMENDED
save.image("PSAdata15gom1gam2gom3.RData") 


load("PSAdata15gom1gam2gom3.RData")

##########################################################################################
##========================================================================================
## PSA: REMOVING THE DRAWS WITH COMPUTATIONAL DIFFICULTIES 
##========================================================================================
##########################################################################################

### identify the error ids for either RFC or FC

errid_<-sort(unique(c(PSASMSARFCwitherrid[[1]],PSASMSAFCwitherrid[[1]])))
length(errid_) #102 : no. of errors
temp<-seq(1,1000,1)
temp2<-temp[-errid_] # position of the 898 usuable draws

temp<-seq(1,1000,1)
temp<-temp[-PSASMSARFCwitherrid[[1]]] # currently 961 draws for rfc 

##### remove the appropriate draws from rfc  to obtain the 898 usuable draws 
temp3<-rep(NA,length(errid_) )
for (i in 1:length(errid_))
  if(length(which(temp==errid_[i]))>0 ) 
    temp3[i]<-which(temp==errid_[i])

temp3<- temp3[!is.na(temp3)]
PSASMSARFC<-PSASMSARFCwitherrid[[2]][-temp3]

##### obtain the usuable draws for fc 
temp<-seq(1,1000,1)
temp<-temp[-PSASMSAFCwitherrid[[1]]] 

temp3<-rep(NA,length(errid_) )
for (i in 1:length(errid_))
  if(length(which(temp==errid_[i]))>0 ) 
    temp3[i]<-which(temp==errid_[i])

temp3<- temp3[!is.na(temp3)]
PSASMSAFC<-PSASMSAFCwitherrid[[2]][-temp3]

##########################################################################################
##========================================================================================
## PSA: MEAN LIFE YEARS 
##========================================================================================
##########################################################################################

### RFC treatment arm 
### mean life years in initial progression-free state 
(PSARFCmeanLYPFSdis<-PSAmeanLY(object=PSASMSARFC, state=1,instate=TRUE,discounted=TRUE, rate=0.035))

### mean life years in progression
(PSARFCmeanLYprogdis<-PSAmeanLY(object=PSASMSARFC,state=2,instate=TRUE,discounted=TRUE, rate=0.035))

### FC treatment arm 
### mean life years in initial progression-free state 
(PSAFCmeanLYPFSdis<-PSAmeanLY(object=PSASMSAFC, state=1,instate=TRUE,discounted=TRUE, rate=0.035))

### mean life years in progression
(PSAFCmeanLYprogdis<-PSAmeanLY(object=PSASMSAFC, state=2,instate=TRUE,discounted=TRUE, rate=0.035))



##########################################################################################
##========================================================================================
## PSA: MEAN QALYs
##========================================================================================
##########################################################################################

set.seed(12345)
utility1=rbeta(1000,800,200)
set.seed(12345)
utility2=rbeta(1000,600,400)

### RFC treatment arm 
### mean QALYs in initial progression-free state 
PSARFCQALYPFSdis<-PSAQALY(object=PSASMSARFC, utility=utility1[1:898], state=1,discounted=TRUE, dis1yronwards=TRUE, rate=0.035)
mean(PSARFCQALYPFSdis)

### mean QALYs in progression
PSARFCQALYprogdis<-PSAQALY(object=PSASMSARFC, utility=utility2[1:898],state=2,discounted=TRUE, rate=0.035)
mean(PSARFCQALYprogdis)

### FC treatment arm 
### mean QALYs in initial progression-free state 
PSAFCQALYPFSdis<-PSAQALY(object=PSASMSAFC, utility=utility1[1:898],state=1,discounted=TRUE, rate=0.035)
mean(PSAFCQALYPFSdis)

### mean QALYs in progression
PSAFCQALYprogdis<-PSAQALY(object=PSASMSAFC, utility=utility2[1:898],state=2,discounted=TRUE, rate=0.035)
mean(PSAFCQALYprogdis)

  
##########################################################################################
##========================================================================================
## PSA: INCORPORATING COSTS AND MEASURING COST-EFFECTIVENESS 
##========================================================================================
##########################################################################################

### costs probabilistic
 
nruns<-1000     
set.seed(12345)   
cost_PFS_support<- rpert(nruns,14,28,42)
     
temp1<-84+5179/((PSARFCmeanLYprogdis+PSAFCmeanLYprogdis)/2*12)
cost_prog_support_dis_<- rpert(nruns,0.5*temp1,temp1,1.5*temp1)

cost_admin_oral <- rpert(nruns, 174, 280, 482)
cost_admin_complex <- rpert(nruns, 210, 430, 795)
cost_bmt<-rpert(nruns, 34318.25, 47565.05, 54646.47)
cost_transfusion<-rpert(nruns, 173.84, 289.73, 405.62)
cost_blood_unit<-rpert(nruns, 96.67, 161.11, 225.26)
     
cost_prog_RFC_dis<-cost_prog_support_dis_*12*PSARFCmeanLYprogdis
cost_prog_FC_dis<-cost_prog_support_dis_*12*PSAFCmeanLYprogdis

cost_PFS_RFC_dis<-10113+1222+2776+1109+21+1109+
  cost_PFS_support*12*PSARFCmeanLYPFSdis+
  cost_bmt*5/402+(cost_transfusion*318+cost_blood_unit*1025)/402

cost_PFS_FC_dis<-2790+1115+22+1115+cost_PFS_support*12*PSAFCmeanLYPFSdis+
  cost_bmt*3/396+(cost_transfusion*269+cost_blood_unit*762)/396

total_cost_RFC_dis<-cost_prog_RFC_dis+cost_PFS_RFC_dis
total_cost_FC_dis<-cost_prog_FC_dis+cost_PFS_FC_dis

mean_cost_prog_RFC_dis<-mean(cost_prog_RFC_dis)
mean_cost_prog_FC_dis<-mean(cost_prog_FC_dis)

mean_cost_PFS_RFC_dis<-mean(cost_PFS_RFC_dis)
mean_cost_PFS_FC_dis<-mean(cost_PFS_FC_dis)

mean_total_cost_RFC_dis<-mean(total_cost_RFC_dis)
mean_total_cost_FC_dis<-mean(total_cost_FC_dis)

mean_cost_prog_RFC_dis
mean_cost_prog_FC_dis

mean_cost_PFS_RFC_dis
mean_cost_PFS_FC_dis

mean_total_cost_RFC_dis
mean_total_cost_FC_dis


##########################################################################################
##========================================================================================
## PSA: COST-EFFECTIVENESS PLANE 
##========================================================================================
##########################################################################################

### incremental discounted qaly     

incQALY<-PSARFCQALYPFSdis+PSARFCQALYprogdis - 
                 PSAFCQALYPFSdis-PSAFCQALYprogdis

### incremental costs
incCost<-total_cost_RFC_dis - total_cost_FC_dis    

incCost<-incCost[1:898]

min(incQALY)
min(incCost)

max(incQALY)
max(incCost)


### plot of cost-effectiveness plane 

pdf('CEplane.pdf', width=8, height=8)
CEplane(x=incQALY, y=incCost,xlower=-2, xupper=2,
        ylower=0, yupper=12000, ICER=NA)  
abline(0,30000)
text(0.75,7000, "WTP=£30,000")
dev.off()

##########################################################################################
##========================================================================================
## PSA: COST-EFFECTIVENESS ACCEPTABILITY CURVE
##========================================================================================
##########################################################################################

### total discounted qaly for each treatment group

QALY_RFC_dis<-PSARFCQALYPFSdis+PSARFCQALYprogdis
QALY_FC_dis<-PSAFCQALYPFSdis+PSAFCQALYprogdis

### cost-effectiveness acceptibility curves 

pdf('CEAC.pdf', width=8, height=8)
CEAC(QALY1=QALY_RFC_dis,QALY2= QALY_FC_dis,
     cost1=total_cost_RFC_dis[1:898], cost2=total_cost_FC_dis[1:898], secondcurve=FALSE,nruns=898) 
dev.off()

##########################################################################################     


     