---
title: "A Tutorial for Time Dependent Markov Model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{A Tutorial for Time Dependent Markov Model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This vignette is based on our tutorial for time-dependent Markov models in R published in *Medical Decision Making*: 

*Alarid-Escudero, F., Krijkamp, E., Enns, E. A., Yang, A., Hunink, M. M., Pechlivanoglou, P., & Jalal, H. (2023). A tutorial on time-dependent cohort state-transition models in r using a cost-effectiveness analysis example. Medical Decision Making, 43(1), 21-41.*

First, we load `twig` library
```{r setup}
library(twig)
```

# Twig 

## model definition
Then, we define the model's structure `twig`
```{r}
n_cycles <- 25 # number of cycles

mytwig <- twig() + 
  decisions(names = c(StandardOfCare, StrategyA, StrategyB, StrategyAB)) + # define decisions
  states(names = c(H, S1, S2, D), # Markov state names
         init_probs = c(1,0,0,0), # everyone starts at H
         max_cycles = c(1,n_cycles, 1, 1)) + # the cohort can stay in S1 for n_cycles
  event(name = death_event,  # first event is death
        options = c(yes,none), # which 2 options
        probs = c(pDie, leftover), # probability function name and its complement
        transitions = c(D, second_event)) + # if death occurs go to D, otherwise, go to the next event (second_event)
  event(name = second_event, # the second event
        options = c(recover, getsick, progress, none), # has 4 options
        probs = c(pRecover, pGetSick, pProgress, leftover), # and 3 named probabilities and a complement 
        transitions = c(H, S1, S2, stay)) + # resulting in transitions to H, S1, S2 or else staying in the original state
  payoffs(names = c(cost, utility), # payoff names
          discount_rates = c(0.03, 0.03)) # payoff discount rates

```

## DecisionTwig
In [DecisionTwig](https://www.dashlab.ca/projects/decision_twig/), this Markov model looks like this:
![Markov](../man/figures/markov_twig.png){width=700px}

Next, we create a data frame of random samples from the model input parameters' distributions
``` {r}

n_sims <- 1000

# Create the data.table with n_sim rows of random samples
params <- data.frame(
  r_HS1         = rbeta(n_sims, 2, 10),             # Transition rate with beta distribution
  r_S1H         = rbeta(n_sims, 5, 5),              # Another transition rate with a different shape
  hr_S1         = rlnorm(n_sims, log(3), 0.2),      # Hazard ratio, log-normal to allow skewness
  hr_S2         = rlnorm(n_sims, log(10), 0.2),     # Higher hazard ratio, same distribution

  hr_S1S2_trtB  = rbeta(n_sims, 6, 4),              # Hazard ratio under treatment with beta distribution

  r_S1S2_scale  = rgamma(n_sims, shape = 2, rate = 25), # Scale parameter, gamma distribution
  r_S1S2_shape  = rgamma(n_sims, shape = 3, rate = 3),  # Shape parameter, gamma distribution

  c_H           = rnorm(n_sims, mean = 2000, sd = 50),   # Annual cost, slight variation for simulation
  c_S1          = rnorm(n_sims, mean = 4000, sd = 100),  # Higher annual cost, slightly varied
  c_S2          = rnorm(n_sims, mean = 15000, sd = 500), # Large cost with moderate variation
  c_D           = 0,                                        # Constant, no variation
  c_trtA        = rnorm(n_sims, mean = 12000, sd = 200), # Cost of treatment A with small variation
  c_trtB        = rnorm(n_sims, mean = 13000, sd = 200), # Cost of treatment B

  u_H           = rbeta(n_sims, 10, 1),                  # Utility close to 1 for Healthy
  u_S1          = rbeta(n_sims, 7.5, 2.5),               # Utility less than Healthy, beta distribution
  u_S2          = rbeta(n_sims, 5, 5),                   # Utility for Sicker
  u_D           = 0,                                        # Utility for Dead is constant
  u_trtA        = rbeta(n_sims, 9.5, 1),                 # Utility with treatment A, close to Healthy

  du_HS1        = rnorm(n_sims, mean = 0.01, sd = 0.005), # Disutility with slight variation
  ic_HS1        = rnorm(n_sims, mean = 1000, sd = 100),   # Cost increase with transition
  ic_D          = rnorm(n_sims, mean = 2000, sd = 100),    # Cost increase when dying
  p0_H          = rbeta(n_sims, 1, 9)                   # Initial probability of being Healthy
)

``` 


# Probability and Payoff Functions
Then, we define the probability and payoff functions used in the `twig` above:

## Probability of recovery (pRecover)
Below, we define the probability of recovery `pRecover` which returns the probability of recovery based on the state and the recovery rate `r_S1H`. However, it is important to first explain logical conditions in R which we will use as shorthand conditions replacing `ifelse` statements in the vectorized function definitions that we use in `twig`.

### Using Logical Conditions in R 

In R, a common way to assign values based on conditions is to use logical conditions. For example, the following function `pRecover` assigns a recovery rate `r_S1H` if the state is "S1" and 0 otherwise using a logical condition: 

``` {r}
pRecover <- function(state, r_S1H){
  rRecover <- r_S1H * (state=="S1") # Assigns r_S1H if in S1, otherwise 0
  rate2prob(rRecover)
}
```

The shorthand notation  `(state == "S1")` is a concise and efficient way to create a logical vector that evaluates to 1 (`TRUE`) when state is "S1" and 0 (`FALSE`) otherwise. This is particularly useful and easier to read than using ifelse statements which can become difficult to read especially when mutliple conditions are nested.  In this example an equivalent function using `ifelse` statement will look like this 

``` {r}
pRecover_ifelse <- function(state, r_S1H){
  rRecover <- ifelse(state == "S1", r_S1H, 0)  # Assigns r_S1H if in S1, otherwise 0
  rate2prob(rRecover)
}
``` 

Although, both of these functions are equivalent in `twig`, the shorthand notation using conditional logic is more concise and easier to read and we will use it throughout this example and the rest of the vignettes.

## Probability of getting sick (pGetSick)
Only those who are H can get sick
``` {r}
pGetSick <- function(state, r_HS1){
  rGetSick <- r_HS1 * (state=="H")
  rate2prob(rGetSick)
}
```

##  Probability of progressing (pProgress)
This depends on the state (only those who are in S1 can progress), the decision (different rates for different decisions), and the cycle_in_state (number of cycles spent in a state - tunnel state)
``` {r}
pProgress <- function(state, decision, cycle_in_state,
                      hr_S1S2_trtB, r_S1S2_scale, r_S1S2_shape){
  
  hr_S1S2 <- hr_S1S2_trtB ^ (decision %in% c("StrategyB", "StrategyAB")) # hazard rate of progression for B or 1 otherwise

  r_S1S2_tunnels <- ((cycle_in_state*r_S1S2_scale)^r_S1S2_shape - 
    ((cycle_in_state - 1)*r_S1S2_scale)^r_S1S2_shape) # hazard rate based on cycle_in_state (tunnel) which follows a weibull distribution

  # only those who are at S1 can progress
  rProgress <- r_S1S2_tunnels * (state=="S1") * hr_S1S2 

  rate2prob(rProgress)
}
```


##  Probability of dying (pDie)
Probabilty of dying depends on age. So, we define age-specific mortality data from the 2015 US Life Tables for ages 24 through 1000. Then, we define the `pDie` as a function of the state (different states have different rates of death), and the cycles since the simulation start which reflects the cohort's age.
``` {r}
v_r_mort_by_age <- vector <- c(
  0.000979, 0.001014, 0.000999, 0.001070, 0.001087, 0.001162, 0.001167, 0.001213, 0.001289, 
  0.001331, 0.001375, 0.001420, 0.001490, 0.001550, 0.001616, 0.001657, 0.001747, 0.001902, 
  0.002052, 0.002173, 0.002395, 0.002559, 0.002807, 0.003023, 0.003349, 0.003712, 0.004085, 
  0.004490, 0.004905, 0.005364, 0.005806, 0.006253, 0.006775, 0.007395, 0.007895, 0.008418, 
  0.008974, 0.009666, 0.010456, 0.011384, 0.011838, 0.012667, 0.013593, 0.014700, 0.015732, 
  0.017340, 0.018758, 0.020967, 0.022917, 0.024913, 0.026767, 0.029707, 0.032412, 0.035982, 
  0.039238, 0.043595, 0.048727, 0.053735, 0.059911, 0.066618, 0.074051, 0.082190, 0.090754, 
  0.103968, 0.115093, 0.124341, 0.137872, 0.154177, 0.172393, 0.194100, 0.212654, 0.243752, 
  0.259087, 0.287781, 0.316429, 0.339149
)

# death depends on the state and age.
pDie <- function(state, cycle,
                 hr_S1, hr_S2){
  r_HD <- v_r_mort_by_age[cycle] # get age-specific mortality 
  rDie <- r_HD * (state=="H") +  # baseline mortality if healthy
          r_HD*hr_S1 * (state=="S1") +  # multiplied by a hazard rate if S1 or 
          r_HD*hr_S2 * (state=="S2") # S2
          # else 0
  rate2prob(rDie)
}
```

## Cost
Cost is a function of the state, decision and whether either event have occured. The events capture transition costs. 
``` {r}
cost <- function(state, decision, second_event, death_event, 
                 ic_HS1, ic_D, c_trtA, c_trtB, 
                 c_H, c_S1, c_S2, c_D){
  # cost of decision is only applied if the state is either S1 or S2
  trans_cost_getting_sick <- ic_HS1 * (second_event=="getsick") # increase in cost when transitioning from Healthy to Sick
  trans_cost_dying <- ic_D * (death_event=="yes") # increase in cost when dying
  
  c_decision <- (state %in% c("S1","S2")) * (
      c_trtA * (decision=="StrategyA") +
      c_trtB * (decision=="StrategyB") +
      (c_trtA + c_trtB) * (decision=="StrategyAB") 
  )
  
  # cost of the state is a function of the state
  c_state <- c_H * (state=="H") + 
             c_S1 * (state=="S1") + 
             c_S2 * (state=="S2") + 
             c_D * (state=="D") 
  # combine both
  return(c_decision + c_state + trans_cost_getting_sick + trans_cost_dying)
}
```

## Utility
Similarly, utility depends on the state, decision and only if the cohort gotsick to apply a transition utility discount for those who make that transition.

``` {r}
utility <- function(state, decision, second_event,
                    du_HS1, u_H, u_trtA, u_S1, u_S2, u_D){

  trans_util_getting_sick <- -du_HS1 * (second_event=="getsick") # apply a utility discount for those who get sick.

  # calcualte state utilities. note that S1 will have utility u_trtA if the decision involves A, and another utility if the decision does not involve A.
  u_state <- u_H * (state=="H") + 
             u_trtA * (state=="S1" & decision %in% c("StrategyA", "StrategyAB")) +
             u_S1 * (state=="S1" & decision %out% c("StrategyA", "StrategyAB")) + 
             u_S2 * (state=="S2") + 
             u_D * (state=="D") 

  # combine the two utilities.
  return(u_state + trans_util_getting_sick)
}

```

# Running the `twig`
We can run the model and check the summary results under `results$mean_ev` of average costs and utilities across all probabilistic analyses.  The results also returns the individual simulation results under `results$sim_ev`

``` {r}
results <- run_twig(twig_obj = mytwig, params = params, n_cycles = n_cycles, progress_bar = FALSE)

results$mean_ev
```

## Parallelization
You can run the same command with parallelization, by setting `parallel = TRUE`. This will speed up the simulation by running each simulation in parallel. [It is commented out here because testing of the `twig` package only runs on a single core.]

``` {r}
# results <- run_twig(twig_obj = mytwig, params = params, n_cycles = n_cycles, parallel = TRUE)

```

## Incremental Cost-Effectiveness Ratio (ICER)
using the `calculate_icers` function adapted from the dampack package, we can retrieve ICER table by passing the payoffs summary table.
``` {R}
calculate_icers(results$mean_ev)
```

## Cost-Effectiveness Acceptibility Curve (CEAC)
Similarly, we can produce the CEAC curve from the simulation results and specifying a series of willingness to pay (WTP) values.

```{r fig.width=7, fig.height=5}
plot_ceac(results$sim_ev, wtp_range = seq(0, 100000, by = 1000))
```

## Summary
This example illustrated the following features of `twig`
* Markov model
* probabilistic anslysis 
* Multiple decisions
* Multiple states
* Multiple sequential events
* tunnel states 
* cycle-dependency
* dependency on prior events in a cycle
* transition payoffs
* using parallel computations
* expected payoffs across probabilistic analyses
* ICER
* CEAC

In addition, various intermediate objects and computations can be returned by enabling `verbose`. To reduce the size of the returned results, this will only use the first PSA sample. For details on these objects check the documentations.

