Kernel alignment for identifying objective criteria from brain MEG recordings in schizophrenia

Kernel Alignment for Identifying Objective Criteria from Brain MEG Recordings in Schizophrenia Mario Mart´ın1, Javier B´ejar1, Gennaro Esp´osito1, Neus Catal`a2, Ulises Cort´es1, Ferran Vi˜nas3, Josep Tarrag´o3, Emilio Rojo3, and Rafal Nowak4 1Knowledge Engineering and Machine Learning Group. Universitat Polit`ecnica de Catalunya 2TALP Research Center. Universitat Polit`ecnica de Catalunya 3Hospital Benito Menni CASM 4Cl´ınica Tekn´on Abstract The current wide access to data from different neuroimaging techniques has permitted to obtain data to explore the possibility of finding objective criteria that can be used for diagnostic purposes. In order to decide which features of the data are relevant for the diagnostic task, we present in this paper a simple method for feature selection based on kernel alignment with the ideal kernel in support vector machines (SVM). The method presented shows state-of-the-art performance while being more efficient than other methods for feature selection in SVM. It is also less prone to overfitting due to the properties of the alignment measure. All these abilities are essential in neuroimaging study, where the number of features representing recordings is usually very large compared with the number of recordings. The method has been applied to a dataset in order to determine objective criteria for the diagnosis of schizophrenia. The dataset analyzed has been obtained from multichannel magnetoencephalogram (MEG) recordings, corresponding to the recordings during the performance of a mismatch negativity (MMN) auditory task by a set of schizophrenia patients and a control group. All signal frequency bands are analyzed (from δ(1-4Hz) to high frequency γ(60-200Hz)) and the signal correlations among the different sensors for these frequencies are used as features. 1 Introduction According to [Institute of Mental Health, 2013], schizophrenia affects approximately 1% of the world population. Despite the huge development of several techniques to analyze brain’s structure and function (i.e. EEG, fMRI, MEG, MRI, PET, SPECT), no specific biomarkers are available for a diagnostic purpose. While there is vast evidence of anatomical changes and abnormalities in neurotransmitter systems in schizophrenia, psychiatry research in recent years has focused in functional disconnectivity of cortical circuits as a core feature of the illness. Abnormal oscillatory activity throughout different neural populations and its aberrant synchronization may account for the symptoms observed in schizophrenia. Thus, defining specific patterns of synchrony in individuals with schizophrenia should allow us to make a reliable diagnosis. Early pioneering work using quantitative EEG laid the foundation for the application of discrimination and classification of electrophysiological brain patterns in health and disease. [John et al., 1977] introduced the concept of neurometrics, to provide quantitative information about brain activity related to anatomical integrity, developmental maturation, and mediation of sensory, perceptual, and cognitive processes. However, as already indicated by [Hinkley et al., 2009], there are significant advantages of using MEG over EEG to study neural processes, especially when modeling how activity within a cortical field can influence and interact with activity in other parts of the brain. As well as EEG, MEG is able to reconstruct neural activity on the order of milliseconds. However, MEG sensors allows investigators to access oscillatory neural activity in higher frequency ranges (e.g. α, β, γ) than those attainable in both fMRI/PET. 1 In addition, the sampling frequency of data acquisition in MEG (generally greater than 100Hz) is not limited by electrode impedance (as in EEG), permitting the examination of ultra-high frequency brain activity in this modality. Moreover, volume conduction artifacts commonly found in other imaging modalities are significantly reduced in MEG, as structures such as the skull and CSF do not interfere with the propagation of the magnetic fields ([Leahy et al., 1998]). Hence, MEG arises as a useful method to study how neural oscillations behave and relate to clinical manifestations in schizophrenia. 1.1 Previous work Several methods using measures of brain activity and brain anatomy have been suggested for discriminating schizophrenia patients from healthy subjects. These methods are based on features extracted from different types of signals. We will focus in the works that use MEG approach to characterize subjects. In [Georgopoulos et al., 2007] a classifier for separating several conditions (six conditions including schizophrenia) from healthy subjects is described. The dataset used was obtained from around one minute of MEG recordings from 142 subjects without engaging any specific task during the recording. This dataset was obtained by prewhitening the signals, computing the cross-correlation of all the pairs of sensors and normalizing the feature values. An initial set of 30,628 features described each recording, so the authors decided to perform feature selection using a genetic algorithm. Linear Discriminant Analysis (LDA) on the complete dataset was used to assess the fitness function of each feature set in the genetic algorithm. In the end, they reported that a set of 16 features, when used for building LDA classifiers, returned 78.9% percent average accuracy for all conditions (without specifying the accuracy for the different conditions, specially for the schizophrenic condition we are studying) using the leave one out (LOO) testing procedure, and 86.4% accuracy for 10 runs of 80/20% random split cross-validation. The genetic algorithm was parameterized to only consider up to 20 features in the feature set. This reduces the overfitting effect that is prone to happen in small datasets described with a lot of features. This paper shows that MEG recordings carry relevant information for condition discrimination. However, from the methodological point of view, the set of 16 features selected to build the classifiers was obtained applying the genetic algorithm and LDA procedure to the complete set of 142 individuals (not to the 80% of individuals from the training set in the crossvalidation step), which biases the results of the validation. In different papers [Ince et al., 2007, Ince et al., 2008, Ince et al., 2009] describe their analysis of MEG oscillatory patterns when subjects engaged in a working memory task during the recordings. They studied a number of other works and reported that a higher accuracy can be obtained from features extracted from recordings of functional brain activity (77 – 94%) than from resting state brain activity (67 – 76%). The goal of their study was to discriminate schizophrenic patients from control subjects. The dataset had 15 patients and 23 control subjects. From the MEG recordings, event related synchronizations (ERS) and event related desynchronizations (ERD) patterns were extracted, accounting for an increase and decrease of amplitude in the rhythmic activity in different frequency bands and time points. Data is filtered between 1 and 48Hz and separated in 8 not standard frequency bands: 4 with 4Hz wide and 4 with 8Hz wide. The patterns obtained were used as features for training different machine learning classifiers. Specifically, decision trees ([Ince et al., 2007]), support vector machines (SVM) ([Ince et al., 2008]) and LDA ([Ince et al., 2009]) were used. For the last two methods, feature selection strategies were also used to reduce the number of features of the final classifier. In particular, for LDA the ROC curve for each feature was used to filter the set of attributes, and for SVM the recursive feature elimination algorithm (RFE [Guyon et al., 2002]) was applied. The highest average accuracy using the LOO testing procedure was reported for the LDA classifier, being of 94.5%. It is important to notice that, when selecting a classifier method, the results in LOO help to find the best parameters of the learning algorithm, but they are not a true indication of the generalization ability of the classifier. An unseen validation test set has to be used in order to measure the accuracy on unseen data. So, results are not indicative of percentage of success in clinical application. In all three cases, data was obtained assuming a set of hypothesis on the schizophrenic conditions, for instance that it can be discriminated in ERS/ERD episodes, or the separation in bands. [Escudero et al., 2013] used resting state MEG recordings to obtain a classifier of schizophrenic patients. The dataset was composed by 15 patients and 17 control subjects. The sensors were divided in five anatomical meaningful areas and the average of the power spectral density, Median Frequency, Spectral Entropy, and 2 Relative Power of the sensors of each area was computed and used as features. Logistic Regression (LR) was used for building the classifier. It is reported that subset feature selection was used to reduce the number of features, but the methodology used was not specified. The mean accuracy reported for the final classifier was of 71,3%. Again, accuracy is computed as the average of a 5-fold cross-validation testing procedure, not on a validation data set. 1.2 Plan of the paper In the current paper, we suggest how to take profit of MEG information for diagnostic purposes, and how to identify some MEG characteristics in the brain relevant to discriminate individuals among a sample of patients with chronic schizophrenia (stable compensate or acute exacerbation state) and healthy controls. In order to achieve these goals, we propose the use of state-of-the-art classification procedures from the machine learning area such as support vector machines (SVM) together with a proposal for feature selection to use with support vector machines that shows some advantages over other methods in terms of reduction of overfitting, speed and need of extra set of samples for stopping the feature selection procedure. The paper is organized as follows: In §2 we describe the materials and the methods to prepare data used in this research. In §3 we describe how Support Vector Machines (SVM) can be used to discriminate between patients and control individuals in the data gathered for this task. In the same section we explain how we set the parameters for the SVM and the initial results obtained. In §4 we describe our novel procedure for feature selection in support vector machines. In §4.1 we analyze and discuss the results obtained using this new feature selection procedure. In §5 we present the validation of the approach. Finally, in §6 we discuss the experimental results and give our conclusions and propose our future lines of research. 2 The dataset A total of 30 subjects were recruited over a 12 month period: 20 patients meeting DSM-IV diagnostic criteria for a chronic psychotic disorder (schizophrenia or schizo-affective disorder) and 10 healthy control subjects. All patients were receiving care at Granollers Mental Health Center from Vall`es Oriental Area (Barcelona) and their diagnosis was confirmed by clinical psychiatrist in charge at least 2 years before recruitment, claiming a longitudinal diagnosis. Their psychopathological status was assessed on the day of the experiment by a psychiatrist by means of the Positive and Negative Syndrome Scale (PANSS). On the other hand, healthy volunteers were picked from the local area, after screening by the structured Clinical Interview for DSM-IV Non-Patient version (SCID-NP) to exclude any history of DSM-IV Axis I/II diagnosis or substance abuse, with an additional exclusion criterion of any first or second-degree biological relatives with psychotic disorder diagnosis. In order to obtain a good representativeness of the heterogeneity of chronic psychotic patients, we aimed to recruit a sample including patients through different stages of the illness (acute vs stable patients). Therefore, the sample was divided in: 1. Ten patients with a chronic psychotic disorder diagnosis, suffering an acute psychotic episode at the moment of the study, scoring PANSS >27 (acute patients). 2. Ten patients with a chronic psychotic disorder with no current relevant positive psychotic symptomatology at the moment of the study (PANSS <14) (stable patients). 3. Ten healthy subjects (control group). Exclusion criteria for all participants included neurological diseases, substance dependence, a history of head injury or a full-scale IQ estimation of less than 70. All subjects proved to be right-handed according to modified version of the Edinburgh Handedness Questionnaire ([Oldfield, 1971]), which asks subjects to demonstrate hand use on various actions. Prior to the study, all subjects were fully informed and gave written consent to participate. This study was approved the Ethics commission of Clinical Research (CEIC) of Hospitals of Hermanas Hospitalarias in Barcelona in accordance with the Declaration of Helsinki protocols. 3 Table 1: Demographic and clinical data of study groups. Stable Disorder Acute Exacerbation Disorder Healthy Control Gender 9 Male/ 1 Female 8 Male/ 0 Female 9 Male/ 1 Female Age (years) 36.7 (±10.7) 37.5 (±13.3) 36.7 (±10.7) Education (years) 10 (±1.7) 9.60 (±2,2) 11 (±1.8) Distribution Diagnosis axis I 7 Schizophrenia disorder 6 Schizophrenia disorder Not applicable (DMS-IV) 3 SchizoAffective disorder 2 SchizoAffective disorder Duration of illness (years) 13.8 (±7.5) 17.5 (±13.3) Not applicable PANSS-P 10.5 (±2.6) 30.3 (±4) Not applicable PANSS-NN 15.4 (±4.8) 19.9 (±3) Not applicable All subjects were Spanish natives or bilingual Catalan/Spanish speakers and they were matched by gender, age and educational degree through the three groups. To guarantee the homogeneity of the groups regarding these characteristics, subjects were recruited according to a sequential design, starting by the acute psychotic patients, following by the stable group, and finally enrolling the control group. Demographic and clinical data for all subjects are presented in Table 1. At the time of MEG recordings acquisition, all patients were using atypical antipsychotic medication. Two patients from the acute psychotic group could not complete the MEG recording due to psychomotor agitation during the tasks, so their recordings were not considered for the purpose of the study. The data used in the experiments corresponds to MEG recordings acquired with a 148-channel whole-head magnetometer (MAGNES 2500WH) in a magnetically shielded room (Vacuumschmelze GmbH, Germany) at the Teknon Hospital in Barcelona, Spain. These were acquired while the subjects were lying on a patient bed with their head inside the helmet-like MEG device. Activity was recorded during an auditory oddball paradigm to assess the mismatch negativity (MMN) response. Acoustic stimuli were generated using the E-Prime software (Psychology Software Tools, USA) and were presented binaurally through Etymotic ER-30 (Etymotic Research, Inc. USA) non-magnetic earpieces. Subjects were instructed to relax and ignore the auditory stimuli during the task. The oddball paradigm was based on sequences of two tones, each starting on a pseudorandom basis, either with a 1000-Hz standard (p=0.8) tone followed by a 800-Hz deviant (p=0.2) tone. Tone duration was 100 ms, with rise and fall times of 10 ms, and the inter-stimulus interval was 350 ms. The recording frequency was 678.19Hz (band-pass 0.1-250Hz). Five position coils were attached to the forehead and to the periauricular points in order to determine the position of the head and to track any head movement during the recording. Data sets in which the relative position of the head changed by >0.5cm throughout the recording session were discarded from further analysis. For each subject, the headshape including the forehead, the nose, and the location of the sensor position coils were digitized using a digitizer wand (Polhemus Fastrak, Polhemus Inc., USA). From the MEG data recordings, six different frequency bands were extracted, using an 8th-order Butterworth digital IIR filter: δ[1-4 Hz], θ[4-8 Hz], α[8-13 Hz], β[13-30 Hz], lower γ [30-60 Hz], and upper γ [60-200 Hz]. The upper γ band was included given that several studies indicate that synchronizations at this frequency range can be linked to schizophrenia (see for example [Uhlhaas et al., 2011]). Studies without band separation were also performed. 2.1 Data Preparation One objective of this study is to find out if there is a set of features, in terms of MEG recording of the brain, for determining an objective, efficient and safe way whether an individual has schizophrenia or not. This set of attributes has to be easy to interpret in terms of brain areas and their relationships. In order to explore this hypothesis, first we must extract from the raw data a suitable group of features that allow the separation of subjects suffering schizophrenia from other ones. We want these features to represent the connectivity of different areas of the brain. They also have to capture the global behavior of the signals independently of their magnitude. Different feature extraction methods that extract synchronization information can be used for MEG signals, being the most common those oriented to the frequency domain like coherence ([Srinivasan et al., 2007]) or to the time domain like 4 s s s Figure 1: Signals from sensors MEG148, MEG91 and MEG61 (from left to right) with 1000 time points length. Correlation is higher for more synchronized signals (corr(MEG148,MEG91)=0.66, corr(MEG148,MEG61)=- 0.008). MEG148 and MEG91 correspond to right side temporal sensors and MEG61 corresponds to an occipital sensor. covariance ([Shenoy et al., 2006]). Given that one of the goals of was to study separately different frequency bands, we judged more adequate to extract the characteristics from the time domain. Among the different alternatives, the decision was to use the Pearson product-moment correlation coefficient because it is a simple way to measure global similarity among the different sensors and normalizes the contribution of each pair of sensors. The coefficients capture this effect as a coincidence of the statistical distribution of the signals as can be seen in figure 1. Subject MEG Recording Auditory Task Band Separation Feature Extraction Correlation Coefficient Dataset . . . . Subject 1 Subject n Window Segmentation Window Correlation Coefficient Averaging Figure 2: Data preparation process. We wanted to use relatively large time periods of the signal (around a minute), so the synchronizations among different areas are captured by the features, capturing a fingerprint of the brain activity. Also, It was needed to reduce the variability of the measures and to obtain characteristics that were not affected by artifacts. Thus, an additional preprocess was performed, consisting in the segmentation of the 10 minute recordings in several windows of one minute to obtain measures at different time segments of the experiment. We proceeded to average the correlation coefficient matrices of these windows to smooth the behavior along the time, obtaining this way a better estimate of the mean behavior of the sensor synchronizations. This procedure also reduces the effects of possible artifacts during the recording, allowing working directly with the raw data. Data preparation for a specific frequency band was performed in the following way (see figure 2): 1. Extract a specific band to study from a file of a subject. 2. Separate the total length of the signal in k= 10 windows of the same size. 3. Compute the correlation coefficient among all pairs of sensors for each window ρxy =covxy σxσy =Pi(xi−¯x)(yi−¯y) pPi(xi−¯x)2pPi(yi−¯y)2 5 4. Average the kcorrelation matrices of all windows building the global correlation matrix for the whole signal In the resulting dataset, separated for each band, each subject file is described by a number of continuous attributes in the range [−1,1] as n2−n 2measures, where nis the number of sensors used in the recordings. For the learning task a random training subset of the data was selected, consisting in 7 control subjects, 7 compensated patients and 6 non compensated patients, for a total of 20 examples. On this data we will study the bands and parameters for the learning algorithm using the standard leave-one-out procedure. As we are concerned with overfitting, we will finally validate the approach selected using the remaining 8 recordings: 3 control subjects, 3 compensated patients and 2 non compensated patients. Given that the number of individuals in training and testing is very small, we explore bootstrapping methods for validation and to obtain a valid p-value for the statistical point of view (see section 5.1). 3 Experimental procedure As stated in §2, preliminary experiments are directed to test whether the dataset has enough information to separate control individuals from the ones with schizophrenia diagnostic. Machine learning literature offers a plenty of tools to perform this task. Nevertheless, the problem under study shows a key feature which reduce possible choices. On one hand, recordings in the data set are described by correlations among 146 (of 148, due to two reported malfunctioning sensors) pairs of sensors entailing 10,585 correlations. On the other hand, we only have 20 recordings. So, we have an extremely small number of recordings described by a lot of features. We have to choose a method that, while showing a state-of-the-art accuracy in most benchmarks, it is also able to deal with such kind of dataset. Support Vector Machines (SVMs) seem a promising approach because they have shown the ability to deal with large dimension spaces compared with the number of instances in several domains, for instance in document classification ([Joachims, 1998]). 3.1 Classification using Support Vector Machines In this section we define the notation that we will be using to describe SVMs and related expressions. Suppose having a set of ndata points xi∈Rdalong with each point’s classification can take on one of two possible values yi∈ {±1}. The linear SVM is defined considering the best hyperplane hx,wi+b= 0 separating the points in two different classes. It is often useful to consider the SVMs in its dual formulation: min α 1 2Pn i,j=1 yiyjK(xixj)αiαj−Pn j=1 αj(1) s.t. Pn i=1 αiyi= 0 0 ≤αi≤C where Kis the selected kernel and C is a parameter usually determined experimentally specifying the trade-off between the empirical error and the complexity term. Once this problem has been solved the coefficients (αi, b) can be used to classify the test point xbased on separating two classes with the hyperplane trough f(x) = sign( n X j=1 αjyjK(xj,x) + b) (2) which classifies xas either +1 or −1. Interpreting the Equation 2, one understands that the dual variables αirepresent the importance of a particular point (Support Vectors) within the model. A high value of αiassociated with a particular support vector xiindicating that the kernel value calculated using this support vector will have more influence on the final value for f(x). We will consider in this paper the αvalues of the resulting solution of the SVM, because it is know that the number of supports (αdifferent of zero) can be used as a bound to the empirical LOO error (see next section). 6 Table 2: Accuracy, Recall and Precision in separation from control to patient diagnosed for each band. Last column shows number of support vectors and proportion of support vectors. See §3.3 for a complete discussion. Band Accu. Recall Prec. Sups. (perc.) Full Band 80 57 80 13 (0.68) δ90 86 86 13 (0.66) θ50 14 20 12 (0.61) α60 29 40 14 (0.75) β55 29 33 15 (0.78) lower γ 80 57 80 14 (0.74) upper γ 80 57 80 13 (0.68) 3.2 Parameters selection and initial results SVMs performance depends on the selected kernel, its parameters and, eventually, on some data preparation. Within our experiments using MEG data, the only pre-process is a standard data normalization, common in the learning of SVMs. On the one hand, selecting a kernel can be done experimentally or using some heuristics, and usually requires some cross-validation technique, to study the possible kernels and its parameters. In order to avoid overfitting in our small data set, we limit our choices by using a linear kernel, which is the simplest approach and does not require further parameters to be tuned. Intuitively, it should be enough for our discriminative task for our dataset, due to the large number of features used to describe the data when compared with the number of recordings. Although SVMs often produce effective solutions for balanced datasets, they could be sensitive to unbalanced datasets and might produce sub-optimal models. A simple method to overcome this problem consists in modifying the SVM objective function by assigning two different misclassification costs, such that the cost of misclassifying training instances of the minority class becomes higher than the penalty for misclassifying instances of the majority class ([Cawley and Talbot, 2001]). A rule of thumb working reasonably well for classification tasks is to set penalty C−for the negative class (assuming that negative class is the majority class) and C+=fcC−for the positive class, where fcis equal to the number of negative examples in the training set divided by the number of positive ones. The effect of this procedure is equivalent to building a SVM with an oversampled minority class to match the number of examples of both classes. Our dataset is imbalanced, having 7 control and 13 patient individuals. We set C−= 1 for negative (patient) and C+= (13/7)C−for positive (control) examples. We tested different values of the C−parameter in the interval [0.1,∞] essentially leaving unchanged the obtained results. Finally, in order to evaluate the classification performance we apply a standard Leave one out (LOO) procedure which gives an almost unbiased estimate of the expected generalization error. LOO consists in removing one element from the training data set, constructing the decision rule on the basis of the remaining training data and then testing on the removed element. The procedure is repeated ntimes (where nis the number of instances of the data set), each time leaving a different individual for testing. Our data set, with 20 examples, allow to repeat this procedure 20 times. Results of the 20 runs were averaged. The LOO error (ELOO) is related with the number of supports in SVMs. It is known that the proportion of support vectors with respect to the total number of instances is a theoretical bound of the LOO error. That is, ELOO ≤nSV nwhere nSV is the number of support vectors (see for instance Remark 7.57 in [Shawe-Taylor and Cristianini, 2004]). As a consequence SVMs with a reduced set of support vectors show a lower leave on out error, which means that we expect from such machines a higher accuracy on data unseen in the learning process (better generalization). 3.3 Classification results Results obtained for the correlation computation on the original data, and on data after filtering each band, are reported in table 2. These results cannot be understood only in terms of accuracy. Note that, being this an imbalanced dataset, a classification method saying that all individuals are schizophrenic leads to an 7 accuracy of 65%. Hence, to avoid this kind of misinterpretation, fairer evaluation scores on the smaller class (in this the control case) are also displayed: •Recall: number of control individuals correctly classified with respect to the total of control individuals. It estimates the probability of detecting the control cases. •Precision: number of the cases classified as control that were actually control instances with respect to the number of cases classified by the system as control. Complementary to recall, it estimates the probability of success in determining the class of the cases classified as control. Table 2 shows that using the correlation of sensors. Considering the full spectrum band, we obtain SVMs with an accuracy of 80% to distinguishing between control and patients individuals. It seems that the correlation of sensors carries information useful for the discrimination task. Separating the signal in different frequency bands shows that θ,αand βdo not help to discriminate control from patients. On the other hand, the δand γbands have a higher accuracy, specially in the case of the δband, that recognizes almost all patients and only fails in 2 cases from 20. However, from the point of view of the robustness of the results, it is worth noting that all classifiers show a high rate of support vectors with respect to the number of examples, most of them around 70% (see last column in table). We have mentioned in section 3.1 that LOO error (ELOO) is theoretically bounded by proportion of supports in the SVM. The higher the proportion of supports vectors needed to build the SVM is, the higher are the expectations on error on unseen cases. So, such machines could show in the worst case unreliable predictions on test data. This is a sign that we may have ended with overfitted SVMs for our data set. Given the small number of examples that we have available and the large number of features describing them, this is an issue that we have to check carefully. Also, being encouraging results, the problem within this approach is that they are not easy to interpret in terms of the visualization of areas of the brain possibly related to schizophrenia. The SVM describes the hyperplane separating positive and negative examples with a set of 10,585 weights in the primal form that are not easily interpretable. Nevertheless, most of them may not be necessary to build an effective hyperplane and we will try to simplify the dataset by using feature selection. 4 Feature Selection by incrementing alignment Fortunately, the interpretation of the classification results is actually possible by studying the effect of feature selection on the dataset. Feature selection methods reduce the feature set in the data collection by removing noisy or redundant features, and thus allowing a better understanding of the classification procedure, which in medical datasets is desirable. Several methods have been proposed for this task in the SVM framework. [Weston et al., 2000] proposes to find, via gradient descent, those features which minimize the leave-one-out error bounds, that depend on the radius which includes all vectors and the margin of the learned SVM. [Guyon et al., 2002] presents the method Recursive Feature Elimination (RFE) that from the whole set of features and the learned SVM, removes the feature that shows the smaller weight win primal after learning (and so decreases less the margin). Features are removed iteratively until a given threshold, in terms of weight, is achieved. The main problem with all these approaches is their computational cost. A feature is only detected after learning takes place. This means that many SVMs have to be learned with the associated cost. Also, those methods rely on cross-validation to adjust the parameters for stopping the removing of features. We propose a different approach that, while taking advantage of the SVM peculiarities, does not require the previous building of SVMs and has also a natural threshold for stopping the feature removal. This approach is based on a measure known as empirical kernel-target alignment ([Cristianini et al., 2002]). SVMs find the maximum margin hyperplane separating data with different labels. To do that, they only need the dot product of the data in the original space or in the feature space. For this, they use the kernel function and the kernel matrix over the data set. A way to measure how a kernel helps to separate the data is by comparing it with the so called ideal kernel:T=yy0. This is a n×nmatrix that by definition presents Ti,j = +1 when xiand xjhave the same label, and −1 otherwise. 8 Empirical kernel alignment measures the fitness of a kernel Kto training labels [Cristianini et al., 2002]: A(K, T ) = hK, T iF phK, KiFhT, TiF (3) where hM, NiFis the Frobenius product of matrices Mand N. Alignment is a measure ranging in [0,1] (the higher the better), having a number of convenient properties, specially (1) it can be evaluated before learning and (2) it is a sharply concentrated statistic around its expected value (stable to different splits of the training data and not prone to overfitting). The last feature is specially relevant in our case because, given the small number of examples and large number of features in our dataset, we want to avoid an overfitting of the results. Alignment has been used for a priori kernel selection (choosing the best kernel from a set of candidates) and for kernel adaption ([Lanckriet et al., 2004]). We will use this a priori model selection procedure to find a suitable set of features, facilitating the learning process. We define kernel KFfrom the set of features Fas the kernel matrix obtained by applying the selected kernel function on the data represented only with features in set F. The idea is to implement feature selection using the empirical kernel-target alignment, using the fact that choosing one set of features or another is equivalent to choosing one kernel or another. Hence, we select the set of features Fdefining the kernel KFwith the highest empirical kernel-target alignment. In order to deal with unbalanced data sets, we used a modified version of alignment ([Kandola et al., 2002]) in which labels to calculate the ideal kernel are changed from +1 to +1/m and from −1 to −1/n, where m and nare the number of positive and negative examples respectively. The learning of the SVM is done with the original set of labels. For feature selection, we implemented an iterative procedure that, starting from the empty set, incrementally adds the feature increasing the most the alignment with the ideal kernel. This stops when no increment in alignment is found. The advantages over other approaches as RFE ([Guyon et al., 2002]) are (1) that our method avoids the costs of building SVMs to consider relevance of features, and (2) that the addition of features naturally stops when alignment is not increased, that is, no validation dataset is required to determine when to stop. However, the cost of the method is still expensive. The cost of computing the linear kernel KFis O(n2f), where nis the number of examples and fis the number of features in F. The kernel has to be calculated for each set of features before computing its alignment. In addition, the number of features fin our dataset is very large, so removing the fterm is desirable. A careful examination of the problem shows that there is an efficient incremental way to compute the kernel after adding a new feature from the same kernel without that feature. This reduces the cost of computing the new kernel for our incremental algorithm from O(n2f) to O(n2). For the linear kernel, the following equation shows how to quickly calculate the new kernel from the old kernel when only adding a new feature i: KF∪{i}(x, y) = X j∈F∪{i} xjyj=X j∈F xjyj+xiyi=KF(x, y) + xiyi with constant cost for one pair of examples. Hence computing the kernel KF∪{i}from KFonly has a cost of O(n2). We will take advantage of this fact in our algorithm. This procedure can be easily extended for non linear kernels, like polynomial or the RBF kernel. In this way, we can reduce also to O(n2) the cost of computing the alignment A(KF∪{i},T) from KF. For large data sets this could still be a problem. However, alignment is a sharply concentrated measure ([Cristianini et al., 2002]) and could be accurately estimated from a sample of the whole training set. With this in mind, the method of incremental Feature Selection by Alignment can be detailed as shown in algorithm 1. The algorithm starts from the empty set of features F. At each iteration it always looks for the feature which increases the alignment the most. This feature is added to Fand the procedure is repeated until no more features remain in the set of candidate features Cincreasing the alignment. Note that after selecting one feature, that feature is removed from set C. The computational time complexity of line 7 is n2. Thus, loop of lines 6-8 has complexity n2|C| ≤ n2f, and thus, the outer loop (lines 5-16) has complexity n2ffa(where fais the final number of features of set F). So, the final computational complexity of the feature selection procedure is then O(n2ffa). 9 [Cristianini et al., 2002] Cristianini, N., Shawe-Taylor, J., Elisseeff, A., and Kandola, J. (2002). On KernelTarget Alignment. In Dietterich, T. G., Becker, S., and Ghahramani, Z., editors, Advances in Neural Information Processing Systems 14. MIT Press. [Escudero et al., 2013] Escudero, J., Ifeachorl, E., Fern´andez, A., L´opez-Ibor, J. J., and Hornero, R. (2013). Changes in the meg background activity in patients with positive symptoms of schizophrenia: spectral analysis and impact of age. Physiological measurement, 34(2):265. [Georgopoulos et al., 2007] Georgopoulos, A. P., Karageorgiou, E., and et al. (2007). Synchronous neural interactions assessed by magnetoencephalography: a functional biomarker for brain disorders. J. Neural Eng., 4:349–355. [Guyon et al., 2002] Guyon, I., J Weston, a. S. B., and Vapnik, V. (2002). Gene selection for cancer classification using support vector machines. Mach. Learn., 46(1-3):389–422. [Harvey et al., 2012] Harvey, P. D., Heaton, R. K., Jr., W. T. C., Green, M. F., Gold, J. M., and Schoenbaum, M. (2012). Diagnosis of schizophrenia: Consistency across information sources and stability of the condition. Schiz. Research, 140(1–3):9 – 14. [Hinkley et al., 2009] Hinkley, L. B. N., Owen, J. P., Fisher, M., Findlay, A. M., Vinogradov, S., and Nagarajan, S. S. (2009). Cognitive Impairments in Schizophrenia as Assessed Through Activation and Connectivity Measures of Magnetoencephalography (MEG) Data. Front Hum Neurosci., 3(73). [Ince et al., 2008] Ince, N., Goksu, F., Pellizzer, G., Tewfik, A., and Stephane, M. (2008). Selection of spectrotemporal patterns in multichannel meg with support vector machines for schizophrenia classification. In 30th Annual International Conference of the IEEE on Engineering in Medicine and Biology, pages 3554– 3557. [Ince et al., 2007] Ince, N., Stephane, M., Tewfik, A., Pellizzer, G., and McClannahan, K. (2007). Schizophrenia classification using working memory meg erd/ers patterns. In 3rd International IEEE/EMBS Conference on Neural Engineering, pages 457–460. [Ince et al., 2009] Ince, N. F., Pellizzer, G., Tewfik, A. H., Nelson, K., Leuthold, A., McClannahan, K., and Stephane, M. (2009). Classification of schizophrenia with spectro-temporo-spatial MEG patterns in working memory. Clinical Neurophysiology, 120:1123–1134. [Institute of Mental Health, 2013] Institute of Mental Health, N. (2013). [Joachims, 1998] Joachims, T. (1998). Text Categorization with Suport Vector Machines: Learning with Many Relevant Features. In ECML ’98: Proceedings of the 10tth European Conference on Machine Learning, pages 137–142. Springer-Verlag. [John et al., 1977] John, E. R., Karmel, B., Corning, W. C., Easton, P., Brown, D., Ahn, H., John, M., Harmony, T., Prichep, L., Toro, A., Gerson, I., Barllet, F., Thatcher, F., Kaye, H., and Valdes, P. (1977). Neurometrics. Science, 196(4297):1393–1410. [Kandola et al., 2002] Kandola, J., Shawe-Taylor, J., and Cristianini, N. (2002). On the Extensions of Kernel Alignment. Technical Report NC-TR-02-120, NeuroCOLT. [Lanckriet et al., 2004] Lanckriet, G. R. G., Cristianini, N., Bartlett, P., Ghaoui, L. E., and Jordan, M. I. (2004). Learning the kernel matrix with semidefinite programming. J. Mach. Learn. Res., 5:27–72. [Leahy et al., 1998] Leahy, R., Mosher, J., Spencer, M., Huang, M., and Lewine, J. (1998). A study of dipole localization accuracy for MEG and EEG using a human skull phantom. Electroencephalogr Clin Neurophysiol., 107(2):159–73. [Oldfield, 1971] Oldfield, R. (1971). The assessment and analysis of handedness: the Edinburgh Inventory. Neuropsychologia, 9:97–113. [Shawe-Taylor and Cristianini, 2004] Shawe-Taylor, J. and Cristianini, N. (2004). Kernel Methods for Pattern Analysis. Cambridge University Press, New York, NY, USA. 16 [Shenoy et al., 2006] Shenoy, P., Krauledat, M., Blankertz, B., Rao, R. P. N., and M¨uller, K.-R. (2006). Towards adaptive classification for bci. Journal of Neural Engineering, 3(1):R13. [Srinivasan et al., 2007] Srinivasan, R., Winter, W. R., Ding, J., and Nunez, P. L. (2007). EEG and MEG coherence: Measures of functional connectivity at distinct spatial scales of neocortical dynamics. Journal of Neuroscience Methods, 166(1):41 – 52. [Uhlhaas et al., 2011] Uhlhaas, P. J., Pipa, G., Neuenschwander, S., Wibral, M., and Singer, W. (2011). A new look at gamma? high-(>60 hz) γ-band activity in cortical networks: function, mechanisms and impairment. Progress in biophysics and molecular biology, 105(1):14–28. [Weston et al., 2000] Weston, J., Mukherjee, S., Chapelle, O., Pontil, M., Poggio, T., and Vapnik, V. (2000). Feature selection for SVMs. In NIPS, pages 668–674. 17