.Rhistory

t1 <- Sys.time() - p1  # calculate time to run the analusis by extracting p from the current system time
################ probabilistic analysis method 0.1 ###########################################
####### Generate random values inside a Markov function - use sapply  #################################
p2 <- Sys.time()          # save system time
set.seed(1)               # set a seed to be able to reproduce the same results
Markov <- function(i) {   # Markov model function - for different values of index i calculate the Markov model
p.HD <- rbeta(1, 20, 980)                 # probability from healthy to dead
p.HS <- rbeta(1, 50, 950)                 # probability from healthy to sick
p.SD <- rbeta(1, 100,900)                 # probability from sick to dead
c.H  <- rnorm(1, 400, 50)                 # cost of being in healthy state
c.S  <- rnorm(1, 100, 80)                 # cost of being in sick state
c.D  <- 0                                 # cost of being dead
u.H  <- rnorm(1, 0.8, 0.02)               # utility of being in healthy state
u.S  <- rnorm(1, 0.5, 0.02)               # utility of being in the sick state
u.D  <- 0                                 # utility of being dead
m.P <- matrix(c(1 - p.HD - p.HS, p.HS, p.HD, # transition matrix
0,           1 - p.SD, p.SD,
0,               0,       1),
byrow = TRUE, nrow = n.s,
dimnames = list(v.n, v.n))
m.TR <- matrix(NA, nrow = n.t + 1, ncol = n.s,
dimnames = list(0:n.t, v.n))   # create Markov trace
m.TR[1, ] <- c(1, 0, 0)                       # initialize Markov trace
for (t in 1:n.t) {                            # throughout the number of cycles
m.TR[t + 1, ] <- m.TR[t, ] %*% m.P          # estimate the Markov trace for cycle t + 1
} # close the for loop for time
v.c <- c(c.H, c.S, c.D)         # vector of costs
v.e <- c(u.H, u.S, u.D)         # vector of utilities
TC <- t(m.TR %*% v.c) %*% v.dw  # total cost for ith simulation
TE <- t(m.TR %*% v.e) %*% v.dw  # total QALYs for ith simulation
c(TC, TE)                       # output of the function TC and TE
} # close the Markov function
Results <- sapply(1:n.sim, Markov)   # Run the Markov function for n.sim times and save the results in the variable Results
TC2     <- Results[1, ] # save costs separately
TE2     <- Results[2, ] # save QALYs separately
t2 <- Sys.time() - p2  # calculate time to run the analysis by extracting p from the current system time
##############################################################################################
################ probabilistic analysis method 1.0 ###########################################
##############################################################################################
################ generate values outside of the for loop #####################################
p3 <- Sys.time()         # save system time
set.seed(1)              # set a seed to be able to reproduce the same results
p.HD <- rbeta(n.sim, 20,  980)                # probability from healthy to dead
p.HS <- rbeta(n.sim, 50,  950)                # probability from healthy to sick
p.SD <- rbeta(n.sim, 100, 900)                # probability from sick to dead
c.H  <- rnorm(n.sim, 400, 50)                 # cost of being in healthy state
c.S  <- rnorm(n.sim, 100, 80)                 # cost of being in sick state
c.D  <- 0
u.H  <- rnorm(n.sim, 0.8, 0.02)               # utility of being in healthy state
u.S  <- rnorm(n.sim, 0.5, 0.02)               # utility of being in the sick state
u.D  <- 0                                     # utility of being dead
TC3 <- TE3 <- vector("numeric", n.sim)        # initiate te vector to store costs and effects
m.TR <- array(0, dim = c(n.t + 1, n.s, n.sim)) # initiate the Markov trace
for (i in 1: n.sim){  # start the loop for every simulation
m.P <- matrix(c(1 - p.HD[i] - p.HS[i], p.HS[i], p.HD[i], # transition matrix
0,                 1 - p.SD[i], p.SD[i],
0,                           0,       1),
byrow = TRUE, nrow = n.s, ncol = n.s,
dimnames = list(v.n, v.n))
m.TR[1, , i] <- c(1, 0, 0)                    # initialize Markov trace
for (t in 1:n.t) {                            # throughout the number of cycles
m.TR[t + 1, , i] <- m.TR[t , , i] %*% m.P   # estimate the Markov trace for cycle t + 1
} # close the for loop for time
v.c <- c(c.H[i], c.S[i], c.D)    # vector of costs
v.e <- c(u.H[i], u.S[i], u.D)    # vector of utilities
TC3[i] <- t(m.TR[, , i] %*% v.c) %*% v.dw  # calculate total discounted cost for ith simulation
TE3[i] <- t(m.TR[, , i] %*% v.e) %*% v.dw  # calculate total discounted QALYs for ith simulation
} # close function
t3 <- Sys.time() - p3  # calculate time to run the analusis by extracting p from the current system time
################ probabilistic analysis method 1.1 ###########################################
################ generate values outside of the Markov function - use sapply #################
p4 <- Sys.time()         # save system time
set.seed(1)              # set a seed to be able to reproduce the same results
p.HD <- rbeta(n.sim, 20,  980)        # probability from healthy to dead
p.HS <- rbeta(n.sim, 50,  950)        # probability from healthy to sick
p.SD <- rbeta(n.sim, 100, 900)        # probability from sick to dead
c.H  <- rnorm(n.sim, 400, 50)         # cost of being in healthy state
c.S  <- rnorm(n.sim, 100, 80)         # cost of being in sick state
c.D  <- 0
u.H  <- rnorm(n.sim, 0.8, 0.02)       # utility of being in healthy state
u.S  <- rnorm(n.sim, 0.5, 0.02)       # utility of being in the sick state
u.D  <- 0
Markov2 <- function(i) {     # Markov model function - for different values of index i calculate the Markov model
m.P <- matrix(c(1 - p.HD[i] - p.HS[i],     p.HS[i], p.HD[i], # Transition matrix
0, 1 - p.SD[i], p.SD[i],
0,           0,      1),
byrow = T, nrow = n.s,
dimnames = list(v.n, v.n))
m.TR <- matrix(NA, nrow = n.t + 1, ncol = n.s,
dimnames = list(0:n.t, v.n))   # create Markov trace
m.TR[1, ] <- c(1, 0, 0)                       # initialize Markov trace
for (t in 1:n.t) {                            # throughout the number of cycles
m.TR[t + 1, ] <- m.TR[t, ] %*% m.P          # estimate the Markov trace for cycle t + 1
} # close the for loop for time
v.c <- c(c.H[i], c.S[i], c.D)   # vector of costs
v.e <- c(u.H[i], u.S[i], u.D)   # vector of utilities
TC4 <- t(m.TR %*% v.c) %*% v.dw  # total discounted cost for ith simulation
TE4 <- t(m.TR %*% v.e) %*% v.dw  # total discounted QALYs for ith simulation
c(TC4, TE4)                      # output of the function TC and TE
} # close the function
Results4 <- sapply(1:n.sim, Markov2) # apply the Markov2 functions for n.Sim times and save the results in the variable Results4
TC4      <- Results4[1, ] # save costs seperately
TE4      <- Results4[2, ] # save QALYs seperately
t4 <- Sys.time() - p4  # calculate time to run the analusis by extracting p from the current system time
################ probabilistic analysis method 1.1 ###########################################
#######    probabilistic analysis using difference equation method    ##################################
p5 <- Sys.time()          # save system time
set.seed(1)               # set a seed to be able to reproduce the same results
m.H <- m.S <- m.D <- matrix(0, nrow =  n.sim, ncol = n.t + 1,
dimnames = list(1:n.sim, 0:n.t))  #initialize state matrices
m.H[, 1] <- 1  # assign every proportion across the simulations to the Healthy state
for (t in 1:n.t) {
m.H[, t + 1] <-  m.H[, t] * (1 - p.HD - p.HS)                     # calculate the prop of cohort in healthy state at time t + 1
m.S[, t + 1] <-  m.S[, t] * (1 - p.SD)       + m.H[, t] * p.HS    # calculate the prop of cohort in sick state at time t + 1
m.D[, t + 1] <-  m.D[, t] + m.H[, t] * p.HD  + m.S[, t] * p.SD    # calculate the prop of cohort in dead state at time t + 1
} # close the loop for time
TE5 <- (u.H * m.H + u.S * m.S + u.D * m.D) %*% v.dw    # total discounted cost for all simulations
TC5 <- (c.H * m.H + c.S * m.S + c.D * m.D) %*% v.dw    # total discounted QALYs for all simulations
mean(TC5)   # calculate the mean discounted costs
mean(TE5)   # calculate the mean discounted QALYs
t5 <- Sys.time() - p5  # calculate time to run the analysis by extracting p from the current system time
TC <- cbind(TC1, TC2, TC3, TC4, TC5)  # create a variable with the total costs for the different methods
TE <- cbind(TE1, TE2, TE3, TE4, TE5)  # create a variable with the total utilities for the different methods
tt <- c(t1, t2, t3, t4, t5)           # count the time it took via each method
colMeans(TE)                          # calculate the column mean of the TE variable
tt                                    # print the time it took for each model to run
TE1 <- data.frame(TE1)
TC1 <- data.frame(TC1)
ScatterCE(strategies = "treatment", m.e = TE1, m.c = TC1) # cost-effectiveness plane
gen_psa <- function(n.sim = 1000, seed = 1){
set.seed(seed)              # set a seed to be able to reproduce the same results
df.psa <- data.frame(
p.HD <- rbeta(n.sim, 20,  980)                # probability from healthy to dead
p.HS <- rbeta(n.sim, 50,  950)                # probability from healthy to sick
p.SD <- rbeta(n.sim, 100, 900)                # probability from sick to dead
c.H  <- rnorm(n.sim, 400, 50)                 # cost of being in healthy state
c.S  <- rnorm(n.sim, 100, 80)                 # cost of being in sick state
c.D  <- 0
u.H  <- rnorm(n.sim, 0.8, 0.02)               # utility of being in healthy state
u.S  <- rnorm(n.sim, 0.5, 0.02)               # utility of being in the sick state
u.D  <- 0                                     # utility of being dead
)
return(df.psa)
}
gen_psa <- function(n.sim = 1000, seed = 1){
set.seed(seed)              # set a seed to be able to reproduce the same results
df.psa <- data.frame(
p.HD <- rbeta(n.sim, 20,  980),                # probability from healthy to dead
p.HS <- rbeta(n.sim, 50,  950),                # probability from healthy to sick
p.SD <- rbeta(n.sim, 100, 900),                # probability from sick to dead
c.H  <- rnorm(n.sim, 400, 50) ,                # cost of being in healthy state
c.S  <- rnorm(n.sim, 100, 80) ,                # cost of being in sick state
c.D  <- 0,
u.H  <- rnorm(n.sim, 0.8, 0.02),               # utility of being in healthy state
u.S  <- rnorm(n.sim, 0.5, 0.02),               # utility of being in the sick state
u.D  <- 0                                     # utility of being dead
)
return(df.psa)
}
gen_psa(10)
gen_psa <- function(n.sim = 1000, seed = 1){
set.seed(seed)              # set a seed to be able to reproduce the same results
df.psa <- data.frame(
p.HD = rbeta(n.sim, 20,  980),                # probability from healthy to dead
p.HS = rbeta(n.sim, 50,  950),                # probability from healthy to sick
p.SD = rbeta(n.sim, 100, 900),                # probability from sick to dead
c.H  = rnorm(n.sim, 400, 50) ,                # cost of being in healthy state
c.S  = rnorm(n.sim, 100, 80) ,                # cost of being in sick state
c.D  = 0,
u.H  = rnorm(n.sim, 0.8, 0.02),               # utility of being in healthy state
u.S  = rnorm(n.sim, 0.5, 0.02),               # utility of being in the sick state
u.D  = 0                                     # utility of being dead
)
return(df.psa)
}
gen_psa(10)
# wrap the code 3-state Markov model in a function
MM.3state <- function(params) {
with(as.list(params),
{
v.dw <- 1 / (1 + d.r) ^ (0:n.t) # calculate discount weight for each cycle based on discount rate d.r
# create a transition matrix
m.P  <- matrix(c(1 - p.HD - p.HS, p.HS, p.HD,
0,            1- p.SD, p.SD,
0,                  0,    1),
nrow = n.s, ncol = n.s, byrow = TRUE,
dimnames = list(v.n, v.n)) # name the columns and rows of the transition matrix
m.TR <- matrix(NA, nrow = n.t + 1 , ncol = n.s,
dimnames = list(0:n.t, v.n))     # create Markov trace (n.t + 1 because R doesn't understand  Cycle 0)
m.TR[1, ] <- c(1, 0, 0)                         # initialize Markov trace
####### PROCESS ###########################################
for (t in 1:n.t){                              # throughout the number of cycles
m.TR[t + 1, ] <- m.TR[t, ] %*% m.P           # estimate the Markov trace for cycle the next cycle (t + 1)
}
####### OUTPUT  ###########################################
# mean cost and QALYs per cycle
v.c <- m.TR %*% c(c.H, c.S, c.D)  # calculate expected costs by multiplying m.TR with the cost vector for the different health states
v.u <- m.TR %*% c(u.H, u.S, u.D)  # calculate expected QALYs by multiplying m.TR with the utilities for the different health states
# discounted QALYs and costs
TC <- t(v.c) %*% v.dw    # Discount costs by multiplying the cost vector with discount weights (v.dw)
TE <- t(v.u) %*% v.dw    # Discount QALYS by multiplying the QALYs vector with discount weights (v.dw)
results <- c(TC, TE)
return(results)
}
)
}
MM.3state(gen_psa(1))
MM.3state(gen_psa(1))
MM.3state(gen_psa(2))
df.psa <- gen_psa(1000)
n.sim <- 1000
n.sim <- 1000
df.psa <- gen_psa(nsim)
n.sim <- 1000
df.psa <- gen_psa(n.sim = nsim)
df.psa <- gen_psa(n.sim = n.sim)
df.out <- matrix(NaN, nrow = n.sim, ncol = 2)
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
# Try it
gen_psa(10) # Works!
# wrap the code 3-state Markov model in a function
MM.3state <- function(params) {
with(as.list(params),
{
v.dw <- 1 / (1 + d.r) ^ (0:n.t) # calculate discount weight for each cycle based on discount rate d.r
# create a transition matrix
m.P  <- matrix(c(1 - p.HD - p.HS, p.HS, p.HD,
0,            1- p.SD, p.SD,
0,                  0,    1),
nrow = n.s, ncol = n.s, byrow = TRUE,
dimnames = list(v.n, v.n)) # name the columns and rows of the transition matrix
m.TR <- matrix(NA, nrow = n.t + 1 , ncol = n.s,
dimnames = list(0:n.t, v.n))     # create Markov trace (n.t + 1 because R doesn't understand  Cycle 0)
m.TR[1, ] <- c(1, 0, 0)                         # initialize Markov trace
####### PROCESS ###########################################
for (t in 1:n.t){                              # throughout the number of cycles
m.TR[t + 1, ] <- m.TR[t, ] %*% m.P           # estimate the Markov trace for cycle the next cycle (t + 1)
}
####### OUTPUT  ###########################################
# mean cost and QALYs per cycle
v.c <- m.TR %*% c(c.H, c.S, c.D)  # calculate expected costs by multiplying m.TR with the cost vector for the different health states
v.u <- m.TR %*% c(u.H, u.S, u.D)  # calculate expected QALYs by multiplying m.TR with the utilities for the different health states
# discounted QALYs and costs
TC <- t(v.c) %*% v.dw    # Discount costs by multiplying the cost vector with discount weights (v.dw)
TE <- t(v.u) %*% v.dw    # Discount QALYS by multiplying the QALYs vector with discount weights (v.dw)
results <- c(TC, TE)
return(results)
}
)
}
p
p <- Sys.time()   # save system time
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
t.psa <- Sys.time() - p  # calculate time to run the analusis by extracting p from the current system time
t.psa
ScatterCE(strategies = "treatment", m.e = TE, m.c = TC) # cost-effectiveness plane
ScatterCE
TE
ScatterCE(strategies = "treatment", m.e = df.out[, "TE"], m.c = df.out[, "TC"]) # cost-effectiveness plane
## Initialize matrix of outcomes
df.out <- matrix(NaN, nrow = n.sim, ncol = 2)
colnames(df.out) <- c("TE", "TC")
## Run PSA
p <- Sys.time()   # save system time
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
t.psa <- Sys.time() - p  # calculate time to run the analusis by extracting p from the current system time
t.psa
p <- Sys.time()   # save system time
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
t.psa <- Sys.time() - p  # calculate time to run the analusis by extracting p from the current system time
t.psa
ScatterCE(strategies = "treatment", m.e = df.out[, "TE"], m.c = df.out[, "TC"]) # cost-effectiveness plane
df.out[, "TE"]
df.out[, "TC"]
## Initialize matrix of outcomes
df.out <- as.data.framd(matrix)(NaN, nrow = n.sim, ncol = 2)
colnames(df.out) <- c("TE", "TC")
## Run PSA
p <- Sys.time()   # save system time
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
t.psa <- Sys.time() - p  # calculate time to run the analusis by extracting p from the current system time
t.psa
ScatterCE(strategies = "treatment", m.e = df.out[, "TE"], m.c = df.out[, "TC"]) # cost-effectiveness plane
str(df.out)
str(df.out[, "TE"])
ScatterCE(strategies = "treatment", m.e = matrix(df.out[, "TE"]), m.c = matrix(df.out[, "TC"])) # cost-effectiveness plane
ScatterCE(strategies = "treatment", m.e = data.frame(df.out[, "TE"]), m.c = data.frame(df.out[, "TC"])) # cost-effectiveness plane
rm(list = ls())  # Delete everything that is in R's memory
setwd("/Users/fernando/Google Drive/Projects/DARTH/Courses/2018 Merck/Participants's material/Day 2/") # set the working directory
source("Functions/PSA_functions.R") # load custom made functions for PSA
######################### INPUT PARAMETERS ####################################
v.n   <- c("Healthy", "Sick", "Dead")     # state names
n.s   <- length(v.n)                      # number of states
n.t   <- 60                               # number of cycles
n.sim <- 10000                            # number of simulations
d.r   <- 0.03                             # discount rate
v.dw  <- 1 / (1 + d.r) ^ (0:n.t)          # calculate discount weight for each cycle based on discount rate d.r
## Function that generates random sample for PSA
gen_psa <- function(n.sim = 1000, seed = 1){
set.seed(seed)              # set a seed to be able to reproduce the same results
df.psa <- data.frame(
p.HD = rbeta(n.sim, 20,  980),                # probability from healthy to dead
p.HS = rbeta(n.sim, 50,  950),                # probability from healthy to sick
p.SD = rbeta(n.sim, 100, 900),                # probability from sick to dead
c.H  = rnorm(n.sim, 400, 50) ,                # cost of being in healthy state
c.S  = rnorm(n.sim, 100, 80) ,                # cost of being in sick state
c.D  = 0,
u.H  = rnorm(n.sim, 0.8, 0.02),               # utility of being in healthy state
u.S  = rnorm(n.sim, 0.5, 0.02),               # utility of being in the sick state
u.D  = 0                                     # utility of being dead
)
return(df.psa)
}
# Try it
gen_psa(10) # Works!
# wrap the code 3-state Markov model in a function
MM.3state <- function(params) {
with(as.list(params),
{
v.dw <- 1 / (1 + d.r) ^ (0:n.t) # calculate discount weight for each cycle based on discount rate d.r
# create a transition matrix
m.P  <- matrix(c(1 - p.HD - p.HS, p.HS, p.HD,
0,            1- p.SD, p.SD,
0,                  0,    1),
nrow = n.s, ncol = n.s, byrow = TRUE,
dimnames = list(v.n, v.n)) # name the columns and rows of the transition matrix
m.TR <- matrix(NA, nrow = n.t + 1 , ncol = n.s,
dimnames = list(0:n.t, v.n))     # create Markov trace (n.t + 1 because R doesn't understand  Cycle 0)
m.TR[1, ] <- c(1, 0, 0)                         # initialize Markov trace
####### PROCESS ###########################################
for (t in 1:n.t){                              # throughout the number of cycles
m.TR[t + 1, ] <- m.TR[t, ] %*% m.P           # estimate the Markov trace for cycle the next cycle (t + 1)
}
####### OUTPUT  ###########################################
# mean cost and QALYs per cycle
v.c <- m.TR %*% c(c.H, c.S, c.D)  # calculate expected costs by multiplying m.TR with the cost vector for the different health states
v.u <- m.TR %*% c(u.H, u.S, u.D)  # calculate expected QALYs by multiplying m.TR with the utilities for the different health states
# discounted QALYs and costs
TC <- t(v.c) %*% v.dw    # Discount costs by multiplying the cost vector with discount weights (v.dw)
TE <- t(v.u) %*% v.dw    # Discount QALYS by multiplying the QALYs vector with discount weights (v.dw)
results <- c(TC, TE)
return(results)
}
)
}
## Draw random sample for PSA
df.psa <- gen_psa(n.sim = n.sim)
## Initialize matrix of outcomes
df.out <- matrix(NaN, nrow = n.sim, ncol = 2)
colnames(df.out) <- c("TE", "TC")
## Run PSA
p <- Sys.time()   # save system time
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
t.psa <- Sys.time() - p  # calculate time to run the analusis by extracting p from the current system time
t.psa
ScatterCE(strategies = "treatment", m.e = data.frame(df.out[, "TE"]), m.c = data.frame(df.out[, "TC"])) # cost-effectiveness plane
rm(list = ls())  # Delete everything that is in R's memory
setwd("/Users/fernando/Google Drive/Projects/DARTH/Courses/2018 Merck/Participants's material/Day 2/") # set the working directory
source("Functions/PSA_functions.R") # load custom made functions for PSA
######################### INPUT PARAMETERS ####################################
v.n   <- c("Healthy", "Sick", "Dead")     # state names
n.s   <- length(v.n)                      # number of states
n.t   <- 60                               # number of cycles
n.sim <- 10000                            # number of simulations
d.r   <- 0.03                             # discount rate
v.dw  <- 1 / (1 + d.r) ^ (0:n.t)          # calculate discount weight for each cycle based on discount rate d.r
####### probabilistic analysis ###########################################
####### Create function that generate random samples ###################################
## Function that generates random sample for PSA
gen_psa <- function(n.sim = 1000, seed = 1){
set.seed(seed)              # set a seed to be able to reproduce the same results
df.psa <- data.frame(
p.HD = rbeta(n.sim, 20,  980),                # probability from healthy to dead
p.HS = rbeta(n.sim, 50,  950),                # probability from healthy to sick
p.SD = rbeta(n.sim, 100, 900),                # probability from sick to dead
c.H  = rnorm(n.sim, 400, 50) ,                # cost of being in healthy state
c.S  = rnorm(n.sim, 100, 80) ,                # cost of being in sick state
c.D  = 0,
u.H  = rnorm(n.sim, 0.8, 0.02),               # utility of being in healthy state
u.S  = rnorm(n.sim, 0.5, 0.02),               # utility of being in the sick state
u.D  = 0                                     # utility of being dead
)
return(df.psa)
}
# Try it
gen_psa(10) # Works!
# wrap the code 3-state Markov model in a function
MM.3state <- function(params) {
with(as.list(params),
{
v.dw <- 1 / (1 + d.r) ^ (0:n.t) # calculate discount weight for each cycle based on discount rate d.r
# create a transition matrix
m.P  <- matrix(c(1 - p.HD - p.HS, p.HS, p.HD,
0,            1- p.SD, p.SD,
0,                  0,    1),
nrow = n.s, ncol = n.s, byrow = TRUE,
dimnames = list(v.n, v.n)) # name the columns and rows of the transition matrix
m.TR <- matrix(NA, nrow = n.t + 1 , ncol = n.s,
dimnames = list(0:n.t, v.n))     # create Markov trace (n.t + 1 because R doesn't understand  Cycle 0)
m.TR[1, ] <- c(1, 0, 0)                         # initialize Markov trace
####### PROCESS ###########################################
for (t in 1:n.t){                              # throughout the number of cycles
m.TR[t + 1, ] <- m.TR[t, ] %*% m.P           # estimate the Markov trace for cycle the next cycle (t + 1)
}
####### OUTPUT  ###########################################
# mean cost and QALYs per cycle
v.c <- m.TR %*% c(c.H, c.S, c.D)  # calculate expected costs by multiplying m.TR with the cost vector for the different health states
v.u <- m.TR %*% c(u.H, u.S, u.D)  # calculate expected QALYs by multiplying m.TR with the utilities for the different health states
# discounted QALYs and costs
TC <- t(v.c) %*% v.dw    # Discount costs by multiplying the cost vector with discount weights (v.dw)
TE <- t(v.u) %*% v.dw    # Discount QALYS by multiplying the QALYs vector with discount weights (v.dw)
results <- c(TC, TE)
return(results)
}
)
}
## Draw random sample for PSA
df.psa <- gen_psa(n.sim = n.sim)
## Initialize matrix of outcomes
df.out <- matrix(NaN, nrow = n.sim, ncol = 2)
colnames(df.out) <- c("TE", "TC")
## Run PSA
p <- Sys.time()   # save system time
for(i in 1:n.sim){
df.out[i,] <- MM.3state(df.psa[i, ])
cat('\r', paste(round(i/n.sim * 100), "% done", sep = " "))       # display the progress of the simulation
}
t.psa <- Sys.time() - p  # calculate time to run the analusis by extracting p from the current system time
t.psa
ScatterCE(strategies = "treatment", m.e = data.frame(df.out[, "TE"]), m.c = data.frame(df.out[, "TC"])) # cost-effectiveness plane
#install.packages(c("reshape2"), dependencies = TRUE)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
source("../Functions/PSA_functions.R") # load custom made functions for PSA
## Clean everything from workspace
rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
source("../Functions/PSA_functions.R")
source("../Functions/DARTH_functions.R")
## Load PSA file
# Read the `.csv` simulation file into `R`.
psa <- read.csv("PSA/PSA.csv",header=TRUE)
## Create initial parameters
# Label our strategies
Strategies <- c("Chemo", "Radio", "Surgery")
## Load PSA file
# Read the `.csv` simulation file into `R`.
psa <- read.csv("PSA.csv",header=TRUE)
## Create initial parameters
# Label our strategies
Strategies <- c("Chemo", "Radio", "Surgery")
Strategies
#Determine the number of strategies
numDepVars <- length(Strategies)
numDepVars
## Create matrices for outcomes and parameters
Outcomes <- psa[, 2:(numDepVars*2+1)]
head(Outcomes)
#######################################
#### Displaying results from a PSA ####
#######################################
### Define matrices with costs and effects
m.c <- Outcomes[, c(1, 3, 5)]
m.e <- Outcomes[, c(2, 4, 6)]
### Define a vector of WTP values
v.wtp <- c(1, seq(5000, 150000, length.out = 31))
### CE Plane and Scatter Plot ####
library(plyr) # To summarize data and generate new data frames from it
library(ggplot2)
ce.mat <- cbind(Strategies = 1:3,
Costs = colMeans(m.c),
Effectviness = colMeans(m.e))
ce.front <- getFrontier(ce.mat, plot = F)
ce.df <- data.frame(Strategy = Strategies,
Cost = colMeans(m.c),
Effectiveness = colMeans(m.e))
plotFrontier(CEmat = ce.df, frontier = ce.front)
### With user-written function & ggplot2
ScatterCE(strategies = Strategies,
m.c = m.c,
m.e = m.e)
#### CEA Curves (CEAC) & Frontier (CEAF) ####
ceaf(v.wtp = v.wtp,
strategies = Strategies,
m.e = m.e,
m.c = m.c)
#### Expected Loss Curves (ELC) ####
elc(v.wtp = v.wtp,
strategies = Strategies,
m.e = m.e,
m.c = m.c)
#### Expected Loss Curves (ELC) ####
evpi(v.wtp = v.wtp,
m.e = m.e,
m.c = m.c)