reference: http://www.stat.wisc.edu/~bates/UseR2008/WorkshopD.pdf from page 120.

Now I want to fit a generalized logistic mixed model, where $\beta$ is fixed effect, $\theta$ is variance covariance parameter of random effects, $u$ is the random effects, $\tilde{u}$ is the empirical random effects. The likelihood function is $$L(\theta, \beta \mid y) = \int p(y \mid X,\beta, Z, u) \; p(u \mid \theta) \, du.$$

I do not understand how fixed effects $\beta$ and variance parameters of random effects $\theta$ are estimated. According to the slides, first they use penalized iteratively reweighted leases squares method to estimate $\tilde{u}(y|\theta,\beta)$, then use Laplace approximation and maximum likelihood method to maximum deviance $d(\beta,\theta|y)$. The deviance is a function of $y,\tilde{u}$, then a function of $y,\beta,\theta$.

From my perspective, $\tilde{u}$ is not a closed form. How can we mamimize a function where the function includes a PIRLS algorithm?


In general, classical maximum likelihood estimation of the parameters in generalized linear mixed effects models requires a combination of numerical integration and numerical optimization. The numerical integration is requires because the integral in the specification of the marginal likelihood function, as you wrote above, does not have a closed-form solution.

For this numerical integration the adaptive Gaussian quadrature is considered one of the best approaches to approximate the involved integrals. The adaptive part of the rule has to do with appropriately centering and scaling the integrands. This step requires locating the mode of the conditional likelihood $$L(u \mid y; \beta, \theta) = p(y \mid u, \beta) \; p(u \mid \theta).$$ In this step PIRLS is used. Then when you have approximated the integral, you optimize the marginal log-likelihood $\ell(\beta, \theta) = \sum_i \log L_i(\beta, \theta \mid y_i, X_i, Z_i)$ you wrote above.

The degree of the approximation in the adaptive Gaussian quadrature rule is controlled by the number of quadrature points. The Laplase approximation is equivalent to 1 quadrature point (and this is why it often does not perform that well, especially for binary data and Poisson data will low expected counts).

  • $\begingroup$ Thank you! I am still confused about one point. When PIRLS is used, $\beta$ and $\theta$ are unknown because they are what we want to estimate in the second step. But if the input is unknown, how do I use PIRLS? To be specific, how do I determine the number of steps of PIRLS? $\endgroup$ – zqin Mar 22 at 17:22
  • $\begingroup$ The procedure works iteratively, you start with initial values for $\beta$ and $\theta$ and for these you calculate the modes $\hat u$ using PIRLS. Then given these modes you do the Laplace approximation and you update the values for $\beta$ and $\theta$. And you continue alternating these two steps until convergence. $\endgroup$ – Dimitris Rizopoulos Mar 22 at 18:07
  • $\begingroup$ Thank you, professor! By the way, the procedure looks like an EM algorithm. Estimating $\tilde{u}$ and Laplace approximation, i.e. the numerical integration is the E-step. And maximizing log-likelihood is the M-step. Is this correct? $\endgroup$ – zqin Mar 22 at 18:40
  • $\begingroup$ Even though it seems like it, this is not an EM algorithm. $\endgroup$ – Dimitris Rizopoulos Mar 23 at 13:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.