[Feature]: Implement calculation of effective degrees of freedom

[Rick-Methot-NOAA]

Describe the solution you would like.

A generic approach to implementing effective degrees of freedom has been developed by Monahan, Thorson, Ianelli. See Jan 13, 2025 email.

Describe alternatives you have considered

current estimates of variance of random effects will benefit from this code

Statistical validity, if applicable

Thorson (2024) https://doi.org/10.1002/ecy.4327

Describe if this is needed for a management application

no

Additional context

No response

[Cole-Monnahan-NOAA]

@Rick-Methot-NOAA I have a prototype of this working for an Alaska Pcod example. There are two challenges from a user standpoint.

First, the user must create a version of their model that has no data. The most straightforward way to do this would seem to be setting all data lambda to 0 in the control file. Can you confirm this is the right approach and provide any guidance on how to efficiently do that? For this cod model I get this:

cbind(mod0$likelihoods_used, mod1$likelihoods_used)
                           values lambdas       values lambdas
TOTAL                 2.43416e+02      NA -15.27930000      NA
Catch                 7.20670e-11      NA   0.00000000      NA
Equil_catch           1.71963e-04      NA   0.00000000      NA
Survey               -5.95181e+01      NA   0.00000000      NA
Length_comp           2.22433e+02      NA   0.00000000      NA
Age_comp              5.56440e+01      NA   0.00000000      NA
Recruitment           4.78685e+00       1 -17.15630000       1
InitEQ_Regime         1.86537e+00       1   1.87267000       1
Forecast_Recruitment  0.00000e+00       1   0.00000000       1
Parm_priors           0.00000e+00       1   0.00000000       1
Parm_softbounds       4.27402e-03      NA   0.00427121      NA
Parm_devs             1.82003e+01       1   0.00000000       1
Crash_Pen             0.00000e+00       1   0.00000000       1

Where the second two columns are the model without data. Does this seem right? I'm surprised the Parm_devs are turned off b/c I did not specify a lambda for "12=parm_dev; "

Second, the user must specify which parameters are impacted by some kind of penalty and thus informed by more than the data. Penalty could be a hyperdistribution e.g. on annual deviations (sigmaR and recdevs), or a Bayesian prior or simple parameter penalty. The user needs a way, in R, to subset those parameters. Here is some pilot code to give you a sense of this:

mod0 <- r4ss::SS_output(dir='Model_24.1/')
mod1 <- r4ss::SS_output(dir='Model_24.1_nodata/')
par_names <- mod0$parameters |>
  dplyr::filter(!is.na(Active_Cnt)) |>
  dplyr::pull(Label)

## declare which parameters are random effects or have priors 
randoms <- "Early_InitAge|Main_RecrDev|rwalk|DEVadd"
lrandom <- grepl(randoms, par_names)
par_names2[lrandom]
# rename some of the vectors 
par_names2 <- case_when(
  grepl('InitAge_', par_names) ~ 'InitAge',
  grepl('DEVadd', par_names) ~ 'Size_DblN_ascnd_se_survey',
  grepl('RecrDev', par_names) ~ 'RecrDev',
  grepl('VonBert_K_Fem_GP_1_DEV_MR_rwalk_bnd', par_names) ~ 'VonB_k',
  grepl('L_at_Amin_Fem_GP_1_DEV_MR_rwalk_bnd', par_names) ~'VonB_Lmin'
)

Can you give some guidance on how a user could check/verify that they have all the parameters in there. I think the devs are obvious but the ones with priors are not.

[iantaylor-NOAA]

@Cole-Monnahan-NOAA,
Thanks for working on this. Two quick notes on ongoing stuff in r4ss that my be able to help with the two challenges you identify.

First, this issue opened today by @e-perl-NOAA r4ss/r4ss#1002 could be a user-friendly approach to removing all the data from the likelihood. I think lambdas would be the most straightforward way to implement that function, but it could also be done with negative fleets (which provide expected values for each observation) or negative years (which should be the same as lambda = 0).

Second, I've just created a pull request for some table-making functions, including one that creates a table of information on all the parameters (basically reformatting the parameters table output by SS3). The that function uses Pr_type == "dev" (in https://github.com/r4ss/r4ss/blob/tables/R/table_pars.R#L80-L81) to identify all the deviation parameters parameter. The determination of which sigma applies to which parameter needs to be improved, but I think that Pr_type != "No_prior" should identify parameters with either an informative prior or a deviation penalty.

[Rick-Methot-NOAA]

I think it may be possible to co-opt the code that does depletion_fleet to also do the needed calculations. See in SS_readcontrol.tpl: SS_Label_Info_4.14.1 #Adjust the phases to negative if beyond turn_off_phase

With depletion_fleet, SS3 adds 1 to the phase in which each parameter becomes active except for the one depletion_fleet which keeps lambda=1 for the "survey" of that fleet. So, in phase 1 there is olye that one datatype for that one fleet that has a positive lambda.
I think we could so something similar to get everything turned off. In fact, as is you might run a model with lambda=0 for a depletion fleet, then read that par file back in with depletion fleet now turned on. Just brainstorming here....

Regarding parmdevs -logL. If all other data types have lambda = 0, then there are no data to tell the model that any dev has a value different than 0.0, so the logL for the parmdevs will be zero.

[Rick-Methot-NOAA]

There already is a concept of depletion_type implemented if a survey is depletion_fleet. So, I think we could create a depletion_type==3 that does what is need for the df calculations.

  if (Svy_units(f) == 34) //  special code for depletion, so prepare to adjust phases and lambdas
  {
    echoinput << "# survey: " << f << " " << fleetname(f) << " is a depletion fleet" << endl;
    depletion_fleet = f;
    depletion_type = Q_setup(f, 2);
    if (depletion_type == 0)
      echoinput << "link_info=0; add 1 to phases of all parms; only R0 active in new phase 1 (same as 3.24 logic)" << endl;
    if (depletion_type == 1)
      echoinput << "link_info=1  only R0 active in phase 1; then exit;  useful for data-limited draws of other fixed parameter" << endl;
    if (depletion_type == 2)
      echoinput << "link_info=2  no phase adjustments, can be used when profiling on fixed R0" << endl;

[Cole-Monnahan-NOAA]

@iantaylor-NOAA I tried negative years but SS3 catches that as an error (rightfully so). I did not try negative fleets. Perhaps that would work.

@Rick-Methot-NOAA We evaluate this no-data model from the .par file from the full model (-maxfn 0). So the dev penalties should be non-zero b/c the devs are non-zero. My question about the lambdas was whether there was a way to tell by the output that it's working correctly. I have never used them before. Does this seem right for a model with 2 fleets?

# 10=recrdev; 11=parm_prior; 12=parm_dev; 13=CrashPen; 14=Morphcomp; 15=Tag-comp; 16=Tag-negbin; 17=F_ballpark; 18=initEQregime
#like_comp fleet  phase  value  sizefreq_method
 1 1 1 0 0
 1 2 1 0 0
 2 1 1 0 1
 2 2 1 0 1
 3 1 1 0 1
 3 2 1 0 1
 4 1 1 0 1
 4 2 1 0 1
 5 1 1 0 1
 5 2 1 0 1
# 6 1 1 0 1
# 6 2 1 0 1
 7 1 1 0 1
 7 2 1 0 1
 8 1 1 0 1
 8 2 1 0 1
 9 1 1 0 1
 9 2 1 0 1

Perhaps it's easier for me to turn over my prototype and let you guys play with it and sort out the best way forward?

[Rick-Methot-NOAA]

lambdas do not turn off the data. They do multiple the -logL by lambda. This means that the -log(s) terms are also being multiplied by lambda and that is probably not what you want.