#' creates the JAGS code for the model
#'
#' @examples
#' bmod <- createJAGScode()
#' @export
createJAGScode <- function(fixedObsError=NULL, smoltData=TRUE, SRtype="BH3"){
if(is.null(fixedObsError)){
obsTau <- "
smoltObsTau ~ dgamma(0.001,0.001)
escObsTau ~ dgamma(0.001,0.001)"
} else if(length(fixedObsError)==2) {
obsTau <- paste(c("\n smoltObsTau"," escObsTau"),fixedObsError^(-2),sep=" <- ", collapse="\n")
} else {
stop("Error: fixedObsError should be NULL or a vector with two values")
}
if(SRtype=="BH3"){
srFuncTxt <- "spawners[i]/(1/prod[stock[i]]^3 + spawners[i]^3/cap[stock[i]]^3)^(1/3)"
}else if(SRtype=="BH"){
srFuncTxt <- "spawners[i]/(1/prod[stock[i]] + spawners[i]/cap[stock[i]])"
}else if(SRtype=="HS"){
srFuncTxt <- "min(spawners[i]*prod[stock[i]],cap[stock[i]])"
}else{
stop("ERROR: unidentified SRtype")
}
if(smoltData){
OCtext1 <- "
logit(oceanSurv[i]) <- oceanSurvL[i] + yearEffectOS[year[i]]
oceanSurvL[i] ~ dnorm(oceanSurvPopL[stock[i]], oceanSurvPopTau[stock[i]])"
OCtext2 <- "
oceanSurvPopL[pop] ~ dnorm(oceanSurvMu, oceanSurvTau)
oceanSurvPopTau[pop] ~ dgamma(0.001,0.001)"
OCtext3 <-"
# ocean survival
oceanSurvMu ~ dnorm(oceanSurvMuPrior[1],oceanSurvMuPrior[2])
oceanSurvSD ~ dt(0,1,1)T(0,) # half cauchy with var=tau=sd=1
oceanSurvTau <- pow(oceanSurvSD,-2)
"
smoltText <- "
# smolt data
for(i in 1:Nsmolt){ # smolt trap count (spring)
smoltObs[i] ~ dlnorm(log(smolt[smoltInd[i]]),smoltObsTau)
}
"
} else {
OCtext1 <- ""
OCtext2 <- ""
OCtext3 <- ""
smoltText <- ""
}
# Create a text string with the JAGS model specification
mod1 <- paste("
model
{
###################################
########### PROCESS MODEL #########
###################################
for(i in 1:N){ # iterate over all years and populations
smolt[i] ~ dlnorm(log(muSmolt[i]), SRresidTau[stock[i]])
muSmolt[i] <- ",srFuncTxt," *
exp(yearEffect[year[i]])
escapement[i] <- smolt[i] * oceanSurv[i] * (1 - HR[i])
",
OCtext1,"
}
# offset escapement to produce spawners
# spawners are in calendar years and escapement are in brood years
for(pop in 1:Npops){
for(i in (escStartInd[pop]+3):escStopInd[pop]){
spawnersWild[i] <- escapement[i-3]
spawners[i] <- spawnersWild[i] / (1-pHOS[i])
}
# fill in the missing years with a vague prior
for(i in escStartInd[pop]:(escStartInd[pop]+2)){
spawnersWild[i] ~ dlnorm(0,0.0001)
spawners[i] <- spawnersWild[i] / (1-pHOS[i])
}
}
### population specific priors ###
for(pop in 1:Npops){
prod[pop] ~ dlnorm(logProdMu, logProdTau)
cap[pop] ~ dlnorm(log(capSlope*habVar[pop]),logCapTau)
SRresidTau[pop] ~ dgamma(0.001,0.001)
SRresidSD[pop] <- 1.0/sqrt(SRresidTau[pop])
",
OCtext2,"
}
### Hyper-priors (i.e. priors describing the distributions of parameters that vary by population)
# productivity
logProdMu ~ dnorm(prodMuPrior[1],prodMuPrior[2])
logProdSD ~ dt(0,1,1)T(0,) # half cauchy with var=tau=sd=1 # dnorm(prodSDPrior[1],prodSDPrior[2])
logProdTau <- 1.0/(logProdSD*logProdSD)
# capacity
capSlope ~ dlnorm(capSlopePrior[1],capSlopePrior[2])
logCapSD ~ dt(0,1,1)T(0,) # half cauchy with var=tau=sd=1 # dunif(0,10) #
logCapTau <- 1.0/(logCapSD*logCapSD)
",
OCtext3,"
### year effects (SR resids), smolt residuals and ocean survival common to all populations
for(k in 1:Nyears){
yearEffectTmp[k] ~ dnorm(0,yearEffectTau)
yearEffect[k] <- yearEffectTmp[k]-mean(yearEffectTmp)
yearEffectTmpOS[k] ~ dnorm(0,yearEffectTauOS)
yearEffectOS[k] <- yearEffectTmpOS[k]-mean(yearEffectTmpOS)
}
yearEffectTau ~ dgamma(0.001,0.001)
yearEffectTauOS ~ dgamma(0.001,0.001)
#######################################
########### OBSERVATION MODEL #########
#######################################
",
smoltText,"
# escapement data
for(i in 1:Nesc){
escapementObs[i] ~ dlnorm(log(spawners[escInd[i]]),escObsTau)
}
",obsTau,"
}
",sep="")
mod1
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.