I’ve fit a nonlinear model to data that measures how several species are affected by competition with one another. The model converges, and the diagnostics— \hat{R}, ESS, posterior predictive check—all look good. I now want to calculate two *quantities of interest* for each pair of species using the parameters estimated by the model. Say, for the sake of simplicity, that my model has two parameters \alpha_j and \beta_j estimated for each species, j. The *quantities* are calculated on each species pair j,k using simple arithmetic operations such as, \beta_j -\sqrt{\frac{\beta_k-\alpha_k}{\alpha_j}}.

The obvious first pass is to use the estimated population level effects (i.e., the central tendency of each parameter’s posterior distribution), but of course the value of having an estimate of the full posterior distribution is that I can propagate the uncertainty associated with each of the parameters to understand its impact on the *quantities of interest*.

I understand that Stan’s posterior samples are autocorrelated, and I’m concerned about how autocorrelation in the posterior samples might bias the calculation of the *quantities of interest*. Suppose I want to calculate these *quantities* using a subset, i, of the total number of posterior samples, N.

I can imagine several different approaches for selecting this subset, i:

- Use
*all*posterior samples (i = N). The disadvantage here is that any visualization of the*quantities of interest*would have substantial overplotting, compounded by the need to expand the axis ranges to encompass rare values far from the central tendency. The upshot is a figure that’s difficult to parse. - Use sequential samples, N_1,N_2,...,N_i. Depending on how many post-warm-up iterations were done on each chain, it would be possible that all samples come from the same chain.
- Draw i samples from N with equal probability across all chains.
- Stratify sampling by drawing every n-th sample. Effectively this would be a thinning approach. n could be chosen to space samples as far apart from one another as possible, or could be drawn from a uniform random distribution with lower & upper bounds chosen to prevent sequential or near sequential samples from being drawn.

**From a conceptual perspective, is there any reason to suspect that these different approaches to subsetting the posterior draws would give different answers when the draws were used in further calculations?** If yes, which (if any) would be the most robust strategy? Is there a principled way to use model diagnostics (e.g., ESS vs. N) to guide the decision?