How to make correct predictions of probabilities and their uncertainty in Bayesian logistic regression?

Question

In the context of Bayesian logistic regression, outcomes $y$ are binary (discrete) and covariates $X$ are given. We assume in particular:

$$ p(y | x, \theta) = Bin(n,\theta) \newcommand{\new}{{\rm new}}$$

where we model $\theta$ as:

$$\theta = {\rm logit}^{-1}(X\beta) $$

so

$$ p(y | x, \beta) = Bin(n,{\rm logit}^{-1}(X\beta)) \newcommand{\new}{{\rm new}}$$

with regression parameters $\beta$.

Using a prior distribution $p(\alpha)$ for $\beta$ we come to the posterior:

$$p(\beta|y,X;\alpha) \propto p(y|X,\beta) p(\alpha)$$

where $p(y|X,\beta)$ the conditional likelihood.

We can now make two types of predictions for a new observation $X_{\new}$, one for $y_{\new} | X_{\new}, X, y \in \{0,1\} $ and one for $p(y_{\new} = 1| X_{\new},X, y) \in (0,1)$.

My question is: How can one make these predictions and quantify the uncertainty in the predictions given the data from the sample $(X,y)$?

I believe the posterior predictive distribution (ppd) is $$p(y_{new}|X,y,X_{new}) = \int p(y_{\new}|X_{\new},X,\beta) p(\beta|y,X) d\beta.$$

We can thus sample from the ppd $D$ times leading to $d=1,...,D$ samples $\tilde{y}_{\new,d}$ from $p(y_{\new}|X,\theta)$. Thus it is clear how to generate predictions for $y_{\new}$, but it is not clear to me how to generate a prediction for $p(y_{\new} = 1| X_{\new},y)$ and quantify its uncertainty.

I believe it is not central to my question, which prior to choose for $\beta$. However, a normal prior would be okay, so that samples from the posterior can be generated by MCMC or Laplace approximation.

Answer 1

I'm going to start by cleaning up a little notational lack of precision in the first part of the question, then go on to the meat of it.

First, if $y \in \{0,1\}$, $y$ is, more precisely than stated in the question, distributed Bernoulli$(\theta)$, not Binomial$(n,\theta)$ for some arbitrary $n$. The two distributions are of course the same when $n=1$, but best to get rid of $n$ altogether if you don't use it subsequently.

On to the question! What the posterior distribution $p(\beta|\cdot)$ allows us to do is calculate an MCMC-based posterior for $\theta$ by sampling from the posterior distribution of $\beta$, then multiplying each sample $\beta$ by $X_{\text{new}}$ and calculating the inverse logit of the result (as in your second equation above.) This gives us a vector of samples of $\theta$ from the posterior distribution of $\theta$. Thanks to the Bernoulli assumption, these samples are samples from the posterior distribution of the expected value of $y_{\text{new}}$ as well as from the posterior distribution of the probability that $y_{\text{new}} = 1 | X_{\text{new}}, y, X$.

Note that we are not using the posterior predictive distribution to do this, just the posterior distribution itself. Your confusion stems from the fact that you tried to work with the posterior predictive distribution of $y_{\text{new}}$, which, for the problem you faced, was a step too far.

The samples of $\theta$ thus generated can be used in the usual MCMC ways in place of the full posterior distribution for $\theta$ or as inputs to a smoothing function of some sort (for calculating a smooth version of the posterior.) Either way, the sample-based posterior distribution contains the information you need to find (estimates of) the expected value of $p(y_{\text{new}} = 1| X_{\text{new}}, y, X)$, (estimates of) quantiles of same, etc. For example, the posterior mean of $p(y_{\text{new}} = 1| X_{\text{new}}, y, X)$ is just the mean of the sampled $\theta$ values, the 5th percentile is just the 5th percentile of the sampled $\theta$ values, etc.