1 Overview

Probabilistic sensitivity analysis (PSA) is used to quantify the impact of parameter uncertainty on the uncertainty of model outputs. PSA is typically performed via a simulation approach whereby the model parameters are randomly sampled from suitable probability distributions and the entire model is simulated for each random draw of the parameters.

In this example, we extend the deterministic simple Markov cohort model to incorporate PSA. We will continue to rely primarily on Base R but will use the hesim package to help make the code more readable.

library("rcea")
library("hesim")
library("data.table")
library("magrittr")
library("ggplot2")

2 Model parameters

2.1 Transition probabilities for SOC

The probability distribution used for transition probabilities will depend on the underlying data. In this case, we assume that summary level data is available on transitions from the Healthy state (n = 900), Sick state (n = 900), and Sicker state (n = 800). The transitions from each state to the other 4 states can be modeled using a Dirichlet distribution (see Appendix).

transitions_soc <- matrix(
  c(848, 150, 0,   2,
    500, 389, 105,  6,
    0,   0,   784, 16,
    0,   0,   0,   23),
  nrow = 4, byrow = TRUE)
state_names <- c("H", "S1", "S2", "D")
colnames(transitions_soc) <- rownames(transitions_soc) <- tolower(state_names)

2.2 Relative risk

We estimate treatment effects in terms of the log relative risk since it is approximately normally distributed; that is, the relative risk follows a lognormal distribution. The mean is given by the log of the point estimate of the relative risk, or \(log(0.8)\).

Since academic studies often report 95% confidence intervals for a parameter, but not its standard error, we will assume that is the case. Specifically, let the lower bound be \(log(0.71)\) and the upper bound be \(log(0.91)\). The standard error is then given by \((log(0.91) - log(0.71))/2z\) where \(z = \Phi^{-1}(0.975)\approx 1.96\). See the Appendix for details.

2.3 Costs

2.3.1 Medical costs

Medical costs are of often assumed to follow a gamma distribution because it can be used to model right skewed distributions. The gamma distribution is defined in terms of a shape and scale (or rate) parameter. However, these parameters can be derived form the mean and standard deviation of a distribution using the method of moments. Mean costs follow those in the deterministic example (H = $2,000, S1 = $4,000, S2 = $2,000), D = \(0\)). The standard deviation in each state is assumed to be equal to the mean.

2.3.2 Treatment costs

Treatment costs are fixed and equal to the costs in the deterministic example.

2.4 Utility

The utility associated with each state is assumed to follow a beta distribution, which is bounded between 0 and 1. The beta distribution is defined in terms of 2 shape parameters, but like the gamma distribution , these can be derived using the method of moments. We assume that the mean (standard error) of utility is estimated to be \(1\) \((0.0)\), \(0.75\) \((0.03)\), \(0.5\) \((0.05)\), and \(0\) \((0.0)\) in state H, S1, S2, and D, respectively.

2.5 Combining the parameters

All model parameters (transition probabilities, relative risk, costs, and utility) can be stored in a list for use in the simulation.

params <- list(
  alpha_soc = transitions_soc,
  lrr_mean = log(.8),
  lrr_lower = log(.71),
  lrr_upper = log(.9),
  c_medical = c(H = 2000, S1 = 4000, S2 = 15000, D = 0),
  c_soc = 2000,
  c_new = 12000,
  u_mean = c(H = 1, S1 = .75, S2 = 0.5, D = 0),
  u_se = c(H = 0, S1 = 0.03, S2 = 0.05, D = 0.0)
)

3 Simulation

The simulation proceeds by (i) randomly sampling the parameters from the probability distributions specified above and (ii) running the Markov model for each draw of the parameters. The result is a draw from the probability distribution of each of the model outputs of interest (i.e, state probabilities, QALYs, and costs).

3.1 Sampling the parameters

While Base R can certainly be used to draw samples of the parameters, the functions hesim::define_rng() and hesim::eval_rng() simplify this process and make the code more readable. Any random number generation function can be used inside the define_rng() block; the only rule is that returned parameter draws must be returned as a list. However, hesim comes with a number of helpful probability distribution functions (type ?rng_distributions for details) to make your life easier.

rng_def <- define_rng({
  lrr_se <- (lrr_upper - lrr_lower)/(2 * qnorm(.975)) # Local object 
  # not returned
  list( # Parameters to return
    p_soc = dirichlet_rng(alpha_soc),
    rr_new = lognormal_rng(lrr_mean, lrr_se),
    c_medical = gamma_rng(mean = c_medical, sd = c_medical),
    c_soc = c_soc,
    c_new = c_new,
    u = beta_rng(mean = u_mean, sd = u_se)
  )
}, n = 1000)
params_rng <- eval_rng(rng_def, params = params)
attr(params_rng, "n") <- rng_def$n
names(params_rng)
## [1] "p_soc"     "rr_new"    "c_medical" "c_soc"     "c_new"     "u"
head(as.matrix(params_rng$p_soc))
##            h_h      h_s1 h_s2         h_d      s1_h     s1_s1     s1_s2
## [1,] 0.8383229 0.1593326    0 0.002344545 0.4811612 0.3956502 0.1146471
## [2,] 0.8501859 0.1478452    0 0.001968814 0.4669832 0.4240769 0.1047213
## [3,] 0.8407569 0.1559604    0 0.003282725 0.5066133 0.3735621 0.1126371
## [4,] 0.8586012 0.1399014    0 0.001497375 0.5017253 0.3801266 0.1107963
## [5,] 0.8152249 0.1825510    0 0.002224100 0.4787484 0.3995763 0.1136843
## [6,] 0.8533603 0.1449687    0 0.001671050 0.5012932 0.3877513 0.1053192
##             s1_d s2_h s2_s1     s2_s2       s2_d d_h d_s1 d_s2 d_d
## [1,] 0.008541637    0     0 0.9787120 0.02128797   0    0    0   1
## [2,] 0.004218618    0     0 0.9819925 0.01800754   0    0    0   1
## [3,] 0.007187433    0     0 0.9837840 0.01621598   0    0    0   1
## [4,] 0.007351737    0     0 0.9765105 0.02348948   0    0    0   1
## [5,] 0.007990968    0     0 0.9794340 0.02056597   0    0    0   1
## [6,] 0.005636229    0     0 0.9778453 0.02215469   0    0    0   1

3.2 Simulating the Markov model

Once samples of the parameters have been drawn, the Markov model can be simulated for each draw.

3.2.1 Input data

One way that a Markov simulation can be generalized is to store “input data” in an object such as a data frame. This input data might consist of different treatment strategies, patients and/or subgroups, or even information about the health states themselves. For instance, if we were simulating different subgroups we might store the age and sex associated with the subgroup which could, in turn, be used as covariates in a statistical model. In this simple example the data will just consist of the names of the two treatment strategies.

data <- data.frame(
  strategy = c("New", "SOC")
)
head(data)
##   strategy
## 1      New
## 2      SOC

3.2.2 Running the simulation

It is a good idea to modularize your R code into functions. This can make your code more readable, maintainable, and reusable. We will work toward a sim_model() function that runs the entire simulation. It will be comprised of three smaller functions: sim_stateprobs(), compute_qalys(), and compute_costs().

The sim_stateprobs() function will simulate health state probabilities for each model cycle (i.e., the Markov trace) for a given treatment strategy and parameter sample. It takes the arguments:

  • p0: the transition probability matrix for SOC,
  • rr: the relative risk)
  • strategy: the name of the treatment strategy as defined in data above.
  • n_cycles: The number of cycles to simulate the model for.

To make the code more readable, we will use the function hesim::tpmatrix() which makes it easy to define a transition probability matrix. The symbol C denotes that a given element is the complement of all other elements in that row, ensuring that the probabilities sum to 1.

sim_stateprobs <- function(p0, rr, strategy, n_cycles){
  rr <- ifelse(strategy == "New", rr, 1)
  p <- tpmatrix(
    C,       p0$h_s1 * rr,  p0$h_s2 * rr,  p0$h_d * rr,
    p0$s1_h, C,             p0$s1_s2 * rr, p0$s1_d * rr,
    p0$s2_h, p0$s2_s1,      C,             p0$s2_d * rr,
    0,       0,             0,             1
  )
  x <- sim_markov_chain(x0 = c(1, 0, 0, 0),
                        p = matrix(as.matrix(p), ncol = 4, byrow = TRUE),
                        n_cycles = n_cycles)
  return(x)
}

To compute (discounted) QALYs for a given Markov trace, we will use a very simple function that is identical to the compute_qalys function used in the deterministic example.

compute_qalys <- function(x, utility, dr = .03){
  n_cycles <- nrow(x) - 1
  pv(x %*% utility, dr, 0:n_cycles)
}

Similarly, our cost function is nearly identical to the compute_costs() function from the deterministic example.

compute_costs <- function(x, costs_medical, costs_treat, dr = .03){
  n_cycles <- nrow(x) - 1
  costs_treat <- c(rep(costs_treat, 3), 0)
  costs <- cbind(
    pv(x %*% costs_medical, dr, 0:n_cycles),
    pv(x %*% costs_treat, dr, 0:n_cycles)
  )
  colnames(costs) <- c("dcost_med", "dcost_treat")
  return(costs)
}

Now that we have created the building blocks for the simulation, we can create our main function to simulate the entire model. 3% discount rates are set by default for costs and QALYs. Most of the function arguments should be self explanatory, but two are worth explaining:

  • params_rng: The output of eval_rng() above.
  • data: The data object defined above

The first part of the function creates an array to store the output. The array is a series of matrices each with n_cycles rows and columns for each output (i.e., state probabilities for the four health states, QALYs, treatment costs, and medical costs). There is one matrix for each parameter sample for the PSA and treatment strategy.

The second part of the function simulates state probabilities (with sim_stateprobs()), QALYs (with compute_qalys()), and costs (with compute_costs()) for each parameter sample and treatment strategy. The number of parameter samples and the names of the treatment strategies are saved as attributes (i.e., metadata) to the array which will be used below to convert the array to a data frame.

sim_model <- function(params_rng, data, n_cycles = 85,
                      dr_qalys = .03, dr_costs = .03){
  # Initialize array of matrices
  n_samples <- attr(params_rng, "n")
  n_strategies <- nrow(data)
  out <- array(NA, dim = c(n_cycles + 1, 7, n_samples * n_strategies))
  dimnames(out) <- list(NULL,
                        c("H", "S1", "S2", "D",
                          "dqalys", "dcosts_med", "dcosts_treat"),
                        NULL)

  # Run the simulation
  i <- 1
  for (s in 1:n_samples){ # Start PSA loop
    for (k in 1:n_strategies) { # Start treatment strategy loop
      x <- sim_stateprobs(p0 = params_rng$p_soc[s, ],
                          rr = params_rng$rr_new[s],
                          strategy = data$strategy[k],
                          n_cycles = n_cycles)
      dqalys <- compute_qalys(x, utility = unlist(params_rng$u[s]),
                              dr = dr_qalys)
      dcosts <- compute_costs(x,
                              costs_medical = unlist(params_rng$c_medical[s]),
                              costs_treat = ifelse(data$strategy[k] == "SOC",
                                                   params_rng$c_soc,
                                                   params_rng$c_new),
                              dr = dr_costs)
      out[, , i] <- cbind(x, dqalys, dcosts)
      i <- i + 1
    } # End treatment strategy loop
  } # End PSA loop

  # Store metadata and return
  attr(out, "n_samples") <- n_samples
  attr(out, "strategies") <- data$strategy
  return(out)
}

Now that we’ve written the function, lets simulate the model with argument defaults (85 model cycles and 3% discount rates). Recall that each array is a matrix.

sim_out <- sim_model(params_rng, data = data)
head(sim_out[, , 1])
##              H        S1          S2           D    dqalys dcosts_med
## [1,] 1.0000000 0.0000000 0.000000000 0.000000000 1.0000000   3034.401
## [2,] 0.8811154 0.1177325 0.000000000 0.001152140 0.9449958   3483.180
## [3,] 0.8349810 0.1529612 0.009056469 0.003001344 0.9040508   3550.256
## [4,] 0.8118710 0.1622592 0.020722824 0.005147011 0.8683033   3502.118
## [5,] 0.7961377 0.1634260 0.032975537 0.007460803 0.8350131   3420.179
## [6,] 0.7828557 0.1620615 0.045182644 0.009900094 0.8033110   3329.309
##      dcosts_treat
## [1,]     12000.00
## [2,]     11637.06
## [3,]     11277.20
## [4,]     10925.18
## [5,]     10582.30
## [6,]     10248.83

Although arrays are computationally efficient objects for storing data, they aren’t often the most useful for summarizing data. We will write to short functions to convert a 3D array to a data.table (with ID columns for the parameter sample and treatment strategy) so that we can summarize outcomes for each parameter sample and treatment strategy very quickly. (Note that other packages such as dplyr could also be used but we prefer data.table for simulations because of its speed).

# rbind an array of matrices into a single matrix
rbind_array <- function(x){
  n_rows <- dim(x)[3] * dim(x)[1]
  x_mat <- matrix(c(aperm(x, perm = c(2, 1, 3))),
                  nrow = n_rows, byrow = TRUE)
  colnames(x_mat) <- dimnames(x)[[2]]
  return(x_mat)
}

# Convert the array into a long dataframe with ID columns
array_to_dt <- function(x){
  id_df <- expand.grid(cycle = 0:(dim(x)[1] - 1),
                       strategy = attr(x, "strategies"),
                       sample = 1:attr(x, "n_samples"))
  x_mat <- rbind_array(x)
  return(as.data.table(cbind(id_df, x_mat)))
}

sim_out <- array_to_dt(sim_out)
head(sim_out)
##    cycle strategy sample         H        S1          S2           D    dqalys
## 1:     0      New      1 1.0000000 0.0000000 0.000000000 0.000000000 1.0000000
## 2:     1      New      1 0.8811154 0.1177325 0.000000000 0.001152140 0.9449958
## 3:     2      New      1 0.8349810 0.1529612 0.009056469 0.003001344 0.9040508
## 4:     3      New      1 0.8118710 0.1622592 0.020722824 0.005147011 0.8683033
## 5:     4      New      1 0.7961377 0.1634260 0.032975537 0.007460803 0.8350131
## 6:     5      New      1 0.7828557 0.1620615 0.045182644 0.009900094 0.8033110
##    dcosts_med dcosts_treat
## 1:   3034.401     12000.00
## 2:   3483.180     11637.06
## 3:   3550.256     11277.20
## 4:   3502.118     10925.18
## 5:   3420.179     10582.30
## 6:   3329.309     10248.83

4 Cost-effectiveness analysis

A cost-effectiveness analysis (CEA) can be performed using the simulation output. A PSA is typically used to represent decision uncertainty using the distribution of (discounted) QALYs and (discounted) total costs. As such, will compute mean discounted QALYs and discounted costs by parameter sample and treatment strategy. As in the previous tutorial, we assume that transitions occur immediately and therefore exclude costs and QALYs measured at the start of the first model cycle.

ce_out <- sim_out[cycle != 0,
                  .(dqalys = sum(dqalys),
                    dcosts = sum(dcosts_med) + sum(dcosts_treat)),
                  by = c("sample", "strategy")]
ce_out
##       sample strategy   dqalys   dcosts
##    1:      1      New 23.76522 446737.9
##    2:      1      SOC 20.75515 164378.8
##    3:      2      New 23.36452 379221.4
##    4:      2      SOC 21.44274 112713.1
##    5:      3      New 23.45367 479096.7
##   ---                                  
## 1996:    998      SOC 21.04841 165033.7
## 1997:    999      New 23.31300 438144.9
## 1998:    999      SOC 21.72958 147964.9
## 1999:   1000      New 23.08725 383215.4
## 2000:   1000      SOC 21.14566 114027.8

In the two treatment strategy case its simpler to “widen” the data so that we can easily compute incremental QALYs and incremental costs.

ce_out_wider <- dcast(ce_out, sample ~ strategy,
                      value.var = c("dqalys", "dcosts"))
ce_out_wider
##       sample dqalys_New dqalys_SOC dcosts_New dcosts_SOC
##    1:      1   23.76522   20.75515   446737.9  164378.76
##    2:      2   23.36452   21.44274   379221.4  112713.14
##    3:      3   23.45367   21.15014   479096.7  206900.16
##    4:      4   23.22162   21.21390   371816.7  105464.51
##    5:      5   22.30074   20.67493   361265.5   93663.15
##   ---                                                   
##  996:    996   23.51433   20.90575   365356.6   90349.23
##  997:    997   22.56368   20.16734   433310.5  172249.36
##  998:    998   22.92116   21.04841   430276.9  165033.73
##  999:    999   23.31300   21.72958   438144.9  147964.87
## 1000:   1000   23.08725   21.14566   383215.4  114027.77

The ICER is computed by taking means across all parameter samples.

ce_out_wider[, idcosts := dcosts_New - dcosts_SOC]
ce_out_wider[, idqalys := dqalys_New - dqalys_SOC]
ce_out_wider[, .(icer = mean(idcosts)/mean(idqalys))]
##        icer
## 1: 124582.4

We can also use a cost-effectiveness plane to visually represent uncertainty in incremental QALYs and incremental costs. The dotted line represents willingness to pay for a QALY (in this example, $100,000) and points below the line are cost-effective while points above the line are not.

A more formal and thorough treatment of CEA and decision uncertainty is given in the CEA tutorial.

5 Appendix

5.1 Dirichlet distribution

The multinomial distribution is a discrete probability distribution for the number of successes for each of k mutually exclusive categories in n trials. The probabilities of the categories are given by \(\pi_1,\ldots, \pi_k\) with \(\sum_{j=1}^k \pi_j=1\) and each \(\pi_j\) defined on \([0,1]\). The Dirichlet distribution is parameterized by the concentration parameters \(\alpha_1,\ldots, \alpha_k\) with \(\alpha_j > 0\). Letting \(x_1,\ldots, x_k\) denote the number of successes in each category, the prior distribution and likelihood are,

\[ \begin{aligned} p(\pi_1,\ldots,\pi_k |\alpha_1,\ldots, \alpha_k) = \text{Dirichlet}(\alpha_1,\ldots,\alpha_k) \\ p(x_1,\ldots,x_k | \pi_1,\ldots,\pi_k) = \text{Multin}(n, \pi_1,\ldots,\pi_k). \end{aligned} \]

5.2 Confidence intervals and standard errors

Let \(\theta\) be a normally distributed random variable, \(\sigma\) be the standard deviation of \(\theta\), and \(z = \Phi^{-1}(1 - \alpha/2)\) for a given confidence level \(\alpha\). A confidence interval for \(\theta\) is then given by \((\theta - z\sigma, \theta + z\sigma)\) with width given by \(2z\sigma\). If lower and upper limits are given by \(\theta_L\) and \(\theta_U\), \(\sigma\) can be solved with

\[ \sigma = \frac{\theta_U - \theta_L}{2z}. \]

A 95% confidence interval is evaluated with \(z = \Phi^{-1}(.975)\).