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*^{1}
(Exercise 2.5) which itself is based on a model described by Chancellor
*et al*^{2}. The model
compares a combination therapy of Lamivudine/Zidovudine versus
Zidovudine monotherapy in people with HIV infection.

The variables used in the model are all numerical constants, and are defined as follows.

```
# transition counts
<- 1251L
nAA <- 350L
nAB <- 116L
nAC <- 17L
nAD <- 731L
nBB <- 512L
nBC <- 15L
nBD <- 1312L
nCC <- 437L
nCD # Healthcare system costs
<- 1701.0 # direct medical costs associated with state A
dmca <- 1774.0 # direct medical costs associated with state B
dmcb <- 6948.0 # direct medical costs associated with state C
dmcc <- 1055.0 # Community care costs associated with state A
ccca <- 1278.0 # Community care costs associated with state B
cccb <- 2059.0 # Community care costs associated with state C
cccc # Drug costs
<- 2278.0 # zidovudine drug cost
cAZT <- 2087.0 # lamivudine drug cost
cLam # Treatment effect
<- 0.509
RR # Discount rates
<- 6.0 # annual discount rate, costs (%)
cDR <- 0.0 # annual discount rate, benefits (%) oDR
```

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 (class `Transition`

) may
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.

```
# create Markov states for monotherapy (zidovudine only)
<- MarkovState$new("A", cost = dmca + ccca + cAZT)
sAm <- MarkovState$new("B", cost = dmcb + cccb + cAZT)
sBm <- MarkovState$new("C", cost = dmcc + cccc + cAZT)
sCm <- MarkovState$new("D", cost = 0.0, utility = 0.0)
sDm # create transitions
<- Transition$new(sAm, sAm)
tAAm <- Transition$new(sAm, sBm)
tABm <- Transition$new(sAm, sCm)
tACm <- Transition$new(sAm, sDm)
tADm <- Transition$new(sBm, sBm)
tBBm <- Transition$new(sBm, sCm)
tBCm <- Transition$new(sBm, sDm)
tBDm <- Transition$new(sCm, sCm)
tCCm <- Transition$new(sCm, sDm)
tCDm <- Transition$new(sDm, sDm)
tDDm # construct the model
<- SemiMarkovModel$new(
m.mono V = list(sAm, sBm, sCm, sDm),
E = list(tAAm, tABm, tACm, tADm, tBBm, tBCm, tBDm, tCCm, tCDm, tDDm),
discount.cost = cDR / 100.0,
discount.utility = oDR / 100.0
)
```

Briggs *et al*^{1}
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.

```
<- nAA + nAB + nAC + nAD
nA <- nBB + nBC + nBD
nB <- nCC + nCD
nC <- matrix(
Pt 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")
)
)$set_probabilities(Pt) m.mono
```

More usually, fully observed transition counts are converted into
transition rates (rather than probabilities), as described by Welton and
Ades^{3}. 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 another at the end is directly related to the per-cycle
probability. As Miller and Homan^{4}, Welton and Ades^{3}, Jones *et al*^{5} and others note, conversion between
rates and probabilities for multi-state Markov models is
non-trivial^{5} and care is needed
when modellers calculate probabilities from published rates for use in
`SemiMarkoModel`

s.

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. The Markov model for monotherapy is
shown in Figure 1.

The states in the model can be tabulated with the function
`tabulate_states`

. For the monotherapy model, the states are
tabulated below. The cost of each state includes the annual cost of AZT
(Zidovudine).

Name | Cost |
---|---|

A | 5034 |

B | 5330 |

C | 11285 |

D | 0 |

The per-cycle transition probabilities, which are the cells of the
Markov transition matrix, can be extracted from the model via the
function `transition_probabilities`

. For the monotherapy
model, the transition matrix is shown below. This is consistent with the
Table 1 of Chancellor *et al*^{2}.

A | B | C | D | |
---|---|---|---|---|

A | 0.7215 | 0.2018 | 0.0669 | 0.009804 |

B | 0 | 0.5811 | 0.407 | 0.01192 |

C | 0 | 0 | 0.7501 | 0.2499 |

D | 0 | 0 | 0 | 1 |

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, without half cycle correction, with discount.
In addition, the proportion of patients alive at each cycle (the Markov
trace) is added to the table. The populations and discounted costs are
consistent with Briggs *et al*, Table 2.3^{1}, and the QALY column is consistent
with Table 2.4 (without half cycle correction). No discount was applied
to the utilities.

```
# create starting populations
<- 1000L
N <- c(A = N, B = 0L, C = 0L, D = 0L)
populations $reset(populations)
m.mono# run 20 cycles
<- m.mono$cycles(ncycles = 20L, hcc.pop = FALSE, hcc.cost = FALSE) MT.mono
```

Years | A | B | C | D | Cost | QALY |
---|---|---|---|---|---|---|

0 | 1000 | 0 | 0 | 0 | 0 | 0 |

1 | 721 | 202 | 67 | 10 | 5153 | 0.99 |

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

7 | 102 | 114 | 341 | 444 | 3301 | 0.556 |

8 | 73 | 87 | 309 | 531 | 2708 | 0.469 |

9 | 53 | 65 | 272 | 610 | 2179 | 0.39 |

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

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 |

The estimated life years is approximated by summing the proportions
of patients left alive at each cycle (Briggs *et al*^{1}, 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
44663 GBP.

For combination therapy, a similar model was created, with adjusted
costs and transition rates. Following Briggs *et al*^{1} 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
<- RR*nAB/nA
pAB <- RR*nAC/nC
pAC <- RR*nAD/nA
pAD <- RR*nBC/nB
pBC <- RR*nBD/nB
pBD <- RR*nCD/nC
pCD # annual transition probability matrix
<- matrix(
Ptc 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")
)
)# create Markov states for combination therapy
<- MarkovState$new("A", cost = dmca + ccca + cAZT + cLam)
sAc <- MarkovState$new("B", cost = dmcb + cccb + cAZT + cLam)
sBc <- MarkovState$new("C", cost = dmcc + cccc + cAZT + cLam)
sCc <- MarkovState$new("D", cost = 0.0, utility = 0.0)
sDc # create transitions
<- Transition$new(sAc, sAc)
tAAc <- Transition$new(sAc, sBc)
tABc <- Transition$new(sAc, sCc)
tACc <- Transition$new(sAc, sDc)
tADc <- Transition$new(sBc, sBc)
tBBc <- Transition$new(sBc, sCc)
tBCc <- Transition$new(sBc, sDc)
tBDc <- Transition$new(sCc, sCc)
tCCc <- Transition$new(sCc, sDc)
tCDc <- Transition$new(sDc, sDc)
tDDc # construct the model
<- SemiMarkovModel$new(
m.comb V = list(sAc, sBc, sCc, sDc),
E = list(tAAc, tABc, tACc, tADc, tBBc, tBCc, tBDc, tCCc, tCDc, tDDc),
discount.cost = cDR / 100.0,
discount.utility = oDR / 100.0
)# set the probabilities
$set_probabilities(Ptc) m.comb
```

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

A | B | C | D | |
---|---|---|---|---|

A | 0.8585 | 0.1027 | 0.03376 | 0.00499 |

B | 0 | 0.7868 | 0.2072 | 0.006069 |

C | 0 | 0 | 0.8728 | 0.1272 |

D | 0 | 0 | 0 | 1 |

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 8 years. The
`reset`

function is used to set the cycle number and elapsed
time of the new run of the mono model.

```
# run combination therapy model for 2 years
<- c("A" = N, "B" = 0L, "C" = 0L, "D" = 0L)
populations $reset(populations)
m.comb# run 2 cycles
<- m.comb$cycles(2L, hcc.pop = FALSE, hcc.cost = FALSE)
MT.comb # feed populations into mono model & reset cycle counter and time
<- m.comb$get_populations()
populations $reset(
m.mono
populations, icycle = 2L,
elapsed = as.difftime(365.25*2.0, units = "days")
)# and run model for next 18 years
<- rbind(
MT.comb $cycles(ncycles = 18L, hcc.pop = FALSE, hcc.cost = FALSE)
MT.comb, m.mono )
```

The Markov trace for combination therapy is as follows:

Years | A | B | C | D | Cost | QALY |
---|---|---|---|---|---|---|

0 | 1000 | 0 | 0 | 0 | 0 | 0 |

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

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

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 |

The ICER is calculated by running both models and calculating the
incremental cost per life year gained. 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 6264
GBP/QALY. This is consistent with the result obtained by Briggs *et
al*^{1} (6276 GBP/QALY),
within rounding error.

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. *Med Decis
Making* **25,** 633–645 (2005).

4.

Miller, D. K. & Homan, S. M. Determining
Transition Probabilities: Confusion and
Suggestions. *Med Decis Making*
**14,** 52–58 (1994).

5.

Jones, E., Epstein, D. & García-Mochón, L.
A
Procedure for Deriving Formulas to
Convert Transition Rates to Probabilities for
Multistate Markov Models. *Med Decis Making*
**37,** 779–789 (2017).