In generalized mixed-effects model, after fixed effects and variance covariance matrix are fitted, how are empirical random effects calculated?

Question

For example, I would like to fit a logistic mixed-effects model.

This article fitting glmm talks about how to fit fixed effects as well as variance covariance matrix of random effects. Theoretically we can use maximum likelihood estimator like this.

But I do not know how to compute the empirical random effects for all groups. In linear mixed-effects model, we can calculate EBLUPS(empirical best linear unbiased predictors) of random effects easily. In generalized linear model, the only method I know is to calculate the posterior distribution of random effects for each group, i.e. $(b_i|x_{it},\beta, y_{it})$, and then calculate the maximum likelihood estimator for each group.

Is there any better method?

In addition, why don't people maximize the joint distribution of $(y_{it}, \beta, b_i, \Sigma_b)$? In this way we can get all the estimates together.

Answer 1

Indeed, when fitting mixed models with maximum likelihood, estimates for the random effects are obtained in a second step using emprical Bayes methodology. In your example, suppose that $\hat \theta = (\hat \beta^\top, \mbox{vech}(\hat \Sigma_b)^\top)^\top$ denote the maximum likelihood estimates for the fixed effects $\beta$ and the unique element of the covariance matrix $\Sigma_b$. Then, the empirical Bayes estimates of the random effects are the modes of the posterior distribution $$p(b_i \mid y_i; \hat\theta) = \frac{p(y_i \mid b_i; \hat \beta) \; p(b_i; \hat \Sigma_b)}{p(y_i; \hat\theta)}, $$ which are typically found by numerically maximizing the log of the numerator, i.e., $$\hat b_i = \arg\max_b \{\log p(y_i \mid b; \hat \beta) + \log p(b; \hat \Sigma_b)\}.$$

A reason why in classical maximum likelihood you do not maximize the log likelihood $\log p(y, b, \beta, \Sigma_b)$ w.r.t. all $(b, \beta, \Sigma_b)$ is because you have the requirement that the dimension of the parameter space is bounded and does not increase with the sample size $n$. If you would treat the random effects $b_i$ as parameters, you have that as $n$ increases, the dimensionality of the parameter space would also increase. In the literature there have been some alternative estimation approaches proposed that do that, such as H-likelihood; however, they have not (yet) displaced the classical maximum likelihood as the primary estimation approach for mixed models.