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")
)A stateprobs object.
Discount rate.
Method used to integrate state values when computing
costs or QALYs. Options are trapz (the default) for the trapezoid rule,
riemann_left left for a left Riemann sum, and
riemann_right right for a right Riemann sum.
Name of the column containing values of the outcome. Default
is "value".
Name of the column indicating the outcome corresponding
to each model. Only used if models is a list. Default is "outcome".
Currently unused.
An object or list of objects of class StateVals used to
model state values. When using sim_qalys(), this should be
a single model for utility. With sim_costs(), a list of models should be
used with one model for each cost category. Finally, with sim_ev(),
this may either be a single model or a list of models. May also be NULL,
in which case length of stay is computed based on the state probabilities
contained in object.
If TRUE, then life-years are simulated in addition to
QALYs.
sim_ev() returns a data.table with the following columns:
A random sample from the PSA.
The treatment strategy ID.
The patient ID.
The subgroup ID.
The health state ID.
The rate used to discount costs.
The outcome corresponding to each model in models.
Only included if models is a list.
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
#> <num> <int> <int> <int> <int> <num> <num>
#> 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
## We can also pass a list of models
sim_ev(stprobs, models = list(`Outcome 1` = outmod))
#> sample strategy_id patient_id grp_id state_id dr outcome value
#> <num> <int> <int> <int> <int> <num> <char> <num>
#> 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
#> <num> <int> <int> <int> <int> <num> <char> <num>
#> 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
#> <num> <int> <int> <int> <num> <num>
#> 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