---
title: "Elementary Semi-Markov Model (Chancellor 1997)"
subtitle: "Monotherapy versus combination therapy for HIV"
author: "Andrew J. Sims"
date: "May 2021"
output: 
  rmarkdown::html_vignette:
    fig_width: 7
    fig_height: 5
    fig_caption: true
    df_print: kable
vignette: >
  %\VignetteIndexEntry{Elementary Semi-Markov Model (Chancellor 1997)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: "REFERENCES.bib"
csl: "nature-no-et-al.csl"
---

```{r}
#| purl = FALSE,
#| include = FALSE
knitr::opts_chunk$set(
  echo = FALSE,
  fig.keep = "last",
  fig.align = "center",
  collapse = TRUE,
  comment = "#>"
)
```

```{r}
#| purl = FALSE
#nolint start
```

```{r}
library(rdecision)
```

```{r}
#| purl = FALSE
#nolint end
```

# Introduction
This vignette is an example of an elementary semi-Markov model using 
the `rdecision` package. It is based on the example given by 
Briggs *et al* [-@briggs2006] (Exercise 2.5) which itself is based on a
model described by Chancellor *et al* [-@chancellor1997]. The model compares 
a combination therapy of Lamivudine/Zidovudine versus Zidovudine monotherapy 
in people with HIV infection.

# Creating the model

## Model structure
The model is constructed by forming a graph, with each state as a 
node and each transition as an edge. Nodes of class `MarkovState` and edges
of class `Transition` have various properties whose values reflect the
variables of the model (costs, rates etc.). Because the model is intended to 
evaluate survival, the utility of states A, B and C are set to 1 (by default) 
and state D to zero. Thus the incremental quality adjusted life years gained per
cycle is equivalent to the survival function. Because the structure of the
model is identical for monotherapy and combination therapy, we will use the 
same model throughout. For this reason, the costs of occupancy of each state
and the costs of making transitions between states are set to zero when the
model is created, and will be changed each time the model is run.

```{r}
#| echo = TRUE
# create Markov states
sA <- MarkovState$new("A")
sB <- MarkovState$new("B")
sC <- MarkovState$new("C")
sD <- MarkovState$new("D", utility = 0.0)
# create transitions
tAA <- Transition$new(sA, sA)
tAB <- Transition$new(sA, sB)
tAC <- Transition$new(sA, sC)
tAD <- Transition$new(sA, sD)
tBB <- Transition$new(sB, sB)
tBC <- Transition$new(sB, sC)
tBD <- Transition$new(sB, sD)
tCC <- Transition$new(sC, sC)
tCD <- Transition$new(sC, sD)
tDD <- Transition$new(sD, sD)
# set discount rates
cDR <- 6.0 # annual discount rate, costs (%)
oDR <- 0.0 # annual discount rate, benefits (%)
# construct the model
m <- SemiMarkovModel$new(
  V = list(sA, sB, sC, sD),
  E = list(tAA, tAB, tAC, tAD, tBB, tBC, tBD, tCC, tCD, tDD),
  discount.cost = cDR / 100.0,
  discount.utility = oDR / 100.0
)
```

## Costs and discounts
The costs and discount rates used in the model (1995 rates) are numerical
constants, and are defined as follows. 
```{r}
#| echo = TRUE
# drug costs
cAZT <- 2278.0 # zidovudine drug cost
cLam <- 2087.0 # lamivudine drug cost

# direct medical and community costs
dmca <- 1701.0 # direct medical costs associated with state A
dmcb <- 1774.0 # direct medical costs associated with state B
dmcc <- 6948.0 # direct medical costs associated with state C
ccca <- 1055.0 # Community care costs associated with state A
cccb <- 1278.0 # Community care costs associated with state B
cccc <- 2059.0 # Community care costs associated with state C

# occupancy costs with monotherapy
cAm <- dmca + ccca + cAZT
cBm <- dmcb + cccb + cAZT
cCm <- dmcc + cccc + cAZT

# occupancy costs with combination therapy
cAc <- dmca + ccca + cAZT + cLam
cBc <- dmcb + cccb + cAZT + cLam
cCc <- dmcc + cccc + cAZT + cLam
```

## Treatment effect
The treatment effect was estimated by Chancellor *et al* [-@chancellor1997]
via a meta-analysis, and is defined as follows:
```{r}
#| echo = TRUE
RR <- 0.509
```

## Transition rates and probabilities
Briggs *et al* [-@briggs2006] interpreted the observed transition counts 
in 1 year as transition probabilities by dividing counts by the total 
transitions observed from each state. With this assumption, the annual 
(per-cycle) transition probabilities are calculated as follows and applied
to the model via the `set_probabilities` function.
```{r}
#| echo = TRUE
# transition counts
nAA <- 1251L
nAB <- 350L
nAC <- 116L
nAD <- 17L
nBB <- 731L
nBC <- 512L
nBD <- 15L
nCC <- 1312L
nCD <- 437L
# create transition matrix
nA <- nAA + nAB + nAC + nAD
nB <- nBB + nBC + nBD
nC <- nCC + nCD
Ptm <- matrix(
  c(nAA / nA, nAB / nA, nAC / nA, nAD / nA,
    0.0, nBB / nB, nBC / nB, nBD / nB,
    0.0,      0.0, nCC / nC, nCD / nC,
    0.0,      0.0,      0.0,      1.0),
  nrow = 4L, byrow = TRUE,
  dimnames = list(
    source = c("A", "B", "C", "D"), target = c("A", "B", "C", "D")
  )
)
```

```{r}
#| echo = FALSE,
#| purl = FALSE
# test that monotherapy transition matrix agrees with Briggs Table 2.2
local({
  E <- matrix(
    c(0.721, 0.202, 0.067, 0.010,
      0.000, 0.581, 0.407, 0.012,
      0.000, 0.000, 0.750, 0.250,
      0.000, 0.000, 0.000, 1.000),   # typo in book (D,D) = 1!
    byrow = TRUE,
    nrow = 4L,
    dimnames = list(
      source = c("A", "B", "C", "D"), target = c("A", "B", "C", "D")
    )
  )
  stopifnot(all.equal(Ptm, E, tolerance = 0.01, scale = 1.0))
})
```

More usually, fully observed transition counts are converted into 
transition rates, rather than probabilities, as described by Welton and 
Ades [-@welton2005]. This is because counting events and measuring total time
at risk includes individuals who make more than one transition during the
observation time, and can lead to rates with values which exceed 1. In contrast,
the difference between a census of the number of individuals in each state at
the start of the interval and a census at the end is directly related to the
per-cycle probability. As Miller and Homan [-@miller1994], Welton and Ades
[-@welton2005], O'Mahony *et al* [-@omahony2015], Jones *et al* [-@jones2017]
and others note, conversion between rates and probabilities for multi-state
Markov models is non-trivial and care is needed when modellers calculate
probabilities from published rates for use in `SemiMarkoModel`s.

# Checking the model

## Diagram
A representation of the model in DOT format ([Graphviz](https://graphviz.org))
can be created using the `as_DOT` function of `SemiMarkovModel`. The function
returns a character vector which can be saved in a file (`.gv` extension) for
visualization with the `dot` tool of Graphviz, or plotted directly in R via
the `DiagrammeR` package. Alternatively, the graph can be saved in the
graph modelling language (GML) format, and imported into the `igraph` package
as a graph. This method offers more options for adjusting the appearance of
the model. The Markov model is shown in the figure below.

```{r}
#| purl = TRUE,
#| fig.cap = "Markov model for comparison of HIV therapy.
#|            A: 200 < cd4 < 500, B: cd4 < 200, C: AIDS, D: Death.",
local({
  # create an igraph object
  gml <- m$as_gml()
  gmlfile <- tempfile(fileext = ".gml")
  writeLines(gml, con = gmlfile)
  ig <- igraph::read_graph(gmlfile, format = "gml")
  # define vertex positions
  yv <- c(A = 1.0, B = 1.0 / 3.0, C = -1.0 / 3.0, D = -1.0)
  # set vertex positions
  layout <- matrix(
    data = c(
      0L, 0L, 0L, 0L,
      vapply(X = igraph::V(ig), FUN.VALUE = 1.0, FUN = function(v) {
        lbl <- igraph::vertex_attr(ig, "label", v)
        return(yv[[lbl]])
      })
    ),
    byrow = FALSE,
    ncol = 2L
  )
  # define edge curvatures
  cm <- matrix(
    data = 0.0, nrow = 4L, ncol = 4L,
    dimnames = list(LETTERS[seq(4L)], LETTERS[seq(4L)])
  )
  cm[["A", "D"]] <- 1.5
  cm[["A", "C"]] <- 1.0
  cm[["B", "D"]] <- -1.0
  # set edge curvatures
  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(1L, 1L, 1L, 1L),
      mar = c(1L, 1L, 1L, 1L),
      xpd = NA
    ),
    code = {
      plot(
        ig,
        rescale = FALSE, asp = 0L,
        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
      )
    }
  )
})
```

## Per-cycle transition probabilities
The per-cycle transition probabilities are the cells of the Markov
transition matrix. For the monotherapy model, the transition matrix is
shown below. This is consistent with the Table 1 of 
Chancellor *et al* [-@chancellor1997].
```{r}
#| purl = TRUE
with(data = as.data.frame(Ptm), expr = {
  data.frame(
    A = round(A, digits = 3L),
    B = round(B, digits = 3L),
    C = round(C, digits = 3L),
    D = round(D, digits = 3L),
    row.names = row.names(Ptm),
    stringsAsFactors = FALSE
  )
})
```

# Running the model
Model function `cycle` applies one cycle of a Markov model to a defined 
starting population in each state. It returns a table with one row per state, 
and each row containing several columns, including the population at the end of
the state and the cost of occupancy of states, normalized by the number of 
patients in the cohort, with discounting applied.

Multiple cycles are run by feeding the state populations at the end of
one cycle into the next. Function `cycles` does this and returns a data frame 
with one row per cycle, and each row containing the state populations and the
aggregated cost of occupancy for all states, with discounting applied. This is
done below for the first 20 cycles of the model for monotherapy, with discount.
For convenience, and future use with probabilistic sensitivity analysis, a
function, `run_mono` is used to wrap up the steps needed to run 20 cycles of 
the model for monotherapy. The arguments to the function are the transition
probability matrix, the occupancy costs for states A, B, and C, and logical
variables which determine whether to apply half-cycle correction to the state
populations, costs and QALYs returned in the Markov trace.
```{r}
#| echo = TRUE
# function to run model for 20 years of monotherapy
run_mono <- function(Ptm, cAm, cBm, cCm, hcc = FALSE) {
  # create starting populations
  N <- 1000L
  populations <- c(A = N, B = 0L, C = 0L, D = 0L)
  m$reset(populations)
  # set costs
  sA$set_cost(cAm)
  sB$set_cost(cBm)
  sC$set_cost(cCm)
  # set transition probabilities
  m$set_probabilities(Ptm)
  # run 20 cycles
  tr <- m$cycles(
    ncycles = 20L, hcc.pop = hcc, hcc.cost = FALSE, hcc.QALY = hcc
  )
  return(tr)
}
```

> Coding note: In function `run_mono`, the occupancy costs for states A, B
  and C are set via calls to function `set_cost()` which is associated with
  a `MarkovState` object. Although these are set *after* the state objects
  `sA`, `sB` and `sC` have been added to model `m`, the updated costs are used
  when the model is cycled. This is because R's R6 objects, such as
  Markov states and transitions, are passed by reference. That is, if an
  R6 object such as a `MarkovState` changes, any other object that refers to 
  it, such as a `SemiMarkovModel` will see the changes. This behaviour is
  different from regular R variable types, such as numeric variables, which
  are passed by value; that is, a copy of them is created within the function
  to which they are passed, and any change to the original would not apply
  to the copy.

The model is run by calling the new function, with appropriate arguments. The
cumulative cost and life years are calculated by summing the appropriate
columns from the Markov trace, as follows:
```{r}
#| echo = TRUE
MT.mono <- run_mono(Ptm, cAm, cBm, cCm)
el.mono <- sum(MT.mono$QALY)
cost.mono <- sum(MT.mono$Cost)
```

```{r}
#| echo = FALSE,
#| purl = FALSE
# test that monotherapy QALY and cost agrees with Briggs tables 2.3, 2.4
local({
  stopifnot(
    all.equal(el.mono, 7.996, tolerance = 0.005, scale = 1.0),
    all.equal(cost.mono, 44663.0, tolerance = 100.0, scale = 1.0)
  )
})
```

The populations and discounted costs are consistent with Briggs *et al*,
Table 2.3 [-@briggs2006], and the QALY column is consistent with Table 2.4
(without half cycle correction). No discount was applied to the utilities.
```{r}
#| purl = TRUE
with(data = MT.mono, expr = {
  data.frame(
    Years = Years,
    A = round(A, digits = 0L),
    B = round(B, digits = 0L),
    C = round(C, digits = 0L),
    D = round(D, digits = 0L),
    Cost = round(Cost, digits = 0L),
    QALY = round(QALY, digits = 3L),
    stringsAsFactors = FALSE
  )
})
```

# Model results

## Monotherapy
The estimated life years is approximated by summing the proportions of patients
left alive at each cycle (Briggs *et al* [@briggs2006], Exercise 2.5). This is 
an approximation because it ignores the population who remain alive after 
21 years, and assumes all deaths occurred at the start of each cycle. For
monotherapy the expected life gained is `r round(el.mono, 3L)` years at a cost
of `r gbp(cost.mono)` GBP.

## Combination therapy
For combination therapy, a similar model was created, with adjusted costs and 
transition rates. Following Briggs *et al* [@briggs2006] the treatment effect
was applied to the probabilities, and these were used as inputs to the model.
More usually, treatment effects are applied to rates, rather than probabilities.
```{r}
#| echo = TRUE
# annual probabilities modified by treatment effect
pAB <- RR * nAB / nA
pAC <- RR * nAC / nC
pAD <- RR * nAD / nA
pBC <- RR * nBC / nB
pBD <- RR * nBD / nB
pCD <- RR * nCD / nC
# annual transition probability matrix
Ptc <- matrix(
  c(1.0 - pAB - pAC - pAD,               pAB,         pAC, pAD,
    0.0, (1.0 - pBC - pBD),         pBC, pBD,
    0.0,               0.0, (1.0 - pCD), pCD,
    0.0,               0.0,         0.0, 1.0),
  nrow = 4L, byrow = TRUE,
  dimnames = list(
    source = c("A", "B", "C", "D"), target = c("A", "B", "C", "D")
  )
)
```

```{r}
#| echo = FALSE,
#| purl = FALSE
# test that combo therapy transition matrix agrees with Briggs Table 2.2
local({
  E <- matrix(
    c(0.858, 0.103, 0.034, 0.005,
      0.000, 0.787, 0.207, 0.006,
      0.000, 0.000, 0.873, 0.127,
      0.000, 0.000, 0.000, 1.000),
    byrow = TRUE,
    nrow = 4L,
    dimnames = list(
      source = c("A", "B", "C", "D"), target = c("A", "B", "C", "D")
    )
  )
  stopifnot(all.equal(Ptc, E, tolerance = 0.01, scale = 1.0))
})
```

The resulting per-cycle transition matrix for the combination therapy is as
follows:
```{r}
#| purl = TRUE
with(data = as.data.frame(Ptc), expr = {
  data.frame(
    A = round(A, digits = 3L),
    B = round(B, digits = 3L),
    C = round(C, digits = 3L),
    D = round(D, digits = 3L),
    row.names = row.names(Ptc),
    stringsAsFactors = FALSE
  )
})
```

In this model, lamivudine is given for the first 2 years, with 
the treatment effect assumed to persist for the same period. The
state populations and cycle numbers are retained by the model between 
calls to `cycle` or `cycles` and can be retrieved by calling `get_populations`.
In this example, the combination therapy model is run for 2 cycles, then the
population is used to continue with the monotherapy model for the remaining
18 years. The `reset` function is used to set the cycle number and elapsed
time of the new run of the mono model. As before, function `run_comb` is created
to wrap up these steps, so they can be used repeatedly for different values of
the model variables.
```{r}
#| echo = TRUE
# function to run model for 2 years of combination therapy and 18 of monotherapy
run_comb <- function(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc, hcc = FALSE) {
  # set populations
  N <- 1000L
  populations <- c(A = N, B = 0L, C = 0L, D = 0L)
  m$reset(populations)
  # set the transition probabilities accounting for treatment effect
  m$set_probabilities(Ptc)
  # set the costs including those for the additional drug
  sA$set_cost(cAc)
  sB$set_cost(cBc)
  sC$set_cost(cCc)
  # run first 2 yearly cycles with additional drug costs and tx effect
  tr <- m$cycles(2L, hcc.pop = hcc, hcc.cost = FALSE, hcc.QALY = hcc)
  # save the state populations after 2 years
  populations <- m$get_populations()
  # revert probabilities to those without treatment effect
  m$set_probabilities(Ptm)
  # revert costs to those without the extra drug
  sA$set_cost(cAm)
  sB$set_cost(cBm)
  sC$set_cost(cCm)
  # restart the model with populations from first 2 years with extra drug
  m$reset(
    populations,
    icycle = 2L,
    elapsed = as.difftime(365.25 * 2.0, units = "days")
  )
  # run for next 18 years, combining the traces
  tr <- rbind(
    tr,
    m$cycles(ncycles = 18L, hcc.pop = hcc, hcc.cost = FALSE, hcc.QALY = hcc)
  )
  # return the trace
  return(tr)
}
```

The model is run by calling the new function, with appropriate arguments,
as follows. The incremental cost effectiveness ratio (ICER) is also calculated,
as the ratio of the incremental cost to the incremental life years of the 
combination therapy compared with monotherapy.
```{r}
#| echo = TRUE
MT.comb <- run_comb(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc)
el.comb <- sum(MT.comb$QALY)
cost.comb <- sum(MT.comb$Cost)
icer <- (cost.comb - cost.mono) / (el.comb - el.mono)
```

```{r}
#| echo = FALSE,
#| purl = FALSE
# test that combo therapy QALY, cost and ICER agree with Briggs ex 2.5
local({
  stopifnot(
    all.equal(el.comb, 8.937, tolerance = 0.02, scale = 1.0),
    all.equal(cost.comb, 50602.0, tolerance = 100.0, scale = 1.0),
    all.equal(icer, 6276.0, tolerance = 20.0, scale = 1.0)
  )
})
```

The Markov trace for combination therapy is as follows:
```{r}
#| purl = TRUE
with(data = MT.comb, expr = {
  data.frame(
    Years = Years,
    A = round(A, digits = 0L),
    B = round(B, digits = 0L),
    C = round(C, digits = 0L),
    D = round(D, digits = 0L),
    Cost = round(Cost, digits = 0L),
    QALY = round(QALY, digits = 3L),
    stringsAsFactors = FALSE
  )
})
```

## Comparison of treatments
Over the 20 year time horizon, the expected life
years gained for monotherapy was `r round(el.mono, 3L)` years at a total cost
per patient of 
`r gbp(cost.mono)` GBP. The expected life years gained with combination therapy
for two years was `r round(el.comb, 3L)` at a total cost per patient of 
`r gbp(cost.comb)` GBP. The incremental change in life years was 
`r round(el.comb - el.mono, 3L)` years at an incremental cost of 
`r gbp(cost.comb - cost.mono)` GBP, giving an ICER of `r gbp(icer)` GBP/QALY.
This is consistent with the result obtained by Briggs *et al* [-@briggs2006]
(6276 GBP/QALY), within rounding error.

## Results with half-cycle correction
With half-cycle correction applied to the state populations, the model can
be recalculated as follows.
```{r}
#| echo = TRUE
MT.mono.hcc <- run_mono(Ptm, cAm, cBm, cCm, hcc = TRUE)
el.mono.hcc <- sum(MT.mono.hcc$QALY)
cost.mono.hcc <- sum(MT.mono.hcc$Cost)
MT.comb.hcc <- run_comb(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc, hcc = TRUE)
el.comb.hcc <- sum(MT.comb.hcc$QALY)
cost.comb.hcc <- sum(MT.comb.hcc$Cost)
icer.hcc <- (cost.comb.hcc - cost.mono.hcc) / (el.comb.hcc - el.mono.hcc)
```

```{r}
#| echo = FALSE,
#| purl = FALSE
# test that model with HCC agrees with Briggs ex 2.5
local({
  stopifnot(
    all.equal(el.mono.hcc, 8.475, tolerance = 0.03, scale = 1.0),
    all.equal(cost.mono.hcc, 44663.0, tolerance = 100.0, scale = 1.0),
    all.equal(el.comb.hcc, 9.42, tolerance = 0.02, scale = 1.0),
    all.equal(cost.comb.hcc, 50602.0, tolerance = 100.0, scale = 1.0),
    all.equal(
      icer.hcc, (50602.0 - 44663.0) / (9.42 - 8.475),
      tolerance = 20.0, scale = 1.0
    )
  )
})
```

Over the 20 year time horizon, the expected life
years gained for monotherapy was `r round(el.mono.hcc, 3L)` years at a total
cost per patient of `r gbp(cost.mono.hcc)` GBP. The expected life years gained
with combination therapy for two years was `r round(el.comb.hcc, 3L)` at a
total cost per patient of `r gbp(cost.comb.hcc)` GBP. The incremental change in
life years was `r round(el.comb.hcc - el.mono.hcc, 3L)` years at an incremental
cost of `r gbp(cost.comb.hcc - cost.mono.hcc)` GBP, giving an ICER of
`r gbp(icer.hcc)` GBP/QALY.

# Probabilistic sensitivity analysis
In their Exercise 4.7, Briggs *et al* [-@briggs2006] extended the original model
to account for uncertainty in the estimates of the values of the model
variables. In this section, the exercise is replicated in `rdecision`, using
the same assumptions.

## Costs
Although it is possible to sample from uncertainty distributions using the
functions in R standard package `stats` (e.g., `rbeta`), `rdecision` introduces
the notion of a `ModVar`, which is an object that can represent a model variable
with an uncertainty distribution. Many of the class methods in `redecision` will
accept a `ModVar` as alternative to a numerical value as an argument, and will
automatically sample from its uncertainty distribution. 

The model costs are represented as `ModVar`s of various types, as follows. The
state occupancy costs for both models involve a summation of other 
variables. Package `rdecision` introduces a form of `ModVar` that is defined
as a mathematical expression (an `ExprModVar`) potentially involving `ModVar`s.
The uncertainty distribution of `cAm`, for example, is complex, because it is a 
sum of two Gamma-distributed variables and a scalar, but `rdecision` takes care
of this when `cAm` is sampled.
```{r}
#| echo = TRUE
# direct medical and community costs (modelled as gamma distributions)
dmca <- GammaModVar$new("dmca", "GBP", shape = 1.0, scale = 1701.0)
dmcb <- GammaModVar$new("dmcb", "GBP", shape = 1.0, scale = 1774.0)
dmcc <- GammaModVar$new("dmcc", "GBP", shape = 1.0, scale = 6948.0)
ccca <- GammaModVar$new("ccca", "GBP", shape = 1.0, scale = 1055.0)
cccb <- GammaModVar$new("cccb", "GBP", shape = 1.0, scale = 1278.0)
cccc <- GammaModVar$new("cccc", "GBP", shape = 1.0, scale = 2059.0)

# occupancy costs with monotherapy
cAm <- ExprModVar$new("cA", "GBP", rlang::quo(dmca + ccca + cAZT))
cBm <- ExprModVar$new("cB", "GBP", rlang::quo(dmcb + cccb + cAZT))
cCm <- ExprModVar$new("cC", "GBP", rlang::quo(dmcc + cccc + cAZT))

# occupancy costs with combination therapy
cAc <- ExprModVar$new("cAc", "GBP", rlang::quo(dmca + ccca + cAZT + cLam))
cBc <- ExprModVar$new("cBc", "GBP", rlang::quo(dmcb + cccb + cAZT + cLam))
cCc <- ExprModVar$new("cCc", "GBP", rlang::quo(dmcc + cccc + cAZT + cLam))
```

## Treatment effect
The treatment effect is also represented by a `ModVar` whose uncertainty follows
a log normal distribution.
```{r}
#| echo = TRUE
RR <- LogNormModVar$new(
  "Tx effect", "RR", p1 = 0.509, p2 = (0.710 - 0.365) / (2.0 * 1.96), "LN7"
)
```

## Transition matrix
The following function generates a transition probability matrix from observed
counts, using Dirichlet distributions, as described by Briggs *et al*. This
could be achieved using the R `stats` function `rgamma`, but `rdecision` offers
the `DirichletDistribition` class for convenience, which is used here.
```{r}
#| echo = TRUE
# function to generate a probabilistic transition matrix
pt_prob <- function() {
  # create Dirichlet distributions for conditional probabilities
  DA <- DirichletDistribution$new(c(1251L, 350L, 116L, 17L)) # from A # nolint
  DB <- DirichletDistribution$new(c(731L, 512L, 15L))  # from B # nolint
  DC <- DirichletDistribution$new(c(1312L, 437L)) # from C # nolint
  # sample from the Dirichlet distributions
  DA$sample()
  DB$sample()
  DC$sample()
  # create the transition matrix
  Pt <- matrix(
    c(DA$r(), c(0.0, DB$r()), c(0.0, 0.0, DC$r()), c(0.0, 0.0, 0.0, 1.0)),
    byrow = TRUE,
    nrow = 4L,
    dimnames = list(
      source = c("A", "B", "C", "D"), target = c("A", "B", "C", "D")
    )
  )
  return(Pt)
}
```

## Running the PSA
The following code runs 1000 iterations of the model. At each run, the model
variables are sampled from their uncertainty distributions, the transition
matrix is sampled from count data, and the treatment effect is applied. 
Functions `run_mono` and `run_comb` are used to generate Markov traces for
each form of therapy, and the incremental costs, life years and ICER for
each run are saved in a matrix.
```{r}
#| echo = TRUE
# create matrix to hold the incremental costs and life years for each run
psa <- matrix(
  data = NA_real_, nrow = 1000L, ncol = 5L,
  dimnames = list(
    NULL, c("el.mono", "cost.mono", "el.comb", "cost.comb", "icer")
  )
)

# run the model repeatedly
for (irun in seq_len(nrow(psa))) {

  # sample variables from their uncertainty distributions
  cAm$set("random")
  cBm$set("random")
  cCm$set("random")
  cAc$set("random")
  cBc$set("random")
  cCc$set("random")
  RR$set("random")

  # sample the probability transition matrix from observed counts
  Ptm <- pt_prob()

  # run monotherapy model
  MT.mono <- run_mono(Ptm, cAm, cBm, cCm, hcc = TRUE)
  el.mono <- sum(MT.mono$QALY)
  cost.mono <- sum(MT.mono$Cost)
  psa[[irun, "el.mono"]] <- el.mono
  psa[[irun, "cost.mono"]] <- cost.mono

  # create Pt for combination therapy (Briggs applied the RR to the transition
  # probabilities - not recommended, but done here for reproducibility).
  Ptc <- Ptm
  for (i in 1L:4L) {
    for (j in 1L:4L) {
      Ptc[[i, j]] <- ifelse(i == j, NA, RR$get() * Ptc[[i, j]])
    }
    Ptc[i, which(is.na(Ptc[i, ]))] <- 1.0 - sum(Ptc[i, ], na.rm = TRUE)
  }

  # run combination therapy model
  MT.comb <- run_comb(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc, hcc = TRUE)
  el.comb <- sum(MT.comb$QALY)
  cost.comb <- sum(MT.comb$Cost)
  psa[[irun, "el.comb"]] <- el.comb
  psa[[irun, "cost.comb"]] <- cost.comb

  # calculate the icer
  psa[[irun, "icer"]] <- (cost.comb - cost.mono) / (el.comb - el.mono)
}
```

> Coding note: The state occupancy costs `cAm`, `cBm` etc. are now 
  `ModVar`s, rather than numeric variables as they were in the deterministic
  model. However, they can still be passed as arguments to
  `MarkovState$set_cost()`, via the arguments to helper functions `run_mono`
  and `run_comb`, and `rdecision` will manage them appropriately, without
  changing any other code. Documentation for functions in `rdecision` explains
  where this is supported by the package.

## Results
The mean (95% confidence interval) for the cost of monotherapy was 
`r gbp(mean(psa[, "cost.mono"]))` 
(`r gbp(quantile(psa[, "cost.mono"], probs = 0.025))` to 
`r gbp(quantile(psa[, "cost.mono"], probs = 0.975))`) GBP,
and the mean (95% CI) cost for combination therapy was
`r gbp(mean(psa[, "cost.comb"]))` 
(`r gbp(quantile(psa[, "cost.comb"], probs = 0.025))` to 
`r gbp(quantile(psa[, "cost.comb"], probs = 0.975))`) GBP. The life years
gained for monotherapy was
`r round(mean(psa[, "el.mono"]), 3L)` 
(`r round(quantile(psa[, "el.mono"], probs = 0.025), 3L)` to 
`r round(quantile(psa[, "el.mono"], probs = 0.975), 3L)`), and the life
years gained for combination therapy was 
`r round(mean(psa[, "el.comb"]), 3L)` 
(`r round(quantile(psa[, "el.comb"], probs = 0.025), 3L)` to 
`r round(quantile(psa[, "el.comb"], probs = 0.975), 3L)`). The mean ICER was
`r gbp(mean(psa[, "icer"]))` GBP/QALY with 95% confidence interval
`r gbp(quantile(psa[, "icer"], probs = 0.025))` to 
`r gbp(quantile(psa[, "icer"], probs = 0.975))` GBP/QALY.

```{r}
#| echo = FALSE,
#| purl = FALSE
# retrieve data set with individual run results from Briggs
data(BriggsEx47, package = "rdecision")
```

From 1000 simulations using an Excel version of the model by Briggs *et al*,
the corresponding values were as follows.
The mean (95% confidence interval) for the cost of monotherapy was 
`r gbp(mean(BriggsEx47[, "Mono.Cost"]))` 
(`r gbp(quantile(BriggsEx47[, "Mono.Cost"], probs = 0.025))` to 
`r gbp(quantile(BriggsEx47[, "Mono.Cost"], probs = 0.975))`) GBP,
and the mean (95% CI) cost for combination therapy was
`r gbp(mean(BriggsEx47[, "Comb.Cost"]))` 
(`r gbp(quantile(BriggsEx47[, "Comb.Cost"], probs = 0.025))` to 
`r gbp(quantile(BriggsEx47[, "Comb.Cost"], probs = 0.975))`) GBP. The life years
gained for monotherapy was
`r round(mean(BriggsEx47[, "Mono.LYs"]), 3L)` 
(`r round(quantile(BriggsEx47[, "Mono.LYs"], probs = 0.025), 3L)` to 
`r round(quantile(BriggsEx47[, "Mono.LYs"], probs = 0.975), 3L)`), and the life
years gained for combination therapy was 
`r round(mean(BriggsEx47[, "Comb.LYs"]), 3L)` 
(`r round(quantile(BriggsEx47[, "Comb.LYs"], probs = 0.025), 3L)` to 
`r round(quantile(BriggsEx47[, "Comb.LYs"], probs = 0.975), 3L)`). The mean ICER
was `r gbp(mean(BriggsEx47[, "ICER"]))` GBP/QALY with 95% confidence interval
`r gbp(quantile(BriggsEx47[, "ICER"], probs = 0.025))` to 
`r gbp(quantile(BriggsEx47[, "ICER"], probs = 0.975))` GBP/QALY.
  
# References