#| purl = FALSE, #| include = FALSE # read vignette source chunks from corresponding testthat script knitr::read_chunk(file.path("..", "tests", "testthat", "test-model-AZT.R")) # read vignette build utility functions knitr::read_chunk(file.path("vutils.R"))
#| purl = FALSE, #| include = FALSE knitr::opts_chunk$set( echo = FALSE, fig.keep = "last", fig.align = "center", collapse = TRUE, comment = "#>" )
#| gbp
#| purl = FALSE #nolint start
library(rdecision)
#| purl = FALSE #nolint end
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.
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.
#| model, #| echo = TRUE
The costs and discount rates used in the model (1995 rates) are numerical constants, and are defined as follows.
#| costs-det, #| echo = TRUE
The treatment effect was estimated by Chancellor et al [-@chancellor1997] via a meta-analysis, and is defined as follows:
#| txeffect-det, #| echo = TRUE
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.
#| pt-mono-det, #| echo = TRUE
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], Jones et al [-@jones2017] and others note, conversion between
rates and probabilities for multi-state Markov models is non-trivial
[@jones2017] 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.
#| fig.cap = "Markov model for comparison of HIV therapy. #| A: 200 < cd4 < 500, B: cd4 < 200, C: AIDS, D: Death.", #| purl = FALSE, #| eval = FALSE DiagrammeR::grViz(m$as_DOT())
#| gv2png
#| fig.cap = "Markov model for comparison of HIV therapy. #| A: 200 < cd4 < 500, B: cd4 < 200, C: AIDS, D: Death.", # images created with dot are more compact than DiagrammeR. pngfile <- gv2png( dot = m$as_DOT(rankdir = "TB", width = 7.0, height = 7.0) ) knitr::include_graphics(pngfile)
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].
#| purl = FALSE pander::pander(Ptm, emphasize.rownames = FALSE, justify = "lcccc")
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.
#| run-mono, #| echo = TRUE
Coding note: In function
run_mono
, the occupancy costs for states A, B and C are set via calls to functionset_cost()
which is associated with aMarkovState
object. Although these are set after the state objectssA
,sB
andsC
have been added to modelm
, 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 aMarkovState
changes, any other object that refers to it, such as aSemiMarkovModel
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:
#| mono-det-results, #| echo = TRUE
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.
#| purl = FALSE keep <- c("Years", "A", "B", "C", "D", "Cost", "QALY") pander::pander(MT.mono[, keep], row.names = FALSE, justify = "rrrrrrr", round = c(2L, 0L, 0L, 0L, 0L, 0L, 3L))
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.
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.
#| pt-comb-det, #| echo = TRUE
The resulting per-cycle transition matrix for the combination therapy is as follows:
#| purl = FALSE pander::pander(Ptc, emphasize.rownames = FALSE, justify = "lcccc")
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.
#| run-combo, #| echo = TRUE
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.
#| comb-det-results, #| echo = TRUE
The Markov trace for combination therapy is as follows:
#| purl = FALSE keep <- c("Years", "A", "B", "C", "D", "Cost", "QALY") pander::pander(MT.comb[, keep], row.names = FALSE, justify = "rrrrrrr", round = c(2L, 0L, 0L, 0L, 0L, 0L, 3L))
Over the 20 year time horizon, the expected life
years gained for monotherapy was r round(el.mono,3)
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.
With half-cycle correction applied to the state populations, the model can be recalculated as follows.
#| hcc-det, #| echo = TRUE
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.
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.
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.
#| costs-psa, #| echo = TRUE
The treatment effect is also represented by a ModVar
whose uncertainty follows
a log normal distribution.
#| txeffect-psa, #| echo = TRUE
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.
#| pt-psa, #| echo = TRUE
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.
#| run-psa, #| echo = TRUE
Coding note: The state occupancy costs
cAm
,cBm
etc. are nowModVar
s, rather than numeric variables as they were in the deterministic model. However, they can still be passed as arguments toMarkovState$set_cost()
, via the arguments to helper functionsrun_mono
andrun_comb
, andrdecision
will manage them appropriately, without changing any other code. Documentation for functions inrdecision
explains where this is supported by the package.
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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.