Continuous time state transition models (CTSTMs) are used to simulate trajectories for patients between mutually exclusive health states. Transitions between health states \(r\) and \(s\) for patient \(i\) with treatment \(k\) at time \(t\) are governed by hazard functions, \(\lambda_{rs}(t|x_{ik})\), that can depend on covariates \(x_{ik}\).
Different assumptions can be made about the time scales used to determine the hazards. In a “clock forward” (i.e., Markov) model, time \(t\) refers to time since entering the initial health state. Conversely, in a “clock reset” (i.e., semi-Markov) model, time \(t\) refers to time since entering the current state \(r\), meaning that time resets to 0 each time a patient enters a new state.
While state occupancy probabilities in “clock forward” models can be
estimated analytically using the Aalen-Johansen
estimator, state occupancy probabilities in “clock reset” models can
only be computed in a general fashion using individual patient
simulation (IPS). hesim
can simulate either “clock
forward”, “clock reset” models, or combinations of the two via IPS.
Discounted costs and quality-adjusted life-years (QALYs) are computed using the continuous time present value of a flow of state values—\(q_{hik}(t)\) for utility and \(c_{m, hik}(t)\) for the \(m\)th cost category—that depend on the health state of a patient on a given treatment strategy at a particular point in time. Discounted QALYs and costs given a model starting at time \(0\) with a time horizon of \(T\) are then given by,
\[ \begin{aligned} QALYs_{hik} &= \int_{0}^{T} q_{hik}(t)e^{-rt}dt, \\ Costs_{m, hik} &= \int_{0}^{T} c_{m, hik}(t)e^{-rt}dt, \end{aligned} \]
where \(r\) is the discount rate. Note that unlike in cohort models, costs and QALYs can depend on time since entering a health state in individual-level models; that is, costs and QALYs can also be “clock-reset” or “clock-forward”.
This example will demonstrate the use of IPS to simulate a clock
forward model. To facilitate comparison to a cohort approach, we will
revisit the total hip replacement (THR) example from the time inhomogeneous Markov
cohort modeling vignette. The following R
packages will
be used during the analysis.
library("hesim")
library("data.table")
library("kableExtra")
library("flexsurv")
library("ggplot2")
The 5 health states—(i) primary THR, (ii) successful primary, (iii) revision THR, (iv) successful revision, and (v) death—-are displayed again for convenience.
We set up the model for two treatment strategies. We will follow the Decision Modeling for Health Economic Evaluation textbook and simulate results for a 60 year-old female. 1,000 patients will be simulated to ensure that results are reasonably stable.
# Treatment strategies
<- data.table(
strategies strategy_id = 1:2,
strategy_name = c("Standard prosthesis", "New prosthesis")
)<- nrow(strategies)
n_strategies
# Patients
<- 1000
n_patients <- data.table(
patients patient_id = 1:n_patients,
gender = "Female",
age = 60
)
# States
<- data.table( # Non-death health states
states state_id = 1:4,
state_name = c("Primary THR", "Sucessful primary", "Revision THR",
"Successful revision")
) <- nrow(states)
n_states
# "hesim data"
<- hesim_data(strategies = strategies,
hesim_dat patients = patients,
states = states)
print(hesim_dat)
## $strategies
## strategy_id strategy_name
## 1: 1 Standard prosthesis
## 2: 2 New prosthesis
##
## $patients
## patient_id gender age
## 1: 1 Female 60
## 2: 2 Female 60
## 3: 3 Female 60
## 4: 4 Female 60
## 5: 5 Female 60
## ---
## 996: 996 Female 60
## 997: 997 Female 60
## 998: 998 Female 60
## 999: 999 Female 60
## 1000: 1000 Female 60
##
## $states
## state_id state_name
## 1: 1 Primary THR
## 2: 2 Sucessful primary
## 3: 3 Revision THR
## 4: 4 Successful revision
##
## attr(,"class")
## [1] "hesim_data"
get_labels()
is used to obtain nice labels for plots and
summary tables.
<- get_labels(hesim_dat)
labs print(labs)
## $strategy_id
## Standard prosthesis New prosthesis
## 1 2
##
## $state_id
## Primary THR Sucessful primary Revision THR Successful revision
## 1 2 3 4
## Death
## 5
The possible transitions to and from each state are summarized with a matrix in an individual patient simulation.
<- rbind(c(NA, 1, NA, NA, 2),
tmat c(NA, NA, 3, NA, 4),
c(NA, NA, NA, 5, 6),
c(NA, NA, 7, NA, 8),
c(NA, NA, NA, NA, NA))
colnames(tmat) <- rownames(tmat) <- names(labs$state_id)
%>%
tmat kable() %>%
kable_styling()
Primary THR | Sucessful primary | Revision THR | Successful revision | Death | |
---|---|---|---|---|---|
Primary THR | NA | 1 | NA | NA | 2 |
Sucessful primary | NA | NA | 3 | NA | 4 |
Revision THR | NA | NA | NA | 5 | 6 |
Successful revision | NA | NA | 7 | NA | 8 |
Death | NA | NA | NA | NA | NA |
The model for the first transition determines the time from the original surgery until it is deemed “successful”. In the Markov model, we used model cycles that were 1-year long meaning that this process took 1-year. Since there are no model cycles in a CTSTM, we can be more flexible. Let’s assume that time to recovery (TTR) from surgery takes on average 6 months and follows an exponential distribution.
<- 2 # Time to recovery rate implies mean time of 1/2 years ttrrPTHR
The second transition depends on the operative mortality rate following primary THR (omrPTHR). As in the cohort model we will assume the probability of death (within 1-year) is drawn from a beta distributions.
# 2 out of 100 patients receiving primary THR died
<- 2
omrPTHR_shape1 <- 98 omrPTHR_shape2
Below we will convert this probability to a rate so that is can be simulated using an exponential distribution. We write a simple function to do this.
<- function(p, t = 1){
prob_to_rate -log(1 - p))/t
( }
Transition 3 depends on the revision rate (rr
) for a
prosthesis, which will again be modeled using a proportional hazards
Weibull model. The coefficients and variance-covariance matrix are
displayed below.
print(rr_coef)
## lngamma cons age male np1
## 0.3740968 -5.4909350 -0.0367022 0.7685360 -1.3444740
print(rr_vcov)
## lngamma cons age male np1
## lngamma 0.002251512 -0.00569100 2.800000e-08 0.00000510 0.0002590
## cons -0.005691000 0.04321908 -7.830000e-04 -0.00724700 -0.0006420
## age 0.000000028 -0.00078300 2.715661e-05 0.00003300 -0.0001110
## male 0.000005100 -0.00724700 3.300000e-05 0.01189539 0.0001840
## np1 0.000259000 -0.00064200 -1.110000e-04 0.00018400 0.1463686
The fourth transition is modeled using the death rate following recovery from the THR (i.e., while in the successful primary state). We will again use age and sex specific mortality rates.
Since we are running the simulation for a 60-year old female, the mortality rate at time \(0\) is \(0.0067\). We are using a clock-forward model, so the hazard depends on time since the start of the model. We can therefore model mortality over time with a piecewise exponential distribution with rates that vary over time. For a 60-year female, the rates will change at 5, 15, and 25 years.
<- c(.0067, .0193, .0535, .1548)
mr <- c(0, 5, 15, 25) mr_times
Transition 5 is similar to transition 1 in that it is a function of time to recovery from a THR. While we could again model this process with an exponential distribution, we will instead assume this time is “fixed” for illustration purposes. Specifically, let’s assume it takes one year for every patient.
<- 1 # There is no rate, the time is fixed ttrRTHR
Following the prior example, the sixth transition depends on the operative mortality rate following revision THR (omrRTHR) and the overall mortality rate (mr). The omrRTHR is assumed to follow the same distribution as the omrPTHR.
These probabilities will be converted to rates below and then added to overall mortality rate (mr) at times 0, 5, 15, and 25 years. Transition 6 will consequently be modeled using a piecewise exponential distribution.
The re-revision rate (rrr) is used to model the seventh transition. Like the omrPTHR, the probability is modeled with a beta distribution and will be converted to a rate below.
# 4 out of 100 patients with a successful revision needed another procedure r
<- 4
omrRTHR_shape1 <- 96 omrRTHR_shape2
The final transition is time to death following a successful revision. We model it with the same model used for transition 4; that is, we assume mortality rates following a successful revision THR are the same as mortality rates following a successful primary THR.
The transition-specific estimates can be combined to create a multi-state model, which is comprised of survival models for each transition. Since we will perform a probabilistic sensitivity analysis (PSA), we will sample the parameters from the distributions described above. 500 iterations will be used.
<- 500 n_samples
Survival models are regression models, so the parameters are
regression coefficients. The coefficients of each model must be a matrix
or data frame (rows for parameter samples and columns for covariates).
Since we often use intercept only models, let’s write a simple function
to convert a vector to a data.table
with a single column
for the intercept.
<- function(v, n = NULL){
vec_to_dt if (length(v) == 1) v <- rep(v, n_samples)
<- data.table(v)
dt colnames(dt) <- "cons"
return(dt)
}
Our new helper function can then be used within
define_rng()
to sample distributions of the
coefficients.
<- define_rng({
transmod_coef_def <- prob_to_rate(beta_rng(shape1 = omrPTHR_shape1,
omrPTHR shape2 = omrPTHR_shape2))
<- fixed(mr)
mr <- omrPTHR + mr
mr_omrPTHR colnames(mr) <- colnames(mr_omrPTHR) <- paste0("time_", mr_times)
<- multi_normal_rng(mu = rr_coef, Sigma = rr_vcov)
rr <- prob_to_rate(beta_rng(shape1 = 4, shape2 = 96))
rrr
list(
log_omrPTHR = vec_to_dt(log(omrPTHR)),
log_mr = log(mr),
log_ttrrPTHR = vec_to_dt(log(ttrrPTHR)),
log_mr_omrPTHR = log(mr_omrPTHR),
rr_shape = vec_to_dt(rr$lngamma),
rr_scale = rr[, -1,],
log_rrr = vec_to_dt(log(rrr))
)n = n_samples)
}, <- eval_rng(transmod_coef_def) transmod_coef
To get a better understanding of the output, lets summarize the coefficients.
summary(transmod_coef)
## param mean sd 2.5% 97.5%
## 1: log_omrPTHR_cons -4.10477726 0.779672371 -5.76310845 -2.88705363
## 2: log_mr_time_0 -5.00564775 0.000000000 -5.00564775 -5.00564775
## 3: log_mr_time_5 -3.94765018 0.000000000 -3.94765018 -3.94765018
## 4: log_mr_time_15 -2.92807363 0.000000000 -2.92807363 -2.92807363
## 5: log_mr_time_25 -1.86562132 0.000000000 -1.86562132 -1.86562132
## 6: log_ttrrPTHR_cons 0.69314718 0.000000000 0.69314718 0.69314718
## 7: log_mr_omrPTHR_time_0 -3.70117065 0.510505368 -4.62115646 -2.77354455
## 8: log_mr_omrPTHR_time_5 -3.26042174 0.337918149 -3.79684654 -2.58972850
## 9: log_mr_omrPTHR_time_15 -2.60918341 0.185042100 -2.87101444 -2.21420249
## 10: log_mr_omrPTHR_time_25 -1.73939478 0.081331045 -1.84553093 -1.55807582
## 11: rr_shape_cons 0.37255535 0.046646488 0.28768668 0.46684078
## 12: rr_scale_cons -5.48054126 0.206073052 -5.86057729 -5.10467604
## 13: rr_scale_age -0.03689076 0.005386752 -0.04699827 -0.02647981
## 14: rr_scale_male 0.76739930 0.107506239 0.56978939 0.97699170
## 15: rr_scale_np1 -1.33841288 0.383069512 -2.12174005 -0.62716192
## 16: log_rrr_cons -3.30694248 0.540256920 -4.48356733 -2.42015321
It is also helpful to take a look at a few of the sampled coefficient data tables.
head(transmod_coef$log_omrPTHR)
## cons
## 1: -3.225382
## 2: -4.166867
## 3: -4.217561
## 4: -4.518022
## 5: -3.915191
## 6: -4.319199
head(transmod_coef$rr_scale)
## cons age male np1
## 1: -5.370253 -0.04343824 0.5988311 -2.0592318
## 2: -5.710184 -0.02507501 0.6306379 -0.3491912
## 3: -5.789766 -0.03138916 0.5626969 -1.6720223
## 4: -5.447601 -0.03706895 0.6568232 -1.5348107
## 5: -4.978150 -0.05180358 0.8589875 -1.7743936
## 6: -5.627855 -0.02914147 0.7767687 -1.2065555
In addition to the coefficients, a complete parameterization of a
transition model requires specification of the survival distributions
and potentially auxiliary information (e.g., the times at which rates
change in a piecewise exponential model). The parameters are stored in a
params_surv_list()
object, which is a list of
params_surv()
objects.
To create the params_surv()
objects, it’s helpful to
first write another short helper function that can convert a single data
table storing the rates for each time in a piecewise exponential to a
list of regression coefficients data tables.
<- function(x) {
as_dt_list lapply(as.list(x), vec_to_dt)
}
The params_surv_list()
object is then created as
follows.
<- params_surv_list(
transmod_params # 1. Primary THR:Successful primary (1:2)
params_surv(coefs = list(rate = transmod_coef$log_ttrrPTHR),
dist = "exp"),
# 2. Primary THR:Death (1:5)
params_surv(coefs = list(rate = transmod_coef$log_omrPTHR),
dist = "exp"),
# 3. Successful primary:Revision THR (2:3)
params_surv(coefs = list(shape = transmod_coef$rr_shape,
scale = transmod_coef$rr_scale),
dist = "weibullPH"),
# 4. Successful primary:Death (2:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp"),
# 5. Revision THR:Successful revision (3:4)
params_surv(coefs = list(est = vec_to_dt(ttrRTHR)),
dist = "fixed"),
# 6. Revision THR:Death (3:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr_omrPTHR),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp"),
# 7. Successful revision:Revision THR (4:3)
params_surv(coefs = list(rate = transmod_coef$log_rrr),
dist = "exp"),
# 8. Successful revision:Death (4:5)
params_surv(coefs = as_dt_list(transmod_coef$log_mr),
aux = list(time = c(0, 5, 15, 25)),
dist = "pwexp")
)
Costs and utilities are unchanged from the cohort model. The mean (standard error) of utility are estimated be \(0.85\) (\(0.03\)), \(0.30\) (\(0.03\)), and \(0.75\) (\(0.04\)) in the successful primary, revision, and successful revision health states, respectively.
<- stateval_tbl(
utility_tbl data.table(state_id = states$state_id,
mean = c(0, .85, .3, .75),
se = c(0, .03, .03, .04)),
dist = "beta"
)head(utility_tbl)
## state_id mean se
## 1: 1 0.00 0.00
## 2: 2 0.85 0.03
## 3: 3 0.30 0.03
## 4: 4 0.75 0.04
The standard prosthesis costs \(£394\) while the new prosthesis costs \(£579\). Both are assumed to be known with certainty.
Since the model assumes that there are no ongoing medical costs, the only remaining cost is the cost of the revision procedure, which was estimated to have a mean of \(£5,294\) and standard error of \(£1,487\).
<- stateval_tbl(
drugcost_tbl data.table(strategy_id = rep(strategies$strategy_id, each = n_states),
state_id = rep(states$state_id, times = n_strategies),
est = c(394, 0, 0, 0,
579, 0, 0, 0)),
dist = "fixed"
)
<- stateval_tbl(
medcost_tbl data.table(state_id = states$state_id,
mean = c(0, 0, 5294, 0),
se = c(0, 0, 1487, 0)),
dist = "gamma",
)
The economic model consists of a model for disease progression and models for assigning utility and cost values to health states.
The transition model is a function of the parameters of the
multi-state model and input data. The input data in this case is a
dataset consisting of one row for each treatment strategy and patient.
It can be generated using expand.hesim_data()
.
<- expand(hesim_dat, by = c("strategies", "patients"))
transmod_data head(transmod_data)
## strategy_id patient_id strategy_name gender age
## 1: 1 1 Standard prosthesis Female 60
## 2: 1 2 Standard prosthesis Female 60
## 3: 1 3 Standard prosthesis Female 60
## 4: 1 4 Standard prosthesis Female 60
## 5: 1 5 Standard prosthesis Female 60
## 6: 1 6 Standard prosthesis Female 60
Since we are using a Weibull regression model to simulate time until a revision hip replacement (Successful primary to Revision THR), we need to create a constant term and covariates for male sex and whether the new prosthesis was used.
:= 1]
transmod_data[, cons := ifelse(gender == "Male", 1, 0)]
transmod_data[, male := ifelse(strategy_name == "New prosthesis", 1, 0)] transmod_data[, np1
The full transition model is created using the transition parameters,
the input data, the transition matrix, and the starting age of the
patients. Since we are using a clock-forward model we use the
clock = "forward"
option.
# Transition model
<- create_IndivCtstmTrans(transmod_params,
transmod input_data = transmod_data,
trans_mat = tmat,
clock = "forward",
start_age = patients$age)
Models based on predicted means (see tparams_mean()
) can
be created directly from the utility and cost tables using since they do
not include covariates and therefore do not require input data.
# Utility
<- create_StateVals(utility_tbl, n = transmod_coef_def$n,
utilitymod hesim_data = hesim_dat)
# Costs
<- create_StateVals(drugcost_tbl, n = transmod_coef_def$n,
drugcostmod method = "starting", hesim_data = hesim_dat)
<- create_StateVals(medcost_tbl, n = transmod_coef_def$n,
medcostmod hesim_data = hesim_dat)
<- list(Drug = drugcostmod,
costmods Medical = medcostmod)
The full economic model is constructed by combining the disease, utility, and cost models.
<- IndivCtstm$new(trans_model = transmod,
econmod utility_model = utilitymod,
cost_models = costmods)
Disease progression is simulated using the
$sim_disease()
method. Unique trajectories are simulated
for each patient, treatment strategy, and PSA sample. Patients
transition from
an old health state that was entered at
time time_start
to
a new health state at time
time_stop
. We will simulate the model for 60 years, or
equivalently, until they reach a maximum age of 120. As shown by the run
time (in seconds), hesim
is quite fast.
<- proc.time()
ptm $sim_disease(max_t = 60, max_age = 120)
econmodproc.time() - ptm
## user system elapsed
## 0.857 0.039 0.897
head(econmod$disprog_)
## sample strategy_id patient_id grp_id from to final time_start time_stop
## 1: 1 1 1 1 1 2 0 0.00000000 0.05815911
## 2: 1 1 1 1 2 5 1 0.05815911 3.75205054
## 3: 1 1 2 1 1 2 0 0.00000000 0.31728230
## 4: 1 1 2 1 2 5 1 0.31728230 32.39675219
## 5: 1 1 3 1 1 2 0 0.00000000 1.64967667
## 6: 1 1 3 1 2 5 1 1.64967667 27.15122529
The simulated patient trajectories can be summarized by computing the probability of being in each health state over time (by averaging across the simulated patients). Here, we will compute state occupancy probabilities at yearly intervals. They are generally quite similar to those in the cohort model.
$sim_stateprobs(t = 0:60) econmod
Following the cohort model, we simulate costs and QALYs with 6% and
1.5% discount rates, respectively. While we are currently assuming that
costs and QALYs are constant over time, we could have let them depend on
time since entering a health state by setting
time_reset = TRUE
in the StateVals
objects.
$sim_qalys(dr = .015)
econmod$sim_costs(dr = .06) econmod
Mean costs and QALYs for each PSA sample are computed with the
$summarize()
method and summary results are produced with
summary.ce()
. Costs are very similar to the cohort model.
Conversely, QALYs are slightly higher in the IPS since we assumed that
patients remained in the primary THR (with utility = 0) state for longer
than in the cohort model. This difference not only suggests that results
may be sensitive to this assumption, but also shows that the flexibility
of individual-level models may make them more realistic than cohort
models.
<- econmod$summarize()
ce_sim summary(ce_sim, labels = labs) %>%
format()
## Discount rate Outcome Standard prosthesis New prosthesis
## 1: 0.015 QALYs 15.38 (14.16, 16.38) 15.42 (14.06, 16.49)
## 2: 0.060 Costs: Drug 394 (394, 394) 579 (579, 579)
## 3: 0.060 Costs: Medical 122 (51, 227) 34 (7, 86)
## 4: 0.060 Costs: total 516 (445, 621) 613 (586, 665)
Since we performed a PSA, decision uncertainty could be formally quantified. Since we cover that in detail in the cost-effectiveness analysis vignette, we will simply compute the incremental cost-effectiveness ratio (ICER) here.
<- cea_pw(ce_sim, comparator = 1, dr_qalys = 0.015, dr_costs = .06,
cea_pw_out k = seq(0, 25000, 500))
icer(cea_pw_out, labels = labs) %>%
format(digits_qalys = 3)
## Outcome New prosthesis
## 1: Incremental QALYs 0.038 (-0.560, 0.608)
## 2: Incremental costs 97 (4, 159)
## 3: Incremental NMB 1,780 (-28,042, 30,282)
## 4: ICER 2,587