Stable isotopes are more and more used in ecology studies. They represent a very efficient way to study food webs and analyses are now cheap enough to allow their use by anyone.

Fig: Isotopic views of food webs in the Everglades (McCutchan et al., 2003)

I am particularly interested in source partitioning, ie in assessing the relative contribution of different sources in the diet of consumers. This was until recently achieved using two sources mixing models and most of the papers were referring to Phillips and Gregg (2001). Nowadays everybody is using Bayesian mixing models and some very user-friendly R packages have been released. We could cite MixSIAR that even come with a GUI.

But all the packages and models are designed to perform the source partitioning of consumers using the raw isotopes values for consumers and the mean + sd values for sources. That makes sense because it is what almost everybody wants to do.

In my case, as I want to use data from previously published studies for a meta-analysis, this is an issue. Indeed, as studies use different models and different fractionation values for the same organisms, I need to recalculate the source partitioning for all studies with the same model and same fractionation values. But most of the studies only reports mean + sd isotope values for consumers. So I had to start from scratch and make a new Bayesian mixing model that allows me to include as a prior the sd of consumers.

I am not a statistician so I asked Andrew Parnell for help. Indeed, he developed simmr, which is already able to compute source partitioning from a unique consumer value but without taking sd into account.

With his precious help we obtained this with y_mean the mean isotope values for the consumer and sigma_obs the standard error for the consumer isotope values:

model {
# Likelihood
for (j in 1:J) {
for (i in 1:N) {
y_mean[i,j] ~ dnorm(inprod(p*q[,j], s_mean[,j]+c_mean[,j]) / inprod(p,q[,j]), 1/var_y[j])
}
var_y[j] <- inprod(pow(p*q[,j],2),pow(s_sd[,j],2)+pow(c_sd[,j],2))/pow(inprod(p,q[,j]),2)
+ pow(sigma_obs[,j],2) + pow(sigma[j],2)
}
# Prior on sigma
for(j in 1:J) { sigma[j] ~ dunif(0,0.01) }
# CLR prior on p
p[1:K] <- expf/sum(expf)
for(k in 1:K) {
expf[k] <- exp(f[k])
f[k] ~ dnorm(mu_f[k],1/pow(sigma_f[k],2))
}
# Prior on mu
for(k in 1:K) {
mu_f[k] ~ dnorm(0,1)
sigma_f[k] ~ dgamma(2,1)
}
}

You just have to run it with rjags. I will try to put an example it on my GIT with an Rnotebook soon.

If you are completely new in this area, I strongly recommend looking here for Andrew’s amazing courses….