Why does recruitment show a StdDev in forecast years when using fixed geometric mean?

[Marta Cousido]

Dear SS Experts,

I am working on projections in Stock Synthesis and have a question about how recruitment uncertainty is handled in the forecast.

In our current assessment,  we decided not to use the stock–recruitment (SR) relationship to estimate recruitment for the intermediate and forecast years. Instead, we fixed recruitment to the geometric mean of the last 10 years, as the SR-based estimate for 2026 was unrealistically high compared to recent estimates.

Specifically, in the forecast file we used option 2 (Multiplier on virgin recruitment) for the forecast recruitment adjustment section.

When reviewing the Report file, we notice that recruitment in these years still has a  StdDev estimate. This confuses us, since recruitment is not being estimated. 

Could someone clarify how  SS computes a standard deviation for recruitment under this setup?

Any clarification on how to correctly interpret this StdDev in projections would be very helpful.

Thanks in advance!

Marta Cousido

[Ian Taylor - NOAA Federal]

Hi Marta,

Rick will know better, but if the model works as I expect, the StdDev for the forecast recruitment should equal the StdDev of Recr_Virgin derived quantity (or scaled by any multiplier you use). The uncertainty in Recr_Virgin should be a direct function of the estimated uncertainty in the SR_LN(R0) parameter. Therefore, if that parameter is fixed, I would expect the forecast recruitments to have negligible or 0 uncertainty.

Let me know if the estimates don't match this description.

-Ian

--
You received this message because you are subscribed to the Google Groups "SS3 - Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ss3-forum+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/ss3-forum/de24bd79-9d7d-471e-915a-8cb0db33b08en%40googlegroups.com.

[Richard Methot - NOAA Federal]

Hello Marta,

There are several factors to consider here that will affect the variance of the forecast recruitments.

1.  as Ian wrote, Virgin recruitment is an estimated quantity with variance, so a multiple of that value will also have variance.

2.  Check this setting in the control file:   0 #_forecast_recruitment phase (incl. late recr) (0 value resets to maxphase+1).  You should change it to -1 if you want the forecast recruitment devs to not be estimated. 

3.  There are other options available to automate the setting of forecast recruitments to the mean of the last 10 years.

a.   use the forecast file specification that allows you to specify a range of years over which to average recruitment.  This option was changed with 3.30.23, so let us know if not clearly described in manual and forecast.ss_new.

or

b.  put a block on the parameter for SR_regime over the range of years:  endyr-X to 9999

and let this block offset be estimated.  This will cause the base recruitment to shift beginning in the year endyr-X and persist throughout the forecast.

Rick

[Marta Cousido]

Dear Ian and Rick,

Thanks very much for your quick answers.

In our SS model, R0 is estimated and therefore we obtain uncertainty associated with Recr_Virgin (Value= 457611, StdDev = 21306.6).

For the forecast, we are using option 2 (multiplier on virgin recruitment), with a multiplier of 0.43. Therefore, the forecast recruitment becomes:

Recr_2026 = 196407, StdDev = 118198.

Then, we expect that the StdDev of the forecast recruitment should be:

21306.6 × 0.43 ≈ 9161,

But this value is much smaller than the reported StdDev (118198).

Therefore, I suspect that I am misunderstanding how SS computes the uncertainty of forecast recruitment under option 2, and that something more complex than a simple scaling of the Recr_Virgin uncertainty is occurring internally.

On the other hand, regarding the forecast recruitment deviations, I also became confused when revisiting the control file. There we set the following line:

0 #_forecast_recruitment phase (incl. late recr)

From the manual description, I interpreted this as allowing forecast recruitment deviations to be estimated after convergence of the rest of parameters. However, when checking the Report.ss output, the forecast years show dev = 0 and no forecast recruitment deviations appear to be estimated. Are we missing something?

Regarding Rick’s option 3a, we are using version 3.30.18 because it is the version used during the benchmark development of our SS model. I checked the manual for 3.30.18  version and the option to use the mean recruitment over a range of last years is clearly described there.

We originally selected option 2 because if we use a multiplier on virgin recruitment we can fix the future recruitment at the geometric mean of recent recruitments, which is commonly used instead of the arithmetic mean. However, since there are no extreme recruitment estimates during the averaging period in our case, using the arithmetic mean through may be reasonable too.

Still, we would like very much to better understand how the forecast recruitment StdDev is computed under option 2.

Thanks again for your help.

[Richard Methot - NOAA Federal]

Marta,

You are correct that the calculation of std.dev. on forecast recruitment is more complex, but that complexity is entirely in the ADMB approach to getting std.dev. on derived quantities; it is not an explicit component of the SS3 code.

When you use the "0" option for forecast recruitment devs, it is normal for the devs to all have a value of 0.0 because there are no data in the forecast to inform that model that any dev should have a value other than 0.0.  However, the forecast devs will still get a std.dev. and that standard dev will be very close to sigmaR because that is the penalty that is keeping those devs close to 0.0.  Ian and I explain this feature in our 2011 paper on recruitment bias:  Can. J. Fish. Aquat. Sci. 68: 1744–1760 (2011) doi:10.1139/F2011-092.

That paper also explores the implications of using arithmetic vs geometric mean recruitment.  SS3 is designed to achieve a high degree of internal consistency when using arithmetic mean recruitment.  The MSY and BMSY values calculated in the reference points section are based on arithmetric mean recruitment.  

We encourage you to update your model to the most recent version to take advantage of some bug fixes and new features.

RIck

[Ian Taylor - NOAA Federal]

Hi All,

In case it's helpful, I calculated the uncertainty of a forecast deviation with forecast option 2 as a function of the estimated uncertainty in the unfished equilibrium recruitment and the uncertainty in the forecast deviation (which equals sigmaR) and got a really close match to the estimated uncertainty of the forecast recruitment.

I used the "simple_small" model included in the r4ss package which can be found at https://github.com/r4ss/r4ss/tree/main/inst/extdata/simple_small but also if you have the package installed you should have the files at file.path(system.file("extdata", package = "r4ss"), "simple_small"). 

I set the forecast file to have the following settings to match Marta's described setup:

2 # Forecast base recruitment:  0= spawn_recr; 1=mult*spawn_recr_fxn; 2=mult*VirginRecr; 3=deprecated; 4=mult*mean_over_yr_range

# for option 4, set phase for fore_recr_devs to -1 in control to get constant mean in MCMC, else devs will be applied

0.43 # multiplier on base recruitment

Then I calculated the uncertainty as follows, where the external calculation of the SD resulted in a value of 1955.483 compared to a model estimate of 1952.5, a difference of less than 0.2%.

# read model output into R

model <- SS_output("inst/extdata/simple_small_forecast_recdev2/")

# mean and SD of unfished equilibrium recruitment

(R0_mean <- model$derived_quants["Recr_Virgin", "Value"])

# [1] 7535.61

(R0_sd <- model$derived_quants["Recr_Virgin", "StdDev"])

# [1] 418.526

# mean and SD of exp(dev) for a forecast year

(exp_dev_mean <- exp(model$parameters["ForeRecr_2032", "Value"])) # should equal exp(0) = 1

# [1] 1

(exp_dev_sd <- model$parameters["ForeRecr_2032", "Parm_StDev"]) # should equal sigmaR

# [1] 0.6

# calculate SD of the product of R0*exp(dev) based on formula at

# https://en.wikipedia.org/wiki/Distribution_of_the_product_of_two_random_variables#Variance_of_the_product_of_independent_random_variables

# with 0.43 multiplier applied

0.43 *

  sqrt(

    (R0_sd^2 + R0_mean^2) *

      (exp_dev_sd^2 + exp_dev_mean^2) -

      (R0_mean^2) * (exp_dev_mean^2)

  )

# [1] 1955.483

# uncertainty in forecast recruitment as derived quantity

model$derived_quants["Recr_2032", "StdDev"]

# [1] 1952.5

--

You received this message because you are subscribed to the Google Groups "SS3 - Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ss3-forum+...@googlegroups.com.

To view this discussion visit https://groups.google.com/d/msgid/ss3-forum/8803df1e-17d2-472e-bd34-3a3f27f7aa55n%40googlegroups.com.

[Marta Cousido]

Hi Rick and Ian,

Thank you very much for your support and patience helping me to understand what was going on with my forecast.

With your latest explanations, I now understand the variance computation and the reasons of its formulation properly. Using Ian’s code, I have been able to obtain nearly the same variance (less than 1% difference compared to the SS-reported value).

I will follow your recommendation and adapt my model files to run it on the most recent SS version.

Thanks again for your help.

Best regards,
Marta