Cost-effectiveness calculation

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))
#> $`1`
#> $`1`$mc_cost
#> [1] 36.31012 36.82432 36.33503 36.25995 37.00919 36.80207
#> 
#> $`1`$mc_health
#> [1] 5.585985e-06 5.370821e-06 4.646076e-06 5.355858e-06 5.362683e-06
#> [6] 5.718115e-06
#> 
#> $`1`$n_tb_screen_all
#>   sim       status  n
#> 1   1 disease-free  3
#> 2   1           tb 17
#> 3   2 disease-free  2
#> 4   2           tb 18
#> 5   3 disease-free  2
#> 6   3           tb 18
#> 
#> $`1`$n_tb_screen_uk
#>   sim       status  n
#> 1   1 disease-free  2
#> 2   1           tb  8
#> 3   2 disease-free  1
#> 4   2           tb  9
#> 5   3 disease-free  0
#> 6   3           tb 10
#> 
#> $`1`$p_LTBI_to_effectiveTx
#> (350,1e+05] 
#>    0.108864 
#> 
#> $`1`$subset_pop
#>    LTBI_pre tests  positive    startTx completeTx      cured LTBI_post
#> 1 0.2868171   0.6 0.1488349 0.08930096 0.05358058 0.03122406 0.2555931
#> 
#> $`1`$call
#> FUN(parameters = X[[i]], N.mc = 10, n.uk_tb = 10, n.exit_tb = 10, 
#>     cost_dectree = "osNode_cost_2009.Rds")
#> 
#> $`1`$N.mc
#> [1] 10
#> 
#> 
#> $`2`
#> $`2`$mc_cost
#> [1] 44.46451 45.99899 44.63913 43.93595 44.20433 44.21356
#> 
#> $`2`$mc_health
#> [1] 7.521093e-06 6.090351e-06 6.261289e-06 6.969095e-06 7.260007e-06
#> [6] 7.248490e-06
#> 
#> $`2`$n_tb_screen_all
#>   sim       status  n
#> 1   1 disease-free  3
#> 2   1           tb 17
#> 3   2 disease-free  3
#> 4   2           tb 17
#> 5   3 disease-free  3
#> 6   3           tb 17
#> 
#> $`2`$n_tb_screen_uk
#>   sim       status n
#> 1   1 disease-free 2
#> 2   1           tb 8
#> 3   2 disease-free 2
#> 4   2           tb 8
#> 5   3 disease-free 2
#> 6   3           tb 8
#> 
#> $`2`$p_LTBI_to_effectiveTx
#> (350,1e+05] 
#>    0.145152 
#> 
#> $`2`$subset_pop
#>    LTBI_pre tests  positive   startTx completeTx      cured LTBI_post
#> 1 0.2868171   0.6 0.1488349 0.1190679 0.07144077 0.04163208 0.2451851
#> 
#> $`2`$call
#> FUN(parameters = X[[i]], N.mc = 10, n.uk_tb = 10, n.exit_tb = 10, 
#>     cost_dectree = "osNode_cost_2009.Rds")
#> 
#> $`2`$N.mc
#> [1] 10
#> 
#> 
#> $`3`
#> $`3`$mc_cost
#> [1] 51.95845 53.77997 53.33703 52.46245 52.78314 52.13228
#> 
#> $`3`$mc_health
#> [1] 8.524595e-06 7.996882e-06 9.289761e-06 9.261866e-06 8.941128e-06
#> [6] 7.988450e-06
#> 
#> $`3`$n_tb_screen_all
#>   sim       status  n
#> 1   1 disease-free  3
#> 2   1           tb 17
#> 3   2 disease-free  4
#> 4   2           tb 16
#> 5   3 disease-free  3
#> 6   3           tb 17
#> 
#> $`3`$n_tb_screen_uk
#>   sim       status n
#> 1   1 disease-free 1
#> 2   1           tb 9
#> 3   2 disease-free 1
#> 4   2           tb 9
#> 5   3 disease-free 1
#> 6   3           tb 9
#> 
#> $`3`$p_LTBI_to_effectiveTx
#> (350,1e+05] 
#>     0.18144 
#> 
#> $`3`$subset_pop
#>    LTBI_pre tests  positive   startTx completeTx     cured LTBI_post
#> 1 0.2868171   0.6 0.1488349 0.1488349 0.08930096 0.0520401  0.234777
#> 
#> $`3`$call
#> FUN(parameters = X[[i]], N.mc = 10, n.uk_tb = 10, n.exit_tb = 10, 
#>     cost_dectree = "osNode_cost_2009.Rds")
#> 
#> $`3`$N.mc
#> [1] 10

to:

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

lapply(scenario_res, function(x) lapply(x, head))
#> $mc_cost
#> $mc_cost$`1`
#> [1] 36.31012 36.82432 36.33503 36.25995 37.00919 36.80207
#> 
#> $mc_cost$`2`
#> [1] 44.46451 45.99899 44.63913 43.93595 44.20433 44.21356
#> 
#> $mc_cost$`3`
#> [1] 51.95845 53.77997 53.33703 52.46245 52.78314 52.13228
#> 
#> 
#> $mc_health
#> $mc_health$`1`
#> [1] 5.585985e-06 5.370821e-06 4.646076e-06 5.355858e-06 5.362683e-06
#> [6] 5.718115e-06
#> 
#> $mc_health$`2`
#> [1] 7.521093e-06 6.090351e-06 6.261289e-06 6.969095e-06 7.260007e-06
#> [6] 7.248490e-06
#> 
#> $mc_health$`3`
#> [1] 8.524595e-06 7.996882e-06 9.289761e-06 9.261866e-06 8.941128e-06
#> [6] 7.988450e-06
#> 
#> 
#> $n_tb_screen_all
#> $n_tb_screen_all$`1`
#>   sim       status  n
#> 1   1 disease-free  3
#> 2   1           tb 17
#> 3   2 disease-free  2
#> 4   2           tb 18
#> 5   3 disease-free  2
#> 6   3           tb 18
#> 
#> $n_tb_screen_all$`2`
#>   sim       status  n
#> 1   1 disease-free  3
#> 2   1           tb 17
#> 3   2 disease-free  3
#> 4   2           tb 17
#> 5   3 disease-free  3
#> 6   3           tb 17
#> 
#> $n_tb_screen_all$`3`
#>   sim       status  n
#> 1   1 disease-free  3
#> 2   1           tb 17
#> 3   2 disease-free  4
#> 4   2           tb 16
#> 5   3 disease-free  3
#> 6   3           tb 17
#> 
#> 
#> $n_tb_screen_uk
#> $n_tb_screen_uk$`1`
#>   sim       status  n
#> 1   1 disease-free  2
#> 2   1           tb  8
#> 3   2 disease-free  1
#> 4   2           tb  9
#> 5   3 disease-free  0
#> 6   3           tb 10
#> 
#> $n_tb_screen_uk$`2`
#>   sim       status n
#> 1   1 disease-free 2
#> 2   1           tb 8
#> 3   2 disease-free 2
#> 4   2           tb 8
#> 5   3 disease-free 2
#> 6   3           tb 8
#> 
#> $n_tb_screen_uk$`3`
#>   sim       status n
#> 1   1 disease-free 1
#> 2   1           tb 9
#> 3   2 disease-free 1
#> 4   2           tb 9
#> 5   3 disease-free 1
#> 6   3           tb 9
#> 
#> 
#> $p_LTBI_to_effectiveTx
#> $p_LTBI_to_effectiveTx$`1`
#> (350,1e+05] 
#>    0.108864 
#> 
#> $p_LTBI_to_effectiveTx$`2`
#> (350,1e+05] 
#>    0.145152 
#> 
#> $p_LTBI_to_effectiveTx$`3`
#> (350,1e+05] 
#>     0.18144 
#> 
#> 
#> $subset_pop
#> $subset_pop$`1`
#>    LTBI_pre tests  positive    startTx completeTx      cured LTBI_post
#> 1 0.2868171   0.6 0.1488349 0.08930096 0.05358058 0.03122406 0.2555931
#> 
#> $subset_pop$`2`
#>    LTBI_pre tests  positive   startTx completeTx      cured LTBI_post
#> 1 0.2868171   0.6 0.1488349 0.1190679 0.07144077 0.04163208 0.2451851
#> 
#> $subset_pop$`3`
#>    LTBI_pre tests  positive   startTx completeTx     cured LTBI_post
#> 1 0.2868171   0.6 0.1488349 0.1488349 0.08930096 0.0520401  0.234777
#> 
#> 
#> $call
#> $call$`1`
#> FUN(parameters = X[[i]], N.mc = 10, n.uk_tb = 10, n.exit_tb = 10, 
#>     cost_dectree = "osNode_cost_2009.Rds")
#> 
#> $call$`2`
#> FUN(parameters = X[[i]], N.mc = 10, n.uk_tb = 10, n.exit_tb = 10, 
#>     cost_dectree = "osNode_cost_2009.Rds")
#> 
#> $call$`3`
#> FUN(parameters = X[[i]], N.mc = 10, n.uk_tb = 10, n.exit_tb = 10, 
#>     cost_dectree = "osNode_cost_2009.Rds")
#> 
#> 
#> $N.mc
#> $N.mc$`1`
#> [1] 10
#> 
#> $N.mc$`2`
#> [1] 10
#> 
#> $N.mc$`3`
#> [1] 10

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
#> $QALY.statusquo
#> $QALY.statusquo[[1]]
#>  [1] 14031.38 14031.38 14031.38 14031.38 14031.38 14031.38 14031.38
#>  [8] 14031.38 14031.38 14031.38
#> 
#> $QALY.statusquo[[2]]
#>  [1] 14031.38 14031.38 14031.38 14031.38 14031.38 14031.38 14031.38
#>  [8] 14031.38 14031.38 14031.38
#> 
#> $QALY.statusquo[[3]]
#>  [1] 14031.38 14031.38 14031.38 14031.38 14031.38 14031.38 14031.38
#>  [8] 14031.38 14031.38 14031.38
#> 
#> 
#> $QALY.screened
#> $QALY.screened[[1]]
#>  [1] 14032.31 14031.98 14031.98 14031.98 14032.00 14032.00 14032.62
#>  [8] 14032.28 14032.31 14031.98
#> 
#> $QALY.screened[[2]]
#>  [1] 14032.31 14032.29 14032.29 14033.23 14031.66 14032.31 14032.92
#>  [8] 14031.69 14031.38 14032.62
#> 
#> $QALY.screened[[3]]
#>  [1] 14032.29 14032.63 14032.31 14032.60 14031.98 14033.86 14031.38
#>  [8] 14032.91 14032.58 14032.00
#> 
#> 
#> $E_cost_screened
#> $E_cost_screened[[1]]
#> [1] 3903219
#> 
#> $E_cost_screened[[2]]
#> [1] 3994908
#> 
#> $E_cost_screened[[3]]
#> [1] 4397815
#> 
#> 
#> $cost.screened_person
#> $cost.screened_person[[1]]
#>  [1] 2892.155 3694.782 5375.227 2237.092 4736.823 3761.111 3073.991
#>  [8] 3802.369 5141.063 4317.578
#> 
#> $cost.screened_person[[2]]
#>  [1] 5907.969 3591.027 3023.731 3648.177 4356.866 3169.036 2444.561
#>  [8] 4059.441 5727.177 4021.097
#> 
#> $cost.screened_person[[3]]
#>  [1] 4167.226 6662.840 3822.676 4525.689 5375.961 4764.232 4329.864
#>  [8] 3010.001 4440.565 2879.095
#> 
#> 
#> $cost.statusquo_person
#> $cost.statusquo_person[[1]]
#>  [1] 2900.635 3700.551 5388.829 2241.034 4748.710 3770.705 3083.310
#>  [8] 3812.376 5149.930 4326.714
#> 
#> $cost.statusquo_person[[2]]
#>  [1] 5929.015 3602.837 3031.512 3672.614 4363.494 3176.638 2452.708
#>  [8] 4065.407 5727.177 4045.260
#> 
#> $cost.statusquo_person[[3]]
#>  [1] 4181.493 6690.999 3831.732 4546.937 5386.166 4811.517 4329.864
#>  [8] 3023.645 4462.383 2888.022
#> 
#> 
#> $cost_incur
#> $cost_incur[[1]]
#>  [1]  -8480.538  -5769.061 -13602.230  -3941.800 -11886.553  -9594.122
#>  [7]  -9318.991 -10007.683  -8867.694  -9136.480
#> 
#> $cost_incur[[2]]
#>  [1] -21045.602 -11809.604  -7781.309 -24437.684  -6627.897  -7602.122
#>  [7]  -8146.484  -5966.299      0.000 -24162.497
#> 
#> $cost_incur[[3]]
#>  [1] -14266.718 -28158.966  -9056.345 -21248.178 -10204.341 -47284.430
#>  [7]      0.000 -13643.529 -21818.627  -8927.014
#> 
#> 
#> $cost.statusquo
#> $cost.statusquo[[1]]
#>  [1] 2900635 3700551 5388829 2241034 4748710 3770705 3083310 3812376
#>  [9] 5149930 4326714
#> 
#> $cost.statusquo[[2]]
#>  [1] 5929015 3602837 3031512 3672614 4363494 3176638 2452708 4065407
#>  [9] 5727177 4045260
#> 
#> $cost.statusquo[[3]]
#>  [1] 4181493 6690999 3831732 4546937 5386166 4811517 4329864 3023645
#>  [9] 4462383 2888022
#> 
#> 
#> $cost.screened
#> $cost.screened[[1]]
#>  [1] 2892155 3694782 5375227 2237092 4736823 3761111 3073991 3802369
#>  [9] 5141063 4317578
#> 
#> $cost.screened[[2]]
#>  [1] 5907969 3591027 3023731 3648177 4356866 3169036 2444561 4059441
#>  [9] 5727177 4021097
#> 
#> $cost.screened[[3]]
#>  [1] 4167226 6662840 3822676 4525689 5375961 4764232 4329864 3010001
#>  [9] 4440565 2879095
#> 
#> 
#> $E_QALY_screened
#> $E_QALY_screened[[1]]
#> [1] 14032.14
#> 
#> $E_QALY_screened[[2]]
#> [1] 14032.27
#> 
#> $E_QALY_screened[[3]]
#> [1] 14032.45
#> 
#> 
#> $QALY.screened_person
#> $QALY.screened_person[[1]]
#>  [1] 14.03231 14.03198 14.03198 14.03198 14.03200 14.03200 14.03262
#>  [8] 14.03228 14.03231 14.03198
#> 
#> $QALY.screened_person[[2]]
#>  [1] 14.03231 14.03229 14.03229 14.03323 14.03166 14.03231 14.03292
#>  [8] 14.03169 14.03138 14.03262
#> 
#> $QALY.screened_person[[3]]
#>  [1] 14.03229 14.03263 14.03231 14.03260 14.03198 14.03386 14.03138
#>  [8] 14.03291 14.03258 14.03200
#> 
#> 
#> $QALY.statusquo_person
#> $QALY.statusquo_person[[1]]
#>  [1] 14.03138 14.03138 14.03138 14.03138 14.03138 14.03138 14.03138
#>  [8] 14.03138 14.03138 14.03138
#> 
#> $QALY.statusquo_person[[2]]
#>  [1] 14.03138 14.03138 14.03138 14.03138 14.03138 14.03138 14.03138
#>  [8] 14.03138 14.03138 14.03138
#> 
#> $QALY.statusquo_person[[3]]
#>  [1] 14.03138 14.03138 14.03138 14.03138 14.03138 14.03138 14.03138
#>  [8] 14.03138 14.03138 14.03138
#> 
#> 
#> $QALYgain
#> $QALYgain[[1]]
#>  [1] 0.9396 0.6012 0.6012 0.6012 0.6264 0.6264 1.2456 0.9072 0.9396 0.6012
#> 
#> $QALYgain[[2]]
#>  [1] 0.9396 0.9144 0.9144 1.8540 0.2880 0.9396 1.5408 0.3132 0.0000 1.2456
#> 
#> $QALYgain[[3]]
#>  [1] 0.9144 1.2528 0.9396 1.2276 0.6012 2.4804 0.0000 1.5336 1.2024 0.6264
#> 
#> 
#> $cost_incur_person
#> $cost_incur_person[[1]]
#>  [1]  -8.480538  -5.769061 -13.602230  -3.941800 -11.886553  -9.594122
#>  [7]  -9.318991 -10.007683  -8.867694  -9.136480
#> 
#> $cost_incur_person[[2]]
#>  [1] -21.045602 -11.809604  -7.781309 -24.437684  -6.627897  -7.602122
#>  [7]  -8.146484  -5.966299   0.000000 -24.162497
#> 
#> $cost_incur_person[[3]]
#>  [1] -14.266718 -28.158966  -9.056345 -21.248178 -10.204341 -47.284430
#>  [7]   0.000000 -13.643529 -21.818627  -8.927014
#> 
#> 
#> $E_cost_incur
#> $E_cost_incur[[1]]
#> [1] -9060.515
#> 
#> $E_cost_incur[[2]]
#> [1] -11757.95
#> 
#> $E_cost_incur[[3]]
#> [1] -17460.81
#> 
#> 
#> $E_cost_incur_person
#> $E_cost_incur_person[[1]]
#> [1] -9.060515
#> 
#> $E_cost_incur_person[[2]]
#> [1] -11.75795
#> 
#> $E_cost_incur_person[[3]]
#> [1] -17.46081
#> 
#> 
#> $QALYgain_person
#> $QALYgain_person[[1]]
#>  [1] 0.0009396 0.0006012 0.0006012 0.0006012 0.0006264 0.0006264 0.0012456
#>  [8] 0.0009072 0.0009396 0.0006012
#> 
#> $QALYgain_person[[2]]
#>  [1] 0.0009396 0.0009144 0.0009144 0.0018540 0.0002880 0.0009396 0.0015408
#>  [8] 0.0003132 0.0000000 0.0012456
#> 
#> $QALYgain_person[[3]]
#>  [1] 0.0009144 0.0012528 0.0009396 0.0012276 0.0006012 0.0024804 0.0000000
#>  [8] 0.0015336 0.0012024 0.0006264
#> 
#> 
#> $E_QALYgain
#> $E_QALYgain[[1]]
#> [1] 0.76896
#> 
#> $E_QALYgain[[2]]
#> [1] 0.89496
#> 
#> $E_QALYgain[[3]]
#> [1] 1.07784
#> 
#> 
#> $E_QALYgain_person
#> $E_QALYgain_person[[1]]
#> [1] 0.00076896
#> 
#> $E_QALYgain_person[[2]]
#> [1] 0.00089496
#> 
#> $E_QALYgain_person[[3]]
#> [1] 0.00107784

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

Nathan Green

2018-04-04

Introduction

Analysis preparation

Calculate cost-effectiveness statistics for each scenario

Contents