How to run simmr
The general structure for running simmr is as
follows:
- Call
simmr_loadon the data to get it into the right format - Plot the data in isotope space (‘iso-space’) using
plot - Run the mixing model with
simmr_mcmcorsimmr_ffvb - Check the model converged with
summary - Check the model fit is calibrated with
posterior_predictive - Explore the results with
plotandsummary, andprior_viz. If you have multiple groups and want to compare output between them, use the `compare_groups’ function
For the next part of this document, we concentrate on simple examples without grouping structure.
Step 1: Getting the data into simmr
simmr requires at minimum 3 input objects; the consumers
or mixtures, the source means, and the source
standard deviations. Optionally, you can also add correction data
(also called trophic enrichment factors, TEFs) represented again as
means and standard deviations, and concentration dependence values. The
easiest way to get data into simmr is to create an Excel file, as shown
in the included vignette('quick_start) guide. Alternatively
you can copy and past your data, comma separated, as below:
mix <- matrix(c(
-10.13, -10.72, -11.39, -11.18, -10.81, -10.7, -10.54,
-10.48, -9.93, -9.37, 11.59, 11.01, 10.59, 10.97, 11.52, 11.89,
11.73, 10.89, 11.05, 12.3
), ncol = 2, nrow = 10)
colnames(mix) <- c("d13C", "d15N")
s_names <- c("Zostera", "Grass", "U.lactuca", "Enteromorpha")
s_means <- matrix(c(-14, -15.1, -11.03, -14.44, 3.06, 7.05, 13.72, 5.96), ncol = 2, nrow = 4)
s_sds <- matrix(c(0.48, 0.38, 0.48, 0.43, 0.46, 0.39, 0.42, 0.48), ncol = 2, nrow = 4)
c_means <- matrix(c(2.63, 1.59, 3.41, 3.04, 3.28, 2.34, 2.14, 2.36), ncol = 2, nrow = 4)
c_sds <- matrix(c(0.41, 0.44, 0.34, 0.46, 0.46, 0.48, 0.46, 0.66), ncol = 2, nrow = 4)
conc <- matrix(c(0.02, 0.1, 0.12, 0.04, 0.02, 0.1, 0.09, 0.05), ncol = 2, nrow = 4)The mix object above contains the stable isotopic data
for the consumers. The data should be listed as the consumer values for
the first isotope, followed by the consumer values for the second
isotope and so on. The matrix function turns this into a
matrix (a rectangle of numbers) with 2 columns. The first column
contains the data for isotope 1, and the second the data for isotope 2.
Any number of isotopes and observations can be used. It is recommended
but not necessary to give the mixtures column names representing the
isotopes to which each column corresponds.
The source names are provided in the s_names object, and
the source means and standard deviations in s_means and
s_sds. These latter objects must also be matrices, where
the number of rows is the number of sources, and the number of columns
the number of isotopes. In each case, the data are included by listing
the values for the first isotope, then the second isotope, and so
on.
The correction data is stored in c_means and
c_sds. Again this should be a matrix of the same dimension
as s_means and s_sds. Finally the
concentration dependencies (i.e. the elemental concentration values) are
included as conc.
Some data sets are also included in simmr for quick
access to examples. See data(package = "simmr") for the
list. They can all be accessed via,
e.g. data("geese_data").
To load the data into simmr, use:
simmr_in <- simmr_load(
mixtures = mix,
source_names = s_names,
source_means = s_means,
source_sds = s_sds,
correction_means = c_means,
correction_sds = c_sds,
concentration_means = conc
)Remember that the correction_means,
correction_sds, and concentration_means are
optional.
Step 2: Plotting the data in iso-space
We can now plot the raw isotopic data with:
plot(simmr_in)
This will produce a biplot with the isotope that is in the first column on the x-axis, and the isotope in the second column on the y-axis. You can make the plot slightly nicer with some extra arguments:
plot(simmr_in,
xlab = expression(paste(delta^13, "C (\u2030)",
sep = ""
)),
ylab = expression(paste(delta^15, "N (\u2030)",
sep = ""
)),
title = "Isospace plot of example data"
)
See the help file help(plot.simmr_input) for more
options on the plotting commands, including the ability to plot
different tracers/isotopes when there are more than 2 isotopes.
If all the mixtures lie inside the mixing polygon defined by the
sources, then the data are acceptable for running simmr.
See Philips et al 2015, Canadian Journal of Zoology for more details on
when data are suitable for running through a mixing model.
Step 3: Running simmr
The next step is to actually run the model. There are two options to choose from here. The code to run the model using a Markov chain Monte Carlo (MCMC) algorithm is as follows:
simmr_out <- simmr_mcmc(simmr_in)This command takes the object simmr_in we created
earlier and uses it as input for the model. It tells simmr
to store the output from the model run in an object called
simmr_out.
Alternatively, it can be run using a Fixed Form Variational Bayes (FFVB) algorithm, as follows:
simmr_out_ffvb <- simmr_ffvb(simmr_in)The model should take less than a minute to run, though this will depend on the speed of the computer you are using. Other data sets might take slightly longer or shorter depending on the number of sources, isotopes, and observations. The progress of the model is displayed on the command line window, which shows the percentage complete.
Step 4: Checking the algorithm converged
Markov chain Monte Carlo (MCMC) works by repeatedly guessing the
values of the dietary proportions and find those values which fit the
data best. The initial guesses are usually poor and are discarded as
part of an initial phase known as the burn-in. Subsequent iterations are
then stored and used for the posterior distribution; the best
estimates of the dietary proportions given the data and the model.
Because it can take many thousands of iterations to move away from the
initial guesses, convergence diagnostics can be created to
check the model has run properly. In simmr this is done
with:
summary(simmr_out, type = "diagnostics")##
## Summary for 1
## Gelman diagnostics - these values should all be close to 1.
## If not, try a longer run of simmr_mcmc.
## deviance Zostera Grass U.lactuca Enteromorpha sd[d13C]
## 1 1 1 1 1 1
## sd[d15N]
## 1If the model run has converged properly the values should be close to
1. If they are above 1.1, we recommend a longer run. See
help(simmr_mcmc) for how to do this.
Fixed Form Variational Bayes (FFVB) is an optimisation based technique. It doesn’t require this diagnostic function.
Step 5: Checking the model fit
You can check the fit of the model with a posterior predictive check. This is similar to a fitted values plot in a linear regression. If the data points (denoted by the plot as \(y\)) broadly lie in the fitted value intervals (denoted \(y_rep\); the default is a 50% interval) then the model is fitting well:
post_pred <- posterior_predictive(simmr_out)
print(post_pred)The output includes a table showing which observations lie outside
the posterior predictive and the proportion doing so, which should
approximately match the proportion specified in the
posterior_predictive function (default 50%).
Step 6: Exploring the results
All SIMMs use informative (usually generalist) prior distributions as a default. These functions work whether the model has been run through MCMC or FFVB. You can plot the priors and the posteriors with:
prior_viz(simmr_out)
simmr produces both textual and graphical summaries of
the model run. Starting with the textual summaries, we can get tables of
the means, standard deviations and credible intervals (the Bayesian
equivalent of a confidence interval) with:
summary(simmr_out, type = "statistics")##
## Summary for 1
## mean sd
## deviance 36.200 3.067
## Zostera 0.217 0.129
## Grass 0.288 0.078
## U.lactuca 0.276 0.032
## Enteromorpha 0.218 0.120
## sd[d13C] 0.545 0.223
## sd[d15N] 0.431 0.227summary(simmr_out, type = "quantiles")##
## Summary for 1
## 2.5% 25% 50% 75% 97.5%
## deviance 32.238 33.957 35.581 37.799 43.634
## Zostera 0.029 0.111 0.200 0.304 0.492
## Grass 0.146 0.233 0.283 0.340 0.454
## U.lactuca 0.216 0.255 0.275 0.296 0.345
## Enteromorpha 0.032 0.122 0.206 0.303 0.469
## sd[d13C] 0.186 0.395 0.514 0.665 1.051
## sd[d15N] 0.051 0.283 0.408 0.558 0.956These suggest that the dietary proportions for this model are quite uncertain. However we can see that the credible interval for U.lactuca is the narrowest, running from approximately 20% to 35% of the diet. The reason this one is the narrowest can be seen from the isospace plot - this source is the most clearly separated from the others.
simmr can also produce histograms, boxplots, density
plots, and matrix plots of the output. Starting with the density
plot:
plot(simmr_out, type = "density")
We can see that Zostera and Enteromorpha are poorly constrained in comparison to Grass and especially U.lactuca. Again this is unsurprising since the isospace plot indicated that these were the two most clearly separated sources.
The most useful output plot is the matrix plot:
plot(simmr_out, type = "matrix")## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
This shows the source histograms on the diagonal, contour plots of the relationship between the sources on the upper diagonal, and the correlation between the sources on the lower diagonal. Large negative correlations indicate that the model cannot discern between the two sources; they may lie close together in iso-space. Large positive correlations are also possible when mixture data lie in a polygon consisting of multiple competing sources. Here the largest negative correlation is between Zostera and Enteromorpha. This is because they lie closest together in isospace. In general, high correlations (negative or positive) are indicative of the model being unable to determine which food sources are being consumed, and are an unavoidable part of stable isotope mixing models.
If you want to compare the dietary proportions between two different
sources, you can use the compare_sources function. This
takes two or more sources and compares the dietary proportions with an
optional plot. For example:
compare_sources(simmr_out,
source_names = c("Zostera", "U.lactuca")
)## Prob ( proportion of Zostera > proportion of U.lactuca ) = 0.328
This produces a direct probability that the dietary proportion for the first source is bigger than that of the second. If you want to compare more than two sources, specify them with:
compare_sources(simmr_out,
source_names = c(
"Zostera",
"U.lactuca",
"Enteromorpha"
)
)## Most popular orderings are as follows:
## Probability
## Zostera > U.lactuca > Enteromorpha 0.2994
## Enteromorpha > U.lactuca > Zostera 0.2767
## U.lactuca > Enteromorpha > Zostera 0.2244
## U.lactuca > Zostera > Enteromorpha 0.1706
## Zostera > Enteromorpha > U.lactuca 0.0164
## Enteromorpha > Zostera > U.lactuca 0.0125
For further information and options on comparing sources, see
help(compare_sources).
