markov-models.Rmd

# Markov Models

Introduction

## Standard Matrix Formulation

## Simulation

#### 1. A simple decision tree {#C1}

This example is taken from @Hazen2014.
The problem is concerned with a competing risk cancer and AIDS decision tree.
We will assume discrete time of single years.
An individual starts in the `Well` state.
They can transition into `Dead`, `Cancer & AIDS`, `Cancer`, `AIDS` or remain in the `Well` state.

Define the transition probabilities:

* Die from other causes: $\delta_0 = 0.001182$
* Die from recurent prostate cancer: $\delta_c = 0.025$
* Die from AIDS: $\delta_a = 0.080$
* Cancer recurs: $\beta_c = 0.0027$
* Develop AIDS: $\beta_a = 0.0083$

Each state has an associated utility or benefit (quality factor in @Hazen2014) accrued by spending one cycle in each state.
Define the state utilities:

* `Well`: $R_w=$ 1.0
* `Cancer`: $R_c=$ 0.60
* `AIDS`: $R_a=$ 0.50
* `Cancer & AIDS`: $R_{ca}=$ 0.30
* `Dead`: $R_d=$ 0

Note that we will not include discounting.

<br>
<style>
div.blue { background-color:#e6f0ff; border-radius: 2px; padding: 5px;}
</style>
<div class = "blue">
C1. Define a (single year) decision tree and calculate the expected quality-adjusted value.
</div>

<br><br>

#### 2. Markov-cycle tree {#C2}

A Markov-cycle tree was introduced by @Hollenberg1984 and is a representation of a Markov process in which the possible events taking place during each cycle are represented by a probability tree.
This is one way of simplifying determining probabilities from multiple paths.

The diagram for the Markov-cycle tree of the example in @Hazen2014 is given below (note that the order of the states is different on the left-hand side and right-hand side).

```{r echo=FALSE, out.width="75%"}

suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(png))

include_graphics("figs/markov_cycle_tree.png")
```


The terminal state are now root or source states, meaning the process returns to the left-hand side to be repeated.

<br>
<style>
div.blue { background-color:#e6f0ff; border-radius: 2px; padding: 5px;}
</style>
<div class = "blue">
C2. Extend the model of C1 for multiple cycles and thus create a Markov-cycle tree. Calculate the mean quality-adjusted lifetime of 90.473.
</div>
<br><br>

#### 3. One-cycle Markov-cycle tree {#C3}

We can rearrange the Markov-cycle tree to closer resemble a Markov model by collapsing the branches into a single cycle and simply combining the probabilities.

In the below figure

* The numbers above each branch are the one-cycle transition probabilities
* The numbers pointing at nodes and names are the mean quality-adjusted durations accrued through $n$ cycles.
* The numbers in brackets are the mean quality-adjusted durations at the start of the cycle.

So for the below figure, the right-most numbers are the mean quality-adjusted durations for cycle 2, the left-most numbers are the mean quality-adjusted durations for cycle 3 and the numbers in brackets are the mean quality-adjusted durations for cycle 1.
@Hazen2014 steps through this calculation in detail.

```{r echo=FALSE, out.width="65%"}
include_graphics("figs/onecycle_markovcycletree.png")
```

<br>
<style>
div.blue { background-color:#e6f0ff; border-radius: 2px; padding: 5px;}
</style>
<div class = "blue">
C3. Modify the model of C2 to create a one-cycle Markov-cycle tree. Calculate the mean quality-adjusted lifetime.
</div>
<br><br>

#### 4. Discrete-time Markov model {#C4}

Clearly, the Markov-cycle tree can also be represented as a discrete-time Markov model.
The transition probabilities can be calculated by combining relevant path probabilities from the decision tree as done for the one-cycle Markov-cycle tree.
The model is shown below (note that death is not shows for simplicity).

```{r, echo=FALSE, out.width="75%"}
suppressPackageStartupMessages(library(heemod))
suppressPackageStartupMessages(library(diagram))

mat_dim <- define_transition(
  state_names = c('Well', 'Cancer', 'Cancer &\n AIDS', "AIDS"),
  a,b,c,d, 
  0,e,f,0,
  0,0,g,0,
  0,0,h,i)

plot(mat_dim)
```

<br>
<style>
div.blue { background-color:#e6f0ff; border-radius: 2px; padding: 5px;}
</style>
<div class = "blue">
C4. Create the equivalent discrete-time Markov model to the one-cycle Markov-cycle tree. Calculate cumulative proportion of patient cycles in each state and take product with health utilities for each respectively to obtain the mean quality-adjusted lifetime.
</div>
<br><br>

#### 5. Roll back Markov-cycle tree {#C5}

A neat strength is that we can calculate the mean quality-adjusted lifetime using the one-cycle Markov-cycle tree representation without calculating the cumulative proportion of time of patient cycles in each health state.
This is done by rolling back using the recursive equation ([value iteration](https://en.wikipedia.org/wiki/Markov_decision_process#Value_iteration)):

$$
V_n(i) = R(i) + \sum_j p_{ij} V_{n-1}(j)
$$
where $V_n(i)$ are the values at node $i$ at step $n$, in our case the mean quality-adjusted lifetime.

<br>
<style>
div.blue { background-color:#e6f0ff; border-radius: 2px; padding: 5px;}
</style>
<div class = "blue">
C5. Calculate the mean quality-adjusted lifetime using the one-cycle Markov-cycle tree and value iteration.
</div>
<br><br>

#### 6. (BONUS CHALLENGE): Roll back stochastic tree {#C6}

So far we have only considered discrete time.
The Markov-cycle tree representation can be extended to continuous time as a _stochastic tree_ 
(see @Hazen2014 for details).
Probabilities are now replaced by rates.
This change is represented by zigzag lines in the diagrams.
This is clearly a more compact representation.

We can calculate mean quality-adjusted lifetime in an analogous way to the discrete-time case by rolling back using the recursive equation:

$$
V(S) = \frac{R(i)}{\sum_j \lambda_j} + \sum_j p_j V(S_j)
$$
The new model diagram is given below.

```{r echo=FALSE, out.width="75%"}
include_graphics("figs/stochastic_tree.png")
```

The rates for state transitions are:

* `Cancer`: $\lambda_c = 0.03250$/year
* `AIDS`: $\lambda_a = 0.10$/year
* `Dead from Cancer`: $\mu_c = 0.3081$/year
* `Dead from AIDS`: $\mu_a = 0.9970$/year
* `Dead other`: $\mu_o = 0.014191$/year

<br>
<style>
div.blue { background-color:#e6f0ff; border-radius: 2px; padding: 5px;}
</style>
<div class = "blue">
C6. Create the stochastic tree model and calculate the mean quality-adjusted lifetime using value iteration.
</div>
<br><br>


&nbsp;

### References





Load supporting packages.

```{r}
library(purrr)
library(knitr)
```

#### C1. A simple decision tree

Define the separate monthly transition probabilities.

```{r}
delta0 <- 0.001182
deltac <- 0.025
deltaa <- 0.08
betac  <- 0.0027
betaa  <- 0.0083
```

and the utilities of being in each state.

```{r}
sutils <- c(dead = 0, CA = 0.3, cancer = 0.6, AIDS = 0.5, well = 1)
```

Note that the order of states is `dead`, `cancer & AIDS`, `cancer`, `AIDS`, `well`.

Define decision trees in terms of the structure, values and probabilities.
We'll use the `tribble` function just because it allows us to specify a matrix by rows rather than columns.

```{r define_tree_probs}
library(tibble)

tree_probs <- list()

# unique state outcomes each branch
tree_probs$well <-
  tribble(~rowname,  ~dead, ~ndead,  ~recurc, ~ncancer, ~CA,  ~cancer, ~AIDS, ~well,
          "well0",   delta0,1-delta0,NA,      NA,       NA,   NA,      NA,    NA,
          "ndead",   NA,    NA,      betac,   1-betac,  NA,   NA,      NA,    NA,
          "recurc",  NA,    NA,      NA,      NA,       betaa,1-betaa, NA,    NA,
          "ncancer", NA,    NA,      NA,      NA,       NA,   NA,      betaa, 1-betaa) %>% 
  column_to_rownames()

tree_probs$cancer <-
  tribble(~rowname, ~dead, ~ndead,  ~diec, ~survc,  ~CA,  ~cancer,
          "cancer0",delta0,1-delta0,NA,    NA,      NA,   NA,
          "ndead",  NA,    NA,      deltac,1-deltac,NA,   NA,
          "survc",  NA,    NA,      NA,    NA,      betaa,1-betaa) %>% 
  column_to_rownames()

tree_probs$AIDS <-
  tribble(~rowname,~dead, ~ndead,  ~diea, ~surva,  ~CA,  ~AIDS,
          "AIDS0", delta0,1-delta0,NA,    NA,      NA,   NA,
          "ndead", NA,    NA,      deltaa,1-deltaa,NA,   NA,
          "surva", NA,    NA,      NA,    NA,      betac,1-betac) %>% 
  column_to_rownames()

tree_probs$CA <-
  tribble(~rowname,~dead, ~ndead,  ~diec, ~survc,  ~diea, ~CA,
          "CA0",   delta0,1-delta0,NA,    NA,      NA,    NA,
          "ndead", NA,    NA,      deltac,1-deltac,NA,    NA,
          "survc", NA,    NA,      NA,    NA,      deltaa,1-deltaa) %>% 
  column_to_rownames()

tree_probs$dead <- 1

tree_probs
```

```{r define_tree_vals, eval = FALSE}
tree_vals <- list()

tree_vals$well <-
  tribble(~rowname, ~dead, ~ndead, ~cancer, ~ncancer, ~CA, ~CnotA, ~AIDS, ~well,
          "well0", 0,0,0,0,0,0,0,0,
          "ndead", 0,0,0,0,0,0,0,0,
          "cancer", 0,0,0,0,0.3,0.6,0,0,
          "ncancer", 0,0,0,0,0,0,0.5,1) %>% 
  column_to_rownames()

tree_vals$cancer <-
  tribble(~rowname, ~dead, ~ndead, ~cancer, ~ncancer, ~CA, ~CnotA, ~AIDS, ~well,
          "well", 0,0,0,0,0,0,0,0,
          "ndead", 0,0,0,0,0,0,0,0,
          "cancer", 0,0,0,0,0.3,0.6,0,0,
          "ncancer", 0,0,0,0,0,0,0.5,1) %>% 
  column_to_rownames()

tree_vals$AIDS <-
  tribble(~rowname, ~dead, ~ndead, ~cancer, ~ncancer, ~CA, ~CnotA, ~AIDS, ~well,
          "well", 0,0,0,0,0,0,0,0,
          "ndead", 0,0,0,0,0,0,0,0,
          "cancer", 0,0,0,0,0.3,0.6,0,0,
          "ncancer", 0,0,0,0,0,0,0.5,1) %>% 
  column_to_rownames()

tree_vals$CA <-
  tribble(~rowname, ~dead, ~ndead, ~cancer, ~ncancer, ~CA, ~CnotA, ~AIDS, ~well,
          "well", 0,0,0,0,0,0,0,0,
          "ndead", 0,0,0,0,0,0,0,0,
          "cancer", 0,0,0,0,0.3,0.6,0,0,
          "ncancer", 0,0,0,0,0,0,0.5,1) %>% 
  column_to_rownames()

```


#### C2. Extend C1 for multiple cycles

Assuming a binomial tree we can forward simulate for a synthetic cohort. This is a brute force approach and is potentially time-consuming.

```{r}
cohort <- list()
n_cohort <- 1000
death_states <- c("diea", "diec", "dead")

for (i in seq_len(n_cohort)) {
  
  traj_s <- NULL
  traj_u <- NULL
  state_name <- "well"
  
  while (!state_name %in% death_states) {
    
    p <- tree_probs[[state_name]]
    binp <- p[state_name, !is.na(p[state_name, ])] #partial match
    
    while (nrow(binp) > 0) {
      
      state_name <- 
        if (runif(1) < binp[1]) names(binp)[1] else names(binp)[2]
      
      binp <- p[state_name == rownames(p), !is.na(p[state_name, ])]
    }
    
    traj_s <- c(traj_s, state_name)
    traj_u <- c(traj_u, sutils[state_name])
  }
  
  cohort[[i]] <- traj_u 
}
```


An example trajectory
```{r}
cohort[[1]]
```

The mean summary statistics should be close to the true expected value.
However, it appears to be pretty noisy and even for fairly large values (in terms of run time) it can be off by one or two.

```{r}
mean(map(cohort, sum, na.rm = TRUE) %>% unlist())
```


```{r eval= FALSE, echo = FALSE}
##TODO

vals <- tree_vals$well
p <- tree_probs$well

num_from_nodes <- nrow(vals)
num_to_nodes <- ncol(vals)

#'to' values
c_hat <- colSums(vals, na.rm = TRUE)

for (i in num_from_nodes:1) {
  
  total <- 0
  for (j in 1:num_to_nodes) {
    
    if (!is.na(p[i,j])) {
      
      total <- total + p[i,j]*c_hat[j]
    }
  }
  c_hat[i] <- total + c_hat[i]
  
}
```



#### C3. Markov-cycle tree

Given the following transition matrix

$$
\tiny
\left(
\begin{matrix}
1 & 0 & 0 & 0 & 0\\
(1 - \delta_0)\delta_c + (1-\delta_0)(1-\delta_c)\delta_a & (1-\delta_0)(1-\delta_c(1-\delta_a) & 0 & 0 & 0\\
\delta_0 + (1-\delta_0)\delta_c & (1-\delta_0)(1-\delta_c)\beta_a & (1-\delta_0)(1-\delta_c)(1-\beta_a) & 0 & 0\\
\delta_0 + (1-\delta_0)\delta_a & (1-\delta_0)\beta_c(1-\delta_a) & 0 & (1-\delta_0)(1-\beta_c)(1-\delta_a) & 0\\
\delta_0 & (1-\delta_0)\beta_c\beta_a & (1-\delta_0)\beta_c(1-\beta_a) & (1-\delta_0)(1-\beta_c)\beta_a & (1-\delta_0)(1-\beta_c)(1-\beta_a)
\end{matrix}
\right)
$$


Then define the transition matrix object 

```{r define_markov_cycle_tree}
p <- list()

p$dead <- c(1,0,0,0,0)

p$CA <-
  c(delta0 + (1-delta0)*deltac + (1-delta0)*(1-deltac)*deltaa, (1-delta0)*(1-deltac)*(1-deltaa),0,0,0)

p$cancer <-
  c(delta0 + (1-delta0)*deltac, (1-delta0)*(1-deltac)*betaa, (1-delta0)*(1-deltac)*(1-betaa),0,0)

p$AIDS <-
  c(delta0 + (1-delta0)*deltaa, (1-delta0)*betac*(1-deltaa), 0, (1-delta0)*(1-betac)*(1-deltaa), 0)

p$well <-
  c(delta0, (1-delta0)*betac*betaa, (1-delta0)*betac*(1-betaa), (1-delta0)*(1-betac)*betaa,
    (1-delta0)*(1-betac)*(1-betaa))

trans <- do.call(rbind, p)
```

Combine the tree data all together into a single list using a function.

```{r}
create_tree <- function(trans, utils) {
  
  if (!all(rowSums(trans) == 1)) stop("probabilities don't sum to one")
  if (nrow(trans) != ncol(trans)) stop("not square matrix")
  if (nrow(trans) != length(utils)) stop("utils length doesnt match transition matrix dimensions")
  
  colnames(trans) <- rownames(trans)
  names(utils) <- rownames(trans)
  
  list(trans = trans,
       utils = utils)
}

my_tree <- create_tree(trans, sutils)
```

Check the input data.
```{r}
str(my_tree)
kable(my_tree$trans, digits = 3)
```


Now we're ready to do the forward cycle.
Basically, using the same approach as above for the separate probabilities, we simulate individuals and then take an average.

```{r}
cohort <- list()
n_cohort <- 1000
p <- my_tree$trans

for (i in seq_len(n_cohort)) {
  
  traj_s <- NULL
  traj_u <- NULL
  state_name <- "well"
  
  while (state_name != "dead") {
    
    res <- rmultinom(n = 1, size = 1, prob = p[state_name, ])
    state_name <- rownames(res)[res[,1] == 1]
    traj_s <- c(traj_s, state_name)
    traj_u <- c(traj_u, sutils[state_name])
  }
  
  cohort[[i]] <- traj_u 
}
```

Here's an example trajectory.
```{r}
cohort[[1]]
```

The expected value is

```{r}
mean(map(cohort, sum, na.rm = TRUE) %>% unlist())
```


#### C4. Regular Markov model

After initialising values for the calculation, for each cycle the probability of being in each of the states is calculated using the transition matrix and the previous cycle state occupancy probabilities.
Similarly, the utilities associated with being in each state are calculated for each cycle.

```{r}
run_model <- function(tree,
                      probs,
                      n_cycles = 1000) {
  
  if (!is.matrix(probs))
    probs <- matrix(probs, nrow = 1)
  
  qalys <- NULL
  costs <- NULL
  
  for (i in seq_len(n_cycles)) {
    
    probs <- rbind(probs, probs[i, ] %*% tree$trans)
    qalys <- rbind(qalys, probs[i, ]*tree$utils)
  }
  
  list(probs = probs,
       qalys = qalys)
}
```

By summing over all cycles we obtain the total utilities for each state.
The total sum is the expected QALYs value for an individual starting in state `well` until `dead`.

```{r}
init_pop <- c(0,0,0,0,1)
res <- run_model(my_tree, init_pop)
colSums(res$qalys)
```

The total expected QALYs are therefore
```{r}
sum(res$qalys)
```


#### C5. Roll back Markov-cycle tree

Let's write a recursive function to do the value iteration.
Because this is more general than a simple binary tree we need to sum over the number of to-nodes at each step.
Also, we need to limit the number of recursions since this function would run until we get a stack overflow error.
In the original paper, Hazen (1992) gives a table for 1, 2, 3, 10, 100, 1000 cylces to show convergence.
Here we inlude a `limit` argument which exits the function call after a certain tree depth is reached.

```{r value_iteration, eval = TRUE}
value_iteration <- function(n,        # starting node number
                            R,        # `reward` (utility/quality factor)
                            p,        # transition probabilty matrix
                            cycle,    # exits recursion after limit tree depth
                            limit = 100) {
  
  ## sub NAs so returns for absorbing state
  p[p == 0] <- NA
  p[p == 1] <- NA

  to_node <- which(!is.na(p[n, ]))

  if (length(to_node) == 0 || cycle == limit) {
    return(R[n])
  }
  else {
    Vsum <- 0
    
    for (i in seq_along(to_node)) {
      Vsum <- Vsum + p[n, to_node[i]]*value_iteration(to_node[i], R, p,
                                                      cycle = cycle + 1, limit = limit)
    }
    
    return(R[n] + Vsum)
  }
}
```



```{r eval = FALSE, echo = FALSE}
# Test on a small tree
value_iteration(
  n = 1,
  R = c(0, 4, 3),
  p = matrix(c(NA, 1/2, 1/3,
               NA, NA, NA,
               NA, NA, NA), ncol = 3, byrow = TRUE),
  cycle = 1)
```

If we run this for the cycles in Table 2 in Hazen (1992), omitting the 1000 cycle because it takes too long to run, then we get the following same values.


```{r, eval = TRUE}
map_dbl(c(1,2,3,10,100),
        function(x) value_iteration(n = 5,
                                    R = my_tree$utils, p = my_tree$trans,
                                    cycle = 1, limit = x))
```

```{r, eval=FALSE, echo=FALSE}

# https://www.radford.edu/~nokie/classes/360/dp-memoized.html
##TODO:
#

val_iter <- function(n, p, R) {
  
  V <- matrix(0, nrow = n, ncol = 5)
  V[1, ] <- R
  V[2, ] <- R + p*t(R)
  
  return(val_iter_aux(n, 5, p, R, V))  
}

val_iter_aux <- function(n, i, p, R, V) {
  
  if (V[n, i] == 0) {
    
    for (j in 1:5){
      
      V[n, i] <- V[n, i] + p[i, j]*val_iter_aux(n - 1, j, p, R, V)
    }
  }
  
  return(V[n, i])
}

```
The problem with this representation is that it grows exponenitally with the number of cycles.
One way to speed things up is to do some of the calculation up-front so that we only do it once.
We can achieve this by nesting a second function inside of the initial as follows.

```{r}
value_iteration2 <- function(n,        # starting node number
                             R,        # `reward` (utility/quality factor)
                             p,        # transition probabilty matrix
                             cycle,    # exits recursion after limit tree depth
                             limit = 100) {
  
  ## sub NAs so returns for absorbing state
  p[p == 0] <- NA
  p[p == 1] <- NA

  to_node <- map(seq(nrow(p)), function(i){ which(!is.na(p[i,])) } )
 
  v_iter <- function(n, cycle) {

    if (length(to_node) == 0 || cycle == limit) {
      return(R[n])
    }
    else {
      Vsum <- 0   
      for (i in seq_along(to_node[[n]])) {
        Vsum <- Vsum + p[n, to_node[[n]][i]]*v_iter(to_node[[n]][i],cycle = cycle + 1)
      }
      return(R[n] + Vsum)
    }
  }
  
  return(v_iter(n, cycle))
}

```

Although we still have the original main problem this does appreciably improve run time.

```{r, eval=FALSE}
microbenchmark::microbenchmark(
  value_iteration(n = 5,
                  R = my_tree$utils,
                  p = my_tree$trans,
                  cycle = 1, limit = 200),
  value_iteration2(n = 5,
                   R = my_tree$utils,
                   p = my_tree$trans,
                   cycle = 1, limit = 200), times = 1)
```

```{r exact_solution, eval=F, echo=F}
compute_v <- function(n,        # starting node number
                      R,        # `reward` (utility/quality factor)
                      p) {      # transition probabilty matrix       
  v <- numeric(n)
  v[1] <- R[1]
  if (n < 2) return(v)
  for(i in seq(2,n)) {
    vsum <- R[i]
    for (j in seq(i-1)) {
       vsum <- vsum + p[i,j]*v[j]
    }
    v[i] <- vsum / (1-p[i,i])
  }
  return(v)
}

# could solve explicitly in matrix formulation
# V = (I - P)^{-1} R
# where I - P must be invertible i.e. non-zero determinant

```