R/extractPredictiveData.R

In AndrewLJackson/tRophicPosition: Bayesian Trophic Position Estimation with Stable Isotopes

#' Function to extract raw data from posterior predictive model-checking
#'  distributions
#'
#' This function parses a data frame where it was stored the predicted data
#' generated by the function \code{\link{TPmodel}}, and returns a list with
#' all the predictive posterior data monitored.
#'
#' @param df data frame, variable where it was saved the output from a posterior
#' predictive model-checking \code{\link{TPmodel}} run.
#' @param get string, could be either "dNcPred" and "dCcPred" (for consumer data
#' ), "dNb1Pred" and "dCb1Pred" (for baseline 1 data), or "dNb2Pred" and
#' "dNC2Pred" (for baseline 2 data).
#' @param all logical, if TRUE it will return a combined list of all monitored
#' data.
#'
#' @return a list of the data generated from a posterior predictive
#' model-checking procedure, if all is TRUE, a vector is returned instead.
#' @export
#'
#' @examples
#' \dontrun{
#' # Generate isotope data
#' isotopeData <- generateTPData(n.obsB = 45,
#'                               sd.dNb1 = 1, sd.dCb1 = 1,
#'                               dCb1 = -34, dCb2 = -24,
#'                               n.obsC = 60, dCc = -29)
#'
#' # Define a one baseline model
#' model.string <- jagsOneBaseline()
#'
#' # Initialize the model
#' model_d <- TPmodel(data = isotopeData,
#'                    model.string = model.string,
#'                    n.chains = 2,
#'                    n.adapt = 5000)
#'
#' # Generate posterior samples of TP, alpha and dNcPred
#' # dNcPred stands for the predicted data dNc (nitrogen values of consumer)
#' samples_d <- posteriorTP(model_d,
#'                          variable.names = c("TP", "alpha", "dNcPred"),
#'                          burnin = 5000,
#'                          n.iter = 5000,
#'                          thin = 10)
#'
#' # Extract posterior predictive data
#' dNcPred <- extractPredictiveData(samples_predicted, get = "dNcPred")
#'
#' # Calculate residuals
#' dNcPred_res <- sweep(dNcPred, 2, isotopeData$dNc, "-")
#'
#' # Combine all residuals
#' dNcPred_resall <- as.numeric(do.call(rbind, dNcPred_res))
#'
#' # Plot a sample of them
#' plot(sample(dNcPred_resall, 10000), xlab = "Index",
#'      ylab = "Residuals")
#'
#' # Extract posterior predictive data combined
#' dNcPred_all <- extractPredictiveData(samples_predicted, get = "dNcPred",
#'                                      all = TRUE)
#'
#' # Plot a histogram of observed data and a density function of predicted data
#' # We need ggplot2 installed previously
#' ggplot2::ggplot() +
#'   ggplot2::geom_histogram(ggplot2::aes(x = a$dNc, y = ..density..),
#'                           binwidth = 0.3, fill = "grey", color = "black") +
#'   ggplot2::geom_density(ggplot2::aes(x = dNcPred_all), color = "blue") +
#'   ggplot2::xlab(expression(paste(delta^{15}, "N (\u2030)"))) +
#'   ggplot2::theme_bw()
#' }
#'
extractPredictiveData <- function (df = NULL, get = NULL, all = FALSE) {

  combined <- as.data.frame(coda::mcmc(do.call(rbind, df)))

  predicted <- combined[data.table::like(names(combined), get)]

  if (isTRUE(all))
    return(as.numeric(do.call(rbind, predicted)))
  else return(predicted)

}