variance vs. standard deviation in beta regression

Question

I have modeled a relationship using beta-regression in the betareg package for R. I can predict values and associated variances of those predictions using the predict function in R. I'd like to convert these variances to standard deviations. This seems striaghtforward, just take the square root of the variance. However, there is an issue with this. First, the standard deviation is larger than the variance in beta regerssion, since the variance is always less than 1, and the square-root of a value less than 1 is always larger than the initial value. Second, and more concerning, when I convert to standard deviation by taking the square root of the variance, I get confidence intervals that overlap 0 and 1. This is a problem, because by its very nature beta-regression is designed to model values on the interval (0,1). Uncertainties beyond 0 and 1 are nonsensical.

How can I go about calculating standard deviation from a variance estimate from a beta-regression model, such that the standard deviation envelope will never be less than 0 or greater than 1?

Answer 1

I think there are two potential sources of confusion here: (1) What the variance pertains to. (2) What kind of intervals are computed.

1. The variance is the predicted variance of the response and not the variance of the predicted mean. Thus, it does not correspond to predict(..., se.fit = TRUE) in an lm() regression. Instead, it would correspond to the residual variance which is assumed to be constant in lm() regressions but not in beta regressions.

2. As the variance pertains to the response (see 1) you should not use +/- 2 times the standard deviation to obtain a confidence interval. A normal approximation makes no sense here. Instead, it would be better to look at the, say, 2.5% and 97.5% quantiles of the predicted beta distribution. You can obtain these in betareg with predict(..., type = "quantile", at = c(0.025, 0.975)).

Finally, if you are indeed looking for the equivalent of predict.lm(..., se.fit = TRUE) then betareg currently has no infrastructure for this. An alternative would be to bootstrap the predictions yourself. Or you could look at what the package lsmeans offers for betareg objects.