# Elementary Semi-Markov Model (Chancellor 1997)

### Monotherapy versus combination therapy for HIV

#### Andrew J. Sims

#### May 2021

# 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 al1 (Exercise 2.5) which itself is based on a model described by Chancellor et al2. 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.

    # 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.

    # 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 al2 via a meta-analysis, and is defined as follows:

    RR <- 0.509

## Transition rates and probabilities

Briggs et al1 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.

    # 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")
      )
    )

More usually, fully observed transition counts are converted into transition rates, rather than probabilities, as described by Welton and Ades3. 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 Homan4, Welton and Ades3, O’Mahony et al5, Jones et al6 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) 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.

Markov model for comparison of HIV therapy. A: 200 < cd4 < 500, B: cd4 < 200, C: AIDS, D: Death.

## 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 al2.

| A | B | C | D
--- | --- | --- | --- | ---
A | 0.721 | 0.202 | 0.067 | 0.010
B | 0.000 | 0.581 | 0.407 | 0.012
C | 0.000 | 0.000 | 0.750 | 0.250
D | 0.000 | 0.000 | 0.000 | 1.000

# 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.

    # 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:

    MT.mono <- run_mono(Ptm, cAm, cBm, cCm)
    el.mono <- sum(MT.mono$QALY)
    cost.mono <- sum(MT.mono$Cost)

The populations and discounted costs are consistent with Briggs et al, Table 2.31, and the QALY column is consistent with Table 2.4 (without half cycle correction). No discount was applied to the utilities.

Years | A | B | C | D | Cost | QALY
--- | --- | --- | --- | --- | --- | ---
0 | 1000 | 0 | 0 | 0 | 0 | 0.000
1 | 721 | 202 | 67 | 10 | 5153 | 0.990
2 | 520 | 263 | 181 | 36 | 5393 | 0.964
3 | 376 | 258 | 277 | 89 | 5368 | 0.911
4 | 271 | 226 | 338 | 165 | 5055 | 0.835
5 | 195 | 186 | 364 | 255 | 4541 | 0.745
6 | 141 | 147 | 361 | 350 | 3929 | 0.650
7 | 102 | 114 | 341 | 444 | 3301 | 0.556
8 | 73 | 87 | 309 | 531 | 2708 | 0.469
9 | 53 | 65 | 272 | 610 | 2179 | 0.390
10 | 38 | 49 | 234 | 679 | 1727 | 0.321
11 | 28 | 36 | 198 | 739 | 1350 | 0.261
12 | 20 | 26 | 165 | 789 | 1045 | 0.211
13 | 14 | 19 | 136 | 830 | 801 | 0.170
14 | 10 | 14 | 111 | 865 | 609 | 0.135
15 | 7 | 10 | 90 | 893 | 460 | 0.107
16 | 5 | 8 | 72 | 915 | 346 | 0.085
17 | 4 | 5 | 57 | 933 | 258 | 0.067
18 | 3 | 4 | 45 | 948 | 192 | 0.052
19 | 2 | 3 | 36 | 959 | 142 | 0.041
20 | 1 | 2 | 28 | 968 | 105 | 0.032

# Model results

## Monotherapy

The estimated life years is approximated by summing the proportions of patients left alive at each cycle (Briggs et al1, 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 7.991 years at a cost of 44,663 GBP.

## Combination therapy

For combination therapy, a similar model was created, with adjusted costs and transition rates. Following Briggs et al1 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.

    # 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")
      )
    )

The resulting per-cycle transition matrix for the combination therapy is as follows:

| A | B | C | D
--- | --- | --- | --- | ---
A | 0.859 | 0.103 | 0.034 | 0.005
B | 0.000 | 0.787 | 0.207 | 0.006
C | 0.000 | 0.000 | 0.873 | 0.127
D | 0.000 | 0.000 | 0.000 | 1.000

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.

    # 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.

    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)

The Markov trace for combination therapy is as follows:

Years | A | B | C | D | Cost | QALY
--- | --- | --- | --- | --- | --- | ---
0 | 1000 | 0 | 0 | 0 | 0 | 0.000
1 | 859 | 103 | 34 | 5 | 6912 | 0.995
2 | 737 | 169 | 80 | 14 | 6736 | 0.986
3 | 532 | 247 | 178 | 43 | 5039 | 0.957
4 | 384 | 251 | 270 | 96 | 4998 | 0.904
5 | 277 | 223 | 330 | 170 | 4713 | 0.830
6 | 200 | 186 | 357 | 258 | 4245 | 0.742
7 | 144 | 148 | 357 | 351 | 3684 | 0.649
8 | 104 | 115 | 337 | 443 | 3102 | 0.557
9 | 75 | 88 | 307 | 530 | 2551 | 0.470
10 | 54 | 66 | 271 | 609 | 2057 | 0.391
11 | 39 | 49 | 234 | 678 | 1633 | 0.322
12 | 28 | 37 | 198 | 737 | 1279 | 0.263
13 | 20 | 27 | 165 | 787 | 990 | 0.213
14 | 15 | 20 | 136 | 829 | 760 | 0.171
15 | 11 | 14 | 111 | 864 | 579 | 0.136
16 | 8 | 11 | 90 | 892 | 437 | 0.108
17 | 6 | 8 | 72 | 914 | 329 | 0.086
18 | 4 | 6 | 58 | 933 | 246 | 0.067
19 | 3 | 4 | 46 | 947 | 183 | 0.053
20 | 2 | 3 | 36 | 959 | 136 | 0.041

## Comparison of treatments

Over the 20 year time horizon, the expected life years gained for monotherapy was 7.991 years at a total cost per patient of 44,663 GBP. The expected life years gained with combination therapy for two years was 8.94 at a total cost per patient of 50,607 GBP. The incremental change in life years was 0.949 years at an incremental cost of 5,944 GBP, giving an ICER of 6,264 GBP/QALY. This is consistent with the result obtained by Briggs et al1 (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.

    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)

Over the 20 year time horizon, the expected life years gained for monotherapy was 8.475 years at a total cost per patient of 44,663 GBP. The expected life years gained with combination therapy for two years was 9.419 at a total cost per patient of 50,607 GBP. The incremental change in life years was 0.944 years at an incremental cost of 5,944 GBP, giving an ICER of 6,295 GBP/QALY.

# Probabilistic sensitivity analysis

In their Exercise 4.7, Briggs et al1 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.

    # 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.

    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.

    # 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.

    # 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 44,943 (22,084 to 90,008) GBP, and the mean (95% CI) cost for combination therapy was 50,792 (27,990 to 93,175) GBP. The life years gained for monotherapy was 8.495 (8.115 to 8.913), and the life years gained for combination therapy was 9.438 (8.91 to 9.913). The mean ICER was 6,396 GBP/QALY with 95% confidence interval 3,091 to 11,325 GBP/QALY.

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 45,924 (22,627 to 98,857) GBP, and the mean (95% CI) cost for combination therapy was 51,759 (28,714 to 103,447) GBP. The life years gained for monotherapy was 8.488 (8.08 to 8.89), and the life years gained for combination therapy was 9.413 (8.87 to 9.93). The mean ICER was 6,516 GBP/QALY with 95% confidence interval 2,692 to 11,319 GBP/QALY.

# References

1.
Briggs, A., Claxton, K. & Sculpher, M. Decision modelling for health economic evaluation. (Oxford University Press, 2006).
2.
Chancellor, J. V., Hill, A. M., Sabin, C. A., Simpson, K. N. & Youle, M. Modelling the cost effectiveness of Lamivudine/Zidovudine combination therapy in HIV infection. Pharmacoeconomics 12, 54–66 (1997).
3.
Welton, N. J. & Ades, A. E. Estimation of Markov Chain Transition Probabilities and Rates from Fully and Partially Observed Data: Uncertainty Propagation, Evidence Synthesis, and Model Calibration. Medical Decision Making 25, 633–645 (2005).
4.
Miller, D. K. & Homan, S. M. Determining Transition Probabilities: Confusion and Suggestions. Medical Decision Making 14, 52–58 (1994).
5.
O’Mahony, J. F., Newall, A. T. & Van Rosmalen, J. Dealing with Time in Health Economic Evaluation: Methodological Issues and Recommendations for Practice. PharmacoEconomics 33, 1255–1268 (2015).
6.
Jones, E., Epstein, D. & García-Mochón, L. A Procedure for Deriving Formulas to Convert Transition Rates to Probabilities for Multistate Markov Models. Medical Decision Making 37, 779–789 (2017).