Bayesian methods using JAGS

Currently, we use Stan and Julia's Turing and also exploring the utility of Google's TensorFlow (via Greta) and R.

NOTE: The following is legacy code. Bugs is still very useful and interesting. But do consider Greta and Turing.

The use of WinBUGS is a powerful means of leveraging Bayesian methods. Its use is however a bit complicated in a non-MSWindows environment. I use the excellent and opensource: JAGS (Just Another Gibbs Sampler) created by Martyn Plummer, which is almost completely compatible with the BUGS syntax used in WinBUGS and ported to Linux, MAC and MSWindows. A direct connection to R is provided by the R-library, rjags. This means that a lot of things are more simply and flexibly done with JAGS. Install JAGS and the R-library rjags ( install.packages( "rjags" ) ).

The following is an example of how to use it to do some parameter fitting to a simple population dynamic model. For more details about the model, see the current version of the snow crab assessment. You will need to have some familiarity with Bayesian methods to understand what is being done below. Sorry, I will work on an introductory level tutorial when I get the chance.

source( "http://bio-r.googlecode.com/files/functions.surplus.production.r" )

source( "http://bio-r.googlecode.com/files/functions.jags.r" )

require(rjags)

# load a few useful modules (JAGS is modular in design)

rjags::load.module("dic")

# rjags::load.module("glm")

# load some example data

load( url( "http://bio-r.googlecode.com/files/L.bugs.test.rdata" ) ) # load L

load( url( "http://bio-r.googlecode.com/files/B.bugs.test.rdata" ) ) # load B

BI = B # rename to avoid confusion as "B" is used below

abundance.index = B

# makes sure data are from overlapping years

oo = which(rownames(L) %in% rownames(BI) )

landings = L[ oo, ]

# parameter list

sb = list(

b0x = c(0.8, 0.6, 0.1), # biomass scaled by carrying capacity

q0x = c(1, 1, 1), # catchability

r0x = c(1, 1, 1), # intrinsic rate of increase

K0x = c(4, 50, 1 ), # carrying capacity

O = abundance.index , # observed index of abundance

C = landings , # catches

sigma = 1/4,

er = 0.2, # target exploitation rate .. for projections

nProj = 5 # TACS for prediction scenarios

)

# here is the core of the model (a observation and process error model)

# It is stored for convenience in "functions.surplus.production.r"

cat ( biomass.logistic.hierarchical.jags )

# set up the model for a run and start a burn-in

m = jags.model( file=jags.setup( biomass.logistic.hierarchical.jags ), data=sb, n.chains=3, n.adapt=20000 )

# choose objects to monitor

tomonitor = c("B", "r", "K", "q", "MSY", "BMSY", "FMSY", "XMSY",

"sd.o", "sd.p", "sd.r", "sd.K", "sd.q", "TAC" )

# check coefficients of the model to make sure starting parameters are reasonable

coef(m)

# pDIC :

m.dic = dic.samples(m, n.iter=10000)

# Mean deviance: 9.24 penalty 189 Penalized deviance: 198

# check for convergence of MCMC

y = jags.samples(m, variable.names=tomonitor, n.iter=20000, thin=10)

gelman.plot(y[["r"]])

gelman.plot(y[["K"]])

# 10000 samples seems enough but a few more to be sure

update(m, n.iter=50000 )

# autocorrelation thinning

y = coda.samples(m, variable.names=c("K", "r", "q"), n.iter=20000, thin=10) # sample from posterior

autocorr.plot(y)

# thin the data

thin=400

niter=1000*thin

# sample from posterior

y = jags.samples(m, variable.names=tomonitor, n.iter=niter, thin=thin)

# extract data and do some summary statistics

B = jags.extract( y, "B", mean )

K = jags.extract( y, "K", mean )

# instantaneous fishing mortality estimates

summary(y$F, median)

# a few plots

ntacs = sb$nProj

yrs0 = as.numeric( as.character( rownames(BI) ) )

yrs = c( yrs0, (max(yrs0)+c(1:ntacs) ) )

yrs.last = max(yrs0) + 0.5

ndata = length(yrs0)

hdat = 1:ndata

graphics.off()

layout( matrix(c(1,2,3), 3, 1 ))

par(mar = c(5, 4, 0, 2))

# timeseries of biomass

for (i in 1:3) {

yran = range(c(0, B[,i], BI[,i]/ q[i], BI[,i] ))*1.05

plot( yrs, B[,i] , type="l", lwd=5, col="darkgray", ylim=yran, xlab="", ylab="" )

points( yrs0, BI[,i]/ q[i], pch=20, col="black", cex=2 )

abline (v=yrs.last , lwd=2, lty="dashed" )

if (i==2) title( ylab="Fishable biomass (kt)" )

if (i==3) title( xlab="Year" )

points( yrs0, BI[,i], pch=20, col="darkgreen" )

lines ( yrs0, BI[,i], lwd=2, col="darkgreen", lty="dashed" )

}

# Harvest Control plots

graphics.off()

layout( matrix(c(1,2,3), 3, 1 ))

par(mar = c(5, 4, 0, 2))

require(car)

for (i in 1:3 ) {

# x11()

plot( B[hdat,i], F[hdat,i], type="b", xlim=c(0, K[i] * 1.25 ),

ylim=c(0, 0.7), col="darkorange", cex=0.8, lwd=2, xlab="", ylab="", pch=20 )

text( B[hdat,i], F[hdat,i], labels=yrs0, pos=3, cex= 0.8 )

nn = as.matrix( cbind( Bx=as.vector( y$B[ndata,i,,] ), Fx = as.vector( y$F[ndata,i,,] ) ))

ellipse.2d(nn[,1], nn[,2], pv=0.1, sc=150)

if (i==3) title( xlab="Fishable biomass (kt)" )

if (i==2) title( ylab="Fishing mortality" )

Fmsy = FMSY[,i]

Fref = FMSY[,i] * 0.2

Bmsy = BMSY[,i]

Bref = K[,i] * 0.2

BK = K[,i]

Fhistorical = mean( F[hdat,i], na.rm=T )

Bhistorical = mean( B[hdat,i], na.rm=T )

yl = 0.6

# abline (h=Fref, lty="dotted")

# text( 0.25, Fref, "Target\n (0.2 * FMSY) ", pos=1 )

abline (h=Fmsy, lty="dotted")

text( 0.05*K[i], Fmsy, "FMSY", pos=3 )

abline (h=Fhistorical, lty="dashed")

text( 0.05*K[i], Fhistorical, "Historical\n mean", pos=1, lwd=2 )

# abline (v=Bref, lty="dotted")

# text( Bref-0.2, 0.25, "Lower biomass reference point\n (LBRP = 0.2 * BMSY)" , srt=90, pos=3)

abline (v=Bmsy, lty="dotted")

text( Bmsy-0.01*K[i], yl, "BMSY" , srt=90, pos=3)

abline (v=BK, lty="dotted")

text( BK-0.01*K[i], yl, "K" , srt=90, pos=3)

abline (v=Bhistorical, lty="dashed")

text( Bhistorical-0.01*K[i], yl, "Historical\n mean" , srt=90, pos=3, lwd=2)

}

# densities of K

graphics.off()

layout( matrix(c(1,2,3), 3, 1 ))

par(mar = c(5, 4, 0, 2))

for (i in 1:3) plot(density( y$K[i,,] ), main="")

qs = apply( y$K[,,], 1, quantile, probs=c(0.025, 0.5, 0.975) )