vignettes/03-markov-cohort-hesim.Rmd
03-markov-cohort-hesim.RmdThis 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.
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)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 25Operations 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.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.7717218Costs 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.87Cost-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.
## 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