Model Validation and Verification: Test the Question, Equations and Code
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
A hospital team must tell its drug committee whether to trust a colleague's model of 1,000 patients with a chronic illness. The model moves people between Well, Ill and Dead once a year for 4 years and reports 3,655.875 life-years. Two questions come first. Is it built right, so the arithmetic does what was intended? That is verification. Is it the right model for these patients and this decision? That is validation. Four checks with answers known in advance test the first question. Every row of chances adds up to 1, and 1,000 people stay 1,000. The second year worked by hand gives the same 640, 247.5 and 112.5 as the model. A weight of 1 for every living state gives a quality-weighted total equal to the life-years. With death switched off, nobody dies. A table filled the wrong way fails at once and shrinks the cohort from 1,000 to 409.6. Validation needs clinicians, sources and outside data.
A colleague's model, 1,000 patients and one number
A hospital's health technology assessment team must advise its drug committee next week. A colleague has built a model of 1,000 patients with a chronic illness, and it reports 3,655.875 life-years for the group over the next 4 years: the years lived, counting one year for each person alive at the start of each year. The committee will compare that total with the same model run for a new treatment. The numbers are illustrative, chosen to teach.
The model has three health states: Well, Ill and Dead. All 1,000 patients start in Well. Each year a person in Well stays well with probability 0.80, falls ill with 0.15 and dies with 0.05. A person in Ill stays ill with 0.85 and dies with 0.15, nobody leaves Dead, and nobody gets well again once ill, an assumption to check with clinicians.
This kind of model has a name. A Markov model follows a group of people between a few health states, one cycle at a time; here a cycle is 1 year. In the first cycle, of the 1,000 people in Well, 800 stay well, 150 fall ill and 50 die; the full lesson is in Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time.
Before anyone uses the 3,655.875 life-years, the team must check the model. This article corrects a common habit: treating a model that runs without an error message and gives a plausible number as a model that has been checked.
Two questions: built right, and the right model
The first question is about the equations and the code. If 1,000 people enter the model, 1,000 must still be counted in every year, alive or dead, and the first year must split them into 800, 150 and 50. Checking that the equations and code do what the builder intended is called verification.
The second question is about the patients and the world. A model can keep every person and still be wrong, because the death chance of 0.05 came from younger patients, or because ill patients really can recover. Judging whether the model is credible enough for this decision and these patients is called validation.
Verification can be finished inside the model, because each of its checks has an answer known in advance. Validation needs knowledge from outside the model: clinicians, the sources of the inputs, other models and data the model never saw. The figure below puts the two questions side by side for the colleague's model.
First check: nobody appears or vanishes
The colleague keeps the yearly chances in a small table, shown below. Each row is where a person is now, and each column is where that person is a year later. People in a model can only move between states; they cannot appear from nowhere or vanish.
The colleague's table of yearly chances
| Now (row), a year later (column) | Well | Ill | Dead | Row sum |
|---|---|---|---|---|
| Well | 0.80 | 0.15 | 0.05 | 1 |
| Ill | 0 | 0.85 | 0.15 | 1 |
| Dead | 0 | 0 | 1 | 1 |
Two parts with known answers
This first check has two parts with known answers. Each row must add up to 1, because a person in Well must end the year somewhere: 0.80 + 0.15 + 0.05 = 1, for Ill 0 + 0.85 + 0.15 = 1 and for Dead 0 + 0 + 1 = 1. The total of Well, Ill and Dead must also stay at 1,000 people at every year boundary, the start (boundary 0) and the end of each of the 4 years. After the first year the split is 800, 150 and 50; after year 4 it is 409.6, 337.2 and 253.2, worked out in the example below, and again adds up to 1,000 people.
The counts carry decimals because a cohort model tracks expected numbers, averages over many possible groups. The table of chances is called the transition matrix. The dead stay in the count because Dead is an absorbing state, a state nobody leaves once in it. Checking the totals is called a conservation check.
Check the rows and the total in R
R code mass-balance.R
# Illustrative colleague's model: 1,000 people in Well, Ill and Dead, four yearly cycles.
# 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
P <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (i in seq_len(cycles)) trace[i + 1, ] <- trace[i, ] %*% P
cat("Transition row sums:", rowSums(P), "\n")
cat("Mass at boundaries 0 to 4 (people):", rowSums(trace), "\n")
cat("Maximum conservation error (people):", max(abs(rowSums(trace) - sum(start))), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")Output of the run mass-balance.log
Transition row sums: 1 1 1
Mass at boundaries 0 to 4 (people): 1000 1000 1000 1000 1000
Maximum conservation error (people): 0
After year 4 Well Ill Dead (people): 409.6 337.2 253.2Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
Second check, worked by hand: the cohort year by year
A second route to the same numbers catches a table read the wrong way round. Write each arrow of the model as one multiplication and add up what arrives in each state.
-
First year
\[ 1000 \times (0.80,\ 0.15,\ 0.05) = (800,\ 150,\ 50) \]
Start with 1,000 people in Well. Well keeps 1,000 × 0.80 = 800 people, sends 1,000 × 0.15 = 150 to Ill and 1,000 × 0.05 = 50 to Dead.
-
Second year, Well
\[ 800 \times 0.80 = 640 \]
Only people already in Well can stay there, because nobody recovers from Ill: 800 × 0.80 = 640 people stay well.
-
Second year, Ill
\[ 120 + 127.5 = 247.5 \]
Ill receives 800 × 0.15 = 120 people from Well and keeps 150 × 0.85 = 127.5 of its own. The 127.5 is an expected count, so it keeps its decimal.
-
Second year, Dead
\[ 50 + 40 + 22.5 = 112.5 \]
Dead keeps its 50, gains 800 × 0.05 = 40 from Well and 150 × 0.15 = 22.5 from Ill.
-
Add up the second year
\[ 640 + 247.5 + 112.5 = 1000 \]
The hand route gives 640, 247.5 and 112.5 people, still 1,000 people in all, so nobody was lost in the second year either.
-
Third year, the same way
\[ 96 + 210.375 = 306.375 \]
Well keeps 640 × 0.80 = 512 people. Ill receives 640 × 0.15 = 96 and keeps 247.5 × 0.85 = 210.375, so it holds 306.375. Dead adds 640 × 0.05 = 32 and 247.5 × 0.15 = 37.125 to its 112.5, so it holds 181.625.
-
Fourth year
\[ 1000 - 409.6 - 337.2 = 253.2 \]
Well keeps 512 × 0.80 = 409.6 and Ill holds 512 × 0.15 + 306.375 × 0.85 = 337.2 after rounding. Dead holds the rest, 1,000 - 409.6 - 337.2 = 253.2 people.
Result: The matrix calculation gives the same counts, a gap of 0 people in the second year and the same 409.6, 337.2 and 253.2 after year 4. If the two routes disagreed, one of them would be reading the table the wrong way round.
The name: an independent calculation
Working a small part of the model out a second, separate way is called an independent calculation. It catches a matrix read in the wrong direction or a wrong index. It reuses the same inputs, so it cannot tell whether 0.15 is the right yearly chance of falling ill.
Check the second year by hand in R
R code expanded-cycle.R
# The second year written out by hand, one arrow at a time, then compared with the matrix.
# 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)
P <- matrix(transition, nrow = length(start), byrow = TRUE)
year_one <- start[1] * P[1, ]
ill_parts <- c(year_one[1] * P[1, 2], year_one[2] * P[2, 2])
dead_parts <- c(year_one[3], year_one[1] * P[1, 3], year_one[2] * P[2, 3])
year_two <- c(year_one[1] * P[1, 1], sum(ill_parts), sum(dead_parts))
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - as.vector(start %*% P %*% P))), "\n")Output of the run expanded-cycle.log
First year by hand Well Ill Dead (people): 800 150 50
Second-year Ill: from Well, stayed Ill (people): 120 127.5
Second-year Dead: already, from Well, from Ill (people): 50 40 22.5
Second year by hand Well Ill Dead (people): 640 247.5 112.5
Largest gap to the matrix (people): 0Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
Check the third and fourth years in R
R code later-years.R
# The third and fourth years the same way, starting from the second year's 640, 247.5 and 112.5.
# 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)
P <- matrix(transition, nrow = length(start), byrow = TRUE)
year_two <- as.vector(start %*% P %*% P)
ill_three <- c(year_two[1] * P[1, 2], year_two[2] * P[2, 2])
dead_three <- c(year_two[1] * P[1, 3], year_two[2] * P[2, 3])
year_three <- as.vector(year_two %*% P)
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", as.vector(year_three %*% P)), "\n")Output of the run later-years.log
Third-year Ill: from Well, stayed Ill (people): 96 210.375
Third-year new deaths: from Well, from Ill (people): 32 37.125
Third year Well Ill Dead (people): 512 306.375 181.625
Fourth year Well Ill Dead (people): 409.6 337.2 253.2Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
Third check: a case where the answer is known
Each person alive at the start of a year counts as one year lived. The people alive at the start of years 1 to 4 are 1,000, 950, 887.5 and 818.375 (the 950 is 800 + 150, the 887.5 is 640 + 247.5 and the 818.375 is 512 + 306.375), so the total is 3,655.875 life-years.
Economic models usually weight each year by its quality, because a year ill is worth less to a patient than a year well. 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.
Now give Well and Ill a weight of 1, and Dead a weight of 0. Every year alive then counts in full, so the QALY total must equal the 3,655.875 life-years exactly. The two totals are 3,655.875 and 3,655.875, a gap of 0. A setting like this, where the right answer is known before the model runs, is called a limiting case.
Compare life-years and weight-1 QALYs in R
R code utility-one.R
# Count life-years, then weighted years with weight 1 for Well and Ill: the two totals must match.
# Each yearly cycle counts one year for every person alive at its start.
# 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
weight_one <- c(1, 1, 0)
P <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (i in seq_len(cycles)) trace[i + 1, ] <- trace[i, ] %*% P
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(sum(alive_start), weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - sum(alive_start), "\n")Output of the run utility-one.log
Alive at the start of years 1 to 4 (people): 1000 950 887.5 818.375
Life-years, weight-1 years: 3655.875 3655.875
Weighted minus life-years: 0Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
Count both totals the same way
The equality holds only when both totals are counted the same way, and time is the easiest place to break that.
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.
Check the two recaps in R
R code weights.R
# The two recaps the weight-1 check leans on: a weighted year, and a future amount valued today.
# Writes: nothing to disk; the results print in the console.
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")Output of the run weights.log
10 years at weight 0.70, gain: 7.0 3.0
Yearly rate back to today (%): 3
1,834 THB due in year 3, worth today (THB): 1,678Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
A mismatch the check catches
Suppose the colleague discounts the QALYs at 3% a year but leaves the life-years as counted. Each year's count is then divided by 1.03 once for each year away: 1,000 / 1.03 = 970.874, 950 / 1.03² = 895.466, 887.5 / 1.03³ = 812.188 and 818.375 / 1.03⁴ = 727.116. These add up to 3,405.644, so the check fails by 3,655.875 - 3,405.644 = 250.231 even though every weight is right. The cause is the mismatched treatment of time, so treat both totals the same way; the weight-1 setting itself is a test of the code, not a claim that ill patients have full health.
Plant the time mismatch in R
R code one-side.R
# Plant a mismatch: shrink the weight-1 years for time at 3% a year, leave the life-years as counted.
# 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
yearly_factor <- 1.03
P <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (i in seq_len(cycles)) trace[i + 1, ] <- trace[i, ] %*% P
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
each_year <- alive_start / yearly_factor^seq_len(cycles)
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
shrunk <- sum(each_year)
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, sum(alive_start) - shrunk)), "\n")Output of the run one-side.log
Each year shrunk (years): 970.874 895.466 812.188 727.116
Shrunk weight-1 years, gap: 3405.644 250.231Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
Fourth check: switch death off
If nobody can die, all 1,000 people must still be alive after 4 years. To switch death off, move each death chance back into staying put: the Well row becomes 0.85, 0.15 and 0, and the Ill row becomes 0, 1 and 0. The Dead row stays as it is.
The living count then stays at 1,000 people in every year and the deaths stay at 0. Simply deleting the death chances would be a different error: the Well row would no longer add up to 1, and people would vanish.
This is a limiting case too. It catches deaths sent to the wrong state, or people sent out of Dead, but it cannot say whether 0.05 is the right yearly death chance for these patients.
Switch death off in R
R code zero-mortality.R
# 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
P_zero <- matrix(transition, nrow = length(start), byrow = TRUE)
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (i in seq_len(cycles)) trace_zero[i + 1, ] <- trace_zero[i, ] %*% P_zero
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")Output of the run zero-mortality.log
Zero-mortality living people: 1000 1000 1000 1000 1000
Zero-mortality deaths: 0 0 0 0 0Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
A broken model fails the checks
Now break the colleague's model on purpose. The nine chances are typed in as one list, row after row: 0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1. Fill the table column by column instead, and the from and to states swap: the rows become 0.80, 0 and 0 (sum 0.8), then 0.15, 0.85 and 0 (sum 1), then 0.05, 0.15 and 1 (sum 1.2).
Well now keeps 0.80 of its people each year and sends nobody anywhere, so the first year leaves 1,000 × 0.80 = 800 people, all in Well. The cohort then shrinks to 800 × 0.80 = 640, then 512 and 409.6 people over four years. The code runs without an error message, yet the row and conservation checks fail at once.
The 409.6 matches the Well count of the intended model after year 4 only by coincidence, because in both Well keeps 0.80 of its people a year: 1,000 × 0.80⁴. Here, though, it is everyone left, because everyone else has vanished.
Break the matrix on purpose in R
R code reversed-matrix.R
# Break the model on purpose: fill the matrix by column instead of by row.
# Rows should be "from" states and columns "to" states; byrow = FALSE swaps them.
start <- c(1000,0,0)
transition <- c(0.80,0.15,0.05,0,0.85,0.15,0,0,1)
cycles <- 4
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (i in seq_len(cycles)) trace_bad[i + 1, ] <- trace_bad[i, ] %*% P_reversed
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")Output of the run reversed-matrix.log
Reversed-matrix row sums: 0.8 1 1.2
Reversed-matrix mass at boundaries 0 to 4 (people): 1000 800 640 512 409.6
Reversed-matrix first cycle Well Ill Dead (people): 800 0 0Requires R ≥ 4.5.2 · tibble 3.3.1 · dplyr 1.2.1 · works in R 4.6.1 · Download model-validation-and-verification.R
Full tidy script
# Model Validation and Verification: Test the Question, Equations and Code
# Check a colleague's cohort model of 1,000 people before a committee uses its life-years.
# ตรวจแบบจำลองกลุ่มประชากร 1,000 คนของเพื่อนร่วมงาน ก่อนที่คณะกรรมการจะใช้ปีชีวิตจากแบบจำลองนั้น
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# Writes: no files; console results only.
# tibble 3.3.1 (stable)
# dplyr 1.2.1 (stable)
# install.packages(c("tibble", "dplyr"))
# update.packages(oldPkgs = c("tibble", "dplyr"))
# remotes::install_version("tibble", "3.3.1")
# remotes::install_version("dplyr", "1.2.1")
stopifnot(getRversion() >= "4.5.2")
library(tibble) # Keep the yearly counts in a named table.
library(dplyr) # Add the totals and the conservation gap to each row.
stopifnot(as.character(utils::packageVersion("tibble")) == "3.3.1")
stopifnot(as.character(utils::packageVersion("dplyr")) == "1.2.1")
#%%@note 1. Write the illustrative inputs once
# States are Well, Ill and Dead. Rows are where people are now, columns where they go next.
# สถานะคือ Well, Ill และ Dead แถวคือที่ที่คนอยู่ตอนนี้ คอลัมน์คือที่ที่จะไปต่อ
start <- c(1000, 0, 0)
transition <- c(0.80, 0.15, 0.05, 0, 0.85, 0.15, 0, 0, 1)
cycles <- 4
weight_one <- c(1, 1, 0)
years_in_state <- 10
state_weight <- 0.70
discount_rate <- 0.03
yearly_factor <- 1.03
one_decimal <- function(value) sprintf("%.1f", value)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
#%%@note 2. Follow the cohort and check that nobody appears or vanishes
# Each row of the table is one year boundary; the gap must stay at zero people.
# แต่ละแถวของตารางคือขอบปีหนึ่ง ผลต่างต้องเป็นศูนย์คนตลอด
transition_matrix <- matrix(transition, nrow = length(start), byrow = TRUE)
trace <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace[cycle_index + 1, ] <- trace[cycle_index, ] %*% transition_matrix
}
validation_table <- tibble::tibble(
boundary = seq_len(nrow(trace)) - 1, # Cohort model: year boundary.
well = trace[, 1], # Cohort model: people in Well.
ill = trace[, 2], # Cohort model: people in Ill.
dead = trace[, 3] # Cohort model: people in Dead.
) |>
dplyr::mutate(
total = well + ill + dead, # Cohort model: everyone, the dead included.
error = total - sum(start) # Cohort model: people gained or lost.
)
stopifnot(max(abs(validation_table$error)) < 1e-9)
#%%@note 3. Write the second year out by hand, then the third
# One product per arrow, so a flipped matrix would give different counts.
# หนึ่งผลคูณต่อหนึ่งลูกศร ถ้าเมทริกซ์กลับด้านจะได้จำนวนคนต่างไป
year_one <- start[1] * transition_matrix[1, ]
ill_parts <- c(
year_one[1] * transition_matrix[1, 2],
year_one[2] * transition_matrix[2, 2]
)
dead_parts <- c(
year_one[3],
year_one[1] * transition_matrix[1, 3],
year_one[2] * transition_matrix[2, 3]
)
year_two <- c(year_one[1] * transition_matrix[1, 1], sum(ill_parts), sum(dead_parts))
ill_three <- c(year_two[1] * transition_matrix[1, 2], year_two[2] * transition_matrix[2, 2])
dead_three <- c(year_two[1] * transition_matrix[1, 3], year_two[2] * transition_matrix[2, 3])
year_three <- as.vector(year_two %*% transition_matrix)
year_four <- as.vector(year_three %*% transition_matrix)
#%%@note 4. Count life-years and weight-1 years on the same clock
# Each yearly cycle counts one year for every person alive at its start.
# คนที่ยังมีชีวิตเมื่อต้นปีนับคนละหนึ่งปี
alive_start <- rowSums(trace[seq_len(cycles), 1:2])
life_years <- sum(alive_start)
weighted_years <- sum(trace[seq_len(cycles), ] %*% weight_one)
quality_years <- years_in_state * state_weight
gain <- years_in_state - quality_years
stopifnot(abs(yearly_factor - (1 + discount_rate)) < 1e-12)
each_year <- alive_start / yearly_factor^seq_len(cycles)
shrunk <- sum(each_year)
#%%@note 5. Switch death off without losing anyone
# Each living state's death chance goes back to staying in that state.
# โอกาสเสียชีวิตของแต่ละสถานะที่มีชีวิตย้ายกลับไปเป็นโอกาสอยู่สถานะเดิม
P_zero <- transition_matrix
alive <- seq_len(length(start) - 1)
P_zero[alive, alive] <- P_zero[alive, alive] + diag(P_zero[alive, length(start)])
P_zero[alive, length(start)] <- 0
trace_zero <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_zero[cycle_index + 1, ] <- trace_zero[cycle_index, ] %*% P_zero
}
#%%@note 6. Break the model on purpose
# Filling the matrix by column swaps the from and to states.
# การเติมเมทริกซ์ตามคอลัมน์สลับสถานะต้นทางกับปลายทาง
P_reversed <- matrix(transition, nrow = length(start), byrow = FALSE)
trace_bad <- matrix(rep(start, cycles + 1), ncol = length(start), byrow = TRUE)
for (cycle_index in seq_len(cycles)) {
trace_bad[cycle_index + 1, ] <- trace_bad[cycle_index, ] %*% P_reversed
}
#%%@note 7. Numbers the article quotes
# Print in the order of the short boxes in the article.
# พิมพ์ตามลำดับกล่องสั้นในบทความ
cat("Transition row sums:", rowSums(transition_matrix), "\n")
cat("Mass at boundaries 0 to 4 (people):", validation_table$total, "\n")
cat("Maximum conservation error (people):", max(abs(validation_table$error)), "\n")
cat("After year 4 Well Ill Dead (people):", sprintf("%.1f", trace[cycles + 1, ]), "\n")
cat("First year by hand Well Ill Dead (people):", year_one, "\n")
cat("Second-year Ill: from Well, stayed Ill (people):", ill_parts, "\n")
cat("Second-year Dead: already, from Well, from Ill (people):", dead_parts, "\n")
cat("Second year by hand Well Ill Dead (people):", year_two, "\n")
cat("Largest gap to the matrix (people):", max(abs(year_two - trace[3, ])), "\n")
cat("Third-year Ill: from Well, stayed Ill (people):", ill_three, "\n")
cat("Third-year new deaths: from Well, from Ill (people):", dead_three, "\n")
cat("Third year Well Ill Dead (people):", year_three, "\n")
cat("Fourth year Well Ill Dead (people):", sprintf("%.1f", year_four), "\n")
cat("Alive at the start of years 1 to 4 (people):", alive_start, "\n")
cat("Life-years, weight-1 years:", sprintf("%.3f", c(life_years, weighted_years)), "\n")
cat("Weighted minus life-years:", weighted_years - life_years, "\n")
cat("10 years at weight 0.70, gain:", one_decimal(c(quality_years, gain)), "\n")
cat("Yearly rate back to today (%):", 100 * discount_rate, "\n")
cat("1,834 THB due in year 3, worth today (THB):", big(1834 / (1 + discount_rate)^3), "\n")
cat("Each year shrunk (years):", sprintf("%.3f", each_year), "\n")
cat("Shrunk weight-1 years, gap:", sprintf("%.3f", c(shrunk, life_years - shrunk)), "\n")
cat("Zero-mortality living people:", rowSums(trace_zero[, alive]), "\n")
cat("Zero-mortality deaths:", trace_zero[, length(start)], "\n")
cat("Reversed-matrix row sums:", rowSums(P_reversed), "\n")
cat("Reversed-matrix mass at boundaries 0 to 4 (people):", rowSums(trace_bad), "\n")
cat("Reversed-matrix first cycle Well Ill Dead (people):", trace_bad[2, ], "\n")
Validation: is it the right model for these patients?
All four checks pass, and none of them says the model fits the hospital's patients. Start with the pathway: in this model nobody recovers from Ill. Show the three states to the clinicians and ask whether a real patient could take a path the model forbids. Whether the structure matches the illness is called conceptual validity.
Then check each input against its source. The 0.05 must be a yearly chance of death for patients like these, not a risk over a longer period or a figure from a younger trial group. Whether the numbers fit the population, the time unit and the question is called input validity; the series article on finding inputs in the literature, Finding Model Inputs in the Literature: Match the Estimate to the Question, shows how to trace such a number to its source.
This model predicts 253.2 deaths per 1,000 people by year 4, and two further checks need something outside the model to test it. Cross-model comparison runs another team's model with the same inputs: if it finds far more or far fewer deaths, one of the two models has an error to find. External validation sets the 253.2 against data the model never used, such as a registry of similar patients. Predictive validation does the same with data collected after the model was built, such as patients followed from then on.
Two published checklists help. AdViSHE, from Vemer and colleagues, lists validation questions for a model's users and decision makers across these areas. TECH-VER, from Büyükkaramikli and colleagues, lists technical verification tests like the four above. Neither makes a model credible by being ticked; use each as a record of what was checked and what was found.
Eddy and colleagues link validation with transparency: a model others can rerun and inspect is easier to check, but being reproducible does not make its predictions accurate. Jalal and colleagues describe R as a good home for such open models.
The validation record for the colleague's model
| Check | Expected before running | Result found | What stays open |
|---|---|---|---|
| Rows and total | Each row sums to 1; 1,000 people at every year boundary | 1, 1 and 1; 1,000 at boundaries 0 to 4; gap 0 | Whether 0.05 is the right death chance |
| Years by hand | One product per arrow gives the matrix's counts | 640, 247.5 and 112.5 both ways; gap 0 | The hand route reuses the same inputs |
| Weight 1 | QALYs equal life-years when both are counted the same way | 3,655.875 and 3,655.875; gap 0 | The real quality weights of ill patients |
| No deaths | 1,000 alive and 0 dead in every year | 1,000 alive and 0 dead at boundaries 0 to 4 | The real yearly death chances |
| Outside data | Model: 253.2 deaths per 1,000 people by year 4 | No outside data yet | Open until a registry of similar patients is compared |
Keep a validation record that can stay open
For each check, write down the claim it tests, the expected result, the actual result, where the evidence is and what remains open. Write the expected result before running the check, so a surprising result cannot quietly become the new expectation. Add the model version and the reviewer's name.
The four verification checks are closed. The external comparison stays open, because no registry data are in hand to set against the 253.2 deaths.
A failed check is a finding to explain, not a reason to adjust inputs until it passes. Record its cause, its fix and what it changed in the result. The team can then tell the committee what the 3,655.875 life-years have been tested for, and what they have not.
Terms to keep distinct
- Verification (การตรวจสอบแบบจำลอง)
- Checking that the equations and code do what the builder intended; here, that 1,000 people stay 1,000.
- Validation (การประเมินความสมเหตุสมผล)
- Judging whether the model is credible enough for this decision and these patients; here, whether 0.05 fits them.
- Transition matrix (เมทริกซ์การย้ายสถานะ)
- The table of yearly chances, one row for where a person is now and one column for where they go; each row sums to 1.
- Absorbing state (สถานะดูดซับ)
- A state nobody leaves once in it; Dead in this model.
- Independent calculation (การคำนวณซ้ำอีกทางหนึ่ง)
- Working part of the model out a second, separate way; here, the years by hand against the matrix.
- Limiting case (กรณีขอบเขต)
- A setting where the right answer is known in advance, such as weight 1 (3,655.875 on both sides) or no deaths.
- Cross-model comparison (การเปรียบเทียบข้ามแบบจำลอง)
- Running a model another team built with the same inputs and comparing the answers.
- External validation (การประเมินความสมเหตุสมผลภายนอก)
- Comparing a prediction, such as 253.2 deaths per 1,000 people by year 4, with data the model never used.
Connect to the clinical foundations
How a cohort model like the colleague's turns evidence into a decision is the subject of the foundation article on decision trees and cohort models. How sure the evidence behind an input such as the 0.05 is, and how to choose it, is the subject of the foundation article on GRADE certainty and choosing effectiveness inputs.
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
- Eddy David M., Hollingworth William, Caro J. Jaime, Tsevat Joel, McDonald Kathryn M., Wong John B.. Model Transparency and Validation: A Report of the ISPOR-SMDM Modeling Good Research Practices Task Force-7. Value in Health. 2012;15(6):843-850. doi:10.1016/j.jval.2012.04.012. https://api.crossref.org/works/10.1016/j.jval.2012.04.012
- Vemer P., Corro Ramos I., van Voorn G. A. K., Al M. J., Feenstra T. L.. AdViSHE: A Validation-Assessment Tool of Health-Economic Models for Decision Makers and Model Users. PharmacoEconomics. 2016;34(4):349-361. doi:10.1007/s40273-015-0327-2. https://api.crossref.org/works/10.1007/s40273-015-0327-2
- Büyükkaramikli Nasuh C., Rutten-van Mölken Maureen P. M. H., Severens Johan L., Al Maiwenn. TECH-VER: A Verification Checklist to Reduce Errors in Models and Improve Their Credibility. PharmacoEconomics. 2019;37(11):1391-1408. doi:10.1007/s40273-019-00844-y. https://api.crossref.org/works/10.1007/s40273-019-00844-y
- 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://api.crossref.org/works/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
- 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 (this article)
- Expected Value of Perfect Information: What Could Better Evidence Be Worth?
Key takeaways
- Ask two questions before using a model's number: is it built right (verification), and is it the right model for these patients (validation)?
- Verify with answers known in advance: every row sums to 1, the 1,000 people stay 1,000, and the years worked by hand give the same 640, 247.5 and 112.5 as the matrix.
- With a weight of 1 for living states, QALYs must equal the 3,655.875 life-years; shrink only one side for time and the check fails by 250.231.
- A table filled by column runs without an error message yet shrinks the cohort from 1,000 to 409.6: only the checks catch it.
- Validation needs clinicians, sources and outside data; keep it open in the record until the 253.2 predicted deaths are compared.
Related in the wiki: Decision Trees, Markov Cohorts and Probabilistic Sensitivity Analysis: How a Model Turns Evidence into a Decision · Beyond the Evidence Pyramid: GRADE Certainty and Choosing Effectiveness Inputs for Economic Models
0