Introduction

This vignette is an example of modelling a decision tree using the rdecision package, with probabilistic sensitivity analysis (PSA). It is based on the model reported by Jenks et al¹ in which a transparent dressing used to secure vascular catheters (Tegaderm CHG) was compared with a standard dressing.

Model variables

Source variables

Eleven source variables were used in the model. The choice of variables, their distributions and their parameters are taken from Table 4 of Jenks et al¹, with the following additional information:

The baseline catheter-related blood stream infection (CRBSI) rate was modelled as a Gamma distribution fitted by the method of moments to a mean of 1.48 (per 1000 catheter days) and a standard deviation of 0.12 (per 1000 catheter days).
The baseline local site infection (LSI) rate was modelled as a Gamma distribution fitted by the method of moments to a mean of 0.14 (per 1000 catheter days) and an assumed standard deviation of 0.5 (/1000 catheter days).
The baseline rate of dermatitis was modelled as a Beta distribution, based on one observed case in a trial of 476 catheter uses.
The effect size of Tegaderm, expressed as the hazard ratio of Tegaderm compared with standard dressings for CRBSI and LSI, and the relative risk of Tegaderm compared with standard dressings for dermatitis, was modelled in each case using a log normal distribution. This was fitted to a sample mean and sample standard deviation on the natural scale, by using the “LN7” parametrization of LogNormModVar.
The probabilities of CRBSI and LSI for standard dressings (p) were modified by the hazard ratio r for Tegaderm using the form p * r. This is an approximation which holds only for very small rates.
Relative risks were also applied as multipliers. This is an approximation which holds only for very small rates.

The model variables were constructed as follows:

# baseline risk
r.CRBSI <- GammaModVar$new(
  "Baseline CRBSI rate",  "/1000 catheter days",
  shape = (1.48 ^ 2L) / (0.12 ^ 2L),
  scale = (0.12 ^ 2L) / 1.48
)
r.LSI <- GammaModVar$new(
  "Baseline LSI rate", "/1000 catheter days",
  shape = (0.14 ^ 2L) / (0.5 ^ 2L),
  scale = (0.5 ^ 2L) / 0.14
)
r.Dermatitis <- BetaModVar$new(
  "Baseline dermatitis risk", "/catheter", alpha = 1L, beta = 475L
)
# relative effectiveness
hr.CRBSI <- LogNormModVar$new(
  "Tegaderm CRBSI HR", "HR",
  p1 = 0.402, p2 = (0.868 - 0.186) / (2L * 1.96), param = "LN7"
)
hr.LSI <- LogNormModVar$new(
  "Tegaderm LSI HR", "HR",
  p1 = 0.402, p2 = (0.868 - 0.186) / (2L * 1.96), param = "LN7"
)
rr.Dermatitis <- LogNormModVar$new(
  "Tegaderm Dermatitis RR", "RR", p1 = 1.0, p2 = 0.5, param = "LN7"
)
# cost variables
c.CRBSI <- GammaModVar$new(
  "CRBSI cost", "GBP",
  shape = (9900.0 ^ 2L) / (3000.0 ^ 2L),
  scale = (3000.0 ^ 2L) / 9900.0
)
c.LSI <- GammaModVar$new(
  "LSI cost", "GBP",
  shape = (100.0 ^ 2L) / (30.0 ^ 2L),
  scale = (30.0 ^ 2L) / 100.0
)
c.Dermatitis <- GammaModVar$new(
  "Dermatitis cost", "GBP",
  shape = (6.0 ^ 2L) / (3.0 ^ 2L),
  scale = (3.0 ^ 2L) / 6.0
)
# number of dressings and days with catheter
n.dressings <- GammaModVar$new(
  "No. dressings", "dressings",
  shape = (3.0 ^ 2L) / (2.0 ^ 2L),
  scale = (2.0 ^ 2L) / 3.0
)
n.cathdays <- GammaModVar$new(
  "No. days with catheter", "days",
  shape = (10.0 ^ 2L) / (5.0 ^ 2L),
  scale = (5.0 ^ 2L) / 10.0
)

Model variable expressions

Variables in the model may be included in the decision tree via mathematical expressions, which involve model variables and are themselves model variables. Forms of expression involving R functions and multiple model variables are supported, provided they conform to R syntax. The following code creates the model variable expressions to be used as values in the decision tree edges.

p.CRBSI.S <- ExprModVar$new(
  "P(CRBSI | standard dressing)", "P",
  rlang::quo(r.CRBSI * n.cathdays / 1000.0)
)
p.CRBSI.T <- ExprModVar$new(
  "P(CRBSI|Tegaderm)", "P",
  rlang::quo(p.CRBSI.S * hr.CRBSI)
)
p.LSI.S <- ExprModVar$new(
  "P(LSI | Standard)", "/patient",
  rlang::quo(r.LSI * n.cathdays / 1000.0)
)
p.LSI.T <- ExprModVar$new(
  "P(LSI | Tegaderm)", "P", rlang::quo(p.LSI.S * hr.LSI)
)
p.Dermatitis.S <- ExprModVar$new(
  "P(dermatitis | standard dressing)", "P",
  rlang::quo(r.Dermatitis)
)
p.Dermatitis.T <- ExprModVar$new(
  "P(dermatitis | Tegaderm)", "P",
  rlang::quo(p.Dermatitis.S * rr.Dermatitis)
)
c.Tegaderm <- ExprModVar$new(
  "Tegaderm CHG cost", "GBP", rlang::quo(6.26 * n.dressings)
)
c.Standard <- ExprModVar$new(
  "Standard dressing cost", "GBP", rlang::quo(1.54 * n.dressings)
)

The decision tree

Constructing the tree

The following code constructs the decision tree based on Figure 2 of Jenks et al¹. In the formulation used by rdecision, the decision tree is constructed from sets of decision, chance and leaf nodes and from edges (actions and reactions). Leaf nodes are synonymous with pathways in Briggs’ terminology². The time horizon is not stated explicitly in the model, and is assumed to be 7 days. It was implied that the time horizon was ICU stay plus some follow-up, and the costs reflect those incurred in that period, so the assumption of 7 days does not affect the rdecision implementation of the model.

The tree is somewhat more complex than Figure 2 of Jenks et al because it allows for patients to have more than one adverse event (AE) during their stay (whereas their Figure 2 implies that only one event per patient is possible). The rates of AE were estimated independently, and allow for multiple events, (figure 1).

In rdecision, if the probability associated with one of the reactions from any chance node is set to missing (NA), it will be computed before each evaluation of the tree to ensure that the probabilities sum to unity.

# create decision tree
th <- as.difftime(7L, units = "days")
# standard dressing
t01 <- LeafNode$new("t01", interval = th)
t02 <- LeafNode$new("t02", interval = th)
c01 <- ChanceNode$new()
e01 <- Reaction$new(
  c01, t01, p = p.Dermatitis.S, cost = c.Dermatitis, label = "Dermatitis"
)
e02 <- Reaction$new(
  c01, t02, p = NA_real_, cost = 0.0, label = "No dermatitis"
)
t03 <- LeafNode$new("t03", interval = th)
t04 <- LeafNode$new("t04", interval = th)
c02 <- ChanceNode$new()
e03 <- Reaction$new(
  c02, t03, p = p.Dermatitis.S, cost = c.Dermatitis, label = "Dermatitis"
)
e04 <- Reaction$new(
  c02, t04, p = NA_real_, cost = 0.0, label = "No dermatitis"
)
c03 <- ChanceNode$new()
e05 <- Reaction$new(c03, c01, p = p.LSI.S, cost = c.LSI, label = "LSI")
e06 <- Reaction$new(c03, c02, p = NA_real_, cost = 0.0, label = "No LSI")
t11 <- LeafNode$new("t11", interval = th)
t12 <- LeafNode$new("t12", interval = th)
c11 <- ChanceNode$new()
e11 <- Reaction$new(
  c11, t11, p = p.Dermatitis.S, cost = c.Dermatitis, label = "Dermatitis"
)
e12 <- Reaction$new(
  c11, t12, p = NA_real_, cost = 0.0, label = "No Dermatitis"
)
t13 <- LeafNode$new("t13", interval = th)
t14 <- LeafNode$new("t14", interval = th)
c12 <- ChanceNode$new()
e13 <- Reaction$new(
  c12, t13, p = p.Dermatitis.S, cost = c.Dermatitis, label = "Dermatitis"
)
e14 <- Reaction$new(
  c12, t14, p = NA_real_, cost = 0.0, label = "No dermatitis"
)
c13 <- ChanceNode$new()
e15 <- Reaction$new(c13, c11, p = p.LSI.S, cost = c.LSI, label = "LSI")
e16 <- Reaction$new(c13, c12, p = NA_real_, cost = 0.0, label = "No LSI")
c23 <- ChanceNode$new()
e21 <- Reaction$new(c23, c03, p = p.CRBSI.S, cost = c.CRBSI, label = "CRBSI")
e22 <- Reaction$new(c23, c13, p = NA_real_, cost = 0.0, label = "No CRBSI")

# Tegaderm branch
t31 <- LeafNode$new("t31", interval = th)
t32 <- LeafNode$new("t32", interval = th)
c31 <- ChanceNode$new()
e31 <- Reaction$new(
  c31, t31, p = p.Dermatitis.T, cost = c.Dermatitis, label = "Dermatitis"
)
e32 <- Reaction$new(
  c31, t32, p = NA_real_, cost = 0.0, label = "no dermatitis"
)
t33 <- LeafNode$new("t33", interval = th)
t34 <- LeafNode$new("t34", interval = th)
c32 <- ChanceNode$new()
e33 <- Reaction$new(
  c32, t33, p = p.Dermatitis.T, cost = c.Dermatitis, label = "Dermatitis"
)
e34 <- Reaction$new(
  c32, t34, p = NA_real_, cost = 0.0, label = "No dermatitis"
)
c33 <- ChanceNode$new()
e35 <- Reaction$new(c33, c31, p = p.LSI.T, cost = c.LSI, label = "LSI")
e36 <- Reaction$new(c33, c32, p = NA_real_, cost = 0.0, label = "No LSI")
t41 <- LeafNode$new("t41", interval = th)
t42 <- LeafNode$new("t42", interval = th)
c41 <- ChanceNode$new()
e41 <- Reaction$new(
  c41, t41, p = p.Dermatitis.T, cost = c.Dermatitis, label = "Dermatitis"
)
e42 <- Reaction$new(
  c41, t42, p = NA_real_, cost = 0.0, label = "No dermatitis"
)
t43 <- LeafNode$new("t43", interval = th)
t44 <- LeafNode$new("t44", interval = th)
c42 <- ChanceNode$new()
e43 <- Reaction$new(
  c42, t43, p = p.Dermatitis.T, cost = c.Dermatitis, label = "Dermatitis"
)
e44 <- Reaction$new(
  c42, t44, p = NA_real_, cost = 0.0, label = "No dermatitis"
)
c43 <- ChanceNode$new()
e45 <- Reaction$new(c43, c41, p = p.LSI.T, cost = c.LSI, label = "LSI")
e46 <- Reaction$new(c43, c42, p = NA_real_, cost = 0.0, label = "No LSI")
c53 <- ChanceNode$new()
e51 <- Reaction$new(c53, c43, p = p.CRBSI.T, cost = c.CRBSI, label = "CRBSI")
e52 <- Reaction$new(c53, c33, p = NA_real_, cost = 0.0, label = "no CRBSI")

# decision node
d1 <- DecisionNode$new("d1")
e9 <- Action$new(d1, c23, label = "Standard", cost = c.Standard)
e10 <- Action$new(d1, c53, label = "Tegaderm", cost = c.Tegaderm)

# create decision tree
V <- list(
  d1,
  c01, c02, c03, c11, c12, c13, c23, c31, c32, c33, c41, c42, c43, c53,
  t01, t02, t03, t04, t11, t12, t13, t14, t31, t32, t33, t34,
  t41, t42, t43, t44
)
E <- list(
  e01, e02, e03, e04, e05, e06, e11, e12, e13, e14, e15, e16, e21, e22,
  e31, e32, e33, e34, e35, e36, e41, e42, e43, e44, e45, e46, e51, e52,
  e9, e10
)
DT <- DecisionTree$new(V, E)

Tree diagram

The draw method of a DecisionTree object creates a graphical representation of the tree, as follows.

DT$draw(border = TRUE)

Figure 1. Decision tree for the Tegaderm model

Summary of model variables in the tree

The model variables which will be associated with actions, reactions and leaf nodes can be tabulated using the method modvar_table. This returns a data frame describing each variable, its description, units and uncertainty distribution. Variables inheriting from type ModVar will be included in the tabulation unless explicitly excluded, regular numeric values will not be listed. In the Tegaderm model, the input model variables are in the following table, including those constructed from expressions.

#>                          Description                   Distribution
#> 1                           LSI cost                   Ga(11.111,9)
#> 2                  P(LSI | Tegaderm)               p.LSI.S * hr.LSI
#> 3                  P(LSI | Standard)        r.LSI * n.cathdays/1000
#> 4                  Baseline LSI rate                Ga(0.078,1.786)
#> 5             No. days with catheter                      Ga(4,2.5)
#> 6                    Tegaderm LSI HR               LN(-0.997,0.414)
#> 7                    Dermatitis cost                      Ga(4,1.5)
#> 8  P(dermatitis | standard dressing)                   r.Dermatitis
#> 9           Baseline dermatitis risk                      Be(1,475)
#> 10                        CRBSI cost              Ga(10.89,909.091)
#> 11                 P(CRBSI|Tegaderm)           p.CRBSI.S * hr.CRBSI
#> 12      P(CRBSI | standard dressing)      r.CRBSI * n.cathdays/1000
#> 13               Baseline CRBSI rate               Ga(152.111,0.01)
#> 14                 Tegaderm CRBSI HR               LN(-0.997,0.414)
#> 15          P(dermatitis | Tegaderm) p.Dermatitis.S * rr.Dermatitis
#> 16            Tegaderm Dermatitis RR               LN(-0.112,0.472)
#> 17            Standard dressing cost             1.54 * n.dressings
#> 18                     No. dressings                 Ga(2.25,1.333)
#> 19                 Tegaderm CHG cost             6.26 * n.dressings

The units, point estimates, lower 95% and upper 9% confidence intervals are are obtained from the same call, in the remaining columns.

#>                                    Variable     Mean  LowerCI   UpperCI
#> 1                             LSI cost, GBP  100.000   50.124   166.811
#> 2                      P(LSI | Tegaderm), P    0.000    0.000     0.004
#> 3               P(LSI | Standard), /patient    0.002    0.000     0.015
#> 4    Baseline LSI rate, /1000 catheter days    0.140    0.000     1.466
#> 5              No. days with catheter, days   10.000    2.725    21.918
#> 6                       Tegaderm LSI HR, HR    0.402    0.164     0.831
#> 7                      Dermatitis cost, GBP    6.000    1.635    13.151
#> 8      P(dermatitis | standard dressing), P    0.002    0.000     0.008
#> 9       Baseline dermatitis risk, /catheter    0.002    0.000     0.008
#> 10                          CRBSI cost, GBP 9900.000 4921.749 16588.599
#> 11                     P(CRBSI|Tegaderm), P    0.006    0.001     0.017
#> 12          P(CRBSI | standard dressing), P    0.015    0.004     0.033
#> 13 Baseline CRBSI rate, /1000 catheter days    1.480    1.254     1.724
#> 14                    Tegaderm CRBSI HR, HR    0.402    0.164     0.831
#> 15              P(dermatitis | Tegaderm), P    0.002    0.000     0.009
#> 16               Tegaderm Dermatitis RR, RR    1.000    0.354     2.258
#> 17              Standard dressing cost, GBP    4.560    0.733    11.467
#> 18                 No. dressings, dressings    3.000    0.433     7.999
#> 19                   Tegaderm CHG cost, GBP   18.780    2.348    48.878

Running the model

Base case

The following code runs a single model scenario, using the evaluate method of a decision node to evaluate each pathway from the decision node, shown in the table. This model did not consider utility, and the columns associated with utility are removed.

RES <- DT$evaluate()

#>   Run       d1   Cost
#> 1   1 Standard 151.29
#> 2   1 Tegaderm  77.75

Univariate sensitivity analysis

The sensitivity of the decision tree results to each source model variable, varied independently of the others, is demonstrated by a tornado diagram. The method tornado can be used to generate such a plot (and also provides a tabulated version of the values used in the plot).

Tornado diagrams compare outcomes for two interventions, labelled as index and ref. In a decision tree, an intervention is defined as a strategy for traversing the tree, expressed as a list of the action edges emanating from each decision node. In trees with a single decision node, the index and ref parameters may be expressed as a single action edge. Source variables are varied over their 95% confidence limits (figure 2).

to <- DT$tornado(index = e10, ref = e9, draw = TRUE)

Figure 2. Tornado diagram for the Tegaderm model

The object returned from method tornado (to) is a data frame which includes the values of the cost difference when each model variable is univariately at the limits of its 95% confidence interval, as follows:

#>                                    Variable      LL       UL Min.CostDiff
#> 1              No. days with catheter, days    2.72    21.92         9.74
#> 2                           CRBSI cost, GBP 4921.75 16588.60        29.48
#> 3                     Tegaderm CRBSI HR, HR    0.16     0.83       108.45
#> 4                  No. dressings, dressings    0.43     8.00        85.66
#> 5  Baseline CRBSI rate, /1000 catheter days    1.25     1.72        60.17
#> 6    Baseline LSI rate, /1000 catheter days    0.00     1.47        73.46
#> 7                             LSI cost, GBP   50.12   166.81        73.50
#> 8                       Tegaderm LSI HR, HR    0.16     0.83        73.58
#> 9                Tegaderm Dermatitis RR, RR    0.35     2.26        73.55
#> 10                     Dermatitis cost, GBP    1.63    13.15        73.54
#> 11      Baseline dermatitis risk, /catheter    0.00     0.01        73.54
#>    Max.CostDiff
#> 1        178.07
#> 2        132.74
#> 3         10.67
#> 4         49.95
#> 5         88.00
#> 6         74.34
#> 7         73.60
#> 8         73.48
#> 9         73.53
#> 10        73.54
#> 11        73.54

Probabilistic sensitivity analysis

Multivariate probabilistic sensitivity analysis is supported through the use of sampling model variables. The same call, with extra parameters, is used to run the PSA and save the results in a data frame. Additionally, the cost difference is computed for each run of the model, as follows:

N <- 1000L
psa <- DT$evaluate(setvars = "random", by = "run", N = N)
psa[, "Difference"] <- psa[, "Cost.Standard"] - psa[, "Cost.Tegaderm"]

The first few runs of PSA are as follows; the by = "run" option reshapes the table to give one row per simulation, rather than one row per run, per strategy.

#>    Run Cost.Tegaderm Cost.Standard Cost.Difference
#> 1    1        152.54        147.35           -5.18
#> 2    2         99.11        285.56          186.46
#> 3    3         84.31        188.53          104.22
#> 4    4         69.28        112.19           42.91
#> 5    5         54.70        134.01           79.32
#> 6    6         72.99        160.49           87.50
#> 7    7        117.75        110.88           -6.87
#> 8    8         99.77        144.00           44.23
#> 9    9        102.35        140.22           37.87
#> 10  10        112.45        236.65          124.20

From PSA (1000 runs), the mean cost of treatment with Tegaderm was 79.69 GBP, the mean cost of treatment with standard dressings was 150.76 GBP and the mean cost saving was 71.07 GBP. The 95% confidence interval for cost saving was -8.86 GBP to 214.10 GBP; the standard deviation of the cost saving was 58.26 GBP. Overall, 95.7% of runs found that Tegaderm was cost saving. These results replicate those reported by Jenks et al (saving of 72.90 GBP, 97.8% cases cost saving; mean cost of standard dressing 151.29 GBP, mean cost of Tegaderm 77.75 GBP).

Scenario - low baseline rate of CRBSI

Jenks et al modelled an additional scenario, in which the baseline rate of CRBSI was 0.3 per 1000 catheter days (modelled as a Gamma distribution fitted to a sample mean of 0.3 and a sample 95% confidence interval of 0.2 to 0.6). A way to achieve this in rdecision is to replace the model variable for the baseline rate of CRBSI, and any other model variables that depend on it via expressions, and then reconstruct the model, as follows.

r.CRBSI <- GammaModVar$new(
  "Baseline CRBSI rate",  "/1000 catheter days",
  shape = (0.30 ^ 2L) / (0.102 ^ 2L),
  scale = (0.102 ^ 2L) / 0.30
)
p.CRBSI.S <- ExprModVar$new(
  "P(CRBSI | standard dressing)", "P",
  rlang::quo(r.CRBSI * n.cathdays / 1000.0)
)
p.CRBSI.T <- ExprModVar$new(
  "P(CRBSI|Tegaderm)", "P",
  rlang::quo(p.CRBSI.S * hr.CRBSI)
)
e21 <- Reaction$new(c23, c03, p = p.CRBSI.S, cost = c.CRBSI, label = "CRBSI")
e22 <- Reaction$new(c23, c13, p = NA_real_, cost = 0.0, label = "No CRBSI")
e51 <- Reaction$new(c53, c43, p = p.CRBSI.T, cost = c.CRBSI, label = "CRBSI")
e52 <- Reaction$new(c53, c33, p = NA_real_, cost = 0.0, label = "no CRBSI")
E <- list(
  e01, e02, e03, e04, e05, e06, e11, e12, e13, e14, e15, e16, e21, e22,
  e31, e32, e33, e34, e35, e36, e41, e42, e43, e44, e45, e46, e51, e52,
  e9, e10
)
DT <- DecisionTree$new(V, E)

The model for this scenario was run under PSA, as for the base case:

N <- 1000L
psa <- DT$evaluate(setvars = "random", by = "run", N = N)
psa[, "Difference"] <- psa[, "Cost.Standard"] - psa[, "Cost.Tegaderm"]

From PSA (1000 runs), the mean cost of treatment with Tegaderm was 30.12 GBP, the mean cost of treatment with standard dressings was 34.22 GBP and the mean cost saving was 4.10 GBP. The 95% confidence interval for cost saving was -25.43 GBP to 43.51 GBP; the standard deviation of the cost saving was 17.16 GBP. Overall, 57.7% of runs found that Tegaderm was cost saving. These results replicate those reported by Jenks et al (saving of 3.56 GBP, 57.9% cases cost saving; mean cost of standard dressing 34.47 GBP, mean cost of Tegaderm 30.79 GBP).

Two threshold analyses were reported for this scenario. This can be achieved in rdecision by using the threshold method of the decision tree. Firstly, the threshold hazard ratio of a CRBSI with Tegaderm versus a CRBSI with a standard dressing was varied in the range 0.1 to 0.9, as follows:

hr_threshold <- DT$threshold(
  index = list(e10),
  ref = list(e9),
  outcome = "saving",
  mvd = "Tegaderm CRBSI HR",
  a = 0.1,
  b = 0.9,
  tol = 0.01
)

This gave a threshold value of 0.53, above which Tegaderm became cost incurring (the reported threshold was 0.53). Secondly, the cost of each CRBSI was varied between 0 GBP and 9900 GBP to find the threshold of cost saving, as follows:

c_crbsi_threshold <- DT$threshold(
  index = list(e10),
  ref = list(e9),
  outcome = "saving",
  mvd = "CRBSI cost",
  a = 0.0,
  b = 9900.0,
  tol = 10.0
)

This gave a threshold value of 7,840.72 GBP, below which Tegaderm became cost incurring (the reported threshold was 8000 GBP).

References

Jenks, M., Craig, J., Green, W., Hewitt, N., Arber, M. & Sims, A. J. Tegaderm CHG IV Securement Dressing for Central Venous and Arterial Catheter Insertion Sites: A NICE Medical Technology Guidance. Applied Health Economics and Health Policy 14, 135–149 (2016).

Briggs, A., Claxton, K. & Sculpher, M. Decision modelling for health economic evaluation. (Oxford University Press, 2006).

Decision tree with PSA (Jenks 2016)