How to approach GLMs using data with beta distribution in R?

Question

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.

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).

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.

Answer 1

The Beta regression family is not a generalised linear model (in any of the strict senses) and so glm can't give you the maximum likelihood estimators (not with the pre-defined family options and not with user-defined ones either).

You can get a reasonable approximation to the Beta model from a glm with family=quasi(link="logit", variance="mu(1-mu)"), since the variance of a Beta with mean $\mu=\alpha/(\alpha+\beta)$ is $$\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}=\mu(1-\mu)\phi$$ with $\phi=1/(\alpha+\beta+1)$ being a dispersion parameter. This is obviously only useful for estimating trends, not for prediction intervals, since there's no parametric model.

In fact, even modelling $\mathrm{logit}\,Y$ with linear regression and a sandwich variance estimator works pretty well if you are interested in trends (you need the sandwich estimator because the variance of $\mathrm{logit}\, Y$ increases as $Y$ moves away from 1/2).

But there is in fact software for maximum likelihood fitting (eg betareg), so you might as well use that instead of forcing it into glm.

All these assume the outcome is within some prespecified interval: by default $[0,1]$ but you can transform the data if the interval is something else. That's pretty fundamental to these models because the main challenge they are dealing with is the hard limits at both ends of the interval.