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 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
<- MarkovState$new("A")
sA <- MarkovState$new("B")
sB <- MarkovState$new("C")
sC <- MarkovState$new("D", utility = 0.0)
sD # create transitions
<- Transition$new(sA, sA)
tAA <- Transition$new(sA, sB)
tAB <- Transition$new(sA, sC)
tAC <- Transition$new(sA, sD)
tAD <- Transition$new(sB, sB)
tBB <- Transition$new(sB, sC)
tBC <- Transition$new(sB, sD)
tBD <- Transition$new(sC, sC)
tCC <- Transition$new(sC, sD)
tCD <- Transition$new(sD, sD)
tDD # set discount rates
<- 6.0 # annual discount rate, costs (%)
cDR <- 0.0 # annual discount rate, benefits (%)
oDR # construct the model
<- SemiMarkovModel$new(
m 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
)
```

The costs and discount rates used in the model (1995 rates) are numerical constants, and are defined as follows.

```
# drug costs
<- 2278.0 # zidovudine drug cost
cAZT <- 2087.0 # lamivudine drug cost
cLam
# direct medical and community 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
# occupancy costs with monotherapy
<- dmca + ccca + cAZT
cAm <- dmcb + cccb + cAZT
cBm <- dmcc + cccc + cAZT
cCm
# occupancy costs with combination therapy
<- dmca + ccca + cAZT + cLam
cAc <- dmcb + cccb + cAZT + cLam
cBc <- dmcc + cccc + cAZT + cLam cCc
```

The treatment effect was estimated by Chancellor *et al*^{2} via a meta-analysis, and is defined
as follows:

`<- 0.509 RR `

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.

```
# transition counts
<- 1251L
nAA <- 350L
nAB <- 116L
nAC <- 17L
nAD <- 731L
nBB <- 512L
nBC <- 15L
nBD <- 1312L
nCC <- 437L
nCD # create transition matrix
<- nAA + nAB + nAC + nAD
nA <- nBB + nBC + nBD
nB <- nCC + nCD
nC <- matrix(
Ptm 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
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 a census 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 is shown in the figure
below.

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*^{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, 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
<- function(Ptm, cAm, cBm, cCm, hcc = FALSE) {
run_mono # create starting populations
<- 1000L
N <- c(A = N, B = 0L, C = 0L, D = 0L)
populations $reset(populations)
m# set costs
$set_cost(cAm)
sA$set_cost(cBm)
sB$set_cost(cCm)
sC# set transition probabilities
$set_probabilities(Ptm)
m# run 20 cycles
<- m$cycles(
tr 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 setafterthe 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:

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

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.

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
44,663 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")
) )
```

The resulting 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 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
<- function(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc, hcc = FALSE) {
run_comb # set populations
<- 1000L
N <- c(A = N, B = 0L, C = 0L, D = 0L)
populations $reset(populations)
m# set the transition probabilities accounting for treatment effect
$set_probabilities(Ptc)
m# set the costs including those for the additional drug
$set_cost(cAc)
sA$set_cost(cBc)
sB$set_cost(cCc)
sC# run first 2 yearly cycles with additional drug costs and tx effect
<- m$cycles(2L, hcc.pop = hcc, hcc.cost = FALSE, hcc.QALY = hcc)
tr # save the state populations after 2 years
<- m$get_populations()
populations # revert probabilities to those without treatment effect
$set_probabilities(Ptm)
m# revert costs to those without the extra drug
$set_cost(cAm)
sA$set_cost(cBm)
sB$set_cost(cCm)
sC# restart the model with populations from first 2 years with extra drug
$reset(
m
populations,icycle = 2L,
elapsed = as.difftime(365.25 * 2.0, units = "days")
)# run for next 18 years, combining the traces
<- rbind(
tr
tr,$cycles(ncycles = 18L, hcc.pop = hcc, hcc.cost = FALSE, hcc.QALY = hcc)
m
)# 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.

```
<- run_comb(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc)
MT.comb <- sum(MT.comb$QALY)
el.comb <- sum(MT.comb$Cost)
cost.comb <- (cost.comb - cost.mono) / (el.comb - el.mono) icer
```

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 |

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 al*^{1} (6276 GBP/QALY), within rounding
error.

With half-cycle correction applied to the state populations, the model can be recalculated as follows.

```
<- run_mono(Ptm, cAm, cBm, cCm, hcc = TRUE)
MT.mono.hcc <- sum(MT.mono.hcc$QALY)
el.mono.hcc <- sum(MT.mono.hcc$Cost)
cost.mono.hcc <- run_comb(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc, hcc = TRUE)
MT.comb.hcc <- sum(MT.comb.hcc$QALY)
el.comb.hcc <- sum(MT.comb.hcc$Cost)
cost.comb.hcc <- (cost.comb.hcc - cost.mono.hcc) / (el.comb.hcc - el.mono.hcc) icer.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.

In their Exercise 4.7, Briggs *et al*^{1} 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.

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)
<- GammaModVar$new("dmca", "GBP", shape = 1.0, scale = 1701.0)
dmca <- GammaModVar$new("dmcb", "GBP", shape = 1.0, scale = 1774.0)
dmcb <- GammaModVar$new("dmcc", "GBP", shape = 1.0, scale = 6948.0)
dmcc <- GammaModVar$new("ccca", "GBP", shape = 1.0, scale = 1055.0)
ccca <- GammaModVar$new("cccb", "GBP", shape = 1.0, scale = 1278.0)
cccb <- GammaModVar$new("cccc", "GBP", shape = 1.0, scale = 2059.0)
cccc
# occupancy costs with monotherapy
<- ExprModVar$new("cA", "GBP", rlang::quo(dmca + ccca + cAZT))
cAm <- ExprModVar$new("cB", "GBP", rlang::quo(dmcb + cccb + cAZT))
cBm <- ExprModVar$new("cC", "GBP", rlang::quo(dmcc + cccc + cAZT))
cCm
# occupancy costs with combination therapy
<- ExprModVar$new("cAc", "GBP", rlang::quo(dmca + ccca + cAZT + cLam))
cAc <- ExprModVar$new("cBc", "GBP", rlang::quo(dmcb + cccb + cAZT + cLam))
cBc <- ExprModVar$new("cCc", "GBP", rlang::quo(dmcc + cccc + cAZT + cLam)) cCc
```

The treatment effect is also represented by a `ModVar`

whose uncertainty follows a log normal distribution.

```
<- LogNormModVar$new(
RR "Tx effect", "RR", p1 = 0.509, p2 = (0.710 - 0.365) / (2.0 * 1.96), "LN7"
)
```

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
<- function() {
pt_prob # create Dirichlet distributions for conditional probabilities
<- DirichletDistribution$new(c(1251L, 350L, 116L, 17L)) # from A # nolint
DA <- DirichletDistribution$new(c(731L, 512L, 15L)) # from B # nolint
DB <- DirichletDistribution$new(c(1312L, 437L)) # from C # nolint
DC # sample from the Dirichlet distributions
$sample()
DA$sample()
DB$sample()
DC# create the transition matrix
<- matrix(
Pt 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)
}
```

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
<- matrix(
psa 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
$set("random")
cAm$set("random")
cBm$set("random")
cCm$set("random")
cAc$set("random")
cBc$set("random")
cCc$set("random")
RR
# sample the probability transition matrix from observed counts
<- pt_prob()
Ptm
# run monotherapy model
<- run_mono(Ptm, cAm, cBm, cCm, hcc = TRUE)
MT.mono <- sum(MT.mono$QALY)
el.mono <- sum(MT.mono$Cost)
cost.mono "el.mono"]] <- el.mono
psa[[irun, "cost.mono"]] <- cost.mono
psa[[irun,
# create Pt for combination therapy (Briggs applied the RR to the transition
# probabilities - not recommended, but done here for reproducibility).
<- Ptm
Ptc for (i in 1L:4L) {
for (j in 1L:4L) {
<- ifelse(i == j, NA, RR$get() * Ptc[[i, j]])
Ptc[[i, j]]
}which(is.na(Ptc[i, ]))] <- 1.0 - sum(Ptc[i, ], na.rm = TRUE)
Ptc[i,
}
# run combination therapy model
<- run_comb(Ptm, cAm, cBm, cCm, Ptc, cAc, cBc, cCc, hcc = TRUE)
MT.comb <- sum(MT.comb$QALY)
el.comb <- sum(MT.comb$Cost)
cost.comb "el.comb"]] <- el.comb
psa[[irun, "cost.comb"]] <- cost.comb
psa[[irun,
# calculate the icer
"icer"]] <- (cost.comb - cost.mono) / (el.comb - el.mono)
psa[[irun, }
```

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.

The mean (95% confidence interval) for the cost of monotherapy was 45,007 (22,208 to 92,719) GBP, and the mean (95% CI) cost for combination therapy was 50,841 (27,858 to 96,642) GBP. The life years gained for monotherapy was 8.488 (8.105 to 8.928), and the life years gained for combination therapy was 9.429 (8.928 to 9.962). The mean ICER was 6,388 GBP/QALY with 95% confidence interval 2,863 to 11,089 GBP/QALY.

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