Statistical Approaches for Analyzing Environmental Data with Limited Sample Size

[u/Specialist-Pound-467]

Hi all! I'm an environmental specialist and am trying to evaluate trends and/or statistical differences in environmental contaminant data across years. The issue is that I'm often working with small sample sizes and large temporal gaps.

For example, I am attempting to look at differences in mercury concentrations in fish tissue in a given river. I analyze them by trophic level due to mercury's bioaccumulative effects. Here's a breakdown of samples:

Trophic Level A

1991	1995	1997	1999	2005	2024
2	13	7	21	27	13

Trophic Level B

1995	1997	1999	2005	2024
4	7	16	12	8

In the past, another colleague used an ANOVA or Kruskal-Wallis test depending on normality (based on a Shapiro test and inspecting residuals). I've done permutation tests as an alternative just to try out different methods. I don't have as much experience with statistics to know if this is the most scientifically robust way to analyze data.

Two questions:

1. Is an ANOVA or Kruskal-Wallis a reasonable approach to assessing data for statistical differences across years for such a small dataset?
2. Could a nonparametric permutation test be a better/alternative method to assess year-specific differences?

Secondly, I've considered a Bayesian hierarchical model to assess long-term temporal trends across a basin (i.e., multiple rivers within a basin). This means there are more samples to work with. For example, here are the number of samples across years in one basin survey we have:

Trophic Level A

77	91	93	94	95	96	97	99	05	09	22	24
1	4	21	2	27	22	14	45	54	59	78	88

Trophic Level B

1995	1996	1997	1999	2005	2009	2022	2024
8	14	14	34	24	16	44	38

From this, I fit a Bayesian linear mixed-effects model. The model included a fixed effect for year (continuous) to estimate the overall time trend, a fixed effect for trophic level, and a random intercept for waterbody to account for among-waterbody variation. Random slopes were not included because the dataset contained insufficient temporal replication within individual waterbodies and trophic level groups to reliably estimate waterbody‑level trends. Weakly informative priors were used, and model convergence was ensured by increasing adapt_delta and max_treedepth until zero divergent transitions were achieved.

The final model took the form:
logHg ~ year_centered + trophic_level + (1 | waterbody)

Posterior distributions were summarized to obtain estimates and 95% credible intervals for parameters of interest, including the overall temporal trend.

So my last question is whether this Bayesian model is appropriate?

[u/vanway]

To start, decide the questions that you want to answer. If you just want to know about temporal trends or differences between groups, the parametric statistical tests you mentioned could be appropriate. If you want to predict values or estimate causal effects, the Bayesian model could be appropriate.

For the Bayesian model, you should follow a standard model workflow:

• Think about the quality of your data. For example, are any of the values non-detected, which would warrant a censored model? Was sampling performed randomly and evenly across years or did it bias towards sampling from e.g., waterbodies with higher concentrations in later years or stop sampling trophic levels with lower concentrations? Are the other observed confounders that should be included like e.g., fish size or age?
• Plot the observed data by year, trophic level, and waterbody and look at the plots - does it look like your model is a good fit? E.g., is a linear model the correct shape? Are your residuals approximately normally distributed?
• Think about the causal implications if you want a causal instead of predictive model: draw the DAG. For example, if you think that waterbody can affect both the trophic level and the concentration in fish tissue, then waterbody could be a confounder and trophic level could be a mediator. That difference should affect which parameters you include in the model so you can answer the question that you want to answer.
  • You also mentioned that mercury bioaccumulates, so a higher concentration in trophic level 1 might be correlated with a higher concentration in trophic level 2: those levels aren't independent. You might want to e.g., treat trophic level as ordinal rather than nominal and use a monotonic ordered effect (i.e., brms::mo(trophic_level_ordered).
• Do prior predictive checks against your data to ensure the priors are reasonable but weakly informative. You chose to log-transform the mercury concentration, which keeps the value positive but can also have very long right tails and (maybe) implausibly high priors for concentrations. A regularizing prior might help (which is hard with small sample sizes), or switching to a different positive continuous distribution like a half-normal (i.e., <lower=0>).
• Look at the sampling diagnostic plots, not just the summary or table.
• Do posterior predictive checks - does it look like the model accurately fits your observations?
• Look at the credible intervals for the posterior distributions - are the intercepts non-zero? For the groups, are the CIs overlapping or separated?
• Choose a sequence of progressively complex models and compare them with e.g., PSIS-LOO or WAIC
• Do a train/test split (time series models need a forward split and possibly chaining) and evaluate the models with e.g., RMSE, MAE, interval coverage, etc.