What does it mean that BLUP is unbiased, given a linear two-level model?

Question

Suppose we have the following mixed effects model for observation $Y_{ij}$ of pupil $i$ in school $j$:

$Y_{ij}=b_0 + u_j + e_{ij}$

Here, $b_0$ is a fixed parameter for the "grand mean", $u_j$ is the random effect of school $j$ and $e_{ij}$ is the residual term for pupil $i$ in school $j$. Once the fixed and random effects have been estimated, the BLUP for each school can be calculated, which shows shrinkage towards the grand mean. I think this entails that the BLUPs of the schools are all biased towards the grand mean; that is, as a prediction of the unknown true $Y$ mean of a given school, the BLUP is biased. However, the U in "BLUP" means "Unbiased". So, my question is: how must I understand this unbiasedness of the BLUP?

My guess is that it means that one must calculate the BLUP for all schools in a sample, and then calculate the difference between the true school means (could be simulated) and the BLUPs. The mean of these differences should approach zero, if the sample is large enough and the true school means are normally distributed, just like the error terms. But is this indeed what the unbiasedness of the BLUP means?

Answer 1

Your guess is a very good one! My go-to reference for these types of questions is the Multilevel and Longitudinal Modeling in Stata book by Sophia Rabe-Hesketh and Anders Skrondal. In chapter 2, they address the unbiasedness part of the BLUP.

In particular, the prediction error for the empirical Bayes prediction, "defined as $\hat{u_i}^{EB} - u_i$ between the prediction and the truth has zero mean over repeated samples of $u_i$ and $e_{ij}$ (or repeated samples of clusters and units from clusters) when model parameters are treated as fixed and known." The model parameters in this case are $\hat\beta$, $\hat\psi$, and $\hat\theta$, corresponding to, respectively, the estimates of the grand mean (fixed intercept) and variances of the random intercept and residual.

Empirical Bayes predictions also have the smallest possible variance (for given model parameters) versus other methods such as maximum likelihood estimates. The latter however, do not employ shrinkage, which leads to a decreased conditional bias relative to the EB predictions. But the unbiasedness in BLUP comes from the small prediction-error variance.

Answer 2

After reading Erik Ruzek's answer, the following addition makes sense, I hope. I will use Erik's notation $u^{EB}_j$ for the BLUP of the term $u_j$. So, we have symbol $u_j$ for school $j$ 's true deviation from the grand mean, and we have $u^{EB}_j$ for the BLUP of this deviation, which is also called "empirical Bayes estimate" of $ u_j$.

The unbiasedness of the BLUP can then be formulated as

$E(u^{EB}_j - u_j) = 0$

where the expectation is taken over all schools. This can also be written as:

$E(u^{EB}_j) - E(u_j) = 0$

and since $E(u_j)=0$ by assumption, we finally have

$E(u^{EB}_j) = 0$

or: the mean of all schools' BLUP (or emperical BAYES) deviations from grand mean $b_0$ is 0.