A Fast Monte Carlo Newton-Raphson Algorithm to Estimate Generalized Linear Mixed Models with Dense Covariance

TL;DR

This paper tackles the computational bottleneck of estimating Generalized Linear Mixed Models with dense covariance structures by introducing a fast Monte Carlo maximum likelihood (MCML) method built on a stochastic Newton-Raphson framework. Random effects are sampled via an efficient importance-sampling scheme, with fixed- and covariance-parameter updates carried out through Monte Carlo-augmented gradients and Hessians, enabling full maximum likelihood estimation without covariance approximations. The authors develop practical stopping criteria and Monte Carlo-sample-sizing rules, and demonstrate substantial speedups using GPU hardware while preserving competitive estimator performance relative to INLA across Poisson and Binomial GLMMs with spatial processes. They further validate the approach on a large real-world dataset, showing feasible runtimes and improved scalability, and discuss limitations and future directions for uncertainty quantification and more extensive benchmarking.

Abstract

Estimation of Generalised linear mixed models (GLMM) including spatial Gaussian process models is often considered computationally impractical for even moderately sized datasets. In this article, we propose a fast Monte Carlo maximum likelihood (MCML) algorithm for the estimation of GLMMs. The algorithm is a stochastic Newton-Raphson method, which approximates the expected Hessian and gradient of the log-likelihood by drawing samples of the random effects. We propose a new stopping criterion for efficient termination and preventing long runs of sampling in the stationary post-convergence phase of the algorithm and discuss Monte Carlo sample size choice. We run a series of simulation comparisons of spatial statistical models alongside the popular integrated nested Laplacian approximation method and demonstrate potential for similar or improved estimator performance and reduced running times. We also consider scaling of the algorithms to large datasets and demonstrate a greater than 100-fold reduction in running times using modern GPU hardware to illustrate the feasibility of full maximum likelihood methods with big spatial datasets.

Figures (4)

• Figure 1: Z-statistic and Bayes Factor (black lines) and running means (red dotted line) for a one-sided test of $\mu_{\Delta L}^{(t)} = 0$. The Bayes Factor uses the prior odds based on $t_0 = 30$
• Figure 2: Values of the log Bayes Factor for different values of $\pi_0$ and $p$ with thresholds for different values of the Bayes Factor
• Figure 3: Estimator distributions for $\beta_0$, $\beta_1$ and log of $\tau^2$ and $\lambda$ for the Poisson GLMM with parameter values reported in the 1st row of Table \ref{['tab:poisres']} with true values indicated by the dashed lines
• Figure 4: Predicted log relative risk and prevalence of of onchocerciasis using MCML and INLA.