revision update

Author: Deming-Yang

README.md

@@ -14,14 +14,14 @@ Model study 3 is the **Case study**, which is used to estimate possible 87Sr/86S
 The "data" folder contains both data used in this study, and published data from [Wooller et al., (2021)](https://www.science.org/doi/abs/10.1126/science.abg1134) for the case studies. 
 
 ## Software requirements
-The code is developed in R (ver. 4.0.5) calling the standalone JAGS (Just Another Gibbs Sampler) program, which is required before the code can be run. The JAGS program can be downloaded via [this link](https://sourceforge.net/projects/mcmc-jags/). Please make sure to download the version that is appropriate for your operating system.
+The code was developed in R (ver. 4.0.5) calling the standalone JAGS (Just Another Gibbs Sampler) program, which is required before the code can be run. The JAGS program can be downloaded via [this link](https://sourceforge.net/projects/mcmc-jags/). Please make sure to download the version that is appropriate for your operating system.
 
 To call JAGS from RStudio, the R packages "rjags" and "R2jags" are also required.
 
 Other R packages specified in the file headers help to visualize data.
 
 ## Related Publication
-This version of the code is for peer review only, and is subject to major updates.
+Yang et al., (in revision), BITS: a Bayesian Isotope Turnover and Sampling model for strontium isotopes in proboscideans and its potential utility in movement ecology, Methods in Ecology and Evolution.
 
 ## How to cite this software
 Yang, D. (2023), BITS: a Bayesian Isotope Turnover and Sampling model. https://github.com/SPATIAL-Lab/BITS

code/2 Driver Sr turnover.R

@@ -10,10 +10,6 @@ library(EnvStats)
 
 source("code/1 Helper functions.R")
 
-plot.col<-viridis(7)
-
-setwd("C:/Users/ydmag/Google Drive/U of U/Elephant movement/Sr-in-ivory")
-
 ####The one pool model is only used in sensitivity test in the supporting information#####
 
 misha.raw <- read.csv("data/Misha raw.csv")
@@ -45,7 +41,7 @@ plot(n.avg.misha.50.dist, n.avg.misha.50.sr, main="50 pt average",
 lines(n.avg.misha.50.dist, n.avg.misha.50.sr)
 hist(n.sd.misha.50.sr)
 
-##################2 pool turnover with 50 pt average (excursion rmv)################
+################## 2 pool turnover with 50 pt average (excursion rmv) ################
 ###prep data###
 ind.remv.50 <- which(n.avg.misha.50.dist > 17000 & n.avg.misha.50.dist < 17500 & n.avg.misha.50.sr < 0.7075)
 ind.remv.50 #36 37 38 39
@@ -57,7 +53,7 @@ n.avg.misha.50.sr.rmv <- n.avg.misha.50.sr[index.50.anom.remv]
 n.sd.misha.50.sr.rmv <- n.sd.misha.50.sr[index.50.anom.remv]
 
 
-#######This is the version used in the article, constraint on parameter b######
+####### This is the version used in the article, constraint on parameter b ######
 R.sd.mea <- n.sd.misha.50.sr.rmv
 dist.mea <- n.avg.misha.50.dist.rmv
 R.mea <- n.avg.misha.50.sr.rmv
@@ -161,5 +157,3 @@ MCMC.CI.pool.ratio$CI_high
 post.misha.pc2p3.Rin.89<- MCMC.CI.bound(post.misha.pc2p3$BUGSoutput$sims.list$Rin, 0.89)
 #calculate the mean of intake after switch
 intake.af <- mean(post.misha.pc2p3.Rin.89[[1]][90:750])
-
-

code/3 Driver Sr intake model.R

@@ -8,11 +8,6 @@ library(MASS)
 library(viridisLite)
 library(EnvStats)
 
-
-plot.col<-viridis(7)
-
-setwd("C:/Users/ydmag/Google Drive/U of U/Elephant movement/Sr-in-ivory")
-
 #adding the model component of food and water mixture as intake#
 intake <- read.csv("data/intake.csv")
 intake$Sr_conc
@@ -113,6 +108,7 @@ save(post.misha.intake, file = "out/post.misha.intake.RData")
 
 post.misha.intake$BUGSoutput$summary
 
+#save MAP estimate and HDI
 w.intake.map <- map_estimate(post.misha.intake$BUGSoutput$sims.list$w.contrib[,1])
 w.intake.hid <- hdi(post.misha.intake$BUGSoutput$sims.list$w.contrib[,1],0.89)
 w.intake.hid[2]

code/4 Driver inversion fidelity test.R

@@ -10,11 +10,7 @@ library(EnvStats)
 
 source("code/1 Helper functions.R")
 
-plot.col<-viridis(7)
-
-setwd("C:/Users/ydmag/Google Drive/U of U/Elephant movement/Sr-in-ivory")
-
-############################# inversion of Misha using params#######
+############################# inversion of Misha using parameters #######
 R.sd.mea <- n.sd.misha.50.sr.rmv
 dist.mea <- n.avg.misha.50.dist.rmv
 R.mea <- n.avg.misha.50.sr.rmv

code/5 Driver inversion case study mammoth.R

@@ -10,19 +10,15 @@ library(EnvStats)
 
 source("code/1 Helper functions.R")
 
-plot.col<-viridis(7)
-
-setwd("C:/Users/ydmag/Google Drive/U of U/Elephant movement/Sr-in-ivory")
-
-####invterting laser abalition data from Wooller et al., 2021####
+#### invterting laser abalition data from Wooller et al., 2021 ####
 Wooller <- read.csv("data/Wooller_Data_S3.csv")
 
 Wooller.sr <- Wooller$Sr_Seg01
 
 #dist is in cm, convert to micron
 wooller.micron <- Wooller$Dist_Seg01*10000
 
-#foward model simulating micromill results (500 micron band)
+#forward model simulating micromill results (500 micron band)
 mm.bwidth <- 500 #microns
 index.wooller.dist<- ceiling(wooller.micron/mm.bwidth) #this is about the same as averaging per 100 data points!
 
@@ -45,7 +41,7 @@ for(i in 1:mm.sim.n){
 
 plot(mm.sim.avg.Wooller.dist, mm.sim.avg.Wooller.sr,type="l")
 
-######ivory extension rate Wooller et al 2021####
+######tusk dentine extension rate Wooller et al 2021####
 wooller.COSr<-read.csv("data/Wooller_isotope_data.csv")
 
 wooller.rate <- rep(NA,27)
@@ -80,7 +76,7 @@ plot(sub.mm.sim.avg.dist, sub.mm.sim.avg.sr, type="l",
 #back to raw data, which is in the dist
 Wooller.sub.raw <- subset(Wooller, (wooller.micron> min(sub.mm.sim.avg.dist -250)) & (wooller.micron< 250 + max(sub.mm.sim.avg.dist)))
 
-########################inversion with precision 1e-7########
+######################## inversion with precision 1e-7 ########
 Ivo.rate.mean <- mean.wooller.rate #microns per day
 Ivo.rate.sd <- sd.wooller.rate
 

code/7 Driver sensitivity tests.R

@@ -10,10 +10,9 @@ library(EnvStats)
 
 source("code/1 Helper functions.R")
 
-plot.col<-viridis(7)
-plot.col.6<-inferno(6)
-
-setwd("C:/Users/ydmag/Google Drive/U of U/Elephant movement/Sr-in-ivory")
+#set plot colors
+plot.col <- viridis(7)
+plot.col.6 <- inferno(6)
 
 #####################################################################################
 ################################      Section 1      ################################
@@ -702,7 +701,7 @@ lines(density(post.misha.inv2p.tsrwca$BUGSoutput$sims.list$a),col="red")
 plot(density(post.misha.inv2p.tsrwca$BUGSoutput$sims.list$Rin.m.cps.ac)) #autocorrelation is weak
 
 ###################################################################################################
-##############case study: Mammoth##############
+############## case study: Mammoth ##############
 ######Mammoth inversion with normal per step error: tsrw ########
 
 Ivo.rate.mean <- mean.wooller.rate #microns per day

code/9 Plots supporting information.R

@@ -11,8 +11,6 @@ library(zoo)
 
 plot.col<-viridis(7)
 
-setwd("C:/Users/ydmag/Google Drive/U of U/Elephant movement/Sr-in-ivory")
-
 source("code/1 Helper functions.R")
 
 ####supplementary figures######