brianstock / MixSIAR

A framework for Bayesian mixing models in R:
http://brianstock.github.io/MixSIAR/
94 stars 75 forks source link

How do I aggregate / combine sources a posteriori? #104

Closed brianstock closed 7 years ago

brianstock commented 7 years ago

Brian, have you written the code for posterior combination of sources after MixSIAR run? I would like to combine source proportions posteriorly. Would you mind to writing code for it if you have not done before. Thank you very much.

Kind regards, Hari

brianstock commented 7 years ago

Hi Hari,

Yes, I did write code to do that for another paper. Just in a single-case (get it to work for this dataset only) way, not in a nice, generalized, well-commented function that will work for any dataset. Would be a nice function to add to the package one of these days...

This analysis had 2 tracers and 6 sources, which we aggregated into 2 categories: "hard-shelled" vs. "soft-bodied" prey. This figure shows the aggregated prior and posterior distributions: https://github.com/brianstock/mantis_shrimp_diet/blob/master/Fig4_prior_posteriors_hard_bw.pdf. "Coral" and "Seagrass" are different habitats, can use the example for either. All you have to do is figure out which sources belong to each category and add them together for each MCMC draw.

# modelobj is your jags model output
attach.jags(modelobj)

# in this ex, cat1 is sources 3, 4, and 6. cat2 is sources 1, 2, and 5
# if you have mix by 1 factor, draws are in p.fac1[draw,fac,src]. E.g. you have mix by Season, p.fac1[,1,] has Season 1 draws)
agg.posterior <- data.frame(cat1 = p.fac1[,1,3]+p.fac1[,1,4]+p.fac1[,1,6], cat2 = p.fac1[,1,1]+p.fac1[,1,2]+p.fac1[,1,5])

# if you have no fixed/random effects, the draws are in p.global[draw,src]
agg.posterior <- data.frame(cat1 = p.global[,3]+p.global[,4]+p.global[,6], cat2 = p.global[,1]+p.global[,2]+p.global[,5])

Here is the full code: https://github.com/brianstock/mantis_shrimp_diet/blob/master/Fig4_prior_posteriors.r. See also: #53.