Skip to contents

Healthy–Sick–Sicker Microsimulation

Wael Mohammed

2025-10-08

Source: vignettes/sickSickerMicrosimPack.Rmd
sickSickerMicrosimPack.Rmd

Introduction

This vignette illustrates how to use the functions provided by the sickSickerMicrosimPack package to simulate a cohort of individuals progressing through a four–state health model. The health states include Healthy (H), Sick (S1), Sicker (S2) and Dead (D). Transition probabilities can depend on the time spent in a state as well as on individual covariates such as age and sex. Costs and utilities are accumulated each cycle and may be discounted.

For more detailed guidance on creating R packages, refer to the online version of the second edition of the R Packages book (Wickham and Bryan 2023) and the Posit package development cheatsheet. These resources provide comprehensive overviews of best practices in packaging, documentation, testing and continuous integration.

Plotting the model diagram

The following diagram summarises the order of operations within the simulation from the specification of inputs through to summarising discounted results.

Transition probability matrix

At each cycle the transition probabilities are updated for every individual using the update_probsV() function. The example below constructs a simple vector of occupied states and times and then computes the corresponding transition matrix.

v_states_names <- c("H", "S1", "S2", "D")
v_occupied_state <- c("H", "S1", "S2", "D")
v_time_in_state <- c(1,2,3,1)
l_trans_probs <- list(
  p_HD = 0.005,
  p_HS1 = 0.15,
  p_S1H = 0.5,
  p_S1S2 = 0.105,
  p_S1D = 1 - exp(-3 * (-log(1-0.005))),
  p_S2D = 1 - exp(-10 * (-log(1-0.005))),
  rp_S1 = 0.2,
  rp_S2 = 0.29
)

m_probs <- update_probsV(
  v_states_names = v_states_names, 
  v_occupied_state = v_occupied_state,
  l_trans_probs = l_trans_probs,
  v_time_in_state = v_time_in_state
)

knitr::kable(
  x = m_probs,
  caption = "Transition probabilities for a four‑individual cohort"
)
Transition probabilities for a four‑individual cohort
H S1 S2 D
H 0.845 0.1500000 0.0000000 0.0050000
S1 0.500 0.3741674 0.1050000 0.0208326
S2 0.000 0.0000000 0.9105244 0.0894756
D 0.000 0.0000000 0.0000000 1.0000000

Running a small simulation

The following code runs a small simulation for ten individuals over five cycles. Individual characteristics are sampled and passed to run_microSimV() together with cost and utility parameters. The average total costs and QALYs per individual are returned.

set.seed(123)
num_i <- 10
num_cycles <- 5
v_start <- rep("H", num_i)
m_feats <- cbind(age = rnorm(num_i, 50, 3), sex = sample(c(0,1), num_i, replace = TRUE))
v_states_costs <- c(H = 2000, S1 = 4000, S2 = 15000, D = 0)
v_cost_coeffs  <- c(age = 11.5, sex = 300)
v_states_utilities <- c(H = 1, S1 = 0.75, S2 = 0.5, D = 0)
v_util_coeffs <- c(age = -0.0018, sex = -0.015)
v_util_t_decs <- c(S1 = -0.0015, S2 = -0.0020)
l_trans_probs <- list(p_HD = 0.005, p_HS1 = 0.15, p_S1H = 0.5,
                      p_S1S2 = 0.105, p_S1D = 0.01488751,
                      p_S2D = 0.04877166, rp_S1 = 0.2, rp_S2 = 0.29)
res <- run_microSimV(
  v_starting_states = v_start,
  num_i = num_i,
  num_cycles = num_cycles,
  m_indi_features = m_feats,
  v_states_names = v_states_names,
  v_states_costs = v_states_costs,
  v_cost_coeffs = v_cost_coeffs,
  v_states_utilities = v_states_utilities,
  v_util_coeffs = v_util_coeffs,
  v_util_t_decs = v_util_t_decs,
  l_trans_probs = l_trans_probs,
  discount_rate_costs = 0.03,
  discount_rate_QALYs = 0.015,
  cycle_length = 1,
  starting_seed = 1
)

res$mean_costs
#> [1] 14428.73
res$mean_qalys
#> [1] 5.167632

The methodology behind this package and its application to cost–effectiveness modelling is described in detail in our tutorial on packaging cost‑effectiveness models in R (Smith, Mohammed, and Schneider 2023).

References

Smith, Robert, Wael Mohammed, and Paul Schneider. 2023. “Packaging Cost-Effectiveness Models in r: A Tutorial.” Wellcome Open Research 8.
Wickham, Hadley, and Jennifer Bryan. 2023. R Packages, 2nd Edition. O’Reilly Media.