Simulate expected values as a function of simulated state occupancy probabilities, with simulation of costs and quality-adjusted life-years (QALYs) as particular use cases.

# S3 method for stateprobs sim_ev( object, models = NULL, dr = 0.03, integrate_method = c("trapz", "riemann_left", "riemann_right"), value_name = "value", outcome_name = "outcome", ... ) sim_qalys( object, model, dr = 0.03, integrate_method = c("trapz", "riemann_left", "riemann_right"), lys = TRUE ) sim_costs( object, models, dr = 0.03, integrate_method = c("trapz", "riemann_left", "riemann_right") )

object | A |
---|---|

dr | Discount rate. |

integrate_method | Method used to integrate state values when computing
costs or QALYs. Options are |

value_name | Name of the column containing values of the outcome. Default
is |

outcome_name | Name of the column indicating the outcome corresponding
to each model. Only used if |

... | Currently unused. |

model, models | An object or list of objects of class |

lys | If |

`sim_ev()`

returns a `data.table`

with the following columns:

- sample
A random sample from the PSA.

- strategy_id
The treatment strategy ID.

- patient_id
The patient ID.

- grp_id
The subgroup ID.

- state_id
The health state ID.

- dr
The rate used to discount costs.

- outcome
The outcome corresponding to each model in

`models`

. Only included if`models`

is a list.- value
The expected value.

The names of the `outcome`

and `value`

columns may be changed with the
`value_name`

and `outcome_name`

arguments. `sim_costs()`

and `sim_qalys()`

return similar objects, that are of class `costs`

and `qalys`

, respectively.

Expected values in cohort models (i.e., those implemented with
the `CohortDtstm`

and `Psm`

classes) are mean outcomes for patients comprising
the cohort. The method used to simulate expected values depends on the
`$method`

field in the `StateVals`

object(s) stored in `model(s)`

. If
`$method = "starting"`

, then state values represent a one-time value that
occurs at time 0.

The more common use case is `$method = "wlos"`

, or a "weighted length of stay".
That is, expected values for each health state can be thought of as state values
weighted by the time a patient spends in each state (and discounted by a
discount factor that depends on the discount rate `dr`

). The
precise computation proceeds in four steps. In the first step, the probability
of being in each health state at each discrete time point is simulated
(this is the output contained in the `stateprobs`

object). Second, a
`StateVals`

model is used to predict state values at each time point.
Third an expected value at each time point is computed by multiplying the
state probability, the state value, and the discount factor. Fourth, the
expected values at each time point are summed across all time points.

The summation in the fourth step can be thought of as a discrete approximation
of an integral. In particular, the limits of integration can be partitioned
into time intervals, with each interval containing a start and an end.
The `integrate_method`

argument determines the approach used
for this approximation:

A left Riemann sum (

`integrate_method = "riemann_left"`

) uses expected values at the start of each time interval.A right Riemann sum (

`integrate_method = "riemann_right"`

) uses expected values at the end of each time interval.The trapezoid rule (

`integrate_method = "trapz"`

) averages expected values at the start and end of each time interval. (This will generally be the most accurate and is recommended.)

Mathematical details are provided in the reference within the "References" section below.

The ID variables in the state value models in `models`

must be
consistent with the ID variables contained in `object`

. In particular,
the `models`

should predict state values for each non-absorbing health state
in `object`

; that is, the number of health states modeled with the
`models`

should equal the number of health states in `object`

less the number
of absorbing states.

The absorbing states are saved as an attribute named `absorbing`

to
`stateprobs`

objects. When simulating state probabilities with a
`CohortDtstmTrans`

object, the absorbing state is determined by the
`absorbing`

field in the class; in a `Psm`

(or with
`sim_stateprobs.survival()`

), the absorbing state is always equal to the
final health state.

Incerti and Jansen (2021). See Section 2.1 for mathematical details.

State probabilities can be simulated using the
`$sim_stateprobs()`

methods from either the `CohortDtstmTrans`

(or `CohortDtstm`

) or `Psm`

classes. State probabilities can also be
computed directly from survival curves with the generic method
`sim_stateprobs.survival()`

.

Costs and QALYs are typically computed within the `R6`

model classes
using the `$sim_costs()`

and `$sim_qalys()`

methods. For instance, see the
documentation and examples for the `CohortDtstm`

and `Psm`

classes.
The `sim_qalys()`

and `sim_costs()`

functions are exported to give users
additional flexibility when creating their own modeling pipelines.
`sim_ev()`

may be useful for computing outcomes other than costs or QALYs.

`costs`

and `qalys`

objects can be passed to `summarize_ce()`

to
create a cost-effectiveness object for performing a cost-effectiveness analysis
with `cea()`

. Although note that typically the `$summarize()`

method
belonging to the `CohortDtstm`

or `Psm`

classes would be used instead.

Use the `IndivCtstm`

class to simulate costs and QALYs with an individual
continuous-time state transition model.

# We need (i) a state probability object and (ii) a model for state values ## We should start by setting up our decision problem hesim_dat <- hesim_data(strategies = data.frame(strategy_id = 1:2), patients = data.frame(patient_id = 1:3), states = data.frame(state_id = 1)) input_data <- expand(hesim_dat, by = c("strategies", "patients")) ## (i) Simulate a state probability object tpmat_id <- tpmatrix_id(input_data, n_samples = 2) p_12 <- ifelse(tpmat_id$strategy_id == 1, .15, .1) tpmat <- tpmatrix( C, p_12, 0, 1 ) transmod <- CohortDtstmTrans$new(params = tparams_transprobs(tpmat, tpmat_id)) stprobs <- transmod$sim_stateprobs(n_cycles = 3) ## Construct model for state values outcome_tbl <- stateval_tbl( data.frame( state_id = 1, est = 5000 ), dist = "fixed" ) outmod <- create_StateVals(outcome_tbl, n = 2, hesim_data = hesim_dat) # We can then simulate expected values ## The generic expected values function sim_ev(stprobs, models = outmod)#> sample strategy_id patient_id grp_id state_id dr value #> 1: 1 1 1 1 1 0.03 11429.69 #> 2: 1 1 2 1 1 0.03 11429.69 #> 3: 1 1 3 1 1 0.03 11429.69 #> 4: 1 2 1 1 1 0.03 12346.79 #> 5: 1 2 2 1 1 0.03 12346.79 #> 6: 1 2 3 1 1 0.03 12346.79 #> 7: 2 1 1 1 1 0.03 11429.69 #> 8: 2 1 2 1 1 0.03 11429.69 #> 9: 2 1 3 1 1 0.03 11429.69 #> 10: 2 2 1 1 1 0.03 12346.79 #> 11: 2 2 2 1 1 0.03 12346.79 #> 12: 2 2 3 1 1 0.03 12346.79#> sample strategy_id patient_id grp_id state_id dr outcome value #> 1: 1 1 1 1 1 0.03 Outcome 1 11429.69 #> 2: 1 1 2 1 1 0.03 Outcome 1 11429.69 #> 3: 1 1 3 1 1 0.03 Outcome 1 11429.69 #> 4: 1 2 1 1 1 0.03 Outcome 1 12346.79 #> 5: 1 2 2 1 1 0.03 Outcome 1 12346.79 #> 6: 1 2 3 1 1 0.03 Outcome 1 12346.79 #> 7: 2 1 1 1 1 0.03 Outcome 1 11429.69 #> 8: 2 1 2 1 1 0.03 Outcome 1 11429.69 #> 9: 2 1 3 1 1 0.03 Outcome 1 11429.69 #> 10: 2 2 1 1 1 0.03 Outcome 1 12346.79 #> 11: 2 2 2 1 1 0.03 Outcome 1 12346.79 #> 12: 2 2 3 1 1 0.03 Outcome 1 12346.79## Suppose the outcome were a cost category. Then we might ## prefer the following: sim_costs(stprobs, models = list(drug = outmod))#> sample strategy_id patient_id grp_id state_id dr category costs #> 1: 1 1 1 1 1 0.03 drug 11429.69 #> 2: 1 1 2 1 1 0.03 drug 11429.69 #> 3: 1 1 3 1 1 0.03 drug 11429.69 #> 4: 1 2 1 1 1 0.03 drug 12346.79 #> 5: 1 2 2 1 1 0.03 drug 12346.79 #> 6: 1 2 3 1 1 0.03 drug 12346.79 #> 7: 2 1 1 1 1 0.03 drug 11429.69 #> 8: 2 1 2 1 1 0.03 drug 11429.69 #> 9: 2 1 3 1 1 0.03 drug 11429.69 #> 10: 2 2 1 1 1 0.03 drug 12346.79 #> 11: 2 2 2 1 1 0.03 drug 12346.79 #> 12: 2 2 3 1 1 0.03 drug 12346.79## Length of stay is computed if there is no state value model sim_ev(stprobs)#> sample strategy_id patient_id grp_id state_id value #> 1: 1 1 1 1 1 2.285938 #> 2: 1 1 2 1 1 2.285938 #> 3: 1 1 3 1 1 2.285938 #> 4: 1 2 1 1 1 2.469358 #> 5: 1 2 2 1 1 2.469358 #> 6: 1 2 3 1 1 2.469358 #> 7: 2 1 1 1 1 2.285938 #> 8: 2 1 2 1 1 2.285938 #> 9: 2 1 3 1 1 2.285938 #> 10: 2 2 1 1 1 2.469358 #> 11: 2 2 2 1 1 2.469358 #> 12: 2 2 3 1 1 2.469358