Diagnostic of an EpiModel Module

Introduction

In this tutorial we explore an easy way to diagnose an EpiModel module. We assume that the module does not crash but produces weird results that we want to inspect directly. If the module is crashing the model, use options(error = recover) to enter debug mode when the crash happens.

For this tutorial, we will use an EpiModelHIV model and we will debug the position module.

Setup

First we load the libraries as usual and define the param and init object.

library("EpiModelHIV")
# pkgload::load_all("../EpiModelHIV-p.git/DoxyPEP")

epistats <- readRDS("data/intermediate/estimates/epistats-local.rds")
netstats <- readRDS("data/intermediate/estimates/netstats-local.rds")
est      <- readRDS("data/intermediate/estimates/netest-local.rds")

param <- param.net(
  data.frame.params   = read.csv("data/input/params.csv"),
  netstats            = netstats,
  epistats            = epistats
)

init <- init_msm(
  prev.ugc = 0.1,
  prev.rct = 0.1,
  prev.rgc = 0.1,
  prev.uct = 0.1
)

The main difference occurs in control, where we set the argument raw.output to TRUE.

  control <- control_msm(
    nsteps              = 100,
    nsims               = 1,
    ncores              = 1,
    cumulative.edgelist = TRUE,
    truncate.el.cuml    = 0,
    verbose             = TRUE,
    raw.output          = TRUE ### NEW ###
  )

This prevents netsim from post processing the simulation at the end of the run. Leaving us with a list of raw dat objects. These are the object used internally by the modules.

Running the simulation

The simulation is ran as usual with a call to netsim. However, here we save the result as dat_list for clarity and create a dat object by getting the first element of dat_list.

dat_list <- netsim(est, param, init, control)
dat <- dat_list[[1]]

At this point, dat is exactly the same thing as the dat object on the last step of the netsim simulation. We save it in order to bypass these steps next time.

saveRDS(dat, "dat.rds")

Getting ready to debug

First, we load dat from the saved file and allow the simulation to run for longer. To do so we add 100 to the dat$control$nsteps.

dat <- readRDS("dat.rds")
dat$control$nsteps <- dat$control$nsteps + 100

Then we run a few more steps of simulation on dat using the EpiModel internal function netsim_run_nsteps. Notice that we use ::: and not :: here. This allow us to access unexported functions.

EpiModel:::netsim_run_nsteps(dat, 30, 1)

This stage is meant to generate a different state of dat each time. Very useful when debugging a stochastic error.

Note: we cannot simply save the state of dat with dat_sav <- dat and then restore it with dat <- dat_sav. That's because it contains C object that will be altered. Reloading the "dat.rds" file resets it correctly though.

Running the simulation up to the module to debug

Now we will want to put our dat object exactly in the state it would be before entering the module we want to debug.

First, we set the at time object to it's correct value.

at <- get_current_timestep(dat) + 1

Now we check how many modules needs to be run before we reach the module to debug. The _modules_list elements to dat contains a named list of modules.

names(dat[["_modules_list"]])
#>
#>  [1] "aging.FUN"        "departure.FUN"    "arrival.FUN"
#>  [4] "partident.FUN"    "hivtest.FUN"      "hivtx.FUN"
#>  [7] "hivprogress.FUN"  "hivvl.FUN"        "resim_nets.FUN"
#> [10] "summary_nets.FUN" "acts.FUN"         "condoms.FUN"
#> [13] "position.FUN"     "prep.FUN"         "hivtrans.FUN"
#> [15] "stitrans.FUN"     "stirecov.FUN"     "stiscreen.FUN"
#> [18] "stitx.FUN"        "prev.FUN"         "cleanup.FUN"

Here we see that position.FUN, our module of interest is the 13th. Therefore, we need to run all modules from 1 to 12.

for (mod in dat[["_modules_list"]][1:12]) {
  dat <- mod(dat, at)
}

Note: applying each module sequentially to dat is actually what netsim does at each timestep.

Debugging the position module

At this point, dat is in the state it should before entering the 13th module (postion.FUN).

Now we simply paste the module's code and run it interactively, inspecting each step as we need to find our error.

# position_msm <- function(dat, at) {
  al <- dat[["temp"]][["al"]]
  if (nrow(al) == 0) {
    return(dat)
  }

  # Attributes
  role.class <- get_attr(dat, "role.class")
  ins.quot   <- get_attr(dat, "ins.quot")

  # Parameters

  ## Process
  p1.role.class <- role.class[al[, "p1"]]
  p2.role.class <- role.class[al[, "p2"]]

  ins <- rep(NA, length(p1.role.class))
  ins[which(p1.role.class == 0)] <- 1
  ins[which(p1.role.class == 1)] <- 0
  ins[which(p2.role.class == 0)] <- 0
  ins[which(p2.role.class == 1)] <- 1

  # Versatile MSM
  vv <- which(p1.role.class == 2 & p2.role.class == 2)
  p1.ins.prob <- ins.quot[al[, 1][vv]] /
                 (ins.quot[al[, 1][vv]] + ins.quot[al[, 2][vv]])
  p1.ins <- runif(length(vv)) < p1.ins.prob
  ins[vv[p1.ins]] <- 1
  ins[vv[!p1.ins]] <- 0

  ## Output
  dat[["temp"]][["al"]] <- cbind(al, ins)

  return(dat)
# }