  There are proposals that extend the classical generalized additive models
(GAMs) to accommodate high-dimensional data ($p>>n$) using group sparse
regularization. However, the sparse regularization may induce excess shrinkage
when estimating smooth functions, damaging predictive performance. Moreover,
most of these GAMs consider an "all-in-all-out" approach for functional
selection, rendering them difficult to answer if nonlinear effects are
necessary. While some Bayesian models can address these shortcomings, using
Markov chain Monte Carlo algorithms for model fitting creates a new challenge,
scalability. Hence, we propose Bayesian hierarchical generalized additive
models as a solution: we consider the smoothing penalty for proper shrinkage of
curve interpolation via reparameterization. A novel two-part spike-and-slab
LASSO prior for smooth functions is developed to address the sparsity of
signals while providing extra flexibility to select the linear or nonlinear
components of smooth functions. A scalable and deterministic algorithm,
EM-Coordinate Descent, is implemented in an open-source R package BHAM.
Simulation studies and metabolomics data analyses demonstrate improved
predictive and computational performance against state-of-the-art models.
Functional selection performance suggests trade-offs exist regarding the effect
hierarchy assumption.
