Skip to content

Commit

Permalink
special case last year of recdevs
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChristineStawitz-NOAA
ChristineStawitz-NOAA committed Nov 17, 2023
1 parent 91fefa8 commit 2614278
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 3 deletions.
2 changes: 1 addition & 1 deletion inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,7 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
}
}

recruitment->log_recruit_devs.resize(this->log_devs.size()+1);
recruitment->log_recruit_devs.resize(this->log_devs.size());
if (this->estimate_log_devs) {
for (size_t i = 0; i < recruitment->log_recruit_devs.size(); i++) {
recruitment->log_recruit_devs[i] = this->log_devs[i];
Expand Down
7 changes: 7 additions & 0 deletions inst/include/population_dynamics/population/population.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -391,13 +391,20 @@ struct Population : public fims_model_object::FIMSObject<Type> {
<< this->recruitment->evaluate(
this->spawning_biomass[year - 1], phi0)
<< std::endl;
if(i_dev==this->nyears){
this->numbers_at_age[i_age_year] =
this->recruitment->evaluate(this->spawning_biomass[year - 1], phi0);
/*the final year of the time series has no data to inform recruitment devs,
so this value is set to the mean recruitment.*/
} else{
this->numbers_at_age[i_age_year] =
this->recruitment->evaluate(this->spawning_biomass[year - 1], phi0) *
/*the log_recruit_dev vector does not include a value for year == 0
and is of length nyears - 1 where the first position of the vector
corresponds to the second year of the time series.*/
fims_math::exp(this->recruitment->log_recruit_devs[i_dev-1]);
this->expected_recruitment[year] = this->numbers_at_age[i_age_year];
}
POPULATION_LOG << " numbers at age at index i_age_year " << i_age_year << " is "
<< this->numbers_at_age[i_age_year] << std::endl;
}
Expand Down
4 changes: 2 additions & 2 deletions tests/testthat/test-fims-estimation.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -691,12 +691,12 @@ test_that("run FIMS in a for loop", {
fims$CreateTMBModel()
parameters <- list(p = fims$get_fixed())
obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")

opt <- with(obj, optim(par, fn, gr,
method = "BFGS",
control = list(maxit = 1000000, reltol = 1e-15)
))

report <- obj$report(obj$par)
expect_false(is.null(report))

Expand Down

15 comments on commit 2614278

@Cole-Monnahan-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm late to this party so maybe misunderstanding something here, but the last year has information from the hyperdistribution and we want that uncertainty to be in the model and eventually projections. So setting to the mean seems not right to me.

@ChristineStawitz-NOAA
Copy link
Contributor Author

@ChristineStawitz-NOAA ChristineStawitz-NOAA commented on 2614278 Nov 17, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the input - it's not the last year, it's nyears + 1. Does that still apply? I need to follow up with the group that worked on PR #514 because they removed the first year of the model comparison recdev vector - it's possible everything should have shifted over one rather than removing the first year?

@Cole-Monnahan-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Conceptually any dev should have the hyperdistribution without any data informing it. That can be the initial years, terminal ones, or projection ones. I'm sure there are special cases where you'd want to maybe use a different sigmaR for projections, or fix certain devs, but the latter is best done through mapping IMO, and will accomplish the same.

@iantaylor-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Andrea-Havron-NOAA may be able to describe it better, but my vague memory of the issue in #514 was that age-0 in the initial age structure and the first recdev were applying to the same year and age so one needed to be removed.
Our feeling was that all this would change in M2 when we moved to more flexible range of years for data, estimated recdevs, etc., so at this point the goal was to just get the model to work.

@Cole-Monnahan-NOAA, I totally agree that we can have uninformed devs in the future. If there's further discussion on this, we could make it easier to find by adding to the "Investigate constrain deviations in recruitment" issue #364.

@Andrea-Havron-NOAA
Copy link
Collaborator

@Andrea-Havron-NOAA Andrea-Havron-NOAA commented on 2614278 Nov 17, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dropped off the first year of recdevs as no data were informing it and leaving it in the model as estimated was resulting in gradients that weren't passing tests (higher than 0.0001) in test-fims-estimation.R. After dropping the first observation, gradients became orders of magnitude smaller. I was hoping this was a temporary fix until we explored more correct ways to implement the recruitment devs under constraints or as random effects in M2. I wasn't running into the same issues happening in this branch. Should we revert the changes I made to the log_recruit_dev vector given the problems happening in the fims_vector branch?

@Cole-Monnahan-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That interesting. I'm surprised they don't go to zero after a step or two given that they're independent, and normal so quadratic where is when Newton's method works the best. Usually there's data later on when the fish are still alive so maybe it was just really poorly informed? In that case there can be a lot of correlation too so maybe that's why the optimizer had a hard time? We could get around these gradient issues by taking a step or two using this algorithm: MLE=MLE-(inverse Hessian)*(gradient).

@ChristineStawitz-NOAA
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Cole-Monnahan-NOAA right now recdevs are annual fixed effects, so there is no hyper distribution. We definitely need to re-do this in M2 when we turn on random effects anyway. This fix is now passing all checks, so I'd prefer that we keep it this way rather than revert log_recruit_dev and have to troubleshoot more gradient issues, unless the gradient issues are pointing towards a different problem.

@Cole-Monnahan-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine to leave as is to keep things moving. Just to highlight the point that the hyperdistribution N(0,sigmaR) affects all devs, regardless if they have data informing them, regardless of if using penalized ML or random effect integration.

@ChristineStawitz-NOAA
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But if there's no data informing the final element of the vector because it's for an unobserved final year, wouldn't that value just be drawn from the hyperdistribution rather than informing what sigmaR + other devs are?

@ChristineStawitz-NOAA
Copy link
Contributor Author

@ChristineStawitz-NOAA ChristineStawitz-NOAA commented on 2614278 Nov 18, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if we add the first element of the recdevs vector back, but remove initial NAA at age of recruitment as a parameter instead?

@msupernaw
Copy link
Contributor

@msupernaw msupernaw commented on 2614278 Nov 18, 2023 via email

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Cole-Monnahan-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChristineStawitz-NOAA Yes correct it will just take the hyperdistribution N(0,sigmaR). But that's what we want, because that uncertainty will be propagated into the NAA and throughout the model. Early devs will affect initial NAA and the last ones will affect the projections (if that's desired). This is pretty common practice, at least in my region, and with SS3. For instance take the following example from Pacific hake 2021 assessment Figure 27. You can see the hyperdistribution at the start and end of the vector.

@msupernaw Also correct that you'd get an estimate of 0 if estimated, and have a gradient of zero if you left off rec_devs[0] from the likelihood calculation. But you'd also have a zero second derivative and hence a non positive definite Hessian for that component and thus the whole model. It also would have no uncertainty and thus underestimate uncertain in early NAA.

I think conceptually we should do ~N(0,sigmaR) for all devs as a baseline, and let the user make adjustments from there through e.g. mapping.

image

@msupernaw
Copy link
Contributor

@msupernaw msupernaw commented on 2614278 Nov 21, 2023 via email

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Cole-Monnahan-NOAA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we'll want some flexibility in the initial NAA configurations. In this case if initial_numbers[0] is a freely estimated fixed effect then recruitment_deviations[0] will just take the hyperdistrubtion but be negatively correlated with the fixed effect. I may be lost in the indexing here. Is recruitment[0] in the first year of the model? If so then I don't see the logic of having initial_numbers[0] in that statement, since that is the age structure at the beginning of the first year, and recruitment happens later in the year given those fish.

My main point is that we should think statistically about our recdevs. The N(0,sigmaR) formulation also implicitly assumes that they're exchangeable random effects, and IMO we should strive to meet those assumptions.

I'm done with my busy assessment season so can be more involved and happy to chat about this stuff in person if it's helpful.

@Andrea-Havron-NOAA
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Cole-Monnahan-NOAA, we will be restructuring recruitment as part of M2. I think we will want to follow a similar structure to WHAM/SAM with a random effect on logN opposed to a recruitment dev structure as the former is sparse and the latter is not. We also need to think if we want to include an option to turn off the random effect and use a constraint instead and what this might look like. I would love your input in the next development cycle.

Please sign in to comment.