Testing informative priors in lieu of age-frequency data

[ebuhle]

@janohlberger, I did some simulation-testing to explore the prospects for using priors on hyperparameters to compensate for sparse or absent age-frequency observations. The basic setup is nearly identical to the single-population example in the iteroparity.R vignette code: a 30-year time series, maiden ages 3-7, Beverton-Holt dynamics, and plausible parameters for WA coastal steelhead. When there are no missing values in S_obs and n_age_obs has an annual sample size of 10-100, this is what you get with default priors:

Note that rho_p here is one arbitrarily selected off-diagonal element of the 4x4 interannual age-structure correlation matrix R_p; as is often the case, there seems to be very little information on these parameters in the data. The estimated states fitted to observed spawner age frequencies by return year look like this:

Next I did the same thing, but with all n_age_obs set to zero. I don't know if I really want such an extreme case to be the example in the vignette, but it's relevant to the WA steelhead problem and helps to bracket the utility of hyperpriors as a substitute for data. These are the results with default priors:

Clearly the age- and kelt survival-related parameters are totally noninformed by data. The spawner-recruit parameters have shifted toward a more linear fit, with Rmax biased high and constrained on the right only by the prior, and alpha biased low but with a long right tail corresponding to a small probability of hockey stick-like B-H fits. That said, in this case the default priors do a pretty good job of keeping the S-R parameters in a reasonable range. To see how much better we might do with informative priors, I set the Dirichlet prior on mu_p and the beta prior on mu_SS, respectively, to have the true mean and a prior sample size of 100 -- if not technically a "best-case" scenario, then certainly much better than one would be likely to do in practice. These are the results:

Counterintuitively, constraining the mean maiden age structure and kelt survival in a small region around their true values actually makes the alpha and Rmax estimates more biased, not less! The fitted S-R curve has moved further toward an almost-linear relationship with a low slope and arbitrarily high asymptote. At this point you might wonder whether the vague default priors on sigma_p and sigma_SS, which allow far greater interannual fluctuations than the true values, are part of the problem. To test this, I once again went with a maximalist solution to bracket the outcomes: I temporarily edited and recompiled the Stan model with no interannual variation in these parameters, so all p are set to mu_p and all s_SS to mu_SS. (Ignore the posteriors for the hyper-SDs in this figure; they're just sampled from their priors and not used in the process model.)

The spawner age structure estimates are more precise, because now the fluctuations are driven solely by cohort strength. (I realized belatedly it would've been helpful to also show the true states in these age time series plots, so I don't know whether the estimates are more accurate.) As far as the key parameters are concerned, though, setting the age-related hyper-variances to zero only makes things worse; the bias in alpha and Rmax has slightly increased.

In practice, I imagine users will seldom have independent, empirical estimates of the interannual SDs of log-ratio cohort age structure or logit kelt survival. More likely, they would use priors on these hyperparameters to subjectively fine-tune the results so there's not "too much" fluctuation. The limiting version of this would be to make the corresponding states time-invariant, as I've done here. At least in this example, setting sigma_p and sigma_SS to zero, which also happens to be much closer to their true values than a random draw from their default priors, did not help to reduce the bias in the key management parameters. I would expect these biases to be even greater, and more robust to the priors on age-related hyperparameters, in the case where F_rate is low and thus the time series does not explore values of S near the origin (cf. the first example in the user-priors.R vignette code). In either case, it seems like a more direct and more effective solution would be to put informative priors on alpha (e.g., based on meta-analysis) and/or Rmax (e.g., based on habitat suitability models).

Also worth noting, Zabel and Levin (2002) showed that time-invariant imputation of age structure leads to biased estimates of S-R parameters. In the context of their run-reconstruction regression approach, those biases overestimated the strength of density dependence, similar to other kinds of observation errors. The specific effect might be different in the state-space context, but I would still expect to see some bias when "smoothly constraining" the age-structure states toward time-invariance via the prior. The last set of results above is consistent with that idea.

An alternative approach would be to model the data-poor populations hierarchically along with the data-richer ones. I would expect to get better results that way than by constraining the hyperparameters in a single-population model, primarily because accounting for the shared process errors in recruitment and kelt survival -- i.e., realized states, which are not amenable to external priors -- typically leads to tighter estimates of the population-specific parameters. The multi-population example in iteroparity.R shows this in a simulated data set with randomly missing S_obs and n_age_obs. I could easily tweak that example so that one or more populations have no age-frequency data at all. Maybe I'll try that next.

[janohlberger]

Thanks for running these @ebuhle

Looks like more informative age priors will only help to some extent and that constraining sigma_p won't necessarily be a good idea. Forcing time-invariant age structure isn't a 'solution' either, as you pointed out, so in cases where it is practical, I'd also favor the hierarchical approach./ As discussed before, there are limitations to that though - thinking about the management/policy process.

I ran some tests with age structure priors by thinning the Quillayute data. Compared to the default Dirichlet prior, using an informative prior on age comp didn't actually change the SR parameters much in that case. Here are a couple of plots for spawners and recruits as well as estimated age structure using the default Dirichlet prior:

IPM-spawners-and-recruits.pdf
IPM-age-proportions.pdf

Here are the same plots when using an informative prior on age structure that roughly corresponds to the estimated age structure when using the full age dataset:

IPM-spawners-and-recruits.pdf
IPM-age-proportions.pdf

[ebuhle]

Thanks for that example @janohlberger, very interesting. How "thin" are the age-frequency data in this case? Looks like there might be no observations at all, as in my simulated data? Regardless, there seems to be some consistency in our results insofar as informative hyperpriors on mean maiden age-at-return (and kelt survival) in a single-population model don't drastically change the estimates of alpha and Rmax, and any such changes don't necessarily reduce bias when the true values are known.

Here's a demo of the spatially hierarchical approach with 8 populations, 20-year time series, and the same life history and similar parameter values to my previous example. All priors are left at their defaults. In the base case, 10% of S_obs and an independent 10% of n_age_obs are missing at random. Here are the posteriors of the population-level S-R parameters when populations are modeled separately vs. hierarchically:

In this particular random realization of the data, the true hyper-means mu_alpha and mu_Rmax are actually in the tails of their respective posteriors, and you can see that reflected in some of the partially pooled estimates as well. Even so, the latter are still quite accurate in almost all cases, and are so much more precise than the independent estimates that RMSE would undoubtedly be lower for the hierarchical model.

Now here's the same thing, but population H has no age-frequency observations. As expected, its unpooled estimates are markedly worse. The partially pooled estimates are a little worse too, but they are much more robust:

Finally, same thing but with a higher proportion of data-poor populations: now E-H have no age-frequency observations. It's kind of amazing how well the partially pooled estimates hold up. In particular, the estimates for the relatively data-rich populations are barely degraded at all by being combined with an equal number of data-poor ones, and notably they are still much better (in terms of bias + variance) than when the data-rich populations are modeled independently. Stein's paradox, and all that.

I wonder if a formalized simulation study along the lines of the examples in this thread would have any influence on the management / policy considerations. If the options for a data-poor population are (1) no external information, poor estimates; (2) information from distant populations (e.g., literature / meta-analysis) incorporated via the prior, likely still poor estimates; or (3) information from nearby closely related populations incorporated via data, likely improved estimates, perhaps managers would feel differently? Anyway, we can discuss this tomorrow.

[janohlberger]

That's great @ebuhle! Nice illustration of the advantage of the partially pooled model and the hierarchical approach more generally. Though I think the main issue from a management/policy perspective is skepticism about the underlying data (in other systems), such that we'd probably want to incorporate varying levels of observation error, and possibly bias, in spawner abundance estimates. I like the idea of developing this into a more formalized simulation approach, and to have that available when interacting at that science-policy interface. Let's discuss later today!