🔍 检索绿色通道
     1.JCR-IF查询      2.CiteScore-IF查询     3.SCI检索     4.中科院分区查询     5.SSCI检索     6.AHCI检索     7.ESCI检索(新增)     8.EI检索     9.EI检索备用     10.ISTP/CPCI检索备用(新增)      11.梅斯期刊智能查询     12.CSCD检索     13.CSSCI检索     14.北大核心检索      15.国自然官方查询     16.国自然梅斯查询     17.国自然科学网查询     18.论文被检索情况查询方法

SIBER(SIBER package)in R:多生态系统间同位素生态位差异比较


(本文于 2016-5-30 22:50 首发于 “科学网”)

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


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

1
2
3
4
5
6
rm(list=ls())
graphics.off()
library(SIBER)
mydata <- read.csv(fname, 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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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

1
2
3
4
5
6
7
8
9
10
11
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# 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

(>看完记得五星好评哦亲<)

分享到:
    友情链接
    e教学: 蓝墨云班课     传课网     慕课网     作业帮
    BBS: 小木虫     丁香园     科学网bbs     零点花园     数据狗论坛     ResearchGate
    预印本: 中国科技论文在线(教育部)     ChinaXiv(中科院)     bioRxiv(冷泉港)     arXiv(康奈尔大学)
    科研资讯: 科学网   科学之家     51science     Publons     J-STAGE     日本の研究
    信息查询: 中国研究生招生信息网     学信网(学历查询)     中国人事考试网     更多>>
    网站建设: Hexo     极简图床     sm.ms     html-unicode     UTF-8转码     颜色代码
    关于本站: 关于我们      Logo含义      广告合作      免责声明      联系我们      RSS订阅