(本文于2016-10-5 20:29 首发于 “科学网”)

在稳定同位素领域,20世纪70 80 年代开始,利用碳氮稳定同位素来研究食物网结构成为有力的技术手段,以弥补传统的胃含物方法中较难分析的空胃率较高的生物的缺陷。

在食物网研究中,主要包括食物来源和营养级的分析。其中,食物来源分析中严格来讲,1种同位素只能对应2种食物来源的分析,n种同位素相对应n+1种食物来源的分析。

问题是,例如有5种食物来源,只用碳、氮2种同位素进行食物来源分析就很困难,20世纪后半叶,普遍采用的方法是Isosource软件进行计算。但是这种方法:①只使用饵料生物同位素的平均值,未考虑标准差的影响;②只使用固定浓缩系数(TEF)值(e.g., 碳氮同位素分别使用1‰ 和3.4‰),也未考虑其标准差的影响等(Parnell AC, Inger R, Bearhop S, et al.,2010)。

本文简要介绍基于贝叶斯算法的最新方法: siar(stable isotope analysis in R)在R软件上的操作过程。

补充R语言基础知识:

上述图片中的 蓝色字体内容均为R代码,并非csv文件里面的内容。上述图片包括3个独立的.csv文件,分别为sources.csv,corrections.csv和consumers.csv文件。而不是都输入在一个.csv文件里。

①上述图片最左一列(1~3,1~19)为excel或csv软件的侧栏编号,请勿写入excel或csv文件的第一列。即,本例中:csv文件的第一列标题为Source或Group。

②R语言软件自身缺点所致,在csv文件数据等所有内容最好不要用宋体等中文字体,建议同意采用新罗马字体,否则容易识别出错或乱码。

③ 写好的csv文件在保存时,请务必选择 CSV(MS-DOS)或苹果电脑选择 CSV(Macintosh)后缀保存,而不要选择其他的CSV后缀保存。1

2

3

4

5

6library(siar)

graphics.off()

sources

consumers

corrections

concs

step2: 作图如下

1

2

3model1

siarplotdata(model1,iso=c(2,1),leg = 1)R中直接输出图之后,将图片copy to clipboard asMetafile 到 ppt中,然后选择该图片并点击ppt工具栏中的“取消组合”,可以对图例等内容进行自由编辑。补充内容之 如何计算各饵料贡献率的标准偏差(Sd值)等参数? 上述模型“model1”建立之后,其实所有贡献率结果已经出来了,我们需要先看一下结果里面到底暗藏了什么玄机,在R软件中输入代码如下1head(model1$output) 得到结果如下图(此处图略) 此时,我们就已经知道model1$output各列代表什么内容了,接下来代码十分重要,务必仔细观察理解并应用!!!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

33

34

35summary(model1$output[,1])

sd(model1$output[,1])

#上述两步即可得到我们想要的 src1对应consumer1的贡献率的各个参数(mean,sd等)数据了,同样地,

summary(model1$output[,2])

sd(model1$output[,2]) # src 2 对应 consumer 1

summary(model1$output[,3])

sd(model1$output[,3]) # src 3 对应 consumer 1

#----

summary(model1$output[,6])

sd(model1$output[,6])

#以上为 src1 对应 consumer2 的贡献率的各个参数(mean,sd等)数据

#注意① 请参照 head(model1$ooutput中的表头内容,相应地更改上述代码中括号内[,#]的数字)

#注意② 和①相同,src 种类不同,上述中括号内的数字会发生变化,请注意!

summary(model1$output[,7])

sd(model1$output[,7]) # src 2 对应 consumer 2

summary(model1$output[,8])

sd(model1$output[,8]) # src 3 对应 consumer 2

#---

summary(model1$output[,11])

sd(model1$output[,11]) # src 1 对应 consumer 3

summary(model1$output[,12])

sd(model1$output[,12]) # src 2 对应 consumer 3

summary(model1$output[,13])

sd(model1$output[,13]) # src 3 对应 consumer 3

# etc... model1$output 数据集中的 “SD1G1,SD2G1” 等其他参数 表达的意思是什么,请大家自主解决这个问题,我就不一一答复了!

step3:作饵料来源所占捕食者食物中的比例的分布图

1

2

3

4siarhistograms(model1)

# 然后在R中输入 1:4, Selection:2,即为上图,其中4表示group4,2表示将所有饵料生物占比在同一张图中表示;同理,

siarhistograms(model1)

# 然后在R中输入 1:1, Selection:2,即为group1的食物来源占比分布图;以此类推。

step4:分布图形式2,以捕食者的group为变量

1

2siarproportionbygroupplot(model1,grp=4)

# 代码中的grp=4表示group4;以此类推。 上述代码运行后,会让你填一些数字后按回车,第一步说明 1:1 for Group1,第二步说明Selection:2 for all source proportion together on one graph.新增内容 2017/03/26 updated

自定义箱线图之利用SIAR计算不同饵料来源对消费者的贡献比例

上述 Step4 我们已经绘制了所有饵料对消费者(以group4为例)的贡献率箱线图,SIAR的默认设置中

siarproportionbygroupplot()函数只绘制贡献率的50%,75%和95%这三个范围,

想要自定义不同的的范围,只需要将上述Step4中的代码1

2siarproportionbygroupplot(model1,grp=4)

# 代码中的grp=4表示group4;以此类推。 更改为如下代码:1

2siarproportionbygroupplot(model1,grp=4,probs = c(95,75,50,25,5))

# 代码中的grp=4表示group4;以此类推。 即可。 附上SIAR工具包原作者AndrewL Jackson对siarproportionbygroupplot()函数的源代码请仔细参考 https://github.com/AndrewLJackson/siar/blob/master/R/siarproportionbygroupplot.R

# step5:分布图形式3,以饵料生物src的种类为变量 #

代码为:

> siarproportionbysourceplot(model,grp=2)

然后在R中按回车键几个显示上图。代码中的grp=2表示饵料生物src2,以此类推。

#step6:导出数据化的结果,见下图左半部分#

代码如下:

1

2siarhdrs(model1)

#src1G1表示第1种饵料生物对group1的捕食者食物来源的贡献比例结果,包括占比的平均值和95%置信区间等。以此类推。

以上。

# siar Q & A #

#Q[1]#

梅老师,您好:

看了你关于SIAR求饵料贡献率的代码,有几个问题想问一下,能不能麻烦您帮我解决一下呢?

1.您所给的数据能具体解释一下是怎么来的吗?是几种饵料呢?(我是初学者,对数据的编辑整理还不是很明白,麻烦您知道一下呢)如果饵料只测了一次而没有标准差该怎么办呢?那个corrections代表什么意思呢?

2.如果我们要测贡献率,那那个直方图该怎么解释呢?

万分感谢~

A: 您好

1. 饵料最好5种,本文例子中为3种:src 1,src2,src3

2. 饵料对应捕食者数目,某一种捕食者比如一种鲫鱼,它也应该有很多个平行样本才有意义,实际情况中这个平行可以是多个采样点

3. corrections表示src对应consumer的同位素的浓缩系数

4. 饵料贡献率直接使用#step6中的数据即可。#step2下面的图中,src的数据是已经加上浓缩系数校正后的值,consumer为原始值。

#Q[2]#

梅老师,您好:

感谢您百忙之中抽空解答我的问题,您说corrections表示src对应consumer的同位素的浓缩系数,我看到您所给的值是一样的,15N是3.54,13C是1.63,能否请问一下您这个浓缩系数值是怎么取值的吗?对应的浓缩系数里面的Sd标准差又是如何得到的呢?

A:

1. 本文只是举例而已,并非实际数值。

2. 不同捕食者对应相同饵料实际情况下浓缩系数肯定不完全相同。

3. 同种捕食者同种饵料,在环境条件不同的情况下浓缩系数也不完全相同,抑或同种捕食者在不同生长时期的实验室饲养情况下浓缩系数也会有差异,就会有sd值了。

4. 以上为已有各个实验参数的理想条件下的sair分析,但是,比如未知特定捕食者对应的饵料的浓缩系数的情况下,用其它已知的相近的捕食者的浓缩系数代替,或者用其它综述论文的浓缩系数平均值代替也一定程度上认为是无可厚非的。

#Q[3]#

梅老师,您好,又要麻烦您啦,我想问一下那个消费者在统计数据时是否需要加上分馏值呢?还是直接用所测数据呢?

A:

捕食者和饵料的同位素值都是直接使用原始值,不需要加上分馏值,因为分馏值是直接在corrections参数中设置。

#Q[4]#

另外,想请问一下cons

A:

concs表示饵料碳氮元素含量(不是同位素值哦!) ,可以使用饵料碳氮元素含量(或百分比)得到加权校正的 SIAR 模型。

当某个饵料的某元素(如氮元素)含量(百分比)显著…..低于或高于其他饵料时,e.g., 两饵料N含量为 12%,而第3个饵料N含量为 1%,推荐使用浓度加权;

当各饵料的碳氮元素含量没有显著变化时则不需要使用浓度加权,即conc

一般情况conc默认设置为0。具体情况具体分析。

更多内容请参考 Donald L. Phillips · Paul L. Koch. Incorporating concentration dependence in stable isotope mixing models. Oecologia (2002) 130:114–125.

#Q[5]#zhao0415:

尊敬的梅老师:

您好!按照您的说明准备了数据,格式如

|sources | Mean15N | sd15N | Mean13C | sd13C

1 | src1 | 6.488984 | 1.4594632 | -11.17023 | 1.2149562

2 | src2 | 4.43216 | 2.2680709 | -30.87984 | 0.6413182

3 | src3 | 11.192613 | 1.1124385 | -11.1709 | 1.9593306

采用命令:sources

Error in make.names(col.names, unique = TRUE) : 多字节字符串1有错

此外: Warning message:

In read.table(file = file, header = header, sep = sep, quote = quote, :

incomplete final line found by readTableHeader on ‘C:/siar/re2.csv’

请问是什么原因呢?谢谢您!

A:

您好,Sd等字母大小写改成和例子完全一样之后再试试看吧

A: 数据等csv所有内容输入格式请用新罗马字体,不要用中文宋体等;另外,文件保存格式为.CSV[MS-DOS]或.CSV[Macintosh],不要保存成其他csv后缀格式。

#Q[6]# 您好。想咨询您在做稳定同位素食性分析时,在R中运行siar包时候,

如果consumer数据分组只有一组时,siarplotdata出来的图横坐标的值并不是N15的同位素值,而是group。

只要分组大于一,图像就会再次正常。想问老师有没有遇到这样的情况,这是不是安装包的bug?怎么解决这样的问题。

A:您好。当consumers数据的Group只有一组时(即只有一种消费者),在运行sair之前,直接将 consumers.csv文件里的”Group”这一列删除掉,即可在siarplotdata时导出正常的图了(e.g., x轴为d13CPI,y轴为d15NPI,见下图)!

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

① 本文仅代表作者个人观点,不代表任何其它立场,欢迎交流合作!^_^ 本文字数统计:2.9k 字

r语言 siar 代码_sair方法及R操作入门——稳定同位素之食物来源比例分析相关推荐

  1. r语言 siar 代码_平滑转换自回归(STAR)模型的R语言编程实现详解

    总体讲,ST(Smooth transition)模型,这块的code混乱的比较狗血,文献中的做法也是千奇百怪, 单变量的有 :STAR, LSTAR, ESTAR, 多变量的有:STVAR, LST ...

  2. ++代码实现 感知机的原理_决策树ID3原理及R语言python代码实现(西瓜书)

    决策树ID3原理及R语言python代码实现(西瓜书) 摘要: 决策树是机器学习中一种非常常见的分类与回归方法,可以认为是if-else结构的规则.分类决策树是由节点和有向边组成的树形结构,节点表示特 ...

  3. R语言编程的高效方法

    R语言编程的高效方法 R语言编程的高效方法 高效编程之细节知识点 并行计算 1.查看核数 2.可以使用并行计算的场景 3.Parellel 包 R语言编程的高效方法 学习资料来源: datacamp ...

  4. r语言 发送邮件html,r语言读取数据的方法

    R 对于基于 SQL 语言的关系型数据库有良好的支持,这些数据库既有商业数据库 Oracle.Microsoft SQL Server.IBM DB2 等,也包含在 GNUGeneral Public ...

  5. R语言统计代码运行耗时实战:计算代码运行时间、使用proc.time函数计算代码运行时间

    R语言统计代码运行耗时实战:计算代码运行时间.使用proc.time函数计算代码运行时间 目录

  6. R语言创建新变量方法

    R语言创建新变量方法 在数据分析中,可能需要对数据进行求和.求均值等处理,并且将处理后的数据重新保存到原来的数据框中,这里提供三种可供解决的方法: 数据框$变量名 <- 表达式.若原数据框中含有 ...

  7. R语言实现常用多重比较方法

    R语言实现常用多重比较方法 在单因素方差分析ANOVA中,如果该因素影响比较显著,那么需要进一步利用多重比较方法比较该因素不同水平的影响,确定不同水平下该因素的影响是否显著.常见的多重比较方法主要有两 ...

  8. R语言时间序列代码整理

    R语言时间序列代码整理 时间序列: data<-read.csv('/Users/ji_fanyang/Desktop/test2.csv') myseries<- ts(data[1], ...

  9. R语言编写代码的方式

    大部分的代码都是基于控制台的交互式操作,但是对于那些要重复好多次的程序片段,将其保存为一段R程序文件是一个不错的选择.通常,R程序以ASCII格式保存,扩展名为".R". 可以在类 ...

最新文章

  1. Java Spring里getBean方法的实现
  2. Java数据库连接池--DBCP浅析
  3. #leetcode刷题之路35-搜索插入位置
  4. 央视消息 | 没考驾照的人可能要恭喜了!
  5. H264帧的分析sps pps
  6. Android 系统优化(35)---Android 中如何计算 App 的启动时间?
  7. idea打开hierarchy面板
  8. linux档案与文件的的压缩与打包
  9. 谷歌:开源捐赠需分成,否则下架!
  10. presto distinct用法
  11. atitit 业务 触发器原理. 与事件原理 docx
  12. 《软件项目管理案例教程》应考概念要点及重点整理
  13. tekton pipelineresource资源
  14. c语言程序设计判断是否是素数,C语言编程输入一个数判断是否为素数(质数)...
  15. 年度盘点 | 2019品牌营销围绕的三个关键词
  16. python中group是啥_python中group和groups函数的区别
  17. 几个内核函数:filp_open、filp_read、IS_ERR、ERR_PTR、PTR_ERR
  18. 第十二周项目二----操作用邻接表存储的图
  19. 百度地图JS API 支持HTTPS
  20. 关于java中输出流flush()方法

热门文章

  1. Win10+PL2030HXA或Prolific USB-to-Serial Comm Port驱动错误或者串口无法打开问题
  2. 谓词推进 PUSHED PREDICATE
  3. 课程学习——数字频率计的设计
  4. DOTA2是暴雪的吗
  5. 青海黑枸杞形象塑造之路看地理标志产品柴达木枸杞品牌形象推广思路
  6. 信息系统项目管理师高频考点(第八章)
  7. 时域反射计(TDR)原理与应用
  8. 小学计算机社团活动简报,多彩社团活动,幸福校园生活——单集镇新河小学开展社团活动简讯...
  9. 高压放大器在工作中的应用实例有哪些
  10. 开发博客系统SpringBoot+Vue教程