การประเมินความสมเหตุสมผลและตรวจสอบแบบจำลอง: ทดสอบคำถาม สมการ และโค้ด

# 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")

# 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")

# 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")

# 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")

# 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")

# 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")

# 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")

# 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")