Plug-In Parameter Estimates in Empirical Bayes Prediction

Why Parameter Substitution Matters

An empirical Bayes prediction for a group effect is conditional on both the group’s observations and the model’s population-level parameters. In practice, those parameters are not usually known: the fixed effects and variance components are estimated from the data. A plug-in prediction uses those fitted values in the conditional distribution of the group effect. Thus, it is not just a calculation based on the group’s own data; it also depends on how the fitted model describes the population and the variation among groups.

The key distinction is between predicting with specified parameter values and predicting after fitting those values. Once estimates are substituted, the resulting prediction is a function of the observed data and the fitted parameters. The calculation can be exact for a Gaussian model conditional on the fitted values, even though the parameter values themselves were estimated.

Definition: Let $y_j$ denote the observations for group $j$, let $b_j$ be its random effect, and let $\theta$ collect the fixed effects and variance components. The empirical Bayes prediction using specified parameters is the conditional mean $E_{\theta}(b_j\mid y_j)$. The plug-in empirical Bayes prediction is $$ \hat b_j = E_{\hat\theta}(b_j\mid y_j), $$ where $\hat\theta$ denotes fitted parameter estimates. The notation means that the conditional prediction formula is evaluated at the estimates; it does not assert that $\hat\theta$ is known without error.

Gaussian Random-Intercept Model

Consider a continuous outcome $y_{ij}$ for observation $i$ in group $j$. In a Gaussian random-intercept model,

y_{ij} = \mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta} + b_j + \varepsilon_{ij}, \qquad b_j \sim N(0,\tau^2), \qquad \varepsilon_{ij} \sim N(0,\sigma^2).

Here, $\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}$ is the fixed-effects mean, $b_j$ is the group-specific deviation from that mean, $\tau^2$ is the random-intercept variance, and $\sigma^2$ is the residual variance. The standard assumptions for the calculation below are that the random effects are independent across groups; the residuals are independent across observations and groups; the random effects and residuals are mutually independent; and both have the stated normal distributions. The fixed-effects design values are treated as observed.

For clarity, first suppose the fixed-effects mean is a common intercept $\beta_0$ for all observations in a group. Write $\bar y_j$ for the group mean and $n_j$ for its size. Conditional on $b_j$, the group mean has distribution $N(\beta_0+b_j,\sigma^2/n_j)$. Its information about $b_j$ therefore depends on the group mean after subtracting the fixed-effects mean, not on the group mean alone.

Deriving the Plug-In Formula

With specified parameters, the prior density for $b_j$ is proportional to $\exp\{-b_j^2/(2\tau^2)\}$. The likelihood contribution from the group mean is proportional to $\exp\{-(\bar y_j-\beta_0-b_j)^2/(2\sigma^2/n_j)\}$. Multiplying these gives a normal conditional density for $b_j$. To find its mean, collect the terms involving $b_j$ in the exponent:

-\frac{1}{2}\left[ \frac{b_j^2}{\tau^2} + \frac{(\bar y_j-\beta_0-b_j)^2}{\sigma^2/n_j} \right].

The coefficient of $b_j^2$ inside the brackets is $1/\tau^2+n_j/\sigma^2$, and the linear coefficient is $-2b_j n_j(\bar y_j-\beta_0)/\sigma^2$. Completing the square gives conditional variance $V_j$ and conditional mean $m_j$:

V_j = \left(\frac{1}{\tau^2}+\frac{n_j}{\sigma^2}\right)^{-1}, \qquad m_j = V_j\frac{n_j}{\sigma^2}(\bar y_j-\beta_0) = \frac{n_j\tau^2}{\sigma^2+n_j\tau^2}(\bar y_j-\beta_0).

Because this conditional distribution is normal, its mean is also its mode. Substituting estimates for the parameters yields the plug-in prediction:

\hat b_j = \frac{n_j\hat\tau^2}{\hat\sigma^2+n_j\hat\tau^2} (\bar y_j-\hat\beta_0).

The factor multiplying the residual group mean is the shrinkage weight. It lies between zero and one when the estimated variance components are positive. The group’s mean is first centered at the fitted fixed-effects mean; the resulting deviation is then pulled toward zero. Both parts of this calculation use fitted parameters: $\hat\beta_0$ determines the centering, while $\hat\tau^2$ and $\hat\sigma^2$ determine the weight.

Worked Example: Substituting Fitted Parameters

Suppose a hypothetical Gaussian random-intercept analysis gives $\hat\beta_0=10.0$ units, $\hat\tau^2=4.00$ square units, and $\hat\sigma^2=4.00$ square units. One group has $n_j=4$ observations and mean $\bar y_j=12.0$ units. Calculate its plug-in prediction.

The group’s observed mean is $12.0-10.0=2.00$ units above the fitted fixed-effects mean. The plug-in shrinkage weight is $$ \hat\lambda_j = \frac{4(4.00)}{4.00+4(4.00)} = \frac{16.0}{20.0} =0.800. $$

Therefore, $$ \hat b_j = 0.800(12.0-10.0)=1.60 \text{ units}. $$ The predicted group-specific mean is $10.0+1.60=11.6$ units. It is below the group’s observed mean of $12.0$ units because the prediction combines the group’s data with the fitted distribution of group effects. This is a prediction of the group’s latent deviation, not a claim that its observed mean must equal $11.6$.

How Fitted Fixed Effects Enter

With covariates, the same principle applies, but the group’s data must be compared with its fitted fixed-effects mean. In the random-intercept model, suppose each observation has predictor $x_{ij}$ and fixed-effects mean $\beta_0+\beta_1x_{ij}$. Define the average fitted fixed-effects mean in group $j$ as $\hat\beta_0+\hat\beta_1\bar x_j$, where $\bar x_j$ is the group’s average predictor. Under the model’s independent Gaussian residual assumption, the group-mean residual is $\bar y_j-\hat\beta_0-\hat\beta_1\bar x_j$. The plug-in prediction multiplies this residual by the corresponding weight.

Worked Example: A Fitted Slope Changes the Prediction

Consider a hypothetical group with $n_j=3$, $\bar y_j=8.40$ outcome units, and $\bar x_j=2.00$ predictor units. Suppose $\hat\beta_0=5.00$ outcome units, $\hat\beta_1=1.20$ outcome units per predictor unit, $\hat\tau^2=1.00$ outcome-units squared, and $\hat\sigma^2=3.00$ outcome-units squared. Calculate the prediction, then see what changes if the fitted slope is $1.00$ instead, with the other fitted values held fixed.

Using the stated slope, the fitted fixed-effects mean for this group is $$ \hat\beta_0+\hat\beta_1\bar x_j =5.00+(1.20)(2.00)=7.40 \text{ outcome units}. $$ The group-mean residual is $8.40-7.40=1.00$ outcome unit. Its shrinkage weight is $$ \hat\lambda_j = \frac{3(1.00)}{3.00+3(1.00)} =\frac{3.00}{6.00}=0.500. $$ Thus, $\hat b_j=0.500(1.00)=0.500$ outcome unit.

If the fitted slope is instead $1.00$ outcome unit per predictor unit, the fitted mean becomes $5.00+(1.00)(2.00)=7.00$ outcome units. The residual is now $8.40-7.00=1.40$ outcome units. With the same variance estimates, the weight remains $0.500$, so the prediction becomes $0.500(1.40)=0.700$ outcome unit. This comparison isolates the fixed-effects contribution: changing the fitted slope changes the centered group residual and therefore the predicted intercept, even though the group’s observations and shrinkage weight have not changed.

Using the Formula Across Groups

A fitted model uses common parameter estimates across groups, but each group can have a different mean, size, or covariate pattern. Consequently, its plug-in prediction need not be the same as another group’s, even when both predictions use the same fitted variance components. The next example applies one set of estimates to two groups and keeps the units explicit.

Worked Example: Comparing Two Group Predictions

Suppose a hypothetical model has $\hat\beta_0=20.0$ mg, $\hat\tau^2=9.00$ mg$^2$, and $\hat\sigma^2=16.0$ mg$^2$. Group A has $n_A=2$ and $\bar y_A=24.0$ mg; Group B has $n_B=8$ and $\bar y_B=21.0$ mg. Assume the Gaussian random-intercept model and its independence assumptions are appropriate. Calculate both plug-in predictions.

For Group A, the weight and prediction are $$ \hat\lambda_A = \frac{2(9.00)}{16.0+2(9.00)} =\frac{18.0}{34.0}=0.5294, \qquad \hat b_A=0.5294(24.0-20.0)=2.118 \text{ mg}. $$

For Group B, $$ \hat\lambda_B = \frac{8(9.00)}{16.0+8(9.00)} =\frac{72.0}{88.0}=0.8182, \qquad \hat b_B=0.8182(21.0-20.0)=0.8182 \text{ mg}. $$

Group A has a larger observed deviation from the fitted intercept, but its predicted effect is only about $2.12$ mg. Group B’s observed deviation is $1.00$ mg and its prediction is about $0.818$ mg. The larger group receives the larger weight because its mean has less residual sampling variation under the model. These are predictions of group deviations; adding $\hat\beta_0$ gives the corresponding fitted group-specific means, $22.1$ mg for A and $20.8$ mg for B.

Beyond the Gaussian Case

The same plug-in definition applies in a non-Gaussian GLMM, but the conditional distribution of a random effect is generally not normal and its mean need not have a closed-form expression. For a group with data $y_j$, conditional likelihood $L_j(b;\boldsymbol{\beta})$, and random-effects density $f(b;\boldsymbol{\psi})$, the plug-in posterior mean is

\hat b_j = \frac{\int b\,L_j(b;\hat{\boldsymbol{\beta}})f(b;\hat{\boldsymbol{\psi}})\,db} {\int L_j(b;\hat{\boldsymbol{\beta}})f(b;\hat{\boldsymbol{\psi}})\,db}.

Here, $\boldsymbol{\beta}$ denotes the fixed effects and $\boldsymbol{\psi}$ the random-effects parameters. The estimated fixed effects enter through the conditional likelihood, while the estimated random-effects parameters enter through the distribution assigned to $b$. Numerical integration or approximation may be needed. In a Gaussian model, the conditional mean and mode coincide; in a non-Gaussian model they can differ, so the target of the reported prediction must be specified. This is distinct from a population-averaged fitted mean, which averages over the random-effects distribution rather than predicting a particular group’s effect.

Diagnostics, Alternatives, and Limitations

A plug-in prediction inherits problems in either the fitted mean structure or the fitted variance structure. Before interpreting a group prediction, check that the model’s conditional distribution and dependence assumptions are plausible for the data. For a Gaussian model, residual plots can reveal strong departures from the assumed residual pattern, and group-level patterns can indicate that the fixed-effects structure is inadequate. Also check that the model fit converged and that parameter estimates are finite and interpretable. If estimates are unstable or the fitted model is poorly specified, the arithmetic may be correct while the prediction remains misleading.

The main alternatives answer different questions. A raw group mean or raw deviation from a fitted mean does not incorporate the random-effects distribution in the same way as an empirical Bayes prediction. A conditional mode is a common alternative to the posterior mean; for the Gaussian random-intercept model shown above, they are identical, but that equality is not general. A population-averaged prediction targets an average over groups, rather than a particular group’s latent deviation. The choice should follow the scientific target.

The plug-in calculation has a specific limitation: it evaluates the conditional prediction at fitted parameters. It does not by itself represent uncertainty in those estimated parameters. Accordingly, a plug-in prediction should not be described as if the parameter estimates were known constants. The prediction is also model-dependent: different fixed-effects specifications or random-effects distributions can produce different group predictions for the same observed data.

Common pitfall: Treating $\bar y_j$ as the group’s residual from the model leads to $\hat b_j=\hat\lambda_j\bar y_j$. This ignores the fitted fixed-effects mean. The correct Gaussian random-intercept calculation centers first: $\hat b_j=\hat\lambda_j(\bar y_j-\hat\beta_0)$, or, with covariates, uses the group mean minus its fitted fixed-effects mean. A second mistake is to call the resulting prediction an observed group mean; the predicted group-specific mean is the fixed-effects mean plus the predicted random effect.

Check Your Understanding

Answer each question using the definitions and calculations above. Solutions follow each question.

A group has $n_j=5$, $\bar y_j=13$ units, $\hat\beta_0=11$ units, $\hat\tau^2=2$ units$^2$, and $\hat\sigma^2=8$ units$^2$. What is its plug-in prediction? Solution: The weight is $5(2)/(8+5(2))=10/18=0.5556$. The centered mean is $13-11=2$ units, so $\hat b_j=0.5556(2)=1.111$ units.
In the previous calculation, what would be wrong with using $0.5556(13)$? Solution: It uses the group mean itself rather than its deviation from the fitted intercept. The prediction must be based on $13-11$, not $13$.
For the common-intercept Gaussian model, which quantities determine the shrinkage weight? Solution: The group size $n_j$, fitted random-intercept variance $\hat\tau^2$, and fitted residual variance $\hat\sigma^2$ determine $\hat\lambda_j=n_j\hat\tau^2/(\hat\sigma^2+n_j\hat\tau^2)$. The fixed-effects estimate determines the centered residual, not this weight.
Why can changing a fitted slope alter a group’s prediction even if its observations and fitted variance components stay the same? Solution: The slope contributes to the fitted fixed-effects mean. Changing it changes the group’s residual from that mean; multiplying the new residual by the unchanged weight gives a different predicted random intercept.
In the Gaussian model, are the conditional mean and conditional mode of $b_j$ different prediction rules? Solution: No. The conditional distribution is normal, whose mean and mode coincide. This equality should not be assumed for a non-Gaussian conditional distribution.