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)}). This way works fine.

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!