simmr: quick start guide

Step 1: install simmr

Use:

install.packages("simmr")

then

library(simmr)

Step 2: load in the data

Some geese isotope data is included with this package. Find where it is with:

system.file("extdata", "geese_data.xls", package = "simmr")

Load into R with:

library(readxl)
path <- system.file("extdata", "geese_data.xls", package = "simmr")
geese_data <- lapply(excel_sheets(path), read_excel, path = path)

If you want to see what the original Excel sheet looks like you can run system(paste('open',path)).

We can now separate out the data into parts

targets <- geese_data[[1]]
sources <- geese_data[[2]]
TEFs <- geese_data[[3]]
concdep <- geese_data[[4]]

Note that if you don’t have TEFs or concentration dependence you can set these all to the value 0 or just leave them blank in the step below.

Step 3: load the data into simmr

geese_simmr <- simmr_load(
  mixtures = targets[, 1:2],
  source_names = sources$Sources,
  source_means = sources[, 2:3],
  source_sds = sources[, 4:5],
  correction_means = TEFs[, 2:3],
  correction_sds = TEFs[, 4:5],
  concentration_means = concdep[, 2:3],
  group = as.factor(paste("Day", targets$Time))
)

Step 4: plot the data

plot(geese_simmr, group = 1:8)

## Warning in geom_errorbar(..., orientation = orientation): Ignoring unknown
## aesthetics: height

Step 5: run through simmr and check convergence

geese_simmr_out <- simmr_mcmc(geese_simmr)
summary(geese_simmr_out,
  type = "diagnostics",
  group = 1
)

Check that the model fitted well:

posterior_predictive(geese_simmr_out, group = 5)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 40
##    Unobserved stochastic nodes: 46
##    Total graph size: 198
## 
## Initializing model

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the bayesplot package.
##   Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Step 6: look at the output

Look at the influence of the prior:

prior_viz(geese_simmr_out)

Look at the histogram of the dietary proportions:

plot(geese_simmr_out, type = "histogram")

compare_groups(geese_simmr_out,
  groups = 1:4,
  source_name = "Enteromorpha"
)

## Most popular orderings are as follows:

##                                     Probability
## Day 428 > Day 124 > Day 398 > Day 1      0.2194
## Day 428 > Day 398 > Day 124 > Day 1      0.1758
## Day 428 > Day 124 > Day 1 > Day 398      0.1536
## Day 428 > Day 398 > Day 1 > Day 124      0.0997
## Day 428 > Day 1 > Day 124 > Day 398      0.0903
## Day 428 > Day 1 > Day 398 > Day 124      0.0622
## Day 398 > Day 428 > Day 124 > Day 1      0.0392
## Day 124 > Day 428 > Day 398 > Day 1      0.0369
## Day 124 > Day 428 > Day 1 > Day 398      0.0303
## Day 398 > Day 428 > Day 1 > Day 124      0.0228
## Day 398 > Day 124 > Day 428 > Day 1      0.0139
## Day 124 > Day 398 > Day 428 > Day 1      0.0128
## Day 1 > Day 428 > Day 398 > Day 124      0.0111
## Day 1 > Day 428 > Day 124 > Day 398      0.0108
## Day 124 > Day 1 > Day 428 > Day 398      0.0069
## Day 1 > Day 398 > Day 428 > Day 124      0.0033
## Day 1 > Day 124 > Day 428 > Day 398      0.0031
## Day 398 > Day 1 > Day 428 > Day 124      0.0022
## Day 124 > Day 1 > Day 398 > Day 428      0.0014
## Day 398 > Day 124 > Day 1 > Day 428      0.0014
## Day 1 > Day 124 > Day 398 > Day 428      0.0011
## Day 1 > Day 398 > Day 124 > Day 428      0.0008
## Day 124 > Day 398 > Day 1 > Day 428      0.0006
## Day 398 > Day 1 > Day 124 > Day 428      0.0003

For the many more options available to run and analyse output, see the main vignette via vignette('simmr')