# SIBER (SIAR package) in R：物种生态位重复度计算

(本文于 2016-5-27 00:10 首发于 “科学网”)
2018/02/01 updated 2017/10/05 updated 2017/3/24 updated 温馨提示：请确保您的 R软件（及R相关软件）安装的系统盘（e.g., C盘），这点十分重要！

SIBER: Stable Isotope Bayesian Ellipses in R

“SIBER” (SIAR package)这个SIAR套件的SIBER分析，主要是用来比较一个生态系统内不同种群之间的同位素生态位竞争的具体百分比值，等等。

Fits bivariate ellipses to stable isotope data using Bayesian inference with the aim being to describe and compare their isotopic niche.

rm(list=ls())
setwd("E:/")
graphics.off()
# 我的数据集包含 3 列: group,x,y.
# x, y 表示 d13C 和 d15N 的值.


library(siar)
attach(mydata)
# now loop through the data and calculate the ellipses.
ngroups<-length(unique(group))

#split the isotope data based on group.
spx<-split(x,group)
spy<-split(y,group)

#create some empty vectors for recording our metrics.
SEA<-numeric(ngroups)
SEAc<-numeric(ngroups)
TA<-numeric(ngroups)

# plot the raw data.
# dev.new() #maybe no need show the plot in a new window.
plot(x,y,col=group,type="p")
legend("topright",legend=as.character(paste("Group",unique(group))),pch=19,col=1:length(unique(group)))


Now loop over each group and plot the small sample size corrected ellipse and the convex hull for each.

for(j in unique(group)){

#fit a standard ellipse to the data.
SE<-standard.ellipse(spx[[j]],spy[[j]],steps=1)

#extract the estimated SEA and SEAc from this object.
SEA[j]<-SE$SEA SEAc[j]<-SE$SEAc

#plot the stand ellipse with d.f.=2(i.e.SEAc).
#these are plotted here as thick solid lines.
lines(SE$xSEAc,SE$ySEAc,col=j,lty=1,lwd=3)   ###

#also,for comparison we can fit and plot the convex hull.
#the convex hull is plotted as dotted thin lines.
# calculate the convex hull for the jth group's isotope values.
# held in the objects created using split()called spx and spy.
CH<-convexhull(spx[[j]],spy[[j]])

#extract the area of the convex hull from this object.
TA[j]<-CH$TA #plot the convex hull lines(CH$xcoords,CH\$ycoords,lwd=1,lty=3)     ###
}


print the area metrics to screen for comparison.

#NB if you are working with real data rather than simulated then you wont be able to calculate the population SEA(pop.SEA).

#if you do this enough times or for enough groups you will easily see the bias in SEA as an estimate of pop.SEA as compared to SEAc which is unbiased.

#both measures are equally variable.

print(cbind(SEA,SEAc,TA))


So for we have fitted the standard ellipses based on frequentist methods and calculated the relevant metrics(SEA and SEAc).

Now we turn our attention to producing a Bayesian estimate of the standard ellipse and its area SEA.B.

reps<-10^4  #循环计算次数.

# Generate the Bayesian estimates for the SEA for each group using the utility function siber.ellipses.

SEA.B<-siber.ellipses(x,y,group,R=reps)
summary(SEA.B)    #N.B. SEA.B数据集包含  10^4 行的内容.


plot out some of the data and results.

Compare two ellipses for significant differences in SEA.

To test whether Group 1 SEA is smaller than Group 2…

You need to calculate the Proportion of G1 ellipses that are less than G2.

To calculate the overlap between two ellipses you can use the following code

One might then wish to calculate the proportion overlap; athough one then runs into a choice as to what the demoninator will be in the equation. You could for instance calculate

the proportion of A that overlaps with B,

the proporiton of B that overlaps with A,

or the proportion of A and B that overlap with each other.

You can also cacluate the overlap between two of the convex hulls, or indeed any polygon using the code that underlies the overlap() function.

And calculate the overlap using the function in spatstat package.

Q & A

For this case from siar package, data only need 3 columns as iso1,iso2,group.

<已有 次阅读>

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

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

（>看完记得五星好评哦亲<）

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