R/evpi_estimators.R

Defines functions plot_evpi_threshold estimate_evpi_evppi_diff_threshold estimate_evpi_evppi estimate_evpi find_nmb_allparams get_all_parameter_list estimate_evppi find_nmb_setofparams get_parameter_list

Documented in estimate_evpi estimate_evpi_evppi estimate_evpi_evppi_diff_threshold estimate_evppi find_nmb_allparams find_nmb_setofparams get_all_parameter_list get_parameter_list plot_evpi_threshold

######################################################################
#' Function to estimate the expected value of partial perfect information
#' @param params_interest the main parameter of interest
#' @param value_params_interest value of the parameter of interest
#' @param names_params_needed names of needed parameters  from param
#' matrix returned
#' @param names_params_model names of parameters in the model
#' @param params_passed parameters passed while running the markov model
#' @param param_file all parameters required to run the model,provided with
#' name of the parameter, distribution, and parameters that define the
#' probability distribution
#' @param colnames_paramdistr col names where the parameter distribution
#' is defined
#' @return current_param_list list of parameters
#' @details
#' this function gets all the parameters except the parameter of interest
#' if the parameter is fixed, just read from file or a distribution,
#' then gets it from a distribution, or if it to be calculated, just
#' give it back as it is. If the parameter is not found, it will be taken
#' same that as that of in the model parameter
#' @examples
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0), method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#'  param_interest <- "tp_well_dis_co"
#'  colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#'  "Param2_value")
#'  value_param_interest <- 0.2
#'  names_params_needed <- colnames(list_markov[1, ]$param_matrix)
#'  names_params_model <- names(list_markov[1, ]$list_param_values)
#'  params_passed <- list_markov[1, ]$list_param_values
#' parameters <- get_parameter_list(param_interest, value_param_interest,
#' names_params_needed, names_params_model, params_passed,
#' param_file, colnames_paramdistr)
#' @export
get_parameter_list <- function(params_interest, value_params_interest,
                               names_params_needed, names_params_model,
                               params_passed,
                               param_file, colnames_paramdistr) {
  list_checks <- list(params_interest, value_params_interest,
                  names_params_needed, names_params_model,
                  params_passed, param_file, colnames_paramdistr)

  results <- sapply(list_checks, packDAMipd::check_null_na)
  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")

  param_distr_data <- utils::read.csv(param_file, row.names = NULL)
  param_colno <- IPDFileCheck::get_columnno_fornames(param_distr_data,
                                                     "parameter")
  if (length(param_colno) > 1) {
    stop("Parameter column should be only one")
  }

  param_colname <- colnames(param_distr_data)[param_colno]
  distr_col_no <- IPDFileCheck::get_columnno_fornames(param_distr_data,
                                                     "distribution")
  if (length(distr_col_no) > 1) {
    stop("Distribution column should be only one")
  }

  current_param_needed <- c()
  current_param_names <- c()
  for (j in 2:length(names_params_needed)) {
    # fixing the param_interest , so that we don't need to sample
    if (!(names_params_needed[j] %in%  params_interest)) {
      # sample others - but identify those require calculations
      #if the selected parameters is needed to run the model
      if (names_params_needed[j] %in% names_params_model) {
        #if the params_passed for the current parameter considered
        # is a numerical value, that has been either read from a file or
        # directly assigned. It might not be calculated. so those can be
        # found from the parameter file
        res <-
          suppressWarnings(as.numeric(params_passed[names_params_needed[j]]))
        if (!is.na(res)) {
          # check if the parameter is defined using a distribution or actually
          # fixed (for death rates)
          row_cor_param <- param_distr_data[param_distr_data[param_colname] ==
                                              names_params_needed[j], ]
          if (nrow(row_cor_param) > 1) {
            stop("One row should correspond to the paramter in parameter file")
          }
          if (nrow(row_cor_param) == 0) {
            curr_parm_value <- params_passed[names_params_needed[j]]
          }else{
            distr_col_val_param <- row_cor_param[[distr_col_no]]
            # if the entry corresponding to distribution is null, na or empty
            # keep this param fixed
            if (packDAMipd::check_null_na(distr_col_val_param) < 0 |
                distr_col_val_param == "") {
              curr_parm_value <-
                packDAMipd::get_parameter_read(names_params_needed[j],
                                               param_file)
            }else{
              curr_parm_value <-
              packDAMipd::get_parameter_def_distribution(names_params_needed[j],
                                                           param_file,
                                                          colnames_paramdistr)
            }
          }
        }else{
          curr_parm_value <- params_passed[names_params_needed[j]]
        }
        curr_parm <- names_params_needed[j]
      }

    }else{
      curr_parm_value <-
        value_params_interest[which(names_params_needed[j] == params_interest)]
      curr_parm <-
        params_interest[which(names_params_needed[j] == params_interest)]
    }
    current_param_needed <- append(current_param_needed, curr_parm_value)
    current_param_names <- append(current_param_names, curr_parm)
  }
  current_param_list <- current_param_needed
  names(current_param_list) <- current_param_names
  return(current_param_list)
}
######################################################################
#' Function to get NMB for a set of parameter
#' This is the inner loop of the two stage Monte Carlo process
#' @param params_interest the main parameter of interest
#' @param value_params_interest value for the parameter of interest
#' @param param_file all parameters required to run the model,provided with
#' name of parameter, distribution and parameters that define the probability
#' distribution
#' @param colnames_paramdistr column names where parameter distribution and
#' values are defined
#' @param list_markov list of markov models to estimate the NMB
#' @param threshold threshold values of WTP
#' @param comparator optional parameter if need ICER as the output
#' @return icer_nmb ICER and NMBs for the set of parameters and given
#' strategies
#' @details
#' For the given set of parameters find the nmb from the packDAMipd function
#' calculate_icer_nmb
#' @export
#' @examples
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0), method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#'  param_interest <- "tp_well_dis_co"
#'  colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#'  "Param2_value")
#'  value_param_interest <- 0.2
#'  threshold <- 2000
#'  res <- find_nmb_setofparams(param_interest, value_param_interest,
#'  param_file, colnames_paramdistr, list_markov, threshold)
find_nmb_setofparams <- function(params_interest, value_params_interest,
                                 param_file, colnames_paramdistr,
                                 list_markov, threshold, comparator = NULL) {
  list_checks <- list(params_interest, value_params_interest,
                      threshold, param_file, colnames_paramdistr)
  results <- sapply(list_checks, packDAMipd::check_null_na)
  if (sum(results) != 0)
    stop("Some of the parameters may by NULL or NA - please check")

  if (is.null(list_markov)) {
    stop("Markov models can not be NULL")
  }
  new_list_markov <- list()
  for (i in 1:nrow(list_markov)) {
    this <- list_markov[i, ]
    # names of strategies extracted - control, intervention etc
    # parameters needed get it from param_matrix
    names_params_needed <- colnames(list_markov[i, ]$param_matrix)
    names_params_model <- names(list_markov[i, ]$list_param_values)
    params_passed <- list_markov[i, ]$list_param_values
    current_param_list <- get_parameter_list(params_interest,
                                             value_params_interest,
                                             names_params_needed,
                                             names_params_model,
                                             params_passed, param_file,
                                             colnames_paramdistr)

    this_markov  <- packDAMipd::markov_model(current_strategy = this$strategy,
                                             cycles = this$cycles,
                                             initial_state = this$initial_state,
                                             discount = this$discount,
                                             parameter_values = current_param_list,
                                             half_cycle_correction =
                                               this$half_cycle_correction,
                                             state_cost_only_prevalent =
                                               this$state_cost_only_prevalent,
                                             state_util_only_prevalent =
                                               this$state_util_only_prevalent,
                                             method = this$method,
                                             startup_cost = this$startup_cost,
                                             startup_util = this$startup_util)

    new_list_markov[[length(new_list_markov) + 1]] <- this_markov
  }
  new_list_markov <- packDAMipd::combine_markov(new_list_markov)
  icer_nmb <- packDAMipd::calculate_icer_nmb(new_list_markov, threshold,
                                             comparator)
  return(icer_nmb)
}

######################################################################
#' Function to estimate the expected value of partial perfect information
#' This is the outer loop of the two stage Monte Carlo process
#' @param params_interest the parameters of interest
#' @param param_file all parameters required to run the model,provided with
#' name of parameter, distribution and parameters that define the probability
#' distribution
#' @param colnames_paramdistr column names where parameter distribution
#' and values are defined
#' @param list_markov list of markov models to estimate the NMB
#' @param threshold threshold values of WTP
#' @param outer_iterations number of iterations for outer loop
#' @param inner_iterations number of iterations for inner loop
#' @param comparator optional parameter if need ICER as the output
#' @return result result after evppi calculation, including three forms
#' of evppi
#' @source https://www.sciencedirect.com/science/article/pii/S1098301510605888
#' @examples
#'  colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#'  "Param2_value")
#'  inner_iterations <- 5
#'  outer_iterations <- 2
#'  threshold <- 2000
#'  param_interest <- "tp_well_dis_co"
#'  value_param_interest <- 0.2
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0), method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#'  param_interest <- "tp_well_dis_co"
#'  res <- estimate_evppi(param_interest, param_file,
#'  colnames_paramdistr, list_markov, threshold, outer_iterations,
#'  inner_iterations)
#' @export
estimate_evppi <- function(params_interest, param_file, colnames_paramdistr,
                           list_markov, threshold, outer_iterations,
                           inner_iterations,
                           comparator = NULL) {

  list_checks <- list(params_interest, param_file, colnames_paramdistr,
                      threshold, outer_iterations, inner_iterations)

  results <- sapply(list_checks, packDAMipd::check_null_na)

  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")

  if (is.null(list_markov)) {
    stop("Markov models can not be NULL")
  }

  cols_needed <- nrow(list_markov)
  all_mean_thetac_nmb_thetai <- matrix(0, nrow = outer_iterations,
                                       ncol = cols_needed)
  all_mean_thetac_maxT_nmbs_thetai <- matrix(0, nrow = outer_iterations,
                                             ncol = 1)
  all_maxT_mean_thetac_maxT_nmbs_thetai <- matrix(0, nrow = outer_iterations,
                                                  ncol = 1)

  #step 9
  for (outer_loop in 1:outer_iterations) {
    #Step 1
    value_params_interest <- c()
    for (m in seq_len(length(params_interest))) {
      value_params_interest <- cbind(value_params_interest,
                                     packDAMipd::get_parameter_def_distribution(
                                       params_interest[m], param_file, colnames_paramdistr))
    }


    all_nmbs_thetai <- matrix(0, nrow = inner_iterations, ncol = cols_needed)
    all_max_nmb_iteration <- list()
    #step  5
    for (inner_loop in 1:inner_iterations) {
      # step 2,3
      this_result <- find_nmb_setofparams(params_interest,
                                          value_params_interest,
                                          param_file, colnames_paramdistr,
                                          list_markov, threshold, comparator)
      nmbs_thetai <- (as.numeric(this_result$NMB))
      # store them
      all_nmbs_thetai[inner_loop, 1:cols_needed] <- nmbs_thetai
      # step 4
      maxT_nmbs_thetai <- max(as.numeric(nmbs_thetai))
      # store
      all_max_nmb_iteration[length(all_max_nmb_iteration) + 1] <-
        maxT_nmbs_thetai
    }
    # step 6
    mean_thetac_nmb_thetai <- colSums(all_nmbs_thetai) / inner_iterations
    # step 7
    mean_thetac_maxT_nmbs_thetai  <- mean(unlist(all_max_nmb_iteration))
    # step 8
    maxT_mean_thetac_maxT_nmbs_thetai <- max(mean_thetac_nmb_thetai)
    #store them
    all_mean_thetac_nmb_thetai[outer_loop, 1:cols_needed] <-
      mean_thetac_nmb_thetai
    all_mean_thetac_maxT_nmbs_thetai[outer_loop] <-
      mean_thetac_maxT_nmbs_thetai
    all_maxT_mean_thetac_maxT_nmbs_thetai[outer_loop] <-
      maxT_mean_thetac_maxT_nmbs_thetai
  }
  #step 10
  mean_thetai_mean_thetac_nmb_thetai <-
    colSums(all_mean_thetac_nmb_thetai) / outer_iterations
  mean_thetai_mean_thetac_maxT_nmb_thetai <-
    mean(all_mean_thetac_maxT_nmbs_thetai)
  mean_thetai_maxT_mean_thetac_nmbs_thetai <-
    mean(all_maxT_mean_thetac_maxT_nmbs_thetai)
  evppi <- mean_thetai_mean_thetac_maxT_nmb_thetai -
    mean_thetai_maxT_mean_thetac_nmbs_thetai
  result <- structure(list(
    max_mean_thetai_mean_thetac_nmb =
      max(mean_thetai_mean_thetac_nmb_thetai),
    mean_thetai_mean_thetac_nmb = mean_thetai_mean_thetac_nmb_thetai,
    mean_thetai_mean_thetac_maxT_nmb = mean_thetai_mean_thetac_maxT_nmb_thetai,
    mean_thetai_maxT_mean_thetac_nmbs =
      mean_thetai_maxT_mean_thetac_nmbs_thetai,
    Chilcott_evppi = evppi
  ))
  return(result)
}


######################################################################
#' Function to estimate the expected value of partial perfect information
#' @param names_params_needed names of needed parameters  from param matrix
#' returned
#' @param names_params_model names of parameters in the model
#' @param params_passed parameters passed while running the markov model
#' @param param_file all parameters required to run the model,provided
#' with name of parameter, distribution and parameters that define the
#' probability distribution
#' @param colnames_paramdistr col names where the parameter distribution is
#' defined
#' @return current_parameter_list set of all parameters
#' @keywords internal
#' @details
#' this function gets all the parameters except the parameter of interest
#' if they parameter is fixed, just read from file or a distribution,
#' then gets it from a distribution, or if it to be calculated, just
#' give it back as it is
#' @export
#' @examples
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0), method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#' param_file <- system.file("extdata", "table_param.csv", package = "packEVPI")
#' colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#'  "Param2_value")
#'  names_params_needed <- colnames(list_markov[1, ]$param_matrix)
#'  names_params_model <- names(list_markov[1, ]$list_param_values)
#'  params_passed <- list_markov[1, ]$list_param_values
#'  parameters <- get_all_parameter_list(names_params_needed, names_params_model,
#'  params_passed, param_file, colnames_paramdistr)
get_all_parameter_list <- function(names_params_needed, names_params_model,
                                   params_passed,
                                   param_file, colnames_paramdistr) {
  list_checks <- list(names_params_needed, names_params_model,
                      params_passed, param_file, colnames_paramdistr)
  results <- sapply(list_checks, packDAMipd::check_null_na)
  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")


  param_distr_data <- utils::read.csv(param_file, row.names = NULL)
  param_colno <- IPDFileCheck::get_columnno_fornames(param_distr_data,
                                                     "parameter")
  if (length(param_colno) > 1) {
    stop("Parameter column should be only one")
  }
  param_colname <- colnames(param_distr_data)[param_colno]
  distr_col_no <- IPDFileCheck::get_columnno_fornames(param_distr_data,
                                                      "distribution")
  if (length(distr_col_no) > 1) {
    stop("Distribution column should be only one")
  }
  current_param_needed <- c()
  current_param_names <- c()
  for (j in 2:length(names_params_needed)) {
    # sample all - but identify those require calculations
    #if the selected parameters is needed to run the model
    if (names_params_needed[j] %in% names_params_model) {
       #if the params_passed for the current parameter considered
      # is a numerical value, that has been either read from a file
      # or directly assigned. It might not be calculated. so those can be
      # found from the parameter file
      res <-
        suppressWarnings(as.numeric(params_passed[names_params_needed[j]]))
      if (!is.na(res)) {
        # check if the parameter is defined using a distribution or actually
        # fixed (for death rates)
        row_cor_param <- param_distr_data[param_distr_data[param_colname] ==
                                            names_params_needed[j], ]
        if (nrow(row_cor_param) > 1) {
          stop("One row should correspond to the paramter in parameter file")
        }
        if (nrow(row_cor_param) == 0) {
          curr_parm_value <- params_passed[names_params_needed[j]]
        }else{
          distr_col_val_param <- row_cor_param[[distr_col_no]]
          # if the entry corresponding to distribution is null, na or
          # empty keep this param fixed
          if (packDAMipd::check_null_na(distr_col_val_param) < 0 |
              distr_col_val_param == "") {
            curr_parm_value <-
              packDAMipd::get_parameter_read(names_params_needed[j],
                                             param_file)
          }else{
            curr_parm_value <-
              packDAMipd::get_parameter_def_distribution(names_params_needed[j],
                                                         param_file,
                                                         colnames_paramdistr)
          }
        }
      }else{
        curr_parm_value <- params_passed[names_params_needed[j]]
      }
      curr_parm <- names_params_needed[j]
    }
    current_param_needed <- append(current_param_needed, curr_parm_value)
    current_param_names <- append(current_param_names, curr_parm)
  }
  current_param_list <- current_param_needed
  names(current_param_list) <- current_param_names
  return(current_param_list)
}
######################################################################
#' Function to get NMB for a set of parameter
#' This is the inner loop of the two stage Monte Carlo process
#' @param param_file all parameters required to run the model,provided with
#' name of parameter, distribution and parameters that define the probability
#' distribution
#' @param colnames_paramdistr column names where parameter distribution and
#' values are defined
#' @param list_markov list of markov models to estimate the NMB
#' @param threshold threshold values of WTP
#' @param comparator optional parameter if need ICER as the output
#' @return icer_nmb ICER and NMB using all sampled parameters
#' @keywords internal
#' @details
#' For the given set of parameter find the nmb from the packDAMipd function
#' calculate_icer_nmb
#' @examples
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' threshold <- 2000
#' colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#' "Param2_value")
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0), method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#' nmbs <- find_nmb_allparams(param_file, colnames_paramdistr,
#' list_markov, threshold)
#' @export
find_nmb_allparams <- function(param_file, colnames_paramdistr,
                               list_markov, threshold, comparator = NULL) {

  list_checks <- list(param_file, colnames_paramdistr,
                      threshold)
  results <- sapply(list_checks, packDAMipd::check_null_na)
  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")

  if (is.null(list_markov)) {
    stop("Markov models can not be NULL")
  }
  new_list_markov <- list()
  for (i in 1:nrow(list_markov)) {
    this <- list_markov[i, ]
    # names of strategies extracted - control, intervention etc
    # parameters needed get it from param_matrix
    names_params_needed <- colnames(list_markov[i, ]$param_matrix)
    names_params_model <- names(list_markov[i, ]$list_param_values)
    params_passed <- list_markov[i, ]$list_param_values
    current_param_list <- get_all_parameter_list(names_params_needed,
                                                 names_params_model,
                                                 params_passed, param_file,
                                                 colnames_paramdistr)
    this_markov  <- packDAMipd::markov_model(current_strategy = this$strategy,
                                             cycles = this$cycles,
                                             initial_state = this$initial_state,
                                             discount = this$discount,
                                             parameter_values = current_param_list,
                                             half_cycle_correction =
                                               this$half_cycle_correction,
                                             state_cost_only_prevalent =
                                               this$state_cost_only_prevalent,
                                             state_util_only_prevalent =
                                               this$state_util_only_prevalent,
                                             method = this$method,
                                             startup_cost = this$startup_cost,
                                             startup_util = this$startup_util,
                                             no_cycles_year = this$no_cycles_year,
                                             start_discount =this$start_discount )
    new_list_markov[[length(new_list_markov) + 1]] <- this_markov
  }
  new_list_markov <- packDAMipd::combine_markov(new_list_markov)
  icer_nmb <- packDAMipd::calculate_icer_nmb(new_list_markov, threshold)
  return(icer_nmb)
}
######################################################################
#' Function to estimate the expected value of partial perfect information
#' This is the outer loop of the two stage Monte Carlo process
#' @param param_file all parameters required to run the model,provided
#' with name of parameter, distribution and parameters that define the
#' probability distribution
#' @param colnames_paramdistr column names where parameter distribution and
#' values are defined
#' @param list_markov list of markov models to estimate the NMB
#' @param threshold threshold values of WTP
#' @param iterations number of iterations needed
#' @param comparator optional parameter if need ICER as the output
#' @return result result after evpi calculation
#' @source https://www.sciencedirect.com/science/article/pii/S1098301510605888
#' @export
#' @examples
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#' "Param2_value")
#' threshold  <-  20000
#' iterations <- 1
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_file <- system.file("extdata", "table_param.csv",
#' package = "packEVPI")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0), method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#' estimate_evpi(param_file,colnames_paramdistr,list_markov, threshold,
#' iterations)
estimate_evpi <- function(param_file, colnames_paramdistr,
                          list_markov, threshold, iterations,
                          comparator = NULL) {

  list_checks <- list(param_file, colnames_paramdistr,
                      threshold, iterations)

  results <- sapply(list_checks, packDAMipd::check_null_na)
  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")
  if (is.null(list_markov)) {
    stop("Markov models can not be NULL")
  }
  names_strategy <- list()
  for (i in 1:nrow(list_markov)) {
    names_strategy <- append(names_strategy, list_markov[[i]]$name_strategy)
  }
  names_strategy <- unlist(names_strategy)


  cols_needed <- nrow(list_markov)
  all_nmbs_thetai <- matrix(0, nrow = iterations, ncol = cols_needed)
  all_max_nmb_iteration <- list()
  for (loop in 1:iterations) {
    this_result <- find_nmb_allparams(param_file, colnames_paramdistr,
                                      list_markov, threshold, comparator)
    nmbs_thetai <- (as.numeric(this_result$NMB))
    # store them
    all_nmbs_thetai[loop, 1:cols_needed] <- nmbs_thetai
    # step 4
    maxT_nmbs_thetai <- max(as.numeric(nmbs_thetai))
    # store
    all_max_nmb_iteration[length(all_max_nmb_iteration) + 1] <-
      maxT_nmbs_thetai

  }
  mean_nmb_thetai <- colSums(all_nmbs_thetai) / iterations
  mean_maxT_nmbs_thetai  <- mean(unlist(all_max_nmb_iteration))
  maxT_mean_nmbs_thetai <- max(mean_nmb_thetai)
  col_with_max_mean_nmbs <- which(mean_nmb_thetai == maxT_mean_nmbs_thetai)
  trt_with_max_mean_nmbs <- names_strategy[col_with_max_mean_nmbs]
  evpi <- mean_maxT_nmbs_thetai - maxT_mean_nmbs_thetai
  result <- structure(list(
    mean_maxT_nmbs = mean_maxT_nmbs_thetai,
    maxT_mean_nmbs = maxT_mean_nmbs_thetai,
    treatment_max_mean_nmbs = trt_with_max_mean_nmbs,
    evpi = evpi
  ))
  return(result)
}

######################################################################
#' Function to estimate the expected value of partial perfect information
#' This is the outer loop of the two stage Monte Carlo process
#' @param params_interest the main parameter of interest
#' @param param_file all parameters required to run the model,provided with
#' name of the parameter, distribution and parameters that define the
#' probability distribution
#' @param colnames_paramdistr column names where parameter distribution
#' and values are defined
#' @param list_markov list of markov models to estimate the NMB
#' @param threshold threshold values of WTP
#' @param outer_iterations number of iterations for outer loop
#' @param inner_iterations number of iterations for inner loop
#' @param comparator optional parameter if need ICER as the output
#' @return result evpi and evppi for a single param
#' @source https://www.sciencedirect.com/science/article/pii/S1098301510605888
#' @export
#' @examples
#' param_file <- system.file("extdata", "table_param.csv", package = "packEVPI")
#' colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#'  "Param2_value")
#'  params_interest <- "tp_well_dis_co"
#'  outer_iterations <- 1
#'  inner_iterations <- 2
#' threshold <- 2000
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#'  utility = "utility_dis_co")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat,
#'  c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#'  "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#'  param_list <- packDAMipd::define_parameters(
#'  tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#'  param_file),
#'  tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#'  param_file),
#'  tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#'  tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'   discount = c(0, 0), method = "half cycle correction", param_list)
#'   list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#'  res <- estimate_evpi_evppi(params_interest, param_file, colnames_paramdistr,
#'  list_markov, threshold, outer_iterations, inner_iterations)
estimate_evpi_evppi <- function(params_interest, param_file,
                                colnames_paramdistr, list_markov,
                                threshold, outer_iterations,
                                inner_iterations,
                                comparator = NULL) {

  list_checks <- list(params_interest, param_file, colnames_paramdistr,
                      threshold, outer_iterations, inner_iterations)
  results <- sapply(list_checks, packDAMipd::check_null_na)

  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")

  if (is.null(list_markov)) {
    stop("Markov models can not be NULL")
  }

  results_evpi <- estimate_evpi(param_file, colnames_paramdistr,
                                list_markov, threshold, outer_iterations,
                                comparator)
  results_evppi <- estimate_evppi(params_interest, param_file,
                                  colnames_paramdistr, list_markov,
                                  threshold, outer_iterations,
                                  inner_iterations, comparator)

  Chilcott_evppi <- results_evppi$Chilcott_evppi
  Brennan_evppi <- results_evppi$mean_thetai_maxT_mean_thetac_nmbs -
    results_evpi$maxT_mean_nmbs
  Claxton_evppi <- results_evppi$mean_thetai_maxT_mean_thetac_nmbs -
    results_evppi$max_mean_thetai_mean_thetac_nmb

  result <- structure(list(
    evpi = results_evpi,
    evppi_2level = results_evppi,
    Chilcott_evppi = Chilcott_evppi,
    Brennan_evppi = Brennan_evppi,
    Claxton_evppi = Claxton_evppi
  ))
  return(result)
}

######################################################################
#' Function to estimate the expected value of partial perfect information
#' This is the outer loop of the two stage Monte Carlo process
#' @param parameters_of_interest parameters of interest
#' @param param_file all parameters required to run the model,provided
#' with name of the parameter, distribution and parameters that define
#' the probability distribution
#' @param colnames_paramdistr column names where parameter distribution and
#' values are defined
#' @param list_markov list of markov models to estimate the NMB
#' @param threshold_values threshold values of WTP
#' @param outer_iterations number of iterations for outer loop
#' @param inner_iterations number of iterations for inner loop
#' @param comparator optional parameter if need ICER as the output
#' @return result result for evpi and evppi for different values of threshold
#' @source https://www.sciencedirect.com/science/article/pii/S1098301510605888
#' @export
#' @examples
#' param_file <- system.file("extdata", "table_param.csv", package = "packEVPI")
#' colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#' "Param2_value")
#' parameters_of_interest <- "tp_well_dis_co"
#' outer_iterations <- 1
#' inner_iterations <- 2
#' threshold <- 2000
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#'  utility = "utility_dis_co")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat,
#'  c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#'  "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#'  param_list <- packDAMipd::define_parameters(
#'  tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#'  param_file),
#'  tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#'  param_file),
#'  tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#'  tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#'  tp_dead_dead =  1,
#'  cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#'  cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#'  cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#'  cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#'  utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#'  param_file),
#'  utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#'  param_file),
#'  tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#'  tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#'  tp_dis_dis_co = "1-( tp_dis_dead)",
#'  tp_dis_dis_in = "1-( tp_dis_dead)")
#'  this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#'  disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#'  utility = "utility_dis_in")
#'  dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#'  tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#'  colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#'  tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in", "tp_well_dead", "tp_dis_dis_in","tp_dis_dead",
#'  "tp_dead_dead"), colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'   discount = c(0, 0), method = "half cycle correction", param_list)
#'   list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#' res <- estimate_evpi_evppi_diff_threshold(parameters_of_interest, param_file,
#' colnames_paramdistr,list_markov, threshold, outer_iterations,
#' inner_iterations)
estimate_evpi_evppi_diff_threshold <- function(parameters_of_interest,
                                               param_file,
                                               colnames_paramdistr,
                                               list_markov, threshold_values,
                                               outer_iterations = NULL,
                                               inner_iterations= NULL,
                                               comparator = NULL) {
  list_checks <- list(parameters_of_interest, param_file, colnames_paramdistr,
                      threshold_values)
  results <- sapply(list_checks, packDAMipd::check_null_na)

  if (sum(results) != 0)
    stop("Some of the parameters may by null or na - please check")

  if (is.null(list_markov)) {
    stop("Markov models can not be NULL")
  }
  if (is.null(outer_iterations))
    outer_iterations <- 500
  if (is.null(inner_iterations))
    inner_iterations <- 1000

  all_evpis <- as.data.frame(0, nrow = length(threshold_values), ncol = 4)
  num_trt_strategies <- nrow(list_markov)

  all_evppis <- as.data.frame(0, nrow = length(threshold_values),
                              ncol = num_trt_strategies + 3)
  named_evppis <- as.data.frame(0, nrow = length(threshold_values), ncol = 3)

  number_threshold <- length(threshold_values)
  for (jj in 1:number_threshold) {
    this_threshold <- threshold_values[jj]
    res_a_param <- estimate_evpi_evppi(parameters_of_interest,
                                       param_file,
                                       colnames_paramdistr,
                                       list_markov,
                                       this_threshold,
                                       outer_iterations,
                                       inner_iterations,
                                       comparator)
    all_evpis[jj, 1] <- as.numeric(res_a_param$evpi[1])
    all_evpis[jj, 2] <- as.numeric(res_a_param$evpi[2])
    all_evpis[jj, 3] <- (res_a_param$evpi[3])
    all_evpis[jj, 4] <- as.numeric(res_a_param$evpi[4])
    all_evppis[jj, 1] <- as.numeric(res_a_param$evppi[1])
    names_strategy <- list()
    for (i in 1:num_trt_strategies) {
      all_evppis[jj, i + 1] <- as.numeric(res_a_param$evppi[[2]][i])
      names_strategy <- append(names_strategy,
                               list_markov[i, ]$strategy$name_strategy)
    }
    names_strategy <- unlist(names_strategy)
    all_evppis[jj, num_trt_strategies + 2] <-
      as.numeric(res_a_param$evppi[3])
    all_evppis[jj, num_trt_strategies + 3] <-
      as.numeric(res_a_param$evppi[4])

    named_evppis[jj, 1] <- as.numeric(res_a_param$Chilcott_evppi)
    named_evppis[jj, 2] <- as.numeric(res_a_param$Brennan_evppi)
    named_evppis[jj, 3] <- as.numeric(res_a_param$Claxton_evppi)

  }
  all_evpis[["thresholds"]] <- threshold_values
  all_evppis[["thresholds"]] <- threshold_values
  named_evppis[["thresholds"]] <- threshold_values
  colnames(all_evpis) <- c("Mean max NMB", "Max Mean NMB",
                           "Treatment with Max Mean NMB",
                           "EVPI", "Thresholds")
  temp_col_names <- list()
  for (i in 1:num_trt_strategies) {
    temp_col_names <- append(temp_col_names,
                             paste("Mean Thetai Mean Thetac NMB ",
                                   names_strategy[i]))
  }
  temp_col_names <- unlist(temp_col_names)
  colnames(all_evppis) <- c("Max Mean Thetai Mean Thetac NMB",
                            temp_col_names,
                            "Mean Thetai Mean Thetac Max NMB",
                            "Mean Thetai Max Mean Thetac NMB",
                            "Thresholds")
  colnames(named_evppis) <- c("Chilcott", "Brennan", "Claxton", "Thresholds")
  result <- list(Parameter = parameters_of_interest,
                 EVPIs = all_evpis,
                 EVPPIs = all_evppis,
                 ThreeEVPPIs = named_evppis)

  return(result)
}

######################################################################
#' Function to estimate the expected value of partial perfect information
#' This is the outer loop of the two stage Monte Carlo process
#' @param result_evpi_evppi result from estimation of evpi and evppi
#' @return plot the plot
#' @export
#' @examples
#' param_file <- system.file("extdata", "table_param.csv", package = "packEVPI")
#' well <- packDAMipd::health_state("well", cost = "cost_well_co", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_co",
#' utility = "utility_dis_co")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat,
#' c("tp_well_well_co","tp_well_dis_co","tp_well_dead", "tp_dis_dis_co",
#' "tp_dis_dead", "tp_dead_dead"),colnames(tmat))
#' health_states <- packDAMipd::combine_state(well, disabled, dead)
#' this.strategy <- packDAMipd::strategy(tm, health_states, "control")
#' param_list <- packDAMipd::define_parameters(
#' tp_well_dis_co = packDAMipd::get_parameter_read("tp_well_dis_co",
#' param_file),
#' tp_well_dis_in =  packDAMipd::get_parameter_read("tp_well_dis_in",
#' param_file),
#' tp_well_dead =  packDAMipd::get_parameter_read("tp_well_dead", param_file),
#' tp_dis_dead =  packDAMipd::get_parameter_read("tp_dis_dead", param_file),
#' tp_dead_dead =  1,
#' cost_well_co =  packDAMipd::get_parameter_read("cost_well_co", param_file),
#' cost_well_in =  packDAMipd::get_parameter_read("cost_well_in", param_file),
#' cost_dis_co =  packDAMipd::get_parameter_read("cost_dis_co", param_file),
#' cost_dis_in =  packDAMipd::get_parameter_read("cost_dis_in", param_file),
#' utility_dis_co =  packDAMipd::get_parameter_read("utility_dis_co",
#' param_file),
#' utility_dis_in =  packDAMipd::get_parameter_read("utility_dis_in",
#' param_file),
#' tp_well_well_co = "1-(tp_well_dis_co + tp_well_dead)",
#' tp_well_well_in = "1-(tp_well_dis_in + tp_well_dead)",
#' tp_dis_dis_co = "1-( tp_dis_dead)",
#' tp_dis_dis_in = "1-( tp_dis_dead)")
#' this_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#' discount = c(0, 0), method = "half cycle correction", param_list)
#' well <- packDAMipd::health_state("well", cost = "cost_well_in", utility = 1)
#' disabled <- packDAMipd::health_state("disabled", cost = "cost_dis_in",
#' utility = "utility_dis_in")
#' dead <- packDAMipd::health_state("dead", cost = 0, utility = 0)
#' tmat <- rbind(c(1, 2, 3), c(NA, 4, 5), c(NA, NA, 6))
#' colnames(tmat) <- rownames(tmat) <- c("well", "disabled", "dead")
#' tm <- packDAMipd::populate_transition_matrix(3, tmat, c("tp_well_well_in",
#'  "tp_well_dis_in","tp_well_dead", "tp_dis_dis_in","tp_dis_dead", "tp_dead_dead"),
#'  colnames(tmat))
#'  health_states <- packDAMipd::combine_state(well, disabled, dead)
#'  this.strategy <- packDAMipd::strategy(tm, health_states, "intervention")
#'  sec_markov <- packDAMipd::markov_model(this.strategy, 24, c(1000, 0, 0),
#'  discount = c(0, 0),method = "half cycle correction", param_list)
#'  list_markov <- packDAMipd::combine_markov(list(this_markov, sec_markov))
#'  parameter_of_interest <- "tp_well_dis_co"
#'  colnames_paramdistr  <- c("Param1_name", "Param1_value", "Param2_name",
#'   "Param2_value")
#'   threshold_values <- c(5000, 10000, 150000, 20000)
#'   res <- estimate_evpi_evppi_diff_threshold(parameter_of_interest, param_file,
#'   colnames_paramdistr, list_markov, threshold_values, outer_iterations = 3,
#'   inner_iterations = 5)
#'   plot_evpi_threshold(res)
plot_evpi_threshold <- function(result_evpi_evppi) {
  if (is.null(result_evpi_evppi)) {
    stop("result from EVPI /EVPPI calculation can not be NULL")
  }
  required_names <- c("Parameter", "EVPIs", "EVPPIs", "ThreeEVPPIs")
  current_names <- names(result_evpi_evppi)
  if (sum(current_names %in% required_names) != length(required_names)) {
    stop("Error - Given result do not contain all the required information")
  }
  EVPIs <- result_evpi_evppi$EVPIs
  EVPPIs <- result_evpi_evppi$ThreeEVPPIs

  name_file_plot <- paste0("EVPI for threshold values for ",
                           result_evpi_evppi$Parameter, ".pdf", sep = "")
  grDevices::pdf(name_file_plot)
  p <- ggplot2::ggplot(EVPIs, ggplot2::aes_(x = ~Thresholds,
                                            y = ~EVPI)) +
    ggplot2::geom_line() +
    ggplot2::labs(x = "Threshold") +
    ggplot2::labs(y = "EVPI") +
    ggplot2::theme(legend.title = ggplot2::element_blank())
  graphics::plot(p) # plot result
  grDevices::dev.off()

  name_file_plot <- paste0("EVPPI for threshold values for ",
                           result_evpi_evppi$Parameter, ".pdf", sep = "")
  grDevices::pdf(name_file_plot)

  EVPPIs_melted <- reshape2::melt(EVPPIs, id.var = "Thresholds")
  p <- ggplot2::ggplot(EVPPIs_melted, ggplot2::aes_(x = ~Thresholds,
                                                    y = ~value, color = ~variable)) +
    ggplot2::geom_line() +
    ggplot2::labs(x = "Threshold") +
    ggplot2::labs(y = "EVPPI") +
    ggplot2::theme(legend.title = ggplot2::element_blank())
  graphics::plot(p) # plot result
  grDevices::dev.off()
}
sheejamk/packEVPI documentation built on April 7, 2023, 8:48 a.m.