Nothing
# Script to construct, in rdecision, the model described by Dong and Buxton (Int
# J Health Tech Assessment 2006;22:191-202) for computer-assisted total knee
# replacement, and to check that its results agree with those presented in the
# paper.
#
# The checks are done as part of the testthat framework, ensuring that
# changes in the package code which unintentionally result in deviations
# from the expected results of the model are identified.
#
# Code to construct and run the model is contained within labelled knitr code
# chunks and do 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 state-names ---------------------------------------------------------
states <- c(
A = "TKR operation for knee problems",
B = "TKR with serious complications",
C = "TKR with minor complications",
D = "Normal health after primary TKR",
E = "Complex revision",
F = "Simple revision",
G = "Other treatments",
H = "Normal health after TKR revision",
I = "Death"
)
## @knitr utils-point ----------------------------------------------------------
utility_A <- 0.72
utility_B <- 0.35
utility_C <- 0.66
utility_D <- 0.78
utility_E <- 0.51
utility_F <- 0.66
utility_G <- 0.72
utility_H <- 0.68
utility_I <- 0.00
## @knitr costs-point ----------------------------------------------------------
cost_A <- 5197.0
cost_B <- 0.0
cost_C <- 0.0
cost_D <- 0.0
cost_E <- 7326.0
cost_F <- 6234.0
cost_G <- 2844.0
cost_H <- 0.0
cost_I <- 0.0
cost_CAS <- 235.0
## @knitr SMM-point ------------------------------------------------------------
# Markov states
sA <- MarkovState$new(states["A"], utility = utility_A)
sB <- MarkovState$new(states["B"], utility = utility_B)
sC <- MarkovState$new(states["C"], utility = utility_C)
sD <- MarkovState$new(states["D"], utility = utility_D)
sE <- MarkovState$new(states["E"], utility = utility_E)
sF <- MarkovState$new(states["F"], utility = utility_F)
sG <- MarkovState$new(states["G"], utility = utility_G)
sH <- MarkovState$new(states["H"], utility = utility_H)
sI <- MarkovState$new(states["I"], utility = utility_I)
States <- list(sA, sB, sC, sD, sE, sF, sG, sH, sI)
# Transitions
tAD <- Transition$new(sA, sD, cost = cost_A)
tAC <- Transition$new(sA, sC, cost = cost_A)
tAB <- Transition$new(sA, sB, cost = cost_A)
tBC <- Transition$new(sB, sC)
tBE <- Transition$new(sB, sE, cost = cost_E)
tBF <- Transition$new(sB, sF, cost = cost_F)
tBG <- Transition$new(sB, sG, cost = cost_G)
tCB <- Transition$new(sC, sB)
tCD <- Transition$new(sC, sD)
tCF <- Transition$new(sC, sF, cost = cost_F)
tCG <- Transition$new(sC, sG, cost = cost_G)
tCC <- Transition$new(sC, sC)
tDC <- Transition$new(sD, sC)
tDB <- Transition$new(sD, sB)
tDD <- Transition$new(sD, sD)
tEB <- Transition$new(sE, sB)
tEH <- Transition$new(sE, sH)
tFB <- Transition$new(sF, sB)
tFC <- Transition$new(sF, sC)
tFG <- Transition$new(sF, sG, cost = cost_G)
tFH <- Transition$new(sF, sH)
tGB <- Transition$new(sG, sB)
tGC <- Transition$new(sG, sC)
tGF <- Transition$new(sG, sF, cost = cost_F)
tGD <- Transition$new(sG, sD)
tHE <- Transition$new(sH, sE, cost = cost_E)
tHF <- Transition$new(sH, sF, cost = cost_F)
tHH <- Transition$new(sH, sH)
tBI <- Transition$new(sB, sI)
tCI <- Transition$new(sC, sI)
tDI <- Transition$new(sD, sI)
tEI <- Transition$new(sE, sI)
tFI <- Transition$new(sF, sI)
tGI <- Transition$new(sG, sI)
tHI <- Transition$new(sH, sI)
tII <- Transition$new(sI, sI)
# Transitions incorporating CAS cost
tAD_CAS <- Transition$new(sA, sD, cost = cost_A + cost_CAS)
tAC_CAS <- Transition$new(sA, sC, cost = cost_A + cost_CAS)
tAB_CAS <- Transition$new(sA, sB, cost = cost_A + cost_CAS)
tBE_CAS <- Transition$new(sB, sE, cost = cost_E + cost_CAS)
tBF_CAS <- Transition$new(sB, sF, cost = cost_F + cost_CAS)
tCF_CAS <- Transition$new(sC, sF, cost = cost_F + cost_CAS)
tGF_CAS <- Transition$new(sG, sF, cost = cost_F + cost_CAS)
tHE_CAS <- Transition$new(sH, sE, cost = cost_E + cost_CAS)
tHF_CAS <- Transition$new(sH, sF, cost = cost_F + cost_CAS)
Transitions_base <- list(
tAD, tAC, tAB, tBC, tBE, tBF, tBG, tCB, tCD, tCF, tCG, tCC,
tDC, tDB, tDD, tEB, tEH, tFB, tFC, tFG, tFH, tGB, tGC, tGF,
tGD, tHE, tHF, tHH, tBI, tCI, tDI, tEI, tFI, tGI, tHI, tII
)
Transitions_CAS <- list(
tAD_CAS, tAC_CAS, tAB_CAS, tBC, tBE_CAS, tBF_CAS, tBG, tCB, tCD, tCF_CAS,
tCG, tCC, tDC, tDB, tDD, tEB, tEH, tFB, tFC, tFG, tFH, tGB, tGC, tGF_CAS,
tGD, tHE_CAS, tHF_CAS, tHH, tBI, tCI, tDI, tEI, tFI, tGI, tHI, tII
)
## @knitr SMM-def-point --------------------------------------------------------
SMM_base <- SemiMarkovModel$new(
V = States, E = Transitions_base,
tcycle = as.difftime(365.25 / 12L, units = "days"), # cycles in months
discount.cost = 0.035,
discount.utility = 0.035
)
SMM_CAS <- SemiMarkovModel$new(
V = States, E = Transitions_CAS,
tcycle = as.difftime(365.25 / 12L, units = "days"), # cycles in months
discount.cost = 0.035,
discount.utility = 0.035
)
## @knitr transitions-point ---------------------------------------------------
# Death
p_death_all <- 0.00341
p_death_primary <- 0.00046
p_death_revision <- 0.00151
# Transitions
p_AtoB <- 0.01495
p_AtoC <- 0.04285
p_AtoD <- NA # alternatively: 1 - p_AtoB - p_AtoC
p_BtoC <- 0.01385
p_BtoE <- 0.02469
p_BtoF <- 0.00523
p_BtoI <- p_death_all + p_death_primary
p_BtoG <- NA # alternatively: 1 - p_BtoC - p_BtoE - p_BtoF - p_BtoI
p_CtoB <- 0.00921
p_CtoC <- 0.02505
p_CtoF <- 0.00250
p_CtoG <- 0.01701
p_CtoI <- p_death_all + p_death_primary
p_CtoD <- NA # alternatively: 1 - p_CtoB - p_CtoC - p_CtoF - p_CtoG - p_CtoI
p_DtoB <- 0.00921
p_DtoC <- 0.01385
p_DtoI <- p_death_all + p_death_primary
p_DtoD <- NA # alternatively: 1 - p_DtoB - p_DtoC - p_DtoI
p_EtoB <- 0.02545
p_EtoI <- p_death_all + p_death_revision
p_EtoH <- NA # alternatively: 1 - p_EtoB - p_EtoI
p_FtoB <- 0.01590
p_FtoC <- 0.00816
p_FtoG <- 0.01701
p_FtoI <- p_death_all + p_death_revision
p_FtoH <- NA # alternatively: 1 - p_FtoB - p_FtoC - p_FtoG - p_FtoI
p_GtoB <- 0.00921
p_GtoC <- 0.01385
p_GtoF <- 0.00250
p_GtoI <- p_death_all + p_death_primary
p_GtoD <- NA # alternatively: 1 - p_GtoB - p_GtoC - p_GtoF - p_GtoI
p_HtoE <- 0.02003
p_HtoF <- 0.01038
p_HtoI <- p_death_all + p_death_revision
p_HtoH <- NA # alternatively: 1 - p_HtoE - p_HtoF - p_HtoI
p_ItoI <- 1.0
# Set transition probabilities
Pt <- matrix(c(
0L, p_AtoB, p_AtoC, p_AtoD, 0L, 0L, 0L, 0L, 0L,
0L, 0L, p_BtoC, 0L, p_BtoE, p_BtoF, p_BtoG, 0L, p_BtoI,
0L, p_CtoB, p_CtoC, p_CtoD, 0L, p_CtoF, p_CtoG, 0L, p_CtoI,
0L, p_DtoB, p_DtoC, p_DtoD, 0L, 0L, 0L, 0L, p_DtoI,
0L, p_EtoB, 0L, 0L, 0L, 0L, 0L, p_EtoH, p_EtoI,
0L, p_FtoB, p_FtoC, 0L, 0L, 0L, p_FtoG, p_FtoH, p_FtoI,
0L, p_GtoB, p_GtoC, p_GtoD, 0L, p_GtoF, 0L, 0L, p_GtoI,
0L, 0L, 0L, 0L, p_HtoE, p_HtoF, 0L, p_HtoH, p_HtoI,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, p_ItoI
), nrow = 9L, byrow = TRUE, dimnames = list(
source = states[LETTERS[1L:9L]], target = states[LETTERS[1L:9L]]
))
SMM_base$set_probabilities(Pt)
## @knitr tx-effect -----------------------------------------------------------
txeffect <- function(Pt, rr) {
# copy transition matrix
pr <- Pt
# calculate reduced probability of transition to state B
derisk <- states[["B"]]
dpb <- pr[, derisk] * rr
# reduce the probability of making a transition to state B and increase the
# probability of making a transition elsewhere by the same amount
uprisks <- c("D", "G", "D", "D", "H", "H", "D", "H", "I")
for (i in 1L:9L) {
s <- states[[i]]
pr[[s, derisk]] <- pr[[s, derisk]] - dpb[[i]]
uprisk <- states[[uprisks[[i]]]]
pr[[s, uprisk]] <- pr[[s, uprisk]] + dpb[[i]]
}
return(pr)
}
## @knitr CAS-transitions-point ------------------------------------------------
# apply CAS_effect to the transition matrix
CAS_effect <- 0.34
Pt_CAS <- txeffect(Pt, CAS_effect)
SMM_CAS$set_probabilities(Pt_CAS)
## @knitr cycle-point ----------------------------------------------------------
# create starting populations
N <- 1000L
populations <- c(N, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)
names(populations) <- states
# run 120 one-month cycles for both models
SMM_base$reset(populations)
SMM_base_10years <- SMM_base$cycles(
ncycles = 120L, hcc.pop = FALSE, hcc.cost = FALSE
)
SMM_CAS$reset(populations)
SMM_CAS_10years <- SMM_CAS$cycles(
ncycles = 120L, hcc.pop = FALSE, hcc.cost = FALSE
)
## @knitr as-table4 ----------------------------------------------------------
# convert a Markov trace (matrix) with monthly cycles to an annual summary
# matrix with cumulative values, as per Dong and Buxton, Table 4
as_table4 <- function(m_trace) {
# calculate cumulative event counts
m_cum <- apply(m_trace, 2L, cumsum)
# create annual summary table
m_t4 <- matrix(
data = c(
m_trace[, "Years"],
m_cum[, "TKR with serious complications"] / 10L,
m_cum[, "TKR with minor complications"] / 10L,
m_cum[, "Complex revision"] / 10L,
m_cum[, "Simple revision"] / 10L,
m_trace[, "Death"] / 10L,
m_cum[, "Cost"] * 1000L,
m_cum[, "QALY"] * 1000L
),
dimnames = list(
NULL,
c(
"Year",
"Cumulative serious complication (%)",
"Cumulative minor complication (%)",
"Cumulative complex revision (%)",
"Cumulative simple revision (%)",
"Cumulative death (%)",
"Discounted costs (£)",
"Discounted QALYs"
)
),
nrow = 121L, ncol = 8L
)
# return in table 4 format
yearly <- (1L:10L) * 12L + 1L
return(m_t4[yearly, ])
}
## @knitr cs-table-4 ----------------------------------------------------------
t4_CS <- as_table4(SMM_base_10years)
## @knitr ---------------------------------------------------------------------
test_that("Table 4 for conventional surgery is replicated", {
expect_identical(nrow(t4_CS), 10L)
expect_intol(
t4_CS[[10L, "Cumulative serious complication (%)"]], 87.32, tolerance = 0.2
)
expect_intol(
t4_CS[[10L, "Cumulative minor complication (%)"]], 135.87, tolerance = 0.2
)
expect_intol(
t4_CS[[10L, "Cumulative complex revision (%)"]], 5.13, tolerance = 0.2
)
expect_intol(
t4_CS[[10L, "Cumulative simple revision (%)"]], 2.55, tolerance = 0.2
)
expect_intol(
t4_CS[[10L, "Cumulative death (%)"]], 37.10, tolerance = 0.2
)
expect_intol(
t4_CS[[10L, "Discounted costs (£)"]], 7791288.9, tolerance = 100000.0
)
expect_intol(
t4_CS[[10L, "Discounted QALYs"]], 5526.4, tolerance = 250.0
)
})
## @knitr cas-table-4 ----------------------------------------------------------
t4_CAS <- as_table4(SMM_CAS_10years)
## @knitr ---------------------------------------------------------------------
test_that("Table 4 for computer-assisted surgery is replicated", {
expect_identical(nrow(t4_CAS), 10L)
expect_intol(
t4_CAS[[10L, "Cumulative serious complication (%)"]], 58.14, tolerance = 0.2
)
expect_intol(
t4_CAS[[10L, "Cumulative minor complication (%)"]], 136.52, tolerance = 0.2
)
expect_intol(
t4_CAS[[10L, "Cumulative complex revision (%)"]], 3.56, tolerance = 0.2
)
expect_intol(
t4_CAS[[10L, "Cumulative simple revision (%)"]], 1.89, tolerance = 0.2
)
expect_intol(
t4_CAS[[10L, "Cumulative death (%)"]], 37.06, tolerance = 0.2
)
expect_intol(
t4_CAS[[10L, "Discounted costs (£)"]], 7208671.4, tolerance = 100000.0
)
expect_intol(
t4_CAS[10L, "Discounted QALYs"], 5541.2, tolerance = 250.0
)
})
## @knitr cea -----------------------------------------------------------------
dcost <- t4_CAS[[10L, "Discounted costs (£)"]] / 1000L -
t4_CS[[10L, "Discounted costs (£)"]] / 1000L
dutil <- t4_CAS[[10L, "Discounted QALYs"]] / 1000L -
t4_CS[[10L, "Discounted QALYs"]] / 1000L
## @knitr ---------------------------------------------------------------------
test_that("Deterministic CEA matches published value", {
expect_intol(dcost, -583.0, tolerance = 50.0)
expect_intol(dutil, 0.0164, tolerance = 0.005)
})
## @knitr utilities-var ------------------------------------------------------
utility_A_nu <- 0.42
utility_B_nu <- 0.80
utility_C_nu <- 0.32
utility_D_nu <- 0.34
utility_E_nu <- 0.55
utility_F_nu <- 0.57
utility_G_nu <- 0.34
utility_H_nu <- 0.38
## @knitr utilities-beta -----------------------------------------------------
utility_A_beta <- BetaModVar$new(
"Utility of state A", "",
utility_A * utility_A_nu, (1.0 - utility_A) * utility_A_nu
)
utility_B_beta <- BetaModVar$new(
"Utility of state B", "",
utility_B * utility_B_nu, (1.0 - utility_B) * utility_B_nu
)
utility_C_beta <- BetaModVar$new(
"Utility of state C", "",
utility_C * utility_C_nu, (1.0 - utility_C) * utility_C_nu
)
utility_D_beta <- BetaModVar$new(
"Utility of state D", "",
utility_D * utility_D_nu, (1.0 - utility_D) * utility_D_nu
)
utility_E_beta <- BetaModVar$new(
"Utility of state E", "",
utility_E * utility_E_nu, (1.0 - utility_E) * utility_E_nu
)
utility_F_beta <- BetaModVar$new(
"Utility of state F", "",
utility_F * utility_F_nu, (1.0 - utility_F) * utility_F_nu
)
utility_G_beta <- BetaModVar$new(
"Utility of state G", "",
utility_G * utility_G_nu, (1.0 - utility_G) * utility_G_nu
)
utility_H_beta <- BetaModVar$new(
"Utility of state H", "",
utility_H * utility_H_nu, (1.0 - utility_H) * utility_H_nu
)
## @knitr costs-gamma --------------------------------------------------------
cost_A_stdev <- (6217L - 4218L) / (2L * 1.96)
cost_E_stdev <- (11307L - 5086L) / (2L * 1.96)
cost_F_stdev <- (7972L - 5043L) / (2L * 1.96)
cost_G_stdev <- (5579L - 1428L) / (2L * 1.96)
cost_A_gamma <- GammaModVar$new(
"Cost of state A", "", cost_A ^ 2L / cost_A_stdev ^ 2L,
cost_A_stdev ^ 2L / cost_A
)
cost_E_gamma <- GammaModVar$new(
"Cost of state E", "", cost_E ^ 2L / cost_E_stdev ^ 2L,
cost_E_stdev ^ 2L / cost_E
)
cost_F_gamma <- GammaModVar$new(
"Cost of state F", "", cost_F ^ 2L / cost_F_stdev ^ 2L,
cost_F_stdev ^ 2L / cost_F
)
cost_G_gamma <- GammaModVar$new(
"Cost of state G", "", cost_G ^ 2L / cost_G_stdev ^ 2L,
cost_G_stdev ^ 2L / cost_G
)
## @knitr cas-cost-gamma ------------------------------------------------------
cost_CAS_stdev <- 4L * cost_A_stdev / cost_A * cost_CAS
cost_CAS_gamma <- GammaModVar$new(
"Cost of CAS", "", cost_CAS ^ 2L / cost_CAS_stdev ^ 2L,
cost_CAS_stdev ^ 2L / cost_CAS
)
cost_A_CAS <- ExprModVar$new(
"Cost of state A with CAS", "", rlang::quo(cost_A_gamma + cost_CAS_gamma)
)
cost_E_CAS <- ExprModVar$new(
"Cost of state E with CAS", "", rlang::quo(cost_E_gamma + cost_CAS_gamma)
)
cost_F_CAS <- ExprModVar$new(
"Cost of state F with CAS", "", rlang::quo(cost_F_gamma + cost_CAS_gamma)
)
## @knitr cas-lognorm ---------------------------------------------------------
CAS_effect_mean <- 0.34
CAS_effect_sd <- 1.25
CAS_effect_lognorm <- LogNormModVar$new(
"Effect of CAS", "", log(CAS_effect_mean), log(CAS_effect_sd)
)
## @knitr SMM-PSA -------------------------------------------------------------
# Markov states as Beta distributions
sA_PSA <- MarkovState$new(states["A"], utility = utility_A_beta)
sB_PSA <- MarkovState$new(states["B"], utility = utility_B_beta)
sC_PSA <- MarkovState$new(states["C"], utility = utility_C_beta)
sD_PSA <- MarkovState$new(states["D"], utility = utility_D_beta)
sE_PSA <- MarkovState$new(states["E"], utility = utility_E_beta)
sF_PSA <- MarkovState$new(states["F"], utility = utility_F_beta)
sG_PSA <- MarkovState$new(states["G"], utility = utility_G_beta)
sH_PSA <- MarkovState$new(states["H"], utility = utility_H_beta)
# state I has no uncertainty associated with it, so a probabilistic
# representation is not required
States_PSA <- list(
sA_PSA, sB_PSA, sC_PSA, sD_PSA, sE_PSA, sF_PSA, sG_PSA, sH_PSA, sI
)
# Transition costs as Gamma distributions
tAD_PSA <- Transition$new(sA_PSA, sD_PSA, cost = cost_A_gamma)
tAC_PSA <- Transition$new(sA_PSA, sC_PSA, cost = cost_A_gamma)
tAB_PSA <- Transition$new(sA_PSA, sB_PSA, cost = cost_A_gamma)
tBC_PSA <- Transition$new(sB_PSA, sC_PSA)
tBE_PSA <- Transition$new(sB_PSA, sE_PSA, cost = cost_E_gamma)
tBF_PSA <- Transition$new(sB_PSA, sF_PSA, cost = cost_F_gamma)
tBG_PSA <- Transition$new(sB_PSA, sG_PSA, cost = cost_G_gamma)
tCB_PSA <- Transition$new(sC_PSA, sB_PSA)
tCD_PSA <- Transition$new(sC_PSA, sD_PSA)
tCF_PSA <- Transition$new(sC_PSA, sF_PSA, cost = cost_F_gamma)
tCG_PSA <- Transition$new(sC_PSA, sG_PSA, cost = cost_G_gamma)
tCC_PSA <- Transition$new(sC_PSA, sC_PSA)
tDC_PSA <- Transition$new(sD_PSA, sC_PSA)
tDB_PSA <- Transition$new(sD_PSA, sB_PSA)
tDD_PSA <- Transition$new(sD_PSA, sD_PSA)
tEB_PSA <- Transition$new(sE_PSA, sB_PSA)
tEH_PSA <- Transition$new(sE_PSA, sH_PSA)
tFB_PSA <- Transition$new(sF_PSA, sB_PSA)
tFC_PSA <- Transition$new(sF_PSA, sC_PSA)
tFG_PSA <- Transition$new(sF_PSA, sG_PSA, cost = cost_G_gamma)
tFH_PSA <- Transition$new(sF_PSA, sH_PSA)
tGB_PSA <- Transition$new(sG_PSA, sB_PSA)
tGC_PSA <- Transition$new(sG_PSA, sC_PSA)
tGF_PSA <- Transition$new(sG_PSA, sF_PSA, cost = cost_F_gamma)
tGD_PSA <- Transition$new(sG_PSA, sD_PSA)
tHE_PSA <- Transition$new(sH_PSA, sE_PSA, cost = cost_E_gamma)
tHF_PSA <- Transition$new(sH_PSA, sF_PSA, cost = cost_F_gamma)
tHH_PSA <- Transition$new(sH_PSA, sH_PSA)
tBI_PSA <- Transition$new(sB_PSA, sI)
tCI_PSA <- Transition$new(sC_PSA, sI)
tDI_PSA <- Transition$new(sD_PSA, sI)
tEI_PSA <- Transition$new(sE_PSA, sI)
tFI_PSA <- Transition$new(sF_PSA, sI)
tGI_PSA <- Transition$new(sG_PSA, sI)
tHI_PSA <- Transition$new(sH_PSA, sI)
# transition I to I also has no probabilistic elements
# Transitions incorporating CAS cost
tAD_CAS_PSA <- Transition$new(sA_PSA, sD_PSA, cost = cost_A_CAS)
tAC_CAS_PSA <- Transition$new(sA_PSA, sC_PSA, cost = cost_A_CAS)
tAB_CAS_PSA <- Transition$new(sA_PSA, sB_PSA, cost = cost_A_CAS)
tBE_CAS_PSA <- Transition$new(sB_PSA, sE_PSA, cost = cost_E_CAS)
tBF_CAS_PSA <- Transition$new(sB_PSA, sF_PSA, cost = cost_F_CAS)
tCF_CAS_PSA <- Transition$new(sC_PSA, sF_PSA, cost = cost_F_CAS)
tGF_CAS_PSA <- Transition$new(sG_PSA, sF_PSA, cost = cost_F_CAS)
tHE_CAS_PSA <- Transition$new(sH_PSA, sE_PSA, cost = cost_E_CAS)
tHF_CAS_PSA <- Transition$new(sH_PSA, sF_PSA, cost = cost_F_CAS)
Transitions_base_PSA <- list(
tAD_PSA, tAC_PSA, tAB_PSA,
tBC_PSA, tBE_PSA, tBF_PSA, tBG_PSA, tBI_PSA,
tCB_PSA, tCD_PSA, tCF_PSA, tCG_PSA, tCC_PSA, tCI_PSA,
tDC_PSA, tDB_PSA, tDD_PSA, tDI_PSA,
tEB_PSA, tEH_PSA, tEI_PSA,
tFB_PSA, tFC_PSA, tFG_PSA, tFH_PSA, tFI_PSA,
tGB_PSA, tGC_PSA, tGF_PSA, tGD_PSA, tGI_PSA,
tHE_PSA, tHF_PSA, tHH_PSA, tHI_PSA,
tII
)
Transitions_CAS_PSA <- list(
tAD_CAS_PSA, tAC_CAS_PSA, tAB_CAS_PSA,
tBC_PSA, tBE_CAS_PSA, tBF_CAS_PSA, tBG_PSA, tBI_PSA,
tCB_PSA, tCD_PSA, tCF_CAS_PSA, tCG_PSA, tCC_PSA, tCI_PSA,
tDC_PSA, tDB_PSA, tDD_PSA, tDI_PSA,
tEB_PSA, tEH_PSA, tEI_PSA,
tFB_PSA, tFC_PSA, tFG_PSA, tFH_PSA, tFI_PSA,
tGB_PSA, tGC_PSA, tGF_CAS_PSA, tGD_PSA, tGI_PSA,
tHE_CAS_PSA, tHF_CAS_PSA, tHH_PSA, tHI_PSA,
tII
)
SMM_base_PSA <- SemiMarkovModel$new(
V = States_PSA, E = Transitions_base_PSA,
tcycle = as.difftime(365.25 / 12L, units = "days"), # cycles in months
discount.cost = 0.035,
discount.utility = 0.035
)
SMM_CAS_PSA <- SemiMarkovModel$new(
V = States_PSA, E = Transitions_CAS_PSA,
tcycle = as.difftime(365.25 / 12L, units = "days"), # cycles in months
discount.cost = 0.035,
discount.utility = 0.035
)
## @knitr fun-dirichlet -------------------------------------------------------
dirichletify <- function(Pt, population = 1L) {
# check argument
stopifnot(
is.matrix(Pt),
is.numeric(Pt),
nrow(Pt) == ncol(Pt),
all(dimnames(Pt)[[2L]] == dimnames(Pt)[[1L]])
)
nNA <- rowSums(is.na(Pt))
sumP <- rowSums(Pt, na.rm = TRUE)
# check if a count or proportion representation
is_count <- any(Pt > 1.0, na.rm = TRUE)
if (is_count) {
# counts cannot have NAs
stopifnot(all(nNA == 0L))
# normalise into proportions
Pt <- Pt / rowSums(Pt)
} else {
# proportions can have NA values, but only 1/row
stopifnot(all(nNA <= 1L))
# store information about which variable was set as NA
whichNA <- apply(Pt, 1L, function(row) which(is.na(row)))
# populate missing values to a total of 1
for (row in names(nNA)[nNA > 0L]) {
Pt[row, whichNA[[row]]] <- 1.0 - sumP[row]
}
}
# build state-wise Dirichlet distributions
for (r in seq_len(nrow(Pt))) {
non0 <- which(Pt[r, ] != 0.0)
# if multiple outgoing transitions are possible, model as Dirichlet
if (length(non0) > 1L) {
dist <- DirichletDistribution$new(Pt[r, non0] * population)
dist$sample() # randomise
Pt[r, non0] <- dist$r()
}
# if only 1 transition is possible, leave as given originally
}
return(Pt)
}
## @knitr cycle-PSA -----------------------------------------------------------
nruns <- 250L
t4_CS_PSA <- array(
dim = c(10L, ncol(t4_CS), nruns),
dimnames = list(NULL, colnames(t4_CS), NULL)
)
t4_CAS_PSA <- array(
dim = c(10L, ncol(t4_CS), nruns),
dimnames = list(NULL, colnames(t4_CAS), NULL)
)
for (run in seq_len(nruns)) {
# reset populations
SMM_base_PSA$reset(populations)
SMM_CAS_PSA$reset(populations)
# find unique modvars and randomise them
mv <- unique(c(
SMM_base_PSA$modvars(),
SMM_CAS_PSA$modvars(),
CAS_effect_lognorm
))
for (m in mv) {
m$set("random")
}
# set the transition matrix, applying the CAS effect for CAS model
pt_cs <- dirichletify(Pt, population = 1000L)
SMM_base_PSA$set_probabilities(pt_cs)
pt_cas <- txeffect(pt_cs, CAS_effect_lognorm$get())
SMM_CAS_PSA$set_probabilities(pt_cas)
# cycle the CS model
mtrace <- SMM_base_PSA$cycles(
ncycles = 120L, hcc.pop = FALSE, hcc.cost = FALSE
)
t4 <- as_table4(as.matrix(mtrace))
t4_CS_PSA[, , run] <- t4
# cycle the CAS model
mtrace <- SMM_CAS_PSA$cycles(
ncycles = 120L, hcc.pop = FALSE, hcc.cost = FALSE
)
t4 <- as_table4(as.matrix(mtrace))
t4_CAS_PSA[, , run] <- t4
}
## @knitr ---------------------------------------------------------------------
test_that("Table 4 for conventional surgery is replicated by mean of PSA", {
# skip on CRAN
skip_on_cran()
# Test that the range of the PSA estimates includes the expected value for
# seven outcomes. This is a crude test of functionality. Statistical tests,
# such as a one sample t-test, are not much help here because the shape of
# the distribution of model results may be non-normal, and the type 1 error
# rate is likely to be higher than suggested by the p-value.
fields <- c(
"Cumulative serious complication (%)",
"Cumulative minor complication (%)",
"Cumulative complex revision (%)",
"Cumulative simple revision (%)",
"Cumulative death (%)",
"Discounted costs (£)",
"Discounted QALYs"
)
for (f in fields) {
r <- range(t4_CS_PSA[10L, f, ])
e <- t4_CS[[10L, f]]
expect_between(e, r[[1L]], r[[2L]])
}
})
## @knitr ---------------------------------------------------------------------
test_that("Table 4 for computer-assisted surgery is replicated by mean PSA", {
# skip on CRAN and set wide tolerance
skip_on_cran()
# Test that the range of the PSA estimates includes the expected value for
# seven outcomes. This is a crude test of functionality. Statistical tests,
# such as a one sample t-test, are not much help here because the shape of
# the distribution of model results may be non-normal, and the type 1 error
# rate is likely to be higher than suggested by the p-value.
fields <- c(
"Cumulative serious complication (%)",
"Cumulative minor complication (%)",
"Cumulative complex revision (%)",
"Cumulative simple revision (%)",
"Cumulative death (%)",
"Discounted costs (£)",
"Discounted QALYs"
)
for (f in fields) {
r <- range(t4_CAS_PSA[10L, f, ])
e <- t4_CAS[[10L, f]]
expect_between(e, r[[1L]], r[[2L]])
}
})
## @knitr t5-CS-PSA ----------------------------------------------------------
fields <- c(
"Cumulative serious complication (%)",
"Cumulative complex revision (%)",
"Cumulative simple revision (%)"
)
t5_CS <- matrix(
data = NA_real_, nrow = length(fields), ncol = 4L,
dimnames = list(fields, c("Mean", "SD", "Q2.5", "Q97.5"))
)
for (f in fields) {
t5_CS[[f, "Mean"]] <- mean(t4_CS_PSA[10L, f, ])
t5_CS[[f, "SD"]] <- sd(t4_CS_PSA[10L, f, ])
t5_CS[[f, "Q2.5"]] <- quantile(t4_CS_PSA[10L, f, ], probs = 0.025)
t5_CS[[f, "Q97.5"]] <- quantile(t4_CS_PSA[10L, f, ], probs = 0.975)
}
## @knitr t5-CAS-PSA ----------------------------------------------------------
fields <- c(
"Cumulative serious complication (%)",
"Cumulative complex revision (%)",
"Cumulative simple revision (%)"
)
t5_CAS <- matrix(
data = NA_real_, nrow = length(fields), ncol = 4L,
dimnames = list(fields, c("Mean", "SD", "Q2.5", "Q97.5"))
)
for (f in fields) {
t5_CAS[[f, "Mean"]] <- mean(t4_CAS_PSA[10L, f, ])
t5_CAS[[f, "SD"]] <- sd(t4_CAS_PSA[10L, f, ])
t5_CAS[[f, "Q2.5"]] <- quantile(t4_CAS_PSA[10L, f, ], probs = 0.025)
t5_CAS[[f, "Q97.5"]] <- quantile(t4_CAS_PSA[10L, f, ], probs = 0.975)
}
## @knitr cea-psa --------------------------------------------------------------
dcost_psa <- t4_CAS_PSA[10L, "Discounted costs (£)", ] / 1000L -
t4_CS_PSA[10L, "Discounted costs (£)", ] / 1000L
dutil_psa <- t4_CAS_PSA[10L, "Discounted QALYs", ] / 1000L -
t4_CS_PSA[10L, "Discounted QALYs", ] / 1000L
icer_psa <- dcost_psa / dutil_psa
## @knitr ---------------------------------------------------------------------
test_that("PSA CEA matches deterministic calculation", {
r <- range(dcost_psa)
expect_between(dcost, r[[1L]], r[[2L]])
r <- range(dutil_psa)
expect_between(dutil, r[[1L]], r[[2L]])
})
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.