Probabilistic Sensitivity Analysis: From Parameter Distributions to CEAC

อ่านฉบับภาษาไทย (Thai version)

Abstract

A health fund must decide whether to pay for option C. The best estimates say C costs 150,000 THB more per patient than usual care A and adds 1.00 year in full health. The fund pays up to 160,000 THB for such a year, so C comes out 10,000 THB ahead. But the trial behind C was small, so its risk ratio, 0.56, could lie anywhere from 0.25 to 1.26, including harm. This article carries that doubt into the decision. Each uncertain input gets a curve of plausible values that suits it, whether a ratio, a rate, a cost or a row of chances. The model then runs 1,000 times with fresh draws. C should come out ahead, with a net gain above 0, in about 599 of them: a probability of 0.60, give or take 3.0 percentage points. Repeating the count at every price gives a curve the fund can read. Short R boxes reproduce every number.

Three steps for the fund's choice between option C and usual care A. First, each uncertain input gets a curve: the risk ratio 0.56 (0.25 to 1.26) never falls below 0, and the response rate 0.30 stays between 0 and 1 on its curve. Then the model runs 1,000 times, and C should come out ahead in about 599 of them. Last, the count is repeated at every price, giving 0.50 at 150,000 THB and 0.60 at 160,000 THB. — Read the three panels from the top: each input becomes a curve, the curves feed 1,000 runs, and the runs give the share for C at every price, so the fund sees how sure it can be. Every number is worked out in the sections below.

Is option C worth 150,000 THB more?

A health fund's committee must decide whether to pay for option C, a new treatment, instead of usual care A. The fund's model gives its best estimates: C costs 150,000 THB more per patient and adds 1.00 year in full health. The fund pays up to 160,000 THB for one such year; we call this ceiling the fund's price, and later we vary it. Valued that way, C brings 160,000 × 1.00 = 160,000 THB of health for 150,000 THB of extra cost, a net gain of 10,000 THB per patient.

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. In these units C adds 1.00 QALY for 150,000 THB more, and the fund pays up to 160,000 THB per QALY.

The trouble is that every input of the model is an estimate. The trial behind C counted 8 patients with the bad outcome (an event) among 66 treated patients and 13 among 60 controls, so the risks were 12.12% and 21.67% in the two groups. That gives a risk ratio of 0.56, but its 95% confidence interval runs from 0.25 to 1.26, and a ratio above 1 would mean C does harm.

Other inputs are uncertain too, and their numbers here are illustrative. The share of patients who respond to C is estimated at 0.30, and a course of care costs 10,000 THB on average. Of 100 healthy patients followed for a year, 80 stayed healthy, 15 fell sick and 5 died. The committee asks how sure it can be that C is worth 150,000 THB more.

This article answers in three steps. First, give each uncertain input a curve that fits what it is. Second, run the model many times, each time with fresh values from every curve. Last, count how often C comes out ahead, meaning a net gain above 0, at 160,000 THB and at every other price.

The case in numbers

R code case.R

# The fund's question in numbers: option C against usual care A, per patient.
# Writes: nothing to disk; the results print in the console.
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")

Output of the run case.log

Extra cost THB: 150,000 | extra health: 1.00 | price THB: 160,000 | net gain THB: 10,000
Trial risks: 12.12% 21.67% | risk ratio: 0.56
Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health: 7.0 3.0

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

One estimate, many plausible values

The response rate of 0.30, with a standard error of 0.05, is not one value but a range of values that could be true, a little above or below 0.30, some more plausible than others. Drawn as a curve, the range is high where values are plausible and low where they are not. Statisticians call such a curve a distribution. Each run of the model will use one value drawn from each input's curve.

The curve describes how sure we are about an average, not how much patients differ. The standard error of 0.05 measures how far the estimated rate could be from the true rate. Each patient still either responds or does not; the curve is about the rate.

The shape must fit what the input is. A ratio cannot fall below 0, a rate must stay between 0 and 1, a cost cannot be negative, and a row of chances must add up to 1 in every row. The next four sections build one curve of each kind from the case's numbers.

The risk ratio: a symmetric bell on the log scale

Start with the trial's risk ratio, 0.56 with a 95% interval from 0.25 to 1.26. The interval is lopsided: 0.25 lies 0.31 below 0.56, but 1.26 lies 0.70 above it. A symmetric curve around 0.56 with that width would also give some values below 0, which no ratio can take.

Taking natural logarithms fixes both problems. The logs of 0.25, 0.56 and 1.26 are -1.386, -0.580 and 0.231, so on this scale the interval is even on both sides of its centre. The plausible values form a symmetric bell: values near the centre are common, values far from it are rare. In such a bell, 95% of the values lie within 1.96 spreads of the centre.

That rule gives the spread. On this scale the interval runs from -1.386 to 0.231, and half of that width is 0.809. Dividing 0.809 by 1.96 gives a spread of 0.413 on the log scale. Check it the other way round: -0.580 − 1.96 × 0.413 = -1.389 and -0.580 + 1.96 × 0.413 = 0.229, which round to -1.39 and 0.23.

This symmetric bell is called a Normal distribution. About 95% of its values fall within 1.96 standard deviations of its mean. Here the mean is -0.58 and the standard deviation is 0.41, the spread just found, so the middle 95% runs from -1.39 to 0.23.

Now convert each value back to a ratio by raising e to its power. The ends become 0.249 and 1.257, the trial's 0.25 and 1.26 once rounded. Every converted value is above 0, the values bunch below 1 and trail off past 1.26, and 8.0% of them lie above 1, where C does harm.

This curve is called a lognormal distribution: the curve of a quantity whose log follows a Normal distribution. It is the usual choice for an uncertain ratio, because it never goes below 0 and it reproduces a lopsided interval.

⤢ click to enlarge

Look at the left panel first: on the log scale the bell is even on both sides of its centre, -0.58, at the dashed line. The right panel shows the same values as ratios, so the curve starts at 0, peaks below 1 and has a long right side toward 1.26. The 0.56 marks the middle value of the curve, and the peak sits a little to its left.

From the trial's interval to the curve

R code log-scale.R

# The trial's risk ratio and its 95% interval, moved to the log scale and back.
# Writes: nothing to disk; the results print in the console.
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * pnorm(logs[2] / spread)), "\n")

Output of the run log-scale.log

Below | above 0.56: 0.31 0.70 | logs: -1.386 -0.580 0.231
Half width: 0.809 | / 1.96, the spread: 0.413 | rounded centre and spread: -0.58 0.41
Centre -/+ 1.96 spreads: -1.389 0.229 | rounded: -1.39 0.23
Back as ratios: 0.249 1.257 | rounded: 0.25 1.26
Draws above 1 (%): 8.0

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

The response rate: a curve that stays between 0 and 1

The response rate is a probability, so every plausible value must lie between 0 and 1. A bell would work near 0.30, but a rate of 0.05 with the same standard error would spill below 0. The curve needs walls at 0 and at 1.

Think of the evidence as a number of patients' worth of information. A rate seen in n patients has a standard error whose square is the rate × (1 − rate) divided by n. For our rate, 0.30 × 0.70 = 0.21 (0.70 is 1 − 0.30) and 0.05 × 0.05 = 0.0025, so n = 0.21 / 0.0025 = 84 patients.

The curve used here divides by its two counts plus one, so the two counts add up to one fewer: 84 − 1 = 83 patients' worth. Of those, 0.30 × 83 = 24.9 responded and 0.70 × 83 = 58.1 did not, and 24.9 / 83 gives back the rate of 0.30 we started from.

A curve built from two such counts is called a Beta distribution. Here its two shapes, 24.9 and 58.1, are those counts. It is a hump that never leaves the range from 0 to 1, so it fits a risk, a response rate or a utility, and here its middle 95% runs from 0.21 to 0.40; that 0.21 matches 0.30 × 0.70 above only by coincidence.

The cost of a course: a curve that stays above zero

A mean cost cannot be negative, and a few expensive patients can pull it up more than cheap ones can pull it down. Its curve should start at 0 and have a long right tail. The mean cost of a course here is 10,000 THB, with a standard error of 2,000 THB.

Think of the cost of a course as 25 small pieces of about 400 THB each. Together they cost 25 × 400 = 10,000 THB, our mean. Because each piece varies, the total varies by the square root of 25, which is 5, times 400 THB, giving 2,000 THB, our standard error. So the two numbers 25 and 400 give back both of our estimates.

To find them from the mean and the standard error, work backwards. The first number, the shape, is the mean divided by the standard error, squared: 10,000 / 2,000 = 5, and 5 × 5 gives 25 as the shape. The second, the scale, is the standard error squared, divided by the mean: 2,000 × 2,000 = 4,000,000, and 4,000,000 / 10,000 = 400 THB.

This curve with a long right tail is called a Gamma distribution, here with shape 25 and scale 400 THB. Its middle 95% runs from 6,471 to 14,284 THB: 3,529 THB below the mean but 4,284 THB above it, the longer right side. Its peak sits at 9,600 THB, a little below the mean.

⤢ click to enlarge

Compare the two panels. On the left, the rate's curve peaks near 0.30 and never leaves the range from 0 to 1 at either end. On the right, the cost's curve starts at 0 and reaches further above its mean of 10,000 THB than below it, so a few dear courses pull the average up.

From a rate and its standard error to two shapes

R code beta.R

# The response rate: an estimate of 0.30 with a standard error of 0.05.
# Writes: nothing to disk; the results print in the console.
rate_mean <- 0.30
rate_se <- 0.05
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")

Output of the run beta.log

0.30 x 0.70: 0.21 | 0.05 x 0.05: 0.0025
Divided: 84 | minus 1: 83
Shapes, responders | non-responders: 24.9 58.1 | their mean: 0.3
Middle 95%: 0.21 0.4

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

From a cost and its standard error to shape and scale

R code gamma.R

# The mean cost of a course of care: 10,000 THB with a standard error of 2,000 THB.
# Writes: nothing to disk; the results print in the console.
cost_mean <- 10000
cost_se <- 2000
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
cost_ends <- qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")

Output of the run gamma.log

Mean / standard error: 5 | squared, the shape: 25
2,000 x 2,000: 4,000,000 | / 10,000, the scale THB: 400 | shape x scale: 10,000
Middle 95% THB: 6,471 14,284 | peak THB: 9,600
Below | above mean THB: 3,529 4,284

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

A row of chances that must add up to 1

The fund's model also follows patients from year to year. A Markov model follows a group of people between a few health states, one cycle at a time. Of 1,000 healthy people, 800 stay healthy, 150 fall sick and 50 die in the first cycle; the full lesson is in Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time. Those moves come from one row of chances, 0.80, 0.15 and 0.05, counted as 80, 15 and 5 of 100 healthy patients.

The three chances share out everyone in the healthy state, so they must add up to 1. Suppose each chance is drawn from its own Beta curve, built from its count (80, 15 or 5) and the rest of the 100 patients. Nothing links the three draws, so their sum wanders: in three rows drawn that way it came out at 0.911, 0.940 and 1.115. The first two rows lose patients, and the last one creates them.

The fix is to draw one positive number for each count, 80, 15 and 5, from a Gamma curve whose shape is that count, and divide each by the total of the three. Three rows drawn that way came out as 0.819, 0.149 and 0.033; 0.842, 0.078 and 0.080; and 0.805, 0.166 and 0.029 for healthy, sick and dead. The three shares split one total, so each row sums to exactly 1 before the digits are rounded, and over many rows they centre on 0.80, 0.15 and 0.05, the counted shares.

Drawing a whole row this way is called a Dirichlet distribution. It does for three or more outcomes what the Beta distribution does for two: with counts 80, 15 and 5, its rows average 0.80, 0.15 and 0.05, and every row sums to 1.

⤢ click to enlarge

Compare the two sides: all three bars on the left reach the dashed line at 1, while on the right the first two bars stop at 0.911 and 0.940 and the third overshoots to 1.115. A model fed the right-hand rows would lose or invent patients.

R code row.R

# Yearly chances for a healthy patient: of 100, 80 stay healthy, 15 fall sick and 5 die.
# Writes: nothing to disk; the results print in the console.
counts <- c(80, 15, 5)
cohort <- 1000
set.seed(20261028)
together <- matrix(rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")

Output of the run row.log

Average row: 0.80 0.15 0.05 | 1,000 healthy, one year on: 800 150 50
Rows drawn together: 0.819 0.149 0.033 0.842 0.078 0.080 0.805 0.166 0.029 | sums: 1 1 1
Rows drawn one share at a time: 0.763 0.112 0.036 0.730 0.138 0.072 0.819 0.208 0.088
Their sums: 0.911 0.940 1.115

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

Run the model many times

The committee's question needs every input to be uncertain at once, because a poor ratio (near or above 1) and a dear course can strike in the same run. So the model is run many times. In each run it draws one value from every curve, works out C and A from that same draw, and records C's extra QALY and extra cost.

In the fund's full model, all four curves would feed each run, moving C's health, its cost and the patients' yearly steps. Building that whole model is a lesson of its own, so here we use only the summary its runs would give. Across runs, C's extra QALY averages 1.00 with a standard deviation of 0.20, and its extra cost averages 150,000 THB with a standard deviation of 24,000 THB. These two spreads are illustrative: the four curves taught the shapes, and the arithmetic from here on needs only these two summaries.

In each run, C's net gain at 160,000 THB per QALY is 160,000 × its extra QALY, minus its extra cost. The worked example below shows where the average net gain and its spread come from, and how often the net gain is above 0.

Worked example: how often is C ahead?

Each run gives C an extra QALY drawn around 1.00 (spread 0.20) and an extra cost drawn around 150,000 THB (spread 24,000 THB). Its net gain is 160,000 × extra QALY, minus extra cost.

The average net gain

\[ 160000 \times 1.00 - 150000 = 10000 \]

Put the averages in: on average C is 10,000 THB ahead.
The spread from each side

\[ 160000 \times 0.20 = 32000 \]

The health side varies by 32,000 THB and the cost side by 24,000 THB.
Spreads add as squares

\[ 32 \times 32 + 24 \times 24 = 1024 + 576 = 1600 = 40 \times 40 \]

When the two sides vary independently, a high health draw is sometimes met by a high cost draw, so the swings partly cancel and the total swings less than the two spreads added. Their squares add instead. The formula works in thousands of THB (32 is 32,000 THB and 24 is 24,000 THB), so its 40 is a spread of 40,000 THB for the net gain.
How far above zero the average sits

\[ 10000 / 40000 = 0.25 \]

The average sits a quarter of a spread above zero.
The share above zero

\[ 0.599 \times 1000 = 599 \]

The net gain, the difference of two bells, is itself a Normal bell. Zero sits 0.25 spreads below its average, so the share of values above zero is 0.599, read from a Normal table or with pnorm(0.25) in R.

Result: C is ahead in about 599 of 1,000 runs: a probability of 0.60 that it is worth its extra cost at 160,000 THB per QALY.

The name for it

Rerunning a model many times, each time with every uncertain input drawn from its curve, is called probabilistic sensitivity analysis (PSA). Its answer is a share of runs rather than a single best result: here the bell expects about 599 of 1,000 runs to favour C. The cloud of runs below shows the same answer as dots.

⤢ click to enlarge

Each dot is one run: across is its extra QALY, up is its extra cost. The line marks where the extra cost equals 160,000 THB × the extra QALY, so green dots below it are runs where C's extra health is worth more than its extra cost, and grey dots above it are runs where A is ahead. A little more than half of the dots are green, the 0.60 of the worked example.

R code net-gain.R

# Net gain of C at 160,000 THB per year of full health, with both of its results uncertain.
# Writes: nothing to disk; the results print in the console.
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- pnorm(gain_mean / gain_sd)
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")

Output of the run net-gain.log

Average net gain THB: 10000 | spread parts THB: 32000 24000
In thousands: 32 24 | squared: 1024 576 | sum: 1600 | root: 40 | spread THB: 40000
Average / spread: 0.25 | share above zero: 0.599 | rounded: 0.6 | runs of 1,000: 599

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

How many runs are enough?

The bell gave 599 runs out of 1,000, but a real set of runs is random. Our own 1,000 runs, drawn from a fixed starting point so anyone can repeat them, put C ahead in 592 runs. Their net gains averaged 9,236 THB with a spread of 39,057 THB, close to the 10,000 THB and 40,000 THB of the worked example.

Early on, the count jumps around. The share was 0.400 after 10 runs, 0.510 after 100 runs and 0.592 after 1,000 runs. More runs make it settle, but some wobble always remains.

The size of that wobble follows from the share and the number of runs. For a share of 0.60, the product 0.60 × 0.40 = 0.24 measures how mixed the runs are. Divide it by 1,000 runs and take the square root, and the standard error of the share is 0.0155 here. Multiplying by 1.96 gives about 3.0 percentage points, so another set of 1,000 runs would most likely land between 0.57 and 0.63 as well.

Answering a question by drawing random inputs many times and counting the results is called Monte Carlo simulation. The wobble that is left is called the Monte Carlo error: here the share's standard error is 0.0155, which moves it about 3.0 percentage points either way at 60%, and more runs shrink it. Whole decision models can be run this way in R (Jalal and colleagues).

⤢ click to enlarge

Follow the blue line from left to right: the share jumps from 0.400 after 10 runs to 0.510 after 100 runs, then settles near the dashed 0.599 line, inside the green band of 0.57 to 0.63 that marks the leftover wobble. So 1,000 runs give an answer good to about 3.0 percentage points.

1,000 runs and the wobble that remains

R code runs.R

# Run the stand-in model 1,000 times and count how often C comes out ahead.
# Writes: nothing to disk; the results print in the console.
set.seed(20261028)
draws <- 1000
gain_draws <- 160000 * rnorm(draws, 1.00, 0.20) - rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")

Output of the run runs.log

Runs ahead of 1,000: 592 | average, spread THB: 9236 39057
Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs: 0.400 0.440 0.480 0.510 0.592 0.584 0.592
Share (%): 60 | 0.60 x 0.40: 0.24
One standard error: 0.0155 | x 1.96, in points: 3.0 | range: 0.57 0.63

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

Repeat the count at every price

The fund's 160,000 THB per QALY is a policy choice, and committees argue about it. So repeat the count at other prices. At each price the average net gain and its spread change, and so does the share of runs with C ahead.

At 100,000 THB per QALY, C's average net gain is 100,000 × 1.00 − 150,000 = -50,000 THB with a spread of 31,241 THB, found the same way as in the worked example. The average sits at -1.60 spreads from zero, so C is ahead in only 0.055 of runs. At 150,000 THB the average is 0 and C is ahead in 0.500 of runs. At 200,000 THB and 300,000 THB the share climbs to 0.858 and 0.990.

Plotting that share against every price draws a curve. This curve is called a cost-effectiveness acceptability curve (CEAC): it shows, for each price the fund might pay per QALY, the share of runs in which C is the better buy. Read at 160,000 THB, it gives 0.599, about 0.60, and A's curve, its mirror image, gives 0.401, so the two add up to 1.

The two curves cross at 0.50 at 150,000 THB. That is C's own cost per extra QALY (150,000 THB / 1.00 = 150,000 THB per QALY), so at that price C's average net gain is 0 THB. Below that price A is ahead in most runs, and above it C is. Even at 300,000 THB, C is ahead in 0.990 of runs, not in all of them, because in a few runs its extra health is small or its extra cost is high.

⤢ click to enlarge

Read the blue line from left to right: C's share rises past 0.50 at 150,000 THB to 0.60 at 160,000 THB, while A's dashed line falls as its mirror. The pale band follows the option with the higher average net gain, so it jumps from A to C at 150,000 THB.

R code curve.R

# Repeat the count at five prices per year of full health.
# Writes: nothing to disk; the results print in the console.
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - pnorm(centre / width)), "\n")

Output of the run curve.log

Prices THB: 100,000 150,000 160,000 200,000 300,000
Average net gain THB: -50,000 0 10,000 50,000 150,000
Spread THB: 31,241 38,419 40,000 46,648 64,622
Average / spread: -1.60 0.00 0.25 1.07 2.32
Share of runs where C is ahead: 0.055 0.500 0.599 0.858 0.990
Share where A is ahead: 0.945 0.500 0.401 0.142 0.010

Requires R ≥ 4.5.2 · ggplot2 4.0.3 · works in R 4.6.1 · Download probabilistic-sensitivity-analysis-distributions-and-ceac.R

Full tidy script

# Uncertain inputs as curves, many runs of a stand-in model, and the share of wins at each price
# Give each uncertain input of the fund's decision a curve, rerun a stand-in model 1,000 times.
# ให้ข้อมูลนำเข้าที่ไม่แน่นอนแต่ละตัวมีเส้นโค้งของตัวเอง แล้วรันแบบจำลองแทน 1,000 ครั้ง
# Requires: R >= 4.5.2
# Tested: works in R 4.6.1
# ggplot2 4.0.3 (stable)
# install.packages(c("ggplot2"))
# update.packages(oldPkgs = c("ggplot2"))
# remotes::install_version("ggplot2", "4.0.3")
library(ggplot2)  # Display the shares across prices without saving files.
stopifnot(getRversion() >= "4.5.2")

#%%@note 1. Define the illustrative inputs
# Every input is synthetic and written out once, in the order the article uses it.
# ค่าตั้งต้นทั้งหมดเป็นข้อมูลสมมติ เขียนไว้ครั้งเดียวตามลำดับที่บทความใช้
extra_cost <- 150000
extra_health <- 1.00
price <- 160000
events <- c(8, 13)
patients <- c(66, 60)
big <- function(value) formatC(value, format = "f", digits = 0, big.mark = ",")
ratio <- c(0.25, 0.56, 1.26)
z95 <- 1.96
rate_mean <- 0.30
rate_se <- 0.05
cost_mean <- 10000
cost_se <- 2000
counts <- c(80, 15, 5)
cohort <- 1000
health_avg_sd <- c(1.00, 0.20)
cost_avg_sd <- c(150000, 24000)
wtp <- 160000
draws <- 1000
prices <- c(100000, 150000, 160000, 200000, 300000)
health_c <- c(1.00, 0.20)
cost_c <- c(150000, 24000)

#%%@note 2. Move the trial's ratio to the log scale
# On the log scale the interval is symmetric, so its half width over 1.96 gives the spread.
# บนสเกลลอการิทึม ช่วงความเชื่อมั่นสมมาตร ครึ่งความกว้างหารด้วย 1.96 จึงได้ค่าการกระจาย
logs <- log(ratio)
half_width <- (logs[3] - logs[1]) / 2
spread <- half_width / z95
ends <- logs[2] + c(-1, 1) * z95 * spread

#%%@note 3. Turn a mean and a standard error into curve shapes
# A rate needs a curve between 0 and 1, a cost needs one above zero.
# อัตราต้องใช้เส้นโค้งระหว่าง 0 ถึง 1 ส่วนต้นทุนต้องใช้เส้นโค้งที่อยู่เหนือศูนย์
se_squared <- rate_se^2
patients_worth <- rate_mean * (1 - rate_mean) / se_squared - 1
shapes <- c(rate_mean, 1 - rate_mean) * patients_worth
cost_shape <- (cost_mean / cost_se)^2
cost_scale <- cost_se^2 / cost_mean
cost_ends <- stats::qgamma(c(0.025, 0.975), cost_shape, scale = cost_scale)

#%%@note 4. Draw whole rows of yearly chances
# Rows drawn together always add up to 1, rows drawn one share at a time do not.
# แถวที่สุ่มพร้อมกันรวมได้ 1 เสมอ ส่วนแถวที่สุ่มทีละค่าไม่รวมได้ 1
set.seed(20261028)
together <- matrix(stats::rgamma(9, shape = rep(counts, 3)), nrow = 3)
together <- together / rep(colSums(together), each = 3)
apart <- matrix(stats::rbeta(9, rep(counts, 3), rep(sum(counts) - counts, 3)), nrow = 3)

#%%@note 5. Combine the two uncertain results into one net gain
# Net gain is price x extra health - extra cost, and independent spreads add as squares.
# ผลได้สุทธิคือราคา x สุขภาพที่เพิ่ม - ต้นทุนที่เพิ่ม ค่าการกระจายที่อิสระต่อกันรวมกันแบบกำลังสอง
gain_mean <- wtp * health_avg_sd[1] - cost_avg_sd[1]
parts <- c(wtp * health_avg_sd[2], cost_avg_sd[2]) / 1000
gain_sd <- sqrt(sum(parts^2)) * 1000
above <- stats::pnorm(gain_mean / gain_sd)

#%%@note 6. Run the stand-in model 1,000 times
# Each run draws extra health and extra cost afresh, then records whether C comes out ahead.
# แต่ละครั้งสุ่มสุขภาพที่เพิ่มและต้นทุนที่เพิ่มใหม่ แล้วบันทึกว่า C ชนะหรือไม่
set.seed(20261028)
gain_draws <- 160000 * stats::rnorm(draws, 1.00, 0.20) - stats::rnorm(draws, 150000, 24000)
running <- cumsum(gain_draws > 0) / seq_len(draws)
mc_se <- sqrt(0.60 * 0.40 / draws)

#%%@note 7. Repeat the count at five prices
# The share of runs where C is ahead, price by price.
# สัดส่วนของการรันที่ C ชนะ ทีละระดับราคา
centre <- prices * health_c[1] - cost_c[1]
width <- sqrt((prices * health_c[2])^2 + cost_c[2]^2)
curve_table <- data.frame(price = prices, share_c = stats::pnorm(centre / width))
print(ggplot2::ggplot(curve_table, ggplot2::aes(x = price, y = share_c)) +
  ggplot2::geom_line() + ggplot2::geom_point())

#%%@note 8. Numbers the article quotes
# Print exactly the short boxes in their reading order.
# พิมพ์ผลลัพธ์ตรงกับกล่องตัวอย่างตามลำดับในบทความ
cat("Extra cost THB:", big(extra_cost), "| extra health:", sprintf("%.2f", extra_health),
  "| price THB:", big(price), "| net gain THB:", big(price * extra_health - extra_cost), "\n")
cat("Trial risks:", paste0(sprintf("%.2f", 100 * events / patients), "%"), "| risk ratio:",
  round((events[1] / patients[1]) / (events[2] / patients[2]), 2), "\n")
cat("Weight 1 = full health, 0 = death | 10 years at 0.70 | gain from full health:",
  sprintf("%.1f", 10 * c(0.70, 1 - 0.70)), "\n")
cat("Below | above 0.56:", sprintf("%.2f", diff(ratio)), "| logs:", sprintf("%.3f", logs), "\n")
cat("Half width:", sprintf("%.3f", half_width), "| / 1.96, the spread:", sprintf("%.3f", spread),
  "| rounded centre and spread:", sprintf("%.2f", c(logs[2], spread)), "\n")
cat("Centre -/+ 1.96 spreads:", sprintf("%.3f", ends), "| rounded:", sprintf("%.2f", ends), "\n")
cat("Back as ratios:", sprintf("%.3f", exp(ends)), "| rounded:", sprintf("%.2f", exp(ends)), "\n")
cat("Draws above 1 (%):", sprintf("%.1f", 100 * stats::pnorm(logs[2] / spread)), "\n")
cat("0.30 x 0.70:", rate_mean * (1 - rate_mean), "| 0.05 x 0.05:",
  formatC(se_squared, format = "f", digits = 4), "\n")
cat("Divided:", rate_mean * (1 - rate_mean) / se_squared, "| minus 1:", patients_worth, "\n")
cat("Shapes, responders | non-responders:", shapes, "| their mean:", shapes[1] / sum(shapes), "\n")
cat("Middle 95%:", round(stats::qbeta(c(0.025, 0.975), shapes[1], shapes[2]), 2), "\n")
cat("Mean / standard error:", cost_mean / cost_se, "| squared, the shape:", cost_shape, "\n")
cat("2,000 x 2,000:", big(cost_se^2), "| / 10,000, the scale THB:", cost_scale,
  "| shape x scale:", big(cost_shape * cost_scale), "\n")
cat("Middle 95% THB:", big(cost_ends), "| peak THB:", big((cost_shape - 1) * cost_scale), "\n")
cat("Below | above mean THB:", big(cost_mean - cost_ends[1]), big(cost_ends[2] - cost_mean), "\n")
cat("Average row:", sprintf("%.2f", counts / sum(counts)), "| 1,000 healthy, one year on:",
  cohort * counts / sum(counts), "\n")
cat("Rows drawn together:", sprintf("%.3f", together), "| sums:", colSums(together), "\n")
cat("Rows drawn one share at a time:", sprintf("%.3f", apart), "\n")
cat("Their sums:", sprintf("%.3f", colSums(apart)), "\n")
cat("Average net gain THB:", gain_mean, "| spread parts THB:", 1000 * parts, "\n")
cat("In thousands:", parts, "| squared:", parts^2, "| sum:", sum(parts^2), "| root:",
  sqrt(sum(parts^2)), "| spread THB:", gain_sd, "\n")
cat("Average / spread:", gain_mean / gain_sd, "| share above zero:", round(above, 3), "| rounded:",
  round(above, 2), "| runs of 1,000:", round(1000 * above), "\n")
cat("Runs ahead of 1,000:", sum(gain_draws > 0), "| average, spread THB:", round(mean(gain_draws)),
  round(stats::sd(gain_draws)), "\n")
cat("Share ahead after 10, 25, 50, 100, 250, 500, 1000 runs:",
  sprintf("%.3f", running[c(10, 25, 50, 100, 250, 500, 1000)]), "\n")
cat("Share (%):", 60, "| 0.60 x 0.40:", 0.60 * 0.40, "\n")
cat("One standard error:", formatC(mc_se, format = "f", digits = 4), "| x 1.96, in points:",
  sprintf("%.1f", 196 * mc_se), "| range:", sprintf("%.2f", 0.60 + c(-1, 1) * 1.96 * mc_se), "\n")
cat("Prices THB:", formatC(prices, format = "d", big.mark = ","), "\n")
cat("Average net gain THB:", formatC(centre, format = "d", big.mark = ","), "\n")
cat("Spread THB:", formatC(round(width), format = "d", big.mark = ","), "\n")
cat("Average / spread:", sprintf("%.2f", centre / width), "\n")
cat("Share of runs where C is ahead:", sprintf("%.3f", stats::pnorm(centre / width)), "\n")
cat("Share where A is ahead:", sprintf("%.3f", 1 - stats::pnorm(centre / width)), "\n")

Run this in R to get the same numbers as the article.

The curve tells the fund how sure it can be; it does not pick the option. Choose the option with the higher average net gain at the fund's price: A below 150,000 THB and C above it. Drawing only that option's share at each price gives the cost-effectiveness acceptability frontier (CEAF), the pale band in the figure.

With three or more options, the option that is ahead most often can differ from the one with the highest average net gain, because it may win by a little and lose by a lot. Choose by the average, and report the share beside it (Briggs, Claxton and Sculpher).

Choose a view. In the first, raise a standard error and watch the rate's and the cost's curves widen; in the second, compare rows drawn together with rows drawn one share at a time. In the third, move the price and read the shares for A and C, which is 0.599 for C at 160,000 THB.

What the fund can say

The committee can now give a fuller answer. At 160,000 THB per QALY, C is the better choice on average, with a net gain of 10,000 THB per patient, and it is ahead in about 0.60 of runs, give or take 3.0 percentage points. That leaves about 0.40 of runs in which the health C adds is worth less than the extra 150,000 THB it costs.

The answer is only as good as the curves. Each curve must come from the evidence behind its input (Briggs and colleagues). Inputs that rise and fall together, such as two costs from the same hospitals, must be drawn together, as the row of chances was, not one at a time. Here we treated extra health and extra cost as independent.

A share of 0.60 also raises a new question: would it be worth paying for more research before deciding? The next article in this series puts a value on removing this uncertainty (Fenwick and colleagues).

Terms as this article uses them

Distribution (การแจกแจง): A curve of the values an uncertain input could take, high where values are plausible.
Normal distribution (การแจกแจงแบบปกติ): The symmetric bell; about 95% of its values lie within 1.96 standard deviations of the mean.
Lognormal distribution (การแจกแจงแบบล็อกนอร์มอล): The curve of a ratio whose log is a Normal bell; it never goes below 0 and has a longer right side.
Beta distribution (การแจกแจงแบบเบตา): A hump between 0 and 1 for an uncertain probability, set by two counts such as 24.9 and 58.1.
Gamma distribution (การแจกแจงแบบแกมมา): A curve above 0 with a long right tail, for an uncertain cost, set by a shape and a scale such as 25 and 400 THB.
Dirichlet distribution (การแจกแจงแบบดีริชเลต์): A way to draw a whole row of chances at once, so every drawn row adds up to 1.
Probabilistic sensitivity analysis (PSA) (การวิเคราะห์ความไวเชิงความน่าจะเป็น): Rerunning the model many times with every uncertain input drawn from its curve, and reporting the share of runs.
Monte Carlo simulation (การจำลองแบบมอนติคาร์โล): Answering a question by drawing random inputs many times and counting; the wobble left is the Monte Carlo error.
Cost-effectiveness acceptability curve (CEAC) (เส้นโค้งการยอมรับความคุ้มค่า): The share of runs in which an option is the better buy, plotted against the price paid per QALY.
Acceptability frontier (CEAF) (ขอบเขตการยอมรับความคุ้มค่า): The share of runs for the option with the higher average net gain, price by price.

Connect to the foundation articles

The foundation article on how decision trees and cohort models turn evidence into a decision shows where this analysis sits in a full model. The one on chance, bias and confounding separates the play of chance, which these curves describe, from bias, which no number of runs removes.

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

Briggs Andrew H., Weinstein Milton C., Fenwick Elisabeth A.L., Karnon Jonathan, Sculpher Mark J., Paltiel A. David. Model Parameter Estimation and Uncertainty: A Report of the ISPOR-SMDM Modeling Good Research Practices Task Force-6. Value in Health. 2012;15(6):835-842. doi:10.1016/j.jval.2012.04.014. https://api.crossref.org/works/10.1016/j.jval.2012.04.014
Andrew Briggs, Karl Claxton, Mark Sculpher. Decision Modelling for Health Economic Evaluation (Handbooks for Health Economic Evaluation). Oxford University Press; 2006. ISBN 9780198526629. https://openlibrary.org/isbn/9780198526629.json
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
Fenwick Elisabeth, Steuten Lotte, Knies Saskia, Ghabri Salah, Basu Anirban, Murray James F., Koffijberg Hendrik (Erik), Strong Mark, Sanders Schmidler Gillian D., Rothery Claire. Value of Information Analysis for Research Decisions-An Introduction: Report 1 of the ISPOR Value of Information Analysis Emerging Good Practices Task Force. Value in Health. 2020;23(2):139-150. doi:10.1016/j.jval.2020.01.001. https://api.crossref.org/works/10.1016/j.jval.2020.01.001

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 (this article)
Model Validation and Verification: Test the Question, Equations and Code
Expected Value of Perfect Information: What Could Better Evidence Be Worth?

Key takeaways

Give each uncertain input a curve that fits it. A ratio gets a bell on the log scale, a rate a hump between 0 and 1, a cost a curve above 0, and a row of chances a whole-row draw.
Rerun the model with fresh draws from every curve: at 160,000 THB per QALY, C is ahead in about 599 of 1,000 runs, a probability of 0.60.
Report the wobble: with 1,000 runs, a share of 0.60 is good to about 3.0 percentage points either way.
Repeat the count at every price: the curve reads 0.50 at 150,000 THB and 0.60 at 160,000 THB. Choose the option with the higher average, which is C above 150,000 THB.

Probabilistic Sensitivity Analysis: From Parameter Distributions to CEAC

อ่านฉบับภาษาไทย (Thai version)

Abstract

Is option C worth 150,000 THB more?

The case in numbers

One estimate, many plausible values

The risk ratio: a symmetric bell on the log scale

From the trial's interval to the curve

The response rate: a curve that stays between 0 and 1

The cost of a course: a curve that stays above zero

From a rate and its standard error to two shapes

From a cost and its standard error to shape and scale

A row of chances that must add up to 1

Rows drawn together and one share at a time

Run the model many times

Worked example: how often is C ahead?

The name for it

The net gain's average, spread and share above zero

How many runs are enough?

1,000 runs and the wobble that remains

Repeat the count at every price

The share of runs at five prices

Choose by the average, report the share

What the fund can say

Terms as this article uses them

Connect to the foundation articles

Try it in the app

References

The series, in learning order

Key takeaways

Comments