Determining an optimum and distribution range for a given species abundance vs an environmental gradient data

Question

This is a fairly complex question so I will attempt to ask it in a fairly basic manner.

I have data on the abundance of 99 different species of estuarine macroinvertebrate species and the sediment mud content (0 - 100 %) in which each observation was obtained. I have a total of 1402 observations for each species (i.e. a massive dataset).

Here is a subset of the raw data for one species to give you an idea of the data I'm working with (if I had 10 reputation points I'd upload a plot of real raw data):

Abundance: 10,14,10,3,3,3,3,4,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,6,6,6,0,0,0,0,12,0,0,0,34,0,0
Mud %:     0.9,4,2,10,13,14,6,5,5,7,22,27,34,37,47,58,54,70,54,80,90,65,56,7,8,34,67,54,32,1,57,45,49,4,78,65,45,35

The primary aim of my research is to determine an "optimum mud % range" (e.g. 15 - 45 %) and "distribution mud % range" (e.g. 0 - 80 %) for each of the 99 invertebrate species.

As you can see the abundance data for the above species contains a significant number of zero values. Although this significantly skews any sort of model that I run on the data (i.e. GLM, GAM), even if I model the non-zero data only, the model for certain species does not fit the data at all well.

So, my question is: what would be the best, most robust way to determine an "optimum" and "distribution" mud range for each species, given that responses vary significantly between species?

---

Just to clarify - the above data is a hypothetical example to give you an idea of how messy the abundance (that is count) data can be for a given species.

Regarding the poisson regression approach: I'm considering conducting a two-step GLM or GAM approach for each species; Step (1) uses logistic regression to model the "probability of presence" for a given species over the mud gradient - using presence/absence data. This obviously takes into account the zero counts; and Step (2) models the "maximum abundance" over the mud gradient - using presence only count data. This step gives me an idea of the species typical response to mud where they DO occur. What are your thoughts on this approach?

I have R code for both steps for one particular species. Heres the code:

     ## BINARY

aa1<-glm(Bin~Mud,dist=binomial,data=Antho)
xmin <- ceiling(min(Antho$Mud))
    xmax <- floor(max(Antho$Mud))
Mudnew <- seq(from=xmin, to=xmax, by=0.1)
pred.dat <- data.frame(Mudnew)
names(pred.dat) <- "Mud"
pred.aa1 <- data.frame(predict.glm(aa1, pred.dat, se.fit=TRUE, type="response"))
pred.aa1.comb <- data.frame(pred.dat, pred.aa1)
names(pred.aa1.comb)
plot(fit ~ Mud, data=pred.aa1.comb, type="l", lwd=2, col=1, ylab="Probability of presence", xlab="Mud content (%)", ylim=c(0,1))


## Maximum abundance

 aa2<-glm(Maxabund~Mud,family=Gamma,data=antho)
 xmin <- ceiling(min(antho$Mud))
     xmax <- floor(max(antho$Mud))
 Mudnew <- seq(from=xmin, to=xmax, by=0.1)
 pred.dat <- data.frame(Mudnew)
 names(pred.dat) <- "Mud"
 pred.aa2 <- data.frame(predict.glm(aa2, pred.dat, se.fit=TRUE, type="response"))
 pred.aa2.comb <- data.frame(pred.dat, pred.aa2)
 names(pred.aa2.comb)
 plot(fit ~ Mud, data=pred.aa2.comb, type="l", lwd=2, col=1, ylab="Maximum abundance per 0.0132 m2", xlab="Mud content (%)")
 AIC(aa2)

My question is: for step (2); will the model code need to be altered (i.e. family=) depending on the shape of each species abundance data, if so, would I just need to inspect a scatter plot of the raw presence only abundance data to confirm the use of a certain function? and how would the code be written for a certain species exhibiting a certain response/functional form?

Answer 1

This seems as much ecological as statistical, not that that's a problem: it's just that you need people with both kinds of expertise to answer fully.

From my limited experience with similar data, I'd regard it as typical that many species are just awkward for your purpose. At each site, many species will be elusive as far as your sample of material is concerned, so you will get many zeros. I don't think you can guarantee a clear result for the optimum from each species. The sample data you posted (38 mud and 37 abundances, by my count) certainly look messy, but take it to be nature's fault, not yours.

Much depends on what your "abundance" variable is precisely. If it's a count, then I would do a Poisson regression for each species on mud and mud$^2$, which guarantees (a) fitting a curve with a turning point (b) fitting as it were on a logarithmic scale (c) respect for zeros too. That's only a start but elementary calculus can be brought in to check that the turning point is a maximum and where any maxima occur (within the range of the data or not?). (I take the zeros to be part of the information.)

If your abundance variable is something else, you need something else, but the basic approach carries over.

I don't think you could guarantee that every species will converge to a solution for a Poisson regression.

Is "distribution range" precisely defined? The simplest definition would be the minimum and maximum mud % for which species had been observed.

Some small points on terminology:

• While 1402 observations is undoubtedly a lot of work if someone is sieving sediment and checking for beasts it won't qualify as "massive" by many standards here.

• The word "significant" in statistics usually implies a significance test. It's not clear what you mean by "a significant number of zero values" or "this significantly skews". Words like "large" and "greatly" are fine for informal use and avoid the ambiguity here.

• "Robust" is too open-ended to create a precise impression. Again there are technical and non-technical senses; if you have a particular technical sense in mind, it's best to spell it out.