maximum likelihood method for generalized mixed effects model 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?
 A: 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).
