Dear all,
I am trying to get an estimate of uncertainty surrounding a single predicted value from a beta regression model (this is similar to a logistic glm - in that it involves a link function and linear predictor - but it uses the beta distribution rather than discrete binomial). For example:
library(betareg)
data("GasolineYield")
test.model<-betareg(yield~gravity+temp,data=GasolineYield)
...and I would like a measure of uncertainty for a single predicted value, such as:
predicted.value<-predict(test.model,newdata=GasolineYield[20,],type="response")
Ideally, this interval estimate would be a confidence interval (more specifically, a prediction confidence interval) or, failing that, a standard error.
However, using predict.betareg {betareg} there is no equivalent to the ‘se.fit=T’ or ‘interval=”confidence”’ arguments that are present in predict.lm(). I note, also, that there is no ‘interval=”confidence”’ in predict.glm. Am I missing something or is this an unresolved issue for beta regression models (and perhaps glms as well)?
It is possible, however, to get the “fitted variances of response” (quoted from ?predict.betareg). For example:
var1<-predict(test.model,newdata=GasolineYield[20,],type="variance")
This seems to be the variance as got from the mean-variance relationship for the beta distribution (though I’ve been unsuccessful in being able to unpack the source code tar.gz file from CRAN in Windows to see this):
precis.param<- predict(test.model,newdata=GasolineYield[20,],type="precision")
#This is from the description of the beta distribution in ?VGAM::betaff:
var2<-predicted.value*(1-predicted.value)/(1+precis.param)
var1==var2
...but I’m not sure what can be done with this value. One can’t simply take the sq rt to get a “standard error” (of sorts), can they? I don’t think, for example, this is taking much account of the uncertainty in the regression parameter estimates.
Even if I was able to get a SE interval, how might I go about converting this to a confidence interval? Would it be correct to use a quantile from the standard normal distribution (i.e. CI<-SE*qnorm(p=0.975)) or, rather, from the beta distribution? I think I’ve seen the former being suggested in the context of non-normal glms (e.g. binomial errors) but the latter would make more sense to me here because, for example, the distribution would be suitably asymmetric near the bounds of the response, i.e. zero and 1.
I found the relevant parameterisation of the beta distribution (involving mu and a ‘precision parameter’, already estimated in the betareg model) in the package ‘gamlss.dist’, so this could be done with:
library(gamlss.dist)
CI<-SE * qBE(p=0.975,mu=predicted.value,sigma=sqrt(var1))
#...where SE has already been calculated somehow!
Thanks in advance for any help/ideas/suggestions! As you can see, I am out of my depth here!
Many thanks,
Oliver

	[[alternative HTML version deleted]]

