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