In this tutorial we use a continuous time state transition model (CTSTM) and relax many of the assumptions made in cohort discrete time state transition models (DTSTMs). First, since the model is in continuous time we do not require model cycles. Second, we estimate the parameters of the health state transitions using a multi-state model so that the simulation model is completely integrated with an underlying statistical model. Third, we use individual patient simulation (IPS) to simulate a semi-Markov model, meaning that (unlike in a Markov model) transitions cannot depend on prior history.

To illustrate, we simplify the sick-sicker model so that it only contains three health states and modify the states—*Stable*, *Progression*, and *Dead*—to mimic an oncology application where patients transition from stable disease to progression to death. There are three transitions: (1) *Stable* to *Progression*, (2) *Stable* to *Dead*, and (3) *Progression* to *Dead*.

The following packages and settings will be used for the analysis. Note that while individual-level simulations can be computationally intensive, they run very quickly in `hesim`

because they are implemented fully in `C++`

under the hood. You can learn more by looking at the `hesim`

multi-state modeling vignette.

Multi-state models are generalizations of survival models to more than 2 states that estimate hazard functions for each possible transition. Flexible survival models can be used for each transition ensuring that the rates at which patients transition between states are consistent with the available data.

Different assumptions can be made about the time scales used to determine the hazards. In a *clock forward* (i.e., *Markov*) model, transition hazards depends on time since entering the initial health state. Conversely, in a *clock reset* (i.e., *semi-Markov*) model, the hazards are a function of time since entering the current state (i.e., time resets to 0 each time a patient enters a new state).

Importantly, state occupancy probabilities in clock-reset models can only be simulated in a general fashion using IPS (although tunnel states can be used in a cohort model to approximate a semi-Markov process). In an IPS multiple patients are simulated and the expected value of outcomes is computed by averaging across patients. It is important to simulate enough patients so that the expected values are stable and not subject to significant Monte Carlo error.

We suggest the tutorial by Putter et al. for additional details on multi-state modeling.

The transitions of a multi-state model in `hesim`

must be characterized by a matrix where each element denotes a transition from a row to a column. If a transition is not possible it is marked with `NA`

; otherwise, an integer denotes the transition number.

tmat <- rbind( c(NA, 1, 2), c(NA, NA, 3), c(NA, NA, NA) ) colnames(tmat) <- rownames(tmat) <- c("Stable", "Progression", "Dead") print(tmat)

```
## Stable Progression Dead
## Stable NA 1 2
## Progression NA NA 3
## Dead NA NA NA
```

As in the cohort model, we must specify the target population and treatment strategies of interest. Unlike the cohort model we simulate 1,000 patients which should be enough to produce reasonably stable results. We also explicitly define the health states, which we will use to model utility and costs.

n_patients <- 1000 patients <- data.table( patient_id = 1:n_patients, age = rnorm(n_patients, mean = 45, sd = 7), female = rbinom(n_patients, size = 1, prob = .51) ) states <- data.table( state_id = c(1, 2), state_name = c("Stable", "Progression") # Non-death health states ) n_states <- nrow(states) strategies <- data.frame( strategy_id = 1:2, strategy_name = c("SOC", "New") ) n_strategies <- nrow(strategies) hesim_dat <- hesim_data( strategies = strategies, patients = patients, states = states ) print(hesim_dat)

```
## $strategies
## strategy_id strategy_name
## 1 1 SOC
## 2 2 New
##
## $patients
## patient_id age female
## 1: 1 42.71774 0
## 2: 2 48.86723 0
## 3: 3 40.27539 1
## 4: 4 46.50052 1
## 5: 5 47.17538 1
## ---
## 996: 996 43.22279 0
## 997: 997 54.22166 0
## 998: 998 37.27172 0
## 999: 999 45.20616 0
## 1000: 1000 39.15981 1
##
## $states
## state_id state_name
## 1: 1 Stable
## 2: 2 Progression
##
## attr(,"class")
## [1] "hesim_data"
```

In the cohort examples, we used parameter estimates derived from external sources. An alternative approach is to estimate the parameters directly with `R`

. We will take that approach here for the multi-state model and estimate the transition specific hazards using `flexsurv`

. Conversely, we will continue to use mean values for utility and costs.

Multi-state data consists of one row for each possible transition from a given health state where only one transition is observed and all others are censored. We created the function `sim_mstate_data()`

which leverages `hesim`

to simulate a hypothetical cohort of 2000 patients for a stable-progression-death model.

data <- rcea::sim_mstate3_data(n = 2000) data[patient_id %in% c(1, 2)]

```
## from to female age patient_id final time_start time_stop status
## 1: 1 2 1 54.37365 1 0 0.000000 7.042086 1
## 2: 1 3 1 54.37365 1 0 0.000000 7.042086 0
## 3: 2 3 1 54.37365 1 1 7.042086 15.000000 0
## 4: 1 2 1 58.93437 2 0 0.000000 7.128476 1
## 5: 1 3 1 58.93437 2 0 0.000000 7.128476 0
## 6: 2 3 1 58.93437 2 1 7.128476 11.038201 1
## transition_id strategy_id strategy_name time
## 1: 1 2 New 7.042086
## 2: 2 2 New 7.042086
## 3: 3 2 New 7.957914
## 4: 1 1 SOC 7.128476
## 5: 2 1 SOC 7.128476
## 6: 3 1 SOC 3.909726
```

Multi-state models can be fit by estimating a joint survival model with interaction terms for different transition (see the `hesim`

introductory vignette for more details) or by fitting separate survival models for each transition. We will take the latter approach and and fit 3 parametric Weibull models. (Note that in an applied example you should compare the performance of multiple distributions. Type `?params_surv`

to view the survival distributions supported by `hesim`

.)

n_trans <- max(tmat, na.rm = TRUE) # Number of transitions wei_fits <- vector(length = n_trans, mode = "list") for (i in 1:length(wei_fits)){ wei_fits[[i]] <- flexsurvreg( Surv(time, status) ~ strategy_name + female, data = data, subset = (transition_id == i) , dist = "weibullPH") } wei_fits <- flexsurvreg_list(wei_fits)

In the cohort model, we assigned utility and cost values to health states using `hesim::define_tparams()`

within the `hesim::define_model()`

framework. An alternative (and more direct) approach is to use a `hesim::stateval_tbl`

. We will continue to assume that the utility associated with each state follows a beta distribution.

utility_tbl <- stateval_tbl( data.table(state_id = states$state_id, mean = c(.8, .6), se = c(0.02, .05) ), dist = "beta", hesim_data = hesim_dat) print(utility_tbl)

```
## state_id mean se
## 1: 1 0.8 0.02
## 2: 2 0.6 0.05
```

Similarly, we continue to assume that medical costs vary by health state and follow a gamma distribution.

medcost_tbl <- stateval_tbl( data.table(state_id = states$state_id, mean = c(2000, 9500), se = c(2000, 9500) ), dist = "gamma", hesim_data = hesim_dat) print(medcost_tbl)

```
## state_id mean se
## 1: 1 2000 2000
## 2: 2 9500 9500
```

Treatment costs remain fixed as in the cohort model, but we now allow costs to vary over time in the progression state. Specifically, we assume that costs for the new treatment are $12,000 for the first 3 months in the progression state and then $10,000 thereafter. This approximate a common scenario in oncology where the new treatment is given in combination with chemotherapy but the chemotherapy is only given for a fixed number of cycles. This would not be possible in a cohort model without creating a tunnel state approximating the 3 months of chemotherapy in the progression state.

n_times <- 2 drugcost_tbl <- stateval_tbl( data.table( strategy_id = rep(strategies$strategy_id, each = n_states * n_times), state_id = rep(rep(states$state_id, each = n_strategies), n_times), time_start = rep(c(0, 3/12), n_states * n_strategies), est = c(rep(2000, 4), # Costs are always the same with SOC 12000, 12000, 12000, 10000 # Costs with new drop after 3 months in progression state ) ), dist = "fixed", hesim_data = hesim_dat) print(drugcost_tbl)

```
## strategy_id state_id time_id time_start time_stop est
## 1: 1 1 1 0.00 0.25 2000
## 2: 1 1 2 0.25 Inf 2000
## 3: 1 2 1 0.00 0.25 2000
## 4: 1 2 2 0.25 Inf 2000
## 5: 2 1 1 0.00 0.25 12000
## 6: 2 1 2 0.25 Inf 12000
## 7: 2 2 1 0.00 0.25 12000
## 8: 2 2 2 0.25 Inf 10000
```

Economic models in `hesim`

can be constructed in one of two ways: first, by defining the model using `hesim::define_model()`

, and secondly, with explicit statistical models. We used the first approach for the cohort model but will use the latter here.

500 iterations are used for the probabilistic sensitivity analysis. (In an applied application you should check the sensitivity of the results to this number and may want to increase it.)

n_samples <- 500

Health state transition models in `hesim`

are a function of input data (i.e., covariates) and a fitted multi-state model (or a parameter object such as `hesim::params_surv()`

). Since we fit separate models for each transition, the input data consists of one observation for each treatment strategy and patient combination (joint models consist of one observation for each treatment strategy, patient, and transition combination). The data can be created easily by using the `hesim::expand()`

function to expand the `hesim_data`

object created above.

```
## strategy_id patient_id strategy_name age female
## 1: 1 1 SOC 42.71774 0
## 2: 1 2 SOC 48.86723 0
## 3: 1 3 SOC 40.27539 1
## 4: 1 4 SOC 46.50052 1
## 5: 1 5 SOC 47.17538 1
## 6: 1 6 SOC 53.21776 0
```

In addition to the input data and model, we must also specify the viable transitions, the clock (“reset” or “forward”), and the starting age of each simulated patient (so that we can control the maximum age they can live to and ensure that they do not live to unrealistically high ages).

transmod <- create_IndivCtstmTrans(wei_fits, transmod_data, trans_mat = tmat, n = n_samples, clock = "reset", start_age = patients$age)

The utility and cost models are based on predicted means (see `hesim::tparams_mean()`

), so they do not include covariates and therefore do not require input data. The cost and utility models can therefore be constructed directly from the utility and cost tables.

Since drug costs depend on time in each state, we must use the `time_reset = TRUE`

option so that time resets when a patient enters a new state (e.g., when entering the progression state); otherwise the time intervals specified would depend on time since the start of the model. The distinction is analogous to the distinction between clock reset and clock forward multi-state models of disease progression.

# Utility utilitymod <- create_StateVals(utility_tbl, n = n_samples) # Costs drugcostmod <- create_StateVals(drugcost_tbl, n = n_samples, time_reset = TRUE) medcostmod <- create_StateVals(medcost_tbl, n = n_samples) costmods <- list(Drug = drugcostmod, Medical = medcostmod)

In the cohort model we initialized the economic model from a `model_def`

object (the output of `define_model()`

) and the separate the transition, utility, and cost models were constructed under the hood. We do that directly here.

econmod <- IndivCtstm$new(trans_model = transmod, utility_model = utilitymod, cost_models = costmods)

Disease progression is simulated using the `$sim_disease()`

method, which simulates unique trajectories through the multi-state model 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`

.

econmod$sim_disease(max_age = 100) 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.000000 2.791395
## 2: 1 1 1 1 2 3 1 2.791395 10.896754
## 3: 1 1 2 1 1 2 0 0.000000 3.445021
## 4: 1 1 2 1 2 3 1 3.445021 10.762468
## 5: 1 1 3 1 1 2 0 0.000000 7.344869
## 6: 1 1 3 1 2 3 1 7.344869 13.312544
```

The IPS can be summarized by computing state occupancy probabilities at different time points. We do that using `$sim_stateprobs()`

.

econmod$sim_stateprobs(t = seq(0, 30 , 1/12))

Like the cohort model, we can compute quality-adjusted life-years (QALYs) after simulating disease progression using the `$sim_qalys()`

method. Let’s compute QALYs for different discount rates, which might be useful for summarizing the impact of the new treatment on QALYs.

```
## sample strategy_id grp_id state_id dr qalys lys
## 1: 1 1 1 1 0 3.186059 4.136118
## 2: 1 1 1 2 0 2.947315 5.912673
## 3: 1 2 1 1 0 4.056467 5.266075
## 4: 1 2 1 2 0 3.158765 6.336868
## 5: 2 1 1 1 0 3.114647 4.001087
## 6: 2 1 1 2 0 2.821441 5.925158
```

Costs for each category are simulated in a similar fashion.

```
## sample strategy_id grp_id state_id dr category costs
## 1: 1 1 1 1 0.03 Drug 7671.951
## 2: 1 1 1 2 0.03 Drug 9410.919
## 3: 1 2 1 1 0.03 Drug 57454.111
## 4: 1 2 1 2 0.03 Drug 48586.545
## 5: 2 1 1 1 0.03 Drug 7430.945
## 6: 2 1 1 2 0.03 Drug 9452.131
```

As in the cohort model, cost-effectiveness analyses can be performed seamlessly from the simulation output.

ce_sim <- econmod$summarize() cea_pw_out <- cea_pw(ce_sim, comparator = 1, dr_qalys = .03, dr_costs = .03, k = seq(0, 25000, 500)) icer_tbl(cea_pw_out, colnames = strategies$strategy_name)

```
## SOC New
## Incremental QALYs "-" "0.88 (0.62, 1.12)"
## Incremental costs "-" "93,344 (86,585, 104,307)"
## Incremental NMB "-" "-49,564 (-60,742, -39,568)"
## ICER "-" "106,607"
## Conclusion "-" "Not cost-effective"
```

Let’s save the cost-effectiveness object so that we can use it in the CEA tutorial.

saveRDS(ce_sim, "ce_sim.rds")