Some test code for checking that the community (and group) labels as strings works.

library(SIBER)

library(viridis)
palette(viridis(5))


mydata <- read.csv("../inst/extdata/test.community.names.csv", 
                   header = TRUE)

test <- createSiberObject(mydata)

community.hulls.args <- list(col = 1, lty = 1, lwd = 1)
group.ellipses.args  <- list(n = 100, p.interval = 0.40, lty = 1, lwd = 2)
group.hull.args      <- list(lty = 2, col = "grey20")



par(mfrow=c(1,1))
plotSiberObject(test,
                ax.pad = 2, 
                hulls = T, community.hulls.args, 
                ellipses = T, group.ellipses.args,
                group.hulls = T, group.hull.args,
                bty = "L",
                iso.order = c(1,2),
                xlab = expression({delta}^13*C~'permille'),
                ylab = expression({delta}^15*N~'permille')
)

# Calculate sumamry statistics for each group: TA, SEA and SEAc
group.ML <- groupMetricsML(test)
print(group.ML)

# Calculate the various Layman metrics on each of the communities.
community.ML <- communityMetricsML(test) 
print(community.ML)

## ------------------------------------------------------------------------

# 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, parms, priors)


## ------------------------------------------------------------------------

# The posterior estimates of the ellipses for each group can be used to
# calculate the SEA.B for each group.
SEA.B <- siberEllipses(ellipses.posterior)

siberDensityPlot(SEA.B, xticklabels = colnames(group.ML), 
                xlab = c("Community | Group"),
                ylab = expression("Standard Ellipse Area " ('permille' ^2) ),
                bty = "L",
                las = 1,
                main = "SIBER ellipses on each group",
                ylims = c(0, 70)
                )

# Add red x's for the ML estimated SEA-c
points(1:ncol(SEA.B), group.ML[3,], col="red", pch = "x", lwd = 2)


AndrewLJackson/SIBER documentation built on June 7, 2024, 3:21 a.m.