Load the required packages:
library(SIDER) library(MixSIAR) library(dplyr) library(ggplot2) # knitr::opts_chunk$set(eval = FALSE)
N.B. to run this example you will need to download and install the jags
software which is separate to R and separate to the R package rjags
. This is required for MixSIAR to run. You can find it at http://mcmc-jags.sourceforge.net
In this example we illustrate how to include a SIDER analysis into a workflow for a mixing model anlaysis using MixSIAR. The section Use SIDER to estimate TDFs (discrimination factors) illustrates how to run SIDER and then use the estimated TDFs in a subsequent mixing model analysis using MixSIAR. The rest of the document is not required for such an analysis, and instead it serves as a comparison of alternative methods, that we would like to see consigned ultimately to history!
Throughout this example we run relatively short runs of SIDER which results in estimates that do not pass the convergence tests. This is done to keep the build time for the package vignettes to something reasonable. As per the introduction vignette, we recommend longer runs, which may take up to an hour or even more, to produce the definitive values you would used subsequently in any analysis.
We will take a single pack from the Wolves example in MixSIAR Semmens et al 2009 which is based on a study by Darimont et al 2008. The SIDER dataset does not contain a TDF for wolf hair (the tissue type sampled in that study) nor indeed any wolf tissue sample. For the original analysis by Darimont et al 2008 and subsequently used by Semmens et al 2009 to present a new heirarchical SIMM, the authors used TDFs from fox hair a nearest species neighbour approximation using data from Roth & Hobson 2000. We compare the effects on the SIMM estimates of diet using this fox TDF with a more recently derived TDF for wolf hair that is not included in SIDER owing to it not being a directly observed TDF from a feeding experiment and the TDFs estimated using SIDER. Using a nearest species neighbour is a commonly adopted approach by researchers when their species does not have a published TDF. We will argue later that researchers should use SIDER to esimate their TDF in nearly all cases even if species-specific data are available.
The instructions to run MixSIAR are taken from the manual for that package, and further information and explanation is included in the MixSIAR package documentation.
First we load the mixture data from a *csv file. Here we are using only a single pack from the mainland population of wolves who were found to consumre more deer than marine mammals or salmon.
mix.filename <- "../inst/extdata/wolves_consumer_mainland.csv" # Load the mixture/consumer data mix <- load_mix_data(filename = mix.filename, iso_names = c("d13C","d15N"), factors = NULL, fac_random = NULL, fac_nested = NULL, cont_effects = NULL)
Now read in the source data.
# Replace the system.file call with the path to your file # source.filename <- system.file("extdata", "wolves_sources.csv", package = "SIDER") source_filename <- "../inst/extdata/wolves_sources.csv" # Load the source data source <- load_source_data(filename = source_filename, source_factors = NULL, conc_dep = FALSE, data_type = "means", mix)
Construct the jags model for this analysis using write_JAGS_model()
.
# Write the JAGS model file model_filename <- "MixSIAR_model.txt" # Name of the JAGS model file resid_err <- FALSE process_err <- TRUE # this writes a jags model as a *txt file to your working directory # which is used later to implement the model. write_JAGS_model(model_filename, resid_err, process_err, mix, source)
The nearest species to the wolf in the dataset with known TDFs is the fox (Roth & Hobson 2000). Following the procedure adopted by both Darimont et al 2008 and Semmens et al 2009 we use only the means and set their variances to 0.
# extract the fox data using dplyr::filter() fox_data <- SIDER::isotope_data %>% filter(species == "Vulpes_vulpes", tissue == "hair") # create the fox_tdfs list and populate it with matrices for the means and sd fox_tdfs <- list() # matrix of means fox_tdfs$mu <- matrix(0, ncol = 2, nrow = 3, dimnames = list(source$source_names, mix$MU_names)) # matrix of standard deviations fox_tdfs$sig2 <- matrix(0, ncol = 2, nrow = 3, dimnames = list(source$source_names, mix$SIG_names)) # add in the means for the fox tdfs that are in the SIDER dataset fox_tdfs$mu[,"Meand13C"] <- fox_data$delta13C fox_tdfs$mu[,"Meand15N"] <- fox_data$delta15N jags_fox <- run_model(run="short", mix, source, fox_tdfs, model_filename, alpha.prior = 1, resid_err, process_err)
Although not from a controlled feeding experiment, Derbridge et al (2015) inferred TDFs for wolves of $\Delta^{13}\text{C} = 1.97$ ‰ and $\Delta^{15}\text{N} = 3.04$ ‰ using a modified Stable Isotope Mixing Model in which the TDFs themselves were estimated.
derbridge_tdfs <- fox_tdfs derbridge_tdfs$mu[,"Meand13C"] <- 1.97 derbridge_tdfs$mu[,"Meand15N"] <- 3.04 derbridge_tdfs$sig2[,"SDd13C"] <- 0.70 derbridge_tdfs$sig2[,"SDd15N"] <- 0.31 jags_derbridge <- run_model(run="short", mix, source, derbridge_tdfs, model_filename, alpha.prior = 1, resid_err, process_err)
The wolf isotope samples are described fully in Darimont, C.T., Papquet, P.C. and Reimchen, T.E. (2009), Landscape heterogeneity and marine subsidy generate exensive intrapopulation niche diversity in a large terrestrial vertebrate. Journal of Animal Ecology, 78: 126-133. doi.
# Read in the data SIDER_data <- scrumpSider(iso.data = "all") ##all trees combined_trees <- scrumpSider(tree = "all") # Checking the data for the species we want to estimate TEF for (Canis lupus) new_data_test <- recipeSider(species = "Canis_lupus", habitat = "terrestrial", taxonomic.class = "mammalia", tissue = "hair", diet.type = "carnivore", tree = combined_trees) # prepare the carbon model tdf_data_c <- prepareSider(new_data_test, isotope_data, combined_trees, "carbon") # prepare the nitrogen model tdf_data_n <- prepareSider(new_data_test, isotope_data, combined_trees, "nitrogen") # formulae for both formula_c <- delta13C ~ diet.type + habitat formula_n <- delta15N ~ diet.type + habitat # common random structure for both random_terms <- ( ~ animal + species + tissue) prior <- list(R = list(V = 1, nu=0.002), G = list(G1=list(V = 1, nu=0.002), G2=list(V = 1, nu=0.002), G3=list(V = 1, nu=0.002))) # model run settings # *Might* want these to be longer runs. nitt <- c(120000) burnin <- c(20000) thin <- c(50) parameters <- c(nitt, thin, burnin) n_chains <- c(2) # convergence settings convergence = c(1.1) ESS = c(1000) TDF_est_c <- imputeSider(mulTree.data = tdf_data_c, formula = formula_c, random.terms = random_terms, prior = prior, output = "wolves_sider_c_run", parameters = parameters, chains = n_chains, convergence = convergence, ESS = ESS) TDF_est_n <- imputeSider(mulTree.data = tdf_data_n, formula = formula_n, random.terms = random_terms, prior = prior, output = "wolves_sider_n_run", parameters = parameters, chains = n_chains, convergence = convergence, ESS = ESS) ###Now lets have a look at the files imputeSider has saved to the current working directory # list.files(pattern = "test_c_run")
Take a look at the summaries of the SIDER runs. First $\Delta^{13}$C.
# Calculate summary statistics of the posterior. # Specifically, the mean and standard deviation would be # taken from here and used in a mixing model analysis using # MixSIAR, MixSIR or SIAR for example. summary(TDF_est_c$tdf_global) # Credible intervals and the mode of the posterior are obtained # using the hdrcde package hdrcde::hdr(TDF_est_c$tdf_global, prob = c(50, 95, 99)) # You can also create density plots of the posterior coda::densplot(TDF_est_c$tdf_global)
Second, $\Delta^{15}$N.
# Calculate summary statistics of the posterior. # Specifically, the mean and standard deviation would be # taken from here and used in a mixing model analysis using # MixSIAR, MixSIR or SIAR for example. summary(TDF_est_n$tdf_global) # Credible intervals and the mode of the posterior are obtained # using the hdrcde package hdrcde::hdr(TDF_est_n$tdf_global, prob = c(50, 95, 99)) # You can also create density plots of the posterior coda::densplot(TDF_est_n$tdf_global)
Now incorporate the new TDF data in a mixsiar model.
# means and variances of the SIDER estimates mu_c <- mean(TDF_est_c$tdf_global) mu_n <- mean(TDF_est_n$tdf_global) sig2_c <- var(TDF_est_c$tdf_global) sig2_n <- var(TDF_est_n$tdf_global) # construct the list object containing the mean and sd of the TDF data # to pass to MixSIAR. This is a list of length 2, each containing a matrix # with row.names of the sources, and specific column names that match the # source file. The code below extracts these from the loaded source and mixture # objects created early in the first section of this document via # load_mix_data() and load_source_data() # construct the list object sider_tdf <- list() # create the matrix of mean TDFs sider_tdf$mu <- matrix(c(rep(mu_c, 3), rep(mu_n, 3)), ncol = 2, nrow = 3, dimnames = list(source$source_names, mix$MU_names)) # create the matrix of variances for the TDFs sider_tdf$sig2 <- matrix(c(rep(sig2_c, 3), rep(sig2_n, 3)), ncol = 2, nrow = 3, dimnames = list(source$source_names, mix$SIG_names))
And finally run mixsiar
jags_sider <- run_model(run="short", mix, source, sider_tdf, model_filename, alpha.prior = 1, resid_err, process_err)
cat("Model 1: use the fox TDFs\n") knitr::kable(jags_fox$BUGSoutput$summary) cat("Model 2: Inferred wolf TDFs from Derbridge paper\n") knitr::kable(jags_derbridge$BUGSoutput$summary) cat("Model 3ehad(post: SIDER derived TDFs\n") knitr::kable(jags_sider$BUGSoutput$summary)
Bundle all the posteriors together using dplyr for plotting with ggplot. Here we focus on the proportion of deer in the diet of this wolf pack, being the largest constituent of the diet. As we can see, adding ad hoc uncertainty to the TDFs used in the original model widens the estimated proportion of deer in the diet and slightly lowers the mean estiamte (c.f. "Original" with "Added Variation" models). Using the fox as the nearest species in the SIDER dataset along with some ad hoc variation has little effect in this instance since their TDFs are similar to those used in the original model. The TDFs inferred by Derbridge et al (2015) differ from the "original" model primarily in having a higher $\Delta^{15}\text{N}$ of 3.04 ‰ compared with 2.60 ‰ and consequently the estimated proportion of deer in the diet is higher. Although the SIDER estimated $\Delta^{15}\text{N}$ is higher still at 3.32 (the $\Delta^{13}\text{C}$ is more similar) the additional uncertainty on these estimates widens and flattens the estimated proportion of deer in the diet, and also brings the esimate away from the boundary at 1, and closer to the prior which in this case is $\approx 0.33$.
We did not include the Derbridge et al (2015) TDFs in the SIDER dataset as they did not conform to our criteria as they are not derived directly from a controlled feeding trial, and instead are inferred indirectly using a modified mixing model approach. In any case, we would argue that unless you have a TDF for your specific population of consumers, and you are satisfied that they best represent their physiology and diet, then SIDER should be used to estimate the TDFs for your study. An additional option is to add the extra data to the SIDER dataset yourself and re-run the analysis to impute a wolf TDF which will now draw additional information from the wolf data in the dataset already. We provided an associated vignette illustrating how to this in SIDER.
post_df <- dplyr::bind_rows(list(as.data.frame(jags_fox$BUGSoutput$sims.matrix), as.data.frame(jags_derbridge$BUGSoutput$sims.matrix), as.data.frame(jags_sider$BUGSoutput$sims.matrix)) , .id = "Model") # change the names to match the sources names(post_df) <- c("Model", "Deviance", source$source_names) # convert to numeric and then factor post_df$Model <- as.numeric(post_df$Model) post_df$Model <- factor(post_df$Model, labels = c("Fox", "Derbridge", "SIDER")) # plot the Deer. p1 <- ggplot(post_df, mapping = aes(x = Model, y = Deer)) + geom_violin(draw_quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975)) + theme_classic() print(p1)
Figure 3. Violin plot of the estimate proportions of Deer in the diet for the 3 alternative models. Horizontal lines show the median at the middle, and then the interquartile range and 95% credible intervals.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.