`vignettes/markov-inhomogeneous-indiv.Rmd`

`markov-inhomogeneous-indiv.Rmd`

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.

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
strategies <- data.table(
strategy_id = 1:2,
strategy_name = c("Standard prosthesis", "New prosthesis")
)
n_strategies <- nrow(strategies)
# Patients
n_patients <- 1000
patients <- data.table(
patient_id = 1:n_patients,
gender = "Female",
age = 60
)
# States
states <- data.table( # Non-death health states
state_id = 1:4,
state_name = c("Primary THR", "Sucessful primary", "Revision THR",
"Successful revision")
)
n_states <- nrow(states)
# "hesim data"
hesim_dat <- hesim_data(strategies = strategies,
patients = patients,
states = states)
print(hesim_dat)
```

```
## $strategies
## strategy_id strategy_name
## <int> <char>
## 1: 1 Standard prosthesis
## 2: 2 New prosthesis
##
## $patients
## patient_id gender age
## <int> <char> <num>
## 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
## <int> <char>
## 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.

```
labs <- get_labels(hesim_dat)
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.

```
tmat <- rbind(c(NA, 1, NA, NA, 2),
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.

`ttrrPTHR <- 2 # Time to recovery rate implies mean time of 1/2 years`

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
omrPTHR_shape1 <- 2
omrPTHR_shape2 <- 98
```

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.

```
prob_to_rate <- function(p, t = 1){
(-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.

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.

`ttrRTHR <- 1 # There is no rate, the time is fixed`

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
omrRTHR_shape1 <- 4
omrRTHR_shape2 <- 96
```

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.

`n_samples <- 500`

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.

```
vec_to_dt <- function(v, n = NULL){
if (length(v) == 1) v <- rep(v, n_samples)
dt <- data.table(v)
colnames(dt) <- "cons"
return(dt)
}
```

Our new helper function can then be used within
`define_rng()`

to sample distributions of the
coefficients.

```
transmod_coef_def <- define_rng({
omrPTHR <- prob_to_rate(beta_rng(shape1 = omrPTHR_shape1,
shape2 = omrPTHR_shape2))
mr <- fixed(mr)
mr_omrPTHR <- omrPTHR + mr
colnames(mr) <- colnames(mr_omrPTHR) <- paste0("time_", mr_times)
rr <- multi_normal_rng(mu = rr_coef, Sigma = rr_vcov)
rrr <- prob_to_rate(beta_rng(shape1 = 4, shape2 = 96))
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)
transmod_coef <- eval_rng(transmod_coef_def)
```

To get a better understanding of the output, lets summarize the coefficients.

`summary(transmod_coef)`

```
## param mean sd 2.5% 97.5%
## <char> <num> <num> <num> <num>
## 1: log_omrPTHR_cons -4.15928533 0.80187587 -6.00448991 -2.89731992
## 2: log_mr_time_0 -5.00564775 0.00000000 -5.00564775 -5.00564775
## 3: log_mr_time_5 -3.94765018 0.00000000 -3.94765018 -3.94765018
## 4: log_mr_time_15 -2.92807363 0.00000000 -2.92807363 -2.92807363
## 5: log_mr_time_25 -1.86562132 0.00000000 -1.86562132 -1.86562132
## 6: log_ttrrPTHR_cons 0.69314718 0.00000000 0.69314718 0.69314718
## 7: log_mr_omrPTHR_time_0 -3.73487676 0.50742315 -4.69207087 -2.78269602
## 8: log_mr_omrPTHR_time_5 -3.28260989 0.33128812 -3.82732859 -2.59732805
## 9: log_mr_omrPTHR_time_15 -2.62136052 0.17900219 -2.88298069 -2.21940665
## 10: log_mr_omrPTHR_time_25 -1.74477637 0.07777364 -1.84980587 -1.56076725
## 11: rr_shape_cons 0.37408141 0.04417159 0.29728783 0.45321402
## 12: rr_scale_cons -5.49129351 0.20292983 -5.88147118 -5.12320765
## 13: rr_scale_age -0.03664679 0.00504719 -0.04621155 -0.02716836
## 14: rr_scale_male 0.76408610 0.11371270 0.55668887 0.98231746
## 15: rr_scale_np1 -1.37312553 0.38143367 -2.10721241 -0.60024772
## 16: log_rrr_cons -3.35450357 0.53554700 -4.44843134 -2.48397327
```

It is also helpful to take a look at a few of the sampled coefficient data tables.

`head(transmod_coef$log_omrPTHR)`

```
## cons
## <num>
## 1: -3.411590
## 2: -4.303414
## 3: -4.168389
## 4: -4.870063
## 5: -4.014536
## 6: -4.502533
```

`head(transmod_coef$rr_scale)`

```
## cons age male np1
## <num> <num> <num> <num>
## 1: -5.485059 -0.04167433 0.8650169 -1.5566642
## 2: -5.328609 -0.04400885 0.8981182 -1.4424299
## 3: -5.608818 -0.03818738 0.8393678 -1.5624932
## 4: -5.567563 -0.03424712 0.6938261 -1.5819033
## 5: -5.466357 -0.03934928 0.9023543 -0.5785861
## 6: -5.602355 -0.03126114 0.6852290 -1.0493641
```

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.

The `params_surv_list()`

object is then created as
follows.

```
transmod_params <- params_surv_list(
# 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.

```
utility_tbl <- stateval_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
## <int> <num> <num>
## 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\).

```
drugcost_tbl <- stateval_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"
)
medcost_tbl <- stateval_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()`

.

```
## strategy_id patient_id strategy_name gender age
## <int> <int> <char> <char> <num>
## 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.

```
transmod_data[, cons := 1]
transmod_data[, male := ifelse(gender == "Male", 1, 0)]
transmod_data[, np1 := ifelse(strategy_name == "New prosthesis", 1, 0)]
```

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
transmod <- create_IndivCtstmTrans(transmod_params,
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
utilitymod <- create_StateVals(utility_tbl, n = transmod_coef_def$n,
hesim_data = hesim_dat)
# Costs
drugcostmod <- create_StateVals(drugcost_tbl, n = transmod_coef_def$n,
method = "starting", hesim_data = hesim_dat)
medcostmod <- create_StateVals(medcost_tbl, n = transmod_coef_def$n,
hesim_data = hesim_dat)
costmods <- list(Drug = drugcostmod,
Medical = medcostmod)
```

The full economic model is constructed by combining the disease, utility, and cost models.

```
econmod <- IndivCtstm$new(trans_model = transmod,
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.

```
## user system elapsed
## 1.666 0.100 1.767
```

`head(econmod$disprog_)`

```
## sample strategy_id patient_id grp_id from to final time_start
## <num> <int> <int> <int> <num> <num> <int> <num>
## 1: 1 1 1 1 1 2 0 0.00000000
## 2: 1 1 1 1 2 5 1 0.06989725
## 3: 1 1 2 1 1 2 0 0.00000000
## 4: 1 1 2 1 2 5 1 0.27152781
## 5: 1 1 3 1 1 2 0 0.00000000
## 6: 1 1 3 1 2 3 0 0.01568356
## time_stop
## <num>
## 1: 0.06989725
## 2: 18.99991466
## 3: 0.27152781
## 4: 28.87905653
## 5: 0.01568356
## 6: 20.97806345
```

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.

`econmod$sim_stateprobs(t = 0:60)`

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.

```
econmod$sim_qalys(dr = .015)
econmod$sim_costs(dr = .06)
```

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.

```
## Discount rate Outcome Standard prosthesis New prosthesis
## <num> <fctr> <char> <char>
## 1: 0.015 QALYs 15.38 (14.24, 16.46) 15.42 (14.23, 16.49)
## 2: 0.060 Costs: Drug 394 (394, 394) 579 (579, 579)
## 3: 0.060 Costs: Medical 120 (50, 236) 33 (6, 83)
## 4: 0.060 Costs: total 514 (444, 630) 612 (585, 662)
```

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_out <- cea_pw(ce_sim, comparator = 1, dr_qalys = 0.015, dr_costs = .06,
k = seq(0, 25000, 500))
icer(cea_pw_out, labels = labs) %>%
format(digits_qalys = 3)
```

```
## Outcome New prosthesis
## <fctr> <char>
## 1: Incremental QALYs 0.041 (-0.514, 0.618)
## 2: Incremental costs 97 (6, 161)
## 3: Incremental NMB 1,976 (-25,751, 30,853)
## 4: ICER 2,347
```