Sick-sicker model

In this example we consider a simple scenario where a person can occupy one of four states: healthy, sick, sicker and dead. At each iteration, we compute the transition probabilities in the population and change the state accordingly. In addition, there are costs and quality adjusted life years associated with each state.

For more information and details about the model, please refer to the original article: Krijkamp EM, Alarid-Escudero F, Enns EA, Jalal HJ, Hunink MGM, Pechlivanoglou P. Microsimulation modeling for health decision sciences using R: A tutorial. Med Decis Making. 2018;38(3):400–22.

The code used for comparative performance purposes can be found on GitHub: https://github.com/DARTH-git/Microsimulation-tutorial

R-code

Here we provide the R-code for implementing the model in the Sima framework. We also see that the performance of Sima in this example is superior.

##
## Implementation of the Sick-Sicker model that appears in:
##
## Krijkamp EM, Alarid-Escudero F, Enns EA, Jalal HJ, Hunink MGM, Pechlivanoglou P. 
##   Microsimulation modeling for health decision sciences using R: 
##   A tutorial. Med Decis Making. 2018;38(3):400–22.
##

library(Sima)
library(data.table)
library(dqrng)
rm(list = ls())

# Vectorized multinomial distribution
rmult <- function(prob) {
    dims <- dim(prob)
    k <- dims[1]
    n <- dims[2]
    for (i in 2:k) {
        prob[i,] <- prob[i,] + prob[i-1,]
    }
    u <- rep(dqrng::dqrunif(n), rep(k, n))
    return(as.integer(1 + colSums(u > prob)))
}

# Initialize the population for the simulation
# -- state : The possible states for an individual (healthy, sick, sicker, dead)
# -- cost  : The initial cost associated with each individual
# -- treatment : Treatment status, 0 = not treated, 1 = treated
# -- QALY : Initial quality adjusted life years
initializer <- function(n) {
    return(data.table(
        state = rep(1L, n),
        cost = rep(2000.0, n),
        treatment = rep(0L, n),
        QALY = rep(0, n)
    ))
}

health_status <- ManipulationEvent$new(
    name = "Health Status",
    description = "Transitions between states (healthy, sick, sicker, dead)",
    parameters = list(
        hd   = 0.005, 
        hs1  = 0.15, 
        s1h  = 0.5,
        s1s2 = 0.105,
        s1d  = 1 - exp(-3  * -log(1 - 0.005)),
        s2d  = 1 - exp(-10 * -log(1 - 0.005))
    ),
    derived = list(
        H  = expression(c(1 - hd - hs1, s1h, 0, 0)),
        S1 = expression(c(hs1, 1 - s1h - s1s2 - s1d, 0, 0)),
        S2 = expression(c(0, s1s2, 1 - s2d, 0)),
        D  = expression(c(hd, s1d, s2d, 1))
    ),
    allocated = list(
        prob = expression(matrix(0, 3, status[ , .N]))
    ),
    mechanism = expression({
        state <- status$state
        prob[1,] <- H[state]
        prob[2,] <- S1[state]
        prob[3,] <- S2[state]
        status[ , state := rmult(prob)]
    })
)

costs <- ManipulationEvent$new(
    name = "Costs",
    description = "Individual costs based on health and treatment",
    parameters = list(
        state_cost = c(2000.0, 4000.0, 15000.0, 0.0),
        treat_cost = 12000.0
    ),
    mechanism = expression({
        status[ , cost := state_cost[state] + (state == 2L) * treatment * treat_cost]
    })
)

quality <- ManipulationEvent$new(
    name = "QALY",
    description = "Quality-adjusted life years",
    parameters = list(
        state_qaly = c(1.0, 0, 0.5, 0.0),
        S1_qaly = 0.75,
        treat_qaly = 0.95
    ),
    mechanism = expression({
        status[ , QALY := state_qaly[state] + (state == 2L) * 
                  (treatment * treat_qaly + (1 - treatment) * S1_qaly)]
    })
)

treatment_intervention <- ManipulationEvent$new(
    name = "Treatment intervention",
    description = "Assign treatment for everyone",
    mechanism = expression({
        status[ , treatment := 1L]
    })
)

sim <- Simulator$new(
    initializer = initializer,
    man_events = list(health_status, costs, quality),
    seeds = 0,
    keep_init = TRUE,
    n = 1e5
)

# A monitor to save the entire history
sim$initialize_monitor(
    aggregators = list(function(dt) data.table::copy(dt)),
    intervals = 1
)

# Run, reset, intervene, run again
system.time({
    sim$run(30, show_progress = FALSE)
    sim$reset()
    sim$intervene(man_events = list(treatment_intervention))
    sim$run(30, show_progress = FALSE)
})
# user  system elapsed 
# 1.49    0.14    1.64 

# Comparison with the microsimulation R package
system.time({
  devtools::source_url("https://raw.githubusercontent.com/DARTH-git/
                       Microsimulation-tutorial/master/Appendix%20D_online_supp.R")
})

# Measured time
comp.time
# Time difference of 5.889181 secs