vignettes/03-markov-cohort-hesim.Rmd
03-markov-cohort-hesim.Rmd
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:
This analysis can be performed using the hesim
package alone.
library("hesim")
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"
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) )
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)
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()
.
## 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:
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 ) ) })
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)
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.8505559
## 3: 1 1 1 1 1 2 0.7980505
## 4: 1 1 1 1 1 3 0.7722735
## 5: 1 1 1 1 1 4 0.7540321
## 6: 1 1 1 1 1 5 0.7380763
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.
## sample strategy_id patient_id grp_id state_id dr qalys lys
## 1: 1 1 1 1 1 0.03 15.651309 15.651309
## 2: 1 1 1 1 2 0.03 2.980147 3.937807
## 3: 1 1 1 1 3 0.03 2.992457 5.878751
## 4: 1 2 1 1 1 0.03 18.300481 18.300481
## 5: 1 2 1 1 2 0.03 2.864980 3.785632
## 6: 1 2 1 1 3 0.03 2.471218 4.854765
## sample strategy_id patient_id grp_id state_id dr category costs
## 1: 1 1 1 1 1 0.03 treatment 31302.618
## 2: 1 1 1 1 2 0.03 treatment 7875.615
## 3: 1 1 1 1 3 0.03 treatment 11757.503
## 4: 1 2 1 1 1 0.03 treatment 219605.766
## 5: 1 2 1 1 2 0.03 treatment 45427.587
## 6: 1 2 1 1 3 0.03 treatment 58257.185
Cost-effectiveness analyses (CEAs) can be performed directly from the simulation output with hesim
. Other R packages such as BCEA could also be considered. Here, we will consider a pairwise comparison between the new treatment and SOC with the hesim::cea_pw()
function.
ce_sim <- econmod$summarize() cea_pw_out <- cea_pw(ce_sim, comparator = 1, dr_qalys = 0.03, dr_costs = 0.03, k = seq(0, 25000, 500))
Although cea_pw()
allows users to summarize output from a PSA we will just create an ICER table using means for now. A complete analysis is provided in the CEA tutorial.
icer_tbl(cea_pw_out, colnames = strategies$strategy_name)
## SOC New
## Incremental QALYs "-" "2.04 (0.93, 3.04)"
## Incremental costs "-" "258,119 (203,020, 287,151)"
## Incremental NMB "-" "-156,284 (-216,660, -85,946)"
## ICER "-" "126,734"
## Conclusion "-" "Not cost-effective"