original/R commands.R
In HealthEconomicsHackathon/ceatutorial: What the Package Does (One Line, Title Case)
##########################################################################################
##========================================================================================
## 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(caToolsq)
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=<a3>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()
##########################################################################################
HealthEconomicsHackathon/ceatutorial documentation built on Nov. 7, 2019, 12:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##########################################################################################
##========================================================================================
## 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(caToolsq)
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=<a3>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()
##########################################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.