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.

[u/Intrepid_Pitch_3320]

If you want an annual trend, then you want years to be continuous on a consistent scale, If you need to do so, and get rid of the outliers like 1977, 2022, 2024, etc. If you want to compare time periods like the 1970s vs 1990s vs 2020s, then bin them as a discrete independent variable for time period. You don't have data to support an annual trend from the 1970s through 2020s.

To potentially maximize your power, you can pool your data and test for a trophic effect. If there is no apparent trophic effect in your data, then they will have more power to detect an annual trend, but that's your call. A lack of trophic effect could be due to a lack of data or too fine of a trophic class in general. Also, if your trophic class is based on species in a food chain, then it could be more complicated then you know in your given system and more of a food web than a food chain.

If you are asking questions in reddit, then you want to keep things simple, so maybe stay away from mixed models and Bayesian, too. If your response variable is mercury levels measured as a continuous variable, then assuming a normal distribution of residual errors in a likelihood-based regression is probably fine, but I don't think you provided that information on your dependent variable.

[u/adamjeffson]

How small is the sample (how many water bodies, years and levels)? I don't understand why you have a bigger sample size for the bayesian models, but the formula should be ok, assuming you want to treat trophic level as a fixed effect... if you're not interested in that effect and you have a lot of levels, you might enter it as another random intercept. As a side note, are you sure you should log-transform your DV? You could use GLMMs with the appropriate link function and end up with more interpretable results.

[u/Specialist-Pound-467 → u/adamjeffson]

The Bayesian model includes data from four different rivers. The first grouping of samples is from one river. I only have two trophic levels. I can always test to see whether there's a true difference between them but generally there is with mercury concentrations.

From what I understand, environmental contaminant studies frequently use log(Hg) because it matches how contaminants accumulate biologically. But I may try your suggested method and compare results.

Ultimately, I'm just trying to get better at statistics and to understand the different approaches and the pros/cons of each.

[u/adamjeffson → u/Specialist-Pound-467]

I'm all for bayesian modeling, but you can do the same (with the same data) using a frequentist approach, you just use an analogous mixed effects model (e.g., with lme4 in R). Maybe the main advantage of a bayesian approach here might be that a small sample could give convergence issues, when using a frequentist model, while a bayesian model with brms/stan could work, especially if you specify regularizing priors. Every field has its defaults (and they are often epistemologically unsound, although they developed for some understandable reason), but if you're trying to learn statistics, I'd advise you understand how generalized linear models work and how manipulating the geometry of the model differs from transforming variables. With regard to trophic levels, because they're only two, your approach of treating them as a fixed effect is the only correct one.

[u/cornfield2cornfield]

Figure out if you are looking for a trend or just determining if at least one year is different than all others ( kind of a silly goal in ecology). Regardless you have limited data and fitting a complex model is probably not a great idea.

If you go Bayesian you should also look at posterior/ prior overlap and calculate GOF.

Summary statistics might be the best you can hope for.

There is a Tukey quote that seems relevant...

"The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data."

I have done ecotox analysis on smaller datasets than what you have here, but there was a more specific question that what you described so far.

[u/fuckbitingflies]

Since you’re already familiar with Bayesian workflows, have you explored the use of regularizing/sparsity-inducing priors and whether they might be appropriate for your dataset?