10 Recruitment
Collins 1892.
Growth, mortality, and recruitment are the three fundamental processes that regulate fish populations. Previous chapters have addressed the first two; here we consider methods for studying recruitment. Recruitment is usually considered to be the annual process of young fish entering the part of the population that is vulnerable to fishing (i.e., the fishable stock). Thus the timing of recruitment depends on fishery selectivity pattern(s), which tend to be gradual rather than knife-edged. One practical approach for defining age at recruitment is to choose the first age at which fishery catches are non-negligible. For simplicity in this chapter, we will simulate recruitment occurring at age one so there is a one-year lag from the time of spawning to recruitment.
What happens between spawning and recruitment is of interest to the fish ecologist and the fishery manager. Typical factors affecting recruitment include environmental conditions, prey availability, and predation. It may seem obvious to assume that spawning stock size affects recruitment, but it is not always easy to detect such a relationship (Hilborn and Walters 1992). The relationship may not be evident if there is little contrast in spawning stock size, high variability in recruitment, too short a time series, or low precision for estimates of spawning stock and recruitment. Characterizing the strength of the relationship can provide valuable guidance to the fishery manager. For example, management will need to be more restrictive in cases where recruitment drops consistently as spawning stock declines. In other cases, a more aggressive harvest policy may be possible because recruitment appears to be relatively stable over a wide range of spawning stock sizes. The risk in any case is of fishing the stock down to a point where recruitment declines sharply. This is termed recruitment overfishing, and it can be very difficult to rebuild a population after that occurs.
The process of recruitment starts with eggs, but egg production is difficult to estimate directly so it is common to use spawning stock biomass (\(B_S\)) as a proxy. One of the most commonly used curves relating spawning stock size and recruitment is the Beverton-Holt equation: \(R=1/(\alpha +\beta/B_S)\). This curve approaches zero as \(B_S\) approaches zero and approaches 1/\(\alpha\) as \(B_S\) approaches infinity. It is often a good descriptor of observed stock-recruitment data in that recruitment is low when spawning stock size is close to zero but appears to be without trend at higher spawning stock sizes.
The best way to explore the behavior of the Beverton-Holt curve is to try different values for the two parameters:
SpBio <- seq(from=10, to=200, by=5)
a <- 0.001
b <- 0.1
R <- 1/(a+b/SpBio)
plot(SpBio, R)In this example code, spawning stock biomass is generated with the seq() (sequence) function over a wide but arbitrary range. The Beverton-Holt parameters are also arbitrary but are chosen to allow for gradually increasing recruitment toward an asymptote of 1000. Experiment with different parameter values to find a set that would maintain an asymptote of 1000 but would produce a steeper curve, with about 90% of asymptotic recruitment at a spawning stock level of 100.
10.1 Fitting a stock-recruitment curve
We begin by fitting a Beverton-Holt model to simulated stock-recruitment data. Section 10.2 describes an approach for comparing full and reduced models, to explore the conditions under which an underlying relationship can be detected.
Our simulation uses a full age-structured model to generate the stock-recruitment data. First, we choose the structure of the population matrix (number of ages and years). Next, we choose arbitrary annual fishing mortality rates. By starting the population at a high level and fishing intensively, we ensure that there is sufficient contrast (Hilborn and Walters 1992) in spawning stock size to be able to detect the underlying relationship. Age-specific fishing mortality rates are obtained by defining the fishery selectivity pattern. We use a logistic curve and arbitrarily choose parameters so that 50% selectivity occurs at age two. (It is useful to plot the selectivity curve (code commented out) while varying the two selectivity parameters.) Obtaining year- and age- specific fishing mortality rates as a product of a year-specific measure of fishing intensity and an age-specific vulnerability to fishing is sometimes referred to as a “separable” model.
rm(list=ls()) # Clear Environment
A <- 10 # Arbitrary. If changed, redo maturity and meanwt vectors
Y <- 12 # Arbitrary. If changed, redo AnnualF vector
AnnualF <- c(0.59, 0.41, 0.74, 0.91, 0.86, 0.74, 1.07, 0.9, 0.87, 1.1, 0.93)
# Arbitrary, increasing trend to create contrast in SpBio. Redo if Y changes
AnnualM <- 0.2
# F by age, using separable model
Fishery_k <- 1 # Slope for logistic function
Fishery_a_50 <- 2 # Age at 0.5 selectivity
SelF <- array(data=NA, dim=A) # Fishery selectivity pattern
F_a <- array(data=NA, dim=c(A,(Y-1))) # Matrix for F by age and year
Z_a <- array(data=NA, dim=c(A,(Y-1))) # Matrix for Z by age and year
S_a <- array(data=NA, dim=c(A,(Y-1))) # Matrix for survival rate by age and year
for (a in 1:A){
SelF[a] <- 1/(1+exp(-Fishery_k*(a-Fishery_a_50)))
} #a
#plot(SelF)
for (y in 1:(Y-1)) {
for (a in 1:A){
F_a[a,y] <- AnnualF[y] * SelF[a]
Z_a[a,y] <- F_a[a,y]+AnnualM
S_a[a,y] <- exp(-Z_a[a,y])
} #a
} #y
N <- array(data=NA, dim=c(A, Y))
BH_alpha <- 1/1E5 # Parameter defining Beverton-Holt asymptotic recruitment
BH_beta <- 0.1 # Arbitrary
SD_rec <- 0.25 # Low level of lognormal error in recruitment
# Set up year-1 vector based on asymptotic recruitment and year-1 survival rates
Mat <- c(0, 0, 0.2, 0.4, 0.8, 1, 1, 1, 1, 1) # Arbitrary maturity schedule. Redo if max age (A) changes
Wgt <- c(0.01, 0.05, 0.10, 0.20, 0.40, 0.62, 0.80, 1.01, 1.30, 1.56) # Mean wt. Redo if A changes
SpBio <- array(data=NA, dim=(Y-1)) # Spawning biomass vector
# Year-1 N (arbitrarily) from asymptotic recruitment, year-1 survival.
N[1,1] <- 1/BH_alpha*rlnorm(n=1, 0, SD_rec)
for (a in 2:A){
N[a,1] <- rbinom(n=1, trunc(N[(a-1),1]), S_a[(a-1),1])
}#a
for (y in 2:Y){
SpBio[y-1] <- sum(N[,y-1]*Wgt*Mat) # Spawning stock biomass
N[1,y] <- 1/(BH_alpha+BH_beta/SpBio[y-1])*rlnorm(n=1, 0, SD_rec)
for (a in 2:A){
N[a,y] <- rbinom(1, trunc(N[(a-1), (y-1)]), S_a[(a-1),(y-1)])
} #a
} #y
#plot(SpBio,N[1,2:Y], ylab="Recruits", xlab="Spawning stock biomass")
# Generate survey data on recruitment and spawning stock
Survey_q <- 1E-3 # Arbitrary catchability coefficient for survey
SD_survey <- 0.2 # Arbitrary ln-scale SD
Exp_ln_B <- log(Survey_q*SpBio)
Survey_B <- rlnorm(n=Y, meanlog=Exp_ln_B, sdlog=SD_survey)
#plot(SpBio, Survey_B[1:Y-1])
Exp_ln_R <- log(Survey_q*N[1,2:Y])
Survey_R <- c(NA,rlnorm(n=(Y-1), meanlog=Exp_ln_R, sdlog=SD_survey))
#plot(N[1,], Survey_R)The next section of code generates the population matrix. The Beverton-Holt parameters are arbitrarily chosen and can be varied to explore different shapes for the curve. Next, we generate the population vector for year one, using approximate equilibrium values. Expected recruitment is set equal to the curve’s asymptote, with lognormal error term used to add random variation; for example, due to environmental conditions. The lognormal distribution is well suited to recruitment data because of the potential for occasional extreme values (Hilborn and Walters 1992). The chosen standard deviation (0.25) is on the lower end of variation in recruitment (Hightower and Grossman 1985). The rest of the year-1 vector is obtained from age-1 recruitment using year-1 mortality rates (\(N_{a,1}=S_{a-1,1}*N_{a-1,1}\)). A binomial function introduces stochastic variation in survival. Population size at the previous age is truncated because the binomial function requires a whole number for the size of the experiment.
The population vector for the remaining years can be filled in once the year-1 vector is determined. For example, recruitment in year 2 depends on spawning stock biomass in year 1. The age-specific vectors for maturity schedule and weight-at-age are arbitrary. Numbers of fish age 2 and older are obtained as survivors from the prior year’s vector; i.e., \(N_{a,y}=S_{a-1,y-1}*N_{a-1,y-1}\), with random variation in survival from a binomial distribution. The final step is to obtain survey observations for recruitment and spawning stock biomass. The random survey values are drawn from a lognormal distribution, using an arbitrary catchability coefficient and ln-scale standard deviation for error. The survey values are referred to as indices but could be estimates of the absolute magnitude of recruitment and spawning stock size. There is not a recruitment survey value for year 1 because the simulated recruitment was not based on spawning stock size.
The JAGS code uses uninformative prior distributions for the two Beverton-Holt parameters and the lognormal error term. Expected recruitment is predicted on ln-scale in order to match the assumed lognormal error. The recruitment index is predicted (R_hat) for use in plotting (see below).
# Load necessary library
library(rjags)
library(R2jags)
sink("StockRecruit.txt")
cat("
model {
# Model parameters: Beverton-Holt alpha and beta, SD for fitted curve
# Priors
BH_alpha.est ~ dunif(0, 10)
BH_beta.est ~ dunif(0, 10)
SD_rec.est ~ dunif(0, 10)
tau <- pow(SD_rec.est, -2)
# Likelihood
for(y in 2:Y) {
lnR_hat[y] <- log(1/(BH_alpha.est+BH_beta.est/Survey_B[y-1]))
Survey_R[y] ~ dlnorm(lnR_hat[y],tau)
R_hat[y] <- exp(lnR_hat[y])
} #y
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("Y", "Survey_R", "Survey_B")
# Initial values
y <- Survey_R[2:Y]
x <- Survey_B[1:(Y-1)]
nls_init <- nls(y ~ 1/(BH_a_start+BH_b_start/x),
start = list(BH_a_start = 0.001, BH_b_start = 0.001),
algorithm="port",
lower=c(1E-6, 0), upper=c(1000,1000))
#summary(nls_init)
nls_save <- coef(nls_init)
#y_hat <- 1/(nls_save[1]+nls_save[2]/x)
#plot(x,y)
#points(x,y_hat, col="red", type="l")
jags.inits <- function(){ list(BH_alpha.est=nls_save[1],
BH_beta.est=nls_save[2],
SD_rec.est=runif(n=1, min=0, max=0.75))}
model.file <- 'StockRecruit.txt'
# Parameters monitored
jags.params <- c("BH_alpha.est", "BH_beta.est", "SD_rec.est",
"R_hat")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 10000,
model.file)
print(jagsfit)
plot(jagsfit)
BH.DIC <- jagsfit$BUGSoutput$DIC # Used in next section
min.y <- min(Survey_R[2:Y],jagsfit$BUGSoutput$mean$R_hat)
max.y <- max(Survey_R[2:Y],jagsfit$BUGSoutput$mean$R_hat)
df.xy <- data.frame(Survey_B[1:(Y-1)], Survey_R[2:Y],jagsfit$BUGSoutput$mean$R_hat)
df.xy <- df.xy[order(df.xy[,1]),] # Sort by survey index for plotting
plot(df.xy[,1],df.xy[,2], ylab="Recruitment index",
xlab="Spawning Stock Index", ylim=c(min.y, max.y))
points(df.xy[,1], df.xy[,3], type="l", col="red")Randomly chosen starting values for the two Beverton-Holt parameters sometimes caused JAGS to crash, so nls(), the R function for nonlinear least-squares, was used to generate good initial values. When unconstrained, nls() occasionally produced a negative estimate of the first parameter (1/asymptotic recruitment), so the “port” algorithm was used to allow bounds to be specified. The nls() code illustrates the compactness of using R functions, compared to a Bayesian analysis, but also the “black box” nature of the analysis.
The JAGS results tend to be reliable for the default settings. The plotted estimates of the recruitment indices capture the decrease due to declining spawning stock size. Estimates of the Beverton-Holt alpha parameter are typically close to 0.01 (BH_alpha/Survey_q), because we fit the curve to a recruitment index (true recruitment * Survey_q) rather than to absolute levels. Estimates of the Beverton-Holt beta parameter tend to be close to the true value. Estimates of the standard deviation for variation about the curve tend to be higher than the true value of 0.25, because the points vary due not only to true recruitment variation but also variation due to survey sampling.
As always, it is very helpful to view the points and fitted curve. We use a data frame to sort the spawning stock, recruitment, and predicted recruitment indicies (link), in order to have a smooth fitted curve (i.e., continuously increasing spawning stock values). Try multiple runs to see how variation in recruitment affects the pattern of points and fitted curve.
10.2 Does spawning stock size affect recruitment?
Section 10.1 illustrated the process of model fitting but did not answer the question of whether there is a detectable relationship between the spawning stock and recruitment indices. To address that question, we return to the approach shown in Section 4.3.6; i.e., comparing DIC scores for candidate models. Here the two candidates are the full Beverton-Holt model and a reduced model with constant recruitment.
We append JAGS code for the constant recruitment case at the end of the previous code so that both models are applied to the same spawning stock and recruitment data. Convergence is extremely rapid for this simpler case. The difference in DIC scores between the full (Beverton-Holt) and reduced (constant recruitment) models is printed to the Console, but it is also useful to compare pD, deviance, and DIC scores.
# Constant recruitment fit
# JAGS code ##################################
sink("StockRecruit.txt")
cat("
model {
# Model parameters: Beverton-Holt alpha, SD for fitted curve
# Priors
BH_alpha.est ~ dunif(0, 10)
SD_rec.est ~ dunif(0, 10)
tau <- pow(SD_rec.est, -2)
# Likelihood
for(y in 2:Y) {
lnR_hat[y] <- log(1/BH_alpha.est)
Survey_R[y] ~ dlnorm(lnR_hat[y],tau)
R_hat[y] <- exp(lnR_hat[y])
} #y
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("Y", "Survey_R")
jags.inits <- function(){ list(BH_alpha.est=nls_save[1],
SD_rec.est=runif(n=1, min=0, max=0.75))}
model.file <- 'StockRecruit.txt'
# Parameters monitored
jags.params <- c("BH_alpha.est", "SD_rec.est", "R_hat")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 10000,
model.file)
print(jagsfit)
plot(jagsfit)
ConstR.DIC <- jagsfit$BUGSoutput$DIC
ConstR.DIC - BH.DICWith the default settings, the DIC score is typically higher (indicating a poorer fit) for the constant recruitment model. Next, try varying some of the simulation settings. The pattern for annual fishing mortality could be altered to produce less of a decline over time in spawning stock (less contrast). Higher levels of recruitment variation are not uncommon [e.g., 0.5-0.75; Hightower and Grossman (1985)], and will make the relationship harder to detect. For a lognormal variance of 0.5 or 0.75, the DIC score may actually be lower for the constant recruitment model because of a smaller penalty for model complexity (smaller estimate of pD). Other potential settings to vary include the length of the time series and the precision of the survey data. Be sure to look back at the scatter plot and fitted Beverton-Holt curve, to get a sense of how the constant recruitment model would compare.
10.3 Exercises
Revise the code for Section 10.1 to fit a Ricker model (\(R=\alpha*B_S*exp(-\beta*B_S)\)). Experiment with the Ricker parameters to produce maximum recruitment at an intermediate observed spawning stock size.
The expected value for \(\alpha\) (0.01 for our default simulation settings) is the true underlying value divided by the survey catchability, because we are fitting the model using survey indices rather than estimates of absolute abundance. Compare estimates of \(\alpha\) for the Beverton-Holt and constant recruitment models, using the approach from Section 10.2. Explain any pattern that you detect.