|||
同步更新于自媒体,完整版请查阅:新位置
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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 07:03
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社