Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
A district clinic looks after 1,000 people with diabetes who have no kidney disease yet. The health office wants a four-year forecast: how many will have kidney disease after four years, how many will have died, and what the group will cost and gain in health. This article builds that forecast one step at a time. First, name three states, healthy, sick and dead, and write the chance of each move in one year. Then move the whole group a year at a time: after the first year 800 are still healthy, 150 are sick and 50 have died, and the total stays the same. Last, decide which count earns each year. Counting only the end of year one gives 0.890 years of full health per person; averaging its start and its end gives 0.945, because people leave the healthy state all through the year. Short R boxes reproduce every number.
A clinic, 1,000 people and a four-year forecast
A district clinic looks after 1,000 people with type 2 diabetes who have no kidney disease yet. The district health office plans its budget four years ahead and asks the clinic manager three questions.
First, how many of these people will have kidney disease after four years, and how many will have died? Second, what will the clinic spend on them? Third, how much healthy life will the group have?
The clinic's records give the yearly chances, synthetic here for teaching. Each year, of the people without kidney disease, a share of 0.15 develop it and a share of 0.05 die. Of the people with kidney disease, a share of 0.15 die each year, and nobody is cured.
A yearly check costs the clinic 100 THB per person per year without kidney disease and 400 THB per person per year with it. Multiplying 1,000 people by four years cannot answer the questions, because the group changes every year: some people move from healthy to sick, and some die. The model has to move the people first and count their cost and health second.
Name the states and the yearly chances
Start by naming the places a person can be at any moment. Here there are three: Healthy, meaning diabetes without kidney disease, Sick, meaning diabetes with kidney disease, and Dead. Each person is in exactly one of them at a time, and these places are the model's health states. The two short names are this article's own: Healthy does not mean cured of diabetes, and Sick means kidney disease has been added.
The model moves in fixed steps of one year, and each step is called a cycle. Over one cycle, a healthy person stays healthy with a chance of 0.80, becomes sick with 0.15 and dies with 0.05. A sick person stays sick with 0.85 and dies with 0.15. A dead person stays dead: Dead has no way out, so it is called an absorbing state.
A transition probability is the chance of moving between two states within one cycle, converted from a rate so that time stays in the calculation. A rate of 0.12 per year over a cycle of 0.25 years gives 0.0296, not the shortcut 0.0300; the full lesson is in Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation. Here one cycle is one year and the clinic's chances are already yearly, so no conversion is needed.
Write the chances as a table. Each row is where people are now, and each column is where they are one cycle later. A row shares out everyone who starts in that state, so its chances must add up to 1: 0.80 + 0.15 + 0.05 for Healthy, 0.85 + 0.15 for Sick. A row that adds up to more or less than 1 creates or loses people.
The yearly chances, from each state (row) to each state (column)
| From | To Healthy | To Sick | To Dead |
|---|---|---|---|
| Healthy | 0.80 | 0.15 | 0.05 |
| Sick | 0 | 0.85 | 0.15 |
| Dead | 0 | 0 | 1 |
Move the whole group one year at a time
All 1,000 people start healthy. During the first year, the chances in the Healthy row share them out: 1,000 × 0.80 = 800 stay healthy, 1,000 × 0.15 = 150 become sick and 1,000 × 0.05 = 50 die. The three counts add up to 1,000 again, so nobody has been lost.
The second year starts from those counts, and now two groups move at once. Sick receives people from two places: new cases from the 800 healthy people, and the sick people who stay sick. The worked example below adds them up.
This way of following people has a name. A Markov model follows a group of people as they move between a few health states, one cycle at a time. Where people go next depends only on where they are now, not on how they got there. Because it follows one whole group, the cohort, rather than one person at a time, it is a Markov cohort model.
Worked example: the second year
Year two starts from 800 healthy, 150 sick and 50 dead. Each count is shared out along its own row of chances, and each state then adds up what it receives.
-
Healthy keeps a share of 0.80
\[ 800 \times 0.80 = 640 \]
Nobody comes back from Sick, so Healthy only loses people: 640 remain.
-
Sick gains new cases and keeps most of its own
\[ 800 \times 0.15 + 150 \times 0.85 = 120 + 127.5 = 247.5 \]
120 healthy people become sick this year, and 127.5 of last year's 150 stay sick.
-
Dead keeps everyone and gains from both
\[ 50 + 800 \times 0.05 + 150 \times 0.15 = 50 + 40 + 22.5 = 112.5 \]
The 50 who died in year one stay dead, and 40 healthy and 22.5 sick people die in year two.
-
Check the total
\[ 640 + 247.5 + 112.5 = 1000 \]
The total is still 1,000 people, so the year has moved people without losing or creating any.
Result: After two years the clinic has 640 healthy, 247.5 sick and 112.5 dead. A count such as 247.5 is an average over many groups like this one, not half a person.
The cohort trace
Repeat the same step and the group keeps moving. The table of counts at the end of every year is called the cohort trace. Down the Healthy column the count falls by a factor of 0.80 each year, from 1,000 to 409.6 after four years. Sick rises to 337.219 and Dead only grows, to 253.181, while every row still adds up to 1,000 people.
The cohort trace: people in each state at the end of each year
| End of year | Healthy | Sick | Dead | Total |
|---|---|---|---|---|
| 0 (start) | 1,000 | 0 | 0 | 1,000 |
| 1 | 800 | 150 | 50 | 1,000 |
| 2 | 640 | 247.5 | 112.5 | 1,000 |
| 3 | 512 | 306.375 | 181.625 | 1,000 |
| 4 | 409.6 | 337.219 | 253.181 | 1,000 |
Check the trace in R
R code trace.R
# 1,000 healthy people move between Healthy, Sick and Dead, one year per cycle.
# Writes: nothing to disk; the results print in the console.
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
moves <- matrix(transition, nrow = 3, byrow = TRUE)
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
flows <- trace[2, ] * moves
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")Output of the run trace.log
Row sums of the chances: 1 1 1 | totals, years 0 to 4: 1000 1000 1000 1000 1000
Year 2 flows from Healthy | Sick | Dead: 640 120 40 | 0 127.5 22.5 | 0 0 50
Healthy, years 0 to 4: 1000 800 640 512 409.6
Sick, years 0 to 4: 0 150 247.5 306.375 337.219
Dead, years 0 to 4: 0 50 112.5 181.625 253.181Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · works in R 4.6.1 · Download building-a-markov-cohort-model-step-by-step.R
Full tidy script
# Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
# Follow 1,000 people through Healthy, Sick and Dead for four years, then count health and cost.
# ติดตามคน 1,000 คนผ่านสถานะสุขภาพดี ป่วย และเสียชีวิต สี่รอบรายปี แล้วนับสุขภาพและต้นทุน
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the states and their yearly values as an explicit table.
library(dplyr) # Add the per-cycle values to that table.
library(ggplot2) # Display the trace without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
cycle_length <- 1
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
states <- tibble::tibble(state = c("Healthy", "Sick", "Dead"), health = health_value,
cost = yearly_cost) |>
dplyr::mutate(health_per_cycle = health * cycle_length, cost_per_cycle = cost * cycle_length)
moves <- matrix(transition, nrow = 3, byrow = TRUE)
stopifnot(all(abs(rowSums(moves) - 1) < 1e-9))
#%%@note 2. Move the cohort one year at a time
# Multiply the counts by the chances each year, and check that the total stays at 1,000.
# คูณจำนวนคนด้วยโอกาสย้ายสถานะทีละปี และตรวจว่ายอดรวมยังเท่ากับ 1,000
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
stopifnot(all(abs(rowSums(trace) - sum(start)) < 1e-9))
flows <- trace[2, ] * moves
people <- trace
#%%@note 3. Count each year at its start, at its end and as the average
# The end count of one year is the start count of the next, so the same trace gives all three.
# จำนวนตอนท้ายปีหนึ่งคือจำนวนตอนต้นปีถัดไป trace เดียวจึงให้การนับได้ทั้งสามแบบ
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
#%%@note 4. Bring later years back to today
# Divide each year's value by 1.03 once for every year that has passed by its end.
# หารค่าของแต่ละปีด้วย 1.03 หนึ่งครั้งต่อทุกปีที่ผ่านไปเมื่อสิ้นปีนั้น
factors <- (1 + rate)^(1:4)
rewards <- cbind(health = states$health_per_cycle, cost = states$cost_per_cycle)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- cbind(start = colSums(ways$start / factors), end = colSums(ways$end / factors),
average = colSums(ways$average / factors))
#%%@note 5. Rework the three foundation examples
# A yearly rate becomes a chance per cycle, a weight turns years into health, a later cost shrinks.
# อัตรารายปีกลายเป็นโอกาสต่อรอบ ค่าน้ำหนักเปลี่ยนปีเป็นสุขภาพ และต้นทุนในอนาคตมีค่าน้อยลงในวันนี้
cycle_chance <- signif(1 - exp(-rate_per_year * cycle_years), 3)
shortcut <- formatC(rate_per_year * cycle_years, format = "f", digits = 4)
weighted_years <- formatC(years_in_state * c(state_value, 1 - state_value), format = "f",
digits = 1)
#%%@note 6. Display the trace
# Plot the number of people in each state at the end of each year, no graphics file is saved.
# แสดงจำนวนคนในแต่ละสถานะเมื่อสิ้นแต่ละปี โดยไม่บันทึกไฟล์ภาพ
trace_table <- tibble::tibble(year = rep(0:cycles, times = 3),
state = rep(states$state, each = cycles + 1), people = as.vector(trace))
print(ggplot2::ggplot(trace_table, ggplot2::aes(x = year, y = people, colour = state)) +
ggplot2::geom_line())
#%%@note 7. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")
cat("Rate 0.12, 0.25 years, chance | shortcut:", cycle_chance, shortcut, "\n")
cat("Ten years at 0.70 | gain from full health:", weighted_years, "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")
Which count earns the year?
The manager's last two questions need a value for each year a person spends in each state, in money and in health.
A QALY, a quality-adjusted life year, is a year of life weighted by its quality, from 1 for full health to 0 for death. Over 10 years, a state valued at 0.70 counts as 7.0 QALYs, so restoring full health gains 3.0; the full lesson is in QALY, DALY and the Utilities Behind Them: What Is Actually Being Measured.
For teaching, the clinic counts a year without kidney disease as full health, 1, a year with kidney disease as 0.6, and a year dead as 0 QALYs. With the yearly costs of 100 THB and 400 THB, each state now carries two values per person per year.
What one person's full year in each state is worth
| State | Cost (THB per year) | Health (QALYs per year) |
|---|---|---|
| Healthy | 100 | 1 |
| Sick | 400 | 0.6 |
| Dead | 0 | 0 |
Year one, counted three ways
Now take year one. At its start all 1,000 people are healthy, which counts as 1,000 × 1 = 1,000 QALYs for the group. At its end 800 are healthy and 150 sick, which counts as 800 × 1 + 150 × 0.6 = 890 QALYs. Neither is right: the people who left Healthy did not all leave on the first day of the year or on the last.
If people leave evenly through the year, the group's average count over the year lies halfway between its two counts. So average them: (1,000 + 890) / 2 = 945 QALYs, or 0.945 per person, instead of 0.890 from the end count or 1.000 from the start. Costs follow the same rule: 1,000 × 100 = 100,000 THB counted at the start and 800 × 100 + 150 × 400 = 140,000 THB at the end, so (100,000 + 140,000) / 2 = 120,000 THB on average.
This averaging has a name. Half-cycle correction is averaging the cohort counts at the start and the end of each cycle, as if people moved halfway through it (Briggs and colleagues). It changes what each year earns, not where anyone goes: the trace above is untouched. Some reports call the same calculation the trapezoidal method, because the shaded area under the blue line in the figure below is a trapezoid.
Check year one in R
R code cycle-one.R
# Year one's cost and health, counted at the start, at the end and as the average of the two.
# Writes: nothing to disk; the results print in the console.
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")Output of the run cycle-one.log
Year 1 health, start | end | average: 1000 890 945
Year 1 health per person: 1.000 0.890 0.945
Year 1 cost THB, start | end | average: 100,000 140,000 120,000Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · works in R 4.6.1 · Download building-a-markov-cohort-model-step-by-step.R
Full tidy script
# Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
# Follow 1,000 people through Healthy, Sick and Dead for four years, then count health and cost.
# ติดตามคน 1,000 คนผ่านสถานะสุขภาพดี ป่วย และเสียชีวิต สี่รอบรายปี แล้วนับสุขภาพและต้นทุน
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the states and their yearly values as an explicit table.
library(dplyr) # Add the per-cycle values to that table.
library(ggplot2) # Display the trace without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
cycle_length <- 1
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
states <- tibble::tibble(state = c("Healthy", "Sick", "Dead"), health = health_value,
cost = yearly_cost) |>
dplyr::mutate(health_per_cycle = health * cycle_length, cost_per_cycle = cost * cycle_length)
moves <- matrix(transition, nrow = 3, byrow = TRUE)
stopifnot(all(abs(rowSums(moves) - 1) < 1e-9))
#%%@note 2. Move the cohort one year at a time
# Multiply the counts by the chances each year, and check that the total stays at 1,000.
# คูณจำนวนคนด้วยโอกาสย้ายสถานะทีละปี และตรวจว่ายอดรวมยังเท่ากับ 1,000
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
stopifnot(all(abs(rowSums(trace) - sum(start)) < 1e-9))
flows <- trace[2, ] * moves
people <- trace
#%%@note 3. Count each year at its start, at its end and as the average
# The end count of one year is the start count of the next, so the same trace gives all three.
# จำนวนตอนท้ายปีหนึ่งคือจำนวนตอนต้นปีถัดไป trace เดียวจึงให้การนับได้ทั้งสามแบบ
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
#%%@note 4. Bring later years back to today
# Divide each year's value by 1.03 once for every year that has passed by its end.
# หารค่าของแต่ละปีด้วย 1.03 หนึ่งครั้งต่อทุกปีที่ผ่านไปเมื่อสิ้นปีนั้น
factors <- (1 + rate)^(1:4)
rewards <- cbind(health = states$health_per_cycle, cost = states$cost_per_cycle)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- cbind(start = colSums(ways$start / factors), end = colSums(ways$end / factors),
average = colSums(ways$average / factors))
#%%@note 5. Rework the three foundation examples
# A yearly rate becomes a chance per cycle, a weight turns years into health, a later cost shrinks.
# อัตรารายปีกลายเป็นโอกาสต่อรอบ ค่าน้ำหนักเปลี่ยนปีเป็นสุขภาพ และต้นทุนในอนาคตมีค่าน้อยลงในวันนี้
cycle_chance <- signif(1 - exp(-rate_per_year * cycle_years), 3)
shortcut <- formatC(rate_per_year * cycle_years, format = "f", digits = 4)
weighted_years <- formatC(years_in_state * c(state_value, 1 - state_value), format = "f",
digits = 1)
#%%@note 6. Display the trace
# Plot the number of people in each state at the end of each year, no graphics file is saved.
# แสดงจำนวนคนในแต่ละสถานะเมื่อสิ้นแต่ละปี โดยไม่บันทึกไฟล์ภาพ
trace_table <- tibble::tibble(year = rep(0:cycles, times = 3),
state = rep(states$state, each = cycles + 1), people = as.vector(trace))
print(ggplot2::ggplot(trace_table, ggplot2::aes(x = year, y = people, colour = state)) +
ggplot2::geom_line())
#%%@note 7. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")
cat("Rate 0.12, 0.25 years, chance | shortcut:", cycle_chance, shortcut, "\n")
cat("Ten years at 0.70 | gain from full health:", weighted_years, "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")
Four years, counted the same way
Each later year works the same way, and the trace already holds every count needed. The end of one year is the start of the next, so year two runs from 890 to 640 × 1 + 247.5 × 0.6 = 788.5 QALYs and averages (890 + 788.5) / 2 = 839.25. Years three and four repeat the step with the next two rows of the trace, and the table puts the three ways of counting side by side.
The group's QALYs in each year, three ways of counting (not yet brought to today's value)
| Year | Counted at the start | Counted at the end | Average of the two |
|---|---|---|---|
| 1 | 1,000 | 890 | 945 |
| 2 | 890 | 788.5 | 839.25 |
| 3 | 788.5 | 695.825 | 742.16 |
| 4 | 695.825 | 611.931 | 653.88 |
| Four years | 3,374.325 | 2,986.256 | 3,180.29 |
What the columns show
Read down the columns. The start column counts every year as if nobody had left yet, and the end column as if everybody who left had gone on the first day. Over four years the start count gives 3,374.325 QALYs and the end count 2,986.256; the average, 3,180.29, sits halfway between them.
The two edge counts differ a lot because the group changes fast within each year. If the counts barely changed within a year, all three columns would agree. When they change quickly, a shorter cycle, such as three months, is safer than any averaging rule.
Check the four years in R
R code four-years.R
# Four years of health from the same trace, counted at the start, at the end and as the average.
# Writes: nothing to disk; the results print in the console.
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
health_value <- c(1, 0.6, 0)
people <- matrix(c(1000, 0, 0), nrow = 5, ncol = 3, byrow = TRUE)
for (year in 1:4) people[year + 1, ] <- people[year, ] %*% matrix(transition, 3, byrow = TRUE)
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")Output of the run four-years.log
Health by year, start: 1000 890 788.5 695.825 | total: 3374.325
Health by year, end: 890 788.5 695.825 611.931 | total: 2986.256
Health by year, average: 945 839.25 742.16 653.88 | total: 3180.29Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · works in R 4.6.1 · Download building-a-markov-cohort-model-step-by-step.R
Full tidy script
# Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
# Follow 1,000 people through Healthy, Sick and Dead for four years, then count health and cost.
# ติดตามคน 1,000 คนผ่านสถานะสุขภาพดี ป่วย และเสียชีวิต สี่รอบรายปี แล้วนับสุขภาพและต้นทุน
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the states and their yearly values as an explicit table.
library(dplyr) # Add the per-cycle values to that table.
library(ggplot2) # Display the trace without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
cycle_length <- 1
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
states <- tibble::tibble(state = c("Healthy", "Sick", "Dead"), health = health_value,
cost = yearly_cost) |>
dplyr::mutate(health_per_cycle = health * cycle_length, cost_per_cycle = cost * cycle_length)
moves <- matrix(transition, nrow = 3, byrow = TRUE)
stopifnot(all(abs(rowSums(moves) - 1) < 1e-9))
#%%@note 2. Move the cohort one year at a time
# Multiply the counts by the chances each year, and check that the total stays at 1,000.
# คูณจำนวนคนด้วยโอกาสย้ายสถานะทีละปี และตรวจว่ายอดรวมยังเท่ากับ 1,000
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
stopifnot(all(abs(rowSums(trace) - sum(start)) < 1e-9))
flows <- trace[2, ] * moves
people <- trace
#%%@note 3. Count each year at its start, at its end and as the average
# The end count of one year is the start count of the next, so the same trace gives all three.
# จำนวนตอนท้ายปีหนึ่งคือจำนวนตอนต้นปีถัดไป trace เดียวจึงให้การนับได้ทั้งสามแบบ
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
#%%@note 4. Bring later years back to today
# Divide each year's value by 1.03 once for every year that has passed by its end.
# หารค่าของแต่ละปีด้วย 1.03 หนึ่งครั้งต่อทุกปีที่ผ่านไปเมื่อสิ้นปีนั้น
factors <- (1 + rate)^(1:4)
rewards <- cbind(health = states$health_per_cycle, cost = states$cost_per_cycle)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- cbind(start = colSums(ways$start / factors), end = colSums(ways$end / factors),
average = colSums(ways$average / factors))
#%%@note 5. Rework the three foundation examples
# A yearly rate becomes a chance per cycle, a weight turns years into health, a later cost shrinks.
# อัตรารายปีกลายเป็นโอกาสต่อรอบ ค่าน้ำหนักเปลี่ยนปีเป็นสุขภาพ และต้นทุนในอนาคตมีค่าน้อยลงในวันนี้
cycle_chance <- signif(1 - exp(-rate_per_year * cycle_years), 3)
shortcut <- formatC(rate_per_year * cycle_years, format = "f", digits = 4)
weighted_years <- formatC(years_in_state * c(state_value, 1 - state_value), format = "f",
digits = 1)
#%%@note 6. Display the trace
# Plot the number of people in each state at the end of each year, no graphics file is saved.
# แสดงจำนวนคนในแต่ละสถานะเมื่อสิ้นแต่ละปี โดยไม่บันทึกไฟล์ภาพ
trace_table <- tibble::tibble(year = rep(0:cycles, times = 3),
state = rep(states$state, each = cycles + 1), people = as.vector(trace))
print(ggplot2::ggplot(trace_table, ggplot2::aes(x = year, y = people, colour = state)) +
ggplot2::geom_line())
#%%@note 7. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")
cat("Rate 0.12, 0.25 years, chance | shortcut:", cycle_chance, shortcut, "\n")
cat("Ten years at 0.70 | gain from full health:", weighted_years, "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")
Bring later years back to today
One more step remains before the forecast is complete. Discounting is converting a future cost or health gain into its value today, shrinking it by the discount rate for each year away. At 3% a cost of 1,834 THB due in year 3 is worth 1,678 THB today; the full lesson is in Whose Cost, in Which Year? Perspective, Opportunity Cost, Inflation and Discounting.
The clinic uses 3% a year for both costs and QALYs, and counts each year's value at the end of that year. Year one's value is divided by 1.03, year two's by 1.03 × 1.03 = 1.0609, year three's by 1.0927 and year four's by 1.1255. For the averaged count that gives 945 / 1.03 = 917.5, 839.25 / 1.0609 = 791.1, 742.16 / 1.0927 = 679.2 and 653.88 / 1.1255 = 581.0 QALYs today. Added before rounding, the four give 2,968.695 QALYs.
Check the yearly factors in R
R code today.R
# Each year's averaged health brought back to today at 3% a year, counted at the end of the year.
# Writes: nothing to disk; the results print in the console.
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
factors <- (1 + rate)^(1:4)
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")Output of the run today.log
Rate: 3% | factors: 1.0300 1.0609 1.0927 1.1255
Average health by year, today's value: 917.5 791.1 679.2 581 | total: 2968.695Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · works in R 4.6.1 · Download building-a-markov-cohort-model-step-by-step.R
Full tidy script
# Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
# Follow 1,000 people through Healthy, Sick and Dead for four years, then count health and cost.
# ติดตามคน 1,000 คนผ่านสถานะสุขภาพดี ป่วย และเสียชีวิต สี่รอบรายปี แล้วนับสุขภาพและต้นทุน
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the states and their yearly values as an explicit table.
library(dplyr) # Add the per-cycle values to that table.
library(ggplot2) # Display the trace without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
cycle_length <- 1
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
states <- tibble::tibble(state = c("Healthy", "Sick", "Dead"), health = health_value,
cost = yearly_cost) |>
dplyr::mutate(health_per_cycle = health * cycle_length, cost_per_cycle = cost * cycle_length)
moves <- matrix(transition, nrow = 3, byrow = TRUE)
stopifnot(all(abs(rowSums(moves) - 1) < 1e-9))
#%%@note 2. Move the cohort one year at a time
# Multiply the counts by the chances each year, and check that the total stays at 1,000.
# คูณจำนวนคนด้วยโอกาสย้ายสถานะทีละปี และตรวจว่ายอดรวมยังเท่ากับ 1,000
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
stopifnot(all(abs(rowSums(trace) - sum(start)) < 1e-9))
flows <- trace[2, ] * moves
people <- trace
#%%@note 3. Count each year at its start, at its end and as the average
# The end count of one year is the start count of the next, so the same trace gives all three.
# จำนวนตอนท้ายปีหนึ่งคือจำนวนตอนต้นปีถัดไป trace เดียวจึงให้การนับได้ทั้งสามแบบ
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
#%%@note 4. Bring later years back to today
# Divide each year's value by 1.03 once for every year that has passed by its end.
# หารค่าของแต่ละปีด้วย 1.03 หนึ่งครั้งต่อทุกปีที่ผ่านไปเมื่อสิ้นปีนั้น
factors <- (1 + rate)^(1:4)
rewards <- cbind(health = states$health_per_cycle, cost = states$cost_per_cycle)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- cbind(start = colSums(ways$start / factors), end = colSums(ways$end / factors),
average = colSums(ways$average / factors))
#%%@note 5. Rework the three foundation examples
# A yearly rate becomes a chance per cycle, a weight turns years into health, a later cost shrinks.
# อัตรารายปีกลายเป็นโอกาสต่อรอบ ค่าน้ำหนักเปลี่ยนปีเป็นสุขภาพ และต้นทุนในอนาคตมีค่าน้อยลงในวันนี้
cycle_chance <- signif(1 - exp(-rate_per_year * cycle_years), 3)
shortcut <- formatC(rate_per_year * cycle_years, format = "f", digits = 4)
weighted_years <- formatC(years_in_state * c(state_value, 1 - state_value), format = "f",
digits = 1)
#%%@note 6. Display the trace
# Plot the number of people in each state at the end of each year, no graphics file is saved.
# แสดงจำนวนคนในแต่ละสถานะเมื่อสิ้นแต่ละปี โดยไม่บันทึกไฟล์ภาพ
trace_table <- tibble::tibble(year = rep(0:cycles, times = 3),
state = rep(states$state, each = cycles + 1), people = as.vector(trace))
print(ggplot2::ggplot(trace_table, ggplot2::aes(x = year, y = people, colour = state)) +
ggplot2::geom_line())
#%%@note 7. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")
cat("Rate 0.12, 0.25 years, chance | shortcut:", cycle_chance, shortcut, "\n")
cat("Ten years at 0.70 | gain from full health:", weighted_years, "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")
Four-year totals
Do the same for the start and end counts, and for cost, and the three ways of counting give these four-year totals at today's value. The average again sits between the two edge counts, for cost as well as for health. The end count gives the highest total cost, because it counts each year with the later head count, when more people have moved into the costly Sick state.
Four years at today's value, for the whole group
| Way of counting each year | Cost (THB) | QALYs |
|---|---|---|
| At the start | 532,594 | 3,149.605 |
| At the end | 604,810 | 2,787.786 |
| Average (half-cycle correction) | 568,702 | 2,968.695 |
The forecast
So the manager can answer. Counting each year at its middle, the 1,000 people will cost the clinic 568,702 THB in today's money over the four years. They will live 2,968.695 QALYs, which is 569 THB and 2.969 QALYs per person. The trace answers the first question: after four years the expected counts are 337.219 people alive with kidney disease and 253.181 deaths in all, averages over many groups like this one.
Check the four-year totals in R
R code totals.R
# Four years of health and cost at today's value, for each way of counting the year.
# Writes: nothing to disk; the results print in the console.
rate <- 0.03
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
rewards <- cbind(health = c(1, 0.6, 0), cost = c(100, 400, 0))
people <- matrix(c(1000, 0, 0), nrow = 5, ncol = 3, byrow = TRUE)
for (year in 1:4) people[year + 1, ] <- people[year, ] %*% matrix(transition, 3, byrow = TRUE)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- sapply(ways, function(rows) colSums(rows / (1 + rate)^(1:4)))
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")Output of the run totals.log
Four-year health today, start | end | average: 3149.605 2787.786 2968.695
Four-year cost today THB, start | end | average: 532594 604810 568702
Per person, average, health | cost THB: 2.969 569Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · works in R 4.6.1 · Download building-a-markov-cohort-model-step-by-step.R
Full tidy script
# Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
# Follow 1,000 people through Healthy, Sick and Dead for four years, then count health and cost.
# ติดตามคน 1,000 คนผ่านสถานะสุขภาพดี ป่วย และเสียชีวิต สี่รอบรายปี แล้วนับสุขภาพและต้นทุน
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the states and their yearly values as an explicit table.
library(dplyr) # Add the per-cycle values to that table.
library(ggplot2) # Display the trace without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
cycle_length <- 1
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
states <- tibble::tibble(state = c("Healthy", "Sick", "Dead"), health = health_value,
cost = yearly_cost) |>
dplyr::mutate(health_per_cycle = health * cycle_length, cost_per_cycle = cost * cycle_length)
moves <- matrix(transition, nrow = 3, byrow = TRUE)
stopifnot(all(abs(rowSums(moves) - 1) < 1e-9))
#%%@note 2. Move the cohort one year at a time
# Multiply the counts by the chances each year, and check that the total stays at 1,000.
# คูณจำนวนคนด้วยโอกาสย้ายสถานะทีละปี และตรวจว่ายอดรวมยังเท่ากับ 1,000
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
stopifnot(all(abs(rowSums(trace) - sum(start)) < 1e-9))
flows <- trace[2, ] * moves
people <- trace
#%%@note 3. Count each year at its start, at its end and as the average
# The end count of one year is the start count of the next, so the same trace gives all three.
# จำนวนตอนท้ายปีหนึ่งคือจำนวนตอนต้นปีถัดไป trace เดียวจึงให้การนับได้ทั้งสามแบบ
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
#%%@note 4. Bring later years back to today
# Divide each year's value by 1.03 once for every year that has passed by its end.
# หารค่าของแต่ละปีด้วย 1.03 หนึ่งครั้งต่อทุกปีที่ผ่านไปเมื่อสิ้นปีนั้น
factors <- (1 + rate)^(1:4)
rewards <- cbind(health = states$health_per_cycle, cost = states$cost_per_cycle)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- cbind(start = colSums(ways$start / factors), end = colSums(ways$end / factors),
average = colSums(ways$average / factors))
#%%@note 5. Rework the three foundation examples
# A yearly rate becomes a chance per cycle, a weight turns years into health, a later cost shrinks.
# อัตรารายปีกลายเป็นโอกาสต่อรอบ ค่าน้ำหนักเปลี่ยนปีเป็นสุขภาพ และต้นทุนในอนาคตมีค่าน้อยลงในวันนี้
cycle_chance <- signif(1 - exp(-rate_per_year * cycle_years), 3)
shortcut <- formatC(rate_per_year * cycle_years, format = "f", digits = 4)
weighted_years <- formatC(years_in_state * c(state_value, 1 - state_value), format = "f",
digits = 1)
#%%@note 6. Display the trace
# Plot the number of people in each state at the end of each year, no graphics file is saved.
# แสดงจำนวนคนในแต่ละสถานะเมื่อสิ้นแต่ละปี โดยไม่บันทึกไฟล์ภาพ
trace_table <- tibble::tibble(year = rep(0:cycles, times = 3),
state = rep(states$state, each = cycles + 1), people = as.vector(trace))
print(ggplot2::ggplot(trace_table, ggplot2::aes(x = year, y = people, colour = state)) +
ggplot2::geom_line())
#%%@note 7. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")
cat("Rate 0.12, 0.25 years, chance | shortcut:", cycle_chance, shortcut, "\n")
cat("Ten years at 0.70 | gain from full health:", weighted_years, "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")
Check the three borrowed examples in R
R code foundations.R
# The three foundation examples behind the recaps: a yearly rate, a health weight and a later cost.
# Writes: nothing to disk; the results print in the console.
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
rate <- 0.03
cat("Rate 0.12, 0.25 years, chance | shortcut:", signif(1 - exp(-rate_per_year * cycle_years), 3),
formatC(rate_per_year * cycle_years, format = "f", digits = 4), "\n")
cat("Ten years at 0.70 | gain from full health:",
formatC(years_in_state * c(state_value, 1 - state_value), format = "f", digits = 1), "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")Output of the run foundations.log
Rate 0.12, 0.25 years, chance | shortcut: 0.0296 0.0300
Ten years at 0.70 | gain from full health: 7.0 3.0
Cost due in year 3 at 3%, today THB: 1678Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · ggplot2 4.0.3 · works in R 4.6.1 · Download building-a-markov-cohort-model-step-by-step.R
Full tidy script
# Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time
# Follow 1,000 people through Healthy, Sick and Dead for four years, then count health and cost.
# ติดตามคน 1,000 คนผ่านสถานะสุขภาพดี ป่วย และเสียชีวิต สี่รอบรายปี แล้วนับสุขภาพและต้นทุน
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# ggplot2 4.0.3 (stable)
# install.packages(c("tibble", "dplyr", "ggplot2"))
# update.packages(oldPkgs = c("tibble", "dplyr", "ggplot2"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
# remotes::install_version("ggplot2", "4.0.3")
library(tibble) # Write the states and their yearly values as an explicit table.
library(dplyr) # Add the per-cycle values to that table.
library(ggplot2) # Display the trace without saving files.
stopifnot(getRversion() >= "4.5.2")
#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
first_start <- c(1000, 0, 0)
first_end <- c(800, 150, 50)
health_value <- c(1, 0.6, 0)
yearly_cost <- c(100, 400, 0)
rate <- 0.03
average_by_year <- c(945, 839.25, 742.1625, 653.878125)
rate_per_year <- 0.12
cycle_years <- 0.25
years_in_state <- 10
state_value <- 0.70
future_cost <- 1834
cycle_length <- 1
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
states <- tibble::tibble(state = c("Healthy", "Sick", "Dead"), health = health_value,
cost = yearly_cost) |>
dplyr::mutate(health_per_cycle = health * cycle_length, cost_per_cycle = cost * cycle_length)
moves <- matrix(transition, nrow = 3, byrow = TRUE)
stopifnot(all(abs(rowSums(moves) - 1) < 1e-9))
#%%@note 2. Move the cohort one year at a time
# Multiply the counts by the chances each year, and check that the total stays at 1,000.
# คูณจำนวนคนด้วยโอกาสย้ายสถานะทีละปี และตรวจว่ายอดรวมยังเท่ากับ 1,000
trace <- matrix(start, nrow = cycles + 1, ncol = 3, byrow = TRUE)
for (year in seq_len(cycles)) trace[year + 1, ] <- trace[year, ] %*% moves
stopifnot(all(abs(rowSums(trace) - sum(start)) < 1e-9))
flows <- trace[2, ] * moves
people <- trace
#%%@note 3. Count each year at its start, at its end and as the average
# The end count of one year is the start count of the next, so the same trace gives all three.
# จำนวนตอนท้ายปีหนึ่งคือจำนวนตอนต้นปีถัดไป trace เดียวจึงให้การนับได้ทั้งสามแบบ
counts <- rbind(first_start, first_end)
health <- c(counts %*% health_value, mean(counts %*% health_value))
spend <- c(counts %*% yearly_cost, mean(counts %*% yearly_cost))
health_edges <- people %*% health_value
at_start <- health_edges[1:4]
at_end <- health_edges[2:5]
average <- (at_start + at_end) / 2
#%%@note 4. Bring later years back to today
# Divide each year's value by 1.03 once for every year that has passed by its end.
# หารค่าของแต่ละปีด้วย 1.03 หนึ่งครั้งต่อทุกปีที่ผ่านไปเมื่อสิ้นปีนั้น
factors <- (1 + rate)^(1:4)
rewards <- cbind(health = states$health_per_cycle, cost = states$cost_per_cycle)
edges <- people %*% rewards
ways <- list(start = edges[1:4, ], end = edges[2:5, ], average = (edges[1:4, ] + edges[2:5, ]) / 2)
today <- cbind(start = colSums(ways$start / factors), end = colSums(ways$end / factors),
average = colSums(ways$average / factors))
#%%@note 5. Rework the three foundation examples
# A yearly rate becomes a chance per cycle, a weight turns years into health, a later cost shrinks.
# อัตรารายปีกลายเป็นโอกาสต่อรอบ ค่าน้ำหนักเปลี่ยนปีเป็นสุขภาพ และต้นทุนในอนาคตมีค่าน้อยลงในวันนี้
cycle_chance <- signif(1 - exp(-rate_per_year * cycle_years), 3)
shortcut <- formatC(rate_per_year * cycle_years, format = "f", digits = 4)
weighted_years <- formatC(years_in_state * c(state_value, 1 - state_value), format = "f",
digits = 1)
#%%@note 6. Display the trace
# Plot the number of people in each state at the end of each year, no graphics file is saved.
# แสดงจำนวนคนในแต่ละสถานะเมื่อสิ้นแต่ละปี โดยไม่บันทึกไฟล์ภาพ
trace_table <- tibble::tibble(year = rep(0:cycles, times = 3),
state = rep(states$state, each = cycles + 1), people = as.vector(trace))
print(ggplot2::ggplot(trace_table, ggplot2::aes(x = year, y = people, colour = state)) +
ggplot2::geom_line())
#%%@note 7. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Row sums of the chances:", rowSums(moves), "| totals, years 0 to 4:", rowSums(trace), "\n")
cat("Year 2 flows from Healthy | Sick | Dead:", flows[1, ], "|", flows[2, ], "|", flows[3, ], "\n")
cat("Healthy, years 0 to 4:", round(trace[, 1], 3), "\n")
cat("Sick, years 0 to 4:", round(trace[, 2], 3), "\n")
cat("Dead, years 0 to 4:", round(trace[, 3], 3), "\n")
cat("Year 1 health, start | end | average:", health, "\n")
cat("Year 1 health per person:", sprintf("%.3f", health / 1000), "\n")
cat("Year 1 cost THB, start | end | average:", big(spend), "\n")
cat("Health by year, start:", at_start, "| total:", sum(at_start), "\n")
cat("Health by year, end:", round(at_end, 3), "| total:", round(sum(at_end), 3), "\n")
cat("Health by year, average:", round(average, 2), "| total:", round(sum(average), 2), "\n")
cat("Rate:", paste0(100 * rate, "%"), "| factors:",
formatC(factors, format = "f", digits = 4), "\n")
cat("Average health by year, today's value:", round(average_by_year / factors, 1),
"| total:", round(sum(average_by_year / factors), 3), "\n")
cat("Four-year health today, start | end | average:", sprintf("%.3f", today["health", ]), "\n")
cat("Four-year cost today THB, start | end | average:", round(today["cost", ]), "\n")
cat("Per person, average, health | cost THB:", sprintf("%.3f", today["health", "average"] / 1000),
round(today["cost", "average"] / 1000), "\n")
cat("Rate 0.12, 0.25 years, chance | shortcut:", cycle_chance, shortcut, "\n")
cat("Ten years at 0.70 | gain from full health:", weighted_years, "\n")
cat("Cost due in year 3 at 3%, today THB:", round(future_cost / (1 + rate)^3), "\n")
Checks before you trust the trace
Three checks catch most arithmetic mistakes in a model like this. Every row of chances adds up to 1. Every year's total is still 1,000 people. Dead never shrinks, because nobody leaves it.
A fourth check is clinical, not arithmetic. The model remembers only where a person is now, so someone who became sick last year faces the same 0.15 chance of dying as someone who has been sick for ten years. If that chance depends on how long a person has been sick, split Sick into a first-year state and a later state, so the model can remember the time (Siebert and colleagues).
Write the timing rules beside the results. Here one cycle is one year, and the chances move people once per cycle. Each year counts as the average of its start and its end, and its value is brought to today at its end.
Two analysts with the same chances but different counting rules report different totals, as the table of four-year totals shows. The short boxes use only base R, the plain approach that Jalal and colleagues describe for health decision models; the downloadable script adds tibble, dplyr and ggplot2 to draw the trace.
Terms as this model uses them
- Health state (สถานะสุขภาพ)
- One of the few places a person can be at any moment: Healthy, Sick or Dead here.
- Cycle (รอบ)
- The fixed time step the model moves in; one year here.
- Transition probability (ความน่าจะเป็นของการเปลี่ยนสถานะ)
- The chance of moving from one state to another within one cycle; each row of them adds up to 1.
- Absorbing state (สถานะดูดกลืน)
- A state nobody leaves once they enter it, such as Dead.
- Markov model (แบบจำลองมาร์คอฟ)
- A model that follows a group between a few health states one cycle at a time; where people go next depends only on where they are now.
- Cohort trace (ตารางจำนวนคนรายรอบ)
- The table of how many people are in each state at the end of each cycle.
- Half-cycle correction (การปรับครึ่งรอบ)
- Counting each cycle as the average of its start and end counts, as if people moved halfway through it.
Connect to the foundation articles
The foundation article on how decision trees and cohort models turn evidence into a decision shows where this model sits among the other model types. The one on QALYs and the weights behind them shows where a weight such as 0.6 comes from. The one on costing and today's value works through the 3% rate.
Try it in the app
Open our Markov model app to change the inputs and watch the result update. It needs a login and a Pro or Max plan.
The app starts from its own example values, so its numbers will not match this article.
References
- Siebert Uwe, Alagoz Oguzhan, Bayoumi Ahmed M., Jahn Beate, Owens Douglas K., Cohen David J., Kuntz Karen M.. State-Transition Modeling: A Report of the ISPOR-SMDM Modeling Good Research Practices Task Force-3. Value in Health. 2012;15(6):812-820. doi:10.1016/j.jval.2012.06.014. https://doi.org/10.1016/j.jval.2012.06.014
- Andrew Briggs, Karl Claxton, Mark Sculpher. Decision Modelling for Health Economic Evaluation (Handbooks for Health Economic Evaluation). Oxford University Press; 2006. ISBN 9780198526629. https://openlibrary.org/isbn/9780198526629
- Jalal Hawre, Pechlivanoglou Petros, Krijkamp Eline, Alarid-Escudero Fernando, Enns Eva, Hunink M. G. Myriam. An Overview of R in Health Decision Sciences. Medical Decision Making. 2017;37(7):735-746. doi:10.1177/0272989x16686559. https://doi.org/10.1177/0272989x16686559
The series, in learning order
- Health Economic Modelling Series Guide: From Evidence to a Decision
- Choosing a Decision-Analytic Model: Let the Clinical Pathway Decide
- Auditable Model Inputs: Build a Parameter Table You Can Trace
- Finding Model Inputs in the Literature: Match the Estimate to the Question
- Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation
- Building a Decision Tree Step by Step: Follow Outcomes, Costs and Health
- Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time (this article)
- Budget Impact Analysis in Practice: Count People Before Counting Money
- ICER and Net Monetary Benefit Across Strategies: Find the Relevant Frontier
- Probabilistic Sensitivity Analysis: From Parameter Distributions to CEAC
- Model Validation and Verification: Test the Question, Equations and Code
- Expected Value of Perfect Information: What Could Better Evidence Be Worth?
Key takeaways
- Name the states first, then the chance of each move in one cycle; each row of chances must add up to 1.
- Move the whole group one cycle at a time: 1,000 healthy people become 800 healthy, 150 sick and 50 dead, and the total stays 1,000 people every year.
- Count each cycle at its middle: averaging year one's start and end gives 0.945 QALYs per person, not 0.890 from the end count.
- Write the timing rules beside the totals, because the same trace gives 2,787.786 to 3,149.605 QALYs at today's value, depending on how each year is counted.
Related in the wiki: Decision Trees, Markov Cohorts and Probabilistic Sensitivity Analysis: How a Model Turns Evidence into a Decision · QALY, DALY and the Utilities Behind Them: What Is Actually Being Measured
0