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?