R/StockRecCheck.R

In szuwalski/GeMS: Generic Management Strategy Evaluation

#' Stock-Recruitment Check
#' 
#' Pre-MSE, used to determine stock-recruitment relationship (Average, Beverton-Holt or Ricker), likely for use somewhere eventually. Currently uses \code{SRstarts} function that has not previously been specified. Probably related to \code{\link{calcREF}}
#' 
#' @param rec Vector of recruitment values
#' @param spbio Vector of spawning biomass values (must be of same length as \code{rec})
#' @param plotSR Logical; Plot stock-recruitment relationship
#' 
#' @return Recruitment function to be used 
#'
StockRecCheck<-function(rec,spbio,plotSR=F)
{
Stx		<-srStarts(rec~spbio,type="BevertonHolt",param=2)
x		<-c(Stx$a,Stx$Rp,1)
outs		<-optim(x,BevHolt,spbio=spbio)
counter	<-1
while(outs$convergence!=0 & counter<20)
{
 x		<-outs$par
 outs		<-optim(x,BevHolt,spbio=spbio)
 counter	<-counter+1
}

Stx1<-srStarts(rec~spbio,type="Ricker",param=2)
x1		<-c(Stx1$a,Stx1$b,1)
outs1		<-optim(x1,Ricker)
x1		<-outs1$par
outs1		<-optim(x1,Ricker)
x1		<-outs1$par
outs1		<-optim(x1,Ricker)

x2		<-c(mean(rec),sd(rec))
outs2		<-optim(x2,StrLine)

AvgAIC	<-outs2$value
RickerAIC	<-outs1$value
BHaic		<-outs$value

AvgPar	<-outs2$par
RickerPars	<-outs1$par
BHpars	<-outs$par

#===check parameters
AICR		<-6+2*RickerAIC
AICB		<-6+2*BHaic
AICA		<-4+2*AvgAIC

#==report decision about stock recruit relationship
SRchoice	<-"Average"
returnPar	<-AvgPar

if(AICR<AICA-2 & AICR<AICB)
{
SRchoice	<-"Rick"
returnPar	<-RickerPars
}

if(AICB<AICA-2 & AICB<AICR)
{
SRchoice	<-"BevH"
returnPar	<-BHpars
}

if(plotSR==T)
 {
  predSpbio<-seq(1,max(spbio),max(spbio)/40)
  alpha<-abs(BHpars[1])
  beta<-abs(BHpars[2])
  plotPred<-alpha*predSpbio/(1+(predSpbio/beta))
  plot(rec~spbio,ylim=c(0,max(rec)),xlim=c(0,max(spbio)))
  lines(plotPred~predSpbio)

  alpha<-RickerPars[1]
  beta<-RickerPars[2]
  plotPred<-predSpbio*exp(alpha-beta*predSpbio)
  lines(plotPred~predSpbio,lty=2,col=2)


  abline(h=AvgPar[1],lty=4,col=4)

  legend("topright",bty='n',col=c(1,2,4),lty=c(1,2,4),legend=c(paste("BH ",round(AICB,2)),
   paste("Rick ", round(AICR,2)),paste("Avg ",round(AICA,2))))
 }
list(SRchoice,returnPar)
}