R/posterior_predict.pandemicEstimated.R

In PandemicLP: Long Term Prediction for Epidemic and Pandemic Data

# Environment created to exchange objects between package functions
pandemic_environment = new.env()
#pandemic_environment$fullPred = list()

#' Draw from the posterior predictive distribution for pandemic data
#'
#' The posterior predictive distribution is the distribution of the outcome
#' implied by the model after using the observed data to update our beliefs
#' about the unknown parameters in the model. Simulating data from the posterior
#' predictive distribution using the observed predictors is useful for checking
#' the fit of the model. Drawing from the posterior predictive distribution at
#' interesting values of the predictors also lets us visualize how a manipulation
#' of a predictor affects (a function of) the outcome(s). With new observations of
#' predictor variables we can use the posterior predictive distribution to generate
#' predicted outcomes.
#' @method posterior_predict pandemicEstimated
#' @param object An object of class \code{pandemicEstimated} created by function \code{\link{pandemic_model}}.
#' @param horizonLong How far into the future the long-term prediction is desired.
#' @param horizonShort How far into the future the short-term prediction is desired.
#' @param ... Currently unused.
#' @references
#' CovidLP Team, 2020. CovidLP: Short and Long-term Prediction for COVID-19. Departamento de Estatistica. UFMG,
#' Brazil. URL: \url{http://est.ufmg.br/covidlp/home/en/}
#' @return An object of class \code{pandemicPredicted}. It includes the sampled predictive distribution
#' the model used to predict, which is the same as the one used to estimate the data. This object can be used
#' directly into the plot function and contains the following elements:
#' \item{\code{predictive_Long}}{
#'   A \code{M x horizonLong} matrix with the full sample of the predictive distribution
#'   for the long-term prediction, where M is the sample size.
#'   The prediction is for daily new cases.
#'   }
#'   \item{\code{predictive_Short}}{
#'   A \code{M x horizonShort} matrix with the full sample of the predictive distribution
#'   for the short-term prediction, where M is the sample size.
#'   The prediction is for daily cumulative cases.
#'   }
#'   \item{\code{data}}{
#'   The data passed on from the \code{\link{pandemicEstimated-objects}} under the element \code{Y$data}.
#'   }
#'   \item{\code{location}}{
#'   A string with the name of the location.
#'   }
#'   \item{\code{cases_type}}{
#'   A string with either "confirmed" or "deaths" to represent the type of data that has been fitted and predicted.
#'   }
#'   \item{\code{pastMu}}{
#'   The fitted means of the data for the observed data points.
#'   }
#'   \item{\code{futMu}}{
#'   The predicted means of the data for the predicted data points.
#'   }
#' Function \code{\link{pandemic_stats}} provides a few useful statistics based on the predictions.
#'
#' @seealso \code{\link{pandemic_model}}, \code{\link{pandemic_stats}} and \code{\link{plot.pandemicPredicted}}. Details
#' about the models behind the calculations can be seen in \code{\link{models}}.
#' @examples
#' \dontrun{
#' dataMG = load_covid("Brazil","MG")
#' estimMG = pandemic_model(dataMG)
#' predMG = posterior_predict(estimMG)
#' predMG}
#' @importFrom rstantools posterior_predict
#' @importFrom methods slot
#' @export
posterior_predict.pandemicEstimated = function(object, horizonLong = 500, horizonShort = 14, ...){

  if (!is(object, "pandemicEstimated")) stop("Please use the output of the pandemic_model() function.")
  if (horizonShort <= 0) stop("Horizons must be positive.")
  if (horizonLong < horizonShort) stop("Long-term horizon may not be lesser than short term horizon.")
  #if (horizonLong > 1000) stop("Long-term horizon larger than 1000 is not currently supported.")

  chains <- as.data.frame(object$fit)
  pop <- object$Y$population
  longPred <- max(1000, horizonLong)

  M = nrow(chains) ## Total iterations
  NA_replacement = 2 * object$Y$population ## Set a NA replacement

  finalTime = sum(grepl("mu", names(chains))) ## How many mu's

  # generate points from the marginal predictive distribution
  if (length(object$seasonal_effect)){
    ss_code = c(seasonal_code(object$Y$data$date, object$seasonal_effect), 0, 0)
    s_code = list(s1 = rep(ss_code[1], object$n_waves),
                  s2 = rep(ss_code[2], object$n_waves),
                  s3 = rep(ss_code[3], object$n_waves))
  }
     else s_code = list(s1 = rep(0, object$n_waves),
                        s2 = rep(0, object$n_waves),
                        s3 = rep(0, object$n_waves))
  pred = generatePredictedPoints_pandemicC(M, chains, longPred, NA_replacement,
                                           object$model_name, finalTime,
                                           object$n_waves, s_code)
  methods::slot(object$fit,"sim")$fullPred = list() # For internal use
  methods::slot(object$fit,"sim")$fullPred$thousandLongPred = pred$yL
  if (object$cases.type == "confirmed")
    methods::slot(object$fit,"sim")$fullPred$thousandShortPred = pred$yS + object$Y$data$cases[nrow(object$Y$data)]
  else
    methods::slot(object$fit,"sim")$fullPred$thousandShortPred = pred$yS + object$Y$data$deaths[nrow(object$Y$data)]
  methods::slot(object$fit,"sim")$fullPred$thousandMus = pred$mu
  y.futL = as.matrix(methods::slot(object$fit,"sim")$fullPred$thousandLongPred[,1:horizonLong])
  y.futS = as.matrix(methods::slot(object$fit,"sim")$fullPred$thousandShortPred[,1:horizonShort])
  errorCheck = which(methods::slot(object$fit,"sim")$fullPred$thousandShortPred[,ncol(methods::slot(object$fit,"sim")$fullPred$thousandShortPred)] > pop)
  if (length(errorCheck)){
    message(paste0(length(errorCheck)," samples were removed from the prediction due to unrealistic results."))
    methods::slot(object$fit,"sim")$fullPred$thousandShortPred = methods::slot(object$fit,"sim")$fullPred$thousandShortPred[-errorCheck,]
    methods::slot(object$fit,"sim")$fullPred$thousandLongPred = methods::slot(object$fit,"sim")$fullPred$thousandLongPred[-errorCheck,]
    y.futL = y.futL[-errorCheck,]
    y.futS = y.futS[-errorCheck,]
  }

  output <- list(predictive_Long = y.futL, predictive_Short = y.futS,
                 data = object$Y$data, location = object$Y$name, cases_type = object$cases.type,
                 pastMu = as.data.frame(object$fit)[grep("mu",names(chains))],
                 futMu = as.matrix(pred$mu[,1:horizonLong]),fit = object$fit,seasonal_effect = object$seasonal_effect,
                 errors = errorCheck)

  class(output) = "pandemicPredicted"
  return(output)
}