Healthy-Sick-Dead Model

Introduction

This document runs through the PANDA framework for the basic healthy-sick-dead model.

Read in Excel File

library(tidyverse)
library(readxl)
library(here)
library(randtoolbox)
library(glue)
library(gt)
library(gtExtras) # remotes::install_github("jthomasmock/gtExtras")
sub_in_expression <- function(x) map2_chr(names(dist_lut),dist_lut, ~ifelse(grepl(glue('^{.x}'),x),gsub(glue('^{.x}'),.y,x),NA)) %>% na.omit() %>% as.vector()

gen_wcc <- function (n_cycles, method = c("Simpson1/3", "half-cycle", "none"))
{
  # Note this is DARTH CODE!!
  if (n_cycles <= 0) {
    stop("Number of cycles should be positive")
  }
  method <- match.arg(method)
  n_cycles <- as.integer(n_cycles)
  if (method == "Simpson1/3") {
    v_cycles <- seq(1, n_cycles + 1)
    v_wcc <- ((v_cycles%%2) == 0) * (2/3) + ((v_cycles%%2) !=
                                               0) * (4/3)
    v_wcc[1] <- v_wcc[n_cycles + 1] <- 1/3
  }
  if (method == "half-cycle") {
    v_wcc <- rep(1, n_cycles + 1)
    v_wcc[1] <- v_wcc[n_cycles + 1] <- 0.5
  }
  if (method == "none") {
    v_wcc <- rep(1, n_cycles + 1)
  }
  return(v_wcc)
}

#############################################
# Read in model structure and parameter info
#############################################

input_file <- normalizePath(here("healthy-sick-dead/healthy-sick-dead.xlsx"))
params_raw <- readxl::read_xlsx(input_file,sheet="parameters")
rates_raw <- readxl::read_xlsx(input_file,sheet = "transition_rates")
probs_raw <- readxl::read_xlsx(input_file,sheet="transition_probabilities")
payoffs_raw <- readxl::read_xlsx(input_file, sheet = "payoffs")

###########################
## Metaparameters
###########################
params_lut <- params_raw %>% select(param,description) %>% deframe()
psa_fmla_lut <- params_raw %>% filter(!is.na(distribution)) %>% select(param,distribution) %>% deframe()

v_names_states <- unique(probs_raw$from); v_names_states

[1] "Healthy" "Sick"    "Dead"   

v_n_states <- length(v_names_states)
v_names_str <- unique(rates_raw$strategy) %>% na.omit() %>% as.vector(); v_names_str

[1] "A" "B"

n_strategies <- length(v_names_str)
# THIS NEEDS TO BE DONE BASED ON THE FORMULA n_cycles <- params_raw %>% filter(param=="n_cycles") %>% pull(baseline) ; n_cycles
n_sim <- params_raw %>% filter(param == "n_sim") %>% pull(baseline); n_sim

[1] 1000

Parse Parameters

# Parameters with scalar values
dist_lut <- list("gamma" = "qgamma",
                 "lognormal" = "qlnorm",
                 "beta" = "qbeta",
                 "normal" = "qnorm",
                 "unif" = "qunif")

# Scalar parameters
params_sc <-
  params_raw  %>%
  filter(is.na(formula)) %>%
  select(param,baseline) %>%
  na.omit() %>%
  deframe() %>%
  as.list()

# parameters with uncertainty distributions
params_psa_ <-
  params_raw  %>%
  filter(!is.na(distribution)) %>%
  select(param,distribution) %>%
  na.omit() %>%
  deframe() %>%
  as.list()

psa_ <- halton(n=params_sc$n_sim, dim = length(params_psa_))
colnames(psa_) <- unlist(params_psa_)
colnames(psa_) <- lapply(colnames(psa_),function(x) sub_in_expression(x))
params_psa <-
  colnames(psa_) %>%
  map(~(glue("{gsub(')$','',.x)}, p={psa_[,.x]})")))  %>%
  set_names(names(params_psa_)) %>%
  bind_cols() %>% #glimpse()
  rowwise() %>%
  mutate_all(~eval(str2expression(.))) %>%
  ungroup() %>%
  as.list()




psa_sum <- 
  params_psa  %>% 
  bind_cols() %>% 
  na.omit() %>% 
  gather(param,value) %>% 
  group_by(param) %>% 
  summarise(
    n = n(),
    mean = mean(value),
    plt_val = list(value),
    .groups = "drop"
  ) %>% 
  mutate(baseline = params_sc[param]) %>% 
  unnest(cols=c(baseline)) %>% 
  arrange(param, desc(plt_val)) %>% 
  mutate(diff = abs((mean - baseline)/baseline)) %>% 
  mutate(desc = params_lut[param]) %>% 
  mutate(dist = psa_fmla_lut[param]) %>% 
  unnest(cols = c(dist)) %>% 
  select(param,desc,baseline,mean,dist,plt_val,n,diff)

Warning in cond & tbl$locname != "title": longer object length is not a multiple
of shorter object length

Warning in cond & tbl$locname != "subtitle": longer object length is not a
multiple of shorter object length

param	desc	baseline	mean	dist	plt_val	n
c_S	Annual cost of sick (S)	1,000	997	gamma(shape = 44.4, scale = 22.5)		1000
c_trtA	Cost of treatment A	25.0	24.8	gamma(shape = 12.5, scale = 2)		1000
c_trtB	Cost of treatment B	1,000	993	gamma(shape = 12, scale = 83.3)		1000
hr_S	Hazard ratio of death in S vs. H	3.00	3.00	lognormal(meanlog = log(3), sdlog = 0.01)		1000
hr_SB_trtB	Reduction in rate of disease progression (HS) as hazard ratio (HR)	0.950	0.950	lognormal(meanlog = log(0.95), sdlog = 0.02)		1000
r_HD	Annual mortality rate (H to D)	0.00600	0.00599	gamma(shape = 60, rate = 10000)		1000
r_HS	Annual disease onset transition rate	0.150	0.150	gamma(shape = 30,rate = 200)		1000
u_H	Utility weight of healthy (H)	1.00	0.985	beta(shape1 = 200, shape2 = 3)		1000
u_S	Utility weight of sick (S)	0.750	0.742	beta(shape1 = 130, shape2 = 45)		1000
Note: Highlighted rows have an absolute percent difference between the mean PSA value and the baseline value greater than 0.01