In n8thangreen/LTBIscreeningproject: LTBI screening cost-effectiveness analysis

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

Introduction

This document demonstrates how we calculate the cost and QALY statistics.

Analysis preparation

Load packages

library(treeSimR)
library(LTBIscreeningproject)
library(dplyr)
library(data.tree)
library(purrr)

We first need to read-in the decision tree output data and the individual cohort data.

dectree_res <- readRDS("decision_tree_output.Rds")
cohort <- readRDS("cohort.Rds")

# data("cost_effectivness_params", package = "LTBIscreeningproject")

It is sometimes arkward to use the decision tree output in the format with the scenario number as the top level. It may be more useful to have each output type at the top level instead, so we transpose the data.

From:

lapply(dectree_res, FUN = function(x) lapply(x, head))

to:

# convert from scenario-wise to parameter-wise format
scenario_res <-
  dectree_res %>%
  purrr::transpose()

lapply(scenario_res, function(x) lapply(x, head))

We need to extract some specific values and initialise the variables for the cost-effectiveness statistics.

Number of TB cases in status-quo scenario.

n_uk_tb <- dectree_res$`1`$call$n.uk_tb
n_exit_tb <- dectree_res$`1`$call$n.exit_tb
n_all_tb <- n_exit_tb + n_uk_tb

Numbers of individuals who avoid getting TB due to screening.

avoid_all_tb <-
  scenario_res$n_tb_screen_all %>%
  purrr::map(function(x) dplyr::filter(x, status == "disease-free")) %>%
  map(function(x) x$n)

avoid_uk_tb  <-
  scenario_res$n_tb_screen_uk %>%
  purrr::map(function(x) dplyr::filter(x, status == "disease-free")) %>%
  map(function(x) x$n)

avoid_tb <-
  mapply(function(x,y) cbind(x, y),
         avoid_all_tb,
         avoid_uk_tb,
         SIMPLIFY = FALSE) %>% 
  map(`colnames<-`, c('all','uk'))

rm(avoid_all_tb,
   avoid_uk_tb)

Define some variables we'll need later.

n.scenarios <- length(dectree_res)
N.mc <- dectree_res$`1`$call$N.mc

interv_cost <- vector(length = N.mc, mode = "list")
interv_QALY <- vector(length = N.mc, mode = "list")
stats_scenario <- vector(length = n.scenarios, mode = "list")

pop_year <-  1000 ##TODO: get real value

We also calculate the expected statistics for reproducability and comparison with the randomly generated values.

mean_cost.aTB_TxDx <-
  unit_cost$aTB_TxDx %>% 
  means_distributions() %>%
  sum()

mean_num_sec_inf <-
  NUM_SECONDARY_INF %>%
  means_distributions() %>%
  unlist()

Extract cost-effectiveness variables from main individual-level data set and create some more expected statistics using these.

costeff_cohort <-
  cohort %>%
  select(cfr,
         QALY_statusquo,
         QALY_diseasefree,
         QALY_cured,
         QALY_fatality,
         uk_notif_discounts,
         all_notif_discounts,
         uk_secondary_inf_discounts,
         all_secondary_inf_discounts) %>% 
  mutate(E_cost_sec_inf = mean_num_sec_inf * mean_cost.aTB_TxDx * all_secondary_inf_discounts,
         E_cost_statusquo = (all_notif_discounts * mean_cost.aTB_TxDx) + E_cost_sec_inf,
         E_QALY_statusquo = (cfr * QALY_fatality) + ((1 - cfr) * QALY_cured)) %>%
  stats::na.omit()

Calculate cost-effectiveness statistics for each scenario

We iterate over all scenarios and all Monte-Carlo samples. The scenario_cost and scenario_QALY functions calculate the respective values for each. We then transpose the returned values from these functions and calculate the final statistics in a single list.

interv_scenario_cost <- partial(scenario_cost,
                                endpoint = interv$ENDPOINT_cost,
                                unit_cost.aTB_TxDx = unit_cost$aTB_TxDx,
                                num_2nd_inf = NUM_SECONDARY_INF,
                                costeff_cohort = costeff_cohort,
                                total_tb = n_all_tb)

interv_scenario_QALY <- partial(scenario_QALY,
                                total = n_all_tb,
                                costeff_data = costeff_cohort)

for (s in seq_len(n.scenarios)) {
  for (i in seq_len(N.mc)) {

    # set.seed(12345)

    num_avoided <- avoid_tb[[s]][i, ]

    interv_cost[[i]] <- interv_scenario_cost(num_avoided)
    interv_QALY[[i]] <- interv_scenario_QALY(num_avoided)
  }

  stats_scenario[[s]] <- costeff_stats(scenario_dat = dectree_res[[s]],
                                       interv_QALY = interv_QALY,
                                       interv_cost = interv_cost,
                                       pop_year = pop_year)
}

aTB_CE_stats <- stats_scenario %>% purrr::transpose()

aTB_CE_stats

save(aTB_CE_stats,
     file = "aTB_CE_stats.RData")