MixSIAR tutorial

MixSIAR 案例教程

MixSIAR,使用贝叶斯模型,依据指纹因子数据,估计物源对混合物的贡献比例,尤其是食物网中的食物来源问题。目前,MixSIAR在cran和github上的版本均为v3.1.12(2020-10-20发布),但是github上的,多了5个函数(函数及作用见案例5)。

本教程共包括若干水土保持学案例,以及11个官方案例(其中前6个来自github的tutorials,后5个补充自官方手册,但案例顺序略有改变)。

涉及资源包括MixSIAR包以及JAGS软件:

  1. CRAN:   https://cran.r-project.org/web/packages/MixSIAR/index.html
  2. github:   https://github.com/brianstock/MixSIAR
  3. JAGS:   https://sourceforge.net/projects/mcmc-jags/files/JAGS/
  4. RDocumentation:   https://www.rdocumentation.org/packages/MixSIAR/versions/3.1.12
  5. Official tutorials:   http://brianstock.github.io/MixSIAR/articles/index.html
  6. 官方手册:MixSIAR安装文件夹下的mixsiar_manual_small.pdf
  7. Stock B C, et al. PeerJ 2018:   https://peerj.com/articles/5096/

总述

MixSIAR经验总结

  1. 必要数据包括三个,混合物、物源、辨别因子

  2. 对于混合物:

    • 三种数据中,需第一个读取
    • 数据结构:必有指纹因子,可有协变量、连续效应变量,需在读取时指定有无、效应类型、嵌套关系等,全部为原始数据(raw形式)
    • 函数load_mix_data,共6个参数:
      • filename:char(文件路径和文件名),指定混合物文件
      • iso_names:char(列名),指定指纹因子
      • factors:NULL(无)/char(列名),指定协变量,并保留输入顺序作为后续变量名识别用
      • fac_random:NULL(无随机、无固定)/T(随机)/F(固定),认为变量是对总体的随机影响,则用随机变量;如果更关心特定部分小总体,则应该选择固定效应
      • fac_nested:NULL(无嵌套)/T(嵌套于其他因子)/F(未嵌套于其他因子)
      • cont_effects:NULL(无)/char(列名),连续效应变量
  3. 对于物源:

    • 数据结构,分为rawmeans两类:
      • raw:第一列名必须为”Sources”(注意结尾有s),变量名与mix中相同
      • means:第一列名必须为”Sources”,则全部指纹因子列名必须有前缀Mean和SD(如geese中的Meand15N、SDd15N),且必须有一列名为”n”标识出样本数量。
      • 此外,如有浓度依赖变量,则数据中需要有Conc*列,*为全部指纹因子(如geese中的浓度依赖列名为Concd15N、Concd13C)
      • 此外,模型的物源顺序按首字母排列,这可能与物源在数据中出现的顺序不同,要特别注意,最好用source$source_names查看
    • 函数load_source_data,共5个参数(不需再定义指纹因子):
      • filename:char(文件路径和文件名),指定物源文件
      • source_factors:NULL(无)/char(列名),效应因子,可以为mix和sources数据中共有的固定、随机、连续效应
      • conc_dep:T(有)/F(无),浓度依赖,数据如有,列名规则固定,只与指纹因子有关,因此这里只为T或F
      • data_typemeans(均值型)/raw(原始型)
      • mix:list,已读取的混合物数据
  4. 对于辨别因子数据

    • 数据结构:第一列为source,后几列名为所有指纹因子的Mean*和SD*;行为每个源地的对应值,每个源仅出现一次
    • 辨别因子项比较简单,仅需要依赖mix进行读取
    • 辨别因子项仅在绘制数据图、计算凸包面积和求解模型时使用
  5. 绘图:

    • 绘制先验分布时,plot_prior函数的alpha.piror,向量形式下应与源地数量相同,否则将重定向为1
    • 连续效应图,由plot_continuous_var函数绘制,包括连续变量的自身、最大值、中位数、最小值的密度函数图
    • 建议将output_JAGS中的heidel设置为默认的FALSE,若为TRUE运行时有可能报错
  6. 模型:

    • 先验信息,要求sum(alpha)==length(alpha)mean(alpha)==1,可以通过乘以长度、除以和的方式实现
    • Rstudio2022.07.2(576)版本,搭配R4.2.2,运行output_JAGS结束后,不能顺利生成png文件(文件大小为0kb,需要重启rstudio),多数情况下不能顺利生成model_result_5_diagnostics.pdf。
    • 如上问题,可以不通过鼠标关闭绘图窗口,而是使用graphics.off()关闭,则不用重启也能完成绘图
    • 运行完run_model后保存模型,再做后续分析时,通常结果会报错!
  7. 结果:

    • 最重要的文件为生成的txt:diagnostics和summary,前者记录了模型是否收敛,后者记录了最终贡献比率
    • 仅当存在随机效应时,才有posterior_density_SD图
    • github版多了5个函数,是对output_JAGS的划分,并且部分可以输出对象以便修改,这5个函数分别为:
      • output_diagnostics,对应output_JAGS的result_5
      • output_pairs, 对应output_JAGS的result_3
      • output_posteriors, 对应output_JAGS的result_2
      • output_stats, 对应output_JAGS的result_1
      • output_xy, 对应output_JAGS的result_4
  8. 其他

    • cran版本为3.1.12比github版本3.1.1还少了5个函数,这些函数是对output_JAGS的拆分
  9. 英文含义:

    • JAGS:Just Another Gibbs Sampler,一款贝叶斯网络建模工具,采用Gibbs采样算法
    • ILR:Isometric Log-Ratio transformation,等距对数比变换
    • DIC:deviance information criterion,偏差信息准则,值越小说明模型对数据更拟合
    • epsilon:multiplicative error term,乘数项误差(误差形式为残差乘以过程误差时才有这项)
    • LOO:(leave-one-out cross-validation),排除一个数据的交叉验证
    • WAIC:(widely applicable information criterion),广泛适当信息准则
    • sigma:标准差,当模型中有随机效应时,output_JAGS或output_posteriors将给出每个效应的标准差密度图,用以比较哪些因素影响捕食贡献
    • TDF:trophic discrimination factor (TDF). TDF is the amount that a consumer’s tissue biotracer values are modified (enriched/depleted) after consuming a source. If tracers are conservative, then set TDF = 0 (ex. essential fatty acids, fatty acid profile data, element concentrations)
  10. MCMC设置: