Data
We set up the model for two treatment strategies and one representative patient.
library("hesim")
library("data.table")
strategies <- data.table(strategy_id = 1:2,
strategy_name = c("Monotherapy", "Combination therapy"))
patients <- data.table(patient_id = 1)
hesim_dat <- hesim_data(strategies = strategies,
patients = patients)
print(hesim_dat)
## $strategies
## strategy_id strategy_name
## 1: 1 Monotherapy
## 2: 2 Combination therapy
##
## $patients
## patient_id
## 1: 1
##
## attr(,"class")
## [1] "hesim_data"
The model will be simulated for each treatment strategy and patient combination. Such data containing treatment strategy and patient characteristics can be generated using expand().
data <- expand(hesim_dat, by = c("strategies", "patients"))
head(data)
## strategy_id patient_id strategy_name
## 1: 1 1 Monotherapy
## 2: 2 1 Combination therapy
Parameters
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
)
Random number generation
To perform a 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 matrix. For instance, we sample medical costs for each (non-death) health state from gamma distributions and they a stored in a matrix (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 matrix, with parameters of the matrix 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)
The model
The model is defined using define_model() function. The names of transition probability matrix, utility, and cost objects are specified, as are the names of the health states.
mod_def <- define_model(tparams_def = tparams_def,
rng_def = rng_def,
params = params)
