Load the required packages:
library(SIDER) library(dplyr) library(ggplot2) # knitr::opts_chunk$set(eval = FALSE)
In this example we will illustrate how you can add additional data for known TDFs to the dataset and run SIDER. We also compare the results obtained from SIDER should the new species be added with the case when it was omitted entirely. In compiling the dataset for SIDER we applied strict rules on the data, including that the TDFs be from a strictly controlled feeding diet. In many cases this excluded data from studies where the TDFs were inferred indirectly from an experiemental manipulation such as when Derbridge et al (2015) inferred TDFs for wolves (Canis lupus) using a modified Stable Isotope Mixing Model in which the TDFs themselves were estimated. Here we will use SIDER to estimate the TDF for wolf hair when wolves are absent from the dataset, and when we include the Derbridge et al (2015) data.
It may seem like over-kill to run SIDER for a species for which you have data, and the obvious question is "why would you do this?". It comes down to how much you trust the TDFs for your species. If you are confident that they most fairly and acccurately reflect the consumers in your mixing model, then you should use them. If however they are from a different population, or at a different time, or they lack the number of replicates etc... then in all likelihood you are not as certain of them as you might be. By adding the data you do have to SIDER and then running SIDER with that same species and tissue missing, then you are getting an estimate that is more informed by the dataset total. Of course it will be weighted more heavily towards the same species data in the dataset owing to the phylogenetic correlation structure in the models, but you will also get estimates of the uncertainty that might be lacking if your data comprise only the mean estimates.
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 make take up to an hour or even more, to produce the definitive values you would used subsequently in any analysis.
##save the entire dataset as an object 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 (Meles meles) data_with_imputed_species <- recipeSider(species = "Canis_lupus", habitat = "terrestrial", taxonomic.class = "mammalia", tissue = "hair", diet.type = "carnivore", tree = combined_trees) # prepare the data for the carbon tdf model tdf_data_c <- prepareSider(data_with_imputed_species, isotope_data, combined_trees, "carbon") # prepare the data for the nitrogen tdf model tdf_data_n <- prepareSider(data_with_imputed_species, isotope_data, combined_trees, "nitrogen") # Define the formulae for both models formula_c <- delta13C ~ diet.type + habitat formula_n <- delta15N ~ diet.type + habitat # Define the random structure which is common to both models 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 and more chains. 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) # estimated TDFs for wolves that are not in (i.e. out) of # the dataset TDF_est_out_c <- imputeSider(mulTree.data = tdf_data_c, formula = formula_c, random.terms = random_terms, prior = prior, output = "wolves_out_sider_c_run", parameters = parameters, chains = n_chains, convergence = convergence, ESS = ESS) TDF_est_out_n <- imputeSider(mulTree.data = tdf_data_n, formula = formula_n, random.terms = random_terms, prior = prior, output = "wolves_out_sider_n_run", parameters = parameters, chains = n_chains, convergence = convergence, ESS = ESS)
Much of the required data has already been read in and objects created in the code above, so all that remains is that we add the wolf data from Derbridge et al (2015) and re-run the model.
# add the wolf data by binding the new data onto the included SIDER dataset derbridge_wolf_tdf_data <- rbind(SIDER_data, data.frame(species = "Canis_lupus", habitat = "terrestrial", taxonomic.class = "mammalia", tissue = "hair", diet.type = "carnivore", source.iso.13C = NA, source.iso.15N = NA, delta13C = 1.97, delta15N = 3.04, citation = c("Derbridge et al (2015)"))) # check it worked tail(derbridge_wolf_tdf_data) # Check the data for the species we want to estimate TEF for (Canis lupus) derbridge_data_with_imputed_species <- recipeSider(species = "Canis_lupus", habitat = "terrestrial", taxonomic.class = "mammalia", tissue = "hair", diet.type = "carnivore", tree = combined_trees) # prepare the carbon model derbridge_tdf_data_c <- prepareSider(derbridge_data_with_imputed_species, isotope_data, combined_trees, "carbon") # prepare the nitrogen model derbridge_tdf_data_n <- prepareSider(derbridge_data_with_imputed_species, isotope_data, combined_trees, "nitrogen") # estimated TDFs for wolves that are in the dataset TDF_est_in_c <- imputeSider(mulTree.data = derbridge_tdf_data_c, formula = formula_c, random.terms = random_terms, prior = prior, output = "wolves_in_sider_c_run", parameters = parameters, chains = n_chains, convergence = convergence, ESS = ESS) TDF_est_in_n <- imputeSider(mulTree.data = derbridge_tdf_data_n, formula = formula_n, random.terms = random_terms, prior = prior, output = "wolves_in_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")
In this example, there is very little difference in the two models. We can see this in the summary statistics and in the figure below we focus on the means and standard devations of the estimated TDFs since these are the statistics that would be passed onwards to a stable isotope mixing model (SIMM) analysis in packages such as MixSIAR, simmr or siar. See the accompanying vignette for an example pipeline on including these estimates in such an analysis.
# bundle the global estimates of TDFs from each of the 4 models into a # wide format data.frame. tdf_wide <- data.frame(DC_in = as.numeric(TDF_est_in_c$tdf_global), DC_out = as.numeric(TDF_est_out_c$tdf_global), DN_in = as.numeric(TDF_est_in_n$tdf_global), DN_out = as.numeric(TDF_est_out_n$tdf_global) ) # get some summary statistics of these estimates summary(tdf_wide) # tidy to long format for nice plotting tdf_long <- tidyr::gather(tdf_wide, model, TDF) # calculate means and standard deviations for plotting in ggplot smrys <- tdf_long %>% group_by(model) %>% summarise(mean = mean(TDF), sd = sd(TDF)) # a data.frame of the derbridge wolf TDF estimates to add to the figure # for comparison der_points = data.frame(x1 = 1, x2 = 2, x3 = 3, x4 = 4, y1 = tail(derbridge_wolf_tdf_data$delta13C,1), y2 = tail(derbridge_wolf_tdf_data$delta13C,1), y3 = tail(derbridge_wolf_tdf_data$delta15N,1), y4 = tail(derbridge_wolf_tdf_data$delta15N,1)) # make a nice ggplot of error bars showing the mean and standard deviations. ggplot(smrys, aes(model, mean)) + geom_errorbar(data = smrys, mapping = aes(ymin = mean - sd, ymax = mean + sd), size = 1, width = 0.3) + geom_point(size = 3) + scale_y_continuous(name="TDF") + geom_segment(data = der_points, mapping = aes(x = x1, y = y1, xend = x2, yend = y2), col = "red", lty = 2) + geom_segment(data = der_points, mapping = aes(x = x3, y = y3, xend = x4, yend = y4), col = "red", lty = 2) + theme_classic()
Figure 1. The means and standard deviations for carbon and nitrogen TDFs for two models: one run with the Derbridge data included in the model, and one with it left out of the model. The red horizontal lines indicate the means of the actual Derbridge et al 2015 data for carbon and nitrogen.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.