Some test code for checking that the group labels as strings works.

library(SIBER)

mydata <- read.csv(file.path("../inst/extdata/test.group.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.95, 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)

Test overlap calculations

Maximum likelihood

testML <- maxLikOverlap("1.C", "1.D", 
                        test, 
                        p.interval = 0.95)

Bayesian overlap

testBayes <- bayesianOverlap("1.C", "1.D", 
                             ellipses.posterior, 
                             p.interval = 0.95, 
                             draws = 1000)

summary(testBayes)


AndrewLJackson/SIBER documentation built on Oct. 21, 2023, 8:09 a.m.