## ----packages, echo=FALSE,message=FALSE,warning=FALSE---------------------------------------------
library(tidyverse)
library(ggplot2)
library(heemod)
library(diagram)


## ---- echo=FALSE----------------------------------------------------------------------------------

## Define transition matrix for first example

mm_simp_heemod <- define_transition(
  state_names = c("healthy", "ill", "dead"),
    .97, .02, 0.01,
    0,    0.96, 0.04,
    0,    0,   1
  )

# display transition matrix
mm_simp_heemod

# plot corresponding diagram
# showing transitions with probabilities
plot(mm_simp_heemod)


# define costs and outcomes associated with each state
healthy_state = define_state(
  cost = 40,
  QALYs = 0.95
)
ill_state = define_state(
  cost = 1000,
  QALYs = 0.6
)
dead_state = define_state(
  cost = 0,
  QALYs = 0
)

## the strategy combines the transition matrix with the information about each state

strat_simp = define_strategy(
  transition = mm_simp_heemod,
  healthy = healthy_state, #[<name in transition matrix> = <defined state>]
  ill=ill_state, 
  dead=dead_state
)

# We can now run this model with the information in the strategy
# for a number of cycles
# In our model cost is one number, but it can be split into 
# `cost_drugs` and `cost_health` to separate the cost of treatment
# from the general healthcare cost (we will see this in the next example)

# we also define the effect (in our case in QALYs, which is what we've used
# when defining the states)

res_simp = run_model(
  sim = strat_simp,
  cycles = 50,
  cost = cost,
  effect= QALYs
)
## See the first few rows of the counts in each state
head(get_counts(res_simp))

# You could make a data frame of all the counts using
countsA = get_counts(res_simp)
# and for example subset this by state
countsA_ill = countsA[countsA$state_names == "ill",]

## see the first few rows of the costs and QALYs for each state
head(countsA)
head(countsA_ill)
# Plot the populations of each state over time
plot(res_simp)
summary(res_simp)


## Section 4 - comparing two strategies


# Define transition matrix for drug A
tm_drugA = define_transition(
  state_names = c("mild", "severe", "dead"),
  0.75, 0.23, 0.02,
  0, 0.9, 0.1, 
  0, 0, 1
)
tm_drugA

# Define transition matrix for drug B
# using relative risk

rr = 0.7  ## Our relative risk
tm_drugB = define_transition(
  state_names = c("mild", "severe", "dead"),
  C, 0.23*rr, 0.02*rr,
  0, C, 0.1*rr,
  0, 0, 1
)

tm_drugB


# setting up some initial cost variables
cost_drugA = 2000
cost_drugB = 2700

# Define states
# splitting drug cost according to strategy using 'dispatch_strategy'
# applying discounting to both cost_health and cost_drug

state_mild= define_state(
  # health and drugs costs decrease by a factor of 0.03 (3%) each year
  cost_health = discount(1500, 0.03),
  cost_drugs = discount(
    dispatch_strategy(
      drugA = cost_drugA, 
      drugB = cost_drugB), 
    0.03),
  cost_total = cost_health + cost_drugs,
  QALYs = 0.9  # QALYs is our choice of effect here. 
  ## We could use a different measure, so long as we linked this to the 
  ## 'effect' argument in run_model.
)

state_severe = define_state(
  cost_health = discount(3000, 0.03),
  cost_drugs = discount(
    dispatch_strategy(
      drugA = cost_drugA, 
      drugB = cost_drugB), 
    0.03),
  cost_total = cost_health + cost_drugs,
  QALYs = 0.6
)

state_dead = define_state(
  cost_health = 0,
  cost_drugs = 0,
  cost_total = 0,
  QALYs = 0
)


# Set up ach strategy

strat_drugA = define_strategy(
  transition = tm_drugA,
  mild = state_mild,
  severe = state_severe,
  dead = state_dead
)

strat_drugB = define_strategy(
  transition = tm_drugB,
  mild = state_mild,
  severe = state_severe,
  dead = state_dead
)

# Run the model
# the names 'drugA' and 'drugB' here match to the names used in 'dispatch_strategy'

res_drugsAB <- run_model(
  drugA = strat_drugA,
  drugB = strat_drugB,
  cycles = 50,
  cost = cost_total,
  effect = QALYs
)

## View and plot results

summary(res_drugsAB)

plot(res_drugsAB, type = "counts", panel = "by_strategy")

plot(res_drugsAB, type = "counts", panel = "by_state")

plot(res_drugsAB, type = "values", panel = "by_value",
     free_y = TRUE)