#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.