MixSIAR贝叶斯混合模型-氮素来源-R语言实现流程

MixSIAR贝叶斯混合模型-氮素来源-R语言实现流程

建议用windows环境，Mac环境R的版本只有最新版可以下载，解决error很麻烦

环境的安装

R、R Studio
直接去官网下载（R语言建议下载R4.2.3，最新和太久的版本都容易报错）
https://mirrors.tuna.tsinghua.edu.cn/CRAN/
选择download for windows，进入下载页面，选择Previous releases，找4.2.3版本

Rstudio的版本可以下载最新

安装吉布斯采样器JAGS：https://sourceforge.net/projects/mcmc-jags/files/
直接下载安装

把上面的都下载完成之后，建议先打开R，第一次打开会让你选择源，选择国内的源对之后安装package很方便

然后关闭R，打开R studio，在Tools菜单里Install package安装mixsiar包，最简便的方法

MixSIAR里有一些实例代码，可以直接拿过来修改使用，下面是我使用的代码：

需要先准备三个csv文件，source、consumer、discrimination，放在一个文件夹下、
source： 来源的数据可以根据自己研究区去查看相关论文数据得到

consumer：是自己的数据

discirmination：和source文件的行列一样，去掉n那一列，一般把值全部设为0即可

# Geese example (concentration dependence)getwd() #得到目前的路径setwd("  ") #设置当前路径library(MixSIAR)# Load mix datamix <- load_mix_data(filename="nitrate_consumer202345.csv",                     iso_names=c("d18O","d15N"), #和自己consumer文件的顺序对应上                     factors="id",                     fac_random=FALSE,                     fac_nested=FALSE,                     cont_effects=NULL)# Load source datasource <- load_source_data(filename="nitrate_sources-1.csv",                           source_factors=NULL,                           conc_dep=FALSE,                           data_type="means",                           mix)# Load discrimination/TDF datadiscr <- load_discr_data(filename="nitrate_discrimination-1.csv", mix)# Make isospace plotplot_data(filename="isospace_plot",          plot_save_pdf=TRUE,          plot_save_png=TRUE,          mix,source,discr)# Calculate standardized convex hull arealibrary(splancs)if(mix$n.iso==2) calc_area(source=source,mix=mix,discr=discr)# Plot your prior#plot_prior(alpha.prior=1,source)# Define model structure and write JAGS model filemodel_filename <- "MixSIAR_model.txt"resid_err <- FALSEprocess_err <- TRUEwrite_JAGS_model(model_filename, resid_err, process_err, mix, source)# Run the JAGS model ("test" first, then "short")# jags.1 <- run_model(run="test", mix, source, discr, model_filename, alpha.prior=1,resid_err,process_err)jags.1 <- run_model(run="normal", mix, source, discr, model_filename, alpha.prior=1)#测试代码的时候选择test模式，要出数据选择normal或者long#设置最大输出行数，当输出过多数据的时候需要设置options(max.print=10000)# Process diagnostics, summary stats, and posterior plotsoutput_JAGS(jags.1, mix, source)

最后可以讲鼠标光标放在第一行，逐行的run，便于解决报错和观察输出结果，最后的输出结果制作氮素来源贡献图时可以选择Mean均值那一列数据。