R/utils.R

Defines functions round_and_format nice_country_name cap_predictions ts_to_incid overlaps assign_epidemic_phase daily_to_weekly extract_predictions_qntls pool_predictions_weighted

Documented in cap_predictions overlaps ts_to_incid

##' @author Sangeeta Bhatia
##' @export
pool_predictions_weighted <- function(outputs, weights, nsim = 10000) {
  models <- names(weights)
  ## Sample model with weights
  n_1 <- sample(
    x = names(weights), nsim, replace = TRUE, prob = weights
  )
  n_1 <- table(n_1)
  if (! all(models %in% names(n_1))) {
    idx <- which(! models %in% names(n_1))
    n_1[[models[idx]]] <- 0
  }
  message("Number of times models picked ")
  message(paste(n_1, collapse = "\n"))

  out <- purrr::imap(
    outputs,
    function(output, model) {
      apply(output, 2, function(y) sample(y, size = n_1[[model]]))
    }
  )
  out <- Reduce("rbind", out)
  out
}

######################################################################
## Extract model quantiles
## Each model output is a list with several elements,
## the key one is Predictions.
## This is a list of data.frames with one data.frame for each country
## For each country, we have a list of 2 components, corresponding
## to the 2 serial intervals being considered.
## it is this last list (country, 2 components) that is passed to
## this function.
##' @export
extract_predictions_qntls <- function(y, prob = c(0.025, 0.25, 0.5, 0.75, 0.975)) {
  names(y) <- paste0("si_", seq_along(y))
  out <- purrr::map_dfr(
    y,
    function(y_si) {
      out2 <- t(
        apply(y_si, 2, stats::quantile, prob = prob, na.rm = TRUE)
      )
      out2 <- as.data.frame(out2)
      out2 <- tibble::rownames_to_column(out2, var = "date")
      out2
    },
    .id = "si"
  )
  out
}

##' @export
daily_to_weekly <- function(y, prob = c(0.025, 0.25, 0.5, 0.75, 0.975)) {
  names(y) <- paste0("si_", seq_along(y))
  out <- purrr::map_dfr(
    y,
    function(y_si) {
      weekly <- rowSums(y_si)
      weekly <- stats::quantile(
        weekly,
        prob = prob,
        na.rm = TRUE
      )
      weekly_df <- as.data.frame(weekly)
      ## This is not the last date for which predictions are
      ## available, but the last date for which observations are
      ## available.
      weekly_df$week_ending <- as.Date(colnames(y_si)[1]) - 1

      weekly_df <- tibble::rownames_to_column(
        weekly_df,
        var = "quantile"
      )

      weekly_df <- tidyr::spread(
        weekly_df, key = quantile, value = weekly
      )
      weekly_df
    },
    .id = "si"
  )
  out
}

##' @export
assign_epidemic_phase <- function(rt) {

  rt$phase <- dplyr::case_when(
    rt$`97.5%` < 1 ~ "decline",
    (rt$`97.5%` - rt$`2.5%` > 1)  ~ "unclear",
    (rt$`2.5%` > 1 &
     ((rt$`97.5%` - rt$`2.5%`) < 1))  ~ "growing",
    (rt$`2.5%` < 1 &
     ((rt$`97.5%` - rt$`2.5%`) < 1))  ~ "stable/growing slowly"
  )
  rt
}

##' Check if two intervals overlap
##'
##'
##'
##' @param x1 a numeric vector of length 2
##' @param x2 a numeric vector of length 2
##' @param digits integer indicating the number of decimal places
##' to be used
##' @return TRUE if x1 and x2 overlap
##' @author Sangeeta Bhatia
##' @export
overlaps <- function(x1, x2, digits) {

  low1 <- round(min(x1), digits = digits)
  high1 <- round(max(x1), digits = digits)
  interval1 <- seq(low1, high1, by = 10^(-digits))

  low2 <- round(min(x2), digits = digits)
  high2 <- round(max(x2), digits = digits)
  interval2 <- seq(low2, high2, by = 10^(-digits))

  common <- intersect(interval1, interval2)
  ## if the two intervals only overlap on the edge, we want to say
  ## they don't overlap.
  overlap <- length(common) > 1

  overlap
}



##' Convert time-series to incidence object
##'
##' incidence package accepts a linelist like object (list of dates)
##' and converts them into an incid object, which is needed by
##' projections package. To use projections package with a time-series
##' we convert time-series object into incid object.
##'
##' @param ts
##' @param date_col
##' @param case_col
##' @return
##' @author Sangeeta Bhatia
##' @export
ts_to_incid <- function(ts, date_col, case_col) {

  first_date <- min(ts[[date_col]])
  last_date <- max(ts[[date_col]])
  x <- tidyr::uncount(ts, weights = ts[[case_col]])
  out <- incidence::incidence(
    x[[date_col]],
    first_date = first_date,
    last_date = last_date
  )
  out
}

##' Cap predictions to twice the observed for visualisation
##'
##'
##'
##' @param pred
##' @return data.frame with capped predictions
##' @author Sangeeta Bhatia
##' @export
cap_predictions <- function(pred) {

  x <- split(pred, pred$country)
  purrr::map_dfr(x, function(y) {
    ymax <- 2 * ceiling(max(y$deaths, na.rm = TRUE) / 10) * 10
    y$`50%`[y$`50%` > ymax] <- NA
    dplyr::mutate_if(y, is.numeric, ~ ifelse(.x > ymax, ymax, .x))
  }
 )
}


##' @export
nice_country_name <- function(x) snakecase::to_title_case(as.character(x))

##' @export
round_and_format <- function(x, digits = 2) {
  format(round(x, digits), nsmall = digits)
}
mrc-ide/rincewind documentation built on Oct. 25, 2020, 12:32 p.m.