In n8thangreen/LTBIscreeningproject: LTBI screening cost-effectiveness analysis

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

options(width = 240)

This document demonstrates running the screening, diagnosis and treatment decision tree model.

First, we need to attach the packages we'll need. treeSimR is written by us too so this will need to be installed from GitHub.

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

Load in data

data("intervention_constants", package = "LTBIscreeningproject")
data("cost_effectivness_params", package = "LTBIscreeningproject")
data("scenario_parameters", package = "LTBIscreeningproject")
load("model_input_cohort.RData")

Intervention specific individual-level data

According to the parameter values we specified in interv we first need to prepare the cohort data set to only include those individuals we're interested in analysing.

Only use a single year cohort

cohort <- dplyr::filter(IMPUTED_sample,
                        issdt_year == interv$year_cohort)

EWNI stay long enough to be screened

cohort <- dplyr::filter(cohort,
                        date_exit_uk1_issdt.years >= interv$min_screen_length_of_stay)

Only keep those inidiviuals that are screened before some other event

if (interv$screen_with_delay) {
  cohort <- dplyr::filter(cohort,
                          screen == 1)
}

Do we include student in the screening programme?

if (interv$no_students) {
  cohort <- dplyr::filter(cohort,
                          visatype2 != "Students")}

Include individuals from 'higher' incidence countries only

cohort <- dplyr::filter(cohort,
                        who_prev_cat_Pareek2011 %in% interv$incidence_grps_screen)

Count numbers of TB cases

n.exit_tb <-
  cohort %>%
  dplyr::filter(exituk_tb) %>%
  dplyr::count()

n.uk_tb <-
  cohort %>%
  dplyr::filter(uk_tb) %>%
  dplyr::count()


num_all_tb_cost <-
  if (interv$ENDPOINT_cost == "exit uk") {
    n.uk_tb
  } else if (interv$ENDPOINT_cost == "death") {
    n.uk_tb + n.exit_tb}

num_all_tb_QALY <-
  if (interv$ENDPOINT_QALY == "exit uk") {
    n.uk_tb
  } else if (interv$ENDPOINT_QALY == "death") {
    n.uk_tb + n.exit_tb}

saveRDS(cohort, file = "cohort.Rds")

Prepare decision tree

Load input files. We have defined the decision tree using YAML syntax. We have a separate trees for cost and health. They have the same structure and probabilities but different values on the branches. We also create our own class.

osNode.cost.fileName <- system.file("data", "LTBI_dtree-cost-symptoms.yaml",
                                    package = "LTBIscreeningproject")

osNode.health.fileName <- system.file("data", "LTBI_dtree-QALYloss-symptoms.yaml",
                                      package = "LTBIscreeningproject")

costeff.cost <- treeSimR::costeffectiveness_tree(yaml_tree = osNode.cost.fileName)
osNode.cost <- costeff.cost$osNode

costeff.health <- treeSimR::costeffectiveness_tree(yaml_tree = osNode.health.fileName)
osNode.health <- costeff.health$osNode

costeff.cost <- treeSimR::costeffectiveness_tree(yaml_tree = "LTBI_dtree-cost-symptoms.yaml")
print(costeff.cost)

Create look-up tables using the cohort data to give proportion in each incidence in country of origin and the probability of LTBI for each of these.

who_levels <- c("(0,50]", "(50,150]", "(150,250]", "(250,350]", "(350,1e+05]")

p_incid_grp <- miscUtilities::prop_table(cohort$who_prev_cat_Pareek2011)

pLatentTB.who <- data.frame(who_prev_cat_Pareek2011 = names(p_incid_grp),
                            LTBI = c(0.03, 0.13, 0.2, 0.3, 0.3))

Then insert these probabilities in to the decision tree

for (i in seq_along(who_levels)) {

  osNode.cost$Set(p = p_incid_grp[i],
                  filterFun = function(x) x$name == who_levels[i])

  osNode.health$Set(p = p_incid_grp[i],
                    filterFun = function(x) x$name == who_levels[i])
}

for (i in who_levels) {

  pLTBI <- subset(pLatentTB.who,
                  who_prev_cat_Pareek2011 == i,
                  select = LTBI)

  osNode.cost$Set(p = pLTBI,
                  filterFun = function(x) x$pathString == miscUtilities::pastef("LTBI screening cost", i, "LTBI"))
  osNode.health$Set(p = pLTBI,
                    filterFun = function(x) x$pathString == miscUtilities::pastef("LTBI screening QALY loss", i, "LTBI"))

  osNode.cost$Set(p = 1 - pLTBI,
                  filterFun = function(x) x$pathString == miscUtilities::pastef("LTBI screening cost", i, "non-LTBI"))
  osNode.health$Set(p = 1 - pLTBI,
                    filterFun = function(x) x$pathString == miscUtilities::pastef("LTBI screening QALY loss", i, "non-LTBI"))
}

Decision tree output

Now we can run the decision tree model using the scenario parameter values, active TB cases, and cost and health decision tree objects.

A single scenario:

decision_tree_cluster(parameters = scenario_parameters[[1]],
                      N.mc = interv$N.mc,
                      n.uk_tb = 10, #as.numeric(n.uk_tb),
                      n.exit_tb = 10, #as.numeric(n.exit_tb),
                      cost_dectree = "osNode_cost_2009.Rds",
                      health_dectree = "osNode_health_2009.Rds")

Multiple scenarios:

lapply(scenario_parameters[1:3],
       decision_tree_cluster,
       N.mc = 1,
       n.uk_tb = 10, #as.numeric(n.uk_tb),
       n.exit_tb = 10, #as.numeric(n.exit_tb),
       cost_dectree = "osNode_cost_2009.Rds",
       health_dectree = "osNode_health_2009.Rds")

lapply(scenario_parameters[1:3],
       decision_tree_cluster,
       N.mc = 10,
       n.uk_tb = 10, #as.numeric(n.uk_tb),
       n.exit_tb = 10, #as.numeric(n.exit_tb),
       cost_dectree = "osNode_cost_2009.Rds",
       health_dectree = "osNode_health_2009.Rds") %>% 
  saveRDS(file = "decision_tree_output.Rds")