---
title: "Simple Markov Models (Homogeneous)"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{Simple Markov Models (Homogeneous)}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r, echo=FALSE, include=FALSE}
library(heemod)
library(ggplot2)
```
The most simple Markov models in health economic evaluation are models were transition probabilities between states do not change with time. Those are called *homogeneous* or *time-homogeneous* Markov models.
## Model description
In this example we will model the cost effectiveness of lamivudine/zidovudine combination therapy in HIV infection ([Chancellor, 1997](https://pubmed.ncbi.nlm.nih.gov/10169387/) further described in [Decision Modelling for Health Economic Evaluation](https://global.oup.com/academic/product/decision-modelling-for-health-economic-evaluation-9780198526629/), page 32. For the sake of simplicity we will not reproduce exactly the analysis from the book. See vignette `vignette("i-reproduction", "heemod")` for an exact reproduction of the analysis.
This model aims to compare costs and utilities of two treatment strategies, *monotherapy* and *combined therapy*.
Four states are described, from best to worst health-wise:
* __A__: CD4 cells > 200 and < 500 cells/mm3;
* __B__: CD4 < 200 cells/mm3, non-AIDS;
* __C__: AIDS;
* __D__: Death.
## Transition probabilities
Transition probabilities for the monotherapy study group are rather simple to implement with `define_transition()`:
```{r}
mat_mono <- define_transition(
.721, .202, .067, .010,
0, .581, .407, .012,
0, 0, .750, .250,
0, 0, 0, 1
)
mat_mono
```
The combined therapy group has its transition probabilities multiplied by $rr = 0.509$, the relative risk of event for the population treated by combined therapy. Since $rr < 1$, the combined therapy group has less chance to transition to worst health states.
The probabilities to stay in the same state are equal to $1 - \sum P_{trans}$ where $P_{trans}$ are the probabilities to change to another state (because all transition probabilities from a given state must sum to 1).
We use the alias `C` as a convenient way to specify the probability complement, equal to $1 - \sum P_{trans}$.
```{r}
rr <- .509
mat_comb <- define_transition(
C, .202*rr, .067*rr, .010*rr,
0, C, .407*rr, .012*rr,
0, 0, C, .250*rr,
0, 0, 0, 1
)
mat_comb
```
We can plot the transition matrix for the monotherapy group:
```{r, fig.width = 6, fig.height=6, fig.align='center'}
plot(mat_mono)
```
And the combined therapy group:
```{r, fig.width = 6, fig.height=6, fig.align='center'}
plot(mat_comb)
```
## State values
The costs of lamivudine and zidovudine are defined:
```{r}
cost_zido <- 2278
cost_lami <- 2086
```
In addition to drugs costs (called `cost_drugs` in the model), each state is associated to healthcare costs (called `cost_health`). Cost are discounted at a 6% rate with the `discount` function.
Efficacy in this study is measured in terms of life expectancy (called `life_year` in the model). Each state thus has a value of 1 life year per year, except death who has a value of 0. Life-years are not discounted in this example.
Only `cost_drug` differs between the monotherapy and the combined therapy treatment groups, the function `dispatch_strategy()` can be used to account for that. For example state A can be defined with `define_state()`:
```{r}
state_A <- define_state(
cost_health = discount(2756, .06),
cost_drugs = discount(dispatch_strategy(
mono = cost_zido,
comb = cost_zido + cost_lami
), .06),
cost_total = cost_health + cost_drugs,
life_year = 1
)
state_A
```
The other states for the monotherapy treatment group can be specified in the same way:
```{r}
state_B <- define_state(
cost_health = discount(3052, .06),
cost_drugs = discount(dispatch_strategy(
mono = cost_zido,
comb = cost_zido + cost_lami
), .06),
cost_total = cost_health + cost_drugs,
life_year = 1
)
state_C <- define_state(
cost_health = discount(9007, .06),
cost_drugs = discount(dispatch_strategy(
mono = cost_zido,
comb = cost_zido + cost_lami
), .06),
cost_total = cost_health + cost_drugs,
life_year = 1
)
state_D <- define_state(
cost_health = 0,
cost_drugs = 0,
cost_total = 0,
life_year = 0
)
```
## Strategy definitions
Strategies can now be defined by combining a transition matrix and a state list with `define_strategy()`:
```{r}
strat_mono <- define_strategy(
transition = mat_mono,
state_A,
state_B,
state_C,
state_D
)
strat_mono
```
For the combined therapy model:
```{r}
strat_comb <- define_strategy(
transition = mat_comb,
state_A,
state_B,
state_C,
state_D
)
```
## Running the model
Both strategies can then be combined in a model and run for 50 years with `run_model()`. Strategies are given names (`mono` and `comb`) in order to facilitate result interpretation.
```{r}
res_mod <- run_model(
mono = strat_mono,
comb = strat_comb,
cycles = 50,
cost = cost_total,
effect = life_year
)
```
By default models are run for 1000 persons starting in the first state (here state **A**).
## Result interpretation
Strategy values can then be compared with `summary()` (optionally net monetary benefits can be calculated with the `threshold` option):
```{r}
summary(res_mod,
threshold = c(1000, 5000, 6000, 1e4))
```
The incremental cost-effectiveness ratio of the combined therapy strategy is thus £`r round(summary(res_mod)$res_comp$.icer[2])` per life-year gained.
The counts per state can be plotted by model:
```{r, fig.align='center', fig.width=6, fig.height=6, message=FALSE}
plot(res_mod, type = "counts", panel = "by_strategy") +
xlab("Time") +
theme_bw() +
scale_color_brewer(
name = "State",
palette = "Set1"
)
```
Or by state:
```{r, fig.align='center', fig.width=6, fig.height=8, message=FALSE}
plot(res_mod, type = "counts", panel = "by_state") +
xlab("Time") +
theme_bw() +
scale_color_brewer(
name = "Strategy",
palette = "Set1"
)
```
The values can also be represented:
```{r, fig.align='center', fig.width=6, fig.height=8, message=FALSE}
plot(res_mod, type = "values", panel = "by_value",
free_y = TRUE) +
xlab("Time") +
theme_bw() +
scale_color_brewer(
name = "Strategy",
palette = "Set1"
)
```
Note that classic `ggplot2` syntax can be used to modify plot appearance.