---
title: "Replicating FRH"
output: html_notebook
editor_options: 
  chunk_output_type: console
bibliography: ../references.bib  
---

# Introduction

The objective of this notebook is to employ our approach to modeling disability-adjusted life years (DALYs) in discrete time Markov cohort models to replicate the examples in @rushby2001 (FRH) and @larson2013calculating. 

```{r setup, echo=TRUE, message=FALSE, warning=FALSE}
#| echo: true
#| message: false
#| warning: false
library(tidyverse)
library(MASS)
library(expm)
library(knitr)
library(kableExtra)
options(scipen = 5) 
transpose <- purrr::transpose
select <- dplyr::select
options(knitr.kable.NA = '')

gen_wcc <- function (n_cycles, method = c("Simpson1/3", "half-cycle", "none")) 
{
    if (n_cycles <= 0) {
        stop("Number of cycles should be positive")
    }
    method <- match.arg(method)
    n_cycles <- as.integer(n_cycles)
    if (method == "Simpson1/3") {
        v_cycles <- seq(1, n_cycles + 1)
        v_wcc <- ((v_cycles%%2) == 0) * (2/3) + ((v_cycles%%2) != 
                                                     0) * (4/3)
        v_wcc[1] <- v_wcc[n_cycles + 1] <- 1/3
    }
    if (method == "half-cycle") {
        v_wcc <- rep(1, n_cycles + 1)
        v_wcc[1] <- v_wcc[n_cycles + 1] <- 0.5
    }
    if (method == "none") {
        v_wcc <- rep(1, n_cycles + 1)
    }
    return(v_wcc)
}
```

FRH consider a woman who develops bipolar depression at age 35, lives for 10 years with the disorder (disability weight = 0.60), and then dies prematurely at age 45. Remaining life expectancy at age 45 is 34.73 years. 

Based on this information, we parameterize a discrete time Markov model in which 35-year old individuals start of sick. The probability of death is then held to zero for 10 years, at which time the probability of death from disease is set to 1.0. We structure our model based on Approach 1, which defines a separate disease-related death transition. 

```{r}
data.frame(state = c("Healthy","Sick","Dead-Other Cause","Dead-Disease"), name = c("H","S1","DOC","DS")) %>% 
  kable(col.names=c("Health State","Health State Label")) %>% 
  kable_styling()
```

# Parameterize 
```{r parameterize}
params_ <- list(
    # Treatment Strategies
    v_tx_names = c("FRH"),      # treatment names
    n_tx = 1, # number of treatment strategies
    
    cycle_correction = "half-cycle",
    
    v_tr_names = c("H","S1"), # transient health states
    v_ab_names = c("DOC","DS"), # absorbing health states
    n_states = 4, # total number of health states
    
    horizon = 50,  # TK 400 # model time horizon (in years)  
    r_v_disc_h  = 0.03,     # annual discount rate for health outcomes
    r_v_disc_c = 0.03,     # annual discount rate for cost outcomes
    Delta_t = 1,      # time step (1 = yearly, 1/12 = monthly, etc.)
    age0 = 35,         # age at baseline
    v_s0T = c(0,1,0,0), # initial state occupancy  
                      # c(1,0,0,0,0) means the modeled cohort starts off healthy
    
    r_HS1 = 0,   # disease onset rate
    r_S1H = 0,    # recovery rate
    r_HD = 0,  # TK .002 # background mortality rate
    
    u_H = 1,           # Healthy utility weight
    u_S1 = 0.4,       # Sick utility weight
    u_D = 0,           # Death utility weight

    dw_S1 = 0.6,      # Sick disability weight
    
    df_ExR =  # Reference life table from GBD
          tibble::tribble(
              ~Age, ~Life.Expectancy,
              0L,      0,
              1L,      0,
              5L,      0,
              10L,     0,
              15L,     0,
              20L,     0,
              25L,     0,
              30L,     0,
              35L,     0,
              40L,     0,
              45L,    34.73 , #FRH2 
              50L,      0,
              55L,      10, # FRH1
              60L,      0,
              65L,      0,
              70L,      0,
              75L,      0,
              80L,      0,
              85L,      0,
              90L,      0,
              95L,      0
          )
)

params <- 
    with(params_,{
        modifyList(params_,list(
            v_names_states = c(v_tr_names, v_ab_names), # health state names
            omega = horizon/Delta_t,  # Total number of cycles
            r_v_disc_h_Delta_t = r_v_disc_h * Delta_t,  # Cycle discount rate: health outcomes
            r_v_disc_c_Delta_t = r_v_disc_c * Delta_t,  # Cycle discount rate: cost outcomes
            ages = (0:(horizon/Delta_t))*Delta_t + age0,  # Age in each cycle
             # Approximation function for reference life table life expectancies:
            f_ExR = function(x) x %>% map_dbl(~(if(.x==45)  34.73 else 0))
        ))
    })

params$ages_trace <- params$ages
params$ages <- params$ages[-length(params$ages)]
```

```{r transition-matrices}
params1 <- with(params,modifyList(params,list(
    # Natural History Transition Rate Matrix
    m_R = 
      ages %>% map(~({
        mR_FRH = 
          matrix(c(
          0,0,0,0, 
          0,0,0,0, 
          0,0,0,0,          
          0,0,0,0),
          nrow = n_states, 
          ncol = n_states,
          byrow=TRUE, 
          dimnames = list(c(v_tr_names,v_ab_names),
                          c(v_tr_names,v_ab_names)
          ))
        
        array(c(as.vector(mR_FRH)), 
              dim = c(length(v_tr_names)+ length(v_ab_names),length(v_tr_names)+ length(v_ab_names),length(v_tx_names)),
          dimnames = list(c(v_tr_names,v_ab_names),c(v_tr_names,v_ab_names),v_tx_names)) %>% 
            apply(.,3,function(x) x, simplify=FALSE) 
        
      }))
    )))

params1 <- with(params1,modifyList(params1,list(
  m_R_ = m_R %>% transpose()
)))
params1$m_R = params1$m_R_
params1[["m_R_"]] <- NULL
names(params1$m_R) = params1$v_tx_names

params1 <- with(params1,modifyList(params1,list(
  m_R_ = v_tx_names %>% map(~({
    p_ <- m_R[[.x]]
    a <- ages
    map2(a,p_,~({
      if (.x==35) {
        .y[1,2]  = 100
        .y[1,1] = -100
      }
      if (.x==44) {
        .y[2,4] = 100
        .y[2,2] = -100
      }
      if (.x>44) {
        .y[1,3] = 100
        .y[1,1] = -100
        .y[2,3] = 100
        .y[2,2] = -100
      }
      .y
    }))
  }))
)))
params1$m_R = params1$m_R_
params1[["m_R_"]] <- NULL
names(params1$m_R) = params1$v_tx_names

params1 <- with(params1,modifyList(params1,list(
    m_P = m_R  %>% map(~({
      mR_ = .x
      mR_ %>% map(~({
              expm(.x * Delta_t)
         }))
      }))
)))

params1 <- with(params1,modifyList(params1,list(
  m_P_ = m_P %>% map(~({
    tmp <- .x
    init <- diag(nrow=nrow(tmp[[1]]),ncol=ncol(tmp[[1]]))
    # FRH replication move
    #init[1,1] = 0
    #init[1,2] = 1
    dimnames(init) = dimnames(tmp[[1]])
    append(.x,list(init),after=0)
  }))
)))
params1$m_P = params1$m_P_
params1[["m_P_"]] <- NULL
names(params1$m_P) = params1$v_tx_names
```

# Outcomes

We then use the resulting Markov trace to calculate outcomes. We estimate years of life with disability (YLD), years of life lost to disease (YLL) and DALY outcomes using a 3% discount rate. This mirrors the FRH DALY(3, 0) approach in @rushby2001 and @larson2013calculating.

First, let's construct and view the Markov trace:

```{r trace}
trace1 <- 
    with(params1, {
        m_P %>% map( ~ ({
            P = .x
            occ <- v_s0T
            P %>% map(~({
              occ <<- occ %*% .x
            })) %>% 
            map(~(data.frame(.x))) %>% 
            bind_rows()
        }))
    })  
trace1 %>% as.data.frame() %>% 
  mutate(age = params1$ages_trace) %>% 
  select(age,everything()) %>% 
  kable() %>% kable_styling()

```

## Continuous Time Discounting Formula

Our first option for discounting is to use the continuous time discounting factor, $e^{-r \Delta_t t}$ where $r$ is the annual discount rate, $\Delta_t$ is the cycle length (a value of 1 indicates an annual cycle) and $t$ is the cycle in question. This option maintains the continuous time discounting appraoch adopted by the Global Burden of Disease study for DALY outcomes---however the formula itself is designed to accomodate a series of discrete payoffs in continuous time (e.g., a payoff of 0.6 at the beggining/end/middle of each year living with disease). 

By comparison, the GBD adopts a continuous time discounting appaoch whereby disability weights accrue as a "flow" of payoffs over time. We can convert a cycle payoff to this "flow" scale by applying a continuous time adjustment factor,

$$
\frac{1}{r}(1-e^{-rt})
$$ 

This factor is multiplied by the disability weight to define the YLD payoff vector for each cycle. 

```{r}
v_disc_h =  # Continuous time discounting
  exp(-params$r_v_disc_h_Delta_t  * 0:(params$omega))

# YLD
yld_ = with(params1,(matrix(c(0,
              dw_S1 * Delta_t * (1/r_v_disc_h_Delta_t) * (1 - exp(-r_v_disc_h_Delta_t)) ,
              0,
              0),
            dimnames = list(c(
                c(v_tr_names,v_ab_names)
            ), c("DW")))
))
yld_ <- 
  with(params1,{
    v_tx_names %>% map(~({
         yld_
    }))
  }) %>% 
  set_names(params1$v_tx_names)

YLDt <- 
   with(params1, {
    v_tx_names %>% map( ~ ({
      P = m_P[[.x]]
      occ <- v_s0T
      d <- yld_[[.x]]
      out <- 0
      P %>% map(~({
        occ <<- occ %*% .x 
        out <<-  occ %*% d
      })) %>% 
        map(~(data.frame(.x))) %>% 
        bind_rows()
    }))
  })  %>% 
  set_names(params1$v_tx_names)
YLD = YLDt %>% map(~sum(.x*  v_disc_h * gen_wcc(params1$omega, method = "none")))

# YLL
new_deaths_from_disease <- 
    map(trace1,~({
        c(0,diff(.x[,"DS"]))
    })) 

remaining_life_expectancy <- 
    with(params1,(1/r_v_disc_h) * (1 - exp(-r_v_disc_h * f_ExR(ages_trace))))
    
YLLt <- 
    new_deaths_from_disease %>% map(~(.x * remaining_life_expectancy ))

YLL <- 
    YLLt %>% map(~(sum(.x * v_disc_h * gen_wcc(params1$omega,method = params1$cycle_correction))))

DALY <- 
    map2(YLL,YLD,~(.x + .y))
DALY
```

This results in a DALY value that exactly matches the example using the GBD equations in @rushby2001. 

## Adapting the Discrete Time Discounting Formula

A more common approach is to use the discrete time discounting formula $\frac{1}{(1+r)^t}$. This formula also assumes a series of discrete payoffs, though in discrete time. 

We can replicate the GBD continuous time discounting approach using the discrete time formula by converting the discount rate $r$ using $r' = exp(r)-1$ @larson2013calculating. We then plug $r'$ into the discrete time discounting formula and simplify to obtain a discounting factor of $(\frac{1}{1+e^{r}})^t-1=\frac{1}{1+e^{rt}-1}=\frac{1}{e^{rt}}=e^{-rt}$, which is exactly the discounting factor used above. It follows trivially, then, that we obtain the same answer as above:

```{r}
v_disc_h = 
  1/exp(params1$r_v_disc_c*(0:(params$omega)))

# YLD
yld_ = with(params1,(matrix(c(0,
              dw_S1 * Delta_t * (1/r_v_disc_h_Delta_t) * (1 - exp(-r_v_disc_h_Delta_t)) ,
              0,
              0),
            dimnames = list(c(
                c(v_tr_names,v_ab_names)
            ), c("DW")))
))
yld_ <- 
  with(params1,{
    v_tx_names %>% map(~({
         yld_
    }))
  }) %>% 
  set_names(params1$v_tx_names)

YLDt <- 
   with(params1, {
    v_tx_names %>% map( ~ ({
      P = m_P[[.x]]
      occ <- v_s0T
      d <- yld_[[.x]]
      out <- 0
      P %>% map(~({
        occ <<- occ %*% .x 
        out <<-  occ %*% d
      })) %>% 
        map(~(data.frame(.x))) %>% 
        bind_rows()
    }))
  })  %>% 
  set_names(params1$v_tx_names)
YLD = YLDt %>% map(~sum(.x*  v_disc_h * gen_wcc(params1$omega, method = "none")))

# YLL
new_deaths_from_disease <- 
    map(trace1,~({
        c(0,diff(.x[,"DS"]))
    })) 

remaining_life_expectancy <- 
    with(params1,(1/r_v_disc_h) * (1 - exp(-r_v_disc_h * f_ExR(ages_trace))))
    
YLLt <- 
    new_deaths_from_disease %>% map(~(.x * remaining_life_expectancy ))

YLL <- 
    YLLt %>% map(~(sum(.x * v_disc_h * gen_wcc(params1$omega,method = params1$cycle_correction))))

DALY <- 
    map2(YLL,YLD,~(.x + .y))
DALY
```


```{r}

```

