In n8thangreen/LTBIdiagTST: Latent TB infection cost-effectiveness analysis

knitr::opts_chunk$set(echo = TRUE)

We will replicate the Excel TB TST-T-SPOT.TB cost-effectiveness model. (File name latesttree.xls.)

Set-up

library(readr)
library(dplyr)
library(tibble)
library(reshape2)
library(assertthat)
library(CEdecisiontree)
library(purrr)

Load in the data.

load(here::here("data", "params.RData"))
load(here::here("data", "trees.RData"))

This included the tree structure variables, the cost and probability arrays, the mapping arrays that inform which nodes have which label with have what cost.

To demonstrate, let us look at the TST and QFT scenario. The decision tree is defined in terms of parents and children in a list.

head(TST_IGRA_fwd_tree, 5)

In order to assign values to the tree we first transform this to a (sparse) matrix format.

probs <- child_list_to_transmat(TST_IGRA_fwd_tree)
print(as_tibble(probs), n = 100)

empty_transmat <-
  as_tibble(matrix(NA_real_,
                   nrow = nrow(probs), ncol = ncol(probs)),
            .name_repair = "minimal")

Next, we specify the labels for each of the edge (or correspondingly to node). We need to do this separately for the probabilities and costs.

pname_from_to <- TST_QFT_fwd_pname_from_to
pname_from_to

cname_from_to <- TST_QFT_fwd_cname_from_to
cname_from_to

Insert probabilities into decision tree

Now that we've set-up the framework for the decision tree, we can assign the input data to it. The probability data is in the form of a list (this is useful for when we want to sample from a distribution later). Lets transform to an array.

label_probs_long <-
  as_tibble(label_probs) %>%
  melt(value.name = "prob",
       variable.name = "name")

Insert the appropriate probabilities by converting the transition matrix to long format, matching branches to labels, matching labels to probabilities and then filling in missing probabilities so that pairs of branches sum to one.

probs_new <-
  probs %>%
  transmat_to_long() %>%
  match_branch_to_label(pname_from_to) %>%
  match_branchlabel_to_prob(label_probs_long) %>%
  fill_complementary_probs()

probs_new

Finally, we insert these new probabilities in to the decision tree.

probs <- insert_to_probmat(dat = probs_new,
                           mat = empty_transmat)
head(probs)

Insert costs into decision tree

We essentially do the same thing now for costs.

label_cost_long <-
  as_tibble(label_costs) %>%
  melt(value.name = "cost",
       variable.name = "name")

head(label_cost_long)

Join the cost names and their associated branches in to a single array.

costs_names <-
  merge(cname_from_to, label_cost_long,
        by = "name", all.x = TRUE) %>%
  mutate(from = as.numeric(as.character(from)),
         to = as.numeric(as.character(to)))
costs_names

Finally, we insert these costs in to the decision tree.

costs <- insert_to_costmat(dat = costs_names,
                           mat = empty_transmat)
head(costs)

Run model

See the CEdecisiontree package for how to use the dectree_expected_value() function. Here we provide the matrix format arguments.

TST_model <-
  define_model(transmat = list(prob = probs,
                               vals = costs))

res <-
  dectree_expected_values(TST_model)

res[1] + label_costs$`TB special nurse visit`

# 196.9457 in xls

Create a single cost-effectiveness decision tree data frame.

all_long <-
  merge(costs_names, probs_new,
        all = TRUE, by = c("from", "to"),
        suffixes = c(".cost", ".prob")) %>%
  rename(vals = cost)

head(all_long)

Probabilistic sensitivity analysis

We first define the decision tree. The difference to previous trees is that we now use the list-column feature to define distributions rather than point values.

Combine the distributions with the original tree specification data frame and remove the point values.

input_SA <- 
  all_long %>% 
  select(-prob, -vals) %>% 
  dplyr::left_join(label_probs_distns, by = "name.prob") %>%
  dplyr::left_join(label_costs_distns, by = "name.cost") %>%
  as_tibble()

input_SA

We can now loop over this tree and generate samples of values for the given distributions. Lets use the sample_distributions() function for this.

tree_dat_sa <- list()

for (i in 1:50) {

  tree_dat_sa[[i]] <-
    define_model(
      dat_long = 
        data.frame(
          from = input_SA$from,
          to   = input_SA$to,
          prob = lapply(input_SA$prob, sample_distributions) %>% unlist(),
          vals = lapply(input_SA$vals, sample_distributions) %>% unlist()), fill_probs = TRUE)
}

This results in a list of trees. Now it is straightforward to map over each of these trees to obtain the total expected values

res <-
  suppressMessages(
    map(tree_dat_sa, .f = dectree_expected_values))
head(res)

hist(map_dbl(res, 1), breaks = 20)

Pathway joint probabilities

We calculate the terminal state total probabilities for each LTBI status to use as inputs for the Markov model.

Define the model.

TST_model_long <- 
  define_model(dat_long = all_long)

For each MM starting state.

LTBI_complete_Tx <- 
  branch_joint_probs(
    TST_model_long,
    nodes = 12) %>% 
  map_dbl(prod) %>% 
  sum()

LTBI_incomplete_Tx <- 
  branch_joint_probs(
    TST_model_long,
    nodes = c(10,48)) %>% 
  map_dbl(prod) %>% 
  sum()

LTBI_no_Tx <- 
  branch_joint_probs(
    TST_model_long,
    nodes = c(14,16,18,20,22,24)) %>% 
  map_dbl(prod) %>% 
  sum()

no_LTBI <- 
  branch_joint_probs(
    TST_model_long,
    nodes = c(33,35,37,39,41,43,45,47,49)) %>% 
  map_dbl(prod) %>% 
  sum()

Check that they sum to one.

all.equal(1 - (LTBI_complete_Tx + LTBI_incomplete_Tx + LTBI_no_Tx), no_LTBI)

Now we can do the same for the SA data to give a distribution of terminal state subpopulations.

terminal_pop <- function(df) {

  res <- list()

  res$completeTx <- 
    branch_joint_probs(
      df,
      nodes = 12) %>% 
    map_dbl(prod) %>% 
    sum()

  res$incompleteTx <- 
    branch_joint_probs(
      df,
      nodes = c(10,48)) %>% 
    map_dbl(prod) %>% 
    sum()

  res$noTx <- 
    branch_joint_probs(
      df,
      nodes = c(14,16,18,20,22,24)) %>% 
    map_dbl(prod) %>% 
    sum()

  res$noLTBI <- 
    branch_joint_probs(
      df,
      nodes = c(33,35,37,39,41,43,45,47,49)) %>% 
    map_dbl(prod) %>% 
    sum()

  return(res)
}

Check

terminal_pop(all_long)

Now we can map over all of the sample trees. Note we restrict to the non-NA branches. This is something to fix properly in the future.

init_states <-
  map_df(tree_dat_sa,
         .f = function(x) terminal_pop(x[1:48, ]))

save(init_states, file = "data/init_states.RData")

Plot the results with the original scenario indicated by the vertical red line.

par(mfrow = c(2,2))

hist(init_states$LTBI_complete_Tx, breaks = 25)
abline(v = LTBI_complete_Tx, col = "red")

hist(init_states$LTBI_incomplete_Tx, breaks = 25)
abline(v = LTBI_incomplete_Tx, col = "red")

hist(init_states$LTBI_no_Tx, breaks = 25)
abline(v = LTBI_no_Tx, col = "red")

hist(init_states$no_LTBI, breaks = 25)
abline(v = no_LTBI, col = "red")