---
title: 'Markov Sick-Sicker model in R - Excerise'
subtitle: "With simulation-time dependence"
author: "The DARTH workgroup"
output:
  html_document: default
  pdf_document: default
keep_tex: yes
self_contained: no
---

Developed by the Decision Analysis in R for Technologies in Health (DARTH) workgroup.
If you use this code as a basis for your own work, please cite our papers (see suggestions below) and in your acknowledgments for both the code and publications please add:
**The authors thank the Decision Analysis in R for Technologies in Health (DARTH) workgroup for providing the coding framework in R.**


This code implements a simulation-time-dependent Sick-Sicker cSTM model to  conduct a CEA of two strategies:
- Standard of Care (SoC): best available care for the patients with the disease. This scenario reflects the natural history of the disease progression.
- Strategy AB: This strategy combines treatment A and treatment B. The strategy treats both those Sick and Sicker. However, only for those Sick the treatment has an effect. For Sick individuals the disease progression is reduced, and individuals in the Sick state have an improved quality of life.

\newpage

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, eval = TRUE)
options(scipen = 999)  # disable scientific notation
rm(list = ls())        # clear memory (removes all the variables from the workspace)
```

Change `eval` to `TRUE` if you want to knit this document.

# Exercise: Construct a Markov Model of the Sick-Sicker Disease

In this exercise, we will model a hypothetical disease that affects individuals with an average age of 25 years and results in increased mortality, increased healthcare costs, and reduced quality of life. The disease has two levels; affected individuals initially become sick but can subsequently progress and become sicker. Two alternative strategies exist for this hypothetical disease: Standard of Care (SoC) and a treatment strategy - Strategy AB. Under the treatment strategy, individuals in the sick and sicker states are treated until they recover (only if sick; individuals in the sicker state cannot recover) or die. The cost of the treatment is additive to the baseline healthcare costs of being sick or sicker. The treatment improves quality of life for those individuals who are sick but has no impact on the quality of life of those who are sicker. Unfortunately, it is not possible to reliably differentiate between people in the sick and sicker states, so treatment cannot be targeted to only those in the sick state. You are asked to evaluate the cost-effectiveness of the treatment.

To model this disease, we will rely on a state-transition cohort model, called the Sick-Sicker model, first described by Enns et al. The Sick-Sicker model consists of four health states: Healthy (H), two disease states, Sick (S1) and Sicker (S2), and Dead (D) (Figure 1). All individuals start in the Healthy state. Over time, healthy individuals may develop the disease and can progress to S1. Individuals in S1 can recover (return to state H), progress further to S2 or die. Individuals in S2 cannot recover (i.e. cannot transition to either S1 or H). Individuals in H have a baseline probability of death; individuals in S1 and S2 experience increased mortality compared to those in the H state, given in terms of hazard ratios. These ratios are used to calculate the probabilities of dying when in S1 and S2. 

```{r, echo = F, warning = F, message = F, out.width = '100%', fig.cap = 'Schematic representation of the Sick-Sicker model', fig.align = 'center'}
if (!require(knitr)) install.packages('knitr')
# include_graphics("../figures/sick_sicker_diagram.png")
```

\newpage

**Table I: Input parameters**

|           **Parameter**            |  **R name** |   **Value**   |
|:-----------------------------------|:------------|:-------------:|
| Cycle length                       | `cycle_length`          | 1  year           |
| Age at baseline                    | `n_age_init`            | 25 years old      |
| Maximum age of follow-up           | `n_age_max`             | 100 years old     |
| Names of health states             | `v_names_states`        | H, S1, S2, D      |
| Names of cycles (time horizon)     | `n_cycles`              | (n_age_max - n_age_init) / cycle_length |
| Annual discount rate (costs/QALYs) | `d_c`/ `d_e`       | 3%                      |
| - Rate of becoming S1 when H       | `r_HS1`     |  0.15          |
| - Rate of becoming H when S1       | `r_S1H`     |  0.5           |
| - Rate of becoming S2 when S1      | `r_S1S2`    |  0.105         |
| Annual mortality                   |             |                |
| - All-cause mortality (H to D)     | `v_r_HD`      | `HMD - info below`   |
| - Hazard ratio of death in S1 vs H | `hr_S1`       |  3          |
| - Hazard ratio of death in S2 vs H | `hr_S2`       |  10         |
| - Hazard ratio of becoming Sicker when Sick under Strategy AB  | `hr_S1S2_trtAB`       |  0.6         |
| Annual costs                       |             |               |
| - Healthy individuals              | `c_H`       |  $2,000       |
| - Sick individuals in S1           | `c_S1`      |  $4,000       |
| - Sick individuals in S2           | `c_S2`      |  $15,000      |
| - Dead individuals                 | `c_D`       |  $0           |
| - Additional costs of sick individuals treated with Strategy AB in S1 or S2       | `c_trtAB` | $25,000 |
| Utility weights                    |             |               |
| - Healthy individuals              | `u_H`       |  1.00         |
| - Sick individuals in S1           | `u_S1`      |  0.75         |
| - Sick individuals in S2           | `u_S2`      |  0.50         |
| - Dead individuals                 | `u_D`       |  0.00         |
| - Utility for individuals treated with Strategy AB in S1 | `u_trtAB` |  0.95        |

*Note: To calculate the probability of dying from S1 and S2, use the hazard ratios provided. To do so, first multiply the rate of dying from healthy by the appropriate hazard ratio; finally, convert this rate back to a probability. Recall that you can convert between rates and probabilities using the following formulas: $r = -log(1-p)$ and $p = 1-e^{(-rt)}$. The package `darthtools` also has the functions `prob_to_rate` and `rate_to_prob` that might be of use to you.

* 1. *HMD* : This is the Human mortality database. We provide you a files called `HMD_USA_Mx.csv` with the age specific mortality rates. The data has the columns X, Year, Age, Female, Male, Total and OpenInterval. We will use the overall age specific mortality rate from the Total column.

### Exercise - set-up

There are quite some steps you need to take in order to create a Markov model reflecting this case. Start by making sure that you have all the course material structured in the suggested way. There should be an R projects in the folder, a file ending on `.Rproj`, double click on this file to open the Rstudio environment. Working from a Rproject avoid struggles with working directories. 

While in you Markov Sick-Sicker R project, open this `markov_sick-sicker_time_exercise_template.Rmd`. This is a template that guides you step by step through the exercise. In addition it helps you to load the packages and the necessary data. The next chunk load the packages. 

# 01 Load packages

```{r, warning = F, message = F}
if (!require('pacman')) install.packages('pacman'); library(pacman) # use this package to conveniently install other packages
# load (install if required) packages from CRAN
p_load("dplyr", "tidyr", "reshape2", "devtools", "scales", "ellipse", "ggplot2", "ggrepel", "gridExtra", "lazyeval", "igraph", "truncnorm", "ggraph", "reshape2", "patchwork", "knitr", "stringr", "diagram", "dampack")                                               
# load (install if required) packages from GitHub
# install_github("DARTH-git/darthtools", force = TRUE) #Uncomment if there is a newer version
p_load_gh("DARTH-git/darthtools")
```

# 02 Load functions

```{r}
# all functions are in the darthtools package
```

# 03 Model input

```{r}
## General setup 
cycle_length <- 1/12   # cycle length equal to one year (use 1/12 for monthly)
n_age_init   <- 25  # age at baseline
n_age_max    <- 100 # maximum age of follow up
n_cycles     <- (n_age_max - n_age_init)/cycle_length # time horizon, number of cycles
# Age labels 
v_age_names  <- paste(rep(n_age_init:(n_age_max - 1), each = 1/cycle_length), 
                      1:(1/cycle_length), 
                      sep = ".")
# the 4 health states of the model:
v_names_states <- c("H",  # Healthy (H)
                    "S1", # Sick (S1)
                    "S2", # Sicker (S2)
                    "D")  # Dead (D)
                                           
n_states <- length(v_names_states)   # number of health states 

### Discounting factors 
# Comment out/ for the relevant country
# USA
#d_c <- d_e <- 0.03  # annual discount rate for costs and QALY
# Canada
#d_c <- d_e <- 0.015 # annual discount rate for costs 

# NL
d_c <- 0.03         # annual discount rate for costs 
d_e <- 0.015        # annual discount rate for QALYs


### Strategies 
v_names_str <- c("Standard of Care", # store the strategy names
                 "Strategy AB") 
n_str       <- length(v_names_str)   # number of strategies

## Within-cycle correction (WCC) using Simpson's 1/3 rule 
v_wcc  <- darthtools::gen_wcc(n_cycles = n_cycles, method = "Simpson1/3")

### Transition rates (annual), and hazard ratios (HRs) 
r_HS1  <- 0.15  # constant annual rate of becoming Sick when Healthy
r_S1H  <- 0.5   # constant annual rate of becoming Healthy when Sick
r_S1S2 <- 0.105 # constant annual rate of becoming Sicker when Sick
hr_S1  <- 3     # hazard ratio of death in Sick vs Healthy 
hr_S2  <- 10    # hazard ratio of death in Sicker vs Healthy 

### Effectiveness of treatment AB 
hr_S1S2_trtAB <- 0.6  # hazard ratio of becoming Sicker when Sick under treatment AB

## Age-dependent mortality rates 
lt_usa <- read.csv("HMD_USA_Mx.csv")
# Extract age-specific all-cause mortality for ages in model time horizon
v_r_mort_by_age <- lt_usa %>% 
  dplyr::filter(Age >= n_age_init & Age < n_age_max) %>%
  dplyr::select(Total) %>%
  as.matrix()

### State rewards 
#### Costs 
c_H     <- 2000  # annual cost of being Healthy
c_S1    <- 4000  # annual cost of being Sick
c_S2    <- 15000 # annual cost of being Sicker
c_D     <- 0     # annual cost of being dead
c_trtAB <- 25000 # annual cost of receiving treatment AB
#### Utilities 
u_H     <- 1     # annual utility of being Healthy
u_S1    <- 0.75  # annual utility of being Sick
u_S2    <- 0.5   # annual utility of being Sicker
u_D     <- 0     # annual utility of being dead
u_trtAB <- 0.95  # annual utility when receiving treatment AB

### Discount weight for costs and effects 
v_dwc   <- 1 / ((1 + (d_c * cycle_length)) ^ (0:n_cycles))
v_dwe   <- 1 / ((1 + (d_e * cycle_length)) ^ (0:n_cycles))

# Process model inputs 
## Age-specific transition rates to the Dead state for all cycles 
v_r_HD_age  <- rep(v_r_mort_by_age, each = 1/cycle_length)
# Name age-specific mortality vector 
names(v_r_HD_age) <- v_age_names

# compute mortality rates
v_r_S1D_age  <- v_r_HD_age * hr_S1 # Age-specific mortality rate in the Sick state 
v_r_S2D_age  <- v_r_HD_age * hr_S2 # Age-specific mortality rate in the Sicker state 
# transform rates to probabilities adjusting by cycle length
p_HS1       <- rate_to_prob(r = r_HS1,  time = cycle_length) # constant annual probability of becoming Sick when Healthy 
p_S1H       <- rate_to_prob(r = r_S1H,  time = cycle_length) # constant annual probability of becoming Healthy when Sick 
p_S1S2      <- rate_to_prob(r = r_S1S2, time = cycle_length) # constant annual probability of becoming Sicker when Sick 
v_p_HD_age   <- rate_to_prob(v_r_HD_age,  time = cycle_length) # Age-specific mortality risk in the Healthy state 
v_p_S1D_age  <- rate_to_prob(v_r_S1D_age, time = cycle_length) # Age-specific mortality risk in the Sick state

v_p_S2D_age  <- rate_to_prob(v_r_S2D_age, time = cycle_length) # Age-specific mortality risk in the Sicker state


## Annual transition probability of becoming Sicker when Sick for treatment AB 
# Apply hazard ratio to rate to obtain transition rate of becoming Sicker when Sick for treatment AB
r_S1S2_trtAB <- r_S1S2 * hr_S1S2_trtAB
# Transform rate to probability to become Sicker when Sick under treatment AB 
# adjusting by cycle length 
p_S1S2_trtAB <- rate_to_prob(r = r_S1S2_trtAB, time = cycle_length)
```

### Exercise - Task 1

Build the Markov model in `R` for Standard of Care (SoC) and Strategy AB including age-dependent mortality and do the following:
(1) Run the model
(2) Plot the cohort trace from the model
(3) Compute state rewards and expected outcomes (total utilities and costs)

# 04 Construct state-transition models

## 04.1 Initial state vector

```{r}
# All starting healthy
v_m_init <- c(H = 1, S1 = 0, S2 = 0, D = 0) # initial state vector
v_m_init
```

## 04.2 Initialize cohort traces

```{r}
### Initialize cohort trace under SoC 
m_M_SoC <- matrix(NA, 
              nrow = (n_cycles + 1), ncol = n_states, 
              dimnames = list(0:n_cycles, v_names_states))
# Store the initial state vector in the first row of the cohort trace
m_M_SoC[1, ] <- v_m_init

### Initialize cohort trace for strategy AB 
# Structure and initial states are the same as for SoC
m_M_strAB <- m_M_SoC # Strategy AB
```

## 04.3 Create transition probability array

```{r, eval = TRUE}
## Create transition probability arrays for strategy SoC 
### Initialize transition probability array for strategy SoC 
a_P_SoC <- array(0,
                 dim  = c(n_states, n_states, n_cycles),
                 dimnames = list(v_names_states, 
                                 v_names_states, 
                                 0:(n_cycles - 1)))
### Fill in array
## From H
# your turn

## From S1
# your turn

## From S2
# your turn

## From D
# your turn


### Initialize transition probability array for strategy AB 
# your turn

# Update only transition probabilities from S1 involving p_S1S2
# your turn

## Check if transition probability arrays are valid 
### Check that transition probabilities are [0, 1] 
check_transition_probability(a_P_SoC,   verbose = TRUE)
check_transition_probability(a_P_trtAB, verbose = TRUE)


### Check that all rows for each slice of the array sum to 1 
check_sum_of_transition_array(a_P_SoC,   n_states = n_states, n_cycles = n_cycles, verbose = TRUE)
check_sum_of_transition_array(a_P_trtAB, n_states = n_states, n_cycles = n_cycles, verbose = TRUE)
```

# 05 Run Markov model

```{r}
# Iterative solution of age-dependent cSTM
# your turn

## Store the cohort traces in a list 
# your turn
```

# 06 Plot Outputs

## 06.1 Plot the cohort trace for strategies SoC and AB

```{r}
# your turn
```
ADD AN DESCRIPTION OF WHAT YOU SEE AND IF THAT MAKES SENSE TO YOU 

# 07 State Rewards 

```{r}
## Scale by the cycle length 
# Vector of state utilities under strategy SoC
# your turn

# Vector of state costs under strategy SoC
# your turn

# Vector of state utilities under strategy AB
# your turn

# Vector of state costs under strategy AB
# your turn

## Store state rewards 
# Store the vectors of state utilities for each strategy in a list 
# your turn

# Store the vectors of state cost for each strategy in a list 
# your turn

# assign strategy names to matching items in the lists
# your turn

```

# 08 Compute expected outcomes 

```{r}
# Create empty vectors to store total utilities and costs 
# your turn

## Loop through each strategy and calculate total utilities and costs 
# your turn

```
ADD AN INTERPRETATION OF YOUR RESULTS

### Exercise - Task 2

Estimate the cost-effectiveness of Strategy AB vs SoC.

# 09 Cost-effectiveness analysis (CEA) 

```{r}
## Incremental cost-effectiveness ratios (ICERs) 
# your turn

```
ADD AN INTERPRETATION OF YOUR RESULTS

### Exercise - Task 3

Create a well-formatted cost-effectiveness table with all results of interest and plot the cost-effectiveness frontier.

```{r}
## CEA table in proper format 
# your turn

```

ADD AN INTERPRETATION OF YOUR RESULTS

```{r}
## CEA frontier 
# your turn
```
ADD AN INTERPRETATION OF YOUR RESULTS

We kindly request you to add the following Acknowledgement paragraph to your further work where DARTH code formed the basis. We also like to remind you that you can add other sources of reference to this paragraph to acknowledge code you got from others. 

# Acknowlegdement

For this work we made use of the template developed by the Decision Analysis in R for Technologies in Health (DARTH) workgroup: <http://darthworkgroup.com>.

The notation of our code is based on the following provided framework and coding convention: Alarid-Escudero, F., Krijkamp, E., Pechlivanoglou, P. et al. A Need for Change! A Coding Framework for Improving Transparency in Decision Modeling. PharmacoEconomics 37, 1329–1339 (2019). <https://doi.org/10.1007/s40273-019-00837-x>.

Other work from DARTH can be found on the website: <http://darthworkgroup.com/publications/>

# Copyright for assignment work

Copyright 2017, THE HOSPITAL FOR SICK CHILDREN AND THE COLLABORATING INSTITUTIONS.All rights reserved in Canada, the United States and worldwide. Copyright, trademarks, trade names and any and all associated intellectual property are exclusively owned by THE HOSPITAL FOR Sick CHILDREN and the collaborating  institutions. These materials may be used, reproduced, modified, distributed and adapted with proper attribution.
