R/simmr_mcmc_tdf.R
simmr_mcmc_tdf.Rd
This function runs a slightly different version of the main
simmr_mcmc
function with the key difference that it estimates
the correction factors (sometimes called trophic enrichment or trophic
discrimination factors; TEFs/TDFs) for a given set of dietary proportions.
simmr_mcmc_tdf(
simmr_in,
p = matrix(rep(1/simmr_in$n_sources, simmr_in$n_sources), ncol = simmr_in$n_sources,
nrow = simmr_in$n_obs, byrow = TRUE),
prior_control = list(c_mean_est = rep(2, simmr_in$n_tracers), c_sd_est = rep(2,
simmr_in$n_tracers)),
mcmc_control = list(iter = 10000, burn = 1000, thin = 10, n.chain = 4)
)
An object created via the function simmr_load
The known dietary proportions for the feeding study. Dietary proportions should be given per individual (even if they are all identical)
A list of values including arguments named means
and sd
which represent the prior means and standard deviations of the
correction factors. These can usually be left at their default values unless
you wish to include to include prior information on them.
A list of values including arguments named iter
(number of iterations), burn
(size of burn-in), thin
(amount
of thinning), and n.chain
(number of MCMC chains).
An object of class simmr_tdf
with two named top-level
components:
The simmr_input
object given to the
simmr_mcmc
function
A set of MCMC chains of class
mcmc.list
from the coda package. These can be analysed using
summary.simmr_output_tdf
.
The idea is that this code can be used for feeding studies where an organism is fed a known proportional diet with a view to estimating the correction factors to be used in a later stable isotope mixing model when the organisms are observed in the field.
The main argument of the function is an object created from
simmr_load
which contains mixture data on a number of tracers
and food source means and standard deviations. Any correction factors
included in this object will be ignored. The known dietary proportions should be provided for each individual (i.e. should be a matrix with the same number of rows as mix
). It is advisable to have multiple different dietary proportion values as part of the feeding experimental design
The output of the function is a posterior distribution on the correction
factors for each food source. Just like the output from
simmr_mcmc
, this should be checked for convergence. Examples
are included below to help assist with this check and further plots
If, after running simmr_mcmc_tdf
the convergence diagnostics in
summary.simmr_output_tdf
are not satisfactory, the values of
iter
, burn
and thin
in mcmc_control
should be
increased by e.g. a factor of 10.
Andrew C. Parnell, Donald L. Phillips, Stuart Bearhop, Brice X. Semmens, Eric J. Ward, Jonathan W. Moore, Andrew L. Jackson, Jonathan Grey, David J. Kelly, and Richard Inger. Bayesian stable isotope mixing models. Environmetrics, 24(6):387–399, 2013.
Andrew C Parnell, Richard Inger, Stuart Bearhop, and Andrew L Jackson. Source partitioning using stable isotopes: coping with too much variation. PLoS ONE, 5(3):5, 2010.
simmr_load
for creating objects suitable for this
function, simmr_mcmc
for estimating dietary proportions,
plot.simmr_input
for creating isospace plots,
summary.simmr_output_tdf
for summarising output
if (FALSE) {
## Example of estimating TDFs for a simple system with known dietary proportions
# Data set 1: 10 obs on 2 isos, 4 sources, with tefs and concdep
# Assume p = c(0.25, 0.25, 0.25, 0.25)
# The data
data(simmr_data_1)
# Load into simmr
simmr_tdf <- with(
simmr_data_1,
simmr_load(
mixtures = mixtures,
source_names = source_names,
source_means = source_means,
source_sds = source_sds,
correction_means = correction_means,
correction_sds = correction_sds,
concentration_means = concentration_means
)
)
# Plot
plot(simmr_tdf)
# MCMC run
simmr_tdf_out <- simmr_mcmc_tdf(simmr_tdf,
p = matrix(
rep(
1 / simmr_tdf$n_sources,
simmr_tdf$n_sources
),
ncol = simmr_tdf$n_sources,
nrow = simmr_tdf$n_obs,
byrow = TRUE
)
)
# Summary
summary(simmr_tdf_out, type = "diagnostics")
summary(simmr_tdf_out, type = "quantiles")
# Now put these corrections back into the model and check the
# iso-space plots and dietary output
simmr_tdf_2 <- with(
simmr_data_1,
simmr_load(
mixtures = mixtures,
source_names = source_names,
source_means = source_means,
source_sds = source_sds,
correction_means = simmr_tdf_out$c_mean_est,
correction_sds = simmr_tdf_out$c_sd_est,
concentration_means = concentration_means
)
)
# Plot with corrections now
plot(simmr_tdf_2)
simmr_tdf_2_out <- simmr_mcmc(simmr_tdf_2)
summary(simmr_tdf_2_out, type = "diagnostics")
plot(simmr_tdf_2_out, type = "boxplot")
}