Several programs are provided in this page to demonstrate how to use Bayesian estimation methods via MCMC techniques to fit fish population dynamics models. Both R code and WinBUGS code are provided. See Quantitative Fish Dynamics by Quinn and Deriso (1999 ) for an exclusive exposure of fish population dynamics modeling.
von Betalanffy growth model: BUGS source code
Fish body growth equation
#*************************************************************************************## BUGS Code for fitting von Beterlanfy growth equation## Author: Dr. Zhenming Su## IFR, DNR and University of Michigan# 2005 - 2013 (C)##*************************************************************************************Ludwig von Bertalanffy Length growth modelmodel { for (i in 1:N) { Lmm[i] ~ dnorm(eta[i], prec) eta[i] <- Linf * (1-exp(-k * (Age[i] - t0))) } # Prior for the error precision (=1/Sigma^2) prec ~ dgamma(1.0E-3, 1.0E-3) sigma2 <- 1 /prec Linf ~ dunif(0, 10000) k ~ dunif(0, 10) t0 ~ dunif(-100, 100) } Data# Northern rockfish male data (from Quinn and Deriso 1999)list(N =183, Age = c(6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8,8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10,10, 10, 11, 11, 11, 11, 11, 11, 11, 12, 12, 13, 13, 13, 13, 13,13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16,16, 16, 16, 16, 17, 17, 17, 18, 18, 20, 23, 24, 24, 24, 25, 26,27, 27, 27, 29, 29, 33, 40, 6, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8,8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10,10, 10, 10, 11, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15,15, 15, 16, 16, 16, 16, 16, 17, 20, 22, 22, 22, 23, 23, 23, 23,23, 23, 23, 24, 24, 24, 24, 25, 25, 27, 28, 29, 31, 32, 32, 32,37, 40, 41, 43), Lmm = c(220, 270, 290, 290, 280, 290, 290, 310, 290, 300,310, 320, 330, 280, 290, 310, 310, 320, 290, 300, 300, 320, 290,310, 310, 310, 320, 340, 350, 300, 300, 310, 320, 320, 320, 320,330, 310, 340, 340, 310, 320, 330, 310, 320, 330, 330, 350, 330,300, 330, 350, 360, 360, 320, 330, 330, 370, 370, 330, 340, 350,350, 330, 340, 330, 330, 360, 370, 340, 340, 370, 340, 350, 360,320, 370, 340, 350, 350, 320, 360, 360, 340, 250, 260, 290, 270,290, 290, 300, 320, 270, 280, 310, 270, 280, 290, 300, 300, 310,310, 310, 310, 320, 300, 310, 330, 330, 270, 290, 290, 300, 310,340, 310, 330, 310, 340, 330, 340, 330, 340, 360, 360, 360, 360,370, 370, 370, 380, 310, 340, 350, 340, 340, 340, 360, 360, 370,370, 380, 330, 350, 360, 360, 360, 380, 340, 370, 360, 360, 400,350, 400, 350, 390, 360, 370, 360, 390, 370, 370, 370, 380, 380,400, 360, 390, 360, 370, 370, 430, 400, 410, 360, 410, 390, 410,390, 390, 410, 400))Intial values:list(Linf = 376.9, k = 0.16, t0 = -1.39, prec=1)list(Linf = 376.9, k = 0.16, t0 = -1.39, sigma=1)Results#--------------------------------------------------------------------------------------# Prior for the error precision (=1/Sigma^2)#-------------------------------------------------------------------------------------- mean sd MC_error val2.5pc median val97.5pc start sample Linf 377.4 4.751 0.06568 368.9 377.1 387.5 1001 14980 k 0.1621 0.02239 3.386E-4 0.121 0.1611 0.2084 1001 14980 sigma2 402.9 42.96 0.3365 327.2 400.2 496.0 1001 14980 t0 -1.578 1.209 0.01949 -4.241 -1.474 0.453 1001 14980DIC Dbar Dhat DIC pD Lmm 1616.0 1613.0 1619.0 2.788total 1616.0 1613.0 1619.0 2.788Minimum deviance1612.0Correlation kLinf -0.8295 kt0 0.9539 Linft0 -0.7216von Betalanffy growth model: R source code for a hybrid Metropolis-Gibbs sampler
R fish body growth MCMC code
#*************************************************************************************## R Code for a Bayesian MCMC sampler for fitting the von Beterlanfy growth equation## Author: Dr. Zhenming Su## suz@michigan.gov# IFR, DNR and University of Michigan# 2005-2013 (C)##*************************************************************************************model <- function (yobs, xobs, para, sigma2){ Lobs <- yobs nd <- length(yobs) age <- xobs Linf <- exp(para[1]) k <- exp(para[2]) t0 <- para[3] Lpred = Linf * (1-exp(-k * (age - t0))) ssq <- sum((Lobs - Lpred)^2, na.rm = T) log.post <- (-nd*log(sigma2)/2 - ssq/(sigma2*2)) list(ssq = ssq, log.post = log.post)}sigma2.update <- function(yobs, xobs, theta, sigma2){ ssq <- model(yobs, xobs, theta, sigma2)$ssq nd <- length(yobs) sigma2 <- 1 / rgamma(1, shape=(0.001 + nd / 2), scale = 1/(0.001 + ssq / 2))}update.theta.Metropolis <- function(yobs, xobs, theta, sigma2, log.post.old, var.jump.matrix){ # ------------------------------------------------------------------------- # Generate candidate random point by a mvnormal jumping distribution # ------------------------------------------------------------------------- theta.star <- mvrnorm(1, theta, var.jump.matrix); # rmvnorm log.post.new <- model(yobs, xobs, theta.star, sigma2)$log.post r.candid.prob <- exp(log.post.new - log.post.old); p.jump <- 0.0 if (r.candid.prob > runif(1)) { theta <- theta.star log.post.old <- log.post.new p.jump <- min(r.candid.prob,1) } list(theta=theta, p.jump=p.jump, log.post.old=log.post.old)}Metropolis.Algorithm <- function(yobs, xobs, theta, sigma2, var.jump.matrix, tune.const, NoIter=10000){ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Implementing a hybrid Gibbs-Metropolis sampler # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ np <- length(theta) var.jump.matrix <- (tune.const / sqrt(np))^2 * var.jump.matrix log.post.old <- model(yobs, xobs, theta, sigma2)$log.post #************************ # * Start the MC chain" * # ************************ res <- matrix(NA, NoIter, np+1, byrow=TRUE) colnames(res) <- c(names(theta),"sigma2") p.jump <- numeric(NoIter) for (i in 1:NoIter) { sigma2 <- sigma2.update(yobs, xobs, theta, sigma2) result <- update.theta.Metropolis(yobs, xobs, theta, sigma2, log.post.old, var.jump.matrix) theta <- result$theta p.jump[i] <- result$p.jump log.post.old <- result$log.post.old res[i,] <- c(theta, sigma2) } list(theta = res, p.jump = p.jump)}# Northern rockfish male data (from Quinn and Deriso 1999)lvbdata <- data.frame(age = c(6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8,8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10,10, 10, 11, 11, 11, 11, 11, 11, 11, 12, 12, 13, 13, 13, 13, 13,13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16,16, 16, 16, 16, 17, 17, 17, 18, 18, 20, 23, 24, 24, 24, 25, 26,27, 27, 27, 29, 29, 33, 40, 6, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8,8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10,10, 10, 10, 11, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,13, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15,15, 15, 16, 16, 16, 16, 16, 17, 20, 22, 22, 22, 23, 23, 23, 23,23, 23, 23, 24, 24, 24, 24, 25, 25, 27, 28, 29, 31, 32, 32, 32,37, 40, 41, 43), Lobs = c(220, 270, 290, 290, 280, 290, 290, 310, 290, 300,310, 320, 330, 280, 290, 310, 310, 320, 290, 300, 300, 320, 290,310, 310, 310, 320, 340, 350, 300, 300, 310, 320, 320, 320, 320,330, 310, 340, 340, 310, 320, 330, 310, 320, 330, 330, 350, 330,300, 330, 350, 360, 360, 320, 330, 330, 370, 370, 330, 340, 350,350, 330, 340, 330, 330, 360, 370, 340, 340, 370, 340, 350, 360,320, 370, 340, 350, 350, 320, 360, 360, 340, 250, 260, 290, 270,290, 290, 300, 320, 270, 280, 310, 270, 280, 290, 300, 300, 310,310, 310, 310, 320, 300, 310, 330, 330, 270, 290, 290, 300, 310,340, 310, 330, 310, 340, 330, 340, 330, 340, 360, 360, 360, 360,370, 370, 370, 380, 310, 340, 350, 340, 340, 340, 360, 360, 370,370, 380, 330, 350, 360, 360, 360, 380, 340, 370, 360, 360, 400,350, 400, 350, 390, 360, 370, 360, 390, 370, 370, 370, 380, 380,400, 360, 390, 360, 370, 370, 430, 400, 410, 360, 410, 390, 410,390, 390, 410, 400))attach(lvbdata)require(mvtnorm)require(MASS)#Parameter estimates by R's nonlinear least squares function: nls# Estimate Std. Error t value Pr(>|t|) #Linf 376.98625 4.35601 86.544 < 2e-16 ***#k 0.16245 0.02076 7.826 4.15e-13 ***#t0 -1.39175 1.07431 -1.295 0.197 #---#Residual standard error: 19.96 on 180 degrees of freedom#Correlation of Parameter Estimates:# Linf k#k -0.8136 #t0 -0.6769 0.968var.jump.matrix <-structure(c(0.000156188435323806, -0.00146906132939133, -0.0112326270158006, -0.00146906132939133, 0.019752582278855, 0.170063983784394, -0.0112326270158006, 0.170063983784394, 1.55449294132603), .Dim = as.integer(c(3,3)), .Dimnames = list(c("Linf", "k", "t0"), c("Linf", "k", "t0")))# some parameters are transformed to the real line# here are log(Linf), log(k), t0, sigma^2theta <- c(5.9, -1.8, -1.39)names(theta) <- c("Linf", "k", "t0")res <- Metropolis.Algorithm(lvbdata$Lobs, lvbdata$age, theta, 403.4, var.jump.matrix, tune.const = 2.4, NoIter = 10000)# transfer parameters to the original scaletheta <- data.frame(Linf = exp(res$theta[,1]), k = exp(res$theta[,2]), t0 = res$theta[,3], Sigma2 = res$theta[,4])mean(res$p.jump)var.jump.matrix <- var(res$theta)var.jump.matrix <- var.jump.matrix[1:3,1:3]# Summariesmean(theta)sd(theta)cor(theta)# scatter plotsX11()pairs(theta)# trace plotsX11()par(mfrow = c(3,1))for (j in 1 : 3) plot(theta[,j], type ="l", ylab = names(theta)[j])# autocorrelationfor (j in 1:3) print(acf(theta[,j], plot = FALSE))X11()par(mfrow = c(2,2))for (j in 1 : 3) acf(theta[,j], plot = TRUE, main = names(theta)[j])detach(lvbdata)
Stock-recruitment model
WinBUGS source code: standard Ricker model
R source code: stock-recruitment analysis by the Gibbs sampler
Surplus Production model
R source code: Hybrid Metropolis-Gibbs sampler for surplus production model
Most fish population dynamics models are nonlinear, and fisheries data are short and of time series nature. The nonlinear nature of the models can pose great difficulties to any estimation methods due to possible high correlations among model parameters. Also data at hand may not provide sufficient information for the estimation of some parameters (e.g., natural mortality). In Bayesian modeling, informative priors are needed for those not-very-well- or un-identified parameters to facilitate estimation and MCMC convergence. Even with informative priors, it may be still quite painful to use WinBUGS for Bayesian estimation of these nonlinear fish population dynamics models due to high cross-correlation among parameters. Re-parameterization and block updating the set of parameters with high correlation may ease the above problems.
For most models except Generalized Linear Model (GLM) (for the fixed-effect parameters), WinBUGS uses single-step (site) Gibbs, Metropolis Hasting and slice sampling algorithms. This can cause very slow and very inefficient sampling for models with strongly correlated parameters. In this case, we may need to write our own sampler that does block updating. The OpenBUGS program is adding more power in this aspect to block update groups of parameters that are likely to be correlated. But I think this is still a on-going process.