Overview

This tutorial repeats the probabilistic sensitivity analysis (PSA) of the Markov cohort model simulation performed in the previous tutorial using hesim. We utilize the cohort discrete time state transition model (cDTSTM) class, which is another name for a (time-homogeneous or time-inhomogeneous) Markov cohort model.

More information about hesim can be found by visiting the package website. We recommend reading the “Articles”—starting with the “Introduction to hesim—to learn more. Economic models can, in general, be simulated with the following steps:

  1. Model setup: Specify the treatment strategies, target population(s), and model structure.
  2. Parameters: Estimate or define the parameters of the economic model.
  3. Simulation:
    1. Construction of model: Create an economic model—consisting of separate statistical models for disease progression, costs, and utilities—that simulate outcomes as a function of input data (derived from Step 1) and parameters (from Step 2).
    2. Simulation of outcomes: Simulate outcomes (disease progression, costs, and quality-adjusted life-years (QALYs)) using the model constructed in Step 3.

This analysis can be performed using the hesim package alone.

Model setup

Before beginning an analysis, it is necessary to define the treatment strategies of interest and the the target population of interest. We continue to run an analysis for two treatment strategies (“SOC” and the “New” treatment) and one representative 25-year old patient.

strategies <- data.frame(
  strategy_id = 1:2,
  strategy_name = c("SOC", "New")
)
patients <- data.frame(
  patient_id = 1,
  age = 25
)
hesim_dat <- hesim_data(
  strategies = strategies,
  patients = patients
)
print(hesim_dat)
## $strategies
##   strategy_id strategy_name
## 1           1           SOC
## 2           2           New
## 
## $patients
##   patient_id age
## 1          1  25
## 
## attr(,"class")
## [1] "hesim_data"
labs <- get_labels(hesim_dat)

Model parameters

We will use the same list of parameters as in the previous tutorial.

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),
  c_soc = 2000,
  c_new = 12000,
  u_mean = c(H = 1, S1 = .75, S2 = 0.5),
  u_se = c(H = 0, S1 = 0.03, S2 = 0.05)
)

Random number generation

The random number generation process for the PSA is once again defined using define_rng().

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)

Transformed parameters

Economic models in hesim consist of separate statistical models used to simulate disease progression, costs, and utility. Like any prediction model, each statistical model simulates outcomes as a function of input data and model parameters. In most statistical models, the underlying parameters (see hesim::params) are transformed (see hesim::tparams) into more relevant parameters for prediction (e.g., the coefficients and covariates from a Weibull regression are used to predict shape and scale parameters).

There are multiple ways to “transform” parameters in hesim. In this example, we will use mathematical expressions via hesim::define_tparams() (designed for Markov cohort models), whereas the next tutorial on multi-state modeling will estimate parameters (see hesim::params_surv) by fitting parametric survival models.

Transformed parameters are modeled here as function of treatment strategies and patients, which can be generated using hesim::expand().

input_data <- hesim::expand(hesim_dat, by = c("strategies", "patients"))
head(input_data)
##    strategy_id patient_id strategy_name age
## 1:           1          1           SOC  25
## 2:           2          1           New  25

Operations in a define_tparams() block are performed using the columns from the input data. To maximize computational efficiency, all operations are vectorized across the rows in the data and PSA samples. Parameters not included in a transformed parameter function are assumed constant across patients and treatment strategies. A list must be returned with the following possible elements:

  • tpmatrix: The transition probability matrix used to simulate transition probabilities in the economic model.
  • utility: The utility values to attach to states and used to simulate QALYs in the economic model. Either a vector (in which case utility is the same in each health state) or a matrix with a column for each (non-death) health state.
  • costs: A named list of costs for each category used to simulate costs in the economic model. Each element of the list must be in the same format as utility.
tparams_def <- define_tparams({
  ## The treatment effect (relative risk) is transformed so that it varies by 
  ## strategies (SOC is the reference strategy)
  rr <- ifelse(strategy_name == "SOC", 1, rr_new)
  
  list(
    tpmatrix = tpmatrix(
      C,          p_soc$h_s1 * rr, p_soc$h_s2 * rr,  p_soc$h_d * rr,
      p_soc$s1_h, C,               p_soc$s1_s2 * rr, p_soc$s1_d * rr,
      p_soc$s2_h, p_soc$s2_s1,     C,                p_soc$s2_d * rr,
      0,          0,               0,                1
    ),
    utility = u,
    costs = list(
      treatment = ifelse(strategy_name == "SOC", c_soc, c_new),
      medical = c_medical
    )
  )
})

Simulation

Construct the model

The economic model is defined by using define_model() to combine the underlying parameters with the expressions for random number generation and parameter transformation.

mod_def <- define_model(tparams_def = tparams_def, 
                        rng_def = rng_def, 
                        params = params)

An economic model (of class CohortDtstm) can then be created from the defined model and input data using the generic function create_CohortDtstm(). The economic model is an R6 object consisting of a transition model for simulating transition probabilities with sim_stateprobs(), a utility model for simulating QALYs with sim_qalys(), and a set of cost models (for each cost category) for simulating costs with sim_costs().

econmod <- create_CohortDtstm(mod_def, input_data)

Simulating outcomes

Health state probabilities

State occupancy probabilities are generated by simulating the discrete time Markov chain.

econmod$sim_stateprobs(n_cycles = 85)
head(econmod$stateprobs_)
##    sample strategy_id patient_id grp_id state_id t      prob
## 1:      1           1          1      1        1 0 1.0000000
## 2:      1           1          1      1        1 1 0.8703459
## 3:      1           1          1      1        1 2 0.8250446
## 4:      1           1          1      1        1 3 0.8024778
## 5:      1           1          1      1        1 4 0.7861723
## 6:      1           1          1      1        1 5 0.7717218

Costs and QALYs

Costs and QALYs are computed by integrating the previously simulated state probabilities. To maintain consistency with our prior analyses, we measure costs and QALYs at the start of each model cycle using a left Riemann sum. We use the option lys = TRUE so that life-years are simulated in addition to QALYs.

econmod$sim_qalys(dr = 0.03, lys = TRUE, integrate_method = "riemann_right")
head(econmod$qalys_)
##    sample strategy_id patient_id grp_id state_id   dr     qalys       lys
## 1:      1           1          1      1        1 0.03 16.980147 16.980147
## 2:      1           1          1      1        2 0.03  2.574968  3.543340
## 3:      1           1          1      1        3 0.03  2.958624  6.016953
## 4:      1           2          1      1        1 0.03 19.525390 19.525390
## 5:      1           2          1      1        2 0.03  2.435161  3.350955
## 6:      1           2          1      1        3 0.03  2.372067  4.824072
econmod$sim_costs(dr = 0.03, integrate_method = "riemann_right")
head(econmod$costs_)
##    sample strategy_id patient_id grp_id state_id   dr  category     costs
## 1:      1           1          1      1        1 0.03 treatment  33960.29
## 2:      1           1          1      1        2 0.03 treatment   7086.68
## 3:      1           1          1      1        3 0.03 treatment  12033.91
## 4:      1           2          1      1        1 0.03 treatment 234304.68
## 5:      1           2          1      1        2 0.03 treatment  40211.46
## 6:      1           2          1      1        3 0.03 treatment  57888.87

Cost-effectiveness analysis

Cost-effectiveness analyses (CEAs) can be performed directly from the simulation output with hesim. Other R packages such as BCEA could also be considered. To perform the CEA, we will first need to sum costs and QALYs across health states.

ce_sim <- econmod$summarize()
saveRDS(ce_sim,  file = "markov-cohort-hesim-ce_sim.rds")
saveRDS(hesim_dat,  file = "markov-cohort-hesim_data.rds")

The analysis can then be performed directly from this output using the hesim::cea() and hesim::cea_pw() functions. The former considers all treatment strategies simultaneously while the latter makes pairwise comparisons of each treatment strategies to a comparator (in our example the SOC). A complete analysis will be provided in the CEA tutorial using the output saved above; for now, we perform a simple pairwise analysis:

cea_pw_out <- cea_pw(ce_sim, comparator = 1, 
                     dr_qalys = 0.03, dr_costs = 0.03,
                     k = seq(0, 25000, 500))

While cea_pw() allows users to summarize output from a PSA we will just create an ICER table using means for now.

format(icer(cea_pw_out, k = 50000, labels = labs))
##              Outcome                          New
## 1: Incremental QALYs            2.05 (1.05, 3.07)
## 2: Incremental costs   257,393 (200,711, 288,833)
## 3:   Incremental NMB -155,092 (-210,411, -75,154)
## 4:              ICER                      125,802