---
title: "Markov model with probabilistic sensitivity analysis"
subtitle: "Computer-assisted total knee replacement"
author: "Paola Cognigni and Andrew Sims"
date: "February 2023"
bibliography: "REFERENCES.bib"
csl: "nature-no-et-al.csl"
output:
rmarkdown::html_vignette:
fig_width: 7
fig_height: 5
fig_caption: true
df_print: kable
vignette: >
%\VignetteIndexEntry{Markov model with probabilistic sensitivity analysis}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r}
#| purl = FALSE,
#| include = FALSE
knitr::opts_chunk$set(
collapse = TRUE,
echo = FALSE,
comment = "#>"
)
```
```{r}
#| purl = FALSE
#nolint start
```
```{r}
library(rdecision)
```
```{r}
#| purl = FALSE
#nolint end
```
The article by Dong and Buxton, 2006 [-@dong2006] describes a Markov model for
the assessment of a new Computer-Assisted Surgery (CAS) technology for Total
Knee Replacement (TKR). This vignette follows their model to calculate the
cost-effectiveness of conventional and computer-assisted TKR surgery using
`rdecision`.
The authors present both a point estimate and a probabilistic interpretation.
Both approaches are shown to highlight the similarities and differences when
using `rdecision`.
# Model description
## States
The patient journey is represented by 9 possible states, with each patient
assigned to exactly one state at any time:
```{r}
#| echo = TRUE
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"
)
```
## Variables
To assess cost effectiveness, the model must incorporate utilities and costs
for each state.
**Utilities** apply in a continuous manner: the longer the time a patient
spends in a particular state, the longer a particular utility value applies.
This is represented in `rdecision` by assigning the utility to a particular
`MarkovState`. This ensures that they are correctly scaled to different cycle
lengths, and that the annual discount rate is applied correctly.
```{r}
#| echo = TRUE
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
```
**Costs**, in this case, represent one-time expenses associated with the
delivery of a surgical procedure, such as a TKR operation or a revision. As
such, they do not depend on the time the patient spends in a particular state,
and are represented by a cost assigned to a `Transition`. Ongoing costs that
apply as long as the patient remains in a state, for example medication, would
instead be assigned to the relevant `MarkovState`. The costs associated with
revisions/further treatments (states **E**, **F** and **G**) can be assigned to
the transitions from any other state into these states. For the cost of the
primary TKR operation (state **A**), the Markov model is designed with no
incoming transitions to the state, so the cost of this operation can be
included into the model by linking it to the transitions from **A** to other
states. When applying costs to transitions, it's important to consider the
possible transitions and ensure that no "double-counting" occurs (for example
by transitions that return to the same state).
```{r}
#| echo = TRUE
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
```
# Point estimate
The initial deterministic analysis uses numerical best estimates for each
parameter (cost, utility, transition probability etc.). This results in a
model with an exact solution, but no estimate of the uncertainty, or of the
impact of individual parameters on the uncertainty of the outcome.
## Markov model
The Markov state transition model specifies which transitions are permitted
between states. Each of the 9 states is represented in `rdecision` by a
`MarkovState` object, with a description, a state cost and a state utility.
Each allowed transition is represented by a `Transition` object, defined by
the source state, the target state, and a transition cost. In all cases, cost
defaults to 0.0 if not set, and utility to 1.0.
In their model, Dong and Buxton incorporate the cost of the CAS technology as
an additional cost attached to interventions - in this case, this would affect
`cost_A`, `cost_E` and `cost_F` (but not `cost_G` which represents non-surgical
procedures). Utilities are unchanged and therefore the same Markov state
definitions can be used by both models.
```{r}
#| echo = TRUE
# 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
)
```
To fully define a semi Markov model, we assign the States, Transitions, duration
of a cycle (as a `difftime` object) and, optionally, annual `discount.cost`
and `discount.utility`.
```{r}
#| echo = TRUE
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
)
```
```{r}
#| echo = TRUE
with(data = SMM_base$tabulate_states(), expr = {
data.frame(
Name = Name,
Utility = Utility,
stringsAsFactors = FALSE
)
})
```
The `SemiMarkovModel` object can be represented graphically in the DOT format
using the `as_DOT()` method for drawing with the `dot` tool from Graphviz, or
with the `DiagrammeR` package. Alternatively it can be exported with the
`as_gml()` method for drawing with the `igraph` package or other package
which imports GML graph representations.
```{r}
#| fig.cap = "Markov state transition model for total knee replacement (TKR)",
#| echo = FALSE,
#| fig.height = 7,
#| purl = TRUE
local({
# create an igraph object
gml <- SMM_base$as_gml()
gmlfile <- tempfile(fileext = ".gml")
writeLines(gml, con = gmlfile)
ig <- igraph::read_graph(gmlfile, format = "gml")
# match layout to Dong & Buxton, fig 1
vxy <- matrix(
data = c(
+0.00, -0.75, +0.00, +0.75, -1.00, -0.25, +0.50, -0.75, +1.25,
+1.00, +0.05, +0.50, +0.00, -0.50, -0.50, -0.50, -1.00, -0.75
),
ncol = 2L,
dimnames = list(states, c("x", "y"))
)
layout <- matrix(
data = c(
vapply(X = igraph::V(ig), FUN.VALUE = 1.0, FUN = function(v) {
lbl <- igraph::vertex_attr(ig, "label", v)
return(vxy[[lbl, "x"]])
}),
vapply(X = igraph::V(ig), FUN.VALUE = 1.0, FUN = function(v) {
lbl <- igraph::vertex_attr(ig, "label", v)
return(vxy[[lbl, "y"]])
})
),
byrow = FALSE,
ncol = 2L
)
# vertex widths
vwidth <- vapply(X = igraph::V(ig), FUN.VALUE = 1.0, FUN = function(v) {
lbl <- igraph::vertex_attr(ig, "label", v)
return(nchar(lbl) * 4.0 + 2.0)
})
# edge curvatures
cm <- matrix(
data = 0.0, nrow = length(states), ncol = length(states),
dimnames = list(states, states)
)
cm[[
"TKR operation for knee problems", "TKR with serious complications"
]] <- -0.7
cm[[
"TKR operation for knee problems", "Normal health after primary TKR"
]] <- +1.4
cm[["Complex revision", "Death"]] <- -0.3
cm[["Simple revision", "Death"]] <- -0.2
cm[["Other treatments", "Death"]] <- -0.1
cm[["TKR with minor complications", "Death"]] <- +2.5
cm[["TKR with serious complications", "Death"]] <- +0.4
curves <- vapply(X = igraph::E(ig), FUN.VALUE = 1.0, FUN = function(e) {
# find source and target labels
trg <- igraph::head_of(ig, e)
trgl <- igraph::vertex_attr(ig, name = "label", index = trg)
src <- igraph::tail_of(ig, e)
srcl <- igraph::vertex_attr(ig, name = "label", index = src)
cr <- cm[[srcl, trgl]]
return(cr)
})
# plot the igraph
withr::with_par(
new = list(
oma = c(0L, 0L, 0L, 0L),
mar = c(1L, 6L, 1L, 6L),
xpd = NA
),
code = {
plot(
ig,
rescale = FALSE, asp = 0L,
vertex.shape = "rectangle", vertex.size = vwidth,
vertex.color = "white", vertex.label.color = "black",
edge.color = "black", edge.curved = curves,
edge.arrow.size = 0.75,
frame = FALSE,
layout = layout,
loop.size = 0.8
)
}
)
})
```
## Transition probabilities
For each allowed transition between a pair of states, probabilities are
reported in Table 2. These are added to the model using the
`set_probabilities` method, which takes for input a numerical matrix of size
N(states) x N(states).
For death outcomes, the probabilities are reported for *primary TKR*,
*TKR revision* and *all reasons*: for each state, the relevant probability of
death (transition to state **I**) will therefore be death(all reasons) +
death(primary TKR *or* TKR revision). Note that Table 5 reports the death
rates for TKR-specific death (primary or revision), which can be recovered
by introducing separate Markov states in the model to represent the 3 possible
causes of death, each with relevant transitions.
The total sum of all outgoing transitions for each state must be 1, so one
transition per state should be calculated as 1 - sum(all other transitions).
This is also verified by the `set_probabilities` method that will return an
error if the probability matrix does not have row-wise sums of 1. Alternatively,
exactly one cell in each row can be left undefined as `NA` and it will
automatically be populated to ensure a row-wise total of 1.
```{r}
#| echo = TRUE
# 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)
```
The CAS technology affects patient outcomes by reducing the probability of
serious complications by about 34%. This applies to all transition
probabilities into state **B**, which correspond to column 2 in the transition
matrix. Note that this also requires adjusting the row-wise sum to 1, by
increasing the probability of making a transition to another state. It is
assumed the states with increased transition probabilities are those defined as
`NA` in matrix `Pt`. Function `txeffect` adjusts the transition probability
matrix, given a relative risk.
```{r}
#| echo = TRUE
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)
}
```
The effect is applied to the transition matrix for conventional surgery to
get the transition matrix for computer-assisted surgery, as follows.
```{r}
#| echo = TRUE
# apply CAS_effect to the transition matrix
CAS_effect <- 0.34
Pt_CAS <- txeffect(Pt, CAS_effect)
SMM_CAS$set_probabilities(Pt_CAS)
```
## Running the model
The simulation described in the paper starts with 1000 patients and progresses
for 10 years, or 120 one-month cycles. At the start of the model, all patients
are in state **A**. The initial state of the simulation can be entered into
the model by applying the `reset` function with a named array describing the
number of simulated patients for each named state (in this case, 1000 in state
**A** and none in the other states).
The `cycles` function simulates the progress of the initial population through
the specified number of cycles, and, optionally, can apply half-cycle
correction to population and cost (`hcc.pop` and `hcc.cost`). In this example,
half-cycle correction is not used.
```{r}
#| echo = TRUE
# 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
)
```
Running the `cycles` function returns a data frame with the population in each
state, the progress of the simulation (in cycles and years), the total cost
(including both transition and occupancy costs), and the calculated QALY. Costs
and QALYs are discounted by the rates specified in the model definition and are
normalised by the initial population.
The table below displays the first few cycles of the model, showing the
transition of patients from state **A** at cycle 0 (starting state) through the
3 allowed transitions to states **B**, **C** and **D**, and subsequently to
further states depending on their assigned transition probabilities. Each cycle
accrues a cost (the cost of TKR surgery is attached here to cycle 1, as it was
assigned to the transitions out of **A**) and this, with the utilities
associated with each state, can be used to calculate the QALY.
```{r}
data.frame(
Cycle = SMM_base_10years[0L : 5L, "Cycle"],
Years = round(SMM_base_10years[0L : 5L, "Years"], digits = 2L),
Cost = gbp(SMM_base_10years[0L : 5L, "Cost"], p = TRUE, char = FALSE),
QALY = round(SMM_base_10years[0L : 5L, "QALY"], digits = 3L),
A = round(SMM_base_10years[0L : 5L, states[["A"]]], digits = 2L),
B = round(SMM_base_10years[0L : 5L, states[["B"]]], digits = 2L),
C = round(SMM_base_10years[0L : 5L, states[["C"]]], digits = 2L),
D = round(SMM_base_10years[0L : 5L, states[["D"]]], digits = 2L),
E = round(SMM_base_10years[0L : 5L, states[["E"]]], digits = 2L),
F = round(SMM_base_10years[0L : 5L, states[["F"]]], digits = 2L),
G = round(SMM_base_10years[0L : 5L, states[["G"]]], digits = 2L),
H = round(SMM_base_10years[0L : 5L, states[["H"]]], digits = 2L),
I = round(SMM_base_10years[0L : 5L, states[["I"]]], digits = 2L),
stringsAsFactors = FALSE
)
```
This progression can also be visualised using the base graphics function
`barplot`:
```{r}
#| fig.cap = "State occupancy progression for conventional surgery.",
#| purl = TRUE
plot_occupancy <- function(SMM) {
withr::with_par(
new = list(mar = c(5L, 4L, 1L, 14L) + 0.1),
code = {
states <- colnames(SMM)[
!colnames(SMM) %in% c("Cost", "QALY", "Cycle", "Years")
]
pal <- sample(hcl.colors(length(states), palette = "Spectral"))
barplot(
t(as.matrix(SMM[, states])),
xlab = "Cycle", ylab = "States", names.arg = SMM[, "Cycle"],
col = pal, border = NA, space = 0L
)
legend(
"topright", legend = states, fill = pal, bty = "n",
inset = c(-0.75, 0.0), xpd = TRUE
)
}
)
}
plot_occupancy(SMM_base_10years)
```
Dong and Buxton [-@dong2006] report results yearly (corresponding to cycles 12,
24, 36 etc.) with some cumulative values. The `cycles` function of
`SemiMarkovModel` returns a Markov trace (state occupancies at the end of each
cycle, with per-cycle costs and QALYs) in the form of a data frame. The Markov
trace can be converted to their Table 4 form via the function `as_table4`.
```{r}
#| echo = TRUE
# 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"],
round(m_cum[, "TKR with serious complications"] / 10L, digits = 2L),
round(m_cum[, "TKR with minor complications"] / 10L, digits = 2L),
round(m_cum[, "Complex revision"] / 10L, digits = 2L),
round(m_cum[, "Simple revision"] / 10L, digits = 2L),
round(m_trace[, "Death"] / 10L, digits = 2L),
round(m_cum[, "Cost"] * 1000L, digits = 2L),
round(m_cum[, "QALY"] * 1000L, digits = 2L)
),
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, ])
}
```
Replication of their Table 4 is given in the next two sections. Note that these
figures **do not match exactly** the published figures for costs and QALYs. This
might be a different application of the discount rate, or rounding errors
associated with the representation of currency in Excel. Small differences
in the CAS calculations can also probably be traced to rounding of the 34%
CAS effect figure.
### Conventional surgery
```{r}
#| echo = TRUE
t4_CS <- as_table4(SMM_base_10years)
```
```{r}
#| purl = FALSE
# test that Table 4 for conventional surgery is replicated"
stopifnot(
all.equal(nrow(t4_CS), 10L),
all.equal(
t4_CS[[10L, "Cumulative serious complication (%)"]], 87.32,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CS[[10L, "Cumulative minor complication (%)"]], 135.87,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CS[[10L, "Cumulative complex revision (%)"]], 5.13,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CS[[10L, "Cumulative simple revision (%)"]], 2.55,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CS[[10L, "Cumulative death (%)"]], 37.10,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CS[[10L, "Discounted costs (£)"]], 7791288.9,
tolerance = 100000.0, scale = 1.0
),
all.equal(
t4_CS[[10L, "Discounted QALYs"]], 5526.4,
tolerance = 250.0, scale = 1.0
)
)
```
```{r}
as.data.frame(t4_CS)
```
### Computer-assisted surgery
```{r}
#| echo = TRUE
t4_CAS <- as_table4(SMM_CAS_10years)
```
```{r}
#| purl = FALSE
# test that Table 4 for computer-assisted surgery is replicated
stopifnot(
all.equal(nrow(t4_CAS), 10L),
all.equal(
t4_CAS[[10L, "Cumulative serious complication (%)"]], 58.14,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CAS[[10L, "Cumulative minor complication (%)"]], 136.52,
tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CAS[[10L, "Cumulative complex revision (%)"]], 3.56, tolerance = 0.2,
scale = 1.0
),
all.equal(
t4_CAS[[10L, "Cumulative simple revision (%)"]], 1.89, tolerance = 0.2,
scale = 1.0
),
all.equal(
t4_CAS[[10L, "Cumulative death (%)"]], 37.06, tolerance = 0.2, scale = 1.0
),
all.equal(
t4_CAS[[10L, "Discounted costs (£)"]], 7208671.4, tolerance = 100000.0,
scale = 1.0
),
all.equal(
t4_CAS[[10L, "Discounted QALYs"]], 5541.2, tolerance = 250.0
)
)
```
```{r}
#| purl = FALSE
as.data.frame(t4_CAS)
```
```{r}
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
```
```{r}
#| purl = FALSE
# test that Deterministic CEA matches published value
local({
stopifnot(
all.equal(dcost, -583.0, tolerance = 50.0, scale = 1.0),
all.equal(dutil, 0.0164, tolerance = 0.005, scale = 1.0)
)
})
```
# Probabilistic Sensitivity Analysis
PSA adds a randomisation-based process to estimate not only the central value
of the estimate, but also its uncertainty. This is done by representing
individual parameters not as single numerical values but as distributions,
from which a range of possible values can be extracted at random. The choice
of distribution type and parametrisation depend on the type of variable and the
values reported in the literature or observed in a relevant patient cohort.
## Parameter distributions
### Utilities
These are constrained to a value between 0 and 1, and have been described with
a Beta distribution:
> A Beta function, with the mean equal to the point estimate and a high
variance to reflect the uncertainty, was used to generate a random utility
for each Markov state
The specific parameters for the distributions are not reported in the paper,
where instead a point estimate and a range (95% CI) are reported for each
utility distribution.
The Beta distribution can be parametrised using mean $\mu$ and sample size
$\nu$ (which relate to the more common parameters as $\alpha = \mu * \nu$
and $\beta = (1 - \mu) * \nu)$, using the point estimate as mean and assessing
numerically the sample size until the confidence intervals match the reported
range. Not all the ranges provided map to a valid Beta distribution with the
given mean, suggesting that the authors used a different method to derive their
Beta distributions. Note also that the authors report that
> To ensure a plausible relationship between the utilities of states “Normal
health after primary TKR,” “TKR with minor complications,” and “TKR with
serious complications,” the process was structured so that the hierarchical
relationship was retained but the differences between the utilities were
randomly drawn from the distribution.
but there is no description of the method used and/or how the resulting values
are constrained to the [0, 1] interval.
```{r}
#| echo = TRUE
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
```
These parameters can now be used to describe the distributions that each
utility value should be drawn from. This is represented in `rdecision` by a
`ModVar` variable, in this case specifically the `BetaModVar` child class. This
type of `ModVar` requires the `alpha` and `beta` parameters, which can be
calculated according to the equations above. Note that for state **I** (Death),
there is no need to define a `ModVar`, as the utility of this state is not
associated with any uncertainty.
```{r}
#| echo = TRUE
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
)
```
### Costs
These can take any non-negative value, and have been described with a Gamma
distribution:
> Variance for the Gamma distribution was estimated from the ranges. A Gamma
function was used to generate a random cost for each Markov state.
These distributions are also reported as means and 95% CI. In this case, the
means of the Gamma distributions are sufficiently far from 0 that the
approximation to a Gaussian distribution is acceptable to estimate the value
of the standard deviation.
The `GammaModVar` class requires the parameters `shape` and `scale`, which
relate to mean $\mu$ and variance $\sigma^2$ as $shape = \mu^2 / \sigma^2$
and $scale = \sigma^2 / \mu$.
```{r}
#| echo = TRUE
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
)
```
The cost of CAS is also modelled with a Gamma distribution, but in this case
the authors estimate the distribution parameters directly:
> We assumed that the ratio of its standard deviation over mean was four times
as large as that of conventional primary TKR
As a result, the total cost of a surgical procedure is the sum of the surgery
cost and the CAS cost, each of which is independently drawn from a defined
gamma distribution. This can be easily represented in `rdecision` using an
`ExprModVar`: a model variable that uses other model variables as operands.
The expression needs to be quoted using the `rlang::quo` command in order to
be stored rather than executed straight away.
```{r}
#| echo = TRUE
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)
)
```
The CAS effect in this model is a reduction in the transition probabilities
to state **B**, and was modelled using a lognormal distribution:
> Lognormal function was used to generate a random “effect of CAS.” The
variance of the “effect of CAS” was estimated from the clinical trials
Because a lognormally distributed value is always positive, the resulting
transition probabilities will be reduced by a non-zero fraction (note however
that values > 100% are possible, which would result an invalid negative
transition probability).
While it is possible to reconstruct the exact parameters of a lognormal
distribution given a mean and 95% CI, the article does not report a range for
this variable. Because exact replication is therefore impossible, an arbitrary
parametrisation is used here that is unlikely to yield an illegal value.
```{r}
#| echo = TRUE
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)
)
```
## Markov model
The same structure as the model used for the deterministic solution can be
used again, but in this case the costs and utilities will be represented by
distributions.
```{r}
#| echo = TRUE
# 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
)
```
The current model includes a number of model variables, that can be displayed
using the `modvar_table` method. This tabulation can also be used to compare
the 95% CI with the values reported in Table 3.
```{r}
#| echo = TRUE
with(data = SMM_base_PSA$modvar_table(), expr = {
data.frame(
Description = Description,
Distribution = Distribution,
Mean = round(Mean, digits = 5L),
SD = round(SD, digits = 5L),
Q2.5 = round(Q2.5, digits = 5L),
Q97.5 = round(Q97.5, digits = 5L),
stringsAsFactors = FALSE
)
})
```
Note that the list of `ModVar`s associated with the CAS model includes not only
the variables that were explicitly assigned to States and Transitions, such as
*Utility of state A*, but also the implicit `ModVar`s underlying `ExprModVar`s,
such as *Cost of CAS*, which is used to calculate *Cost of state A with CAS*.
```{r}
#| echo = TRUE
with(data = SMM_CAS_PSA$modvar_table(), expr = {
data.frame(
Description = Description,
Distribution = Distribution,
Mean = round(Mean, digits = 5L),
SD = round(SD, digits = 5L),
Q2.5 = round(Q2.5, digits = 5L),
Q97.5 = round(Q97.5, digits = 5L),
stringsAsFactors = FALSE
)
})
```
## Transition probabilities
> To generate a logical multi-Markov state probabilistic transition matrix from
the initial point estimates, the Dirichlet distribution was used (11). We
estimated a count for each Markov state by the transition probabilities (we
assumed that the total counts equalled 1,000). We used random number and
Gamma distribution formulae to generate a one-parameter (standard) Gamma
distribution for each cell of the transition matrix. The one-parameter Gamma
distribution in Excel was obtained by setting the first (alpha) parameter
equal to the estimated count and the second (beta) parameter equal to 1. The
final step was to “normalize” the realizations from the Gamma distribution
back to the 0–1 scale by dividing each realization through by the
corresponding row total.
The complex description above is an approximation of a Dirichlet distribution
when constrained by the functionalities of Excel. In `rdecision`, variables
drawn from a multinomial probabilistic distribution can be defined using
in-built functions:
* For each state, construct a `DirichletDistribution` from the known counts or
proportions for each exiting transition. If using proportions, a population
should also specified to correctly assign variance to the probabilistic
distribution. In this example, the authors assumed a population of 1000.
* For each state-specific distribution, draw a random realisation of the
multivariate distribution (`sample()` method) and populate the allowed
transitions with numerical values (`r()` method). This is performed by
the `dirichletify` function.
* Repeat the process at every cycle to sample different outgoing transition
probabilities from the distribution.
```{r}
#| echo = TRUE
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)
}
```
## PSA runs
The following code executes 250 runs of the model for conventional surgery and
for computer-assisted surgery, each of 120 cycles. For each run the set of
model variables are sampled from their uncertainty distributions. The Markov
traces for each run, in matrix form, are converted to the format used by Dong
et al in their Table 4, and added to 3D arrays which store the traces from all
the runs, separately for conventional and computer-assisted surgery.
The transition probability matrix is sampled from an assumed matrix of counts
using the `dirichletify` function.
In the case of CAS, there is an additional step in the simulation process due
to the fact that the transition probabilities into state **B** (column 2 in the
transition matrix) are modified by a multiplier, *Effect of CAS*, which is
itself drawn from a distribution. In this case, the relevant `ModVar` needs to be
randomised at every cycle, and applied to the transition matrix.
This a "paired" simulation; that is, we are interested in the
difference between computer-assisted surgery and conventional surgery,
taking account of the variation in inputs. In this model, many variables are
common to both models, and we make sure to randomise all variables once, before
running each model, every cycle. This is analogous, in statistics, to a paired
t-test rather than a two-sample t-test. When there are many common variables,
the results of the two models will be correlated.
```{r}
#| echo = TRUE
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
}
```
## Results
The means of the cohort simulation results should be very similar, but not
necessarily identical, to those in Table 4 (and replicated analytically above).
Each run of this script is expected to yield a slightly different result.
The addition of probabilistic elements can be used to describe not only the
means (which should approximate the deterministic cohort model results), but
also the distributions of the estimates, giving the uncertainties associated
with the two models.
In the following two sections, the mean of all runs are presented as per Table 4
of Dong & Buxton [-@dong2006] and the uncertainties are presented as per their
Table 5, with the addition of the upper and lower confidence limits of the mean.
### Conventional surgery
```{r}
#| echo = TRUE
local({
t4_CS_PSA_mean <- apply(t4_CS_PSA, MARGIN = c(1L, 2L), FUN = mean)
t4_CS_PSA_mean <- apply(
t4_CS_PSA_mean, MARGIN = c(1L, 2L), FUN = round, digits = 2L
)
as.data.frame(t4_CS_PSA_mean)
})
```
```{r}
#| echo = TRUE
local({
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"]] <- round(mean(t4_CS_PSA[10L, f, ]), digits = 2L)
t5_CS[[f, "SD"]] <- round(sd(t4_CS_PSA[10L, f, ]), digits = 2L)
t5_CS[[f, "Q2.5"]] <- round(
quantile(t4_CS_PSA[10L, f, ], probs = 0.025), digits = 2L
)
t5_CS[[f, "Q97.5"]] <- round(
quantile(t4_CS_PSA[10L, f, ], probs = 0.975), digits = 2L
)
}
as.data.frame(t5_CS)
})
```
### Computer-assisted surgery
```{r}
#| echo = TRUE
local({
t4_CAS_PSA_mean <- apply(t4_CAS_PSA, MARGIN = c(1L, 2L), FUN = mean)
t4_CAS_PSA_mean <- apply(
t4_CAS_PSA_mean, MARGIN = c(1L, 2L), FUN = round, digits = 2L
)
as.data.frame(t4_CAS_PSA_mean)
})
```
```{r}
#| echo = TRUE
local({
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"]] <- round(mean(t4_CAS_PSA[10L, f, ]), digits = 2L)
t5_CAS[[f, "SD"]] <- round(sd(t4_CAS_PSA[10L, f, ]), digits = 2L)
t5_CAS[[f, "Q2.5"]] <- round(
quantile(t4_CAS_PSA[10L, f, ], probs = 0.025), digits = 2L
)
t5_CAS[[f, "Q97.5"]] <- round(
quantile(t4_CAS_PSA[10L, f, ], probs = 0.975), digits = 2L
)
}
as.data.frame(t5_CAS)
})
```
Note that the means correlate well with those reported in Table 5, but the
assertion that "[CAS] had lower variances for each indicator" is not reliably
replicated for every variable. This is likely to be due to the fact that the
parameters for the distribution of the *Effect of CAS* could not be extracted,
and the values chosen here may assign excessive uncertainty to this variable.
### Cost effectiveness analysis
```{r}
#| echo = TRUE
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
```
```{r}
#| purl = FALSE
# test that PSA CEA matches deterministic calculation
local({
rc <- range(dcost_psa)
ru <- range(dutil_psa)
stopifnot(
dcost >= rc[[1L]],
dcost <= rc[[2L]],
dutil >= ru[[1L]],
dutil <= ru[[2L]]
)
})
```
From PSA, after 10 years the mean (95% CI) cost change per patient was
`r round(mean(dcost_psa), digits = 2L)`
(`r round(quantile(dcost_psa, probs = 0.025), digits = 2L)` to
`r round(quantile(dcost_psa, probs = 0.975), digits = 2L)`) £,
and the mean change in utility was
`r round(mean(dutil_psa), digits = 2L)`
(`r round(quantile(dutil_psa, probs = 0.025), digits = 2L)` to
`r round(quantile(dutil_psa, probs = 0.975), digits = 2L)`) QALY, compared
with the deterministic values of -£583 and +0.0148 respectively, as reported in
the paper. The incremental cost effectiveness ratio was less than
£30,000 per QALY for
`r round(100 * sum(icer_psa < 30000) / length(icer_psa), 2L)`% of runs. The
distribution of results on the cost effectiveness plane is shown in the
figure below (for comparison with figure 2 of the paper).
```{r}
#| fig.cap = "Cost-effectiveness plane. QALY = Quality-adjusted life year",
#| purl = TRUE
withr::with_par(
new = list(mar = c(1L, 2L, 2L, 1L)),
code = {
plot(
dcost_psa ~ dutil_psa,
pch = 20L,
xlim = c(-0.1, 0.15), ylim = c(-5000.0, 1000.0),
axes = FALSE,
xlab = "", ylab = ""
)
axis(side = 1L, pos = 0.0, las = 1L)
axis(side = 2L, pos = 0.0, las = 1L)
mtext(
text = expression(paste(Delta, "QALYs")),
side = 2L, adj = 5L / 6L
)
mtext(
text = expression(paste(Delta, "Cost (£)")),
side = 3L, adj = 2L / 5L
)
}
)
```