1 Overview

Individualized cost-effectiveness analysis (iCEA) evaluates the cost-effectiveness of treatments at the individual (or subgroup) level. This has two major implications:

The hesim package help facilitate iCEA by providing a number of functions for analyzing subgroup level health and cost outcomes from simulation models that quantify parameter uncertainty using probabilistic sensitivity analysis (PSA). These functions take simulation output and generate measures commonly used for technology assessment including:

The rest of this document provides an overview of iCEA and how it can be conducted using hesim. The perspective is Bayesian in nature in that it is concerned with estimating the entire distribution of outcomes rather than just expected values (Baio 2012; Baio and Dawid 2015). It also stresses that both optimal treatments and the cost-effectiveness of those treatments vary across individuals (Basu and Meltzer 2007; Espinoza et al. 2014).

2 Net monetary benefits

Decision analysis provides a formal framework for making treatment decisions based on the utility that a therapy provides to a patient population. Decisions are typically made using a net benefit approach grounded in expected utility theory. The optimal treatment strategy is the one that maximizes expected NMBs where expected NMBs are calculated by averaging over the patient population and uncertain parameters \(\theta\). For a given subgroup \(g\) and parameter set \(\theta\), NMBs are computed as the difference between the monetized health gains from an intervention less costs, or,

\[ \begin{aligned} NMB_g(j,\theta) = e_{gj}\cdot k- c_{gj}, \end{aligned} \]

where \(e_{gj}\) and \(c_{gj}\) are measures of clinical effectiveness (e.g. QALYs) and costs in subgroup \(g\) using treatment \(j\) respectively, and \(k\) is a decision makers willingness to pay per unit of clinical effectiveness. The optimal treatment for a given subgroup is the one that maximizes expected NMBs,

\[ \begin{aligned} j^{*}_g = \text{argmax}_j E_{\theta} \left[NMB_g(j,\theta)\right]. \end{aligned} \]

In practice, new interventions are usually compared to a standard treatment often referred to as the comparator. In these cases, a new treatment in a given subgroup is preferred to the comparator if the expected incremental net monetary benefit (INMB) of the new treatment is positive; that is, treatment 1 is preferred to treatment 0 in subgroup \(g\) if \(E_\theta \left[INMB_g\right] > 0\) where the INMB in a particular subgroup is given by

\[ \begin{aligned} INMB_g(\theta) = NMB_g(j = 1, \theta) - NMB_g(j = 0, \theta). \end{aligned} \] Equivalently, treatment \(1\) is preferred to treatment \(0\) in subgroup \(g\) if the incremental cost-effectiveness ratio (ICER) is greater than the willingness to pay threshold \(k\),

\[ \begin{aligned} k > \frac{c_{g1} - c_{s0}}{e_{g1} - e_{g0}} = ICER_g. \end{aligned} \]

3 Probabilistic sensitivity analysis

Expected NMBs are expected values and ignore parameter uncertainty. This implies that NMBs are uncertain and that optimal treatment strategies may be selected incorrectly. This uncertainty can be quantified using PSA, which uses Bayesian and quasi-Bayesian techniques to estimate the distribution of NMBs given the distribution of the parameters for each treatment strategy

Since the joint distribution of the model parameters cannot be derived analytically (except in the simplest cases), the distribution of \(\theta\) is approximated by simulating the parameters from their joint posterior distribution and calculating relevant quantities of interest as a function of the simulated parameters. For each treatment strategy and subgroup, PSA therefore produces \(n\) random draws from the posterior distribution of clinical effectiveness and costs,

\[ \begin{aligned} e_{gj} &= [e_{gj}^1, e_{gj}^2, \dots, e_{gj}^n] \\ c_{gj} &= [c_{gj}^1, c_{gj}^2, \dots, c_{gj}^n]. \end{aligned} \]

Below we simulate costs and QALYs for three treatment strategies and two subgroups (in a real world analysis, this output would be derived from a detailed health-economic simulation model). Strategy 1 is the current standard of care; it is the cheapest therapy, but also the least efficacious. Strategies 2 and 3 are equally costly, but Strategy 2 is more effective in subgroup 1 while Strategy 3 is more effective in subgroup 2.

set.seed(131)
n_samples <- 1000

# cost
c <- vector(mode = "list", length = 6)
names(c) <- c("Strategy 1, Grp 1", "Strategy 1, Grp 2", "Strategy 2, Grp 1",
              "Strategy 2, Grp 2", "Strategy 3, Grp 1", "Strategy 3, Grp 2")
c[[1]] <- rlnorm(n_samples, 2, .1)
c[[2]] <- rlnorm(n_samples, 2, .1)
c[[3]] <- rlnorm(n_samples, 11, .15)
c[[4]] <- rlnorm(n_samples, 11, .15)
c[[5]] <- rlnorm(n_samples, 11, .15)
c[[6]] <- rlnorm(n_samples, 11, .15)

# effectiveness
e <- c
e[[1]] <- rnorm(n_samples, 8, .2)
e[[2]] <- rnorm(n_samples, 8, .2)
e[[3]] <- rnorm(n_samples, 10, .8)
e[[4]] <- rnorm(n_samples, 10.5, .8)
e[[5]] <- rnorm(n_samples, 8.5, .6)
e[[6]] <- rnorm(n_samples, 11, .6)

# cost and effectiveness by strategy and simulation
library("data.table")
ce <- data.table(sample = rep(seq(n_samples), length(e)),
                 strategy = rep(paste0("Strategy ", seq(1, 3)), 
                                each = n_samples * 2),
                 grp = rep(rep(c("Group 1", "Group 2"),
                               each = n_samples), 3),
                 cost = do.call("c", c), qalys = do.call("c", e))
head(ce)
##    sample   strategy     grp     cost    qalys
## 1:      1 Strategy 1 Group 1 6.808115 7.984456
## 2:      2 Strategy 1 Group 1 6.997936 7.984085
## 3:      3 Strategy 1 Group 1 8.183247 8.221359
## 4:      4 Strategy 1 Group 1 7.423647 8.155250
## 5:      5 Strategy 1 Group 1 7.063846 8.292916
## 6:      6 Strategy 1 Group 1 7.234418 8.066187

For any given willingness to pay \(k\), expected NMBs can be calculated by strategy, subgroup, and simulation number. For example, with \(k=150,000\), a reasonable estimate of the value of a life-year in the United States, Strategy 2 provides the highest expected NMBs in subgroup 2 while Strategy 3 provides the highest expected NMBs in subgroup 2.

ce <- ce[, nmb := 150000 * qalys - cost]
enmb <- ce[, .(enmb = mean(nmb)), by = c("strategy", "grp")]
enmb <- dcast(enmb, strategy ~ grp, value.var = "enmb")
print(enmb)
##      strategy Group 1 Group 2
## 1: Strategy 1 1200672 1200323
## 2: Strategy 2 1447770 1509467
## 3: Strategy 3 1215793 1589370

A number of measures have been proposed in the health economics literature to summarize the PSA. Below we describe the most common measures, which can be calculated using the functions icea and icea_pw. The icea function summarizes results by taking into account each treatment strategy in the analysis, while the function icea_pw summarizes “pairwise” results in which each treatment is compared to a comparator.

We can use the icea function to summarize results from our data.table object of simulated output for a range of willingness to pay values,

library("hesim")
ktop <- 200000
icea <-  icea(ce, k = seq(0, ktop, 500), sample = "sample", strategy = "strategy",
                 grp = "grp", e = "qalys", c = "cost")

The most important input in icea is the data.table object (x) containing columns for simulation number (sim), treatment strategy (strategy), subgroup (grp), clinical effectiveness (e), and costs (c). Users specify the names of the relevant columns in their output table as strings. The other relevant parameter is \(k\), which is a range of willingness to pay values to use for estimating NMBs.

Likewise, we can use icea_pw to summarize the PSA when directly comparing the two treatment strategies (Strategy 2 and Strategy 3) to the comparator (Strategy 1).

icea_pw <-  icea_pw(ce,  k = seq(0, ktop, 500), comparator = "Strategy 1",
                       sample = "sample", strategy = "strategy", grp = "grp",
                    e = "qalys", c = "cost")

The same inputs are used as in icea except users must specify the name of the comparator strategy.

3.1 Probability most cost-effective

A useful summary measure for quantifying uncertainty is the probability that each treatment strategy is the most cost effective. For a particular subgroup, this is estimated from simulation output as the proportion of simulation draws that each strategy has the highest NMB. For example, consider a random sample of 10 draws from the PSA simulation output and suppose \(k\) is again equal to \(150,000\).

random_rows <- sample(1:n_samples, 10)
nmb_dt <- dcast(ce[sample %in% random_rows & grp == "Group 2"], 
                sample ~ strategy, value.var = "nmb")
setnames(nmb_dt, colnames(nmb_dt), c("sample", "nmb1", "nmb2", "nmb3"))
nmb_dt <- nmb_dt[, maxj := apply(nmb_dt[, .(nmb1, nmb2, nmb3)], 1, which.max)]
nmb_dt <- nmb_dt[, maxj := factor(maxj, levels = c(1, 2, 3))]
sample nmb1 nmb2 nmb3 maxj
227 1205465 1402045 1629226 3
353 1193747 1647590 1480650 2
373 1200018 1607944 1623927 3
445 1166575 1578200 1439357 2
497 1222893 1506394 1625707 3
555 1202491 1405106 1457402 3
701 1232418 1405196 1532673 3
770 1193758 1672907 1414983 2
907 1178634 1534693 1666771 3
992 1157191 1528678 1685557 3 
mce <- prop.table(table(nmb_dt$maxj))
print(mce)
## 
##   1   2   3 
## 0.0 0.3 0.7

In this example, treatments 1, 2, and 3 have the highest NMBs a fraction 0, 0.3, and 0.7 of the time respectively. The icea function performs this same calculations for a range of values of \(k\) and all nsims random draws of the simulation output. The output is a tidy data.table which facilitates plotting with ggplot.

library("ggplot2")
library("scales")
theme_set(theme_minimal())
ggplot2::ggplot(icea$mce, aes(x = k, y = prob, col = factor(strategy))) +
  geom_line() + facet_wrap(~grp) + xlab("Willingness to pay") +
  ylab("Probability most cost-effective") +
  scale_x_continuous(breaks = seq(0, ktop, 100000), label = scales::dollar) +
  theme(legend.position = "bottom") + scale_colour_discrete(name = "Strategy")

In group 1, Strategy 2 provides the greatest NMBs with high probability for almost all reasonable values of k. In group 2, the results are less certain, although Strategy 3 provides the greatest NMBs with a higher probability than Strategy 2.

3.2 Cost-effectiveness acceptability frontier (CEAF)

One drawback of the previous plot is that the probability of being cost-effective cannot be used to determine the optimal treatment option. Instead, if a decision-makers objective is to maximize health gain, then decisions should be based on the expected net monetary benefit (Barton, Briggs, and Fenwick 2008). The cost-effectiveness acceptability frontier (CEAF), which plots the probability that the optimal treatment strategy (i.e., the strategy with the highest expected NMB) is cost-effective, is appropriate in this context.

A CEAF curve can be easily created by using the best column to subset to the treatment strategy with the highest expected NMB for each willingness to pay value and group.

ggplot2::ggplot(icea$mce[best == 1], aes(x = k, y = prob, col = strategy)) +
  geom_line() + facet_wrap(~grp) + xlab("Willingness to pay") +
  ylab("Probability most cost-effective") +
  scale_x_continuous(breaks = seq(0, ktop, 100000), label = scales::dollar) +
  theme(legend.position = "bottom") + scale_colour_discrete(name = "Strategy")

3.3 Value of perfect information

One draw back of the previous measure is that it ignores the magnitude of cost or QALY gains. A measure which combines the probability of being most effective with the magnitude of the expected NMB is the expected value of perfect information (EVPI). Intuitively, the EVPI provides an estimate of the amount that a decision maker would be willing to pay to collect additional data and completely eliminate uncertainty. Mathematically, the EVPI is defined as the difference between the maximum expected NMB given perfect information and the maximum expected NMB given current information. In other words, we calculate the NMB for the optimal treatment strategy for each random draw of the parameters and compare that to the NMB for the treatment strategy that is optimal when averaging across all parameters. Mathematically, the EVPI for subgroup \(g\) is,

\[ \begin{aligned} EVPI_g &= E_\theta \left[max_j NMB_g(j, \theta)\right] - max_j E_\theta \left [ NMB_g(j, \theta)\right]. \\ \end{aligned} \]

To illustrate consider the same random sample of 10 draws from our simulation output used above.

strategymax_g2 <- which.max(enmb[[3]])
nmb_dt <- nmb_dt[, nmbpi := apply(nmb_dt[, .(nmb1, nmb2, nmb3)], 1, max)]
nmb_dt <- nmb_dt[, nmbci := nmb_dt[[strategymax_g2 + 1]]]
kable(nmb_dt, digits = 0, format = "html")
sample nmb1 nmb2 nmb3 maxj nmbpi nmbci
227 1205465 1402045 1629226 3 1629226 1629226
353 1193747 1647590 1480650 2 1647590 1480650
373 1200018 1607944 1623927 3 1623927 1623927
445 1166575 1578200 1439357 2 1578200 1439357
497 1222893 1506394 1625707 3 1625707 1625707
555 1202491 1405106 1457402 3 1457402 1457402
701 1232418 1405196 1532673 3 1532673 1532673
770 1193758 1672907 1414983 2 1672907 1414983
907 1178634 1534693 1666771 3 1666771 1666771
992 1157191 1528678 1685557 3 1685557 1685557


To calculate EVPI, we average NMBs given current information and NMBs given perfect information across simulation draws.

enmbpi <- mean(nmb_dt$nmbpi)
enmbci <- mean(nmb_dt$nmbci)
print(enmbpi)
## [1] 1611996
print(enmbci)
## [1] 1555625
print(enmbpi - enmbci)
## [1] 56370.66

The icea function performs this same calculation across all simulation draws from the PSA and for a number of values of willingness to pay values \(k\). A plot by group of the the EVPI for different values of \(k\) is shown below. The kinks in the plot represent values of \(k\) where the optimal strategy changes.

ggplot2::ggplot(icea$evpi, aes(x = k, y = evpi)) +
  geom_line() + facet_wrap(~grp) + xlab("Willingness to pay") +
  ylab("Expected value of perfect information") +
  scale_x_continuous(breaks = seq(0, ktop, 100000), label = scales::dollar) +
  scale_y_continuous(label = scales::dollar) +
  theme(legend.position = "bottom")