A Framework for Cost-Effectiveness Analyses

DecisionTree <- function(data, N = NA, ...) {
  call_obj <- match.call()
  structure(list(data = data,
                 N = N),
            call = call_obj,
            class = c("DecisionTree", "Model"))
}

MarkovModel <- function(init_probs = NA,
                        trans_matrix = NA,
                        cost_matrix = NA,
                        q_matrix = NA,
                        mapping = NA,
                        N = NA,
                        n_cycles = 2, ...) {
  
  call_obj <- match.call()
  structure(list(init_probs   = init_probs,
                 trans_matrix = trans_matrix,
                 cost_matrix  = cost_matrix,
                 q_matrix     = q_matrix,
                 mapping      = mapping,
                 N            = N,
                 n_cycles     = n_cycles),
            call = call_obj,
            class = c("MarkovModel", "Model"))
}

CombinedModel <- function(...) {
  models <- list(...)
  structure(models, class = "CombinedModel")
}

update_model <- function(model, result) {
  UseMethod("update_model")
}

update_model.MarkovModel <- function(model, result) {
  model$init_probs <- map_decision_to_markov(result, model$mapping)
  model
}

update_model.DecisionTree <- function(model, result) {
  model$data <- map_markov_to_decision(model, result)
  model
}

run_model.CombinedModel <- function(model_chain) {
  result <- list()
  
  for (i in seq_along(model_chain)) {
    current_model <- model_chain[[i]]
    
    if (i > 1) {
      # Update model via S3 dispatch
      current_model <- update_model(current_model, result[[i - 1]])
    }
    result[[i]] <- run_model(current_model)
  }
  
  structure(result,
            class = c("output", class(model_chain)))
}

run_model.DecisionTree <- function(model) {
  ## transform from model?
  
  # run model
  path_results <- calculate_pathways(model)
  
  ## transform from decision tree?
  
  # ce values
  res <- 
    path_results %>%
    group_by(treatment) %>%
    summarise(
      expected_cost = sum(terminal_prob * path_cost),
      expected_eff = sum(terminal_prob * path_eff)
    )
  
  structure(c(
    res,
    path_results = list(path_results)),
    class = c("output", class(model)))
}

run_model.MarkovModel <- function(model) {
  ## transform from model
  
  # run model
  out <- markov_model(start_pop = model$init_probs,
                      p_matrix = model$trans_matrix,
                      state_c_matrix = model$cost_matrix,
                      state_q_matrix = model$q_matrix,
                      n_cycles = model$n_cycles)
  
  ## transform from markov model
  
  res <- list()
  res$terminal <- out$final_state
  res$pop <- out$pop
  # ce values
  res$expected_cost <- out$total_costs
  res$expected_eff <- out$total_QALYs
  
  structure(res,
            class = c("output", class(model)))
}

decision_tree <- tribble(
  ~treatment, ~from,                ~to,              ~prob, ~cost,  ~eff,
  # Initial population split (true underlying state for Model A)
  "A",        "root",               "Healthy",        0.8,    0,      0, # Higher prob of starting healthy
  "A",        "root",               "Pre_clinical_Disease", 0.15,  0, 0,
  "A",        "root",               "Disease",        0.05,   0,      0,
  # Direct transitions to death states
  "A",        "root",               "Death_from_Disease", 0.0,    0,  0,
  "A",        "root",               "Other_Cause_Death",  0.0,    0,  0,
  "A",        "root",               "Treated",        0.0,    0,  0,
  # starting "Healthy" (true health status)
  "A",        "Healthy",            "FP",             0.1,   200,     0,    # False Positive
  "A",        "Healthy",            "TN",             0.9,   50,      0,    # True Negative
  # starting "Pre-clinical Disease" (true health status)
  # Assuming diagnostic test operates differently or has different outcomes
  "A",        "Pre_clinical_Disease", "TP_Preclinical", 0.7,   200,     0,  # True Positive for pre-clinical
  "A",        "Pre_clinical_Disease", "FN_Preclinical", 0.3,   50,      0,  # False Negative for pre-clinical
  # starting "Disease" (true health status)
  # Assuming diagnostic test operates differently or has different outcomes
  "A",        "Disease",            "TP_Disease",     1,   200,     0,     # True Positive for active disease
  "A",        "Disease",            "FN_Disease",     0,   50,      0,     # False Negative for active disease
)

# Mapping from decision tree terminal nodes to Markov states
# same for all treatments
mapping <- c(Treated = "Treated",
             FP = "Treated",
             TN = "Healthy",
             TP_Preclinical = "Treated",
             FN_Preclinical = "Pre_clinical_Disease",
             TP_Disease = "Treated",
             FN_Disease = "Disease",
             Death_from_Disease = "Death_from_Disease",
             Other_Cause_Death = "Other_Cause_Death")

trans_prob_mat <- array(
  c(
    # FROM Healthy (H)
    0.90,  # To Healthy
    0.05,  # To Pre-clinical Disease
    0.00,  # To Disease (assuming progression via pre-clinical)
    0.00,  # To Treated (assuming treatment only after diagnosis)
    0.00,  # To Death from Disease (assuming no direct death from healthy)
    0.05,  # To Other-Cause Death
    
    # FROM Pre-clinical Disease (P)
    0.00, # To Healthy (assuming no spontaneous recovery to healthy from pre-clinical)
    0.65, # To Pre-clinical Disease (staying untreated)
    0.20, # To Disease (progressing to active disease)
    0.00, # To Treated (initiating treatment at pre-clinical stage, if applicable)
    0.10, # To Death from Disease (due to progression during pre-clinical stage)
    0.05, # To Other-Cause Death
    
    # FROM Disease (D)
    0.00, # To Healthy
    0.00, # To Pre-clinical Disease (unlikely to regress to pre-clinical)
    0.60, # To Disease (remaining in active disease, untreated)
    0.20, # To Treated (initiating treatment for active disease)
    0.20, # To Death from Disease (progression of active disease)
    0.00, # To Other-Cause Death
    
    # FROM Treated (T) - Adjusted for new "Disease" state
    0.00, # To Healthy
    0.00, # To Pre-clinical Disease (unlikely from treated state)
    0.00, # To Disease (treatment failure/relapse to active disease)
    0.95, # To Treated (remaining treated and stable)
    0.00, # To Death from Disease (despite treatment)
    0.05, # To Other-Cause Death
    
    # FROM Death from Disease (DD) - Absorbing state
    0.00, 0.00, 0.00, 0.00, 1.00, 0.00, # All to Death from Disease
    
    # FROM Other-Cause Death (OCD) - Absorbing state
    0.00, 0.00, 0.00, 0.00, 0.00, 1.00  # All to Other-Cause Death
  ),
  dim = c(6, 6, 1),
  dimnames = list(to = states,
                  from = states,
                  scenario = "A")
) |> 
  aperm(c(2, 1, 3))  # rearrange dimensions

cost_mat <- array(
  c(0,     # Healthy
    0,     # Pre-clinical Disease
    0,     # Disease
    100,   # Treated
    0,     # Death from Disease
    0),    # Other-Cause Death
  dim = c(1, 6, 1), # 1 row, 6 states, 1 scenario
  dimnames = list(NULL, state = states, scenario = "A"))
q_mat <- array(
  c(1,      # Healthy
    0.8,    # Pre-clinical Disease
    0.5,    # Disease
    0.9,    # Treated
    0,      # Death from Disease
    0),     # Other-Cause Death
  dim = c(1, 6, 1), # 1 row, 6 states, 1 scenario
  dimnames = list(NULL, state = states, scenario = "A"))

dt <- DecisionTree(decision_tree)

mm10 <- MarkovModel(trans_matrix = trans_prob_mat,
                    cost_matrix = cost_mat,
                    q_matrix = q_mat,
                    init_probs = p_init,
                    mapping = mapping,
                    n_cycles = 10)

pair_model10 <- CombinedModel(dt, mm10)

screening_submodels10 <- rep(pair_model10, 2)
screening_model10 <- do.call(CombinedModel, args = screening_submodels10)
screening_results10 <- run_model(screening_model10)