A decision tree is a decision model that represents all possible pathways through sequences of events (nodes), which can be under the experimenter’s control (decisions) or not (chances). A decision tree can be represented visually according to a standardised grammar:
Nodes are linked by edges:
rdecision
builds a Decision Tree model by defining these
elements and their relationships. For example, consider the fictitious
and idealized decision problem, introduced in the package README file,
of choosing between providing two forms of lifestyle advice, offered to
people with vascular disease, which reduce the risk of needing an
interventional procedure. The cost to a healthcare provider of the
interventional procedure (e.g., inserting a stent) is 5000 GBP; the cost
of providing the current form of lifestyle advice, an appointment with a
dietician (“diet”), is 50 GBP and the cost of providing an alternative
form, attendance at an exercise programme (“exercise”), is 750 GBP. If
an advice programme is successful, there is no need for an
interventional procedure. These costs can be defined as scalar
variables, as follows:
<- 50.0
cost_diet <- 750.0
cost_exercise <- 5000.0 cost_stent
The model for this fictional scenario can be defined by the following elements:
<- DecisionNode$new("Programme") decision_node
<- ChanceNode$new("Outcome")
chance_node_diet <- ChanceNode$new("Outcome") chance_node_exercise
<- LeafNode$new("No intervention")
leaf_node_diet_no_stent <- LeafNode$new("Intervention")
leaf_node_diet_stent <- LeafNode$new("No intervention")
leaf_node_exercise_no_stent <- LeafNode$new("Intervention") leaf_node_exercise_stent
These nodes can then be wired into a decision tree graph by defining the edges that link pairs of nodes as actions or reactions.
These are the two programmes being tested. The cost of each action, as described in the example, is embedded into the action definition.
<- Action$new(
action_diet cost = cost_diet, label = "Diet"
decision_node, chance_node_diet,
)<- Action$new(
action_exercise cost = cost_exercise, label = "Exercise"
decision_node, chance_node_exercise, )
These are the possible outcomes of each programme (success or failure), with their relevant probabilities.
To continue our fictional example, in a small trial of the “diet” programme, 12 out of 68 patients (17.6%) avoided having an interventional procedure within one year, and in a separate small trial of the “exercise” programme 18 out of 58 patients (31.0%) avoided the interventional procedure within one year (it is assumed that the baseline characteristics in the two trials were comparable).
<- 12L
s_diet <- 56L
f_diet <- 18L
s_exercise <- 40L f_exercise
Epidemiologically, we can interpret the trial results in terms of the incidence proportions of having an adverse event (needing a stent) within one year, for each treatment strategy.
<- f_diet / (s_diet + f_diet)
ip_diet <- f_exercise / (s_exercise + f_exercise)
ip_exercise <- 1.0 / (ip_diet - ip_exercise) nnt
The incidence proportions are 0.82 and 0.69 for diet and exercise, respectively, noting that we define a programme failure as the need to insert a stent. The number needed to treat is the reciprocal of the difference in incidence proportions; 7.47 people must be allocated to the exercise programme rather than the diet programme to save one adverse event.
These trial results can be represented as probabilities of outcome
success (p_diet
, p_exercise
) derived from the
incidence proportions of the trial results. The probabilities of outcome
failure are denoted with “q” prefixes.
<- 1.0 - ip_diet
p_diet <- 1.0 - ip_exercise
p_exercise <- 1.0 - p_diet
q_diet <- 1.0 - p_exercise q_exercise
These probabilities, as well as the cost associated with each outcome, can then be embedded into the reaction definition.
<- Reaction$new(
reaction_diet_success
chance_node_diet, leaf_node_diet_no_stent,p = p_diet, cost = 0.0, label = "Success"
)
<- Reaction$new(
reaction_diet_failure
chance_node_diet, leaf_node_diet_stent,p = q_diet, cost = cost_stent, label = "Failure"
)
<- Reaction$new(
reaction_exercise_success
chance_node_exercise, leaf_node_exercise_no_stent,p = p_exercise, cost = 0.0, label = "Success"
)
<- Reaction$new(
reaction_exercise_failure
chance_node_exercise, leaf_node_exercise_stent,p = q_exercise, cost = cost_stent, label = "Failure"
)
When all the elements are defined and satisfy the restrictions of a
Decision Tree (see the documentation for the DecisionTree
class for details), the whole model can be built:
<- DecisionTree$new(
dt V = list(
decision_node,
chance_node_diet,
chance_node_exercise,
leaf_node_diet_no_stent,
leaf_node_diet_stent,
leaf_node_exercise_no_stent,
leaf_node_exercise_stent
),E = list(
action_diet,
action_exercise,
reaction_diet_success,
reaction_diet_failure,
reaction_exercise_success,
reaction_exercise_failure
) )
rdecision
includes a draw
method to
generate a diagram of a defined Decision Tree.
$draw() dt
As a decision model, a Decision Tree takes into account the costs,
probabilities and utilities encountered as each strategy is traversed
from left to right. In this example, only two strategies (Diet or
Exercise) exist in the model and can be compared using the
evaluate()
method.
<- dt$evaluate() rs
Run | Programme | Probability | Cost | Benefit | Utility | QALY |
---|---|---|---|---|---|---|
1 | Diet | 1 | 4167.65 | 0 | 1 | 1 |
1 | Exercise | 1 | 4198.28 | 0 | 1 | 1 |
From the evaluation of the two strategies, it is apparent that the Diet strategy has a marginally lower net cost by 30.63 GBP.
Because this example is structurally simple, we can verify the results by direct calculation. The net cost per patient of the diet programme is the cost of delivering the advice (50.00 GBP) plus the cost of inserting a stent (5,000.00 GBP) multiplied by the proportion who require a stent, or the failure rate of the programme, 0.82, equal to 4,167.65 GBP. By a similar argument, the net cost per patient of the exercise programme is 4,198.28 GBP, as the model predicts. If 7.47 patients are required to change from the diet programme to the exercise programme to save one stent, the incremental increase in cost of delivering the advice is the number needed to treat multiplied by the difference in the cost of the programmes, 750.00 GBP - 50.00 GBP, or 5,228.79 GBP. Because this is greater than the cost saved by avoiding one stent (5,000.00 GBP), we can see that the additional net cost per patient of delivering the programme is the difference between these costs, divided by the number needed to treat, or 30.63 GBP, as the model predicts.
Note that this approach aggregates multiple paths that belong to the
same strategy (for example, the Success and Failure paths of the Diet
strategy). The option by = "path"
can be used to evaluate
each path separately.
<- dt$evaluate(by = "path") rp
Run | Programme | Leaf | Probability | Cost | Benefit | Utility | QALY |
---|---|---|---|---|---|---|---|
1 | Diet | Intervention | 0.824 | 4158.82 | 0 | 0.824 | 0.824 |
1 | Diet | No.intervention | 0.176 | 8.82 | 0 | 0.176 | 0.176 |
1 | Exercise | Intervention | 0.690 | 3965.52 | 0 | 0.690 | 0.690 |
1 | Exercise | No.intervention | 0.310 | 232.76 | 0 | 0.310 | 0.310 |
Cost is not the only consideration that can be modelled using a Decision Tree. Suppose that requiring an intervention reduces the quality of life of patients, associated with attending pre-operative appointments, pain and discomfort of the procedure, and adverse events. This is estimated to be associated with a mean disutility of 0.05 for those who receive a stent, assumed to persist over 1 year.
To incorporate this into the model, we can set the utility of the two leaf nodes which are associated with having a stent. Because we are changing the property of two nodes in the tree, and not changing the tree structure, we do not have to rebuild the tree.
<- 0.05
du_stent $set_utility(1.0 - du_stent)
leaf_node_diet_stent$set_utility(1.0 - du_stent)
leaf_node_exercise_stent<- dt$evaluate() rs
Run | Programme | Probability | Cost | Benefit | Utility | QALY |
---|---|---|---|---|---|---|
1 | Diet | 1 | 4167.65 | 0 | 0.959 | 0.959 |
1 | Exercise | 1 | 4198.28 | 0 | 0.966 | 0.966 |
In this case, while the Diet strategy is preferred from a cost
perspective, the utility of the Exercise strategy is superior.
rdecision
also calculates Quality-adjusted life-years
(QALYs) taking into account the time horizon of the model (in this case,
the default of one year was used, and therefore QALYs correspond to the
Utility values). From these figures, the Incremental cost-effectiveness
ratio (ICER) can be easily calculated:
<- rs[[which(rs[, "Programme"] == "Exercise"), "Cost"]] -
delta_c which(rs[, "Programme"] == "Diet"), "Cost"]]
rs[[<- rs[[which(rs[, "Programme"] == "Exercise"), "Utility"]] -
delta_u which(rs[, "Programme"] == "Diet"), "Utility"]]
rs[[<- delta_c / delta_u icer
resulting in a cost of 4575.76 GBP per QALY gained in choosing the more effective Exercise strategy over the cheaper Diet strategy.
This can be verified by direct calculation, by dividing the difference in net costs of the two programmes (30.63 GBP) by the increase in QALYs due to stents saved (0.05 multiplied by the difference in success rates of the programme, (0.31 - 0.176), or 0.007 QALYs), an ICER of 4,575.76 GBP / QALY.
The model shown above uses a fixed value for each parameter, resulting in a single point estimate for each model result. However, parameters may be affected by uncertainty: for example, the success probability of each strategy is extracted from a small trial of few patients. This uncertainty can be incorporated into the Decision Tree model by representing individual parameters with a statistical distribution, then repeating the evaluation of the model multiple times with each run randomly drawing parameters from these defined distributions.
In rdecision
, model variables that are described by a
distribution are represented by ModVar
objects. Many
commonly used distributions, such as the Normal, Log-Normal, Gamma and
Beta distributions are included in the package, and additional
distributions can be easily implemented from the generic
ModVar
class. Additionally, model variables that are
calculated from other r probabilistic variables using an expression can
be represented as ExprModVar
objects.
Fixed costs can be left as numerical values, or also be represented
by ModVar
s which ensures that they are included in variable
tabulations.
<- ConstModVar$new("Cost of diet programme", "GBP", 50.0)
cost_diet <- ConstModVar$new("Cost of exercise programme", "GBP", 750.0)
cost_exercise <- ConstModVar$new("Cost of stent intervention", "GBP", 5000.0) cost_stent
In our simplified example, the probability of success of each strategy should include the uncertainty associated with the small sample that they are based on. This can be represented statistically by a Beta distribution, a probability distribution constrained to the interval [0, 1]. A Beta distribution that captures the results of the trials can be defined by the alpha (observed successes) and beta (observed failures) parameters.
<- BetaModVar$new(
p_diet alpha = s_diet, beta = f_diet, description = "P(diet)", units = ""
)<- BetaModVar$new(
p_exercise alpha = s_exercise, beta = f_exercise, description = "P(exercise)", units = ""
)
<- ExprModVar$new(
q_diet ::quo(1.0 - p_diet), description = "1 - P(diet)", units = ""
rlang
)<- ExprModVar$new(
q_exercise ::quo(1.0 - p_exercise), description = "1 - P(exercise)", units = ""
rlang )
These distributions describe the probability of success of each
strategy; by the constraints of a Decision Tree, the sum of all
probabilities associated with a chance node must be 1, so the
probability of failure should be calculated as 1 - p(Success). This can
be represented by an ExprModVar
.
The newly defined ModVars
can be incorporated into the
Decision Tree model using the same grammar as the non-probabilistic
model. Because the actions and reactions are objects already included in
the tree, we can change their properties using set_
calls
and those new properties will be used when the tree is evaluated.
$set_cost(cost_diet)
action_diet$set_cost(cost_exercise) action_exercise
$set_probability(p_diet)
reaction_diet_success
$set_probability(q_diet)
reaction_diet_failure$set_cost(cost_stent)
reaction_diet_failure
$set_probability(p_exercise)
reaction_exercise_success
$set_probability(q_exercise)
reaction_exercise_failure$set_cost(cost_stent) reaction_exercise_failure
All the probabilistic variables included in the model can be
tabulated using the modvar_table()
method, which details
the distribution definition and some useful parameters, such as mean, SD
and 95% CI.
Description | Units | Distribution | Mean | E | SD | Q2.5 | Q97.5 | Est |
---|---|---|---|---|---|---|---|---|
Cost of stent intervention | GBP | Const(5000) | 5000.000 | 5000.000 | 0.000 | 5000.000 | 5000.000 | FALSE |
1 - P(diet) | 1 - p_diet | 0.824 | 0.824 | 0.046 | 0.722 | 0.902 | TRUE | |
P(diet) | Be(12,56) | 0.176 | 0.176 | 0.046 | 0.096 | 0.275 | FALSE | |
Cost of exercise programme | GBP | Const(750) | 750.000 | 750.000 | 0.000 | 750.000 | 750.000 | FALSE |
Cost of diet programme | GBP | Const(50) | 50.000 | 50.000 | 0.000 | 50.000 | 50.000 | FALSE |
P(exercise) | Be(18,40) | 0.310 | 0.310 | 0.060 | 0.199 | 0.434 | FALSE | |
1 - P(exercise) | 1 - p_exercise | 0.690 | 0.691 | 0.060 | 0.565 | 0.798 | TRUE |
A call to the evaluate()
method with the default
settings uses the expected (mean) value of each variable, and so
replicates the point estimate above.
<- dt$evaluate() rs
Run | Programme | Probability | Cost | Benefit | Utility | QALY |
---|---|---|---|---|---|---|
1 | Diet | 1 | 4167.65 | 0 | 0.96 | 0.96 |
1 | Exercise | 1 | 4198.28 | 0 | 0.97 | 0.97 |
However, because each variable is described by a distribution, it is now possible to explore the range of possible values consistent with the model. For example, a lower and upper bound can be estimated by setting each variable to its 2.5-th or 97.5-th percentile:
<- dt$evaluate(setvars = "q2.5")
rs_025 <- dt$evaluate(setvars = "q97.5") rs_975
The costs for each choice when all variables are at their upper and lower confidence levels are as follows:
Q2.5 | Q97.5 | |
---|---|---|
Diet | 4569.40 | 3676.00 |
Exercise | 4754.75 | 3579.87 |
To sample the possible outcomes in a completely probabilistic way,
the setvar = "random"
option can be used, which draws a
random value from the distribution of each variable. Repeating this
process a sufficiently large number of times builds a collection of
results compatible with the model definition, which can then be used to
calculate ranges and confidence intervals of the estimated values.
<- 1000L
N <- dt$evaluate(setvars = "random", by = "run", N = N) rs
The estimates of cost for each intervention can be plotted as follows:
plot(
"Cost.Diet"],
rs[, "Cost.Exercise"],
rs[, pch = 20L,
xlab = "Cost of diet (GBP)", ylab = "Cost of exercise (GBP)",
main = paste(N, "simulations of vascular disease prevention model")
)abline(a = 0.0, b = 1.0, col = "red")
A tabular summary is as follows:
Cost.Diet | Cost.Exercise | |
---|---|---|
Min. :3351 | Min. :3181 | |
1st Qu.:4032 | 1st Qu.:4008 | |
Median :4188 | Median :4229 | |
Mean :4169 | Mean :4207 | |
3rd Qu.:4329 | 3rd Qu.:4422 | |
Max. :4742 | Max. :5018 |
The variables can be further manipulated, for example calculating the difference in cost between the two strategies for each run of the randomised model:
"Difference"] <- rs[, "Cost.Diet"] - rs[, "Cost.Exercise"]
rs[, <- quantile(rs[, "Difference"], c(0.025, 0.975)) CI
hist(
"Difference"], 100L, main = "Distribution of saving",
rs[, xlab = "Saving (GBP)"
)
::kable(
knitr: 10L, c(1L, 3L, 8L, 12L)], digits = 2L,
rs[1L row.names = FALSE
)
Run | Cost.Diet | Cost.Exercise | Difference |
---|---|---|---|
1 | 4294.88 | 4086.98 | 207.90 |
2 | 4251.22 | 4032.02 | 219.20 |
3 | 3882.10 | 3880.43 | 1.67 |
4 | 3697.33 | 4243.93 | -546.59 |
5 | 4221.99 | 4422.15 | -200.16 |
6 | 4201.93 | 4271.73 | -69.80 |
7 | 4212.78 | 4236.26 | -23.49 |
8 | 3788.91 | 4184.11 | -395.19 |
9 | 4087.15 | 3829.22 | 257.92 |
10 | 3698.31 | 4390.85 | -692.54 |
Plotting the distribution of the difference of the two costs reveals that, in this model, the uncertainties in the input parameters are large enough that either strategy could have a lower net cost, within a 95% confidence interval [-771.19, 749.44].
rdecision
provides a threshold
method to
compare two strategies and identify, for a given variable, the value at
which one strategy becomes cost saving over the other:
<- dt$threshold(
cost_threshold index = list(action_exercise),
ref = list(action_diet),
outcome = "saving",
mvd = cost_exercise$description(),
a = 0.0, b = 5000.0, tol = 0.1
)
<- dt$threshold(
success_threshold index = list(action_exercise),
ref = list(action_diet),
outcome = "saving",
mvd = p_exercise$description(),
a = 0.0, b = 1.0, tol = 0.001
)
By univariate threshold analysis, the exercise program will be cost saving when its cost of delivery is less than 719.38 GBP or when its success rate is greater than 31.7%.
These can be verified by direct calculation. The cost of delivering the exercise programme at which its net cost equals the net cost of the diet programme is when the difference between the two programme delivery costs multiplied by the number needed to treat becomes equal to the cost saved by avoiding one stent. This is 719.37, in agreement with the model. The threshold success rate for the exercise programme is when the number needed to treat is reduced such that the net cost of the two programmes is equal, i.e., to 7.14, from which we can calculate the programme success rate threshold as 0.32, in agreement with the model.