同步更新于自媒体，完整版请查阅：新位置

https://cran.r-project.org/web/packages/SIBER/vignettes/Introduction-to-SIBER.html

SIBER (SIBER package)是用来比较“两个生态系统之间”的同位素差异（overlap percentage）。

#--------step1. Setting up our R session for this demonstration-------------------------

rm(list=ls())

graphics.off()

library(SIBER)

mydata <- read.csv(fname, header=T)

siber.example <- createSiberObject(mydata)

#--------------------------step2.Plotting the raw data--------------------------------------

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(siber.example,

ax.pad = 2,

hulls = F, 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~'u2030'),

ylab = expression({delta}^15*N~'u2030')

)

#------------------step3.Summary statistics and custom graphic additions------------------

par(mfrow=c(1,1))

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")

plotSiberObject(siber.example,

ax.pad = 2,

hulls = F, community.hulls.args,

ellipses = F, group.ellipses.args,

group.hulls = F, group.hull.args,

bty = "L",

iso.order = c(1,2),

xlab=expression({delta}^13*C~'u2030'),

ylab=expression({delta}^15*N~'u2030'),

cex = 0.5

)

group.ML <- groupMetricsML(siber.example)

print(group.ML)

plotGroupEllipses(siber.example, n = 100, p.interval = 0.95,

lty = 1, lwd = 2)

plotGroupEllipses(siber.example, n = 100, p.interval = 0.95, ci.mean = T,

lty = 1, lwd = 2)

plotSiberObject(siber.example,

ax.pad = 2,

hulls = T, community.hulls.args,

ellipses = F, group.ellipses.args,

group.hulls = F, group.hull.args,

bty = "L",

iso.order = c(1,2),

xlab=expression({delta}^13*C~'u2030'),

ylab=expression({delta}^15*N~'u2030'),

cex = 0.5

)

plotGroupEllipses(siber.example, n = 100, p.interval = 0.95,

ci.mean = T, lty = 1, lwd = 2)

# Calculate the various Layman metrics on each of the communities.

community.ML <- communityMetricsML(siber.example)

print(community.ML)

#-------------------step4.Fitting the Bayesian models to the data-------------------------

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

priors <- list()

priors$R <- 1 * diag(2)

priors$k <- 2

priors$tau.mu <- 1.0E-3

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

#------------------step5.Comparing among groups using the Standard Ellipse Area-----------

SEA.B <- siberEllipses(ellipses.posterior)

siberDensityPlot(SEA.B, xticklabels = colnames(group.ML),

xlab = c("Community | Group"),

ylab = expression("Standard Ellipse Area " ('\u2030' ^2) ),

bty = "L",

las = 1,

main = "SIBER ellipses on each group"

)

points(1:ncol(SEA.B), group.ML[3,], col="red", pch = "x", lwd = 2)

cr.p <- c(0.95, 0.99) # vector of quantiles

SEA.B.credibles <- lapply(

as.data.frame(SEA.B),

function(x,...){tmp<-hdrcde::hdr(x)$hdr},

prob = cr.p)

SEA.B.modes <- lapply(

as.data.frame(SEA.B),

function(x,...){tmp<-hdrcde::hdr(x)$mode},

prob = cr.p, all.modes=T)

#--------------step6.Comparison of entire communities using the Layman metrics------------

mu.post <- extractPosteriorMeans(siber.example, ellipses.posterior)

layman.B <- bayesianLayman(mu.post)

siberDensityPlot(layman.B[[1]], xticklabels = colnames(layman.B[[1]]),

bty="L", ylim = c(0,20))

comm1.layman.ml <- laymanMetrics(siber.example$ML.mu[[1]][1,1,],

siber.example$ML.mu[[1]][1,2,]

)

points(1:6, comm1.layman.ml$metrics, col = "red", pch = "x", lwd = 2)

siberDensityPlot(layman.B[[2]], xticklabels = colnames(layman.B[[2]]),

bty="L", ylim = c(0,20))

相关博文：

在R中运行        SIAR：http://blog.sciencenet.cn/blog-651374-1006886.html  

在R中运行      SIBER（SIBER package）： http://blog.sciencenet.cn/blog-651374-981332.html  

在R中运行      SIBER（SIAR package） http://blog.sciencenet.cn/blog-651374-979988.html