Evolutionary Generalized Radial Basis Function neural networks for improving prediction accuracy in gene classification using feature selection

Evolutionary Generalized Radial Basis Function neural networks for improving prediction accuracy in gene classification using feature selection Francisco Fernández-Navarro , César Hervás-Martínez , Roberto Ruiz , Jose C. Riquelme Keywords: Generalized Radial Basis Function Generalized Gaussian Distribution Evolutionary algorithm Gene classification Feature selection a b s t r a c t Radial Basis Function Neural Networks (RBFNNs) have been successfully employed in several function approximation and pattern recognition problems. The use of different RBFs in RBFNN has been reported in the literature and here the study centres on the use of the Generalized Radial Basis Function Neural Networks (GRBFNNs). An interesting property of the GRBF is that it can continuously and smoothly repro-duce different RBFs by changing a real parameter . In addition, the mixed use of different RBF shapes in only one RBFNN is allowed. Generalized Radial Basis Function (GRBF) is based on Generalized Gaussian Distribution (GGD), which adds a shape parameter, , to standard Gaussian Distribution. Moreover, this paper describes a hybrid approach, Hybrid Algorithm (HA), which combines evolutionary and gradient-based learning methods to estimate the architecture, weights and node topology of GRBFNN classifiers. The feasibility and benefits of the approach are demonstrated by means of six gene microarray classi-fication problems taken from bioinformatic and biomedical domains. Three filters were applied: Fast Correlation-Based Filter (FCBF), Best Incremental Ranked Subset (BIRS), and Best Agglomerative Ranked Subset (BARS); this was done in order to identify salient expression genes from among the thousands of genes in microarray data that can directly contribute to determining the class membership of each pattern. After different gene subsets were obtained, the proposed methodology was performed using the selected gene subsets as new input variables. The results confirm that the GRBFNN classifier leads to a promising improvement in accuracy. 1. Introduction In traditional RBFNN, the Gaussian function is selected as the network activation function [1,2], although other functional forms have been used for the RBFNNs including some types of thin-plate spline functions, multi-quadratic functions and sigmoidal functions [3]. Nevertheless, there are still some problems in standard Gaussian RBFNN. First, if, the underlying curve representing training pat-terns is nearly constant in a specific interval, it is difficult to utilize a Gaussian function to approximate this constant valued function unless its width tends to infinity. In this case, a RBFNN would be an inefficient model to approximate constant valued functions. Second, in high-dimensional space, all pairwise distance between patterns seem to be very similar, i.e., the distances to nearest and furthest neighbours look nearly identical. Therefore, the widely used Gaussian kernel and Euclidean distance are not necessarily appropriate functions to quantify similarity in high dimensional spaces because distances in this kind of problem are concentrated and the Gaussian kernel loses its interpretation in terms of locality around its centre [4,5]. Despite this, the use of Euclidean distance in high-dimensional space is not questioned in the machine learning community since it corresponds to distance as we define it in our three-dimensional world. Third, several papers have included this as future experimentation to be performed with more general RBF models that achieve a good compromise between low training effort and flexible modelling capabilities [6–9]. In order to take care of these problems, a new activation function is presented in this paper. The proposed RBF is based on Generalized Gaussian Distribution (GGD) [10]. GGD adds a shape parameter to the normal distribution. In the same way that GGD adds a  shape parameter to Gaussian Distribution, the novel RBF proposed, called Generalized Radial Basis Function (GRBF), also adds a shape parameter to standard Gaussian RBF (SRBF). The GRBF allows better matching between the shape of the kernel and the distribution of the distances, since the parameter provokes concavity or convexity around the point where the distance is the radii of the kernel, r. On the other hand, although a great number of algorithms have been developed to estimate the parameters of a RBFNN for a fixed topology, most of them are hill-climbing procedures, which usually fall in a local optimum. Evolutionary algorithms (EAs) have proved to be very effective and robust search methods for locating zones in the search space where good solutions can be found, even if this space is large and contains multiple local optimums. The application of the EA to optimize RBFNN is justified by a smaller overall computational cost, compared to the methods of trial and error, and their robustness as opposed to constructive/pruning methods. In general EAs are less efficient than local search techniques in finding the local optimum, so it is convenient to allow the EA to select initial solutions in good areas of the search space, and to locate local optimum in these areas afterwards. For that reason, a Hybrid Algorithm (HA) was employed to estimate the parameters of the model proposed. The HA combines both a global search procedure, EA, and a local improvement procedure based on gradient descent, the iRprop+ algorithm. Thus, the EA optimizes RBFNN parameters so that RBFNN parameters are located in an area of global optimum and, the iRprop+ algorithm refines its parameters to improve the results reported by the EA. The performance of the proposed methodology was evaluated in six well-known deoxyribonucleic acid (DNA) microarray classification problems. DNA microarray allows relative levels of ribonucleic acid (RNA) or messenger RNA (mRNA) abundance to be determined in a set of tissues or cell populations for thousands of genes simultaneously. The importance of the use of Artificial Neural Networks (ANNs) in the classification of microarray gene expression [11] as an alternative to other techniques was stated in several research works [12,13] due to their flexibility and the high degree of accuracy in their fit to experimental data. These datasets were selected to justify the use of this kernel model in the classification of high dimensionality problems. Furthermore, other soft computing techniques have been implemented to address this problem [14–16]. The motivation for applying feature selection (FS) techniques has shifted from being optional to becoming a real prerequisite for model building. A typical microarray dataset may contain thousands of genes but only a small number of samples (often less than two hundred). Theoretically, having more genes should give us more discriminating power. However, this can cause several problems: increased computational complexity and cost; too many redundant or irrelevant genes; and estimation degradation in the classification error. There are two ways to group feature selection algorithms, depending on the evaluation measure chosen: one, according to the model used (filter or wrapper) and two, according to the way in which the features are evaluated (individually or by subsets). The filter model evaluates features according to heuristics based on overall data characteristics, notwithstanding the classification method applied, whereas the wrapper uses the behaviour of the target classification algorithm as the feature evaluation criterion. Based on the generation procedure, FS can be divided into individual feature ranking (FR) and feature subset selection (FSS) [17,18]. FR measures feature-class relevance, then ranks features by their scores and selects the top-ranked ones. In contrast, FSS attempts to find a set of features that performs well. It integrates the metrics for measuring feature-class relevance and feature–feature interactions. A hybrid model was proposed to handle large datasets to take advantage of the above two approaches (FR, FSS). These methods decouple relevance analysis and redundancy analysis, and have proven to be more effective than ranking methods and more efficient than subset evaluation methods in many traditional high-dimensional datasets. In this framework, FCBF (Fast Correlation-Based Filter) [19], BIRS (Best Incremental Ranked −2 −1.5 −1 −0.5 0 0.5 11.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 τvalues τ= 0.5 τ= 1.0 τ= 2.0 τ= 4.0 τ= 8.0 τ= 16.0 Fig. 1. Probability density function of the Generalized Gaussian Distribution (GGD) with different values of ,c= 0 and r= 1. Subset) [20] and BARS (Best Agglomerative Ranked Subset) [21] are the methods proposed to obtain relevant features and to remove redundancy. These features are considered input variables in the network models (GRBFNN) that we propose in this paper One of the major advantages of the proposed method is the reduced number of features and GRBFs included in the final expression, since the HA reduces its complexity by pruning connections and hidden nodes. This can result in a better interpretability of the model, which is especially important when dealing with real problems. Therefore, using the proposed approach, the feature selection is performed in two stages: firstly, in preprocessing by means of the feature selector and secondly, in the HA by pruning connections. This paper is organized as follows: Section 2 formally presents the GRBF model considered in this work and the main characteristics of the algorithm used for training the model. Section 3 introduces the feature selection algorithms used in this paper. Section 4 describes the experiments carried out and discusses the results obtained. Finally, Section 5 completes the paper with the main conclusions and future directions suggested by this study. 2. Classification method 2.1. Generalized Gaussian Distribution Although Gaussian Distribution has a principal role in statistical applications, the analysis of real data often leads to rejecting the hypothesisthatdatahavebeengeneratedbynormaldistribution.In these circumstances the adoption of more flexible models allowing therepresentationofdatageneratedby distributions inanarea near the Gaussian one may be appropriate. In particular, models which embed Gaussian distribution as a special case are of great interest, since their use permits deviations to be dealt with from normality, while preserving the possibility to test the adequacy of Gaussian distribution to the data. Generalized Gaussian Distribution (GGD) adds only one additional parameter to Gaussian distribution, the shape parameter . The GGD can approximate a large class of statistical distributions by modifying this parameter, for instance: the Gaussian distribution is obtained for = 2, the Laplacian distribution for = 1, and by making →0we can obtain a distribution close to uniform distribution (Fig. 1). The analytical equation for the −2 −1.5 −1 −0.5 0 0.5 11.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 τvalues τ= 0.5 τ= 1.0 τ= 2.0 τ= 4.0 τ= 8.0 τ= 16.0 Fig. 2. Radial unit activation in one-dimensional space with c= 0 and r= 1 for the Generalized RBF (GRBF) with different values of . probability density function of GGD, for a unidimensional input variable is given by pj(x;i, cj, rj)=j 2rj(1/j)exp −x−cjj rj j, j =1, . . . , J (1) where Jis the number of classes in the problem, cj,rjand j> 0 are the mean, the scale or width and a shape parameters, respectively, and j> 0 of the ith class-conditional distribution, respectively. (z)=∞ 0z−1ezdz, for z> 0. The scale parameter rjthat expresses the width of the distribution is related to the normal standard deviation by the equation: rj=j(1/j) (3/j)(2) Based on the GGD probability distribution, we define a novel RBF, by removing the constraints of a probability function, that can generalize to the Standard Gaussian RBF (SRBF) by adding a new parameter which can relax or contract the basis functions. In this way, the Generalized Radial Basis Function (GRBF) is defined using the following expression, for a K-dimensional input space: (x)=exp −x−cjj rj j(3) where Kis the number of inputs, xi= (x1i,. . .,xKi) is the vector of measurements, rjthe width of the GRBF, cj= (cj1,. . .,cjK) the centre and jthe exponent of the jth GRBF. Fig. 2 presents the radial unit activation for the GRBF for different values of . From Fig. 2, another observation can be made: the shape of is approximately rectangular for high values of . This implies that the GRBF kernel should be a good candidate for approximating constant functions at specific intervals. Moreover, as shown in Fig. 2, GRBF has a unique maximum at radial symmetry, and a local support property which complies with the fundamental properties of RBF used in Artificial Neural Networks (ANNs). Finally, due to the  parameter of the GRBF kernel, concavity or convexity is provoked around the point where the distance is the radii of the kernel, so the GRBF demonstrates appropriate kernel functions to quantify similarity in high dimensional spaces. When a normal distribution is 0 2 4 6 8 10 12 14 16 18 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 q f(q,d) d=1 d=2 d=20 d=100 d=200 (a) 0 2 4 6 8 10 12 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 q f(q,d = 2) f(q,d = 100) GRBF (r = 10.6, τ= 13) SRBF (r = 10.6, τ= 2) SRBF (r = 1.6, τ= 2) (b) Fig. 3. Probability density of a point from a normal distribution to be at distance q. assumed, the probability density function to find a point at distance qfrom the mean of the distribution is given by: f(q, d)=qd−1 2(d/2)−1×e−q2/2 (d/2)(=1),(4) where dis the input data dimension. For one dimension (d= 1), f(q, d) is maximum at the mean, but when the dimension grows, f(q,d) moves away from the mean (see Fig. 3a). Therefore, the probability of finding patterns near the mean, when the dimension is high, is almost zero. In conclusion, SRBF are no longer local in higher dimensions, and those models that have been seen as sums of local kernels do not behave as such in high dimensions. Fig. 3b shows that the GRBF provides better matching to the patterns in high dimensional spaces, because the parameter allows modification of the shape of the GRBF curvature. As can be observed in Fig. 3b when d= 100, the SRBF needs a high value of rto model these patterns. When the SRBF has a high value of r, it has a slightly pronounced curvature which causes the SRBF to assign very similar grades of membership to the patterns located far from the cluster centre. This effect can be observed in Fig. 3b: when d= 2, the SRBF shows its ability to fit distance distribution, assigning membership values Fig. 4. Hybrid Algorithm (HA) framework. in the interval [0, 1]; however when d= 100, the SRBF assigns membership values in the interval [0.27–0.57], while the GRBF assigns membership values in the interval [0–1]. In our opinion and based also on the experimental results (Section 4.4), this justifies our considering GRBF to be an appropriate kernel to quantify similarity in high dimensional spaces. 2.2. Base classifier: probabilistic Generalized Radial Basis Functions In a classification problem, measurements xi,i= 1, 2, . . .,K, of a single individual (or object) are taken, and the individuals are to be classified into one of the Jclasses based on these measurements. A training sample D={(xn,yn) ; n= 1, 2, . . .,N}is available, where xn= (x1n,. . .,xkn) is the random vector of measurements taking values in ⊂RK, and ynis the class level of the nth individual, where the common technique of representing class levels using a “1-of-J” encoding vector is adopted, y= (y(1),y(1),. . .,y(J)), and the Correctly Classified Rate or accuracy of the classifier is defined by C=(1/N)N n=1I(C(xn)=yn), where I(.) is the zero-one lost function. In order to tackle this classification problem, the outputs of the GRBF model have been interpreted from the point of view of probability through the use of the softmax activation function, which is given by: gl(x,␪l)=exp fl(x,␪l) J j=1exp fj(x,␪j), l =1,2, . . . , J (5) where Jis the number of classes in the problem, fj(x,␪l) is the output of the joutput neuron for pattern xand gl(x,␪l) is the probability a pattern xhas of belonging to class j. The model to estimate the function fl(x,␪l) is defined by the following equation: fl(x,␪l)=ˇl 0+ M  j=1 ˇl jexp −x−cjj rj j, l =1,2, . . . , J (6) where Mis the number of GRBFs or number of nodes in the hidden layer. Using the softmax activation function presented in Eq. (5), the class predicted by the NN corresponds to the node in the output layer with the greatest output value. The function used to evaluate a GRBF Neural Network (GRBFNN) is the function of cross-entropy error and it is given by the following expression: l(␪)= N  n=1 [− J  l=1 y(l) nfl(xn,␪l)+log J  l=1 expfl(xn,␪l)] (7) where ␪= (␪1,. . .,␪J). The error surface associated with the model is very convoluted with numerous local optima and the Hessian matrix of the error function l(␪) is, in general, indefinite. Moreover, the optimal number of basis functions in the model (i.e. the number of hidden nodes in the neural network) is unknown. Thus, we estimate the parameters ␪by means of a Hybrid Algorithm (HA). 2.3. Hybrid Algorithm The proposed Hybrid Algorithm (HA) is composed by two stages. In the first stage, an Evolutionary Algorithm is used as a global stochastic search algorithm which generates candidate RBFNNs. In the second stage, the iRprop+ algorithm performs a local optimization procedure to the best RBFNN individual of the last generation. Fig. 4 describes the procedure to estimate the parameters of GRBFNN. The basic framework of the EA is the following: the search begins with an initial population of RBFNNs and, in each iteration, the population is updated using a population-update algorithm which evolves both its structure and weights. The population is subject to operations of replication and mutation. The main characteristics of the algorithm are the following: 1. Representation of the individuals. The algorithm evolves architectures and connection weights simultaneously, each individual being a fully specified RBFNN. The neural networks are represented using an object-oriented approach and the algorithm deals directly with the RBFNN phenotype. Each connection is specified by a binary value indicating if the connection exists, and the real value representing its weights. 2. Error and fitness functions.We consider l(␪) (Eq. 7) as the error function of an individual gin the population. The fitness measure needed for evaluating the individuals is a strictly decreasing transformation of the error function l(␪) given by A(␪) = (1/1 + l(␪)), where 0 < A(␪)≤1. 3. Initialization of the population. The initial population is generated trying to obtain RBFNNs with the maximum possible fitness. First, 5000 random RBFNNs are generated. The centres of the radial units are firstly defined by the k-means algorithm for different values of k, where k∈[Mmin,Mmax], Mmin and Mmax being the minimum and maximum numbers of hidden nodes allowed for any RBFNN model in the HA. The widths of the RBFNNs are initialized to the geometric mean of the distance to the two nearest neighbourhood and the parameter to 2, since when = 2 the GRBF reduces to the standard Gaussian RBF (SRBF). A random value in the [−I,I] interval is assigned for the weights between the hidden layer and the output layer. The individuals obtained areevaluatedusingthefitnessfunctionandtheinitialpopulation is finally obtained by selecting the best 500 RBFNNs. 4. Parametric and structural mutations. Parametric mutation consists of a simulated annealing algorithm [22]. Structural mutation implies a modification in the structure of the RBFNNs and allows the exploration of different regions in the search space, helping to keep the diversity of the population. There are four different structural mutations: hidden node addition, hidden node deletion, connection addition and connection deletion. These four mutations are applied sequentially to each network. More information about the genetic operators proposed can be seen in [23]. It is important to note the structural and parametric mutations of : •Structural mutation: If the structural mutator adds a new node in the RBFNN, the parameter is assigned to 2, since when = 2 the GRBF reproduces to the Gaussian RBF. •Parametric mutation: The parameter is updated by adding a εvalue, where ε∈[−0.25, 0.25], since the modification of the GRBF is very sensitive to the variation (Fig. 2). 5. iRprop+Local Optimizer. The local optimization algorithm used in our paper is the iRprop+[24] optimization method. The iRprop+ is believed to be a fast and robust learning algorithm. This algorithm applies a backtracking strategy (i.e. it decides whether to take a step back along a weight direction or not by means of a heuristic). In the methodology proposed, we run the EA and then apply the local optimization algorithm to the best solution obtained by the EA in the last generation. Further details about the adaptation of the iRprop+ local improvement procedure to the softmax activation function can be seen in Appendix A of the paper. Since the accuracy of the GRBFNN model was evaluated with bioinformatic datasets and a typical microarray dataset may contain thousands of genes, applying the FS techniques are a prerequisite for building the GRBFNN model in this context. The methodology proposed is shown in detail in Fig. 5. 3. Hybrid-generation feature selection 3.1. Introduction The limitations of both the approaches introduced in Section 1(FR and FSS) in high-dimensional spaces, clearly suggest the need for a hybrid model. The three methods used in this work can be labelled as this kind of framework, Hybrid-Generation Feature Selection. In feature subset selection, it is a fact that two types of features are generally perceived as being unnecessary: features that are irrelevant to the target concept, and features that are redundant due to other features. The purpose of a feature subset algorithm is to identify relevant features according to a definition of relevance. However, the notion of relevance in machine learning has not yet been rigorously defined by common agreement [25]. On the other hand, notions of feature redundancy are normally in terms of feature correlation. It is widely accepted that two features are redundant to each other if their values are completely correlated. There are two widely used types of measures for the correlation between two variables: linear and non-linear. In the linear, the Pearson correlation coefficient is used, and in the case of nonlinear, many measures are based on the concept of entropy, or the measurement of the uncertainty of a random variable. Symmetrical uncertainty (SU) [26] is frequently used, defined as SU(x,y)=2IG(x|y) H(x)+H(y) where H(x)= −K iP(xi)log2(P(xi)) is the entropy of a variable x and IG(x|y) = H(x)−H(x|y) is the information gain from xprovided by y. Both of them are between pairs of variables. However, it may not be as straightforward in determining feature redundancy when one is correlated with a set of features. CFS [27] is one of well-known techniques to rank the relevance of features by measuring correlation between features and classes and between features and other features. The heart of the CFS (Correlation-based Feature Selection) algorithm contains a heuristic for evaluating the worth or merit of a subset of features. This heuristic takes into account the usefulness of individual features for predicting the class label, along with the level of intercorrelation among them. The hypothesis on which the heuristic is based is: Good feature subsets contain features that are highly correlated with the class, yet uncorrelated with one another. MeritS=k×rcf k+k×(k−1) ×rff where MeritSis the heuristic of a feature subset Scontaining kfeatures, rcf the average feature-class correlation, and rff the average feature–feature intercorrelation. For discrete class problems, CFS first discretizes numeric features and then uses symmetrical uncertainty to estimate the degree of association between the discretized features. Due to the high computational cost in a high dimensional domain, we discard the wrapper approach, and the three methods use correlation concepts as relevance and redundancy criteria. 3.2. Fast Correlation-Based Filter (FCBF) Aiming to achieve high efficiency, FCBF calculates SU–correlation between any Fifeature and class Cgenerating a list in descending order, and heuristically decides a Fifeature to be relevant if it is highly correlated with class C, i.e., if SUi,c>ı, where ıis a relevance threshold which can be determined by users. The relevant features selected are then subject to redundancy analysis. Similarly, FCBF evaluates the SU-correlation between individual features for redundancy analysis based on an approximate Markov blanket concept. For two relevant Fiand Fj(i/=j) features, Fjcan be eliminated if SUi,c≥SUj,cand SUi,j≥SUj,c. The iteration starts from the first element in the ranking and continues as follows. For all the remaining features, if Fihappens to form an approximate Markov blanket for Fj,Fjwill be removed from the list. After one round of filtering features based on Fi, the algorithm will take the remaining feature right next to Fias the new reference to repeat the filtering process. The algorithm stops when no more features can be eliminated. Fig. 6 describes the FCBF method. 3.3. Best Incremental Ranked Subset (BIRS) BIRS considers that relevance and redundancy concepts are included in the following “incremental usefulness” definition: given a sample of data, an evaluation measure L, a feature space Fig. 5. Flow diagram of the GRBF method. Fand a feature subset S(S⊆F), the feature Fiis incrementally useful to Lwith respect to Sif the evaluation of the hypothesis that Lproduces using the group of features {Fi}∪Sis better than the evaluation achieved using just the subset of features S. i.e., if Fiis not incrementally useful to Lwith respect to S, then the evaluation value given the subset Sis equal or better than the subset evaluation result known {Fi}∪S. It suggests that Figives no information beyond what is already in S, therefore, Ficould be removed safely, or in this case, Fiwould not be added to S. However, since the computational complexity to determine all possible interactions among features is very high (mainly in high-dimensional domains), BIRS considers using a guided search in preference to an ordered list of attributes. BIRS deals with incremental ranked usefulness in order to devise an approach to explicitly identify relevant features and not take into account redundant features. The idea is to choose the Fifeature from a ranked list one by one in the following way: firstly, the features are ranked according to some evaluation measure (SUi,c); secondly, BIRS deals with the list of features once, crossing the ranking from the beginning to the last ranked feature. The evaluation results using CFS with the first feature in the list are obtained and it is marked as selected. The result is obtained again with the first and second features. The second will be marked as selected depending on whether the evaluation obtained is statistically significantly better. The process is repeated until the last feature on the ranked list is reached. Finally, the algorithm returns the best subset found, and Fig. 6. FCBF algorithm. it will not contain irrelevant or redundant features. Fig. 7 describes the BIRS method. 3.4. Best Agglomerative Ranked Subset (BARS) BARS is called agglomerative due to the way it constructs the final subset of selected features. The method begins by generating a ranking. Then, pairs of features are obtained with the ranking’s first features, in combination with each one of the remaining features on the list. The pairs of features are ranked according to the value of the evaluation, and the process is repeated, that is, the subsets made up by the first sets on the new list are compared with the rest of the sets. The process continues until the final subset obtained. Step one generates a feature ranking ranging from best to worst according to a correlation measure (CFS). Next, a list of solutions is generated, in such a way that a solution for each individual feature is created and the same ranking order is maintained. The agglomerative search consists of making a subset of relevant features by joining subsets with a lower number of features. With every iteration a new list of solutions from the previous structure is generated. Each candidate set, made by joining two sets from the previous list of solutions, will become part of the next list of solutions if, when the subset evaluator (CFS) is applied to it, it gives back a higher measure value than the one obtained with the best (or first) subset from the previous list of solutions. To prevent the algorithm from becoming prohibitively time consuming, new sets of features are generated by joining the first sets to the remaining previous list of solutions. That is, the first set on the list is joined to the second set, next the first set is joined to the third set, and so on until the end of the list. Next, the second set of the list is joined to the third set, the second set and the fourth set, and so on until the last set on the list. This process of combining a set of features with the rest of the sets on the list is carried out with the best kfeature sets from the previous list of solutions. The process ends when only one feature subset is left, or when combining the subsets no longer causes an improvement. At the end, the algorithm returns the best positioned feature subset of all the subsets evaluated. Fig. 8 describes a sequence of numbered steps of the BARS method. 4. Experiments This section presents the experimental results and analysis of GRBF models on 6 public microarray datasets with high dimensionality/small sample size. At the beginning, the datasets and several machine learning algorithms used in this analysis are briefly described. Subsequently, experimental results are given and discussed with respect to various aspects. 4.1. Microarray data To validate the effectiveness of our method, a series of experiments were performed on 6 publicly available gene microarray datasets (Table 1). These datasets were taken from bioinformatic and biomedical domains. They are often used to validate the performance of the classifier and gene selector. Due to high dimensionality and small sample size, gene selection is an essential prerequisite for further data analysis. Brief descriptions of them are given in continuation. Breast consists of 97 samples collected from breast cancer patients. 46 of them are from patients labelled as relapse, the rest of the 51 samples are from patients who remain healthy from the disease and are regarded as non-relapse. Each sample is described by 24,481 genes. CNS (Central Nervous System) is derived from patient samples in embryonal tumours of the central nervous system. The total number of genes to be tested is 7129 and the number of samples is 60. There are two types of samples in the dataset, where 21 are survivors (who survive the treatment) and 39 are failures (who succumbed to the disease). Colon uses Affymetrix oligonucleotide arrays to monitor expression levels of over 6500 human genes from 40 tumour and 22 normal colon tissue samples. The 2000 genes with the highest minimal intensity across the 62 tissues were used in this analysis. Leukaemia refers to the primary disorders of bone marrow. This dataset contains 72 samples with malignant neoplasms of haematopoietic stem cells, of which 47 are acute lymphobastic leukaemia (ALL) and 25 acute myeloid leukaemia (AML). The total number of genes to be tested is 7129. Lung has 12,600 genes in 203 samples. The 203 samples consist of 139 lung adenocarcinomas (AD), 21 squamous (SQ) cell carcinoma cases, 20 pulmonary carcinoid (COID) tumours and 6 small cell lung cancer cases (SCLC), as well as 17 normal lung (NL) samples. GCM contains 190 samples. These samples are divided into 14 varieties of tumour. The expression levels of 16,063 genes are reported. Fig. 7. BIRS algorithm. Table 1 Characteristics of the six datasets used for the experiments: feature selection type (FS), number of instances (Size), number of Real (R), Binary (B) and Nominal (N) input variables, total number of inputs (# In), number of classes (# Out), per-class distribution of the instances (Distribution), minimum and maximum number of hidden nodes used for each dataset ([Mmin,Mmax]) and the number of generations (# Gen). Dataset Source FS Size R B N # In # Out Distribution [Mmin,Mmax]Gen Breast Van’t Veer et al. [40] BARS 97 183 – – 183 2 (46,51) [1, 3] 100 BIRS 261 – – 261 FCBF 493 – – 493 CNS Pomeroy et al. [41] BARS 60 187 – – 187 2 (21,39) [1, 3] 10 BIRS 206 – – 206 FCBF 170 – – 170 Colon Alon et al. [42] BARS 62 58 – – 58 2 (40,22) [1, 3] 10 BIRS 93 – – 93 FCBF 59 – – 59 Leukaemia Golub et al. [43] BARS 72 225 – – 225 2 (42,25) [1, 3] 50 BIRS 240 – – 240 FCBF 203 – – 203 Lung Bhattacharjee et al. [44] BARS 203 237 – – 237 5 (139,17,6,21,20) [5, 8] 100 BIRS 263 – – 263 FCBF 250 – – 250 GCM Ramaswamy et al. [45] BARS 253 311 – – 311 14 (11,10,11,11,22, [25, 28] 400 BIRS 288 – – 288 11,10,10,30,11, FCBF 264 – – 264 11,11,11,20) In these 6 microarray datasets, all gene expression values are numeric. For convenience sake, they were standardized before our experiments, that is, for each gene represented, the mean and standard deviation were zero and one, respectively, after the standardized operation had been performed. 4.2. Alternative statistical and artificial intelligence methods used for comparison purposes Different state-of-the-art statistical and artificial intelligence algorithms have been implemented for comparison purposes. Specifically, the results of the following algorithms have been compared to the GRBF method presented in this paper: 1. A Gaussian Radial Basis Function Network (RBFN) [28], deriving the centres and width of hidden units using k-means, and combining the outputs obtained from the hidden layer using logistic regression. 2. The MultiLogistic (MLogistic) algorithm. It is a method for building a multinomial logistic regression model with a ridge estimator to guard against overfitting by penalizing large coefficients [29]. 3. The SimpleLogistic (SLogistic) algorithm. It is based on applying LogitBoost algorithm with simple regression functions and determining the optimum number of iterations by a five fold cross-validation [30]. 4. The C4.5 classification tree inducer [31]. 5. The Naive Bayes standard learning algorithm (NaiveBayes) [28]. 6. The Logistic Model Tree (LMT) [30] classifier. 7. The IB1 classifier [32]. It uses a simple distance measure to find the training instance closest to the given test instance, and predicts the same class as this training instance. 8. The Support Vector Machine (SVM) classifier [33] with RBF kernels. These algorithms have been selected for comparison since they are some of the best performing algorithms in recent literature on classification problems. Many of these approaches have also been tested before in the classification problem on microarray gene expression. The detailed description and some previous results of these methods can be found in [34,28,30]. Fig. 8. BARS algorithm. Table 2 Comparison of the proposed method to other probabilistic methods: mean and standard deviation (SD) of the accuracy results (CG(%)) from 30 executions, mean accuracy (CG( %)) and mean ranking (R). Dataset FS Method (CG(%)) RBFN MLogistic SLogistic C4.5 NaiveBayes LMT IB1 SVM GRBF Result Result Result Result Result Result Result Result Mean ±SD Breast BARS 80.00 80.00 72.00 72.00 52.00 72.00 76.00 80.00 83.06 ±2.07 BIRS 72.00 76.00 64.00 64.00 76.00 64.00 80.00 84.00 82.26 ±5.32 FCBF 80.00 84.00 84.00 64.00 80.00 84.00 76.00 76.00 86.66 ±3.53 CNS BARS 80.00 73.33 66.66 66.66 66.66 66.66 73.33 66.67 76.88 ±3.38 BIRS 80.00 93.33 60.00 73.33 80.00 60.00 86.66 66.67 82.00 ±3.97 FCBF 86.66 100.00 80.00 60.00 86.66 80.00 93.33 66.67 96.44 ±5.46 Colon BARS 93.75 93.75 100.00 81.25 93.75 100.00 87.75 62.50 100.00 ±0.00 BIRS 81.25 68.75 67.50 87.50 81.25 87.50 81.25 62.50 85.83 ±4.32 FCBF 87.50 75.00 81.25 75.00 81.25 75.00 81.25 62.50 86.04 ±3.91 Leukaemia BARS 94.44 100.00 88.89 83.33 100.00 88.89 94.44 66.67 100.00 ±0.00 BIRS 94.44 100.00 94.44 83.33 100.00 94.44 100.00 66.67 100.00 ±0.00 FCBF 94.44 94.44 83.33 83.33 100.00 83.33 100.00 66.67 100.00 ±0.00 Lung BARS 96.07 96.07 98.03 94.11 96.07 98.03 98.03 98.03 99.64 ±0.98 BIRS 94.11 90.19 98.03 94.11 98.03 98.03 98.03 98.03 98.45 ±0.67 FCBF 94.11 94.11 98.03 74.50 94.11 98.03 94.11 94.11 97.50 ±0.98 GCM BARS 75.00 73.07 63.49 57.69 75.00 71.15 59.61 75.00 82.32 ±4.30 BIRS 76.92 78.84 75.00 57.69 78.84 75.00 75.00 76.92 79.87 ±3.45 FCBF 82.00 80.76 71.15 48.07 71.15 67.30 69.23 80.76 79.75 ±4.90 CG( %) 85.70 86.20 80.32 73.32 83.93 81.29 84.66 75.02 89.81 RCG4.41 4.22 5.69 7.69 4.66 5.58 4.69 6.05 1.97 The best result is in bold face and the second best result in italics. 4.3. Experimental design The evaluation of the different models has been performed using two different measures: Correctly Classified Rate (CCR) or accuracy and Root Mean Square Error (RMSE). CCR represents threshold metrics and RMSE a rank metric. RMSE is a metric corresponding to the expected value of the squared error loss or quadratic loss. RMSE is a frequently used measurement of the differences between values predicted by a model or an estimator, and the values actually observed in what is being modelled or estimated. All the parameters used in the HA (Section 2.3) except the maximum and minimum number of RBFs in the hidden layer ([Mmin, Mmax]) and the number of generations (# Gen) have the same values in all problems analysed below (Table 1). The connections between hidden and output layer are initialized in the [−5, 5] interval (i.e. [−I,I] = [−5, 5]). The size of the population is N= 500. For the structural mutation, the number of nodes that can be added or removed is within the [1,2] interval, and the number of connections to add or delete in the hidden and the output layers during structural mutations is within the [1,7] interval. For the selection of the SVM hyperparameters (regularization parameter, C, and width of the Gaussian functions, ), a grid search algorithm has been applied with a ten-fold cross-validation, using the following ranges: C∈{2−5, 2−3,. . ., 215}and ∈{2−15, 2−13, ..., 23}. The experimental design was conducted using a holdout cross validation procedure with 3n/4 instances for the training dataset and n/4 instances for the generalization dataset. In order to evaluate the stability of the methods, the evolutionary algorithm is run 30 times. The HA and the model proposed was implemented in JAVA. We also used “libsvm” [35] to obtain the results of the SVM method, and WEKA1to obtain the results of the remaining methods. 1http://www.cs.waikato.ac.nz/ml/weka/ 4.4. Comparison of the GRBF model with other classifiers In this subsection, the GRBF model is compared to other base line classifiers. The purpose of this section is to show the improvement in accuracy in the classification problem of the microarray gene expression. Tables 2 and 3 show the mean and the standard deviation of the correct classification rate (CG) and the Root Mean Square Error (RMSEG) in the generalization set for each dataset and the RBFN, MLogistic, SLogistic, C4.5, NaiveBayes, LMT, IB1, SVM and GRBF methods. Based on the mean CGand RMSEG, the ranking of each method in each dataset (R= 1 for the best performing method and R= 9 for the worst one) is obtained and the mean accuracy and RMSE (CGand RMSEG) and the mean ranking (RCGand RRMSEG) are also included in Tables 2 and 3. From the analysis of the results, it can be concluded, from a purely descriptive point of view, that the GRBF method obtained the best results for ten datasets in CGand for fourteen datasets in RMSEG. Furthermore, the GRBF method yield the best mean (CG=89.81 %) and ranking (RCG=1.97) in CG, and, taking RMSEG into account, the GRBF got the best performance in both measures (RMSEG=0.18, RRMSEG=1.75). To determine the statistical significance of the rank differences observed for each method in the different datasets, we have carried out a non-parametric Friedman test [36] with the ranking of CGand RMSEGof the best models as the test variables (since a previous evaluation of the CGand RMSEGvalues results in rejecting the normalityand the equality of the variances’hypothesis).The test shows that the effect of the method used for classification is statistically significant at a significance level of 5%, as the confidence interval is C0= (0, F0.05 = 2.00) and the F-distribution statistical values are F*= 8.19 /∈C0for CGand F*= 9.19 /∈C0for RMSEG. Consequently, we reject the null-hypothesis stating that all algorithms perform equally in mean ranking. Based on this rejection, the Nemenyi post hoc test is used to compare all the classifiers to each other. This test considers that the performance of any two classifiers is deemed significantly