siberMVN throws error after adding 3 more lake communities...ideas?

[megredgar]

Hello!

I am trying to run code to look at the isotopic food webs for five different reservoir communities. My code works when only running two specific communities, but the addition of the three other reservoirs causes issues when running JAGS (Error in update.jags(model, n.iter, ...) : RUNTIME ERROR: Failed to get Cholesky decomposition of R).

*This error only happens when attempting to run the data for isotopic food webs... there are no issues when running the model to assess the isotopic niche.

I've already tried reducing the number of iterations and numbering communities/groups instead of having them as text.

Any help or tips would be greatly appreciated!

Thank you!

Here's my code for the model:

parms <- list() parms$n.iter <- 1 10^4 parms$n.burnin <- 1 10^3 parms$n.thin <- 10
parms$n.chains <- 2

define the priors priors <- list() priors$R <- 1 * diag(2) priors$k <- 2 priors$tau.mu <- 1.0E-3

library(rjags)

ellipses.posterior <- siberMVN(siber.reservoirs, parms, priors)

[AndrewLJackson]

hi. Ive tested SIBER with more than two communities and it works. see code below. This error suggests at least one of your groups is causing problems for the estimation of the covariance matrix. Check the output from groupMetricsML() which might indicate if the "standard" way to calculate a covariance matrix is working on all your groups. Low sample size (n<3) and/or data that are replicates or perfectly collinear might cause problems.

best wishes Andrew

---

title: "Test SIBER on more than two communities" output: html_notebook

library(SIBER)

Simulate some data

test_data <- generateSiberData(n.groups = 3, 
                               n.communities = 3,
                               n.obs = 30,
                               mu.range = c(-1,1,-1,1),
                               wishSigmaScale = 1)

Create the SIBER object

# create the SIBER class object
test_siber_obj <- createSiberObject(test_data)

# test it by calculating max likelihood summary statistics
test_groupML <- groupMetricsML(test_siber_obj)

print(test_groupML)

Fit the Bayesian ellipses

# options for running jags
parms <- list()
parms$n.iter <- 2 * 10^4   # number of iterations to run the model for
parms$n.burnin <- 1 * 10^3 # discard the first set of values
parms$n.thin <- 10     # thin the posterior by this many
parms$n.chains <- 2        # run this many chains

# define the priors
priors <- list()
priors$R <- 1 * diag(2)
priors$k <- 2
priors$tau.mu <- 1.0E-3

# fit the ellipses which uses an Inverse Wishart prior
# on the covariance matrix Sigma, and a vague normal prior on the 
# means. Fitting is via the JAGS method.
ellipses.posterior <- siberMVN(test_siber_obj, 
                               parms, 
                               priors)

Calculate the community level "layman" metrics.

# extract the posterior means
mu.post <- extractPosteriorMeans(test_siber_obj,
                                 ellipses.posterior)

# calculate the corresponding distribution of layman metrics
layman.B <- bayesianLayman(mu.post)

Plot the first community to check its working

# --------------------------------------
# Visualise the first community
# --------------------------------------
siberDensityPlot(layman.B[[1]], xticklabels = colnames(layman.B[[1]]), 
                bty="L", ylim = c(0,3))

# add the ML estimates (if you want). Extract the correct means 
# from the appropriate array held within the overall array of means.
comm1.layman.ml <- laymanMetrics(
  test_siber_obj$ML.mu[[1]][1,1,],
  test_siber_obj$ML.mu[[1]][1,2,]
                                 )
points(1:6, comm1.layman.ml$metrics, col = "red", pch = "x", lwd = 2)

[megredgar]

@AndrewLJackson Thank you so much! This helped me figure out the solution to my problem. When I just have two communities, sample sizes less than 3 don't cause issues. But, with more than 2 communities, sample sizes less than 3 cause problems.

Many thanks! Meg