I want to make predictions for a bayesian multilevel which basically looks like: $y_{it} = \alpha_{i} + x_{it}\beta$.

I was told that I could make predictions by using (in the case that $y_{it}$ is normally dsitributed): $\hat{y}_{it} = E(y_{it}) = \frac{1}{M}\sum_{m = 1}^{M}(\alpha_i^{(m)} + x_{it}\beta^{(m)})$.

However, I do not fully understand how this follows from the formula to make predictions: $p(y_{N+x}) = \int_{-\infty}^{+\infty} p(y_{N+x}|\theta)p(\theta|y)d\theta$.

I searched for papers all over the internet but I could not find one which explains this. Could someone explain the relationship to me or knows a paper where this is explained.

Thanks in advance!


Lets go ahead and unpack your second equation first. This is commonly called the posterior predictive distribution, and represents your beliefs about the $N+x$th datapoint conditioned on the data you have already seen. The first term $p(y_{N+x} | \alpha, \beta)$ is the model likelihood, and can be loosely interpreted as describing how probable a y value is, for a specific value of $\alpha, \beta$. Since you have already specified that y is normal, we have that $y_i|\alpha_i, \beta \sim N(\alpha_i + x_{it}\beta, \sigma^2)$. The second term $p(\alpha, \beta|y)$ is the posterior, and represents our beliefs about $\alpha, \beta$ after we have observed our data.

Now, because we have some uncertainty about the true values of our parameters, we don't want to simply use one value of our $\alpha, \beta$ to determine our predictions. Instead, we want to take into account all possible values of $\alpha, \beta$, and compute the prediction for each one, while weighting by how probable these values actually are. This is what corresponds to that integral that we see - $p(y|\alpha, \beta)$ represents the "prediction" part (sort of), while $p(\alpha, \beta | y)$ represents how probable these given values of $\alpha, \beta$ are.

In practice, the posterior is intractable, so we have to resort to sampling methods and represent our posterior with a finite number of samples. Thus our sample posterior will be discrete, and can be represented as

$$ \hat{p}(\alpha, \beta | y) = \frac{1}{M}$$

Because of the finite number of our samples, the predictive distribution equation then turns from an integral into a sum:

$$ \hat{p}(Y_{i,N+x} | y) = \sum_{m=1}^M p(Y_{i,N+x}|\alpha_i^m, \beta^m)\hat{p}(\alpha^m, \beta^m | y) = \sum_{m=1}^M \frac{p(Y_{i,N+x} | \alpha_i^m, \beta^m)}{M}$$

Thus, we have an approximation to our posterior predictive distribution. But, in this case, we don't just want a distribution, we want a specific point value prediction. So, lets take the expected value of this distribution with respect to $Y_{N+x}$ and use that as our prediction. Recall that the expected value of a normal is simply its mean: that means that our expected value of $y_{it}$ for a given value of $\alpha_i, \beta$ is simply $\alpha_i + x_{it}\beta$. Therefore, putting all of the pieces together, we have

$$\hat{Y_{i,N+x}} = \frac{1}{M}\sum_{m=1}^M \alpha_{i} + x_{i,N+x}\beta$$

You can see that this is an approximation to the expected value of the posterior predictive distribution: we are making a prediction for each given value of $\alpha, \beta$, while weighting by how probable those values actually are.


You asked for more information about why $ \hat{p}(\alpha, \beta | y) = \frac{1}{M}$. $ \hat{p}(\alpha, \beta | y)$ is whats called an empirical distribution: it is comprised of a finite sample of draws from some underlying true distribution. Since we want the posterior to be an actual probability distribution, we have the constraint that $\sum_{m=1}^M \hat{p}(\alpha, \beta | y) = 1$. The most common way to deal with this constraint is to set $\hat{p}(\alpha, \beta | y) = \frac{1}{M}$ for all values. As for the prior and likelihood, they are incorporated via the locations of the samples, not the probabilities. Lets consider a simple example.

a <- rnorm(100, 0, 1)
points(a, rep(0, len(a)))


Note that each point is given probability mass $\frac{1}{100}$, but it is the location of the points that encodes the true density. The $\frac{1}{100}$ is just there to ensure that the probability distribution sums to 1.

  • $\begingroup$ Thanks, great explanation. The only thing I don't fully understand is why we can approximate the posterior distribution $\hat{p}(\alpha, \beta|y) = \frac{1}{M}$, since it should be proportional to prior times likelihood, right? Could you maybe explain this a little bit more? $\endgroup$ – Activation Feb 11 at 12:04
  • $\begingroup$ I've added an edit expanding on that question. Hope it helps $\endgroup$ – snickerdoodles777 Feb 11 at 16:20
  • $\begingroup$ One final question/remark. Shouldn't the restriction $\sum_{m=1}^{M} \hat{p}(\alpha, \beta | y) = 1$ be $\sum_{m=1}^{M} \hat{p}(\alpha^{(m)}, \beta^{(m)} | y) = 1$? $\endgroup$ – Activation Feb 12 at 8:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.