Nothing
# Script to construct, in rdecision, the model reported by Jenks et al
# (App Health Econ 2016;14:135-149) for a securement dressing for central
# venous and arterial catheters.
#
# The checks are done as part of the testthat framework, which ensures that
# any changes in the package code which unintentionally result in deviations
# from the reported results of the model are identified.
#
# Code to construct and run the model is contained within labelled knitr code
# chunks and does not contain test expectations, so can be used by a vignette.
# Unlabelled code chunks may contain testthat expectations and should be
# ignored by a vignette.
## @knitr variables -----------------------------------------------------------
# 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
)
## @knitr ---------------------------------------------------------------------
test_that("variables have expected values", {
# baseline CRBSI
expect_intol(r.CRBSI$mean(), 1.48, 0.02)
q <- r.CRBSI$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 1.28, 0.05)
expect_intol(q[[2L]], 1.75, 0.05)
# baseline LSI
expect_intol(r.LSI$mean(), 0.14, 0.01)
q <- r.LSI$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 0.0, 0.05)
# baseline dermatitis
expect_intol(r.Dermatitis$mean(), 1L / 476L, 0.0001)
q <- r.Dermatitis$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 0.000, 0.005)
expect_intol(q[[2L]], 0.010, 0.005)
# HR of CRBSI for Tegaderm
expect_intol(hr.CRBSI$mean(), 0.402, 0.010)
q <- hr.CRBSI$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 0.186, 0.05)
expect_intol(q[[2L]], 0.868, 0.05)
# HR of LSI for Tegaderm
expect_intol(hr.LSI$mean(), 0.402, 0.010)
q <- hr.LSI$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 0.186, 0.05)
expect_intol(q[[2L]], 0.868, 0.05)
# RR of dermatitis
expect_intol(rr.Dermatitis$mean(), 1.0, 0.010)
q <- rr.Dermatitis$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 0.35, 0.05)
expect_intol(q[[2L]], 2.26, 0.05)
# cost of CRBSI
expect_intol(c.CRBSI$mean(), 9900.0, 10.0)
q <- c.CRBSI$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 4921.0, 10.0)
expect_intol(q[[2L]], 16589.0, 10.0)
# cost of LSI
expect_intol(c.LSI$mean(), 100.0, 10.0)
q <- c.LSI$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 50.1, 1.0)
expect_intol(q[[2L]], 166.8, 1.0)
# cost of dermatitis
expect_intol(c.Dermatitis$mean(), 6.0, 0.1)
q <- c.Dermatitis$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 1.64, 0.1)
expect_intol(q[[2L]], 13.1, 0.1)
# number of dressings
expect_intol(n.dressings$mean(), 3.0, 0.1)
q <- n.dressings$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 0.4, 0.1)
expect_intol(q[[2L]], 8.0, 0.1)
# number of catheter days
expect_intol(n.cathdays$mean(), 10.0, 0.1)
q <- n.cathdays$quantile(probs = c(0.025, 0.975))
expect_intol(q[[1L]], 2.7, 1.0)
expect_intol(q[[2L]], 21.9, 1.0)
})
## @knitr expressions ---------------------------------------------------------
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)
)
## @knitr tree ----------------------------------------------------------------
# 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)
## @knitr ---------------------------------------------------------------------
test_that("model variables are as expected", {
mv <- DT$modvars()
expect_length(mv, 19L)
MVT <- DT$modvar_table()
expect_identical(nrow(MVT), 19L)
expect_identical(sum(MVT$Est), 8L)
MVT <- DT$modvar_table(FALSE)
expect_identical(nrow(MVT), 11L)
})
## @knitr basecase ------------------------------------------------------------
RES <- DT$evaluate()
## @knitr --------------------------------------------------------------------
test_that("EAC base case agrees with direct calculation", {
# values from Table 4
r_crbsi <- 1.48
r_lsi <- 0.14
r_derm <- 0.0021
hr_crbsi_teg <- 0.402
hr_lsi_teg <- 0.402
rr_derm_teg <- 1.0
c_crbsi <- 9900.0
c_derm <- 6.0
c_lsi <- 100.0
n_cdays <- 10.0
n_dress <- 3L
c_teg <- 6.26
c_std <- 1.54
# probabilities
p_crbsi_std <- r_crbsi * (n_cdays / 1000.0)
p_lsi_std <- r_lsi * (n_cdays / 1000.0)
p_derm_std <- r_derm
p_crbsi_teg <- p_crbsi_std * hr_crbsi_teg
p_lsi_teg <- p_lsi_std * hr_lsi_teg
p_derm_teg <- rr_derm_teg * p_derm_std
# component costs
cdress_std <- c_std * n_dress
cdress_teg <- c_teg * n_dress
ccrbsi_std <- c_crbsi * p_crbsi_std
ccrbsi_teg <- c_crbsi * p_crbsi_teg
clsi_std <- c_lsi * p_lsi_std
clsi_teg <- c_lsi * p_lsi_teg
cderm_std <- c_derm * p_derm_std
cderm_teg <- c_derm * p_derm_teg
# total costs
c_std <- cdress_std + ccrbsi_std + clsi_std + cderm_std
c_teg <- cdress_teg + ccrbsi_teg + clsi_teg + cderm_teg
# check against the model
expect_intol(RES$Cost[RES$d1 == "Standard"], c_std, 2.00)
expect_intol(RES$Cost[RES$d1 == "Tegaderm"], c_teg, 2.00)
# check against the Excel model
expect_intol(RES$Cost[RES$d1 == "Standard"], 151.29, 2.00)
expect_intol(RES$Cost[RES$d1 == "Tegaderm"], 77.75, 2.00)
})
## @knitr --------------------------------------------------------------------
test_that("illegal arguments to evaluate are rejected", {
expect_error(DT$evaluate(setvars = 42L), class = "setvars_not_character")
expect_error(DT$evaluate(setvars = "mode"), class = "setvars_invalid")
expect_silent(DT$evaluate(setvars = "q2.5"))
expect_silent(DT$evaluate(setvars = "q50"))
expect_silent(DT$evaluate(setvars = "q97.5"))
expect_silent(DT$evaluate(setvars = "current"))
})
## @knitr --------------------------------------------------------------------
test_that("tornado method rejects illegal arguments", {
grDevices::pdf(file = NULL)
# tornado
expect_error(DT$tornado(), class = "missing_strategy")
expect_error(
DT$tornado(index = list(e10), ref = list(e52)),
class = "invalid_strategy"
)
expect_error(
DT$tornado(index = list(e10), ref = list(e9), outcome = "survival"),
class = "invalid_outcome"
)
expect_error(
DT$tornado(index = list(e10), ref = list(e9), draw = 42L),
class = "invalid_draw"
)
expect_error(
DT$tornado(
index = list(e10), ref = list(e9), exclude = 42L,
draw = TRUE
),
class = "exclude_not_list"
)
expect_error(
DT$tornado(
index = list(e10), ref = list(e9), exclude = list("SN1", "SN2", "SNX"),
draw = TRUE
),
class = "exclude_element_not_modvar"
)
expect_silent(
DT$tornado(
index = list(e10), ref = list(e9), exclude = list(hr.CRBSI$description()),
draw = TRUE
)
)
to <- DT$tornado(
index = list(e10), ref = list(e9),
draw = TRUE
)
expect_identical(nrow(to), 11L)
grDevices::dev.off()
})
## @knitr PSA -----------------------------------------------------------------
N <- 1000L
psa <- DT$evaluate(setvars = "random", by = "run", N = N)
psa$Difference <- psa$Cost.Standard - psa$Cost.Tegaderm
## @knitr ---------------------------------------------------------------------
test_that("PSA replicates values from Excel model", {
skip_on_cran()
expect_intol(mean(psa$Cost.Standard), 151.29, 10.00)
expect_intol(mean(psa$Cost.Tegaderm), 77.75, 10.00)
expect_intol(mean(psa$Difference), 72.90, 10.00)
expect_intol(sum(psa$Difference > 0.0) / N, 0.978, tol = 0.05)
})
## @knitr scenario -----------------------------------------------------------
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)
## @knitr --------------------------------------------------------------------
test_that("scenario case agrees with direct calculation", {
# evaluate the scenario as a point estimate
s_sco <- DT$evaluate()
# values from Table 4
r_crbsi <- 0.30
r_lsi <- 0.14
r_derm <- 0.0021
hr_crbsi_teg <- 0.402
hr_lsi_teg <- 0.402
rr_derm_teg <- 1.0
c_crbsi <- 9900.0
c_derm <- 6.0
c_lsi <- 100.0
n_cdays <- 10.0
n_dress <- 3L
c_teg <- 6.26
c_std <- 1.54
# probabilities
p_crbsi_std <- r_crbsi * (n_cdays / 1000.0)
p_lsi_std <- r_lsi * (n_cdays / 1000.0)
p_derm_std <- r_derm
p_crbsi_teg <- p_crbsi_std * hr_crbsi_teg
p_lsi_teg <- p_lsi_std * hr_lsi_teg
p_derm_teg <- rr_derm_teg * p_derm_std
# component costs
cdress_std <- c_std * n_dress
cdress_teg <- c_teg * n_dress
ccrbsi_std <- c_crbsi * p_crbsi_std
ccrbsi_teg <- c_crbsi * p_crbsi_teg
clsi_std <- c_lsi * p_lsi_std
clsi_teg <- c_lsi * p_lsi_teg
cderm_std <- c_derm * p_derm_std
cderm_teg <- c_derm * p_derm_teg
# total costs
c_std <- cdress_std + ccrbsi_std + clsi_std + cderm_std
c_teg <- cdress_teg + ccrbsi_teg + clsi_teg + cderm_teg
# check against the model
expect_intol(s_sco$Cost[s_sco$d1 == "Standard"], c_std, 2.00)
expect_intol(s_sco$Cost[s_sco$d1 == "Tegaderm"], c_teg, 2.00)
# check against the excel model
expect_intol(s_sco$Cost[s_sco$d1 == "Standard"], 34.47, 2.00)
expect_intol(s_sco$Cost[s_sco$d1 == "Tegaderm"], 30.79, 2.00)
})
## @knitr scenario-PSA --------------------------------------------------------
N <- 1000L
psa <- DT$evaluate(setvars = "random", by = "run", N = N)
psa$Difference <- psa$Cost.Standard - psa$Cost.Tegaderm
## @knitr ---------------------------------------------------------------------
test_that("scenario PSA agrees with reported results", {
skip_on_cran()
expect_intol(mean(psa$Cost.Standard), 34.47, 5.00)
expect_intol(mean(psa$Cost.Tegaderm), 30.79, 5.00)
expect_intol(mean(psa$Difference), 3.56, 3.00)
expect_intol(sum(psa$Difference > 0.0) / N, 0.579, tol = 0.2)
})
## @knitr threshold-hr -------------------------------------------------------
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
)
## @knitr --------------------------------------------------------------------
test_that("scenario hazard rate threshold agrees with that reported", {
expect_intol(hr_threshold, 0.53, tol = 0.05)
})
## @knitr threshold-ccrbsi ----------------------------------------------------
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
)
## @knitr --------------------------------------------------------------------
test_that("scenario CRBSI cost threshold agrees with reported value", {
expect_intol(c_crbsi_threshold, 8000.0, tol = 300.0)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.