git-svn-id: https://svn.d4science.research-infrastructures.eu/gcube/trunk/data-analysis/EcologicalEngineSmartExecutor@115361 82a268e6-3cf1-43bd-a215-b396298e98cf
This commit is contained in:
parent
0f5b73a868
commit
de93d14fbf
|
|
@ -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
|
||||
|
|
@ -0,0 +1,435 @@
|
|||
set.seed(999) ## for same random sequence
|
||||
#require(hacks)
|
||||
#13/05/2015
|
||||
#setwd("C:/Users/Ye/Documents/Data poor fisheries/Martell Froese Method/")
|
||||
|
||||
## Read Data for stock, year=yr, catch=ct, and resilience=res. Expects space delimited file with header yr ct and years in integer and catch in real with decimal point
|
||||
## For example
|
||||
## stock res yr ct
|
||||
## cap-icel Medium 1984 1234.32
|
||||
|
||||
## filename <- "RAM_MSY.csv"
|
||||
##filename <- "ICESct2.csv"
|
||||
|
||||
cat("Step 1","\n")
|
||||
TestRUN <- F # if it is true, just run on the test samples, false will go for a formal run!
|
||||
|
||||
filename <- "D20.csv"
|
||||
outfile <- "CatchMSY_Output.csv"
|
||||
outfile2 <- paste("NonProcessedSpecies.csv",sep="")
|
||||
|
||||
#cdat <- read.csv2(filename, header=T, dec=".")
|
||||
cdat1 <- read.csv(filename)
|
||||
cat("\n", "File", filename, "read successfully","\n")
|
||||
|
||||
|
||||
cat("Step 2","\n")
|
||||
if(file.exists("cdat.RData"))
|
||||
{load("cdat.RData")} else
|
||||
{
|
||||
|
||||
dim(cdat1)
|
||||
yrs=1950:2013
|
||||
|
||||
# to set NA as 0
|
||||
cdat1[is.na(cdat1)] <- 0
|
||||
nrow <- length(cdat1[,1])
|
||||
ndatColn <- length(cdat1[1,c(-1:-12)])
|
||||
rownames(cdat1) <- NULL
|
||||
|
||||
cdat <- NULL
|
||||
|
||||
for(i in 1:nrow)
|
||||
#for(i in 1:5)
|
||||
|
||||
{#i=1
|
||||
#a <- ctotal3[i,-1]
|
||||
tmp=data.frame(stock=rep(as.character(cdat1[i,"Stock_ID"]),ndatColn),
|
||||
species=rep(as.character(cdat1[i,"Scientific_name"]),ndatColn),
|
||||
yr=yrs,ct=unlist(c(cdat1[i,-c(1:12)])),
|
||||
res=rep(cdat1[i,"ResilienceIndex"],ndatColn))
|
||||
|
||||
cdat <- rbind(cdat,tmp)
|
||||
#edit(cdat)
|
||||
}
|
||||
save(cdat,file="cdat.RData")
|
||||
}
|
||||
|
||||
StockList=unique(as.character(cdat$stock))
|
||||
|
||||
cat("Step 3","\n")
|
||||
## FUNCTIONS are going to be used subsequently
|
||||
.schaefer <- function(theta)
|
||||
{
|
||||
with(as.list(theta), { ## for all combinations of ri & ki
|
||||
bt=vector()
|
||||
ell = 0 ## initialize ell
|
||||
J=0 #Ye
|
||||
for (j in startbt)
|
||||
{
|
||||
if(ell == 0)
|
||||
{
|
||||
bt[1]=j*k*exp(rnorm(1,0, sigR)) ## set biomass in first year
|
||||
for(i in 1:nyr) ## for all years in the time series
|
||||
{
|
||||
xt=rnorm(1,0, sigR)
|
||||
bt[i+1]=(bt[i]+r*bt[i]*(1-bt[i]/k)-ct[i])*exp(xt)
|
||||
## calculate biomass as function of previous year's biomass plus net production minus catch
|
||||
}
|
||||
|
||||
#Bernoulli likelihood, assign 0 or 1 to each combination of r and k
|
||||
ell = 0
|
||||
if(bt[nyr+1]/k>=lam1 && bt[nyr+1]/k <=lam2 && min(bt) > 0 && max(bt) <=k && bt[which(yr==interyr)]/k>=interbio[1] && bt[which(yr==interyr)]/k<=interbio[2])
|
||||
ell = 1
|
||||
J=j # Ye
|
||||
}
|
||||
}
|
||||
return(list(ell=ell,J=J)) # Ye adding J=J
|
||||
|
||||
|
||||
})
|
||||
}
|
||||
|
||||
sraMSY <-function(theta, N)
|
||||
{
|
||||
#This function conducts the stock reduction
|
||||
#analysis for N trials
|
||||
#args:
|
||||
# theta - a list object containing:
|
||||
# r (lower and upper bounds for r)
|
||||
# k (lower and upper bounds for k)
|
||||
# lambda (limits for current depletion)
|
||||
|
||||
|
||||
with(as.list(theta),
|
||||
{
|
||||
ri = exp(runif(N, log(r[1]), log(r[2]))) ## get N values between r[1] and r[2], assign to ri
|
||||
ki = exp(runif(N, log(k[1]), log(k[2]))) ## get N values between k[1] and k[2], assing to ki
|
||||
itheta=cbind(r=ri,k=ki, lam1=lambda[1],lam2=lambda[2], sigR=sigR)
|
||||
## assign ri, ki, and final biomass range to itheta
|
||||
M = apply(itheta,1,.schaefer) ## call Schaefer function with parameters in itheta
|
||||
i=1:N
|
||||
## prototype objective function
|
||||
get.ell=function(i) M[[i]]$ell
|
||||
ell = sapply(i, get.ell)
|
||||
get.J=function(i) M[[i]]$J # Ye
|
||||
J=sapply(i,get.J) # Ye
|
||||
return(list(r=ri,k=ki, ell=ell, J=J)) # Ye adding J=J
|
||||
})
|
||||
}
|
||||
|
||||
|
||||
getBiomass <- function(r, k, j)
|
||||
{
|
||||
BT <- NULL
|
||||
bt=vector()
|
||||
for (v in 1:length(r))
|
||||
{
|
||||
bt[1]=j[v]*k[v]*exp(rnorm(1,0, sigR)) ## set biomass in first year
|
||||
for(i in 1:nyr) ## for all years in the time series
|
||||
{
|
||||
xt=rnorm(1,0, sigR)
|
||||
bt[i+1]=(bt[i]+r[v]*bt[i]*(1-bt[i]/k[v])-ct[i])*exp(xt)
|
||||
## calculate biomass as function of previous year's biomass plus net production minus catch
|
||||
}
|
||||
BT=rbind(BT, t(t(bt)))
|
||||
}
|
||||
return(BT)
|
||||
}
|
||||
|
||||
## The End of Functions section
|
||||
|
||||
|
||||
cat("Step 4","\n")
|
||||
stockLoop <- StockList
|
||||
# randomly select stocks from randomly selected 5 area codes first
|
||||
if(TestRUN)
|
||||
{
|
||||
set.seed(999)
|
||||
AreaCodeList <- unique(cdat1$AREA_Code)
|
||||
sampledAC <- sample(AreaCodeList,size=5,replace=F)
|
||||
stockLoop <- cdat1[cdat1$AREA_Code %in% sampledAC,c("Stock_ID")]
|
||||
}
|
||||
|
||||
#setup counters
|
||||
counter1 <- 0
|
||||
counter2 <- 0
|
||||
|
||||
cat("Step 4","\n")
|
||||
## Loop through stocks
|
||||
for(stock in stockLoop)
|
||||
{
|
||||
t0<-Sys.time()
|
||||
##stock = "3845" # NB only for test single loop!
|
||||
## make graph file names:
|
||||
b <- with(cdat1,cdat1[Stock_ID == stock,c(1,3,5,12)]) # Stock_ID,AREA_Names,Country,"Species"
|
||||
bb <- do.call(paste,b)
|
||||
|
||||
yr <- cdat$yr[as.character(cdat$stock)==stock]
|
||||
ct <- as.numeric(cdat$ct[as.character(cdat$stock)==stock])/1000 ## assumes that catch is given in tonnes, transforms to '000 tonnes
|
||||
res <- unique(as.character(cdat$res[as.character(cdat$stock)==stock])) ## resilience from FishBase, if needed, enable in PARAMETER SECTION
|
||||
nyr <- length(yr) ## number of years in the time series
|
||||
|
||||
cat("\n","Stock",stock,"\n")
|
||||
flush.console()
|
||||
|
||||
## PARAMETER SECTION
|
||||
mvlen=3
|
||||
ma=function(x,n=mvlen){filter(x,rep(1/n,n),sides=1)}
|
||||
|
||||
## If resilience is to be used, delete ## in rows 1-4 below and set ## in row 5 below
|
||||
start_r <- if(res == "Very low"){c(0.015, 0.1)}else{
|
||||
if(res == "Low") {c(0.05,0.5)}else {
|
||||
if(res == "High") {c(0.6,1.5)}else {c(0.2,1)}
|
||||
}
|
||||
}
|
||||
## Medium, or default if no res is found
|
||||
##start_r <- c(0.5,1.5) ## disable this line if you use resilience
|
||||
start_k <- c(max(ct),50*max(ct)) ## default for upper k e.g. 100 * max catch
|
||||
## startbio <- c(0.8,1) ## assumed biomass range at start of time series, as fraction of k
|
||||
##startbio <- if(ct[1]/max(ct) < 0.5) {c(0.5,0.9)} else {c(0.3,0.6)} ## use for batch processing
|
||||
|
||||
## NB: Yimin's new idea on 20Jan14
|
||||
startbio<- if(mean(ct[1:5])/max(ct) < 0.3) {c(0.6,0.95)} else {
|
||||
if(mean(ct[1:5])/max(ct)>0.3&mean(ct[1:5])/max(ct)<0.6) {c(0.3,0.7)} else {
|
||||
c(0.2,0.6)}}
|
||||
|
||||
interyr <- yr[2] ## interim year within time series for which biomass estimate is available; set to yr[2] if no estimates are available
|
||||
interbio <- c(0, 1) ## biomass range for interim year, as fraction of k; set to 0 and 1 if not available
|
||||
## finalbio <- c(0.8, 0.9) ## biomass range after last catches, as fraction of k
|
||||
## finalbio <- if(ct[nyr]/max(ct) > 0.5) {c(0.3,0.7)} else {c(0.01,0.4)} ## use for batch processing
|
||||
|
||||
## Yimin's new stuff on 10Mar14
|
||||
#######> pre-classification
|
||||
|
||||
pre.clas=ct
|
||||
pre.clas[pre.clas==0]=0.1
|
||||
tx=ma(as.numeric(pre.clas),n=mvlen)
|
||||
Myr=which.max(tx)
|
||||
Maxc=pre.clas[which.max(tx)]
|
||||
|
||||
|
||||
if(Myr==1)startbio=c(0.05,0.6)else
|
||||
{
|
||||
if (ct[1]/Maxc>=0.5) startbio=c(0.4,0.85)
|
||||
else startbio=c(0.65,0.95)
|
||||
}
|
||||
|
||||
if (Myr==length(yr))finalbio=c(.4,.95) else # ie from fully to overexploited
|
||||
{
|
||||
if (tx[length(ct)]/Maxc>=0.5) finalbio=c(.4,.85)
|
||||
else finalbio=c(.05,.6)
|
||||
}
|
||||
|
||||
|
||||
# if (Myr==length(yr))finalbio=c(.5,.9)
|
||||
# #if (Myr<length(yr)){
|
||||
# # if ((tx[length(ct)]/Maxc)>=0.8) finalbio=c(.4,.8) else
|
||||
# # if (tx[length(ct)]/Maxc>0.5) finalbio=c(.3,.7) else finalbio=c(.05,.6)}
|
||||
# # below is the last used (20 Feb)
|
||||
# if (Myr<length(yr))
|
||||
# {
|
||||
# if (tx[length(ct)]/Maxc>0.5) finalbio=c(.2,.8)
|
||||
# else finalbio=c(.05,.6)
|
||||
# }
|
||||
|
||||
##############<
|
||||
n <- 30000 ## number of iterations, e.g. 100000
|
||||
sigR <- 0.0 ## process error; 0 if deterministic model; 0.05 reasonable value? 0.2 is too high
|
||||
|
||||
startbt <- seq(startbio[1], startbio[2], by = 0.05) ## apply range of start biomass in steps of 0.05
|
||||
parbound <- list(r = start_r, k = start_k, lambda = finalbio, sigR)
|
||||
|
||||
cat("Last year =",max(yr),", last catch =",1000*ct[nyr],"\n")
|
||||
cat("Resilience =",res,"\n")
|
||||
cat("Process error =", sigR,"\n")
|
||||
cat("Assumed initial biomass (B/k) =", startbio[1],"-", startbio[2], " k","\n")
|
||||
cat("Assumed intermediate biomass (B/k) in", interyr, " =", interbio[1],"-",interbio[2]," k","\n")
|
||||
cat("Assumed final biomass (B/k) =", parbound$lambda[1],"-",parbound$lambda[2]," k","\n")
|
||||
cat("Initial bounds for r =", parbound$r[1], "-", parbound$r[2],"\n")
|
||||
cat("Initial bounds for k =", format(1000*parbound$k[1], digits=3), "-", format(1000*parbound$k[2],digits=3),"\n")
|
||||
|
||||
flush.console()
|
||||
|
||||
## MAIN
|
||||
|
||||
R1 = sraMSY(parbound, n)
|
||||
|
||||
## Get statistics on r, k, MSY and determine new bounds for r and k
|
||||
r1 <- R1$r[R1$ell==1]
|
||||
k1 <- R1$k[R1$ell==1]
|
||||
j1 <- R1$J[R1$ell==1] # Ye
|
||||
msy1 <- r1*k1/4
|
||||
mean_msy1 <- exp(mean(log(msy1)))
|
||||
max_k1a <- min(k1[r1<1.1*parbound$r[1]]) ## smallest k1 near initial lower bound of r
|
||||
max_k1b <- max(k1[r1*k1/4<mean_msy1]) ## largest k1 that gives mean MSY
|
||||
max_k1 <- if(max_k1a < max_k1b) {max_k1a} else {max_k1b}
|
||||
|
||||
if(length(r1)<10)
|
||||
{
|
||||
cat("Too few (", length(r1), ") possible r-k combinations,
|
||||
check input parameters","\n")
|
||||
appendPar <- ifelse(counter1==0,F,T)
|
||||
colnamePar <- ifelse(counter1==0,T,F)
|
||||
|
||||
NoModellingSpe <- as.data.frame(cbind(stock,length(r1),b))
|
||||
names(NoModellingSpe) <- c("Stock","No_of_r1",names(b))
|
||||
write.table(NoModellingSpe,file=outfile2,
|
||||
append = appendPar, row.names = FALSE,
|
||||
col.names=colnamePar,sep=",")
|
||||
flush.console()
|
||||
counter1 <- counter1 + 1
|
||||
}
|
||||
|
||||
if(length(r1)>=10)
|
||||
{
|
||||
## set new upper bound of r to 1.2 max r1
|
||||
parbound$r[2] <- 1.2*max(r1)
|
||||
## set new lower bound for k to 0.9 min k1 and upper bound to max_k1
|
||||
parbound$k <- c(0.9 * min(k1), max_k1)
|
||||
|
||||
cat("First MSY =", format(1000*mean_msy1, digits=3),"\n")
|
||||
cat("First r =", format(exp(mean(log(r1))), digits=3),"\n")
|
||||
cat("New upper bound for r =", format(parbound$r[2],digits=2),"\n")
|
||||
cat("New range for k =", format(1000*parbound$k[1], digits=3), "-", format(1000*parbound$k[2],digits=3),"\n")
|
||||
|
||||
## Repeat analysis with new r-k bounds
|
||||
R1 = sraMSY(parbound, n)
|
||||
|
||||
## Get statistics on r, k and msy
|
||||
r = R1$r[R1$ell==1]
|
||||
k = R1$k[R1$ell==1]
|
||||
j = R1$J[R1$ell==1] # Ye
|
||||
msy = r * k / 4
|
||||
mean_ln_msy = mean(log(msy))
|
||||
|
||||
##############################################################
|
||||
##> Ye
|
||||
# BT=0
|
||||
|
||||
##
|
||||
R2<-getBiomass(r, k, j)
|
||||
|
||||
#R2<-R2[-1,]
|
||||
runs<-rep(1:length(r), each=nyr+1)
|
||||
years=rep(yr[1]:(yr[length(yr)]+1),length=length(r)*(length(yr)+1))
|
||||
|
||||
runs=t(runs)
|
||||
years=t(years)
|
||||
stock_id=rep(stock,length(runs))
|
||||
R3<-cbind(as.numeric(runs), as.numeric(years), stock_id, as.numeric(R2) )
|
||||
|
||||
## changed this, as otherwise biomass is the level of the factor below
|
||||
R4<-data.frame(R3, stringsAsFactors=FALSE)
|
||||
names(R4)<-c("Run", "Year", "Stock","Biomass")
|
||||
|
||||
Bmsy_x<-k*0.5
|
||||
Run<-c(1:length(r))
|
||||
BMSY<-cbind(Run, Bmsy_x)
|
||||
R5<-merge(R4, BMSY, by="Run", all.x=T, all.y=F)
|
||||
R5$B_Bmsy<-as.numeric(paste(R5$Biomass))/R5$Bmsy_x
|
||||
|
||||
### B/Bmsy calculated for all feasible combinations of r,K,B0
|
||||
R6<-aggregate(log(B_Bmsy)~as.numeric(Year)+Stock, data=R5,
|
||||
FUN=function(z){c(mean=mean(z),sd=sd(z),upr=exp(quantile(z, p=0.975)),
|
||||
lwr=exp(quantile(z, p=0.025)), lwrQ=exp(quantile(z, p=0.25)),
|
||||
uprQ=exp(quantile(z, p=0.75)))}) # from directly calculated from R5 becasue B_Bmsy has a lognormal dist
|
||||
|
||||
R6<-data.frame(cbind(R6[,1:2],R6[,3][,1],R6[,3][,2],R6[,3][,3],R6[,3][,4],R6[,3][,5], R6[,3][,6]))
|
||||
names(R6)<-c("Year", "Stock", "BoverBmsy", "BoverBmsySD","BoverBmsyUpper","BoverBmsyLower","BoverBmsylwrQ","BoverBmsyuprQ")
|
||||
##remove last entry as it is 1 greater than number of years
|
||||
## removed final year here for ease of dataframe output below
|
||||
R6<-R6[-length(R6),]
|
||||
## geometric mean
|
||||
GM_B_Bmsy<-exp(R6$BoverBmsy)
|
||||
GM_B_BmsySD=R6$BoverBmsySD #add
|
||||
## arithmetic mean
|
||||
M_B_Bmsy<-exp(R6$BoverBmsy+R6$BoverBmsySD^2/2)
|
||||
|
||||
### r,k, and MSY
|
||||
|
||||
#del GM_B_Bmsy=c(rep(0,(min(yr)-1940)),GM_B_Bmsy)
|
||||
#del GM_B_BmsySD=c(rep(0,(min(yr)-1940)),GM_B_BmsySD) ######
|
||||
#del M_B_Bmsy=c(rep(0,(min(yr)-1940)),M_B_Bmsy)
|
||||
#del yr1=seq(1940,max(yr))
|
||||
|
||||
yr1=yr #add
|
||||
|
||||
stockInfo <- with(cdat1,cdat1[Stock_ID==stock,1:12])
|
||||
temp=c(startbio[1],startbio[2],finalbio[1],finalbio[2],res,
|
||||
mean(log(r)),sd(log(r)),mean(log(k)),sd(log(k)),mean(log(msy)),
|
||||
sd(log(msy)),sigR,min(yr),max(yr),max(ct),length(r),GM_B_Bmsy,GM_B_BmsySD,M_B_Bmsy)
|
||||
|
||||
#add, adding "GM_B_BmsySD" in the line above
|
||||
|
||||
output=as.data.frame(matrix(temp,nrow=1))
|
||||
output <- cbind(stockInfo,output)
|
||||
names(output) <- c(names(cdat1)[1:12],"startbio[1]","startbio[2]","finalbio[1]","finalbio[2]",
|
||||
"res","mean(log(r))","sd(log(r))","mean(log(k))","sd(log(k))",
|
||||
"mean(log(msy))","sd(log(msy))","sigR","min(yr)","max(yr)","max(ct)",
|
||||
"length(r)",paste("GM_B_msy",yr1,sep="_"),paste("GM_B_msySD",yr1,sep="_"),paste("M_B_Bmsy",yr1,sep="_"))
|
||||
|
||||
#add, adding "paste("GM_B_msySD",yr1,sep="_")"in the line above
|
||||
|
||||
######< Ye
|
||||
########################################################
|
||||
|
||||
## plot MSY over catch data
|
||||
pdf(paste(bb,"graph.pdf",sep="_"))
|
||||
|
||||
par(mfcol=c(2,3))
|
||||
plot(yr, ct, type="l", ylim = c(0, max(ct)), xlab = "Year",
|
||||
ylab = "Catch (1000 t)",main = paste("StockID",stock,sep=":"))
|
||||
abline(h=exp(mean(log(msy))),col="red", lwd=2)
|
||||
abline(h=exp(mean_ln_msy - 2 * sd(log(msy))),col="red")
|
||||
abline(h=exp(mean_ln_msy + 2 * sd(log(msy))),col="red")
|
||||
|
||||
hist(r, freq=F, xlim=c(0, 1.2 * max(r)), main = "")
|
||||
abline(v=exp(mean(log(r))),col="red",lwd=2)
|
||||
abline(v=exp(mean(log(r))-2*sd(log(r))),col="red")
|
||||
abline(v=exp(mean(log(r))+2*sd(log(r))),col="red")
|
||||
|
||||
plot(r1, k1, xlim = start_r, ylim = start_k, xlab="r", ylab="k (1000t)")
|
||||
|
||||
hist(k, freq=F, xlim=c(0, 1.2 * max(k)), xlab="k (1000t)", main = "")
|
||||
abline(v=exp(mean(log(k))),col="red", lwd=2)
|
||||
abline(v=exp(mean(log(k))-2*sd(log(k))),col="red")
|
||||
abline(v=exp(mean(log(k))+2*sd(log(k))),col="red")
|
||||
|
||||
|
||||
plot(log(r), log(k),xlab="ln(r)",ylab="ln(k)")
|
||||
abline(v=mean(log(r)))
|
||||
abline(h=mean(log(k)))
|
||||
abline(mean(log(msy))+log(4),-1, col="red",lwd=2)
|
||||
abline(mean(log(msy))-2*sd(log(msy))+log(4),-1, col="red")
|
||||
abline(mean(log(msy))+2*sd(log(msy))+log(4),-1, col="red")
|
||||
|
||||
hist(msy, freq=F, xlim=c(0, 1.2 * max(msy)), xlab="MSY (1000t)",main = "")
|
||||
abline(v=exp(mean(log(msy))),col="red", lwd=2)
|
||||
abline(v=exp(mean_ln_msy - 2 * sd(log(msy))),col="red")
|
||||
abline(v=exp(mean_ln_msy + 2 * sd(log(msy))),col="red")
|
||||
|
||||
graphics.off()
|
||||
|
||||
cat("Possible combinations = ", length(r),"\n")
|
||||
cat("geom. mean r =", format(exp(mean(log(r))),digits=3), "\n")
|
||||
cat("r +/- 2 SD =", format(exp(mean(log(r))-2*sd(log(r))),digits=3),"-",format(exp(mean(log(r))+2*sd(log(r))),digits=3), "\n")
|
||||
cat("geom. mean k =", format(1000*exp(mean(log(k))),digits=3), "\n")
|
||||
cat("k +/- 2 SD =", format(1000*exp(mean(log(k))-2*sd(log(k))),digits=3),"-",format(1000*exp(mean(log(k))+2*sd(log(k))),digits=3), "\n")
|
||||
cat("geom. mean MSY =", format(1000*exp(mean(log(msy))),digits=3),"\n")
|
||||
cat("MSY +/- 2 SD =", format(1000*exp(mean_ln_msy - 2 * sd(log(msy))),digits=3), "-", format(1000*exp(mean_ln_msy + 2 * sd(log(msy))),digits=3), "\n")
|
||||
|
||||
## Write results into outfile, in append mode (no header in file, existing files will be continued)
|
||||
## output = data.frame(stock, sigR, startbio[1], startbio[2], interbio[1], interbio[2], finalbio[1], finalbio[2], min(yr), max(yr), res, max(ct), ct[1], ct[nyr], length(r), exp(mean(log(r))), sd(log(r)), min(r), quantile(r,0.05), quantile(r,0.25), median(r), quantile(r,0.75), quantile(r,0.95), max(r), exp(mean(log(k))), sd(log(k)), min(k), quantile(k, 0.05), quantile(k, 0.25), median(k), quantile(k, 0.75), quantile(k, 0.95), max(k), exp(mean(log(msy))), sd(log(msy)), min(msy), quantile(msy, 0.05), quantile(msy, 0.25), median(msy), quantile(msy, 0.75), quantile(msy, 0.95), max(msy))
|
||||
|
||||
#write.table(output, file = outfile, append = TRUE, sep = ";", dec = ".", row.names = FALSE, col.names = FALSE)
|
||||
appendPar <- ifelse(counter2==0,F,T)
|
||||
colnamePar <- ifelse(counter2==0,T,F)
|
||||
write.table(output, file = outfile, append = appendPar, sep = ",", dec = ".",
|
||||
row.names = FALSE, col.names = colnamePar)
|
||||
|
||||
counter2 <- counter2 + 1
|
||||
|
||||
}
|
||||
cat("Elapsed: ",Sys.time()-t0," \n")
|
||||
} ## End of stock loop, get next stock or exit
|
||||
|
|
@ -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>
|
||||
Binary file not shown.
|
|
@ -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.
|
|
@ -0,0 +1,5 @@
|
|||
#!/bin/sh
|
||||
# LWR
|
||||
cd $1
|
||||
|
||||
java -Xmx1024M -classpath ./:./c3p0-0.9.1.2.jar:./common-configuration-scanner-1.0.1-SNAPSHOT.jar:./common-encryption-1.0.1-3.5.0.jar:./common-gcore-resources-1.2.0-3.5.0.jar:./common-gcore-stubs-1.2.0-3.5.0.jar:./common-scope-1.2.1-SNAPSHOT.jar:./common-scope-maps-1.0.2-3.5.0.jar:./commons-collections-3.1.jar:./commons-io-1.2.jar:./discovery-client-1.0.1-3.5.0.jar:./dom4j-1.6.1.jar:./ecological-engine-1.8.1-SNAPSHOT.jar:./EcologicalEngineExecutor-1.6.4-SNAPSHOT.jar:./hibernate3.jar:./ic-client-1.0.1-3.5.0.jar:./jaxen-1.1.2.jar:./jta-1.1.jar:./log4j-1.2.16.jar:./mongo-java-driver-2.12.4.jar:./postgresql-8.4-702.jdbc4.jar:./slf4j-api-1.6.0.jar:./slf4j-log4j12-1.6.0.jar:./storage-manager-core-2.1.3-3.6.0.jar:./storage-manager-wrapper-2.1.0-3.5.0.jar:./xalan-2.6.0.jar:./xpp3_min-1.1.4c.jar:./xstream-1.3.1.jar:./YASMEEN-matcher-1.2.0.1.jar:./YASMEEN-parser-1.2.0.jar org.gcube.dataanalysis.executor.nodes.algorithms.LWR $2 execution.output
|
||||
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.
Loading…
Reference in New Issue