R/epiestim_config.R

In RussellXu/rtestimate: Estimation Of Reproductive Number

#' This function was contributed by Rich Fitzjohn. It modifies default arguments
#' using user-provided values. The argument 'strict' triggers and error
#' behaviour: if strict==TRUE: all new values need to be part of the defaults.
#'
#' @param defaults
#' @param x
#' @param strict
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
modify_defaults <- function(defaults, x, strict = TRUE) {
  extra <- setdiff(names(x), names(defaults))
  if (strict && (length(extra) > 0L)) {
    stop("Additional invalid options: ", paste(extra, collapse=", "))
  }
  utils::modifyList(defaults, x, keep.null = TRUE) # keep.null is needed here
}


#' Set and check parameter settings of epiestim estimation
#'
#' @param ...
#' @param incid
#' @param method
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
make_config <- function(..., incid = NULL,
                        method = c("non_parametric_si",
                                   "parametric_si",
                                   )){

  config <- list(...)
  if (length(config) == 1L && is.list(config[[1]])) {
    config <- config[[1]]
  }

  ## SET DEFAULTS
  defaults <- list(t_start = NULL,
                   t_end = NULL,
                   n1 = 500,
                   n2 = 50,
                   mean_si = NULL,
                   std_si = NULL,
                   std_mean_si = NULL,
                   min_mean_si = NULL,
                   max_mean_si = NULL,
                   std_std_si = NULL,
                   min_std_si = NULL,
                   max_std_si = NULL,
                   si_distr = NULL,
                   si_parametric_distr = NULL,
                   seed = NULL,
                   mean_prior = 5,
                   std_prior = 5,
                   cv_posterior = 0.3)

  ## MODIFY CONFIG WITH ARGUMENTS ##
  config <- modify_defaults(defaults, config)

  ## checking and processing incid
  if(!is.null(incid))
  {
    incid <- process_I(incid)
    T <- nrow(incid)

    ## filling in / checking t_start and t_end
    if(is.null(config$t_start) || is.null(config$t_end))
    {
      msg <- "Default config will estimate R on weekly sliding windows.
    To change this change the t_start and t_end arguments. "
      message(msg)
      config$t_start <- seq(2, T-6)
      config$t_end <- seq(8, T)
    }else
    {
      check_times(config$t_start, config$t_end, T)
    }
  }

  class(config) <- "estimate_R_config"
  return(config)

}


#' Use incidence create a single column for the Infected person counts.
#'
#' @param incid
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
process_I <- function(incid) {
  if (inherits(incid, "incidence")) {
    I_inc   <- incid
    incid   <- as.data.frame(I_inc)
    incid$I <- rowSums(incidence::get_counts(I_inc))
  }
  vector_I        <- FALSE
  single_col_df_I <- FALSE
  if (is.vector(incid)) {
    vector_I <- TRUE
  } else if (is.data.frame(incid)) {
    if (ncol(incid) == 1) {
      single_col_df_I <- TRUE
    }
  }
  if (vector_I | single_col_df_I) {
    if (single_col_df_I) {
      I_tmp <- incid[[1]]
    } else {
      I_tmp <- incid
    }
    incid      <- data.frame(local = I_tmp, imported = rep(0, length(I_tmp)))
    I_init     <- sum(incid[1, ])
    incid[1, ] <- c(0, I_init)
  } else {
    if (!is.data.frame(incid) |
        (!("I" %in% names(incid)) &
         !all(c("local", "imported") %in% names(incid)))) {
      stop("incid must be a vector or a dataframe with either i) a column
           called 'I', or ii) 2 columns called 'local' and 'imported'.")
    }
    if (("I" %in% names(incid)) &
        !all(c("local", "imported") %in% names(incid))) {
      incid$local    <- incid$I
      incid$local[1] <- 0
      incid$imported <- c(incid$I[1], rep(0, nrow(incid) - 1))
    }
    if (incid$local[1] > 0) {
      warning("incid$local[1] is >0 but must be 0, as all cases on the first
              time step are assumed imported. This is corrected automatically
              by cases being transferred to incid$imported.")
      I_init <- sum(incid[1, c("local", "imported")])
      incid[1, c("local", "imported")] <- c(0, I_init)
    }
  }

  incid[which(is.na(incid))] <- 0
  date_col <- names(incid) == "dates"
  if (any(date_col)) {
    if (any(incid[, !date_col] < 0)) {
      stop("incid must contain only non negative integer values.")
    }
  } else {
    if (any(incid < 0)) {
      stop("incid must contain only non negative integer values.")
    }
  }

  return(incid)
}


#' Check SI model distribution
#'
#' @param si_distr
#' @param sumToOne
#' @param method
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
check_si_distr <- function(si_distr, sumToOne = c("error", "warning"),
                           method = "non_parametric_si")
  ## this only produces warnings and errors, does not return anything
{
  sumToOne <- match.arg(sumToOne)
  if (is.null(si_distr)) {
    stop(paste0("si_distr argument is missing but is required for method ",
                method, "."))
  }
  if (!is.vector(si_distr)) {
    stop("si_distr must be a vector.")
  }
  if (si_distr[1] != 0) {
    stop("si_distr should be so that si_distr[1] = 0.")
  }
  if (any(si_distr < 0)) {
    stop("si_distr must be a positive vector.")
  }
  if (abs(sum(si_distr) - 1) > 0.01) {
    if (sumToOne == "error") {
      stop("si_distr must sum to 1.")
    }
    else if (sumToOne == "warning") {
      warning("si_distr does not sum to 1.")
    }
  }
}

#' incidence's date must be an object of class date or numeric
#'
#' @param incid
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
check_dates <- function(incid) {
  dates <- incid$dates
  if (class(dates) != "Date" & class(dates) != "numeric") {
    stop("incid$dates must be an object of class date or numeric.")
  } else {
    if (unique(diff(dates)) != 1) {
      stop("incid$dates must contain dates which are all in a row.")
    } else {
      return(dates)
    }
  }
}

#' This only produces warnings and errors, does not return anything
#'
#' @param t_start
#' @param t_end
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
check_times <- function(t_start, t_end, T)
{
  if (!is.vector(t_start)) {
    stop("t_start must be a vector.")
  }
  if (!is.vector(t_end)) {
    stop("t_end must be a vector.")
  }
  if (length(t_start) != length(t_end)) {
    stop("t_start and t_end must have the same length.")
  }
  if (any(t_start > t_end)) {
    stop("t_start[i] must be <= t_end[i] for all i.")
  }
  if (any(t_start < 2 | t_start > T | t_start %% 1 != 0)) {
    stop("t_start must be a vector of integers between 2 and the number of
         timesteps in incid.")
  }
  if (any(t_end < 2 | t_end > T | t_end %% 1 != 0)) {
    stop("t_end must be a vector of integers between 2 and the number of
         timesteps in incid.")
  }
}


#' Epiestim config set up
#'
#' @param config
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
process_config <- function(config) {
  if (!("mean_prior" %in% names(config))) {
    config$mean_prior <- 5
  }

  if (!("std_prior" %in% names(config))) {
    config$std_prior <- 5
  }

  if (config$mean_prior <= 0) {
    stop("config$mean_prior must be >0.")
  }
  if (config$std_prior <= 0) {
    stop("config$std_prior must be >0.")
  }

  if (!("cv_posterior" %in% names(config))) {
    config$cv_posterior <- 0.3
  }

  return(config)
}


#' Epiestim config check
#'
#' @param config
#' @param method
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
check_config <- function(config, method) {
  if (method == "non_parametric_si") {
    check_si_distr(config$si_distr, method = method)
  }
  if (method == "parametric_si") {
    if (is.null(config$mean_si)) {
      stop("method parametric_si requires to specify the config$mean_si
           argument.")
    }
    if (is.null(config$std_si)) {
      stop("method parametric_si requires to specify the config$std_si
           argument.")
    }
    if (config$mean_si <= 1) {
      stop("method parametric_si requires a value >1 for config$mean_si.")
    }
    if (config$std_si <= 0) {
      stop("method parametric_si requires a >0 value for config$std_si.")
    }
  }
  if (method == "uncertain_si") {
    if (is.null(config$mean_si)) {
      stop("method uncertain_si requires to specify the config$mean_si
           argument.")
    }
    if (is.null(config$std_si)) {
      stop("method uncertain_si requires to specify the config$std_si
           argument.")
    }
    if (is.null(config$n1)) {
      stop("method uncertain_si requires to specify the config$n1 argument.")
    }
    if (is.null(config$n2)) {
      stop("method uncertain_si requires to specify the config$n2 argument.")
    }
    if (is.null(config$std_mean_si)) {
      stop("method uncertain_si requires to specify the config$std_mean_si
           argument.")
    }
    if (is.null(config$min_mean_si)) {
      stop("method uncertain_si requires to specify the config$min_mean_si
           argument.")
    }
    if (is.null(config$max_mean_si)) {
      stop("method uncertain_si requires to specify the config$max_mean_si
           argument.")
    }
    if (is.null(config$std_std_si)) {
      stop("method uncertain_si requires to specify the config$std_std_si
           argument.")
    }
    if (is.null(config$min_std_si)) {
      stop("method uncertain_si requires to specify the config$min_std_si
           argument.")
    }
    if (is.null(config$max_std_si)) {
      stop("method uncertain_si requires to specify the config$max_std_si
           argument.")
    }
    if (config$mean_si <= 0) {
      stop("method uncertain_si requires a >0 value for config$mean_si.")
    }
    if (config$std_si <= 0) {
      stop("method uncertain_si requires a >0 value for config$std_si.")
    }
    if (config$n2 <= 0 || config$n2 %% 1 != 0) {
      stop("method uncertain_si requires a >0 integer value for config$n2.")
    }
    if (config$n1 <= 0 || config$n1 %% 1 != 0) {
      stop("method uncertain_si requires a >0 integer value for config$n1.")
    }
    if (config$std_mean_si <= 0) {
      stop("method uncertain_si requires a >0 value for config$std_mean_si.")
    }
    if (config$min_mean_si < 1) {
      stop("method uncertain_si requires a value >=1 for config$min_mean_si.")
    }
    if (config$max_mean_si < config$mean_si) {
      stop("method uncertain_si requires that config$max_mean_si >=
           config$mean_si.")
    }
    if (config$mean_si < config$min_mean_si) {
      stop("method uncertain_si requires that config$mean_si >=
           config$min_mean_si.")
    }
    if (signif(config$max_mean_si - config$mean_si, 3) != signif(config$mean_si -
                                                                 config$min_mean_si, 3)) {
      warning("The distribution you chose for the mean SI is not centered around
              the mean.")
    }
    if (config$std_std_si <= 0) {
      stop("method uncertain_si requires a >0 value for config$std_std_si.")
    }
    if (config$min_std_si <= 0) {
      stop("method uncertain_si requires a >0 value for config$min_std_si.")
    }
    if (config$max_std_si < config$std_si) {
      stop("method uncertain_si requires that config$max_std_si >=
           config$std_si.")
    }
    if (config$std_si < config$min_std_si) {
      stop("method uncertain_si requires that config$std_si >=
           config$min_std_si.")
    }
    if (signif(config$max_std_si - config$std_si, 3) != signif(config$std_si -
                                                               config$min_std_si, 3)) {
      warning("The distribution you chose for the std of the SI is not centered
              around the mean.")
    }
  }
  if (config$cv_posterior < 0) {
    stop("config$cv_posterior must be >0.")
  }
}


#' vnapply function
#'
#' @param X
#' @param FUN
#' @param ...
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
vnapply <- function(X, FUN, ...) {
  vapply(X, FUN, numeric(1), ...)
}


#' Overall infectivity
#'
#' @param incid
#' @param si_distr
#'
#' @return
#' @export
#'
#' @examples
#'
#' @family EpiEstim
overall_infectivity <- function(incid, si_distr) {
  incid <- process_I(incid)
  T <- nrow(incid)
  check_si_distr(si_distr, "warning")
  lambda <- vector()
  lambda[1] <- NA
  for (t in seq(2, T))
    lambda[t] <- sum(si_distr[seq_len(t)] *
                       rowSums(incid[seq(t, 1), c("local", "imported")]),
                     na.rm = TRUE)
  return(lambda)
}


#' Discretized Generation Time Distribution Assuming A Shifted Gamma distribution
#'
#' @param k
#' @param mu
#' @param sigma
#'
#' @family EpiEstim

discr_si <- function(k, mu, sigma)
{
  if (sigma < 0) {
    stop("sigma must be >=0.")
  }
  if (mu <= 1) {
    stop("mu must be >1")
  }
  if (any(k < 0)) {
    stop("all values in k must be >=0.")
  }

  a <- ((mu - 1) / sigma)^2
  b <- sigma^2 / (mu - 1)

  cdf_gamma <- function(k, a, b) stats::pgamma(k, shape = a, scale = b)

  res <- k * cdf_gamma(k, a, b) +
    (k - 2) * cdf_gamma(k - 2, a, b) - 2 * (k - 1) * cdf_gamma(k - 1, a, b)
  res <- res + a * b * (2 * cdf_gamma(k - 1, a + 1, b) -
                          cdf_gamma(k - 2, a + 1, b) - cdf_gamma(k, a + 1, b))
  res <- vnapply(res, function(e) max(0, e))

  return(res)
}