In rBatt/fishyWarnings: How should management respond to warning signals of haddock collapse?

The examples in some of the fishyWarnings examples are fairly thorough simulations + analyses geared toward the purpose of the package, which of course is to see if early warnings can help avoid a haddock collapse.

Rebuilding Haddock

First, simulate the time series. Haddock starts at a low biomass (the 1993 value of 13.8 kt). Then we simulate forward 50 years under 6 different harvest rates.

library(fishyWarnings)
nY <- 50
# 13.8 kt is from pg 2922, second column, penultimate paragraph; 
# 13.8 is the starting value in 1993
makeB <- function(){c(13.8, rep(NA, nY-1))}
Year <- seq(1993, length.out=nY)
qE <- c(0, 0.06, 0.12, 0.18, 0.21, 0.24)
Bvec <- replicate(length(qE), makeB())
for(j in 1:length(qE)){
    for(i in 2:nrow(Bvec)){
        Bvec[i,j] <- Bstep(Bvec[i-1,j], qE[j], sdU=0.1)
    }
}
ltys <- c("solid", "dotted", "dashed", "longdash", "twodash", "dotdash")
ylim <- range(Bvec, na.rm=TRUE)

With the simulation complete, we can plot the results, recreating Figure 5A:

plot(Year, Bvec[,1], type='l', 
    lty=ltys[1], ylim=ylim, 
    ylab="Biomass (kt)", xlab="time", main="STH model"
)
for(j in 2:ncol(Bvec)){
    lines(Year, Bvec[,j], lty=ltys[j])
}
legend("topleft", paste("F", qE, sep="="), lty=ltys)

As you can see, depending on the harvest rate, the stock may or may not rebuild.

Collapse Haddock, Find Warnings

In this case we are going to start with a high biomass of haddock, and increase harvest rate as time passes.

In each year, we will calculate an early warning signal. The early warning signals used here are 1) standard deviation, 2) first order autocorrelation (AR(1)), and 3) redShift, which is a measure of spectral slope.

According to Spencer and Collie (1997), the critical point for collapse (going from high stable point to low) is at $qE = 0.36$. As the time series evolves and harvest intensifies and approach 0.36, the early warning indicators should beging to "fire off". AR(1) and SD should be elevated, and the redShift signal should decline.

Note: I am still playing with the redShift indicator, but the idea is that the variance spectrum should become stronger and should shift to relatively more power at lower frequencies. The statistic reported as redShift in these figures is the slope of log(power) vs 1/log(freq) (where freq is stats::spectrum()$freq, which is probably the inverse of what you would expect). Thus, a negative value for the spectrum slope means that more power is at lower frequencies, and decreases in this slope should mean a relative increase in low-frequency variability. So far, this is not the pattern that I have observed, so something may be wrong.

First, set up some options regarding simulation length and harvest value gradient.

# ==================
# = Set up Options =
# ==================
nYear <- 5E2
Year <- 1:nYear
tripVal <- c(0.36, 0.21)
nBurn <- 100
sdU <- 0.05

# set up qE to increase harvest until collapse
qE_range <- c(tripVal[1]-0.2, tripVal[1]+0.05)

# goal here is to have harvest rate precisely hit certain 
# values of qE during an integer year,
# while retaining linear change throughout
qEPoints <- c(qE_range[1], tripVal[2], tripVal[1], qE_range[2])
qEPoints <- qEPoints[qEPoints>=min(qE_range) & qEPoints<=max(qE_range)]
seq_arg0 <- c(as.list(range(qEPoints)), list(length.out=nYear))
yrProbs <- (qEPoints-min(qEPoints))/diff(range(qEPoints))
yrPoints <- quantile(Year, yrProbs)
qE <- approx(x=yrPoints, y=qEPoints, xout=Year)$y

Method 1: Approach equilibrium per qE

For each value of qE (harvest rate), simulate {r nBurn} years of data, calculating the warning indicators for that entire time series (it's important to note that qE is constant here). The idea is that the system is more-or-less at equilibrium, which contrasts with what we'll do next (Method 2).

# ===========================================
# = Simulate Biomass and Calculate Warnings =
# ===========================================
# For each year, simulate nBurn years
# B in the top-level year is the mean of the nBurn sims
# Early warnings calculated from the nBurn, so no rolling window
warnMat <- matrix(NA, ncol=3, nrow=nYear, dimnames=list(NULL, c("sd","ac1","redShift")))
stats <- colnames(warnMat)
Bmu <- numeric(nYear)
B0 <- 400
rsl <- structure(vector('list', nYear), class="rslist", .Names=Year) # for the plot.rslist method I wrote
for(i in 1:nYear){
    B0 <- Burn(B0, n=100, qE=qE[i], sdU=0, accumulate=FALSE)
    B_t <- Burn(B0, n=nBurn, qE=qE[i], sdU=sdU, accumulate=TRUE)
    for(s in 1:ncol(warnMat)){
        if(stats[s]=="redShift"){
            rsl[[i]] <- ews(B_t, stat=stats[s])
            warnMat[i,s] <- rsl[[i]]
        }else{
            warnMat[i,s] <- ews(B_t, stat=stats[s])
        }
    }
    Bmu[i] <- mean(B_t)
}

# =============================
# = Plot Biomass and Warnings =
# =============================
par(mfrow=c(2,2), mar=c(1.75,1.75,0.5,1.75), ps=8, cex=1, mgp=c(0.75,0.15,0), tcl=-0.15)
plot(Year, Bmu, type='l')
abline(v=yrPoints, lty='dashed')
text_y <- min(Bmu)+yrProbs*diff(range(Bmu))
text(yrPoints, y=text_y, label=paste0("qE=", qEPoints), pos=c(4,4,2,2))

plot(Year, warnMat[,"sd"], col="forestgreen", type="l", ylab="Standard Deviation")
abline(v=yrPoints, lty='dashed')
plot(Year, warnMat[,"ac1"], col="blue", type='l', ylab="AR(1)")
abline(v=yrPoints, lty='dashed')
plot(rsl, xaxs='r')
par(new=TRUE)
plot(Year, warnMat[,"redShift"], 
    col=adjustcolor("white",0.5), 
    lwd=3, type='l', ylab="", xaxt='n', yaxt='n'
)
lines(Year, warnMat[,"redShift"], col=adjustcolor("black",0.5), lwd=0.5)
axis(side=4)
abline(v=yrPoints, lty='dashed')

Method 2: Let biomass evolve as qE changes slowly

In this approach biomass is not iterated forward nBurn time steps per qE value. Instead, biomass and harvest rate both change dynamically. This is a more realistic scenario, because the early warnings are calculated in a rolling window and biomass may never be at equilibrium b/c qE keeps changing (chasing a moving target).

Simulate time series and calculate rolling window statistics.

# ==============================
# = Calculate time series of B =
# ==============================
B0 <- 400
B0 <- Burn(B0, n=nBurn, qE=qE[1], sdU=0)
Bvec <- c(B0, rep(NA, nYear-1))
for(i in 2:nYear){
    Bvec[i] <- Bstep(B=Bvec[i-1], qE=qE[i], sdU=sdU)
}

# ==================================
# = Calculate rolling window stats =
# ==================================
sdVec <- rollWarn(Bvec, stat='sd', win=min(50, nYear/5))
ac1Vec <- rollWarn(Bvec, stat='ac1', win=min(50, nYear/5))
redVec0 <- rollWarn(Bvec, stat="redShift", win=min(50, nYear/5)) # holds the full spectrum
redList <- attributes(redVec0)$rsl # this gets the full spectrum formatted for plotting
redVec <- sapply(redVec0, function(x)x[[1]]) # this gets just the spectral slope

Plot the early warnings:

# =======================
# = Plot B and Warnings =
# =======================
par(mfrow=c(2,2), mar=c(1.75,1.75,0.5,1.75), ps=8, cex=1, mgp=c(0.75,0.15,0), tcl=-0.15)
plot(Year, Bvec, type='l')
abline(v=yrPoints, lty='dashed')
text_y <- min(Bvec)+yrProbs*diff(range(Bvec))
text(yrPoints, y=text_y, label=paste0("qE=", qEPoints), pos=c(4,4,2,2))

plot(Year, sdVec, col="forestgreen", type="l", ylab="Standard Deviation")
abline(v=yrPoints, lty='dashed')
plot(Year,ac1Vec, col="blue", type='l', ylab="AR(1)")
abline(v=yrPoints, lty='dashed')
plot(redList, xaxs='r')
par(new=TRUE)
plot(Year, redVec, col=adjustcolor("white",0.5), 
    lwd=3, type='l', ylab="", xaxt='n', yaxt='n'
)
lines(Year, redVec, col=adjustcolor("black",0.5), lwd=0.5)
axis(side=4)
mtext("redShift", side=4, line=0.75)
abline(v=yrPoints, lty='dashed')

Warning of Collapse Conclusion

As you can see, there are definite changes in the early warning statistics near the time that the biomass collapses. However, the standard deviation increases rather drastically right before the shift in some cases, and the redShift is not behaving like I'd expect it to (though, the metric I've cooked up is likely imperfect). Furthermore, it seems odd to me that the shift is not occurring exactly at $qE=0.36$; this is when I'd expect the stock to collapse. It may be that the parameterization I'm using here is slightly different from that described in the Spencer and Collie (1997) text, and I just didn't notice it (though I took care to try to select the parameter estimates reported in the paper).

rBatt/fishyWarnings documentation built on May 26, 2019, 7:45 p.m.