R1 offset, R regime, and historical fishing mortality parameters - Public Forums - Stock Synthesis

I definitely don't have all the answers, but I'll post a few things here and hopefully Rick can correct anything incorrect that I've said in addition to filling in the gaps that I leave, although he's on leave this week, so that might take a little while.

Starting with the old R1 offset, my understanding is that for an annual model, the equilibria should be equal to

Unfished equilibrium (associated with startyr - 2):

for a = 0: N_{a,startyr-2} = R0

for a = 1, 2, ..., Amax - 1: N_{a,startyr-2} = N_{a,startyr-2} * exp(-M) [exponential decay from R0]

for a = Amax: more complex stuff happens to account for Amax as plus group.

Initial equilibrium (associated with startyr - 1):

If R1_offset = 0 then the numbers at age should exactly match the unfished equilibrium. If R1_offset is non-zero, then R1 = R0*exp(R1_offset). Otherwise, the rest is the same.

for a = 0: N_{a,startyr-1} = R1

for a = 1, 2, ..., Amax - 1: N_{a,startyr-1} = N_{a,startyr-1} * exp(-M) [that is, exponential decay from R1 instead of R0]

for a = Amax: more complex stuff happens to account for Amax as plus group.

Time series (startyr and beyond):

for a = 0: N_{a,y} = f(R0, SB_y, h, etc.) * exp(-0.5*b_y*sigmaR^2 + recdev_y),

where f() is the stock-recruit function, b_y is the bias adjustment for year y and recdev_y is the recruitment deviation for year y.

for a = 1, 2, ..., Amax - 1: N_{a, y} = N_{a-1, y-1} * exp(-Z_{a-1,y-1})

for a = Amax: N{a,y} = N_{a-1, y-1} * exp(-Z_{a-1,y-1}) + N_{a, y-1} * exp(-Z_{a,y-1})

What's different with 3.30 and the SR_regime option?:

For any given year, including startyr - 1, the R0 is replaced by R0*exp(SR_regime[y]). So to mimic the R1_offset, you need to put a block on SR_regime for y = startyr - 1. But if SR_regime has some change during your main time series, that change will be filtered through the stock-recruit relationship into an impact on the numbers at age for whatever years are impacted. This can also apply to the forecast. A block on SR_regime is also easier than some of the old dummy environmental variables that were created in the past to adjust recruitment up and down for long periods.

Hopefully this makes some sense and I haven't made too many mistakes.

-Ian

This is great information to have on R1 offset and R regime. Thank you, Ian. It would be great to know, also, the initial conditions equations when initial F is used (at the same time as R regime and/or R1 offset). 

Cheers,

Carey

Here is some additional information about the initial equilibrium calculations in SS.  Note that the dynamics described by Ian above are also described in section 1.1 of the appendix to Methot&Wetzel (2012) Technical Description of SS.

1.  recruitment penalty:  This was not included in the technical description document.  In 3.24, the R1_offset was penalized to keep it from being too large.   This was done by calculating an approximation of the average age of fish:

//  ProgLabel_17.2  calc an ave_age for the first gp as a scaling factor in logL for initial recruitment (R1) deviation
        if(Do_AveAge==0)
        {
          Do_AveAge=1;
          ave_age = 1.0/natM(1,gg,nages/2)-0.5;
        }
 

Then using this in objective function:

      recr_like += 0.5 * square( log(R1/R1_exp) / (sigmaR/ave_age) );
 

In 3.30.12, this penalty is now separated into a separate component of the total -logL:

      regime_like = 0.5 * square( log(R1/R1_exp) / (sigmaR/ave_age) );
and it can now be modified by a separate lambda (lambda change #18).  So you can eliminate this penalty by giving it a lambda=0.0.

2.  Program flow:  The attached file contains the population dynamics code for SS3.30.12.  It has three main FUNCTIONs:  Initial Equilibrium,  Time_Series, Do_Equil_Calc where Initial Equilibrium calls Do_Equil_Calc first with Z=M and R=R0, then with Z=M+init_F and R=R1.   The Do_Equil_Calc function is used repeatedly in SS.  FIrst for the initial equilibrium, then in benchmark and whereever equilibrium calculations are needed.

3.  Plus Group:  Do_Equil_Calc processes the population dynamics on an age by age basis out to 3X maxage, as described in Methot&Wetzel.  Whereas a simple age-structured model could simply get the analytical calculation for the numbers-at-age in the plus group,  SS includes more complex possibilities such as multiple areas and hermaphroditism, so the analytical solution is not universally possible.  So as SS does the Do_Equil_Calc calculations out to 3X maxage, it takes into account all aspects of the fleet, area, biological structure that occur during the time series.

However, the body size of fish in the plus group for the initial equilibrium is more difficult if body size is calculated from growth equations within SS and fish do not attain L_infinity by the time they reach maxage.  In the growth module of SS, it goes out to 3x maxage to calculate mean size-at-age and then assigns the mean size at maxage to equal a weighted mean of size at maxage out to 3x maxage.  The weighting should approximate the mortality, but mortality is not know at the time of the growth calculations, so approximations are available in the growth setup.  Then if you have time-varying growth active, SS will update the mean size at maxage each year according to the weighted mean of the mean size of fish already in the plus group and the mean size of fish just growing into the plus group.  So if M+F during the time series is greater than the M+F in the initial equilibrium calculation, you will see the mean size in the plus group decline over time.

It would be more correct if SS did this updating of the mean size in the plus group even when time-varying growth was not active.  Alas, simply making that happen would make SS always run as slowly as it does when growth is time-varying (because the age-length key needs constant updating), so a more surgical approach will be needed if we attempt to make this improvement in the SS approach. 

4.  Spawner-recruitment in the initial equilibrium:   There is a switch in SS that either sets R1 equal to some specific recruitment level then calculates the initial equilibrium scaled by that R1, or that first uses M+F and Do_Equil_Calc to calculate spawners per recruit, then uses that spawners per recruit to calculate the equilibrium recruitment that would come from that.  Originally, SS had only the first option.  Then when more applications started using late start years and high initial catch levels, it was more necessary to admit that the initial equilibrium F had already begun to affect the equilibrium recruitment level.  Note that the second option does identical calculations as are done in benchmark to calculate MSY.  The corollary of this is that initial equilibrium catch cannot exceed MSY, and that as initial equilibrium catch approaches MSY, then initial equilibrium catch is not longer monotonic with initial F.  So beware.

5.  Testing:  I have created test cases in which I first found MSY and Fmsy, then set up an example in which initial equilibrium catch and each year's catch was set to MSY, and the forecast was set to use Fmsy.  In this case, the population and its age-structure remained identical from the initial equilibrium through each year of the time series and the forecast and in the calculated MSY benchmark.  I will try to find the example and post it.

Rick