Estimate random effects for a new individual with a linear mixed effects model

Question

Consider repeated observations $\mathcal{Y} = (y_{i,j})_{i,j}$ obtained for $p$ individuals ($1 \leq i \leq p$), at different time points $t_{i,j}$ $(1 \leq j \leq n_i$). The "random slope and intercept" model writes:

$$ y_{i,j} = \left( \beta_0 + b_{i,0} \right) + \left( \beta_1 + b_{i,1} \right) t_{i,j} + \varepsilon_{i,j}, $$

where $\beta = \begin{bmatrix} \beta_0 & \beta_1 \end{bmatrix}^{\top}$ denote the fixed effects of the model and $$b_i = \begin{bmatrix} b_{i,0} & b_{i,1} \end{bmatrix}^{\top} \sim \mathcal{N}\left( 0, \mathbf{D} \right), \quad b_i \perp\kern-5pt\perp \varepsilon_i,$$ denote the random effects, $\varepsilon_{i,j} \sim \mathcal{N}\left( 0, \sigma^2 \right)$. Let $\theta = \left( \beta, \mathrm{vech}\left( \mathbf{D} \right), \sigma^2 \right)$ denote the model parameters.

Given $\mathcal{Y}$, one can obtain an estimator $\hat{\theta}$ of $\theta$ by maximizing the model likelihood (or restricted likelihood). Now, say that we have some data $\mathbf{y}_{\mathrm{new}}^{\ast} = \left( y_{\mathrm{new},1}, \ldots, y_{\mathrm{new}, n^{\ast}}\right)$ for a new individual.

We want to estimate the trajectory (i.e., the straight line) of this new individual. To do that, we only need to estimate his random effects $\mathbf{b}_{\mathrm{new}}$. How do we do that?

1. One could get $\mathbf{b}_{\mathrm{new}}$ from the posterior $p\left( \mathbf{b}_{\mathrm{new}} \mid \mathbf{y}_{\mathrm{new}}, \hat{\theta} \right)$. Unless I am mistaken, this is what D. Rizopoulos proposed in his answer to a similar question. Using the Bayes rule, we get:

$$ p\left( \mathbf{b} \mid \mathbf{y}_{\mathrm{new}}, \hat{\theta} \right) \propto p\left( \mathbf{y}_{\mathrm{new}} \mid \mathbf{b}, \hat{\theta} \right) p\left( \mathbf{b} \mid \hat{\theta} \right), $$

and we could have:

$$ \mathbf{b}_{\mathrm{new}} \in \mathop{\mathrm{argmax}} \limits_{\mathbf{b}} p\left( \mathbf{y}_{\mathrm{new}} \mid \mathbf{b}, \hat{\theta} \right) p\left( \mathbf{b} \mid \hat{\theta} \right), $$

which would yield, unless I am mistaken, the BLUP (Best Linear Unbiased Predictor) of this new individual's random effects.

2. Would it make sense to estimate $\mathbf{b}_{\mathrm{new}}$ by maximizing the following instead?

$$ \int p\left( \mathbf{b} \mid \mathbf{y}_{\mathrm{new}}, \theta \right) p\left( \theta \mid \mathcal{Y} \right) \, d\theta, $$

which would be $\mathbb{E}_{p\left( \theta \mid \mathcal{Y} \right)}\left[ p\left( \mathbf{b} \mid \mathbf{y}_{\mathrm{new}}, \theta \right) \right]$. I am not sure this makes sense but I was thinking of something similar to the posterior predictive distribution.

Answer 1

The difference between the two approaches is actually a difference between an Empirical Bayes versus a fully Bayesian approach to estimate the same thing.

If you fit the mixed model using maximum likelihood, then you typically follow the first option, whereas, under the fully Bayesian approach from which you take posterior samples also for $\theta$, you would go for the second one. I would not expect to see big differences with regard to point estimates between the two. But if you are also interested in the variance of these estimates, then the second approach accounts for the uncertainty in estimating $\theta$ whereas the first one does not.

You could still use the second approach to get $\mathbf{b}_{\mathrm{new}}$ even you fit the model with maximum likelihood by approximating the posterior distribution of the parameters $[\theta \mid \mathcal{Y}]$ with a multivariate normal distribution with mean the maximum likelihood estimates (MLEs) and covariance matrix the variance-covariance matrix of the MLEs.