Reporting metrics

[kalebentley]

I'd like to initiate this thread to summarize and discuss the current list of metrics that are generated for Columbia River chum salmon "populations" (see below) and reported to various portals including WDFW's adult and juvenile estimate databases (SaSI and JMS, respectively, and summarized via SCoRE), Coordinated Assessment (CA), and project-level reports (e.g., Grays River Chum Salmon Enhancement Program).

Although worthy of a larger discussion, I've mentioned to @Hillsont that we may be able to replace the current adult "run-reconstruction" process with estimates that are generated from the IPM. However, this requires that the IPM generate the metrics we currently report. As @Hillsont said, we don't want to create a "Frankenstein" reporting process and so it would be great if we can get the IPM to generate everything we want.

Below is a list of the current reporting metrics for Columbia River chum salmon that @Hillsont, @BradGarnerWDFW, and I pulled together late last week.

Estimates of abundance

• Total adult/spawner abundance by return year
• Adult/spawner abundance by composition (totals and percentages)
  • Origin (i.e., natural/program -- see list below);
  • NOSA (Natural-origin Spawner Abundance by TRT populations)
  • Age (total age);
  • Sex (male, female)
• Brood Table (Adult Abundance by age population and brood year/age
• Juvenile abundance for monitoring sites

Other annual metrics (estimates/data needed in parentheses)

• Mining Rate (Broodstock and spawners by populations and return year)
• Progeny per Parent (Adult-to-adult; R/S) by origin
  • Hatchery = Broodstock, Recruits – adults by brood year
  • Natural = Spawners and corresponding recruits 3-6 years later
  • Duncan Channels = same as Natural
• pHOS (total spawners & hatchery spawners)
• pNOB/pHOB (broodstock by origin)
• PNI (pNOB & pHOS)
• SARs by origin (ignoring strays for now)
  • Hatchery = Released hatchery fry and corresponding HORs (adults) 3-6 years later
  • Natural = estimated NOR fry and corresponding NORs (adults) 3-6 years later
  • Duncan = estimated NOR fry and corresponding NORs (adults) 3-6 years later

The above metrics will need to be generated by year (duh!) and "reporting population". There are two general groups of "reporting populations". The first population grouping is the same as the populations in our IPM datasets (there are currently 12). The second population grouping is based on lower Columbia River TRT designations (see pg. 26 of LCR Recovery Plan). For populations that are currently in the IPM, this includes Grays/Chinook (partial estimate as no estimates from Chinook), Washougal/I-205 (not sure if this is considered partial), and Lower Gorge (all mainstem & WA tributaries combined). We also report on the Upper Gorge population using counts at Bonneville Dam along with some metric of fall backs (@Hillsont will have to help clarify this one) but perhaps not important since Upper Gorge data are not in our model currently (but perhaps could be?).

As indicated above, some of the metrics need to be generated for both natural origin and "program" fish. Here's a list of the existing programs for chum salmon in the Lower Columbia (NOTE: I've denoted the method used to identify these fish in parentheses but not important for this thread):

• Duncan Hatchery (thermal 01 - 13, PBT 14-current)
• Lewis Hatchery (PBT all)
• Grays Hatchery (PBT brood, thermal all juveniles)
• Big Creek Hatchery - Oregon - (PBT brood, thermal all juveniles)
• Duncan Channels (strontium 03-13, PBT)
• Petersen RSI (PBT brood) <- not currently part of IPM but should discuss integrating at some point

To move this discussion along, @ebuhle will you help me designate each metric that I've listed above into one of the three categories?

1. already generated by IPM
2. currently not generated but could "easily" be added (e.g., with a few lines of code in the generated quantities block),
3. currently not generated and can't without modifying the model (i.e., would require a non-trivial amount of work)

[ebuhle]

Thanks for this detailed list @kalebentley, this is great. I agree there's a larger conversation to be had about reporting; e.g., although I'm not super familiar with any of the data portals you mention, I wouldn't think we'd want to report IPM output in place of the "raw" data or sample-based estimates that are needed as inputs to modeling. For now I'll stick to categorizing these metrics vis a vis the current IPM.

Estimates of abundance

* Total adult/spawner abundance by return year

Done. Estimates available as posterior distribution of states ("smoothing distribution") or posterior predictive distribution (PPD) of replicate data (i.e., includes observation noise).

* Adult/spawner abundance by composition (totals and percentages)
  - Origin (i.e., natural/program -- see list below);

Done, with limitations (see, e.g., here). Currently we estimate natural vs. hatchery, except hatchery includes Duncan Channel. Program-specific estimates will come from modeling hatcheries and straying, which is nontrivial but on the agenda.

  - NOSA (Natural-origin Spawner Abundance by TRT populations)

Easy. Aggregating by TRT populations is no problem except insofar as we don't have data for all the subpopulations.

  - Age (total age);

Done.

  - Sex (male, female)

Done, serendipitously!

* Brood Table (Adult Abundance by age population and brood year/age

Done.

* Juvenile abundance for monitoring sites

Done.

Other annual metrics (estimates/data needed in parentheses)

* Mining Rate (Broodstock and spawners by populations and return year)

Done. In fact, this is precisely the B_rate parameter that salmonIPM estimates as a computational kludge, so it would be great if it had practical utility.

* Progeny per Parent (Adult-to-adult; R/S) by origin

  * Hatchery = Broodstock, Recruits – adults by brood year

Nontrivial. Will require modeling hatcheries / straying. Also, I had imagined the hatchery life cycle going from smolt release to adult return, so this warrants further discussion if we need hatchery R/S.

  * Natural = Spawners and corresponding recruits 3-6 years later
  * Duncan Channels = same as Natural

Easy. Not yet calculated as such, but it's trivial (ignoring straying for now).

* pHOS (total spawners & hatchery spawners)

Done, with qualifications about Duncan Channel mentioned above.

* pNOB/pHOB (broodstock by origin)

Nontrivial? I've been assuming only NORs are taken as broodstock, although now I realize I haven't actually checked that against the data. In principle I think it wouldn't be too difficult to estimate a separate annual B_rate for natural- and hatchery-origin spawners, informed by the origin composition of bio_data samples with location != disposition, but I've never considered it before.

Never mind, I realized this refers to broodstock origin-composition aggregated by recipient hatchery, not by donor population. So it will require modeling hatcheries / straying, although see note above about starting the hatchery life cycle with M vs. S. (That said, 21% of the B_take_obs from Grays MS is hatchery-origin, so the assumption of NOR-only broodstock could be problematic. All other pops range from 1-9%. Ugh.)

* PNI (pNOB & pHOS)

Nontrivial. (See above.)

* SARs by origin  ** (ignoring strays for now) **

  * Hatchery = Released hatchery fry and corresponding HORs (adults) 3-6 years later

Nontrivial. Will need model of hatcheries / straying.

  * Natural = estimated NOR fry and corresponding NORs (adults) 3-6 years later
  * Duncan = estimated NOR fry and corresponding NORs (adults) 3-6 years later

Done.

We also report on the Upper Gorge population using counts at Bonneville Dam along with some metric of fall backs (@Hillsont will have to help clarify this one) but perhaps not important since Upper Gorge data are not in our model currently (but perhaps could be?).

This is very interesting. I've been wondering about the Upper Gorge population since I first saw your PPT slides. Maybe this should go on the list of potential populations to include that we talked about recently?

• Petersen RSI (PBT brood) <- not currently part of IPM but should discuss integrating at some point

Likewise, I don't know what this is but it sounds interesting.

• Big Creek Hatchery - Oregon - (PBT brood, thermal all juveniles)

I'll just note that there were a total of 30 recoveries from Big Creek Hatchery (out of 1605 total hatchery recoveries), so I doubt we'll be able to do much here.

[ebuhle]

Hey Eric,

Would you be able to (easily) generate & share a plot of fry per female spawner for the Hamilton Springs spawning channel and/or an estimate of fry capacity (median, 95% CIs)?

In short, we are running the fry trap and are interested in generating a prediction of fry outmigrants. The estimate of female spawners in fall 2021 was 577 (not sure the uncertainty).

Right now, I am generating an inseason prediction of outmigants using observed daily trapping data and fitting a log-normal outmigration timing/distribution (using MLL) in Excel (using Solver) by varying the total abundance and peak outmigration day (mean & SD). It’s mostly fine for our purposes but given it is earlier in the outmigration the fit isn’t “good” without providing some strong constraints…

Anyhow, again, this isn’t super important so again if this my request isn’t easy then no worries.

Kale

Hi @kalebentley, yes, it would be very easy to generate a plot (or a numerical summary, which sounds like it would be more directly useful if I'm understanding you) of fry per female for each brood year in the data set and/or Mmax. But if you're specifically looking for a forecast of either the 2022 fry outmigration or fry per female for the 2021 brood year, it would be necessary to add a row to the data for 2022 and re-fit the model. More importantly, that forecast would be much more precise if the 2021 spawner data (abundance, age and sex composition) were included; currently* those are all missing. Maybe that's not what you're after, though?

For the moment, I can say that the posterior median and 95% credible interval of fry capacity for Hamilton Channel in the Ricker model (in units of fish, after expanding Mmax by the habitat length of 0.4 km) are 866603 (492735, 6327817).

*The good news is that I did eventually get the model to fit the updated data. It still generates a bunch of NaN or Inf values in multiple variables during initialization for reasons I don't understand, but eventually the sampler is able to initialize itself and get moving. I suspect this may not be the last we'll see of this issue, but for now the results look fine.

[Hillsont]

All,

If I'm following this correctly the IPM is expanding 577 females to a potential egg deposition of 2166507. That would equal ~3,750 eggs per female which is way above what the hatchery data indicates, I started at ~3K but dropped to ~2.75K when I was making outmigrant predictions for Hamilton Springs and Duncan channels. I also noticed the point estimate was outside of the CI (1231836, 15819543), mis-click? My back of the napkin math predicts ~793K outmigrants (577 females at 2,750 eggs per female and 50% egg to fry survival)

Thanks

TH

[kalebentley]

Todd - responding to your post quickly. The "2166507" number was specifically in regards to the total capacity of the Hamilton Spring spawning channel and NOT the forecast of 2022 outmigrants based on the 2021 female spawner abundance. That said, it looks like Eric edited his response and the capacity estimate is 866603 (492735, 6327817), Quite the range!

[ebuhle]

@Hillsont, the figures provided above are for maximum recruitment aka fry "capacity" for Hamilton Channel. As noted, I didn't attempt a prediction of the 2022 fry outmigration yet, pending clarification from @kalebentley and/or 2021 spawner data. And yes, I edited that comment to correct a typo as well as a dumb math error in my initial post.

EDIT: @kalebentley beat me to it!

[ebuhle]

Quite the range!

Yeah, many of these parameters that are positive by definition have pretty skewed posteriors on the log scale, which become extremely skewed on the natural scale. FWIW, I've used 90% credible intervals throughout the reporting on this model; the 90% CI in this case is a somewhat more reasonable (531686, 3540780).

[Hillsont]

Sorry, my bad. I have a habit of skimming and jumping to the end, in this case Kale needing to predict outmigrants in 2022 for permitting issues.

TH

[ebuhle]

No worries. I have a habit of commenting and then immediately realizing I need to edit my comment, which isn't much of a problem on GitHub but can be confusing if you're interacting with Issues through email.

[ebuhle]

@kalebentley @Hillsont Following up on our email thread about power analysis for estimating p_HOS, I wrote out the IPM-based estimates (posterior medians and 95% credible intervals, only for those cases with fit_p_HOS == TRUE) and compared them to frequentist ones (MLEs and 95% confidence intervals based on n_W_obs and n_H_obs, for all cases). The point estimates aren't directly comparable, as the MLE would correspond to the posterior mode, not median; however, with a handful of exceptions at the low end, the uncertainty intervals are nearly identical. This supports my claim that you could do the power analysis using just the origin-composition samples and get basically the same answer you'd get from the IPM, but much faster. Also, a Bayesian "power analysis" using the beta posterior given a uniform beta prior will give you the same answers as the frequentist version shown here. (Substitute Dirichlet for beta if you're interested in the multinomial origin-composition rather than simply p_HOS.)

To @kalebentley's question about quantifying the uncertainty in proportions: yeah, I was going to point out that the CV of the components of a Dirichlet (including beta as a special case) is not constant with respect to the point estimates / underlying probabilities. It's not even approximately constant, as in the lognormal. This means that, for a given total sample size, the CVs of different components (e.g., in this example, p_HOS vs p_NOS = 1 - p_HOS) will differ, and there isn't a single summary statistic for the entire vector of proportions. Intuitively, this happens because proportions closer to 0.5 have more wiggle room than ones close to 0 or 1.

Possible solutions off the top of my head:

1. Compute the CV of each component and base your power analysis on the component with the largest CV.
2. Ditto, but use the average CV (this doesn't seem like a good idea IMO).
3. Transform the proportions from the simplex to the real line, and base your power analysis on the SDs in the transformed space. This could be done by (i) taking the logit of each marginal proportion vs. the sum of the others, or (ii) taking the additive log ratios (sometimes called the "multinomial logit") of all components vs. a reference class, with the obvious choice being Natural spawner. This will stabilize the componentwise variances somewhat, but they will still differ from one another, so you'll be back to (1) or (2).

A final note for now: unless you are transforming the Dirichlet posterior to log-ratio space, it's not necessary to actually sample from it. You can just use the analytical formula for the componentwise variance to compute the CV.

[kalebentley]

@kalebentley @Hillsont Following up on our email thread about power analysis for estimating p_HOS, I wrote out the IPM-based estimates (posterior medians and 95% credible intervals, only for those cases with fit_p_HOS == TRUE) and compared them to frequentist ones (MLEs and 95% confidence intervals based on n_W_obs and n_H_obs, for all cases). The point estimates aren't directly comparable, as the MLE would correspond to the posterior mode, not median; however, with a handful of exceptions at the low end, the uncertainty intervals are nearly identical. This supports my claim that you could do the power analysis using just the origin-composition samples and get basically the same answer you'd get from the IPM, but much faster. Also, a Bayesian "power analysis" using the beta posterior given a uniform beta prior will give you the same answers as the frequentist version shown here. (Substitute Dirichlet for beta if you're interested in the multinomial origin-composition rather than simply p_HOS.)

Cool. Appreciate the demonstration.

To @kalebentley's question about quantifying the uncertainty in proportions: yeah, I was going to point out that the CV of the components of a Dirichlet (including beta as a special case) is not constant with respect to the point estimates / underlying probabilities. It's not even approximately constant, as in the lognormal. This means that, for a given total sample size, the CVs of different components (e.g., in this example, p_HOS vs p_NOS = 1 - p_HOS) will differ, and there isn't a single summary statistic for the entire vector of proportions. Intuitively, this happens because proportions closer to 0.5 have more wiggle room than ones close to 0 or 1.

"...proportions closer to 0.5 having more wiggle room..." may be part of it but what I noticed is that the standard deviations of each multinomial component is the same. Therefore, any component with a relatively high mean will have a smaller CV while any component with a low mean will have a larger CV.

Possible solutions off the top of my head:

1. Compute the CV of each component and base your power analysis on the component with the largest CV.

Assuming the mean-level of pHOS is 5%, it will take approx. 875 samples to achieve a CV <15% (see plot below).

For quick comparision, if we run only 400 samples, here's the results where origin 2 would be the equivalent of pHOS (and origin 1 would be pNOS)

2. Ditto, but use the average CV (this doesn't seem like a good idea IMO).

Hmm...ok.

3. Transform the proportions from the simplex to the real line, and base your power analysis on the SDs in the transformed space. This could be done by (i) taking the logit of each marginal proportion vs. the sum of the others, or (ii) taking the additive log ratios (sometimes called the "multinomial logit") of all components vs. a reference class, with the obvious choice being Natural spawner. This will stabilize the componentwise variances somewhat, but they will still differ from one another, so you'll be back to (1) or (2).

I mostly follow you here. If it's worth pursuing this alternative approach, I may have some specific questions on how to actually implement it.

All said - I still wonder if CV is an appropriate measure of precision for pHOS (or any proportion) given what is highlighted above. I scanned through a few NOAA monitoring guidance documents (e.g., Crawford & Rumsey 2011) to see if there was any clarity on this topic but couldn't find anything in the few minutes I looked. Ultimately, I would think we would be interested in the precision of the abundance estimate for NORs and HORs, which together get you estimates of pNOS and pHOS and the precision is a result of both the precision of the (M-R) abundance estimate and the precision of the composition estimate. So setting a precision level on hatchery abundance is still going to be driven by the precision of the composition estimate.

Hmm...getting back to the original question re: "how many samples". We could set our own threshold. For instance, instead of CV we could set a threshold on range or quantile. For example, range of 95% CI is <0.10. Without having a specific application for the precision, any threshold seems a bit artibrary. Any other thoughts?

A final note for now: unless you are transforming the Dirichlet posterior to log-ratio space, it's not necessary to actually sample from it. You can just use the analytical formula for the componentwise variance to compute the CV.

I follow. I realized this after a little bit of reading yesterday but wasn't sure how to calculate specific quantiles (e.g., lower and upper 95% CIs). Generating the samples in R with the rdirichlet function takes a fraction of a second.

[ebuhle]

Thanks for that numerical illustration, @kalebentley.

"...proportions closer to 0.5 having more wiggle room..." may be part of it but what I noticed is that the standard deviations of each multinomial component is the same. Therefore, any component with a relatively high mean will have a smaller CV while any component with a low mean will have a larger CV.

I was being a bit breezy and imprecise there, so let me (hopefully) clarify this. For a Dirichlet with parameter vector $\boldsymbol{\alpha}$, where $\sum{\boldsymbol{\alpha}} = \alpha_0$, the SD of the $i$-th component $p_i$ is

$$ \frac{1}{\alpha_0} \sqrt{ \frac{\alpha_i (\alpha_0 - \alpha_i)}{\alpha_0 + 1} } . $$

In the special case of 2 components, where the Dirichlet collapses to a beta, $\textrm{SD}(p_1) = \textrm{SD}(p_2)$ because $\alpha_i(\alpha_0 - \alpha_i) = \alpha_1 \alpha_2$ and thus the numerator inside the square root is the same. Or to put it much more simply, $\textrm{SD}(p_1) = \textrm{SD}(1 - p_1) = \textrm{SD}(p_2)$. With more than 2 components, this equality doesn't hold unless all the components of $\boldsymbol{\alpha}$ are equal.

But we're really interested in the CV of $p_i$, given by

$$ \sqrt{ \frac{\alpha_0 - \alpha_i}{\alpha_i(\alpha_0 + 1)} } = \sqrt{ \frac{1}{\alpha_0 + 1} \left( \frac{\alpha_0}{\alpha_i} - 1 \right) } = \sqrt{ \frac{1}{\alpha_0 + 1} \left( \frac{1}{\mathbb{E}(p_i)} - 1 \right) }, $$

which is strictly decreasing in the expectation of $p_i$, as you observed. Bigger proportions will have smaller CVs, all else equal.

Given that fact, I agree that using the CV as a measure of uncertainty for a proportion is a bit odd, and not something I've really seen done before as far as I can recall. On the other hand, I'm not aware of a general summary statistic for the precision of a probability vector that would address your question either. Since we're working with the Dirichlet, the canonical measure of precision is simply $\alpha_0$, the total effective sample size, but that's not very helpful for figuring out what the total sample size ought to be! Pragmatically, the CV still seems useful so long as you stipulate that the goal is to manage the precision of the smallest $p_i$. Your idea of using the width of the credible interval could work too, although it might lead to some counterintuitive results as $p_i$ approaches 0 or 1 (which goes back to my point about the absolute variation being greatest at 0.5).

I follow. I realized this after a little bit of reading yesterday but wasn't sure how to calculate specific quantiles (e.g., lower and upper 95% CIs). Generating the samples in R with the rdirichlet function takes a fraction of a second.

Fair enough; I've certainly got nothing against brute-force simulation! But if you ever do want to compute the posterior quantiles directly, you can use the fact that the marginal distribution of $p_i$ is beta:

ci95 <- qbeta(c(0.025, 0.975), alpha[i], sum(alpha) - alpha[i])

[kalebentley]

Thanks @ebuhle. I appreciate the detailed explanation of the Dirichlet distribution. Super helpful.

Once again, returning to Todd's original question & our problem statement: "how many (otolith) samples do we need to run to estimate origin?” with the focus on potentially decreasing the current target number given increased analysis costs this coming year. Currently, we have 900 otolith (Coast strata populations) and 1200 PBT samples (Cascade and Gorge strata populations plus Grays brood to look for strays into these areas) budgeted to run for origin determination, which right now cost $13 per ololith sample and $30 per PBT sample (not sure what the cost increase is going to be).

Focusing on the Coast strata first, if we want to stick with meeting a precision target of a CV <= 15% for pHOS, which is going to be our smallest pi then we'll need to stick with a target sample size of approx. 900 assuming pHOS is going to be approx. 5%. If we are ok with pHOS estimates having a CV of say 20-25% then we could get away with running approx. half the number of samples. I don't see an issue with this decrease in precision given that pHOS is so slow. Similar to what I highlighted in my previous post, assuming the mean pHOS will be 5%, here's the results you get with 450 samples: And here's what you get with 900 samples: Sure - on average, you met your 15% CV goal with 900 samples but the 95% CIs are practically identical to the 450 sample results.

Curious Todd - what is the increase in otolith sample cost this coming year? Running 450 fewer samples assuming the same cost of $13/sample will save us $5,850. While I get that we want to spend every dollar as wisely as possible, this seems like a pretty small amount overall.

I am sure it is obvious to everyone but these results are appliable to a single report grouping (e.g., Coast strata) with only two groups (e.g., pHOS and pNOS). I am not certain how our sample target of 900 has been distributed among the Coast populations but guessing most have been from the Grays. If we want to make population-level estimates of pHOS (e.g., Grays, Elochoman, Skamokawa) and still meet these same precision targets then we'll need to increase our sample sizes.

As for origin calls for the Cascade & Gorge Stratas and corresponding populations, we could take a closer look at the estimated proportions of Duncan channel, Duncan hatchery, & Lewis hatchery origin fish relative to NORs to see if we could reduce the number of PBT samples. However, as you stated before @Hillsont, "we’ve always run more origin samples for the Cascade and Gorge strata populations, assumption has been lower numbers of project origin adults needing higher sample rate to detect." Scanning through the table of pHOS estimates @ebuhle generated earlier this week, it looks to me that the "pHOS" estimates are pretty darn precise for most years and populations by simplying looking at the upper credible interval he report (which I'm guessing is the 90% quantile). Curious, @ebuhle, do these estimates of pHOS lump together all non-NOR samples (i.e., do they include Duncan Channel origin fish or only Duncan and/or Lewis Hatchery fish)? Regardless, since these are "population-level" estimates (as opposed to Strata-level), I assume the precision will only improve as the data are lumped.

Thoughts on next steps @Hillsont? Do you want to "sit down" together and explore various sampling scenarios to evaluate how it'll affect the precision of our estimates vs. the cost-savings?

[ebuhle]

If we are ok with pHOS estimates having a CV of say 20-25% then we could get away with running approx. half the number of samples.

Given those results, I agree.

Scanning through the table of pHOS estimates @ebuhle generated earlier this week, it looks to me that the "pHOS" estimates are pretty darn precise for most years and populations by simplying looking at the upper credible interval he report (which I'm guessing is the 90% quantile).

97.5th quantile, actually; those are 95% credible intervals and frequentist confidence intervals. Sorry, that detail is buried somewhere upthread but I should've put it in the variable names. And yes, the uncertainty around a small proportion has a long tail, but by 100-200 IID samples you're well within the regime of diminishing returns:

Curious, @ebuhle, do these estimates of pHOS lump together all non-NOR samples (i.e., do they include Duncan Channel origin fish or only Duncan and/or Lewis Hatchery fish)?

Yes they do. In this version of the model, "hatchery-origin" really means "known nonlocal origin".

Regardless, since these are "population-level" estimates (as opposed to Strata-level), I assume the precision will only improve as the data are lumped.

Well ... actually, this raises a very important point that I had totally glossed over, namely that all of this probability math relies on the assumption that samples are IID. We make that assumption within a population and year, even though it's violated in myriad ways, because it's the best we can do. But when moving to the Strata level and combining independent samples from different populations, you'll need to use a hierarchical estimator to account for among-population heterogeneity. Otherwise you'll overestimate the precision of the Strata-level mean. That hierarchical structure is inherent in the IPM, but not in the analytical beta posterior.

[kalebentley]

Hey @ebuhle,

I just had a call with @Hillsont and @BradGarnerWDFW on the topic of reporting metrics. In short, there is a report that is due on May 15th specific to the Duncan "program". Below is a mock-up of the data/estimates from the IPM that Brad and Todd would like to use in the report. I provide some notes on each reporting metric of interest. The summary .csv table should include all available years from our data set/analysis, which I believe is 2002 - 2021.

Let's still plan to generate .csv results files for every metric that we outlined in the thread above (by population and year) sometime relatively soon but this specific Duncan summary is the top near-term priority and ideally it could be generated by Wed., May 8th.

Thanks, Eric!

Kale

[ebuhle]

Hi @kalebentley, thanks for that template. Most of these quantities are straightforward to estimate from the existing IPM.

The proportion of green females is considered known data (p_G_obs in fish_data), not an unknown state, so the IPM has nothing to add beyond the observed sample statistic.

We cannot presently estimate % transported from source pop in general; recall that the current model formulation doesn't link broodstock removals to adult transfers into hatcheries or artificial channels, so this will have to wait for completion of the broodstock-to-smolt phase of the life cycle (#24). Regarding Duncan Channel specifically, adults returning to Duncan Creek are considered natural recruitment to Duncan Channel, which is why there's no Duncan Creek in fish_data$pop:

https://github.com/ebuhle/LCRchumIPM/blob/2782a463b174c2e23deb6b6446c95454746773ef/analysis/R/01_LCRchumIPM_data.R#L81-L89

In other words, the % transported from the source Duncan Creek into Duncan Channel is 100% by definition. But that aside, it's not clear to me how the quantity defined in the table (local recruitment to Duncan Creek as a proportion of total adults in Duncan Channel) represents the proportion of adults transported from their source populations (not only Duncan Creek but also Ives, etc.) to Duncan Channel. We can easily estimate the former quantity, but it would be better described as "% local broodstock" or something like that.

Everything else can be easily done, for Duncan Channel and all other populations, and I'll do that in the next few days.

[ebuhle]

@kalebentley @Hillsont @BradGarnerWDFW, I put together a posterior summary table of (almost) all the reporting metrics discussed in this thread that the IPM in its current form is able to estimate. See /analysis/results/reporting_metrics.csv (click in upper right to download).

All quantities are summarized as posterior median (.med) and 90% credible interval (.05, .95). States (e.g. abundance) are indexed by the calendar year in which the life stage was observed, and transition rates by the beginning (e.g. SAR by outmigration year). This coincides with brood year in the case of spawner abundance and composition, but not in general.

The notation follows the model as described in our report:

• S Spawner abundance.
• E_hat Expected egg deposition. Recall that eggs are not an independent state in our model; egg deposition and egg-to-fry survival are both subsumed in the S-R function and a single process error term. As we've done in the past, we can compute expected egg deposition $\widehat{E_t}$ as a deterministic function of spawners, age-weighted fecundity, sex ratio and non-green female discount.
• M Juvenile outmigrant (fry / smolt) abundance.
• M/E_hat Egg-to-fry survival based on expected eggs. Because of process error, this ratio can exceed 1 if juvenile abundance (the realized state) exceeds expected egg deposition.
• q_3, q_4, q_5 Proportions of 3-, 4- and 5-year-old spawners.
• p_HOS Proportion hatchery-origin spawners.
• p_local Proportion of spawners that recruited locally vs. being translocated. Only defined for Duncan Channel, in which case local recruits return to Duncan Creek. (This is my interpretation of @kalebentley's "% transported from source pop"; let me know if I got it right.)
• q_F Proportion female spawners.
• s_MS SAR (survival from smolts M to spawners S).
• R/S Adult recruits per spawner, indexed by spawning year.

Recall that the current model formulation does not have states representing spawner or smolt abundance in hatchery populations, so the only available quantity of interest is SAR. This will change when we complete the hatchery life cycle (#24) with a process model of hatchery spawner abundance and composition and juvenile production.

I didn't summarize the broodstock mining rate because it's not monitored by default, so the samples are not in my saved model object. I would just have to re-fit the model and save B_rate. But anyway the hatchery broodstock-to-smolt piece (#24) will make this more interesting and less of a nuisance parameter.

[kalebentley]

Awesome. I was about to ask if you could elaborate on how M/E_hat could be >1 but followed your link over to issue #6 and realized I need to spend time re-familiarizing myself with this LONG thread/topic. Appreciate your fast turnaround time getting this reporting metrics file generated (and your post above defining each variable is SUPER helpful). @Hillsont, @BradGarnerWDFW, and I will review this information this week and circle back with you, as needed. Thanks @ebuhle.

[kalebentley]

In other words, the % transported from the source Duncan Creek into Duncan Channel is 100% by definition. But that aside, it's not clear to me how the quantity defined in the table (local recruitment to Duncan Creek as a proportion of total adults in Duncan Channel) represents the proportion of adults transported from their source populations (not only Duncan Creek but also Ives, etc.) to Duncan Channel. We can easily estimate the former quantity, but it would be better described as "% local broodstock" or something like that.

I forgot to add - I looked at the "reporting_metrics" file and the value generated for p_local is exactly what I was looking for in my original template (and fwiw, my column header "% transported from source" was poorly worded).

As an aside, having these metrics in one file is awesome! Out of curiosity, I generated the two plots below which show that the total abundance of spawners in Duncan Channel (top) and the percentage of the total spawners that came from the Duncan trap (bottom), as opposed to being translocated from the mainstem or Hamilton, has steadily increased over the past 20 years.

Interestingly, the percentage of adults that have been assigned as a Duncan Channel origin recruit has not increased but rather decreased over the past 15-ish years. I haven't completely thought this one through but my initial hunch is that it is related to the relative size of the Duncan population compared to larger nearby populations and the overall increase in abundance observed for the Lower Gorge strata/MPG in the past 6-8 years (see final plot below). It will be fun to start digging into these metrics and exploring more!

[ebuhle]

Yeah, I was having some PTSD flashbacks re M/E_hat. It occurred to me that an alternative way to represent egg-to-fry survival would be M_hat/E_hat, where $\widehat{M_t}$ is expected fry abundance, a deterministic function of $\widehat{E_t}$ and the S-R parameters. This way, we would ignore the process error in juvenile production altogether (conditional on the states that determine it). The advantage would be that "expected egg-to-fry survival" cannot exceed 1; the disadvantage is that it ignores information in the data indicating that realized fry abundance $M_t$ departs from its expected value.

Re p_local, I also noticed some interesting patterns. In the data frame below, p_local_obs is the observed sample proportion:

              pop year S_obs S_add_obs p_local_obs p_local.med p_local.05 p_local.95
1  Duncan Channel 2002    65        65        0.00        0.12       0.03       0.40
2  Duncan Channel 2003    54        54        0.00        0.15       0.03       0.46
3  Duncan Channel 2004    69        68        0.01        0.08       0.02       0.24
4  Duncan Channel 2005    46        39        0.15        0.18       0.11       0.33
5  Duncan Channel 2006    54        31        0.43        0.42       0.32       0.51
6  Duncan Channel 2007    38        30        0.21        0.22       0.15       0.32
7  Duncan Channel 2008    42        40        0.05        0.14       0.06       0.27
8  Duncan Channel 2009    51        25        0.51        0.51       0.41       0.59
9  Duncan Channel 2010    76        26        0.66        0.66       0.61       0.74
10 Duncan Channel 2011   144        70        0.51        0.52       0.43       0.62
11 Duncan Channel 2012    50        46        0.08        0.30       0.15       0.49
12 Duncan Channel 2013    88        61        0.31        0.31       0.21       0.43
13 Duncan Channel 2014    96        72        0.25        0.35       0.24       0.55
14 Duncan Channel 2015   182        31        0.83        0.83       0.79       0.86
15 Duncan Channel 2016   268        62        0.77        0.75       0.58       0.79
16 Duncan Channel 2017    75        68        0.09        0.27       0.13       0.45
17 Duncan Channel 2018   192        63        0.67        0.67       0.55       0.72
18 Duncan Channel 2019   106        30        0.72        0.72       0.67       0.78
19 Duncan Channel 2020   123        45        0.63        0.65       0.59       0.75
20 Duncan Channel 2021   530        42        0.92        0.92       0.89       0.93
21 Duncan Channel 2022    NA         0          NA        1.00       1.00       1.00

The sample- and IPM-based point estimates are mostly very close, and the 90% CIs almost always contain the former. The main exceptions are related to boundary issues (the IPM process model will always produce some local recruitment, so p_local can never be identically 0), but especially so in 2002-2004. Looking at spawner abundance in those years, it's evident that the estimated states are biased slightly high relative to the point observations. I think this is driven by shrinkage toward the common ESU trend: Duncan Channel doesn't show the same pattern of starting high, then rapidly descending seen in most other pops. But it also reflects the awkward situation of needing to impute the "unknown" $\tau_S$ in Duncan Channel, which results in unrealistically large wiggle room for the states relative to the observations, as I pointed out in the context of hatchery / channel observation error.