This is a replication file for the manuscript. It can be used to produce an HTML file that displays the results without the text with:

# rmarkdown::render("replication.R", output_format = "html_document")

Replication was performed by the authors on macOS Mojave 10.14.6. See the Session Information below for further details.

Settings and dependencies

Package versioning is managed using the renv package. Packages can be installed by cloning the GitHub repository available at https://github.com/hesim-dev/hesim-manuscript and running renv::restore().

library("data.table")
library("flexsurv")
library("ggplot2")
library("hesim")
library("knitr")
library("magrittr")
library("msm")
opts_chunk$set(
  tidy = FALSE, eval = TRUE, cache = TRUE,
  fig.height = 4, fig.width = 6
)
options(prompt = "> ", continue = "+  ", width = 120, useFancyQuotes = TRUE)
ggplot2::theme_set(ggplot2::theme_bw())
set.seed(102)

Illustrative example

Individual continuous time state transition model

Setup model

Target population, treatment strategies, and model structure

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

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

# hesim data
n_patients <- 1000
patients <- data.table(
  patient_id = 1:n_patients,
  age = rnorm(n_patients, mean = 45, sd = 7),
  female = rbinom(n_patients, size = 1, prob = .51)
)

states <- data.table(
  state_id = c(1, 2),
  state_name = c("Stable", "Progression") # Non-death health states
)

strategies <- data.table(
  strategy_id = 1:3,
  strategy_name = c("SOC", "New 1", "New 2"),
  soc = c(1, 0, 0),
  new1 = c(0, 1, 0),
  new2 = c(0, 0, 1)
)

hesim_dat <- hesim_data(
  strategies = strategies,
  patients = patients,
  states = states
)
print(hesim_dat)

## $strategies
##    strategy_id strategy_name soc new1 new2
## 1:           1           SOC   1    0    0
## 2:           2         New 1   0    1    0
## 3:           3         New 2   0    0    1
## 
## $patients
##       patient_id      age female
##    1:          1 46.26366      1
##    2:          2 50.49314      1
##    3:          3 35.52785      1
##    4:          4 58.88309      0
##    5:          5 53.66930      0
##   ---                           
##  996:        996 54.60555      0
##  997:        997 50.10654      0
##  998:        998 45.02035      0
##  999:        999 45.60573      1
## 1000:       1000 39.45426      0
## 
## $states
##    state_id  state_name
## 1:        1      Stable
## 2:        2 Progression
## 
## attr(,"class")
## [1] "hesim_data"

Labels for ID variables

labs_indiv <- get_labels(hesim_dat)
print(labs_indiv)

## $strategy_id
##   SOC New 1 New 2 
##     1     2     3 
## 
## $state_id
##      Stable Progression       Death 
##           1           2           3

Parameterization

Transition model

# Estimation data
onc3[patient_id %in% c(1, 2)] %>%
  kable()

from	to	strategy_name	age	patient_id	time_start	time_stop	status	transition_id	strategy_id	time
Stable	Progression	New 2	59.85813	1	0.000000	2.420226	1	1	3	2.420226
Stable	Death	New 2	59.85813	1	0.000000	2.420226	0	2	3	2.420226
Progression	Death	New 2	59.85813	1	2.420226	14.620259	1	3	3	12.200032
Stable	Progression	New 2	62.57282	2	0.000000	7.497464	0	1	3	7.497464
Stable	Death	New 2	62.57282	2	0.000000	7.497464	0	2	3	7.497464

# Fit multi-state model
n_trans <- max(tmat, na.rm = TRUE) 
wei_fits <- vector(length = n_trans, mode = "list")
f <- as.formula(Surv(time, status) ~ factor(strategy_name) + female + age)
for (i in 1:length(wei_fits)){
  if (i == 3) f <- update(f, .~.-factor(strategy_name)) 
  wei_fits[[i]] <- flexsurvreg(f, data = onc3, 
                               subset = (transition_id == i),
                               dist = "weibull")
}
wei_fits <- flexsurvreg_list(wei_fits)

Utility model

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

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

Cost models

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

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

# Drug costs
n_times <- 2
n_states <- nrow(states)
n_strategies <- nrow(strategies)
n_i <- n_states * n_times
drugcost_dt <-  data.table(
    strategy_id = rep(strategies$strategy_id, each = n_i),
    state_id = rep(rep(states$state_id, each = n_times), n_strategies),
    time_start = rep(c(0, 3/12), n_states * n_strategies),
    est = c(2000, 2000, 1500, 1200,
           12000, 12000, 1500, 1200,
           15000, 15000, 1500, 1200)
) 

drugcost_tbl <- stateval_tbl(
  drugcost_dt,
  dist = "fixed")
print(drugcost_tbl)

##     strategy_id state_id time_id time_start time_stop   est
##  1:           1        1       1       0.00      0.25  2000
##  2:           1        1       2       0.25       Inf  2000
##  3:           1        2       1       0.00      0.25  1500
##  4:           1        2       2       0.25       Inf  1200
##  5:           2        1       1       0.00      0.25 12000
##  6:           2        1       2       0.25       Inf 12000
##  7:           2        2       1       0.00      0.25  1500
##  8:           2        2       2       0.25       Inf  1200
##  9:           3        1       1       0.00      0.25 15000
## 10:           3        1       2       0.25       Inf 15000
## 11:           3        2       1       0.00      0.25  1500
## 12:           3        2       2       0.25       Inf  1200

Simulation

Constructing the economic model

n_samples <- 1000 # Number of samples for PSA

Transition model

transmod_data <- expand(hesim_dat,
                        by = c("strategies", "patients"))
head(transmod_data)

##    strategy_id patient_id strategy_name soc new1 new2      age female
## 1:           1          1           SOC   1    0    0 46.26366      1
## 2:           1          2           SOC   1    0    0 50.49314      1
## 3:           1          3           SOC   1    0    0 35.52785      1
## 4:           1          4           SOC   1    0    0 58.88309      0
## 5:           1          5           SOC   1    0    0 53.66930      0
## 6:           1          6           SOC   1    0    0 53.40432      1

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

Utility and cost models

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

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

Economic model

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

Simulating outcomes

Individual trajectories through model

ictstm$sim_disease(max_age = 100)
head(ictstm$disprog_)

##    sample strategy_id patient_id grp_id from to final time_start time_stop
## 1:      1           1          1      1    1  3     1   0.000000 11.183651
## 2:      1           1          2      1    1  2     0   0.000000  8.058198
## 3:      1           1          2      1    2  3     1   8.058198 12.794488
## 4:      1           1          3      1    1  2     0   0.000000 15.250551
## 5:      1           1          3      1    2  3     1  15.250551 23.791357
## 6:      1           1          4      1    1  3     1   0.000000  6.047286

State probabilities

ictstm$sim_stateprobs(t = seq(0, 30 , 1/12))
head(ictstm$stateprobs_)

##    sample strategy_id grp_id state_id          t  prob
## 1:      1           1      1        1 0.00000000 1.000
## 2:      1           1      1        1 0.08333333 1.000
## 3:      1           1      1        1 0.16666667 1.000
## 4:      1           1      1        1 0.25000000 0.999
## 5:      1           1      1        1 0.33333333 0.998
## 6:      1           1      1        1 0.41666667 0.996

autoplot(ictstm$stateprobs_, labels = labs_indiv,
         ci = FALSE)

QALYs and costs

# QALYs
ictstm$sim_qalys(dr = c(0,.03))
head(ictstm$qalys_)

##    sample strategy_id grp_id state_id dr    qalys      lys
## 1:      1           1      1        1  0 4.933188 6.374648
## 2:      1           1      1        2  0 3.233055 5.255885
## 3:      1           2      1        1  0 5.604984 7.242740
## 4:      1           2      1        2  0 2.683823 4.363013
## 5:      1           3      1        1  0 6.289275 8.126979
## 6:      1           3      1        2  0 2.915707 4.739981

# Costs
ictstm$sim_costs(dr = .03)
head(ictstm$costs_)

##    sample strategy_id grp_id state_id   dr category      costs
## 1:      1           1      1        1 0.03     Drug  11352.373
## 2:      1           1      1        2 0.03     Drug   4691.177
## 3:      1           2      1        1 0.03     Drug  76292.008
## 4:      1           2      1        2 0.03     Drug   3819.849
## 5:      1           3      1        1 0.03     Drug 105209.810
## 6:      1           3      1        2 0.03     Drug   4042.658

Summarize QALYs and costs

ce_sim_ictstm <- ictstm$summarize()
summary(ce_sim_ictstm, labels = labs_indiv) %>%
  format() %>%
  kable()

Discount rate	Outcome	SOC	New 1	New 2
0.00	QALYs	7.87 (7.12, 8.60)	8.59 (7.81, 9.40)	9.20 (8.42, 10.06)
0.03	QALYs	6.58 (6.03, 7.14)	7.10 (6.52, 7.69)	7.52 (6.96, 8.11)
0.03	Costs: Drug	15,352 (14,431, 16,234)	78,465 (73,854, 83,200)	105,306 (99,228, 111,686)
0.03	Costs: Medical	46,784 (4,518, 144,256)	46,222 (4,543, 137,303)	46,840 (4,746, 143,093)
0.03	Costs: total	62,136 (19,888, 159,037)	124,687 (82,038, 218,484)	152,146 (109,625, 245,279)

Partitioned survival model

Setup model

Target population for cohort model

hesim_dat$patients <- data.table(
  patient_id = 1:4,
  grp_id = c(1, 1, 2, 2),
  patient_wt = rep(1/4, 4),
  age = c(55, 65, 55, 65),
  female = c(1, 1, 0, 0),
  grp_name = rep(c("Female", "Male"), 
                 each = 2)
)

Labels for cohort model

labs_cohort <- get_labels(hesim_dat)

Parameterization

Survival models

Parameters are estimated using both individual patient data and summary data: a model is fit for SOC and relative treatment effects are obtained from the literature.

# Estimation data
onc_pfs_os <- as_pfs_os(onc3, patient_vars = c("patient_id", "female", "age", 
                                               "strategy_name"))
onc_pfs_os[patient_id %in% c(1, 2)]

##    patient_id female      age strategy_name pfs_time pfs_status os_status   os_time
## 1:          1      0 59.85813         New 2 2.420226          1         1 14.620258
## 2:          2      0 62.57282         New 2 7.497464          0         0  7.497464

# Parameters
## Fit models for SOC
onc_pfs_os_soc <- onc_pfs_os[strategy_name == "SOC"]
fit_pfs_soc_wei <- flexsurv::flexsurvreg(
  Surv(pfs_time, pfs_status) ~ female,
  data = onc_pfs_os_soc,
  dist = "weibull")

fit_os_soc_wei <- flexsurvreg(
  Surv(os_time, os_status) ~  female,
  data = onc_pfs_os_soc,
  dist = "weibull")

psmfit_soc_wei <- partsurvfit(
  flexsurvreg_list(pfs = fit_pfs_soc_wei, os = fit_os_soc_wei),
  data = onc_pfs_os_soc
)

## Obtain relative treatment effects
## Fit using the "onc_pfs_os" dataset but assume the estimates are
## obtained from the literature for illustrative purposes
### Create dummies
onc_pfs_os[, new1 := ifelse(strategy_name == "New 1", 1, 0)]
onc_pfs_os[, new2 := ifelse(strategy_name == "New 2", 1, 0)]

### Fit models
fit_pfs_wei <- flexsurv::flexsurvreg(
  Surv(pfs_time, pfs_status) ~ female + new1 + new2,
  data = onc_pfs_os,
  dist = "weibull")

fit_os_wei <- flexsurvreg(
  Surv(os_time, os_status) ~ female + new1 + new2,
  data = onc_pfs_os,
  dist = "weibull")

### Extract coefficients (assumed taken from literature)
params_pfs_wei <- create_params(fit_pfs_wei, n = n_samples)
params_os_wei <- create_params(fit_os_wei, n = n_samples)
coef_pfs_wei <- params_pfs_wei$coefs$scale[, c("new1", "new2")]
coef_os_wei <- params_os_wei$coefs$scale[, c("new1", "new2")]
head(coef_pfs_wei)

##           new1      new2
## [1,] 0.1814189 0.2654431
## [2,] 0.1390351 0.2907570
## [3,] 0.1298460 0.2307868
## [4,] 0.1729702 0.3143110
## [5,] 0.1988663 0.3179721
## [6,] 0.1305311 0.2261630

Costs

# Update drug costs to work with cohort model
drugcost_tbl <- stateval_tbl(
  drugcost_dt[time_start ==  0][, time_start := NULL],
  dist = "fixed")
print(drugcost_tbl)

##    strategy_id state_id   est
## 1:           1        1  2000
## 2:           1        2  1500
## 3:           2        1 12000
## 4:           2        2  1500
## 5:           3        1 15000
## 6:           3        2  1500

Simulation

Constructing the economic model

Survival models

# Input data
survmods_data <- expand(hesim_dat, by = c("strategies", "patients"))
head(survmods_data)

##    strategy_id patient_id strategy_name soc new1 new2 grp_id patient_wt age female grp_name
## 1:           1          1           SOC   1    0    0      1       0.25  55      1   Female
## 2:           1          2           SOC   1    0    0      1       0.25  65      1   Female
## 3:           1          3           SOC   1    0    0      2       0.25  55      0     Male
## 4:           1          4           SOC   1    0    0      2       0.25  65      0     Male
## 5:           2          1         New 1   0    1    0      1       0.25  55      1   Female
## 6:           2          2         New 1   0    1    0      1       0.25  65      1   Female

# Parameters
## Parameters for SOC. We obtain a warning because the model did not
## converge during one of the bootstrap samples; this is allowed because
## max_errors = 5 >= 1.
survmods_params <- create_params(psmfit_soc_wei, n = n_samples, 
                                 uncertainty = "bootstrap", 
                                 max_errors = 5, silent = TRUE)

## Warning in survreg.fit(X, Y, weights, offset, init = init, controlvals = control, : Ran out of iterations and did not
## converge

## Warning in bootstrap.partsurvfit(object, B = n, ...): There were 1 errors when fitting the statistical model during the
## bootstrap procedure.

## "Manually" add relative treatment effects
survmods_params$pfs$coefs$scale <- cbind(survmods_params$pfs$coefs$scale,
                                         coef_pfs_wei)
survmods_params$os$coefs$scale <- cbind(survmods_params$os$coefs$scale,
                                        coef_os_wei)
survmods_params$pfs$coefs$scale[1:2, ]

##      (Intercept)     female      new1      new2
## [1,]    1.784666 -0.1364781 0.1814189 0.2654431
## [2,]    1.776743 -0.1173777 0.1390351 0.2907570

# Combine input data and parameters to create the survival models
survmods_data[, ("(Intercept)") := 1]
survmods <- create_PsmCurves(survmods_params, 
                             input_data = survmods_data)

Utility and cost models

utilitymod <- create_StateVals(utility_tbl, n = n_samples, 
                               hesim_data = hesim_dat)
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(Drug = drugcostmod, Medical = medcostmod)

Economic model

psm <- Psm$new(survival_models = survmods,
               utility_model = utilitymod,
               cost_models = costmods)

Simulating outcomes

Survival

times <- seq(0, 30, by = .1)
psm$sim_survival(t = times)
head(psm$survival_)

##    sample strategy_id patient_id grp_id patient_wt curve   t  survival
## 1:      1           1          1      1       0.25     1 0.0 1.0000000
## 2:      1           1          1      1       0.25     1 0.1 0.9998540
## 3:      1           1          1      1       0.25     1 0.2 0.9993126
## 4:      1           1          1      1       0.25     1 0.3 0.9982992
## 5:      1           1          1      1       0.25     1 0.4 0.9967670
## 6:      1           1          1      1       0.25     1 0.5 0.9946815

autoplot(psm$survival_, 
         labels = c(labs_cohort,
                    list(curve = c("PFS" = 1, "OS" = 2))
                    ),
         ci = TRUE, ci_style = "ribbon")

State probabilities

psm$sim_stateprobs()
psm$stateprobs_[sample == 1 & patient_id == 1 & state_id == 2 & t == 12]

##    sample strategy_id patient_id grp_id patient_wt state_id  t      prob
## 1:      1           1          1      1       0.25        2 12 0.3215069
## 2:      1           2          1      1       0.25        2 12 0.3934037
## 3:      1           3          1      1       0.25        2 12 0.4148810

QALYs and costs

psm$sim_qalys(dr = .03) 
psm$sim_costs(dr = .03)

Summarize QALY and costs

ce_sim_psm <- psm$summarize(by_grp = TRUE)
summary(ce_sim_psm, labels = labs_cohort) %>%
  format() %>%
  kable()

Group	Discount rate	Outcome	SOC	New 1	New 2
Female	0.03	QALYs	6.10 (5.55, 6.63)	6.61 (6.00, 7.22)	6.96 (6.36, 7.60)
Female	0.03	Costs: Drug	15,225 (14,629, 15,900)	66,473 (62,592, 70,188)	88,466 (83,561, 93,548)
Female	0.03	Costs: Medical	52,296 (4,435, 180,774)	52,717 (4,871, 171,715)	53,138 (5,156, 169,727)
Female	0.03	Costs: total	67,521 (19,657, 195,939)	119,189 (70,613, 239,546)	141,604 (92,445, 260,157)
Male	0.03	QALYs	6.51 (5.99, 7.08)	7.05 (6.46, 7.70)	7.42 (6.81, 8.07)
Male	0.03	Costs: Drug	16,247 (15,553, 17,004)	74,027 (69,772, 78,081)	98,691 (93,156, 104,387)
Male	0.03	Costs: Medical	52,677 (4,767, 171,325)	52,709 (5,188, 165,850)	52,874 (5,417, 162,113)
Male	0.03	Costs: total	68,924 (20,902, 187,307)	126,736 (78,810, 238,080)	151,565 (103,556, 262,961)

Cohort discrete time state transition model

Parameterization

Transition model

We will first obtain the underlying parameter estimates. A model is fit for SOC and relative treatment effects are obtained from the literature.

# Estimation data
head(onc3p)

##          state strategy_name female      age patient_id      time strategy_id state_id
## 1:      Stable         New 2      0 59.85813          1  0.000000           3        1
## 2: Progression         New 2      0 59.85813          1  2.420226           3        2
## 3:       Death         New 2      0 59.85813          1 14.620258           3        3
## 4:      Stable         New 2      0 62.57282          2  0.000000           3        1
## 5:      Stable         New 2      0 62.57282          2  7.497464           3        1
## 6:      Stable           SOC      1 61.44379          3  0.000000           1        1

# Fit multi-state model for SOC with msm
qinit <- matrix(0, nrow = 3, ncol = 3)
qinit[1, 2] <-  0.28; qinit[1, 3] <-  0.013; qinit[2, 3] <-  0.10
msm_fit <- msm(state_id ~ time, subject = patient_id, 
               data = onc3p[strategy_name == "SOC"],
               exacttimes = FALSE,
               covariates = ~ age + female,
               qmatrix = qinit, gen.inits = FALSE)

# Compare model fit with msm using exact times to an exponential model
## msm with exact transition times
tmat2 <- tmat; tmat2[is.na(tmat2)] <- 0
msm_efit <- msm(state_id ~ time, subject = patient_id, 
                data = onc3p[strategy_name == "SOC"],
                qmatrix = tmat2, exacttimes = TRUE)

## flexsurv with exponential model
exp_fits <- vector(length = n_trans, mode = "list")
for (i in 1:length(wei_fits)){
  exp_fits[[i]] <- flexsurvreg(Surv(time, status) ~ 1, data = onc3, 
                               subset = (transition_id == i & 
                                         strategy_name == "SOC"),
                               dist = "exponential")
}

## Compare the coefficients
msm_efit

## 
## Call:
## msm(formula = state_id ~ time, subject = patient_id, data = onc3p[strategy_name ==     "SOC"], qmatrix = tmat2, exacttimes = TRUE)
## 
## Maximum likelihood estimates
## 
## Transition intensities
##                           Baseline                   
## Stable - Stable           -0.17427 (-0.1880,-0.16155)
## Stable - Progression       0.12999 ( 0.1191, 0.14191)
## Stable - Dead              0.04428 ( 0.0381, 0.05146)
## Progression - Progression -0.09088 (-0.1027,-0.08040)
## Progression - Dead         0.09088 ( 0.0804, 0.10272)
## 
## -2 * log-likelihood:  6173.881 
## [Note, to obtain old print format, use "printold.msm"]

exp_fits[[1]]$res

##            est      L95%      U95%          se
## rate 0.1299939 0.1190743 0.1419149 0.005819325

exp_fits[[2]]$res

##             est       L95%       U95%         se
## rate 0.04428652 0.03810548 0.05147018 0.00339662

exp_fits[[3]]$res

##            est       L95%      U95%          se
## rate 0.0908797 0.08040199 0.1027228 0.005679981

# Relative treatment effects (i.e., relative risk)
# Assume estimated from external source and simulate distribution
# for PSA
params_rr <- list(
  lrr_12_est = c(soc = log(1), new1 = log(.80), new2 = log(.71)),
  lrr_12_se = c(soc = 0, new1 = .03, new2 = .04),
  lrr_13_est = c(soc = log(1), new1 = log(.90), new2 = log(.85)),
  lrr_13_se = c(soc = 0, new1 = .02, new2 = .03)
)  

params_rr_def <- define_rng({
  list(
    rr_12 = lognormal_rng(lrr_12_est, lrr_12_se),
    rr_13 = lognormal_rng(lrr_13_est, lrr_13_se)
  )}, n = n_samples)
params_rr_rng <- eval_rng(params_rr_def, params_rr)

We now use the underlying parameter estimates and input data to generate the “transformed” parameters.

# Start by predicting transition intensity matrices for each patient
# with SOC. Use a multivariate normal distribution to simulate samples
# for the PSA
transmod_data <- survmods_data
qmat_soc <- qmatrix(msm_fit, 
                    newdata = transmod_data[strategy_name == "SOC"],
                    uncertainty = "normal", n = n_samples)
dim(qmat_soc)

## [1]    3    3 4000

qmat_soc[,, 1]

##            [,1]       [,2]       [,3]
## [1,] -0.3533989  0.3371267 0.01627226
## [2,]  0.0000000 -0.1190753 0.11907530
## [3,]  0.0000000  0.0000000 0.00000000

# Take the matrix exponential to convert the transition intensity matrices to
# transition probability matrices
cycle_len <- 1/4
pmat_soc <- expmat(qmat_soc, t = cycle_len)
pmat_soc[,, 1]

##           [,1]       [,2]       [,3]
## [1,] 0.9154407 0.07945955 0.00509979
## [2,] 0.0000000 0.97066990 0.02933010
## [3,] 0.0000000 0.00000000 1.00000000

# Align the transition matrices with the ID variables
tpmat_id <- tpmatrix_id(transmod_data, n_samples)
head(tpmat_id)

##    sample strategy_id patient_id grp_id patient_wt
## 1:      1           1          1      1       0.25
## 2:      1           1          2      1       0.25
## 3:      1           1          3      2       0.25
## 4:      1           1          4      2       0.25
## 5:      1           2          1      1       0.25
## 6:      1           2          2      1       0.25

# Compute relative risks for each row in the input data and each parameter 
# sample for the PSA
xbeta <- function(x, beta) c(x %*% t(beta))
x_rr <- as.matrix(transmod_data[, .(soc, new1, new2)])
rr <- cbind(xbeta(x_rr, params_rr_rng$rr_12),
            xbeta(x_rr, params_rr_rng$rr_13))
head(rr) # Corresponds to 1st 6 rows in tpmat_id

##          [,1]      [,2]
## [1,] 1.000000 1.0000000
## [2,] 1.000000 1.0000000
## [3,] 1.000000 1.0000000
## [4,] 1.000000 1.0000000
## [5,] 0.833283 0.9039595
## [6,] 0.833283 0.9039595

# Apply the relative risks to SOC to get transition probabilities for each
# parameter sample, treatment strategy, and representative patient
pmat <- apply_rr(pmat_soc, rr = rr,
                 index = list(c(1, 2), c(1, 3)))
tprobs <- tparams_transprobs(pmat, tpmat_id)

Simulation

Constructing the economic model

Transition model

transmod <- CohortDtstmTrans$new(params = tprobs, 
                                 cycle_length = cycle_len)

Economic model: use the utility and cost models from the PSM

cdtstm <- CohortDtstm$new(trans_model = transmod,
                          utility_model = psm$utility_model,
                          cost_models = psm$cost_models)

Simulating outcomes

State probabilities

cdtstm$sim_stateprobs(n_cycles = 30/cycle_len)

QALYs and costs

cdtstm$sim_qalys(dr = .03)
cdtstm$sim_costs(dr = .03)

Summarize QALYs and costs

ce_sim_cdtstm <- cdtstm$summarize()
summary(ce_sim_cdtstm, labels = labs_cohort) %>%
  format() %>%
  kable()

Discount rate	Outcome	SOC	New 1	New 2
0.03	QALYs	5.56 (4.78, 6.26)	5.92 (5.10, 6.67)	6.14 (5.27, 6.84)
0.03	Costs: Drug	13,862 (12,563, 14,711)	46,196 (41,201, 50,104)	60,714 (53,418, 66,584)
0.03	Costs: Medical	62,515 (4,186, 217,838)	62,011 (4,351, 223,642)	61,636 (4,671, 213,904)
0.03	Costs: total	76,378 (18,215, 231,236)	108,207 (50,113, 267,630)	122,350 (64,868, 274,308)

Cost-effectiveness analysis

# "ce" object to use for analysis
ce_sim_ictstm

## $costs
##       category   dr sample strategy_id     costs grp_id
##    1:     Drug 0.03      1           1  16043.55      1
##    2:     Drug 0.03      1           2  80111.86      1
##    3:     Drug 0.03      1           3 109252.47      1
##    4:     Drug 0.03      2           1  14645.16      1
##    5:     Drug 0.03      2           2  77962.49      1
##   ---                                                  
## 8996:    total 0.03    999           2  95024.16      1
## 8997:    total 0.03    999           3 122357.69      1
## 8998:    total 0.03   1000           1  52638.13      1
## 8999:    total 0.03   1000           2 119410.47      1
## 9000:    total 0.03   1000           3 151164.45      1
## 
## $qalys
##         dr sample strategy_id    qalys grp_id
##    1: 0.00      1           1 8.166243      1
##    2: 0.00      1           2 8.288806      1
##    3: 0.00      1           3 9.204982      1
##    4: 0.00      2           1 7.440970      1
##    5: 0.00      2           2 8.414148      1
##   ---                                        
## 5996: 0.03    999           2 7.134710      1
## 5997: 0.03    999           3 7.399217      1
## 5998: 0.03   1000           1 6.435002      1
## 5999: 0.03   1000           2 7.103850      1
## 6000: 0.03   1000           3 7.584523      1
## 
## attr(,"class")
## [1] "ce"

Perform cost-effectiveness analysis

wtp <- seq(0, 250000, 500) # Willingness to pay per QALY
cea_ictstm <- cea(ce_sim_ictstm, dr_costs = .03, dr_qalys = .03, k = wtp)
cea_pw_ictstm <- cea_pw(ce_sim_ictstm, comparator = 1,
                        dr_qalys = .03, dr_costs = .03,
                        k = wtp)

Incremental cost-effectiveness ratio

icer(cea_pw_ictstm, k = 100000, labels = labs_indiv) %>%
  format() %>%
  kable()

Outcome	New 1	New 2
Incremental QALYs	0.52 (0.16, 0.90)	0.93 (0.58, 1.27)
Incremental costs	62,551 (51,685, 70,389)	90,010 (77,342, 100,581)
Incremental NMB	-11,003 (-42,132, 21,247)	3,223 (-28,582, 33,456)
ICER	121,345	96,543

Representing decision uncertainty

Cost-effectiveness plane

plot_ceplane(cea_pw_ictstm, labels = labs_indiv)

Cost-effectiveness acceptability curve

Simultaneous comparison

plot_ceac(cea_ictstm, labels = labs_indiv)

Pairwise comparison

plot_ceac(cea_pw_ictstm, labels = labs_indiv)

Cost-effectiveness acceptability frontier

plot_ceaf(cea_ictstm, labels = labs_indiv)

Expected value of perfect information

plot_evpi(cea_ictstm, labels = labs_indiv)

Session information

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
## [1] msm_1.6.8         magrittr_2.0.1    knitr_1.31        hesim_0.5.0       ggplot2_3.3.3     flexsurv_1.1.1   
## [7] survival_3.2-7    data.table_1.14.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6         highr_0.8          pillar_1.4.7       compiler_4.0.3     RColorBrewer_1.1-2 tools_4.0.3       
##  [7] digest_0.6.27      evaluate_0.14      lifecycle_1.0.0    tibble_3.0.6       gtable_0.3.0       lattice_0.20-41   
## [13] pkgconfig_2.0.3    rlang_0.4.10       Matrix_1.3-2       yaml_2.2.1         mvtnorm_1.1-1      expm_0.999-6      
## [19] xfun_0.21          withr_2.4.1        stringr_1.4.0      dplyr_1.0.4        generics_0.1.0     vctrs_0.3.6       
## [25] grid_4.0.3         tidyselect_1.1.0   deSolve_1.28       mstate_0.3.1       glue_1.4.2         R6_2.5.0          
## [31] rmarkdown_2.7      farver_2.0.3       tidyr_1.1.2        purrr_0.3.4        MASS_7.3-53.1      codetools_0.2-18  
## [37] htmltools_0.5.1.1  scales_1.1.1       ellipsis_0.3.1     splines_4.0.3      colorspace_2.0-0   renv_0.12.5       
## [43] labeling_0.4.2     quadprog_1.5-8     stringi_1.5.3      muhaz_1.2.6.1      munsell_0.5.0      crayon_1.4.1

hesim: Health Economic Simulation Modeling and Decision Analysis

Devin Incerti, Jeroen P Jansen

2021-02-27

Settings and dependencies

Illustrative example

Individual continuous time state transition model

Setup model

Target population, treatment strategies, and model structure

Labels for ID variables

Parameterization

Transition model

Utility model

Cost models

Simulation

Constructing the economic model

Transition model

Utility and cost models

Economic model

Simulating outcomes

Individual trajectories through model

State probabilities

QALYs and costs

Summarize QALYs and costs

Partitioned survival model

Setup model

Target population for cohort model

Labels for cohort model

Parameterization

Survival models

Costs

Simulation

Constructing the economic model

Survival models

Utility and cost models

Economic model

Simulating outcomes

Survival

State probabilities

QALYs and costs

Summarize QALY and costs

Cohort discrete time state transition model

Parameterization

Transition model

Simulation

Constructing the economic model

Transition model

Economic model: use the utility and cost models from the PSM

Simulating outcomes

State probabilities

QALYs and costs

Summarize QALYs and costs

Cost-effectiveness analysis

Perform cost-effectiveness analysis

Incremental cost-effectiveness ratio

Representing decision uncertainty

Cost-effectiveness plane

Cost-effectiveness acceptability curve

Simultaneous comparison

Pairwise comparison

Cost-effectiveness acceptability frontier

Expected value of perfect information

Session information