Time to Code
Introduction to the Breast Cancer Modeling Problem
We will develop a 4-state Markov model to perform an analysis of the expected outcomes for a cohort of 55-year-old women who have undergone tumor excision for localized breast cancer. After surgery, all patients are initially tumor-free (i.e., will start in the “Local” Markov state). Each year, patients in the Local state face a 2% chance that they will experience a recurrence. Of these recurrences, 75% will present as metastatic disease (and patients will enter the “Mets” Markov state), and the remainder present as local disease (and patients will enter the “Recur” Markov state). Patients in the Recur state undergo repeat resection. Following this operation, they face a 6% chance of recurrence each year (assume that all local recurrences - 1st, 2nd, 3rd, etc. - have equivalent prognoses). Of these recurrences, 90% will present as metastatic disease. If the patients do not recur after the initial local recurrence, they remain in the Recur state since the prognosis remains the same each year. Only patients with metastatic disease can die of breast cancer. All women face a probability of dying from other causes based on the U.S. life tables. Women who die (either of breast cancer or other causes) will enter the “Dead” Markov state.
The Markov model described above represents usual care for women with breast cancer. Once you have constructed that Markov model you will then add a second strategy of a hypothetical new treatment. The effect of this new treatment is to reduce the initial risk of cancer recurrence by 50%; however, the treatment also lowers the quality of life of women taking the medication because of side effects. Treatment is only taken by women who have not experienced a recurrence (i.e., those who reside in the Local state).
After comparing the two strategies based on quality-adjusted life expectancy, you will incorporate costs and calculate the incremental cost per quality-adjusted life year gained of treatment compared with usual care.
Utilities
In a standard decision tree, the utility values are assigned at the end of the tree at each of the terminal nodes. In a Markov model, the terminal node specifies what state to go to for the next cycle, rather than the payoff.
There are three types of utilities in a Markov model:
The initial utility (Init) of a state is the utility a “person” accrues because of being in a state when the model starts (applied only in the first cycle).
The incremental utility (Incr) of a state is the utility a “person” accrues each time she or he arrives in that state as the model runs (applied many times).
The final utility (Final) of a state is the utility a “person” accrues if she or he is in that state at the time the model finishes running (applied only in the final cycle).
You will define these utilities using utility variables (uLocal
, uRecur
, and uMets
) in order to be able to calculate quality-adjusted life expectancy, and use a discount variable (r
) that will allow you to discount life expectancy (and quality-adjusted life expectancy). The variable “uLocal
” is the utility for being in the Local state. The variable “uRecur
” is the utility for being in the Recur state. The variable “uMets
” is the utility for being in the Mets state. If these utility values are equal to 1, then the expected value represents life expectancy (the model gives “credit” of 1 full year for spending a year in any “non-dead” state). We will use these to calculate quality-adjusted life expectancy.
Branch | Initial* | Utility Assignments Incremental | Final |
---|---|---|---|
Local | 0.5*uLocal | uLocal/(1+r)^_stage | 0 |
Recur | 0 | uRecur/(1+r)^_stage | 0 |
Mets | 0 | uMets/(1+r)^_stage | 0 |
Dead | 0 | 0 | 0 |
* First year gets credit for half a year - known as the “half-cycle correction.”
Tasks
Setting up the R Scripts:
- Create an R script named
utils.R
. This script will house all the functions we will develop for our analysis. - Create another R script named
CEA.R
(Cost-Effectiveness Analysis). This will be the main script where we will perform our analysis.
- Create an R script named
Initializing the
CEA.R
Script: Define the following variables:Age
,muMets
,pBCA2
,pMets1
,pMets2
,r
,uMets
, anduRecur.
Variable | Value |
---|---|
uRecur | 0.80 |
uMets | 0.40 |
r | 0.03 |
pMets2 | 0.90 |
pMets1 | 0.75 |
pBCA2 | 0.06 |
muMets | 0.50 |
Age | 55 |
These are the variables that you will define globally. In the CEA.R
script, start by copying and pasting the following code block:
##############################################################################
# Author: [Your Name]
# Purpose: Main analysis file for breast cancer treatment simulation
##############################################################################
# Libraries used
library(readxl) # used to hazard tables in from Excel
library(ggplot2) # used to graph results
# Functions for our analysis
source('code/utils.R')
# Mortality tables
<- read_excel("data/mufem.xlsx", sheet = 'FEMALES')
mufem
##############################################################################
# Run the simulation
##############################################################################
# Global variables
<- 0.80
uRecur <- 0.40
uMets <- 0
uDead <- 0.03
r <- 0.90
pMets2 <- 0.75
pMets1 <- 0.06
pBCA2 <- 0.50
muMets <- 55 Age
Save and Commit:
- Save both scripts in your project directory.
- Commit these initial scripts to your Git repository, ensuring that you’re tracking the changes and setting a strong foundation for version control.
Writing some helper functions:
The function for pDie (shown below) converts the annual rates from the life table (in mufem) to annual probabilities. The first cycle (year 1) through the Markov model, stage is equal to 0. Hence, the rate selected from the table is that associated with the starting age of the cohort (55). The second cycle through the model, stage is equal to 1 and the annual rate selected from the mufem table is that associated with 56-year-old women.
<- function(stage) {
pDie <- Age + stage
age_index <- mufem$`annual hazard rate`[age_index]
rate return(1 - exp(-rate))
}
The overall annual mortality rate for women with metastatic disease is equal to the sum of the rates of dying from other causes and the excess disease-specific mortality rate (defined as muMets
). The function for pDieMets
(shown below) converts the sum of these rates to annual probabilities.
<- function(stage) {
pDieMets <- Age + stage
age_index <- mufem$`annual hazard rate`[age_index] + muMets
rate return(1 - exp(-rate))
}
Lets add both of these functions to our utils.R
script.
- Writing main simulation function and defining local variables:
You should have 2 undefined variables now: pBCA1
and uLocal
. Since the value of pBCA1
and uLocal
will be different depending on the strategy, you will define these two variables locally within our simulation function.
“Usual Care”:
uLocal
= 0.95 andpBCA1
= 0.02.“Treat”:
uLocal
= 0.90 andpBCA1
= 0.01.
Add the following function to your utils.R
script.
# Define the Markov simulation function
<- function(initial_population=1, cycles=45, strategy='Usual Care') {
simulate_markov
if(strategy == 'Usual Care'){
<- 0.02
pBCA1 <- 0.95
uLocal
else{
} <- 0.01
pBCA1 <- 0.90
uLocal
}
<- c(uLocal, uRecur, uMets, uDead)
utilities
<- data.frame(Cycle = 1:cycles,
results Local = rep(0, cycles),
Recur = rep(0, cycles),
Mets = rep(0, cycles),
Dead = rep(0, cycles),
cumulative.utility = rep(0, cycles))
<- c(Local = initial_population, Recur = 0, Mets = 0, Dead = 0)
current_population
for (cycle in 1:cycles) {
if(cycle == 1){ #half cycle correction
<- sum(current_population*utilities) * 0.5
utility_value
}
# Transition logic
<- current_population["Local"] * pDie(cycle)
local_death <- current_population["Recur"] * pDie(cycle)
recur_death
<- current_population["Local"] * pBCA1
local_recur <- local_recur * pMets1
local_mets <- local_recur * (1 - pMets1)
local_recur
<- current_population["Recur"] * pBCA2
recur_recur <- recur_recur * pMets2
recur_mets <- recur_recur * (1 - pMets2)
recur_recur
<- current_population["Mets"] * pDieMets(cycle)
mets_death
# Update populations
"Local"] <- current_population["Local"] - local_recur - local_mets - local_death
current_population["Recur"] <- current_population["Recur"] + local_recur - recur_recur - recur_mets - recur_death
current_population["Mets"] <- current_population["Mets"] + local_mets + recur_mets - mets_death
current_population["Dead"] <- current_population["Dead"] + mets_death + local_death + recur_death
current_population[
# Update QALY
<- utility_value + (sum(current_population*utilities)/((1 + r)^cycle))
utility_value
# Store results for this cycle
2:6] <- c(current_population, utility_value)
results[cycle,
}
return(results)
}
- Run the model and save some outputs:
Great now we have everything we need to run this first analysis! Lets add the following code under “Run the simulation” in the CEA.R
script.
##############################################################################
# Run the simulation
##############################################################################
# Global variables
<- 0.80
uRecur <- 0.40
uMets <- 0
uDead <- 0.03
r <- 0.90
pMets2 <- 0.75
pMets1 <- 0.06
pBCA2 <- 0.50
muMets <- 55
Age
<- simulate_markov(strategy = 'Usual Care')
results_usual <- simulate_markov(strategy = 'Treatment') results_treat
If this code runs we are going to add one more section:
##############################################################################
# Visualization
##############################################################################
# Plot the results for both groups
<- ggplot(results_usual, aes(x = Cycle)) +
plot_usual geom_line(aes(y = Local, color = "Local")) +
geom_line(aes(y = Recur, color = "Recur")) +
geom_line(aes(y = Mets, color = "Mets")) +
geom_line(aes(y = Dead, color = "Dead")) +
scale_y_continuous(labels = scales::percent) +
labs(title = "Markov Model Simulation for Usual Care",
y = "Number of Individuals") +
theme_minimal()
ggsave('output/figures/Usual Care Progression.png',
width = 7, height = 4, bg='white')
plot_usual,
<- ggplot(results_treat, aes(x = Cycle)) +
plot_treatment geom_line(aes(y = Local, color = "Local")) +
geom_line(aes(y = Recur, color = "Recur")) +
geom_line(aes(y = Mets, color = "Mets")) +
geom_line(aes(y = Dead, color = "Dead")) +
scale_y_continuous(labels = scales::percent) +
labs(title = "Markov Model Simulation for Treatment",
y = "Number of Individuals") +
theme_minimal()
ggsave('output/figures/Treatment Progression.png',
width = 7, height = 4, bg='white') plot_treatment,
This code generates two figures showing the disease progression of the cohort and exports them straight to our output/figures
folder. This means that everytime we run this script these figures are getting generated. This is good because we want code that can generate all results and figures from scratch everytime!
If you have not commited and pushed code in a while you should do so!!!!
- Let’s add in costs!:
Cost-Effectiveness Analysis
You will now add costs as a 2nd Markov reward for each of the 4 states. Any differences between the two strategies in terms of costs will have to be locally assigned. The cost assignments are given in the following table. You will need to create variables when you enter unrecognized variables.
Cost Assignments
Branch | Initial* | Incremental | Final |
---|---|---|---|
Local | 0.5*(cLocal+cRx) | (cLocal+cRx)/(1+r)^_stage | 0 |
Recur | 0 | cRecur/(1+r)^_stage | 0 |
Mets | 0 | cMets/(1+r)^_stage | 0 |
Dead | 0 | 0 | 0 |
* First year gets credit for half a year - known as the “half-cycle correction.”
The annual costs of breast cancer (cLocal, cRecur, and cMets) will be the same for both strategies
Variable | Value |
---|---|
cRecur | 5000 |
cMets | 20000 |
cLocal | 500 |
Therefore, these variables can be assigned globally. So let’s go ahead and add the following to CEA.R
:
<- 5000
cRecur <- 20000
cMets <- 500 cLocal
To locally assign a cost for cRx (the annual cost of treatment), we will do so in our markov model function we wrote earlier. Here we can either decided to modify the one we had previously written, or we could write a new function. I am going to write a new function and add it utils.R
. We will need to locally assign the cost of treatment (cRx
), which is 0 if “Usual Care” and 1000 if “Treatment”
# Define the Markov simulation function with costs
<- function(initial_population=1, cycles=45, strategy='Usual Care') {
simulate_markov_costs
if(strategy == 'Usual Care'){
<- 0.02
pBCA1 <- 0.95
uLocal <- 0
cRx
else{
} <- 0.01
pBCA1 <- 0.90
uLocal <- 1000
cRx
}
<- c(uLocal, uRecur, uMets, uDead)
utilities <- c(cLocal + cRx, cRecur, cMets, 0)
costs
<- data.frame(Cycle = 1:cycles,
results Local = rep(0, cycles),
Recur = rep(0, cycles),
Mets = rep(0, cycles),
Dead = rep(0, cycles),
cumulative.utility = rep(0, cycles),
cumulative.costs = rep(0,cycles))
<- c(Local = initial_population, Recur = 0, Mets = 0, Dead = 0)
current_population
for (cycle in 1:cycles) {
if(cycle == 1){ #half cycle correction
<- sum(current_population*utilities) * 0.5
utility_value <- sum(current_population*costs) * 0.5
cost_value
}
# Transition logic
<- current_population["Local"] * pDie(cycle)
local_death <- current_population["Recur"] * pDie(cycle)
recur_death
<- current_population["Local"] * pBCA1
local_recur <- local_recur * pMets1
local_mets <- local_recur * (1 - pMets1)
local_recur
<- current_population["Recur"] * pBCA2
recur_recur <- recur_recur * pMets2
recur_mets <- recur_recur * (1 - pMets2)
recur_recur
<- current_population["Mets"] * pDieMets(cycle)
mets_death
# Update populations
"Local"] <- current_population["Local"] - local_recur - local_mets - local_death
current_population["Recur"] <- current_population["Recur"] + local_recur - recur_recur - recur_mets - recur_death
current_population["Mets"] <- current_population["Mets"] + local_mets + recur_mets - mets_death
current_population["Dead"] <- current_population["Dead"] + mets_death + local_death + recur_death
current_population[
# Update QALY
<- utility_value + (sum(current_population*utilities)/((1 + r)^cycle))
utility_value
# Update QALY
<- cost_value + (sum(current_population*costs)/((1 + r)^cycle))
cost_value
# Store results for this cycle
2:7] <- c(current_population, utility_value, cost_value)
results[cycle,
}
return(results)
}
And now we are almost done! Let’s just use this function to calculate our costs and we will export those tables to our outputs/tables
folder. So now in our CEA.R
script let’s add:
##############################################################################
# Export cost tables
##############################################################################
<- 5000
cRecur <- 20000
cMets <- 500
cLocal
<- simulate_markov_costs(strategy = 'Usual Care')
costs_usual <- simulate_markov_costs(strategy = 'Treatment')
costs_treat
write.csv(costs_usual, 'output/tables/Costs Usual Care.csv')
write.csv(costs_treat, 'output/tables/Costs Treatment.csv')
Don’t forget to commit and push your changes!!!
And congratulations! If you made it successfully through this your repository should look something like mine: