Package 'hesim'

Description

Elements of transition probability matrices are multiplied by relative risks and the transition probability matrices are adjusted so that rows sum to 1. Operations are vectorized and each relative risk is multiplied by every transition matrix (stored in 3-dimensional arrays).

Usage

apply_rr(x, rr, index, complement = NULL)

Arguments

Details

This function is useful for applying relative treatment effects measured using relative risks to an existing transition probability matrix. For example, a transition probability matrix for the reference treatment strategy may exist or have been estimated from the data. Relative risks estimated from a meta-analysis or network meta-analysis can then be applied to the reference transition probability matrix. If the number of rows in rr exceeds x, then the arrays in x are recycled to the number of rows in rr, which facilitates the application of relative risks from multiple treatment strategies to a reference treatment.

Value

A 3-dimensional array where each slice contains matrices of the same dimension as each matrix in x and the number of slices is equal to the number of rows in rr.

Examples

p_12 <- c(.7, .5)
p_23 <- c(.1, .2)
x <- as_array3(tpmatrix(
  C, p_12, .1,
  0, C,     p_23,
  0, 0,     1
))

# There are the same number of relative risk rows and transition probability matrices
rr_12 <- runif(2, .8, 1)
rr_13 <- runif(2, .9, 1)
rr <- cbind(rr_12, rr_13)
apply_rr(x, rr, 
         index = list(c(1, 2), c(1, 3)),
         complement = list(c(1, 1), c(2, 2)))
         
# There are more relative risk rows than transition probability matrices
rr_12 <- runif(4, .8, 1)
rr_13 <- runif(4, .9, 1)
rr <- cbind(rr_12, rr_13)
apply_rr(x, rr, 
         index = list(c(1, 2), c(1, 3)),
         complement = list(c(1, 1), c(2, 2)))

Description

Convert a 2-dimensional tabular object where each row stores a flattened square matrix to a 3-dimensional array of square matrices and vice versa. This allows multiple transition matrices to be stored as either tabular objects (e.g., matrices, data frames, etc) or as arrays.

Usage

as_array3(x)

as_tbl2(
  x,
  output = c("data.table", "data.frame", "matrix", "tpmatrix"),
  prefix = "",
  sep = "_"
)

Arguments

Value

For as_array3() a 3-dimensional array of square matrices; for as_tbl2() a 2-dimensional tabular object as specified by output.

See Also

tpmatrix

Examples

p_12 <- c(.7, .6)
pmat <- tpmatrix(
 C, p_12,
 0, 1
)
pmat

as_array3(pmat)
as_array3(as.matrix(pmat))
as_tbl2(as_array3(pmat))
as_tbl2(as_array3(pmat), prefix = "p_", sep = ".")

Description

Convert a multi-state dataset with irreversible transitions containing 3 health states to a dataset with one row per patient and progression-free survival (PFS) and overall survival (OS) time-to-event outcomes.

Usage

as_pfs_os(
  data,
  patient_vars,
  status = "status",
  time_stop = "time_stop",
  transition = "transition_id"
)

Arguments

Value

A data.table with one row per patient containing each variable in patient_vars as well as a time variable and status indicator for both PFS (pfs_status, pfs_time) and OS (os_time, os_status).

Examples

as_pfs_os(onc3, patient_vars = c("patient_id", "strategy_name", "female", "age"))

Description

Creates a data.table that combines the transition probability matrices and ID variables from a tparams_transprobs object. This is often useful for debugging.

Usage

## S3 method for class 'tparams_transprobs'
as.data.table(x, ..., prefix = "prob_", sep = "_", long = FALSE)

Arguments

Value

The output always contains columns for the ID variables and the transition probabilities, but the form depends on on the long argument. If FALSE, then a data.table with one row for each transition probability matrix is returned; otherwise, the data.table contains one row for each transition and columns from (the state being transitioned from) and to (the state being transitioned to) are added.

See Also

tparams_transprobs()

Examples

# Create tparams_transprobs object
hesim_dat <- hesim_data(strategies = data.frame(strategy_id = 1:2),
                        patients = data.frame(patient_id = 1:3))
input_data <- expand(hesim_dat, by = c("strategies", "patients"))    
tpmat_id <- tpmatrix_id(input_data, n_samples = 2)      
p_12 <- runif(nrow(tpmat_id), .6, .7) + 
  .05 * (tpmat_id$strategy_id == 2)
tpmat <- tpmatrix(
  C, p_12,
  0, 1
)
tprobs <- tparams_transprobs(tpmat, tpmat_id)

# Convert to data.table in "wide" format
as.data.table(tprobs)
as.data.table(tprobs, prefix = "")
as.data.table(tprobs, prefix = "", sep = ".")

# Convert to data.table in "long: format
as.data.table(tprobs, long = TRUE)

Description

Quickly plot state probabilities stored in a stateprobs object.

Usage

## S3 method for class 'stateprobs'
autoplot(
  object,
  labels = NULL,
  ci = FALSE,
  prob = 0.95,
  ci_style = c("ribbon", "line"),
  geom_alpha = 0.3,
  ...
)

Arguments

Value

A ggplot object.

Note

If there are multiple patients/groups, then state probabilities are averaged across patients/groups (using the weights in patient_wt if available) prior to plotting.

See Also

Psm for an example.

Description

Quickly plot survival curves stored in a survival object.

Usage

## S3 method for class 'survival'
autoplot(
  object,
  labels = NULL,
  ci = FALSE,
  prob = 0.95,
  ci_style = c("ribbon", "line"),
  geom_alpha = 0.3,
  ...
)

Arguments

Value

A ggplot object.

Note

If there are multiple patients, then survival probabilities are averaged across patients (using the weights in patient_wt if available) prior to plotting.

See Also

Psm for an example.

Description

An object that summarizes simulated measures of clinical effectiveness and costs from a simulation model for use in a cost-effectiveness analysis.

Format

A list containing two elements:

costs

Total (discounted) costs by category.

qalys

(Discounted) quality-adjusted life-years.

Costs

The costs data.table contains the following columns:

category

The cost category.

dr

The discount rate.

sample

A randomly sampled parameter set from the probabilistic sensitivity analysis (PSA)

strategy_id

The treatment strategy ID.

grp_id

An optional column denoting a subgroup. If not included, it is assumed that a single subgroup is being analyzed.

costs

Costs.

Quality-adjusted life-years

The qalys data.table contains the following columns:

dr

The discount rate.

sample

A randomly sampled parameter set from the probabilistic sensitivity analysis (PSA)

strategy_id

The treatment strategy ID.

grp_id

An optional column denoting a subgroup. If not included, it is assumed that a single subgroup is being analyzed.

qalys

Quality-adjusted life-years

Description

Conduct cost-effectiveness analysis (CEA) given output of an economic model; that is, summarize a probabilistic sensitivity analysis (PSA), possibly by subgroup.

• cea() computes the probability that each treatment is most cost-effective, output for a cost-effectiveness acceptability frontier, the expected value of perfect information, and the net monetary benefit for each treatment.

• cea_pw() conducts pairwise CEA by comparing strategies to a comparator. Computed quantities include the incremental cost-effectiveness ratio, the incremental net monetary benefit, output for a cost-effectiveness plane, and output for a cost-effectiveness acceptability curve.

Usage

cea(x, ...)

cea_pw(x, ...)

## Default S3 method:
cea(x, k = seq(0, 2e+05, 500), sample, strategy, grp = NULL, e, c, ...)

## Default S3 method:
cea_pw(
  x,
  k = seq(0, 2e+05, 500),
  comparator,
  sample,
  strategy,
  grp = NULL,
  e,
  c,
  ...
)

## S3 method for class 'ce'
cea(x, k = seq(0, 2e+05, 500), dr_qalys, dr_costs, ...)

## S3 method for class 'ce'
cea_pw(x, k = seq(0, 2e+05, 500), comparator, dr_qalys, dr_costs, ...)

Arguments

Value

cea() returns a list of four data.table elements.

summary

A data.table of the mean, 2.5% quantile, and 97.5% quantile by strategy and group for clinical effectiveness and costs.

mce

The probability that each strategy is the most effective treatment for each group for the range of specified willingness to pay values. In addition, the column best denotes the optimal strategy (i.e., the strategy with the highest expected net monetary benefit), which can be used to plot the cost-effectiveness acceptability frontier (CEAF).

evpi

The expected value of perfect information (EVPI) by group for the range of specified willingness to pay values. The EVPI is computed by subtracting the expected net monetary benefit given current information (i.e., the strategy with the highest expected net monetary benefit) from the expected net monetary benefit given perfect information.

nmb

The mean, 2.5% quantile, and 97.5% quantile of net monetary benefits for the range of specified willingness to pay values.

cea_pw also returns a list of four data.table elements:

summary

A data.table of the mean, 2.5% quantile, and 97.5% quantile by strategy and group for incremental clinical effectiveness and costs.

delta

Incremental effectiveness and incremental cost for each simulated parameter set by strategy and group. Can be used to plot a cost-effectiveness plane.

ceac

Values needed to plot a cost-effectiveness acceptability curve by group. The CEAC plots the probability that each strategy is more cost-effective than the comparator for the specified willingness to pay values.

inmb

The mean, 2.5% quantile, and 97.5% quantile of incremental net monetary benefits for the range of specified willingness to pay values.

Examples

library("data.table")
library("ggplot2")
theme_set(theme_bw())

# Simulation output
n_samples <- 30

sim <- data.table(sample = rep(seq(n_samples), 4),
                  c = c(rlnorm(n_samples, 5, .1), rlnorm(n_samples, 5, .1),
                        rlnorm(n_samples, 11, .1), rlnorm(n_samples, 11, .1)),
                  e = c(rnorm(n_samples, 8, .2), rnorm(n_samples, 8.5, .1),
                        rnorm(n_samples, 11, .6), rnorm(n_samples, 11.5, .6)),
                  strategy_id = rep(1:2, each = n_samples * 2),
                  grp_id = rep(rep(1:2, each = n_samples), 2)
)

# Cost-effectiveness analysis
cea_out <- cea(sim, k = seq(0, 200000, 500), sample = "sample", 
               strategy = "strategy_id", grp = "grp_id", 
               e = "e", c = "c")
names(cea_out)

## Some sample output
## The probability that each strategy is the most cost-effective 
## in each group with a willingness to pay of 20,000
cea_out$mce[k == 20000]

# Pairwise cost-effectiveness analysis
cea_pw_out <-  cea_pw(sim,  k = seq(0, 200000, 500), comparator = 1,
                      sample = "sample", strategy = "strategy_id", 
                      grp = "grp_id", e = "e", c = "c")
names(cea_pw_out)

## Some sample output
## The cost-effectiveness acceptability curve
head(cea_pw_out$ceac[k >= 20000])

# Summarize the incremental cost-effectiveness ratio
labs <- list(strategy_id = c("Strategy 1" = 1, "Strategy 2" = 2),
             grp_id = c("Group 1" = 1, "Group 2" = 2))
format(icer(cea_pw_out, labels = labs))

# Plots
plot_ceplane(cea_pw_out, label = labs)
plot_ceac(cea_out, label = labs)
plot_ceac(cea_pw_out, label = labs)
plot_ceaf(cea_out, label = labs)
plot_evpi(cea_out, label = labs)

Description

Simulate outcomes from a cohort discrete time state transition model.

Format

An R6::R6Class object.

Public fields

Methods

Public methods

• CohortDtstm$new()

• CohortDtstm$sim_stateprobs()

• CohortDtstm$sim_qalys()

• CohortDtstm$sim_costs()

• CohortDtstm$summarize()

• CohortDtstm$clone()

Method new()

Create a new CohortDtstm object.

Usage

CohortDtstm$new(trans_model = NULL, utility_model = NULL, cost_models = NULL)

Arguments

Returns

A new CohortDtstm object.

Method sim_stateprobs()

Simulate health state probabilities using CohortDtstmTrans$sim_stateprobs().

Usage

CohortDtstm$sim_stateprobs(n_cycles)

Arguments

Returns

An instance of self with simulated output of class stateprobs stored in stateprobs_.

Method sim_qalys()

Simulate quality-adjusted life-years (QALYs) as a function of stateprobs_ and utility_model. See sim_qalys() for details.

Usage

CohortDtstm$sim_qalys(
  dr = 0.03,
  integrate_method = c("trapz", "riemann_left", "riemann_right"),
  lys = TRUE
)

Arguments

Returns

An instance of self with simulated output of class qalys stored in qalys_.

Method sim_costs()

Simulate costs as a function of stateprobs_ and cost_models. See sim_costs() for details.

Usage

CohortDtstm$sim_costs(
  dr = 0.03,
  integrate_method = c("trapz", "riemann_left", "riemann_right")
)

Arguments

Returns

An instance of self with simulated output of class costs stored in costs_.

Method summarize()

Summarize costs and QALYs so that cost-effectiveness analysis can be performed. See summarize_ce().

Usage

CohortDtstm$summarize(by_grp = FALSE)

Arguments

Method clone()

The objects of this class are cloneable with this method.

Usage

CohortDtstm$clone(deep = FALSE)

Arguments

References

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.

See Also

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.

Examples

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))

Description

Simulate health state transitions in a cohort discrete time state transition model.

Format

An R6::R6Class object.

Public fields

Active bindings

Methods

Public methods

• CohortDtstmTrans$new()

• CohortDtstmTrans$sim_stateprobs()

• CohortDtstmTrans$clone()

Method new()

Create a new CohortDtstmTrans object.

Usage

CohortDtstmTrans$new(
  params,
  input_data = NULL,
  trans_mat = NULL,
  start_stateprobs = NULL,
  cycle_length = 1,
  absorbing = NULL
)

Arguments

Returns

A new CohortDtstmTrans object.

Method sim_stateprobs()

Simulate probability of being in each health state during each model cycle.

Usage

CohortDtstmTrans$sim_stateprobs(n_cycles)

Arguments

Returns

An object of class stateprobs.

Method clone()

The objects of this class are cloneable with this method.

Usage

CohortDtstmTrans$clone(deep = FALSE)

Arguments

See Also

create_CohortDtstmTrans() creates a CohortDtstmTrans object from either a fitted statistical model or a parameter object. A complete economic model can be implemented with the CohortDtstm class.

Examples

library("msm")
library("data.table")
set.seed(101)

# We consider two examples that have the same treatment strategies and patients.
# One model is parameterized by fitting a multi-state model with the "msm"
# package; in the second model, the parameters are entered "manually" with
# a "params_mlogit_list" object.

# MODEL SETUP
strategies <- data.table(
  strategy_id = c(1, 2, 3),
  strategy_name = c("SOC", "New 1", "New 2")
)
patients <- data.table(patient_id = 1:2)
hesim_dat <- hesim_data(
  strategies = strategies,
  patients = patients
)

# EXAMPLE #1: msm
## Fit multi-state model with panel data via msm
qinit <- rbind(
  c(0, 0.28163, 0.01239),
  c(0, 0, 0.10204),
  c(0, 0, 0)
)
fit <- msm(state_id ~ time, subject = patient_id,
           data = onc3p[patient_id %in% sample(patient_id, 100)],
           covariates = list("1-2" =~ strategy_name),
           qmatrix = qinit)

## Simulation model
transmod_data <- expand(hesim_dat)
transmod <- create_CohortDtstmTrans(fit,
                                    input_data = transmod_data,
                                    cycle_length = 1/2,
                                    fixedpars = 2,
                                    n = 2)
transmod$sim_stateprobs(n_cycles = 2)

# EXAMPLE #2: params_mlogit_list
## Input data
transmod_data[, intercept := 1]
transmod_data[, new1 := ifelse(strategy_name == "New 1", 1, 0)]
transmod_data[, new2 := ifelse(strategy_name == "New 2", 1, 0)]

## Parameters
n <- 10
transmod_params <- params_mlogit_list(
  
  ## Transitions from stable state (stable -> progression, stable -> death)
  stable = params_mlogit(
    coefs = list(
      progression = data.frame(
        intercept = rnorm(n, -0.65, .1),
        new1 = rnorm(n, log(.8), .02),
        new2 = rnorm(n, log(.7, .02))
      ),
      death = data.frame(
        intercept = rnorm(n, -3.75, .1),
        new1 = rep(0, n),
        new2 = rep(0, n)
      )
    )
  ),
  
  ## Transition from progression state (progression -> death)
  progression = params_mlogit(
    coefs = list(
      death = data.frame(
        intercept = rnorm(n, 2.45, .1),
        new1 = rep(0, n),
        new2 = rep(0, n)
      )
    )
  )
)
transmod_params

## Simulation model
tmat <- rbind(c(0, 1, 2),
              c(NA, 0, 1),
              c(NA, NA, NA))
transmod <- create_CohortDtstmTrans(transmod_params, 
                                    input_data = transmod_data,
                                    trans_mat = tmat, cycle_length = 1)
transmod$sim_stateprobs(n_cycles = 2)

Description

An object of class costs returned from methods ⁠$sim_costs()⁠ in model classes that store simulated costs.

Components

A costs object inherits from data.table and contains 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.

category

The cost category (e.g., drug costs, medical costs, etc).

costs

The simulated cost values.

Description

A generic function for creating an object of class CohortDtstm.

Usage

create_CohortDtstm(object, ...)

## S3 method for class 'model_def'
create_CohortDtstm(
  object,
  input_data,
  cost_args = NULL,
  utility_args = NULL,
  ...
)

Arguments

Description

A generic function for creating an object of class CohortDtstmTrans.

Usage

create_CohortDtstmTrans(object, ...)

## S3 method for class 'multinom_list'
create_CohortDtstmTrans(
  object,
  input_data,
  trans_mat,
  n = 1000,
  uncertainty = c("normal", "none"),
  ...
)

## S3 method for class 'msm'
create_CohortDtstmTrans(
  object,
  input_data,
  cycle_length,
  n = 1000,
  uncertainty = c("normal", "none"),
  ...
)

## S3 method for class 'params_mlogit_list'
create_CohortDtstmTrans(object, input_data, trans_mat, ...)

Arguments

Details

Disease models may either be created from a fitted statistical model or from a parameter object. In the case of the former, input_data is a data frame like object that is used to look for variables from the statistical model that are required for simulation. In this sense, input_data is very similar to the newdata argument in most predict() methods (e.g., see predict.lm()). In other words, variables used in the formula of the statistical model must also be in input_data.

In the case of the latter, the columns of input_data must be named in a manner that is consistent with the parameter object. In the typical case (e.g., with params_surv or params_mlogit), the parameter object contains coefficients from a regression model, usually stored as matrix where rows index parameter samples (i.e., for a probabilistic sensitivity analysis) and columns index model terms. In such instances, there must be one column from input_data with the same name as each model term in the coefficient matrix; that is, the columns in input_data are matched with the columns of the coefficient matrices by name. If there are model terms in the coefficient matrices that are not contained in input_data, then an error will be thrown.

See Also

See CohortDtstmTrans for examples.

Description

A generic function for creating an object of class IndivCtstmTrans.

Usage

create_IndivCtstmTrans(object, ...)

## S3 method for class 'flexsurvreg_list'
create_IndivCtstmTrans(
  object,
  input_data,
  trans_mat,
  clock = c("reset", "forward"),
  n = 1000,
  uncertainty = c("normal", "none"),
  ...
)

## S3 method for class 'flexsurvreg'
create_IndivCtstmTrans(
  object,
  input_data,
  trans_mat,
  clock = c("reset", "forward"),
  n = 1000,
  uncertainty = c("normal", "none"),
  ...
)

## S3 method for class 'params_surv'
create_IndivCtstmTrans(
  object,
  input_data,
  trans_mat,
  clock = c("reset", "forward", "mix", "mixt"),
  reset_states = NULL,
  transition_types = NULL,
  ...
)

## S3 method for class 'params_surv_list'
create_IndivCtstmTrans(
  object,
  input_data,
  trans_mat,
  clock = c("reset", "forward", "mix", "mixt"),
  reset_states = NULL,
  transition_types = NULL,
  ...
)

Arguments

Details

Disease models may either be created from a fitted statistical model or from a parameter object. In the case of the former, input_data is a data frame like object that is used to look for variables from the statistical model that are required for simulation. In this sense, input_data is very similar to the newdata argument in most predict() methods (e.g., see predict.lm()). In other words, variables used in the formula of the statistical model must also be in input_data.

In the case of the latter, the columns of input_data must be named in a manner that is consistent with the parameter object. In the typical case (e.g., with params_surv or params_mlogit), the parameter object contains coefficients from a regression model, usually stored as matrix where rows index parameter samples (i.e., for a probabilistic sensitivity analysis) and columns index model terms. In such instances, there must be one column from input_data with the same name as each model term in the coefficient matrix; that is, the columns in input_data are matched with the columns of the coefficient matrices by name. If there are model terms in the coefficient matrices that are not contained in input_data, then an error will be thrown.

Value

Returns an R6::R6Class object of class IndivCtstmTrans.

See Also

See IndivCtstmTrans and IndivCtstm for examples.

Description

create_params is a generic function for creating an object containing parameters from a fitted statistical model. If uncertainty != "none", then random samples from suitable probability distributions are returned.

Usage

create_params(object, ...)

## S3 method for class 'lm'
create_params(object, n = 1000, uncertainty = c("normal", "none"), ...)

## S3 method for class 'multinom'
create_params(object, n = 1000, uncertainty = c("normal", "none"), ...)

## S3 method for class 'multinom_list'
create_params(object, n = 1000, uncertainty = c("normal", "none"), ...)

## S3 method for class 'flexsurvreg'
create_params(object, n = 1000, uncertainty = c("normal", "none"), ...)

## S3 method for class 'flexsurvreg_list'
create_params(object, n = 1000, uncertainty = c("normal", "none"), ...)

## S3 method for class 'partsurvfit'
create_params(
  object,
  n = 1000,
  uncertainty = c("normal", "bootstrap", "none"),
  max_errors = 0,
  silent = FALSE,
  ...
)

Arguments

Value

An object prefixed by params_. Mapping between create_params and the classes of the returned objects are:

• create_params.lm -> params_lm

• create_params.multinom -> params_mlogit

• create_params.multinom_list -> params_mlogit_list

• create_params.flexsurvreg -> params_surv

• create_params.flexsurvreg_list -> params_surv_list

• create_params.partsurvfit -> params_surv_list

See Also

These methods are typically used alongside create_input_mats() to create model objects as a function of input data and a fitted statistical model. For instance, create_PsmCurves() creates the survival model for a partitioned survival model, create_IndivCtstmTrans() creates the transition model for an individual continuous time state transition model, create_CohortDtstmTrans() creates the transition model for a cohort discrete time state transition model, and create_StateVals() creates a health state values model.

Examples

# create_params.lm
fit <- lm(costs ~ female, data = psm4_exdata$costs$medical)
n <- 5
params_lm <- create_params(fit, n = n)
head(params_lm$coefs)
head(params_lm$sigma)

# create_params.flexsurvreg
library("flexsurv")
fit <- flexsurvreg(formula = Surv(futime, fustat) ~ 1, 
                   data = ovarian, dist = "weibull")
n <- 5
params_surv_wei <- create_params(fit, n = n)
print(params_surv_wei$dist)
head(params_surv_wei$coefs)

Description

A generic function for creating a PsmCurves object.

Usage

create_PsmCurves(object, ...)

## S3 method for class 'flexsurvreg_list'
create_PsmCurves(
  object,
  input_data,
  n = 1000,
  uncertainty = c("normal", "bootstrap", "none"),
  est_data = NULL,
  ...
)

## S3 method for class 'params_surv_list'
create_PsmCurves(object, input_data, ...)

Arguments

Details

Disease models may either be created from a fitted statistical model or from a parameter object. In the case of the former, input_data is a data frame like object that is used to look for variables from the statistical model that are required for simulation. In this sense, input_data is very similar to the newdata argument in most predict() methods (e.g., see predict.lm()). In other words, variables used in the formula of the statistical model must also be in input_data.

In the case of the latter, the columns of input_data must be named in a manner that is consistent with the parameter object. In the typical case (e.g., with params_surv or params_mlogit), the parameter object contains coefficients from a regression model, usually stored as matrix where rows index parameter samples (i.e., for a probabilistic sensitivity analysis) and columns index model terms. In such instances, there must be one column from input_data with the same name as each model term in the coefficient matrix; that is, the columns in input_data are matched with the columns of the coefficient matrices by name. If there are model terms in the coefficient matrices that are not contained in input_data, then an error will be thrown.

Value

Returns an R6Class object of class PsmCurves.

See Also

See PsmCurves and Psm for examples. PsmCurves provides an example in which a model is parameterized both with (via create_PsmCurves.flexsurvreg_list()) and without (via create_PsmCurves.params_surv_list()) access to patient-level data. The Psm example shows how state probabilities, costs, and utilities can be computed from predicted survival curves.

Description

create_StateVals() is a generic function for creating an object of class StateVals from a fitted statistical model or a stateval_tbl object.

Usage

create_StateVals(object, ...)

## S3 method for class 'lm'
create_StateVals(
  object,
  input_data = NULL,
  n = 1000,
  uncertainty = c("normal", "none"),
  ...
)

## S3 method for class 'stateval_tbl'
create_StateVals(object, hesim_data = NULL, n = 1000, ...)

Arguments

Details

If object is a stateval_tbl, then a hesim_data object is used to specify treatment strategies, patients, and/or health states not included as columns in the table, or, to match patients in the table to groups. Not required if the table includes one row for each treatment strategy, patient, and health state combination. Patients are matched to groups by specifying both a patient_id and a grp_var column in the patients table.

Value

A StateVals object.

See Also

See StateVals for documentation of the class and additional examples. An example use case for create_StateVals.stateval_tbl() is provided in the stateval_tbl() documentation.

Examples

set.seed(10)

# EXAMPLE FOR `create_statevals.lm()`
## Simple example comparing two treatment strategies where
## medical costs vary by sex and health state

## Setup model
hesim_dat <- hesim_data(
  strategies = data.frame(strategy_id = c(1, 2)),
  patients = data.frame(
    patient_id = c(1, 2),
    female = c(1, 0)
  ),
  states = data.frame(
    state_id = c(1, 2, 3),
    state_name = c("state1", "state2", "state3")
  )
)

## Fit model
medcost_estimation_data <- psm4_exdata$costs$medical
medcost_estimation_data$time5 <- rbinom(nrow(medcost_estimation_data), 
                                        1, .5) # Illustrative time dummy
medcost_fit <- lm(costs ~ female + state_name + time5, 
                  data = medcost_estimation_data)

## Create medical cost model
### Allow medical costs to vary across time in addition to by patient and 
### health state
medcost_times <- time_intervals(
  data.frame(time_start = c(0, 3, 5),
            time5 = c(0, 0, 1)) # Time dummy corresponds to time > 5
)
medcost_input_data <- expand(hesim_dat, 
                             by = c("strategies", "patients", "states"),
                             times = medcost_times)
medcost_model <- create_StateVals(medcost_fit, medcost_input_data,
                                  n = 1)

## Explore predictions from medical cost model
### We can assess predictions at multiple time points
medcost_model$sim(t = c(1, 6), type = "predict")

Description

Create a data table of health state transitions from a transition matrix describing the states and transitions in a multi-state model suitable for use with hesim_data.

Usage

create_trans_dt(trans_mat)

Arguments

Value

Returns a data.table::data.table in tidy format with three columns:

transition_id

Health state transition ID.

from

The starting health state.

to

The health state that will be transitions to.

Examples

tmat <- rbind(c(NA, 1, 2),
              c(NA, NA, 3),
              c(NA, NA, NA))
create_trans_dt(tmat)

Description

A model expression is defined by specifying random number generation functions for a probabilistic sensitivity analysis (PSA) and transformations of the sampled parameters as a function of input_data. The unevaluated expressions are evaluated with eval_model() and used to generate the model inputs needed to create an economic model.

Usage

define_model(tparams_def, rng_def, params = NULL, n_states = NULL)

eval_model(x, input_data)

Arguments

Details

eval_model() evaluates the expressions in an object of class model_def returned by define_model() and is, in turn, used within functions that instantiate economic models (e.g., create_CohortDtstm()). The direct output of eval_model() can also be useful for understanding and debugging model definitions, but it is not used directly for simulation.

Economic models are constructed as a function of input data and parameters:

1. Input data: Objects of class expanded_hesim_data consisting of the treatment strategies and patient population.

2. Parameters: The underlying parameter estimates from the literature are first stored in a list (params argument). Random number generation is then used to sample the parameters from suitable probability distributions for the PSA (rng_def argument). Finally, the sampled parameters are transformed as a function of the input data into values (e.g., elements of a transition probability matrix) used for the simulation (tparams_def argument). The params argument can be omitted if the underlying parameters values are defined inside a define_rng() block.

Value

define_model() returns an object of class model_def, which is a list containing the arguments to the function. eval_model() returns a list containing ID variables identifying parameter samples, treatment strategies, patient cohorts, and time intervals; the values of parameters of the transition probability matrix, utilities, and/or cost categories; the number of health states; and the number of random number generation samples for the PSA.

See Also

define_tparams(), define_rng()

Examples

# Data
library("data.table")
strategies <- data.table(strategy_id = 1:2,
                         strategy_name = c("Monotherapy", "Combination therapy"))
patients <- data.table(patient_id = 1)
hesim_dat <- hesim_data(strategies = strategies,
                       patients = patients)
data <- expand(hesim_dat)

# Model parameters
rng_def <- define_rng({
  alpha <- matrix(c(1251, 350, 116, 17,
                    0, 731, 512, 15,
                    0, 0, 1312, 437,
                    0, 0, 0, 469),
                  nrow = 4, byrow = TRUE)
  rownames(alpha) <- colnames(alpha) <- c("A", "B", "C", "D")
  lrr_mean <- log(.509)
  lrr_se <- (log(.710) - log(.365))/(2 * qnorm(.975))
  
  list(
    p_mono = dirichlet_rng(alpha),
    rr_comb = lognormal_rng(lrr_mean, lrr_se),
    u = 1,
    c_zido = 2278,
    c_lam = 2086.50,
    c_med = gamma_rng(mean = c(A = 2756, B = 3052, C = 9007),
                      sd = c(A = 2756, B = 3052, C = 9007))
  )
}, n = 2)

tparams_def <- define_tparams({
  rr = ifelse(strategy_name == "Monotherapy", 1, rr_comb)
  list(
    tpmatrix = tpmatrix(
      C, p_mono$A_B * rr, p_mono$A_C * rr, p_mono$A_D * rr,
      0, C, p_mono$B_C * rr, p_mono$B_D * rr,
      0, 0, C, p_mono$C_D * rr,
      0, 0, 0, 1),
    utility = u,
    costs = list(
      drug = ifelse(strategy_name == "Monotherapy",
                    c_zido, c_zido + c_lam),
      medical = c_med
    ) 
  )
})

# Simulation
## Define the economic model
model_def <- define_model(
  tparams_def = tparams_def,
  rng_def = rng_def)

### Evaluate the model expression to generate model inputs
### This can be useful for understanding the output of a model expression
eval_model(model_def, data)

## Create an economic model with a factory function
econmod <- create_CohortDtstm(model_def, data)

Description

Random number generation expressions are used to randomly sample model parameters from suitable distributions for probabilistic sensitivity analysis. These functions are typically used when evaluating an object of class model_def defined using define_model().

Usage

define_rng(expr, n = 1, ...)

eval_rng(x, params = NULL, check = TRUE)

Arguments

Details

hesim contains a number of random number generation functions that return parameter samples in convenient formats and do not typically require the number of samples, n, as arguments (see rng_distributions). The random number generation expressions are evaluated using eval_rng() and used within expr in define_rng(). If a multivariate object is returned by eval_rng(), then the rows are random samples and columns are distinct parameters (e.g., costs for each health state, elements of a transition probability matrix).

Value

define_rng() returns an object of class rng_def, which is a list containing the unevaluated random number generation expressions passed to expr, n, and any additional arguments passed to ... . eval_rng() evaluates the rng_def object and returns an eval_rng object containing the evaluated expression.

See Also

Parameters can be conveniently sampled from probability distributions using a number of random number generation functions (see rng_distributions). An economic model can be created with create_CohortDtstm() by using define_rng() (or a previously evaluated eval_rng object) alongside define_tparams() to define a model with define_model(). It can be useful to summarize an evaluated expression with summary.eval_rng().

Examples

params <- list(
  alpha = matrix(c(75, 25, 33, 67), byrow = TRUE, ncol = 2),
  inptcost_mean = c(A = 900, B = 1500, C = 2000),
  outptcost_mean = matrix(c(300, 600, 800,
                            400, 700, 700),
                           ncol = 3, byrow = TRUE)
)
rng_def <- define_rng({
  aecost_mean <- c(500, 800, 1000) # Local object not 
                                   # not returned by eval_rng()
  list( # Sampled values of parameters returned by eval_rng()
    p = dirichlet_rng(alpha), # Default column names
    inptcost = gamma_rng(mean = inptcost_mean, # Column names based on 
                         sd = inptcost_mean),  # named vector
    outptcost = outptcost_mean, # No column names because
                                # outptcost_mean has none.
    aecost = gamma_rng(mean = aecost_mean, # Explicit naming of columns
                       sd = aecost_mean,
                       names = aecost_colnames)
  )
}, n = 2, aecost_colnames = c("A", "B", "C")) # Add aecost_colnames to environment
params_sample <- eval_rng(x = rng_def, params)
summary(params_sample)
params_sample

Description

Transformed parameter expressions are used to transform the parameter values sampled with eval_rng() as a function of input data (treatment strategies and patients) and time intervals. These functions are used when evaluating an object of class model_def defined using define_model(). The transformed parameters are ultimately converted into tparams objects and used to simulate outcomes with an economic model.

Usage

define_tparams(expr, times = NULL, ...)

eval_tparams(x, input_data, rng_params)

Arguments

Details

define_tparams() is evaluated when creating economic models as a function of model_def objects defined with define_model(). Operations are "vectorized" in the sense that they are performed for each unique combination of input_data and params. expr is evaluated in an environment including each variable from input_data, all elements of rng_params, and a variable time containing the values from times. The time variable can be used to create models where parameters vary as a function of time. eval_tparams() is not exported and is only meant for use within eval_model().

Value

define_tparams() returns an object of class tparams_def, which is a list containing the unevaluated "transformation" expressions passed to expr, times, and any additional arguments passed to ... . eval_tparams() evaluates the tparams_def object and should return a list of transformed parameter objects.

See Also

define_model(), define_rng()

Description

An object of class disprog returned from methods ⁠$sim_disease()⁠ in model classes. It contains simulated trajectories through a multi-state model.

Components

A disprog object inherits from data.table and contains the following columns:

sample

A random sample from the PSA.

strategy_id

The treatment strategy ID.

patient_id

The patient ID.

from

The health state ID transitioned from.

to

The health state ID transitioned to.

final

An indicator equal to 1 if a patient is in their final health state during the simulation and 0 otherwise.

time_start

The time at the start of the interval.

time_stop

The time at the end of the interval.

The object also contains size and absorbing attributes. The size attribute is a numeric vector with the elements n_samples, n_strategies, n_patients, and n_states denoting the number of samples, treatment strategies, patients, and health states The absorbing attribute is a numeric vector containing the absorbing health states; i.e., the health states that cannot be transitioned from. Operationally, an absorbing state is a row in a transition matrix (as in the trans_mat field of the IndivCtstmTrans class) with all NAs.

See Also

A disease progression object can be simulated with either the IndivCtstm or IndivCtstmTrans classes.

Description

Create a data table in long format from all combinations of specified tables from an object of class hesim_data and optionally time intervals. See "Details" for an explanation of how the expansion is done.

Usage

## S3 method for class 'hesim_data'
expand(object, by = c("strategies", "patients"), times = NULL)

Arguments

Details

This function is similar to expand.grid(), but works for data frames or data tables. Specifically, it creates a data.table from all combinations of the supplied tables in object and optionally the start of times intervals in times. The supplied tables are determined using the by argument. The resulting dataset is sorted by prioritizing ID variables as follows: (i) strategy_id, (ii) patient_id, (iii) the health-related ID variable (either state_id or transition_id, and (iv) the time intervals from times.

Value

An object of class expanded_hesim_data, which is a data.table with an "id_vars" attribute containing the names of the ID variables in the data table and, if times is not NULL, a time_intervals object derived from times.

Examples

strategies <- data.frame(strategy_id = c(1, 2))
patients <- data.frame(patient_id = seq(1, 3), age = c(65, 50, 75),
                          gender = c("Female", "Female", "Male"))
states <- data.frame(state_id =  seq(1, 3),
                     state_var = c(2, 1, 9))
hesim_dat <- hesim_data(strategies = strategies,
                        patients = patients,
                        states = states)
expand(hesim_dat, by = c("strategies", "patients"))
expand(hesim_dat, by = c("strategies", "patients"),
       times = c(0, 2, 10))

Description

This is a wrapper around msm::MatrixExp() that computes the exponential of multiple square matrices.

Usage

expmat(x, t = 1, ...)

Arguments

Details

This function is most useful when exponentiating transition intensity matrices to produce transition probability matrices. To create transition probability matrices for discrete time state transition models with annual cycles, set t=1. An array of matrices is returned which can be used to create the value element of a tparams_transprobs object. See qmatrix() for an example.

Value

Returns an array of exponentiated matrices. If length(t) > 1, then length(t) arrays are returned for each element in x.

See Also

qmatrix.msm(), qmatrix.data.table()

Description

Draw random samples from a generalized gamma distribution using the parameterization from flexsurv. Written in C++ for speed. Equivalent to flexsurv::rgengamma.

Usage

fast_rgengamma(n, mu = 0, sigma = 1, Q)

Arguments

Value

A vector of random samples from the generalized gamma distribution. The length of the sample is determined by n. The numerical arguments other than n are recycled so that the number of samples is equal to n.

Examples

n <- 1000
m <- 2 ; s <- 1.7; q <- 1
ptm <- proc.time()
r1 <- fast_rgengamma(n, mu = m, sigma = s, Q = q)
proc.time() - ptm
ptm <- proc.time()
library("flexsurv")
r2 <- flexsurv::rgengamma(n, mu = m, sigma = s, Q = q)
proc.time() - ptm
summary(r1)
summary(r2)

Description

Combine flexsurvreg objects into a list.

Usage

flexsurvreg_list(...)

Arguments

Value

An object of class flexsurvreg_list.

Examples

library("flexsurv")
 fit1 <- flexsurv::flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, dist = "weibull")
 fit2 <- flexsurv::flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, dist = "exp")
 fsreg_list <- flexsurvreg_list(wei = fit1, exp = fit2)
 class(fsreg_list)

Description

Get value labels for the ID variables in a hesim_data object and create a list of named vectors that can be passed to formatting and plotting functions. This lets users create nice labels for treatment strategies, subgroups, health states, and/or transitions when presenting results.

Usage

get_labels(
  object,
  strategy = "strategy_name",
  grp = "grp_name",
  state = "state_name",
  transition = "transition_name",
  death_label = "Death"
)

Arguments

Value

A list of named vectors containing the values and labels of variables. The elements of each vector are the values of a variable and the names are the labels. The names of the list are the names of the ID variables.

See Also

hesim_data(), set_labels()

Examples

library("data.table")
strategies <- data.table(
  strategy_id = c(1, 2),
  strategy_name = c("Strategy 1", "Strategy 2")
)
patients <- data.table(
  patient_id = seq(1, 4),
  age = c(50, 55, 60, 65),
  grp_id = c(1, 1, 2, 2),
  grp_name = rep(c("Age 50-59", "Age 60-69"), each = 2)
)
states <- data.table(
  state_id =  seq(1, 2),
  state_name = c("State 1", "State 2")
)
hesim_dat <- hesim_data(
  strategies = strategies,
  patients = patients,
  states = states
)
labs <- get_labels(hesim_dat)
labs

# Pass to set_labels()
d <- data.table(strategy_id = c(1, 1, 2, 2),
                grp_id = c(1, 2, 1, 2))
set_labels(d, labs, new_name = c("strategy_name", "grp_name"))
d

Description

A list of tables required for health economic simulation modeling. This object is used to setup models by defining the treatment strategies, target population, and model structure.

Usage

hesim_data(strategies, patients, states = NULL, transitions = NULL)

Arguments

Value

Returns an object of class hesim_data, which is a list of data tables for health economic simulation modeling.

Note

Each table must either be a data.frame or data.table. All ID variables within each table must be numeric vectors of integers and should be of the form 1,2,...N where N is the number of unique values of the ID variable.

See Also

expand.hesim_data(), get_labels()

Examples

strategies <- data.frame(strategy_id = c(1, 2))
patients <- data.frame(patient_id = seq(1, 3), age = c(65, 50, 75),
                       gender = c("Female", "Female", "Male"))
states <- data.frame(state_id =  seq(1, 3),
                     state_var = c(2, 1, 9))
hesim_dat <- hesim_data(strategies = strategies,
                        patients = patients,
                        states = states)

Description

Generate a tidy table of incremental cost-effectiveness ratios (ICERs) given output from cea_pw() with icer() and format for pretty printing with format.icer().

Usage

icer(x, prob = 0.95, k = 50000, labels = NULL, ...)

## S3 method for class 'icer'
format(
  x,
  digits_qalys = 2,
  digits_costs = 0,
  pivot_from = "strategy",
  drop_grp = TRUE,
  pretty_names = TRUE,
  ...
)

Arguments

Details

Note that icer() will report negative ICERs; however, format() will correctly note whether a treatment strategy is dominated by or dominates the reference treatment.

Value

icer() returns an object of class icer that is a tidy data.table with the following columns:

strategy

The treatment strategy.

grp

The subgroup.

outcome

The outcome metric.

estimate

The point estimate computed as the average across the PSA samples.

lower

The lower limit of the confidence interval.

upper

The upper limit of the confidence interval.

format.icer() formats the table according to the arguments passed.

See Also

cea_pw()

Description

Computes incremental effect for all treatment strategies on outcome variables from a probabilistic sensitivity analysis relative to a comparator.

Usage

incr_effect(x, comparator, sample, strategy, grp = NULL, outcomes)

Arguments

Value

A data.table containing the differences in the values of each variable specified in outcomes between each treatment strategy and the comparator.

Examples

# simulation output
n_samples <- 100
sim <- data.frame(sample = rep(seq(n_samples), 4),
             c = c(rlnorm(n_samples, 5, .1), rlnorm(n_samples, 5, .1),
                    rlnorm(n_samples, 11, .1), rlnorm(n_samples, 11, .1)),
             e = c(rnorm(n_samples, 8, .2), rnorm(n_samples, 8.5, .1),
                   rnorm(n_samples, 11, .6), rnorm(n_samples, 11.5, .6)),
             strategy = rep(paste0("Strategy ", seq(1, 2)),
                           each = n_samples * 2),
             grp = rep(rep(c("Group 1", "Group 2"),
                           each = n_samples), 2)
)
# calculate incremental effect of Strategy 2 relative to Strategy 1 by group
ie <- incr_effect(sim, comparator = "Strategy 1", sample = "sample",
                        strategy = "strategy", grp = "grp", outcomes = c("c", "e"))
head(ie)

Description

Simulate outcomes from an individual-level continuous time state transition model (CTSTM). The class supports "clock-reset" (i.e., semi-Markov), "clock-forward" (i.e., Markov), and mixtures of clock-reset and clock-forward multi-state models as described in IndivCtstmTrans.

Format

An R6::R6Class object.

Public fields

Methods

Public methods

• IndivCtstm$new()

• IndivCtstm$sim_disease()

• IndivCtstm$sim_stateprobs()

• IndivCtstm$sim_qalys()

• IndivCtstm$sim_costs()

• IndivCtstm$summarize()

• IndivCtstm$clone()

Method new()

Create a new IndivCtstm object.

Usage

IndivCtstm$new(trans_model = NULL, utility_model = NULL, cost_models = NULL)

Arguments

Returns

A new IndivCtstm object.

Method sim_disease()

Simulate disease progression (i.e., individual trajectories through a multi-state model) using IndivCtstmTrans$sim_disease().

Usage

IndivCtstm$sim_disease(max_t = 100, max_age = 100, progress = NULL)

Arguments

Returns

An instance of self with simulated output stored in disprog_.

Method sim_stateprobs()

Simulate health state probabilities as a function of time using the simulation output stored in disprog.

Usage

IndivCtstm$sim_stateprobs(t)

Arguments

Returns

An instance of self with simulated output of class stateprobs stored in stateprobs_.

Method sim_qalys()

Simulate quality-adjusted life-years (QALYs) as a function of disprog_ and utility_model.

Usage

IndivCtstm$sim_qalys(
  dr = 0.03,
  type = c("predict", "random"),
  lys = TRUE,
  by_patient = FALSE
)

Arguments

Returns

An instance of self with simulated output of class qalys stored in qalys_.

Method sim_costs()

Simulate costs as a function of disprog_ and cost_models.

Usage

IndivCtstm$sim_costs(
  dr = 0.03,
  type = c("predict", "random"),
  by_patient = FALSE,
  max_t = Inf
)

Arguments

Returns

An instance of self with simulated output of class costs stored in costs_.

Method summarize()

Summarize costs and QALYs so that cost-effectiveness analysis can be performed. See summarize_ce().

Usage

IndivCtstm$summarize(by_grp = FALSE)

Arguments

Method clone()

The objects of this class are cloneable with this method.

Usage

IndivCtstm$clone(deep = FALSE)

Arguments

References

Incerti and Jansen (2021). See Section 2.2 for a mathematical description of an individual-level CTSTM and Section 4.1 for an example in oncology.

See Also

The IndivCtstmTrans documentation describes the class for the transition model and the StateVals documentation describes the class for the cost and utility models. An IndivCtstmTrans object is typically created using create_IndivCtstmTrans().

There are currently two relevant vignettes. vignette("mstate") shows how to parameterize IndivCtstmTrans objects in cases where patient-level data is available by fitting a multi-state models. vignette("markov-inhomogeneous-indiv") shows how the time inhomogeneous Markov cohort model in vignette("markov-inhomogeneous-cohort") can be developed as an individual patient simulation; in doing so, it shows how IndivCtstm models can be used even without access to patient-level data.

Examples

library("flexsurv")

# Treatment strategies, target population, and model structure
strategies <- data.frame(strategy_id = c(1, 2))
patients <- data.frame(patient_id = seq(1, 3),
                       age = c(45, 50, 60),
                       female = c(0, 0, 1))
states <- data.frame(state_id = c(1, 2))
hesim_dat <- hesim_data(strategies = strategies,
                        patients = patients,
                        states = states)

# Parameter estimation
## Multi-state model
tmat <- rbind(c(NA, 1, 2),
              c(3, NA, 4),
              c(NA, NA, NA))
fits <- vector(length = max(tmat, na.rm = TRUE), mode = "list")
surv_dat <- data.frame(mstate3_exdata$transitions)
for (i in 1:length(fits)){
  fits[[i]] <- flexsurvreg(Surv(years, status) ~ factor(strategy_id),
                           data = surv_dat,
                           subset = (trans == i),
                           dist = "weibull")
}
fits <- flexsurvreg_list(fits)

## Utility
utility_tbl <- stateval_tbl(data.frame(state_id = states$state_id,
                                       mean = mstate3_exdata$utility$mean,
                                       se = mstate3_exdata$utility$se),
                            dist = "beta")
## Costs
drugcost_tbl <- stateval_tbl(data.frame(strategy_id = strategies$strategy_id,
                                        est = mstate3_exdata$costs$drugs$costs),
                             dist = "fixed")
medcost_tbl <- stateval_tbl(data.frame(state_id = states$state_id,
                                       mean = mstate3_exdata$costs$medical$mean,
                                       se = mstate3_exdata$costs$medical$se),
                            dist = "gamma")

# Economic model
n_samples = 2

## Construct model
### Transitions
transmod_data <- expand(hesim_dat)
transmod <- create_IndivCtstmTrans(fits, input_data = transmod_data,
                                   trans_mat = tmat,
                                   n = n_samples)

### Utility
utilitymod <- create_StateVals(utility_tbl, n = n_samples, hesim_data = hesim_dat)

### Costs
drugcostmod <- create_StateVals(drugcost_tbl, n = n_samples, hesim_data = hesim_dat)
medcostmod <- create_StateVals(medcost_tbl, n = n_samples, hesim_data = hesim_dat)
costmods <- list(drugs = drugcostmod,
                 medical = medcostmod)

### Combine
ictstm <- IndivCtstm$new(trans_model = transmod,
                         utility_model = utilitymod,
                         cost_models = costmods)


## Simulate outcomes
head(ictstm$sim_disease()$disprog_)
head(ictstm$sim_stateprobs(t = c(0, 5, 10))$stateprobs_[t == 5])
ictstm$sim_qalys(dr = .03)
ictstm$sim_costs(dr = .03)

### Summarize cost-effectiveness
ce <- ictstm$summarize()
head(ce)
format(summary(ce), pivot_from = "strategy")

Description

Simulate health state transitions in an individual-level continuous time state transition model using parameters from a multi-state model.

Format

An R6::R6Class object.

Super class

hesim::CtstmTrans -> IndivCtstmTrans

Public fields

Methods

Public methods

• IndivCtstmTrans$new()

• IndivCtstmTrans$sim_disease()

• IndivCtstmTrans$sim_stateprobs()

• IndivCtstmTrans$check()

• IndivCtstmTrans$clone()

Method new()

Create a new IndivCtstmTrans object.

Usage

IndivCtstmTrans$new(
  params,
  input_data,
  trans_mat,
  start_state = 1,
  start_age = 38,
  death_state = NULL,
  clock = c("reset", "forward", "mix", "mixt"),
  reset_states = NULL,
  transition_types = NULL
)

Arguments

Returns

A new IndivCtstmTrans object.

Method sim_disease()

Simulate disease progression (i.e., individual trajectories through a multi-state model using an individual patient simulation).

Usage

IndivCtstmTrans$sim_disease(max_t = 100, max_age = 100, progress = NULL)

Arguments

Returns

An object of class disprog.

Method sim_stateprobs()

Simulate health state probabilities from a disprog object.

Usage

IndivCtstmTrans$sim_stateprobs(t, disprog = NULL, ...)

Arguments

Returns

An object of class stateprobs.

Method check()

Input validation for class. Checks that fields are the correct type.

Usage

IndivCtstmTrans$check()

Method clone()

The objects of this class are cloneable with this method.

Usage

IndivCtstmTrans$clone(deep = FALSE)

Arguments

See Also

IndivCtstmTrans objects are conveniently created from either fitted models or parameter objects with create_IndivCtstmTrans(). A complete economic model can be implemented with the IndivCtstm class.

Examples

library("flexsurv")

# Simulation data
strategies <- data.frame(strategy_id = c(1, 2, 3))
patients <- data.frame(patient_id = seq(1, 3),
                       age = c(45, 50, 60),
                       female = c(0, 0, 1))

# Multi-state model with transition specific models
tmat <- rbind(c(NA, 1, 2),
              c(NA, NA, 3),
              c(NA, NA, NA))
fits <- vector(length = max(tmat, na.rm = TRUE), mode = "list")
for (i in 1:length(fits)){
  fits[[i]] <- flexsurvreg(Surv(years, status) ~ 1,
                           data = bosms3[bosms3$trans == i, ],
                           dist = "exp")
}
fits <- flexsurvreg_list(fits)

# Simulation model
hesim_dat <- hesim_data(strategies = strategies,
                        patients = patients)
fits_data <- expand(hesim_dat)
transmod <- create_IndivCtstmTrans(fits, input_data = fits_data,
                                   trans_mat = tmat,
                                   n = 2)
head(transmod$hazard(c(1, 2, 3)))
head(transmod$cumhazard(c(1, 2, 3)))

## Simulate disease progression and state probabilities together
transmod$sim_stateprobs(t = c(0, 5, 10))[t == 5]

## Simulate disease progression and state probabilities separately
disprog <- transmod$sim_disease(max_t = 10)
transmod$sim_stateprobs(t = c(0, 5, 10), disprog = disprog)[t == 5]

Description

An object of class input_mats contains input matrices for simulating a statistical model. Consists of (i) input matrices, X, and (ii) metadata used to index each matrix in X.

Once created, an input_mats object can be converted to a data.table::data.table with as.data.table(), which is a helpful way to check that the object is as expected. The print() method summarizes the object and prints it to the console.

More details are provided under "Details" below.

Usage

input_mats(X, ...)

## S3 method for class 'input_mats'
as.data.table(x, ...)

## S3 method for class 'input_mats'
print(x, ...)

Arguments

Details

input_mats objects are used with params objects to simulate disease models, cost models, and/or utility models. Each column of ⁠$X⁠ contains variables from the params object and a given row corresponds to a combination of the ID variables. An input matrix must always have rows for the treatment strategies (strategy_id) and patients (patient_id); it may optionally have rows for health variables (state_id or transition_id) and time intervals (time_id). The rows must be sorted by prioritizing (i) strategy_id, (ii) patient_id, (iii) the health related ID variable (either state_id or transition_id) and (iv) time_id.

While input_mats objects can be created directly with input_mats(), it is rarely a good idea to do so. They are typically created as the input_data field when creating model objects (e.g., with create_IndivCtstmTrans(), create_CohortDtstmTrans(), and create_PsmCurves()). Internally, these functions create the input matrices using create_input_mats() methods, which ensure that they are in the correct format. Users may also use create_input_mats() methods, but there is not usually a good reason to do so.

as.data.table.input_mats() will convert input matrices into a single data.table() that column binds the ID variables and the unique combinations of variables contained in the elements of ⁠$X⁠. print.input_mats() prints a call to as.data.table() and provides additional information about the ID variables.

See Also

See IndivCtstmTrans() and PsmCurves() for examples in which the input_data field of an instance of a model class is an input_mats object.

Examples

library("data.table")

# Input matrices are typically created as part of model objects
# Let's illustrate with a partitioned survival model (PSM)

## Model setup
strategies <- data.frame(strategy_id = c(1, 2),
                         new_strategy = c(0, 1))
patients <- data.frame(patient_id = seq(1, 3),
                       age = c(45, 47, 60),
                       female = c(1, 0, 0),
                       group = factor(c("Good", "Medium", "Poor")))
hesim_dat <- hesim_data(strategies = strategies,
                        patients = patients)

## Create survival models for PSM
### Parameters
n <- 2
survmod_params <- params_surv_list(
  # Progression free survival (PFS) 
  pfs = params_surv(
    coefs = list(
      rate = data.frame(intercept = rnorm(n, log(1/5), 1),
                        new_strategy = rnorm(n, log(.8), 1))
      ),
    dist = "exp"
  ),
  
  # Overall survival (OS)
  os = params_surv(
    coefs = list(
      rate = data.frame(intercept = rnorm(n, log(1/10), 1))
    ),
    dist = "exp"
  )
)

### Input data
survmod_input_data <- expand(hesim_dat)[, intercept := 1]

### Model object
survmod <- create_PsmCurves(survmod_params, input_data = survmod_input_data)

## Inspect input data
survmod$input_data # Print "input_mats" object to console
as.data.table(survmod$input_data) # Convert "input_mats" object to data.table

Description

Compute the parameters shape1 and shape2 of the beta distribution using method of moments given the mean and standard deviation of the random variable of interest.

Usage

mom_beta(mean, sd)

Arguments

Details

If $\mu$ is the mean and $\sigma$ is the standard deviation of the random variable, then the method of moments estimates of the parameters shape1 = $\alpha > 0$ and shape2 = $\beta > 0$ are:

$\alpha = \mu \left(\frac{\mu(1-\mu)}{\sigma^2}-1 \right)$

and

$\beta = (1 - \mu) \left(\frac{\mu(1-\mu)}{\sigma^2}-1 \right)$

Value

A list containing the parameters shape1 and shape2.

Examples

mom_beta(mean = .8, sd = .1)
# The function is vectorized.
mom_beta(mean = c(.6, .8), sd = c(.08, .1))

Description

Compute the shape and scale (or rate) parameters of the gamma distribution using method of moments for the random variable of interest.

Usage

mom_gamma(mean, sd, scale = TRUE)

Arguments

Details

If $\mu$ is the mean and $\sigma$ is the standard deviation of the random variable, then the method of moments estimates of the parameters shape = $\alpha > 0$ and scale = $\theta > 0$ are:

$\theta = \frac{\sigma^2}{\mu}$

and

$\alpha = \frac{\mu}{\theta}$

The inverse of the scale parameter, $\beta = 1/\theta$, is the rate parameter.

Value

If scale = TRUE, then a list containing the parameters shape and scale; otherwise, if scale = FALSE, then a list containing the parameters shape and rate.

Examples

mom_gamma(mean = 10000, sd = 2000)
# The function is vectorized.
mom_gamma(mean = c(8000, 10000), sd = c(1500, 2000))

Description

Example multi-state data for parameterizing a continuous time state transition model. Costs and utility are also included to facilitate cost-effectiveness analysis.

Usage

mstate3_exdata

Format

A list containing the following elements:

transitions

A data frame containing the times at which patient transitions between health states based on the prothr dataset from the mstate package.

costs

A list of data frames. The first data frame contains summary medical cost estimates and the second data frame contains drug cost data.

utility

A data frame of summary utility estimates.

Transitions data

The data frame has the following columns:

strategy_id

Treatment strategy identification number.

patient_id

Patient identification number.

age

Patient age (in years).

female

1 if a patient is female; 0 if male.

from

Starting state.

to

Receiving state.

trans

Transition number.

Tstart

Starting time.

Tstop

Transition time.

years

Elapsed years between Tstart and Tstop.

status

Status variable; 1=transition, 0=censored.

Cost data

The cost list contains two data frames. The first data frame contains data on the drug costs associated with each treatment strategy.

strategy_id

The treatment strategy identification number.

costs

Annualized drug costs.

The second data frame contains summary data on medical costs by health state, and contains the following columns:

state_id

The health state identification number.

mean

Mean costs.

se

Standard error of medical costs.

Utility data

The data frame has the following columns:

state_id

The health state identification number.

mean

Mean utility

se

Standard error of utility

Description

Combine multinom objects into a list.

Usage

multinom_list(...)

Arguments

Value

An object of class multinom_list.

Examples

library("nnet")
 library("data.table")
 trans_data <- data.table(multinom3_exdata$transitions)
 dat_healthy <- trans_data[state_from == "Healthy"]
 fit_healthy <- multinom(state_to ~ strategy_name + female + age_cat + year_cat, 
                          data = dat_healthy)
 dat_sick <- trans_data[state_from == "Sick"]
 dat_sick$state_to <- droplevels(dat_sick$state_to)
 fit_sick <- multinom(state_to ~ strategy_name + female + age_cat + year_cat, 
                      data = dat_sick)
 fits <- multinom_list(healthy = fit_healthy, sick = fit_sick)
 class(fits)

Description

Example discrete time health state transitions data simulated using multinomial logistic regression. Costs and utility are also included to facilitate cost-effectiveness analysis.

Usage

multinom3_exdata

Format

A list containing the following elements:

transitions

A data frame containing patient transitions between health states at discrete time intervals (i.e., on a yearly basis).

costs

A list of data frames. The first data frame contains drug cost data and the second contains summary medical cost estimates.

utility

A data frame of summary utility estimates.

Transitions data

The data frame has the following columns:

patient_id

Patient identification number.

strategy_id

Treatment strategy identification number.

strategy_name

Treatment strategy name.

age

Patient age (in years).

age_cat

A factor variable with 3 age groups: (i) age less than 40, (ii) age at least 40 and less than 60, and (iii) age at least 60.

female

1 if a patient is female; 0 if male.

year

The year since the start of data collection with the first year equal to 1.

state_from

State making a transition from.

state_to

State making a transition to.

year_cat

Factor variable for year with 3 categories: (i) year 3 and below, (ii) year between 3 and 6, and (iii) year 7 and above.

Cost data

The cost list contains two data frames. The first data frame contains data on the drug costs associated with each treatment strategy.

strategy_id

The treatment strategy identification number.

strategy_name

The treatment strategy name.

costs

Annualized drug costs.

The second data frame contains summary data on medical costs by health state, and contains the following columns:

state_id

The health state identification number.

state_name

The name of the health state.

mean

Mean medical costs.

se

Standard error of medical costs.

Utility data

The data frame has the following columns:

state_id

The health state identification number.

state_name

The name of the health state.

mean

Mean utility

se

Standard error of utility.

Description

Simulated 3-state dataset in oncology with three health states (Stable, Progression, and Death) and three possible transitions (Stable -> Progression, Stable -> Death, and Progression -> Death).

Usage

onc3

Format

A data.table with the following columns:

from

Health state making a transition from.

to

Health state making a transition to.

strategy_name

Standard of care (SOC), new treatment 1 (New 1), or new treatment 2 (New 2).

female

1 if a patient is female; 0 if male.

age

Patient age (in years).

patient_id

Patient identification number.

time_start

Starting time.

time_stop

Stopping time.

status

Status indicator: 1=transition, 0=censored.

transition_id

Integer denoting transition: 1 = Stable -> Progression, 2 = Stable -> Death, 3 = Progression -> Death.

strategy_id

Strategy identification number.

time

Elapsed years between time_start and time_stop.

See Also

onc3p

Examples

head(onc3)

Description

The same dataset as onc3 converted into a panel data format in which health states are recorded at a finite series of times.

Usage

onc3p

Format

A data.table with the following columns:

state

The name of the health state (Stable, Progression, and Death).

strategy_name

Standard of care (SOC), new treatment 1 (New 1), or new treatment 2 (New 2).

female

1 if a patient is female; 0 if male.

age

Patient age (in years).

patient_id

Patient identification number.

time

Time that state was recorded.

strategy_id

Strategy identification number.

state_id

The health state identification number.

See Also

onc3

Examples

head(onc3p)

Description

Objects prefixed by "params_" are lists containing the parameters of a statistical model used for simulation modeling. The parameters are used to simulate outcomes as a function of covariates contained in input matrices (input_mats).

See Also

tparams

Description

Create a list containing the parameters of a fitted linear regression model.

Usage

params_lm(coefs, sigma = 1)

Arguments

Details

Fitted linear models are used to predict values, $y$, as a function of covariates, $x$,

$y = x^T\beta + \epsilon.$

Predicted means are given by $x^T\hat{\beta}$ where $\hat{\beta}$ is the vector of estimated regression coefficients. Random samples are obtained by sampling the error term from a normal distribution, $\epsilon \sim N(0, \hat{\sigma}^2)$.

Value

An object of class params_lm, which is a list containing coefs, sigma, and n_samples. n_samples is equal to the number of rows in coefs. The coefs element is always converted into a matrix.

See Also

This parameter object is useful for modeling health state values when values can vary across patients and/or health states as a function of covariates. In many cases it will, however, be simpler, and more flexible to use a stateval_tbl. For an example use case see the documentation for create_StateVals.lm().

Examples

library("MASS")
n <- 2
params <- params_lm(
  coefs = mvrnorm(n, mu = c(.5,.6),
                  Sigma = matrix(c(.05, .01, .01, .05), nrow = 2)),
  sigma <- rgamma(n, shape = .5, rate = 4)
)
summary(params)
params