Calculate the normalized surface area of the source convex hull

Description

calc_area() calculates the normalized surface area of the SOURCE + TDF convex hull, only if there are exactly 2 biotracers.

Usage

calc_area(source, mix, discr)

Arguments

Details

Important detail is that, unlike in Brett (2014), calc_area uses the combined SOURCE + TDF variance to normalize the surface area:

\sqrt{\sigma^2_source + \sigma^2_discr}

This is the variance used in fitting the mixing model.

calc_area() relies on the splancs::areapl() function from the splancs package. If splancs is not installed, a WARNING message will appear.

Value

If source$by_factor = FALSE, calc_area returns a scalar, the normalized surface area of the SOURCE + TDF convex hull

If source$by_factor = TRUE, calc_area returns a vector, where the entries are the normalized surface areas of the convex hull of each source factor level (e.g. source data by 3 Regions, returns a 3-vector of the areas of the Region 1 convex hull, Region 2 convex hull, etc.)

See Also

Brett (2014): https://www.researchgate.net/profile/Michael_Brett/publication/269873625_Resource_polygon_geometry_predicts_Bayesian_stable_isotope_mixing_model_bias/links/549884090cf2519f5a1de635.pdf

Combine sources from a finished MixSIAR model (a posteriori)

Description

combine_sources aggregates the proportions from multiple sources. Proportions are summed across posterior draws, since the source proportions are correlated.

Usage

combine_sources(jags.1, mix, source, alpha.prior = 1, groups)

Arguments

Details

Note: Aggregating sources after running the mixing model (a posteriori) effectively changes the prior weighting on the sources. Aggregating uneven numbers of sources will turn an 'uninformative'/generalist prior into an informative one. Because of this, combine_sources automatically generates a message describing this effect and a figure showing the original prior, the effective/aggregated prior, and what the 'uninformative'/generalist prior would be if sources were instead grouped before running the mixing model (a priori).

Value

combined, a list including:

• combined$post: matrix, posterior draws with new source groupings

• combined$source.new: list, original source list with modified entries for n.sources and source_names

• combined$groups: (input) list, shows original and combined sources

• combined$jags.1: (input) rjags model object

• combined$source.old: (input) list of original source data

• combined$mix: (input) list of original mix data

• combined$prior.old: (input) prior vector on original sources

• combined$prior.new: (output) prior vector on combined sources

See Also

summary_stat and plot_intervals

Examples

## Not run: 
# first run mantis shrimp example
# combine 6 sources into 2 groups of interest (hard-shelled vs. soft-bodied)
#   'hard' = 'clam' + 'crab' + 'snail'           # group 1 = hard-shelled prey
#   'soft' = 'alphworm' + 'brittlestar' + 'fish' # group 2 = soft-bodied prey
combined <- combine_sources(jags.1, mix, source, alpha.prior=alpha, 
                groups=list(hard=c("clam","crab","snail"), soft=c("alphworm","brittlestar","fish")))

# get posterior medians for new source groupings
apply(combined$post, 2, median)
summary_stat(combined, meanSD=FALSE, quantiles=c(.025,.5,.975), savetxt=FALSE)

## End(Not run)

Compare the predictive accuracy of 2 or more MixSIAR models

Description

compare_models uses the 'loo' package to compute LOO (leave-one-out cross-validation) or WAIC (widely applicable information criterion) for 2 of more fit MixSIAR models.

Usage

compare_models(x, loo = TRUE)

Arguments

Details

LOO and WAIC are "methods for estimating pointwise out-of-sample prediction accuracy from a fitted Bayesian model using the log-likelihood evaluated at the posterior simulations of the parameter values". See Vehtari, Gelman, & Gabry (2017). In brief:

• LOO and WAIC are preferred over AIC or DIC

• LOO is more robust than WAIC

• 'loo' estimates standard errors for the difference in LOO/WAIC between two models

• We can calculate the relative support for each model using LOO/WAIC weights

Value

Data frame with the following columns:

• Model: names of x (input list)

• LOOic / WAIC: LOO information criterion or WAIC

• se_LOOic / se_WAIC: standard error of LOOic / WAIC

• dLOOic / dWAIC: difference between each model and the model with lowest LOOic/WAIC. Best model has dLOOic = 0.

• se_dLOOic / se_dWAIC: standard error of the difference between each model and the model with lowest LOOic/WAIC

• weight: relative support for each model, calculated as Akaike weights (p.75 Burnham & Anderson 2002). Interpretation: "an estimate of the probability that the model will make the best predictions on new data, conditional on the set of models considered" (McElreath 2015).

See Also

'loo' package

Vehtari, A, A Gelman, and J Gabry. 2017. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing.

Pages 75-88 in Burnham, KP and DR Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach. Springer Science & Business Media.

Pages 188-201 in McElreath, R. 2016. Statistical rethinking: a Bayesian course with examples in R and Stan. CRC Press.

Examples

## Not run: 
# Model 1 = wolf diet by Region + Pack
mix.1 <- load_mix_data(filename=mix.filename, 
                     iso_names=c("d13C","d15N"), 
                     factors=c("Region","Pack"), 
                     fac_random=c(TRUE,TRUE), 
                     fac_nested=c(FALSE,TRUE), 
                     cont_effects=NULL)
source.1 <- load_source_data(filename=source.filename, source_factors="Region",
                             conc_dep=FALSE, data_type="means", mix.1)
discr.1 <- load_discr_data(filename=discr.filename, mix.1)

# Run Model 1
jags.1 <- run_model(run="test", mix.1, source.1, discr.1, model_filename, 
                    alpha.prior = 1, resid_err=T, process_err=T)
                    
# Model 2 = wolf diet by Region (no Pack)
mix.2 <- load_mix_data(filename=mix.filename, 
                     iso_names=c("d13C","d15N"), 
                     factors=c("Region"), 
                     fac_random=c(TRUE), 
                     fac_nested=c(FALSE), 
                     cont_effects=NULL)
source.2 <- load_source_data(filename=source.filename, source_factors="Region",
                             conc_dep=FALSE, data_type="means", mix.2)
discr.2 <- load_discr_data(filename=discr.filename, mix.2)

# Run Model 2
jags.2 <- run_model(run="test", mix.2, source.2, discr.2, model_filename, 
                    alpha.prior = 1, resid_err=T, process_err=T)
                    
# Compare models 1 and 2 using LOO
compare_models(x=list(jags.1, jags.2), loo=TRUE)

# Compare models 1 and 2 using WAIC
compare_models(x=list(jags.1, jags.2), loo=FALSE)

# Get WAIC for model 1
compare_models(x=list(jags.1), loo=FALSE)

# Create named list of rjags objects to get model names in summary
x <- list(jags.1, jags.2)
names(x) <- c("Region + Pack", "Region")
compare_models(x)

## End(Not run)

Load trophic discrimination factor (TDF) data

Description

load_discr_data loads the trophic discrimination factor (TDF) data. TDF is the amount that a consumer's tissue biotracer values are modified (enriched/depleted) after consuming a source. If tracers are conservative, then set TDF = 0 (ex. essential fatty acids, fatty acid profile data, element concentrations).

Usage

load_discr_data(filename, mix)

Arguments

Value

discr, a list including:

• discr$mu, matrix of discrimination means

• discr$sig2, matrix of discrimination variances

Load mixture data

Description

load_mix_data loads the mixture data file and names the biotracers and any Fixed, Random, or Continuous Effects.

Usage

load_mix_data(
  filename,
  iso_names,
  factors,
  fac_random,
  fac_nested,
  cont_effects
)

Arguments

Value

mix, a list including:

• mix$data: dataframe, raw mix/consumer data (all columns in 'filename'),

• mix$data_iso: matrix, mix/consumer biotracer/isotope values (those specified in 'iso_names'),

• mix$n.iso: integer, number of biotracers/isotopes,

• mix$n.re: integer, number of random effects,

• mix$n.ce: integer, number of continuous effects,

• mix$FAC: list of fixed/random effect values, each of which contains:

  • values: factor, values of the effect for each mix/consumer point

  • levels: numeric vector, total number of values

  • labels: character vector, names for each factor level

  • lookup: numeric vector, if 2 factors and Factor.2 is nested within Factor.1, stores Factor.1 values for each level of Factor.2 (e.g. Wolf Ex has 8 Packs in 3 Regions, and mix$FAC[[2]]$lookup = c(1,1,1,2,2,2,2,3), the Regions each Pack belongs to).

  • re: T/F, is the factor a Random Effect? (FALSE = Fixed Effect)

  • name: character, name of the factor (e.g. "Region")

• mix$CE: list of length n.ce, contains the cont_effects values centered (subtract the mean) and scaled (divide by SD)

• mix$CE_orig: list of length n.ce, contains the original (unscaled) cont_effects values

• mix$CE_center: vector of length n.ce, means of each cont_effects

• mix$CE_scale: vector of length n.ce, SD of each cont_effects

• mix$cont_effects: vector of length n.ce, names of each cont_effects

• mix$MU_names: vector of biotracer/iso MEAN column headings to look for in the source and discrimination files (e.g. 'd13C' in iso_names, 'Meand13C' here)

• mix$SIG_names: vector of biotracer/iso SD column headings to look for in the source and discrimination files (e.g. 'd13C' in iso_names, 'SDd13C' here)

• mix$iso_names: vector of isotope column headings in 'filename' (same as input)

• mix$N: integer, number of mix/consumer data points

• mix$n.fe: integer, number of Fixed Effects

• mix$n.effects: integer, number of Fixed Effects + Random Effects

• mix$factors: vector of length n.effects, names of the Fixed and Random Effects

• mix$fac_random: T/F vector of length n.effects indicating which effects are Random (= TRUE) and Fixed (= FALSE)

• mix$fac_nested: T/F vector of length n.effects indicating which effects are nested within the other, if any

• mix$fere: TRUE if there are 2 Fixed Effects or 1 Fixed Effect and 1 Random Effect, FALSE otherwise. Used by write_JAGS_model.

If no biotracer/isotope columns are specified, a WARNING prompts the user to select 2, 1, or 0.

If more than 2 Fixed/Random Effects are selected, a WARNING prompts the user to select 2, 1, or 0.

If more than 1 Continuous Effect is selected, a WARNING prompts the user to select 1 or 0.

Load source data

Description

load_source_data specifies the source data structure (factors, concentration dependence, data type) and loads the source data file. Sources are sorted alphabetically.

Usage

load_source_data(filename, source_factors = NULL, conc_dep, data_type, mix)

Arguments

Details

WARNING messages check for:

• More than one source factor selected

• Source factor not in mixture data

• Source sample sizes missing or entered incorrectly

• Source SD = 0

Value

source, a list including:

• source$n.sources: integer, number of sources

• source$source_names: vector, source names/labels

• source$S_MU: matrix, source means used for plotting - NOT passed to JAGS. If sources are by factor, then the third column of S_MU will be the factor values (e.g. for 4 sources and 3 Regions: 1 2 3 1 2 3 1 2 3 1 2 3)

• source$S_SIG: matrix, source SDs used for plotting - NOT passed to JAGS. Same structure as S_MU.

• source$S_factor1: factor or NULL, factor values if sources are by factor.

• source$S_factor_levels: scalar or NULL, number of S_factor1 levels if sources are by factor.

• source$conc: matrix or NULL, concentration dependence values for each isotope

• source$MU_array: array of source means, dim(src,iso,f1) or dim(src,iso) if data_type="means", NULL if data_type="raw".

• source$SIG2_array: array of source variances, dim(src,iso,f1) or dim(src,iso) if data_type="means", NULL if data_type="raw".

• source$n_array: vector/matrix of source sample sizes, dim(src,f1) or dim(src) if data_type="means", NULL if data_type="raw".

• source$SOURCE_array: array of source data, dim(src,iso,f1,replicate) or dim(src,iso,replicate) if data_type="raw", NULL if data_type="means".

• source$n.rep: vector/matrix of source sample sizes, dim(src,f1) or dim(src) if data_type="raw", NULL if data_type="means".

• source$by_factor: NA or factor number, are the source data by a Fixed or Random Effect?

• source$data_type: "raw" or "means", same as input.

• source$conc_dep: T/F, same as input.

See Also

load_mix_data and load_discr_data

mixsiar

Description

Added only to pass R CMD check with 0 NOTEs

Process mixing model output from JAGS

Description

output_JAGS processes the mixing model output, prints and saves (in the working directory):

• diagnostics

• summary statistics

• posterior density plots

• pairs plot

• trace/XY plots

Usage

output_JAGS(
  jags.1,
  mix,
  source,
  output_options = list(summary_save = TRUE, summary_name = "summary_statistics",
    sup_post = FALSE, plot_post_save_pdf = TRUE, plot_post_name = "posterior_density",
    sup_pairs = FALSE, plot_pairs_save_pdf = TRUE, plot_pairs_name = "pairs_plot", sup_xy
    = TRUE, plot_xy_save_pdf = FALSE, plot_xy_name = "xy_plot", gelman = TRUE, heidel =
    FALSE, geweke = TRUE, diag_save = TRUE, diag_name = "diagnostics", indiv_effect =
    FALSE, plot_post_save_png = FALSE, plot_pairs_save_png = FALSE, plot_xy_save_png =
    FALSE, diag_save_ggmcmc = TRUE)
)

Arguments

Value

p.both – only if 2 fixed effects OR 1 fixed + 1 random, otherwise NULL).

p.both holds the MCMC chains for the estimated proportions at the different factor levels. Dimensions = [n.draws, f1.levels, f2.levels, n.sources].

Calculated by combining the ilr offsets from global intercept: ilr.both[,f1,f2,src] = ilr.global[,src] + ilr.fac1[,f1,src] + ilr.fac2[,f2,src] And then transforming from ilr- to proportion-space.

Plot proportions by a continuous covariate

Description

plot_continuous_var creates a plot of how the mixture proportions change according to a continuous covariate, as well as plots of the mixture proportions for the individuals with minimum, median, and maximum covariate values. Called by output_JAGS if any continuous effects are in the model.

Usage

plot_continuous_var(
  jags.1,
  mix,
  source,
  output_options,
  alphaCI = 0.05,
  exclude_sources_below = 0.1
)

Arguments

Details

MixSIAR fits a continuous covariate as a linear regression in ILR/transform-space. Two terms are fit for the proportion of each source: an intercept and a slope. The plotted line uses the posterior median estimates of the intercept and slope, and the lines are curved because of the ILR-transform back into proportion-space. The 95% credible intervals are shaded.

If the model contains both a continuous AND a categorical (factor) covariate, MixSIAR fits a different intercept term for each factor level and all levels share the same slope term.

See Also

Francis et al. 2011

Plot biotracer data

Description

plot_data creates plot(s) of the biotracer data and saves the plot(s) to file(s) in the working directory. All 3 required data files must have been loaded by load_mix_data, load_source_data, and load_discr_data. Behavior depends on the number of tracers:

• 1 tracer: calls plot_data_one_iso to create a 1-D plot.

• 2 tracers: calls plot_data_two_iso to create a biplot.

• >2 tracers: calls plot_data_two_iso in a loop to create biplots for each pairwise combination of biotracers.

Usage

plot_data(
  filename,
  plot_save_pdf,
  plot_save_png,
  mix,
  source,
  discr,
  return_obj = FALSE
)

Arguments

Details

An important detail is that plot_data_two_iso and plot_data_one_iso plot the raw mix data and add the TDF to the source data, since this is the polygon that the mixing model uses to determine proportions. The plotted source means are:

\mu_source + \mu_discr

The source error bars are +/- 1 standard deviation, calculated as a combination of source and TDF variances:

\sqrt{\sigma^2_source + \sigma^2_discr}

plot_data looks for 'C', 'N', 'S', and 'O' in the biotracer column headers and assumes they are stable isotopes, labeling the axes with, e.g., expression(paste(delta^13, "C (u2030)",sep="")).

See Also

plot_data_two_iso, plot_data_one_iso

Plot biotracer data (1-D)

Description

plot_data_one_iso creates a 1-D plot of mix and source tracer data and saves the plot to a file in the working directory

Usage

plot_data_one_iso(
  mix,
  source,
  discr,
  filename,
  plot_save_pdf,
  plot_save_png,
  return_obj = FALSE
)

Arguments

Details

An important detail is that plot_data_one_iso plots the raw mix data and adds the TDF to the source data, since this is the polygon that the mixing model uses to determine proportions. The plotted source means are:

\mu_source + \mu_discr

The source error bars are +/- 1 standard deviation, calculated as a combination of source and TDF variances:

\sqrt{\sigma^2_source + \sigma^2_discr}

plot_data_one_iso looks for 'C', 'N', 'S', and 'O' in the biotracer column headers and assumes they are stable isotopes, labeling the axes with, e.g., expression(paste(delta^13, "C (u2030)",sep="")).

See Also

plot_data

Plot biotracer data (2-D)

Description

plot_data_two_iso creates a 2-D plot of mix and source tracer data and saves the plot to a file in the working directory

Usage

plot_data_two_iso(
  isotopes,
  mix,
  source,
  discr,
  filename,
  plot_save_pdf,
  plot_save_png,
  return_obj = FALSE
)

Arguments

Details

An important detail is that plot_data_two_iso plots the raw mix data and adds the TDF to the source data, since this is the polygon that the mixing model uses to determine proportions. The plotted source means are:

\mu_source + \mu_discr

The source error bars are +/- 1 standard deviation, calculated as a combination of source and TDF variances:

\sqrt{\sigma^2_source + \sigma^2_discr}

plot_data_two_iso looks for 'C', 'N', 'S', and 'O' in the biotracer column headers and assumes they are stable isotopes, labeling the axes with, e.g., expression(paste(delta^13, "C (u2030)",sep="")).

See Also

plot_data

Plot posterior uncertainty intervals from a MixSIAR model

Description

plot_intervals plots the posterior interval estimates (quantile-based) from the MCMC draws in a MixSIAR model. Calls bayesplot::mcmc_intervals.

Usage

plot_intervals(
  combined,
  toplot = "p",
  levels = NULL,
  groupby = "factor",
  savepdf = FALSE,
  filename = "post_intervals",
  ...
)

Arguments

See Also

combine_sources and summary_stat

Examples

## Not run: 
# 1. run mantis shrimp example
original <- combine_sources(jags.1, mix, source, alpha, 
                groups=list(alphworm="alphworm",brittlestar="brittlestar",clam="clam",
                           crab="crab",fish="fish",snail="snail"))
# 2. combine 6 sources into 2 groups of interest (hard-shelled vs. soft-bodied)
#   'hard' = 'clam' + 'crab' + 'snail'           # group 1 = hard-shelled prey
#   'soft' = 'alphworm' + 'brittlestar' + 'fish' # group 2 = soft-bodied prey
combined <- combine_sources(jags.1, mix, source, alpha.prior=alpha, 
                groups=list(hard=c("clam","crab","snail"), soft=c("alphworm","brittlestar","fish")))

plot_intervals(combined,toplot="fac1")
plot_intervals(original,toplot="fac1")
plot_intervals(combined,toplot="fac1",levels=1)
plot_intervals(combined,toplot="fac1",levels=2)


## End(Not run)

Plot prior

Description

plot_prior plots your prior on the global diet proportions (p.global) and the uninformative prior side-by-side. Your prior is in red, and the "uninformative"/generalist prior (alpha = 1) in dark grey.

Usage

plot_prior(
  alpha.prior = 1,
  source,
  plot_save_pdf = TRUE,
  plot_save_png = FALSE,
  filename = "prior_plot"
)

Arguments

Run the JAGS model

Description

run_model calls JAGS to run the mixing model created by write_JAGS_model. This happens when the "RUN MODEL" button is clicked in the GUI.

Usage

run_model(
  run,
  mix,
  source,
  discr,
  model_filename,
  alpha.prior = 1,
  resid_err = NULL,
  process_err = NULL
)

Arguments

Value

jags.1, a rjags model object

Note: Tracer values are normalized before running the JAGS model. This allows the same priors to be used regardless of scale of the tracer data, without using the data to select the prior (i.e. by setting the prior mean equal to the sample mean). Normalizing the tracer data does not affect the proportion estimates (p_k), but does affect users seeking to plot the posterior predictive distribution for their data. For each tracer, we calculate the pooled mean and standard deviation of the mix and source data, then subtract the pooled mean and divide by the pooled standard deviation from the mix and source data. For details, see lines 226-269.

Summary statistics from posterior of MixSIAR model

Description

summary_stat prints and saves summary statistics

Usage

summary_stat(
  combined,
  toprint = "all",
  groupby = "factor",
  meanSD = TRUE,
  quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975),
  savetxt = TRUE,
  filename = "summary_statistics"
)

Arguments

See Also

combine_sources and plot_intervals

Examples

## Not run: 
# first run mantis shrimp example
# combine 6 sources into 2 groups of interest (hard-shelled vs. soft-bodied)
#   'hard' = 'clam' + 'crab' + 'snail'           # group 1 = hard-shelled prey
#   'soft' = 'alphworm' + 'brittlestar' + 'fish' # group 2 = soft-bodied prey
combined <- combine_sources(jags.1, mix, source, alpha.prior=alpha, 
                groups=list(hard=c("clam","crab","snail"), soft=c("alphworm","brittlestar","fish")))

summary_stat(combined)
summary_stat(combined, savetxt=FALSE)
summary_stat(combined, meanSD=FALSE)
summary_stat(combined, quantiles=c(.05,.5,.95))
summary_stat(combined, toprint="fac1")
summary_stat(combined, toprint="p")
summary_stat(combined, toprint="global")

## End(Not run)

Write the JAGS model file

Description

write_JAGS_model creates "MixSIAR_model.txt", which is passed to JAGS by run_model when the "RUN MODEL" button is clicked in the GUI. Several model options will have already been specified when loading the mix and source data, but here is where the error term options are selected:

1. Residual * Process (resid_err = TRUE, process_err = TRUE)

2. Residual only (resid_err = TRUE, process_err = FALSE)

3. Process only (resid_err = FALSE, process_err = TRUE)

Usage

write_JAGS_model(
  filename = "MixSIAR_model.txt",
  resid_err = TRUE,
  process_err = TRUE,
  mix,
  source
)

Arguments

Details

WARNING messages are displayed if:

• resid_err = FALSE and process_err = FALSE are both selected.

• N=1 mix data point and did not choose "Process only" error model (MixSIR)

• Fitting each individual mix data point separately as a Fixed Effect, but did not choose "Process only" error model (MixSIR).