This code forms the basis for the state-transition model of the tutorial: ‘A Tutorial on Time-Dependent Cohort State-Transition Models in R using a Cost-Effectiveness Analysis Example’
Authors: - Fernando Alarid-Escudero - Eline Krijkamp - Eva A. Enns - Alan Yang - M.G. Myriam Hunink - Petros Pechlivanoglou - Hawre Jalal Please cite the article when using this code
To program this tutorial we used: R version 4.0.5 (2021-03-31) Platform: 64-bit operating system, x64-based processor Running under: Mac OS 12.2.1 RStudio: Version 1.4.1717 2009-2021 RStudio, Inc
This code implements a simulation-time-dependent Sick-Sicker cSTM
model to conduct a CEA of two strategies: - Standard of Care (SoC): best
available care for the patients with the disease. This scenario reflects
the natural history of the disease
progression. - Strategy AB: This strategy combines treatment A and
treatment B. The disease progression is reduced, and individuals in the
Sick state have an
improved quality of life.
Change eval to TRUE if you want to knit
this document.
rm(list = ls()) # clear memory (removes all the variables from the workspace)if (!require('pacman')) install.packages('pacman'); library(pacman) # use this package to conveniently install other packages
# load (install if required) packages from CRAN
p_load("dplyr", "tidyr", "reshape2", "devtools", "scales", "ellipse", "ggplot2", "ggrepel", "gridExtra", "lazyeval", "igraph", "truncnorm", "ggraph", "reshape2", "patchwork", "knitr", "stringr", "diagram", "dampack")
# load (install if required) packages from GitHub
# install_github("DARTH-git/darthtools", force = TRUE) #Uncomment if there is a newer version
p_load_gh("DARTH-git/darthtools")# all functions are in the darthtools package## General setup
cycle_length <- 1 # cycle length equal to one year (use 1/12 for monthly)
n_age_init <- 25 # age at baseline
n_age_max <- 100 # maximum age of follow up
n_cycles <- (n_age_max - n_age_init)/cycle_length # time horizon, number of cycles
# Age labels
v_age_names <- paste(rep(n_age_init:(n_age_max-1), each = 1/cycle_length),
1:(1/cycle_length),
sep = ".")
# the 4 health states of the model:
v_names_states <- c("H", # Healthy (H)
"S1", # Sick (S1)
"S2", # Sicker (S2)
"D") # Dead (D)
n_states <- length(v_names_states) # number of health states
### Discounting factors
# Comment out/ for the relevant country
# USA
#d_c <- d_e <- 0.03 # annual discount rate for costs and QALY
# Canada
#d_c <- d_e <- 0.015 # annual discount rate for costs
# NL
d_c <- 0.04 # annual discount rate for costs
d_e <- 0.015 # annual discount rate for QALYs
### Strategies
v_names_str <- c("Standard of care", # store the strategy names
"Strategy AB")
n_str <- length(v_names_str) # number of strategies
## Within-cycle correction (WCC) using Simpson's 1/3 rule
v_wcc <- gen_wcc(n_cycles = n_cycles, method = "Simpson1/3")
### Transition rates (annual), and hazard ratios (HRs)
r_HS1 <- 0.15 # constant annual rate of becoming Sick when Healthy
r_S1H <- 0.5 # constant annual rate of becoming Healthy when Sick
r_S1S2 <- 0.105 # constant annual rate of becoming Sicker when Sick
hr_S1 <- 3 # hazard ratio of death in Sick vs Healthy
hr_S2 <- 10 # hazard ratio of death in Sicker vs Healthy
### Effectiveness of treatment AB
hr_S1S2_trtAB <- 0.6 # hazard ratio of becoming Sicker when Sick under treatment AB
## Age-dependent mortality rates
lt_usa_2015 <- read.csv("../data/HMD_USA_Mx_2015.csv")
# Extract age-specific all-cause mortality for ages in model time horizon
v_r_mort_by_age <- lt_usa_2015 %>%
dplyr::filter(Age >= n_age_init & Age < n_age_max) %>%
dplyr::select(Total) %>%
as.matrix()
### State rewards
#### Costs
c_H <- 2000 # annual cost of being Healthy
c_S1 <- 4000 # annual cost of being Sick
c_S2 <- 15000 # annual cost of being Sicker
c_D <- 0 # annual cost of being dead
c_trtAB <- 25000 # annual cost of receiving treatment AB
#### Utilities
u_H <- 1 # annual utility of being Healthy
u_S1 <- 0.75 # annual utility of being Sick
u_S2 <- 0.5 # annual utility of being Sicker
u_D <- 0 # annual utility of being dead
u_trtAB <- 0.95 # annual utility when receiving treatment AB
### Transition rewards
du_HS1 <- 0.01 # disutility when transitioning from Healthy to Sick
ic_HS1 <- 1000 # increase in cost when transitioning from Healthy to Sick
ic_D <- 2000 # increase in cost when dying
### Discount weight for costs and effects
v_dwc <- 1 / ((1 + (d_e * cycle_length)) ^ (0:n_cycles))
v_dwe <- 1 / ((1 + (d_c * cycle_length)) ^ (0:n_cycles))
# Process model inputs
## Age-specific transition rates to the Dead state for all cycles
v_r_HDage <- rep(v_r_mort_by_age, each = 1/cycle_length)
# Name age-specific mortality vector
names(v_r_HDage) <- v_age_names
# compute mortality rates
v_r_S1Dage <- v_r_HDage * hr_S1 # Age-specific mortality rate in the Sick state
v_r_S2Dage <- v_r_HDage * hr_S2 # Age-specific mortality rate in the Sicker state
# transform rates to probabilities adjusting by cycle length
p_HS1 <- rate_to_prob(r = r_HS1, t = cycle_length) # constant annual probability of becoming Sick when Healthy conditional on surviving
p_S1H <- rate_to_prob(r = r_S1H, t = cycle_length) # constant annual probability of becoming Healthy when Sick conditional on surviving
p_S1S2 <- rate_to_prob(r = r_S1S2, t = cycle_length) # constant annual probability of becoming Sicker when Sick conditional on surviving
v_p_HDage <- rate_to_prob(v_r_HDage, t = cycle_length) # Age-specific mortality risk in the Healthy state
v_p_S1Dage <- rate_to_prob(v_r_S1Dage, t = cycle_length) # Age-specific mortality risk in the Sick state
v_p_S2Dage <- rate_to_prob(v_r_S2Dage, t = cycle_length) # Age-specific mortality risk in the Sicker state
## Annual transition probability of becoming Sicker when Sick for treatment AB
# Apply hazard ratio to rate to obtain transition rate of becoming Sicker when Sick for treatment AB
r_S1S2_trtAB <- r_S1S2 * hr_S1S2_trtAB
# Transform rate to probability to become Sicker when Sick under treatment AB
# adjusting by cycle length conditional on surviving
p_S1S2_trtAB <- rate_to_prob(r = r_S1S2_trtAB, t = cycle_length)# All starting healthy
v_m_init <- c(H = 1, S1 = 0, S2 = 0, D = 0) # initial state vector
v_m_init### Initialize cohort trace under SoC
m_M_SoC <- matrix(NA,
nrow = (n_cycles + 1), ncol = n_states,
dimnames = list(0:n_cycles, v_names_states))
# Store the initial state vector in the first row of the cohort trace
m_M_SoC[1, ] <- v_m_init
### Initialize cohort trace for strategy AB
# Structure and initial states are the same as for SoC
m_M_strAB <- m_M_SoC # Strategy AB## Create transition probability arrays for strategy SoC
### Initialize transition probability array for strategy SoC
# All transitions to a non-death state are assumed to be conditional on survival
a_P_SoC <- array(0,
dim = c(n_states, n_states, n_cycles),
dimnames = list(v_names_states,
v_names_states,
0:(n_cycles - 1)))
### Fill in array
## From H
# your turn
a_P_SoC["H", "H", ] <- (1 - v_p_HDage) * (1 - p_HS1)
a_P_SoC["H", "S1", ] <- (1 - v_p_HDage) * p_HS1
a_P_SoC["H", "D", ] <- v_p_HDage
## From S1
# your turn
a_P_SoC["S1", "S1", ] <- (1 - p_S1S2) * (1 -v_p_S1Dage)
a_P_SoC["S1", "S2", ] <-(1 -v_p_S1Dage)* p_S1S2
a_P_SoC["S1", "D", ] <- v_p_S1Dage
## From S2
# your turn
a_P_SoC["S2", "S2", ] <-(1 -v_p_S2Dage)
a_P_SoC["S2", "D", ] <- v_p_S2Dage
## From D
# your turn
a_P_SoC["D", "D", ] <- 1
### Initialize transition probability array for strategy AB
# your turn
a_P_strAB <- a_P_SoC
# Update only transition probabilities from S1 involving p_S1S2
# your turn
a_P_strAB["S1", "S1", ] <- (1 - v_p_S1Dage) * (1 - p_S1S2_trtAB)
a_P_strAB["S1", "S2", ] <- (1 - v_p_S1Dage) * p_S1S2_trtAB
## Check if transition probability arrays are valid
### Check that transition probabilities are [0, 1]
check_transition_probability(a_P_SoC, verbose = TRUE)
check_transition_probability(a_P_strAB, verbose = TRUE)
### Check that all rows for each slice of the array sum to 1
check_sum_of_transition_array(a_P_SoC, n_states = n_states, n_cycles = n_cycles, verbose = TRUE)
check_sum_of_transition_array(a_P_strAB, n_states = n_states, n_cycles = n_cycles, verbose = TRUE)# These arrays will capture transitions from each state to another over time
### Initialize transition dynamics array for strategy SoC
# your turn
# Set first slice of a_A_SoC with the initial state vector in its diagonal
# your turn
### Initialize transition-dynamics array for strategy AB
# Structure and initial states are the same as for SoC
# your turn# Iterative solution of age-dependent cSTM
# your turn
for(t in 1:n_cycles){
# For SoC
m_M_SoC[t + 1, ] <- m_M_SoC[t, ] %*% a_P_SoC[, , t]
# For strategy AB
m_M_strAB[t + 1, ] <- m_M_strAB[t, ] %*% a_P_strAB[, , t]}
## Store the cohort traces in a list
# your turn
l_m_M <- list(SoC = m_M_SoC,
AB = m_M_strAB)
names(l_m_M) <- v_names_str
## Store the transition dynamics array for each strategy in a list
# your turn# your turn
plot_trace(m_M_SoC) ## Scale by the cycle length
# Vector of state utilities under strategy SoC
# your turn
v_u_SoC <- c(H = u_H,
S1 = u_S1,
S2 = u_S2,
D = u_D) * cycle_length
# Vector of state costs under strategy SoC
# your turn
v_c_SoC <- c(H = c_H,
S1 = c_S1,
S2 = c_S2,
D = c_D) * cycle_length
# Vector of state utilities under strategy AB
# your turn
v_u_strAB <- c(H = u_H,
S1 = u_S1,
S2 = u_S2,
D = u_D) * cycle_length
# Vector of state costs under strategy AB
# your turn
v_c_strAB <- c(H = c_H + c_trtAB,
S1 = c_S1,
S2 = c_S2,
D = c_D) * cycle_length
## Store state rewards
# Store the vectors of state utilities for each strategy in a list
# your turn
l_u <- list(SQ = v_u_SoC,
AB = v_u_strAB)
# Store the vectors of state cost for each strategy in a list
# your turn
l_c <- list(SQ = v_c_SoC,
AB = v_c_strAB)
# assign strategy names to matching items in the lists
# your turn
names(l_u) <- names(l_c) <- v_names_str# Create empty vectors to store total utilities and costs
# your turn
v_tot_qaly <- v_tot_cost <- vector(mode = "numeric", length = n_str)
names(v_tot_qaly) <- names(v_tot_cost) <- v_names_str
## Loop through each strategy and calculate total utilities and costs
# your turn
for (i in 1:n_str) {
v_u_str <- l_u[[i]] # select the vector of state utilities for the i-th strategy
v_c_str <- l_c[[i]] # select the vector of state costs for the i-th strategy
### Expected QALYs and costs per cycle
## Vector of QALYs and Costs
# Apply state rewards
v_qaly_str <- l_m_M[[i]] %*% v_u_str # sum the utilities of all states for each cycle
v_cost_str <- l_m_M[[i]] %*% v_c_str # sum the costs of all states for each cycle
### Discounted total expected QALYs and Costs per strategy and apply within-cycle correction if applicable
# QALYs
v_tot_qaly[i] <- t(v_qaly_str) %*% (v_dwe * v_wcc)
# Costs
v_tot_cost[i] <- t(v_cost_str) %*% (v_dwc * v_wcc)}## Incremental cost-effectiveness ratios (ICERs)
# your turn
df_cea <- calculate_icers(cost = v_tot_cost,
effect = v_tot_qaly,
strategies = v_names_str)
df_cea## CEA table in proper format
# your turn
table_cea <- format_table_cea(df_cea)
table_cea## CEA frontier
# your turn
plot(df_cea, label = "all", txtsize = 14) +
expand_limits(x = max(table_cea$QALYs) + 0.1) +
theme(legend.position = c(0.8, 0.3))