## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

# Numeric notation in this document
# With thanks to https://stackoverflow.com/questions/18965637/set-global-thousand-separator-on-knitr/18967590#18967590
knitr::knit_hooks$set(
    inline = function(x) {
        if(!is.numeric(x)){
            x
        } else if(x<100) {
            prettyNum(round(x, 3), big.mark=",")
        } else {
            prettyNum(round(x, 0), big.mark=",")
        }
    }
)

## ----setup, echo=TRUE, message=FALSE------------------------------------------
library(dplyr)
library(ggplot2)
library(lubridate)
library(flexsurv)
library(heemod)
library(tidyr)
library(dynamicpv)

## ----constants----------------------------------------------------------------
# Time constants
days_in_year <- 365.25
days_in_week <- 7
cycle_years <- days_in_week / days_in_year # Duration of a one week cycle in years

# Time horizon (years) and number of cycles
thoz <- 20
Ncycles <- ceiling(thoz/cycle_years)

# Real discounting
disc_year <- 0.03 # Per year
disc_cycle <- (1+disc_year)^(cycle_years) - 1 # Per cycle

# Price inflation
infl_year <- 0.025 # Per year
infl_cycle <- (1+infl_year)^(cycle_years) - 1 # Per cycle

# Nominal discounting
nomdisc_year <- (1+disc_year)*(1+infl_year) - 1
nomdisc_cycle <- (1+nomdisc_year)^(cycle_years) - 1 # Per cycle

## ----heemod-------------------------------------------------------------------
# State names
state_names = c(
  progression_free = "PF",
  progression = "PD",
  death = "Death"
  )

# PFS distribution for SoC with Exp() distribution and mean of 50 weeks
surv_pfs_soc <- heemod::define_surv_dist(
  distribution = "exp",
  rate = 1/50
)

# OS distribution for SoC with Lognorm() distribution, meanlog = 4.5, sdlog = 1
# This implies a mean of exp(4 + 0.5 * 1^2) = exp(4.5) = 90 weeks
surv_os_soc <- heemod::define_surv_dist(
  distribution = "lnorm",
  meanlog = 4,
  sdlog = 1
)

# PFS and OS distributions for new
surv_pfs_new <- heemod::apply_hr(surv_pfs_soc, hr=0.5)
surv_os_new <- heemod::apply_hr(surv_os_soc, hr=0.6)

# Define partitioned survival model, soc
psm_soc <- heemod::define_part_surv(
  pfs = surv_pfs_soc,
  os = surv_os_soc,
  terminal_state = FALSE,
  state_names = state_names
  )

# Define partitioned survival model, soc
psm_new <- heemod::define_part_surv(
  pfs = surv_pfs_new,
  os = surv_os_new,
  terminal_state = FALSE,
  state_names = state_names
  )

# Parameters
params <- heemod::define_parameters(
  # Discount rate
  disc = disc_cycle,
  # Disease management costs
  cman_pf = 80,
  cman_pd = 20,
  # Drug acquisition costs - the SoC regime only uses SoC drug, the New regime only uses New drug
  cdaq_soc = dispatch_strategy(
    soc = 400,
    new = 0
  ),
  cdaq_new = dispatch_strategy(
    soc = 0,
    new = 1500
  ),
  # Drug administration costs
  cadmin = dispatch_strategy(
    soc = 50,
    new = 75
  ),
  # Adverse event risks
  risk_ae = dispatch_strategy(
    soc = 0.08,
    new = 0.1
  ),
  # Adverse event average costs
  uc_ae = 2000,
  # Subsequent treatments
  csubs = dispatch_strategy(
    soc = 1200,
    new = 300
  ),
  # Health state utilities
  u_pf = 0.8,
  u_pd = 0.6,
)

# Define PF states
state_PF <- heemod::define_state(
  # Costs for the state
  cost_daq_soc = discount(cdaq_soc, disc_cycle),
  cost_daq_new = discount(cdaq_new, disc_cycle),
  cost_dadmin = discount(cadmin, disc_cycle),
  cost_dman = discount(cman_pf, disc_cycle),  
  cost_ae = risk_ae * uc_ae,
  cost_subs = 0,
  cost_total = cost_daq_soc + cost_daq_new + cost_dadmin + cost_dman + cost_ae + cost_subs,
  # Health utility, QALYs and life years
  pf_year = discount(cycle_years, disc_cycle),
  life_year = discount(cycle_years, disc_cycle),
  qaly = discount(cycle_years * u_pf, disc_cycle)
  )

# Define PD states
state_PD <- heemod::define_state(
  # Costs for the state
  cost_daq_soc = 0,
  cost_daq_new = 0,
  cost_dadmin = 0,
  cost_dman = discount(cman_pd, disc_cycle),  
  cost_ae = 0,
  cost_subs = discount(csubs, disc_cycle),
  cost_total = cost_daq_soc + cost_daq_new + cost_dadmin + cost_dman + cost_ae + cost_subs,
  # Health utility, QALYs and life years
  pf_year = 0,
  life_year = heemod::discount(cycle_years, disc_cycle),
  qaly = heemod::discount(cycle_years * u_pd, disc_cycle)
  )

# Define Death state
state_Death <- heemod::define_state(
  # Costs are zero
  cost_daq_soc = 0,
  cost_daq_new = 0,
  cost_dadmin = 0,
  cost_dman = 0,
  cost_ae = 0,
  cost_subs = 0,
  cost_total = cost_daq_soc + cost_daq_new + cost_dadmin + cost_dman + cost_ae + cost_subs,
  # Health outcomes are zero
  pf_year = 0,
  life_year = 0,
  qaly = 0,
)

# Define strategy for SoC
strat_soc <- heemod::define_strategy(
    transition = psm_soc,
    "PF" = state_PF,
    "PD" = state_PD,
    "Death" = state_Death
  )

# Define strategy for new
strat_new <- heemod::define_strategy(
  transition = psm_new,
  "PF" = state_PF,
  "PD" = state_PD,
  "Death" = state_Death
)

# Create heemod model
heemodel <- heemod::run_model(
  soc = strat_soc,
  new = strat_new,
  parameters = params,
  cycles = Ncycles,
  cost = cost_total,
  effect = qaly,
  init = c(1, 0, 0),
  method = "life-table"
)

## ----dynpricing---------------------------------------------------------------
# Dates
# Date of calculation = 1 September 2025
doc <- lubridate::ymd("20250901")
# Date of LOE for SoC = 1 January 2028
loe_soc_start <- lubridate::ymd("20280101")
# Maturation of SoC prices by LOE + 1 year, i.e. = 1 January 2029
loe_soc_end <- lubridate::ymd("20290101") 
# Date of LOE for new treatment = 1 January 2031
loe_new_start <- lubridate::ymd("20310101")
# Maturation of new treatment prices by LOE + 1 year, i.e. = 1 January 2032
loe_new_end <- lubridate::ymd("20320101") 

# Effect of LoEs on prices once mature
loe_effect_soc <- 0.7
loe_effect_new <- 0.5

# Calculation of weeks since DoC for LoEs and price maturities
wk_start_soc <- floor((loe_soc_start-doc) / lubridate::dweeks(1))
wk_end_soc <- floor((loe_soc_end-doc) / lubridate::dweeks(1))
wk_start_new <- floor((loe_new_start-doc) / lubridate::dweeks(1))
wk_end_new <- floor((loe_new_end-doc) / lubridate::dweeks(1))

# Price maturity times
wk_mature_soc <- wk_end_soc - wk_start_soc
wk_mature_new <- wk_end_new - wk_start_new

# Create a tibble of price indices of length 2T, then pull out columns as needed
# We only need of length T for now, but need of length 2T for future calculations later
pricetib <- tibble(
  model_time = 1:(2*Ncycles),
  model_year = model_time * cycle_years,
  static = 1,
  geninfl = (1 + infl_cycle)^(model_time - 1),
  loef_soc = pmin(pmax(model_time - wk_start_soc, 0), wk_mature_soc) / wk_mature_soc,
  loef_new = pmin(pmax(model_time - wk_start_new, 0), wk_mature_new) / wk_mature_new,
  dyn_soc = geninfl * (1 - loe_effect_soc * loef_soc),
  dyn_new = geninfl * (1 - loe_effect_new * loef_new)
)

# Price indices required for calculations
prices_oth <- pricetib$geninfl
prices_static <- pricetib$static
prices_dyn_soc <- pricetib$dyn_soc
prices_dyn_new <- pricetib$dyn_new

## ----dynuptake----------------------------------------------------------------
# Time for uptake to occur
uptake_years <- 2

# Uptake vector for non-dynamic uptake
uptake_single <- c(1, rep(0, Ncycles-1))

# Uptake vector for dynamic uptake
uptake_weeks <- round(uptake_years / cycle_years)
share_multi <- c((1:uptake_weeks)/uptake_weeks, rep(1, Ncycles-uptake_weeks))
uptake_multi <- rep(1, Ncycles) * share_multi

## ----static1------------------------------------------------------------------
heemodel

## ----static2------------------------------------------------------------------
# Pull out the payoffs of interest from oncpsm
payoffs <- get_dynfields(
    heemodel = heemodel,
    payoffs = c("cost_daq_new", "cost_daq_soc", "cost_total", "qaly", "life_year"),
    discount = "disc"
    ) |>
    mutate(
      model_years = model_time * cycle_years,
      # Derive costs other than drug acquisition, as at time zero
      cost_nondaq = cost_total - cost_daq_new - cost_daq_soc,
      # ... and at the start of each timestep
      cost_nondaq_rup = cost_total_rup - cost_daq_new_rup - cost_daq_soc_rup
    )

# Create and view dataset for SoC
hemout_soc <- payoffs |> filter(int=="soc")
head(hemout_soc)

# Create and view dataset for new intervention
hemout_new <- payoffs |> filter(int=="new")
head(hemout_new)

## ----scen1--------------------------------------------------------------------
# SOC, costs other than drug acquisition
s1_soc_othcost <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_soc$cost_nondaq_rup,
    prices = prices_oth,
    discrate = nomdisc_cycle
    )

# SOC, drug acquisition costs
s1_soc_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_soc$cost_daq_soc_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# SOC, total costs
s1_soc_cost <- s1_soc_daqcost + s1_soc_othcost

# SOC, QALYs
s1_soc_qaly <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_soc$qaly_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# New intervention, costs other than drug acquisition
s1_new_othcost <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_new$cost_nondaq_rup,
    prices = prices_oth,
    discrate = nomdisc_cycle
    )

# New intervention, drug acquisition costs
s1_new_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_new$cost_daq_new_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# New intervention, total costs
s1_new_cost <- s1_new_daqcost + s1_new_othcost

# New intervention, QALYs
s1_new_qaly <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_new$qaly_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# Incremental cost
s1_icost <- s1_new_cost - s1_soc_cost

summary(s1_icost)

# Incremental QALY
s1_iqaly <- s1_new_qaly - s1_soc_qaly

summary(s1_iqaly)

# ICER
s1_icer <- total(s1_icost) / total(s1_iqaly)
s1_icer

## ----calc_scen2---------------------------------------------------------------
# SOC, costs other than drug acquisition are unchanged
s2_soc_othcost <- s1_soc_othcost

# SoC, drug acquisition costs
s2_soc_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_soc$cost_daq_soc_rup,
    prices = prices_dyn_soc,
    discrate = nomdisc_cycle
    )

# SoC, total costs
s2_soc_cost <- s2_soc_daqcost + s2_soc_othcost

# SoC, QALYs are unchanged
s2_soc_qaly <- s1_soc_qaly

# New intervention, costs other than drug acquisition are unchanged
s2_new_othcost <- s1_new_othcost

# New intervention, drug acquisition costs
s2_new_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_single,
    payoffs = hemout_new$cost_daq_new_rup,
    prices = prices_dyn_new,
    discrate = nomdisc_cycle
    )

# New intervention, total costs
s2_new_cost <- s2_new_daqcost + s2_new_othcost

# New intervention, QALYs are unchanged
s2_new_qaly <- s1_new_qaly

# Incremental cost
s2_icost <- s2_new_cost - s2_soc_cost

summary(s2_icost)

# Incremental QALY
s2_iqaly <- s2_new_qaly - s2_soc_qaly

summary(s2_iqaly)

# ICER
s2_icer <- total(s2_icost) / total(s2_iqaly)
s2_icer

## ----calc_scen3---------------------------------------------------------------
# SOC, costs other than drug acquisition
s3_soc_othcost <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_soc$cost_nondaq_rup,
    prices = prices_oth,
    discrate = nomdisc_cycle
    )

# SoC, drug acquisition costs
s3_soc_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_soc$cost_daq_soc_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# SoC, total costs
s3_soc_cost <- s3_soc_daqcost + s3_soc_othcost

# SoC, QALYs
s3_soc_qaly <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_soc$qaly_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# New intervention, costs other than drug acquisition
s3_new_othcost <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_new$cost_nondaq_rup,
    prices = prices_oth,
    discrate = nomdisc_cycle
    )

# New intervention, drug acquisition costs
s3_new_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_new$cost_daq_new_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# New intervention, total costs
s3_new_cost <- s3_new_daqcost + s3_new_othcost

# New intervention, QALYs
s3_new_qaly <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_new$qaly_rup,
    prices = prices_static,
    discrate = disc_cycle
    )

# Incremental costs
s3_icost <- s3_new_cost - s3_soc_cost

summary(s3_icost)

# Incremental QALYs
s3_iqaly <- s3_new_qaly - s3_soc_qaly

summary(s3_iqaly)

# ICER
s3_icer <- total(s3_icost) / total(s3_iqaly)
s3_icer

## ----calc_scen4---------------------------------------------------------------
# SOC, costs other than drug acquisition are unchanged
s4_soc_othcost <- s3_soc_othcost

# SoC, drug acquisition costs
s4_soc_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_soc$cost_daq_soc_rup,
    prices = prices_dyn_soc,
    discrate = nomdisc_cycle
    )

# SoC, total costs
s4_soc_cost <- s4_soc_daqcost + s4_soc_othcost

# SoC, QALYs are unchanged
s4_soc_qaly <- s3_soc_qaly

# New intervention, costs other than drug acquisition are unchanged
s4_new_othcost <- s3_new_othcost

# New intervention, drug acquisition costs
s4_new_daqcost <- dynamicpv::dynpv(
    uptakes = uptake_multi,
    payoffs = hemout_new$cost_daq_new_rup,
    prices = prices_dyn_new,
    discrate = nomdisc_cycle
    )

# New intervention, total costs
s4_new_cost <- s4_new_daqcost + s4_new_othcost

# New intervention, QALYs
s4_new_qaly <- s3_new_qaly

# Incremental costs
s4_icost <- s4_new_cost - s4_soc_cost

summary(s4_icost)

# Incremental QALYs
s4_iqaly <- s4_new_qaly - s4_soc_qaly

summary(s4_iqaly)

# ICER
s4_icer <- total(s4_icost) / total(s4_iqaly)
s4_icer

## ----future_calc--------------------------------------------------------------
# Times at which to plot ICER
gtimes <- round((0:(2*thoz))/cycle_years/2)
# SOC drug acquisition costs
gc_soc_daq <- gtimes |>
    purrr::map_vec(\(l) mean(futurepv(
        tzero = l,
        payoffs = hemout_soc$cost_daq_soc_rup,
        prices = prices_dyn_soc,
        discrate = nomdisc_cycle
    ))
)
# SOC other costs
gc_soc_oth <- gtimes |>
    purrr::map_vec(\(l) mean(futurepv(
        tzero = l,
        payoffs = hemout_soc$cost_nondaq_rup,
        prices = prices_oth,
        discrate = nomdisc_cycle
    ))
)
# New drug acquisition costs
gc_new_daq <- gtimes |>
    purrr::map_vec(\(l) mean(futurepv(
        tzero = l,
        payoffs = hemout_new$cost_daq_new_rup,
        prices = prices_dyn_new,
        discrate = nomdisc_cycle
    ))
)
# New other costs
gc_new_oth <- gtimes |>
    purrr::map_vec(\(l) mean(futurepv(
        tzero = l,
        payoffs = hemout_new$cost_nondaq_rup,
        prices = prices_oth,
        discrate = nomdisc_cycle
    ))
)
# Combine in a tibble
ds <- tibble(
    # Time in weeks and years
    time_weeks = gtimes,
    time_years = time_weeks * cycle_years,
    # Evaluation date
    evaldate = doc + time_weeks * 7,
    # Price/inflation index
    pinfl = prices_oth[gtimes + 1],
    # Total costs for each intervention
    totcost_new = gc_new_daq + gc_new_oth,
    totcost_soc = gc_soc_daq + gc_soc_oth,
    # Scenario 1/2 QALYs
    qaly_soc = total(s1_soc_qaly),
    qaly_new = total(s1_new_qaly)
  ) |>
  mutate(
    # Incremental cost and QALYs
    icost = totcost_new-totcost_soc,
    iqaly = qaly_new-qaly_soc,
    # Nominal ICER, and real (inflation adjusted) ICER
    Nominal = icost/iqaly,
    Real = Nominal / pinfl
  ) |>
  # Pivot to long so can be used in a graphic
  pivot_longer(
    cols = c("Nominal", "Real"),
    names_to = "Type",
    values_to = "ICER"
  )

## ----calc_future--------------------------------------------------------------
# Plot real and nominal present value over time
ggplot(ds,
  aes(x = evaldate, y = ICER, color=Type)) +
  geom_line() +
  labs(x = "Evaluation date") +
  geom_hline(yintercept = ds$ICER[1], linetype='dotted') +
  geom_vline(xintercept = loe_new_start, linetype='dashed') +
  geom_vline(xintercept = loe_soc_start, linetype='dashed') +
  scale_y_continuous(
    labels = scales::comma,
    limits=c(0, 150000)
    )