Nice data - ugly residuals - too many puzzle pieces

Question

Dear statistics fans,

I am absolutely clueless in an analysis and hope you can help because I have nobody else to ask right now and got totally lost in four statistics books and a zillion tutorials over the past five weeks.

And here is an overview of what I am trying to do (also see toy graph for visual aid):

• I try modelling a continuous response variable as function of a whole bunch of predictors along a land use gradient (ie a succession of typical anthropogenic land use classes like agriculture, meadow... to unused, protected landscape)
• I got observation plots in five land uses (A-E in the graph) along that land use gradient, replicated in two different ecosystems; the ecosystems got equal numbers of observations, only some land uses are underrepresented/unbalanced
• there is some spatial autocorrelation as within ecosystems and land uses the oberservation plots are clumped in the landscaped
• the two ecosystems generally have a different baseline level in the response var (orange vs darkred solid line), but the overall pattern is the same; only when plotting histograms of the response vars their distribution always looks bimodal (red solid line in small inset)
• I would like to include the two ecosystems as a fixed factor (same for the five land uses) as I am interested in their difference, not seeking generalization across these
• the predictors along the land use gradient are of two types: i) the ones that follow the same trend all across the gradient, some increasing, others decreasing (blue dashed lines), and ii) the ones that only have a non-zero trend across half of the gradient and are often zeros or close to zero in the other half of the gradient (green dotted lines)

What has happened so far:

• a predecessor of mine tried modelling the whole gradient with linear models, but got no significant results because the response var is decreasing to either end
• hence I set out and separated the dataset in two halves (using the land use ‘C’ in both sets) and tried modelling this with linear and linear mixed models; it is a pity though, as that separation also decreases the number of observations drastically for each subset
• also, this complicates matters in so far, that in order to interprete the gradient as a whole I was advised to use the same set of predictors in both models, even so some of them got a lot of zeros or close to zero in 2 out of 3 land use locations of a half gradient
• anyways, I tried starting with data visualizations and plotting, did anovas between land uses that were promising
• followed by extensive correlation checks and PCAs of all predictors, compiled some nice predictor sets that have no correlation greater |0.75| (although I had to drop some nice predictor vars in the process, unfortunately)
• I started out with lm() and stepAIC(), then added interactions for a few hand-picked ecologically meaningful term combinations
• then moved over to lmer() to include the nestedness of land use within ecosystem, but then I got confused if that is correct, as that means both factors are introduced in the random part of the model, but I wanted to keep them as fixed factors ... ?
• then I tried lme() ans stepAIC() and also doodled around with a variance structure to allow for different variances per ecosystem and land use but I ain’t sure what I am doing there and if I do it correctly, as the same factors I used the variance structure for are also my main fixed effects of interest
• then I took the book by Zuur & Ieno (Mixed effects models and extensions in ecology with R, 2009) an tried to work through their stepwise model construction process with gls() and lme() but checking residual plots of the models is not really fun, I believe partly because of the zeros in predictors, and some oddballs like in one particular ecosystem-land use combination having all residuals turning to zero variance
• I then tried standardizing and transforming the hell out of predictors and even response variables, but nothing really helps much, plus I bet you people know a way nicer method than that

My Goal:

• get together a data analysis in R that adequately deals with the data I have (n<100 observations), preferably picturing the whole gradient as it is an ecologically meaningful real-life gradient
• find a proper way of identifying the most important predictors and their effect sizes on my response variable in such a way that the outcomes are still interpreteable in ecologically meaningful terms (e.g. no four-way-interactions please etc)
• I need to get this done with my limited understanding of statistical mathematics, yet I am eager to learn and understand new things (just not all at once)

I would be really grateful if one of you could hint me to the right kind of analysis, as right now I ain’t even sure what I should look for and in which order I should implement a zillion bits and pieces of advice into one stringent analysis. I feel like I get a lot of suggestions, yet lack the understanding to meaningfully add them together and am often not sure if adding more things just messes everything else up.

Thank you so very much in advance!

Answer 1

With fewer than 100 observations, you are limited in what you can accomplish. The risk in this type of modeling is overfitting the data, matching your model to the peculiarities of the data sample instead of elucidating true underlying associations between predictors and outcome.

For the unpenalized regression models you seem mostly to be using, you typically need on the order of 15 observations per coefficient that you will estimate. That would restrict you to 6 or 7 unpenalized coefficients, one of which you have already assigned to ecosystem. That count of coefficients includes not only the specific predictors you have in mind, but also any interaction terms needed to evaluate differences in predictor associations with outcome between the ecosystems.

In your work so far, you have exacerbated that problem with data-driven attempts to transform the data and stepwise predictor selection. Any use of the outcomes to choose the predictors or their transformations will, at the least, make p-values and the like unreliable. Stepwise or other methods for automated predictor selection are poor choices.

What you need is a model that can flexibly fit the data in a way that lets the data tell you the functional shapes of the association between predictors and outcome, but done in a way that penalizes the magnitudes of coefficients to avoid overfitting the data at hand. A generalized additive model (GAM) can be a good choice for this type of data.

There's a reasonably simple explanation of the principles in Section 7.7 of An Introduction to Statistical Learning. The mgcv package is a popular implementation in R. This site has nearly 900 questions tagged with generalized-additive-model, and several individuals with experience in such models who regularly contribute to this site.