# Defines the jags model to fit the two baselines trophic position full model
#
# Takes some parameters and returns a jags model object as a
# character string for passing to \code{\link[rjags]{jags.model}}. Although
# it is possible to use a number of predefined or customized
# distributions (see
# \href{https://sourceforge.net/projects/mcmc-jags/files/Manuals/}
# {JAGS documentation}),
# it is likelly that most of the time
# you will be using a normal distribution. This is the default option (i.e.
# when the function is called without arguments) and it is like this:
# "mu ~ dnorm(0, 0.0001)". In this case, a prior of normally distributed mu is
# defined, with a mean 0, and a standard deviation of 0.0001. This is a normal
# distributed prior, although uninformative. You might want to change the mean
# and/or the standard deviation according to your previously knowledge of the
# system you are working on. As well as the prior for mu, JAGS uses "tau",
# which is the precision. Precision is a deterministic function (instead of the
# distributional "~"), and it is calculated as "tau <- power(sigma, -2)", thus
# you have to define as well sigma, which stands for the standard deviation.
#
# @param muCb1 a distribution defining prior for mean (mu) for C of baseline 1.
# @param sigmaCb1 a distribution defining sigma (std dev) for C of baseline 1.
# @param muNb1 a distribution defining prior for mean (mu) for N of baseline 1.
# @param sigmaNb1 a distribution defining sigma (std dev) for N of baseline 1.
# @param muCb2 a distribution defining prior for mean (mu) for C of baseline 2.
# @param sigmaCb2 a distribution defining sigma (std dev) for C of baseline 2.
# @param muNb2 a distribution defining prior for mean (mu) for N of baseline 2.
# @param sigmaNb2 a distribution defining sigma (std dev) for N of baseline 2.
# @param alpha a distribution defining alpha (mixing model between 2 sources).
# @param sigmaCc a distribution defining sigma (std dev) for C of consumer.
# @param TP a distribution defining prior of trophic position.
# @param sigmaNc a distribution defining sigma (std dev) for N of consumer.
# @param muDeltaN a distribution defining prior for the mean (mu) of
# deltaN. deltaN stands for trophic enrichment factor of Nitrogen.
# @param sigmaDeltaN a value defining sigma (std dev) for the mean (mu) of
# deltaN.
# @param lambda an integer indicating the trophic position of the baseline.
# @param muDeltaC a distribution defining prior for the mean (mu) of
# deltaC. deltaC stands for trophic enrichment factor of Carbon
# @param sigmaDeltaC a value defining sigma (std dev) for the mean (mu) of
# deltaC.
#
# @return A jags model as a character string
#
jagsTwoBaselinesPost <- function (muCb1 = NULL,
sigmaCb1 = NULL,
muNb1 = NULL,
sigmaNb1 = NULL,
muCb2 = NULL,
sigmaCb2 = NULL,
muNb2 = NULL,
sigmaNb2 = NULL,
alpha = NULL,
sigmaCc = NULL,
TP = NULL,
sigmaNc = NULL,
muDeltaN = NULL,
sigmaDeltaN = NULL,
muDeltaC = NULL,
sigmaDeltaC = NULL,
lambda = NULL)
{
# ----------------------------------------------------------------------------
# JAGS code for fitting Inverse Wishart version of SIBER to two groups
# ----------------------------------------------------------------------------
modelString <- "
model {
# -----------------------------------------------------------------------
# First we define all the likelihood functions
# Likelihood for dC of baseline 1 (normally distributed)
for (i in 1:length(dCb1)){
dCb1[i] ~ dnorm(muCb1, tauCb1)
}
# The same is repeated for dN of baseline 1
for (i in 1:length(dNb1)){
dNb1[i] ~ dnorm(muNb1, tauNb1)
}
# And is also repeated for dN of baseline 2
for (i in 1:length(dNb2)){
dNb2[i] ~ dnorm(muNb2, tauNb2)
}
# The same is repeated for dC of baseline 2
for (i in 1:length(dCb2)){
dCb2[i] ~ dnorm(muCb2, tauCb2)
}
# Likelihood for deltaN
for (j in 1:length(deltaN)){
deltaN[j] ~ dnorm(muDeltaN, tauDeltaN)
}
# Likelihood for deltaC
for (j in 1:length(deltaC)){
deltaC[j] ~ dnorm(muDeltaC, tauDeltaC)
}
# ----------------------------------------------------------------------------
#And now we are ready to calculate the trophic position
# ----------------------------------------------------------------------------
# Likelihood for dC of consumer (dCc) is a simple mixing model of
# the dC of the two baselines
# dCc is modelled as having a normal distribution
# with mean calculated with the two baselines weighted by alpha
for (i in 1:length(dCc)) {
dCc[i] ~ dnorm(muCb2 - (muDeltaC * TP) - (alpha * (muCb2 - muCb1)), tauCc)
}
# ----------------------------------------------------------------------------
# Likelihood for the nitrogen data in the consumer uses the estimated
# proportion of baseline 1 and 2 in the consumer to inform trophic position.
for (i in 1:length(dNc)){
dNc[i] ~ dnorm(muDeltaN * (TP - lambda) + muNb1*alpha + muNb2 * (1 - alpha),
tauNc)
}"
# ----------------------------------------------------------------------------
# Priors
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# Priors for dCb1
if (is.null(muCb1)) {
newString <- "muCb1 ~ dnorm(0, 0.0001)"
} else {
newString <- paste("muCb1 ~", toString(muCb1))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaCb1)) {
newString <- "tauCb1 <- pow(sigmaCb1, -2)
sigmaCb1 ~ dunif(0, 100)"
} else {
newString <- "tauCb1 <- pow(sigmaCb1, -2)"
newString2 <- paste("sigmaCb1 ~", toString(sigmaCb1))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors for dNb1
if (is.null(muNb1)) {
newString <- "muNb1 ~ dnorm(0, 0.0001)"
} else {
newString <- paste("muNb1 ~", toString(muNb1))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaNb1)) {
newString <- "tauNb1 <- pow(sigmaNb1, -2)
sigmaNb1 ~ dunif(0, 100)"
} else {
newString <- "tauNb1 <- pow(sigmaNb1, -2)"
newString2 <- paste("sigmaNb1 ~", toString(sigmaNb1))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors for dCb2
if (is.null(muCb2)) {
newString <- "muCb2 ~ dnorm(0, 0.0001)"
} else {
newString <- paste("muCb2 ~", toString(muCb2))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaCb2)) {
newString <- "tauCb2 <- pow(sigmaCb2, -2)
sigmaCb2 ~ dunif(0, 100)"
} else {
newString <- "tauCb2 <- pow(sigmaCb2, -2)"
newString2 <- paste("sigmaCb2 ~", toString(sigmaCb2))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors for dNb2
if (is.null(muNb2)) {
newString <- "muNb2 ~ dnorm(0, 0.0001)"
} else {
newString <- paste("muNb2 ~", toString(muNb2))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaNb2)) {
newString <- "tauNb2 <- pow(sigmaNb2, -2)
sigmaNb2 ~ dunif(0, 100)"
} else {
newString <- "tauNb2 <- pow(sigmaNb2, -2)"
newString2 <- paste("sigmaNb2 ~", toString(sigmaNb2))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors on the carbon mixing model
if (is.null(alpha)) {
newString <- "alpha ~ dbeta(1,1)"
} else {
newString <- paste("alpha ~", toString(alpha))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaCc)) {
newString <- "tauCc <- pow(sigmaCc, -2)
sigmaCc ~ dunif(0, 100)"
} else {
newString <- "tauCc <- pow(sigmaCc, -2)"
newString2 <- paste("sigmaCc ~", toString(sigmaCc))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors on the dN in the consumer
if (is.null(TP)) {
newString <- "TP ~ dunif(lambda, 10)"
} else {
newString <- paste("TP ~", toString(TP))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaNc)) {
newString <- "tauNc <- pow(sigmaNc, -2)
sigmaNc ~ dunif(0, 100)"
} else {
newString <- "tauNc <- pow(sigmaNc, -2)"
newString2 <- paste("sigmaNc ~", toString(sigmaNc))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors on the deltaN (Nitrogen trophic enrichment factor)
if (is.null(muDeltaN)) {
newString <- "muDeltaN ~ dnorm(0, 0.0001)"
} else {
newString <- paste("muDeltaN ~", toString(muDeltaN))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaDeltaN)) {
newString <- "tauDeltaN <- pow(sigmaDeltaN, -2)
sigmaDeltaN ~ dunif(0, 100)"
} else {
newString <- "tauDeltaN <- pow(sigmaDeltaN, -2)"
newString2 <- paste("sigmaDeltaN ~", toString(sigmaDeltaN))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
# ----------------------------------------------------------------------------
# Priors on the deltaC (Carbon trophic enrichment factor)
if (is.null(muDeltaC)) {
newString <- "muDeltaC ~ dnorm(0, 0.0001)"
} else {
newString <- paste("muDeltaC ~", toString(muDeltaC))
}
modelString <- paste (modelString, newString, sep = "\n")
if (is.null(sigmaDeltaC)) {
newString <- "tauDeltaC <- pow(sigmaDeltaC, -2)
sigmaDeltaC ~ dunif(0, 100)"
} else {
newString <- "tauDeltaC <- pow(sigmaDeltaC, -2)"
newString2 <- paste("sigmaDeltaC ~", toString(sigmaDeltaC))
newString <- paste(newString, newString2, sep = "\n")
}
modelString <- paste (modelString, newString, sep = "\n")
#constants
if (is.null(lambda)) {
newString <- "lambda <- 2"
} else {
newString <- paste("lambda <- ", toString(lambda))
}
modelString <- paste (modelString, newString, sep = "\n")
newString <- "}" # end of jags model script
modelString <- paste (modelString, newString, sep = "\n")
return(modelString)
} # end of function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.