数据资源下载地址：https://download.csdn.net/download/weixin_68126662/89101333
R代码参考：
#############################################################################
###
### GLBH0046: Advanced Economic Evaluation
###
### Practical 1: Resource use and costs
###
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 1-2")
# Load the data (using the 'haven' library, which reads Stata data files directly)
library(haven)
data<- read_dta('evar.dta')
# Familiarise with dataset
names(data) # variable names
print(data) # visualise data (first 10 rows)
summary(data) # basic descritptive stats for all the variables
### Question 3
#define unit costs
device_evar <- 5500
device_or <- 700
drug_evar <- 600
drug_or <- 400
oheads <- 2
staff_evar <- 12
staff_or <- 9
cc <- 1400
gm <- 260
out <- 140
gp <- 60
nurse <- 20
icare <- 160
prod <- 65
# Consumables - assign medical device and drug costs to individuals by arm
device_cost <- ifelse(data$arm==1, device_evar, device_or)
drug_cost <- ifelse(data$arm==1, drug_evar, drug_or)
data$cons_cost <- device_cost + drug_cost
# Theatre costs - combine time (in min) spent in theatre with unit costs
overhead_cost <- data$theatre_los*oheads
staff_cost <- ifelse(data$arm==1, data$theatre_los*staff_evar, data$theatre_los*staff_or)
data$theatre_cost <- overhead_cost + staff_cost
# Critical care cost - combine bed-days spent in critical care with unit cost
data$cc_cost <-data$cc_los*cc
# General medical ward cost - combine bed-days spent in general medical ward with unit cost
data$gm_cost <-data$gm_los*gm
# Outpatient cost - combine outpatient visits with unit cost
data$out_cost <-data$out_los*out
# GP cost - combine GP and nurse visits with unit cost
gp_cost <- data$gp_visit*gp
nurse_cost <- data$nurse_visit*nurse
data$primary_cost <- gp_cost + nurse_cost
# Summarise costs across components using with either 'tapply' or 'aggregate' functions
#tapply(data$theatre_cost, data$arm, function(x) c(mean(x), sd(x)))
tc<-with(data, aggregate(data$theatre_cost ~ data$arm, FUN=function(x) c(mean=mean(x),SD=sd(x) )))
cc<-with(data, aggregate(data$cc_cost ~ data$arm, FUN=function(x) c(mean=mean(x),SD=sd(x) )))
gc<-with(data, aggregate(data$gm_cost ~ data$arm, FUN=function(x) c(mean=mean(x),SD=sd(x) )))
oc<-with(data, aggregate(data$out_cost ~ data$arm, FUN=function(x) c(mean=mean(x),SD=sd(x) )))
pc<-with(data, aggregate(data$primary_cost ~ data$arm, FUN=function(x) c(mean=mean(x),SD=sd(x) )))
#you may wish to combine results in a table
table <- rbind(t(tc), t(cc)[2:3,], t(gc)[2:3,], t(oc)[2:3,], t(pc)[2:3,])
table
# More compact using 'dplyr' library (advanced)
library(dplyr)
data %>%
group_by(arm) %>%
summarise(across(
.cols = theatre_cost:primary_cost,
.fns = list(Mean = mean, SD = sd), na.rm = TRUE,
.names = "{col}_{fn}"
))
### Question 5
# Calculate total cost
data$total_cost <- data$cons_cost + data$theatre_cost + data$cc_cost +
data$gm_cost + data$out_cost + data$primary_cost
# Report mean difference (and 95% CI) in total cost for EVAR and Open Repair
tapply(data$total_cost, data$arm, function(x) c(mean(x), sd(x)))
lm(data$total_cost ~ data$arm)
t.test(data$total_cost ~ data$arm, conf.level = 0.95)
### Question 6
# Informal care cost - combine number days of informal care with unit cost
data$inf_cost <-data$carer*icare
# Productivity losses - combine number days off work with unit cost
data$work_cost <-data$off_work*prod
# Calculate total cost under societal perspective
data$total_cost_soc <- data$total_cost + data$inf_cost + data$work_cost
# Report mean difference (and 95% CI) in total cost for EVAR and Open Repair
tapply(data$total_cost_soc, data$arm, function(x) c(mean(x), sd(x)))
lm(data$total_cost_soc ~ data$arm)
t.test(data$total_cost_soc ~ data$arm, conf.level = 0.95)
### Question 7
# plot distribution (histogram) of total healthcare costs by trial arm
library(ggplot2)
ggplot(data, aes(x = total_cost)) +
geom_histogram(fill = "white", colour = "black") +
facet_grid(data$arm ~ .)
# Alternative forms to address the skewness
# Parametric
# GLM, considering a distribution that better fits the skewed costs (e.g. Gamma)
gamma.mod<-glm(data$total_cost ~ data$arm, family=Gamma(link="identity"))
cbind(coef(gamma.mod), confint(gamma.mod))
qqnorm(residuals(gamma.mod), main="")
qqline(residuals(gamma.mod))
# Non-parametric
# Non-parametric bootstrap using 'boot' package
library(boot)
# create function that estimates parameter of interest
inc.cost <- function(data, indices) {
d <- data[indices,] # allows boot to sample from the original data
lm <- lm(total_cost ~ arm, data=d)
inc.cost<-coef(lm)[2]
return(inc.cost)
}
# run R=1000 bootstrap replications
set.seed(90210) # set seed to ensure same results after re-runs
b <- boot(data=data, statistic=inc.cost, R=1000)
plot(b)
mean(b$t) # incremental cost is just the average of the Bootstrap replicates
# Report 95% CI using 'boot.ci'
# It takes the 2.5% and 97.5% percentiles of the distribution of Bootstrap replicates
# Here we consider bias-corrected and accelerated (BCa) percentiles (more appropriate for skewed data)
boot.ci(b, type="bca")
# save the dataset (cost data will be needed for next practical)
write_dta(data,'costs.dta')
#############################################################################
###
### GLBH0046: Advanced Economic Evaluation
###
### Practical 2: QALYs
###
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 1-2")
# Load the data (using the 'haven' library, which reads Stata data files directly)
library(haven)
data<- read_dta('evar.dta')
### QUESTION 2
library(eq5d)
library(ggplot2)
#Estimate eq5d score at 1 month
data$eq5d_1mo <- eq5d(scores= data.frame(
MO=c(data$q1_1mo),
SC=c(data$q2_1mo),
UA=c(data$q3_1mo),
PD=c(data$q4_1mo),
AD=c(data$q5_1mo)),
country="UK", version="3L", type="TTO", ignore.invalid=TRUE) #ignore NAs
ggplot(data, aes(x = eq5d_1mo)) + # histogram of EQ-5D score at 1 month
geom_histogram(fill = "white", colour = "black") +
facet_grid(arm ~ .)
#Estimate eq5d score at 3 months
data$eq5d_3mo <- eq5d(scores= data.frame(
MO=c(data$q1_3mo),
SC=c(data$q2_3mo),
UA=c(data$q3_3mo),
PD=c(data$q4_3mo),
AD=c(data$q5_3mo)),
country="UK", version="3L", type="TTO", ignore.invalid=TRUE) #ignore NAs
#Estimate eq5d score at 12 months
data$eq5d_12mo <- eq5d(scores= data.frame(
MO=c(data$q1_12mo),
SC=c(data$q2_12mo),
UA=c(data$q3_12mo),
PD=c(data$q4_12mo),
AD=c(data$q5_12mo)),
country="UK", version="3L", type="TTO", ignore.invalid=TRUE) #ignore NAs
### QUESTION 3
# Summarise EQ-5D-3L for survivors (i.e. ignore NAs for deceased patients)
# 1 month
tapply(data$eq5d_1mo[data$death30days==0], data$arm[data$death30days==0], function(x) c(mean(x), sd(x)))
lm1<-lm(data$eq5d_1mo ~ data$arm)
cbind(coef(lm1), confint(lm1))
# 3 months
tapply(data$eq5d_3mo[data$death3mnth==0], data$arm[data$death3mnth==0], function(x) c(mean(x), sd(x)))
lm3<-lm(data$eq5d_3mo ~ data$arm)
cbind(coef(lm3), confint(lm3))
# 12 months
tapply(data$eq5d_12mo[data$death1year==0], data$arm[data$death1year==0], function(x) c(mean(x), sd(x)))
lm12<-lm(data$eq5d_12mo ~ data$arm)
cbind(coef(lm12), confint(lm12))
### QUESTION 4
# 1 month
# Report number of deceased per treatment arm
tapply(data$death30days, data$arm, function(x) c(sum(x), sum(x)/length(x)))
# run logistic regression and report odds-ratio (exponentiate results)
logist1 <- glm(data$death30days ~ data$arm, data = data, family = "binomial"(link="logit"))
exp(cbind(coef(logist1), confint(logist1)))
# 3 months
tapply(data$death3mnth, data$arm, function(x) c(sum(x), sum(x)/length(x)))
logist3 <- glm(data$death3mnth ~ data$arm, data = data, family = "binomial"(link="logit"))
exp(cbind(coef(logist3), confint(logist3)))
# 12 months
tapply(data$death1year, data$arm, function(x) c(sum(x), sum(x)/length(x)))
logist12 <- glm(data$death1year ~ data$arm, data = data, family = "binomial"(link="logit"))
exp(cbind(coef(logist12), confint(logist12)))
# Plot survival (Kaplan-Meier) curves
library(survival)
library(survminer)
data$ttd[data$ttd>365]<-365.25 # focus on time till death or censoring (ttd) within 12 months
kmfit <- survfit( Surv(data$ttd, data$death1year==1) ~ data$arm)
# plot survival curves using 'ggsurvplot'
ggsurvplot(
kmfit, # survfit object with calculated KM statistics
data = data,
risk.table = TRUE, # show no. patients at risk
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for point estimates of survival curves.
xlim = c(0,366), # edit X axis
)
### QUESTION 5
# Calculate QALYs for individuals who survive up to 12 months (assuming eq-5d is zero at baseline)
data$qaly <- 0.125*data$eq5d_1mo + (5.5/12)*data$eq5d_3mo + 0.375*data$eq5d_12mo
# Replace QALYs=0 for individuals who die between baseline and 1st month
data$qaly[data$death30days==1] <- 0
# Calculate QALYs for individuals who die between 1 and 3 months
data$qaly[data$death30days==0 & data$death3mnth==1] <- (data$ttd[data$death30days==0 & data$death3mnth==1]/365)*0.5*data$eq5d_1mo[data$death30days==0 & data$death3mnth==1]
# Calculate QALYs for individuals who die between 3 and 12 months
data$qaly[data$death3mnth==0 & data$death1year==1] <- 0.125*data$eq5d_1mo[data$death3mnth==0 & data$death1year==1] +
(data$ttd[data$death3mnth==0 & data$death1year==1]/365-(1/12))*0.5*data$eq5d_3mo[data$death3mnth==0 & data$death1year==1]
# Plot QALYs
ggplot(data, aes(x = qaly)) +
geom_histogram(fill = "white", colour = "black") +
facet_grid(arm ~ .)
### QUESTION 6
# Incremental QALY at 12 months
tapply(data$qaly, data$arm, function(x) c(mean(x), sd(x)))
# estimate 95% CI via Bootstrap
library(boot)
inc.qaly <- function(data, indices) {
d <- data[indices,] # allows boot to sample from the original data
lm <- lm(qaly ~ arm, data=d)
inc.qaly<-coef(lm)[2]
return(inc.qaly)
}
set.seed(90210) # set seed to ensure same results after re-runs
b <- boot(data=data, statistic=inc.qaly, R=1000)
mean(b$t) # incremental qaly
boot.ci(b, type="bca")
### QUESTION 7
# Cost-effectiveness at 12 months
# add total cost variable to main dataset (make sure both datasets are sorted by 'patientid')
data.cost<-read_dta('costs.dta')
data$total_cost<-data.cost$total_cost
# create function that extracts both incremental cost and QALY
cea <- function(data, indices) {
d <- data[indices,] # allows boot to sample from the original data
inc.cost<-coef(lm(total_cost ~ arm, data=d))[2]
inc.qaly<-coef(lm(qaly ~ arm, data=d))[2]
return(c(inc.qaly,inc.cost))
}
# run bootstrap
set.seed(90210) # set seed to ensure same results after re-runs
b <- boot(data=data, statistic=cea, R=1000)
mean(b$t[,2])/mean(b$t[,1]) # ICER
mean(b$t[,1]*30000-b$t[,2]) # INB (assuming a threshold of ?30,000 per QALY)
boot.cea<-cbind(b$t[,1], b$t[,2]) # bootstrap sample of incremental QALY and cost
# cost-effectiveness plane
par(las=1)
plot(boot.cea, main=" ", xlab="", ylab="",
xlim=c(-.1,.2), ylim=c(-6000,6000),pch=19, cex=0.7, col="darkslategray", axes=F )
points(x=mean(b$t[,1]), y=mean(b$t[,2]), col='red', cex=1, pch=16)
axis(1, pos=0, cex.axis=1)
axis(2, pos=0, cex.axis=1)
abline(0,30000, lty='dotted')
text(0.15, 3500, "WTP ?30,000", cex=1)
mtext(side=1, "Incremental QALY", cex=1)
mtext(side=2, "Incremental cost", cex=1, las=0)
#############################################################################
###
### GLBH0046: Advanced Economic Evaluation
###
### Practical 3: Markov model set-up
###
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 4-7")
### QUESTION 6
### 6.1 - create vectors for key inputs
# treatments
treat <- c("EVAR", "OR") # treatments
n_treat <- length(treat) # no. treatments
# health states (E: event-free survival, C: CVD event, D: Death)
states <- c("E", "C", "D")
n_states <- length(states) # no. states
# cohort
n <- "size" # size of the cohort
# cycles
n_cycles <- 10 # no. cycles (*to be modified)
### 6.2 costs for each health state
c_E <- "cost_E" # cost in E state
c_C <- "cost_C" # cost in C state
c_D <- "cost_D" # cost in D state
c_evar <- "cost_evar" # Cost of EVAR treatment (one-off)
c_or <- "cost_or" # Cost of Open Repair treatment (one-off)
### 6.2.1 costs matrix across states and treatment groups
c_matrix <- matrix(
c(c_E, c_C, c_D,
c_E, c_C, c_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
### 6.3 utilities for each health state
u_E <- "utility_E" # utility in E state
u_C <- "utility_C" # utility in C state
u_D <- "utility_D" # utility in D state
### 6.3.1 utilities matrix across states and treatment groups
u_matrix <- matrix(
c(u_E, u_C, u_D,
u_E, u_C, u_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
### 6.4 transition probabilities between health states
p_EC <- "prob_EC" # transition probability from E to C state
p_CE <- "prob_CE" # transition probability from C to E state
p_ED <- "prob_ED" # transition probability from E to D state
p_CD <- "prob_CD" # transition probability from C to D state
effect<- "treat_effect" # effect of treatment (reduces progression)
# 6.4.1 create matrix with all transition rates per treatment arm
# initialise transition matrix
p_mat <- array(data = NA,
dim = c(n_states, n_states, n_treat),
dimnames = list(from = states,
to = states,
treat))
# Fill in the matrix (Open Repair)
# From E
p_mat["E", "E", "OR"] <- "1 - p_EC - p_ED"
p_mat["E", "C", "OR"] <- "p_EC"
p_mat["E", "D", "OR"] <- "p_ED"
# From C
p_mat["C", "E", "OR"] <- "p_CE"
p_mat["C", "C", "OR"] <- "1 - p_CD - p_CE"
p_mat["C", "D", "OR"] <- "p_CD"
# From D
p_mat["D", "E", "OR"] <- 0
p_mat["D", "C", "OR"] <- 0
p_mat["D", "D", "OR"] <- 1
# Fill in the matrix (EVAR)
# From E
p_mat["E", "E", "EVAR"] <- "1 - p_EC*(1-effect) - p_ED"
p_mat["E", "C", "EVAR"] <- "p_EC*(1-effect)"
p_mat["E", "D", "EVAR"] <- "p_ED"
# From C
p_mat["C", "E", "EVAR"] <- "p_CE"
p_mat["C", "C", "EVAR"] <- "1 - p_CD - p_CE"
p_mat["C", "D", "EVAR"] <- "p_CD"
# From D
p_mat["D", "E", "EVAR"] <- 0
p_mat["D", "C", "EVAR"] <- 0
p_mat["D", "D", "EVAR"] <- 1
### 6.5 Create empty arrays to store model output
# Matrix to store numbers of individuals in each health state across cycles
pop <- array(data = NA,
dim = c(n_states, n_cycles+1, n_treat),
dimnames = list(state = states,
cycle = NULL,
treatment = treat))
# first cycle
pop["E", cycle = 1, ] <- n
pop["C", cycle = 1, ] <- 0
pop["D", cycle = 1, ] <- 0
# Matrix to store cost-effectiveness endpoints
costs <- LYs <- QALYs <- array(data = NA,
dim = c(n_treat, n_cycles+1),
dimnames = list(treatment = treat,
cycle = NULL))
# vector to store total life years (LYs), QALYs and costs
total_LYs <- total_QALYs <- total_costs <- setNames(c(NA, NA), treat)
### 6.6 set up loops to run the model
# outer loop - loop over treatment groups
for (i in 1:n_treat) {
# inner loop - loop over cycles
for (j in 2:n_cycles) { # starts in cycle 2 as cycle 1 is already defined above
pop[, cycle = j, treatment = i] <-
pop[, cycle = j - 1, treatment = i] %*% p_mat[, , treatment = i]
}
}
#############################################################################
###
### GLBH0046: Advanced Economic Evaluation
###
### Practical 4: Populate Markov model
###
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 4-7")
### QUESTION 5
# Populate Markov model
# treatments
treat <- c("EVAR", "OR") # treatments
n_treat <- length(treat) # no. treatments
# health states (E: event-free survival, C: CVD event, D: Death)
states <- c("E", "C", "D")
n_states <- length(states) # no. states
# cohort
n <- 1000 # size of the cohort
age0 <- 55 # initial age of cohort
# cycles
n_cycles <- 43 # no. cycles
# Discount rate (same for costs and QALYs)
dr <- 0.035
# costs for each health state
c_E <- 500 # cost in E state
c_C <- 2500 # cost in C state
c_D <- 0 # cost in D state
c_evar <- 18000 # Cost of EVAR treatment (one-off)
c_or <- 15500 # Cost of Open Repair treatment (one-off)
# costs matrix across states and treatment groups
c_matrix <- matrix(
c(c_E, c_C, c_D,
c_E, c_C, c_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# utilities for each health state
u_E <- 0.80 # utility in E state
u_C <- 0.60 # utility in C state
u_D <- 0 # utility in D state
# utilities matrix across states and treatment groups
u_matrix <- matrix(
c(u_E, u_C, u_D,
u_E, u_C, u_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# annual transition probabilities between health states
p_EC <- 0.05 # transition probability from E to C state
p_CE <- 0.22 # transition probability from C to E state
p_ED <- 0.04 # transition probability from E to D state
p_CD <- 0.14 # transition probability from C to D state
effect<- 0.1 # effect of treatment (reduces progression)
# 6.4.1 create matrix with all transition rates per treatment arm
# initialise transition matrix
p_mat <- array(data = NA,
dim = c(n_states, n_states, n_treat),
dimnames = list(from = states,
to = states,
treat))
# Fill in the matrix (Open Repair)
# From E
p_mat["E", "E", "OR"] <- 1 - p_EC - p_ED
p_mat["E", "C", "OR"] <- p_EC
p_mat["E", "D", "OR"] <- p_ED
# From C
p_mat["C", "E", "OR"] <- p_CE
p_mat["C", "C", "OR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "OR"] <- p_CD
# From D
p_mat["D", "E", "OR"] <- 0
p_mat["D", "C", "OR"] <- 0
p_mat["D", "D", "OR"] <- 1
# Fill in the matrix (EVAR)
# From E
p_mat["E", "E", "EVAR"] <- 1 - p_EC*(1-effect) - p_ED
p_mat["E", "C", "EVAR"] <- p_EC*(1-effect)
p_mat["E", "D", "EVAR"] <- p_ED
# From C
p_mat["C", "E", "EVAR"] <- p_CE
p_mat["C", "C", "EVAR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "EVAR"] <- p_CD
# From D
p_mat["D", "E", "EVAR"] <- 0
p_mat["D", "C", "EVAR"] <- 0
p_mat["D", "D", "EVAR"] <- 1
# Create empty matrices to store model output
# Matrix to store numbers of individuals in each health state across cycles
pop <- array(data = NA,
dim = c(n_states, n_cycles+1, n_treat),
dimnames = list(state = states,
cycle = NULL,
treatment = treat))
# first cycle
pop["E", cycle = 1, ] <- n
pop["C", cycle = 1, ] <- 0
pop["D", cycle = 1, ] <- 0
# Matrix to store cost-effectiveness endpoints
costs <- LYs <- QALYs <- array(data = NA,
dim = c(n_treat, n_cycles+1),
dimnames = list(treatment = treat,
cycle = NULL))
# vector to store total life years (LYs), QALYs and costs
total_LYs <- total_QALYs <- total_costs <- setNames(c(NA, NA), treat)
### Run the model
# outer loop - loop over treatment groups
for (i in 1:n_treat) {
# inner loop - loop over cycles
for (j in 1:n_cycles) {
pop[, cycle = j+1, treatment = i] <- # starts in cycle 2 as cycle 1 is already defined above
pop[, cycle = j, treatment = i] %*% p_mat[, , treatment = i]
} # end of inner loop
# Life-years - 1 if person alive (i.e states E and C), 0 otherwise (state D)
LYs[i, ] <- (c(1,1,0) %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles) # cycle 1 = Year 0
# Quality-adjusted life-years (LYs adjusted by the quality of life)
QALYs[i, ] <- (u_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
# Costs
costs[i, ] <-(c_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
total_LYs[i] <- sum(LYs[treatment = i, -1]) # -1 means 'exclude Year 0'
total_QALYs[i] <- sum(QALYs[treatment = i, -1])
total_costs[i] <- sum(costs[treatment = i, -1])
} # end of outer loop
# Add one-off cost (per patient) of the intervention
total_costs<-c(total_costs["EVAR"] + c_evar*n, total_costs["OR"] + c_or*n)
### Question 6
# incremental LY
inc_ly <- total_LYs["EVAR"] - total_LYs["OR"]
# incremental QALY
inc_qaly <- total_QALYs["EVAR"] - total_QALYs["OR"]
# Incremental costs
inc_cost <- total_costs["EVAR"] - total_costs["OR"]
# Incremental net benefit, assuming a WTP threshold of �30,000 per QALY4
inb<-inc_qaly*30000-inc_cost
# Incremental cost-effectiveness ratio
icer <- inc_cost/inc_qaly
total_LYs
total_QALYs
total_costs
inc_ly
inc_qaly
inc_cost
inb
Icer
#############################################################################
###
### GLBH0046: Advanced Economic Evaluation
###
### Practical 5: Deterministic sensitivity analysis
###
#############################################################################
#############################################################################
#
### Question 4
#
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 3-6")
# Re-run Markov model with age0=75 and n_cycles=23
# treatments
treat <- c("EVAR", "OR") # treatments
n_treat <- length(treat) # no. treatments
# health states (E: event-free survival, C: CVD event, D: Death)
states <- c("E", "C", "D")
n_states <- length(states) # no. states
# cohort
n <- 1000 # size of the cohort
age0 <- 75 # initial age of cohort
# cycles
n_cycles <- 23 # no. cycles
# Discount rate (same for costs and QALYs)
dr <- 0.035
# costs for each health state
c_E <- 500 # cost in E state
c_C <- 2500 # cost in C state
c_D <- 0 # cost in D state
c_evar <- 18000 # Cost of EVAR treatment (one-off)
c_or <- 15500 # Cost of Open Repair treatment (one-off)
# costs matrix across states and treatment groups
c_matrix <- matrix(
c(c_E, c_C, c_D,
c_E, c_C, c_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# utilities for each health state
u_E <- 0.80 # utility in E state
u_C <- 0.60 # utility in C state
u_D <- 0 # utility in D state
# utilities matrix across states and treatment groups
u_matrix <- matrix(
c(u_E, u_C, u_D,
u_E, u_C, u_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# annual transition probabilities between health states
p_EC <- 0.05 # transition probability from E to C state
p_CE <- 0.22 # transition probability from C to E state
p_ED <- 0.04 # transition probability from E to D state
p_CD <- 0.14 # transition probability from C to D state
effect<- 0.1 # effect of treatment (reduces progression)
# 6.4.1 create matrix with all transition rates per treatment arm
# initialise transition matrix
p_mat <- array(data = NA,
dim = c(n_states, n_states, n_treat),
dimnames = list(from = states,
to = states,
treat))
# Fill in the matrix (Open Repair)
# From E
p_mat["E", "E", "OR"] <- 1 - p_EC - p_ED
p_mat["E", "C", "OR"] <- p_EC
p_mat["E", "D", "OR"] <- p_ED
# From C
p_mat["C", "E", "OR"] <- p_CE
p_mat["C", "C", "OR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "OR"] <- p_CD
# From D
p_mat["D", "E", "OR"] <- 0
p_mat["D", "C", "OR"] <- 0
p_mat["D", "D", "OR"] <- 1
# Fill in the matrix (EVAR)
# From E
p_mat["E", "E", "EVAR"] <- 1 - p_EC*(1-effect) - p_ED
p_mat["E", "C", "EVAR"] <- p_EC*(1-effect)
p_mat["E", "D", "EVAR"] <- p_ED
# From C
p_mat["C", "E", "EVAR"] <- p_CE
p_mat["C", "C", "EVAR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "EVAR"] <- p_CD
# From D
p_mat["D", "E", "EVAR"] <- 0
p_mat["D", "C", "EVAR"] <- 0
p_mat["D", "D", "EVAR"] <- 1
# Create empty matrices to store model output
# Matrix to store numbers of individuals in each health state across cycles
pop <- array(data = NA,
dim = c(n_states, n_cycles+1, n_treat),
dimnames = list(state = states,
cycle = NULL,
treatment = treat))
# first cycle
pop["E", cycle = 1, ] <- n
pop["C", cycle = 1, ] <- 0
pop["D", cycle = 1, ] <- 0
# Matrix to store cost-effectiveness endpoints
costs <- LYs <- QALYs <- array(data = NA,
dim = c(n_treat, n_cycles+1),
dimnames = list(treatment = treat,
cycle = NULL))
# vector to store total life years (LYs), QALYs and costs
total_LYs <- total_QALYs <- total_costs <- setNames(c(NA, NA), treat)
### Run the model
# outer loop - loop over treatment groups
for (i in 1:n_treat) {
# inner loop - loop over cycles
for (j in 1:n_cycles) {
pop[, cycle = j+1, treatment = i] <- # starts in cycle 2 as cycle 1 is already defined above
pop[, cycle = j, treatment = i] %*% p_mat[, , treatment = i]
} # end of inner loop
# Life-years - 1 if person alive (i.e states E and C), 0 otherwise (state D)
LYs[i, ] <- (c(1,1,0) %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles) # cycle 1 = Year 0
# Quality-adjusted life-years (LYs adjusted by the quality of life)
QALYs[i, ] <- (u_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
# Costs
costs[i, ] <-(c_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
total_LYs[i] <- sum(LYs[treatment = i, -1]) # -1 means 'exclude Year 0'
total_QALYs[i] <- sum(QALYs[treatment = i, -1])
total_costs[i] <- sum(costs[treatment = i, -1])
} # end of outer loop
# Add one-off cost (per patient) of the intervention
total_costs<-c(total_costs["EVAR"] + c_evar*n, total_costs["OR"] + c_or*n)
# incremental LY
inc_ly <- total_LYs["EVAR"] - total_LYs["OR"]
# incremental QALY
inc_qaly <- total_QALYs["EVAR"] - total_QALYs["OR"]
# Incremental costs
inc_cost <- total_costs["EVAR"] - total_costs["OR"]
# Incremental net benefit at �30,000/QALY threshold
inb <- inc_qaly*30000-inc_cost
# Incremental cost-effectiveness ratio
icer <- inc_cost/inc_qaly
total_LYs
total_QALYs
total_costs
inc_ly
inc_qaly
inc_cost
inb
icer
#############################################################################
#
### Question 5
#
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 3-6")
#############################################################################
### Step 1: initialise model set-up as before until transitions matrix
# Time-dependent transition probabilities
# treatments
treat <- c("EVAR", "OR") # treatments
n_treat <- length(treat) # no. treatments
# health states (E: event-free survival, C: CVD event, D: Death)
states <- c("E", "C", "D")
n_states <- length(states) # no. states
# cohort
n <- 1000 # size of the cohort
age0 <- 55 # initial age of cohort
# cycles
n_cycles <- 43 # no. cycles
# Discount rate (same for costs and QALYs)
dr <- 0.035
# costs for each health state
c_E <- 500 # cost in E state
c_C <- 2500 # cost in C state
c_D <- 0 # cost in D state
c_evar <- 18000 # Cost of EVAR treatment (one-off)
c_or <- 15500 # Cost of Open Repair treatment (one-off)
# costs matrix across states and treatment groups
c_matrix <- matrix(
c(c_E, c_C, c_D,
c_E, c_C, c_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# utilities for each health state
u_E <- 0.80 # utility in E state
u_C <- 0.60 # utility in C state
u_D <- 0 # utility in D state
# utilities matrix across states and treatment groups
u_matrix <- matrix(
c(u_E, u_C, u_D,
u_E, u_C, u_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
#############################################################################
### Step 2: Write function to create transition probabilities that change with age cohort
# initialise matrix as before
p_mat <- array(data = NA,
dim = c(n_states, n_states, n_treat),
dimnames = list(from = states,
to = states,
treat))
# Create 'lookup' function to return a matrix that includes age-related transition probabilities
p_mat_cycle <- function(p_mat, age, cycle,
p_EC = 0.05,
p_CE = 0.22,
excess=0.10, # CHD-related excess mortality
effect = 0.1) {
# p_ED will vary according to age cohort 55-64, 65-74, 75-84, 85+
p_ED_lookup <-
c("(54,64]" = 0.008, # round bracket means 'exclusive'; square bracket means 'inclusive'
"(64,74]" = 0.019,
"(74,84]" = 0.056,
"(84,98]" = 0.203)
# in each cycle, p_ED value will be picked according to age cohort
age_group <- cut(age, breaks = c(54,64,74,84,98))
p_ED <- p_ED_lookup[age_group]
p_CD <- p_ED + excess
# Fill in the matrix (Open Repair)
# p_EC and p_CE will be constant across cycles
# From E
p_mat["E", "E", "OR"] <- 1 - p_EC - p_ED
p_mat["E", "C", "OR"] <- p_EC
p_mat["E", "D", "OR"] <- p_ED
# From C
p_mat["C", "E", "OR"] <- p_CE
p_mat["C", "C", "OR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "OR"] <- p_CD
# From D
p_mat["D", "E", "OR"] <- 0
p_mat["D", "C", "OR"] <- 0
p_mat["D", "D", "OR"] <- 1
# Fill in the matrix (EVAR)
# p_EC and p_CE will be constant across cycles
# From E
p_mat["E", "E", "EVAR"] <- 1 - p_EC*(1-effect) - p_ED
p_mat["E", "C", "EVAR"] <- p_EC*(1-effect)
p_mat["E", "D", "EVAR"] <- p_ED
# From C
p_mat["C", "E", "EVAR"] <- p_CE
p_mat["C", "C", "EVAR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "EVAR"] <- p_CD
# From D
p_mat["D", "E", "EVAR"] <- 0
p_mat["D", "C", "EVAR"] <- 0
p_mat["D", "D", "EVAR"] <- 1
return(p_mat)
} # end of lookup function
#############################################################################
### Step 3: Create empty matrices to store model outputs (as before)
# Matrix to store numbers of individuals in each health state across cycles
pop <- array(data = NA,
dim = c(n_states, n_cycles+1, n_treat),
dimnames = list(state = states,
cycle = NULL,
treatment = treat))
# first cycle
pop["E", cycle = 1, ] <- n
pop["C", cycle = 1, ] <- 0
pop["D", cycle = 1, ] <- 0
# Matrix to store cost-effectiveness endpoints
costs <- LYs <- QALYs <- array(data = NA,
dim = c(n_treat, n_cycles+1),
dimnames = list(treatment = treat,
cycle = NULL))
# vector to store total life years (LYs), QALYs and costs
total_LYs <- total_QALYs <- total_costs <- setNames(c(NA, NA), treat)
#############################################################################
### Step 4: Re-run the model by replacing p_mat by p_mat_cycle across each cycle
# outer loop - loop over treatment groups
for (i in 1:n_treat) {
# set the age 'clock' starting at age0
age<-age0
# inner loop - loop over cycles
for (j in 1:n_cycles) {
# in each cycle replace p_mat by p_mat_cycle that includes age-related transition probabilities
p_mat <- p_mat_cycle(p_mat, age, j)
pop[, cycle = j+1, treatment = i] <- # starts in cycle 2 as cycle 1 is already defined above
pop[, cycle = j, treatment = i] %*% p_mat[, , treatment = i]
# clock 'ticks' before next cycle
age <- age + 1
} # end of inner loop
# Life-years - 1 if person alive (i.e states E and C), 0 otherwise (state D)
LYs[i, ] <- (c(1,1,0) %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles) # cycle 1 = Year 0
# Quality-adjusted life-years (LYs adjusted by the quality of life)
QALYs[i, ] <- (u_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
# Costs
costs[i, ] <-(c_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
total_LYs[i] <- sum(LYs[treatment = i, -1]) # -1 means 'exclude Year 0'
total_QALYs[i] <- sum(QALYs[treatment = i, -1])
total_costs[i] <- sum(costs[treatment = i, -1])
} # end of outer loop
# Add one-off cost (per patient) of the intervention
total_costs<-c(total_costs["EVAR"] + c_evar*n, total_costs["OR"] + c_or*n)
# incremental LY
inc_ly <- total_LYs["EVAR"] - total_LYs["OR"]
# incremental QALY
inc_qaly <- total_QALYs["EVAR"] - total_QALYs["OR"]
# Incremental costs
inc_cost <- total_costs["EVAR"] - total_costs["OR"]
# Incremental net benefit at �30,000/QALY threshold
inb <- inc_qaly*30000-inc_cost
# Incremental cost-effectiveness ratio
icer <- inc_cost/inc_qaly
total_LYs
total_QALYs
total_costs
inc_ly
inc_qaly
inc_cost
inb
icer
#############################################################################
#
### Question 6
#
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 3-6")
#############################################################################
### Step 1: Create function that will be used to re-run the model for each parameter value
### chd_mod below is a function of all parameters which values we wish to vary in the
### deterministic sensitivity analysis (listed in Table 2).
### all other parameters are fixed (e.g. age0, dr, c_D, u_D, etc.)
chd_mod <- function(all_params) { # all_params - set of all parameters
with(as.list(all_params), {
# treatments
treat <- c("EVAR", "OR") # treatments
n_treat <- length(treat) # no. treatments
# health states (E: event-free survival, C: CVD event, D: Death)
states <- c("E", "C", "D")
n_states <- length(states) # no. states
# cohort
n <- 1000 # size of the cohort
age0 <- 55 # initial age of cohort
# cycles
n_cycles <- 43 # no. cycles
# Discount rate (same for costs and QALYs)
dr <- 0.035
# costs for each health state
c_E <- c_E # cost in E state
c_C <- c_C # cost in C state
c_D <- 0 # cost in D state
c_evar <- c_evar # Cost of EVAR treatment (one-off)
c_or <- c_or # Cost of Open Repair treatment (one-off)
# costs matrix across states and treatment groups
c_matrix <- matrix(
c(c_E, c_C, c_D,
c_E, c_C, c_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# utilities for each health state
u_E <- u_E # utility in E state
u_C <- u_C # utility in C state
u_D <- 0 # utility in D state
# utilities matrix across states and treatment groups
u_matrix <- matrix(
c(u_E, u_C, u_D,
u_E, u_C, u_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# annual transition probabilities between health states
p_EC <- p_EC # transition probability from E to C state
p_CE <- p_CE # transition probability from C to E state
p_ED <- p_ED # transition probability from E to D state
p_CD <- p_CD # transition probability from C to D state
effect<- effect # effect of treatment (reduces progression)
# 6.4.1 create matrix with all transition rates per treatment arm
# initialise transition matrix
p_mat <- array(data = NA,
dim = c(n_states, n_states, n_treat),
dimnames = list(from = states,
to = states,
treat))
# Fill in the matrix (Open Repair)
# From E
p_mat["E", "E", "OR"] <- 1 - p_EC - p_ED
p_mat["E", "C", "OR"] <- p_EC
p_mat["E", "D", "OR"] <- p_ED
# From C
p_mat["C", "E", "OR"] <- p_CE
p_mat["C", "C", "OR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "OR"] <- p_CD
# From D
p_mat["D", "E", "OR"] <- 0
p_mat["D", "C", "OR"] <- 0
p_mat["D", "D", "OR"] <- 1
# Fill in the matrix (EVAR)
# From E
p_mat["E", "E", "EVAR"] <- 1 - p_EC*(1-effect) - p_ED
p_mat["E", "C", "EVAR"] <- p_EC*(1-effect)
p_mat["E", "D", "EVAR"] <- p_ED
# From C
p_mat["C", "E", "EVAR"] <- p_CE
p_mat["C", "C", "EVAR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "EVAR"] <- p_CD
# From D
p_mat["D", "E", "EVAR"] <- 0
p_mat["D", "C", "EVAR"] <- 0
p_mat["D", "D", "EVAR"] <- 1
# Create empty matrices to store model output
# Matrix to store numbers of individuals in each health state across cycles
pop <- array(data = NA,
dim = c(n_states, n_cycles+1, n_treat),
dimnames = list(state = states,
cycle = NULL,
treatment = treat))
# first cycle
pop["E", cycle = 1, ] <- n
pop["C", cycle = 1, ] <- 0
pop["D", cycle = 1, ] <- 0
# Matrix to store cost-effectiveness endpoints
costs <- LYs <- QALYs <- array(data = NA,
dim = c(n_treat, n_cycles+1),
dimnames = list(treatment = treat,
cycle = NULL))
# vector to store total life years (LYs), QALYs and costs
total_LYs <- total_QALYs <- total_costs <- setNames(c(NA, NA), treat)
### Run the model
# outer loop - loop over treatment groups
for (i in 1:n_treat) {
# inner loop - loop over cycles
for (j in 1:n_cycles) {
pop[, cycle = j+1, treatment = i] <- # starts in cycle 2 as cycle 1 is already defined above
pop[, cycle = j, treatment = i] %*% p_mat[, , treatment = i]
} # end of inner loop
# Life-years - 1 if person alive (i.e states E and C), 0 otherwise (state D)
LYs[i, ] <- (c(1,1,0) %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles) # cycle 1 = Year 0
# Quality-adjusted life-years (LYs adjusted by the quality of life)
QALYs[i, ] <- (u_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
# Costs
costs[i, ] <-(c_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
total_LYs[i] <- sum(LYs[treatment = i, -1]) # -1 means 'exclude Year 0'
total_QALYs[i] <- sum(QALYs[treatment = i, -1])
total_costs[i] <- sum(costs[treatment = i, -1])
} # end of outer loop
# Add one-off cost (per patient) of the intervention
total_costs<-c(total_costs["EVAR"] + c_evar*n, total_costs["OR"] + c_or*n)
# incremental LY
inc_ly <- total_LYs["EVAR"] - total_LYs["OR"]
# incremental QALY
inc_qaly <- total_QALYs["EVAR"] - total_QALYs["OR"]
# Incremental costs
inc_cost <- total_costs["EVAR"] - total_costs["OR"]
# Incremental net benefit (INB), assuming a �30,000/QALY threshold
inb <- inc_qaly*30000-inc_cost
# Incremental cost-effectiveness ratio (ICER)
icer <- inc_cost/inc_qaly
# create dataframe to store all outputs of interest (cost, QALY, NMB, ICER) by treatment arm
ce <- data.frame(Strategy = treat,
Cost = numeric(n_treat),
QALY = numeric(n_treat),
INB = numeric(n_treat),
ICER = numeric(n_treat),
stringsAsFactors = FALSE)
ce[1, c("Cost", "QALY", "INB", "ICER")] <- c(total_costs["EVAR"], total_QALYs["EVAR"], inb, icer)
ce[2, c("Cost", "QALY", "INB", "ICER")] <- c(total_costs["OR"], total_QALYs["OR"], 0, 0)
return(ce)
})
} # end of CHD model
#############################################################################
### Step 2: Define range of values for parameters (listed in Table 2)
# parameter values used in the base case (reference point)
base <- list(c_evar=18000, c_or=15500, c_E=500, c_C=2500,u_E=0.8, u_C=0.6,
p_ED=0.04, p_CD=0.14, p_EC=0.05, p_CE=0.22, effect=0.1)
# parameter range for the DSA (in this case, lower and upper limits)
param_range <- data.frame(
pars = c("c_evar", "c_or", "c_E", "c_C","u_E", "u_C",
"p_ED", "p_CD", "p_EC", "p_CE", "effect"),
min = c(16000, 14000, 0, 1000, 0.7, 0.5, 0.02, 0.09, 0.03, 0.15, 0.07),
max = c(20000, 17000, 1000, 4000, 0.9, 0.7, 0.06, 0.2, 0.07, 0.3, 0.13))
#############################################################################
### Step 3: run DSA using 'dampack' library
library(dampack)
# 'run_owsa_det' is a function that essentially loops 'chd_mod' over each parameter value
dsa <- run_owsa_det(params_range = param_range, # range of parameters for DSA
params_basecase = base, # base-case parameter values
nsamp = 2, # 2 values (lower and upper limits)
FUN = chd_mod, # function created above
outcomes = c("Cost", "QALY", "INB", "ICER"), # outputs of interest
strategies = c("EVAR", "OR"), # treatment arms
progress = FALSE)
# DSA results are stored in 'dsa$owsa' (stacked by treatment arm)
dsa$owsa_Cost
dsa$owsa_QALY
dsa$owsa_INB
dsa$owsa_ICER
#############################################################################
### Step 4: Report the results
# there are many ways to do this.
# Option 1: report results in a table (less straightforward to interpret)
# Option 2: plot the results, often using a tornado diagram to help visualisation
# Here I will make use of the 'ggplot2' library but others could be used.
library(ggplot2)
library(plyr)
library(dplyr)
library(tidyverse)
# create dataframe with the outputs of interest (here, let's focus on ICER)
dsa.res <- data.frame(
pars = c("c_evar", "c_or", "c_E", "c_C","u_E", "u_C",
"p_ED", "p_CD", "p_EC", "p_CE", "effect"),
min = NA, # ICER results from lower limit values
max = NA, # ICER results from upper limit values
dif = NA) # difference between ICER of lower and upper limits (to sort results by magnitude of CI)
dsa.res[1, 2:4] <- c(dsa$owsa_ICER[1,"outcome_val"],dsa$owsa_ICER[2,"outcome_val"],abs(dsa$owsa_ICER[1,"outcome_val"]-dsa$owsa_ICER[2,"outcome_val"]))
dsa.res[2, 2:4] <- c(dsa$owsa_ICER[5,"outcome_val"],dsa$owsa_ICER[6,"outcome_val"],abs(dsa$owsa_ICER[5,"outcome_val"]-dsa$owsa_ICER[6,"outcome_val"]))
dsa.res[3, 2:4] <- c(dsa$owsa_ICER[9,"outcome_val"],dsa$owsa_ICER[10,"outcome_val"],abs(dsa$owsa_ICER[9,"outcome_val"]-dsa$owsa_ICER[10,"outcome_val"]))
dsa.res[4, 2:4] <- c(dsa$owsa_ICER[13,"outcome_val"],dsa$owsa_ICER[14,"outcome_val"],abs(dsa$owsa_ICER[13,"outcome_val"]-dsa$owsa_ICER[14,"outcome_val"]))
dsa.res[5, 2:4] <- c(dsa$owsa_ICER[17,"outcome_val"],dsa$owsa_ICER[18,"outcome_val"],abs(dsa$owsa_ICER[17,"outcome_val"]-dsa$owsa_ICER[18,"outcome_val"]))
dsa.res[6, 2:4] <- c(dsa$owsa_ICER[21,"outcome_val"],dsa$owsa_ICER[22,"outcome_val"],abs(dsa$owsa_ICER[21,"outcome_val"]-dsa$owsa_ICER[22,"outcome_val"]))
dsa.res[7, 2:4] <- c(dsa$owsa_ICER[25,"outcome_val"],dsa$owsa_ICER[26,"outcome_val"],abs(dsa$owsa_ICER[25,"outcome_val"]-dsa$owsa_ICER[26,"outcome_val"]))
dsa.res[8, 2:4] <- c(dsa$owsa_ICER[29,"outcome_val"],dsa$owsa_ICER[30,"outcome_val"],abs(dsa$owsa_ICER[29,"outcome_val"]-dsa$owsa_ICER[30,"outcome_val"]))
dsa.res[9, 2:4] <- c(dsa$owsa_ICER[33,"outcome_val"],dsa$owsa_ICER[34,"outcome_val"],abs(dsa$owsa_ICER[33,"outcome_val"]-dsa$owsa_ICER[34,"outcome_val"]))
dsa.res[10, 2:4] <- c(dsa$owsa_ICER[37,"outcome_val"],dsa$owsa_ICER[38,"outcome_val"],abs(dsa$owsa_ICER[37,"outcome_val"]-dsa$owsa_ICER[38,"outcome_val"]))
dsa.res[11, 2:4] <- c(dsa$owsa_ICER[41,"outcome_val"],dsa$owsa_ICER[42,"outcome_val"],abs(dsa$owsa_ICER[41,"outcome_val"]-dsa$owsa_ICER[42,"outcome_val"]))
# ICER from base-case
base.case <- 19374
# sort results (from largest to narrowest) of lower-upper ICER difference
order.pars <- dsa.res %>% arrange(dif) %>%
mutate(pars=factor(x=pars, levels=pars)) %>%
select(pars) %>% unlist() %>% levels()
# width of columns in plot (value between 0 and 1)
width <- 0.95
# get data frame in shape for ggplot and geom_rect
tornado <- dsa.res %>%
gather(key='type', value='output.value', min:max) %>% # gather columns min and max into a single column
select(pars, type, output.value, dif) %>% # reorder columns
mutate(pars=factor(pars, levels=order.pars), # create the columns for geom_rect
ymin=pmin(output.value, base.case),
ymax=pmax(output.value, base.case),
xmin=as.numeric(pars)-width/2,
xmax=as.numeric(pars)+width/2)
# create plot
png(width = 960, height = 540)
ggplot() + geom_rect(data = tornado, aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin, fill=type)) +
theme_bw() + theme(axis.title=element_text(size = 16), legend.text=element_text(size=16),
axis.text=element_text(size=16), legend.title = element_blank(), legend.position = 'bottom') +
ylab("ICER") + geom_hline(yintercept = base.case) +
scale_x_continuous(breaks = c(1:length(order.pars)), labels = order.pars) + coord_flip()
dev.off() # plot will be stored in your working directory
#############################################################################
#
### Question 7
#
#############################################################################
# two-way DSA will combine parameter values for 'c_evar' and 'c_or'
two_param_range <- data.frame(pars = c("c_evar", "c_or"),
min = c(16000, 14000),
max = c(20000, 17000))
# 'run_twsa_det' is similar to 'run_owsa_det'
two_dsa <- run_twsa_det(params_range = two_param_range, # range of parameters for two_DSA
params_basecase = base, # base-case parameter values
nsamp = 10, # 10 values (between lower and upper limits)
FUN = chd_mod, # same function created above
outcomes = c("Cost", "QALY", "INB", "ICER"), # outputs of interest
strategies = c("EVAR", "OR"), # treatment arms
progress = FALSE)
# DSA results are stored as previously
twsa_INB<-two_dsa$twsa_INB
# report results using two-way SA plot
plot(twsa_INB, col = "bw", txtsize=16)
#############################################################################
###
### GLBH0046: Advanced Economic Evaluation
###
### Practical 7: Probabilistic sensitivity analysis
###
#############################################################################
rm(list=ls())
# Set working directory and load the dataset
setwd("N:/Manuel Gomes/Teaching/Advanced EE/Practicals/Practicals 3-6")
set.seed(20910) # Set random number generator -> ensures simulations can be reproduced
#############################################################################
### Question 3
#############################################################################
# Plot each distribution in Table 1 and check plausibility
# EVAR costs
c_evar_dist<-rgamma(n=1000, shape=100, scale=180)
hist(c_evar_dist)
# Open Repair costs
c_or_dist<-rgamma(n=1000, shape=100, scale=155)
hist(c_or_dist)
# Costs of E state
c_E_dist<-rgamma(n=1000, shape=100, scale=5)
c_E_dist1<-rgamma(n=1000, shape=25, scale=20)
par(mfrow=c(1,2), cex.lab=1.5, cex.axis=1.5, cex.main=1.5)
hist(c_E_dist, main="SE=1/10 mean")
hist(c_E_dist1, main="SE=1/5 mean")
dev.off()
# Costs of C state
c_C_dist<-rgamma(n=1000, shape=100, scale=25)
dev.off() # reset par() function
hist(c_C_dist)
# Utility of E state
u_E_dist<-rbeta(n=1000, shape1=19.2, shape2=4.8)
hist(u_E_dist)
# Utility of C state
u_C_dist<-rbeta(n=1000, shape1=39.4, shape2=26.3)
hist(u_C_dist)
# transition probability p_ED
p_ED_dist<-rbeta(n=1000, shape1=15.3, shape2=367.7)
hist(p_ED_dist)
# transition probability p_CD
p_CD_dist<-rbeta(n=1000, shape1=42, shape2=258)
hist(p_CD_dist)
# transition probability p_EC
p_EC_dist<-rbeta(n=1000, shape1=23.7, shape2=450.3)
hist(p_EC_dist)
# transition probability p_EC
p_EC_dist<-rbeta(n=1000, shape1=23.7, shape2=450.3)
hist(p_EC_dist)
# Effect parameter
effect_dist<-rlnorm(n=1000, meanlog=log(0.1), sdlog=0.0686)
effect_dist1<-rlnorm(n=1000, meanlog=log(0.1), sdlog=0.1)
par(mfrow=c(1,2), cex.lab=1.5, cex.axis=1.5, cex.main=1.5)
hist(effect_dist, main="sdlog=0.068")
hist(effect_dist1, main="sdlog=0.1")
# transition probability p_CE
p_CE_dist<-rbeta(n=1000, shape1=94.2, shape2=333.8)
dev.off() # reset par() function
hist(p_CE_dist)
#############################################################################
### Question 4
#############################################################################
#############################################################################
### Step 1: Create function that will be used to re-run the model for each PSA value
### chd_mod below is a function of all parameters which values we wish to vary in the
### deterministic sensitivity analysis (listed in Table 2).
### all other parameters are fixed (e.g. age0, dr, c_D, u_D, etc.)
chd_mod <- function(c_evar, c_or, c_E, c_C,u_E, u_C,
p_ED, p_CD, p_EC, p_CE, effect, wtp) {
# treatments
treat <- c("EVAR", "OR") # treatments
n_treat <- length(treat) # no. treatments
# health states (E: event-free survival, C: CVD event, D: Death)
states <- c("E", "C", "D")
n_states <- length(states) # no. states
# cohort
n <- 1000 # size of the cohort
age0 <- 55 # initial age of cohort
# cycles
n_cycles <- 43 # no. cycles
# Discount rate (same for costs and QALYs)
dr <- 0.035
# costs for each health state
c_E <- c_E # cost in E state
c_C <- c_C # cost in C state
c_D <- 0 # cost in D state
c_evar <- c_evar # Cost of EVAR treatment (one-off)
c_or <- c_or # Cost of Open Repair treatment (one-off)
# costs matrix across states and treatment groups
c_matrix <- matrix(
c(c_E, c_C, c_D,
c_E, c_C, c_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# utilities for each health state
u_E <- u_E # utility in E state
u_C <- u_C # utility in C state
u_D <- 0 # utility in D state
# utilities matrix across states and treatment groups
u_matrix <- matrix(
c(u_E, u_C, u_D,
u_E, u_C, u_D),
byrow = TRUE,
nrow = n_treat,
dimnames = list(treat,states))
# annual transition probabilities between health states
p_EC <- p_EC # transition probability from E to C state
p_CE <- p_CE # transition probability from C to E state
p_ED <- p_ED # transition probability from E to D state
p_CD <- p_CD # transition probability from C to D state
effect<- effect # effect of treatment (reduces progression)
# 6.4.1 create matrix with all transition rates per treatment arm
# initialise transition matrix
p_mat <- array(data = NA,
dim = c(n_states, n_states, n_treat),
dimnames = list(from = states,
to = states,
treat))
# Fill in the matrix (Open Repair)
# From E
p_mat["E", "E", "OR"] <- 1 - p_EC - p_ED
p_mat["E", "C", "OR"] <- p_EC
p_mat["E", "D", "OR"] <- p_ED
# From C
p_mat["C", "E", "OR"] <- p_CE
p_mat["C", "C", "OR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "OR"] <- p_CD
# From D
p_mat["D", "E", "OR"] <- 0
p_mat["D", "C", "OR"] <- 0
p_mat["D", "D", "OR"] <- 1
# Fill in the matrix (EVAR)
# From E
p_mat["E", "E", "EVAR"] <- 1 - p_EC*(1-effect) - p_ED
p_mat["E", "C", "EVAR"] <- p_EC*(1-effect)
p_mat["E", "D", "EVAR"] <- p_ED
# From C
p_mat["C", "E", "EVAR"] <- p_CE
p_mat["C", "C", "EVAR"] <- 1 - p_CD - p_CE
p_mat["C", "D", "EVAR"] <- p_CD
# From D
p_mat["D", "E", "EVAR"] <- 0
p_mat["D", "C", "EVAR"] <- 0
p_mat["D", "D", "EVAR"] <- 1
# Create empty matrices to store model output
# Matrix to store numbers of individuals in each health state across cycles
pop <- array(data = NA,
dim = c(n_states, n_cycles+1, n_treat),
dimnames = list(state = states,
cycle = NULL,
treatment = treat))
# first cycle
pop["E", cycle = 1, ] <- n
pop["C", cycle = 1, ] <- 0
pop["D", cycle = 1, ] <- 0
# Matrix to store cost-effectiveness endpoints
costs <- LYs <- QALYs <- array(data = NA,
dim = c(n_treat, n_cycles+1),
dimnames = list(treatment = treat,
cycle = NULL))
# vector to store total life years (LYs), QALYs and costs
total_LYs <- total_QALYs <- total_costs <- setNames(c(NA, NA), treat)
### Run the model
# outer loop - loop over treatment groups
for (i in 1:n_treat) {
# inner loop - loop over cycles
for (j in 1:n_cycles) {
pop[, cycle = j+1, treatment = i] <- # starts in cycle 2 as cycle 1 is already defined above
pop[, cycle = j, treatment = i] %*% p_mat[, , treatment = i]
} # end of inner loop
# Life-years - 1 if person alive (i.e states E and C), 0 otherwise (state D)
LYs[i, ] <- (c(1,1,0) %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles) # cycle 1 = Year 0
# Quality-adjusted life-years (LYs adjusted by the quality of life)
QALYs[i, ] <- (u_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
# Costs
costs[i, ] <-(c_matrix[treatment = i, ] %*% pop[, , treatment = i]) * 1/(1 + dr)^(0:n_cycles)
total_LYs[i] <- sum(LYs[treatment = i, -1]) # -1 means 'exclude Year 0'
total_QALYs[i] <- sum(QALYs[treatment = i, -1])
total_costs[i] <- sum(costs[treatment = i, -1])
} # end of outer loop
# Add one-off cost (per patient) of the intervention
total_costs<-c(total_costs["EVAR"] + c_evar*n, total_costs["OR"] + c_or*n)
# incremental LY
inc_ly <- total_LYs["EVAR"] - total_LYs["OR"]
# incremental QALY
inc_qaly <- total_QALYs["EVAR"] - total_QALYs["OR"]
# Incremental costs
inc_cost <- total_costs["EVAR"] - total_costs["OR"]
# net monetary benefit at �30,000 per QALY - will be needed for value of information analysis (VOI)
nb.evar <- total_QALYs["EVAR"]*wtp - total_costs["EVAR"]
nb.or <- total_QALYs["OR"]*wtp - total_costs["OR"]
# Incremental net benefit (INB)
inb <- inc_qaly*wtp-inc_cost
# Cost-effective: =1 if INB>0, 0 otherwise (could also be calculated as nb.evar-nb.or)
pce <- ifelse(inb>0,1,0)
# Incremental cost-effectiveness ratio (ICER)
icer <- inc_cost/inc_qaly
ce <- c(inc_qaly,inc_cost,nb.evar,nb.or,pce)
return(ce)
} # end of CHD model
#############################################################################
### Step 2: Create an empty matrix to store all the outputs of interest
M<-1000 # define number of PSA runs
# Create an empty matrix to store outputs of interest
psa <- matrix(NA, nrow = M, ncol = 5,
dimnames = list(NULL, c("inc_qaly","inc_cost", "nb.evar", "nb.or", "pce")))
#############################################################################
### Step 3: Re-run chd_mod M=1000 times
for (i in 1:M) {
psa.res <- chd_mod(
c_evar=c_evar_dist[i],
c_or =c_or_dist[i],
c_E =c_E_dist1[i],
c_C =c_C_dist[i],
u_E =u_E_dist[i],
u_C =u_C_dist[i],
p_ED =p_ED_dist[i],
p_CD =p_CD_dist[i],
p_EC =p_EC_dist[i],
p_CE =p_CE_dist[i],
effect=effect_dist1[i],
wtp =30000)
psa [i,] <- psa.res
}
#############################################################################
### Step 4: Plot results in the cost-effectiveness (CE) plane.
n<-1000 # cohort size
wtp<-30000 # willingness-to-pay threshold for QALY gain
base_case<-c(121, 2351145) # incremental QALY and incremental cost from base-case
# plot results in the CE plane
par(las=1)
plot(x = psa[,1]/n, y = psa[,2]/n,
main=" ", xlab="", ylab="",
xlim = c(0, 0.4), ylim = c(-5000, 10000),
pch=19, cex=1.2, col="darkslategray", axes=F )
points(x=base_case[1]/n, y=base_case[2]/n, col='red', cex=1.5, pch=16)
axis(1, pos=0, cex.axis=1.5)
axis(2, pos=0, cex.axis=1.5)
abline(a = 0, b = wtp, lty='dotted', lwd = 2) # Willingness-to-pay threshold
text(0.35, 7500, "WTP - �30,000", cex=1.5)
mtext(side=1, line=1, "Incremental QALY", cex=1.5)
mtext(side=2, line=3, "Incremental cost", cex=1.5, las=0)
#############################################################################
### Question 5
#############################################################################
#############################################################################
### Report cost-effectiveness acceptability curve (CEAC)
### Plot CEAC
wtp <- seq(0,100000, length=21) # Willingness-to-pay (WTP) values for QALY gain
n_wtp <- length(wtp) # no. Of WTP
ceac<-c()
for (j in 1:n_wtp) { # outer loop - across alternative WTP values
for (i in 1:M) { # inner loop - across PSA values
psa.res <- chd_mod(
c_evar=c_evar_dist[i],
c_or =c_or_dist[i],
c_E =c_E_dist1[i],
c_C =c_C_dist[i],
u_E =u_E_dist[i],
u_C =u_C_dist[i],
p_ED =p_ED_dist[i],
p_CD =p_CD_dist[i],
p_EC =p_EC_dist[i],
p_CE =p_CE_dist[i],
effect=effect_dist1[i],
wtp = wtp[j])
psa [i,] <- psa.res
}
ceac[j] <- sum(psa[,"pce"])/M # % PSA runs for which INB>0
}
plot(wtp, ceac, type="l", ylim=c(0,1), main=" ", col = "dark red", cex.lab = 1.5,
xaxt='n', yaxt='n',
ylab="Prob EVAR being cost-effective", xlab="Willingness to pay for a QALY (�, thousands)" )
axis(1,at=c(0,10000,20000,30000,40000,50000,60000,70000,80000,90000,100000), labels=c(0,10,20,30,40,50,60,70,80,90,100),cex.axis=1.5)
axis(2,at=seq(0.1,.9,0.2), labels=seq(0.1,0.9,0.2), cex.axis=1.5, las=1)
abline(v=20000, lty=3, col="grey")
abline(v=30000, lty=3, col="grey")
#############################################################################
### Question 6
#############################################################################
#############################################################################
### Calculate expected value of perfect information (EVPI)
### Plot EVPI
evpi <-c()
for (j in 1:n_wtp) { # outer loop - across alternative WTP values
for (i in 1:M) { # inner loop - across PSA values
psa.res <- chd_mod(
c_evar=c_evar_dist[i],
c_or =c_or_dist[i],
c_E =c_E_dist1[i],
c_C =c_C_dist[i],
u_E =u_E_dist[i],
u_C =u_C_dist[i],
p_ED =p_ED_dist[i],
p_CD =p_CD_dist[i],
p_EC =p_EC_dist[i],
p_CE =p_CE_dist[i],
effect=effect_dist1[i],
wtp = wtp[j])
psa [i,] <- psa.res
}
nb <- psa[,3:4] # Net benefit matrix
# EVPI=E(max NB)- max E(NB)
# E(max NB) - average of the maximum benefit between EVAR or OR across M runs
# max E(NB) - maximum of the average net benefit between EVAR and OR
evpi[j] <- mean(apply(nb, 1, max)) - max(colMeans(nb))
}
par(mar = c(5,7,4,2) + 0.1) # increase axis space
par(mgp=c(3.5,1,0)) # adjust space between label and axis
plot(wtp, evpi, type="l", main=" ", col = "dark red", xaxt='n', yaxt='n', cex.lab = 1.5,
ylab="Population EVPI (�, thousands)", xlab="Willingness to pay for a QALY (�, thousands)" )
axis(1,at=c(0,10000,20000,30000,40000,50000,60000,70000,80000,90000,100000), labels=c(0,10,20,30,40,50,60,70,80,90,100), cex.axis=1.5)
axis(2,at=c(0,1e+5,2e+5,3e+5,4e+5,5e+5,6e+5,7e+5,8e+5,9e+5,1e+6),labels=c(0,100,200,300,400,500,600,700,800,900,1000), cex.axis=1.5, las=1)
abline(h=500000, lty=3, col="grey")