 Solution

For the Gaussian linear mixed model, temporarily regard $\beta ,D,R$ as known, with $R={\sigma }^{2}I$ here and $V=ZD{Z}^T+R$. The joint multivariate normal distribution isBy the conditional multivariate normal distribution formula,A nonsingular Gaussian density has its mode at its mean. Thus the conditional mode of Gaussian random effects equals the conditional expectation. Replacing the unknown parameters by fitted values givesFor positive-definite $D$ and $R$, an equivalent derivation minimizes the negative conditional log density,Differentiation yieldsThe inverse-covariance penalty shrinks furnace intercepts and slopes toward the population effects. For furnace $j$, set $Z_j=(1,t_j)$ and $D_j=\mathrm{diag}({\tau }_0^2,{\tau }_1^2)$; the same formula provides its two-component conditional mode from that furnace's observations, conditional on the fitted common parameters. If a variance component is zero, the covariance form remains meaningful and its corresponding random effect is zero, whereas the formula with $D^{-1}$ needs a limiting interpretation. With known covariance parameters these are best linear unbiased predictors; plug-in variance estimates give their empirical version.