# 多生态系统间同位素生态位差异比较： SIBER（SIBER）in R

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)

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,

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,

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,

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

http://wap.sciencenet.cn/blog-651374-981332.html

## 相关博文

GMT+8, 2021-1-25 19:43