I’m trying to run some models on bee presence with five predictor variables. A snippet of the data is attached, but essentially I measured floral abundance and richness, calculated floral evenness and diversity, measured canopy cover (%), and the dependent variable is the guild/functional group diversity (last column). Elevation band is a factor with 5 levels. Snippet of dataset. Column labels are: ElevationBand, FloralAbund, FloralRich, FloralEven, FloralDiv, Canopy, GuildDiv

I have tested for normality; the response data is not normally distributed even if log transformed.

I tried fitting to Poisson and negative binomial, but it appears that beta fits best. I used descdist(fxndat$GuildDiv, boot=10000) to get the Cullen and Frey graph (attached). Cullen and Frey graph showing beta distribution I’m pretty comfortable with the glm() function with normally distributed data, but I’m having trouble finding documentation on how to approach modeling with the beta distribution as that’s not an option for the “family” argument (or at least, it’s not listed in the glm function documentation).

Is there a different function I should use? Additionally, my diversity indices don’t necessarily fit within the (0,1) parameter that I’ve read for the beta distribution, so I’m not sure I should even use a GLM.

Thank you in advance!