R/dcm.R

In EpiModel: Mathematical Modeling of Infectious Disease Dynamics

#' @title Deterministic Compartmental Models
#'
#' @description Solves deterministic compartmental epidemic models for
#'              infectious disease.
#'
#' @param param Model parameters, as an object of class [param.dcm()].
#' @param init Initial conditions, as an object of class [init.dcm()].
#' @param control Control settings, as an object of class
#'        [control.dcm()].
#'
#' @details
#' The `dcm` function uses the ordinary differential equation solver in
#' the `deSolve` package to model disease as a deterministic compartmental
#' system. The parameterization for these models follows the standard approach
#' in `EpiModel`, with epidemic parameters, initial conditions, and control
#' settings.
#'
#' The `dcm` function performs modeling of both base model types and
#' original models with new structures. Base model types include one-group
#' and two-group models with disease types for Susceptible-Infected (SI),
#' Susceptible-Infected-Recovered (SIR), and Susceptible-Infected-Susceptible
#' (SIS). Both base and original models require the `param`,
#' `init`, and `control` inputs.
#'
#' Solver settings such as tolerances are passed to `deSolve` through the
#' `ode.args` argument of [control.dcm()], and the `solver` argument replaces
#' the solver altogether, for example to restart it at abrupt changes in a
#' model's rates.
#'
#' @return
#' A list of class `dcm` with the following elements:
#'
#'  * **param:** the epidemic parameters passed into the model through
#'        `param`, with additional parameters added as necessary.
#'  * **control:** the control settings passed into the model through
#'        `control`, with additional controls added as necessary.
#'  * **epi:** a list of data frames, one for each epidemiological
#'        output from the model. Outputs for base models always include the
#'        size of each compartment, as well as flows in, out of, and between
#'        compartments.
#'
#'
#' @references
#' Soetaert K, Petzoldt T, Setzer W. Solving Differential Equations in
#' R: Package deSolve. Journal of Statistical Software. 2010; 33(9): 1-25.
#' \doi{10.18637/jss.v033.i09}.
#'
#' @keywords model
#'
#' @seealso Extract the model results with [as.data.frame.dcm()].
#' Summarize the time-specific model results with [summary.dcm()].
#' Plot the model results with [plot.dcm()]. Plot a compartment flow
#' diagram with [comp_plot()].
#'
#' @export
#'
#' @examples
#' ## Example 1: SI Model (One-Group)
#' # Set parameters
#' param <- param.dcm(inf.prob = 0.2, act.rate = 0.25)
#' init <- init.dcm(s.num = 500, i.num = 1)
#' control <- control.dcm(type = "SI", nsteps = 500)
#' mod1 <- dcm(param, init, control)
#' mod1
#' plot(mod1)
#'
#' ## Example 2: SIR Model with Vital Dynamics (One-Group)
#' param <- param.dcm(inf.prob = 0.2, act.rate = 5,
#'                    rec.rate = 1/3, a.rate = 1/90, ds.rate = 1/100,
#'                    di.rate = 1/35, dr.rate = 1/100)
#' init <- init.dcm(s.num = 500, i.num = 1, r.num = 0)
#' control <- control.dcm(type = "SIR", nsteps = 500)
#' mod2 <- dcm(param, init, control)
#' mod2
#' plot(mod2)
#'
#' ## Example 3: SIS Model with act.rate Sensitivity Parameter
#' param <- param.dcm(inf.prob = 0.2, act.rate = seq(0.1, 0.5, 0.1),
#'                    rec.rate = 1/50)
#' init <- init.dcm(s.num = 500, i.num = 1)
#' control <- control.dcm(type = "SIS", nsteps = 500)
#' mod3 <- dcm(param, init, control)
#' mod3
#' plot(mod3)
#'
#' ## Example 4: SI Model with Vital Dynamics (Two-Group)
#' param <- param.dcm(inf.prob = 0.4,  inf.prob.g2 = 0.1,
#'                    act.rate = 0.25, balance = "g1",
#'                    a.rate = 1/100, a.rate.g2 = NA,
#'                    ds.rate = 1/100, ds.rate.g2 = 1/100,
#'                    di.rate = 1/50, di.rate.g2 = 1/50)
#' init <- init.dcm(s.num = 500, i.num = 1,
#'                  s.num.g2 = 500, i.num.g2 = 0)
#' control <- control.dcm(type = "SI", nsteps = 500)
#' mod4 <- dcm(param, init, control)
#' mod4
#' plot(mod4)
#'
#' ## Example 5: SI Model with Varying Initial Conditions
#' param <- param.dcm(inf.prob = 0.2, act.rate = 0.25)
#' init <- init.dcm(s.num = 500, i.num = c(1, 5, 25))
#' control <- control.dcm(type = "SI", nsteps = 500)
#' mod5 <- dcm(param, init, control)
#' mod5
#' plot(mod5)
#'
dcm <- function(param, init, control) {
  check.control.class("dcm", "EpiModel dcm")

  crosscheck.dcm(param, init, control)

  # Model selection ---------------------------------------------------------
  if (is.null(control$new.mod)) {

    if (control$type == "SI") {
      if (param$groups == 1) {
        if (param$vital == FALSE) {
          model <- mod_SI_1g_cl
          init <- c(init, si.flow = 0)
        } else {
          model <- mod_SI_1g_op
          init <- c(init, si.flow = 0, a.flow = 0, ds.flow = 0, di.flow = 0)
        }
      } else {
        if (param$vital == FALSE) {
          model <- mod_SI_2g_cl
          init <- c(init, si.flow = 0, si.flow.g2 = 0)
        } else {
          model <- mod_SI_2g_op
          init <- c(init,
                    si.flow = 0, a.flow = 0, ds.flow = 0, di.flow = 0,
                    si.flow.g2 = 0, a.flow.g2 = 0, ds.flow.g2 = 0,
                    di.flow.g2 = 0)
        }
      }
    }
    if (control$type == "SIR") {
      if (param$groups == 1) {
        if (param$vital == FALSE) {
          model <- mod_SIR_1g_cl
          init <- c(init, si.flow = 0, ir.flow = 0)
        } else {
          model <- mod_SIR_1g_op
          init <- c(init, si.flow = 0, ir.flow = 0, a.flow = 0,
                    ds.flow = 0, di.flow = 0, dr.flow = 0)
        }
      } else {
        if (param$vital == FALSE) {
          model <- mod_SIR_2g_cl
          init <- c(init, si.flow = 0, ir.flow = 0,
                    si.flow.g2 = 0, ir.flow.g2 = 0)
        } else {
          model <- mod_SIR_2g_op
          init <- c(init, si.flow = 0, ir.flow = 0, a.flow = 0,
                    ds.flow = 0, di.flow = 0, dr.flow = 0,
                    si.flow.g2 = 0, ir.flow.g2 = 0, a.flow.g2 = 0,
                    ds.flow.g2 = 0, di.flow.g2 = 0, dr.flow.g2 = 0)
        }
      }
    }
    if (control$type == "SIS") {
      if (param$groups == 1) {
        if (param$vital == FALSE) {
          model <- mod_SIS_1g_cl
          init <- c(init, si.flow = 0, is.flow = 0)
        } else {
          model <- mod_SIS_1g_op
          init <- c(init, si.flow = 0, is.flow = 0,
                    a.flow = 0, ds.flow = 0, di.flow = 0)
        }
      } else {
        if (param$vital == FALSE) {
          model <- mod_SIS_2g_cl
          init <- c(init, si.flow = 0, is.flow = 0,
                    si.flow.g2 = 0, is.flow.g2 = 0)
        } else {
          model <- mod_SIS_2g_op
          init <- c(init, si.flow = 0, is.flow = 0,
                    a.flow = 0, ds.flow = 0, di.flow = 0,
                    si.flow.g2 = 0, is.flow.g2 = 0,
                    a.flow.g2 = 0, ds.flow.g2 = 0, di.flow.g2 = 0)
        }
      }
    }

  } else {
    # Sub in new model
    model <- control$new.mod
  }


  # Print model option ------------------------------------------------------
  if (control$print.mod == TRUE) {
    return(print(model))
  }


  # Sensitivity parameters and initial conditions ----------------------------
  if (control$sens.param == FALSE) {
    control$nruns <- 1
  } else {
    control$nruns <- max(sapply(param, length), sapply(init, length))
  }

  # Validate that all varying vectors have the same length
  if (control$nruns > 1) {
    all.lengths <- c(sapply(param, length), sapply(init, length))
    long.lengths <- all.lengths[all.lengths > 1]
    if (length(unique(long.lengths)) > 1) {
      stop("All varying parameters and initial conditions must have the ",
           "same length")
    }

    longv <- which(sapply(param, length) == control$nruns)
    longvn <- names(longv)
    lim.p <- param[!(names(param) %in% longvn)]

    longi <- which(sapply(init, length) == control$nruns)
    longin <- names(longi)
    lim.i <- init[!(names(init) %in% longin)]
  }


  # Model runs --------------------------------------------------------------
  verbose.dcm(control, type = "startup")
  for (s in 1:control$nruns) {

    ## Sensitivity parameter input
    if (control$nruns > 1) {
      all.p <- as.list(lim.p)
      for (j in seq_along(longvn)) {
        all.p[[longvn[j]]] <- param[[longvn[j]]][s]
      }
    } else {
      all.p <- param
    }

    ## Per-run initial conditions
    if (control$nruns > 1 && length(longin) > 0) {
      all.i <- lim.i
      for (j in seq_along(longin)) {
        all.i[[longin[j]]] <- init[[longin[j]]][s]
      }
    } else {
      all.i <- init
    }
    t0 <- setNames(as.numeric(all.i), names(all.i))

    ## Timesteps
    if (length(control$nsteps) == 1) {
      times <- seq(1, control$nsteps, control$dt)
    } else {
      times <- control$nsteps
    }


    ## Differential equation solvers
    if (!is.null(control$solver)) {
      df <- data.frame(control$solver(y = t0,
                                      times = times,
                                      func = model,
                                      parms = all.p,
                                      control = control))
    } else if (control$dede == FALSE) {
      df <- data.frame(do.call(ode, c(list(y = t0,
                                           times = times,
                                           func = model,
                                           parms = all.p,
                                           method = control$odemethod),
                                      control$ode.args)))
    } else {
      df <- data.frame(do.call(dede, c(list(y = t0,
                                            times = times,
                                            func = model,
                                            parms = all.p),
                                       control$ode.args)))
    }

    ## Recalculate Flows
    isFlow <- grep(".flow", x = names(df))
    if (length(isFlow) > 0) {
      for (j in isFlow) {
        df[, j] <- c(diff(df[, j]), NA)
      }
    }


    # Output ------------------------------------------------------------------

    ## Save to out object
    if (s == 1) {
      out <- saveout.dcm(df, s, param, control)
    } else {
      out <- saveout.dcm(df, s, param, control, out)
    }

    ## Progress
    verbose.dcm(control, type = "progress", s)

  }

  out$init <- init

  class(out) <- "dcm"
  invisible(out)
}

## Match a requested time to the time vector of a solved `dcm` model.
##
## `control$nsteps` may hold either the maximum time or the full vector of
## times, and `dt` may be non-integer, so requested times are matched against
## `control$timesteps` rather than compared against a scalar. The comparison
## allows for the floating point error of a non-integer `dt` sequence, in which
## a time such as 1.3 is not stored exactly. Returns the index into
## `timesteps`.
match_dcm_time <- function(at, timesteps) {
  if (length(at) != 1 || !is.numeric(at) || is.na(at)) {
    stop("Specify at as a single numeric time step")
  }

  tol <- sqrt(.Machine$double.eps) * max(1, abs(at))
  idx <- which(abs(timesteps - at) <= tol)

  if (length(idx) == 0) {
    stop("Specify at as one of the time steps of the model, between ",
         min(timesteps), " and ", max(timesteps))
  }

  return(idx[1])
}


## Describe the time coordinate of a solved `dcm` model for printing. Reports
## the number of time steps, adding the range when the model was solved over
## times other than 1, 2, ..., n, as a non-integer `dt` or an explicit vector of
## times in `nsteps` produces.
format_dcm_timesteps <- function(control) {
  timesteps <- control$timesteps

  if (is.null(timesteps)) {
    return(paste(control$nsteps, collapse = " "))
  }
  if (identical(as.numeric(timesteps), as.numeric(seq_along(timesteps)))) {
    return(paste(length(timesteps)))
  }

  paste0(length(timesteps), " (", min(timesteps), " to ",
         max(timesteps), ")")
}