---
title: "Simple sick sicker model"
author: "Sheeja Manchira Krishnan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simple sick sicker model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ref.bib
---

This document demonstrates developing For a simple markov model explained in 
reference \ref[alarid2020cohort] using packDAMipd. This model is explained in 
STM_01.R file in Cohort-modeling-tutorial provided in DARTH-git github page. 

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = ">", fig.width = 6, fig.height = 4
)
```

```{r setup}
library(packDAMipd)
```
Define health states. Each health state must have a name, cost, and utility. 
The cost and utility are the state values. That means the individuals having 
the state membership will occur the particular cost and utility. If any 
cost/utility due to transitions occur, they need to specify similar to defining 
the transition probabilities.

```{r}
H <- health_state("H", cost = "c_H ", utility = "u_H")
S1 <- health_state("S1", cost = "c_S1", utility = "u_S1")
S2 <- health_state("S2", cost = "c_S2", utility = "u_S2")
D <- health_state("D", cost = "c_D ", utility = "u_D")
```

Define allowed transition probabilities and number them. The below matrix is 
numbered so that the maximum entry in the matrix gives the total number of 
allowed transitions. Column names and row names are just the names of the 
health states."NA" indicates a non allowed transition. 
```{r}
tmat <- rbind(c(1, 2, NA, 3), c(4, 5, 6, 7),c(NA, NA, 8, 9), c(NA, NA, NA, 10))
colnames(tmat) <- rownames(tmat) <- c("H","S1" ,"S2","D")
```

As there is cost and utility during some transitions, we define tmat_cost and tmat_utility similar to defining transition probability. Again only non zero transition cost and transition utility are non NA.

```{r}
tmat_cost <- rbind(c(NA, 1, NA, 2), c(NA, NA, NA, 3), c(NA, NA, NA, 4), 
                   c(NA, NA, NA, NA))
colnames(tmat_cost) <- rownames(tmat_cost) <- c("H", "S1", "S2", "D")
tmat_util <- rbind(c(NA, 1, NA, NA), c(NA, NA, NA, NA), c(NA, NA, NA, NA), 
                   c(NA, NA, NA, NA))
colnames(tmat_util) <- rownames(tmat_util) <- c("H", "S1", "S2", "D")

```

Now we define the transition matrix , transition cost matrix and transition 
utility matrix (the latter two are optional,we use only if there are 
cost/utility during transitions) using the tmat defined earlier, and number of
health states.

```{r}
tm <- populate_transition_matrix(4, tmat, c("p_HH","p_HS1","p_HD","p_S1H",
                                   "p_S1S1", "p_S1S2", "p_S1D",
                                   "p_S2S2","p_S2D","p_DD" ), colnames(tmat))
tm_cost <- transition_cost_util(4, tmat_cost, c("ic_HS1","ic_D","ic_D","ic_D"), colnames(tmat_cost))
tm_util <- transition_cost_util(4, tmat_util, c("du_HS1" ), colnames(tmat_util))
```
Now define the parameter list using define_parameters(). If there is any formula 
or expression to be evaluated during the cycle, it has be given as a string (see below). For example. p_S1D is calculated using a formula with parameters p_HD 
and hr_S1 and is thus given as a string expression.
```{r}
param_list <- define_parameters(p_HD = 0.002, p_HS1  = 0.15, p_S1H = 0.5,
                     p_S1S2  = 0.105,hr_S1   = 3,hr_S2   = 10,
                     p_S1D  = "1 - exp(log(1 - p_HD) * hr_S1)",
                     p_S2D  = "1 - exp(log(1 - p_HD) * hr_S2)",
                     p_HH = "1 - (p_HS1 + p_HD)",
                     p_S1S1 = "1 - (p_S1H + p_S1S2+ p_S1D)",
                     p_S2S2 = "1 - ( p_S2D)",
                     p_DD = 1,
                     c_H   = 2000,c_S1  = 4000,c_S2  = 15000,
                     c_D   = 0, c_Trt = 12000,u_H   = 1,
                     u_S1  = 0.75,u_S2  = 0.5,u_D   = 0,
                     u_Trt = 0.95, du_HS1 = -0.01,ic_HS1 = 1000,ic_D   = 2000 )
```

We define the strategy "Usual care" by using the function 'strategy' with the 
health states combined (using combine_state()), the transition matrix (tm), 
transition cost and utility matrices. Using the defined strategy 'uc_strategy', 
we run the markov model using markov_model() with cycles, initial states, 
initial cost, initial utility values, discount rates and parameter list. The 
results of the markov model (trace matrix) can be plotted using plot_model() 
function.
```{r}
health_states <- combine_state(H,S1,S2,D)
uc_strategy <- strategy(tm, health_states, "Usual care",tm_cost,tm_util)
uc_markov <- markov_model( uc_strategy, 85, c(1, 0,0,0), c(0.03,0.03), 
                           param_list, TRUE, method = "half cycle correction")
# plot_model(uc_markov)
```

Similarly, we define the health states, strategy, and markov model for the "New treatment". Parameter list is the same as that defined for the usual care. 
Costs for the health states in the new treatment strategy differs, thus we need 
to redefine the health states.


```{r}
H <- health_state("H", cost = "c_H ", utility = "u_H")
S1 <- health_state("S1", cost = "c_S1 + c_Trt ",utility = "u_Trt")
S2 <- health_state("S2", cost = "c_S2 + c_Trt",utility = "u_S2")
D <- health_state("D", cost = "c_D ",utility = "u_D")
health_states <- combine_state(H,S1,S2,D)
trt_strategy <- strategy(tm, health_states, "New treatment",tm_cost,tm_util)
trt_markov <- markov_model(current_strategy = trt_strategy, cycles = 85, 
                           initial_state = c(1, 0,0,0), 
                           discount = c(0.03,0.03), parameter_values = 
                             param_list, TRUE, method = 
                             "half cycle correction")

```

Once we have the results form the markov models for both the strategies, we 
combine these two to define list_markov. Finally we obtain incremental cost effectiveness ratio (ICER) and net monetary benefit (NMB) for the cost 
effectiveness analysis using the function calculate_icer_nmb(). The results 
of the cost effectiveness acceptability curve(CEAC) can be plotted using 
plot_nmb_lambda function as shown below. 
For calculating NMB, threshold value  and the name of the comparator should be provided, while for plotting CEAC, range of values of thresholds should be 
provided along with the name of the comparator.

The \ref[alarid2020cohort] reports the results as follows which is the same as obtained using packDAMipd.

The cost of usual care obtained from STM_01.R of  Cohort-modeling-tutorial in DARTH-git  is 156,703 USD while that of the new treatment is  287,912 USD.
The  QALYs for the usual care and the new treatment are 21.556  and 22.289 respectively.
The ICER obtained is 178,792.

```{r}
list_markov <- combine_markov(uc_markov, trt_markov)
icerandnamb = calculate_icer_nmb(list_markov, threshold = 20000, comparator = 
                                   "Usual care")
icerandnamb
```

```{r}
# plot_nmb_lambda(list_markov,threshold_values = c(1000, 2000, 5000, 7000, 10000, 150000,20000), comparator = "Usual care" , currency = "USD")
```
