Maximum likelihood for generalized linear models with nested random effects via high-order, multivariate Laplace approximation

Nested random effects models are often used to represent similar processes occurring in each of many clusters. Suppose that, given cluster-specific ra ndom effects b, the data y are distributed according to f(y\b, theta), whil e b follows a density p(b\theta). Likelihood inference requires maximizatio n of integral f(y\b, theta)p(b\theta)db with respect to theta. Evaluation o f this integral often proves difficult, making likelihood inference difficu lt to obtain. We propose a multivariate Taylor series approximation of the log of the integrand that can be made as accurate as desired if the integra nd and all its partial derivatives with respect to b are continuous in the neighborhood of the posterior mode of b\theta, y. We then apply a Laplace a pproximation to the integral and maximize the approximate integrated likeli hood via Fisher scoring. We develop computational formulas that implement t his approach for two-level generalized linear models with canonical link an d multivariate normal random effects. A comparison with approximations base d on penalized quasi-likelihood, Gauss-Hermite quadrature, and adaptive Gau ss-Hermite quadrature reveals that, for the hierarchical logistic regressio n model under the simulated conditions, the sixth-order Laplace approach is remarkably accurate and computationally fast.