Hierarchical additive interaction modelling with Gaussian process prior and its efficient implementation for multidimensional grid data

Additive Gaussian process (GP) models offer flexible tools for modelling complex non-linear relationships and interaction effects among covariates. While most studies have focused on predictive performance, relatively little attention has been given to identifying the underlying interaction structure, which may be of scientific interest in many applications. In practice, the use of additive GP models in this context has been limited by the cubic computational cost and quadratic storage requirements of GP inference. This paper presents a fast hierarchical additive interaction GP model for multi-dimensional grid data. A hierarchical ANOVA decomposition kernel forms the foundation of our model, which incorporate main and interaction effects under the principle of marginality. Kernel centring ensures identifiability and provides a unique, interpretable decomposition of lower- and higher-order effects. For datasets forming a multi-dimensional grid, efficient implementation is achieved by exploiting the Kronecker product structure of the covariance matrix. Our contribution is to extend Kronecker-based computation to handle any interaction structure within the proposed class of hierarchical additive GP models, whereas previous methods were limited to separable or fully saturated cases. The benefits of the proposed approach are demonstrated through simulation studies and an application to high-frequency nitrogen dioxide concentration data in London.