---
title: "Reproducing Exact Results from DMHEE"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Reproducing Exact Results from DMHEE}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r, echo=FALSE, include=FALSE}
library(heemod)
```

The purpose of this vignette is to present how to reproduce exactly the results from [Decision Modelling for Health Economic Evaluation](https://global.oup.com/academic/product/decision-modelling-for-health-economic-evaluation-9780198526629). While other vignettes such as `vignette("c-homogeneous", "heemod")`, `vignette("d-non-homogeneous", "heemod")` or `vignette("e-probabilistic", "heemod")` are greatly inspired from this book, key elements differ (mostly for the sake of clarity) and thus results differ, sometimes significantly, from the book. Here we show how to exactly reproduce the results with the `heemod` package.

# HIV model

Key differences in DMHEE:

  1. transitions occur at the end of each year,
  2. cost are counted starting from year 1, not year 0,
  3. treatment stops after 2 years,
  4. rounding errors.


It is possible to reproduce 1. and 2. by making transition happen at the end of each year with method = "end". Since with this method the transition occur 1 year after the beginning, costs should be discounted from the first cycle with the argument first = TRUE in discount().
Point 3. is reproduced by making `rr` and `cost_lami` a time changing variable like this `rr = ifelse(model_time <= 2, .509, 1.00)`.

The last point is reproduced by writing transition probabilities as fractions.

```{r}
par_mod <- define_parameters(
  rr = ifelse(model_time <= 2, .509, 1),
  cost_lami = ifelse(model_time <= 2, 2086.5, 0),
  cost_zido = 2278
)

mat_mono <- define_transition(
  1251/1734, 350/1734, 116/1734,  17/1734,
  0,         731/1258, 512/1258,  15/1258,
  0,         0,        1312/1749, 437/1749,
  0,         0,        0,         1.00
)

mat_comb <- define_transition(
  C, 350/1734*rr, 116/1734*rr, 17/1734*rr,
  0, C,           512/1258*rr, 15/1258*rr,
  0, 0,           C,           437/1749*rr,
  0, 0,           0,           1.00
)

A_mono <- define_state(
  cost_health = 2756,
  cost_drugs = cost_zido,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 1
)
B_mono <- define_state(
  cost_health = 3052,
  cost_drugs = cost_zido,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 1
)
C_mono <- define_state(
  cost_health = 9007,
  cost_drugs = cost_zido,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 1
)
D_mono <- define_state(
  cost_health = 0,
  cost_drugs = 0,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 0
)

A_comb <- define_state(
  cost_health = 2756,
  cost_drugs = cost_zido + cost_lami,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 1
)
B_comb <- define_state(
  cost_health = 3052,
  cost_drugs = cost_zido + cost_lami,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 1
)
C_comb <- define_state(
  cost_health = 9007,
  cost_drugs = cost_zido + cost_lami,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 1
)
D_comb <- define_state(
  cost_health = 0,
  cost_drugs = 0,
  cost_total = discount(
    cost_health + cost_drugs, .06, first = T),
  life_year = 0
)

mod_mono <- define_strategy(
  transition = mat_mono,
  A_mono,
  B_mono,
  C_mono,
  D_mono
)
mod_comb <- define_strategy(
  transition = mat_comb,
  A_comb,
  B_comb,
  C_comb,
  D_comb
)

res_mod <- run_model(
  mono = mod_mono,
  comb = mod_comb,
  parameters = par_mod,
  cycles = 20,
  cost = cost_total,
  effect = life_year,
  method = "end",
  init = c(1, 0, 0, 0)
)
summary(res_mod)
```

# THR model

Key difference in DMHEE:

  1. Mortality rates are much higher in the book.
  
This can be corrected by using a user-specified mortality table and then fetch the values with:

```r
look_up(death_prob, age = age, sex = sex, bin = "age")
```

```{r}
# a function to return age-related mortality rate
# given age and sex
death_prob <- data.frame(
  age = rep(seq(35, 85, 10), each = 2),
  sex = rep(1:0, 6),
  value = c(
    1.51e-3, .99e-3, 3.93e-3,
    2.6e-3, 10.9e-3, 6.7e-3,
    31.6e-3, 19.3e-3, 80.1e-3,
    53.5e-3, 187.9e-3, 154.8e-3
  )
)
death_prob

param <- define_parameters(
    age_init = 60,
    sex = 0,
    # age increases with cycles
    age = age_init + model_time,
    
    # operative mortality rates
    omrPTHR = .02,
    omrRTHR = .02,
    
    # re-revision mortality rate
    rrr = .04,
    
    # parameters for calculating primary revision rate
    cons = -5.490935,
    ageC = -.0367022,
    maleC = .768536,
    lambda = exp(cons + ageC * age_init + maleC * sex),
    lngamma = 0.3740968,
    gamma = exp(lngamma),
    lnrrNP1 = -1.344474,
    rrNP1 = exp(lnrrNP1),
   
    # revision probability of primary procedure
    standardRR = 1 - exp(lambda * ((model_time - 1) ^ gamma -
                                     model_time ^ gamma)),
    np1RR = 1 - exp(lambda * rrNP1 * ((model_time - 1) ^ gamma - 
                                        model_time ^ gamma)),
    
    # age-related mortality rate
    sex_cat = ifelse(sex == 0, "FMLE", "MLE"),
    mr = look_up(death_prob, age = age, sex = sex, bin = "age"),
    
    u_successP = .85,
    u_revisionTHR = .30,
    u_successR = .75,
    c_revisionTHR = 5294
)

mat_standard <- define_transition(
    state_names = c(
      "PrimaryTHR",
      "SuccessP",
      "RevisionTHR",
      "SuccessR",
      "Death"
    ),
    0, C, 0,          0, omrPTHR,
    0, C, standardRR, 0, mr,
    0, 0, 0,          C, omrRTHR+mr,
    0, 0, rrr,        C, mr,
    0, 0, 0,          0, 1
)

mat_np1 <- define_transition(
    state_names = c(
      "PrimaryTHR",
      "SuccessP",
      "RevisionTHR",
      "SuccessR",
      "Death"
    ),
    0, C, 0,     0, omrPTHR,
    0, C, np1RR, 0, mr,
    0, 0, 0,     C, omrRTHR+mr,
    0, 0, rrr,   C, mr,
    0, 0, 0,     0, 1
)

mod_standard <- define_strategy(
  transition = mat_standard,
  PrimaryTHR = define_state(
    utility = 0,
    cost = 394
  ),
  SuccessP = define_state(
    utility = discount(u_successP, .015),
    cost = 0
  ),
  RevisionTHR = define_state(
    utility = discount(u_revisionTHR, .015),
    cost = discount(c_revisionTHR, .06)
  ),
  SuccessR = define_state(
    utility = discount(u_successR, .015),
    cost = 0
  ),
  Death = define_state(
    utility = 0,
    cost = 0
  )
)

mod_np1 <- define_strategy(
  transition = mat_np1,
  PrimaryTHR = define_state(
    utility = 0,
    cost = 579
  ),
  SuccessP = define_state(
    utility = discount(u_successP, .015),
    cost = 0
  ),
  RevisionTHR = define_state(
    utility = discount(u_revisionTHR, .015),
    cost = discount(c_revisionTHR, .06)
  ),
  SuccessR = define_state(
    utility = discount(u_successR, .015),
    cost = 0
  ),
  Death = define_state(
    utility = 0,
    cost = 0
  )
)

res_mod <- run_model(
  standard = mod_standard,
  np1 = mod_np1,
  parameters = param,
  cycles = 60,
  cost = cost,
  effect = utility,
  method = "beginning",
  init = c(1, 0, 0, 0, 0)
)
summary(res_mod)
```