Milos Hauskrecht

Richard Pelikan

Michal Valko

James Lyons- Weiler

FEATURE SELECTION/WEIGHTING AND DIMENSIONALITY REDUCTION IN GENOMICS AND PROTEOMICS

Keywords: Bioinformatics, Classification, Feature selection

Finding reliable, meaningful patterns in data with high numbers of attributes can be extremely difficult. Feature selection helps us to decide what attributes or combination of attributes are most important for finding these patterns. In this chapter, we study feature selection methods for building classification models from high-throughput genomic (microarray) and proteomic (mass spectrometry) data sets. Thousands of feature candidates must be analyzed, compared and combined in such data sets. We describe the basics of four different approaches used for feature selection and illustrate their effects on an MS cancer proteomic data set. The closing discussion provides assistance in performing an analysis in high-dimensional genomic and proteomic data.

Introduction

As technology improves, the amount of information we collect about the world increases. Sensor networks collect traffic or weather information in real-time, documents and news articles are distributed and searched on-line, information in medical records is collected and stored in electronic form. All of this information can be mined so that the relations among components of the underlying systems are better understood and their models can be built. Microarray and mass spectrometry (MS) technologies are producing large quantities of genomic and proteomic data relevant for our understanding of the behavior and function of an organism, or characteristics of disease and its dynamics. Thousands of genes are measured in a typical microarray assay; tens of thousands of measurements comprise a mass spectrometry proteomic profile. The high-dimensional nature of the data demands the development of special data analysis procedures that are able to adequately handle such data. The central question of this process becomes the identification of those features (measurements, attributes) that are most relevant for characterizing the system and its behavior. We study this problem in the context of classification tasks where our goal is to find features that discriminate well among classes of samples, such as samples from people with and without a certain disease.

Feature selection is a process that aims to identify a small subset of features from a large number of features collected in the data set. Two closely-related objectives may drive the feature selection process: (1) Building a reliable classification model which discriminates disease from control samples with high accuracy. The model is then applied to early detection and diagnosis of the disease. (2) Biomarker discovery task where a small set of features (genes in DNA microarrays, or peaks in proteomic spectra) that discriminate well between disease and control groups is identified so that the responsible features can be subjected to further laboratory exploration.

In principle, building a good classification model does not require feature selection. However, when the sample size is small in comparison to the number of features, feature selection may be necessary before a classification model can be reliably learned. With a small sample size, the estimates of parameters of the model may become unreliable and may cause overfitting, a phenomenon in which each datum is fit so rigidly that the model lacks flexibility for future data. To avoid overfitting, feature selection is applied to balance the number of features in proportion to the sample size. On the other hand, identification of a small panel of features for biomarker discovery purposes requires a classification model so that the discriminative behavior of the panel can be assessed.

The dimensionality of typical genomic and proteomic data sets one has to analyze surpasses the number of samples collected in typical studies by a large margin. For example, a typical microarray study can consist of up to a hundred samples with thousands of gene-expression measurements. Mass spectrometry (MS) proteomic profiling is less expensive and as a result one can often see data sets with two to three hundred profiles. MS profiles consist of thousands of measurements. Typically, 'peaks' are selected among those measurements, and number in the hundreds. In either case, feature selection becomes important for both the biomarker discovery and interpretive analysis tasks; one has to seek a robust combination of feature selection methods and classification models to assure their reliability and success. Finally, feature selection may be a one-shot process, but typically, it is a search problem where more than one feature subset is evaluated and compared. Since the number of possible feature subsets is exponential in the number of constituent features, efficient feature selection methods are typically sought.

Feature selection methods are typically divided into three main groups: filter, wrapper and embedded methods. Filter methods rank each feature according to some univariate metric, and only the highest ranking features are used; the remaining features are eliminated. Wrapper algorithms [START_REF] Kohavi | The wrapper approach[END_REF] search for the best subset of features.

To assess the quality of a feature set, these methods rely on and interact with a classification algorithm and its ability to discriminate among the classes. The wrapper algorithm treats a classification algorithm as a black box, so any classification method can be combined with the wrapper. Standard optimization techniques (hill climbing, simulated annealing or genetic algorithms) can be used. Embedded methods search among different feature subsets, but unlike wrappers, the process is tied closely to a certain classification model and takes advantage of its characteristics and structure. In addition to feature selection approaches, in which a subset of original features is searched, the dimensionality problem can be often resolved via feature construction. The process of feature construction builds a new set of features by combining multiple existing features with the expectation that their combination improves our chance to discriminate among the classes as compared to the original feature space.

In this chapter, we first introduce the main ideas of four different methods for feature selection and dimensionality reduction and describe some of their representatives in greater depth. Later, we apply the methods to the analysis of one MS proteomic cancer data set. We analyze each method with respect to the quality of features selected and stress differences among the methods. Since our measuring criterion for feature effectiveness is how well it allows us to classify our samples, we compare the methods and their classification accuracy by combining them with a fixed classification method -a linear support vector machine [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. In closing, we analyze the results and give recommendations on the methods.

Basic Concepts Filter Methods

Filter methods perform feature selection in two steps. In the first step, the filter method assesses each feature individually for its potential in discriminating among classes in the data. In the second step, features falling beyond some thresholding criterion are eliminated, and the smaller set of remaining features is used. This score-and-filter approach has been used in many recent publications, due to its relative simplicity. Scoring methods generally focus on measuring the differences between distributions of features. The resulting score is intended to reflect the quality of each feature in terms of its discriminative power. Many scoring criteria exist. For example, in the Fisher score [START_REF] Pavlidis | Gene functional classification from heterogeneous data[END_REF],

V (i) = µ (+) (i) -µ (-) (i) σ 2 (+) (i) + σ 2 (-) (
i) the quality of each feature is expressed in terms of the difference among the empirical means of two distributions, normalized by the sum of their variances. Table 1.1 displays examples of scoring criteria used in bioinformatics literature. Note that some of the scores can be applied directly to continuous quantities, while others require discretization. Scores can be limited to two classes, like the Fisher score, while others, such as the mutual information score, can be used in the presence of 3 or more classes. For the remainder of this chapter, we will assume our scoring metrics deal with binary decisions, where the data either belong to a positive (+) or negative (-) group.

Criteria based on hypothesis testing.

Some of the scoring criteria are related to statistical hypothesis testing and significance of their results. For example, the t-statistic is related to the null hypothesis H 0 under which the two class-conditional distributions (p(x|y = (+)) and p(x|y = (-))) have the identical mean, that is µ (+) = µ (-) . The degree of violation of H 0 is captured by the p-value of the t-statistic with respect to the Student distribution. As a result, features can be

Criterion References

Fisher Score [START_REF] Golub | Molecular classification of cancer: class discovery and class prediction by gene expression monitoring[END_REF][START_REF] Furey | Support vector machine classification and validation of cancer tissue samples using microarray expression data[END_REF] SAM Scoring Criterion [START_REF] Tusher | Significance analysis of microarrays applied to the ionizing radiation response[END_REF][START_REF] Storey | The analysis of gene expression data: methods and software, chapter SAM thresholding and false discovery rates for detecting differential gene expression in DNA microarrays[END_REF] Student t-test [START_REF] Baldi | A bayesian framework for the analysis of microarray expression data: regularized t-test and statistical inferences of gene changes[END_REF][START_REF] Gosser | The probable error of a mean[END_REF] Mutual Information [START_REF] Tzannes | The mutual information principle and applications[END_REF] χ 2 (Chi Square) [START_REF] Chernoff | The use of maximum likelihood estimates in chi2 tests for goodness-of-fit[END_REF][START_REF] Liu | Chi2: Feature selection and discretization of numeric attributes[END_REF] AUC [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (roc) curve[END_REF] J-measure [START_REF] Smyth | An information theoretic approach to rule induction from databases[END_REF] J5 Score [START_REF] Patel | cageda: a web application for the integrated analysis of global gene expression patterns in cancer[END_REF] See Appendix (table A.1) for mathematical definitions of these scores.

ranked using the inverse of their p-value. Similarly, one can rank the features according to the inverse of the p-value of the Wilcoxon ranksum test [START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF], a nonparametric method, testing the null hypothesis that the class-conditional densities of individual features are equal.

Permutation tests.

Any differential scoring metric (statistic) can be incorporated into and evaluated within the hypothesis testing framework via permutation tests. Permutation (or randomization) tests define a class of non-parametric techniques developed in the statistics literature [START_REF] Kendall | The treatment of ties in ranking problems[END_REF][START_REF] Good | Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses[END_REF], that are used to estimate the probability distribution of a statistic under the null (random) hypothesis from the available data. The estimate of the probability distribution of a scoring metric (Fisher score, J-measure, t-score, etc.) under the null condition allows us to estimate the p-value of the score observed in the data, similarly to the t-test or Wilcoxon rank-sum test. From the viewpoint of feature selection, the null hypothesis assumes that the conditional probability distributions for the two classes (y = (+) or (-)) are identical under a feature x, that is, p(x|y = (+)) = p(x|y = (-)); or equivalently, that the data and the labels are independent, p(x, y) = p(x)p(y). The distribution of data under the null hypothesis is generated where N T b ≥T is the number of permutations for which the test statistic T b is better than T ; Return p; } through random permutations (of labels) in the data. The permutation test algorithm is shown in figure 1.1. The main cycle of the algorithm either scans through all possible permutations of labels, or, if this set is too large, a large number B of permutations is generated randomly. With sufficient cycles, the distribution of the test statistic under the null hypothesis can be estimated reliably.

Choosing features based on the score.

Differential scores or their associated p-value scores allow us to rank all feature candidates. However, it is still not clear how many features should be filtered out. The task is easy if we always seek a fixed set of k features. In such a case, the top k features are selected with respect to the ordering imposed by ranking features by their score. However, the quality of these features may vary widely, so selecting the features based solely on the order may cause some poor features to be included in the set. An alternative method is to choose features by introducing a threshold on the value of the score. Unfortunately, not every scoring criterion has an interpretable meaning, so it is unclear how to select an appropriate threshold. The statistic typically used for this purpose is the p-value associated with the hypothesis test. For example, if the p-value threshold is 0.05 then there is a 5% chance the feature is not differentially expressed at the threshold value. Such a setting allows us to control the chance of false positive selections. These are features which appear discriminative by chance.

Feature set selection and controlling false positives.

The high-dimensional nature of biological data sources necessitates that many features (genes or MS-profile peaks) be tested and evaluated simultaneously. Unfortunately, this increases the chance that false positives are selected. To illustrate this, assume we measure the expression of 10,000 independent genes and none of them are differentially expressed. Despite the fact that there is no differential expression, we might expect 100 features to have their p-value smaller than 0.01. An individual feature with p-value 0.01 may appear good in isolation, but may become a suspect if it is selected from thousands of tested features. In such a case, the p-value of the combined set of the top 100 features selected out of 10,000 is quite different. Thus, adjustment of the p-value when performing multiple tests in parallel is necessary.

The Bonferroni correction adjusts the p-value for each individual test by dividing the target p-value for all findings by the number of findings. This assures that the probability of falsely rejecting any null hypotheses is less than or equal to the target p. The limitation of the Bonferroni correction is that it operates under the assumption of independence and as a result it is too conservative if features are correlated. Two alternatives to the Bonferroni correction are offered by: (1) the Familywise Error Rate method (FWER, [START_REF] Westfall | Resamplingbased multiple testing: examples and methods for p-value adjustment[END_REF]) and ( 2) methods for controlling the False Discovery Rate (FDR, [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF][START_REF] Tusher | Significance analysis of microarrays applied to the ionizing radiation response[END_REF]. FWER takes into account the dependence structure among features, which often translates to higher power. [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] suggest to control FDR instead of the p-value. The FDR is defined as the mean of the number of false rejections divided by the total number of rejections. The Significance Analysis of Microarrays (SAM) method [START_REF] Storey | The analysis of gene expression data: methods and software, chapter SAM thresholding and false discovery rates for detecting differential gene expression in DNA microarrays[END_REF] is used as an estimate of the FDR. Depending on the chosen threshold value for the test statistic T , it estimates the expected proportion of false positives on the feature list using a permutation scheme.

Correlation filtering.

To keep the feature set small, the objective is to diversify the features as much as possible. The selected features should be discriminative as well as independent from each other as much as possible. The rationale is that two or more independent features will be able to discriminate the two classes better than any of them individually. Each feature may differentiate different sets of data well, and independence between the features tends to reduce the overlap of the sets. Similarly, highly dependent features tend to favor the same data and thus are less likely to help when both are included in the panel. The extreme case is when the two features are exact duplicates, in which case one feature can be eliminated.

Correlation filters [START_REF] Ross | Systematic variation in gene expression patterns in human cancer cell lines[END_REF][START_REF] Hauskrecht | Feature selection for classification of seldi-tof-ms proteomic profiles[END_REF] try to remove highly correlated features since these are less likely to add new discriminative information [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. Various elimination schemes are used within these filters to reduce the chance of selected features being highly correlated. Typically, correlation filters are used in combination with other differential scoring methods. For example, features can be selected incrementally according to their p-value; the feature to be added next is checked for correlation with previously selected features. If the new feature exceeds some correlation threshold, it is eliminated [START_REF] Hauskrecht | Feature selection for classification of seldi-tof-ms proteomic profiles[END_REF].

Wrapper Methods

Wrapper methods [START_REF] Kohavi | The wrapper approach[END_REF] search for the best feature subset in combination with a fixed classification method. The goodness of a feature subset is determined using internal-validation methods, such as, k-fold or leave-one-out cross-validation [START_REF] Krus | Computer-assisted multicrossvalidation in regression analysis[END_REF]. Since the number of all combinations is exponential in the number of features, the efficiency of the search methods is often critical for its practical acceptance. Different heuristic optimization frameworks have been applied to search for the best subset. These include: forward selection, backward elimination [START_REF] Blum | Selection of relevant features and examples in machine learning[END_REF], hill-climbing, beam search (Russel and Norvig, 1995), and randomized algorithms such as genetic algorithms [START_REF] Koza | Survey of genetic algorithms and genetic programming[END_REF] or simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. In general, these methods explore the search space (subsets of all features) starting with no features, all features, or a random selection of features. For example, the forward selection approach builds a feature set by starting from an empty feature set and incrementally adding the feature that improves the current feature set the most. The procedure stops when no improvement in the feature set quality is possible.

Embedded Methods

Embedded methods incorporate variable selection as part of the model building process. A classic example of an embedded method is CART (Classification and Regression Trees, [START_REF] Breiman | Classification and regression trees[END_REF]).

CART searches the range of each individual feature to find the split that optimally divides the observed data into a more homogeneous groups (with respect to the outcome variable). Beginning with the subsets of the variable that produces the most homogeneous split, each variable is again searched across its range to find the next optimal split. This pro-cess is continued within each new subset until all data are perfectly fit by the resulting tree, or the terminal nodes have a small sample size. The group constituting the majority of data points in each node determines the classification accuracy of the derived terminal nodes. Misclassification error from internal cross-validation can be used to backprune the decision tree and optimize its projected generalization performance on additional independent test examples.

Regularization/Shrinkage methods.

Regularization or shrinkage methods [START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Xing | Feature selection for high-dimensional genomic microarray data[END_REF] offer an alternative way to learn classifications for data sets with large number of features but small sample size. These methods trim the space of features directly during classification. In other words, regularization effectively shuts down (or zeros the influence of) unnecessary features.

Regularization can be incorporated either into the error criterion or directly into the model. Let w be a set of parameters defining a classification model (e.g., the weights of a logistic regression model), and let Error(w, D) be an error function reflecting the fit of the model to data (e.g., least-squares as likelihood-based error). A regularized error function is then defined as:

Error Reg (w, D) = Error(w, D) + λ||w||,
where λ > 0 is a regularization constant, and ||.|| is either the L1 or L2 norm. Intuitively, the regularization term penalizes the model for nonzero weights so the optimization of the new error function drives all unnecessary parameters to 0. Automatic Relevance Determination (ARD, [START_REF] Neal | Assessing relevance determination methods using delve generalization in neural networks and machine learning[END_REF][START_REF] Mackay | Bayesian Methods for Adaptive Models[END_REF]) achieves regularization effects in a slightly different way. The relevance of an individual feature is represented explicitly via model parameters and the values of these parameters are learned through Bayesian methods. In both cases, the output of the learning is a feature-restricted classification model, so features are selected in parallel with model learning.

Support vector machines.

Regularization effects are at work also in one of the most popular classification frameworks these days: the support vector machine (SVM) [START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF][START_REF] Schölkopf | Learning with Kernels[END_REF]. The SVM defines a linear decision boundary (hyperplane) that separates case and control examples. The boundary maximizes the distance (also called margin) in between the two sample groups. The effects of margin optimization are twofold: only a small set of data points (support vectors) are critical for the separation; the dimensions unncessary for separation are penalized. Both of these processes help to fight the prob-lem of model overfit. As a result, the SVM offers a robust classification framework that works very well for situations with a moderately large number of features and relatively small sample sizes.

Feature construction

Better discriminatory performance can be often achieved using features constructed from the original input features. Building a new feature is an opportunity to incorporate domain specific knowledge into the process and hence to improve the quality of features. Nevertheless, a number of generic feature construction methods exist: clustering; linear (affine) projections of the original feature space; as well as more sophisticated space transformations such as wavelet or kernel transforms. In the following, we briefly review three basic feature construction approaches: clustering, PCA and linear discriminative projections.

Clustering.

Clustering groups data components (data points or features) according to their similarity. Every data component is assigned to one of the groups (clusters); components falling into the same cluster are assigned the same value in the new (reduced) representation. Clustering is typically used to identify distinguished sample groups in data [START_REF] Ben-Dor | Tissue classification with gene expression profiles[END_REF][START_REF] Slonim | Class prediction and discovery using gene expression data[END_REF]. In contrast to supervised learning techniques that rely heavily on class label information, clustering is unsupervised and the information about the target groups (classes) is not used. From the dimensionality reduction perspective, a data point is assigned a cluster label which is then used as its representation.

Clustering methods rely on the similarity matrix -a matrix of distances between data components. The similarity matrix can be built using one of the standard distance metrics such as Euclidean, Mahalanobis, Minkowski, etc, but more complex distances based on, for example, functional similarity of genes [START_REF] Speer | Spectral clustering gene ontology terms to group genes by function[END_REF], are possible. Table 1.2 gives a list of some standard distance metrics one may use in clustering.

Clustering algorithms.

The goal of clustering is to optimize intra and inter cluster distances among the components. Two basic clustering algorithms are: k-means clustering [START_REF] Mcqueen | Some methods for classification and analysis of multivariate observations[END_REF][START_REF] Ball | A clustering technique for summarizing multivariate data[END_REF], and hierarchical agglomerative clustering [START_REF] Cormack | A review of classification[END_REF][START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF].

Briefly, the k-means algorithm clusters data into groups by iteratively optimizing positions of cluster centers (means) so that the sum of within-cluster distances (the distances between data points and their cluster centers) is minimized. Initial positions for cluster centers are gen-

(r, s) = (xr -xs)trace(Σ) -1 (xr -xs) ′ Mahalanobis distance d(r, s) = (xr -xs)Σ -1 (xr -xs) ′ City Block metric d(r, s) = n j=1 |xrj -xsj| Minkowski metric d(r, s) = p n j=1 |xrj -xsj| p Cosine distance d(r, s) = 1 - xr x ′ s √ x ′ r xr √ x ′ s x s) Correlation distance d(r, s) = 1 - (xr -xr )(xs-xs) ′ √ (xr -xr )(xr -xr ) ′ √ (xs-xs)(xs-xs) ′ Hamming distance d(r, s) = #(x rj =x sj ) n Jaccard distance d(r, s) = #[(x rj =x sj )∧((x rj =0)∨(x sj =0))] #[(x rj =0)∨(x sj =0)]
x and x ′ denote a column vector and its transpose respectively. xr and xs indicate the r th and s th samples in the data set, respectively.

x rj indicates the j th feature of the r th sample in the data set.

xr indicates the mean of all features in the r th sample in the data set. Σ is the sample covariance matrix.

The symbol # denotes counts; the number of instances satisfying the associated property.

erated randomly or by using heuristics. The algorithm is not guaranteed to converge to the optimal solution. On the other hand, hierarchical agglomerative methods work by combining pairs of data entities (features) or clusters into a hierarchical structure (called a dendrogram). The algorithm starts from unit clusters and merges them greedily (i.e. choosing the merge which most improves the fit of the clusters to the data) into larger clusters using an a priori selected similarity measure.

Probabilistic (soft) clustering.

The k-means and agglomerative clustering methods assign every data point into a single cluster. However, sometimes it may be hard to decide what cluster the point belongs to. In probabilistic (soft) clustering methods, a data point belongs to all clusters, but the strength (weight) of its association with clusters differs by how well it fits cluster descriptions. Typically, the weight has probabilistic meaning and defines a probability with which a data point belongs to a cluster.

To calculate the probability, an underlying probabilistic model must be first fit to the data. Briefly, data are assumed to be generated from k different classes that correspond to clusters. Each class has its own distribution for generating data points. The parameters of these distributions as well as class (cluster) priors are fit (learned) using Expectation-Maximization techniques [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. Once the model parameters are known, the probabilistic weights relating a data point and clusters are posterior probabilities of the point belonging to classes. A classic example of a probabilistc model often used in clustering is the Mixture of Gaussians model [START_REF] Mclachlan | Clustering via normal mixture models[END_REF], where k clusters are modeled using k Gaussian distributions.

Clustering features.

Clustering methods can be applied to group either data points or features in the data. When clustering features, the dimensionality reduction is achieved by selecting a representative feature (typically the feature that is closest to the cluster center, [START_REF] Guyon | An introduction to variable and feature selection[END_REF]), or by aggregating all features within the cluster via averaging to build a new (mean) feature. If we assume k different feature clusters, the original feature space is reduced to a new k-dimensional space. An example method of feature clustering is to cluster features based on intra-correlation, and use the cluster center as a representative. Closely correlated features are not likely to help when separated, so grouping them away from more unrelated features will help diversify the resulting features.

Principal component analysis.

Principal Component Analysis (PCA), [START_REF] Jolliffe | Principal Component Analysis[END_REF]) is a widely used method for reducing the dimensionality of data. PCA finds projections of high dimensional data into a lower dimensional subspace such that the variance retained in the projected data is maximized. Equivalently, PCA gives uncorrelated linear projections of data while minimizing their least square reconstruction error. Additionally, PCA works fully unsupervised; class labels are ignored. PCA can be extended to nonlinear projections using kernel methods [START_REF] Bach | Kernel independent component analysis[END_REF]. Dimensionality reduction methods similar to PCA that let us project high dimensional features into a lower dimensional space include multidimensional scaling (MDS) [START_REF] Cox | Multidimensional Scaling[END_REF] used often for data visualization purposes or independent component analysis (ICA) [START_REF] Jutten | Blind separation of sources, part 1: An adaptive algorithm based on neuromimetic architecture[END_REF].

Discriminative projections.

Principal component analysis identifies affine (linear) projections of data that maximize the variance observed in data. The method operates in a fully unsupervised manner; no knowledge of class labels is used to find the principal projections. The question is whether there is a way to identify linear projections of features such that they optimize the discriminability among the two classes. Techniques which try to achieve this goal include Fisher's linear discriminant (FLD) [START_REF] Duda | Pattern Classification[END_REF], linear discriminant analysis [START_REF] Hastie | The Elements of Statistical Learning[END_REF] and more complex methods like partial least squares (PLS) [START_REF] Denham | Implementing partial least squares[END_REF][START_REF] Dijkstra | Some comments on maximum likelihood and partial least squares methods[END_REF].

Take for example, the linear discriminant analysis model. The model assumes that cases and controls are generated from two Gaussian distributions with means µ (-) , µ (+) and the same covariance matrix Σ. The parameters of the two distributions are estimated from data using the maximum likelihood methods. The decision boundary that is defined by data points that give the same probability for both distributions is a line. The linear projection is defined as:

w = Σ -1 (µ (+) -µ (-) ),
where µ (-) , µ (+) are the means of the two groups and Σ is the covariance for both groups, where p(x|y) ∼ N (µ, Σ).

Advantages and Disadvantages

Each of the aforementioned methods comes with advantages and disadvantages. The following text briefly summarizes them.

Filter methods:

Advantages: Univariate scores are very easy to calculate and thus, filter methods have a short running time. If our goal is a prediction, they often perform well in combinations with more robust classification methods such as the SVM.

Disadvantages: Many differential scoring methods exist, it is unclear which one is best for the data set at hand. The features are analyzed independent of each other. This is a problem if our goal is to identify a small panel of discriminative features (biomarkers). Multivariate relations/dependencies must be incorporated through additional criteria, e.g. correlation filters.

Wrapper methods

Advantages: More comprehensive search of the feature-set space.

The feature set with the best discriminative potential on a fixed classification method is selected.

Disadvantages: Running time is much longer than filter methods; many feature sets need to be analyzed and assessed. In addition, scoring of feature sets is based on internal cross-validation methods which lengthens their running time. The reliability of the estimate of the internal cross-validation error needs to be considered. Low reliability of the internal validation error in combination with a large number of subsets examined can be lethal especially in various greedy search schemes.

Embedded methods:

Advantages: Features and their selection are tuned to a specific model. Learning methods which incorporate aspects of regularization, like the SVM or regularized logistic regression, can learn very good predictive models even in the presence of high-dimensional data. We recommend trying SVM as a first step if the goal is only to build a predictive model.

Disadvantages: Identification of a small set of features may be problematic. Backward feature elimination routines [START_REF] Guyon | An introduction to variable and feature selection[END_REF] can be used to reduce the feature panel to a more reasonable size.

Feature construction methods:

Advantages: May incorporate the domain knowledge which may translate to improved feature sets.

Disadvantages: If features are constructed using one of the out-ofbox methods (e.g., PCA) the new features may be hard to interpret biologically. In addition, many feature construction techniques (e.g. clustering, PCA, ICA) work in an unsupervised mode, so high quality features for discriminatory purposes are not guaranteed.

Case Study: Pancreatic Cancer

To illustrate some of the advantages and disadvantages of feature selection methods, we use a data set of MS proteomic profiles for pancreatic cancer collected at the University of Pittsburgh Cancer Institute (UPCI). Since full feature-selection comparison is very hard to do without a full predictive model that combines both the feature-selection and the classification stages we test feature selection methods in combination with one classification method -the linear support vector machine (SVM) [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. All classification results presented in the following text were obtained by using the repeated random subsampling strategy with 40 different train/test data splits using 70/30 train/test split ratio. The optimization criterion for the SVM method was a zeroone loss function, which focuses on improving classification error instead of sensitivity or specificity. The statistics reported are: ACE (average test classification error), SN (Sentivity) and SP (Specificity) and their standard deviations.

Data and preprocessing

The data set consists of 116 MS profiles, with 57 cancer cases (+ group) and 59 controls, matched according to their smoking history, age, and gender (-group). The data were generated using Ciphergen Biosystems Inc. SELDI-TOF (Surface-Enhanced Laser Desorption/Ionization Time-of-Flight) mass spectrometry. Compounds such as proteins, peptides and nucleic acids for masses of up to 200 thousand Daltons are recorded using this technology. Before applying feature selection techniques the data set was preprocessed using the Proteomic Data Analysis Package (PDAP, [START_REF] Hauskrecht | Feature selection for classification of seldi-tof-ms proteomic profiles[END_REF]). The following preprocessing steps were applied: (1) cuberoot variance stabilization, (2) local minwindow baseline correction, (3) Gaussian kernel smoothing, (4) rangerestricted intensity normalization, and (5) peak-based profile alignment. The quality of all profiles were tested beforehand on raw MS profile readings using total ion current (TIC). None of the profiles differed by more than two standard deviations from the mean TIC, which is our current quality-assurance/quality-control threshold for sample exclusion. After basic preprocessing, peaks in the range of 1500-1650 Daltons were identified and their corresponding intensities were extracted.1 This gave us a data set of 116 samples with 602 peak features.

Filter methods

Basic Filter Methods.

Many univariate scoring metrics that assess the individual quality of features were proposed in the literature. An important question is how the rankings and subsequent feature selection induced by these metrics vary. Table 1.3 shows the number of overlapping features for the top 20 features selected according to four frequently used scoring criteria: correlation, Fisher, t-statistic and Wilcoxon's pvalue measures. The table shows that different scoring metrics may induce rather different feature orders and as a result, different feature panels. It is very hard to argue that any one of them is the best. The quality depends strongly on the classification technique used in the next step, but even there the story is often unclear, and the best method tends to vary among the data sets. Table 1.4 illustrates the results obtained using top 20 choices of four scoring methods from Table 1.3 after we combine them with the linear SVM model. Standard deviations of performance statistics are also given. We see that the best classification error was obtained using the features selected based on the t-statistic score. While our experience is that the t-statistic score performs well on many proteomic data sets, other scoring metrics may often outperform it.

Controlling false positive selections.

A problem with high dimensional data is that some features may appear as good discriminators simply by chance. The problem of false positive identifications of features is critical for the biomarker discovery task. Clearly, a more comprehensive analysis and validation of the feature in the lab may incur a significant monetary cost. While positive feature selections may influence also the generalizations of the predictive model and its classification accuracy, the classification methods are often more robust to handle them and the problem of false positive features is less pressing than for the biomarker discovery applications.

The false positive selection rate can be controlled via p-value on individual features, Bonferroni corrected p-value for the panel of features, or through false discovery rate. Table 1.5 shows the number of features out of 602 original features selected by each of these methods. Assuming that all features are independent and random, we expect to see about 30 false positive features under the simple p-value of 0.05 for each feature. Using this estimate and the fact that we see only 13 features for the p-value of 0.05 would lead us to the conclusion that all of these are likely obtained by chance. The caveat is that when features are dependent and correlated the expected numbers are very different. Indeed, features in this and other proteomic data sets exhibit a large amount of correlation among the features; so the result in the table is indicative of such a dependency. The Bonferroni correction typically leads to a very conservative bound that may be very hard to satisfy. For example, none of the features in our cancer data passed Bonferronicorrected p-value of 0.05. FWER and FDR methods and their thresholds give better estimates of false positive selections and their rates for the real-world data and should be preferred over simple and Bonferronicorrected p-value thresholding.

When selecting features, our objective is to strike the right balance between the number of features, the flexibility they may offer when building multivariate discriminators, and the risk of inclusion of false positive features. The FWER and FDR methods give better control over risks of false positives. However, choosing the optimal thresholds for these techniques is a matter of personal preference. For example, two different approaches can be taken. If the selected features are meant only for use with an automated classification routine, it may be more acceptable to risk selecting false positives, and thusly the threshold can be less stringent. On the other hand, if the selected features are to be investigated more thoroughly (e.g. to analyze them using wetlab techniques), it would be far less acceptable to suggest that false positives are informative features. In this case, the threshold should be set more aggressively.

Correlation filters.

Biological (genomic and proteomic) data sets often exhibit a relatively high number of correlations. The correlations can be introduced by the technology producing the data or they reflect true underlying dependencies among measured species. For example, a peak in a proteomic profile is formed by a collection of correlated measurements, triple or double charged ions cause the same signal to be replicated at different parts of the profiles, and finally some peaks are correlated because they share a common regulatory (or interaction) pathway.

Selecting two features that are near duplicates, even if they are highly discriminative, does not help the classification model and its accuracy. Correlation filtering alleviates the problem by removing features highly correlated with existing features in the panel. Table 1.6 illustrates the number of features one obtains by filtering out correlated features at different maximum allowed absolute correlation (MAC) thresholds from the original 602 features. We note that the amounts of correlates filtered out at higher thresholds are statistically significantly different (at p=0.01) from what one would obtain for independent feature sets. Figure 1.2 illustrates the effect of correlation filtering when it is combined with the univariate feature scoring based on the t-statistic. We see that test errors for smaller feature sets (size 5) are improved if feature panels are decorrelated. However, for larger feature panels the effect of feature decorrelation may vanish since some good features that add some discriminative value to the panel are filtered out. For example, for 20 features in Figure 1.2 the effect of correlation filtering has disappeared and the SVM classifier based on the unrestricted t-statistic score performs better than classifiers with correlation thresholds of 0.75 and 0.5. This illustrates one of the problems of the method, identification of an appropriate MAC threshold. We must note that the effect as seen in Figure 1.2 may be less pronounced on other classification methods or on other data sets, while in some cases correlation thresholds may lead to superior performance. These differing outcomes are the results of tradeoffs of feature quality and overfit processes.

The plain correlation threshold filtering method suffers from a couple of problems. First, an identification of an appropriate correlation threshold in advance is hard. Moreover, for different feature sizes there appears to be a different threshold that works best so switching of thresholds may be appropriate. One solution to this problem is the parallel correlation filtering method [START_REF] Hauskrecht | Feature selection for classification of seldi-tof-ms proteomic profiles[END_REF]) that works at multiple correlation threshold levels in parallel and uses internal cross-validation methods to decide on what feature (correlation level) to select next. The performance of the method is compared to the unresticted t-statistic filter and two correlation filtering methods based on simple MAC thresholds in Figure 1.2. Wrapper methods search for the best subset of features by trying them in combination with a fixed classification method. However, there is a natural tradeoff between the quality of the feature set found, and the time taken to search for it. Table 1.7 displays performance statistics for two search methods: greedy forward selection and simulated annealing. The forward selection approach, also called the greedy approach, adds the feature which improves the set the most. The panel begins empty and is built incrementally, stopping when no improvement in the feature set is possible. Simulated annealing is a randomized algorithm and if it is left to search long enough all possible combinations may be reached and evaluated. Thus, simulated annealing may arrive at a better solution than the greedy method when given enough time. This quality/time tradeoff is captured in the table. The model based on the greedy forward selection method leads to average errors of 0.1750 while simulated annealing approaches 0.1660. To reach the result, 7037.4 feature sets were evaluated on average by forward selection, while simulated annealing was run for 10000 steps on every train/test split.

Wrapper methods

Evaluating a new feature set in any wrapper method is done by internal validation methods, such as k-fold cross-validation or leave-one-out validation. The overhead incurred by the evaluation step contributes to the running time of the algorithm. In general, using more internal splits improves the estimate of the error for each feature set. The price paid for it is an additional increase in the running time. Despite the downfalls, the results obtained from wrapper methods powered by various search heuristics are often quite good, especially when computational time is not an issue.

Embedded methods

Table 1.8 shows the results of three classification methods with embedded feature selection: CART [START_REF] Breiman | Classification and regression trees[END_REF], Regularized logistic regression (RLR) [START_REF] Hastie | The Elements of Statistical Learning[END_REF] and Support Vector Machines (SVM) [START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF]. Each of these methods handles features differently, and consequently leads to different classification accuracies. We see that two of the methods, RLR and SVM, achieved results comparable or better than filter and wrapper methods. While this is not the rule, the linear SVM appears to be a very stable method across a large range of features so we always recommend to try it on the full feature set. Embedded methods may not be optimal, if we want to use them for biomarker discovery, that is, if our objective is to find a small set of original features with a good discriminatory performance. The embedded methods may rely on too many features so a follow-up selection of a smaller subset is necessary. Wrapper methods based on the backward feature elimination [START_REF] Guyon | An introduction to variable and feature selection[END_REF] achieve this by gradually eliminating the features that affect the performance the least.

Feature construction methods

To illustrate feature construction methods we use three unsupervised methods: sample clustering, feature clustering and PCA projections, all aimed to reduce the dimensionality of data. The results of these methods in combination with the linear SVM are in table 1.9. The first entry in the table (sample clustering with Euclidean distance) illustrates the major weakness of clustering methods: the clustering does not give reasoning as to why the data components group together, other that their distance is close, which obviously depends on the choice of the metric. Thus, one has to assure that the distance selected is not arbitrary and makes sense for the data and the prediction task. The result for clustering of features based on the correlation metric also supports this point. There are many feature correlates in the proteomic data set, so grouping the features based on their mutual correlation and replacing the features in each cluster with a feature corresponding to the cluster center tends to eliminate high correlates in the new (reduced) data. This is very similar in spirit to the correlation filtering method. The difference is that the correlation filtering is closely combined with and benefits from the univariate score filtering, while correlation clustering works fully unsupervised.

PCA constructs features using linear projections of complete data. Since PCA arranges projections along uncorrelated axes, it helps to relieve us from identifying feature correlates. As a result, we see an improvement in classification error over some other construction and filtering methods. Note that PCA can be a good "one shot" technique, avoiding necessities like the choice of the number of clusters, k, in k-means clustering, or scoring metric in filtering methods. The effort saved by not choosing parameters is in exchange for knowledge about a targetable panel of biomarkers, but PCA can still be convenient if the only interest is constructing a predictive model.

Summary of analysis results and recommendations

There are multiple feature selection / dimensionality reduction methods one may apply to reduce the feature size of the data and make it 'comparable' to its sample size. Unfortunately, there is no perfect recipe for what method to choose but here are some guidelines.

Having prior information about how features can be related to the prediction task will always help feature selection and its subsequent application. So whenever possible try to use this information. For example, when the biological relevance of features can be ascertained, the potentially irrelevant or obvious features can also be eliminated.

In the presence of no prior information, more generic information can be used for steering feature selection in the right direction. The effect of a feature on the target class and the presence of multivariate dependencies (e.g. correlations) among feature candidates appear to be the most important ones. The importance of a feature is captured by a univariate scoring metric. Dealing with highly correlated features, either by grouping them or eliminating redundancies, can help the selection process by narrowing the choice of features.

Feature selection coupled with more robust classification methods, like SVM, can perform extremely well on all features. Backward feature elimination methods can be applied if we would like to identify a smaller panel of informative features.

The feature selection method applied to data does not have to match a single method. A combination of feature selection methods may be beneficial and may work much better [START_REF] Xing | Feature selection for high-dimensional genomic microarray data[END_REF]. For example, it may help to exclude some features outright with a basic filtering method by removing the lowest-scoring features and apply other methods (e.g. wrapper or PCA methods) only on the remaining features.

Since there are many feature selection methods, one may be tempted to try many of them in combination with a specific classifier and pick the one that gives the best test-set result post-hoc. Note that in such a case the error is biased and does not objectively report on the generalizability of the approach. Model selection methods based, for example, on an internal cross-validation loop should be applied whenever a choice out of many candidates is allowed.

In closing, it is important to note that the selection of the feature selection technique should first be driven by prior knowledge about the data, and then by the primary goal you wish to accomplish by analyzing the data: obtain a small, easy to interpret, feature panel or build a good classification model. Feature selection techniques vary in their complexity and interpretability, and the issues discussed above must be taken into careful consideration.

Conclusions

In this chapter, we have presented four basic approaches to feature selection and dimensionality reduction. Filter, wrapper, and embedded methods work with the available features and choose those which appear important. In slight contrast, feature construction methods build new features which can be more powerful than previous ones. To discuss the entire gamut of feature selection methods would be exhaustive, as researchers must constantly meet their needs of analyzing highdimensional data. The techniques covered here are among the most effective for analyzing genomic and proteomic data, in terms of building predictive models and developing biologically relevant information.

Figure 1 . 1 .

 11 Figure 1.1. Permutation Test Algorithm

Figure 1 . 2 .

 12 Figure 1.2. Effect of correlation filtering on classification errors. Results of correlation filtering on the t-statistic score and SVM are shown.

Table 1 .

 1 1. Examples of Univariate Scoring Criteria for Filter Methods. See Appendix for mathematical definitions of these scores.

Table 1 .

 1 2. Examples of Distance Metrics for Clustering

	Metric

Table 1 .

 1 3. Overlap of top 20 features for four different metrics.

		Correlation	Fisher	t-statistics	Wilcoxon
	Correlation	-	18	12	18
	Fisher	18	-	11	16
	t-statistics	12	11	-	11
	Wilcoxon	18	16	11	-

Table 1 .

 1 4. Results for classifiers based on different feature filtering methods and the linear SVM. Standard deviations are given in parentheses.

		Correlation	Fisher	t-statistics	Wilcoxon
	ACE	0.2500 (0.1178)	0.2188 (0.1075)	0.1743 (0.0684)	0.2611 (0.1091)
	SN	0.8022 (0.0387)	0.8102 (0.1210)	0.8259 (0.0997)	0.7956 (0.1200)
	SP	0.7142 (0.1249)	0.7628 (0.1423)	0.8327 (0.0852)	0.6961 (0.1607)
	PPV	0.7229 (0.1452)	0.7580 (0.1265)	0.8179 (0.0881)	0.7115 (0.1358)

Table 1

 1

	.5. p-value for t-statistics
	original number of features	p < 0.05	Bonferroni p < 0.05	FDR 0.2
	602	13	0	5

Table 1

 1

	.6. Effect of correlation filtering
	threshold	1	0.9	0.8	0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
	number of features 602 460 247 119 52	22	12	9	6	3	1

Table 1 .

 1 7. Wrapper methods with two search algorithms: forward selection and simulated annealing. Standard deviations are given in parentheses.

	602	Greedy	Simulated Annealing
	ACE	0.1750 (0.0668)	0.1660 (0.0603)
	SN	0.8239 (0.1123)	0.8149 (0.1097)
	SP	0.8261 (0.1100)	0.8614 (0.0784)
	# steps	7037.4	10000

Table 1 .

 1 8. Performance statistics for embedded methods. Standard deviations are given in parentheses.

	602	CART	Regularized LR	SVM
	ACE	0.3681 (0.0897)	0.1382 (0.0584)	0.1382 (0.0623)
	SN	0.6321 (0.1888)	0.8619 (0.1026)	0.8536 (0.0913)
	SP	0.6361 (0.2088)	0.8624 (0.0942)	0.8769 (0.0881)
	PPV	0.6270 (0.1442)	0.8561 (0.0933)	0.8617 (0.0952)

Table 1 .

 1 9. Construction methods: Sample clustering using squared Euclidian distance, feature clustering using correlation coefficient, and PCA. Standard deviations are given in parentheses.

	602	sample clustering	feature clustering	PCA projections
	ACE	0.4525 (0.0810)	0.2104 (0.0652)	0.1681 (0.0594)
	SN	0.4721 (0.1604)	0.7932 (0.1426)	0.8223 (0.0984)
	SP	0.6444 (0.1633)	0.7968 (0.0920)	0.8492 (0.0842)

The region below 1500 Daltons is unsuitable for analysis because of known signal reproducibility problems. The region is often referred to as the junk region. On the other hand, signals for higher mass-to-charge-ratios are of lower intensity which makes them hard to separate from the noise. An a priori upper limit is typically set to restrict the search for signal.

Appendix Filtering Methods

Table A.1. Formulae for popular filter scores

Filter Name Formula

Fisher Score score(i)

SAM scoring criterion.

The standard SAM technique is meant to be used in a permutation setting, however, the scoring criteria can still be used for filtering methods.

The correcting constants s(i) and s0 are computed as follows:

s0 = 1 for purposes of simplicity.


Michal Valko
email: michal@cs.pitt.edu

Nuno C Marques

Marco Castellani

Evolutionary Feature Selection for Spiking Neural Network Pattern Classifiers

Keywords: Artificial Intelligence, Neural Networks, Genetic Algorithms, Spiking Neuron Models, JASTAP

This paper presents an application of the biologically realistic JASTAP neural network model to classification tasks. The JASTAP neural network model is presented as an alternative to the basic multi-layer perceptron model. An evolutionary procedure previously applied to the simultaneous solution of feature selection and neural network training on standard multi-layer perceptrons is extended with JASTAP model. Preliminary results on IRIS standard data set give evidence that this extension allows the use of smaller neural networks that can handle noisier data without any degradation in classification accuracy.

I. INTRODUCTION

Classification of patterns is important to many data mining processes. Among possible classifiers, artificial neural network classifiers have proven to be one of the most robust classification systems. Their capability to deal with noisy input patterns, to handle both noisy and continuous value data did prove to be an essential tool for classification and prediction [START_REF] Mitchell | Machine Learning[END_REF].

This paper aims to present an alternative to the basic perceptron model [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF] by using the neural realistic JASTAP model [START_REF] Janco | Modeling of neuronal functions: A neuronlike element with the graded response[END_REF]. JASTAP model implements biologically plausible neural networks in a well parameterized model that can be statistically related [START_REF] Valko | Evolving neural networks for statistical decision theory[END_REF]. Since the usual simplifications for other spiking models [START_REF]Pulsed Neural Networks[END_REF] are not done a priori, we can keep all the flexibility known in biological neurons. On the same time, since parameters can also be predefined in advance (giving rise to more simple models) JASTAP proves as a very powerful new way of representing both model knowledge and experience. This way, only needed features can be activated in the model. A first demonstrative study on the use of this new kind of artificial neural network for classification tasks is presented on IRIS dataset.

Standard backpropagation learning [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF] on JASTAP neural networks (and on simpler spiking models in general) is still an open issue (mainly due to the different weights that need to be adjusted and to cyclic connections). So, a more generic and powerful learning method is used: we have incorporated in JASTAP the evolutionary learning model of FeaSANNT [START_REF] Castellani | A technical report on the evolutionary feature selection for artificial neural network pattern classifiers[END_REF].

FeaSANNT is an evolutionary procedure for simultaneous solution of the two combinatory search problems of feature selection and parameter learning in artificial neural network classifier systems. FeaSANNT was already successfully applied to feature selection on traditional artificial neural network pattern classifiers [START_REF] Castellani | A technical report on the evolutionary feature selection for artificial neural network pattern classifiers[END_REF]. In this paper we extend previous work by applying the JASTAP model to the same type of classification problems.

An additional motivation for this work relies on the need to encode recurrence and context in neural networks. The simple MLP (and the related popular backpropagation learning algorithm [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF]), can not handle cyclic connections. Indeed, the well known usage of MLPs as universal approximators can only be achieved by means of enlarging the hidden layer neurons as needed (e.x. Kürková's Universal Approximators article in [START_REF] Arbib | The Handbook of Brain Theory and Neural Networks[END_REF]). As a result, in problems where context or time representation is needed, repetitive iterations can not be used. Spiking neurons in general and JASTAP model in particular may present a solution for this problem.

The remaining of this paper is organized as follows. Section II presents the JASTAP spiking neural model. Then the learning procedure used is presented in section III. The changes made to adapt this procedure and some advantages and properties specific for JASTAP networks are described in section IV. Experimental results are presented in section V. Finally, the contributions of this work are presented in section VI.

II. JASTAP

A. JASTAP As a Spiking Neuron Model

JASTAP 1 model aims to simulate some biologically realistic functions of neural networks [START_REF] Janco | Modeling of neuronal functions: A neuronlike element with the graded response[END_REF]. JASTAP belongs to the family of spiking models. Formal spiking neuron models in general work with temporal coding with mostly biologically relevant action potentials. Spiking models can also perform complex non-linear classification [START_REF] Bohte | Error-backpropagation in temporally encoded networks of spiking neurons[END_REF]. This can be achieved with less units (spiking neurons) than with classical rate-coded networks (please see section V-C).

In this section we briefly present JASTAP model [START_REF] Janco | Modeling of neuronal functions: A neuronlike element with the graded response[END_REF] and the parameters used for easier adaptation of this model to classification tasks. More details on the JASTAP model can be found in [START_REF] Janco | Modeling of neuronal functions: A neuronlike element with the graded response[END_REF].

B. Model Description

A JASTAP model is an artificial NN, which consists of JASTAP neurons as the basic elements. A neuron is described with:

• SET OF SYNAPSES As it is usual in any neural network model, an JASTAP neuron is interconnected with its environment by one or more synaptic inputs and a single output (axon). The output can be connected with one or more neuron synapses in the network.

For each synapse we consider:

-INPUT -can be internal (connected with axon of other neuron) or external (from the outer environment). -SHAPE OF POSTSYNAPTIC POTENTIAL (PSP) PRO-TOTYPE -Biological neuron waveform evoked by a spike arriving at a synapse is described in JASTAP model by

PSP(t) = k • 1 -e -t t 1 2 • e -2t t 2 (1)
The waveform inter alia (controlled by parameters t 1 and t 2 ) emulates whether the synapse is located on a soma or on a dendritic tree. t 1 and t 2 can vary from synapse to synapse. They determine the potential decay. For example, [START_REF] Redman | The time course of synaptic potentials evoked in cat spinal motoneurones at identified group ia synapses[END_REF] mentions t 1 = 0.3 ms and t 2 = 2.7 ms. As the neuron carries out the time-and-space summation of the input potentials and due to time discretisation during simulation, a more moderate decay can cause more sophisticated information processing. Therefore, we used t 1 up to 5 ms and t 2 up to 15 ms (please see figure 1). -LATENCY -the time delay of the synaptic transmission and the axonal conduction. In other words, synapse is silent for the latency time before it starts to translate the input impulse into an action potential. -SYNAPTIC WEIGHT (SW) -a value from -1, 1 that represents the strength of the synaptic input. An excitatory PSP belongs to a synapse with a positive SW and an inhibitory PSP to a synapse with a negative one. -PLASTIC CHANGES -Depending of the synapse type, SW can be influenced by some mechanisms, for instance Hebbian learning or heterosynaptic presynaptic mechanism. In future this can be used to enable online learning. In this work this feature of the model is not used, because the learning during information processing phase is not our goal.

• INSTANTANEOUS MEMBRANE POTENTIAL (MP) -a quantity within the -1, 1 range, determined as the sum of PSPs limited by the non-linear function (figure 2)

MP(t) = 2 π • atan synapses PSP(t) (2)
• THRESHOLD -θ, a value from 0, 1 -that determines the limit for firing. • SPIKE FREQUENCY -the spike frequency is restricted by the absolute refractory period. This is managed by setting minimum I min and maximum I max inter-spike interval for the firing pattern. The standard value we use is I min = 1 ms for the lowest and I max = 10 ms for the highest value. The actual inter-spike interval I a is determined as:

I a = I max -(I max -I min ) • 2 π • atan MP -θ 1 -MP (3)
The spike frequency condition does not allow the neuron to fire sooner, even if the MP exceeds the threshold.

C. Tasks Where JASTAP Has Already Succeeded JASTAP was primarily designed to model real information flow in human brain. Until now, most experiments with JASTAP intend to prove this behavior. On one hand, [START_REF] Pavlasek | Temporal coding and recognition of uncued temporal patterns in neuronal spike trains: biologically plausible network of coincidence detectors and coordinated time delays[END_REF] shows that the JASTAP networks can handle (even noised) information coded in temporal patterns. On the other hand, [START_REF] Pavlasek | Rate coding: neurobiological network performing detection of the difference between mean spiking rates[END_REF] shows its ability to recognize and distinguish different spike rates. This is particularly interesting since the standard MLP weights can be seen as the mean spiking rate concept.

JASTAP model also presents several capabilities that are related with the Gamma distribution [START_REF] Valko | Evolving neural networks for statistical decision theory[END_REF]. Gamma distribution was chosen for experiments because of biological evidence of its plausibility [START_REF] Koch | Biophysics of Computation: Information Processing in Single Neurons[END_REF]. The results showed that the JASTAP models are able to make decisions about features (mean rate, coefficient of variation) of Gamma distribution even with very few neuron units. Moreover [START_REF] Valko | Evolving neural networks for statistical decision theory[END_REF] explored evolved networks and described at low-level how the actual decisions are made.

III. FEASANNT

A. What Is FeaSANNT?

One of the major problems in the use of ANNs is to train the frequently large set of parameters (that is, the connection weights). Because most of ANN weight training procedures are based on gradient descent of the error surface, they are prone to sub-optimal convergence to local minima. Due to its large number of parameters, JASTAP model is particularly sensible to this problem. Also, since standard backpropagation learning [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF] is still an open issue in spiking neural networks and in JASTAP neural networks in particular, we have applied a more powerful learning model for this problem.

Global search techniques such as evolutionary algorithms (EAs) or simulated annealing are known to produce more robust results in difficult search spaces such as the one presented by JASTAP model. FeaSANNT is an evolutionary procedure for simultaneous solution of the two combinatorial search problems of feature selection and parameter learning for artificial neural network classifier systems [START_REF] Castellani | A technical report on the evolutionary feature selection for artificial neural network pattern classifiers[END_REF]. Because FeaSANNT incorporates the use of an embedded approach,

B. Feature Selection in an Embedded Approach

Feature selection can be regarded as a search problem in the discrete space of the subsets of data attributes. Unfortunately, due to the often large set of attributes and their interactions, selecting the optimal feature vector is usually difficult and time consuming. Once again, the search space is noisy, complex, nondifferentiable, multimodal and deceptive and the results are strongly related to the type of classifier system used. However, as pattern classification is based on the information carried by the feature vector, feature selection has a crucial impact on the classifier accuracy, learning capabilities and size. The problem of defining the feature vector stems from the uncertainty about which set of attributes best characterises the data patterns. Finding the optimal one requires removing unnecessary, conflicting, overlapping and redundant features in order to maximise the classifier accuracy, compactness and learning capabilities. Due to the often large number of attributes, searching exhaustively for the best subset is impractical.

Most of current feature selection algorithms either select variables by ranking them according to statistical, geometrical or information theoretical measures or directly evaluate subsets of features by running the induction algorithm and taking the classification accuracy as fitness measure. While the former approach is faster but liable to unforeseen failure as it does not take into account the classifier system, the latter usually implies extensive trial and error guided by hill climbing of the optimisation landscape [START_REF] Blum | Selection of Relevant Features and Examples in Machine Learning[END_REF]. Moreover, due to the locality of the search method, algorithms following the second approach are likely to be led astray by noise or to get trapped into secondary peaks of performance.

Two main approaches are customary for selection of the input vector of features, the filter approach and the wrapper approach [START_REF] Blum | Selection of Relevant Features and Examples in Machine Learning[END_REF]. The former pre-processes (filters) the initial set of attributes based on desirable properties such as orthogonality or information content while the latter evaluates each candidate feature vector on the learning results of the classifier. Feature extraction routines embedded in the learning algorithm (embedded approach, [START_REF] Blum | Selection of Relevant Features and Examples in Machine Learning[END_REF] were also reported even though their application is so far less common. While the filter approach is prone to unexpected failures since it completely ignores the induction algorithm, the wrapper approach involves a far more severe computational effort that may hinder adequate exploration of the search space. Moreover, since the learning results are often dependent on the initialisation of the ANN weights, the latter approach must take into account possibly inaccurate evaluations of the candidate solutions.

FeaSANNT implements an embedded approach in an evolutionary feature selection paradigm. The search is guided by using GAs for learning how to mask inputs. Including the selection of the input features into the learning algorithm of the classifier avoids the computational overheads of repeating the whole ANN training procedure for every evaluation of an individual.

This paper reports on the extension of FeaSANNT model to handle JASTAP Spiking Neural Networks. FeaSANNT's embedded approach was found as particularly useful when a genetic algorithm was used for evolution as well as for learning. This way we were able to use the global nature of the evolutionary search to avoid being trapped by suboptimal peaks of performance while learning both the best set of features, neuron threshold, synapse weights and synapse latencies.

C. What has been already done?

Previous experiments with FeaSANNT [START_REF] Castellani | A technical report on the evolutionary feature selection for artificial neural network pattern classifiers[END_REF], shown that the simultaneous evolution of the input vector and the ANN weights allows significant saving of computational resources. FeaSANNT algorithm was applied on the popular multi-layer perceptron classifier. Also several experiments were preformed on six real-world numerical data sets that gave accurate and robust learning results. Significant reduction of the initial set of input features was achieved in most of the benchmark problems considered. Examination of the evolution curves revealed selection of the optimal feature vector takes place at an initial stage of the search. FeaSANNT seemed also to compare well with other classification algorithms in the literature. However, the proposed algorithm entails lower computational costs due to the embedded feature selection strategy.

IV. USING JASTAP NN MODEL FOR PATTERN CLASSIFICATION

A. Expected Benefits

Expected benefits of spiking models, like JASTAP for pattern recognition (together with feature selection), include: a) Smaller network structures: As stated previously, when compared to MLP, JASTAP neurons should require a smaller number of units. There is already some evidence that the decrease in number of neuron units can be seen even in the simple XOR problem [START_REF] Bohte | Error-backpropagation in temporally encoded networks of spiking neurons[END_REF] when using Spike Response Model [START_REF] Gerstner | Time structure of the activity in neural network models[END_REF].

b) Noise filtering: In [START_REF] Pavlasek | Temporal coding and recognition of uncued temporal patterns in neuronal spike trains: biologically plausible network of coincidence detectors and coordinated time delays[END_REF] authors point out that JASTAP is able to extract information with background spiking noise. In section V it will shown how much noise JASTAP model can bear when applied to the IRIS dataset. We will be especially interested in processing data with Gamma noise, as this kind of noise that is biologically relevant [START_REF] Koch | Biophysics of Computation: Information Processing in Single Neurons[END_REF]. c) Results amenable to analysis: A smaller number of units and detailed processing (temporal coding instead of rate coding) give us a chance to take an insight to internal operations at the level of synapses. Simply speaking, we will try to decode what is the evolved network actually doing.

B. Encoding Variables into Spike-Time Patterns

Encoding processed data into time patterns is still an open issue even for Neurophysiology. For example, it is not clear if the temporal [START_REF] Singer | Time as coding space?[END_REF] or rate coding representations play a crucial role in information processing [START_REF] Abbott | Neural codes and distributed representations: foundations of neural computation[END_REF].

We found a first mention to this issue in [START_REF] Bohte | Error-backpropagation in temporally encoded networks of spiking neurons[END_REF]. Authors used receptive fields to transform real-world data into patterns. Although this approach can be neurologically plausible, for our (FeaSANNT) purposes it was decided to use a repetitive encoding when the same inter-spike interval representing the data value is repeated as an input to input neuron. When the mask element is applied to a particular input, that neuron just receives an empty input train.

Every input is associated with one input neuron as an external synapse. Spike potentials for input neurons are precalculated from input data in such a way that the inter-spike lengths are linearly dependent on the data value (in the range of 5-15 ms). If no noise is taken into account interspike intervals corresponding to one input train are the same. When we test the noise handling, each interspike interval receives a value chosen from the Gamma distribution. This kind of noise is considered as more biologically relevant (rather than applying noise to a particular value) because it indicates cumulative noise during information processing over several neurons (fig. 3).

Decoding to output can be simple. The following method is used: Some neurons in the network are designated as the output neurons. When any of the output neurons fire, they are taken as a hot-spot network decision. The meaning of the decision is task-dependent. Thereafter the simulation in progress can be stopped, because the subsequent activity has no impact on the final decision.

C. Designing Fitness Function for JASTAP Classification

During the tests preformed the selection of fitness function was a crucial step when evolving (recurrent) JASTAP networks. The fitness function was crucial not only for the number of iterations for evolution (i.e. the time the network took to learn), but also in effectivity of learning.

In the first place, the simple successful ratio based fitness functions gave very poor results. A problem with these functions is that they cannot distinguish between networks with no response from those which are able to fire (even if not correctly). Furthermore local minima are common in these cases. This is the consequence of the emergence of the individuals responding to all inputs in the same way. To avoid such problems it was created a fitness function that favors the following criteria:

• early or non-silent responses on input pattern can be classified as incorrect. • heterogeneousness: A network responding (correctly) to many patterns from various classes is better that than a network that can correctly respond to only one. (With same overall ratio).

• selectivity: we should evaluate with higher fitness networks responding mostly correctly to one class, even if responding randomly (or not responding) to the other classes. Indeed these networks exhibit selectivity to some class are better that networks which respond to all patterns in the same way (this is often case, when a population is initialized with random weights). • For datasets with many classes, individuals responding in several classes should be evaluated with higher fitness.

A linear combination of the several (sub-)fitness criteria stated above was used for the final fitness function. First, by favoring small differences in correct responses:

c 1 / 1 + ε - i∈classes (correct-ratio[i])
Second, heterogeneousness and selectivity is maintained by emphasizing the minimal fitness, i. e. the fitness based on minimal (with respect to classes) correct responses:

c 2 / 1 + ε -min i∈classes (correct-ratio[i]) 1. 2. 3. 4.
3: Example of encoding data to temporal code: first iris-setosa training example (i.e. [5.1, 3.5, 1.4, 0.2] scaled to [7.2, 11.3, 5.7, 5.4]). Each row represents an input to 1 of 4 input neurons. Data are scaled to 5-15 ms and repeated over time period of 300 ms. Every interspike-interval is noised with ±1 ms of Gamma noise.

Finally, to handle and favor ability to respond in multiclassed datasets, the function gives little credit when the network is able to respond correctly in any pair of classes.

c 3 • i =j∈classes min(correct-ratio[i], correct-ratio[j])
Experiments were carried out with the values

[c 1 , c 2 , c 3 ] = [1, 30 , 7]
.

V. RESULTS AND COMPARISONS

A. IRIS Data Set

In order to evaluate results the JASTAP in FeaSANNT framework was tested on IRIS UCI ML [START_REF] Hettich | UCI repository of machine learning databases[END_REF] well-known benchmark problem. Inputs were converted into temporal code as described in section IV-B.

The Iris data set was firstly presented in a statistical study for determining the classification of three sets of flowers. The classification set has three distinct classes of plants and four numeric features (corresponding to the petal and sepal width and length). A sample of 150 examples were collected (50 for each class). Only one class is linearly separable from the remaining ones, because some data points on the other two classes are intersecting each other (however most of data points are fairly distinct). From statistical evidence only the measures of petal width and length are enough for classification.

In FeaSANNT algorithm only two features were selected with an average accuracy of 94.7. By using the four features with an MLP with backpropagation learning 96.2% of accuracy were achieved.

B. Set-Up

1) Chromosomes: Structure is not evolved so the individual consists of chromosome defining values subject to evolution and the binary mask for feature selection. Gray binary coding was used for the chromosome network values. The evolution of neuron thresholds and synapse weights and latencies is crucial for network function. Following results of [START_REF] Valko | Evolving neural networks for statistical decision theory[END_REF], JASTAP is not evolved for PSP's t 1 and t 2 parameters I min and I max defining bounding firing periods.

C. Results

1) Iris data set:

This experiment on IRIS data set aimed at showing the potential of JASTAP due to the smaller number of neuron units needed for good classification. Therefore the hidden layer was totally dropped leaving only 7 neurons in structure (4 for inputs and 3 for outputs). In addition we interconnected neurons in input layer to enable mutual information exchange and recurrence. The results are shown in table IV. BPrule stands for classical BP algorithm [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF]. ANNT and FeaSANNT refer to the EA algorithm either disabled or enabled feature selection module [START_REF] Castellani | A technical report on the evolutionary feature selection for artificial neural network pattern classifiers[END_REF]. Present results are in V: Noised Iris data set FeaSTAP 2 column.

2) Handling the noise: JASTAP model was designed to be as detailed as it is needed to simulate biorealistic functions in reasonable way. With this in mind, we have conducted experiments to test the noise handling during classification (please see subsection IV-B). To keep on being biologically inspired we have used the noise generated from Gamma distribution. The noise value in cortical spike times is about 1-2 ms [START_REF]Pulsed Neural Networks[END_REF] what is about 10 % due to our encoding. Such a noise can be generated from G (α = 25, β = 0.8) Smaller amounts of noise (at the 1 % level) were also tested.

Gamma distribution generators GS * and GKM1 from [START_REF] Fishman | Monte-Carlo -concepts algorithms and applications[END_REF] were used. These generators use a combination of the acceptance-rejection, composition and inverse transform methods.

This results show that noise levels until the 10 % level don't reduce the classification quality. However the problem becomes more difficult to learn on the 10 % noise level. By increasing the noise level beyond this limits classifier accuracy starts being compromised. It should be noted that noise was applied to each inter-spike interval independently (even in the same input train).

VI. CONCLUSIONS

This paper main goal is to show that the basic perceptron model can be replaced by the biologically realistic JASTAP neural network model for classification tasks. In order to do so, an evolutionary procedure for simultaneous solution of feature selection and for JASTAP neural network parameter adjusting was used.

Preliminary experimental results show that not only JASTAP seems to be able to deal with the same learning problems with smaller neural networks, but also JASTAP unique features enable the insertion of artificially generated noise into the learning set without degrading learning performance. Indeed, 2 FeaSTAP means JASTAP in FeaSANNT results on the Iris standard data set [START_REF] Hettich | UCI repository of machine learning databases[END_REF] use smaller neural networks without compromising accuracy. Also, noise was artificial inserted into training data without any degradation in classification accuracy. However, it was noticed that the feature selection was not so significant as in original FeaSANNT work. The JASTAP network seems to use the spike potentials even from irrelevant inputs just to improve precision. This is probably the reason why 100% precision is achieved on the IRIS data-set. Since this data set is known not to be fully separable, this high value is probably due to a statistical abnormality either in sample data or in the evaluation procedure. Anyway, care should be taken on future experiments to avoid opportunistic overlearning strategies.

Although not discussed in the paper, the behavior of the evolved networks for IRIS classification was also studied. Indeed, it was observed that the networks do the classification by emphasizing differences in the features correlating with the classes petal length and width.

A major problem with JASTAP learning model was the computational time it took to learn. This paper is presenting a learning model for a difficult problem: there are several parameters to tune, time is directly simulated and we are performing feature selection. As a result the learning of more complicated networks or problems causes very slow learning (i.e. the number of iterations per unit of time decreases). In future work we hope to overcome some of these problems by improving the time simulation and by restricting more the parameters that should be learned. Running FeaSANNT procedure in parallel should help to reduce learning times.

Also in future work, we hope to use JASTAP neural networks to represent logic models as neural networks. Recent work presents the so called neuro-symbolic networks [START_REF] Garcez | [END_REF]. These systems provide simple ways to represent predicate logic programs as MLP neural networks [START_REF] Hitzler | Logic programs and connectionist networks[END_REF]. Unfortunately, providing simple and adjustable structures for encoding previous knowledge or algorithms in a neural network is a difficult task (e.x. Siegelmann's Neural Automata and Analog Computational Complexity article in [START_REF] Arbib | The Handbook of Brain Theory and Neural Networks[END_REF]). Probably a more compact and modular JASTAP neural network would provide a good solution for this problem. In such a neural network, the architecture and some set of weights and parameters will be fixed for encoding background information. For that the relations of the JASTAP model with Gamma distribution [START_REF] Valko | Evolving neural networks for statistical decision theory[END_REF] should also be studied.

Several other features of JASTAP model can also be used with advantage for classification tasks. For example the ability for plastic changes (already studied in JASTAP) can be very interesting for online adaptation of the neural network classifier to small changes in the learned model. Because the JASTAP model works with temporal code, we argue that it can more easily discover statistical regularities over longer time in input pattern trains and encode them on the output. This could be useful in classification problems where context is an issue (e.g. [START_REF] Marques | Tagging with small training corpora[END_REF], [START_REF] Castellani | Automatic Detection of Meddies through Texture Analysis of Sea Surface Temperature Maps[END_REF]).

t 1 1 : 2 :

 112 = 0.3 ms, t 2 = 2.7 ms t 1 = 5 ms, t 2 = 15 ms Postsynaptic potential (PSP) with different waveform inter alia. Limiting non-linear function both feature selection and NN parameter learning can be achieved at (relatively) reduced computational costs.

  Table II displays used JASTAP setup.2) Evolutionary Algorithm: Selection is made in an elitism-like manner: Recombined individuals are evaluated with the same examples as their parents and then the best part (a half) of all chromosomes (parents and offspring together) is taken as the next population.

	feature	min value max value
	synapse weight	-1	1
	synapse latency	0 ms	40 ms
	synapse waveform t 1 value	5 ms
	synapse waveform t 2 value	15 ms
	neuron threshold	0	1
	minimal neuron firing period	1 ms
	maximal neuron firing period	10 ms
	maximal network response time	300 ms
	II: Parameter boundaries used for evolution
	feature	value
	trials	1
	generations	100-500
	population size	100
	values crossover rate	0.7
	values mutation rate	0.05
	mask mutation rate	0.05

III: Parameter settings of evolutionary algorithm