hesim-dev / hesim

Health economic simulation modeling and decision analysis
https://hesim-dev.github.io/hesim/
62 stars 15 forks source link

Building a model with known transition probabilities for each strategy #103

Closed studentdoctorcoder closed 1 year ago

studentdoctorcoder commented 1 year ago

Hi! I'm trying to use the package to create a cohort with two transition probability matrices using the create_CohortDtstm. I've written the entire code below since I'm not entirely sure what is causing the error. There are 2 transition matrices because I want to test two different treatments, and they each carry with them different probabilities of moving between the 6 states in this model (both treatment groups experience the same 6 states).

library("hesim") library("data.table")

OUTPUT 1: BASE CASE

strategies <- data.table(strategy_id = 1:2, strategy_name = c("SA", "SB")) patients <- data.table(patient_id = 1) states <- data.table(state_id = 1:6, state_name = c("A", "B", "C", "D", "E", "F")) hesim_dat <- hesim_data(strategies = strategies, patients = patients, states = states)

trans_matrix1 <- matrix (c([assume I have a TP here]), byrow = byRow, nrow = 6) colnames(trans_matrix1 <- rownames(trans_matrix1) <- c("A", "B", "C", "D", "E", "F") trans_matrix2 <- matrix (c([assume I have a different TP here]), byrow = byRow, nrow = 6) colnames(trans_matrix2 <- rownames(trans_matrix1) <- c("A", "B", "C", "D", "E", "F")

Base Costs

*Assume I set up all cost variables like: cost1 <- 100 cost2 <- 200

Base Utility

*Assume I set up the utility of each of the 6 states, with each greater than 0 or less than 1 uA <- 0.5

Cost and Utility Maps

cost_map <- c(C=#, C2=#,R=#, RNS=#,SD=#,BM=#) ## costs utility_map <- c(C=#, C2=#, R=#, RNS=#, SD=#, BM=#) #utility

Parameters

params <- list( alpha_matrix1 = trans_matrix1 alpha_matrix2 = trans_matrix2 cost_map_mean = cost_map, utility_map_mean = utility_map, cost1_mean <- cost1 [assume I do the same with the rest of the cost variables] utility1_mean <- utility1 [assume I do the same with the rest of the utility variables] )

Random Number Generation

rng_def <- define_rng({ list ( p_matrix1 = gamma_rng(alpha_matrix1, sd = alpha_matrix1), p_matrix2 = gamma_rng(alpha_matrix2, sd = alpha_matrix2), r_utility_map = gamma_rng(utility_map_mean, sd = utility_map_mean), r_cost1 <- gamma_rng(cost1_mean, sd = cost1_mean), r_utility1 <- gamma_rng(utility1_mean, sd = utility1_mean), ) }, n = 1000)

Transformed Parameters

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

k <- 0 j <- 0 tparams_def <- define_tparams({ if (strategy_name == "SA") { k <- ifelse(state_name == "R", k <- k + 1, 0) list ( tpmatrix = p_matrix1, utility = p_utility_map, costs = list ( all costs listed here) ) ) } if (strategy_name == "SB") { j <- ifelse(state_name == "", j <- j + 1, 0) list ( tpmatrix = p_matrix2, utility = p_utility_map, costs = list ( all costs listed here) ) ) } }, times = c(1, Inf))

mod_def <- define_model(tparams_def = tparams_def, rng_def = rng_def, params = params)

econmod <- create_CohortDtstm(mod_def, input_data)


That code produces this error: Error in switch(class, matrix = as.matrix(z), data.table = as.data.table(z), : EXPR must be a length 1 vector

Does anyone know what's going on? I've looked all over Stack Exchange and cannot find a solution, since I am using this package and its create_CohortDtstm command.

dincerti commented 1 year ago

Hi @studentdoctorcoder, please include a reproducible example that generates the error.

studentdoctorcoder commented 1 year ago

library("hesim") library("data.table")

OUTPUT 1: BASE CASE

strategies <- data.table(strategy_id = 1:2, strategy_name = c("SA", "SB")) patients <- data.table(patient_id = 1) states <- data.table(state_id = 1:6, state_name = c("A", "B", "C", "D", "E", "F")) hesim_dat <- hesim_data(strategies = strategies, patients = patients, states = states) trans_SA <- matrix(c(0.92, 0, 0.04, 0, 0.04, 0,

SBHealthStates <- c("A", "B", "C", "D", "E", "F") byRow <- TRUE trans_SB <- matrix(c(0.92, 0, 0.04, 0, 0.04, 0, 0, 0.92, 0, 0.04, 0.04, 0,

Base Costs

cost1 <- 246.40 cost2 <- 6109.06 SA <- 17490.45 SB <- 12655.02 cost3 <- 397.9 cost4 <- 10225.13 cost5 <- 102.43 cost6 <- 3624.15 cost7 <- 4.90

Base Utility

uc <- 0.83 uncf3 <- 0.1025 uncf4 <- 0.07 ur <- 0.73 ud <- 0.2

Cost and Utility Maps

cost_map <- c(C=0, C2=0,R=0, RNS=0,SD=0,BM=0) ## costs utility_map <- c(C=uc, C2=uc, R=ur, RNS=ur, SD=ud, BM=ud) #utility

Parameters

params <- list( alpha_SA = trans_SA, alpha_SB = trans_SB, cost_map_mean = cost_map, utility_map_mean = utility_map, cost1_mean <- cost1, cost2_mean <- cost2, SA_mean <- SA, SB_mean <- SB, cost3_mean <- cost3, cost4_mean <- cost4, cost5_mean <- cost5, cost6_mean <- cost6, cost7_mean <- cost7, uncf3_mean <- uncf3, uncf4_mean <- uncf4 )

Random Number Generation

rng_def <- define_rng({ list ( p_SA = gamma_rng(alpha_SA, sd = alpha_SA), p_SB = gamma_rng(alpha_SB, sd = alpha_SB), r_utility_map = gamma_rng(utility_map_mean, sd = utility_map_mean), r_cost1 <- gamma_rng(cost1_mean, sd = cost1_mean), r_cost2 <- gamma_rng(cost2_mean, sd = cost2_mean), r_SA <- gamma_rng(SA_mean, sd = SA_mean), r_SB <- gamma_rng(SB_mean, sd = SB_mean), r_cost3 <- gamma_rng(cost3_mean, sd = cost3_mean), r_cost4 <- gamma_rng(cost4_mean, sd = cost4_mean), r_cost5 <- gamma_rng(cost5_mean, sd = cost5_mean), r_cost6 <- gamma_rng(cost6_mean, sd = cost6_mean), r_cost7 <- gamma_rng(cost7_mean, sd = cost7_mean), r_uncf3 <- gamma_rng(uncf3_mean, sd = uncf3_mean), r_uncf4 <- gamma_rng(uncf4_mean, sd = uncf4_mean) ) }, n = 1000)

Transformed Parameters

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

k <- 0 j <- 0 tparams_def <- define_tparams({ if (strategy_name == "SA") { k <- ifelse(state_name == "C", k <- k + 1, 0) list ( tpmatrix = p_SA, utility = p_utility_map, costs = list ( treatment = ifelse(state_name == "A" | time < 2, p_SA, 0), treatmentSB = ifelse(state_name == "C" | k==1, p_SB, 0), ##should be treatmentSB? cost3therapy = p_cost3, cost4zolimab = p_cost4, cost1cost = (p_cost1/3), cost5cost = (p_cost5/3), cost6cost = ifelse (state_name == "E" | state_name == "F", p_cost6, 0), mem1 = ifelse (time <= 1, 5, 0), mem2 = ifelse (time >=2 | time <= 6, 10, 0) ) ) } if (strategy_name == "SB") { j <- ifelse(state_name == "C", j <- j + 1, 0) list ( tpmatrix = p_SB, utility = p_utility_map, costs = list ( treatment = ifelse(state_name == "A" | time < 2, p_SB, 0), treatmentSB = ifelse(state_name == "C" | j==1, p_SA, 0), cost3therapy = p_cost3, cost4zolimab = p_cost4, cost1cost = (p_cost1/3), cost5cost = (p_cost5/3), cost6cost = ifelse (state_name == "E" | state_name == "F", p_cost6, 0), mem1 = ifelse (time <= 1, 5, 0), mem2 = ifelse (time >=2 | time <= 6, 10, 0) ) ) } }, times = c(1, Inf))

mod_def <- define_model(tparams_def = tparams_def, rng_def = rng_def, params = params)

econmod <- create_CohortDtstm(mod_def, input_data)

dincerti commented 1 year ago

@studentdoctorcoder this problem is best tackled without the define_model interface. I instead recommend setting up you transition probabilities more directly via tparams_transprobs.

Below is a complete example. Note that I assume no parameter uncertainty. If you want that, you could draw the transition matrices from Dirichlet distributions and use relevant distributions for the costs and utilities when using stateval_tbl (e.g., beta for utility, gamma for medical costs).

library("hesim")
library("data.table")
library("abind")

# Decision problem and model structure -----------------------------------------
# Note: states E and F are assumed to be death states since they cannot be 
# transitioned out of
strategies <- data.table(
  strategy_id = 1:2, strategy_name = c("SA", "SB")
)
patients <- data.table(patient_id = 1)
states <- data.table(
  state_id = 1:6, state_name = c("A", "B", "C", "D", "E", "F")
)
non_death_states <- states[!state_name %in% c("E", "F")]
hesim_dat <- hesim_data(
  strategies = strategies, patients = patients, states = non_death_states
)
N_SAMPLES <- 1  # No parameter uncertainty

# Transition matrices ----------------------------------------------------------
tpmat_sa <-  matrix(
  c(
    0.92, 0, 0.04, 0, 0.04, 0,
    0, 0.92, 0, 0.04, 0.04, 0,
    0, 0.46, 0, 0, 0.04, 0.5,
    0, 0, 0, 0.46, 0.04, 0.5,
    0, 0, 0, 0, 1, 0,
    0, 0, 0, 0, 0, 1
    ), 
  byrow = TRUE, 
  nrow = 6,
  dimnames = list(states$state_name, states$state_name)
)
tpmat_sb <- matrix(
  c(
    0.92, 0, 0.04, 0, 0.04, 0,
    0, 0.92, 0, 0.04, 0.04, 0,
    0, 0.46, 0, 0, 0.04, 0.5,
    0, 0, 0, 0.46, 0.04, 0.5,
    0, 0, 0, 0, 1, 0,
    0, 0, 0, 0, 0, 1
  ), 
  byrow = TRUE, 
  nrow = 6,
  dimnames = list(states$state_name, states$state_name)
)

transition_matrices <- tparams_transprobs(
  abind(tpmat_sa, tpmat_sb, along = 3),
  tpmatrix_id = tpmatrix_id(expand(hesim_dat), n_samples = N_SAMPLES)
)

# Utility and cost parameters --------------------------------------------------
utility_tbl <- stateval_tbl(
  data.table(
    state_id = non_death_states$state_id,
    est = c(A = 0.9, B = 0.8, C = 0.7, D = 0.6)
  ),
  dist = "fixed"
)

drugcost_tbl <- stateval_tbl(
  data.table(
    strategy_id = strategies$strategy_id,
    est = c(1000, 500)
  ),
  dist = "fixed"
)

medcost_tbl <- stateval_tbl(
  data.table(
    state_id = non_death_states$state_id,
    est = c(A = 500 , B = 600, C = 700, D = 800)
  ),
  dist = "fixed"
)

# Construction of the economic model -------------------------------------------
transmod <- CohortDtstmTrans$new(params = transition_matrices)
utilitymod <- create_StateVals(
  utility_tbl, 
  n = N_SAMPLES, 
  hesim_data = hesim_dat
)
costmods <- list(
  Drug = create_StateVals(drugcost_tbl, n = N_SAMPLES, hesim_data = hesim_dat),
  Medical = create_StateVals(medcost_tbl, n = N_SAMPLES, hesim_data = hesim_dat)
)

econmod <- CohortDtstm$new(
  trans_model = transmod,
  utility_model = utilitymod,
  cost_models = costmods
)

# Simulation of outcomes -------------------------------------------------------
labs <- get_labels(hesim_dat)

# State probabilities
econmod$sim_stateprobs(n_cycles = 60)
autoplot(econmod$stateprobs_, labels = labs)

# QALYs and costs
econmod$sim_qalys()
econmod$sim_costs()

# Decision analysis ------------------------------------------------------------
# Summarize costs and QALYs by treatment strategy
ce <- econmod$summarize()
format(summary(ce, labels = labs))

# Peform cost-effectiveness analysis
cea_out <- cea(ce, dr_qalys = .03, dr_costs = .03)
cea_pw_out <- cea_pw(ce, comparator = 1, dr_qalys = .03, dr_costs = .03)
format(icer(cea_pw_out))