Age-1 anomalies in age-1 leaks from unkriged outputs for age-2+ calculations (for 2019)

[rebeccathomas-NOAA]

I have noticed two things that are related to age-1 numbers, and I am not sure if they are related . However, I will describe them as two separate issues to keep things simple. One of the metrics I look at is the age-1 "leaks" from unkriged, ("transect') reports. This is one of the things that goes into calculation of the age-1 index used by stock assessment. It's cell B43 in sheet 3 of the transect_length_age_biomass_report.xlsx . I have looked at EchoPop outputs both from my running of the data as well as the results from the transect_length_age_biomass_report (e.g. at https://docs.google.com/spreadsheets/d/1Ev91Qbh3iBB6KlRIq78-G4FEpyXcHkWX/edit?gid=958259491#gid=958259491), and there is about a 10% difference between the EchoPro output of ~0.00787 (e.g. in un-kriged_len_age_biomass_table) and the EchoPop output of ~0.00869. Do we know where this difference is coming from?

[brandynlucca]

@rebeccathomas-NOAA

I have traced down the sources for this discrepancy.

There was a change to how proportions were stored in the xarray.DataArray and xarray.Dataset objects in #464 and #474. One consequence of this that appears to have leaked through the testing cracks was how the age-1 number proportions were being computed. Previously, the within-group proportions were used directly as opposed to the overall proportions. However, these are no longer both stored; only the overall proportions are stored within the objects. Consequently, the number proportions computed for age-1 fish within each stratum were not appropriately normalized, which resulted with much lower proportions (and therefore elevated abundance estimates). This has been corrected in #486.

The second is less a bug and more an oversight in the workflow design. The Echopop-EchoPro comparisons were done using the general transect/kriging population estimates files rather than the table reports. Specifically with the snippet concerning the distribution of unaged biomass estimates over age distributions:

https://github.com/OSOceanAcoustics/echopop/blob/a6389c2f05de6a8106f612a71b3ef033b9ac886b/echopop/workflow_examples/year_specific/hake_2019.py#L877-L884

dict_ds_transect_biomass_table[
    "standardized_unaged"
] = apportionment.distribute_unaged_from_aged(
    population_table = dict_ds_transect_biomass_table["unaged"],
    reference_table = dict_ds_transect_biomass_table["aged"],
    stratum_dim = "stratum_ks",
    impute = False 
)

Here, there is some leakage when impute = False due to a mismatch in the presence of unaged fish at a particular $\ell$ that may not have any values in the aged dataset. This results in NaN (due to divide-by-zero). So that results in a loss of biomass which can vary in impact depending on the dataset. So:

>>> dict_ds_transect_biomass_table["standardized_unaged"].sum()
<xarray.DataArray 'biomass' ()> Size: 8B
array(1.5055274e+09)

Setting impute = True and defining just age_bin as the imputation variable via impute_variable = ['age_bin'], i.e.,

dict_ds_transect_biomass_table[
    "standardized_unaged"
] = apportionment.distribute_unaged_from_aged(
    population_table = dict_ds_transect_biomass_table["unaged"],
    reference_table = dict_ds_transect_biomass_table["aged"],
    stratum_dim = "stratum_ks",
    impute = True,
    impute_variable = ['age_bin'],
)

returns the original, pre-standardized biomass:

>>> dict_ds_transect_biomass_table["standardized_unaged"].sum()
<xarray.DataArray 'biomass' ()> Size: 8B
array(1.51563142e+09)

So the imputation helps avoid this leakage.

Once that is merged (with an updated version of at least hake_2019.py, this issue will be closed. Feel free to re-open this issue if there is a related observation due to a bug elsewhere.

[rebeccathomas-NOAA]

Thanks for looking into this, Brandyn! Your solution in issue 483 (age-2 + abundances being off) seems to have done the trick, and the age-2+ abundances match for me now. Thanks!
However, when I run the latest hake_2019 with your version of the data, in the transect_length_age_biomass_report I still get age-1 biomass to be about 10% high at ~0.00869. I am not sure if it is helpful, but I do note that I do get the same result as you here:
dict_ds_transect_biomass_table["standardized_unaged"].sum() <xarray.DataArray 'biomass' ()> Size: 8B array(1.51563142e+09)
Thanks,
RT

[brandynlucca]

The source:

A = df_nasc.groupby(["stratum_ks"])["biomass"].sum()
B = dict_da_weight_proportion["aged"] / dict_da_weight_proportion["aged"].sum(dim = ["length_bin", "age_bin", "sex"])
(A.to_xarray() * B)["proportion"].sel(age_bin = 1).sum() * 1e-9

So currently, Echopop runs the table imputation for the transect-based, age-length distributed biomass estimates. EchoPro did not, and instead replicated the above code logic. @rebeccathomas-NOAA is mulling over the ramifications for which way to implement.

[rebeccathomas-NOAA]

After discussion, I believe that the difference between the unkriged age-1 biomass ("age-1 leaks") between EchoPop and EchoPro is real, but does not need to be "fixed." My understanding is that in Chu's EchoPro code, the unaged fish were not used in the conversion of biomass to biomass at age and length for the unkriged data. However, the unaged fish were used in EchoPop (and imputed to age) when converting biomass to biomass at age for the unkriged results. As the length & age samples are more limited than the length-only samples, we can see some differences between proportions at age when just-length samples are used and imputed to age vs. when the aged samples only are used. Both methods, EchoPro and EchoPop, used aged and unaged fish in the post-kriging distribution of biomass. This change in unkriged biomass proportions-at-age will not propagate through to the kriged biomass-proportions-at-age, as the total biomass is first kriged and then apportioned.

[brandynlucca]

@rebeccathomas-NOAA

So thinking about this a bit more...

Currently, impute_kriged_table uses a nearest neighbor imputation, whereby the proxy length bin, $\ell^\star$, minimizes the absolute length difference from the empty bin, $\ell_0$. So from the set of valid nonempty aged bins, $L_\alpha$, $\ell^\star$ is defined by:

$$ \ell^\star = \min \left( \arg \min_{\ell \in L_\alpha} |\ell - \ell_0| \right). $$

So this finds the $\ell$ from the aged data that is closest to $\ell_0$. In the case of ties, it then takes the minimum of those lengths. So if there are two $\ell$, the larger one is effectively discarded. For example, if an empty 70 cm bin is the exact midpoint between a 68 and 72 cm non-empty bin, we ignore the 72 cm bin entirely and impute the value strictly from the 68 cm fish. This introduces an operational bias due to the tie-breaker criterion that systematically skews the imputed population estimates toward smaller fish. This also makes the end-result more vulnerable to sampling noise (e.g., a very low sample size or random 'dip' in the distribution that just so happens to exist at $\ell^\star$.

Alternatively, a distance-weighted pooling scheme may address these limitations. Instead of arbitrarily selecting a single nearest neighbor and discarding adjacent data, this approach would identify the immediate bounding interval surrounding $\ell_0$ to preserve the distribution on both sides. So from the set of non-empty aged bins $L_\alpha$, we can identify the closest lower and upper neighboring proxies relative to $\ell_0$:

$$ \ell_\mathrm{lower} = \max \{ \ell \in L_\alpha | \ell < \ell_0 \}, $$

$$ \ell_\mathrm{upper} = \min \{ \ell \in L_\alpha | \ell > \ell_0 \}. $$

The absolute 'distances' from $\ell_0$ are therefore:

$$ d_\mathrm{lower} = \ell_0 - \ell_\mathrm{lower}, $$

$$ d_\mathrm{upper} = \ell_\mathrm{upper} - \ell_0. $$

So the conditional proxy set $L^\star \subseteq L_\alpha$ (valid proxy values in the nonempty aged length bins):

$$ L^\star = \begin{cases} \{\ell_{\text{lower}}, \ell_{\text{upper}}\} & \text{if both neighbors exist} \\ \{\ell_{\text{lower}}\} & \text{if } \ell_0 > \max(L_\alpha) \\ \{\ell_{\text{upper}}\} & \text{if } \ell_0 < \min(L_\alpha) \end{cases} $$

We can then use the distances to prevent far away bins from disproportionately influencing the empty bin (i.e., the normalized inverse-distance weights) for each $\ell \in L^\star$:

$$ w_\ell = \begin{cases} \frac{d_{\text{upper}}}{d_{\text{lower}} + d_{\text{upper}}}, & \text{if } \ell = \ell_{\text{lower}} \text{ and both neighbors exist} \\ \frac{d_{\text{lower}}}{d_{\text{lower}} + d_{\text{upper}}}, & \text{if } \ell = \ell_{\text{upper}} \text{ and both neighbors exist} \\ 1, & \text{if } L^\star \text{ contains only a single proxy} \end{cases} $$

So now the initial unallocated unaged biomass at $\ell_0$, $B_0$, can be redistributed across the age-classes $\alpha$ by incorporating the inverse distance weights into a mass-weighted allocation:

$$ B_{\ell_0, \alpha} = B_0 \times \frac{ w_{\text{lower}} B_{\ell_{\text{lower}}, \alpha} + w_{\text{upper}} B_{\ell_{\text{upper}}, \alpha} }{ w_{\text{lower}} \sum_{j} B_{\ell_{\text{lower}}, j} + w_{\text{upper}} \sum_{j} B_{\ell_{\text{upper}}, j} } $$

This joint-weighting approach helps eliminate that operational bias toward smaller fish while still ensuring that the allocated biomass across all $\alpha$ still sums to the original grand total biomass. The IDW can also be removed where ties from the original nearest neighbor approach are instead weighted by their pooled contribution. For instance, say the lower and upper bins are 10 and 40 kg, respectively, at 68 and 72 cm where $\ell = 70 \mathrm{cm}$. Then:

$$ B_{\ell_0, \alpha} = B_0 \times \frac{B_{\ell_{\text{lower}}, \alpha} + B_{\ell_{\text{upper}}, \alpha}}{\sum_{j} B_{\ell_{\text{lower}}, j} + \sum_{j} B_{\ell_{\text{upper}}, j}} $$

which for this example becomes:

$$ B_{70, \alpha} = B_0 \times \frac{ B_{68, \alpha} + B_{72, \alpha} }{ 10 + 40 } = B_0 \times \frac{B_{68, \alpha} + B_{72, \alpha}}{50} $$

So if, say, $\alpha = 3$ where the two equivalent bins are 7 and 30 kg, respectively, then:

$$ B_{70, \alpha = 3} = B_0 \times \frac{7 + 30}{10 + 40} = B_0 \times \frac{37}{50} $$

Or something to that effect. Just some food for thought for shoring up some of the potential issues with how this imputation is being handled.