These are some extra functions I wrote to help process GJAM output.

source("gjamExtras.r")
library(tidyverse)

It’s important to run each model multiple times to make sure the chains converge. Once you have multiple runs of your model, use gelman() to determine which parameter chains have converged and combineChains() to combine the chains from multiple runs into one parameter estimate.

The typical cut off for Gelman indices is 1.2. Higher values indicate that the chains did not converge.

outNames <- c("run1.rdata", "run2.rdata", "run3.rdata")
gel <- gelman(outNames = outNames)

converged <- gel[which(gel$gelman < 1.2),]


# Visualize gelman indices by species
ggplot(gel) +
  geom_hline(yintercept = 0) +
  geom_histogram(aes(x = gelman), binwidth = 0.2) +
  facet_wrap(~ species, scales = "free_x") +
  geom_vline(xintercept = 1.2, linetype = "dashed") +
  theme_bw() +
  theme(strip.background = element_blank()) +
  labs(x = "Gelman-Ruben Index")

# Visualize gelman indices by predictor
ggplot(gel) +
  geom_hline(yintercept = 0) +
  geom_histogram(aes(x = gelman), binwidth = 0.2) +
  facet_wrap(~ predictor, scales = "free_x") +
  geom_vline(xintercept = 1.2, linetype = "dashed") +
  theme_bw() +
  theme(strip.background = element_blank()) +
  labs(x = "Gelman-Ruben Index")

Use plotChains() to see the chains. You may want to plot the chains for a certain species or predictor based on the gelman indices. By default, all predictors for 5 random species are plotted. Use the plotPred, plotSpec, and numChains parameters to select certain species or predictors. This function only plots chains from one output at a time.

plotChains(out, plotPred = "all", plotSpec = "random")

Combine chains from all runs and calculate point estimate and credible interval. credInt can be 95 or 90.

# combine chains from all runs
outNames <- c("OUT/run1.rdata", "OUT/run2.rdata", "OUT/run3.rdata")
b_combined <- combineChains(outNames = outNames, credInt = "95")

Remove parameters with gelman indices greater than 1.2.

b <- inner_join(b_combined, converged, by = c("species", "predictor"))

Plot all the betas using plotBetas(). To color species by trait, include a dataframe with two columns named species and trait.

plotBetas(b)

Finally, plot observed and predicted values for each species. This function only uses one model run.

plotPredObs(out)