R/createJAGScode.R

In MartinLiermann/coastalCohoSS:

#' 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
}