extractPredictiveData: Function to extract raw data from posterior predictive...

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

Function to extract raw data from posterior predictive model-checking distributions

Description

This function parses a data frame where it was stored the predicted data generated by the function TPmodel, and returns a list with all the predictive posterior data monitored.

Usage

extractPredictiveData(df = NULL, get = NULL, all = FALSE)

Arguments

df	data frame, variable where it was saved the output from a posterior predictive model-checking TPmodel run.
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).
all	logical, if TRUE it will return a combined list of all monitored data.

Value

a list of the data generated from a posterior predictive model-checking procedure, if all is TRUE, a vector is returned instead.

Examples

## Not run: 
# 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()

## End(Not run)

clquezada/tRophicPosition documentation built on Jan. 4, 2023, 12:31 p.m.