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 Ive 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 Im not sure what can be done with this value. One cant simply take the
sq rt to get a standard error (of sorts), can they? I dont 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 Ive 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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
