Proper analysis for vernal pool study data

Question

I am trying to analyze a dataset consisting of counts of amphibian egg masses (3 different species) in nine vernal pools over a four-year period (consecutive years, 2014 to 2017). During each count (one per year), data were also collected on pool area, pool depth, pH, and conductivity. Initially, I did a standard ANOVA to determine whether pool (9 of them), year, area, depth, pH or conductivity were predictors of egg mass numbers in the pools. I was told (by journal manuscript reviewers) that ANOVA was not an appropriate analysis (because of the single value for each combination of factors, which I agree with), and they suggested two alternatives:

1. Multiple linear regression or similar model that allows for continuous and categorical variables (new variables could be added, such as percent surrounding forest, distance to nearest pool, any surrounding development, or low/high development pressure), both to be included as opposed to univariate correlations that can mask complex relationships among independent variables.
2. Because the same ponds were sampled repeatedly, samples from the same pond in multiple years are not independent replicates. Recommend a repeated measure ANOVA, averaging all the years for each pool and just analyze the data for pool averages.

Please note that sample sizes are small (9 pools, 4 years), so any analysis may not have much statistical power. But even a finding of no effect would have meaning for the study. Any thoughts on appropriate tests/approaches (these suggested ones or others) would be appreciated.

---

Additional Update:

Thanks for the responses. The Bayesian approach is probably the proper one, but I don't think I'm willing to put in the background work to understand/apply it, at least for this project. I'll likely take a crack at the multilevel (mixed?) model, with one IV (independent variable?) for each model, as Peter Flom suggests. I assume that this means running multiple separate models? Some clarification (and maybe appropriate routine/package in R) would be helpful.

The "hypothesis" here is that at least one of the independent variables affects total egg mass counts for a pool. From some preliminary analyses, I already know there is a significant correlation (Spearman) of pool size (area) with counts (the bigger the pool, the more egg masses), so that's one. I'm just looking to verify that with a model (or models, multilevel/hierarchical or otherwise) with which I can also analyze other variables. Perhaps a single model is not appropriate for these data, as mkt points out.

Basic visualization of the data may be all that's needed, but I've found that most journals these days find that insufficient and want hard stats to back it up.

Answer 1

Another approach is Bayesian canonical polyadic decomposition (BCPD), which is generally considered to be an exploratory method. You would start with the notion that your data array $\mathcal{A}$ can be decomposed as a linear combination of tensor products:

$\mathcal{A} = \sum_{k=1}^r \lambda_k \bigotimes_{j=1}^h \vec X_{k,j}$

See Tamara G. Kolda: "Tensor Decomposition" for a (non-Bayesian) explanation of canonical polyadic decomposition and some of the visualization techniques. For a Bayesian version of those decomposition plots you can either plot the MAP values or make a contour or surface plot of the posterior. For a more detailed review of tensor decomposition, I suggest the landmark paper Kolda & Bader 2009.

The gist of CPD is signal separation according to different combinations of levels of factors. Why I like BCPD is for encoding background information and imputing missing values.

Offhand I don't know of a package that does this, especially with customizable priors. But it can be readily built using probabilistic programming packages (PyMC, Stan, etc). At some point I should package my own code to share.

Answer 2

Your journal reviewers seem ... confused.

Multiple linear regression will not solve the problem of having a single value for each set of factors. In fact, linear regression and ANOVA are mathematically the same. In matrix form, both are given by $Y = XB + e$.

If you average all the years for each pool, then you no longer have repeated measures, and this lowers your N to 9 (the number of pools). Either you do RM-ANOVA or you average for each pool, but you can't do both.

In addition, RM-ANOVA makes unrealistic assumptions, including sphericity. A multilevel model would be better. Also, you have five independent variables. That's too many for the amount of data.

Finally, your last paragraph is troubling. If you have low power, then a finding of a non-significant effect is not very meaningful, because only a very large effect would be significant. HOWEVER, if you find small effects, then that may be meaningful. Small sample size does not bias the estimates, it just makes them imprecise.

I think what I would do in your situation is multilevel models with one IV for each model. But you may need to hire a statistical consultant to help you with this. (Not me, I'm retired).

Answer 3

While I 100% agree with Peter Flom's answer (+1), I personally think you should drop the null hypothesis testing approach for this project. This advice is notwithstanding what you can convince your editor is reasonable.

Instead, I would develop a Bayesian mixed effects model. With a Bayesian analysis you will/should spend substantial amounts of time considering your choice of prior and likelihood, but there are no hard constraints on sample size per se. However, with small sample size you should expect that with weak priors you will have large standard errors. That's fine provided that you are okay with being imprecise if that is all your background knowledge (encoded into priors) and the data can support.

If you have insufficient clarity about what I just suggested due to Bayesian inference being unfamiliar to you, I recommend watching Richard McElreath's Statistical Rethinking lectures and then reading Gelman et al 2020.

---

_{^{(I am not retired and sometimes take small modelling projects on the side. I can be found if you look hard enough.)}}

Answer 4

The multilevel/hierarchical models suggested by Peter Flom and Galen are the best solution, in the technical sense. But given your small sample size (9 points, 4 years, 6-10 predictors), and your apparent lack of a specific goal (such as a parameter to estimate or a hypothesis to evaluate), I'd consider instead a simple descriptive/exploratory analysis.

It's not obvious that the hierarchical model will tell you a whole lot more than some basic visualisation. And the hierarchical models will be much more complex to execute and to interpret if you and your audience are unfamiliar with them.

UPDATE based on the addition to the question:

About journals wanting 'hard stats to back it up': you may be right, but statistics doesn't solve the basic problem here. Any model you could feasibly come up with is going to be a vast oversimplification of the true causal network. Among other things, there's plenty of important unmeasured variables - predation pressure perhaps being the most important. And the weak observational design will not allow a model to provide strong evidence for any causal claims. Some evidence of association based on a model may be useful, but there aren't going to be any 'hard stats' you can point to with confidence. This is why I think a descriptive analysis is a good option - it's simple and the limitations aren't buried under model complexity.

Answer 5

Not my area of expertise (limnology, amphibians, or statistical models), but I seem to recall a friend having a similar problem (except with kidney transplant survival data) and he basically fed his data into a few different machine learning frameworks and then pried the lid off the top of the black box and was able to pick out predictors and levels of significance. Depending on the ML framework, e.g. a Bayesian classifier, might have similar results. I will leave it to others to comment on the post-hoc extraction of measures of association and prediction from ML libraries.

I think there is another, probably totally off-topic for a stats forum, perspective that I am surprised the reviewer didn't mention.

Others commented on the problems w/ having a data set of non-independent (see below) measurements w/ only 9 sites × 4 years, where there are likely significant unmeasured differences, e.g. temperature in the weeks leading up to the sample time, rainfall--and timing is key here, changes in predator species (which may be highly localized), spread of pathogens, or even a change in what brand of pesticide the local farmer is using, etc. would all have potential to confound the crap out of the data.

So I think the first thing would be to back up a bit and just look at the data and see if there is a story to tell. You didn't mention what the data actually showed.

I don't think you can use year as a separate variable either, as I don't recall frogs living all that long...so while it is a separate sample w/ a different set of variables, it is only a rough proxy of comparable conditions between your sites, assuming they are in the same general area and elevation. If they are ponds up in the mountains at different elevations, then they very well could have markedly different temperature and precipitation profiles.

Most ponds in areas w/ sedimentary rock underlying a decent layer of soil are roughly a spherical cap (or a spherical segment if the sediment has truncated the water column to a frustum) has a volume $\mathrm V = \frac{1}{6}\pi h (3a^2 + h^2)$. a=area, h=depth

Since both pH and conductivity (which you didn't mention the range) depend on the volume of water in the pond (as well as runoff) these variables are inherently linked--not just by chance, but intrinsically so.

Finally, you didn't mention the year-to-year, pond-to-pond variability. Obviously, your stochastic noise level is pretty high to begin with. One thing you didn't mention, but what you can actually get your hands on is the actual weather data which you might want to try and figure into your analysis as the temperature, rainfall, and to a lesser extent, wind and cloud cover, will impact the water balance which frogs are rather particular about.

If you know the location of the ponds, you can probably pull up enough GIS data to get an idea of land use--but it won't tell you how much organochlorine and other hormonal disruptive pesticides were used year to year (although it should be less in each, we hope).

I know that isn't the answer to the "what statistical model do I use" question, but might help you answer the "how do I make sense of this data" question.