Prediction intervals for the outcome of a logistic regression with binomial response

Suppose we have a logistic regression model:

\begin{align} P(y=1\vert\mathbf{x}) &= p \\ \log\left(\frac{p}{1-p}\right) &= \boldsymbol{\beta}\mathbf{x} \end{align}

Given a random sample $D=\{\mathbf{X},\mathbf{y}\}$ of size $N$, we can compute confidence intervals for the $\boldsymbol{\beta}$ and correspondingly prediction intervals for $p$, given a certain value $\mathbf{x}^*$ of the predictor vector. This is all very standard and detailed, for example, here.

Suppose instead that I'm interested in a prediction interval for $y$, given $\mathbf{x}^*$. Of course, it doesn't make any sense at all to compute a prediction interval for a single realization of $y$, because $y$ can only take the values 0 and 1, and no value in between. However, if we consider $m$ realizations of $y$ for the same fixed value of $\mathbf{x}^*$ , then this becomes similar (but not identical) to the question of computing a prediction interval for a binomial random variable. This is basically the same situation described by Glen_b in the comments to this answer. Does this question have an answer, apart from the trivial one "use nonparametric bootstrap"?

• can you compute a prediction interval for $log(p / (1-p))$ instead perhaps? – Hugh Perkins Mar 31 '18 at 1:15
• @HughPerkins I think that the issue is how to combine the uncertainty in p with the uncertainty in binomial sampling also given the uncertainty in p. Is there a closed-form solution? – EdM Mar 31 '18 at 2:17
• @EdM you got my point. I wonder if there is a closed form solution or an analytical approximation. – DeltaIV Mar 31 '18 at 12:03
• [offtopic] random idea, it occurs to me that it could be interesting to have a tag like 'open-research-opportunity' for questions like this which/if they are answered in the negative – Hugh Perkins Mar 31 '18 at 16:05

One way this should work without bootstrapping (which in practice may be the fastest thing tho implement), would be:

1. Assume that a normal approximation for the predicted log-odds ($x \hat{\beta}$) plus/minus its standard error works. Any logistic regression software will provide this.
2. The percentiles of this distribution transform to probabilities via the anti-logit.
3. One can find a (mixture of) beta distribution(s) that approximates the predictive distribution for the probability well.
4. The predictive distribution for the outcome is then a (mixture of) beta-binomial distribution(s with the same mixing weights as used in step 3).

Alternatively, one can "just" integrate out the log-odds from the joint predictive of outcome and log-odds, but I believe that will be a complete mess with no closed form solution.

• You could also just directly simulate from the asymptotic multivariate normal for $\beta-\hat{\beta}$, and then form a mixture of binomials over those values. – Glen_b Apr 1 '18 at 6:28
• I like the overall idea, but I'm not sure about the details. For example, "find a (mixture of) beta distribution(s) that approximates the predictive distribution for the probability well", How do you this in practice? Could you add an example? Even a low-dimensional one would suffice. – DeltaIV Apr 2 '18 at 6:07
• About your second alternative, if I understand correctly you suggest to write down the joint posterior distribution of outcomes and log-odds (which is a joint distribution between continuous and discrete variables) and integrate out log-odds. I definitely agree that that's going to be a mess. One could salvage the approach using MCMC to estimate the integral, I think. However, I would need MCMC algorithms which work for discrete-continuous variables...I'm not sure which one do. I think rstan (my tool of choice for MCMC) doesn't. – DeltaIV Apr 2 '18 at 6:11
• I can write this up as something in the form of an answer if you prefer -- I don't mind either way. – Glen_b Apr 3 '18 at 14:58
• @Glen_b I'd really appreciate that. – DeltaIV Apr 6 '18 at 18:25