# Reproducing Exact Results from DMHEE

#### 2020-05-11

The purpose of this vignette is to present how to reproduce exactly the results from Decision Modelling for Health Economic Evaluation. 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 beginning of each year (equivalent to transition happening at the end + ignoring the first year) with `method = "beginning"`. Since with this method the first year is actually the second, costs should be discounted from the start 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(markov_cycle <= 2, .509, 1.00)`.

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

    par_mod <- define_parameters(
      rr = ifelse(markov_cycle <= 2, .509, 1),
      cost_lami = ifelse(markov_cycle <= 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
    )

    ## No named state -> generating names.

    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
    )

    ## No named state -> generating names.

    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
    )

    ## No named state -> generating names.

    mod_comb <- define_strategy(
      transition = mat_comb,
      A_comb,
      B_comb,
      C_comb,
      D_comb
    )

    ## No named state -> generating names.

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

    ## 2 strategies run for 20 cycles.
    ##
    ## Initial state counts:
    ##
    ## A = 1
    ## B = 0
    ## C = 0
    ## D = 0
    ##
    ## Counting method: 'beginning'.
    ##
    ## Values:
    ##
    ##      cost_health cost_drugs cost_total life_year
    ## mono    45541.24   18203.97   44663.45  7.991207
    ## comb    48082.83   24492.28   50601.65  8.937389
    ##
    ## Efficiency frontier:
    ##
    ## mono -> comb
    ##
    ## Differences:
    ##
    ##      Cost Diff. Effect Diff.     ICER Ref.
    ## comb   5938.198    0.9461822 6275.956 mono

# 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:

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

    # 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

    ##    age sex   value
    ## 1   35   1 0.00151
    ## 2   35   0 0.00099
    ## 3   45   1 0.00393
    ## 4   45   0 0.00260
    ## 5   55   1 0.01090
    ## 6   55   0 0.00670
    ## 7   65   1 0.03160
    ## 8   65   0 0.01930
    ## 9   75   1 0.08010
    ## 10  75   0 0.05350
    ## 11  85   1 0.18790
    ## 12  85   0 0.15480

    param <- define_parameters(
        age_init = 60,
        sex = 0,
        # age increases with cycles
        age = age_init + markov_cycle,

        # 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 * ((markov_cycle - 1) ^ gamma -
                                         markov_cycle ^ gamma)),
        np1RR = 1 - exp(lambda * rrNP1 * ((markov_cycle - 1) ^ gamma -
                                            markov_cycle ^ 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 = "end",
      init = c(1, 0, 0, 0, 0)
    )
    summary(res_mod)

    ## 2 strategies run for 60 cycles.
    ##
    ## Initial state counts:
    ##
    ## PrimaryTHR = 1
    ## SuccessP = 0
    ## RevisionTHR = 0
    ## SuccessR = 0
    ## Death = 0
    ##
    ## Counting method: 'end'.
    ##
    ## Values:
    ##
    ##           utility     cost
    ## standard 14.65283 512.4327
    ## np1      14.69734 610.3112
    ##
    ## Efficiency frontier:
    ##
    ## standard -> np1
    ##
    ## Differences:
    ##
    ##     Cost Diff. Effect Diff.     ICER     Ref.
    ## np1   97.87858   0.04451095 2198.977 standard