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

2017-03-15

(本文于 2016-5-30 22:50 首发于 “科学网”) 2018/02/01 updated

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

各原始数据在csv中的输入格式参考下图（相关说明参考这里）：

step1. 调入数据. Setting up our R session for this demonstration

rm(list=ls())
graphics.off()
library(SIBER)
mydata <- read.csv("C:/~/Desktop/mydata-SIBERinSIBER.csv",, header=T)
siber.example <- createSiberObject(mydata)
#Plotting the raw data

step2.初绘散点分布图. Plotting the raw data

Create lists of plotting arguments to be passed onwards to each of the three plotting functions.

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

# this time we will make the points a bit smaller by
# cex = 0.5

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
)

Calculate sumamry statistics for each group: TA, SEA and SEAc

1 2	group.ML <- groupMetricsML(siber.example) print(group.ML)

You can add more ellipses by directly calling plot.group.ellipses()

Add an additional p.interval % prediction ellilpse.

1 2	plotGroupEllipses(siber.example, n = 100, p.interval = 0.95, lty = 1, lwd = 2)

Or you can add the XX% confidence interval around the bivariate means by specifying ci.mean = T along with whatever p.interval you want.

1	plotGroupEllipses(siber.example, n = 100, p.interval = 0.95, ci.mean = T, lty = 1, lwd = 2)

A second plot provides information more suitable to comparing the two communities based on the community-level Layman metrics.

This time we will make the points a bit smaller by cex = 0.5

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
)

Or you can add the XX% confidence interval around the bivariate means by specifying ci.mean = T along with whatever p.interval you want.

1	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.

1 2	community.ML <- communityMetricsML(siber.example) print(community.ML)

step4.贝叶斯模型拟合. Fitting the Bayesian models to the data

# 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(siber.example, parms, priors)

step5.Groups 间的标准椭圆面积的比较.Comparing among groups using the Standard Ellipse Area

# 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 " ('\u2030' ^2) ),
bty = "L",
las = 1,
main = "SIBER ellipses on each group"
)

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

# Calculate some credible intervals

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

# call to hdrcde:hdr using lapply()

SEA.B.credibles <- lapply(
as.data.frame(SEA.B),
function(x,...){tmp<-hdrcde::hdr(x)$hdr},
prob = cr.p)

# do similar to get the modes, taking care to pick up multimodal posterior distributions if present

SEA.B.modes <- lapply(
as.data.frame(SEA.B),
function(x,...){tmp<-hdrcde::hdr(x)$mode},
prob = cr.p, all.modes=T)

step6.基于Layman矩阵的全生态系统之间的比较.Comparison of entire communities using the Layman metrics

# extract the posterior means

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

# calculate the corresponding distribution of layman metrics

layman.B <- bayesianLayman(mu.post)

# Visualise the first community

siberDensityPlot(layman.B[[1]], xticklabels = colnames(layman.B[[1]]),
bty="L", ylim = c(0,20))

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

# Visualise the second community


siberDensityPlot(layman.B[[2]], xticklabels = colnames(layman.B[[2]]),
bty="L", ylim = c(0,20))


# Alternatively, pull out TA from both and aggregate them into a
# single matrix using cbind() and plot them together on one graph.

Q & A

1. 在 SIBER package 中运行 SIBER 程序时, 数据集必须包括 4 列内容： iso1,iso2,group,community.

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

<已有次阅读>

由于本文作者水平有限，文中如有错误之处，欢迎大家批评指正！

① 本文仅代表作者个人观点，不代表任何其它立场，欢迎交流合作！

② 转载与分享请注明：本文源于 http://meiweiping.cn

^_^ 本文字数统计：1.1k 字

梅卫平の阅览室

科研同位素

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