匠人府分享 http://blog.sciencenet.cn/u/meiweipingg

博文

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

已有 5090 次阅读 2016-5-30 22:50 |个人分类:同位素|系统分类:科研笔记| 同位素, niche, SIBER, overlap, RStudio

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


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



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

上一篇:安装运行SIBER(in SIAR)时,如何解决 Error: not a matrix
下一篇:R统计分析之一:两独立样本检验及注意事项

0

该博文允许注册用户评论 请点击登录 评论 (1 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

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

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部