A Partial Least Squares based algorithm for parsimonious variable selection.

BACKGROUND: In genomics, a commonly encountered problem is to extract a subset of variables out of a large set of explanatory variables associated with one or several quantitative or qualitative response variables. An example is to identify associations between codon-usage and phylogeny based definitions of taxonomic groups at different taxonomic levels. Maximum understandability with the smallest number of selected variables, consistency of the selected variables, as well as variation of model performance on test data, are issues to be addressed for such problems.
RESULTS: We present an algorithm balancing the parsimony and the predictive performance of a model. The algorithm is based on variable selection using reduced-rank Partial Least Squares with a regularized elimination. Allowing a marginal decrease in model performance results in a substantial decrease in the number of selected variables. This significantly improves the understandability of the model. Within the approach we have tested and compared three different criteria commonly used in the Partial Least Square modeling paradigm for variable selection; loading weights, regression coefficients and variable importance on projections. The algorithm is applied to a problem of identifying codon variations discriminating different bacterial taxa, which is of particular interest in classifying metagenomics samples. The results are compared with a classical forward selection algorithm, the much used Lasso algorithm as well as Soft-threshold Partial Least Squares variable selection.
CONCLUSIONS: A regularized elimination algorithm based on Partial Least Squares produces results that increase understandability and consistency and reduces the classification error on test data compared to standard approaches.

Background

With the tremendous increase in data collection techniques in modern biology, it has become possible to sample observations on a huge number of genetic, phenotypic and ecological variables simultaneously. It is now much easier to generate immense sets of raw data than to establish relations and provide their biological interpretation [1-3]. Considering cases of supervised statistical learning, huge sets of measured/collected variables are typically used as explanatory variables, all with a potential impact on some response variable, e.g. a phenotype or class label. In many situations we have to deal with data sets having a large number of variables p in comparison to the number of samples n. In such 'large p small n' situations selection of a smaller number of influencing variables is important for increasing the performance of models, to diminish the curse of dimensionality, to speed up the learning process and for interpretation purposes [4,5]. Thus, some kind of variable selection procedure is frequently needed to eliminate unrelated features (noise) for providing a more observant analysis of the relationship between a modest number of explanatory variables and the response. Examples include the selection of gene expression markers for diagnostic purposes, selecting SNP markers for explaining phenotype differences, or as in the example presented here, selecting codon preferences discriminating between different bacterial phyla. The latter is particularly relevant to the classification of samples in metagenomic studies [6]. Multivariate approaches like correspondence analysis and principal component analysis has previously been used to analyze variations in codon usage among genes [7]. However, in order to relate the selection specifically to a response vector, like the phylum assignment, we need a selection based on a supervised learning method.

Partial Least Square (PLS) regression is a supervised method specifically established to address the problem of making good predictions in the 'large p small n' situation, see [8]. PLS in its original form has no implementation of variable selection, since the focus of the method is to find the relevant linear subspace of the explanatory variables, not the variables themselves. However, a very large p and small n can spoil the PLS regression results, as demonstrated by Keles et. al. [9], discovering that the asymptotic consistency of the PLS estimators for univariate responses do not hold, and by [10], who observed a large variation on test set.

Boulesteix has theoretically explored a tight connection between PLS dimension reduction and variable selection [11] and work in this field has existed for many years. Examples are [8,9,11-23]. For an optimum extraction of a set of variables, we need to look for all possible subsets of variables, which is impossible if p is large enough. Normally a set of variables with a reasonable performance is a compromise over the optimal sub set.

In general, variable selection procedures can be categorized [5] into two main groups: filter methods and wrapper methods. Filter methods select variables as a preprocessing step independently of some classifier or prediction model, while wrapper methods are based on some supervised learning approach [12]. Hence, any PLS-based variable selection is a wrapper method. Wrapper methods need some sort of criterion that relies solely on the characteristics of the data as described by [5,12]. One candidate among these criteria is the PLS loading weights, where down-weighting small PLS loading weights is used for variable selection [8,11,13-17,24-27]. A second possibility is to use the magnitude of the PLS regression coefficients for variable selection [18-20,28-34]. Jackknifing and/or bootstrapping on regression coefficients has been utilized to select influencing variables [20,30,31,33,34]. A third commonly used criterion is the Variables Importance on PLS projections (VIP) introduced by Eriksson et. al. [21] and is commonly used in practise [22,31,35-37].

There are several PLS-based wrapper selection algorithms, for example uninformative variable elimination (UVE-PLS) [18], where artificial random variables are added to the data as a reference such that the variable with least performance are eliminated. Iterative PLS (IPLS) adds new variable(s) in the model or remove variables from the model if it improves the model performance [19]. A backward elimination procedure based on leave one variable out is another example [5]. Although wrapper based methods perform well the number of variables selected is still often large [5,12,38], which may make interpretation hard ([23,39,40]).

Among recent advancements in PLS methodology itself we find that Indahl et. al. [41] propose a new data compression method for estimating optimal latent variables classification and regression problems by combining PLS methodology and canonical correlation analysis (CCA), called Canonical Powered PLS (CPPLS). In our work we have adopted this new methodology and proposed a regularized greedy algorithm based on a backward elimination procedure. The focus is on classification problems, but the same concept can be used for prediction problems as well. Our principle idea is to focus on a parsimonious selection, achieved by tolerating a minor performance deviation from any 'optimum' if this gives a substantial decrease in the number of selected variables. This is implemented as a regularization of the search for optimum performance, making the selection less dependent on 'peak performance' and hence more stable. In this respect, the choice of the CPPLS variant is not important, and even the use of non-PLS based methods could in principle be implemented with some minor adjustments. Both loading weights, PLS regression coefficients significance obtained from jackknifing and VIP are explored here for ordering the variables with respect to their importance.

1 Methods

1.1 Model fitting

We consider a classification problem where every object belongs to one out of two possible classes, as indicated by the n × 1 class label vector . From we create the n × 1 numeric response vector by dummy coding, i.e. contains only 0's and 1's. The association between and the n × p predictor matrix is assumed to be explained by the linear model E() = β where β are the p × 1 vector of regression coefficients. The purpose of variable selection is to find a column subset of capable of satisfactory explaining the variations in .

From a modeling perspective, ordinary least square fitting is no option when n

1.2 Canonical Powered PLS (CPPLS) Regression

PLS is an iterative procedure where relation between and is found through the latent variables. The PLS estimate of the regression coefficients for the above given model based on k components can be achieved by

where is the pl × k matrix of -loadings that is summary of X-variables, is the a vector of -loadings i.e. summary of y-variables and is the p × k matrix of loading weights, for details see [8]. Recently, Indahl et. al. [41] proposed a new data compression method for estimating optimal latent variables by combining PLS methodology and canonical correlation analysis (CCA). They introduce a flexible trade-off between the element wise correlations and variances specified by a power parameter γ, ranging from 0 to 1. Defines the loading weights as

where sdenotes the sign of the kcorrelation and Kis a scaling constant assuring unit length w(γ). In this study we restricted γ to lower region (0.001, 0.050) and to upper region (0.950, 0.999). This means we consider combinations of γ for emphasizing either the variance (γ close to 0) or the correlations (γ close to 1). The γ value from above regions that optimizes the canonical correlation is always selected for each component of CPPLS algorithm, see Indahl et. al. [41] for details on the CPPLS algorithm.

Based on the CPPLS estimated regression coefficients we can predict the dummy-variables by

and from the data set we build a classifier using straightforward linear discriminant analysis [42].

1.3 First regularization - model dimension estimate

The CPPLS algorithm assumes that the column space of has a subspace of dimension α containing all information relevant for predicting (known as the relevant subspace) [43]. In order to estimate α we use cross-validation and the performance Pdefined as the fraction of correctly classified observations in a cross-validation procedure, using a components in the CPPLS algorithm.

The cross-validation estimate of α can be found by systematically trying out a range of dimensions a = 1,..., A, compute Pfor each a, and choose as the a where we reach the maximum P. Let us denote this value a*. It is well known that in many cases Pwill be almost equally large for many choices of a. Thus, estimating α by this maximum value is likely to be a rather unstable estimator. To stabilize this estimate we use a regularization based on the principle of parsimony where we search for the smallest possible a whose corresponding performance is not significantly worse than the optimum. If ρis the probability of a correct classification using the a-component model, and ρsimilar for the a*-component model, we test H0 : ρ= ρagainst the alternative H1 : ρ<ρ. In practice Pand Pare estimates of ρand ρ. The smallest a where we cannot reject H0 is our estimate . The testing is done by analyzing the 2 × 2 contingency table of correct and incorrect classifications for the two choices of a, using the McNemar test [44]. This test is appropriate since the model classification at a specific component depends on the model classification at the other components.

This regularization depends on a user-defined rejection level c of the McNemar test. Using a large c (close to 1) means we easily reject H0, and the estimate is often similar to a*. By choosing a smaller c (closer to 0) we get a harder regularization, i.e. a smaller and more stability at the cost of a lower performance.

1.4 Selection criteria

We have implemented and tried out three different criteria for PLS-based variable selection:

1.4.1 Loading weights

Variable j can be eliminated if the relative loading weight, rfor a given PLS component satisfies for some chosen threshold u ∈ [0, 1].

1.4.2 Regression coefficients

Variable j can be eliminated if the corresponding regression coefficient β= 0. Testing H0 : β= 0 against H1: β≠ = 0 can be done by a jackknife t-test. All computations needed have already been done in the cross-validation used for estimating the model dimension α. For each variable we compute the corresponding false discovery rate (q -value) which is based on the p values from jackknifing, and variable j can be eliminated if q>u for some fixed threshold u ∈ [0, 1].

1.4.3 Variable importance on PLS projections (VIP)

VIP for the variable j is defined according to [21] as

where a = 1, 2, ..., a*, wis the loading weight for variable j using a components and , and p2are CPPLS scores, loading weights and y-loadings respectively corresponding to the acomponent. [22] explains the main difference between the regression coefficient βand v. The vweights the contribution of each variable according to the variance explained by each PLS component, i.e. where (/||||)2 represents the importance of the jvariable. Variable j can be eliminated if v 1, see [21,22,36], but a proper threshold between 0.83 and 1.21 can maximize the performance [36].

1.5 Backward elimination

When we have n <