CohortDtstm | R Documentation |
Simulate outcomes from a cohort discrete time state transition model.
An R6::R6Class object.
trans_model
The model for health state transitions. Must be an object
of class CohortDtstmTrans
.
utility_model
The model for health state utility. Must be an object of
class StateVals
.
cost_models
The models used to predict costs by health state.
Must be a list of objects of class StateVals
, where each element of the
list represents a different cost category.
stateprobs_
An object of class stateprobs
simulated using $sim_stateprobs()
.
qalys_
An object of class qalys
simulated using $sim_qalys()
.
costs_
An object of class costs
simulated using $sim_costs()
.
new()
Create a new CohortDtstm
object.
CohortDtstm$new(trans_model = NULL, utility_model = NULL, cost_models = NULL)
trans_model
The trans_model
field.
utility_model
The utility_model
field.
cost_models
The cost_models
field.
A new CohortDtstm
object.
sim_stateprobs()
Simulate health state probabilities using CohortDtstmTrans$sim_stateprobs()
.
CohortDtstm$sim_stateprobs(n_cycles)
n_cycles
The number of model cycles to simulate the model for.
An instance of self
with simulated output of class stateprobs
stored in stateprobs_
.
sim_qalys()
Simulate quality-adjusted life-years (QALYs) as a function of stateprobs_
and
utility_model
. See sim_qalys()
for details.
CohortDtstm$sim_qalys( dr = 0.03, integrate_method = c("trapz", "riemann_left", "riemann_right"), lys = TRUE )
dr
Discount rate.
integrate_method
Method used to integrate state values when computing
costs or QALYs. Options are trapz
for the trapezoid rule,
riemann_left
for a left Riemann sum, and
riemann_right
for a right Riemann sum.
lys
If TRUE
, then life-years are simulated in addition to QALYs.
An instance of self
with simulated output of class qalys stored
in qalys_
.
sim_costs()
Simulate costs as a function of stateprobs_
and cost_models
.
See sim_costs()
for details.
CohortDtstm$sim_costs( dr = 0.03, integrate_method = c("trapz", "riemann_left", "riemann_right") )
dr
Discount rate.
integrate_method
Method used to integrate state values when computing
costs or QALYs. Options are trapz
for the trapezoid rule,
riemann_left
for a left Riemann sum, and
riemann_right
for a right Riemann sum.
An instance of self
with simulated output of class costs stored
in costs_
.
summarize()
Summarize costs and QALYs so that cost-effectiveness analysis can be performed.
See summarize_ce()
.
CohortDtstm$summarize(by_grp = FALSE)
by_grp
If TRUE
, then costs and QALYs are computed by subgroup. If
FALSE
, then costs and QALYs are aggregated across all patients (and subgroups).
clone()
The objects of this class are cloneable with this method.
CohortDtstm$clone(deep = FALSE)
deep
Whether to make a deep clone.
Incerti and Jansen (2021). See Section 2.1 for a description of a cohort DTSTM and details on simulating costs and QALYs from state probabilities. An example in oncology is provided in Section 4.3.
CohortDtstm
objects can be created from model objects as
documented in create_CohortDtstm()
. The CohortDtstmTrans
documentation
describes the class for the transition model and the StateVals
documentation
describes the class for the cost and utility models. A CohortDtstmTrans
object is typically created using create_CohortDtstmTrans()
.
There are currently three relevant vignettes. vignette("markov-cohort")
details a relatively simple Markov model and
vignette("markov-inhomogeneous-cohort")
describes a more complex time
inhomogeneous model in which transition probabilities vary in every model
cycle. The vignette("mlogit")
shows how a transition model can be parameterized
using a multinomial logistic regression model when transition data is collected
at evenly spaced intervals.
library("data.table")
library("ggplot2")
theme_set(theme_bw())
set.seed(102)
# NOTE: This example replicates the "Simple Markov cohort model"
# vignette using a different approach. Here, we explicitly construct
# the transition probabilities "by hand". In the vignette, the transition
# probabilities are defined using expressions (i.e., by using
# `define_model()`). The `define_model()` approach does (more or less) what
# is done here under the hood.
# (0) Model setup
hesim_dat <- hesim_data(
strategies = data.table(
strategy_id = 1:2,
strategy_name = c("Monotherapy", "Combination therapy")
),
patients <- data.table(patient_id = 1),
states = data.table(
state_id = 1:3,
state_name = c("State A", "State B", "State C")
)
)
n_states <- nrow(hesim_dat$states) + 1
labs <- get_labels(hesim_dat)
# (1) Parameters
n_samples <- 10 # Number of samples for PSA
## Transition matrix
### Input data (one transition matrix for each parameter sample,
### treatment strategy, patient, and time interval)
p_id <- tpmatrix_id(expand(hesim_dat, times = c(0, 2)), n_samples)
N <- nrow(p_id)
### Transition matrices (one for each row in p_id)
p <- array(NA, dim = c(n_states, n_states, nrow(p_id)))
#### Baseline risk
trans_mono <- rbind(
c(1251, 350, 116, 17),
c(0, 731, 512, 15),
c(0, 0, 1312, 437),
c(0, 0, 0, 469)
)
mono_ind <- which(p_id$strategy_id == 1 | p_id$time_id == 2)
p[,, mono_ind] <- rdirichlet_mat(n = 2, trans_mono)
#### Apply relative risks
combo_ind <- setdiff(1:nrow(p_id), mono_ind)
lrr_se <- (log(.710) - log(.365))/(2 * qnorm(.975))
rr <- rlnorm(n_samples, meanlog = log(.509), sdlog = lrr_se)
rr_indices <- list( # Indices of transition matrix to apply RR to
c(1, 2), c(1, 3), c(1, 4),
c(2, 3), c(2, 4),
c(3, 4)
)
rr_mat <- matrix(rr, nrow = n_samples, ncol = length(rr_indices))
p[,, combo_ind] <- apply_rr(p[, , mono_ind],
rr = rr_mat,
index = rr_indices)
tp <- tparams_transprobs(p, p_id)
## Utility
utility_tbl <- stateval_tbl(
data.table(
state_id = 1:3,
est = c(1, 1, 1)
),
dist = "fixed"
)
## Costs
drugcost_tbl <- stateval_tbl(
data.table(
strategy_id = c(1, 1, 2, 2),
time_start = c(0, 2, 0, 2),
est = c(2278, 2278, 2278 + 2086.50, 2278)
),
dist = "fixed"
)
dmedcost_tbl <- stateval_tbl(
data.table(
state_id = 1:3,
mean = c(A = 1701, B = 1774, C = 6948),
se = c(A = 1701, B = 1774, C = 6948)
),
dist = "gamma"
)
cmedcost_tbl <- stateval_tbl(
data.table(
state_id = 1:3,
mean = c(A = 1055, B = 1278, C = 2059),
se = c(A = 1055, B = 1278, C = 2059)
),
dist = "gamma"
)
# (2) Simulation
## Constructing the economic model
### Transition probabilities
transmod <- CohortDtstmTrans$new(params = tp)
### Utility
utilitymod <- create_StateVals(utility_tbl,
hesim_data = hesim_dat,
n = n_samples)
### Costs
drugcostmod <- create_StateVals(drugcost_tbl,
hesim_data = hesim_dat,
n = n_samples)
dmedcostmod <- create_StateVals(dmedcost_tbl,
hesim_data = hesim_dat,
n = n_samples)
cmedcostmod <- create_StateVals(cmedcost_tbl,
hesim_data = hesim_dat,
n = n_samples)
costmods <- list(drug = drugcostmod,
direct_medical = dmedcostmod,
community_medical = cmedcostmod)
### Economic model
econmod <- CohortDtstm$new(trans_model = transmod,
utility_model = utilitymod,
cost_models = costmods)
## Simulating outcomes
econmod$sim_stateprobs(n_cycles = 20)
autoplot(econmod$stateprobs_, ci = TRUE, ci_style = "ribbon",
labels = labs)
econmod$sim_qalys(dr = 0, integrate_method = "riemann_right")
econmod$sim_costs(dr = 0.06, integrate_method = "riemann_right")
# (3) Decision analysis
ce_sim <- econmod$summarize()
wtp <- seq(0, 25000, 500)
cea_pw_out <- cea_pw(ce_sim, comparator = 1, dr_qalys = 0, dr_costs = .06,
k = wtp)
format(icer(cea_pw_out))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.