Fit generalized linear mixed models (GLMMs) with normal random
effects using first-order Laplace, fully exponential Laplace (FEL) with
mean-only corrections, and FEL with mean and variance-diagonal corrections in
the E-step of an expectation-maximization (EM) algorithm. The current
development version provides a matrix-based interface (y, X, Z) and
supports binary logit and probit, and Poisson log-link models. An EM
framework is used to update fixed effects, random effects, and a single
variance component tau^2 for G = tau^2 I, with staged approximations
(Laplace -> FEL mean-only -> FEL full) for efficiency and stability. A
pseudo-likelihood engine glmmFEL_pl() implements the working-response /
working-weights linearization approach of Wolfinger and O'Connell (1993)