# Zero-inflated highly skewed predictor variables

I've thoroughly searched this website and multiple others and can't seem to find an answer to my question. This is also my first post so I hope I've followed all the rules. I apologise for the length, it's complicated!

I'm modelling prey fish functional group biomass (my response variable) against 16 predictor benthic habitat variables (in the form of % cover). Although the majority of the benthic variables are fine, some of them are causing me problems.

1. Several of the habitat variables are highly zero inflated and skewed. Out of the 40 sites I surveyed, some benthic variables were only recorded at 3 or 4 sites and were completely absent from the others. As an example Silt was absent from 38 sites, and had only 1% cover at the remaining two sites.
2. There are also some variables with clear outliers with much higher % cover. For example, sponge cover was minimal at all sites with the exception of one where it was really high. Although outliers, there is obviously huge variation in the natural environment so these values are not errors.

My problem is that when I run my GLMs (using glmulti to find the best model due to the high number of predictor variables), what seems to happen is these zero inflated, highly skewed predictor variables or the variables with the outliers are almost always included in the final best model. They often have a significant interaction which seems solely driven by the outlier sites or skewed nature of the variables and seems to hide other actual relationships occurring. I am fairly sure this isn't right but I'm unsure of the most appropriate way to move forward.

An example:

After running glmulti, best model of PlanktivoreBiomass is below. As you can see EnCA is zero inflated and highly skewed yet included, and sponge has a clear outlier which seems to be driving the positive relationship.

PlanktivoreBiomass ~ 1 + Complexity + Cyanobac + EncrustingCoralAlgal + Sponge

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  401.162    127.756   3.140  0.00342 **
Comp        -421.337    185.999  -2.265  0.02979 *
Cyano         14.236      5.540   2.569  0.01460 *
EnCA         -46.886     20.768  -2.258  0.03031 *
Sponge         9.325      2.992   3.117  0.00364 **


[

I've tried two things:

1. removing the outlier sites, rerunning the models then comparing AIC values of the models and z-values of predictor variables. As soon as those outlier sites are removed the variables are no longer included and the best model is more sensible (based on my ecological knowledge). However, I do then lose important information on the other variables from that site.
2. removing all variables recorded at less than 5% of sites before all analysis.

Are either of those two approaches appropriate? Ideally I would like to follow the approach of removing the variables prior to analysis as then I keep the other information from each site. I have documented everything and in my methods I would explain the approach I have taken plus add an appendix showing the impact the variable removal has on the final models (vs leaving them in for example).

I should add that although I present the final best model for each group, I also build all models within 2 AIC of it and average them out to present the top averaged model results also.