Impact of climate change on seabird species off the east coast of Scotland and potential implications for environmental assessments: study

Appendix E: JAGS model code

Productivity model

model{

for(j in 1:nobs.prod){
BS.Fledged[j] ~ dbin(p.prod[j], mbs * BS.Pairs[j])
log(p.prod[j]) <- fixed.prod[j] + re.prod.year[iyear[j]] + re.prod.site[isite[j]] + re.prod.ilre[j]
fixed.prod[j] <- alpha.prod + sum(beta.prod[1:nxp] * x.prod[j,1:nxp])
re.prod.ilre[j] ~ dnorm(0, tau.prod.ilre)
}

for(s in 1:nsites){
re.prod.site[s] ~ dnorm(0, tau.prod.site)
}

for(t in 1:nyears){
re.prod.year[t] ~ dnorm(0, tau.prod.year)
}

alpha.prod ~ dnorm(0, hyper.prod.alpha)

for(k in 1:nxp){
beta.prod[k] ~ dnorm(0, hyper.prod.beta)

}

tau.prod.site ~ dgamma(hyper.prod.site[1], hyper.prod.site[2])
tau.prod.year ~ dgamma(hyper.prod.year[1], hyper.prod.year[2])
tau.prod.ilre ~ dgamma(hyper.prod.ilre[1], hyper.prod.ilre[2])

}

Survival model

model{

for(j in 1:nobs.asurv){

Pairs[j] ~ dpois(Pairs.true[j])
Pairs.true[j] <- Pairs.surv[j] + nrecuits.pairs[j]
Pairs.surv[j] ~ dpois(p.asurv[j] * Pairs.prev[j])
nrecuits.pairs[j] <- step(round(nrecruits[j] / 2))
nrecruits[j] ~ dbin(p.recruit[j], Fledged.mafb[j])
Fledged.mafb[j] <- BS.Fledged.mafb[j] + US.Fledged.mafb[j]
BS.Fledged.mafb[j] ~ dbin(mup.prod[j], mbs * BS.Pairs.mafb[j])
US.Fledged.mafb[j] ~ dbin(mup.prod[j], mbs * (Pairs.mafb[j] - BS.Pairs.mafb[j]))
p.recruit[j] <- pow(p.jsurv[j], afb)
mup.prod[j] ~ dbeta(hyper.mup.prod[1], hyper.mup.prod[2])
log(p.jsurv[j]) <- log(pe.jsurv) + link.jsurv * (log(p.asurv[j]) - log(pe.asurv))
log(p.asurv[j]) <- fixed.asurv[j] + re.asurv.year[iyear[j]] + re.asurv.site[isite[j]]
fixed.asurv[j] <- alpha.asurv + sum(beta.asurv[1:nxa] * x.asurv[j,1:nxa])
}
for(s in 1:nsites){
re.asurv.site[s] ~ dnorm(0, tau.asurv.site)
}
for(t in 1:nyears){
re.asurv.year[t] ~ dnorm(0, tau.asurv.year)
}
alpha.asurv ~ dnorm(0, hyper.asurv.alpha)
for(k in 1:nxa){
beta.asurv[k] ~ dnorm(0, hyper.asurv.beta)
}
tau.asurv.site ~ dgamma(hyper.asurv.site[1], hyper.asurv.site[2])
tau.asurv.year ~ dgamma(hyper.asurv.year[1], hyper.asurv.year[2])
}

Trend model

model{

for(j in 1:nobs.trend){
Pairs[j] ~ dpois(mu[j])
mu[j] <- log(Pairs.prev[j] + 1) + ratmu[j]
ratmu[j] <- fixed.trend[j] + re.trend.year[iyear[j]] + re.trend.site[isite[j]] + re.trend.ilre[j]
fixed.trend[j] <- alpha.trend + sum(beta.trend[1:nxt] * x.trend[j,1:nxt])
re.trend.ilre[j] ~ dnorm(0, tau.trend.ilre)
}

for(s in 1:nsites){
re.trend.site[s] ~ dnorm(0, tau.trend.site)
}

for(t in 1:nyears){
re.trend.year[t] ~ dnorm(0, tau.trend.year)
}

alpha.trend ~ dnorm(0, hyper.trend.alpha)
for(k in 1:nxt){
beta.trend[k] ~ dnorm(0, hyper.trend.beta)
}

tau.trend.site ~ dgamma(hyper.trend.site[1], hyper.trend.site[2])
tau.trend.year ~ dgamma(hyper.trend.year[1], hyper.trend.year[2])
tau.trend.ilre ~ dgamma(hyper.trend.ilre[1], hyper.trend.ilre[2])
}