PaperlayerView exact revision on arXiv

Source-linked AI summary

Spike-and-Slab Priors for Function Selection in Structured Additive Regression Models

Fabian Scheipl, Ludwig Fahrmeir, Thomas Kneib

arXiv:1105.5250v2stat.MEstat.AP

TL;DR

Function selection is challenging in complex regression models with many potentially nonlinear effects. The paper proposes spike-and-slab priors for selecting individual coefficients or coefficient blocks, reporting good simulation performance and improved prediction accuracy while noting diminished performance under very strong concurvity.

  • Problem

    Automatic function selection remains difficult in complex regression models with a large number of potentially nonlinear effects.

  • Method

    The paper introduces a spike-and-slab prior framework that selects individual coefficients or batches representing model terms.

  • Results

    The approach demonstrates good performance in simulations and offers improved prediction accuracy.

  • Takeaways & Limitations

    The framework extends function-selection methods to a broader model class while allowing principled inclusion of uncertainty.

  • Takeaways & Limitations

    Very strong concurvity leads to diminished performance, and automatic function selection is not suitable in some complex settings.

Abstract

from arXiv · show

Structured additive regression provides a general framework for complex Gaussian and non-Gaussian regression models, with predictors comprising arbitrary combinations of nonlinear functions and surfaces, spatial effects, varying coefficients, random effects and further regression terms. The large flexibility of structured additive regression makes function selection a challenging and important task, aiming at (1) selecting the relevant covariates, (2) choosing an appropriate and parsimonious representation of the impact of covariates on the predictor and (3) determining the required interactions. We propose a spike-and-slab prior structure for function selection that allows to include or exclude single coefficients as well as blocks of coefficients representing specific model terms. A novel multiplicative parameter expansion is required to obtain good mixing and convergence properties in a Markov chain Monte Carlo simulation approach and is shown to induce desirable shrinkage properties. In simulation studies and with (real) benchmark classification data, we investigate sensitivity to hyperparameter settings and compare performance to competitors. The flexibility and applicability of our approach are demonstrated in an additive piecewise exponential model with time-varying effects for right-censored survival times of intensive care patients with sepsis. Geoadditive and additive mixed logit model applications are discussed in an extensive appendix.

1. INTRODUCTION

Structured additive regression accommodates diverse covariate effects, but its flexibility makes selecting relevant functions, representations, and interactions difficult. The paper develops a Bayesian spike-and-slab framework with parameter expansion to perform function selection across broad structured additive models.

  • Motivation: Structured additive regression combines additive predictors with nonlinear functions, interactions, spatial effects, and random effects represented through suitable design matrices.Effects may be represented using penalized splines, Gaussian Markov random fields, or Gaussian random effects.
  • Motivation: Function selection seeks parsimonious special cases by identifying negligible effects, selecting covariates, and distinguishing linear, nonlinear, and interaction effects.The motivating applications include selecting categorical covariates, determining linearity, and assessing interactions or spatial effects.
  • Limitations of existing methods: Existing likelihood-ratio tests are limited to Gaussian regression and are unsuitable for automatic selection in complex models with many potentially nonlinear effects.The paper motivates a Bayesian counterpart based on spike-and-slab priors for selecting function variances.
  • Proposed approach: A straightforward spike-and-slab approach creates severe convergence and mixing problems because small block variances and small coefficients reinforce one another.Blockwise samplers may remain near a zero-variance basin of attraction.
  • Proposed approach: The proposed multiplicative parameter expansion yields efficient MCMC, regularization resembling Lq-penalization with q < 1, and selection across non-Gaussian STAR models and diverse effect types.The framework covers splines, spatial and random effects, interactions, and high-dimensional predictors with hundreds of model terms.

2. NMIG PRIORS FOR FUNCTION SELECTION

The paper reparameterizes structured additive model terms into penalized and unpenalized components, enabling function selection for both parts. It then uses NMIG and parameter-expanded NMIG priors to shrink coefficients or coefficient blocks while supporting parsimonious function representations.

  • 2.1. Generic Parameterization: Spectral decomposition separates each function into penalized and unpenalized components with an orthogonal basis representation.The penalized component receives a proper Gaussian representation, while null-space coefficients represent the unpenalized component.
  • 2.1. Generic Parameterization: Reduced-rank computation retains eigenvectors explaining at least .995 of total eigenvalue variation, often reducing a 20-function basis to 8–12 penalized columns.For a cubic P-spline with 20 basis functions, the unpenalized component commonly has one linear-trend column because the constant is absorbed into the global intercept.
  • 2.1. Generic Parameterization: The reparameterization permits function selection on both penalized and unpenalized parts, distinguishing nonlinear effects from linear representations.Separating polynomial null-space terms from non-polynomial penalized terms supports decisions about whether an effect is absent, linear, or nonlinear.
  • 2.2. Parameter-Expanded NMIG Prior: The multiplicative parameterization represents a coefficient block as βj = αjξj, with αj controlling block importance and ξj distributing it across entries.This construction extends coefficient-level shrinkage to blocks of coefficients representing model terms.
  • 2.2. Parameter-Expanded NMIG Prior: NMIG combines a spike at small variance with a diffuse slab, so coefficients assigned to the spike are strongly shrunk toward zero and can be excluded.A Beta prior on the mixture weight can encode sparsity preferences for overparameterized models.
  • 2.3. Shrinkage Properties: The peNMIG prior combines an infinite spike at zero with heavy tails and generalizes shrinkage to multiple coefficient blocks.Its marginal shape is described as close to the horseshoe prior, while preserving block-level sampling and shrinkage properties.

3. EMPIRICAL EVALUATION

The empirical evaluation compares the peNMIG-based approach with existing selection and prediction methods across simulations, binary-classification benchmarks, and a survival application. Results show competitive prediction, sparse and robust selection in many settings, but sensitivity under strong concurvity and for large coefficient blocks in non-Gaussian models.

  • Simulation studies: Simulations cover varied response types, sample sizes, signal-to-noise ratios, covariate correlations, concurvity, and sparsity levels.Competitors include component-wise boosting, ACOSSO, SPAM, double-shrinkage GAMs, HGAM, and Bayesian additive regression trees.
  • Simulation studies: The proposed approach is highly competitive while supporting generalized exponential-family responses and broader model-term classes than most previous suggestions.The paper contrasts its scope with methods restricted to Gaussian responses or univariate smooth functions.
  • Simulation studies: Selection of large coefficient blocks, such as random effects, can be problematic for Poisson and binary responses, with low selection accuracy but no adverse estimation effects.This limitation concerns model selection rather than the resulting estimation performance.
  • Case Study: Survival of Surgical Patients with Severe Sepsis: The application examples and benchmark classification analyses were successfully fitted using a default prior specification.The reported default uses v0 = 0.00025, (aτ, bτ) = (5, 25), and (aw, bw) = (1, 1).
  • Binary Classification Benchmarks: Across 21 binary-classification datasets, the approach usually has more variable performance than mboost but lower median predictive deviances under all prior specifications.Prediction performance is reported as robust to different hyperparameter settings, including settings with large p/n ratios.
  • Binary Classification Benchmarks: On most benchmark datasets, the approach achieves relatively more precise predictions with smaller models than mboost, without clear dependence on basic dataset characteristics.The comparison considers predictive deviance and model sparsity, while prior differences are generally more sensitive to v0 than to the remaining hyperparameters.

4. CONCLUSIONS

The paper proposes a general Bayesian function-selection framework based on spike-and-slab priors and multiplicative parameter expansion for structured additive regression. It supports broad model classes, improves computational behavior, and shows strong empirical performance with relatively low hyperparameter sensitivity.

  • Conclusions: A non-identifiable multiplicative parameter expansion associates selection of coefficient blocks with a scalar scaling factor.The blocks can represent spline bases or random intercepts.
  • Conclusions: The reparameterization alleviates mixing problems that arise in naive implementations of the proposed prior.The paper links this construction to efficient MCMC strategies and desirable shrinkage properties.
  • Conclusions: The peNMIG prior is applicable to various response types, particularly non-Gaussian responses.The framework extends function selection to exponential-family structured additive regression models.
  • Conclusions: Good performance is demonstrated in simulations and applications, with fairly low sensitivity to hyperparameter settings and prior specifications.The paper also reports robustness in complex models and theoretical investigations of peNMIG shrinkage properties.
  • Conclusions: The model class can be extended fairly easily to other latent Gaussian or latent exponential-family models.This supports broader use of the framework beyond the models directly evaluated.
  • Conclusions: The framework extends Bayesian model averaging beyond linear and generalized linear models to exponential-family structured additive regression.It permits uncertainty about term selection and model structure to enter inferential statements.

COMPUTATIONAL DETAILS

The approach is implemented in the R package spikeSlabGAM.

  • COMPUTATIONAL DETAILS: The proposed approach is implemented in the R package spikeSlabGAM.The package is identified as Scheipl (2011b).

B. Problems of the Conventional NMIG Prior when Selecting

The conventional NMIG blockwise sampler mixes poorly when selecting long coefficient blocks, with switching probabilities becoming effectively trapped as block dimension grows. A multiplicative parameter expansion is proposed to remedy this and provide desirable shrinkage.

  • Blockwise Gibbs sampling is ill suited for NMIG posteriors even for moderately large coefficient blocks.
  • For d = 1, switching between spike and slab states remains likely across realistic coefficient magnitudes.
  • For d = 5, inclusion probabilities already resemble a step function, while for d = 20 switching away from inclusion is practically zero for most draws.
  • Mixing of γ is therefore very slow for long subvectors, producing posterior inclusion probabilities near either 0 or 1 that depend strongly on chain initialization.
  • A multiplicative parameter expansion offers a possible remedy and induces desirable shrinkage properties for the resulting estimates.

C. Simulation Results

The simulations evaluate prediction and complexity recovery across Gaussian and Poisson settings, sparsity levels, sample sizes, dependence structures, and competing methods. They use replicated data-generating scenarios with known relevant and irrelevant model terms.

  • The study compares peNMIG with component-wise boosting, ACOSSO for Gaussian responses, and oracle GAM fits using simulated Gaussian and Poisson data.
  • The Gaussian comparisons include ACOSSO, whereas ACOSSO is unavailable for non-Gaussian responses.
  • The data-generating process varies sparsity, covariate dependence, sample size, signal-to-noise ratio, and response distribution.
  • Each setting uses 50 replications, with predictive MSE evaluated on test sets containing 5000 observations.
  • Complexity recovery measures the proportion of model terms correctly included or excluded, counting both true positives and true negatives.

C.1. Gaussian response

For Gaussian responses, peNMIG delivers robust prediction across prior settings and often approaches the oracle model, while complexity recovery depends strongly on the spike variance parameter v0.

  • Predictive performance is robust to prior settings, with especially similar behavior across settings within replications.
  • Median relative prediction MSE stays below 2 for peNMIG across Gaussian scenarios, while boosting and ACOSSO exceed 4 in most settings.
  • In large-sample correlated-covariate cases, ACOSSO and boosting reach relative prediction MSEs above 32 and 64, respectively.
  • For large-sample low-sparsity scenarios, peNMIG approaches the oracle model, with relative prediction MSEs close to one.
  • Estimated inclusion probabilities are highly sensitive to v0 but comparatively robust to (aτ, bτ), with v0 = 0.00025 producing sensitivities consistently above 0.7.
  • PeNMIG has specificity above .97 across settings, supporting better accuracy than mboost in sparse scenarios despite mboost's higher sensitivity.

C.2. Poisson response

For Poisson responses, prediction is robust across prior settings and peNMIG recovers model complexity better than boosting across the evaluated settings. The choice of v0 trades sensitivity against specificity.

  • Predictive performance is very robust against different prior settings for Poisson responses.
  • PeNMIG's predictions are more precise than mboost's, especially for smaller datasets and correlated responses.
  • For low-sparsity correlated settings, peNMIG approaches the oracle GAM, with relative prediction errors mostly between 1 and 1.5 and occasionally below the oracle error.
  • Smaller v0 tends to improve unsparse-setting performance by increasing sensitivity and reducing specificity, whereas larger v0 has the opposite trade-off.
  • Complexity recovery is much better for peNMIG than for boosting across the evaluated settings and priors.
  • In the low-sparsity uncorrelated scenario, mboost includes practically all model terms, resulting in very low specificity.

C.3. Gaussian GAM with concurvity

The simulation evaluates function selection under varying concurvity, signal-to-noise ratios, and scenarios where influential and noise covariates are functionally related. The proposed approach achieves strong prediction and selection performance, while inclusion probabilities become less reliable when effects cannot be disentangled.

  • Simulation design: The simulation varies three concurvity scenarios, SNR values of 1 and 5, and 50 replications for each settings combination.Predictive MSE is evaluated on test sets with 5000 observations.
  • Prediction and selection: The proposed approach dominates prediction accuracy in difficult concurvity settings, with BART and mboost as fairly close competitors.Figure 14 compares prediction MSE across scenarios, signal-to-noise ratios, and concurvity levels.
  • Prediction and selection: The proposed approach outperforms the other methods in selection accuracy, while double shrinkage is a close second in noisy settings but performs much worse in prediction.Selection accuracy measures correctly selected or removed covariates.
  • Concurvity effects: Inclusion probabilities are unreliable for intermediate to strong concurvity when separate effects are difficult to distinguish, although interactions are often selected in such cases.For x3 and x4, inclusion probabilities decrease under intermediate concurvity and recover somewhat under perfect curvilinearity.
  • Concurvity effects: For a noise variable related to an influential covariate, the correct model was identified in 44 of 50 low-SNR and 17 of 50 high-SNR replicates, except under perfect curvilinearity.At perfect curvilinearity, the effects cannot be disentangled.
  • Concurvity effects: With strong concurvity ≥0.6 in the noisy third scenario, the true model was identified at least 24 of 50 times for low SNR and at least 41 of 50 times for high SNR.The scenario contains a noisy version of a spurious covariate with true effects.

C.4. Summary

The generalized additive model simulations show that peNMIG is competitive for estimation and function selection across complex settings. Estimation is robust to hyperparameter configurations and strong concurvity, whereas selection is more sensitive, particularly to v0.

  • Simulation summary: The peNMIG model is very competitive in estimation accuracy and robust to different hyperparameter configurations in complex models with strong concurvity.These simulations extend across generalized additive model settings.
  • Simulation summary: Function selection is more sensitive to hyperparameter configurations, especially v0.Smaller v0 improves the distinction between important and irrelevant terms.
  • Simulation summary: With smaller v0, spikeSlabGAM distinguishes important and irrelevant terms fairly reliably.
  • Simulation summary: peNMIG is competitive with component-wise boosting and clearly dominates other function-selection approaches in the concurvity simulation study.

D. Additional Case Studies

The additional case studies apply blockwise function selection to Munich rental data, combining spatial, smooth, categorical, and interaction terms. The selected models retain predictive accuracy, produce stable term selection, and avoid evident selection-induced attenuation of important effects.

  • Munich rental guide: The Munich rental model contains 594 coefficients across 269 terms, including spatial, smooth, and 265 additional potentially influential covariates.Term selection is challenging because the available covariates include redundant and highly collinear variables.
  • Munich rental guide: The selected additive predictor is dominated by balcony presence, tenancy start date, residential-area quality, attic presence, and playground presence.Terms with posterior inclusion probability greater than 10% are listed in Table 3.
  • Cross-validation: Prediction accuracy is not diminished by selecting only relevant terms, with equivalent performance for the peNMIG and BayesX expert models.The authors attribute this to the adaptive shrinkage properties of the peNMIG prior.
  • Cross-validation: Adding many interaction terms produces no noticeable prediction change in most folds because only one interaction has inclusion probability above 0.1 and its effect is small.
  • Cross-validation: The approach estimates important effects without selection-induced attenuation bias despite using variable-selection priors.Its predictive performance is similar to Bayesian lasso, Bayesian ridge, and conventional NMIG for most folds.
  • Cross-validation: Variable selection is stable across all ten folds for both expert models, while the full model identifies a core set of 23 covariates in at least nine folds.Nine additional covariates exceed the inclusion threshold in at least one fold.

D.1. Case Study: Hymenoptera Venom Allergy

The Hymenoptera venom allergy case study uses a highly complex model to assess nonlinear effects, interactions, center heterogeneity, and selection stability. The analysis finds broadly higher risk for wasp patients, stable important main effects, and convergence challenges that can be mitigated with multiple chains.

  • Model and objectives: The analysis models severe allergic reactions using nonlinear age and tryptase effects, bee-versus-wasp differences, study-center effects, and interactions.The full peNMIG model contains 267 coefficients in 66 terms.
  • Estimated effects: Posterior inclusion probabilities identify the interlocking interaction among CAP class, tryptase, and culprit insect as a joint effect for interpretation.Figure 20 displays effects for terms with P(γ = 1) > .1.
  • Estimated effects: The estimated three-way interaction has substantial uncertainty, but risk is generally higher for wasp patients, with a culprit-insect odds ratio of 1.16 (80%CI: 1-2.43).The increase in risk for wasp patients appears smaller at lower and larger at higher tryptase concentrations.
  • Convergence: Poor mixing causes some single-chain inclusion probabilities to end near zero or one, whereas many parallel chains alleviate this problem.Chains can remain trapped near posterior modes for long periods.
  • Convergence: Although posterior means converge slowly for problematic terms, 10 to 20 chains reliably distinguish important, intermediate, and negligible effects.The authors do not claim that this number fully explores the high-dimensional model space or reliably estimates joint posterior model probabilities.
  • Predictive performance and stability: More complex models show slight decreases in predictive accuracy but still outperform an unregularized GAMM.
  • Predictive performance and stability: Across subsamples, marginal term inclusion probabilities are fairly stable, and as few as 8 chains may provide rough term-importance estimates.All model specifications identify the same subset of important main effects.
Loading 1105.5250v2…