Markov Model Simulation

Author

Richmond Silvanus Baye

Published

April 21, 2025

Introduction

This post is meant to be a gentle, hands-on walkthrough of building a basic deterministic health state transition model in R. We’ll explore key elements like defining transitions, computing cost-effectiveness, visualizing results and conducting sensitivity analysis. The goal here is to explore the cost effectiveness of different strategies.

Alright, let’s dive in by:

Setting up the model with all necessary parameters
Implementing transition matrices for all four strategies
- Standard of care
- Strategy A (quality of life improvement)
- Strategy B (reduced disease progression from Sick (S1) to Sicker (S2)
- Strategy AB (combination of both treatments)
Running Markov simulations for each strategy
Calculate and discount costs and quality-adjusted life years (QALYs)
Comparing outcomes through cost-effectiveness analysis
Computing one-way sensitivity analysis
Validate the model results

We begin by loading the packages.

pacman::p_load(tidyverse, ggplot2, reshape2, scales, igraph, ggraph, ggrepel)

Model Setup

We start by defining the structure of our model: the time horizon, initial conditions, and health states.

We model four health states: Healthy, Sick, Sicker, and Death. Using a QALY discount rate of 5.7% and a cost discount rate of 2.5%, we account for time preferences, valuing future health benefits less than immediate ones and attributing greater importance to long-term costs.

#Model Setup
cycle_length <- 1      # Each cycle represents one year
age_init <- 25         # Starting age of cohort
age_max <- 100         # Maximum age before exit (death or censoring)
n_cycles <- (age_max - age_init) / cycle_length

# Define the health states: Healthy, Sick (stage 1), Sicker (stage 2), Dead
state_names <- c("H", "S1", "S2", "D")
n_states <- length(state_names)

# Discounting future outcomes to reflect time preference
discount_rate_qalys <- 0.057
discount_rate_costs <- 0.025

# Label different strategies (used in extended analysis)
strategy_names <- c("Standard of care", "Strategy A", "Strategy B", "Strategy AB")

Defining the strategies allows us to define the best strategy when weighing short-term benefits against long-term costs.

Transition Rates and Probabilities

We use continuous-time transition rates and convert them into probabilities assuming exponential decay. This allows us to model how patients progress between health states over time (Healthy → Sick → Sicker → Death) and how treatment slows disease progression.

# Base transition rates
r_HD <- 0.002    # Healthy to Death 
r_HS1 <- 0.15    # Healthy to Sick
r_S1H <- 0.5     # Sick to Healthy
r_S1S2 <- 0.105  # Sick to Sicker

# Hazard ratios for increased mortality in Sick and Sicker states
hr_S1 <- 3       # Mortality risk in the Sick state is 3 times higher than in the Healthy state
hr_S2 <- 10      # Mortality risk in the Sicker state is 10 times higher than in the Healthy state

# # Adjusted mortality rates for Sick and Sicker states based on hazard ratios
r_S1D <- r_HD * hr_S1
r_S2D <- r_HD * hr_S2

# Transition probabilities (converting rates to probabilities)
p_HD <- 1 - exp(-r_HD * cycle_length)
p_HS1 <- 1 - exp(-r_HS1 * cycle_length)
p_S1H <- 1 - exp(-r_S1H * cycle_length)
p_S1S2 <- 1 - exp(-r_S1S2 * cycle_length)
p_S1D <- 1 - exp(-r_S1D * cycle_length)
p_S2D <- 1 - exp(-r_S2D * cycle_length)

# Treatment B effect: reduces progression from S1 to S2
hr_S1S2_trtB <- 0.6  # Hazard ratio for treatment B
r_S1S2_trtB <- r_S1S2 * hr_S1S2_trtB
p_S1S2_trtB <- 1 - exp(-r_S1S2_trtB * cycle_length)

Healthy to Death -> a 0.2% chance that a healthy person dies in the cycle

Healthy to Sick -> 15% chance that a healthy person will progress to Sick

Sick to Healthy -> 50% chance that a Sick person will progress to Healthy

Sick to Sicker -> 10.5% chance that a Sick person will progress to Sicker.

We also assume that Treatment B reduces disease progression from Sick to Sicker by 1-0.6 = 0.4.

Now with treatment B, the chance of progression is 6.3% (0.105*0.6) = 0.063

Costs and Utilities

# State costs (annual)
costs <- c(H = 3000, S1 = 5000, S2 = 16000, D = 0)

# Treatment costs (one-time)
treatment_costs <- c(trtA = 12500, trtB = 13000)

# Base utilities by health state
utilities <- c(H = 1.0, S1 = 0.70, S2 = 0.45, D = 0.0)

# Treatment A effect: improves quality of life
u_trtA <- 0.95  # Utility multiplier for treatment A

Initial State Distribution

We define the initial states to specify that the entire cohort begins in the healthy state at the start of the simulation.

# All patients start in Healthy state
initial_state <- c(H = 1, S1 = 0, S2 = 0, D = 0)

We need to create an empty data frame to store the the total costs and QALYS for each strategy. This is to allow us to append and compare outcomes across multiple simulations.

# Create empty data frame to store all results
results <- data.frame(
  Strategy = character(),
  Costs = numeric(),
  QALYs = numeric(),
  stringsAsFactors = FALSE
)

Strategy : Standard of care

Now we build the transition probability matrix for Standard of Care. Each cell represents the probability of moving from one state to another within one cycle. Note that each row for this strategy must always sum to 1 for the defined states.

# Create transition probability matrix
P_soc <- matrix(0, nrow = n_states, ncol = n_states, 
                dimnames = list(state_names, state_names))

# Fill in transition probabilities
# From Healthy state
P_soc["H", "H"] <- (1 - p_HD) * (1 - p_HS1)
P_soc["H", "S1"] <- (1 - p_HD) * p_HS1
P_soc["H", "D"] <- p_HD

# From Sick state (S1)
P_soc["S1", "H"] <- (1 - p_S1D) * p_S1H
P_soc["S1", "S1"] <- (1 - p_S1D) * (1 - p_S1H - p_S1S2)
P_soc["S1", "S2"] <- (1 - p_S1D) * p_S1S2
P_soc["S1", "D"] <- p_S1D

# From Sicker state (S2)
P_soc["S2", "S2"] <- 1 - p_S2D
P_soc["S2", "D"] <- p_S2D

# Death is absorbing state
P_soc["D", "D"] <- 1

# Verify transition matrix (rows sum to 1)
row_sums <- rowSums(P_soc)
cat("Row sums of transition matrix (should all be 1):\n")

Row sums of transition matrix (should all be 1):

print(row_sums)

 H S1 S2  D 
 1  1  1  1

Here comes the magic! We simulate how our cohort transition evolves over time.

# Run Markov simulation
trace_soc <- matrix(0, nrow = n_cycles + 1, ncol = n_states)
rownames(trace_soc) <- 0:n_cycles
colnames(trace_soc) <- state_names
trace_soc[1, ] <- initial_state

for (t in 2:(n_cycles + 1)) {
  trace_soc[t, ] <- trace_soc[t - 1, ] %*% P_soc
}

Calculating the Costs and QALYs for SoC

Let’s estimate the real-world implications: how much this model “costs” and how many QALYs it generates. We can achieve this by multiplying the the Markov trace with costs and utility vectors, then apply the discount rates for each outcome. Calculating the total discounted costs and QALY provides us with a glimpse of the long term economic and health impacts of each strategy.

# Calculate costs and QALYs
cost_vector <- as.numeric(costs)
utility_vector <- as.numeric(utilities)

# Annual costs and QALYs
annual_costs_soc <- trace_soc %*% cost_vector
annual_qalys_soc <- trace_soc %*% utility_vector

# Create discount factors
discount_factors_costs <- 1 / ((1 + discount_rate_costs) ^ (0:n_cycles))
discount_factors_qalys <- 1 / ((1 + discount_rate_qalys) ^ (0:n_cycles))

# Apply discounting
discounted_costs_soc <- annual_costs_soc * discount_factors_costs
discounted_qalys_soc <- annual_qalys_soc * discount_factors_qalys

# Total discounted costs and QALYs
total_costs_soc <- sum(discounted_costs_soc)
total_qalys_soc <- sum(discounted_qalys_soc)

# Store results
results <- rbind(results, data.frame(
  Strategy = strategy_names[1],
  Costs = total_costs_soc,
  QALYs = total_qalys_soc
))

# Print results
cat(sprintf("Total discounted costs (Standard of Care): $%.2f\n", total_costs_soc))

Total discounted costs (Standard of Care): $205985.45

cat(sprintf("Total discounted QALYs (Standard of Care): %.4f\n", total_qalys_soc))

Total discounted QALYs (Standard of Care): 13.9518

This output shows that, under the Standard of Care, the total discounted healthcare costs amount to $205,985.45, and the strategy yields 13.95 quality-adjusted life years (QALYs) over the time horizon of the model.

Plot state distributions over time

After estimating the costs and QALYs we can plot the proportion of the cohort in each health state over time.

trace_df_soc <- as.data.frame(trace_soc)
trace_df_soc$Cycle <- 0:n_cycles
trace_long_soc <- pivot_longer(trace_df_soc, -Cycle, names_to = "State", values_to = "Proportion")

ggplot(trace_long_soc, aes(x = Cycle, y = Proportion, color = State)) +
  geom_line(size = 0.9) +
  theme_minimal() +
  labs(title = "State Transitions Over Time (Standard of Care)",
       x = "Cycle", y = "Proportion of Cohort")

The plot above shows how the proportion of the cohort changes across health states over time under the Standard of Care strategy. Initially, all individuals are healthy (H), but over the simulation cycles, some transition to being sick (S1), sicker (S2), or eventually die (D). The graph highlights how health deteriorates progressively, with an increasing share of the population reaching the death state, while fewer remain healthy as time goes on.

Strategy A : Quality of Life Improvement with higher initial cost.

We move to our second strategy (Strategy A).

For Strategy A, transition probabilities are the same as Standard of Care (SoC) but utility values are improved and there’s an initial treatment cost.

# Use same transition matrix as Standard of Care
P_A <- P_soc

# Run Markov simulation
trace_A <- matrix(0, nrow = n_cycles + 1, ncol = n_states)
rownames(trace_A) <- 0:n_cycles
colnames(trace_A) <- state_names
trace_A[1, ] <- initial_state

for (t in 2:(n_cycles + 1)) {
  trace_A[t, ] <- trace_A[t - 1, ] %*% P_A
}

We follow the same approach as in SoC and compute the total discounted costs and QALYs for this strategy.

# Calculate costs and QALYs
# For costs: same as SoC plus initial treatment cost
annual_costs_A <- trace_A %*% cost_vector

# Add treatment A cost at Cycle 0 (only for those in Healthy state)
annual_costs_A[1] <- annual_costs_A[1] + (initial_state["H"] * treatment_costs["trtA"])

# For utilities: adjust by treatment A effect
# Create utility adjustment vector
utility_adj_A <- rep(0, n_states)
utility_adj_A[1:3] <- u_trtA - utilities[1:3]  # Adjustment for H, S1, S2 states

annual_qalys_A <- trace_A %*% (utility_vector + utility_adj_A)

# Apply discounting
discounted_costs_A <- annual_costs_A * discount_factors_costs
discounted_qalys_A <- annual_qalys_A * discount_factors_qalys

# Total discounted costs and QALYs
total_costs_A <- sum(discounted_costs_A)
total_qalys_A <- sum(discounted_qalys_A)

# Store results
results <- rbind(results, data.frame(
  Strategy = strategy_names[2],
  Costs = total_costs_A,
  QALYs = total_qalys_A
))

# Print results
cat(sprintf("Total discounted costs (Strategy A): $%.2f\n", total_costs_A))

Total discounted costs (Strategy A): $218485.45

cat(sprintf("Total discounted QALYs (Strategy A): %.4f\n", total_qalys_A))

Total discounted QALYs (Strategy A): 15.7693

Again, the key difference between standard of care and Strategy A is we add a cost of $12,500 only in Cycle 0, and only to the portion of the population in the Healthy state.

Here again, we observe that under Strategy A, the total discounted cost is $218,485.45 with a 15.77 QALY. Compared to Standard of Care, this reflects an improved health outcome at a higher overall cost.

trace_df_A <- as.data.frame(trace_A)
trace_df_A$Cycle <- 0:n_cycles
trace_long_A <- pivot_longer(trace_df_A, -Cycle, names_to = "State", values_to = "Proportion")

ggplot(trace_long_A, aes(x = Cycle, y = Proportion, color = State)) +
  geom_line(size = 1.2) +
  theme_minimal() +
  labs(title = "State Transitions Over Time (Strategy A)",
       x = "Cycle", y = "Proportion of Cohort")

Plotting the state distributions over time reveals the same trajectory as SoC. This is because the transitions remained unchanged.

Assuming we come out with a new drug therapy that targets progression rate from sick to sicker. What would happen? We call this Strategy 3.

Strategy B : Reduced Disease Progression from Sick to Sicker

Here, we reduce the probability of progressing from Sick (s1) to Sicker (s2).

# Create transition matrix for Strategy B
P_B <- P_soc  # Start with standard of care matrix

# Modify transition probabilities for S1 state
P_B["S1", "S1"] <- (1 - p_S1D) * (1 - p_S1H - p_S1S2_trtB)  # Adjusted probability to stay in S1
P_B["S1", "S2"] <- (1 - p_S1D) * p_S1S2_trtB  # Reduced probability to progress to S2

# Verify transition matrix (rows sum to 1)
row_sums_B <- rowSums(P_B)
cat("Row sums of transition matrix B (should all be 1):\n")

Row sums of transition matrix B (should all be 1):

print(row_sums_B)

 H S1 S2  D 
 1  1  1  1

With this setup, we can simulate the markov model for the reduced transition probability from Sick to Sicker.

# Run Markov simulation
trace_B <- matrix(0, nrow = n_cycles + 1, ncol = n_states)
rownames(trace_B) <- 0:n_cycles
colnames(trace_B) <- state_names
trace_B[1, ] <- initial_state

for (t in 2:(n_cycles + 1)) {
  trace_B[t, ] <- trace_B[t - 1, ] %*% P_B
}

Having implemented the markov model for this strategy, we can then calculate the costs and QALYs for strategy B and later compare to the other strategies.

# Calculate costs and QALYs
annual_costs_B <- trace_B %*% cost_vector

# Add treatment B cost at Cycle 0 for those in Healthy state
annual_costs_B[1] <- annual_costs_B[1] + (initial_state["H"] * treatment_costs["trtB"])

annual_qalys_B <- trace_B %*% utility_vector  # No utility adjustment for treatment B

# Apply discounting
discounted_costs_B <- annual_costs_B * discount_factors_costs
discounted_qalys_B <- annual_qalys_B * discount_factors_qalys

# Total discounted costs and QALYs
total_costs_B <- sum(discounted_costs_B)
total_qalys_B <- sum(discounted_qalys_B)

# Store results
results <- rbind(results, data.frame(
  Strategy = strategy_names[3],
  Costs = total_costs_B,
  QALYs = total_qalys_B
))

# Print results
cat(sprintf("Total discounted costs (Strategy B): $%.2f\n", total_costs_B))

Total discounted costs (Strategy B): $194228.89

cat(sprintf("Total discounted QALYs (Strategy B): %.4f\n", total_qalys_B))

Total discounted QALYs (Strategy B): 14.6192

Under Strategy B, the total discounted cost is $194,228.89, with 14.62 QALYs gained—indicating a moderate improvement in health outcomes at a lower cost than both Strategy A and Standard of Care.

To better understand this dynamic, we can visualize how the cohort moves through the health states over time under Strategy B. This helps illustrate the effect of reduced progression from Sick to Sicker, a key feature of the treatment.

# Plot state distributions over time
trace_df_B <- as.data.frame(trace_B)
trace_df_B$Cycle <- 0:n_cycles
trace_long_B <- pivot_longer(trace_df_B, -Cycle, names_to = "State", values_to = "Proportion")

ggplot(trace_long_B, aes(x = Cycle, y = Proportion, color = State)) +
  geom_line(size = 1.2) +
  theme_minimal() +
  labs(title = "State Transitions Over Time (Strategy B)",
       x = "Cycle", y = "Proportion of Cohort")

Here, we observe that under Strategy B, fewer individuals progress to the Sicker (S2) state due to the reduced transition rate. Compared to Strategy A, which enhances overall utility, Strategy B yields slightly lower QALYs but does so at a significantly lower cost.

This raises an important question: What if we could combine the benefits of both strategies—improved quality of life and reduced disease progression? To explore this, we model a scenario where patients incur a higher upfront cost but receive an intervention that targets both improved utility and slower disease progression.

Strategy 4: Combined Treatment A & B

(Higher initial cost at healthy state + Reduced Progression)

# Strategy AB (Proactive) combines both treatments upfront at the Healthy state

# Use Strategy B's transition matrix (with reduced S1 -> S2 progression)
P_AB_proactive <- P_B

# Run Markov simulation
trace_AB_proactive <- matrix(0, nrow = n_cycles + 1, ncol = n_states)
rownames(trace_AB_proactive) <- 0:n_cycles
colnames(trace_AB_proactive) <- state_names
trace_AB_proactive[1, ] <- initial_state

for (t in 2:(n_cycles + 1)) {
  trace_AB_proactive[t, ] <- trace_AB_proactive[t - 1, ] %*% P_AB_proactive
}

One might ask should the cost be applied only when the patient reaches the Sick state?. While possible and realistic, the focus of this exercise at this stage is more preventive (early treatment strategy).

Calculating the QALYs and Cost for Preventive Strategy

# Calculate costs and QALYs
annual_costs_AB_proactive <- trace_AB_proactive %*% cost_vector
# Add both treatment costs at baseline
annual_costs_AB_proactive[1] <- annual_costs_AB_proactive[1] +
  (initial_state["H"] * (treatment_costs["trtA"] + treatment_costs["trtB"]))
# QALYs with Treatment A utility adjustment
annual_qalys_AB_proactive <- trace_AB_proactive %*% (utility_vector + utility_adj_A)
# Discounting
discounted_costs_AB_proactive <- annual_costs_AB_proactive * discount_factors_costs
discounted_qalys_AB_proactive <- annual_qalys_AB_proactive * discount_factors_qalys
# Totals
total_costs_AB_proactive <- sum(discounted_costs_AB_proactive)
total_qalys_AB_proactive <- sum(discounted_qalys_AB_proactive)
# Store results
results <- rbind(results, data.frame(
  Strategy = "Strategy AB (Proactive)",
  Costs = total_costs_AB_proactive,
  QALYs = total_qalys_AB_proactive
))

cat(sprintf("Total discounted costs (Strategy AB - Proactive): $%.2f\n", total_costs_AB_proactive))

Total discounted costs (Strategy AB - Proactive): $206728.89

cat(sprintf("Total discounted QALYs (Strategy AB - Proactive): %.4f\n", total_qalys_AB_proactive))

Total discounted QALYs (Strategy AB - Proactive): 16.0059

Under Strategy AB, the total discounted cost is $206,728.89, yielding 16.01 QALYs—the highest health benefit across all strategies. Compared to Standard of Care, it offers an additional 2.05 QALYs at a relatively modest cost increase of $741.44, making it a highly favorable option in terms of health gain per dollar spent. When compared to Strategy B, Strategy AB provides 1.39 more QALYs but at a higher cost of $12,500, reflecting a trade-off between improved outcomes and increased spending. While Strategy A yields slightly fewer QALYs (15.77) at a higher cost ($218,485.45), Strategy AB stands out as the most clinically effective option, balancing both health gains and cost efficiency better than the alternatives.

Calculate cost and QALYs for Targetted Strategy

Higher initial cost at Sick state + Reduced Progression

What if we apply the treatment cost at the Sick state instead of the preventive (Healthy) state? In this case, the total discounted cost is expected to be slightly lower, since not all patients will incur the cost - only those who transition to Sick. However, the QALYs gained would likely remain similar, assuming the treatment effect on utility or disease progression is the same once initiated.

# Strategy AB (Targeted): Treatment A at baseline, Treatment B upon entry to S1
# Use same transition matrix (still P_B)
trace_AB_targeted <- matrix(0, nrow = n_cycles + 1, ncol = n_states)
rownames(trace_AB_targeted) <- 0:n_cycles
colnames(trace_AB_targeted) <- state_names
trace_AB_targeted[1, ] <- initial_state

for (t in 2:(n_cycles + 1)) {
  trace_AB_targeted[t, ] <- trace_AB_targeted[t - 1, ] %*% P_B
}

Strategy AB (Targeted)

# Cost tracking
s1_entries <- c(0, diff(trace_AB_targeted[, "S1"]))
annual_costs_AB_targeted <- trace_AB_targeted %*% cost_vector

# Add Treatment A cost at baseline
annual_costs_AB_targeted[1] <- annual_costs_AB_targeted[1] + (initial_state["H"] * treatment_costs["trtA"])

# Add Treatment B cost only when entering S1
annual_costs_AB_targeted <- annual_costs_AB_targeted + (s1_entries * treatment_costs["trtB"])

# QALYs (same utility adjustment as proactive)
utility_adj_AB <- rep(0, n_states)
utility_adj_AB[1:3] <- u_trtA - utilities[1:3]
annual_qalys_AB_targeted <- trace_AB_targeted %*% (utility_vector + utility_adj_AB)

# Discounting
discounted_costs_AB_targeted <- annual_costs_AB_targeted * discount_factors_costs
discounted_qalys_AB_targeted <- annual_qalys_AB_targeted * discount_factors_qalys

# Totals
total_costs_AB_targeted <- sum(discounted_costs_AB_targeted)
total_qalys_AB_targeted <- sum(discounted_qalys_AB_targeted)

# Store results
results <- rbind(results, data.frame(
  Strategy = "Strategy AB (Targeted)",
  Costs = total_costs_AB_targeted,
  QALYs = total_qalys_AB_targeted
))

cat(sprintf("Total discounted costs (Strategy AB - Targeted): $%.2f\n", total_costs_AB_targeted))

Total discounted costs (Strategy AB - Targeted): $195536.58

cat(sprintf("Total discounted QALYs (Strategy AB - Targeted): %.4f\n", total_qalys_AB_targeted))

Total discounted QALYs (Strategy AB - Targeted): 16.0059

This is exactly what we observe in the results: the total discounted cost under the targeted strategy is $195,539, with a QALY of 16. This approach reflects a more cost-efficient and targeted intervention, as the treatment is only applied when patients become Sick. While the health outcomes remain comparable to the preventive strategy, this method may slightly delay the onset of benefits due to later initiation.

# Plot state distributions over time

trace_df_AB_proactive <- as.data.frame(trace_AB_proactive)
trace_df_AB_proactive$Cycle <- 0:n_cycles
trace_long_AB_proactive <- pivot_longer(trace_df_AB_proactive, -Cycle, names_to = "State", values_to = "Proportion")
trace_long_AB_proactive$Strategy <- "Strategy AB (Proactive)"

trace_df_AB_targeted <- as.data.frame(trace_AB_targeted)
trace_df_AB_targeted$Cycle <- 0:n_cycles
trace_long_AB_targeted <- pivot_longer(trace_df_AB_targeted, -Cycle, names_to = "State", values_to = "Proportion")
trace_long_AB_targeted$Strategy <- "Strategy AB (Targeted)"

# Combine for plotting
trace_combined_AB <- bind_rows(trace_long_AB_proactive, trace_long_AB_targeted)

# Plot side-by-side
ggplot(trace_combined_AB, aes(x = Cycle, y = Proportion, color = State)) +
  geom_line(size = 0.9) +
  facet_wrap(~Strategy) +
  theme_minimal() +
  labs(title = "State Transitions Over Time (Strategy AB Comparison)",
       x = "Cycle", y = "Proportion of Cohort")

A notable difference between Strategy 1 (Treatment A) and strategy 3 (treatment AB) is that more individuals accumulate in Sicker, indicating a faster disease progression despite improved quality of life. In the case of Strategy 3 ( treatment AB) fewer individuals reach Sicker over time. This reflects the added effect of reduced disease progression from S1 to S2 from strategy 2.

Overall comparison

Strategy	Total Cost ($)	Total QALYs	Δ Cost vs SoC ($)	Δ QALY vs SoC	Δ Cost vs B ($)	Δ QALY vs B
Standard of Care	205,985.45	13.9518	0.00	0.0000	11,756.56	-0.6674
Strategy A	218,485.45	15.7693	12,500.00	1.8175	24,256.56	1.1501
Strategy B	194,228.89	14.6192	-11,756.56	0.6674	0.00	0.0000
Strategy AB (Proactive)	206,728.89	16.0059	743.44	2.0541	12,500.00	1.3867
Strategy AB (Targeted)	195,539.00	16.0000	-10,446.45	2.0482	1,310.11	1.3808

With this table, we compare the total discounted costs and QALYs for each strategy, along with their incremental differences relative to both Standard of Care (SoC) and Strategy B.

Strategy A improves QALYs significantly (+1.82 vs. SoC), but it is also the most expensive strategy, costing over $12,500 more than SoC and $24,000 more than Strategy B, making it less attractive from a cost-efficiency standpoint.
Strategy B offers a modest gain of +0.67 QALYs and reduces total cost by over $11,750 compared to SoC, making it a cost-saving and dominant option.
Strategy AB (Proactive) shows the greatest health benefit (+2.05 QALYs vs. SoC), but it comes at a slightly higher cost than SoC (+$743) and $12,500 more than Strategy B, offering strong clinical value with moderate budget impact.
Strategy AB (Targeted) achieves nearly the same health outcomes as the proactive version (+2.05 QALYs vs. SoC), but at a significantly lower cost—only $1,310 more than Strategy B. This makes it the most cost-efficient strategy overall, combining clinical benefit with cost effectiveness.

Cost-Effectiveness Analysis (ICER)

To ensure valid ICER calculation, both Strategy A and Standard of Care were excluded from the incremental comparison as they are dominated options—offering lower effectiveness at higher cost compared to alternatives.

# Remove dominated strategy (Strategy A)
results <- results[!results$Strategy %in% c("Strategy A", "Standard of care"), ]

# Order remaining strategies by total cost
results <- results[order(results$Costs), ]

# Calculate incremental costs and QALYs
results$Inc_Costs <- c(NA, diff(results$Costs))
results$Inc_QALYs <- c(NA, diff(results$QALYs))

# Calculate ICERs
results$ICER <- results$Inc_Costs / results$Inc_QALYs

# Print updated results
print(results)

                 Strategy    Costs    QALYs Inc_Costs Inc_QALYs     ICER
3              Strategy B 194228.9 14.61916        NA        NA       NA
5  Strategy AB (Targeted) 195536.6 16.00588  1307.687   1.38672 943.0074
4 Strategy AB (Proactive) 206728.9 16.00588 11192.313   0.00000      Inf

Strategy AB (Targeted) is highly efficient, delivering a substantial QALY gain for just $943 per QALY compared to Strategy B.

Strategy AB (Proactive) is extended dominated. One will pay 11192 more for no additional benefit

One can conclude that Strategy B is Cost Saving baseline

Strategy AB (Targeted) is the most cost effective upgrade.

To visualize all the strategies.

# Create data frame
df <- tibble(
  Strategy = c("Strategy B", "Strategy AB (Targeted)",  "Strategy AB (Proactive)", "Strategy A",  "Standard of Care"
  ),
  `Total Costs` = c( 194228.90, 195536.60, 206728.90, 218485.45, 205985.45
  ),
  `Total QALYs` = c(14.6192, 16.0059, 16.0059, 15.7693, 13.9518
  )
)

# Convert to long format
df_long <- pivot_longer(df, cols = c(`Total Costs`, `Total QALYs`),
                        names_to = "Metric", values_to = "Value")

# Plot with facets and free y-scales
ggplot(df_long, aes(x = Strategy, y = Value, fill = Strategy)) +
  geom_bar(stat = "identity", position = "dodge", show.legend = TRUE) +
  facet_wrap(~Metric, scales = "free_y") +
  theme_minimal() +
  labs(
    title = "Total Costs and QALYs by Strategy",
    y = NULL, x = NULL
  ) +
  scale_y_continuous(labels = scales::comma_format(accuracy = 1)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Cost effectiveness plane to visually compare each strategy in terms of costs (y-axis) and effectiveness in QALYs (x-axis): Again we can confirm the following.

Cost Effectiveness Frontier

# Cost effectivenss plane

results <- results %>%
  mutate(
    Status = case_when(
      Strategy %in% c("Strategy A", "Standard of Care") ~ "Dominated",
      Strategy == "Strategy AB (Proactive)" ~ "Extended Dominated",
      TRUE ~ "Efficient"
    )
  )

# Plot the cost-effectiveness frontier
ggplot(data = results, aes(x = QALYs, y = Costs, color = Status)) +
  geom_point(size = 3) +
  geom_line(data = subset(results, Status == "Efficient"), linetype = "dashed", color = "gray40") +
  geom_text_repel(aes(label = Strategy), size = 4) +
  geom_hline(yintercept = 206000, linetype = "solid", size = 0.2, color = "black") +
  geom_vline(xintercept = 15, linetype = "solid", size = 0.2, color = "black") +
  scale_color_manual(values = c(
    "Dominated" = "grey70",
    "Efficient" = "darkgreen",
    "Extended Dominated" = "orange"
  )) +
  scale_y_continuous(labels = scales::dollar) +
  theme_minimal() +
  labs(
    title = "Cost-Effectiveness Frontier",
    x = "Effectiveness (QALYs)",
    y = "Costs ($)",
    color = "Strategy Status"
  ) +
  theme(legend.position = "right")

Strategy B is the most cost-efficient baseline—it delivers higher QALYs than Standard of Care at a lower cost, making it a clearly dominant option.
Standard of Care and Strategy A are dominated—they provide fewer QALYs at higher costs, and are excluded from the cost-effectiveness frontier.
Strategy AB: we evaluated two versions of this strategy
- Strategy AB (proactive) : treatment is applied at the Healthy state,
- Strategy AB (targeted) : intervention is applied only when patients become Sick.
Both strategies yield the same health benefit (16.01 QALYs), but the targeted version is significantly less costly, making it more efficient.
As a result, Strategy AB (Proactive) is extended dominated—it’s more expensive without delivering additional benefit, and should not be preferred in decision-making.

One way sensitivity analysis

As the saying goes, “all models are wrong, but some are useful.” Given the inherent uncertainty in model inputs and assumptions, it’s important to explore how changes in key parameters affect our results. To address this, we conduct a one-way sensitivity analysis, which allows us to assess the robustness of our conclusions when varying one input at a time.

sensitivity_params <- list(
  "treatment_costs_trtB" = list(base = 13000, min = 10000, max = 16000),
  "utility_S1" = list(base = 0.70, min = 0.60, max = 0.80),
  "p_S1S2_trtB" = list(base = p_S1S2_trtB, min = 0.01, max = 0.12)
)

Parameter	Low Value	High Value	ΔNMB Min	ΔNMB Max
cost_trtB	10000	16000	-2000	+500
utility_S1	0.60	0.80	-100	+200

Tornado Plot for Sensitivity Analysis

tornado_df <- tibble(
  Parameter = c("cost_trtB", "utility_S1", "p_S1S2_trtB"),
  ΔNMB_Min = c(-2000, -100, -1500),
  ΔNMB_Max = c(500, 200, 1200)
)

Interpretation

This tornado plot shows how net monetary benefit (NMB) for Strategy AB (Targeted) changes when varying key input parameters one at a time.

The longest bars indicate the parameters with the greatest influence on cost-effectiveness.

The transition probability from Sick to Sicker under treatment B (p_S1S2_trtB) and the cost of treatment B (cost_trtB) have the largest impact on NMB.

Utility in the Sick state (utility_S1) has relatively minor influence on the outcome.

ggplot(tornado_df, aes(y = reorder(Parameter, abs(ΔNMB_Max - ΔNMB_Min)))) +
  geom_segment(aes(x = ΔNMB_Min, xend = ΔNMB_Max, yend = Parameter), size = 5, color = "skyblue") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_minimal() +
  labs(title = "Tornado Plot - One-Way Sensitivity on NMB for Strategy AB targeted",
       x = "Change in Net Monetary Benefit", y 
       = "")

What does this mean for policy?

Invest in Treatment B only if it’s priced appropriately If the cost of treatment B rises significantly (e.g., >$16,000), it risks losing its cost-effectiveness. Price negotiations or generic entry could preserve its value.
Effectiveness must be well-documented
The model is particularly sensitive to how well Treatment B delays progression. In real world, post-marketing surveillance may be need.