# 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

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

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)$$.