Decision Trees, Markov Cohorts and Probabilistic Sensitivity Analysis: How a Model Turns Evidence into a Decision

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

Abstract

A district hospital committee has been told that a new preventive treatment, C, gives 6.00 quality-adjusted life years (QALYs) per patient over a lifetime for 250,000 THB. Usual care gives 5.00 for 100,000 THB. The trial behind C lasted one year. This article opens the model that turned that year into a lifetime. A decision tree prices the first year: each outcome's chance times its cost or health, added up. One year alone makes C cost close to a million baht per QALY, because its benefit comes later. So the model follows 1,000 people who move once a year between Well, Sick and Dead until all have died, which gives the committee's numbers. Changing one input at a time shows that the trial's risk ratio decides the answer. Rerunning the model 1,000 times with its two least certain inputs drawn at random puts C ahead in 542 of the 1,000 runs, 54.2%: a share, not one number.

🔍 Click to enlarge Click anywhere to close

Where do lifetime numbers come from?

A district hospital committee must decide whether to fund C, a new preventive treatment for adults with a chronic condition, or keep usual care, A. The analysts' report says that over a lifetime A costs 100,000 THB (Thai baht) per patient and gives 5.00 QALYs. C costs 250,000 THB and gives 6.00. 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 QALY, DALY and the Utilities Behind Them: What Is Actually Being Measured.

So C buys 1.00 extra QALY for 150,000 THB more, a price of 150,000 THB per extra QALY, under the committee's limit of 160,000 THB. The committee's limit, 160,000 THB per QALY, is also called the threshold.

Value the extra QALY at the threshold and subtract the extra cost: 160,000 × 1.00 − 150,000 = 10,000 THB. That is the incremental net monetary benefit (INMB). Positive means the extra cost is worth paying. We call this 10,000 THB C's lead; the full lesson is in the article on the ICER and net benefit.

The trial behind C followed patients for one year. On usual care 13 of 60 had the complication the treatment prevents, 21.67%, against 8 of 66 on C, 12.12%.

So a committee member asks two questions. Where do lifetime numbers come from, when nobody was watched for a lifetime? And the trial's risk ratio of 0.56 could, by its 95% confidence interval, lie anywhere from 0.25 to 1.26: how sure is C's lead of 10,000 THB?

Both questions have the same answer, a model. A common shortcut treats a model's output as one fixed number. This article shows where the number comes from and how far it can move. It uses a decision tree for the first year, one group of patients followed year by year for the lifetime, and the model rerun with its least certain inputs drawn at random.

The first year as a decision tree

Read first on this site. These posts teach in full how to build a decision tree step by step and when a tree is the right structure. This section keeps only the committee's first year.

• Building a Decision Tree Step by Step: Follow Outcomes, Costs and Health
• Choosing a Decision-Analytic Model: Let the Clinical Pathway Decide

Start with what the trial saw: one year. Within it, three things can happen to a patient: they stay well, they have the complication, or they die of another cause. A life table for people of this age shows that 6.5% of them die of other causes each year, whatever the treatment. On usual care 21.67% have the complication, so 71.83% stay well; on C 12.12% have it, so 81.38% stay well.

Each outcome has a cost and a health value for the year. A study of patients' preferences gives each kind of year a quality weight, its utility: 0.80 QALYs for a year lived well and 0.50 for a year with the complication. On C a well year also carries the drug, at its price of 37,800 THB a year, and the hospital's cost list prices a year of complication care at 16,900 THB. A death adds no cost and no health, and costs that both options share are left out, so a well year on usual care costs 0.

To get one patient's expected cost, multiply each outcome's cost by its chance and add them up. On usual care that is 21.67% × 16,900 = 3,662 THB, since the other two branches cost nothing. On C it is 81.38% × 37,800 + 12.12% × 16,900 = 32,810 THB. Health works the same way: 71.83% × 0.80 + 21.67% × 0.50 = 0.683 QALYs on usual care, and 0.712 on C.

This drawing is a decision tree [6]: a square for the choice, a circle where chance takes over, and a payoff at the end of each branch. Working from the ends back to the circles is called rolling back the tree, and its rule fits on one line, where p is a branch's chance and V its payoff, a cost or QALYs:

$$E[V]=\sum_k p_k\,V_k$$

Why one year is not enough

Read the tree as the whole answer and C looks hopeless. In year 1 it costs 29,148 THB more per patient and adds 0.029 QALYs. That is close to a million baht for each extra QALY, more than six times the committee's limit.

The tree is not wrong; it is too short. C pays off through the complications it prevents, and a complication does not end with the year: the patient stays sick, needs care every year after, and dies sooner. A tree would need a new fork for every year of every patient's life, so the drawing would grow with each year it covers. A neater way is to stop following one patient down every branch and follow a whole group instead.

# The first year as a decision tree: each branch's chance times its cost or health, added up.
f <- function(x, d = 0) formatC(x, format = "f", digits = d, big.mark = ",")
p_sick <- c(A = 13 / 60, C = 8 / 66)            # complication within the year, from the trial
p_well <- 1 - p_sick - 0.065                     # 6.5% die of other causes; the rest stay well
cost <- p_well * c(0, 37800) + p_sick * 16900    # C's drug while well; a year of complication care
qaly <- p_well * 0.80 + p_sick * 0.50            # a death adds no cost and no health
cat("Stay well (%): A", f(100 * p_well[1], 2), " C", f(100 * p_well[2], 2), "\n")
cat("Complication (%): A", f(100 * p_sick[1], 2), " C", f(100 * p_sick[2], 2), "\n")
cat("Risk ratio, C over A:", f(p_sick[2] / p_sick[1], 2), "\n")
cat("Expected cost per patient (THB): A", f(cost[1]), " C", f(cost[2]), "\n")
cat("Expected QALYs per patient: A", f(qaly[1], 3), " C", f(qaly[2], 3), "\n")
cat("C minus A: cost (THB)", f(diff(cost)), " QALYs", f(diff(qaly), 3), "\n")

Stay well (%): A 71.83  C 81.38
Complication (%): A 21.67  C 12.12
Risk ratio, C over A: 0.56
Expected cost per patient (THB): A 3,662  C 32,810
Expected QALYs per patient: A 0.683  C 0.712
C minus A: cost (THB) 29,148  QALYs 0.029

One group, three states, one step a year

Read first on this site. This post teaches in full how a share for one step comes from a rate, and why the length of the step stays in the calculation. The full lesson on moving a whole group between states is linked under the model's name below; this section keeps only the committee's group of 1,000.

• Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation

Put 1,000 people in one of three states, Well, Sick or Dead, and once a year move a fixed share of each state to another; the shares are the tree's chances. From Well, 21.67% fall sick on usual care, 6.5% die of another cause and 71.83% stay. From Sick, 13% die each year, the yearly share of deaths in a cohort study of patients followed after the complication, and 87% stay sick. Nobody recovers, and nobody leaves Dead.

These shares are chances per one-year step, not rates. A rate counts events per person-time. A share for one step is derived from it with the length of the step kept in the calculation. The life table's 6.5% a year, for example, matches a death rate of 0.0672 per person-year, because $1-e^{-0.0672}=0.065$.

The same one-year step decides how the risk ratio is used. The trial's shares can be used as they are only because the trial and the step both last one year. For the same reason its risk ratio of 0.56, a ratio of one-year risks, multiplies usual care's one-year share; a hazard ratio would be applied to the rate instead.

After year 1 on usual care the group holds 718.3 Well, 216.7 Sick and 65.0 Dead: the tree's three branches, now as counts. Fractions of a person are fine, because these are expected counts, the average over many such groups. Year 2 starts from those counts, and the steps below show how 216.7 Sick become 344.1.

# 1,000 people start Well and move once a year by the shares in each row; counted at year end.
f <- function(x, d = 1) formatC(x, format = "f", digits = d, big.mark = ",")
p_sick <- 13 / 60                         # usual care: the trial's one-year complication risk
move <- rbind(c(1 - p_sick - 0.065, p_sick, 0.065),  # from Well: stay, fall sick, die
              c(0, 0.87, 0.13),                    # from Sick: stay sick, die
              c(0, 0, 1))                          # Dead stays dead
utility <- c(0.80, 0.50, 0)
year1 <- c(1000, 0, 0) %*% move; year2 <- year1 %*% move
cat("Year 1, Well Sick Dead:", f(year1), "\n")
cat("Year 2, Well Sick Dead:", f(year2), "\n")
cat("Year 2 Sick = still sick + newly sick:", f(year1[2] * 0.87), "+", f(year1[1] * p_sick), "\n")
cat("QALYs per person, years 1, 2:", f(c(year1 %*% utility, year2 %*% utility) / 1000, 3), "\n")
cat("Death rate behind the 6.5% yearly share (per person-year):", f(-log(1 - 0.065), 4), "\n")

Year 1, Well Sick Dead: 718.3 216.7 65.0
Year 2, Well Sick Dead: 516.0 344.1 139.9
Year 2 Sick = still sick + newly sick: 188.5 + 155.6
QALYs per person, years 1, 2: 0.683 0.585
Death rate behind the 6.5% yearly share (per person-year): 0.0672

The name in the title: a Markov cohort

This kind of model carries the name in this article's title: a Markov model follows a group of people between a few health states, one cycle at a time. In the series' own example, 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.

Its defining rule is the one used above: where a person goes next depends only on the state they are in now, not on how they got there [1, 2]. The group it follows is called a cohort and each step a cycle, so the committee's model is a Markov cohort with a one-year cycle.

An overview on this site sets decision trees and these models side by side:

• Decision Trees and Markov Models in Health Economic Evaluation: Principles, Construction, and Applications

# The series' own example of one cycle: 1,000 people start healthy and move by one row of shares.
start <- 1000
row <- c(healthy = 0.80, sick = 0.15, dead = 0.05)   # from Healthy: stay, fall sick, die
cat("People who start healthy:", start, "\n")
cat("Still healthy after the first cycle:", start * row[["healthy"]], "\n")
cat("Fallen sick in the first cycle:", start * row[["sick"]], "\n")
cat("Died in the first cycle:", start * row[["dead"]], "\n")

People who start healthy: 1000
Still healthy after the first cycle: 800
Fallen sick in the first cycle: 150
Died in the first cycle: 50

A lifetime, added up

Keep stepping for 100 years, until everyone has died, and add up the years each person spends in each state. On usual care that is 2.55 years Well and 5.92 years Sick per person; on C it is 4.37 years Well and 5.01 Sick.

A shortcut gives the same years by hand. Someone in Well stays there 71.83 ÷ (21.67 + 6.5) = 2.55 years on average: the yearly chance of staying divided by the yearly chance of leaving. Of those who leave Well, the share 21.67 ÷ (21.67 + 6.5) fall sick; each is counted Sick at the end of that year and then stays 87 ÷ 13 more years, so 21.67 ÷ (21.67 + 6.5) × (1 + 87 ÷ 13) = 5.92. On C the same steps give 81.38 ÷ (12.12 + 6.5) = 4.37 and 5.01.

Now value those years. Usual care gives 2.55 × 0.80 + 5.92 × 0.50 = 5.00 QALYs and costs 5.92 × 16,900 = 100,000 THB to the nearest 1,000. C gives 4.37 × 0.80 + 5.01 × 0.50 = 6.00 QALYs. It costs 4.37 × 37,800 + 5.01 × 16,900 = 250,000 THB to the nearest 1,000.

These are the committee's numbers. C's extra QALY comes from more years spent well, and C costs more because each of those years carries the drug. The lead is 160,000 × 1.00 − 150,000 = 10,000 THB.

# Run each option's group of 1,000 for 100 years, until all have died, and add up the years.
f <- function(x, d = 0) formatC(x, format = "f", digits = d, big.mark = ",")
years_in <- function(p_sick, m = c(1000, 0, 0), total = 0) {   # m: Well, Sick, Dead at the start
  move <- rbind(c(1 - p_sick - 0.065, p_sick, 0.065), c(0, 0.87, 0.13), c(0, 0, 1))
  for (t in 1:100) { m <- m %*% move; total <- total + m }; total / 1000 }  # years per person
a <- years_in(13 / 60); c_ <- years_in(8 / 66)       # A usual care; C, from the trial's risks
qaly <- c(sum(a * c(0.80, 0.50, 0)), sum(c_ * c(0.80, 0.50, 0)))
cost <- c(sum(a * c(0, 16900, 0)), sum(c_ * c(37800, 16900, 0)))
cat("Well and sick years per person: A", f(a[1:2], 2), " C", f(c_[1:2], 2), "\n")
cat("QALYs: A", f(qaly[1], 2), " C", f(qaly[2], 2), "\n")
cat("Cost, nearest 1,000 THB: A", f(round(cost[1], -3)), " C", f(round(cost[2], -3)), "\n")
cat("C minus A: cost", f(round(diff(cost), -3)), " QALYs", f(diff(qaly), 2), "\n")
cat("Lead of C at 160,000 per QALY (THB):", f(round(160000 * diff(qaly) - diff(cost), -3)), "\n")

Well and sick years per person: A 2.55 5.92  C 4.37 5.01
QALYs: A 5.00  C 6.00
Cost, nearest 1,000 THB: A 100,000  C 250,000
C minus A: cost 150,000  QALYs 1.00
Lead of C at 160,000 per QALY (THB): 10,000

What the model forgets

Read first on this site. These posts teach in full when a model must remember how long someone has been sick, and which count earns each year. This section keeps only what both rules do to the committee's numbers.

• Choosing a Decision-Analytic Model: Let the Clinical Pathway Decide
• Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time

The rule that only the present state matters is also the model's main limit. Here a patient in their first year with the complication and one in their tenth both die at the same 13% each year. If the first year after a complication were the dangerous one, the model would need to split Sick in two: a first sick year that everyone leaves after one step, which modellers call a tunnel state, and later sick years. Another way is to follow patients one at a time, each with their own history [2].

Counting everyone at the end of each year is a choice too. A person who dies during year 1 earns nothing for it. C's drug, which all 1,000 start on the first day, is charged only to those still well at the year's end.

The usual refinement counts each year at the average of its start and end, as if people moved halfway through it (called the half-cycle correction) [3]. In the series' own example, a cycle's healthy year is weighted 1.0 and a sick year 0.6 (that example's weights, not the committee's 0.80 and 0.50), and 800 of 1,000 people stay healthy while 150 fall sick in the first cycle.

Counted at the end, that is 0.800 × 1.0 + 0.150 × 0.6 = 0.890 QALYs per person. Counted at the average of start and end, it is (1.000 + 0.800) ÷ 2 × 1.0 + 0.150 ÷ 2 × 0.6 = 0.945. The full lesson is in Building a Markov Cohort Model Step by Step: Follow the Cohort Through Time.

Counting each year at the average of its start and end would have added half a year more in Well for everyone. Both options would gain 0.40 QALYs (half of 0.80), so C's extra stays 1.00 QALY. But C would also pay for half a year more of its drug, 0.5 × 37,800 = 18,900 THB.

C's lead of 10,000 THB would then become 10,000 − 18,900 = −8,900, about 9,000 THB behind. Its cost per extra QALY would rise from 150,000 to 168,900 THB, above the 160,000 limit. So if the committee had counted this way, the verdict would flip.

The committee counted at the year's end, so the lead is 10,000 THB, the cost is 150,000 THB per QALY, and every run in this article keeps that count. The counting rule is an input to report with the result; the explorer above lets you switch to the start of each year or the average of start and end.

# The series' own example, first cycle: 1,000 start healthy; 800 stay healthy, 150 fall sick.
f <- function(x, d) formatC(x, format = "f", digits = d)
w <- c(healthy = 1.0, sick = 0.6)                  # that example's weights, not the committee's
start <- c(1000, 0) / 1000; end <- c(800, 150) / 1000   # shares at the cycle's start and end
cat("Weights of a healthy and a sick year:", f(w, 1), "\n")
cat("Shares at the start of the cycle, healthy and sick:", f(start, 3), "\n")
cat("Shares at the end of the cycle, healthy and sick:", f(end, 3), "\n")
cat("QALYs per person, counted at the end:", f(sum(end * w), 3), "\n")
cat("QALYs per person, average of start and end count:", f(sum((start + end) / 2 * w), 3), "\n")

Weights of a healthy and a sick year: 1.0 0.6
Shares at the start of the cycle, healthy and sick: 1.000 0.000
Shares at the end of the cycle, healthy and sick: 0.800 0.150
QALYs per person, counted at the end: 0.890
QALYs per person, average of start and end count: 0.945

# Count each year at its end, or at the average of its start and end; all 1,000 start Well.
f <- function(x, d = 0) formatC(x, format = "f", digits = d, big.mark = ",")
yrs <- function(p) c((1 - p - 0.065) / (p + 0.065), p / (p + 0.065) / 0.13)  # Well, Sick years
a <- yrs(13 / 60); c_ <- yrs(8 / 66)                          # counted at year ends
dq <- round(sum((c_ - a) * c(0.80, 0.50)), 2)                  # C's extra QALYs, 1.00
dc <- round(c_[1] * 37800 + (c_[2] - a[2]) * 16900, -3); extra <- 0.5 * 37800  # THB
lead <- 160000 * dq - dc - c(0, extra)            # C's lead: year ends, average count
cat("QALYs added to each option by the average count:", f(0.5 * 0.80, 2), "\n")
cat("Drug cost added to C by the average count (THB):", f(extra), "\n")
cat("Lead of C (THB), year ends and average:", f(lead), " nearest 1,000:", f(round(lead, -3)), "\n")
cat("Cost per extra QALY (THB), year ends and average:", f(dc / dq), f((dc + extra) / dq), "\n")

QALYs added to each option by the average count: 0.40
Drug cost added to C by the average count (THB): 18,900
Lead of C (THB), year ends and average: 10,000 -8,900  nearest 1,000: 10,000 -9,000
Cost per extra QALY (THB), year ends and average: 150,000 168,900

Where every input comes from

Each number the model uses comes from a different place, as the table shows, and the model's job is to join them. In this teaching example every input is illustrative. In a real model each one would carry its source, its value and how uncertain it is, in a table like this one.

The two ways out of Well come from different sources, the trial and the life table, yet they share the same patients. Someone who dies of another cause can no longer have the complication. Turning each into a yearly share on its own and adding the two can overstate how many leave Well each year. A fuller model therefore adds the two rates first and then splits the one exit between them, as Rates, Probabilities and Transition Probabilities: Keep Time in the Calculation shows.

The drug is taken only while a patient is well; after a complication, complication care replaces it. Every year also counts the same however far ahead it lies, to match the committee's figures; the costing article shows how future years are usually shrunk.

Check the model before using it

Read first on this site. This post teaches in full how to test a model's equations and code, and how to judge whether it is the right model. This section keeps only the checks on the committee's model.

• Model Validation and Verification: Test the Question, Equations and Code

Before trusting the model, check it in two ways. First, was it built right? Every row of shares adds up to 100%, nobody comes back from Dead, and the first year gives back what it was built from: 216.7 of 1,000 fall sick on usual care, the trial's 21.67%, and 12.12% on C. Modellers call these checks verification.

Second, is it the right model for these patients? The years lived can be compared with a life table, and the model's forecasts with data it was not built on, such as patients followed for longer than the trial. That is validation, from face validity (does the structure make clinical sense?) to external validity (does it predict data it was not built on?) [2].

Which input could flip the answer?

Read first on this site. This post teaches sensitivity analysis for health economic evaluation in full, one input at a time and all inputs together. This section keeps only the inputs of the committee's model.

• Uncertainty and Sensitivity Analysis in Health Economic Evaluation: Principles, Tools, and Application

A lead of 10,000 THB on a lifetime cost of 250,000 THB is small. One way to see how fragile it is: change one input at a time to a low and a high value. Keep the rest at their usual values, and run the model again [4].

At the ends of the trial's interval, a risk ratio of 0.25 gives C a lead of 188,000 THB, and 1.26 leaves it 131,000 THB behind. A drug price of 30,000 or 45,000 THB a year moves the lead to 44,000 or −21,000 THB. A utility of 0.40 or 0.60 for a sick year moves it to 25,000 or −4,000 THB. A sick-care cost of 11,000 or 24,000 THB a year gives 5,000 or 17,000 THB: a dearer complication helps C, because C prevents it.

Drawn as bars, widest at the top, these results make a tornado diagram. The bars show the lead in baht rather than a cost per QALY, because a ratio swings wildly, and can change sign, when the extra QALYs come close to zero. Only the risk ratio swings the answer far to both sides.

The lead is exactly zero at a risk ratio of 0.59, the switching value. C stays ahead only while the true risk ratio is below it, yet the trial cannot rule out anything from 0.25 to 1.26.

# Move one input at a time to a low and a high value; C's lead (THB) at 160,000 THB per QALY.
f <- function(x) formatC(round(x, -3), format = "f", digits = 0, big.mark = ",")
lead <- function(rr = 8 / 66 / (13 / 60), drug = 37800, sick = 16900, u_sick = 0.50) {
  yrs <- function(p) c((1 - p - 0.065) / (p + 0.065), p / (p + 0.065) / 0.13)  # well, sick years
  a <- yrs(13 / 60); c_ <- yrs(13 / 60 * rr)
  160000 * sum((c_ - a) * c(0.80, u_sick)) - c_[1] * drug - (c_[2] - a[2]) * sick }
cat("Usual values:", f(lead()), "\n")
cat("Risk ratio 0.25 and 1.26:", f(lead(rr = 0.25)), "and", f(lead(rr = 1.26)), "\n")
cat("Drug 30,000 and 45,000:", f(lead(drug = 30000)), "and", f(lead(drug = 45000)), "\n")
cat("Sick care 11,000 and 24,000:", f(lead(sick = 11000)), "and", f(lead(sick = 24000)), "\n")
cat("Sick utility 0.40 and 0.60:", f(lead(u_sick = 0.40)), "and", f(lead(u_sick = 0.60)), "\n")
cat("Risk ratio at a zero lead:", round(uniroot(function(r) lead(rr = r), c(0.3, 1))$root, 2), "\n")

Usual values: 10,000
Risk ratio 0.25 and 1.26: 188,000 and -131,000
Drug 30,000 and 45,000: 44,000 and -21,000
Sick care 11,000 and 24,000: 5,000 and 17,000
Sick utility 0.40 and 0.60: 25,000 and -4,000
Risk ratio at a zero lead: 0.59

Run the model 1,000 times

Changing one input at a time misses what happens when several are off at once, and it treats 1.26 as if it were as likely as 0.56. The trial says values near 0.56 are likely and values near 1.26 are possible but rare. So run the model many times, and each time draw the uncertain inputs at random, each value turning up as often as the evidence makes it plausible. A full analysis draws every uncertain input this way [4]; here we draw the two we are least sure of, the risk ratio and the yearly cost of sick care, and hold the rest at their usual values.

For the risk ratio, the draws are centred on 0.56 and spread so that 95 in 100 land inside the trial's own interval, 0.25 to 1.26. About 1 in 40 fall below it, 1 in 40 above, and none below zero, as a risk ratio cannot be negative. For the yearly cost of sick care, the draws centre on 16,900 THB, never fall below zero and stretch further up than down, because a few patients cost far more than the rest. About 19 in 20 land between 11,000 and 24,000 THB.

The first run drew a risk ratio of 0.69 and a sick-care cost of 23,016 THB, and C came out 28,793 THB behind. The second drew 0.36 and 19,712 THB, and C led by 111,378 THB. The third drew 0.59, close to the risk ratio at which C's lead vanishes, and 21,313 THB: C led by only 2,244 THB, helped by the dearer sick care it prevents. In every run the Well row still adds up to 100%, because the share who stay well is whatever is left after the drawn complication share and the 6.5% who die of other causes.

Over 1,000 runs, C comes out ahead in 542 of them, 54.2%. That share answers the committee's second question: on the evidence so far, C is the better buy in a little more than half of the plausible worlds. The runs keep the committee's year-end count.

The method you have just seen has a name and a full lesson on this site. That post teaches in full how each uncertain input gets a curve of plausible values, how many runs are enough and how the runs are summed up. This section kept only the two inputs of the committee's model that we are least sure of.

• Probabilistic Sensitivity Analysis: From Parameter Distributions to CEAC

A share and an average, not one number

The right-hand chart repeats the count of runs with C ahead at every price per QALY the committee might accept, so it rises as the limit rises. Even at the highest prices on its right, it stops short of 100%: in some runs C adds little health, or none, while still costing far more, so no price makes it worthwhile [5].

The decision itself rests on the average lead across the runs, 14,054 THB, which is higher than the 10,000 THB from the usual inputs. The lead bends with the risk ratio: low draws help C more than high draws hurt it. So fund the option with the highest average lead, which here is C at 14,054 THB, and report the share beside it. The average weighs how much C wins or loses by, not only how often.

Never average the 1,000 costs per QALY instead. In runs where C adds almost no health that ratio shoots to huge positive or negative values, while the lead stays an ordinary sum of money.

1,000 runs are a sample too: in 19 of 20 repeats with fresh draws, the share would land within about 3.1 points of 54.2%. Fixing the starting number of the random draws, called the seed, makes the same 1,000 runs come out every time, so anyone can check them.

# Rerun the model 1,000 times; each run draws the risk ratio and the yearly cost of sick care.
f <- function(x, d = 0) formatC(x, format = "f", digits = d, big.mark = ",")
set.seed(2026); n <- 1000
rr <- exp(rnorm(n, log(8 / 66 / (13 / 60)), sqrt(1 / 8 - 1 / 66 + 1 / 13 - 1 / 60)))  # the trial
care <- rgamma(n, shape = 25, scale = 16900 / 25)       # centred on 16,900 THB; never below zero
well_years <- function(p) (1 - p - 0.065) / (p + 0.065)  # years in Well, counted at year ends
sick_years <- function(p) p / (p + 0.065) / 0.13         # share who fall sick x 1 / 0.13 years
p <- 13 / 60 * rr; a <- 13 / 60                           # C's yearly risk in each run; A's
d_well <- well_years(p) - well_years(a); d_sick <- sick_years(p) - sick_years(a)
lead <- 160000 * (0.80 * d_well + 0.50 * d_sick) - 37800 * well_years(p) - care * d_sick
cat("Runs 1-3: risk ratio", f(rr[1:3], 2), " sick care", f(care[1:3]), " lead", f(lead[1:3]), "\n")
w <- mean(lead > 0); cat("C ahead in", sum(lead > 0), "of", n, "runs:", f(100 * w, 1), "%\n")
cat("Average lead:", f(mean(lead)), " give or take:", f(196 * sqrt(w * (1 - w) / n), 1), "points\n")

Runs 1-3: risk ratio 0.69 0.36 0.59  sick care 23,016 19,712 21,313  lead -28,793 111,378 2,244
C ahead in 542 of 1000 runs: 54.2 %
Average lead: 14,054  give or take: 3.1 points

The name in the title: probabilistic sensitivity analysis

This is the method in the title: probabilistic sensitivity analysis (PSA) reruns the model many times with every uncertain input drawn from its distribution. The result is a share of runs rather than one number, and it shows how sensitive the answer is to inputs drawn by chance [4]. The series' own example is a different one from the committee's. There, an extra net benefit that averages 10,000 THB with a standard deviation of 40,000 THB gives about 599 of 1,000 runs in favour of the treatment.

A cost-effectiveness acceptability curve (CEAC) shows, at each threshold, the share of simulation runs in which a treatment is the best buy. In the series' example it reads 0.60 at 160,000 THB per QALY, because about 599 of 1,000 runs favour the treatment. The committee's curve, on the right of the figure above, reads 54.2% there [5]. The same lesson teaches the method and the curve in full: Probabilistic Sensitivity Analysis: From Parameter Distributions to CEAC.

# The series' own example of 1,000 reruns: the extra net benefit changes from run to run.
f <- function(x, d = 0) formatC(x, format = "f", digits = d, big.mark = ",")
gain_mean <- 10000; gain_sd <- 40000; runs <- 1000    # THB, THB, reruns of the model
share <- pnorm(gain_mean / gain_sd)                    # share of runs with an extra gain above 0
cat("Extra net benefit, average (THB):", f(gain_mean), "\n")
cat("Extra net benefit, standard deviation (THB):", f(gain_sd), "\n")
cat("Runs that favour the treatment, of", f(runs), "runs:", round(share * runs), "\n")
cat("Height of the acceptability curve at 160,000 THB per QALY:", f(share, 2), "\n")

Extra net benefit, average (THB): 10,000
Extra net benefit, standard deviation (THB): 40,000
Runs that favour the treatment, of 1,000 runs: 599
Height of the acceptability curve at 160,000 THB per QALY: 0.60

What 1,000 runs cannot fix

Every run uses the same model and the same trial. Suppose the trial was biased, for instance because sicker patients ended up on usual care. Then every draw of the risk ratio is centred on the wrong value, and 1,000 runs give a precise picture of the wrong answer Chance, Bias and Confounding: Why an Association Is Not Yet a Treatment Effect.

If the structure is wrong, for instance if the first sick year is the dangerous one, no amount of drawing fixes it: build the other structure and compare the answers. So report the inputs with their sources, the structure, the one-at-a-time results and the share from the runs together, so that a reader can check each part.

The series, in learning order

1. Health Economics Series Guide: Read the Evidence Before the Costs
2. Incidence, Prevalence and the Right Denominator: The Epidemiology Beneath Cost-of-Illness and Markov Models
3. Design, Model and Estimand: What Each Study Design Can Actually Estimate
4. Risk Ratio, Odds Ratio, Risk Difference and NNT: One 2×2 Table, Four Answers
5. The Hazard Ratio Is Not a Risk Ratio: Proportional Hazards, Restricted Mean Survival Time and Model Inputs
6. Chance, Bias and Confounding: Why an Association Is Not Yet a Treatment Effect
7. Beyond the Evidence Pyramid: GRADE Certainty and Choosing Effectiveness Inputs for Economic Models
8. From Sensitivity to Cost per Case Detected: Predictive Values, Prevalence and the Economics of Screening
9. Whose Cost, in Which Year? Perspective, Opportunity Cost, Inflation and Discounting
10. QALY, DALY and the Utilities Behind Them: What Is Actually Being Measured
11. ICER, Dominance and the Cost-Effectiveness Plane: Reading the Thai Threshold of 160,000 Baht per QALY
12. Decision Trees, Markov Cohorts and Probabilistic Sensitivity Analysis: How a Model Turns Evidence into a Decision (this article)
13. Cost-Effective Is Not Affordable, and Efficient Is Not Fair: Budget Impact and Equity
14. Parameter Sourcing, Reference Cost Lists and Utility Catalogues When Local Data Are Sparse
15. Early Health Technology Assessment, Headroom Pricing and Benefit-Package Decisions
16. Communicating Health Economic Evidence: From Model Outputs to Policy Briefs

Related in the wiki: Incidence, Prevalence and the Right Denominator: The Epidemiology Beneath Cost-of-Illness and Markov Models · The Hazard Ratio Is Not a Risk Ratio: Proportional Hazards, Restricted Mean Survival Time and Model Inputs · Beyond the Evidence Pyramid: GRADE Certainty and Choosing Effectiveness Inputs for Economic Models · ICER, Dominance and the Cost-Effectiveness Plane: Reading the Thai Threshold of 160,000 Baht per QALY · Chance, Bias and Confounding: Why an Association Is Not Yet a Treatment Effect · Whose Cost, in Which Year? Perspective, Opportunity Cost, Inflation and Discounting · QALY, DALY and the Utilities Behind Them: What Is Actually Being Measured