Suppose I fit a Binomial regression and obtain the point estimates and variance-covariance matrix of the regression coefficients. That will allow me to get a CI for the expected proportion of successes in a future experiment, $p$, but I need a CI for the observed proportion. There have been a few related answers posted, including simulation (suppose I don't want to do that) and a link to Krishnamoorthya et al (which does not quite answer my question).
My reasoning is as follows: if we use just the Binomial model, we are forced to assume that $p$ is sampled from Normal distribution (with the corresponding Wald CI) and therefore it is impossible to get CI for the observed proportion in closed form. If we assume that $p$ is sampled from beta distribution, then things are much easier because the count of successes will follow Beta-Binomial distribution. We will have to assume that there is no uncertainty in the estimated beta parameters, $\alpha$ and $\beta$.
There are three questions:
1) A theoretical one: is it ok to use just the point estimates of beta parameters? I know that to construct a CI for future observation in multiple linear regression
$Y = x'\beta + \epsilon, \epsilon \sim N(0, \sigma^2)$
they do that w.r.t. error term variance, $\sigma^2$. I take it (correct me if I am wrong) that the justification is that in practice $\sigma^2$ is estimated with a far greater precision than the regression coefficients and we won't gain much by trying to incorporate the uncertainty of $\sigma^2$. Is a similar justification applicable to the estimated beta parameters, $\alpha$ and $\beta$?
2) What package is better (R: gamlss-bb, betareg, aod?; I also have access to SAS).
3) Given the estimated beta parameters, is there an (approximate) shortcut to get the quantiles (2.5%, 97.5%) for the count of future successes or, better yet, for the proportion of future successes under Beta-Binomial distribution.