D1_decision_tree_DARTH_HVE_example.RmdThis is an example that we use to illustrate decision trees in our DARTH teaching materials. The example involves the treatment of herpes simplex ecncephalopathy (HVE) by either treating everyone without biopsy (Treat), not treating or taking biopsy (DoNotTreat), and Taking biopsy and treat if biopsy posivite (Biopsy).
We will first define the gmod model and the parameters and then we will show how one would code the model direclty in R and compare the outcomes.
rm(list = ls())
library(gmod)
#> Please note that active development of gmod is now moved to twig. Please install twig https://hjalal.github.io/twig/ and review the vignettes for the latest features.gmod model definition p_comp() takes in two arguments, one is the special keyword decision, and the other is a custom one. This is a good example that custom arguments and values can be passed to define functions for the event_mapping and the payoffs. Thus, the return values can be a function of these arguments and it can reduce the amount of typing required if one choses to define multiple functions.
mygmod <- gmod() +
decisions("DoNotTreat", "Treat", "Biopsy") +
#final_outcomes("Death", "HVE_comp", "no_HVE_comp", "OVE_comp", "no_OVE_comp") +
#events("DIE", "HVE","get_comp") +
event(name = "DIE",
scenarios = c(T, F),
probs = c(pDie(decision), Inf),
outcomes = c("Death", "HVE_event")) +
event(name = "HVE_event",
scenarios = c(T, F),
probs = c(p_HVE, Inf),
outcomes = c("get_HVE_comp", "get_OVE_comp")) +
event(name = "get_HVE_comp",
scenarios = c(T, F),
probs = c(p_comp(decision, HVE = TRUE), Inf),
outcomes = c("HVE_comp", "no_HVE_comp")) +
event(name = "get_OVE_comp",
scenarios = c(T, F),
probs = c(p_comp(decision, HVE = FALSE), Inf),
outcomes = c("OVE_comp", "no_OVE_comp")) +
#payoffs(cost = cost(decision, final_outcome, prop_with_event("HVE"=TRUE, decision)),
payoffs(cost = cost(decision, final_outcome),
effectiveness = effectiveness(decision, final_outcome))Define the model parameters.
v_names_str <- c("No Tx", "Tx All", "Biopsy") # names of strategies
n_str <- length(v_names_str) # number of strategies
wtp <- 100000 # willingness to pay threshold
params <- list(
wtp = 100000 , # willingness to pay threshold
# Probabilities,
p_HVE = 0.52 ,# prevalence of HVE
p_HVE_comp = 0.71 ,# complications with untreated HVE
p_OVE_comp = 0.01 ,# complications with untreated OVE
p_HVE_comp_tx = 0.36 ,# complications with treated HVE
p_OVE_comp_tx = 0.20 ,# complications with treated OVE
p_biopsy_death = 0.005 ,# probability of death due to biopsy
# Costs,
c_VE = 1200 ,# cost of viral encephalitis care without complications
c_VE_comp = 9000 ,# cost of viral encephalitis care with complications
c_tx = 9500 ,# cost of treatment
c_biopsy = 25000 ,# cost of brain biopsy
# QALYs,
q_VE = 20 ,# remaining QALYs for those without VE-related complications
q_VE_comp = 19 ,# remaining QALYs for those with VE-related complications
q_loss_biopsy = 0.01 ,# one-time QALY loss due to brain biopsy
q_death_biopsy = 0 # remaining QALYs for those who died during biopsy
)define custom user functions
pDie <- function(decision){
if (decision == "Biopsy") p_biopsy_death else 0
}
p_comp <- function(decision, HVE){
if (decision == "DoNotTreat" & HVE ) return(p_HVE_comp)
if (decision == "DoNotTreat" & !HVE )return(p_OVE_comp)
if (decision == "Treat" & HVE ) return(p_HVE_comp_tx)
if (decision == "Treat" & !HVE ) return(p_OVE_comp_tx)
if (decision == "Biopsy" & HVE) return(p_HVE_comp_tx)
if (decision == "Biopsy" & !HVE) return(p_OVE_comp)
}
c_HVE <- function(decision){
if (decision == "biopsy") c_tx else 0
}
cost <- function(decision, final_outcome){
c_biopsy*(decision=="Biopsy") +
c_tx*(decision=="Treat" | (decision=="Biopsy" & final_outcome %in% c("HVE_comp", "no_HVE_comp"))) +
c_VE_comp*(final_outcome %in% c("HVE_comp", "OVE_comp")) +
c_VE*(final_outcome %in% c("no_HVE_comp", "no_OVE_comp"))
}
effectiveness <- function(decision, final_outcome){
-q_loss_biopsy*(decision=="Biopsy") +
q_VE_comp*(final_outcome %in% c("HVE_comp", "OVE_comp")) +
q_VE*(final_outcome %in% c("no_HVE_comp", "no_OVE_comp"))
}testing some of the functions outside of the model
#p_comp(decision = "Biopsy", HVE = FALSE)
#
#cost("Biopsy", "HVE_comp")
#effectiveness("Biopsy", "Death")
model_struc <- gmod_build(mygmod)
model_struc
#> $decisions
#> [1] "DoNotTreat" "Treat" "Biopsy"
#>
#> $n_decisions
#> [1] 3
#>
#> $n_events
#> [1] 4
#>
#> $final_outcomes
#> [1] "Death" "HVE_comp" "no_HVE_comp" "OVE_comp" "no_OVE_comp"
#>
#> $n_final_outcomes
#> [1] 5
#>
#> $events
#> [1] "DIE" "HVE_event" "get_HVE_comp" "get_OVE_comp"
#>
#> $payoffs
#> $payoffs$cost
#> cost(decision, final_outcome)
#>
#> $payoffs$effectiveness
#> effectiveness(decision, final_outcome)
#>
#>
#> $payoff_names
#> [1] "cost" "effectiveness"
#>
#> $n_payoffs
#> [1] 2
#>
#> $final_outcome_formulae
#> # A tibble: 15 × 10
#> # Groups: decision, final_outcome [15]
#> decision final_outcome path_id probs DIE HVE_event get_HVE_comp
#> <chr> <chr> <dbl> <chr> <chr> <chr> <chr>
#> 1 Biopsy Death 1 (pDie('Biopsy'… TRUE FALSE FALSE
#> 2 Biopsy HVE_comp 2 (1-(pDie('Biop… FALSE TRUE TRUE
#> 3 Biopsy OVE_comp 4 (1-(pDie('Biop… FALSE FALSE FALSE
#> 4 Biopsy no_HVE_comp 3 (1-(pDie('Biop… FALSE TRUE FALSE
#> 5 Biopsy no_OVE_comp 5 (1-(pDie('Biop… FALSE FALSE FALSE
#> 6 DoNotTreat Death 1 (pDie('DoNotTr… TRUE FALSE FALSE
#> 7 DoNotTreat HVE_comp 2 (1-(pDie('DoNo… FALSE TRUE TRUE
#> 8 DoNotTreat OVE_comp 4 (1-(pDie('DoNo… FALSE FALSE FALSE
#> 9 DoNotTreat no_HVE_comp 3 (1-(pDie('DoNo… FALSE TRUE FALSE
#> 10 DoNotTreat no_OVE_comp 5 (1-(pDie('DoNo… FALSE FALSE FALSE
#> 11 Treat Death 1 (pDie('Treat')) TRUE FALSE FALSE
#> 12 Treat HVE_comp 2 (1-(pDie('Trea… FALSE TRUE TRUE
#> 13 Treat OVE_comp 4 (1-(pDie('Trea… FALSE FALSE FALSE
#> 14 Treat no_HVE_comp 3 (1-(pDie('Trea… FALSE TRUE FALSE
#> 15 Treat no_OVE_comp 5 (1-(pDie('Trea… FALSE FALSE FALSE
#> # ℹ 3 more variables: get_OVE_comp <chr>, cost <chr>, effectiveness <chr>
#>
#> $summary_formulae
#> # A tibble: 3 × 3
#> decision cost effectiveness
#> <chr> <chr> <chr>
#> 1 Biopsy (pDie('Biopsy'))*cost('Biopsy', 'Death')+(1-(pDie('B… (pDie('Biops…
#> 2 DoNotTreat (pDie('DoNotTreat'))*cost('DoNotTreat', 'Death')+(1-… (pDie('DoNot…
#> 3 Treat (pDie('Treat'))*cost('Treat', 'Death')+(1-(pDie('Tre… (pDie('Treat…
#>
#> attr(,"class")
#> [1] "gmod_decision"
model_res <- gmod_gen_model_function(model_struc)
#> Warning in gmod_gen_model_function.gmod_decision(model_struc): Model function my_decision_model is generated. It can be run by calling it directly:
#> my_decision_model(params)
model_res
#> $decisions
#> [1] "DoNotTreat" "Treat" "Biopsy"
#>
#> $n_decisions
#> [1] 3
#>
#> $n_events
#> [1] 4
#>
#> $final_outcomes
#> [1] "Death" "HVE_comp" "no_HVE_comp" "OVE_comp" "no_OVE_comp"
#>
#> $n_final_outcomes
#> [1] 5
#>
#> $events
#> [1] "DIE" "HVE_event" "get_HVE_comp" "get_OVE_comp"
#>
#> $payoffs
#> $payoffs$cost
#> cost(decision, final_outcome)
#>
#> $payoffs$effectiveness
#> effectiveness(decision, final_outcome)
#>
#>
#> $payoff_names
#> [1] "cost" "effectiveness"
#>
#> $n_payoffs
#> [1] 2
#>
#> $final_outcome_formulae
#> # A tibble: 15 × 10
#> # Groups: decision, final_outcome [15]
#> decision final_outcome path_id probs DIE HVE_event get_HVE_comp
#> <chr> <chr> <dbl> <chr> <chr> <chr> <chr>
#> 1 Biopsy Death 1 (pDie('Biopsy'… TRUE FALSE FALSE
#> 2 Biopsy HVE_comp 2 (1-(pDie('Biop… FALSE TRUE TRUE
#> 3 Biopsy OVE_comp 4 (1-(pDie('Biop… FALSE FALSE FALSE
#> 4 Biopsy no_HVE_comp 3 (1-(pDie('Biop… FALSE TRUE FALSE
#> 5 Biopsy no_OVE_comp 5 (1-(pDie('Biop… FALSE FALSE FALSE
#> 6 DoNotTreat Death 1 (pDie('DoNotTr… TRUE FALSE FALSE
#> 7 DoNotTreat HVE_comp 2 (1-(pDie('DoNo… FALSE TRUE TRUE
#> 8 DoNotTreat OVE_comp 4 (1-(pDie('DoNo… FALSE FALSE FALSE
#> 9 DoNotTreat no_HVE_comp 3 (1-(pDie('DoNo… FALSE TRUE FALSE
#> 10 DoNotTreat no_OVE_comp 5 (1-(pDie('DoNo… FALSE FALSE FALSE
#> 11 Treat Death 1 (pDie('Treat')) TRUE FALSE FALSE
#> 12 Treat HVE_comp 2 (1-(pDie('Trea… FALSE TRUE TRUE
#> 13 Treat OVE_comp 4 (1-(pDie('Trea… FALSE FALSE FALSE
#> 14 Treat no_HVE_comp 3 (1-(pDie('Trea… FALSE TRUE FALSE
#> 15 Treat no_OVE_comp 5 (1-(pDie('Trea… FALSE FALSE FALSE
#> # ℹ 3 more variables: get_OVE_comp <chr>, cost <chr>, effectiveness <chr>
#>
#> $summary_formulae
#> # A tibble: 3 × 3
#> decision cost effectiveness
#> <chr> <chr> <chr>
#> 1 Biopsy (pDie('Biopsy'))*cost('Biopsy', 'Death')+(1-(pDie('B… (pDie('Biops…
#> 2 DoNotTreat (pDie('DoNotTreat'))*cost('DoNotTreat', 'Death')+(1-… (pDie('DoNot…
#> 3 Treat (pDie('Treat'))*cost('Treat', 'Death')+(1-(pDie('Tre… (pDie('Treat…
#>
#> attr(,"class")
#> [1] "gmod_decision"
my_decision_model(params)
#> cost effectiveness
#> Biopsy 32599.41 19.69896
#> DoNotTreat 4117.20 19.62600
#> Treat 12908.96 19.71680And here is what the code looks like if one wants to code the transitions all in R directly without using gmod
# Create vector of weights for each strategy
v_w_no_tx <- c( p_HVE * p_HVE_comp , # HVE, complications
p_HVE * (1 - p_HVE_comp) , # HVE, no complications
(1 - p_HVE) * p_OVE_comp , # OVE, complications
(1 - p_HVE) * (1 - p_OVE_comp)) # OVE, no complications
v_w_tx <- c( p_HVE * p_HVE_comp_tx , # HVE w/tx, complications
p_HVE * (1 - p_HVE_comp_tx) , # HVE w/tx, no complications
(1 - p_HVE) * p_OVE_comp_tx , # OVE w/tx, complications
(1 - p_HVE) * (1 - p_OVE_comp_tx)) # OVE w/tx, no complications
v_w_biopsy <- c(p_biopsy_death , # biopsy death
# no biopsy death., HVE w/tx, complications
(1-p_biopsy_death) * p_HVE * p_HVE_comp_tx ,
# no biopsy death., HVE w/tx, no complications
(1-p_biopsy_death) * p_HVE * (1-p_HVE_comp_tx) ,
# no biopsy death., OVE, complications
(1-p_biopsy_death) * (1-p_HVE) * p_OVE_comp ,
# no biopsy death., OVE, no complications
(1-p_biopsy_death) * (1-p_HVE) * (1 - p_OVE_comp))
# Create vector of final_outcomes (QALYs) for each strategy
v_qaly_no_tx <- c(q_VE_comp , # HVE, complications
q_VE , # HVE, no complications
q_VE_comp , # OVE, complications
q_VE) # OVE, no complications
v_qaly_tx <- c(q_VE_comp , # HVE, complications
q_VE , # HVE, no complications
q_VE_comp , # OVE, complications
q_VE) # OVE, no complications
v_qaly_biopsy <- -q_loss_biopsy + # loss due to biopsy
c(q_death_biopsy , # biopsy complications
q_VE_comp , # no biopsy comp., HVE w/tx, complications
q_VE , # no biopsy comp., HVE w/tx, no complications
q_VE_comp , # no biopsy comp., OVE, complications
q_VE) # no biopsy comp., OVE, no complications
# Create vector of costs for each strategy
v_cost_no_tx <- c(c_VE_comp , # HVE, complications
c_VE , # HVE, no complications
c_VE_comp , # OVE, complications
c_VE) # OVE, no complications
v_cost_tx <- c_tx + # cost of treatment
c(c_VE_comp , # HVE, complications
c_VE , # HVE, no complications
c_VE_comp , # OVE, complications
c_VE) # OVE, no complications
v_cost_biopsy <- c_biopsy + # cost of biopsy procedure
c(0 , # cost of death (zero)
c_VE_comp + c_tx , # no biopsy comp., HVE w/tx, complications
c_VE + c_tx , # no biopsy comp., HVE w/tx, no complications
c_VE_comp , # no biopsy comp., OVE, complications
c_VE) # no biopsy comp., OVE, no complications
# Calculate total utilities for each strategy
total_qaly_no_tx <- v_w_no_tx %*% v_qaly_no_tx
total_qaly_tx <- v_w_tx %*% v_qaly_tx
total_qaly_biopsy <- v_w_biopsy %*% v_qaly_biopsy
# Calculate total costs for each strategy
total_cost_no_tx <- v_w_no_tx %*% v_cost_no_tx
total_cost_tx <- v_w_tx %*% v_cost_tx
total_cost_biopsy <- v_w_biopsy %*% v_cost_biopsy
# vector of total QALYs
v_total_qaly <- c(total_qaly_no_tx, total_qaly_tx, total_qaly_biopsy)
# vector of total costs
v_total_cost <- c(total_cost_no_tx, total_cost_tx, total_cost_biopsy)
# calculate vector of nmb
v_nmb <- v_total_qaly * wtp - v_total_cost
# Name final_outcomes
names(v_total_qaly) <- v_names_str # names for the elements of the total QALYs vector
names(v_total_cost) <- v_names_str # names for the elements of the total cost vector
names(v_nmb) <- v_names_str # names for the elements of the nmb vector
df_output <- data.frame(Strategy = v_names_str,
Cost = v_total_cost,
Effect = v_total_qaly,
NMB = v_nmb)
# model output
df_output
#> Strategy Cost Effect NMB
#> No Tx No Tx 4117.20 19.62600 1958483
#> Tx All Tx All 12908.96 19.71680 1958771
#> Biopsy Biopsy 32599.41 19.69896 1937297