tests/testthat/example_project/R/run_sickSicker_model.R

In assertHE: Visualisation and Verification of Health Economic Decision Models

#' Run the SickSicker State-Transition Model
#'
#' @description Run the SickSicker Markov model.
#'
#' @param params_ List containing the model parameters. This list should contain
#' the following objects:
#' \describe{
#'  \item{age_init_}{Age at baseline, default = `25`}
#'  \item{age_max_}{Maximum age of follow up, default = `55`}
#'  \item{discount_rate_}{Discount rate for costs and QALYs}
#'  \item{p_HD}{Probability to die when healthy}
#'  \item{p_HS1}{Probability to become sick when healthy}
#'  \item{p_S1H}{Probability to become healthy when sick}
#'  \item{p_S1S2}{Probability to become sicker when sick}
#'  \item{hr_S1}{Hazard ratio of death in sick v healthy}
#'  \item{hr_S2}{Hazard ratio of death in sicker v healthy}
#'  \item{c_H}{Cost of remaining one cycle in the healthy state}
#'  \item{c_S1}{Cost of remaining one cycle in the sick state}
#'  \item{c_S2}{Cost of remaining one cycle in the sicker state}
#'  \item{c_Trt}{Cost of treatment(per cycle)}
#'  \item{c_D}{Cost of being in the death state}
#'  \item{u_H}{Utility when healthy}
#'  \item{u_S1}{Utility when sick}
#'  \item{u_S2}{Utility when sicker}
#'  \item{u_D}{Utility when dead}
#'  \item{u_Trt}{Utility when being treated}
#' }
#'
#' @return A matrix containing the costs and QALYs generated by running the
#' SickSicker model.
#'
#' @family modelrun
#'
#' @export
#' }
run_sickSicker_model <- function(
    params_ = NULL
) {

  ## Run the model using the objects in the params list:

  with(
    data = params_,
    expr = {
      ## Calculate interim parameters:

      # strategy names
      Strategies = c("No Treatment", "Treatment")

      # time horizon, number of cycles
      n_t = age_max_ - age_init_
      # the 4 states of the model: Healthy (H), Sick (S1), Sicker (S2), Dead (D)
      v_n  <- c("H", "S1", "S2", "D")

      # rate of death in healthy
      r_HD    <- - log(1 - p_HD)
      # rate of death in sick
      r_S1D   <- hr_S1 * r_HD
      # rate of death in sicker
      r_S2D   <- hr_S2 * r_HD

      # probability of death in sick
      p_S1D   <- 1 - exp(-r_S1D)
      # probability of death in sicker
      p_S2D   <- 1 - exp(-r_S2D)

      # initial Markov trace
      initial_trace <- c("H" = 1, "S1" = 0, "S2" = 0, "D" = 0)

      # =================================== #
      #  Check the initial trace conforms   #
      # =================================== #
      assertHE::check_init(x = initial_trace)
      # =================================== #
      # =================================== #

      # create vectors for the costs and utility of each treatment group
      c_trt    <- c("H" = c_H, "S1" = c_S1 + c_Trt, "S2" = c_S2 + c_Trt, "D" = c_D)
      c_no_trt <- c("H" = c_H, "S1" = c_S1, "S2" = c_S2, "D" = c_D)
      u_trt    <- c("H" = u_H, "S1" = u_Trt, "S2" = u_S2, "D" = u_D)
      u_no_trt <- c("H" = u_H, "S1" = u_S1, "S2" = u_S2, "D" = u_D)

      ## Calculate the discount weight for each cycle:

      v_dw <- calculate_discounting_weights(
        discount_rate_ = discount_rate_,
        time_horizon_ = n_t,
        first_cycle_ = TRUE)

      ## Define the model's transition matirix:

      m_P <- define_transition_matrix(
        states_nms_ = v_n,
        transition_probs_ = c(
          1 - (p_HS1 + p_HD),                        p_HS1,         0,  p_HD,
          p_S1H,              1 - (p_S1H + p_S1S2 + p_S1D),    p_S1S2, p_S1D,
          0,                                             0, 1 - p_S2D, p_S2D,
          0,                                             0,         0,     1
        )
      )

      # ============================== #
      #  Check the probability matrix  #
      # ============================== #
      assertHE::check_trans_prob_mat(m_P)
      # ============================== #
      # ============================== #

      ## Create the model's Markov trace:

      m_TR <- create_Markov_trace(
        transition_matrix_ = m_P,
        time_horizon_ = n_t,
        states_nms_ = v_n,
        initial_trace_ = initial_trace
      )

      # ============================== #
      #  we can check the markov trace #
      # ============================== #
      assertHE::check_markov_trace(m_TR = m_TR,
                         dead_state = "D",
                         confirm_ok = TRUE)
      # ============================== #
      # ============================== #


      ## Calculate costs:

      v_C_no_trt <- calculate_costs(
        Markov_trace_ = m_TR,
        costs_ = c_no_trt,
        discounting_weights_ = v_dw
      )
      v_C_trt <- calculate_costs(
        Markov_trace_ = m_TR,
        costs_ = c_trt,
        discounting_weights_ = v_dw
      )
      tc_no_trt  <- sum(v_C_no_trt)
      tc_trt     <- sum(v_C_trt)

      ## Calculate QALYs:

      v_E_no_trt <- calculate_QALYs(
        Markov_trace_ = m_TR,
        utilities_ = u_no_trt,
        discounting_weights_ = v_dw
      )
      v_E_trt <- calculate_QALYs(
        Markov_trace_ = m_TR,
        utilities_ = u_trt,
        discounting_weights_ = v_dw
      )
      te_no_trt <- sum(v_E_no_trt)
      te_trt    <- sum(v_E_trt)

      ## Prepare results:
      results <- c(
        "Cost_no_treatment"  = tc_no_trt ,
        "Cost_treatment"     = tc_trt ,
        "QALYs_no_treatment" = te_no_trt,
        "QALYs_treatment"    = te_trt
      )

      return(results)
    }
  )
}