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

1. 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.

2. Initializing the CEA.R Script: Define the following variables: Age, muMets, pBCA2, pMets1, pMets2, r, uMets, and uRecur.

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
mufem <- read_excel("data/mufem.xlsx", sheet = 'FEMALES') 

##############################################################################
# Run the simulation
##############################################################################

# Global variables
uRecur <- 0.80
uMets <- 0.40
uDead <- 0
r <- 0.03
pMets2 <- 0.90
pMets1 <- 0.75
pBCA2 <- 0.06
muMets <- 0.50
Age <- 55

3. 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.

4. 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.

pDie <- function(stage) {
  age_index <- Age + stage
  rate <- mufem$`annual hazard rate`[age_index]
  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.

pDieMets <- function(stage) {
  age_index <- Age + stage
  rate <- mufem$`annual hazard rate`[age_index] + muMets
  return(1 - exp(-rate))
}

Lets add both of these functions to our utils.R script.

5. 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 and pBCA1 = 0.02.

• “Treat”: uLocal = 0.90 and pBCA1 = 0.01.

Add the following function to your utils.R script.

# Define the Markov simulation function
simulate_markov <- function(initial_population=1, cycles=45, strategy='Usual Care') {

if(strategy == 'Usual Care'){
  pBCA1 <- 0.02
  uLocal <- 0.95
  
} else{
  pBCA1 <- 0.01
  uLocal <- 0.90
}

utilities <- c(uLocal, uRecur, uMets, uDead)

results <- data.frame(Cycle = 1:cycles, 
                      Local = rep(0, cycles), 
                      Recur = rep(0, cycles), 
                      Mets = rep(0, cycles),
                      Dead = rep(0, cycles),
                      cumulative.utility = rep(0, cycles))

current_population <- c(Local = initial_population, Recur = 0, Mets = 0, Dead = 0)

for (cycle in 1:cycles) {
  
  if(cycle == 1){ #half cycle correction
    utility_value <- sum(current_population*utilities) * 0.5
  }
  
  # Transition logic
  local_death <- current_population["Local"] * pDie(cycle)
  recur_death <- current_population["Recur"] * pDie(cycle)
  
  local_recur <- current_population["Local"] * pBCA1
  local_mets <- local_recur * pMets1
  local_recur <- local_recur * (1 - pMets1)
  
  recur_recur <- current_population["Recur"] * pBCA2
  recur_mets <- recur_recur * pMets2
  recur_recur <- recur_recur * (1 - pMets2)
  
  mets_death <- current_population["Mets"] * pDieMets(cycle)
  
  # Update populations
  current_population["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
  
  # Update QALY
  utility_value <- utility_value + (sum(current_population*utilities)/((1 + r)^cycle))

  # Store results for this cycle
  results[cycle, 2:6] <- c(current_population, utility_value)
}

return(results)
}

6. 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
uRecur <- 0.80
uMets <- 0.40
uDead <- 0
r <- 0.03
pMets2 <- 0.90
pMets1 <- 0.75
pBCA2 <- 0.06
muMets <- 0.50
Age <- 55

results_usual <- simulate_markov(strategy = 'Usual Care')
results_treat <- simulate_markov(strategy = 'Treatment')

If this code runs we are going to add one more section:

##############################################################################
# Visualization
##############################################################################

# Plot the results for both groups
plot_usual <- ggplot(results_usual, aes(x = Cycle)) +
  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', 
       plot_usual, width = 7, height = 4, bg='white')

plot_treatment <- ggplot(results_treat, aes(x = Cycle)) +
  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', 
       plot_treatment, width = 7, height = 4, bg='white')

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!!!!

7. 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:

cRecur <- 5000
cMets <- 20000
cLocal <- 500

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
simulate_markov_costs <- function(initial_population=1, cycles=45, strategy='Usual Care') {
  
  if(strategy == 'Usual Care'){
    pBCA1 <- 0.02
    uLocal <- 0.95
    cRx <- 0

  } else{
    pBCA1 <- 0.01
    uLocal <- 0.90
    cRx <- 1000
  }
  
  utilities <- c(uLocal, uRecur, uMets, uDead)
  costs <- c(cLocal + cRx, cRecur, cMets, 0)
  
  results <- data.frame(Cycle = 1:cycles, 
                        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))
  
  current_population <- c(Local = initial_population, Recur = 0, Mets = 0, Dead = 0)
  
  for (cycle in 1:cycles) {
    
    if(cycle == 1){ #half cycle correction
      utility_value <- sum(current_population*utilities) * 0.5
      cost_value <- sum(current_population*costs) * 0.5
    }
    
    # Transition logic
    local_death <- current_population["Local"] * pDie(cycle)
    recur_death <- current_population["Recur"] * pDie(cycle)
    
    local_recur <- current_population["Local"] * pBCA1
    local_mets <- local_recur * pMets1
    local_recur <- local_recur * (1 - pMets1)
    
    recur_recur <- current_population["Recur"] * pBCA2
    recur_mets <- recur_recur * pMets2
    recur_recur <- recur_recur * (1 - pMets2)
    
    mets_death <- current_population["Mets"] * pDieMets(cycle)
    
    # Update populations
    current_population["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
    
    # Update QALY
    utility_value <- utility_value + (sum(current_population*utilities)/((1 + r)^cycle))

    # Update QALY
    cost_value <- cost_value + (sum(current_population*costs)/((1 + r)^cycle))
    
    # Store results for this cycle
    results[cycle, 2:7] <- c(current_population, utility_value, cost_value)
  }
  
  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
##############################################################################

cRecur <- 5000
cMets <- 20000
cLocal <- 500

costs_usual <- simulate_markov_costs(strategy = 'Usual Care')
costs_treat <- simulate_markov_costs(strategy = 'Treatment')

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:

https://github.com/jacobjameson/BCP_Jameson_Coots/