Sparse additive Gaussian process with soft interactions

by   Garret Vo, et al.

Additive nonparametric regression models provide an attractive tool for variable selection in high dimensions when the relationship between the response and predictors is complex. They offer greater flexibility compared to parametric non-linear regression models and better interpretability and scalability than the non-parametric regression models. However, achieving sparsity simultaneously in the number of nonparametric components as well as in the variables within each nonparametric component poses a stiff computational challenge. In this article, we develop a novel Bayesian additive regression model using a combination of hard and soft shrinkages to separately control the number of additive components and the variables within each component. An efficient algorithm is developed to select the importance variables and estimate the interaction network. Excellent performance is obtained in simulated and real data examples.



There are no comments yet.


page 14

page 16

page 17

page 24

page 27


Novel Semi-parametric Tobit Additive Regression Models

Regression method has been widely used to explore relationship between d...

Bayesian Variable Selection for Gaussian copula regression models

We develop a novel Bayesian method to select important predictors in reg...

Additivity Assessment in Nonparametric Models Using Ratio of Pseudo Marginal Likelihoods

Nonparametric regression models such as Bayesian Additive Regression Tre...

Causal Inference with the Instrumental Variable Approach and Bayesian Nonparametric Machine Learning

We provide a new flexible framework for inference with the instrumental ...

Flexible, Non-parametric Modeling Using Regularized Neural Networks

Non-parametric regression, such as generalized additive models (GAMs), i...

An Optimal Test for the Additive Model with Discrete or Categorical Predictors

In multivariate nonparametric regression the additive models are very us...

Variable selection in Functional Additive Regression Models

This paper considers the problem of variable selection when some of the ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Variable selection has played a pivotal role in scientific and engineering applications, such as biochemical analysis (Manavalan and Johnson, 1987), bioinformatics (Saeys et al, 2007), and text mining (Kwon et al, 2003)

, among other areas. A significant portion of existing variable selection methods are only applicable to linear parametric models. Despite the linearity and additivity assumption, variable selection in linear regression models has been popular since 1970; refer to Akaike information criterion [AIC;

Akaike (1973)]; Bayesian information criterion [BIC; Schwarz et al (1978)] and Risk inflation criterion [RIC; Foster and George (1994)].

Popular classical sparse-regression methods such as Least absolute shrinkage operator [LASSO; Tibshirani (1996); Efron et al (2004)], and related penalization methods (Fan and Li, 2001; Zou and Hastie, 2005; Zou, 2006; Zhang, 2010) have gained popularity over the last decade due to their simplicity, computational scalability and efficiency in prediction when the underlying relation between the response and the predictors can be adequately described by parametric models. Bayesian methods (Mitchell and Beauchamp, 1988; George and McCulloch, 1993, 1997) with sparsity inducing priors offers greater applicability beyond parametric models and are a convenient alternative when the underlying goal is in inference and uncertainty quantification. However, there is still a limited amount of literature which seriously considers relaxing the linearity assumption, particularly when the dimension of the predictors is high. Moreover, when the focus is on learning the interactions between the variables, parametric models are often restrictive since they require very many parameters to capture the higher-order interaction terms.
      Smoothing based non-additive nonparametric regression methods (Lafferty and Wasserman, 2008; Wahba, 1990; Green and Silverman, 1993; Hastie and Tibshirani, 1990) can accommodate a wide range of relationships between predictors and response leading to excellent predictive performance. Such methods have been adapted for different methods of functional component selection with non-linear interaction terms: component selection and smoothing operator [COSSO; Lin et al (2006)], sparse addictive model [SAMS; Ravikumar et al (2009)] and variable selection using adaptive nonlinear interaction structure in high dimensions [VANISH; Radchenko and James (2010)]. However, when the number of variables is large and their interaction network is complex, modeling each functional component is highly expensive.
      Nonparametric variable selection based on kernel methods are increasingly becoming popular over the last few years. Liu et al (2007)

provided a connection between the least square kernel machine (LKM) and the linear mixed models.

Zou et al (2010); Savitsky et al (2011) introduced Gaussian process with dimension-specific scalings for simultaneous variable selection and prediction. Yang et al (2015) argued that a single Gaussian process with variable bandwidths can achieve optimal rate in estimation when the true number of covariates . However, when the true number of covariates is relatively high, the suitability of using a single Gaussian process is questionable. Moreover, such an approach is not convenient to recover the interaction among variables. Fang et al (2012) used the nonnegative Garotte kernel to select variables and capture interaction. Though these methods can successfully perform variable selection and capture the interaction, non-additive nonparametric models are not sufficiently scalable when the dimension of the relevant predictors is even moderately high. Fang et al (2012) claimed that extensions to additive models may cause over-fitting issues in capturing the interaction between variables (i.e. capture more interacting variables than the ones which are influential).
      To circumvent this bottleneck, Yang et al (2015); Qamar and Tokdar (2014)

introduced the additive Gaussian process with sparsity inducing priors for both the number of components and variables within each component. The additive Gaussian process captures interaction among variables, can scale up to moderately high dimensions and is suitable for low sparse regression functions. However, the use of two component sparsity inducing prior forced them to develop a tedious Markov chain Monte Carlo algorithm to sample from the posterior distribution. In this paper, we propose a novel method, called the additive Gaussian process with soft interactions to overcome this limitation. More specifically, we decompose the unknown regression function

into components, such as , for hard shrinkage parameters , . Each component is assigned a Gaussian process prior. To induce sparsity within each Gaussian process, we introduce an additional level of soft shrinkage parameters. The combination of hard and soft shrinkage priors makes our approach very straightforward to implement and computationally efficient, while retaining all the advantages of the additive Gaussian process proposed by Qamar and Tokdar (2014). We propose a combination of Markov chain Monte Carlo (MCMC) and the Least Angle Regression algorithm (LARS) to select the Gaussian process components and variables within each component.
      The rest of the paper is organized as follows. § 2 presents the additive Gaussian process model. § 3 describes the two-level regularization and the prior specifications. The posterior computation is detailed in § 4 and the variable selection and interaction recovery approach is presented in § 5. The simulation study results are presented in § 6. A couple of real data examples are considered in § 7. We conclude with a discussion in § 8.

2 Additive Gaussian Process

For observed predictor-response pairs , where (i.e. is the sample size and is the dimension of the predictors), an additive nonparametric regression model can be expressed as


The regression function in (1) is a sum of regression functions, with the relative importance of each function controlled by the set of non-negative parameters . Typically the unknown parameter is assumed to be sparse to prevent from over-fitting the data.
      Gaussian process (GP) (Rasmussen, 2006) provides a flexible prior for each of the component functions in . GP defines a prior on the space of all continuous functions, denoted for a fixed function and a positive definite function defined on such that for any finite collection of points , the distribution of is multivariate Gaussian with mean

and variance-covariance matrix

. The choice of the covariance kernel is crucial to ensure the sample path realizations of the Gaussian process are appropriately smooth. A squared exponential covariance kernel

with an Gamma hyperprior assigned to the inverse-bandwidth parameter

ensures optimal estimation of an isotropic regression function (van der Vaart and van Zanten, 2009) even when a single component function is used (). When the dimension of the covariates is high, it is natural to assume that the underlying regression function is not isotropic. In that case, Bhattacharya et al (2014a) showed that a single bandwidth parameter might be inadequate and dimension specific scalings with appropriate shrinkage priors are required to ensure that the posterior distribution can adapt to the unknown dimensionality. However, Yang et al (2015) showed that single Gaussian process might not be appropriate to capture interacting variables and also does not scale well with the true dimension of the predictor space. In that case, an additive Gaussian process is a more effective alternative which also leads to interaction recovery as a bi-product. In this article, we work with the additive representation in (1) with dimension specific scalings (inverse-bandwidth parameters) along dimension for the th Gaussian process component, and .

      We assume that the response vector

in (1) is centered and scaled. Let with


In the next section, we discuss appropriate regularization on and . A shrinkage prior on the facilitates the selection of variables within component and allows adaptive local smoothing. An appropriate regularization on allows to adapt to the degree of additivity in the data without over-fitting.

3 Regularization

A full Bayesian specification will require placing prior distribution on both and . However, such a specification requires tedious posterior sampling algorithms to sample from the posterior distribution as seen in Qamar and Tokdar (2014). Moreover, it is difficult to identify the role of and since one can remove the effect of the th component by either setting to zero or by having . This ambiguous representation causes mixing issues in a full-blown MCMC. To facilitate computation, we adopt a partial Bayesian approach to regularize and . We propose a hybrid-algorithm which is a combination of i) MCMC, to sample conditional on ii) and optimization to estimate conditional on . With this viewpoint, we propose the following regularization on and .

3.1 regularization for

Conditional on , (1) is linear in . Hence we impose regularization on which are updated using least absolute shrinkage and selection operator (LASSO) (Tibshirani, 1996; Tibshirani et al, 1997; Hastie et al, 2005). This enforces sparsity on at each stage of the algorithm, thereby pruning the un-necessary Gaussian process components in .

3.2 Global-local shrinkage for

The parameters controls the effective number of variables within each component. For each , are assumed to be sparse. As opposed to the two component mixture prior on in Qamar and Tokdar (2014), we enforce weak-sparsity using a global-local continuous shrinkage prior which potentially have substantial computational advantages over mixture priors. A rich variety of continuous shrinkage priors being proposed recently (Park and Casella, 2008; Tipping, 2001; Griffin and Brown, 2010; Carvalho et al, 2010, 2009; Bhattacharya et al, 2014b), which can be unified through a global-local (GL) scale mixture representation of Polson and Scott (2010) below,


for each fixed , where and are densities on the positive real line. In (3), controls global shrinkage towards the origin while the local parameters allow local deviations in the degree of shrinkage for each predictor. Special cases include Bayesian lasso (Park and Casella, 2008), relevance vector machine (Tipping, 2001), normal-gamma mixtures (Griffin and Brown, 2010) and the horseshoe (Carvalho et al, 2010, 2009) among others. Motivated by the remarkable performance of horseshoe, we assume both and

to be square-root of half-Cauchy distributions.

4 Hybrid algorithm for prediction, selection and interaction recovery

In this section, we develop a fast algorithm which is a combination of optimization and conditional MCMC to estimate the parameters and for . Conditional on , (1) is linear in and hence we resort to the least angle regression procedure (Efron et al, 2004) with five fold cross validation to estimate . The computation of the lasso solutions is a quadratic programming problem, and can be tackled by standard numerical analysis algorithms.
      Next, we describe the conditional MCMC to sample from and at a new point conditional on the parameters . For two collection of vectors and of size and respectively, denote by the matrix . Let and define and

denote the corresponding matrices. For a random variable

, we denote by the conditional distribution of given the remaining random variables.
      Observe that the algorithm does not necessarily produces samples which are approximately distributed as the true posterior distribution. The combination of optimization and conditional sampling is similar to stochastic EM (Diebolt et al, 1994; Meng and Rubin, 1994) which is employed to avoid computing costly integrals required to find maximum likelihood in latent variable models. One can expect convergence of our algorithm to the true parameters by appealing to the theory of consistency and asymptotic normality of stochastic EM algorithms (Nielsen, 2000).

4.1 MCMC to sample

Conditional on , we cycle through the following steps:

  1. Compute Compute the predictive mean

  2. Compute the predictive variance

  3. Sample .

  4. Compute the predictive

  5. Update by sampling from the distribution

  6. Update by sampling from the distribution


4.1.1 Algorithm to Sample and

In the MCMC algorithm above, the conditional distributions of and are not available in closed form. Therefore, we sample them using Metropolis-Hastings algorithm (Hastings, 1970). In this paper, we give the algorithm for updating only, as the steps for are similar. Assuming that the chain is currently at the iteration , the Metropolis-Hastings algorithm to sample independently for proceeds as following:

  1. Propose .

  2. Compute the Metropolis ratio:

  3. Sample . If then , else .

The proposal variance

is tuned to ensure that the acceptance probability is between 20% - 40%. We also propose a similar Metropolis-Hastings algorithm to sample from the conditional distribution of


5 Variable Selection and Interaction Recovery

In this section, we first state a generic algorithm to select important variables based on the samples of a parameter vector . This algorithm is independent of the prior for and unlike other variable selection algorithms, it requires little tuning parameters making it suitable for practical purposes. The idea is based on finding the most probable set of variables in the posterior median of on . Since the distribution for the number of important variables is more stable and largely unaffected by the Metropolis-Hastings algorithm, we find the mode H of the distribution for the number of important variables. Then, we select the H largest coefficients from the posterior mean of .
      In this algorithm, we use -means algorithm (Bishop, 2006; Han et al, 2011) with at each iteration to form two clusters, corresponding to signal and noise variables respectively. One cluster contains values concentrating around zero, corresponding to the noise variables. The other cluster contains values concentrating away from zeros, corresponding to the signals. At the iteration, the number of non-zero signals is estimated by the smaller cluster size out of the two clusters. We take the mode over all the iterations to obtain the final estimate H for the number of non-zero signals i.e. . The H largest entries of the posterior median of are identified as the non-zero signals.
      We run the algorithm for 5,000 iterations with a burn-in of 2,000 to ensure convergence. Based on the remaining iterates, we apply the algorithm to for each component to select important variables within each for . Using this approach, we can select the important variables within each function. We define the inclusion probability of a variable as the proportion of functions (out of ) which contains that variable. Next, we apply the algorithm to and select the important functions. Let us denote by

the set of active functions. The

probability of interaction between a pair of variables is defined as the proportion of functions within in which the pair appears together. As a result, we can find the interaction between important variables with optimal number of active components. Observe that the inclusion probability and the probability of interaction is not a functional of the posterior distribution and is purely a property of the additive representation. Hence, we do not require the sampler to convergence to the posterior distribution. As illustrated in § 6, these inclusion and the interaction probabilities provide an excellent representation of a variable or an interaction being present or absent in the model.

6 Simulation examples

In this section, we consider eight different simulation settings with 50 replicated datasets each and test the performance of our algorithm with respect to variable selection, interaction recovery, and prediction. To generate the simulated data, we draw , and , where , and . Table 1 summarize the results for the eight different datasets with different combinations of and for both non-interaction and interaction cases, respectively.

Equation for the Dataset
Simulated Dataset Non-interaction Data Interaction Data
1 100 10
2 100 100
3 100 20
4 100 100
Table 1: Summary of interaction simulated datasets.

6.1 Variable Selection

We compute the inclusion probability for each variable in each simulated dataset, then provide the bar plots as in Figures 1-4 below.

(a) Non-interaction Case
(b) Interaction Case
Figure 1: Inclusion Probability for dataset 1.
(a) Non-interaction Case
(b) Interaction Case
Figure 2: Inclusion Probability for dataset 2.
(a) Non-interaction Case
(b) Interaction Case
Figure 3: Inclusion Probability for dataset 3.
(a) Non-interaction Case
(b) Interaction Case
Figure 4: Inclusion Probability for dataset 4.

It is clear from the figures, that the estimated inclusion probability of the active variables are larger than the noise variables. To quantify the performance, we use a threshold of greater than to identify a variable as signal. With these included variables, we compute the false positive rate (FPR), which is the proportion of true signals not detected by our algorithm, and false negative rate (FNR), which is the proportion of false signals detected by our algorithm. Both values are recorded in Table 2 to assess the quantitative performance of our algorithm.

Non-interaction Dataset Interaction Dataset
1 0.0 0.0 0.0 0.0
2 0.0 0.0 0.0 0.0
3 0.0 0.0 0.0 0.05
4 0.0 0.01 0.0 0.01
Table 2: The average false positive (FPR) and false negative (FNR) for replicated datasets

Based on the results in Table 2, it is immediate that the algorithm is very successful in delivering accurate variable selection for both non-interaction and interaction cases.

6.2 Interaction Recovery

In order to capture the interaction network, we compute the probability of interaction between two variables by calculating the proportion of functions in which both the variables jointly appear. With these probability values, we provide the interaction heat map for each dataset for both the non-interaction and interaction cases.

(a) Non-interaction Case
(b) Interaction Case
Figure 5: Interaction heat map for dataset 1.
(a) Non-interaction Case
(b) Interaction Case
Figure 6: Interaction heat map for dataset 2.
(a) Non-interaction Case
(b) Interaction Case
Figure 7: Interaction heat map for dataset 3.
(a) Non-interaction Case
(b) Interaction Case
Figure 8: Interaction heat map for dataset 4.

Based on these interaction heat maps, it is evident that the estimated interaction probabilities for the non-interacting variables are less than the corresponding number for interacting variables. Using a threshold value of 0.5, we discard these non-interacting probability values to construct the interaction network among the variables. The interaction networks for each dataset in both non-interaction and interaction are as follows.

(a) Non-interaction 1
(b) Interaction 1
(c) Non-interaction 2
(d) Interaction 2
Figure 9: Interaction network for dataset 1 and 2, respectively.
(a) Non-interaction 3
(b) Interaction 3
(c) Non-interaction 4
(d) Interaction 4
Figure 10: Interaction network for dataset 3 and 4, respectively.

Based on the interaction network above, our algorithm successfully captures the interaction network in all the datasets for selected variables according to the inclusion probability.

6.3 Predictive performance

We randomly partition each dataset into training (50%) and test (50%) observations. We apply our algorithm on the training data and compare the performance on the test dataset. For the sake of brevity we plot the predicted vs. the observed test observations only for a few cases in Figure 11.

(a) Prediction for Non-interaction 1
(b) Prediction for Non-interaction 3
(c) Prediction for Interaction 1
(d) Prediction for Interaction 3
Figure 11: Prediction versus Response for Simulated Data

From Figure 11, the predicted observations and the true observations fall very closely along the line demonstrating a good predictive performance. We compare our results with Fang et al (2012). However, their additive model was not able to capture higher order interaction and thus have a poor predictive performance compared to our method.

6.4 Comparison with BART

Bayesian Additive Regression Tree [BART; Chipman et al (2010)] is a state of the art method for variable selection in nonparametric regression problems. BART is a Bayesian “sum of tree” framework which fits and infers the data through an iterative back-fitting MCMC algorithm to generate samples from a posterior. Each tree in BART Chipman et al (2010) is constrained by a regularization prior. Hence BART is similar to our method which also resorts to back-fitting MCMC to generate samples from a posterior.
      Since BART is well-known to deliver excellent prediction results, its performance in terms of variable selection and interaction recovery in high-dimensional setting is worth investigating. In this section, we compare our method with BART in all the three aspects: variable selection, interaction recovery and predictive performance. For comparison, with BART, we use the same simulation settings as in Table 1 with all combinations of (, ), where and .
      We use 50 replicated datasets and compute average inclusion probabilities for each variable. Similar to § 6.1, the variable must have average probability value bigger than 0.1 in order to be selected. Then, we compute the false positive and false negative rates for both algorithms as in Table 2. These values are recorded in Table 3.

Our Algorithm BART
Dataset p n Non-interaction Interaction Non-interaction Interaction
1 10 100 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 100 100 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
3 20 100 0.0 0.05 0.0 0.05 1.0 1.0 1.0 1.0
4 100 100 0.0 0.01 0.0 0.01 1.0 1.0 1.0 1.0
1 150 100 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
4 150 100 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
1 200 100 0.01 0.0 0.0 0.0 NA NA NA NA
4 200 100 0.0 0.0 0.0 0.0 NA NA NA NA
Table 3: Comparison between our algorithm and BART for variable selection

In Table 3, the first column indicates which equations are used to generate the data with the respective and values in the second and third column for both non-interaction and interaction cases. For example, if the dataset is 1, the equations to generate the data is and for non-interaction and interaction case, respectively. NA value means that the algorithm cannot run at all for that particular combination of and values.
      According to Table 3, BART performs similar to our algorithm when and . However, as increases, BART fails to perform adequately while our algorithm still performs well even when is less than . When is twice as , BART fails to run while our algorithm provides excellent results in variable selection. Overall, our algorithm performs significantly better than BART in terms of variable selection.

7 Real data analysis

In this section, we demonstrate the performance of our method in two real data sets. We use the Boston housing data and concrete slump test datasets obtained from UCI machine learning repository. Both data have been used extensively in the literature.

7.1 Boston Housing Data

In this section, we use the Boston housing data to compare the performance between BART and our algorithm. The Boston housing data (Harrison and Rubinfeld, 1978)

contains information collected by the United States Census Service on the median value of owner occupied homes in Boston, Massachusetts. The data has 506 number of instances with thirteen continuous variables and one binary variable. The data is split into 451 training and 51 test observations. The description for each variable is summarized in Table


Variables Abbreviation Description
1 CRIM Per capita crime rate
2 ZN Proportion of residential land zoned for lots over 25,000 squared feet
3 INDUS Proportion of non-retail business acres per town

Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)

5 NOX Nitric oxides concentration (parts per 10 million)
6 RM Average number of rooms per dwelling
7 AGE Proportion of owner-occupied units built prior to 1940
8 DIS Weighted distances to five Boston employment centers
9 RAD Index of accessibility to radial highways
10 TAX Full-value property-tax rate per $10,000
11 PTRATIO Pupil-teacher ratio by town
12 B where Bk is the proportion of blacks by town
13 LSTAT Percentage of lower status of the population
14 MEDV Median value of owner-occupied homes in $1000’s
Table 4: Boston housing dataset variable

MEDV is chosen as the response and the remaining variables are included as predictors. We ran our algorithm for 5000 iterations and the prediction result for both algorithms is shown in Figure 12.

(a) Our Algorithm
(b) BART
Figure 12: Prediction versus Response’s for Boston Housing Dataset

Although our algorithm has a comparable prediction error with BART, we argue below that we have a more convincing result in terms of variable selection. We displayed the inclusion probability barplot in Figure 13.

(a) Our Algorithm
(b) BART
Figure 13: Inclusion Probability for the Boston Housing dataset

Savitsky et al (2011) previously analyzed this dataset and selected variables RM, DIS and LSTAT. BART only selected NOX and RM, while oue algorithm selected CRIM, ZN, NOX, DIS, B and LSTAT based on inclusion probabilities greater than or equal to 0.1. Clearly the set of selected variables from our method has more common elements with that of Savitsky et al (2011).

7.2 Concrete Slump Test

In this section we consider an engineering application to compare our algorithm against BART. The concrete slump test dataset records the test results of two executed tests on concrete to study its behavior (Yeh et al, 2008; Yeh, 2007).
      The first test is the concrete-slump test, which measures concrete’s plasticity. Since concrete is a composite material with mixture of water, sand, rocks and cement, the first test determines whether the change in ingredients of concrete is consistent. The first test records the change in the slump height and the flow of water. If there is a change in a slump height, the flow must be adjusted to keep the ingredients in concrete homogeneous to satisfy the structure ingenuity. The second test is the “Compressive Strength Test”, which measures the capacity of a concrete to withstand axially directed pushing forces. The second test records the compressive pressure on the concrete.
      The concrete slump test dataset has 103 instances. The data is split into 53 instances for training and 50 instances for testing. There are seven continuous input variables, which are seven ingredients to make concrete, and three outputs, which are slump height, flow height and compressive pressure. Here we only consider the slump height as the output. The description for each variable and output is summarized in table 5.

Variables Ingredients Unit
1 Cement kg
2 Slag kg
3 Fly ash kg
4 Water kg
5 Super-plasticizer (SP) kg
6 Coarse Aggregation kg
7 Fine Aggregation kg
8 Slump cm
9 Flow cm
10 28-day Compressive Strength Mpa
Table 5: Concrete Slump Test dataset

The predictive performance is illustrated in figure 14.

(a) Our Algorithm
(b) BART
Figure 14: Prediction versus Response’s for Concrete Slump Test Dataset

Similar to the Boston housing dataset, our algorithm performs closely to BART in prediction. Next, we investigate the performances in terms of variable selection. We plot the barplot of the inclusion probability for each variable in Figure 15.

(a) Our Algorithm
(b) BART
Figure 15: Inclusion Probability for the Concrete Slump Test dataset

Yurugi et al (2000) determined that coarse aggregation has a significant impact on the plasticity of a concrete. Since the difference in slump’s height is to measure the plasticity of a concrete, coarse aggregation is a critical variable in the concrete slump test. According to Figure 15, our algorithm selects coarse aggregation as the most important variable unlike BART, which clearly demonstrates the efficacy of our algorithm compared to BART.

8 Conclusion

In this paper, we propose a novel Bayesian nonparametric approach for variable selection and interaction recovery with excellent performance in selection and interaction recovery in both simulated and real datasets. Our method obviates the computation bottleneck in recent unpublished work Qamar and Tokdar (2014) by proposing a simpler regularization involving a combination of hard and soft shrinkage parameters. Moreover, our algorithm is computationally efficient and highly scalable.

Although such sparse additive models are well known to adapt to the underlying true dimension of the covariates (Yang et al, 2015), literature on consistent selection and interaction recovery in the context of nonparametric regression models is missing. As a future work, we propose to investigate consistency of the variable selection and interaction of our method.


  • Akaike (1973) Akaike H (1973) Maximum likelihood identification of gaussian autoregressive moving average models. Biometrika 60(2):255–265
  • Bhattacharya et al (2014a) Bhattacharya A, Pati D, Dunson D (2014a) Anisotropic function estimation using multi-bandwidth Gaussian processes. The Annals of Statistics 42(1):352–381
  • Bhattacharya et al (2014b) Bhattacharya A, Pati D, Pillai NS, Dunson DB (2014b) Dirichlet-laplace priors for optimal shrinkage. Journal of the American Statistical Association (just-accepted):00–00
  • Bishop (2006)

    Bishop CM (2006) Pattern recognition and machine learning. springer

  • Carvalho et al (2009)

    Carvalho C, Polson N, Scott J (2009) Handling sparsity via the horseshoe. In: International Conference on Artificial Intelligence and Statistics, pp 73–80

  • Carvalho et al (2010) Carvalho CM, Polson NG, Scott JG (2010) The horseshoe estimator for sparse signals. Biometrika p asq017
  • Chipman et al (2010) Chipman HA, George EI, McCulloch RE (2010) Bart: Bayesian additive regression trees. The Annals of Applied Statistics pp 266–298
  • Diebolt et al (1994) Diebolt J, Ip E, Olkin I (1994) A stochastic em algorithm for approximating the maximum likelihood estimate. Tech. rep., Technical Report 301, Department of Statistics, Stanford University, California
  • Efron et al (2004) Efron B, Hastie T, Johnstone I, Tibshirani R (2004) Least angle regression. The Annals of statistics 32(2):407–499
  • Fan and Li (2001) Fan J, Li R (2001) Variable selection via nonconcave penalized likelihood and its oracle properties. JASA 96(456):1348–1360, DOI 10.1198/016214501753382273,
  • Fang et al (2012) Fang Z, Kim I, Schaumont P (2012) Flexible variable selection for recovering sparsity in nonadditive nonparametric models. arXiv preprint arXiv:12062696
  • Foster and George (1994) Foster DP, George EI (1994) The risk inflation criterion for multiple regression. The Annals of Statistics pp 1947–1975
  • George and McCulloch (1993) George EI, McCulloch RE (1993) Variable selection via gibbs sampling. Journal of the American Statistical Association 88(423):881–889
  • George and McCulloch (1997) George EI, McCulloch RE (1997) Approaches for Bayesian variable selection. Statistica sinica 7(2):339–373
  • Green and Silverman (1993) Green PJ, Silverman BW (1993) Nonparametric regression and generalized linear models: a roughness penalty approach. CRC Press
  • Griffin and Brown (2010) Griffin J, Brown P (2010) Inference with normal-gamma prior distributions in regression problems. Bayesian Analysis 5(1):171–188
  • Han et al (2011) Han J, Kamber M, Pei J (2011) Data mining: concepts and techniques: concepts and techniques. Elsevier
  • Harrison and Rubinfeld (1978) Harrison D, Rubinfeld DL (1978) Hedonic housing prices and the demand for clean air. Journal of environmental economics and management 5(1):81–102
  • Hastie et al (2005) Hastie T, Tibshirani R, Friedman J, Franklin J (2005) The elements of statistical learning: data mining, inference and prediction. The Mathematical Intelligencer 27(2):83–85
  • Hastie and Tibshirani (1990) Hastie TJ, Tibshirani RJ (1990) Generalized additive models, vol 43. CRC Press
  • Hastings (1970) Hastings WK (1970) Monte carlo sampling methods using markov chains and their applications. Biometrika 57(1):97–109
  • Kwon et al (2003) Kwon OW, Chan K, Hao J, Lee TW (2003) Emotion recognition by speech signals. In: INTERSPEECH, Citeseer
  • Lafferty and Wasserman (2008) Lafferty J, Wasserman L (2008) Rodeo: sparse, greedy nonparametric regression. The Annals of Statistics pp 28–63
  • Lin et al (2006) Lin Y, Zhang HH, et al (2006) Component selection and smoothing in multivariate nonparametric regression. The Annals of Statistics 34(5):2272–2297
  • Liu et al (2007) Liu D, Lin X, Ghosh D (2007) Semiparametric regression of multidimensional genetic pathway data: Least-squares kernel machines and linear mixed models. Biometrics 63(4):1079–1088
  • Manavalan and Johnson (1987) Manavalan P, Johnson WC (1987) Variable selection method improves the prediction of protein secondary structure from circular dichroism spectra. Analytical biochemistry 167(1):76–85
  • Meng and Rubin (1994) Meng XL, Rubin DB (1994) On the global and componentwise rates of convergence of the em algorithm. Linear Algebra and its Applications 199:413–425
  • Mitchell and Beauchamp (1988) Mitchell TJ, Beauchamp JJ (1988) Bayesian variable selection in linear regression. Journal of the American Statistical Association 83(404):1023–1032
  • Nielsen (2000) Nielsen SF (2000) The stochastic em algorithm: estimation and asymptotic results. Bernoulli pp 457–489
  • Park and Casella (2008) Park T, Casella G (2008) The Bayesian lasso. Journal of the American Statistical Association 103(482):681–686
  • Polson and Scott (2010)

    Polson N, Scott J (2010) Shrink globally, act locally: sparse Bayesian regularization and prediction. Bayesian Statistics 9:501–538

  • Qamar and Tokdar (2014) Qamar S, Tokdar ST (2014) Additive gaussian process regression. arXiv preprint arXiv:14117009
  • Radchenko and James (2010) Radchenko P, James GM (2010) Variable selection using adaptive nonlinear interaction structures in high dimensions. Journal of the American Statistical Association 105(492):1541–1553
  • Rasmussen (2006) Rasmussen CE (2006) Gaussian processes for machine learning
  • Ravikumar et al (2009) Ravikumar P, Lafferty J, Liu H, Wasserman L (2009) Sparse additive models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71(5):1009–1030
  • Saeys et al (2007)

    Saeys Y, Inza I, Larrañaga P (2007) A review of feature selection techniques in bioinformatics. bioinformatics 23(19):2507–2517

  • Savitsky et al (2011) Savitsky T, Vannucci M, Sha N (2011) Variable selection for nonparametric gaussian process priors: Models and computational strategies. Statistical science: a review journal of the Institute of Mathematical Statistics 26(1):130
  • Schwarz et al (1978) Schwarz G, et al (1978) Estimating the dimension of a model. The annals of statistics 6(2):461–464
  • Tibshirani (1996) Tibshirani R (1996) Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B (Methodological) pp 267–288
  • Tibshirani et al (1997) Tibshirani R, et al (1997) The lasso method for variable selection in the cox model. Statistics in medicine 16(4):385–395
  • Tipping (2001) Tipping M (2001) Sparse Bayesian learning and the relevance vector machine. The Journal of Machine Learning Research 1:211–244
  • van der Vaart and van Zanten (2009) van der Vaart AW, van Zanten JH (2009) Adaptive bayesian estimation using a gaussian random field with inverse gamma bandwidth. The Annals of Statistics pp 2655–2675
  • Wahba (1990) Wahba G (1990) Spline models for observational data, vol 59. Siam
  • Yang et al (2015) Yang Y, Tokdar ST, et al (2015) Minimax-optimal nonparametric regression in high dimensions. The Annals of Statistics 43(2):652–674
  • Yeh et al (2008) Yeh I, et al (2008) Modeling slump of concrete with fly ash and superplasticizer. Computers and Concrete 5(6):559–572
  • Yeh (2007)

    Yeh IC (2007) Modeling slump flow of concrete using second-order regressions and artificial neural networks. Cement and Concrete Composites 29(6):474–480

  • Yurugi et al (2000) Yurugi M, Sakata N, Iwai M, Sakai G (2000) Mix proportion for highly workable concrete. Proceedings of Concrete pp 579–589
  • Zhang (2010) Zhang CH (2010) Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics pp 894–942
  • Zou et al (2010) Zou F, Huang H, Lee S, Hoeschele I (2010) Nonparametric bayesian variable selection with applications to multiple quantitative trait loci mapping with epistasis and gene–environment interaction. Genetics 186(1):385–394
  • Zou (2006) Zou H (2006) The adaptive lasso and its oracle properties. JASA 101(476):1418–1429, DOI 10.1198/016214506000000735,
  • Zou and Hastie (2005) Zou H, Hastie T (2005) Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society, Series B 67:301–320