2.1. Feature Selection and Classification through GMDH Algorithm

GMDH-type neural network algorithm is a heuristic self-organization method that investigates the relations among the variables. The algorithm defines its structure itself. Ivakhnenko [11] presented the following polynomial—known as the Ivakhnenko polynomial—to construct a high order polynomial:

where m is the number of variables to be regressed in each neuron and a,b,c,… are weights of variables in the polynomial. Here, y is a response variable, xi and xj are the exploratory variables. In this study, only the main effects are included in the model as presented below.

The GMDH algorithm, in general, investigates all pairwise combinations of p exploratory variables. Therefore, m is specified as 2 in Equation (2). For this algorithm, there exist three weights to be estimated in each neuron. The weights are estimated via least square estimation.

In model building and evaluation process, the data are divided into three sets: train (60%), validation (20%), and test (20%) sets. Train set is included in model building. Validation set is used for neuron selection. Test set is utilized to estimate the performance of the methods on unseen data.

The GMDH algorithm can be depicted as follows:

1. Each pairwise combination goes into one neuron.

2. Weights are estimated with train set in each neuron at layer k.

3. The predicted probabilities of train set are estimated in each neuron at layer k.

4. The predicted probabilities of validation set are estimated in each neuron at layer k.

5. The external criterion EC (i.e., mean square error) is calculated using validation set in each neuron at layer k.

6. Selection pressure (α) and the maximum number of neurons to be selected need to be specified.

7. The neurons of which external criteria are smaller than α⋅minEC+1−α⋅maxEC/2 are selected. If the number of selected neurons is larger than the specified maximum number of neurons, the neurons—as many as the specified maximum number of neurons—having smaller external criterion compared to the rest of them are selected.

8. The predicted probabilities of train set obtained from selected neurons become the inputs for the next layer.

9. This process (i) to (viii) continues until the stopping rule is realized.

10. There are three stopping rules to conclude the algorithm. The first one is an increase in minimum external criterion at consecutive layers. Second, the algorithm stops when the specified maximum number of layers is reached. The third one is that the algorithm stops if only one neuron in a layer is selected.

11. At the last layer, only one neuron having minimum EC is selected.

GMDH algorithm is a system of layers where the neurons are present. The number of neurons in each layer is determined by the number of inputs. For example, providing that the number of inputs going into a layer is equal to p, the number of neurons in that layer becomes h=p2, since all pairwise combinations of inputs are considered. This does not mean that all layers include h neurons. For instance, the number of inputs in the input layer defines just the number of neurons in first layer. The number of neurons selected in the first layer determines the number of neurons in second layer. The algorithm organizes the architecture itself. Sample architecture of GMDH algorithm is placed in Figure 1 when there exist three layers and four inputs.

Figure 1

Architecture of group method of data handling (GMDH) algorithm.

In the GMDH architecture shown in Figure 1, there exist four inputs (X1, X2, X3, X4). From these input variables, three of them (X1, X2, X4) are dominating the system. X3 does not have an impact on classification. GMDH algorithm selects these important features having an effect on classification.

2.2. Diverse Classifiers Ensemble Based on GMDH Algorithm

dce-GMDH algorithm is the GMDH algorithm which assemble the well-known classifiers—support vector machines, random forest, naive Bayes, elastic net logistic regression, artificial neural network. These classifiers are available in e1071 [5], randomForest [7], e1071 [5], glmnet [25], nnet [6] packages, respectively. Specifically, these classifiers are available in svm (e1071), randomForest (randomForest), naiveBayes (e1071), cv.glmnet (glmnet), nnet (nnet) functions, respectively. Unlike GMDH algorithm, dce-GMDH algorithm includes base layer. The classifiers are placed at base layer. Predicted probabilities are obtained using all inputs through these classifiers. The predicted probabilities obtained from these classifiers continue their way as inputs of first layer without applying any neuron selection process. The rest of the algorithm is same as GMDH algorithm. The sample architecture of dce-GMDH algorithm is stated in Figure 2.

Figure 2

Architecture of diverse classifiers ensemble based on group method of data handling (dce-GMDH) algorithm.

The dce-GMDH algorithm is a system of layers where the neurons exist. The number of neurons in a base layer is five since the five classifiers are included. The number of neurons in other layers is defined by the number of inputs. The algorithm assembles the most appropriate classifiers by organizing itself.

In the dce-GMDH architecture shown in Figure 2, there exist four inputs (X1, X2, X3, X4). These four inputs enter each neuron at base layer. There exists a different classifier in each neuron at base layer. Predicted probabilities are obtained by utilizing four inputs through the classifiers. These predicted probabilities obtained from these classifiers continue to first layer without applying any neuron selection process. Since five inputs will enter in the first layer, the number of neurons in that layer becomes 52=10. According to external criterion, four neurons are selected and six neurons are eliminated from the network. Since four neurons are selected in the first layer, the number of neurons in the second layer becomes 42=6. This process continues until one of the stopping rules is realized. Also, the algorithm returns which classifiers are assembled.

3.1. Table of the Descriptive Statistics: Table()

Table() produces a table for simple descriptive statistics for a binary response. It returns frequency (percentage) for the variables with class of factor/ordered. Also, this function returns mean ± standard deviation (median, minimum–maximum) or mean ± standard deviation (median, quartile1 - quartile3) for the variables with class of numeric/integer. The option argument is used to return minimum - maximum or quartile1–quartile3 values. When this argument is set to “min-max”, this function returns mean ± standard deviation (median, minimum–maximum). When this argument is set to “Q1–Q3,” this function returns mean ± standard deviation (median, quartile1–quartile3). The percentages can be specified with the percentages argument as row, column, or total percentages. The ndigits argument is a vector of two numbers utilized to specify the number of digits. The first one is used to specify the number of digits for numeric/integer variables. The second one specifies the number of digits for percentages of factor/ordered variables. Default is set to ndigits = c(2,1). There exists output argument to return the output in a specified format (R, LaTeX, HTML). In this example, we use “LaTeX” output. # obtain a table for simple descriptive  statistics for a binary response R> Table (x, y, option = "min-max",  percentages = "column", ndigits =  c(2,1), output = "LaTeX")

Some portion of the output is given below and its table version is presented in Table 1. \begin{table}[ht] \centering \begin{tabular}{rrrrr}   \hline   &  & benign &   & malignant \\    \hline Observations &   & 444 &   & 239 \\    Cl.thickness &   &  &   &  \\    1 &   & 136 (30.6\%) &   &  3 ( 1.3\%) \\    2 &   &  46 (10.4\%) &   &  4 ( 1.7\%) \\    3 &   &  92 (20.7\%) &   & 12 ( 5.0\%) \\    4 &   &  67 (15.1\%) &   & 12 ( 5.0\%) \\    5 &   &  83 (18.7\%) &   & 45 (18.8\%) \\    6 &   &  15 ( 3.4\%) &   & 18 ( 7.5\%) \\    7 &   &   1 ( 0.2\%) &   & 22 ( 9.2\%) \\    8 &   &   4 ( 0.9\%) &   & 40 (16.7\%) \\    9 &   &   0 ( 0.0\%) &   & 14 ( 5.9\%) \\    10 &   &   0 ( 0.0\%) &   & 69 (28.9\%) \\    Cell.size &   &  &   &  \\    1 &   & 369 (83.1\%) &   &  4 ( 1.7\%) \\    2 &   &  37 ( 8.3\%) &   &  8 ( 3.3\%) \\    3 &   &  27 ( 6.1\%) &   & 25 (10.5\%) \\    4 &   &   8 ( 1.8\%) &   & 30 (12.6\%) \\    5 &   &   0 ( 0.0\%) &   & 30 (12.6\%) \\    6 &   &   0 ( 0.0\%) &   & 25 (10.5\%) \\    7 &   &   1 ( 0.2\%) &   & 18 ( 7.5\%) \\    8 &   &   1 ( 0.2\%) &   & 27 (11.3\%) \\    9 &   &   1 ( 0.2\%) &   &  5 ( 2.1\%) \\    10 &   &   0 ( 0.0\%) &   & 67 (28.0\%) \\     \hline \end{tabular} \end{table}

3.2. Feature Selection and Classification through GMDH Algorithm: GMDH()

In this section, we demonstrate GMDH() function for feature selection and classification. It constructs GMDH algorithm, returns summary statistics of GMDH architecture and important variables. First, we randomly divide data into train, validation and test sets, and then call the GMDH() function. The first and second arguments in this function are a matrix of the exploratory variables and a factor in training set, respectively. The third and fourth arguments are a matrix of the exploratory variables and a factor in validation set, respectively. The alpha argument is the selection pressure. The maxlayers argument is the maximum number of layers specified. The maxneurons argument is the maximum number of neurons allowed in the second and the later layers. The exCriterion argument is the external criterion to be used for neuron selection. The verbose argument is utilized to print the output in R console. # change the class of x to a matrix R> x <- data.matrix(x) # the seed number is fixed to 12345 for  reproducibility R> seed <- 12345 # the number of observations R> nobs <- length(y) R> set.seed(seed) # to split train, validation and test  sets # to shuffle data R> indices <- sample(1:nobs) # the number of observations in each set R> ntrain <- round(nobs*0.6,0) R> nvalid <- round(nobs*0.2,0) R> ntest <- nobs-(ntrain+nvalid) # obtain the indices of sets R> train.indices <- sort(indices[1:ntrain]) R> valid.indices <- sort(indices[(ntrain+1) :(ntrain+nvalid)]) R> test.indices <- sort(indices[(ntrain +nvalid+1):nobs]) # obtain train, validation and test sets R> x.train <- x[train.indices,] R> y.train <- y[train.indices] R> x.valid <- x[valid.indices,] R> y.valid <- y[valid.indices] R> x.test <- x[test.indices,] R> y.test <- y[test.indices] R> set.seed(seed) # construct model via GMDH algorithm R> model <- GMDH(x.train, y.train,  x.valid, y.valid, alpha = 0.6,  maxlayers = 10, maxneurons = 15,  exCriterion = "MSE", verbose = TRUE)  Structure :   Layer   Neurons  S. neurons    Min MSE      1        36          15    0.06316      2       105          15    0.05310      3       105          15    0.05188      4       105          15    0.05161      5       105          15    0.05127      6       105          15    0.05110      7       105          15    0.05098      8       105          11    0.05096      9        55          15    0.05096     10       105           1    0.05095  External criterion : Mean Square Error   Feature selection  : 8 out of 9   variables are selected.                       Cl.thickness        Cell.size    Marg.adhesion     Epith.c.size      Bare.nuclei      Bl.cromatin  Normal.nucleoli          Mitoses

Here, the structure includes layer, neurons, s. neurons, and min MSE in the output above. The layer shows the number of layer. The neurons represent the number of neurons in corresponding layer. The s. neurons mean the number of selected neurons. The min MSE respresents the minimum external criterion which is calculated for the neuron gives the minimum external criterion on validation set in the corresponding layer. There exist two options for the external criterion namely, mean square error and mean absolute error.

Eight variables—clump thickness, uniformity of cell size, marginal adhesion, single epithelial cell size, bare nuclei, bland chromatin, normal nucleoli, mitoses—are selected by the algorithm. Minimum external criterion can be plotted across layers (presented in Figure 3) by the following code: R> plot(model)

Predictions for test set can be made after model building process is completed. Test set has 136 observations, but only 10 of them are reported to save space. R> predict(model, x.test, type = "class")   [1] benign benign benign benign benign    [6] benign malignant benign benign    [10] benign Levels: benign malignant R> predict(model, x.test, type =  "probability")             benign   malignant   [1,] 1.000000000 0.000000000   [2,] 0.643870382 0.356129618   [3,] 0.670641964 0.329358036   [4,] 0.974398179 0.025601821   [5,] 0.920988111 0.079011889   [6,] 0.994693987 0.005306013   [7,] 0.436033878 0.563966122   [8,] 0.951034736 0.048965264   [9,] 1.000000000 0.000000000  [10,] 0.994693987 0.005306013

Figure 3

Minimum external criterion across layers (GMDH algorithm).

The GMDH algorithm predicts that the probability of benign for the first and second persons are 100% and 64.4%, respectively. Since the predicted probability of benign is greater than the predicted probability of malignant, these persons are classified as benign.

3.3. Confusion Matrix and Related Statistics: confMat()

The confMat() function produces a confusion matrix for a binary response. It also returns some related statistics. These statistics are accuracy, no information rate, unweighted Kappa statistic, Matthews correlation coefficient, sensitivity, specificity, positive predictive value, negative predictive value, prevalence, balanced accuracy, youden index, detection rate, detection prevalence, precision, recall, and F1 measure. The formulation of these statistics are not stated in this paper, but presented in the manual of GMDH2 package. The positive argument is an optional character string used to specify the positive factor level. The verbose argument is utilized to print the output in R console. # obtain predicted classes for test set R> y.test_pred <- predict(model, x.test,  type = "class") # obtain confusion matrix and some  statistics for test set R> confMat(y.test_pred, y.test, positive  = "malignant") Confusion Matrix and Statistics             reference data        malignant benign   malignant        51      1   benign            5     79      Accuracy             :   0.9559      No Information Rate  :   0.5882      Kappa                :   0.9079      Matthews Corr Coef   :   0.9097      Sensitivity          :   0.9107      Specificity          :   0.9875      Positive Pred Value  :   0.9808      Negative Pred Value  :   0.9405      Prevalence           :   0.4118      Balanced Accuracy    :   0.9491      Youden Index         :   0.8982      Detection Rate       :   0.375      Detection Prevalence :   0.3824      Precision            :   0.9808      Recall               :   0.9107      F1                   :   0.9444       Positive Class       :   malignant 

Accuracy of GMDH algorithm is estimated to be 0.9559. This algorithm classifies 95.59% of persons in a correct class. Also, sensitivity and specificity are calculated as 0.9107 and 0.9875. The algorithm classifies 91.07% of the persons having breast cancer, 98.75% of the persons not having breast cancer.

3.4. Scatter Plots with Classification Labels: cplot2d() and cplot3d()

The cplot2d() and cplot3d() functions provide interactive 2-dimensional (Figure 4a) and 3-dimensional (Figure 4b) scatter plots with classification labels, respectively. These functions originally use the plot_ly function from plotly [28] package. The first two arguments of cplot2d() are the exploratory variables stated in the x and y axes of Figure 4a. The first three arguments of cplot3d() are the exploratory variables placed in the x, y, and z axes of Figure 4b. The ypred and yobs arguments are predicted and observed classes. The colors and symbols arguments are used to specify the colors and symbols of true/false classification labels, respectively. The size of symbols can be changed with the size argument. The names of axes can be changed with the arguments xlab, ylab, zlab, and title. # to produce 2D scatter plot with  classification labels for test set R> cplot2d(x.test[,1], x.test[,2],  y.test_pred, y.test, colors = c("red",  "black"), xlab = "clump thickness", ylab = "uniformity of cell size") # to produce 3D scatter plot with  classification labels for test set R> cplot3d(x.test[,1], x.test[,2],  x.test[,6], y.test_pred, y.test,  colors = c("red", "black"),  xlab = "clump thickness", ylab = "uniformity of cell size",  zlab = "bare nuclei")

Figure 4

Scatter plots with classification labels.

3.5. Diverse Classifiers Ensemble Based on GMDH Algorithm: dceGMDH()

In this part, we demonstrate dceGMDH() function for classification. It constructs dce-GMDH algorithm, returns summary statistics of dce-GMDH architecture and assembled classifiers. Like GMDH() function, the first and second arguments are a matrix of the exploratory variables and a factor in training set, respectively. The third and fourth arguments are a matrix of the exploratory variables and a factor in validation set, respectively. The alpha argument is the selection pressure. The maxlayers argument is the specified maximum number of layers. The maxneurons argument is the maximum number of neurons allowed in the second and later layers. The exCriterion argument is the external criterion to be utilized for neuron selection. The verbose argument is utilized to print the output in R console. Also, there are the arguments for options of classifiers. The svm_options argument is a list for options of svm. The randomForest_options argument is a list for options of randomForest. The naiveBayes_options argument is a list for options of naiveBayes. The cv.glmnet_options argument is a list for options of cv.glmnet (the elastic net mixing parameter is fixed to 0.5 as default). The nnet_options argument is a list for options of nnet. R> set.seed(seed) # construct model via dce-GMDH algorithm R> model <- dceGMDH(x.train, y.train,  x.valid, y.valid, alpha = 0.6, maxlayers  = 10, maxneurons = 15, exCriterion =  "MSE", verbose = TRUE)  Structure :   Layer   Neurons  S. neurons    Min MSE      0         5           5    0.04669      1        10           1    0.04641  External criterion   : Mean Square Error   Classifiers ensemble : 2 out of 5   classifiers are assembled.                    svm  cv.glmnet

In this example, two classifiers—support vector machine and elastic net logistic regression—are assembled by the algorithm. Minimum external criterion can be plotted across layers (presented in Figure 5) by the following line: R> plot(model)

Figure 5

Minimum external criterion across layers (dce-GMDH algorithm).

Predictions for test set can be made after model building process is completed. Test set has 136 observations; therefore, 10 of them are reported to save space. R> predict(model, x.test, type = "class")   [1] benign benign malignant benign    [5] benign benign malignant benign    [9] benign benign     Levels: benign malignant R> predict(model, x.test, type =  "probability")              benign    malignant   [1,] 0.9571287282 4.287127e-02   [2,] 0.8317147956 1.682852e-01   [3,] 0.3400820793 6.599179e-01   [4,] 1.0000000000 0.000000e+00   [5,] 0.9876416020 1.235840e-02   [6,] 1.0000000000 0.000000e+00   [7,] 0.2762650840 7.237349e-01   [8,] 1.0000000000 0.000000e+00   [9,] 1.0000000000 0.000000e+00  [10,] 1.0000000000 0.000000e+00

The dce-GMDH algorithm predicts that the probability of benign for the first and second persons are 95.7% and 83.2%, respectively. Since the predicted probability of benign is greater than the predicted probability of malignant, these persons are classified as benign.

Confusion matrix and related statistics are obtained through the following codes to investigate the performance measures for the test set: # obtain predicted classes for test set R> y.test_pred <- predict(model, x.test,  type = "class") # obtain confusion matrix and some  statistics for test set R> confMat(y.test_pred, y.test, positive  = "malignant") Confusion Matrix and Statistics             reference data        malignant benign   malignant        54      1   benign            2     79      Accuracy             :   0.9779      No Information Rate  :   0.5882      Kappa                :   0.9543      Matthews Corr Coef   :   0.9545      Sensitivity          :   0.9643      Specificity          :   0.9875      Positive Pred Value  :   0.9818      Negative Pred Value  :   0.9753      Prevalence           :   0.4118      Balanced Accuracy    :   0.9759      Youden Index         :   0.9518      Detection Rate       :   0.3971      Detection Prevalence :   0.4044      Precision            :   0.9818      Recall               :   0.9643      F1                   :   0.973       Positive Class       :   malignant 

Accuracy rate of dce-GMDH algorithm is estimated to be 0.9779. This algorithm classifies 97.79% of persons in a correct class.

All in all, using dce-GMDH algorithm increases the classification performance approximately 2% in accuracy compared to GMDH algorithm for this data set.

5.1. Case Study on Real Datasets

In this section, we will implement GMDH2 using several real datasets. The first data set is Wisconsin breast cancer data collected by Wolberg and Mangasarian [26]. The implementation of the package demonstrated on this data set is introduced in Section 3.

The second data set is ionosphere data downloaded from mlbench [27] R package. This data set was collected by a system in Goose Bay, Labrador. Free electrons were targeted in the ionosphere. This data set contains 351 radar returns with 34 attributes. Providing that the radar returns illustrate evidence of some type of structure in the ionosphere, they are called “good” radar returns. Otherwise, their signals pass through the ionosphere and they are called “bad” returns. There exist 225 and 126 observations in each group, respectively.

Finally, we utilized sonar data downloaded from mlbench [27] R package. Sonar signals are classified to be bounced off a metal cylinder or a roughly cylindrical rock. This data set includes 208 sonar signals (111 and 97 observations in each group, respectively) with 60 attributes. Each attribute is in the range 0.0 to 1.0. Each attribute represents the energy within a particular frequency band.

5.2. Implementation of the GMDH2 Package

We perform both GMDH and dce-GMDH classifiers on Wisconsin breast cancer data, Ionosphere data, and Sonar data. In the methods, selection pressure is set to 0.6. The maximum number of neurons is fixed to 15. The maximum number of layers is set to 10. MSE is used as an external criterion in all analyses.

The data set is split into three parts as train, validation and test sets including 60%, 20%, and 20% of all samples, respectively. In this part, the performances of the GMDH classifiers are reported based on the confusion matrices of true and predicted classes for test sets. The seed number is fixed to “12345” for the reproducibility of the results.

5.3. Performance Comparison of GMDH Classifiers

In this section, we discuss the performance of the GMDH classifiers. As it is mentioned in the earlier sections, several measures are considered for the model performances. We reported the results in Table 2 with accuracy, no information rate, Kappa, Matthews correlation coefficient, sensitivity, specificity, positive predictive value, negative predictive value, prevalence, balanced accuracy, youden index, detection rate, detection prevalence, precision, recall, and F1 measure.

According to the measures of accuracy, kappa, Matthews correlation coefficient, positive predictive value (precision), balanced accuracy, youden index and F1, dce-GMDH algorithm is superior to GMDH algorithm in all data sets. With respect to sensitivity (recall) and negative predictive value, dce-GMDH algorithm performs better compared to GMDH algorithm on both breast cancer and sonar data sets, but vice versa is true on ionosphere data set. According to the measure of specificity, dce-GMDH classifier performs as well as GMDH classifier except for ionosphere data, but it performs better than GMDH classifier on ionosphere data.