3.1 Estimates from the literature

Counts of transitions at yearly intervals between each of the 4 health states were available from a longitudinal dataset from Chelsea and Westminister hospital.

trans_mono <- matrix(c(1251, 350, 116, 17,
                       0, 731, 512, 15,
                       0, 0, 1312, 437,
                       0, 0, 0, 469),
                      ncol = 4, nrow = 4, byrow = TRUE)
colnames(trans_mono) <- rownames(trans_mono) <-  c("A", "B", "C", "D")
print(trans_mono)

##      A   B    C   D
## A 1251 350  116  17
## B    0 731  512  15
## C    0   0 1312 437
## D    0   0    0 469

The transition probabilities from state \(r\) to each of the 4 states (State A, State B, State C, and State D) can be modeled using Dirichlet distributions, which are parameterized by a vector of concentration parameters \(\alpha = (\alpha_{r1}, \alpha_{r2}, \alpha_{r3}, \alpha_{r4})\). The Dirichlet distribution is particularly useful for modeling uncertainty in the transition probabilities because it is a conjugate prior for the multinomial distribution. The posterior distribution of \(\alpha\) consequently follows a Dirichlet distribution which can be thought of as the “effective sample size” for each transition (see Appendix).

Treatment effects are defined in terms of the relative risk, which is assumed to reduce the probability of all transitions to a more severe health state by an equal amount. We work with the log relative risk since it is approximately normally distributed. The estimated relative risk of disease progression is 0.509 (95% CI 0.365 to 0.710). We can then calculate the mean and standard error on the log scale.

The cost-effectiveness analysis is based on life-years, meaning that each health state is assigned a utility value of 1. The costs of zidovudine and lamivudine are $2,278 and $2,086.50, respectively. Medical costs are split into direct medical and community medical costs.

The parameter estimates are stored in a list.

params <- list(
  alpha_mono = trans_mono, 
  lrr_mean = log(.509), 
  lrr_lower <- log(.365),
  lrr_upper = log(.710),
  c_dmed_mean = c(A = 1701, B = 1774, C = 6948),
  c_cmed_mean = c(A = 1055, B = 1278, C = 2059),
  c_zido = 2278,
  c_lam = 2086.50,
  u = 1
)

3.2 Random number generation

To perform a probabilistic sensitivity analysis (PSA), it is necessary to draw samples of the parameters from suitable probability distributions. While the probability distributions could have been drawn manually and stored in the params object above, it can be more easily done by using random number generation functions with define_rng().

As noted above, we draw the transition probabilities with monotherapy from a Dirichlet distribution and the relative risk (of a transition to a worse state with combination therapy relative to monotherapy) from a lognormal distribution. Drug costs and utility are assumed to be fixed while medical costs are drawn from gamma distributions.

Sampled parameters values are stored in a list, meaning that samples for each parameter can be stored as either a vector or data.table (or matrix or data.frame). For instance, we sample medical costs for each (non-death) health state from gamma distributions and they a stored in a data.table (one row for each sample and one column for each health state). Likewise, the parameters of each cell in the transition probability matrix are store in a data.table, with parameters of the data.table ordered rowwise across the columns. On the other hand, samples of the relative risk are stored in a vector.

rng_def <- define_rng({
  lrr_se <- (lrr_upper - lrr_lower)/(2 * qnorm(.975)) # Local object 
                                                      # not returned
  list( # Parameters to return
    p_mono = dirichlet_rng(alpha_mono),
    rr_comb = lognormal_rng(lrr_mean, lrr_se),
    c_zido = c_zido,
    c_lam = c_lam,
    c_dmed = gamma_rng(mean = c_dmed_mean, sd = c_dmed_mean),
    c_cmed = gamma_rng(mean = c_cmed_mean, sd = c_cmed_mean),
    u = u
  )
}, n = 1000)

3.3 Transformed parameters

In most statistical models, the underlying parameters are “transformed” into more relevant parameters for prediction. For instance, if a multinomial logistic regression is used parameterize transition probabilities, the coefficients are used to predict the probabilities. While transition, cost, and utility models can be parameterized using supported statistical models with hesim::params objects, they can also, in some cases, be parameterized using the transformed parameters themselves with hesim::tparams objects.

We will use the transformed parameter approach here by using define_tparams() to define tparams objects as a function of the underlying parameters and input data (treatment strategies and patients). We generate the input data using hesim::expand().

input_data <- expand(hesim_dat, by = c("strategies", "patients"))
head(input_data)

##    strategy_id patient_id       strategy_name
## 1:           1          1         Monotherapy
## 2:           2          1 Combination therapy

Operations in a define_tparams() block are performed using the columns from the input data and a list of parameters . 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 elements for the transition probability matrix, costs, and utilities.

An important helper function is tpmatrix(), which allows users to define transition matrices as a function of other parameters. In this example, transition probabilities are defined in terms of the transition matrix for monotherapy (p_mono) and the relative risk (rr). The relative risk is simply 1 for monotherapy (since the corresponding transition matrix is just p_mono) and equal to rr_comb for combination therapy.

The only costs in the model that vary across treatment strategies (rather than just states) are drug costs. The monotherapy treatment strategy include the cost of zidovudine only whereas the combination therapy includes the costs of both zidovudine and lamivudine.

Parameters may depend on time (since the start of the model) by using the package defined time variable. The time variable corresponds to the stopping time of time intervals, which are specified using the times argument in define_tparams(). Since patients are only assumed to use lamivudine for two years, rr_comb and the costs of lamivudine are not applied beyond year 2.

tparams_def <- define_tparams({
  ## The treatment effect (relative risk) is transformed so that it varies by 
  ## strategies and only applies for the first 2 years (Monotherapy is 
  ## the reference strategy). 
  rr <- ifelse(strategy_name == "Monotherapy" | time > 2, 1, rr_comb)
  
  list(
    tpmatrix = tpmatrix(
      C, p_mono$A_B * rr, p_mono$A_C * rr, p_mono$A_D * rr,
      0, C,               p_mono$B_C * rr, p_mono$B_D * rr,
      0, 0,               C,               p_mono$C_D * rr,
      0, 0,               0,               1
    ),
    utility = u,
    costs = list(
        drug = ifelse(strategy_name == "Monotherapy" | time > 2,
                      c_zido, c_zido + c_lam),
        community_medical = c_cmed,
        direct_medical = c_dmed
    )
  )
}, times = c(2, Inf))

4.1 Constructing the economic 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)

A economic model (of class CohortDtstm) can 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 (trans_model) for simulating transition probabilities with sim_stateprobs(), a utility model for simulating quality-adjusted life-years 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)

4.2 Simulating outcomes

econmod$sim_stateprobs(n_cycles = 20)
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.7473092
## 3:      1           1          1      1        1 2 0.5584711
## 4:      1           1          1      1        1 3 0.4173506
## 5:      1           1          1      1        1 4 0.3118899
## 6:      1           1          1      1        1 5 0.2330782

autoplot(econmod$stateprobs_, labels = labs,
         ci = TRUE, ci_style = "ribbon")

econmod$sim_qalys(dr = 0, 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 2.948678 2.948678
## 2:      1           1          1      1        2  0 1.680248 1.680248
## 3:      1           1          1      1        3  0 3.766338 3.766338
## 4:      1           2          1      1        1  0 3.978809 3.978809
## 5:      1           2          1      1        2  0 1.722319 1.722319
## 6:      1           2          1      1        3  0 3.738614 3.738614

econmod$sim_costs(dr = 0.06, 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.06     drug  5407.661
## 2:      1           1          1      1        2 0.06     drug  2866.793
## 3:      1           1          1      1        3 0.06     drug  5392.669
## 4:      1           2          1      1        1 0.06     drug 10438.419
## 5:      1           2          1      1        2 0.06     drug  3206.960
## 6:      1           2          1      1        3 0.06     drug  5228.013