In feralaes/dampack: Package with useful decision-analytic modeling tools
knitr::opts_chunk$set(echo = TRUE)
devtools::load_all()
library(kableExtra) # https://haozhu233.github.io/kableExtra/awesome_table_in_html.html
## chunk will ensure that:
# all the figures generated by the report will be placed in the figs/sub-directory
# all the figures will be 6.5 x 4 inches and centered in the text.
knitr::opts_chunk$set(fig.width=6.5, fig.height=4, fig.align="center")
This document shows an example of how to conduct decision-analytic models in R [@Jalal2017b] using the R package dampack.
The R package dampack has a couple functions that could be used to evaluate Markov models. Specifically, the function CalculateMarkovTrace takes a transition probability matrix (time-homogeneous) or array (time-dependent) and then simulates a cohort of hypothetical individuals, and returns a cohort trace and transition array. The function ExpectedValue takes a transition array and rewards matrix, and calculates expected outcomes. The package can be used to compare competing strategies and calculate their cost-effectiveness. Below, we provide an applied example of these functions.
Model Description
We will develop a 4-state Markov model to quantify the expected outcomes (costs and quality-adjusted life years [QALYs]) for a hypothetical cohort of 55-year-old women who have undergone tumor excision for localized breast cancer (BCA). After surgery, all patients are initially tumor-free (i.e., will start in the "Local" Markov state). Each year, patients in the Local state face a 2% chance that they will experience a recurrence. Of these recurrences, 75% will present as metastatic disease (and patients will enter the "Mets"" Markov state), and the remainder present as local disease (and patients will enter the "Recur"" Markov state). Patients in the Recur state undergo repeat resection. Following this operation, they face a 6% chance of recurrence each year. (Assume that all local recurrences - 1st, 2nd, 3rd, etc. - have equivalent prognoses.) Of these recurrences, 90% will present as metastatic disease. Only patients with metastatic disease can die of breast cancer. Finally, all patients in the local and recur state face a 0.08 annual rate of mortality of dying from non-BCA related causes and patients in the mets state face a 0.5 anual rate of of dying of metastatic cancer in addition to dying from other causes. The state-transition diagram of this 4-state Markov model is shown in Figure \@ref(fig:state-transition-diag) and the Markov tree diagram that shows how transitions occur is shown in Figure \@ref(fig:Markov-tree-diag).
knitr::include_graphics("figs/Markov-Diagram-BRCA.pdf")
knitr::include_graphics("figs/Markov-tree-BRCA.pdf")
The Markov model described above represents usual care for women with breast cancer. Once we construct that Markov model, we will then add a second strategy of a hypothetical new treatment. The effect of this new treatment is to reduce the initial rate of cancer recurrence by 50%; however, the treatment also lowers the quality of life of women (disutility of 0.05) taking the medication because of side effects. Treatment is only taken by women who have not experienced a recurrence (i.e., those who reside in the Local state). The Markov model for the treatment strategy will be structurally identical to that of the usual care strategy.
After comparing the two strategies based on quality-adjusted life expectancy (QALE) and costs, we will incorporate costs and calculate the incremental cost per quality-adjusted life year (QALY) gained of treatment compared with usual care.
Markov structure: Usual care
The cancer model consists of 4-states (Local, Mets, Recur, and Dead). Therefore, we will be creating a 4x4 transition matrix. We start by labeling the states in the transition matrix.
state.names <- c("Local", "Recur", "Mets", "Dead")
We then define the number of states in the model.
n.s <- length(state.names)
From the decision problem, we know that patients in the local state face a 2% chance of recurrence under usual care (defined as p.BCA1). Of these recurrences, 75% will present as metastatic disease (and patients will enter the “Mets” Markov state; defined as p.Mets1), and the remainder present as local disease (and patients will enter the “Recur” Markov state; i.e., 1-p.Mets1). We will create variables to represent these probabilities of the transition matrix. Patients in the recur state face a 6% chance of recurrence each year (defined as p.BCA2). Of these recurrences, 90% will present as metastatic disease (defined as p.Mets2).
All patients in the local and recur state face a 0.076 annual rate of dying from non-BCA-related causes, $\mu_{Die}$ (defined as mu.Die). Patients in the Mets state face a 0.5 anual rate of of dying of metastatic cancer in addition to dying from other causes, $\mu_{DieMets}$ (defined as mu.DieMets).
To compute the annual probability of dying from non-BCA-related causes (defined as p.Die) we transform the annual rate assuming constant exponential rate
$$
p_{Die} = 1-exp(-\mu_{Die})
$$
The overall annual mortality rate for women with metastatic disease is equal to the sum of the rates of dying from other causes and the excess disease-specific mortality rate (defined as mu.DieMets). Therefore, the annual probability of dying in the Mets state (defined as p.DieMets) is converted by summing these rates and convert them to annual probabilities
$$
p_{DieMets} = 1-exp(-(\mu_{Die}+\mu_{DieMets}))
$$
The annual costs of breast cancer in the Local, Recur and Mets states (c.Local, c.Recur, and c.Mets, respectively) will be the same for both strategies. Table \@ref(tab:param-table) shows parameters of the model of usual care.
Table: (#tab:param-table) Table of parameters.
| Parameter | R name | Value |
|:-----------------------------------|:-----------:|:-------------:|
| Time horizon ($n_t$) | n.t | 51 years |
| Annual discount rate (costs/QALYs) | d.c/d.e | 3% |
| Annual probability of recurrence | p.BCA1 | 0.02 |
| after first local tumor diagnosis | | |
| under usual care | | |
| Proportion of metastatic disease | p.Mets1 | 0.75 |
| of recurrences from the Local state| | |
| Annual probability of recurrence | p.BCA2 | 0.06 |
| after > 1 local tumor diagnoses | | |
| Proportion of metastatic disease | p.Mets2 | 0.90 |
| of recurrences from the Recur state| | |
| Annual mortality | | |
| - All-cause mortality rate | mu.Die | 0.08 or age-specific |
| - Cancer specific mortality rate | mu.DieMets| 0.5 |
| Annual costs | | |
| - Cost of new treatment | c.Rx | 1000 |
| - Local state | c.Local | 500 |
| - Recur state | c.Recur | 5000 |
| - Mets state | c.Mets | 20000 |
| Utilities | | |
| - Local state under usual care | u.Local | 0.95 |
| - Recur state | u.Recur | 0.80 |
| - Mets state | u.Mets | 0.40 |
| - Dead state | u.Dead | 0.00 |
We input the parameters as follows:
n.t <- 51
d.c <- 0.03
d.e <- 0.03
p.BCA1 <- 0.02
p.Mets1 <- 0.75
p.BCA2 <- 0.06
p.Mets2 <- 0.90
mu.Die <- 0.08 # Mean hazard
p.Die <- 1 - exp(-mu.Die)
mu.DieMets <- 0.5
p.DieMets <- 1 - exp(-(mu.Die + mu.DieMets))
u.Local <- 0.95
u.Recur <- 0.80
u.Mets <- 0.40
u.Dead <- 0.00
c.Local <- 500
c.Recur <- 5000
c.Mets <- 20000
c.Rx <- 1000
c.Dead <- 0
Transition matrix
To construct the transition matrix for the usual care (UC) strategy, we first need to define a matrix that we will name m.P.UC. We initialize the matrix with NaN so then it can be filled with numeric elements.
# Initialize matrix
m.P.UC <- matrix(NaN,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.P.UC
Each row in m.P.UC represent a health state and each of its elements store the transition probbailities from that particular state. We need to manually store each of these transition probabilities into the matrix.
# Fill in matrix
# From Local
m.P.UC["Local", "Local"] <- (1 - p.Die) * (1 - p.BCA1)
m.P.UC["Local", "Recur"] <- (1 - p.Die) * p.BCA1 * (1 - p.Mets1)
m.P.UC["Local", "Mets"] <- (1 - p.Die) * p.BCA1 * p.Mets1
m.P.UC["Local", "Dead"] <- p.Die
# From Recur
m.P.UC["Recur", "Local"] <- 0
m.P.UC["Recur", "Recur"] <- (1 - p.Die) * ((1 - p.BCA2) + (p.BCA2 * (1 - p.Mets2)))
m.P.UC["Recur", "Mets"] <- (1 - p.Die) * p.BCA2 * p.Mets2
m.P.UC["Recur", "Dead"] <- p.Die
# From Mets
m.P.UC["Mets", "Local"] <- 0
m.P.UC["Mets", "Recur"] <- 0
m.P.UC["Mets", "Mets"] <- 1 - p.DieMets
m.P.UC["Mets", "Dead"] <- p.DieMets
# From Dead
m.P.UC["Dead", "Local"] <- 0
m.P.UC["Dead", "Recur"] <- 0
m.P.UC["Dead", "Mets"] <- 0
m.P.UC["Dead", "Dead"] <- 1
m.P.UC
We assume that all women start in the Local state, therefore the initial state vector defined as `s0'
s0 <- c(1, 0, 0, 0) # Everyone starts in the local state
Run Markov model
All of the components necessary to run the Markov model are now in place. We use the function CalculateMarkovTrace from the dampack package to run the model. The function requires a transition probability matrix, m.P.UC, in this case, the initial state vector (s0), and the number of cycles we want the model to run, n.t. The function produces a list with the cohort trace and the transition array, which we will call m.Markov.UC for the UC strategy.
l.Markov.UC <- CalculateMarkovTrace(M = m.P.UC, p0 = s0, n.cycles = n.t)
Rewards matrices
To calculate life years, quality adjusted life years, and cost, we will specify a reward matrix for each of these outcomes. The elements of each of the columns of these matrices must contain the state rewards of the health state corresponding to that column.
For the life years, for each year spent in any of the non-dead states, the cohort will get one year of life. That is, the elements of all the columns corresponding non-dead states equal to one.
m.LY.UC <- matrix(rep(c(1, 1, 1, 0), (n.s - 1)),
byrow = TRUE,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.LY.UC
For the QALY, we follow the similar procedure for the LY's reward matrix but now each column has the utilities that correspond with that particular state
m.QALY.UC <- matrix(rep(c(u.Local, u.Recur, u.Mets, u.Dead), n.s),
byrow = TRUE,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.QALY.UC
m.C.UC <- matrix(rep(c(c.Local, c.Recur, c.Mets, c.Dead), n.s),
byrow = TRUE,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.C.UC
Compute expected values
Using the ExpectedValue function from the dampack package, we calculate the expected life years of the Markov cohort under usual care. In the expected value function we input the transition array (generated from the CalculateMarkovTrace function), the reward matrix for life years m.LY.UC, discount rate d.e and ask to conduct a half cycle correction.
ev.LY.UC <- ExpectedValue(trans = l.Markov.UC$trans,
rwd = m.LY.UC,
r = d.e, half = TRUE)
ev.LY.UC
To compute the expected QALYs, we follow the same procedure but now using the reward matrix for QALYs m.QALY.UC
ev.QALY.UC <- ExpectedValue(trans = l.Markov.UC$trans,
rwd = m.QALY.UC,
r = d.e, half = T)
ev.QALY.UC
To compute the expected costs, we follow the same procedure but now using the reward matrix for costs m.C.UC
ev.C.UC <- ExpectedValue(trans = l.Markov.UC$trans,
rwd = m.C.UC,
r = d.c, half = T)
ev.C.UC
Markov structure: New treatment
We will now model the effect of a hypothetical new treatment. This new treatment competes with usual care. From the decision problem, we know that the treatment is only taken by women who have not experienced a recurrence (i.e., those who reside in the Local state).
The new treatment reduces the initial rate of cancer recurrence by 50% (i.e., a hazard ratio or recurrence under new treatment, defined as hr.Recur). To compute the probability of cancer recurrence associated with treatment (defined as p.BCA1.Rx), we first transform the probability of cancer recurrence without treatment to a rate. We then multiply the rate of cancer recurrence without treatment by the treatment effect (hr.Recur = 0.50). Finally, we convert the rate back to a probability to represent the probability of cancer recurrence for those who receive treatment.
hr.Recur <- 0.5
p.BCA1.Rx <- 1 - exp(log(1 - p.BCA1) * hr.Recur)
Transition matrix
We will incorporate the effects of treatment in a new transition matrix (defined as m.P.Rx). The treatment transition matrix will incorporate the treatment effects which impact transitions from the local disease state to the disease recurrence state and from the local disease state to the metastatic disease state. All other transitions in the model are identical to those in the usual care transition matrix (Table \@ref(tab:param-table)).
As before, each row in the transition matrix represents a health state and each element of the matrix store the transition probabilities. Ultimately, we will be able to compare the expected outcomes between the usual care and treatment models. We manually enter the transition probabilities into the matrix.
# Initialize matrix
m.P.Rx <- matrix(NaN,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.P.Rx
# Fill in matrix
# From Local
m.P.Rx["Local", "Local"] <- (1 - p.Die) * (1 - p.BCA1.Rx)
m.P.Rx["Local", "Recur"] <- (1 - p.Die) * p.BCA1.Rx * (1 - p.Mets1)
m.P.Rx["Local", "Mets"] <- (1 - p.Die) * p.BCA1.Rx * p.Mets1
m.P.Rx["Local", "Dead"] <- p.Die
# From Recur
m.P.Rx["Recur", "Local"] <- 0
m.P.Rx["Recur", "Recur"] <- (1 - p.Die) * ((1 - p.BCA2) + (p.BCA2 * (1 - p.Mets2)))
m.P.Rx["Recur", "Mets"] <- (1 - p.Die) * p.BCA2 * p.Mets2
m.P.Rx["Recur", "Dead"] <- p.Die
# From Mets
m.P.Rx["Mets", "Local"] <- 0
m.P.Rx["Mets", "Recur"] <- 0
m.P.Rx["Mets", "Mets"] <- 1 - p.DieMets
m.P.Rx["Mets", "Dead"] <- p.DieMets
# From Dead
m.P.Rx["Dead", "Local"] <- 0
m.P.Rx["Dead", "Recur"] <- 0
m.P.Rx["Dead", "Mets"] <- 0
m.P.Rx["Dead", "Dead"] <- 1
m.P.Rx
Run Markov model
l.Markov.Rx <- CalculateMarkovTrace(M = m.P.Rx, p0 = s0, n.cycles = n.t)
Reward matrices
The LY's reward matrix for Rx is the same than for UC.
m.LY.Rx <- m.LY.UC
As noted in the decision problem, treatment lowers the quality of life (i.e., utility decrement) of women taking the medication due to side effects (defined as du.Local). We incorporate the disutility (0.05) associated with treatment into the reward matrix for the Local state of the Rx model.
du.Local <- 0.05
The reward matrix for Rx changes by adding the disutility from the treatment to the Local state.
m.QALY.Rx <- matrix(rep(c(u.Local - du.Local, u.Recur, u.Mets, u.Dead), n.s),
byrow = TRUE,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.QALY.Rx
m.C.Rx <- matrix(rep(c(c.Local + c.Rx, c.Recur, c.Mets, c.Dead), n.s),
byrow = TRUE,
nrow = n.s, ncol = n.s,
dimnames = list(state.names, state.names))
m.C.Rx
Compute expected values
ev.LY.Rx <- ExpectedValue(trans = l.Markov.Rx$trans,
rwd = m.LY.Rx,
r = d.e, half = TRUE)
ev.LY.Rx
To compute the expected QALYs, we follow the same procedure but now using the reward matrix for QALYs m.QALY.Rx
ev.QALY.Rx <- ExpectedValue(trans = l.Markov.Rx$trans,
rwd = m.QALY.Rx,
r = d.e, half = TRUE)
ev.QALY.Rx
To compute the expected costs, we follow the same procedure but now using the reward matrix for costs m.C.Rx
ev.C.Rx <- ExpectedValue(trans = l.Markov.Rx$trans,
rwd = m.C.Rx,
r = d.c, half = TRUE)
ev.C.Rx
Outputs
Table \@ref(tab:Output) shows the outputs of the Markov model for the two strategies.
markov.out <- data.frame(Strategy = c("Usual care", "New Treatment"),
LYs = round(c(ev.LY.UC, ev.LY.Rx), 2),
QALYs = round(c(ev.QALY.UC, ev.QALY.Rx), 2),
Costs = scales::dollar(c(ev.C.UC, ev.C.Rx)))
knitr::kable(markov.out,
booktabs = TRUE,
caption = "Outcomes of Markov model by strategy") %>%
kable_styling(latex_options = "hold_position")
Cost-effectiveness analysis
Finally, in Table \@ref(tab:CEA), we show the results for the CEA of the Markov model.
v.c <- c(ev.C.UC, ev.C.Rx)
v.e <- c(ev.QALY.UC, ev.QALY.Rx)
markov.cea <- data.frame(Strategy = c("Usual care", "New Treatment"),
QALYs = round(v.e, 2),
`Incremental QALYs` = c("-",
round(diff(v.e), 2)),
Costs = scales::dollar(v.c),
"Incremental Costs" = c("-",
scales::dollar(diff(v.c))),
ICER = c("-", paste(scales::dollar(diff(v.c)/diff(v.e)), "/QALY")))
knitr::kable(markov.cea,
booktabs = TRUE,
caption = "CEA of new treatment for breast cancer patients") %>%
kable_styling(latex_options = "hold_position")
References
feralaes/dampack documentation built on May 16, 2019, 12:48 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE) devtools::load_all() library(kableExtra) # https://haozhu233.github.io/kableExtra/awesome_table_in_html.html
## chunk will ensure that: # all the figures generated by the report will be placed in the figs/sub-directory # all the figures will be 6.5 x 4 inches and centered in the text. knitr::opts_chunk$set(fig.width=6.5, fig.height=4, fig.align="center")
This document shows an example of how to conduct decision-analytic models in R [@Jalal2017b] using the R package dampack.
The R package dampack has a couple functions that could be used to evaluate Markov models. Specifically, the function CalculateMarkovTrace takes a transition probability matrix (time-homogeneous) or array (time-dependent) and then simulates a cohort of hypothetical individuals, and returns a cohort trace and transition array. The function ExpectedValue takes a transition array and rewards matrix, and calculates expected outcomes. The package can be used to compare competing strategies and calculate their cost-effectiveness. Below, we provide an applied example of these functions.
Model Description
We will develop a 4-state Markov model to quantify the expected outcomes (costs and quality-adjusted life years [QALYs]) for a hypothetical cohort of 55-year-old women who have undergone tumor excision for localized breast cancer (BCA). After surgery, all patients are initially tumor-free (i.e., will start in the "Local" Markov state). Each year, patients in the Local state face a 2% chance that they will experience a recurrence. Of these recurrences, 75% will present as metastatic disease (and patients will enter the "Mets"" Markov state), and the remainder present as local disease (and patients will enter the "Recur"" Markov state). Patients in the Recur state undergo repeat resection. Following this operation, they face a 6% chance of recurrence each year. (Assume that all local recurrences - 1st, 2nd, 3rd, etc. - have equivalent prognoses.) Of these recurrences, 90% will present as metastatic disease. Only patients with metastatic disease can die of breast cancer. Finally, all patients in the local and recur state face a 0.08 annual rate of mortality of dying from non-BCA related causes and patients in the mets state face a 0.5 anual rate of of dying of metastatic cancer in addition to dying from other causes. The state-transition diagram of this 4-state Markov model is shown in Figure \@ref(fig:state-transition-diag) and the Markov tree diagram that shows how transitions occur is shown in Figure \@ref(fig:Markov-tree-diag).
knitr::include_graphics("figs/Markov-Diagram-BRCA.pdf")
knitr::include_graphics("figs/Markov-tree-BRCA.pdf")
The Markov model described above represents usual care for women with breast cancer. Once we construct that Markov model, we will then add a second strategy of a hypothetical new treatment. The effect of this new treatment is to reduce the initial rate of cancer recurrence by 50%; however, the treatment also lowers the quality of life of women (disutility of 0.05) taking the medication because of side effects. Treatment is only taken by women who have not experienced a recurrence (i.e., those who reside in the Local state). The Markov model for the treatment strategy will be structurally identical to that of the usual care strategy.
After comparing the two strategies based on quality-adjusted life expectancy (QALE) and costs, we will incorporate costs and calculate the incremental cost per quality-adjusted life year (QALY) gained of treatment compared with usual care.
Markov structure: Usual care
The cancer model consists of 4-states (Local, Mets, Recur, and Dead). Therefore, we will be creating a 4x4 transition matrix. We start by labeling the states in the transition matrix.
state.names <- c("Local", "Recur", "Mets", "Dead")
We then define the number of states in the model.
n.s <- length(state.names)
From the decision problem, we know that patients in the local state face a 2% chance of recurrence under usual care (defined as p.BCA1). Of these recurrences, 75% will present as metastatic disease (and patients will enter the “Mets” Markov state; defined as p.Mets1), and the remainder present as local disease (and patients will enter the “Recur” Markov state; i.e., 1-p.Mets1). We will create variables to represent these probabilities of the transition matrix. Patients in the recur state face a 6% chance of recurrence each year (defined as p.BCA2). Of these recurrences, 90% will present as metastatic disease (defined as p.Mets2).
All patients in the local and recur state face a 0.076 annual rate of dying from non-BCA-related causes, $\mu_{Die}$ (defined as mu.Die). Patients in the Mets state face a 0.5 anual rate of of dying of metastatic cancer in addition to dying from other causes, $\mu_{DieMets}$ (defined as mu.DieMets).
To compute the annual probability of dying from non-BCA-related causes (defined as p.Die) we transform the annual rate assuming constant exponential rate
$$
p_{Die} = 1-exp(-\mu_{Die})
$$
The overall annual mortality rate for women with metastatic disease is equal to the sum of the rates of dying from other causes and the excess disease-specific mortality rate (defined as mu.DieMets). Therefore, the annual probability of dying in the Mets state (defined as p.DieMets) is converted by summing these rates and convert them to annual probabilities
$$ p_{DieMets} = 1-exp(-(\mu_{Die}+\mu_{DieMets})) $$
The annual costs of breast cancer in the Local, Recur and Mets states (c.Local, c.Recur, and c.Mets, respectively) will be the same for both strategies. Table \@ref(tab:param-table) shows parameters of the model of usual care.
Table: (#tab:param-table) Table of parameters.
| Parameter | R name | Value |
|:-----------------------------------|:-----------:|:-------------:|
| Time horizon ($n_t$) | n.t | 51 years |
| Annual discount rate (costs/QALYs) | d.c/d.e | 3% |
| Annual probability of recurrence | p.BCA1 | 0.02 |
| after first local tumor diagnosis | | |
| under usual care | | |
| Proportion of metastatic disease | p.Mets1 | 0.75 |
| of recurrences from the Local state| | |
| Annual probability of recurrence | p.BCA2 | 0.06 |
| after > 1 local tumor diagnoses | | |
| Proportion of metastatic disease | p.Mets2 | 0.90 |
| of recurrences from the Recur state| | |
| Annual mortality | | |
| - All-cause mortality rate | mu.Die | 0.08 or age-specific |
| - Cancer specific mortality rate | mu.DieMets| 0.5 |
| Annual costs | | |
| - Cost of new treatment | c.Rx | 1000 |
| - Local state | c.Local | 500 |
| - Recur state | c.Recur | 5000 |
| - Mets state | c.Mets | 20000 |
| Utilities | | |
| - Local state under usual care | u.Local | 0.95 |
| - Recur state | u.Recur | 0.80 |
| - Mets state | u.Mets | 0.40 |
| - Dead state | u.Dead | 0.00 |
We input the parameters as follows:
n.t <- 51 d.c <- 0.03 d.e <- 0.03 p.BCA1 <- 0.02 p.Mets1 <- 0.75 p.BCA2 <- 0.06 p.Mets2 <- 0.90 mu.Die <- 0.08 # Mean hazard p.Die <- 1 - exp(-mu.Die) mu.DieMets <- 0.5 p.DieMets <- 1 - exp(-(mu.Die + mu.DieMets)) u.Local <- 0.95 u.Recur <- 0.80 u.Mets <- 0.40 u.Dead <- 0.00 c.Local <- 500 c.Recur <- 5000 c.Mets <- 20000 c.Rx <- 1000 c.Dead <- 0
Transition matrix
To construct the transition matrix for the usual care (UC) strategy, we first need to define a matrix that we will name m.P.UC. We initialize the matrix with NaN so then it can be filled with numeric elements.
# Initialize matrix m.P.UC <- matrix(NaN, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.P.UC
Each row in m.P.UC represent a health state and each of its elements store the transition probbailities from that particular state. We need to manually store each of these transition probabilities into the matrix.
# Fill in matrix # From Local m.P.UC["Local", "Local"] <- (1 - p.Die) * (1 - p.BCA1) m.P.UC["Local", "Recur"] <- (1 - p.Die) * p.BCA1 * (1 - p.Mets1) m.P.UC["Local", "Mets"] <- (1 - p.Die) * p.BCA1 * p.Mets1 m.P.UC["Local", "Dead"] <- p.Die # From Recur m.P.UC["Recur", "Local"] <- 0 m.P.UC["Recur", "Recur"] <- (1 - p.Die) * ((1 - p.BCA2) + (p.BCA2 * (1 - p.Mets2))) m.P.UC["Recur", "Mets"] <- (1 - p.Die) * p.BCA2 * p.Mets2 m.P.UC["Recur", "Dead"] <- p.Die # From Mets m.P.UC["Mets", "Local"] <- 0 m.P.UC["Mets", "Recur"] <- 0 m.P.UC["Mets", "Mets"] <- 1 - p.DieMets m.P.UC["Mets", "Dead"] <- p.DieMets # From Dead m.P.UC["Dead", "Local"] <- 0 m.P.UC["Dead", "Recur"] <- 0 m.P.UC["Dead", "Mets"] <- 0 m.P.UC["Dead", "Dead"] <- 1 m.P.UC
We assume that all women start in the Local state, therefore the initial state vector defined as `s0'
s0 <- c(1, 0, 0, 0) # Everyone starts in the local state
Run Markov model
All of the components necessary to run the Markov model are now in place. We use the function CalculateMarkovTrace from the dampack package to run the model. The function requires a transition probability matrix, m.P.UC, in this case, the initial state vector (s0), and the number of cycles we want the model to run, n.t. The function produces a list with the cohort trace and the transition array, which we will call m.Markov.UC for the UC strategy.
l.Markov.UC <- CalculateMarkovTrace(M = m.P.UC, p0 = s0, n.cycles = n.t)
Rewards matrices
To calculate life years, quality adjusted life years, and cost, we will specify a reward matrix for each of these outcomes. The elements of each of the columns of these matrices must contain the state rewards of the health state corresponding to that column.
For the life years, for each year spent in any of the non-dead states, the cohort will get one year of life. That is, the elements of all the columns corresponding non-dead states equal to one.
m.LY.UC <- matrix(rep(c(1, 1, 1, 0), (n.s - 1)), byrow = TRUE, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.LY.UC
For the QALY, we follow the similar procedure for the LY's reward matrix but now each column has the utilities that correspond with that particular state
m.QALY.UC <- matrix(rep(c(u.Local, u.Recur, u.Mets, u.Dead), n.s), byrow = TRUE, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.QALY.UC
m.C.UC <- matrix(rep(c(c.Local, c.Recur, c.Mets, c.Dead), n.s), byrow = TRUE, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.C.UC
Compute expected values
Using the ExpectedValue function from the dampack package, we calculate the expected life years of the Markov cohort under usual care. In the expected value function we input the transition array (generated from the CalculateMarkovTrace function), the reward matrix for life years m.LY.UC, discount rate d.e and ask to conduct a half cycle correction.
ev.LY.UC <- ExpectedValue(trans = l.Markov.UC$trans, rwd = m.LY.UC, r = d.e, half = TRUE) ev.LY.UC
To compute the expected QALYs, we follow the same procedure but now using the reward matrix for QALYs m.QALY.UC
ev.QALY.UC <- ExpectedValue(trans = l.Markov.UC$trans, rwd = m.QALY.UC, r = d.e, half = T) ev.QALY.UC
To compute the expected costs, we follow the same procedure but now using the reward matrix for costs m.C.UC
ev.C.UC <- ExpectedValue(trans = l.Markov.UC$trans, rwd = m.C.UC, r = d.c, half = T) ev.C.UC
Markov structure: New treatment
We will now model the effect of a hypothetical new treatment. This new treatment competes with usual care. From the decision problem, we know that the treatment is only taken by women who have not experienced a recurrence (i.e., those who reside in the Local state).
The new treatment reduces the initial rate of cancer recurrence by 50% (i.e., a hazard ratio or recurrence under new treatment, defined as hr.Recur). To compute the probability of cancer recurrence associated with treatment (defined as p.BCA1.Rx), we first transform the probability of cancer recurrence without treatment to a rate. We then multiply the rate of cancer recurrence without treatment by the treatment effect (hr.Recur = 0.50). Finally, we convert the rate back to a probability to represent the probability of cancer recurrence for those who receive treatment.
hr.Recur <- 0.5 p.BCA1.Rx <- 1 - exp(log(1 - p.BCA1) * hr.Recur)
Transition matrix
We will incorporate the effects of treatment in a new transition matrix (defined as m.P.Rx). The treatment transition matrix will incorporate the treatment effects which impact transitions from the local disease state to the disease recurrence state and from the local disease state to the metastatic disease state. All other transitions in the model are identical to those in the usual care transition matrix (Table \@ref(tab:param-table)).
As before, each row in the transition matrix represents a health state and each element of the matrix store the transition probabilities. Ultimately, we will be able to compare the expected outcomes between the usual care and treatment models. We manually enter the transition probabilities into the matrix.
# Initialize matrix m.P.Rx <- matrix(NaN, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.P.Rx
# Fill in matrix # From Local m.P.Rx["Local", "Local"] <- (1 - p.Die) * (1 - p.BCA1.Rx) m.P.Rx["Local", "Recur"] <- (1 - p.Die) * p.BCA1.Rx * (1 - p.Mets1) m.P.Rx["Local", "Mets"] <- (1 - p.Die) * p.BCA1.Rx * p.Mets1 m.P.Rx["Local", "Dead"] <- p.Die # From Recur m.P.Rx["Recur", "Local"] <- 0 m.P.Rx["Recur", "Recur"] <- (1 - p.Die) * ((1 - p.BCA2) + (p.BCA2 * (1 - p.Mets2))) m.P.Rx["Recur", "Mets"] <- (1 - p.Die) * p.BCA2 * p.Mets2 m.P.Rx["Recur", "Dead"] <- p.Die # From Mets m.P.Rx["Mets", "Local"] <- 0 m.P.Rx["Mets", "Recur"] <- 0 m.P.Rx["Mets", "Mets"] <- 1 - p.DieMets m.P.Rx["Mets", "Dead"] <- p.DieMets # From Dead m.P.Rx["Dead", "Local"] <- 0 m.P.Rx["Dead", "Recur"] <- 0 m.P.Rx["Dead", "Mets"] <- 0 m.P.Rx["Dead", "Dead"] <- 1 m.P.Rx
Run Markov model
l.Markov.Rx <- CalculateMarkovTrace(M = m.P.Rx, p0 = s0, n.cycles = n.t)
Reward matrices
The LY's reward matrix for Rx is the same than for UC.
m.LY.Rx <- m.LY.UC
As noted in the decision problem, treatment lowers the quality of life (i.e., utility decrement) of women taking the medication due to side effects (defined as du.Local). We incorporate the disutility (0.05) associated with treatment into the reward matrix for the Local state of the Rx model.
du.Local <- 0.05
The reward matrix for Rx changes by adding the disutility from the treatment to the Local state.
m.QALY.Rx <- matrix(rep(c(u.Local - du.Local, u.Recur, u.Mets, u.Dead), n.s), byrow = TRUE, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.QALY.Rx
m.C.Rx <- matrix(rep(c(c.Local + c.Rx, c.Recur, c.Mets, c.Dead), n.s), byrow = TRUE, nrow = n.s, ncol = n.s, dimnames = list(state.names, state.names)) m.C.Rx
Compute expected values
ev.LY.Rx <- ExpectedValue(trans = l.Markov.Rx$trans, rwd = m.LY.Rx, r = d.e, half = TRUE) ev.LY.Rx
To compute the expected QALYs, we follow the same procedure but now using the reward matrix for QALYs m.QALY.Rx
ev.QALY.Rx <- ExpectedValue(trans = l.Markov.Rx$trans, rwd = m.QALY.Rx, r = d.e, half = TRUE) ev.QALY.Rx
To compute the expected costs, we follow the same procedure but now using the reward matrix for costs m.C.Rx
ev.C.Rx <- ExpectedValue(trans = l.Markov.Rx$trans, rwd = m.C.Rx, r = d.c, half = TRUE) ev.C.Rx
Outputs
Table \@ref(tab:Output) shows the outputs of the Markov model for the two strategies.
markov.out <- data.frame(Strategy = c("Usual care", "New Treatment"), LYs = round(c(ev.LY.UC, ev.LY.Rx), 2), QALYs = round(c(ev.QALY.UC, ev.QALY.Rx), 2), Costs = scales::dollar(c(ev.C.UC, ev.C.Rx))) knitr::kable(markov.out, booktabs = TRUE, caption = "Outcomes of Markov model by strategy") %>% kable_styling(latex_options = "hold_position")
Cost-effectiveness analysis
Finally, in Table \@ref(tab:CEA), we show the results for the CEA of the Markov model.
v.c <- c(ev.C.UC, ev.C.Rx) v.e <- c(ev.QALY.UC, ev.QALY.Rx) markov.cea <- data.frame(Strategy = c("Usual care", "New Treatment"), QALYs = round(v.e, 2), `Incremental QALYs` = c("-", round(diff(v.e), 2)), Costs = scales::dollar(v.c), "Incremental Costs" = c("-", scales::dollar(diff(v.c))), ICER = c("-", paste(scales::dollar(diff(v.c)/diff(v.e)), "/QALY"))) knitr::kable(markov.cea, booktabs = TRUE, caption = "CEA of new treatment for breast cancer patients") %>% kable_styling(latex_options = "hold_position")
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.