A Unified and Computationally Efficient Non-Gaussian Statistical Modeling Framework

Abstract: Datasets that exhibit non-Gaussian characteristics are common in many fields, while the current modeling framework and available software for non-Gaussian models is limited. We introduce Linear Latent Non-Gaussian Models (LLnGMs), a unified and computationally efficient statistical modeling framework that extends a class of latent Gaussian models to allow for latent non-Gaussian processes. The framework unifies several popular models, from simple temporal models to complex spatial-temporal and multivariate models, facilitating natural non-Gaussian extensions. Computationally efficient Bayesian inference, with theoretical guarantees, is developed based on stochastic gradient descent estimation. The R package \texttt{ngme2}, which implements the framework, is presented and demonstrated through a wide range of applications including novel non-Gaussian spatial and spatio-temporal models.

• The paper introduces Linear Latent non-Gaussian Models, a unified operator-based framework that extends latent Gaussian models with generalized hyperbolic noise and supports temporal, spatial, multivariate, and non-separable spatio-temporal applications.
• The paper develops Rao-Blackwellized stochastic-gradient MAP estimation and SGLD posterior sampling using sparse Gaussian updates and GIG Gibbs steps, with five Gibbs samples per iteration often providing stable optimization.
• The paper reports that non-Gaussian models improved predictive scores for grasshopper, kidney-function, and wind-field data, while Gaussian models performed best for transformed precipitation, showing that added complexity is most valuable when heavy tails or skewness remain.

The LLnGM framework

The paper introduces Linear Latent non-Gaussian Models (LLnGMs), a hierarchical modeling framework that extends latent Gaussian models (LGMs) by replacing the Gaussian assumption on both the latent field and the measurement noise with generalized hyperbolic (GH) distributions. The model is specified through an operator equation $$\mathbf{K}(\boldsymbol\theta)\mathbf{W} = \boldsymbol\epsilon^{\mathbf{W}}$$ for the latent process, coupled with a linear observation model $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$. The noise vectors are given normal mean-variance mixtures whose mixing variables $\mathbf{V}$ have independent generalized inverse Gaussian (GIG) components, which induces GH marginals on the noise terms. Because the Gaussian distribution is a limiting case of the GH family, setting $\mathbf{V} = \mathbf{h}$ recovers a standard LGM, so the framework embeds the Gaussian baseline as a nested special case.

A central structural property is that the conditional distributions of the mixing variables remain in the GIG family, and $\mathbf{W}\mid\mathbf{V}, \mathbf{Y}$ is multivariate Gaussian with sparse precision $$\mathbf{K}^\top\mathbf{D}_{\mathbf{V}^W}^{-1}\mathbf{K} + \mathbf{A}^\top\mathbf{D}_{\mathbf{V}^Y}^{-1}\mathbf{A}$$. This closure property under both the modeling flexibility and the computational strategy of the paper.

Unified process specification

The operator matrix $\mathbf{K}(\boldsymbol\theta)$ provides a single representation for diverse model classes. Banded difference operators yield AR($p$), random walk, and Ornstein–Uhlenbeck models; finite element discretization of the Matérn SPDE $(\kappa^2 - \Delta)^{\alpha/2}u = \mathcal{W}$ yields spatial fields; Kronecker products of component operators yield separable spatio-temporal and replicated (random effects) structures; and block-diagonal constructions with a dependence matrix $\mathbf{D}$ yield multivariate Type G Matérn models in the spirit of Bolin and Wallin's multivariate SPDE framework. Two distinct bivariate spatio-temporal compositions are distinguished: an AR(1)-in-time model driven by cross-correlated spatial fields versus a bivariate coregionalization model built from independent space-time components.

The framework also admits non-separable spatio-temporal models via discretized advection–diffusion evolution SPDEs, producing sparse block-lower-bidiagonal operators that preserve computational tractability while encoding transport and directionality that separable constructions cannot represent. The appendix further establishes that although the conditional precision $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$0 does not factorize as a Kronecker product, the unconditional covariance remains second-order separable under variance-mixture non-Gaussianity — a non-obvious property relevant to interpretation of these models.

Notably, the mixture parameters ($$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$1, $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$2, and the GIG shape parameters) may be specified as functions of covariates, enabling non-stationary skewness and heavy-tailedness. The authors state that applications of non-stationary skewness models appear largely unexplored, positioning this capability as a novel contribution of the framework.

Inference methodology

Because the marginal likelihood is intractable, the authors develop MAP estimation via stochastic gradient descent using Fisher's identity to express the gradient through expectations over the augmented posterior. Gibbs sampling exploits the conditional conjugacy: $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$3 is Gaussian with sparse precision, and each $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$4 conditionally follows a GIG distribution. A Rao-Blackwellized gradient estimator, obtained by applying Fisher's identity a second time to integrate out $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$5 analytically, reduces Monte Carlo variance. Its consistency follows from geometric ergodicity of the Gibbs sampler established in companion theoretical work (Awadelkarim et al., 8 Feb 2026); convergence holds almost surely and in $$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$6 by Birkhoff's theorem. In practice, preconditioned SGD with five Gibbs samples per iteration sufficed for stable optimization.

Convergence diagnostics combine the Gelman–Rubin statistic across parallel chains, within-chain checkpoint stabilization criteria, and a cumulative inner-product diagnostic on successive gradient estimates. Full posterior inference proceeds via stochastic gradient Langevin dynamics (SGLD): after optimization converges near the MAP, injected Gaussian noise converts the SGD trajectory into a posterior sampler, avoiding Metropolis–Hastings acceptance steps entirely. Prediction propagates posterior uncertainty by drawing latent variables from the augmented posterior and evaluating the linear predictor at new locations or times.

The estimation approach depends on the ergodicity conditions of the underlying Gibbs sampler holding; while the cited theory covers this class, the practical behavior under extreme parameter regimes is not characterized beyond empirical evidence.

Software implementation

The R package ngme2 implements the framework with an R interface over a C++ backend using Eigen, CHOLMOD, and PARDISO for sparse computations, with OpenMP-based parallel chain estimation. Model specification follows R-INLA-like syntax through f() terms within standard formula interfaces, with constructors such as noise_nig() and noise_gal() for noise components.

Empirical results across four applications

Grasshopper abundance (Montana, 1948–1990). An NIG-driven AR(1) model detected significant positive skewness ($$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$$7, 95% CI 1.75–3.10) and heavy tails ($$\hat\nu = 1.33), revealing asymmetric ecological shocks invisible to the Gaussian model. Cross-validation favored NIG on MAE, CRPS, and sCRPS (e.g., CRPS 0.964 vs. 1.032), though MSE was essentially tied — a point the authors report transparently rather than overstating.</p> <p><strong>Longitudinal kidney function</strong> (Salford cohort; subsample of 500 patients, 13,196 observations). The NIG Matérn-SPDE model outperformed the Gaussian model on all four metrics (MAE 0.131 vs. 0.142; −CRPS 0.100 vs. 0.113), consistent with prior findings on this data type. The estimated shape parameter$$\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$8 indicates pronounced heavy-tailedness in eGFR trajectories.

Climate reanalysis precipitation. Here the results run counter to a non-Gaussian preference: the Gaussian non-stationary model achieved the best scores on all metrics for the 1989 hold-out year, and the authors explicitly attribute this to the cube-root transformation having rendered the data sufficiently close to Gaussianity. This honest negative result usefully delineates when the additional complexity of non-Gaussian latent fields is unwarranted.

Bivariate wind fields (Balkan Peninsula, 1,078 space-time observations). The NIG multivariate spatio-temporal model beat the Gaussian model on MAE, CRPS, and sCRPS, with the Gaussian model retaining a marginal MSE advantage.

Comparison with Stan and TMBStan

On a simulated NIG-AR(1) model with known parameters, the comparison produced sharp contrasts:

Method	Time (s)	Convergence	Noise KLD
ngme2	20.0	Converged	0.011
Stan MCMC	37.7	Poor ($\mathbf{Y} = \mathbf{A}\mathbf{W} + \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon^{\mathbf{Y}}$9)	0.075
TMBStan	360.0	Converged ($\mathbf{V}$0)	2.026

Full HMC in Stan exhibited funnel geometry pathologies: divergent transitions, low BFMI, effective sample sizes as low as 4, and severe underestimation of the process scale ($\mathbf{V}$1 vs. true value 2.00). TMBStan resolved mixing but introduced a more serious failure mode: its Laplace approximation imposes Gaussian geometry on a genuinely non-Gaussian posterior, collapsing $\mathbf{V}$2 to approximately 0.10 and inflating the measurement error to compensate, yielding a noise KLD of 2.026 — roughly two orders of magnitude worse than ngme2's 0.011. These are strong claims about the reliability of Laplace-based marginalization for latent non-Gaussian models, and they are supported by the reported estimates. However, the comparison rests on a single simulation configuration ($\mathbf{V}$3); generalization to other dimensions, noise strengths, or model classes is not demonstrated.

Limitations and open questions

The paper concedes several boundaries of the framework. Only NIG and GAL noises admit efficient inference, since these are the sole convolution-closed GH subclasses, restricting the usable mixture families. Fractional-smoothness SPDEs (via rational approximations), fields on metric graphs, and parameterized observation matrices $\mathbf{V}$4 are identified as extensions not yet incorporated. The climate application shows the framework offers no automatic benefit when preprocessing already induces approximate Gaussianity, so the practical value hinges on genuine latent non-Gaussianity. Whether the SGLD-based uncertainty quantification is calibrated in high-dimensional settings — where only MAP point estimation was benchmarked against reference samplers — remains unverified.

Conclusion

This work consolidates previously fragmented non-Gaussian latent variable modeling — temporal, longitudinal, spatial, multivariate, separable and non-separable spatio-temporal — into a single operator-based formulation with provably convergent, Rao-Blackwellized stochastic gradient inference and full Bayesian posterior sampling via SGLD. The empirical evidence indicates that when latent non-Gaussianity is present, the framework recovers the driving-noise structure substantially better than generic HMC or Laplace-based alternatives at lower cost, while the precipitation analysis appropriately demonstrates that the Gaussian special case should not be abandoned without justification.