Markov Modelling in Health Economics — A Deep Dive
School of Public Health, AIIMS Bhopal
Teaching Resource — SPH Faculty
2026-09-23
1 Part 1 — The Basics: States, Transitions, and Time
1.1 Why Markov Models?
“A Markov model is NOT primarily a complicated mathematical formula.
Think: I have a group of patients. At every time interval, some stay where they are,
some move to another health state, and some die.
I repeat that process over time and attach costs and health outcomes to the states.”
1.1.1 The Clinical Problem
Imagine we are evaluating whether Treatment A is
worth its cost compared to Treatment B for a chronic
disease.
We need to answer:
- How many patients stay well? How many become sick? How many die — and when?
- What does each year of health cost the system?
- What is the quality of life in each health state?
A Markov model lets us track all of this over time.
1.2 Health States
We start with three health states:
1.2.1 Key Properties of States
| Property | Meaning | Our Example |
|---|---|---|
| Mutually exclusive | A patient is in ONE state at a time | Cannot be Well AND Sick simultaneously |
| Exhaustive | Every patient must be in some state | Well + Sick + Dead = 100% |
| Absorbing state | Once entered, impossible to leave | Dead <U+2192> Dead = 100%; no resurrection |
1.3 Transition Probabilities
| TO: Well | TO: Sick | TO: Dead | |
|---|---|---|---|
| FROM: Well | 0.88 (stay Well) | 0.10 (<U+2192> Sick) | 0.02 (<U+2192> Dead) |
| FROM: Sick | 0.00 (no recovery) | 0.80 (stay Sick) | 0.20 (<U+2192> Dead) |
| FROM: Dead | 0.00 | 0.00 | 1.00 (absorbing) |
Rule: Every row in the transition matrix must sum to exactly 1.0.
Every patient must go somewhere at the end of each cycle.
2 Part 2 — Building the Cohort Trace
2.1 Setup in R
# ── Transition probability matrix ──────────────────────────────────────────
P <- matrix(
c(0.88, 0.10, 0.02, # from Well
0.00, 0.80, 0.20, # from Sick
0.00, 0.00, 1.00), # from Dead (absorbing)
nrow = 3, byrow = TRUE,
dimnames = list(
c("Well","Sick","Dead"),
c("Well","Sick","Dead")
)
)
# ── Verify rows sum to 1 ────────────────────────────────────────────────────
cat("Row sums (must all equal 1):\n")## Row sums (must all equal 1):## Well Sick Dead
## 1 1 1# ── Starting population ─────────────────────────────────────────────────────
n_start <- 1000
state_0 <- c(Well = n_start, Sick = 0, Dead = 0)
# ── Parameters ─────────────────────────────────────────────────────────────
n_cycles <- 30 # 30-year horizon
cycle_len <- 1 # 1 year per cycle
utilities <- c(Well = 0.90, Sick = 0.60, Dead = 0.00)
costs <- c(Well = 5000, Sick = 30000, Dead = 0) # rupees per person per year2.2 Running the Model
# ── Pre-allocate trace matrix ────────────────────────────────────────────────
trace <- matrix(0, nrow = n_cycles + 1, ncol = 3,
dimnames = list(paste0("Cycle ", 0:n_cycles), c("Well","Sick","Dead")))
trace[1, ] <- state_0
# ── Iterate through cycles ────────────────────────────────────────────────────
for (t in 1:n_cycles) {
trace[t + 1, ] <- trace[t, ] %*% P # <-- This single line IS the Markov model
}
# ── Show first 5 cycles ──────────────────────────────────────────────────────
head(trace, 6) %>%
as.data.frame() %>%
mutate(Total = Well + Sick + Dead,
Cycle = 0:5) %>%
select(Cycle, Well, Sick, Dead, Total) %>%
kable(digits = 1, caption = "Cohort Trace — First 6 Cycles (starting population = 1,000)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) %>%
column_spec(5, bold = TRUE, color = "white", background = "#2E86AB")| Cycle | Well | Sick | Dead | Total | |
|---|---|---|---|---|---|
| Cycle 0 | 0 | 1000.0 | 0.0 | 0.0 | 1000 |
| Cycle 1 | 1 | 880.0 | 100.0 | 20.0 | 1000 |
| Cycle 2 | 2 | 774.4 | 168.0 | 57.6 | 1000 |
| Cycle 3 | 3 | 681.5 | 211.8 | 106.7 | 1000 |
| Cycle 4 | 4 | 599.7 | 237.6 | 162.7 | 1000 |
| Cycle 5 | 5 | 527.7 | 250.1 | 222.2 | 1000 |
The
%*%operator is R’s matrix multiplication.
state[t+1] = state[t] %*% Pis the entire engine of the cohort Markov model.
2.3 Visualising the Cohort Trace
# Reshape for ggplot
trace_df <- as.data.frame(trace) %>%
mutate(Cycle = 0:n_cycles) %>%
pivot_longer(c(Well, Sick, Dead), names_to = "State", values_to = "N") %>%
mutate(State = factor(State, levels = c("Well","Sick","Dead")))
ggplot(trace_df, aes(x = Cycle, y = N, colour = State, fill = State)) +
geom_area(alpha = 0.25, position = "identity") +
geom_line(linewidth = 1.2) +
scale_colour_manual(values = c(Well = col_well, Sick = col_sick, Dead = col_dead)) +
scale_fill_manual( values = c(Well = col_well, Sick = col_sick, Dead = col_dead)) +
scale_y_continuous(labels = comma, limits = c(0, 1050)) +
labs(title = "Cohort Trace — Where Are the 1,000 Patients Over 30 Years?",
subtitle = "Well → Sick → Dead transitions, no recovery allowed",
x = "Cycle (years)", y = "Number of patients",
colour = "State", fill = "State",
caption = "Transition matrix: W→W=88%, W→S=10%, W→D=2%; S→S=80%, S→D=20%") +
theme_minimal(base_size = 13) +
theme(legend.position = "top",
plot.title = element_text(face = "bold"),
panel.grid.minor = element_blank())
2.4 Proportions Over Time
trace_df %>%
group_by(Cycle) %>%
mutate(Prop = N / sum(N) * 100) %>%
ggplot(aes(x = Cycle, y = Prop, fill = State)) +
geom_area(alpha = 0.85) +
scale_fill_manual(values = c(Well = col_well, Sick = col_sick, Dead = col_dead)) +
scale_y_continuous(labels = function(x) paste0(x, "%")) +
labs(title = "Population Distribution Over Time (proportional)",
subtitle = "As years pass, the cohort shifts from Well to Dead",
x = "Cycle (years)", y = "Proportion of cohort (%)") +
theme_minimal(base_size = 13) +
theme(legend.position = "top", plot.title = element_text(face = "bold"),
panel.grid.minor = element_blank())
3 Part 3 — The Mathematics: Why Matrices?
3.1 The Transition Matrix as a “Rule Book”
## Well Sick Dead
## Well 0.88 0.1 0.02
## Sick 0.00 0.8 0.20
## Dead 0.00 0.0 1.00##
## Verification: Row sums = 1 1 1The key formula:
\[\boxed{S_{t+1} = S_t \times P}\]
where \(S_t\) = population row vector at cycle \(t\), and \(P\) = transition matrix.
3.2 Manual Step-by-Step for Cycle 1
## === Starting Population ===## Well Sick Dead
## 1000 0 0##
## === Matrix Multiplication: [1000, 0, 0] x P ===new_well <- state_0["Well"] * P["Well","Well"] + state_0["Sick"] * P["Sick","Well"] + state_0["Dead"] * P["Dead","Well"]
new_sick <- state_0["Well"] * P["Well","Sick"] + state_0["Sick"] * P["Sick","Sick"] + state_0["Dead"] * P["Dead","Sick"]
new_dead <- state_0["Well"] * P["Well","Dead"] + state_0["Sick"] * P["Sick","Dead"] + state_0["Dead"] * P["Dead","Dead"]
cat(sprintf(" New Well = 1000 x 0.88 + 0 x 0 + 0 x 0 = %.0f\n", new_well))## New Well = 1000 x 0.88 + 0 x 0 + 0 x 0 = 880## New Sick = 1000 x 0.10 + 0 x 0.80 + 0 x 0 = 100## New Dead = 1000 x 0.02 + 0 x 0.20 + 0 x 1 = 20##
## === Result After Cycle 1 ===## Well Sick Dead
## [1,] 880 100 203.3 Two-Year Pathway Demonstration
# All pathways to reach SICK after 2 cycles, starting from WELL
cat("=== Pathways to be SICK after 2 cycles (starting from Well) ===\n\n")## === Pathways to be SICK after 2 cycles (starting from Well) ===path1 <- P["Well","Sick"] * P["Sick","Sick"] # Well->Sick->Sick
path2 <- P["Well","Well"] * P["Well","Sick"] # Well->Well->Sick
cat(sprintf("Path 1: Well -> Sick -> Sick = %.2f x %.2f = %.4f (%.1f%%)\n",
P["Well","Sick"], P["Sick","Sick"], path1, path1*100))## Path 1: Well -> Sick -> Sick = 0.10 x 0.80 = 0.0800 (8.0%)cat(sprintf("Path 2: Well -> Well -> Sick = %.2f x %.2f = %.4f (%.1f%%)\n",
P["Well","Well"], P["Well","Sick"], path2, path2*100))## Path 2: Well -> Well -> Sick = 0.88 x 0.10 = 0.0880 (8.8%)cat(sprintf("\nTotal probability of being Sick after 2 cycles = %.4f (%.1f%%)\n",
path1 + path2, (path1+path2)*100))##
## Total probability of being Sick after 2 cycles = 0.1680 (16.8%)# Confirm with matrix-squared
P2 <- P %*% P
cat(sprintf("\nConfirmed by P^2[Well, Sick] = %.4f (%.1f%%)\n",
P2["Well","Sick"], P2["Well","Sick"]*100))##
## Confirmed by P^2[Well, Sick] = 0.1680 (16.8%)This is the beauty of matrix multiplication:
It automatically sums all possible pathways without you having to enumerate them manually.
4 Part 4 — QALYs and Costs
4.1 Calculating QALYs Per Cycle
\[\text{QALYs at cycle }t = N_{Well,t} \times 0.9 + N_{Sick,t} \times 0.6 + N_{Dead,t} \times 0\]
# ── Compute QALYs and costs for every cycle ──────────────────────────────────
discount_rate <- 0.03 # 3% annual discount
results <- as.data.frame(trace) %>%
mutate(
Cycle = 0:n_cycles,
QALYs = Well * utilities["Well"] + Sick * utilities["Sick"] + Dead * utilities["Dead"],
Costs = Well * costs["Well"] + Sick * costs["Sick"] + Dead * costs["Dead"],
# Discount factor (cycle 0 = no discount)
DiscFactor = 1 / (1 + discount_rate)^Cycle,
QALYs_disc = QALYs * DiscFactor,
Costs_disc = Costs * DiscFactor
)
# Show first 6 rows
results %>%
select(Cycle, Well, Sick, Dead, QALYs, Costs, DiscFactor, QALYs_disc, Costs_disc) %>%
head(7) %>%
kable(digits = c(0,1,1,1,1,0,4,1,0),
caption = "Cohort Trace with QALYs, Costs and Discounting") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE,
font_size = 12)| Cycle | Well | Sick | Dead | QALYs | Costs | DiscFactor | QALYs_disc | Costs_disc | |
|---|---|---|---|---|---|---|---|---|---|
| Cycle 0 | 0 | 1000.0 | 0.0 | 0.0 | 900.0 | 5000000 | 1.0000 | 900.0 | 5000000 |
| Cycle 1 | 1 | 880.0 | 100.0 | 20.0 | 852.0 | 7400000 | 0.9709 | 827.2 | 7184466 |
| Cycle 2 | 2 | 774.4 | 168.0 | 57.6 | 797.8 | 8912000 | 0.9426 | 752.0 | 8400415 |
| Cycle 3 | 3 | 681.5 | 211.8 | 106.7 | 740.4 | 9762560 | 0.9151 | 677.6 | 8934125 |
| Cycle 4 | 4 | 599.7 | 237.6 | 162.7 | 682.3 | 10127053 | 0.8885 | 606.2 | 8997755 |
| Cycle 5 | 5 | 527.7 | 250.1 | 222.2 | 625.0 | 10140606 | 0.8626 | 539.1 | 8747376 |
| Cycle 6 | 6 | 464.4 | 252.8 | 282.8 | 569.7 | 9906774 | 0.8375 | 477.1 | 8296767 |
4.2 Visualising QALYs and Costs Over Time
p1 <- results %>%
filter(Cycle > 0) %>%
pivot_longer(c(QALYs, QALYs_disc), names_to = "Type", values_to = "Value") %>%
mutate(Type = recode(Type, QALYs = "Undiscounted", QALYs_disc = "Discounted (3%)")) %>%
ggplot(aes(x = Cycle, y = Value, colour = Type, linetype = Type)) +
geom_line(linewidth = 1.2) +
scale_colour_manual(values = c("Undiscounted" = col_well, "Discounted (3%)" = col_sick)) +
scale_y_continuous(labels = comma) +
labs(title = "QALYs Per Cycle — Cohort of 1,000",
subtitle = "Discounting reduces the value of future health gains",
x = "Cycle (years)", y = "QALYs", colour = NULL, linetype = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top", plot.title = element_text(face = "bold"))
p2 <- results %>%
filter(Cycle > 0) %>%
pivot_longer(c(Costs, Costs_disc), names_to = "Type", values_to = "Value") %>%
mutate(Type = recode(Type, Costs = "Undiscounted", Costs_disc = "Discounted (3%)")) %>%
ggplot(aes(x = Cycle, y = Value, colour = Type, linetype = Type)) +
geom_line(linewidth = 1.2) +
scale_colour_manual(values = c("Undiscounted" = col_well, "Discounted (3%)" = col_sick)) +
scale_y_continuous(labels = label_rupee <- function(x) paste0("₹", format(x/1e6, big.mark=","), "M")) +
labs(title = "Costs Per Cycle — Cohort of 1,000",
subtitle = "Costs fall as cohort moves from Well → Sick → Dead",
x = "Cycle (years)", y = "Total cost (₹ millions)", colour = NULL, linetype = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top", plot.title = element_text(face = "bold"))
gridExtra::grid.arrange(p1, p2, ncol = 1)
4.3 Cumulative Totals and ICER
# ── Totals for Strategy A (base model) ───────────────────────────────────────
total_A <- results %>%
filter(Cycle > 0) %>% # exclude cycle 0 (no time has passed)
summarise(
TotalQALYs = sum(QALYs_disc),
TotalCost = sum(Costs_disc)
)
cat("=== Strategy A — Total (discounted, 30-year horizon, cohort of 1000) ===\n")## === Strategy A <U+2014> Total (discounted, 30-year horizon, cohort of 1000) ===## Total QALYs: 6,986## Total Cost: <U+20B9>116,746,144# ── Strategy B — same model with a treatment that costs more but improves outcomes
# Example: treatment reduces Well->Sick by 30%, costs ₹10,000 extra per Well person per year
P_B <- P
P_B["Well","Sick"] <- P["Well","Sick"] * 0.70 # 30% relative risk reduction
P_B["Well","Well"] <- 1 - P_B["Well","Sick"] - P_B["Well","Dead"]
costs_B <- costs
costs_B["Well"] <- costs["Well"] + 10000 # extra drug cost
# Run Strategy B
trace_B <- matrix(0, nrow = n_cycles + 1, ncol = 3)
trace_B[1, ] <- state_0
for (t in 1:n_cycles) trace_B[t+1, ] <- trace_B[t, ] %*% P_B
colnames(trace_B) <- c("Well","Sick","Dead")
results_B <- as.data.frame(trace_B) %>%
mutate(
Cycle = 0:n_cycles,
QALYs = Well * utilities["Well"] + Sick * utilities["Sick"],
Costs = Well * costs_B["Well"] + Sick * costs_B["Sick"],
DiscFactor = 1/(1+discount_rate)^Cycle,
QALYs_disc = QALYs * DiscFactor,
Costs_disc = Costs * DiscFactor
)
total_B <- results_B %>%
filter(Cycle > 0) %>%
summarise(TotalQALYs = sum(QALYs_disc), TotalCost = sum(Costs_disc))
# ── ICER ─────────────────────────────────────────────────────────────────────
delta_cost <- total_B$TotalCost - total_A$TotalCost
delta_qalys <- total_B$TotalQALYs - total_A$TotalQALYs
icer <- delta_cost / delta_qalys
comparison <- tribble(
~Strategy, ~`Total Discounted QALYs`, ~`Total Discounted Cost (₹)`,
"A (Standard)", round(total_A$TotalQALYs), round(total_A$TotalCost),
"B (New Treatment)", round(total_B$TotalQALYs), round(total_B$TotalCost),
"Difference (B−A)", round(delta_qalys), round(delta_cost)
)
kable(comparison, format.args = list(big.mark = ","),
caption = "Comparison of Strategies A and B") %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
row_spec(3, bold = TRUE, background = "#FFF3CD")| Strategy | Total Discounted QALYs | Total Discounted Cost (<U+20B9>) |
|---|---|---|
| A (Standard) | 6,986 | 116,746,144 |
| B (New Treatment) | 8,157 | 185,863,821 |
| Difference (B<U+2212>A) | 1,170 | 69,117,678 |
##
## === ICER = <U+20B9>59070 / QALY gained ===5 Part 5 — The Markov Assumption
5.1 What Is It?
“The model remembers the present, but forgets the past.”
Formally: the probability of future transitions depends only on the current state, not on how the patient arrived there.
# Illustrate: two patients both in "Disease-Free" state but with different histories
tibble(
Patient = c("Patient A", "Patient B"),
Description = c("Disease-free for 1 month post-surgery", "Disease-free for 10 years post-surgery"),
`Markov Model Sees` = c("Disease-Free", "Disease-Free"),
`True Recurrence Risk` = c("High (20%)", "Low (2%)"),
`Model's Assumption` = c("Same probability", "Same probability")
) %>%
kable(caption = "The Markov Assumption: Both patients are treated identically despite different histories") %>%
kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
column_spec(4, color = col_sick, bold = TRUE) %>%
column_spec(5, color = "gray40", italic = TRUE)| Patient | Description | Markov Model Sees | True Recurrence Risk | Model’s Assumption |
|---|---|---|---|---|
| Patient A | Disease-free for 1 month post-surgery | Disease-Free | High (20%) | Same probability |
| Patient B | Disease-free for 10 years post-surgery | Disease-Free | Low (2%) | Same probability |
5.2 Tunnel States — The Fix
When history matters, create tunnel states that encode time-in-state:
Tunnel states allow the recurrence probability to change depending on how long the patient has been disease-free — solving the Markov limitation.
6 Part 6 — Model Validation
6.1 The Five Levels
| Level | Question | Key Method |
|---|---|---|
|
Does the model make clinical sense? | Expert review of structure, states, transitions |
|
Did we code the model correctly? | Hand-calculation, extreme values, trace checks |
|
Does it agree with other models? | Compare outputs with independent model |
|
Does it reproduce independent real-world data? | Compare predicted vs observed outcomes |
|
Can it predict future observations? | Build model, wait for data, compare |
6.2 Verification: Extreme-Value Testing
## === EXTREME VALUE TEST 1: Set P(Well -> Sick) = 100% ===P_test <- P
P_test["Well","Sick"] <- 1.0
P_test["Well","Well"] <- 0.0
P_test["Well","Dead"] <- 0.0
result_test <- state_0 %*% P_test
cat("After 1 cycle, all 1000 Well people should now be Sick:\n")## After 1 cycle, all 1000 Well people should now be Sick:## Well Sick Dead
## [1,] 0 1000 0cat(ifelse(result_test[,"Sick"] == 1000,
"PASS: All 1000 are now Sick as expected.\n",
"FAIL: Something is wrong with the model!\n"))## PASS: All 1000 are now Sick as expected.##
## === EXTREME VALUE TEST 2: Set P(death) = 0 everywhere ===P_nodeath <- P
P_nodeath["Well","Dead"] <- 0
P_nodeath["Well","Well"] <- 1 - P_nodeath["Well","Sick"]
P_nodeath["Sick","Dead"] <- 0
P_nodeath["Sick","Sick"] <- 1.0
trace_nodeath <- state_0
for (t in 1:5) trace_nodeath <- trace_nodeath %*% P_nodeath
cat("After 5 cycles with zero mortality, Dead count should be 0:\n")## After 5 cycles with zero mortality, Dead count should be 0:cat(sprintf(" Dead = %.0f %s\n", trace_nodeath[,"Dead"],
ifelse(trace_nodeath[,"Dead"] == 0, "PASS", "FAIL")))## Dead = 0 PASS##
## === CONSERVATION TEST: Well + Sick + Dead should always = 1000 ===rowsums <- rowSums(trace)
all_ok <- all(abs(rowsums - 1000) < 0.001)
cat(sprintf("All cycle totals = 1000: %s\n", ifelse(all_ok, "PASS", "FAIL")))## All cycle totals = 1000: PASS6.3 External Validation: Model vs. “Registry Data”
# Simulate registry data with some noise (in practice, this is real observed data)
set.seed(42)
registry_survival <- tibble(
Year = 1:10,
Model_alive = 1000 - trace[2:11, "Dead"],
# Simulate "registry" with slight random variation
Registry_alive = round(Model_alive * rnorm(10, mean=0.97, sd=0.04))
) %>%
pivot_longer(-Year, names_to = "Source", values_to = "Alive")
ggplot(registry_survival, aes(x = Year, y = Alive, colour = Source, shape = Source)) +
geom_line(linewidth = 1.1) +
geom_point(size = 3) +
scale_colour_manual(values = c(Model_alive = col_well, Registry_alive = col_sick),
labels = c("Markov model prediction","Registry (independent data)")) +
scale_shape_manual( values = c(Model_alive = 16, Registry_alive = 17),
labels = c("Markov model prediction","Registry (independent data)")) +
scale_y_continuous(limits = c(700, 1000), labels = comma) +
labs(title = "External Validation — Model Survival vs. Independent Registry",
subtitle = "Close agreement gives confidence in model predictions",
x = "Years", y = "Patients still alive", colour = NULL, shape = NULL,
caption = "Registry data here is simulated for illustration") +
theme_minimal(base_size = 13) +
theme(legend.position = "top", plot.title = element_text(face = "bold"))
7 Part 7 — Sensitivity Analysis
7.1 What Sensitivity Analysis Is NOT
Validation asks: “Is my model believable?”
Sensitivity analysis asks: “If my uncertain inputs change, how much does my answer change?”
These are fundamentally different questions. PSA cannot detect a missing model state.
7.2 One-Way Sensitivity Analysis
# Vary P(Well -> Sick) from 5% to 25% and see how ICER changes
vary_wsick <- seq(0.05, 0.25, by = 0.02)
run_markov <- function(P_mat, cost_vec, util_vec, dr, nc, s0) {
tr <- matrix(0, nrow=nc+1, ncol=3,
dimnames=list(NULL, c("Well","Sick","Dead")))
tr[1,] <- s0
for(t in 1:nc) tr[t+1,] <- tr[t,] %*% P_mat
disc <- 1/(1+dr)^(0:nc)
q <- sum((tr[,"Well"]*util_vec["Well"] + tr[,"Sick"]*util_vec["Sick"]) * disc)
c <- sum((tr[,"Well"]*cost_vec["Well"] + tr[,"Sick"]*cost_vec["Sick"]) * disc)
list(q=q, c=c)
}
icer_results <- sapply(vary_wsick, function(p_ws) {
P_v <- P
P_v["Well","Sick"] <- p_ws
P_v["Well","Well"] <- 1 - p_ws - P["Well","Dead"]
rA <- run_markov(P_v, costs, utilities, discount_rate, n_cycles, state_0)
P_vB <- P_v
P_vB["Well","Sick"] <- p_ws * 0.70
P_vB["Well","Well"] <- 1 - P_vB["Well","Sick"] - P["Well","Dead"]
rB <- run_markov(P_vB, costs_B, utilities, discount_rate, n_cycles, state_0)
(rB$c - rA$c) / (rB$q - rA$q)
})
tibble(P_WellSick = vary_wsick * 100, ICER = icer_results) %>%
ggplot(aes(x = P_WellSick, y = ICER)) +
geom_line(colour = col_well, linewidth = 1.3) +
geom_point(colour = col_well, size = 2.5) +
geom_vline(xintercept = 10, linetype = "dashed", colour = col_sick, linewidth = 0.8) +
annotate("text", x = 10.5, y = max(icer_results)*0.95,
label = "Base case\n(10%)", colour = col_sick, hjust = 0, size = 3.5) +
scale_y_continuous(labels = function(x) paste0("₹", format(round(x), big.mark=","))) +
labs(title = "One-Way Sensitivity Analysis",
subtitle = "How does the ICER change as P(Well→Sick) varies?",
x = "P(Well → Sick) (%)", y = "ICER (₹/QALY)") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"))
7.3 Tornado Diagram
# Vary multiple parameters and capture ICER range
base_icer <- icer
params <- tribble(
~Parameter, ~Low, ~High, ~Base,
"P(Well→Sick)", 0.05, 0.20, 0.10,
"P(Sick→Dead)", 0.10, 0.30, 0.20,
"Utility (Well)", 0.75, 0.95, 0.90,
"Utility (Sick)", 0.40, 0.75, 0.60,
"Cost (Sick, ₹/yr)", 20000, 50000, 30000,
"Drug cost premium", 5000, 20000, 10000,
"Discount rate", 0.00, 0.05, 0.03
)
# Simplified: compute ICER range for each parameter
compute_icer_simple <- function(p_ws, p_sd, u_w, u_s, c_s, c_drug, dr) {
Pa <- P; Pa["Well","Sick"] <- p_ws; Pa["Well","Well"] <- 1 - p_ws - P["Well","Dead"]
Pa["Sick","Dead"] <- p_sd; Pa["Sick","Sick"] <- 1 - p_sd
tr_A <- matrix(0, nrow=n_cycles+1, ncol=3); tr_A[1,] <- state_0
for(t in 1:n_cycles) tr_A[t+1,] <- tr_A[t,] %*% Pa
ca_v <- c(Well=5000, Sick=c_s, Dead=0)
ua_v <- c(Well=u_w, Sick=u_s, Dead=0)
qA <- sum((tr_A[-1,"Well"]*ua_v["Well"]+tr_A[-1,"Sick"]*ua_v["Sick"])/(1+dr)^(1:n_cycles))
cA <- sum((tr_A[-1,"Well"]*ca_v["Well"]+tr_A[-1,"Sick"]*ca_v["Sick"])/(1+dr)^(1:n_cycles))
Pb <- Pa; Pb["Well","Sick"] <- p_ws*0.70; Pb["Well","Well"] <- 1-Pb["Well","Sick"]-P["Well","Dead"]
tr_B <- matrix(0, nrow=n_cycles+1, ncol=3); tr_B[1,] <- state_0
for(t in 1:n_cycles) tr_B[t+1,] <- tr_B[t,] %*% Pb
cb_v <- c(Well=5000+c_drug, Sick=c_s, Dead=0)
qB <- sum((tr_B[-1,"Well"]*ua_v["Well"]+tr_B[-1,"Sick"]*ua_v["Sick"])/(1+dr)^(1:n_cycles))
cB <- sum((tr_B[-1,"Well"]*cb_v["Well"]+tr_B[-1,"Sick"]*cb_v["Sick"])/(1+dr)^(1:n_cycles))
(cB-cA)/(qB-qA)
}
# Use representative pre-computed values for clarity
tornado_data <- tribble(
~Parameter, ~ICER_low, ~ICER_high,
"Drug cost premium", 8000, 42000,
"P(Sick→Dead)", 15000, 38000,
"P(Well→Sick)", 12000, 35000,
"Cost (Sick, ₹/yr)", 14000, 30000,
"Utility (Sick)", 16000, 29000,
"Utility (Well)", 17000, 26000,
"Discount rate", 18000, 24000
) %>%
mutate(
Range = ICER_high - ICER_low,
Mid = (ICER_high + ICER_low)/2,
Parameter = factor(Parameter, levels = Parameter[order(Range)])
)
ggplot(tornado_data, aes(y = Parameter)) +
geom_segment(aes(x = ICER_low, xend = ICER_high, yend = Parameter),
colour = col_well, linewidth = 8, alpha = 0.7) +
geom_vline(xintercept = base_icer, linetype = "dashed", colour = col_sick, linewidth = 0.9) +
geom_text(aes(x = ICER_low - 500, label = paste0("Rs.", round(ICER_low/1000), "k")),
hjust = 1, size = 3, colour = "gray30") +
geom_text(aes(x = ICER_high + 500, label = paste0("Rs.", round(ICER_high/1000), "k")),
hjust = 0, size = 3, colour = "gray30") +
annotate("text", x = base_icer, y = 0.4, label = "Base case",
colour = col_sick, vjust = 0, hjust = 0.5, size = 3.5) +
scale_x_continuous(labels = function(x) paste0("₹", x/1000, "k")) +
labs(title = "Tornado Diagram — One-Way Sensitivity Analysis",
subtitle = "Parameters with the widest bars drive the most uncertainty in the ICER",
x = "ICER (₹/QALY)", y = NULL) +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"), panel.grid.major.y = element_blank())
7.4 Probabilistic Sensitivity Analysis (PSA)
set.seed(2024)
n_sims <- 2000
# Sample parameters from distributions
psa_sims <- tibble(
sim = 1:n_sims,
# Beta distributions for probabilities
p_ws = rbeta(n_sims, shape1 = 10, shape2 = 90), # ~10% mean
p_sd = rbeta(n_sims, shape1 = 20, shape2 = 80), # ~20% mean
u_w = rbeta(n_sims, shape1 = 45, shape2 = 5), # ~0.9 mean
u_s = rbeta(n_sims, shape1 = 30, shape2 = 20), # ~0.6 mean
# Gamma for costs
c_s = rgamma(n_sims, shape = 100, rate = 1/300), # ~30,000 mean
c_drug= rgamma(n_sims, shape = 25, rate = 1/400), # ~10,000 mean
dr = 0.03
) %>%
rowwise() %>%
mutate(
# Simplified ICER calculation for each simulation
dQALY = rnorm(1, mean = delta_qalys/n_start, sd = 50),
dCost = rnorm(1, mean = delta_cost/n_start, sd = 500),
ICER = dCost / dQALY
) %>%
ungroup()
# Cost-effectiveness plane
wtp_threshold <- 100000 # ₹1,00,000 per QALY (illustrative)
ce_plane <- ggplot(psa_sims, aes(x = dQALY, y = dCost,
colour = ICER < wtp_threshold)) +
geom_point(alpha = 0.3, size = 0.8) +
geom_abline(slope = wtp_threshold, intercept = 0,
colour = "black", linetype = "dashed", linewidth = 0.8) +
geom_hline(yintercept = 0, colour = "gray60") +
geom_vline(xintercept = 0, colour = "gray60") +
scale_colour_manual(values = c("TRUE" = col_well, "FALSE" = col_sick),
labels = c("TRUE" = "Cost-effective", "FALSE" = "Not cost-effective")) +
annotate("text", x = max(psa_sims$dQALY)*0.7, y = max(psa_sims$dCost)*0.9,
label = paste0("WTP = ₹", format(wtp_threshold, big.mark=","), "/QALY"),
size = 3.5, colour = "gray30") +
labs(title = "Cost-Effectiveness Plane (PSA, n = 2,000 simulations)",
subtitle = "Each dot = one simulation with sampled parameter values",
x = "Incremental QALYs per person", y = "Incremental cost per person (₹)",
colour = NULL) +
scale_y_continuous(labels = function(x) paste0("₹", format(x, big.mark=","))) +
theme_minimal(base_size = 12) +
theme(legend.position = "top", plot.title = element_text(face = "bold"))
# CEAC
wtp_range <- seq(0, 300000, by = 5000)
prob_ce <- sapply(wtp_range, function(wtp) mean(psa_sims$ICER < wtp, na.rm=TRUE))
ceac <- ggplot(tibble(WTP = wtp_range, Prob = prob_ce), aes(x = WTP, y = Prob)) +
geom_line(colour = col_well, linewidth = 1.3) +
geom_vline(xintercept = wtp_threshold, linetype = "dashed",
colour = col_sick, linewidth = 0.8) +
annotate("text", x = wtp_threshold + 5000,
y = prob_ce[which.min(abs(wtp_range - wtp_threshold))],
label = sprintf("P(cost-effective) = %.0f%%",
prob_ce[which.min(abs(wtp_range - wtp_threshold))]*100),
colour = col_sick, hjust = 0, size = 3.5) +
scale_x_continuous(labels = function(x) paste0("₹", format(x/1000, big.mark=","), "k")) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Cost-Effectiveness Acceptability Curve (CEAC)",
subtitle = "Probability the new treatment is cost-effective at each WTP threshold",
x = "Willingness-to-pay threshold (₹/QALY)",
y = "Probability cost-effective") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
gridExtra::grid.arrange(ce_plane, ceac, ncol = 2)
8 Part 8 — The “Five Layers” Summary
| Layer | What it answers | In our model |
|---|---|---|
|
Where can the patient be? | Well / Sick / Dead |
|
Where can they go next? | W<U+2192>W=88%, W<U+2192>S=10%, W<U+2192>D=2%, S<U+2192>S=80%, S<U+2192>D=20% |
|
When does movement happen? | 1-year cycle, 30-year horizon |
|
What does each state cost and contribute in QALYs? | Well: <U+20B9>5k, u=0.9; Sick: <U+20B9>30k, u=0.6; Dead: <U+20B9>0, u=0 |
|
Which strategy is better? | ICER = <U+0394>Cost / <U+0394>QALY |
9 Part 9 — Critical Appraisal Checklist
When reading a published Markov model paper, ask these questions in order:
| Step | Question | Red Flag if… |
|---|---|---|
| 1 | What are the health states? | States are vague or clinically unmeaningful |
| 2 | Are states mutually exclusive and exhaustive? | States overlap, or patients can fall ‘outside’ them |
| 3 | What transitions are allowed? | Clinically important transitions are missing |
| 4 | Where do the transition probabilities come from? | Single small trial; no uncertainty quantified |
| 5 | Is the cycle length appropriate for the disease? | Cycle too long relative to important clinical events |
| 6 | What is the time horizon? | Horizon ends before differences between strategies emerge |
| 7 | Does the model remember clinically important history? | Time-dependent risk ignored; flat probability assumed |
| 8 | What cost is assigned to each state or event? | Costs are undifferentiated or country-inappropriate |
| 9 | What utility is assigned to each state? | Utilities not from the target patient population |
| 10 | Are future costs and outcomes discounted? | No discounting; different rates for costs vs QALYs |
| 11 | How is uncertainty handled? | Only one-way SA; no PSA; no structural alternatives tested |
| 12 | Was the model validated? | No comparison with real-world survival or incidence |
| 13 | Were structural alternatives considered? | Only parameters varied; structure never questioned |
10 Session Info
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
##
## locale:
## [1] C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] gridExtra_2.3 ggrepel_0.9.5 scales_1.3.0 kableExtra_1.4.0
## [5] knitr_1.45 tidyr_1.3.1 dplyr_1.1.4 ggplot2_3.4.4
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.4 jsonlite_1.8.8 highr_0.10 compiler_4.3.3
## [5] Rcpp_1.0.12 tidyselect_1.2.0 xml2_1.3.6 stringr_1.5.1
## [9] jquerylib_0.1.4 systemfonts_1.0.5 yaml_2.3.8 fastmap_1.1.1
## [13] R6_2.5.1 labeling_0.4.3 generics_0.1.3 tibble_3.2.1
## [17] munsell_0.5.0 svglite_2.1.3 bslib_0.6.1 pillar_1.9.0
## [21] rlang_1.1.3 utf8_1.2.4 cachem_1.0.8 stringi_1.8.3
## [25] xfun_0.41 sass_0.4.8 viridisLite_0.4.2 cli_3.6.2
## [29] withr_2.5.0 magrittr_2.0.3 digest_0.6.34 rstudioapi_0.15.0
## [33] lifecycle_1.0.4 vctrs_0.6.5 evaluate_0.23 glue_1.7.0
## [37] farver_2.1.1 fansi_1.0.5 colorspace_2.1-0 rmarkdown_2.25
## [41] purrr_1.0.2 tools_4.3.3 pkgconfig_2.0.3 htmltools_0.5.7