Gianpaolo Coro 2015-03-30 15:31:04 +00:00
parent 2a58553145
commit b7b5df763b
33 changed files with 1362 additions and 0 deletions

View File

@ -0,0 +1,696 @@
##--------------------------------------------------------
## CMSY analysis with estimation of total biomass, including Bayesian Schaefer
## written by Rainer Froese with support from Gianpaolo Coro in 2013-2014
## This version adjusts biomass to average biomass over the year
## It also contains the FutureCrash option to improve prediction of final biomass
## Version 21 adds the purple point to indicate the 25th percentile of final biomass
## Version 22 accepts that no biomass or CPUE area available
##--------------------------------------------------------
library(R2jags) # Interface with JAGS
library(coda)
#-----------------------------------------
# Some general settings
#-----------------------------------------
# set.seed(999) # use for comparing results between runs
rm(list=ls(all=TRUE)) # clear previous variables etc
options(digits=3) # displays all numbers with three significant digits as default
graphics.off() # close graphics windows from previous sessions
#-----------------------------------------
# General settings for the analysis
#-----------------------------------------
sigR <- 0.02 # overall process error; 0.05 works reasonable for simulations, 0.02 for real data; 0 if deterministic model
n <- 10000 # initial number of r-k pairs
batch.mode <- T # set to TRUE to suppress graphs
write.output <- T # set to true if table of output is wanted
FutureCrash <- "No"
#-----------------------------------------
# Start output to screen
#-----------------------------------------
cat("-------------------------------------------\n")
cat("Catch-MSY Analysis,", date(),"\n")
cat("-------------------------------------------\n")
#------------------------------------------
# Read data and assign to vectors
#------------------------------------------
# filename_1 <- "AllStocks_Catch4.csv"
# filename_2 <- "AllStocks_ID4.csv"
# filename_1 <- "SimCatch.csv"
# filename_2 <- "SimSpec.csv"
# filename_2 <- "SimSpecWrongS.csv"
# filename_2 <- "SimSpecWrongI.csv"
# filename_2 <- "SimSpecWrongF.csv"
# filename_2 <- "SimSpecWrongH.csv"
# filename_2 <- "SimSpecWrongL.csv"
# filename_1 <- "FishDataLim.csv"
# filename_2 <- "FishDataLimSpec.csv"
filename_1 <- "WKLIFE4Stocks.csv"
filename_2 <- "WKLIFE4ID.csv"
outfile<-"outfile"
outfile.txt <- "outputfile.txt"
cdat <- read.csv(filename_1, header=T, dec=".", stringsAsFactors = FALSE)
cinfo <- read.csv(filename_2, header=T, dec=".", stringsAsFactors = FALSE)
cat("Files", filename_1, ",", filename_2, "read successfully","\n")
# Stocks with total biomass data and catch data from StartYear to EndYear
# stocks <- sort(as.character(cinfo$stock)) # All stocks
stocks<-"HLH_M07"
# select one stock after the other
for(stock in stocks) {
# assign data from cinfo to vectors
res <- as.character(cinfo$Resilience[cinfo$stock==stock])
StartYear <- as.numeric(cinfo$StartYear[cinfo$stock==stock])
EndYear <- as.numeric(cinfo$EndYear[cinfo$stock==stock])
r_low <- as.numeric(cinfo$r_low[cinfo$stock==stock])
r_hi <- as.numeric(cinfo$r_hi[cinfo$stock==stock])
stb_low <- as.numeric(cinfo$stb_low[cinfo$stock==stock])
stb_hi <- as.numeric(cinfo$stb_hi[cinfo$stock==stock])
intyr <- as.numeric(cinfo$intyr[cinfo$stock==stock])
intbio_low <- as.numeric(cinfo$intbio_low[cinfo$stock==stock])
intbio_hi <- as.numeric(cinfo$intbio_hi[cinfo$stock==stock])
endbio_low <- as.numeric(cinfo$endbio_low[cinfo$stock==stock])
endbio_hi <- as.numeric(cinfo$endbio_hi[cinfo$stock==stock])
Btype <- as.character(cinfo$Btype[cinfo$stock==stock])
FutureCrash <- as.character(cinfo$FutureCrash[cinfo$stock==stock])
comment <- as.character(cinfo$comment[cinfo$stock==stock])
# extract data on stock
yr <- as.numeric(cdat$yr[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])
ct <- as.numeric(cdat$ct[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000 ## assumes that catch is given in tonnes, transforms to '000 tonnes
if(Btype=="observed" | Btype=="CPUE" | Btype=="simulated") {
bt <- as.numeric(cdat$TB[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000 ## assumes that biomass is in tonnes, transforms to '000 tonnes
} else {bt <- NA}
nyr <- length(yr) # number of years in the time series
if(Btype!="observed") {bio <- bt}
# change biomass to moving average as assumed by Schaefer (but not for simulations or CPUE)
# for last year use reported bio
if(Btype=="observed") {
ma <- function(x){filter(x,rep(1/2,2),sides=2)}
bio <- ma(bt)
bio[length(bio)] <- bt[length(bt)] }
# initialize vectors for viable r, k, bt
rv.all <- vector()
kv.all <- vector()
btv.all <- matrix(data=vector(),ncol=nyr+1)
#----------------------------------------------------
# Determine initial ranges for parameters and biomass
#----------------------------------------------------
# initial range of r from input file
if(is.na(r_low)==F & is.na(r_hi)==F) {
start_r <- c(r_low,r_hi)
} else {
# initial range of r and CatchMult values based on resilience
if(res == "High") {
start_r <- c(0.6,1.5)} else if(res == "Medium") {
start_r <- c(0.2,0.8)} else if(res == "Low") {
start_r <- c(0.05,0.5)} else { # i.e. res== "Very low"
start_r <- c(0.015,0.1)}
}
# initial range of k values, assuming k will always be larger than max catch
# and max catch will never be smaller than a quarter of MSY
start_k <- c(max(ct),16*max(ct)/start_r[1])
# initial biomass range from input file
if(is.na(stb_low)==F & is.na(stb_hi)==F) {
startbio <- c(stb_low,stb_hi)
} else {
# us low biomass at start as default
startbio <- c(0.1,0.5)
}
MinYear <- yr[which.min(ct)]
MaxYear <- yr[which.max(ct)]
# use year and biomass range for intermediate biomass from input file
if(is.na(intbio_low)==F & is.na(intbio_hi)==F) {
intyr <- intyr
intbio <- c(intbio_low,intbio_hi)
# else if year of minimum catch is at least 3 years away from StartYear and EndYear of series, use min catch
} else if((MinYear - StartYear) > 3 & (EndYear - MinYear) > 3 ) {
# assume that biomass range in year before minimum catch was 0.01 - 0.4
intyr <- MinYear-1
intbio <- c(0.01,0.4)
# else if year of max catch is at least 3 years away from StartYear and EndYear of series, use max catch
} else if((MaxYear - StartYear) > 3 & (EndYear - MaxYear) > 3 ) {
# assume that biomass range in year before maximum catch was 0.3 - 0.9
intyr <- MaxYear-1
intbio <- c(0.3,0.9)
} else {
# assume uninformative range 0-1 in mid-year
intyr <- as.integer(mean(c(StartYear, EndYear)))
intbio <- c(0,1) }
# end of intbio setting
# final biomass range from input file
if(is.na(endbio_low)==F & is.na(endbio_hi)==F) {
endbio <- c(endbio_low,endbio_hi)
} else {
# else use Catch/maxCatch to estimate final biomass
endbio <- if(ct[nyr]/max(ct) > 0.5) {c(0.4,0.8)} else {c(0.01,0.4)}
} # end of final biomass setting
#----------------------------------------------
# MC with Schaefer Function filtering
#----------------------------------------------
Schaefer <- function(ri, ki, startbio, intyr, intbio, endbio, sigR, pt) {
# if stock is not expected to crash within 3 years if last catch continues
if(FutureCrash == "No") {
yr.s <- c(yr,EndYear+1,EndYear+2,EndYear+3)
ct.s <- c(ct,ct[yr==EndYear],ct[yr==EndYear],ct[yr==EndYear])
nyr.s <- length(yr.s)
} else{
yr.s <- yr
ct.s <- ct
nyr.s <- nyr
}
# create vector for initial biomasses
startbt <-seq(from =startbio[1], to=startbio[2], by = (startbio[2]-startbio[1])/10)
# create vectors for viable r, k and bt
rv <- array(-1:-1,dim=c(length(ri)*length(startbt))) #initialize array with -1. The -1 remaining after the process will be removed
kv <- array(-1:-1,dim=c(length(ri)*length(startbt)))
btv <- matrix(data=NA, nrow = (length(ri)*length(startbt)), ncol = nyr+1)
intyr.i <- which(yr.s==intyr) # get index of intermediate year
#loop through r-k pairs
npoints = length(ri)
nstartb = length(startbt)
for(i in 1 : npoints) {
if (i%%1000==0)
cat(".")
# create empty vector for annual biomasses
bt <- vector()
# loop through range of relative start biomasses
for(j in startbt) {
# set initial biomass, including process error
bt[1]=j*ki[i]*exp(rnorm(1,0, sigR)) ## set biomass in first year
#loop through years in catch time series
for(t in 1:nyr.s) { # for all years in the time series
xt=rnorm(1,0, sigR) # set new random process error for every year
# calculate biomass as function of previous year's biomass plus surplus production minus catch
bt[t+1]=(bt[t]+ri[i]*bt[t]*(1-bt[t]/ki[i])-ct.s[t])*exp(xt)
# if biomass < 0.01 k or > 1.1 k, discard r-k pair
if(bt[t+1] < 0.01*ki[i] || bt[t+1] > 1.1*ki[i]) { break } # stop looping through years, go to next upper level
if ((t+1)==intyr.i && (bt[t+1]>(intbio[2]*ki[i]) || bt[t+1]<(intbio[1]*ki[i]))) { break } #intermediate year check
} # end of loop of years
# if last biomass falls without expected ranges goto next r-k pair
if(t < nyr.s || bt[yr.s==EndYear] > (endbio[2]*ki[i]) || bt[yr.s==EndYear] < (endbio[1]*ki[i])) {
next } else {
# store r, k, and bt, plot point, then go to next startbt
rv[((i-1)*nstartb)+j] <- ri[i]
kv[((i-1)*nstartb)+j] <- ki[i]
btv[((i-1)*nstartb)+j,] <- bt[1:(nyr+1)]/ki[i] #substitute a row into the matrix, exclude FutureCrash years
if(pt==T) {points(x=ri[i], y=ki[i], pch=".", cex=2, col="black")
next }
}
} # end of loop of initial biomasses
} # end of loop of r-k pairs
rv=rv[rv!=-1]
kv=kv[kv!=-1]
btv=na.omit(btv) #delete first line
cat("\n")
return(list(rv, kv,btv))
} # end of Schaefer function
#------------------------------------------------------------------
# Uniform sampling of the r-k space
#------------------------------------------------------------------
# get random set of r and k from log space distribution
ri1 = exp(runif(n, log(start_r[1]), log(start_r[2])))
ki1 = exp(runif(n, log(start_k[1]), log(start_k[2])))
#-----------------------------------------------------------------
# Plot data and progress
#-----------------------------------------------------------------
#windows(14,9)
par(mfcol=c(2,3))
# plot catch
plot(x=yr, y=ct, ylim=c(0,1.2*max(ct)), type ="l", bty="l", main=paste(stock,"catch"), xlab="Year",
ylab="Catch", lwd=2)
points(x=yr[which.max(ct)], y=max(ct), col="red", lwd=2)
points(x=yr[which.min(ct)], y=min(ct), col="red", lwd=2)
# plot r-k graph
plot(ri1, ki1, xlim = start_r, ylim = start_k, log="xy", xlab="r", ylab="k", main="Finding viable r-k", pch=".", cex=2, bty="l", col="lightgrey")
#1 - Call MC-Schaefer function to preliminary explore the space without prior information
cat(stock, ": First Monte Carlo filtering of r-k space with ",n," points\n")
MCA <- Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
rv.all <- append(rv.all,MCA[[1]])
kv.all <- append(kv.all,MCA[[2]])
btv.all <- rbind(btv.all,MCA[[3]])
#take viable r and k values
nviablepoints = length(rv.all)
cat("* Found ",nviablepoints," viable points from ",n," samples\n");
#if few points were found then resample and shrink the k log space
if (nviablepoints<=1000){
log.start_k.new <- log(start_k)
max_attempts = 3
current_attempts = 1
while (nviablepoints<=1000 && current_attempts<=max_attempts){
if(nviablepoints > 0) {
log.start_k.new[1] <- mean(c(log.start_k.new[1], min(log(kv.all))))
log.start_k.new[2] <- mean(c(log.start_k.new[2], max(log(kv.all)))) }
n.new=n*current_attempts #add more points
ri1 = exp(runif(n.new, log(start_r[1]), log(start_r[2])))
ki1 = exp(runif(n.new, log.start_k.new[1], log.start_k.new[2]))
cat("Shrinking k space: repeating Monte Carlo in the interval [",exp(log.start_k.new[1]),",",exp(log.start_k.new[2]),"]\n")
cat("Attempt ",current_attempts," of ",max_attempts," with ",n.new," points","\n")
MCA <- Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
rv.all <- append(rv.all,MCA[[1]])
kv.all <- append(kv.all,MCA[[2]])
btv.all <- rbind(btv.all,MCA[[3]])
nviablepoints = length(rv.all) #recalculate viable points
cat("* Found altogether",nviablepoints," viable points \n");
current_attempts=current_attempts+1 #increment the number of attempts
}
}
# If tip of viable r-k pairs is 'thin', do extra sampling there
gm.rv = exp(mean(log(rv.all)))
if(length(rv.all[rv.all > 0.9*start_r[2]]) < 10) {
l.sample.r <- (gm.rv + max(rv.all))/2
cat("Final sampling in the tip area above r =",l.sample.r,"\n")
log.start_k.new <- c(log(0.8*min(kv.all)),log(max(kv.all[rv.all > l.sample.r])))
ri1 = exp(runif(50000, log(l.sample.r), log(start_r[2])))
ki1 = exp(runif(50000, log.start_k.new[1], log.start_k.new[2]))
MCA <- Schaefer(ri=ri1, ki=ki1, startbio=startbio, intyr=intyr, intbio=intbio, endbio=endbio, sigR=sigR, pt=T)
rv.all <- append(rv.all,MCA[[1]])
kv.all <- append(kv.all,MCA[[2]])
btv.all <- rbind(btv.all,MCA[[3]])
nviablepoints = length(rv.all) #recalculate viable points
cat("Found altogether", length(rv.all), "unique viable r-k pairs and biomass trajectories\n")
}
# ------------------------------------------------------------
# Bayesian analysis of catch & biomass with Schaefer model
# ------------------------------------------------------------
if(Btype == "observed" | Btype=="simulated") {
cat("Running Schaefer MCMC analysis....\n")
mcmc.burn <- as.integer(30000)
mcmc.chainLength <- as.integer(60000) # burn-in plus post-burn
mcmc.thin = 10 # to reduce autocorrelation
mcmc.chains = 3 # needs to be at least 2 for DIC
# Parameters to be returned by JAGS
jags.save.params=c('r','k','sigma.b', 'alpha', 'sigma.r') #
# JAGS model
Model = "model{
# to avoid crash due to 0 values
eps<-0.01
# set a quite narrow variation from the expected value
sigma.b <- 1/16
tau.b <- pow(sigma.b,-2)
Bm[1] <- log(alpha*k)
bio[1] ~ dlnorm(Bm[1],tau.b)
for (t in 2:nyr){
bio[t] ~ dlnorm(Bm[t],tau.b)
Bm[t] <- log(max(bio[t-1] + r*bio[t-1]*(1 - (bio[t-1])/k) - ct[t-1], eps))
}
# priors
alpha ~ dunif(0.01,1) # needed for fit of first biomass
#inverse cubic root relationship between the range of viable r and the size of the search space
inverseRangeFactor <- 1/((start_r[2]-start_r[1])^1/3)
# give sigma some variability in the inverse relationship
sigma.r ~ dunif(0.001*inverseRangeFactor,0.02*inverseRangeFactor)
tau.r <- pow(sigma.r,-2)
rm <- log((start_r[1]+start_r[2])/2)
r ~ dlnorm(rm,tau.r)
# search in the k space from the center of the range. Allow high variability
km <- log((start_k[1]+start_k[2])/2)
tau.k <- pow(km,-2)
k ~ dlnorm(km,tau.k)
#end model
}"
# Write JAGS model to file
cat(Model, file="r2jags.bug")
### random seed
set.seed(runif(1,1,500)) # needed in JAGS
### run model
jags_outputs <- jags(data=c('ct','bio','nyr', 'start_r', 'start_k'),
working.directory=NULL, inits=NULL,
parameters.to.save= jags.save.params,
model.file="r2jags.bug", n.chains = mcmc.chains,
n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength,
refresh=mcmc.burn/20, )
# ------------------------------------------------------
# Results from JAGS Schaefer
# ------------------------------------------------------
r_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$r))
k_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$k))
## sigma_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$sigma.b))
alpha_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$alpha))
## sigma.r_out <- as.numeric(mcmc(jags_outputs$BUGSoutput$sims.list$sigma.r))
mean.log.r.jags <- mean(log(r_out))
SD.log.r.jags <- sd(log(r_out))
lcl.log.r.jags <- mean.log.r.jags-1.96*SD.log.r.jags
ucl.log.r.jags <- mean.log.r.jags+1.96*SD.log.r.jags
gm.r.jags <- exp(mean.log.r.jags)
lcl.r.jags <- exp(lcl.log.r.jags)
ucl.r.jags <- exp(ucl.log.r.jags)
mean.log.k.jags <- mean(log(k_out))
SD.log.k.jags <- sd(log(k_out))
lcl.log.k.jags <- mean.log.k.jags-1.96*SD.log.k.jags
ucl.log.k.jags <- mean.log.k.jags+1.96*SD.log.k.jags
gm.k.jags <- exp(mean.log.k.jags)
lcl.k.jags <- exp(lcl.log.k.jags)
ucl.k.jags <- exp(ucl.log.k.jags)
mean.log.MSY.jags<- mean(log(r_out)+log(k_out)-log(4))
SD.log.MSY.jags <- sd(log(r_out)+log(k_out)-log(4))
gm.MSY.jags <- exp(mean.log.MSY.jags)
lcl.MSY.jags <- exp(mean.log.MSY.jags-1.96*SD.log.MSY.jags)
ucl.MSY.jags <- exp(mean.log.MSY.jags+1.96*SD.log.MSY.jags)
} # end of MCMC Schaefer loop
#------------------------------------
# get results from CMSY
#------------------------------------
# get estimate of most probable r as median of mid log.r-classes above cut-off
# get remaining viable log.r and log.k
rem.log.r <- log(rv.all[rv.all > gm.rv])
rem.log.k <- log(kv.all[rv.all>gm.rv])
# get vectors with numbers of r and mid values in about 25 classes
hist.log.r <- hist(x=rem.log.r, breaks=25, plot=F)
log.r.counts <- hist.log.r$counts
log.r.mids <- hist.log.r$mids
# get most probable log.r as mean of mids with counts > 0
log.r.est <- median(log.r.mids[which(log.r.counts > 0)])
lcl.log.r <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.025))
ucl.log.r <- as.numeric(quantile(x=log.r.mids[which(log.r.counts > 0)], 0.975))
r.est <- exp(log.r.est)
lcl.r.est <- exp(lcl.log.r)
ucl.r.est <- exp(ucl.log.r)
# do linear regression of log k ~ log r with slope fixed to -1 (from Schaefer)
reg <- lm(rem.log.k ~ 1 + offset(-1*rem.log.r))
int.reg <- as.numeric(reg[1])
sd.reg <- sd(resid(reg))
se.reg <- summary(reg)$coefficients[2]
# get estimate of log(k) from y where x = log.r.est
log.k.est <- int.reg + (-1) * log.r.est
# get estimates of CL of log.k.est from y +/- SD where x = lcl.log r or ucl.log.r
lcl.log.k <- int.reg + (-1) * ucl.log.r - sd.reg
ucl.log.k <- int.reg + (-1) * lcl.log.r + sd.reg
k.est <- exp(log.k.est)
lcl.k.est <- exp(lcl.log.k)
ucl.k.est <- exp(ucl.log.k)
# get MSY from remaining log r-k pairs
log.MSY.est <- mean(rem.log.r + rem.log.k - log(4))
sd.log.MSY.est <- sd(rem.log.r + rem.log.k - log(4))
lcl.log.MSY.est <- log.MSY.est - 1.96*sd.log.MSY.est
ucl.log.MSY.est <- log.MSY.est + 1.96*sd.log.MSY.est
MSY.est <- exp(log.MSY.est)
lcl.MSY.est <- exp(lcl.log.MSY.est)
ucl.MSY.est <- exp(ucl.log.MSY.est)
# get predicted biomass vectors as median and quantiles of trajectories
median.btv <- apply(btv.all,2, median)
lastyr.bio <- median.btv[length(median.btv)-1]
nextyr.bio <- median.btv[length(median.btv)]
lcl.btv <- apply(btv.all,2, quantile, probs=0.025)
q.btv <- apply(btv.all,2, quantile, probs=0.25)
ucl.btv <- apply(btv.all,2, quantile, probs=0.975)
lcl.lastyr.bio <- lcl.btv[length(lcl.btv)-1]
ucl.lastyr.bio <- ucl.btv[length(lcl.btv)-1]
lcl.nextyr.bio <- lcl.btv[length(lcl.btv)]
ucl.nextyr.bio <- ucl.btv[length(lcl.btv)]
# -----------------------------------------
# Plot results
# -----------------------------------------
# Analysis of viable r-k pairs
plot(x=rv.all, y=kv.all, xlim=start_r,
ylim=c(0.9*min(kv.all, ifelse(Btype == "observed",k_out,NA), na.rm=T), 1.1*max(kv.all)),
pch=16, col="grey",log="xy", bty="l",
xlab="r", ylab="k", main="Analysis of viable r-k")
abline(v=gm.rv, lty="dashed")
# plot points and best estimate from full Schaefer analysis
if(Btype == "observed"|Btype=="simulated") {
# plot r-k pairs from MCMC
points(x=r_out, y=k_out, pch=16,cex=0.5)
# plot best r-k pair from MCMC
points(x=gm.r.jags, y=gm.k.jags, pch=19, col="green")
lines(x=c(lcl.r.jags, ucl.r.jags),y=c(gm.k.jags,gm.k.jags), col="green")
lines(x=c(gm.r.jags,gm.r.jags),y=c(lcl.k.jags, ucl.k.jags), col="green")
}
# if data are from simulation, plot true r and k
if(Btype=="simulated") {
l.stock <- nchar(stock) # get length of sim stock name
r.char <- substr(stock,l.stock-1,l.stock) # get last character of sim stock name
r.sim <- NA # initialize vector for r used in simulation
if(r.char=="_H") {r.sim=1; lcl.r.sim=0.8; ucl.r.sim=1.25} else
if(r.char=="_M") {r.sim=0.5;lcl.r.sim=0.4;ucl.r.sim=0.62} else
if(r.char=="_L") {r.sim=0.25;lcl.r.sim=0.2;ucl.r.sim=0.31} else {r.sim=0.05;lcl.r.sim=0.04;ucl.r.sim=0.062}
# plot true r-k point with error bars
points(x=r.sim, y=1000, pch=19, col="red")
# add +/- 20% error bars
lines(x=c(lcl.r.sim,ucl.r.sim), y=c(1000,1000), col="red")
lines(x=c(r.sim,r.sim), y=c(800,1250), col="red")
}
# plot blue dot for proposed r-k, with 95% CL lines
points(x=r.est, y=k.est, pch=19, col="blue")
lines(x=c(lcl.r.est, ucl.r.est),y=c(k.est,k.est), col="blue")
lines(x=c(r.est,r.est),y=c(lcl.k.est, ucl.k.est), col="blue")
# plot biomass graph
# determine k to use for red line in b/k plot
if(Btype=="simulated") {k2use <- 1000} else
if(Btype == "observed") {k2use <- gm.k.jags} else {k2use <- k.est}
# determine hight of y-axis in plot
max.y <- max(c(bio/k2use,ucl.btv,0.6,startbio[2], intbio[2],endbio[2]),na.rm=T)
plot(x=yr,y=median.btv[1:nyr], lwd=2, xlab="Year", ylab="Relative biomass b/k", type="l",
ylim=c(0,max.y), bty="l", main=paste("Pred. biomass vs ", Btype,sep=""))
lines(x=yr, y=lcl.btv[1:nyr],type="l")
lines(x=yr, y=ucl.btv[1:nyr],type="l")
points(x=EndYear,y=q.btv[yr==EndYear], col="purple", cex=1.5, lwd=2)
abline(h=0.5, lty="dashed")
abline(h=0.25, lty="dotted")
lines(x=c(yr[1],yr[1]), y=startbio, col="blue")
lines(x=c(intyr,intyr), y=intbio, col="blue")
lines(x=c(max(yr),max(yr)), y=endbio, col="blue")
# if observed biomass is available, plot red biomass line
if(Btype == "observed"|Btype=="simulated") {
lines(x=yr, y=bio/k2use,type="l", col="red", lwd=1)
}
# if CPUE data are available, scale to predicted biomass range, plot red biomass line
if(Btype == "CPUE") {
par(new=T) # prepares for new plot on top of previous
plot(x=yr, y=bio, type="l", col="red", lwd=1,
ann=F,axes=F,ylim=c(0,1.2*max(bio, na.rm=T))) # forces this plot on top of previous one
axis(4, col="red", col.axis="red")
}
# plot yield and biomass against equilibrium surplus parabola
max.y <-max(c(ct/MSY.est,ifelse(Btype=="observed"|Btype=="simulated",ct/gm.MSY.jags,NA),1.2),na.rm=T)
# plot parabola
x=seq(from=0,to=2,by=0.001)
y=4*x-(2*x)^2
plot(x=x, y=y, xlim=c(0,1), ylim=c(0,max.y), type="l", bty="l",xlab="Relative biomass b/k",
ylab="Catch / MSY", main="Equilibrium curve")
# plot catch against CMSY biomass estimates
points(x=median.btv[1:nyr], y=ct/MSY.est, pch=16, col="grey")
points(x=q.btv[yr==EndYear],y=ct[yr==EndYear]/MSY.est, col="purple", cex=1.5, lwd=2)
# plot catch against observed biomass or CPUE
if(Btype == "observed"|Btype=="simulated") {
points(x=bio/k2use, y=ct/gm.MSY.jags, pch=16, cex=0.5)
}
# plot exploitation rate u against u.msy
# get u derived from predicted CMSY biomass
u.CMSY <- ct/(median.btv[1:nyr]*k.est)
u.msy.CMSY <- 1-exp(-r.est/2) # # Fmsy from CMSY expressed as exploitation rate
# get u from observed or simulated biomass
if(Btype == "observed"|Btype=="simulated") {
u.bio <- ct/bio
u.msy.bio <- 1-exp(-gm.r.jags/2)
}
# get u from CPUE
if(Btype == "CPUE") {
q=max(median.btv[1:nyr][is.na(bio)==F],na.rm=T)*k.est/max(bio,na.rm=T)
u.CPUE <- ct/(q*bio)
}
# determine upper bound of Y-axis
max.y <- max(c(1.5, 1.2*u.CMSY/u.msy.CMSY,ct[yr==EndYear]/(q.btv[yr==EndYear]*k.est)/u.msy.CMSY,
ifelse(Btype=="observed"|Btype=="simulated",max(u.bio[is.na(u.bio)==F]/u.msy.bio),0),
na.rm=T))
# plot u from CMSY
plot(x=yr,y=u.CMSY/u.msy.CMSY, type="l", bty="l", ylim=c(0,max.y), xlab="Year",
ylab="u / u_msy", main="Exploitation rate")
abline(h=1, lty="dashed")
points(x=EndYear,y=ct[yr==EndYear]/(q.btv[yr==EndYear]*k.est)/u.msy.CMSY, col="purple", cex=1.5, lwd=2)
# plot u from biomass
if(Btype == "observed"|Btype=="simulated") lines(x=yr, y=u.bio/u.msy.bio, col="red")
# plot u from CPUE
if(Btype == "CPUE") {
par(new=T) # prepares for new plot on top of previous
plot(x=yr, y=u.CPUE, type="l", col="red", ylim=c(0, 1.2*max(u.CPUE,na.rm=T)),ann=F,axes=F)
axis(4, col="red", col.axis="red")
}
if(batch.mode == TRUE) {dev.off()} # close plot window
# ------------------------------------------
# print input and results to screen
cat("---------------------------------------\n")
cat("Species:", cinfo$ScientificName[cinfo$stock==stock], "\n")
cat("Name and region:", cinfo$EnglishName[cinfo$stock==stock], ",", cinfo$Name[cinfo$stock==stock], "\n")
cat("Stock:",stock,"\n")
cat("Catch data used from years", min(yr),"-", max(yr), "\n")
cat("Prior initial relative biomass =", startbio[1], "-", startbio[2], "\n")
cat("Prior intermediate rel. biomass=", intbio[1], "-", intbio[2], "in year", intyr, "\n")
cat("Prior final relative biomass =", endbio[1], "-", endbio[2], "\n")
cat("If current catches continue, is the stock likely to crash within 3 years?",FutureCrash,"\n")
cat("Prior range for r =", format(start_r[1],digits=2), "-", format(start_r[2],digits=2),
", prior range for k =", start_k[1], "-", start_k[2],"\n")
# if data are simulated, print true r-k
if(filename_1=="SimCatch.csv") {
cat("True r =", r.sim, "(because input data were simulated with Schaefer model)\n")
cat("True k = 1000 \n")
cat("True MSY =", 1000*r.sim/4,"\n")
cat("True biomass in last year =",bio[length(bio)],"or",bio[length(bio)]/1000,"k \n")
cat("True mean catch / MSY ratio =", mean(ct)/(1000*r.sim/4),"\n")
}
# print results from full Schaefer if available
if(Btype == "observed"|Btype=="simulated") {
cat("Results from Bayesian Schaefer model using catch & biomass (",Btype,")\n")
cat("MSY =", gm.MSY.jags,", 95% CL =", lcl.MSY.jags, "-", ucl.MSY.jags,"\n")
cat("Mean catch / MSY =", mean(ct)/gm.MSY.jags,"\n")
if(Btype != "CPUE") {
cat("r =", gm.r.jags,", 95% CL =", lcl.r.jags, "-", ucl.r.jags,"\n")
cat("k =", gm.k.jags,", 95% CL =", lcl.k.jags, "-", ucl.k.jags,"\n")
}
}
# results of CMSY analysis
cat("Results of CMSY analysis \n")
cat("Altogether", nviablepoints,"unique viable r-k pairs were found \n")
cat(nviablepoints-length(rem.log.r),"r-k pairs above the initial geometric mean of r =", gm.rv, "were analysed\n")
cat("r =", r.est,", 95% CL =", lcl.r.est, "-", ucl.r.est,"\n")
cat("k =", k.est,", 95% CL =", lcl.k.est, "-", ucl.k.est,"\n")
cat("MSY =", MSY.est,", 95% CL =", lcl.MSY.est, "-", ucl.MSY.est,"\n")
cat("Predicted biomass in last year =", lastyr.bio, "2.5th perc =", lcl.lastyr.bio,
"97.5th perc =", ucl.lastyr.bio,"\n")
cat("Predicted biomass in next year =", nextyr.bio, "2.5th perc =", lcl.nextyr.bio,
"97.5th perc =", ucl.nextyr.bio,"\n")
cat("----------------------------------------------------------\n")
## Write some results into outfile
if(write.output == TRUE) {
# write data into csv file
output = data.frame(cinfo$ScientificName[cinfo$stock==stock], stock, StartYear, EndYear, mean(ct)*1000,
ifelse(Btype=="observed"|Btype=="simulate",bio[length(bio)],NA), # last biomass on record
ifelse(Btype == "observed"|Btype=="simulated",gm.MSY.jags,NA), # full Schaefer
ifelse(Btype == "observed"|Btype=="simulated",lcl.MSY.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",ucl.MSY.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",gm.r.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",lcl.r.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",ucl.r.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",gm.k.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",lcl.k.jags,NA),
ifelse(Btype == "observed"|Btype=="simulated",ucl.k.jags,NA),
r.est, lcl.r.est, ucl.r.est, # CMSY r
k.est, lcl.k.est, ucl.k.est, # CMSY k
MSY.est, lcl.MSY.est, ucl.MSY.est, # CMSY r
lastyr.bio, lcl.lastyr.bio, ucl.lastyr.bio, # last year bio
nextyr.bio, lcl.nextyr.bio, ucl.nextyr.bio)# last year + 1 bio
write.table(output, file=outfile, append = T, sep = ",",
dec = ".", row.names = FALSE, col.names = FALSE)
# write some text into text outfile.txt
cat("Species:", cinfo$ScientificName[cinfo$stock==stock], "\n",
"Name:", cinfo$EnglishName[cinfo$stock==stock], "\n",
"Region:", cinfo$Name[cinfo$stock==stock], "\n",
"Stock:",stock,"\n",
"Catch data used from years", min(yr),"-", max(yr),", biomass =", Btype, "\n",
"Prior initial relative biomass =", startbio[1], "-", startbio[2], "\n",
"Prior intermediate rel. biomass=", intbio[1], "-", intbio[2], "in year", intyr, "\n",
"Prior final relative biomass =", endbio[1], "-", endbio[2], "\n",
"Future crash with current catches?", FutureCrash, "\n",
"Prior range for r =", format(start_r[1],digits=2), "-", format(start_r[2],digits=2),
", prior range for k =", start_k[1], "-", start_k[2],"\n",
file=outfile.txt,append=T)
if(filename_1=="SimCatch.csv") {
cat(" True r =", r.sim, "(because input data were simulated with Schaefer model)\n",
"True k = 1000, true MSY =", 1000*r.sim/4,"\n",
"True biomass in last year =",bio[length(bio)],"or",bio[length(bio)]/1000,"k \n",
"True mean catch / MSY ratio =", mean(ct)/(1000*r.sim/4),"\n",
file=outfile.txt,append=T)
}
if(Btype == "observed"|Btype=="simulated") {
cat(" Results from Bayesian Schaefer model using catch & biomass \n",
"r =", gm.r.jags,", 95% CL =", lcl.r.jags, "-", ucl.r.jags,"\n",
"k =", gm.k.jags,", 95% CL =", lcl.k.jags, "-", ucl.k.jags,"\n",
"MSY =", gm.MSY.jags,", 95% CL =", lcl.MSY.jags, "-", ucl.MSY.jags,"\n",
"Mean catch / MSY =", mean(ct)/gm.MSY.jags,"\n",
file=outfile.txt,append=T)
}
cat(" Results of CMSY analysis with altogether", nviablepoints,"unique viable r-k pairs \n",
nviablepoints-length(rem.log.r),"r-k pairs above the initial geometric mean of r =", gm.rv, "were analysed\n",
"r =", r.est,", 95% CL =", lcl.r.est, "-", ucl.r.est,"\n",
"k =", k.est,", 95% CL =", lcl.k.est, "-", ucl.k.est,"\n",
"MSY =", MSY.est,", 95% CL =", lcl.MSY.est, "-", ucl.MSY.est,"\n",
"Predicted biomass last year b/k =", lastyr.bio, "2.5th perc b/k =", lcl.lastyr.bio,
"97.5th perc b/k =", ucl.lastyr.bio,"\n",
"Precautionary 25th percentile b/k =",q.btv[yr==EndYear],"\n",
"----------------------------------------------------------\n",
file=outfile.txt,append=T)
}
} # end of stocks loop

View File

@ -0,0 +1,119 @@
##--------------------------------------------------------
## CMSY analysis with estimation of total biomass, including Bayesian Schaefer
## written by Rainer Froese with support from Gianpaolo Coro in 2013-2014
## This version adjusts biomass to average biomass over the year
## It also contains the FutureCrash option to improve prediction of final biomass
## Version 21 adds the purple point to indicate the 25th percentile of final biomass
## Version 22 accepts that no biomass or CPUE area available
##--------------------------------------------------------
library(R2jags) # Interface with JAGS
library(coda)
#-----------------------------------------
# Some general settings
#-----------------------------------------
# set.seed(999) # use for comparing results between runs
rm(list=ls(all=TRUE)) # clear previous variables etc
options(digits=3) # displays all numbers with three significant digits as default
graphics.off() # close graphics windows from previous sessions
#-----------------------------------------
# General settings for the analysis
#-----------------------------------------
sigR <- 0.02 # overall process error; 0.05 works reasonable for simulations, 0.02 for real data; 0 if deterministic model
n <- 10000 # initial number of r-k pairs
batch.mode <- T # set to TRUE to suppress graphs
write.output <- T # set to true if table of output is wanted
FutureCrash <- "No"
#-----------------------------------------
# Start output to screen
#-----------------------------------------
cat("-------------------------------------------\n")
cat("Catch-MSY Analysis,", date(),"\n")
cat("-------------------------------------------\n")
#------------------------------------------
# Read data and assign to vectors
#------------------------------------------
# filename_1 <- "AllStocks_Catch4.csv"
# filename_2 <- "AllStocks_ID4.csv"
# filename_1 <- "SimCatch.csv"
# filename_2 <- "SimSpec.csv"
# filename_2 <- "SimSpecWrongS.csv"
# filename_2 <- "SimSpecWrongI.csv"
# filename_2 <- "SimSpecWrongF.csv"
# filename_2 <- "SimSpecWrongH.csv"
# filename_2 <- "SimSpecWrongL.csv"
# filename_1 <- "FishDataLim.csv"
# filename_2 <- "FishDataLimSpec.csv"
filename_1 <- "WKLIFE4Stocks.csv"
filename_2 <- "WKLIFE4ID.csv"
outfile<-"outfile"
outfile.txt <- "outputfile.txt"
cdat <- read.csv(filename_1, header=T, dec=".", stringsAsFactors = FALSE)
cinfo <- read.csv(filename_2, header=T, dec=".", stringsAsFactors = FALSE)
cat("Files", filename_1, ",", filename_2, "read successfully","\n")
# Stocks with total biomass data and catch data from StartYear to EndYear
# stocks <- sort(as.character(cinfo$stock)) # All stocks
stocks<-"HLH_M07"
# select one stock after the other
for(stock in stocks) {
# assign data from cinfo to vectors
res <- as.character(cinfo$Resilience[cinfo$stock==stock])
StartYear <- as.numeric(cinfo$StartYear[cinfo$stock==stock])
EndYear <- as.numeric(cinfo$EndYear[cinfo$stock==stock])
r_low <- as.numeric(cinfo$r_low[cinfo$stock==stock])
r_hi <- as.numeric(cinfo$r_hi[cinfo$stock==stock])
stb_low <- as.numeric(cinfo$stb_low[cinfo$stock==stock])
stb_hi <- as.numeric(cinfo$stb_hi[cinfo$stock==stock])
intyr <- as.numeric(cinfo$intyr[cinfo$stock==stock])
intbio_low <- as.numeric(cinfo$intbio_low[cinfo$stock==stock])
intbio_hi <- as.numeric(cinfo$intbio_hi[cinfo$stock==stock])
endbio_low <- as.numeric(cinfo$endbio_low[cinfo$stock==stock])
endbio_hi <- as.numeric(cinfo$endbio_hi[cinfo$stock==stock])
Btype <- as.character(cinfo$Btype[cinfo$stock==stock])
FutureCrash <- as.character(cinfo$FutureCrash[cinfo$stock==stock])
comment <- as.character(cinfo$comment[cinfo$stock==stock])
# extract data on stock
yr <- as.numeric(cdat$yr[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])
ct <- as.numeric(cdat$ct[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000 ## assumes that catch is given in tonnes, transforms to '000 tonnes
if(Btype=="observed" | Btype=="CPUE" | Btype=="simulated") {
bt <- as.numeric(cdat$TB[cdat$stock==stock & cdat$yr >= StartYear & cdat$yr <= EndYear])/1000 ## assumes that biomass is in tonnes, transforms to '000 tonnes
} else {bt <- NA}
nyr <- length(yr) # number of years in the time series
cat("->---------------------------------------
Species: NA
Name and region: NA , NA
Stock: HLH_M07
Catch data used from years 1 - 50
Prior initial relative biomass = 0.5 - 0.9
Prior intermediate rel. biomass= 0.01 - 0.4 in year 25
Prior final relative biomass = 0.4 - 0.8
If current catches continue, is the stock likely to crash within 3 years? No
Prior range for r = 0.2 - 0.8 , prior range for k = 125 - 9965
Results from Bayesian Schaefer model using catch & biomass ( simulated )
MSY = 91.7 , 95% CL = 83.9 - 100
Mean catch / MSY = 0.882
r = 0.425 , 95% CL = 0.374 - 0.483
k = 863 , 95% CL = 783 - 951
Results of CMSY analysis
Altogether 2055 unique viable r-k pairs were found
1142 r-k pairs above the initial geometric mean of r = 0.343 were analysed
r = 0.522 , 95% CL = 0.349 - 0.782
k = 683 , 95% CL = 438 - 1067
MSY = 89.2 , 95% CL = 82.2 - 96.7
Predicted biomass in last year = 0.676 2.5th perc = 0.435 97.5th perc = 0.768
Predicted biomass in next year = 0.673 2.5th perc = 0.433 97.5th perc = 0.758
----------------------------------------------------------
",file=outfile.txt,append=T)
}

View File

@ -0,0 +1,17 @@
<?xml version='1.0' encoding='UTF-8'?>
<hibernate-configuration>
<session-factory>
<property name="connection.driver_class">org.postgresql.Driver</property>
<property name="connection.provider_class">org.hibernate.connection.C3P0ConnectionProvider</property>
<property name="connection.url">jdbc:postgresql://localhost/testdb</property>
<property name="connection.username">gcube</property>
<property name="connection.password">d4science2</property>
<property name="dialect">org.hibernate.dialect.PostgreSQLDialect</property>
<property name="transaction.factory_class">org.hibernate.transaction.JDBCTransactionFactory</property>
<property name="c3p0.timeout">0</property>
<property name="c3p0.max_size">1</property>
<property name="c3p0.max_statements">0</property>
<property name="c3p0.min_size">1</property>
<property name="current_session_context_class">thread</property>
</session-factory>
</hibernate-configuration>

View File

@ -0,0 +1,530 @@
#### R and JAGS code for estimating LWR-parameters from previous studies
#### Meant for updating the ESTIMATE table in FishBase
#### Created by Rainer Froese in March 2013, including JAGS models by James Thorston
#### Modified in June 2013 to include subfamilies
rm(list=ls(all=TRUE)) # remove previous variables and data
options(digits=3) # 3 significant digits as default
library(R2jags) # Interface with JAGS
runif(1) # sets random seed
#### Read in data
DataFile = "RF_LWR2.csv" # RF_LWR4 was extracted from FishBase in June 2013
Data = read.csv(DataFile, header=TRUE)
cat("Start", date(), "\n")
cat("Data file =", DataFile, "\n")
# Get unique, sorted list of Families
Fam.All <- sort(unique(as.character(Data$Family)))
Families <- Fam.All[Fam.All== "Acanthuridae" | Fam.All == "Achiridae"]
OutFile = "LWR_Test1.csv"
JAGSFILE = "dmnorm_0.bug"
# Get unique, sorted list of body shapes
Bshape <- sort(unique(as.character(Data$BodyShapeI)))
#------------------------------------------
# Functions
#------------------------------------------
#---------------------------------------------------------
# Function to get the priors for the respective body shape
#---------------------------------------------------------
Get.BS.pr <- function(BS) {
### Assignment of priors based on available body shape information
# priors derived from 5150 LWR studies in FishBase 02/2013
if (BS == "eel-like") { # eel-like prior for log(a) and b
prior_mean_log10a = -2.99
prior_sd_log10a = 0.175
prior_tau_log10a = 1/prior_sd_log10a^2
prior_mean_b = 3.06
prior_sd_b = 0.0896
prior_tau_b = 1/prior_sd_b^2
} else
if (BS == "elongated") { # elongate prior for log(a) and b
prior_mean_log10a = -2.41
prior_sd_log10a = 0.171
prior_tau_log10a = 1/prior_sd_log10a^2
prior_mean_b = 3.12
prior_sd_b = 0.09
prior_tau_b = 1/prior_sd_b^2
} else
if (BS == "fusiform / normal") { # fusiform / normal prior for log(a) and b
prior_mean_log10a = -1.95
prior_sd_log10a = 0.173
prior_tau_log10a = 1/prior_sd_log10a^2
prior_mean_b = 3.04
prior_sd_b = 0.0857
prior_tau_b = 1/prior_sd_b^2
} else
if (BS == "short and / or deep") { # short and / or deep prior for log(a) and b
prior_mean_log10a = -1.7
prior_sd_log10a = 0.175
prior_tau_log10a = 1/prior_sd_log10a^2
prior_mean_b = 3.01
prior_sd_b = 0.0905
prior_tau_b = 1/prior_sd_b^2
} else
# priors across all shapes, used for missing or other BS
{
prior_mean_log10a = -2.0
prior_sd_log10a = 0.313
prior_tau_log10a = 1/prior_sd_log10a^2
prior_mean_b = 3.04
prior_sd_b = 0.119
prior_tau_b = 1/prior_sd_b^2
}
# Priors for measurement error (= sigma) based on 5150 studies
# given here as shape mu and rate r, for gamma distribution
SD_rObs_log10a = 6520
SD_muObs_log10a = 25076
SD_rObs_b = 6808
SD_muObs_b = 37001
# Priors for between species variability (= sigma) based on 5150 studies for 1821 species
SD_rGS_log10a = 1372
SD_muGS_log10a = 7933
SD_rGS_b = 572
SD_muGS_b = 6498
prior.list <- list(mean_log10a=prior_mean_log10a, sd_log10a=prior_sd_log10a,
tau_log10a=prior_tau_log10a, mean_b=prior_mean_b, sd_b=prior_sd_b,
tau_b=prior_tau_b, SD_rObs_log10a=SD_rObs_log10a, SD_muObs_log10a=SD_muObs_log10a,
SD_rObs_b=SD_rObs_b, SD_muObs_b=SD_muObs_b, SD_rGS_log10a=SD_rGS_log10a,
SD_muGS_log10a=SD_muGS_log10a, SD_rGS_b=SD_rGS_b, SD_muGS_b=SD_muGS_b)
return(prior.list)
}
#--------------------------------------------------------------------
# Function to do a Bayesian analysis including LWR from relatives
#--------------------------------------------------------------------
SpecRelLWR <- function(a, b, wts, GenusSpecies, Nspecies, prior_mean_b, prior_tau_b,
prior_mean_log10a, prior_tau_log10a, SD_rObs_log10a, SD_muObs_log10a,
SD_rObs_b, SD_muObs_b, SD_rGS_log10a, SD_muGS_log10a,
SD_rGS_b, SD_muGS_b){
### Define JAGS model
Model = "
model {
#### Process model -- effects of taxonomy
# given the likelihood distributions and the priors,
# create normal posterior distributions for log10a, b,
# and for the process error (=between species variability sigmaGS)
abTrue[1] ~ dnorm(prior_mean_log10a,prior_tau_log10a)
abTrue[2] ~ dnorm(prior_mean_b,prior_tau_b)
sigmaGSlog10a ~ dgamma( SD_rGS_log10a, SD_muGS_log10a)
sigmaGSb ~ dgamma( SD_rGS_b, SD_muGS_b)
# given the posterior distributions and the process errors,
# establish for every species the expected witin-species
# parameter distributions; no correlation roGS between species
roGS <- 0
tauGenusSpecies[1] <- pow(sigmaGSlog10a,-2)
tauGenusSpecies[2] <- pow(sigmaGSb,-2)
for(k in 1:Nspecies){
abGenusSpecies[k,1] ~ dnorm(abTrue[1],tauGenusSpecies[1])
abGenusSpecies[k,2] ~ dnorm(abTrue[2],tauGenusSpecies[2])
}
### Observation model
## Errors
# given the data and the priors, establish distributions
# for the observation errors sigmaObs
sigmaObslog10a ~ dgamma( SD_rObs_log10a, SD_muObs_log10a)
sigmaObsb ~ dgamma( SD_rObs_b, SD_muObs_b)
# create inverse covariance matrix, with negative parameter correlation roObs
roObs ~ dunif(-0.99,0)
CovObs[1,1] <- pow(sigmaObslog10a,2)
CovObs[2,2] <- pow(sigmaObsb,2)
CovObs[1,2] <- roObs * sigmaObslog10a * sigmaObsb
CovObs[2,1] <- CovObs[1,2]
TauObs[1:2,1:2] <- inverse(CovObs[1:2,1:2])
## likelihood
# given the data, the priors and the covariance,
# create multivariate likelihood distributions for log10(a) and b
for(i in 1:N){
TauObsI[i,1:2,1:2] <- TauObs[1:2,1:2] * pow(Weights[i],2) # weighted precision
ab[i,1:2] ~ dmnorm(abGenusSpecies[GenusSpecies[i],1:2],TauObsI[i,1:2,1:2])
}
}
"
# Write JAGS model
cat(Model, file=JAGSFILE)
# JAGS settings
Nchains = 3 # number of MCMC chains to be used in JAGS
Nburnin = 1e4 # number of burn-in iterations, to be discarded; 1e4 = 10000 iterations for burn-in
Niter = 3e4 # number of iterations after burn-in; 3e4 = 30000 iterations
Nthin = 1e1 # subset of iterations to be used for analysis; 1e1 = every 10th iteration
# Run JAGS: define data to be passed on in DataJags;
# determine parameters to be returned in Param2Save;
# call JAGS with function Jags()
DataJags = list(ab=cbind(log10(a),b), N=length(a), Weights=wts, Nspecies=Nspecies, GenusSpecies=GenusSpecies,
prior_mean_b=prior_mean_b, prior_tau_b=prior_tau_b,
prior_mean_log10a=prior_mean_log10a, prior_tau_log10a=prior_tau_log10a,
SD_rObs_log10a=SD_rObs_log10a, SD_muObs_log10a=SD_muObs_log10a,
SD_rObs_b=SD_rObs_b, SD_muObs_b=SD_muObs_b,
SD_rGS_log10a=SD_rGS_log10a, SD_muGS_log10a=SD_muGS_log10a,
SD_rGS_b=SD_rGS_b, SD_muGS_b=SD_muGS_b)
Params2Save = c("abTrue","abGenusSpecies","sigmaGSlog10a","sigmaGSb","sigmaObslog10a","sigmaObsb","roObs")
Jags <- jags(inits=NULL, model.file=JAGSFILE, working.directory=NULL, data=DataJags,
parameters.to.save=Params2Save, n.chains=Nchains, n.thin=Nthin, n.iter=Niter, n.burnin=Nburnin)
Jags$BUGSoutput # contains the results from the JAGS run
# Analyze output for the relatives
abTrue <- Jags$BUGSoutput$sims.list$abTrue
R_mean_log10a <- mean(abTrue[,1]) # true mean of log10(a)
R_sd_log10a <- sd(abTrue[,1]) # true SE of log10(a)
R_mean_b <- mean(abTrue[,2]) # true mean of b
R_sd_b <- sd(abTrue[,2]) # true SE of b
# Analyze output for the target species
abGenusSpecies <- Jags$BUGSoutput$sims.list$abGenusSpecies
mean_log10a <- mean(abGenusSpecies[,1,1]) # true mean of log10(a) for the first species= target species
sd_log10a <- sd(abGenusSpecies[,1,1]) # true SE of log10(a)
mean_b <- mean(abGenusSpecies[,1,2]) # true mean of b
sd_b <- sd(abGenusSpecies[,1,2]) # true SE of b
mean_sigma_log10a <- mean(Jags$BUGSoutput$sims.list$sigmaObslog10a) # measurement error of log10(a)
sd_sigma_log10a <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigmaObslog10a), 2, sd)
mean_sigma_b <- mean(Jags$BUGSoutput$sims.list$sigmaObsb) # measurement error of b
sd_sigma_b <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigmaObsb), 2, sd)
ro_ab <- mean(Jags$BUGSoutput$sims.list$roObs) # measurement correlation of log10(a),b
out.list <- list(N=length(a), mean_log10a=mean_log10a, sd_log10a=sd_log10a, mean_b=mean_b, sd_b=sd_b,
R_mean_log10a=R_mean_log10a, R_sd_log10a=R_sd_log10a, R_mean_b=R_mean_b, R_sd_b=R_sd_b)
return(out.list)
}
#-----------------------------------------------------------------------------
# Function to do a Bayesian LWR analysis with studies for target species only
#-----------------------------------------------------------------------------
SpecLWR <- function(a, b, wts, prior_mean_b, prior_tau_b,
prior_mean_log10a, prior_tau_log10a, SD_rObs_log10a, SD_muObs_log10a,
SD_rObs_b, SD_muObs_b, SD_rGS_log10a, SD_muGS_log10a,
SD_rGS_b, SD_muGS_b){
# Define JAGS model
Model = "
model {
sigma1 ~ dgamma( SD_rObs_log10a, SD_muObs_log10a) # posterior distribution for measurement error in log10a
sigma2 ~ dgamma( SD_rObs_b, SD_muObs_b) # posterior distribution for measurement error in log10a
ro ~ dunif(-0.99,0) # uniform prior for negative correlation between log10a and b
abTrue[1] ~ dnorm(prior_mean_log10a,prior_tau_log10a) # normal posterior distribution for log10a
abTrue[2] ~ dnorm(prior_mean_b,prior_tau_b) # normal posterior distribution for b
CovObs[1,1] <- pow(sigma1,2)
CovObs[2,2] <- pow(sigma2,2)
CovObs[1,2] <- ro * sigma1 * sigma2
CovObs[2,1] <- CovObs[1,2]
TauObs[1:2,1:2] <- inverse(CovObs[1:2,1:2]) # create inverse covariance matrix
for(i in 1:N){
TauObsI[i,1:2,1:2] <- TauObs[1:2,1:2] * pow(Weights[i],2) # converts prior SD into prior weighted precision
# given the data, the priors and the covariance, create multivariate normal posteriors for log(a) and b
ab[i,1:2] ~ dmnorm(abTrue[1:2],TauObsI[i,1:2,1:2])
}
}
"
# Write JAGS model
cat(Model, file=JAGSFILE)
# JAGS settings
Nchains = 3 # number of MCMC chains to be used in JAGS
Nburnin = 1e4 # number of burn-in runs, to be discarded; 10000 iterations for burn-in
Niter = 3e4 # number of iterations after burn-in; 3e4 = 30000 iterations
Nthin = 1e1 # subset of iterations to be used for analysis; 1e1 = every 10th iteration
# Run JAGS: define data to be passed on in DataJags; determine parameters to be returned in Param2Save; call JAGS with function Jags()
DataJags = list(ab=cbind(log10(a),b), N=length(a), Weights=wts, prior_mean_b=prior_mean_b,
prior_tau_b=prior_tau_b, prior_mean_log10a=prior_mean_log10a, prior_tau_log10a=prior_tau_log10a,
SD_rObs_log10a=SD_rObs_log10a, SD_muObs_log10a=SD_muObs_log10a,
SD_rObs_b=SD_rObs_b, SD_muObs_b=SD_muObs_b)
Params2Save = c("abTrue","sigma1","sigma2","ro")
Jags <- jags(inits=NULL, model.file=JAGSFILE, working.directory=NULL, data=DataJags, parameters.to.save=Params2Save, n.chains=Nchains, n.thin=Nthin, n.iter=Niter, n.burnin=Nburnin)
Jags$BUGSoutput # contains the results from the JAGS run
# Analyze output
abTrue <- Jags$BUGSoutput$sims.list$abTrue
mean_log10a <- mean(abTrue[,1]) # true mean of log10(a)
sd_log10a <- sd(abTrue[,1]) # true SE of log10(a)
mean_b <- mean(abTrue[,2]) # true mean of b
sd_b <- sd(abTrue[,2]) # true SE of b
mean_sigma_log10a <- mean(Jags$BUGSoutput$sims.list$sigma1) # measurement error of log10(a)
sd_sigma_log10a <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigma1), 2, sd)
mean_sigma_b <- mean(Jags$BUGSoutput$sims.list$sigma2) # measurement error of b
sd_sigma_b <- apply(as.matrix(Jags$BUGSoutput$sims.list$sigma2), 2, sd)
ro_ab <- mean(Jags$BUGSoutput$sims.list$ro) # measurement correlation of log10(a),b
out.list <- list(N=length(a), mean_log10a=mean_log10a, sd_log10a=sd_log10a, mean_b=mean_b, sd_b=sd_b)
return(out.list)
} # End of Functions section
#--------------------------------
# Analysis by Family
#--------------------------------
# Do LWR analysis by Family, Subfamily and Body shape, depending on available LWR studies
# for(Fam in "Acanthuridae") {
for(Fam in Families) {
Subfamilies <- sort(unique(Data$Subfamily[Data$Family==Fam]))
for(SF in Subfamilies) {
for(BS in Bshape) {
# get species (SpecCodes) in this Subfamily and with this body shape
SpecCode.SF.BS <- unique(Data$SpecCode[Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS])
# if there are species with this body shape
if(length(SpecCode.SF.BS > 0)) {
# get priors for this body shape
prior <- Get.BS.pr(BS)
# get LWR data for this body shape
b_raw <- Data$b[Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS]
cat("\n")
cat("Family =", Fam, ", Subfamily =", SF, ", Body shape =", BS, ", Species =", length(SpecCode.SF.BS), ", LWR =",
length(b_raw[is.na(b_raw)==F]), "\n")
# if no LWR studies exist for this body shape, assign the respective priors to all species
if(length(b_raw[is.na(b_raw)==F])==0) {
# assign priors to species with no LWR in this Subfamily with this body shape
cat("Assigning overall body shape prior to", length(SpecCode.SF.BS), " species \n")
for(SpC in SpecCode.SF.BS) {
out.prior <- data.frame(Fam, SF, BS, SpC, 0, prior$mean_log10a, prior$sd_log10a, prior$mean_b, prior$sd_b,
paste("all LWR estimates for this BS"))
write.table(out.prior, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F)
}
} else {
# Update priors for this body shape using existing LWR studies
# get LWR data for this Subfamily and body shape
Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0)
wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses)
a <- Data$a[Keep]
b <- Data$b[Keep]
GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep])
# add a first dummy record with prior LWR and low score = 0.3, as pseudo target species
# Name of dummy target species is Dum1 dum1
TargetSpec = paste("Dum1", "dum1")
wts <- c(0.3, wts)
a <- c(10^(prior$mean_log10a), a)
b <- c(prior$mean_b, b)
GenSpec <- c(TargetSpec, GenSpec)
# Relabel GenSpec so that TargetSpec = level 1
OtherSpecies = unique(GenSpec[GenSpec != TargetSpec])
GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies))
Nspecies = nlevels(GenusSpecies) # number of species
# run Bayesian analysis for pseudo target species with Subfamily members
# The resulting R_mean_log10a, R_sd_log10a, R_mean_b, R_sd_b will be used for species without LWR
cat("Updating Subfamily-Bodyshape prior using", Nspecies-1, "species with LWR studies \n")
prior.SFam.BS <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies, prior_mean_b=prior$mean_b,
prior_tau_b=prior$tau_b, prior_mean_log10a=prior$mean_log10a,
prior_tau_log10a=prior$tau_log10a, SD_rObs_log10a=prior$SD_rObs_log10a,
SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b,
SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a,
SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b,
SD_muGS_b=prior$SD_muGS_b)
#------------------------------------------------------------------------------------------
# if there are Genera with >= 5 species with LWR, update body shape priors for these Genera
#------------------------------------------------------------------------------------------
Genera <- unique(as.character(Data$Genus[Keep]))
# create empty list of lists for storage of generic priors
prior.Gen.BS <- rep(list(list()),length(Genera)) # create a list of empty lists
names(prior.Gen.BS) <- Genera # name the list elements according to the Genera
for(Genus in Genera){
# check if Genus contains >= 5 species with LWR data
if(length(unique(Data$SpecCode[Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F &
Data$Score>0 & Data$Genus==Genus]))>=5) {
# run Subfamily analysis with only data for this genus
Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0 &
Data$Genus==Genus)
wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses)
a <- Data$a[Keep]
b <- Data$b[Keep]
GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep])
# add a first dummy record with prior LWR and low score = 0.3, as pseudo target species
# Name of dummy target species is Dum1 dum1
TargetSpec = paste("Dum1", "dum1")
wts <- c(0.3, wts)
a <- c(10^(prior$mean_log10a), a)
b <- c(prior$mean_b, b)
GenSpec <- c(TargetSpec, GenSpec)
# Relabel GenSpec so that TargetSpec = level 1
OtherSpecies = unique(GenSpec[GenSpec != TargetSpec])
GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies))
Nspecies = nlevels(GenusSpecies) # number of species
# run Bayesian analysis for pseudo target species with Genus members
# R_mean_log10a, R_sd_log10a, R_mean_b, R_sd_b will be used for species without LWR
cat("Updating prior for Genus =", Genus, ", with", Nspecies -1, "LWR Species \n")
prior.Gen.BS[[Genus]] <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies,
prior_mean_b=prior.SFam.BS$R_mean_b,
prior_tau_b=1/prior.SFam.BS$R_sd_b^2,
prior_mean_log10a=prior.SFam.BS$R_mean_log10a,
prior_tau_log10a=1/prior.SFam.BS$R_sd_log10a,
SD_rObs_log10a=prior$SD_rObs_log10a,
SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b,
SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a,
SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b,
SD_muGS_b=prior$SD_muGS_b)
}
}
# new Subfamily-BS priors have been generated
# for some genera, new Genus-BS priors have been generated
# ---------------------------------------------------------------------
# Loop through all species in this Subfamily-BS; assign LWR as appropriate
# ---------------------------------------------------------------------
for(SpC in SpecCode.SF.BS) {
Genus <- as.character(unique(Data$Genus[Data$SpecCode==SpC]))
Species <- as.character(unique(Data$Species[Data$SpecCode==SpC]))
TargetSpec = paste(Genus, Species)
LWR <- length(Data$b[Data$SpecCode==SpC & is.na(Data$b)==F & Data$Score>0])
LWRGenspec <- length(unique(Data$SpecCode[Data$BodyShapeI==BS & is.na(Data$b)==F &
Data$Score>0 & Data$Genus==Genus]))
LWRSFamspec <- length(unique(Data$SpecCode[Data$BodyShapeI==BS & is.na(Data$b)==F &
Data$Score>0 & Data$Family==Fam & Data$Subfamily==SF]))
#---------------------------------------------------------
# >= 5 LWR in target species, run single species analysis
if(LWR >= 5) {
# Run analysis with data only for this species
Keep <- which(Data$SpecCode==SpC & is.na(Data$b)==F & Data$Score>0)
wts = Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses)
a = Data$a[Keep]
b = Data$b[Keep]
# determine priors to be used
if(LWRGenspec >= 5) {
prior_mean_b=prior.Gen.BS[[Genus]]$R_mean_b
prior_tau_b=1/prior.Gen.BS[[Genus]]$R_sd_b^2
prior_mean_log10a=prior.Gen.BS[[Genus]]$R_mean_log10a
prior_tau_log10a=1/prior.Gen.BS[[Genus]]$R_sd_log10a^2
} else
if (LWRSFamspec > 0) {
prior_mean_b=prior.SFam.BS$R_mean_b
prior_tau_b=1/prior.SFam.BS$R_sd_b^2
prior_mean_log10a=prior.SFam.BS$R_mean_log10a
prior_tau_log10a=1/prior.SFam.BS$R_sd_log10a^2
} else {
prior_mean_b=prior$mean_b
prior_tau_b=prior$tau_b
prior_mean_log10a=prior$mean_log10a
prior_tau_log10a=prior$tau_log10a
}
cat("Running single species analysis for", TargetSpec, "LWR =", LWR, ", LWR species in Genus=",LWRGenspec,"\n" )
# call function for single species analysis
post <- SpecLWR(a, b, wts, prior_mean_b=prior_mean_b,
prior_tau_b=prior_tau_b, prior_mean_log10a=prior_mean_log10a,
prior_tau_log10a=prior_tau_log10a, SD_rObs_log10a=prior$SD_rObs_log10a,
SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b,
SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a,
SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b,
SD_muGS_b=prior$SD_muGS_b)
out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(post$mean_log10a, digits=3), format(post$sd_log10a, digits=3), format(post$mean_b, disgits=3), format(post$sd_b, digits=3),
paste("LWR estimates for this species"))
write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F)
} else
#--------------------------------------------------------
# 1-4 LWR in target species and >= 5 LWR species in Genus
# run hierarchical analysis for genus members, with Subfamily-BS prior
if(LWR >= 1 & LWRGenspec >=5) {
# run Subfamily analysis with only data for this genus
Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0 &
Data$Genus==Genus)
wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses)
a <- Data$a[Keep]
b <- Data$b[Keep]
GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep])
# Relabel GenSpec so that TargetSpec = level 1
OtherSpecies = unique(GenSpec[GenSpec != TargetSpec])
GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies))
Nspecies = nlevels(GenusSpecies) # number of species
# run Bayesian analysis for target species with Genus members
cat("Running analysis with congeners for", TargetSpec, ", LWR =", LWR,", LWR species in Genus =", LWRGenspec,"\n")
post <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies,
prior_mean_b=prior.SFam.BS$R_mean_b,
prior_tau_b=1/prior.SFam.BS$R_sd_b^2,
prior_mean_log10a=prior.SFam.BS$R_mean_log10a,
prior_tau_log10a=1/prior.SFam.BS$R_sd_log10a^2,
SD_rObs_log10a=prior$SD_rObs_log10a,
SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b,
SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a,
SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b,
SD_muGS_b=prior$SD_muGS_b)
out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(post$mean_log10a, digits=3), format(post$sd_log10a, digits=3), format(post$mean_b, disgits=3), format(post$sd_b, digits=3),
paste("LWR estimates for species & Genus-BS"))
write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F)
} else
#-------------------------------------------------------
# 1-4 LWR in target species and < 5 LWR species in Genus
# run hierarchical analysis for Subfamily members, with bodyshape prior
if(LWR >= 1 & LWRSFamspec > 1) {
# run Subfamily analysis
Keep <- which(Data$Family==Fam & Data$Subfamily==SF & Data$BodyShapeI==BS & is.na(Data$b)==F & Data$Score>0)
wts <- Data$Score[Keep] # Un-normalized weights (so that Cov is comparable among analyses)
a <- Data$a[Keep]
b <- Data$b[Keep]
GenSpec <- paste(Data$Genus[Keep],Data$Species[Keep])
# Relabel GenSpec so that TargetSpec = level 1
OtherSpecies = unique(GenSpec[GenSpec != TargetSpec])
GenusSpecies = factor(GenSpec, levels=c(TargetSpec, OtherSpecies))
Nspecies = nlevels(GenusSpecies) # number of species
# run Bayesian analysis for target species with Subfamily members
cat("Running analysis with Subfamily members for", TargetSpec, ", LWR =", LWR,", LWR species in Subfamily-BS =",
LWRSFamspec, "\n")
post <- SpecRelLWR(a, b, wts, GenusSpecies, Nspecies,
prior_mean_b=prior$mean_b,
prior_tau_b=prior$tau_b,
prior_mean_log10a=prior$mean_log10a,
prior_tau_log10a=prior$tau_log10a,
SD_rObs_log10a=prior$SD_rObs_log10a,
SD_muObs_log10a=prior$SD_muObs_log10a, SD_rObs_b=prior$SD_rObs_b,
SD_muObs_b=prior$SD_muObs_b, SD_rGS_log10a=prior$SD_rGS_log10a,
SD_muGS_log10a=prior$SD_muGS_log10a, SD_rGS_b=prior$SD_rGS_b,
SD_muGS_b=prior$SD_muGS_b)
out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(post$mean_log10a, digits=3), format(post$sd_log10a, digits=3),
format(post$mean_b, disgits=3), format(post$sd_b, digits=3),
paste("LWR estimates for species & Subfamily-BS"))
write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F)
} else
#--------------------------------------------------
# assign Genus-BS priors to target species
if(LWRGenspec >= 5) {
cat("Assign Genus-BS prior for", TargetSpec, "\n")
out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(prior.Gen.BS[[Genus]]$mean_log10a, digits=3),
format(prior.Gen.BS[[Genus]]$sd_log10a, digits=3),
format(prior.Gen.BS[[Genus]]$mean_b, digits=3), format(prior.Gen.BS[[Genus]]$sd_b, digits=3),
paste("LWR estimates for this Genus-BS"))
write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F)
} else {
# -----------------------------------------------
# assign Subfamily-BS priors to target species
cat("Assign Subfamily-BS prior for", TargetSpec,"\n")
out.SpC <- data.frame(Fam, SF, BS, SpC, LWR, format(prior.SFam.BS$mean_log10a, digits=3), format(prior.SFam.BS$sd_log10a, digits=3),
format(prior.SFam.BS$mean_b, digits=3), format(prior.SFam.BS$sd_b, digits=3), paste("LWR estimates for this Subfamily-BS"))
write.table(out.SpC, file=OutFile, append = T, sep=",", dec=".", row.names=F, col.names=F)
}
} # end of species loop for this Subfamily and body shape
} # end of section dealing with Subfamily - body shapes that contain LWR estimates
} # end of section that deals with Subfamily - body shapes that contain species
} # end of body shape section
} # end of Subfamily section
} # end of Family section
cat("End", date(),"\n")

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.