# bayesian predictions in multilevel model for panel data

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.

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.

edit:

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)
plot(density(a))
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.

• 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? – Activation Feb 11 at 12:04
• I've added an edit expanding on that question. Hope it helps – snickerdoodles777 Feb 11 at 16:20
• 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$? – Activation Feb 12 at 8:19