An extensive experimental survey of regression methods

Rúa Jenaro de la Fuente, s/n –Campus Vida – Universidade de Santiago de Compostela -15782 Santiago de Compostela – citius.usc.es © 2018 Elsevier Ltd. This manuscript version is made available under the CC-BY-NC-ND 4.0 license An extensive experimental survey of regression methods M. Fernández-Delgado, M.S. Sirsat, E. Cernadas, S. Alawadi, S. Barro and M. Febrero-Bande Version: accepted article M. Fernández-Delgado, M.S. Sirsat, E. Cernadas, S. Alawadi, S. Barro and M. Febrero-Bande (2019) An extensive experimental survey of regression methods. Neural Networks, 111, 11 - 34. Doi: https://doi.org/10.1016/j.neunet.2018.12.010 How to cite: Copyright information: An extensive experimental survey of regression methods M. Fern´andez-Delgado1,∗, M.S. Sirsat1, E. Cernadas1, S. Alawadi1, S. Barro1, M. Febrero-Bande2 Abstract Regression is a very relevant problem in machine learning, with many different available approaches. The current work presents a comparison of a large collection composed by 77 popular regression models which belong to 19 families: linear and generalized linear models, generalized additive models, least squares, projection methods, LASSO and ridge regression, Bayesian models, Gaussian processes, quantile regression, nearest neighbors, regression trees and rules, random forests, bagging and boosting, neural networks, deep learning and support vector regression. These methods are evaluated using all the regression datasets of the UCI machine learning repository (83 datasets), with some exceptions due to technical reasons. The experimental work identifies several outstanding regression models: the M5 rule-based model with corrections based on nearest neighbors (cubist), the gradient ∗M. Fern´andez-Delgado Email address: [email protected] (M. Fern´andez-Delgado ) URL: https://citius.usc.es/equipo/persoal-adscrito/manuel-fernandez-delgado (M. Fern´andez-Delgado ) 1Centro Singular de Investigaci´on en Tecnolox´ıas da Informaci´on da USC (CiTIUS), University of Santiago de Compostela, Campus Vida, 15782, Santiago de Compostela, Spain. 2Dept. of Statistics, Mathematical Analysis and Optimization, University of Santiago de Compostela, Campus Vida, 15782, Santiago de Compostela, Spain. Preprint submitted to Neural Networks January 7, 2019 boosted machine (gbm), the boosting ensemble of regression trees (bstTree) and the M5 regression tree. Cubist achieves the best squared correlation (R2) in 15.7% of datasets being very near to it, with difference below 0.2 for 89.1% of datasets, and the median of these differences over the dataset collection is very low (0.0192), compared e.g. to the classical linear regression (0.150). However, cubist is slow and fails in several large datasets, while other similar regression models as M5 never fail and its difference to the best R2is below 0.2 for 92.8% of datasets. Other well-performing regresors are the committee of neural networks (avNNet), extremely randomized regression trees (extraTrees, which achieves the best R2in 33.7% of datasets), random forest (rf) and ε-support vector regression (svr), but they are slower and fail in several datasets. The fastest regression model is least angle regression lars, which is 70 and 2,115 times faster than M5 and cubist, respectively. The model which requires least memory is non-negative least squares (nnls), about 2 GB, similarly to cubist, while M5 requires about 8 GB. For 97.6% of datasets there is a regression model among the 10 bests which is very near (difference below 0.1) to the best R2, which increases to 100% allowing differences of 0.2. Therefore, provided that our dataset and model collection are representative enough, the main conclusion of this study is that, for a new regression problem, some model in our top-10 should achieve R2near to the best attainable for that problem. Keywords: Regression, UCI machine learning repository, cubist, M5, gradient boosted machine, extremely randomized regression tree, support vector regression penalized linear regression. 2 1. Introduction1 The objective of this paper is to provide a “road map” for researchers2 who want to solve regression problems and need to know how well work3 the currently available regression methods. In machine learning, regression4 methods are designed to predict continuous numeric outputs where an order5 relation is defined. Regression has been widely studied from the statistics6 field, which provides different approaches to this problem: linear and gen-7 eralized linear regression, least and partial least squares regression (LS and8 PLS), least absolute shrinkage and selection operator (LASSO) and ridge9 regression, multivariate adaptive regression splines (MARS), least angle re-10 gression (LARS), among others. Furthermore, several methods arising from11 the field of machine learning were designed to be universal function aproxi-12 mators, so they can be applied both for classification and regression: neural13 networks, support vector machines, regression trees and rules, bagging and14 boosting ensembles, random forests and others. The current work develops15 an empirical quantitative comparison of a very large collection of regression16 techniques which is intended to provide the reader: 1) a list of the currently17 available regression models, grouped by families of related methods; 2) a18 brief description and list of references about each approach, alongside with19 technical details about its execution such as software implementation, list of20 tunable hyperparameters and recommended values; 3) a ranking of the avail-21 able models according to its performance and speed, identifying the best22 performing approach and the performance level which can be expected for it;23 and 4) the code to run all the regression models considered in this study for24 3 any regression problem3. In this comparison, we use the whole collection of25 regression datasets provided of the UCI machine learning repository (except-26 ing some datasets excluded by technical reasons), which a large collection of27 regression problems, and it should allow to develop a realistic and significant28 evaluation of the regression methods. As we explained in a previous paper29 comparing classifiers [1], provided that the size of the model collection used30 in the current comparison is large enough, we can assume that the best per-31 formance, measured in terms of squared correlation (R2), achieved by some32 regression model for each dataset (denoted as R2 best) is the highest attainable33 performance for that dataset. For a model which achieves a given R2in that34 dataset, the difference ∆ = R2 best −R2, averaged over the dataset collection,35 can be used as an estimation of the expected ∆ for that model and a new36 dataset D, not included in the collection. For the best model Xon the cur-37 rent comparison, it is expected that ∆ &0, i.e., the R2achieved by Xshould38 not be too far from Rbest in average over the data collection. Thus, although39 by the No-Free-Lunch theorem [2] we can not guarantee that Xwill be the40 best model for D, we can expect that Xwill achieve R2> R2 best −∆, so that41 Xwill not be very far from R2 best for dataset D. Consequently, the current42 paper may be useful for researchers who want to know how far a given model43 (e.g. the best model X) will be from the best available performance (which44 is, of course, unknown) for a new dataset. On the other hand, in general the45 best models in the current comparison achieve the best, or very near to the46 best, performances for most datasets in the collection. Therefore, although47 3https://nextcloud.citius.usc.es/index.php/s/Yb8LZQQFrgckjFk (visited December 14, 2018). 4 Xwill not be the best regression model for a new dataset D, we can expect48 that some of the best models in our comparison will achieve the best R2.49 Thus, the current comparison may be also useful to provide to the reader a50 reduced list (e.g., the 10 best performing models of the collection) which is51 expected to include the one which provides the highest available performance52 for a new dataset D.53 The section 2 describes the materials and methods used for this compar-54 ison, which include the list of datasets and regression methods, grouped by55 families. The description of regression models and issues related to their ex-56 ecution (software implementation, number of tunable hyperparameters and57 their values) are included in Appendix B. The section 3 reports the results58 of the experimental work and discusses them globally, by families of regres-59 sion models and by datasets, best model for each dataset, elapsed time and60 memory. Finally, the section 4 compiles the conclusions of this study.61 2. Materials and methods62 This section describes the scope of the current work, defined by the collec-63 tion of datasets used in this comparison (subsection 2.1) and by the regression64 methods that will be compared (subsection 2.2).65 Original UCI name Datasets #patterns #inputs 3D Road network 3Droad 434,874 4/3 Airfoil self-noise airfoil 1,503 5 Air quality air-quality-CO, air-quality-NMHC 1,230 8 air-quality-NO2, air-quality-NOx air-quality-O3 Appliances energy prediction appliances-energy 19,735 28/26 Auto MPG auto-MPG 398 8/23 Automobile automobile 205 26/66 Continued on next page. 5 Table 2 – Continued from previous page. Original UCI name Datasets #patterns #inputs Beijing PM2.5 beijing-pm25 41,758 12 Bike sharing bike-day 731 13/30 bike-hour 17,379 14/42 Blog feedback blog-feedback 60,021 280/13 Buzz in social media buzz-twitter 583,250 77 Combined cycle power plant combined-cycle 9,568 4 Communities & crime com-crime 1,994 122 Communities & crime unnormalized com-crime-unnorm 2,215 124/126 Computer hardware com-hd 209 7 Concrete compressive strength compress-stren 1,030 8 Concrete slump test slump 103 9 slump-comp, slump-flow 7 Condition based maintenance of naval propulsion plants cond-turbine 11,934 13 Conventional and Social Media Movies 14/15 csm1415 231 12/11 Cuff-less blood pressure estimation cuff-less 61,000 3/2 Daily Demand Forecasting Orders daily-demand 60 13/12 Dynamic features of VirusShare Executables dynamic-features 107,856 482/265 Energy efficiency energy-cool, energy-heat 768 8/7 Facebook comment volume facebook-comment 40,949 54/48 Facebook metrics facebook-metrics 500 19 Forestfires forestfires 517 12/39 Gas sensor array under dynamic gas-dynamic-CO 58 438/57 gas mixtures gas-dynamic-methane Geographical original of music geo-lat, geo-long 1,059 116/72 geo-music-lat, geo-music-long 68 GPS trajectories gps-trajectory 163 10 Greenhouse gas observing network greenhouse-net 955,167 15 Housing housing 452 13 Individual household electric power consumption household-consume 2,049,280 6/5 Istanbul stock exchange stock-exchange 536 8 KEGG metabolic reaction network (undirected) KEGG-reaction 65,554 27/25 KEGG metabolic relation network (directed) KEGG-relation 54,413 22/17 Online news popularity online-news 39,644 59/55 Online video characteristics and transcoding time dataset video-transcode 68,784 20/8 Continued on next page. 6 Table 2 – Continued from previous page. Original UCI name Datasets #patterns #inputs Parkinson speech dataset with multiple types of sound recordings park-speech 1,040 26 Parkinson’s telemonitoring park-motor-UPDRS, park-total- UPDRS 5,875 16 PM2.5 Data 5 Chinese Cities pm25-beijing-dongsi 24,237 13 pm25-beijing-dongsihuan 20,166 pm25-beijing-nongzhanguan 24,137 pm25-beijing-us-post 49,579 pm25-chengdu-caotangsi 22,997 pm25-chengdu-shahepu 23,142 pm25-chengdu-us-post 27,368 pm25-guangzhou-city-station 32,351 pm25-guangzhou-5th-middle-school 21,095 pm25-guangzhou-us-post 32,351 pm25-shanghai-jingan 22,099 pm25-shanghai-us-post 31,180 pm25-shanghai-xuhui 23,128 pm25-shenyang-taiyuanji 22,992 pm25-shenyang-us-post 20,452 pm25-shenyang-xiaoheyan 23,202 Physicochemical properties of protein tertiary structure physico-protein 45,730 9 Relative location of CT slices on axial axis CT-slices 53,500 385/355 Servo servo 167 4/15 SML2010 SML2010 4,137 20/18 Stock portfolio stock-abs, stock-annual, stock-excess 252 6 stock-rel, stock-systematic, stock-total Student performance student-mat 395 32/77 student-por 649 32/56 UJIIndoorLoc UJ-lat, UJ-long 21,048 528/373 Yacht hydrodynamics yacht-hydro 308 6 YearPredictionMSD year-prediction 2,000 90 Table 2: Collection of 83 datasets from the UCI repository. Each column reports: original name in the UCI repository; datasets created from the original one; number of patterns (or observations) and inputs, before and after preprocessing. 66 7 Excluded dataset Reason Amazon access samples Huge number of inputs (20,000) Breast cancer Wisconsin (Prognostic) Too few recurrent patterns (47) Cargo 2000 Freight Tracking and Tracing Less than 10 different output values (3) Challenger USA space shuttle O-ring Too few patterns (23) and inputs (3) Condition based maintenance of naval propulsion plants (compress output) Less than 10 different output values (9) Container crane controller Too few patterns (15) DrivFace Less than 10 different output values (4 subjects) Early biomarkers of Parkinsons disease based on natural connected speech Data are not available Educational process mining Inputs and output for regression are not clear ElectricityLoadDiagrams Huge number of inputs (140,256) Fertility Less than 10 different output values (2) Gas sensor array drift dataset at different concentrations Less than 10 different output values (7) Gas sensor array exposed to turbulent gas mixtures Huge number of inputs (150,000) Gas sensor array under flow modulation Less than 10 different output values (4) Geo-Magnetic field and WLAN Data format very complex Improved spiral test using digitized graphics tablet for monitoring parkinsons disease Data are not available Insurance Company Benchmark (COIL 2000) Less than 10 different output values (3) KDC-4007 dataset Collection Less than 10 different output values (8) KDD cup 1998 Format too complex Las Vegas Strip Less than 10 different output values (5) News popularity in multiple social media platforms Data are text instead of numbers Noisy office Format too complex (PNG images) Open university learning analytics Format too complex Paper Reviews Less than 10 different output values (5) Parkinson disease spiral drawings using digitized graphics tablet Less than 10 different output values (3) Skillcraft1 master table Less than 10 different output values (7) Solar flare Less than 10 different output values Tamilnadu electricity board hourly readings Less than 10 different output values (2) Tennis major tournament match statistics Format problems Twin gas sensor arrays Less than 10 different output values (4) UJIIndoorLoc-Mag Output almost constant, format very complex wiki4HE Less than 10 different output values (7) Wine quality (white/red) Less than 10 different output values (7/6) Table 1: List of the UCI regression datasets which were excluded from this study with the reason to be discarded. In datasets with discrete outputs the number of different output (or response) values is between parentheses. 8 Family Regression models Family Regression models Linear regression (LR) 1. lm [11] Least absolute shrinkage and selection operator (LASSO) 21. lasso [12] 2. rlm [13] 22. relaxo [14] Generalized linear regression (PLM) 3. penalized [15] 23. lars [16] 4. enet [12] Ridge 24. ridge [12] 5. glmnet [17] 25. spikeslab [18] 6. glmSAIC [19] 26. foba [20] Additive models (AM) 7. gam [21] Bayesian models (BYM) 27. bayesglm [22] 8. earth [23] 28. brnn [24] Least squares (LS) 9. nnls [25] 29. bMachine [26] 10. krlsRadial [27] Gaussian processes (SGP) 30. gprLin [28] Projection methods (PRJ) 11. spls [29] 31. gprRad [28] 12. simpls [30] 32. gprPol [28] 13. kpls [31] Quantile regression (QTR) 33. rqlasso [32] 14. wkpls [33] 34. rqnc [34] 15. enpls.fs [35] 35. qrnn [36] 16. plsRglm [37] Nearest neighbors (NN) 36. kknn [38] 17. ppr [39] Regression trees (RGT) 37. rpart [40] 18. pcr [41] 38. nodeHarvest [42] 19. icr [43] 39. ctree2 [44] 20. superpc [45] 40. partDSA [46] 41. evtree [47] Table 5: List of regression models and references grouped by families (see Appendix B for a brief description of each model). these cases, although the caret model list specifies that hyperparameter as177 tunable, in the practice only one value is used. The list of hyperparameter178 values which are used in our experiments for a model and dataset is included179 in the file results model dataset.dat, where model and dataset stand for180 the names of the regression model and dataset, respectively, which is placed181 in the directory results/dataset/model implem, where implem may be R,182 C,Python or Matlab. For some models (ridge,rlm,mlpWD,mlpWDml,dnn,183 krlsRad and icr), the value list provided by the getModelInfo function was184 not valid due to several reasons, so we directly specify the hyperparameter185 values used for tuning in the file programs/R/initialize.R. The models186 in Matlab, C++ and Python use pre-specified values, listed in the script187 15 Family Regression models Family Regression models Regression rules (RGR) 42. M5 [48] Boosting (BST) (continued) 60. gbm [49] 43. cubist [50] 61. blackboost [51] 44. SBC [52] 62. xgbTree [53] Random forests (RF) 45. rf [54] 63. xgbLinear [53] 46. Boruta [55] Neural networks (NET) 64. mlWD [56] 47. RRF [57] 65. mlWDml [56] 48. cforest [58] 66. avNNet [6] 49. qrf [59] 67. rbf [56] 50. extraTrees [60] 68. grnn [61] Bagging (BAG) 51. bag [62] 69. elm [63] 52. bagEarth [6] 70. kelm [63] 53. treebag [64] 71. pcaNNet [6] Boosting (BST) 54. rndGLM [65] 72. bdk [66] 55. BstLm [53] Deep learning (DL) 73. dlkeras [67] 56. bstSm [53] 74. dnn [68] 57. bstTree [53] Support vector regression (SVR) 75. svr [9] 58. glmboost [69] 76. svmRad [70] 59. gamboost [69] 77. rvmRad [71] Table 6: Continuation of Table 5. run model dataset.sh, which are the same for all datasets. Tables 5 and188 6 list the the collection of 77 regression models used in this work, grouped189 by families, which are described in Appendix B, specifying the software im-190 plementation (R package or other platforms), their tunable hyperparameters191 and the values used.192 3. Results and discussion193 The experimental work [72] uses the following methodology: for each194 dataset with less than 10,000 patterns, N= 500 random partitions are gen-195 erated, using the 50% of the patterns for training, 25% for validation (in196 hyperparameter tuning) and 25% for test. For each dataset with more than197 10,000 patterns, a 10-fold cross validation is developed, so there are N= 10198 training, validation and test percentages. The rationale is to limit the compu-199 16 tational overhead of 500 trials to smaller datasets, using a lighter methodol-200 ogy (10-fold), although statistically less significant, for larger datasets. Each201 regression model is trained on the training partitions for each combination202 of its hyperparameter values, and it is tested on its corresponding validation203 partition. The performance measures used are the root mean square error204 (RMSE), the squared correlation (R2) and the mean absolute error (MAE).205 We use these three different measures in order to give more significance to206 the results, which in this way can be observed from three different points of207 view, and also in order to evaluate whether they are coherent suggesting sim-208 ilar conclusions. For each combination of hyperparameter values, the average209 RMSE over the validation partitions is calculated, and the combination with210 the lowest average RMSE is selected for testing (quantile regression models211 as rqlasso,rqnc and qrnn are designed to optimize the quantile error, which212 is used instead of RMSE). Finally, the model is trained on the training parti-213 tions using the selected combination of its hyperparameter value and tested214 on the test partitions. The performance measurements are the RMSE, R2 215 and MAE between the true and predicted output values concatenated for the216 Ntest sets. Note that the R2is calculated using the predicted and true out-217 puts for the test patterns, while it is often used to measure the percentage of218 variance explained by the model on the training patterns. Those regression219 models which lack tunable hyperparameters are trained on the training parti-220 tions and tested on the corresponding test partitions, and the average RMSE,221 R2and MAE over the test partitions are the quality measurements. Some222 models which are specially sensitive to collinear inputs are trained, for each223 partition, using only those inputs which are not collinear. Although collinear224 17 inputs have been removed from the dataset in the initial preprocessing, for225 certain partitions some inputs in the training set may be collinear despite of226 being not collinear considering the whole dataset. To avoid the subsequent227 errors, these inputs are discarded for these models. All the continuous inputs228 and the output are standardized to have zero mean and standard deviation229 one, using the mean and deviation calculated in each training partition.230 We run this collection of 77 regression models over 83 datasets, developing231 a total of 6,391 experiments, which were developed on a cluster whose nodes232 are equipped with 64 Intel Xeon E5-2650L processors and 4 GB of RAM233 memory each processor, although those regression models which required234 more memory with large data sets were run using several processors and up235 to 128 GB of RAM memory. Since certain models failed for some datasets,236 we developed a preliminar study to evaluate the datasets according to their237 size, given by its population, and “difficulty”, estimated by the R2achieved238 by the linear regression model (lm). We selected lm because it is a classical239 approach which can be considered as a baseline reference for other models and240 it does not require large time nor memory, so it does not fail in any dataset.241 Figure 1 plots R2 lm for all the datasets vs. their populations Np. According242 to this plot, we divided the datasets into four groups: group SD includes243 20 datasets with Np<5,000 and R2 lm <0.6, i.e., small-difficult datasets;244 group SE includes 23 datasets with R2 lm ≥0.6 and Np<5,000 (small-easy245 datasets); group LD with 33 datasets where R2 lm <0.6 and Np≥5,000246 (large-difficult datasets); and group LE with 7 datasets where R2 lm ≥0.6 and247 Np≥5,000 (large-easy datasets). Table 7 lists the datasets of each group.248 In order to compare the R2values achieved by the regression models over249 18 101102103104105106 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 #patterns R2 lm Figure 1: Values of R2achieved by lm for all the datasets plotted against their populations (in logarithmic scale), dividing datasets in groups small-difficult (SD, lower left quarter of the figure), small-easy (SE, upper left), large-difficult (LD, lower right) and large-easy (LE, upper right). all the datasets, averaging would weight more those datasets with high R2,250 favoring models which perform well in easy datasets and biasing the results.251 In order to do a neutral comparison, the solution is to average over all the252 datasets the model positions in a ranking sorted by decreasing R2, instead of253 directly averaging R2values, because the model positions belong to the same254 range for all the datasets. This is done using the Friedman ranking [73] of255 the R2coefficient, which evaluates the position where each regression model256 is placed, in average over all the datasets, when R2is sorted by decreasing257 values. The R2Friedman ranking of the M= 77 models over the D= 83258 datasets can be calculated as follows. For each dataset d= 1,...,D, the R2 259 values of all the models are sorted decreasingly. For each model m= 1,...,M260 let pmd be its position in dataset d. The Friedman rank Fmof model mis261 defined as Fm=1 DPD d=1 pmd, i.e., the average position of model mover all262 19 the sortings of R2for the different datasets. For example, a model with rank263 5 achieves the 5th highest R2coefficient in average over all the datasets.264 A number of run-time errors happened for certain models and datasets.265 There are more errors in large datasets, because some model implementations266 may not be designed to process large amounts of data. When a model fails267 for a dataset (because it overcomes the maximum allowed time of 48 hours,268 because it requires more than 128 GB of RAM, or due to other reasons), and269 in order to calculate the Friedman ranking, its R2is intended to be zero,270 while its RMSE and MAE are assigned as:271 RMSE = max    max m∈R [RMSEm],v u u t 1 N N X i=1 (ti−¯ ti)2   (1) MAE = max (max m∈R [MAEm],1 N N X i=1 |ti−¯ t|)(2) where Ris the set of models which did not fail in that dataset, tiis the272 true output for test pattern iand Nis the number of test patterns. Besides,273 denoting by kthe test partition to which pattern ibelongs, ¯ tiis the mean of274 the true output values over the patterns in the k-th training partition. The275 rationale behind this is that a regression model which fails behaves as if it276 would predict the mean of the true output values for all the test patterns, so it277 should be the last of the list. For some models, the errors happen only during278 tuning for some partitions, which are not considered to calculate the average279 RMSE corresponding to that combination of hyper-parameter values. When280 a model fails for a given combination of hyper-parameter values and all the281 20 partitions, that combination is not selected for testing. When a model fails282 for all the combinations of hyper-parameter values, or when it fails for some283 test partition, the model is considered that fails for that dataset. Overall,284 the number of experiments where the model failed is 1,205 and represents285 18.85% of the 6,391 experiments.286 Group (#datasets) Datasets SD (20): small-difficult airfoil com-crime-unnorm csm1415 forestfires geo-lat geo-long geo-music- lat geo-music-long gps-trajectory park-speech slump slump-flow stock-abs stock-annual stock-excess stock-rel stock-systematic stock-total student-mat student-por SE (23): small-easy air-quality-CO air-quality-NMHC air-quality-NO2 air-quality-NOx airquality-O3 automobile auto-mpg bike-day com-crime com-hd compress-stren daily-demand energy-cool energy-heat facebook-metrics gas-dynamic-CO gasdynamic-methane housing servo slump-comp SML2010 stock-exchange yachthydro LD (33): large-difficult 3Droad appliances-energy beijing-pm25 blog-feedback buzz-twitter cuff-less dynamic-features facebook-comment greenhouse-net household-consume KEGG-relation online-news park-motor-UPDRS park-total-UPDRS physico-protein pm25-beijing-dongsi pm25-beijing-dongsihuan pm25- beijing-nongzhanguan pm25-beijing-us-post pm25-chengdu-caotangsi pm25- chengdu-shahepu pm25-chengdu-us-post pm25-guangzhou-5th-middle-school pm25-guangzhou-city-station pm25-guangzhou-us-post pm25-shanghai- jingan pm25-shanghai-us-post pm25-shanghai-xuhui pm25-shenyang- taiyuanji pm25-shenyang-us-post pm25-shenyang-xiaoheyan video-transcode year-prediction LE (7): large-easy bike-hour combined-cycle cond-turbine CT-slices KEGG-reaction UJ-lat UJ- long Table 7: Groups of datasets according to its size (small/large) and complexity (easy/difficult). 3.1. Discussion by dataset group287 Table 8 reports the 20 best regression models according to the Friedman288 ranking of R2, RMSE and MAE for the datasets of group SD, which in-289 cludes 20 small-difficult datasets. For each model in the R2ranking, the290 percentage of datasets where it failed is also reported (column %Error). The291 last two columns report the models which achieved the best R2for some292 dataset and the percentage of datasets where this happened. First of all,293 21 penalized achieves the first positions in the three rankings, being the best294 R2for 25% of datasets. ExtraTrees,rf and kelm are the following in the295 R2ranking, although the former achieves much lower positions in RMSE and296 MAE rankings. Specifically, extraTrees achieves the best R2in 40% of the297 datasets, although it fails in 5% of them, so it can be considered less regu-298 lar as penalized. Other good positions in the R2ranking are for qrf and299 bstTree, followed by avNNet,svr and svmRad, Gaussian process (gprRad300 and gprPol), bagEarth and cubist, which achieves the best R2for 10% of301 datasets.302 R2RMSE MAE Best R2 Pos. Model Rank %Error Model Rank Model Rank Model %Best 1 penalized 8.45 0.0 penalized 9.65 penalized 13.40 extraTrees 40.0 2 extraTrees 13.05 5.0 kelm 13.15 svmRad 13.90 penalized 25.0 3 rf 15.35 5.0 gprPol 14.45 svr 15.70 cubist 10.0 4 kelm 19.15 5.0 bagEarth 17.55 kelm 16.25 brnn 10.0 5 qrf 20.75 0.0 svmRad 18.00 bstTree 19.15 rbf 5.0 6 bstTree 21.00 0.0 cforest 18.65 gprPol 19.25 qrf 5.0 7 avNNet 21.25 5.0 bstTree 19.10 cubist 19.65 bagEarth 5.0 8 svr 22.20 10.0 svr 19.35 bagEarth 19.75 — — 9 svmRad 23.15 5.0 enet 21.50 cforest 21.45 — — 10 gprRad 23.20 0.0 BstLm 22.75 qrf 23.35 — — 11 RRF 24.10 20.0 glmboost 23.80 avNNet 23.85 — — 12 bagEarth 24.10 0.0 gbm 24.20 gbm 24.35 — — 13 gprPol 24.35 0.0 foba 24.70 grnn 27.00 — — 14 gbm 24.60 0.0 bMachine 25.30 gprRad 28.35 — — 15 cubist 26.20 5.0 grnn 26.25 extraTrees 29.05 — — 16 ridge 27.85 0.0 spls 27.25 rf 29.15 — — 17 treebag 29.45 0.0 spikeslab 27.25 BstLm 29.45 — — 18 foba 29.55 0.0 rf 27.90 treebag 29.50 — — 19 spls 29.85 0.0 lars 28.35 rqlasso 29.75 — — 20 lars 30.25 0.0 avNNet 28.90 glmboost 29.95 — — Table 8: List of the 20 best regression models according to the Friedman rank of R2, RMSE and MAE for dataset group SD, with 20 small-difficult datasets. The last two columns list the models which achieve the best R2for some dataset, sorted by decreasing number of datasets. Since this group includes only small datasets, most models exhibit low303 error percentages (i.e., most models never or rarely fail on datasets of this304 22 group), although some models with errors achieve good positions, e.g. svr305 (10% of errors), RRF (20%), extraTrees and cubist (5% each one). Besides,306 qrnn and nodeHarvest are very slow and they were shutdown after 48 h. for307 the 20 datasets of this group. Considering memory errors, rndGLM is the308 model which requires more memory, overcoming the memory and time limits309 in 1 and 5 datasets of this group, respectively.310 R2RMSE MAE Best R2 Pos. Model Rank %Error Model Rank Model Rank Model %Best 1 cubist 6.48 0.0 cubist 6.26 cubist 5.65 cubist 21.7 2 avNNet 10.13 4.3 avNNet 10.04 avNNet 10.96 avNNet 13.0 3 bstTree 12.57 0.0 bstTree 12.39 bstTree 13.70 extraTrees 13.0 4 gbm 12.87 0.0 gbm 12.74 ppr 13.91 gbm 8.7 5 bagEarth 14.57 0.0 bagEarth 14.30 gbm 13.96 penalized 8.7 6 ppr 14.65 0.0 ppr 14.52 bagEarth 15.96 bMachine 8.7 7 bMachine 14.96 4.3 bMachine 15.22 bMachine 16.70 kelm 4.3 8 extraTrees 17.13 8.7 earth 18.35 M5 16.83 M5 4.3 9 earth 18.57 0.0 kelm 18.65 qrf 17.43 rf 4.3 10 kelm 18.70 17.4 extraTrees 18.87 extraTrees 18.00 brnn 4.3 11 rf 19.26 4.3 rf 19.65 kelm 18.74 bagEarth 4.3 12 M5 20.30 0.0 M5 19.87 rf 19.70 bstTree 4.3 13 RRF 22.43 8.7 RRF 22.61 earth 21.00 — — 14 qrf 22.48 0.0 qrf 23.43 brnn 22.48 — — 15 brnn 23.65 21.7 brnn 23.61 RRF 22.83 — — 16 gprPol 24.00 4.3 gprPol 24.22 pcaNNet 24.17 — — 17 pcaNNet 25.04 0.0 pcaNNet 25.09 gprPol 25.78 — — 18 dlkeras 27.35 0.0 dlkeras 27.61 rqlasso 26.52 — — 19 Boruta 27.43 17.4 Boruta 27.78 cforest 27.13 — — 20 enet 28.52 0.0 enet 28.09 Boruta 27.70 — — Table 9: List of the 20 best regression models according to the Friedman rank of R2, RMSE and MAE over 23 datasets of group SE (small-easy). Considering small-easy datasets (group SE, 23 datasets, table 9), the311 three rankings are even more coherent than for group SD, sharing the first312 three positions: cubist, which achieves the best R2for 21.7% of datasets,313 avNNet (the best R2for 13% of datasets) and bstTree.Penalized is not314 present in this list (although it is the best in 8.7% of datasets), but gbm315 and bMachine (which are the bests in 8.7% of datasets), bagEarth,ppr,316 23 extraTrees (the best in 13% of datasets), earth and kelm are in positions 4-317 10. Other models with good results are rf,M5 (the best for 4.3% of datasets),318 RRF and qrf.319 R2RMSE MAE Best R2 Pos. Model Rank %Error Model Rank Model Rank Model %Best 1 M5 9.48 0.0 M5 9.39 M5 9.55 extraTrees 48.5 2 cubist 12.39 15.2 gbm 12.61 kknn 12.55 bstTree 12.1 3 gbm 12.48 3.0 cubist 12.70 cubist 12.61 cubist 9.1 4 xgbTree 14.24 6.1 xgbTree 14.15 gbm 13.42 dlkeras 6.1 5 kknn 14.48 12.1 kknn 14.48 xgbTree 14.82 xgbTree 6.1 6 bstTree 15.12 12.1 bstTree 15.27 bstTree 16.00 ppr 3.0 7 blackboost 16.79 0.0 blackboost 17.27 grnn 17.33 kknn 3.0 8 dlkeras 18.36 15.2 pcaNNet 18.33 blackboost 17.70 M5 3.0 9 svr 18.58 27.3 svr 18.45 svr 18.76 rf 3.0 10 pcaNNet 18.76 0.0 dlkeras 18.67 pcaNNet 18.82 qrf 3.0 11 grnn 19.70 18.2 ppr 19.48 dlkeras 19.18 bMachine 3.0 12 ppr 19.73 3.0 grnn 19.52 ppr 19.58 — — 13 qrf 21.33 27.3 qrf 21.27 qrf 20.30 — — 14 svmRad 21.88 24.2 svmRad 21.79 svmRad 20.58 — — 15 earth 22.52 0.0 earth 22.21 extraTrees 22.33 — — 16 extraTrees 23.03 27.3 bag 22.91 bag 22.61 — — 17 bag 23.03 15.2 avNNet 23.42 earth 22.88 — — 18 avNNet 23.52 21.2 extraTrees 23.91 avNNet 23.64 — — 19 bMachine 24.76 24.2 bMachine 24.64 bMachine 25.24 — — 20 cforest 25.79 27.3 cforest 25.91 rpart 26.00 — — Table 10: List of the 20 best regression models according to the Friedman rank of R2, RMSE and MAE for the 33 datasets of group LD (large-difficult). In large-difficult datasets (group LD, 33 datasets, table 10), the M5320 achieves the first positions in the three rankings (although it achieves the best321 R2only in 3% of datasets), followed by cubist (which achieves the best R2 322 and errors in 9.1% and 15.2% of datasets, respectively) and gbm. Other mod-323 els with good performance are xgbTree,knn,bstTree,blackboost,dlkeras324 (15.2% of errors), svr (with errors in 27.3% of datasets) and pcaNNet. The325 high error frequency of several models (either by overcoming limits on mem-326 ory or time) is due to the large size of datasets in this group. ExtraTrees also327 overcomes the maximum time in 27.3% of datasets and, as in groups SD and328 24 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 10 20 30 40 50 60 70 80 90 100 Threshold (θ) for ∆ Percentage of datasets where ∆ = R2 best − R2 > θ cubist (7 errors) gbm (1 error) bstTree (5 errors) M5 (0 errors) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 cubistgbm bstTree M5 avNNet extraTrees qrf pcaNNet rf bMachine ∆=R2 best−R2 Figure 3: Left panel: percentage of datasets where the difference ∆ = R2 best −R2, with R2 achieved by cubist,gbm,bstTree and M5, overcomes a given threshold θ. Right panel: boxplots of the differences R2 best −R2for the 10 best regressors. respectively) causes that blue and black lines fall to zero for θ > 0.4 (outside416 the plot), while green and red lines (gbm and M5, respectively) fall to zero417 at 0.322 and 0.338. Note that gbm fails (achieving R2= 0) only for dataset418 year-prediction, for which by chance R2 best is low (0.338), so ∆ = 0.338 for419 this dataset and the green curve falls to zero at θ= 0.338. If the R2 best were420 higher, the green line would continue to the right without falling to zero,421 similarly to blue and black lines. The right panel of Figure 3 shows the422 boxplots of the differences R2 best −R2for the first 10 models in the global423 ranking. The blue boxes report the 25% and 75% quantiles, while the red424 line inside the box is the median, and the blue points outside the box are the425 outliers. Although cubist,gbm,bstTree,extraTrees and rf exhibit low426 medians, all the models have several outliers (caused by datasets where they427 fail) with high ∆ values, excepting M5, the only one which guarantees low ∆428 values (below 0.322) for all the datasets.429 Figure 4 (left panel) plots R2 best and the R2achieved by M5 and gbm.M5 is430 31 0 10 20 30 40 50 60 70 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dataset R2 gbm M5 Best −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 5 10 15 20 25 δ = R2 M5 − R2 gbm #Datasets δ > 0 for 53% of datasets Σδ > 0 δ = 2.3 Σδ ≤ 0 δ = −1.3 Figure 4: Left panel: R2 gbm (blue), R2 M5(red) and R2 best (green) for each dataset, sorted by decreasing values of R2 gbm values. Right panel: histogram of the difference R2 M5−R2 gbm. near R2 best more often than gbm, and in several cases gbm is clearly below M5,431 but the former rarely outperforms the latter, and in these cases with lower432 difference. The right panel shows the histogram of the difference R2 M5−R2 gbm:433 its values are positive (i.e., M5 outperforms gbm) for 52.4% of the datasets,434 and when they are positive, they are higher (in absolute value) than when435 they are negative, so the sum of positive ∆ values (2.3) outperforms the436 sum of negative values (-1.3). This shows that overall M5 outperforms gbm,437 although the latter is higher in the global ranking (Table 13). Remember438 that cubist,gbm and bstTree fail for some datasets, while M5 never fails.439 In the left panel of Figure 4, the maximum difference R2 best −R2 M5is 0.322440 in dataset student-mat, whose output is discrete with more than 10 values441 (see the left panel of Figure 5), so the dataset was not excluded. The low442 R2 best = 0.3673 for this dataset means that no model, and not only M5, fits443 accurately the true output, and both blue and green points fit equally bad the444 red line. The right panel of the same figure plots the difference R2 best −R2 M5 445 against R2 best. This difference is low for all the datasets (note that the vertical446 32 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 Predicted output True output R2 best=0.3673 (penalized) R2 M5=0.0450 M5 penalized 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 R2 best R2 best − R2 M5 Figure 5: Left panel: true against predicted output for M5 (blue points) and the best regression model (penalized) for dataset student-mat. Right panel: difference R2 best−R2 M5 against R2 best for all the datasets. scale is 0–0.4), being below 0.2 (resp. 0.1) for 92.8% (resp. 60.2%) of the447 datasets.448 10 20 30 40 50 60 70 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dataset R2 lm rlm R2 best 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1−R2 lm R2 M5 R2 M5 R2 best Figure 6: Left: values of R2 lm,R2 rlm and R2 best, sorted by decreasing R2 lm. Right: R2 M5and R2 best against 1 −R2 lm. The left panel of Figure 6 compares R2 lm,R2 rlm and R2 best in all datasets.449 Since both models only differ in the robustness against outliners, the dif-450 ference between them identifies those datasets with outliers. This difference451 overcomes 0.05 only in 6 datasets and its highest value is 0.31, so that dataset452 33 outliers are few and not very relevant. In order to study the behavior of M5453 with the dataset complexity, the right panel shows R2 best and R2 M5against454 1−R2 lm, which measures the difficulty of the regression problem. The differ-455 ence R2 best −R2 M5, instead of raising with 1−R2 lm, achieves the highest values456 for 0.65 <1−R2 lm <0.9. However, in the most difficult datasets, where457 1−R2 lm >0.9, the R2 M5follows very well R2 best, so M5 performs well even for458 hard datasets.459 10 20 30 40 50 60 70 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dataset R2 dlkeras M5 cubist 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 R2 max cubist M5 dlkeras R2 Figure 7: Left: values of R2achieved by dlkeras,M5 and cubist (datasets sorted by decreasing R2 max). Right: boxplots of R2 best and R2achieved by cubist,M5 and dlkeras over all the datasets. Figure 7 (left panel) compares dlkeras to M5 and cubist over all the460 datasets (points with R2= 0 correspond to datasets where cubist or dlkeras461 fail). The cubist model (green line) achieves almost always the highest per-462 formance (in 61 of 83 datasets), while M5 overcomes cubist and dlkeras463 in 18 datasets, and dlkeras only in 4 datasets. In fact, cubist outper-464 forms dlkeras in 71 datasets, while dlkeras outperforms cubist only in 8465 datasets. The difference between dlkeras and M5 is lower (44 and 39 datasets466 favoring M5 and dlkeras, respectively. The right panel of Figure 7 shows the467 34 boxplots of R2 best and R2of cubist,M5 and dlkeras: this last box is clearly468 below cubist and M5, but its median is similar to M5 and lower than cubist.469 The upper box ends of cubist and M5 are near to R2 best, and the median470 of cubist is also very near to R2 best, but the lower box ends of cubist and471 dlkeras are much below R2 best and M5 due to errors. Analyzing the param-472 eter tuning of dlkeras, the largest available size (753= 421,875 neurons in473 three hidden layers) was selected only in 17 of 83 datasets. Therefore, in the474 remaining 66 datasets R2did not increase with larger networks, so they are475 not expected to provide better performances. However, they would spend476 higher computation times, overcoming the maximum allowed time (48 h.)477 more frequently, so dlkeras would achieve more errors than cubist and M5,478 which never fails.479 Small-difficult Small-easy Large-difficult Large-easy Family-model Pos. Family-model Pos. Family-model Pos. Family-model Pos. PLM-penalized 1 RGR-cubist 1 RGR-M5 1 RGR-M5 1 RF-extraTrees 2 NET-avNNet 2 BST-gbm 3 BST-gbm 3 NET-kelm 4 BST-bstTree 3 NN-kknn 5 BAG-bag 4 BST-bstTree 6 BAG-bagEarth 5 DL-dlkeras 8 NET-pcaNNet 7 SVR-svr 8 PRJ-ppr 6 SVR-svr 9 LR-lm 9 SGP-gprRad 10 BYM-bMachine 7 NET-pcaNNet 10 AM-earth 10 RF-extraTrees 8 AM-earth 9 Table 16: Best model of each family within the 10 best positions in the R2Friedman ranking for each dataset group. 3.3. Discussion by family of regression model480 It is interesting to analyze the behavior of the best model of each family.481 Table 16 reports the families with models in the top-10 of the R2ranking for482 each dataset group. The boosting family (BST, with models bstTree and483 35 gbm) and neural networks (NET, models kelm,avNNet and pcaNNet) families,484 are present in all the groups, while regression rules (RGR), with models485 cubist and M5, achieves the first position in three of four groups (small-easy,486 large-difficult, large-easy), and penalized linear regression (PLM) achieves487 the first position (penalized) in small-difficult datasets. Bagging (BAG,488 with models bag and bagEarth) and support vector regression (SVR, svr)489 are present in two groups, while RF (extraTrees), projection methods (PRJ,490 ppr), Gaussian processes (SGP, gprRad), nearest neighbors (NN, kknn), deep491 learning (DL, dlkeras) and linear regression (LR, lm and rlm) are only492 present in just one group.493 Family Best Pos. Family Best Pos. Regression rules cubist 1 Nearest neighbors kknn 16 Boosting gbm 2 Penalized linear models penalized 17 Neural networks avNNet 5 Deep learning dlkeras 18 Random forests extraTrees 6 Ridge foba 27 Bayesian models bMachine 10 Lasso lars 28 Bagging bagEarth 11 Linear regression lm 33 Support vector regression svr 12 Regression trees ctree2 37 Projection methods ppr 13 Gaussian processes gprPol 50 Generalized additive models earth 14 Quantile regression rqlasso 56 Table 17: Best regression model of each family and position in the global ranking. Considering the global ranking, Table 17 reports the families, sorted by494 the position of their best models in Table 13. Only regression rules, boost-495 ing and neural networks are in the top-5, followed by random forests and496 Bayesian models with positions below 10. Most of the remaining families497 have best models which outperform lm (position 33), while regression trees,498 Gaussian processes and quantile regression achieve positions even higher.499 36 0 0.2 0.4 0.6 0.8 1 R2 airfoil gps−trajectory slump−flow com−crime−unnorm slump student−por geo−lat student−mat geo−long geo−music−lat geo−music−long park−speech stock−annual csm1415 stock−excess stock−rel stock−total forestfires stock−systematic stock−abs Best penalized lm (a) Small-difficult (SD). 0 0.2 0.4 0.6 0.8 1 R2 bike−day air−quality−NMHC yacht−hydro com−hd energy−heat SML2010 slump−comp energy−cool facebook−metrics daily−demand air−quality−NO2 air−quality−CO air−quality−O3 compress−stren air−quality−NOx automobile auto−mpg servo housing stock−exchange gas−dynamic−methane gas−dynamic−CO com−crime Best cubist lm (b) Small-easy (SE). 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Data set R2 Best M5 lm (c) Large-difficult (LD). 0 0.2 0.4 0.6 0.8 1 R2 bike−hour cond−turbine UJ−long CT−slices UJ−lat KEGG−reaction combined−cycle Best M5 lm (d) Large-easy (LE). Figure 8: Best R2(in green), R2achieved by the model with the best R2Friedman rank (in blue), and R2achieved by lm (in red) for each group. 37 3.4. Best result for each dataset500 The green line of Figure 4 shows R2 best for the 83 datasets. For 41 of501 them (which represents 49.39%) the R2 best >0.8, so that at least one model502 was able to predict the output with an acceptable accuracy (alternatively,503 these datasets might be considered as “easy” to learn). Besides, for other504 20 datasets (24.09% of the total) the R2 best is between 0.6 and 0.8, which505 is still an acceptable accuracy (datasets with middle difficulty). However,506 R2 best is between 0.2 and 0.6 for 13 datasets, which represents 15.66% (hard507 datasets), while R2 best <0.2 for 9 datasets (10.84%), where the models could508 not learn the regression problem at all. Some datasets are really hard, e.g.509 stock-abs, where R2 best = 0.059. Figure 8 plots the best R2, alongside with510 R2achieved by the best model and by lm for each dataset group. In group511 SD (Figure 8a), the best model (penalized) is near to the best R2except512 for datasets airfoil,gps-trajectory,slump-flow and slump. Since this group513 includes small-difficult datasets, the R2of lm is always below 0.6, but for the514 first five datasets some model achieves higher R2values. The penalized is515 also better than lm for all datasets, although the difference is low for datasets516 after geo-music-long. For group small-easy (SE, Figure 8b), the R2values of517 lm are higher, but the best model (cubist) is always very near to the best518 R2with some difference with respect to lm. In the large-difficult group (LD,519 Figure 8c), the lm values are very low and the best model (M5) is far from520 lm, following the best R2very closely for 12 of 33 datasets with a margin521 of 0.2-0.4 for the remaining 21 datasets. Finally, in group LE (large-easy522 datasets, Figure 8d) the lm is already near the best R2, although the best523 model (M5 again) always achieves the best R2.524 38 3.5. Discussion by elapsed times and memory consumption525 We studied the memory and time required by each regression model over526 all the datasets. Table 18 reports the information of the 20 best models ac-527 cording to the R2Friedman rank in each column: %Best reports the percent-528 age of datasets where they achieve the best R2;%Error reports the percentage529 of datasets where they failed; %ME reports the percentage of memory errors530 caused by overcoming the largest allowed memory (128 GB); Datasets/mem531 reports the number of datasets run on the memory queues with {2i}6 i=1 GB.532 In order to measure the time required by each model, the time spent in533 hyper-parameter tuning is discarded to avoid biasing caused by differences534 among models in the number of hyper-parameters and hyper-parameter val-535 ues. The column %TE reports the percentage of time errors, i.e., datasets536 where the model overcomes the maximum allowed time (48 h.). Although it537 may surprise that some models are not able to finish within 48 h., we must538 consider the size of some datasets (more than 2 millions of patterns, up to539 640 inputs) and the high number of trials (500) for some datasets. Gen-540 erally, high values in the %TE column happen with slow models, specially541 for large datasets. Since some models fail but do not overcome the allowed542 memory nor time, the sum of columns %ME and %TE is not always equal543 to column %Error, e.g. nnls has no memory nor time errors, but %Error544 is 4.8%. The column Time reports the time (in sec.) spent by the model545 for a training+testing trial on dataset compress-stren, whose size might be546 considered “standard”: 1,030 patterns and 8 inputs. The time is set to the547 maximum allowed time for models with errors in this dataset.548 Comparing cubist,gbm,bstTree and M5 in terms of column %Best,549 39 #Datasets/mem(GB) Pos. Model Rank %Best %Error %ME 2-4-8-16-32-64-128 %TE Time 1 cubist 13.92 15.7 8.4 68-6-8-1-0-0-0 8.4 2.47 2 gbm 15.42 2.4 1.2 78-3-2-0-0-0-0 1.2 1.46 3 bstTree 15.84 6.0 6.0 74-5-3-1-0-0-0 6.0 3.84 4 M5 18.20 2.4 0-36-29-18-0-0-0 1.32 5 avNNet 19.23 3.6 14.5 77-3-2-0-0-1-0 14.5 3.20 6 extraTrees 19.61 33.7 19.3 72-9-2-0-0-0-0 20.5 3.62 7 qrf 22.41 2.4 14.5 62-12-4-1-3-1-0 14.5 3.99 8 pcaNNet 23.49 77-3-3-0-0-0-0 1.36 9 rf 23.82 3.6 24.1 69-6-2-0-5-1-0 24.1 3.05 10 bMachine 23.83 3.6 15.7 54-27-2-0-0-0-0 16.8 12.10 11 bagEarth 24.14 2.4 7.2 76-1-4-1-1-0-0 7.2 2.21 12 svr 24.17 27.7 78-5-0-0-0-0-0 27.7 172800 13 ppr 24.57 1.2 4.8 78-3-2-0-0-0-0 4.8 1.24 14 earth 25.52 77-4-2-0-0-0-0 1.46 15 blackboost 25.69 75-1-3-0-3-0-1 1.48 16 kknn 26.24 1.2 6.0 80-2-1-0-0-0-0 6.0 2.41 17 penalized 27.70 8.4 12.0 77-4-2-0-0-0-0 12.0 1.54 18 dlkeras 28.07 3.6 7.2 83-0-0-0-0-0-0 7.2 6.17 19 svmRad 29.14 28.9 77-2-3-0-0-1-0 28.9 172800 20 grnn 29.61 9.6 6.0 47-13-6-5-5-1-1 3.6 0.22 28 lars 33.16 79-2-2-0-0-0-0 0.03 77 qrnn 77.00 100.0 77-4-2-0-0-0-0 100.0 172800 63 rndGLM 51.80 51.8 44.6 0-0-0-12-22-10-2 7.2 2.54 Table 18: List of the 20 first regression models sorted by increasing R2Friedman rank, with the percentage of datasets where each model achieves the best R2(column %Best), percentage datasets with errors (column %Error), percentage of memory errors (column %ME), number of datasets for each memory size (column #Datasets/mem), percentage of time errors (%TE) and training+test time (in sec.) spent for dataset compress-stren (column Time). Empty cells correspond to zero values. cubist achieves often the best R2(15.7% of datasets) followed by bstTree550 (6%) while gbm and M5 tie (2.4%). Cubist and bstTree fail in 8.4% and 6%551 of datasets, respectively, while gbm fails less (the three overcome the allowed552 time) and M5 never fails. None of them overcomes the memory limits, but M5553 requires more memory (4-16 GB), while the others require 2-8 GB. Finally,554 M5 and gbm are faster (1.32 and 1.46 s./trial, respectively), while bstTree555 40 has 41,757 patterns, despite its description in the UCI documentation spec-678 ifies 43,824 because 2067 patterns whose output is missing, so it can not679 be predicted, were removed. The cuff-less dataset has 73,200,000 patterns,680 while its description specifies 12,000. Instead of discarding it, we used the681 first 61,000 patterns. The greenhouse-net dataset has 2,921 files with 327682 patterns per file, which gives 955,167 patterns instead of 2,921 in the dataset683 description. The household-consume dataset has 2,049,280 patterns instead684 of 2,075,259 as listed in the UCI documentation, because 25,979 original pat-685 terns have missing values (labeled as ‘?’) for all the inputs and the output.686 The online-news dataset has 39,644 patterns instead of 39,797 as listed in687 the UCI documentation. For the UJIIndoorLoc datasets, output floor was688 discarded and did not give a separate regression dataset because it has only689 three different values.690 Appendix B. Listing of regression methods691 This appendix describes the regression model used in the current work,692 grouped by families, alongside with their software implementations and val-693 ues of their tunable hyper-parameters. Default values are assumed for all the694 model parameters not cited explicitly.695 I. Linear regression (LR)696 1. lm is the linear regression model implemented by the stats package697 [11]. Collinear inputs exhibit undefined coefficients in the linear regres-698 sion model returned by lm, being discarded by it and by other models699 in the list, as we told above.700 47 2. rlm implements the robust linear model (MASS package), fitted using701 iteratively re-weighted least squares with maximum likelihood type es-702 timation, which is robust to outliers in the output although not in703 inputs [13]. The only hyperparameter is the Ψ function, which can be704 huber (Huber function, which leads to a convex optimization problem),705 hampel and Tukey bisquare, both with local minima. In our experi-706 ments, these functions are selected as the best Ψ for 16%, 82% and 2%,707 respectively, of the datasets.708 II. Penalized linear regression (PLM)709 3. penalized is the penalized linear regression (penalized package), which710 fits generalized linear models with a combination of L1 and L2 penal-711 ties. The L1 penalty, also named LASSO, penalizes the sum of absolute712 values of the coefficients, thus reducing the coefficients of inputs which713 are not relevant, similarly to input selection. The L2 penalty (also714 named ridge) penalizes the sum of squared coefficients, reducing the715 effects of input collinearity. The regression is regularized by weighting716 both penalties [15], whose weights are given by hyperparameters λ1,717 tuned with values 1, 2, 4, 8 and 16, and λ2, with values 1, 2, 4 and 8.718 In our experiments, λ1=λ2= 1 in the 87.9% of the datasets, and only719 in 10 of 83 datasets λ16= 1 or λ26= 1.720 4. enet is the elastic-net regression model (elasticnet package), com-721 puted using the least angle regression - elasticnet (LARS-EN) algorithm722 [12]. Elastic-net provides a model for regularization and input selec-723 tion, grouping together the inputs which are strongly correlated. This724 model is specially useful when the number of inputs is higher than725 48 the number of patterns, as opposed to LASSO models. There are two726 hyperparameters (5 values each one): the quadratic penalty, or regular-727 ization, hyperparameter (λ, with values 0, {10−i}3 1) and the fraction s728 of the L1 norm of the coefficient vector relative to the norm at the full729 least squares solution (the fraction mode is used in the predict.enet730 function, with values 0.05, 0.28, 0.52, 0.76, 1).731 5. glmnet is the LASSO and elastic-net regularization for generalized732 linear models (GLM) implemented in the glmnet package [17]. The733 glmnet model uses penalized maximum likelihood to fit a GLM for the734 LASSO and elastic-net non-convex penalties. The mixing percentage735 αis tuned with 5 values from 0.1 to 1: the value α=1 (resp. <1) cor-736 responds to the LASSO (resp. elastic-net) penalty. The selected value737 for αduring hyperparameter tuning was 0.1 in 41.7% of the datasets.738 The regularization hyperparameter λis also tuned with values 0.00092,739 0.0092 and 0.092.740 6. glmSAIC is the generalized linear model with stepwise feature se-741 lection [19] using the Akaike information criterion and the stepAIC742 function in the MASS package (model glmStepAIC in the caret model743 list).744 III. Additive models (AM)745 7. gam is the generalized additive model (GAM) using splines (mgcv pack-746 age). This model [21] is a GLM whose linear predictor is a sum of747 smooth functions (penalized regression splines) of the covariates. The748 estimation of the spline parameters uses the generalized cross valida-749 tion criterion. The only hyperparameter is select, a boolean flag that750 49 adds an extra penalty term to each function penalizing its wiggliness751 (waving).752 8. earth is the multivariate adaptive regression spline (MARS) in the753 earth package. This method [23] is a hybrid of GAM and regression754 trees (see family XII) which uses an expansion of product spline func-755 tions to model non-linear data and interactions among inputs. The756 spline number and parameters are automatically determined from the757 data using recursive partitioning, and distinguishing between additive758 contributions of each input and interactions among them. The func-759 tions are added iteratively to reduce maximally the residual, until its760 change is too small or a number of iterations is reached. The maximum761 number of terms in the model (nprune) is tuned with 15 values (less762 for some datasets) between 2 and 24.763 IV. Least squares (LS)764 9. nnls is the non-negative least squares regression (nnls package), which765 uses the Lawson-Hanson NNLS method [25] to solve for xthe optimiza-766 tion problem minx|Ax −b|subject to x≥0, where Ais the input data767 matrix, bis the true output and xis the linear predictor.768 10. krlsRad is the radial basis function kernel regularized least squares769 regression (KRLS package), which uses Gaussian radial basis functions770 to learn the best fitting function which minimizes the squared loss of a771 Tikhonov regularization problem [27]. The KRLS method, which cor-772 responds to the krlsRadial in the caret model list, learns a closed773 form function which is so interpretable as ordinary regression models.774 The only hyperparameter is the kernel spread (σ), with 10 values in the775 50 set {10i}2 −7. By default, this method determines the trade-off between776 model fit and complexity, which is defined by the λparameter, by min-777 imizing the sum of squared leave-one-out errors. The getModelInfo778 function only lists one value for λ, despite being listed as a tunable779 hyperparameter in the caret model list.780 V. Projection methods (PRJ)781 11. spls is the sparse partial least squares regression (spls package). This782 method [29] uses sparse linear combinations of the inputs in the dimen-783 sionality reduction of PLS in order to avoid lack of consistency of PLS784 with high dimensional patterns. The hyperparameters are the number785 of latent components (K), with values 1, 2 and 3, and the threshold (η),786 with 7 values from 0.1 to 0.9.787 12. simpls fits a PLS regression model with the simpls method [30], imple-788 mented by the plsr function in the pls package, with method=simpls.789 The PLS method projects the inputs and the output to a new space and790 it searches the direction in the input space which explains the maxi-791 mum output variance. Simpls is particularly useful when there are792 more inputs than patterns and inputs are collinear. It directly calcu-793 lates the PLS factors as linear combinations of the inputs maximizing a794 covariance criterion with orthogonality and normalization constraints.795 The only hyperparameter is the number of components (ncomp) used796 by the simpls model, with values from 1 to min(10,#inputs-1).797 13. kpls is the PLS regression with method=kernelpls [31] in the same798 function and package as simpls, using the same hyperparameter setting799 as simpls with 6 values. This is the model named kernelpls in the800 51 caret model list.801 14. wkpls uses method=widekernelpls [33] for PLS, tuning the number of802 components (ncomp) as simpls also with 10 values (model widekernelpls803 in the caret model list).804 15. enpls.fs is an ensemble of sparse partial least squares (spls, see model805 #12) regression models implemented by the enpls package [35]. The806 getModelInfo function lists only one value for the number of com-807 ponents (maxcomp), while the threshold argument, specified as a hy-808 perparameter by the caret model list, is missing in the enpls.fit809 function.810 16. plsRglm is the partial least squares generalized linear model (plsRglm811 package) with modele=pls-glm-gaussian [37]. The hyperparameters812 are the number of extracted components (nt), tuned with values 1, 2,813 3 and 4, and the input significance level (alpha.pvals.expli), with814 values in the set {10i}2 −2.815 17. ppr performs the projection pursuit regression (stats package), which816 models the output as a sum of averaging functions (mean, median,817 etc.) of linear combinations of the inputs [39]. The coefficients are818 iteratively calculated to minimize a projection pursuit (fitting criterion,819 given by the fraction of unexplained variance which is explained by each820 function) until it falls below a predefined threshold. The only tunable821 hyperparameter is the number of terms of the final model (nterms),822 with values from 1 to 10.823 18. pcr develops principal component regression (pls package), which mod-824 els the output using classical linear regression with coefficients esti-825 52 mated with principal component analysis (PCA), i.e., using the prin-826 cipal components as inputs [41]. It works in three stages: 1) performs827 PCA and selects a subset of the principal components; 2) uses ordinary828 least squares to model the output vector using linear regression on the829 selected components; 3) uses the eigenvectors corresponding to the se-830 lected components in order to calculate the final pcr estimator trans-831 forming the modeled output vector to the original space, and estimates832 the regression coefficients for the original outputs. The number of com-833 ponents (ncomp) is tuned with values from 1 to min(10,#inputs-1).834 19. icr is the independent component regression (caret package). The835 icr fits a linear regression model using independent component analy-836 sis (ICA), implemented by the fastICA package, instead of the original837 inputs [43]. The input data are considered a linear combination of a838 number of independent and non-Gaussian components (sources), so the839 training set matrix is written as the product of the source matrix and a840 linear mixed matrix, which contains the coefficients of the linear com-841 bination. The ICA estimates a “separating” matrix, which multiplied842 by the original data, provides the sources. This matrix must maxi-843 mize the non-Gaussianity of the sources, measured by the neg-entropy.844 The only hyperparameter is the number of independent components845 n.comp, with values from 1 to min(10,#inputs-1).846 20. superpc is the supervised PCA (superpc package). This method [45]847 retains only a subset of the principal components which are correlated848 to the output. The tunable hyperparameters are the number of princi-849 pal components (n.components), tuned with values 1, 2 and 3 (in all850 53 the datasets the value 1 is selected), and the threshold for retaining851 the input scores, with values 0.1 and 0.9.852 VI. Least absolute shrinkage and selection operator (LASSO)10 853 21. lasso performs LASSO regression, using the enet function in the elas-854 ticnet package with λ= 0 to obtain the LASSO solution.855 22. relaxo develops relaxed LASSO (relaxo package), which generalizes856 the LASSO shrinkage method for linear regression [14]. This method857 is designed to overcome the trade-off between speed and convergence858 in the L2-loss function of the regular LASSO, specially for sparse high-859 dimensional patterns. It provides solutions sparser than LASSO with860 better prediction error. The relaxation hyperparameter (φ) is tuned861 with 7 values from 0.1 to 0.9, while the penalty hyperparameter (λ) is862 tuned with 3 data-dependent values.863 23. lars is the least angle regression (lars package), a model selection864 method [16] which is less greedy than the typical forward selection865 methods. It starts with zero coefficients for all the inputs and finds866 the input imost correlated with the output, increasing step-by-step its867 coefficient until another input jhas high correlation with the current868 residual (i.e., the error, or difference between the true and predicted869 outputs). The coefficients of inputs iand jare increased in the equi-870 angular direction between inputs iand juntil some other input kis so871 correlated with the residual as input j. Then, it proceeds in the equi-872 angular direction among i,jand k, which is the “least angle direction”,873 10Due to the high number of models, LASSO models are included in a specific family and not in the penalized linear regression family. 54 and so on until all the coefficients are non-zero (i.e., all the inputs are in874 the model). The lasso and fraction options are specified for training875 and prediction respectively, and the fraction hyperparameter (ratio876 between the L1 norm of the coefficient vector and the norm at the full877 LS solution) is tuned with 10 values between 0.05 and 1 (for 46.7% of878 datasets the selected value of fraction was 1).879 VII. Ridge regression (RIDGE)11 880 24. ridge develops ridge regression (elasticnet package), which intro-881 duces a regularization term, alongside with the squared difference be-882 tween the desired and true outputs, in the function to optimize. This883 term, which evaluates the model complexity (e.g., the matrix norm for884 linear models), is weighted by the penalty or regularization hyperpa-885 rameter (λ). We use the enet function in the elasticnet package,886 already used for model enet, tuning λwith 5 values between 0.01 and887 0.1 (these two values are selected for 50% and 30% of the datasets,888 respectively).889 25. spikeslab implements the spike and slab regression (spikeslab pack-890 age), which computes weighted generalized ridge regression estimators891 using Bayesian spike and lab models [18]. The spikeslab method892 combines filtering for dimensionality reduction, model averaging using893 Bayesian model averaging, and variable selection using the gnet estima-894 tor. The only tunable hyperparameter is the number of selected inputs895 (vars), with the two values listed by the getModelInfo function: 2896 11Similarly to LASSO, ridge regression models are grouped in a separate family instead of the penalized linear regression family. 55 and the number of inputs (both selected with similar frequencies).897 26. foba is the ridge regression with forward, backward and sparse input898 selection [20], implemented in the foba package. We use the adap-899 tive forward-backward greedy version of the method (with the default900 value foba for the type argument of the foba function), which does a901 backward step when the ridge penalized risk increases in less than the902 parameter ν(with value 0.5 by default) multiplied by the ridge penal-903 ized risk reduction in the previous forward step. The hyperparameters904 are regularization for ridge regression (λ), with 10 values between 10−5 905 and 0.1, and the number of selected inputs or sparsity (k) for the pre-906 diction, with two values: 2 and the number of inputs.907 VIII. Bayesian models (BYM)908 27. bayesglm is the Bayesian GLM, implemented by the arm package. It909 uses expectation maximization to update the βcoefficients of the GLM910 at each iteration, using an augmented regression to represent the prior911 information [22]. The coefficients are calculated using a Student-t prior912 distribution.913 28. brnn is the Bayesian regularized neural network (brnn package), a net-914 work with one hidden layer trained using Gauss–Newton optimization.915 The training minimizes a combination of squared error and a regular-916 ization term which uses the squared network weights [24]. The Bayesian917 regularization [80] determines the weights of both terms based on infer-918 ence techniques. This requires an iterative computation of the Hessian919 matrix (or its Gauss–Newton approximation) of the performance with920 respect to the weights and biases until a goal is met or a maximum921 56 ations or the test p-value (0.01 by default) are increased. The only1072 hyperparameter is mtry, tuned as in rf.1073 47. RRF is the regularized random forest (RRF package), which uses reg-1074 ularization for input selection in rf, penalizing the selection of a new1075 input for splitting when its Gini information gain is similar to the in-1076 puts included in the previous splits [57]. The hyperparameters are1077 mtry, with 3 linearly spaced values between 2 and the number of in-1078 puts, and the regularization coefficient (coefReg), with values 0.01 and1079 1, both selected with similar frequencies.1080 48. cforest is a forest ensemble of conditional inference trees [54], each one1081 fitting one bootstrap sample (party package [58]). The only hyperpa-1082 rameter is the number of selected inputs (mtry, with values 2 and the1083 number of inputs) of the conditional trees.1084 49. qrf is the quantile regression forest (quantregForest package [59]),1085 a tree-based ensemble which generalizes RF in order to estimate con-1086 ditional quantile functions. This regression model grows several RFs,1087 storing all the training patterns associated to each node in each tree.1088 For each test pattern, the weight of each training pattern is the average1089 of the weights of all the training patterns in the leaves activated by that1090 pattern in the different trees of the forest. Using these weights, the dis-1091 tribution function of each output value, and the conditional quantiles,1092 are estimated. The only hyperparameter is mtry (tuned with 2 values1093 as cforest). The quantile prediction threshold (argument what in the1094 predict.quantregForest function) is set to 0.5.1095 50. extraTrees is the ensemble of extremely randomized regression trees1096 63 [60] implemented by the extraTrees package. It randomizes the input1097 and cut-point of each split (or node in the tree), using a parameter1098 which tunes the randomization strength. The full training set is used1099 instead of a bootstrap replica. It is expected that explicit randomiza-1100 tion of input and cut-point splittings combined with ensemble averaging1101 should reduce the variance more than other methods. Its hyperparam-1102 eters are the number of inputs randomly selected at each node (mtry,1103 tuned with 2 values as cforest) and the minimum sample size to split1104 a node (numRandomCuts), tuned with integer values from 1 to 10 (the1105 selected value was 1 for 48.3% of the datasets).1106 XV. Bagging ensembles (BAG)1107 51. bag [62] is the bagging ensemble of conditional inference regression1108 trees (see model #39) implemented by the caret package. The output1109 for a test pattern is the average of the outputs over the base regression1110 trees.1111 52. bagEarth is the bagged MARS (caret package), a bagging ensemble1112 of MARS base regression models implemented in the earth package1113 (see model #9), which learns a MARS model with degree=1 for each1114 bootstrap sample. The only hyperparameter is the maximum number1115 of terms (nprune) in the pruned regression model (10 values).1116 53. treebag is the bagged CART, a bagging ensemble of rpart regression1117 base trees (see model #37), implemented by the ipredbagg function1118 in the ipred package [64].1119 XVI. Boosting ensembles (BST)1120 54. rndGLM is a boosting ensemble of GLMs [65] implemented by the1121 64 randomGLM package (also named randomGLM in the caret model list).1122 It uses several bootstrap samples (nBags=100 by default) of the train-1123 ing set, randomly selecting inputs and interaction terms among them1124 depending on the maxInteractionOrder hyperparameter, tuned with1125 values 1, 2 and 3 (selected with frequencies 53.3%, 40% and 6.7%, re-1126 spectively). For each sample, inputs are ranked by its correlation with1127 the output, and a predefined number of them are selected, using forward1128 selection, to create a multivariate GLM. For a test pattern, the pre-1129 dicted value is the average of the GLM outputs. This regression model1130 has very high memory requirements, overcoming the largest available1131 memory (128GB) in 38 datasets, and requiring 128, 64, 32 and 16GB1132 in 2, 10, 22 and 12 datasets, respectively.1133 55. BstLm is the gradient boosting machine with linear base models, im-1134 plemented in the bst package. Gradient boosting optimizes arbitrary1135 differentiable loss functions defining the fitting criteria [53]. Boosting1136 combines weak base regression models into a strong ensemble by it-1137 eratively adding base models, and in each iteration the new model is1138 trained to fit the error (residual) of the previous ensemble. Since the er-1139 ror can be viewed as the negative gradient of the squared error loss func-1140 tion, boosting can be considered a gradient descent method. BstLm1141 uses the bst function with linear base models (argument learner=lm)1142 and Gaussian family, since squared error loss is used. The only hyper-1143 parameter is the number of boosting iterations (mstop), with 10 values1144 from 50 to 500.1145 56. bstSm is the gradient boosting with smoothing spline base regression1146 65 models (learner=sm in the bst function of the same package). The1147 number of boosting iterations (mstop) is tuned with 10 values as BstLm.1148 57. bstTree is the gradient boosting with regression base trees (learner=1149 tree, same function and package as BstLm). The hyperparameters are1150 the number of boosting iterations (mstop, 4 values from 40 to 200) and1151 the maximum depth of nodes in the final tree (maxdepth item in the list1152 of the control.tree argument of the bst function), with integer values1153 between 1 and 5 (this last value is selected in 55% of the datasets).1154 58. glmboost is the gradient boosting ensemble with GLM base regression1155 models (glmboost function in the mboost package), tuning the number1156 of boosting iterations (mstop, 10 values).1157 59. gamboost is the boosted generalized additive model (mboost package),1158 a gradient boosting ensemble of GAM base regression models [69]. The1159 ensemble minimizes a weighted sum of the loss function evaluated at1160 the training patterns by computing its negative gradient. The base re-1161 gression models are component-wise models (P-splines with a B-spline1162 base, by default). The only hyperparameter is the number of initial1163 boosting iterations (mstop), with 10 values from 50 to 500, selecting1164 500 as the best value for 56.7% of the datasets.1165 60. gbm is the generalized boosting regression model (gbm package [49]),1166 named stochastic gradient boosting in the caret model list. The hyper-1167 parameters are the maximum depth of input interactions (interaction.1168 depth), with integer values from 1 to 5 (the last value was selected in1169 48.3% of the datasets), and number of trees for prediction (n.trees),1170 with values from 50 to 250 with step 50 (the value 250 was selected in1171 66 45% of the datasets). We use a Gaussian distribution and shrinkage=1172 0.1 (default values).1173 61. blackboost is the gradient boosting (blackboost function in the mboost1174 package) with conditional inference regression base trees (ctree in the1175 party package, see model #40) and arbitrary loss functions [51]. The1176 only hyperparameter is the maximum tree depth (maxdepth argument1177 in the party::ctree control function, used as tree controls argu-1178 ment of the blackboost function), with integer values from 1 to 5,1179 value selected in 79% of the datasets.1180 62. xgbTree is the extreme gradient boosting [53], using the xgb.train1181 function in the xgboost package with booster=gbtree, root mean1182 squared error as evaluation metric and linear regression as objective1183 function. The hyperparameters are the maximum tree depth (max depth),1184 with values 1, 2 and 3 (max depth=3 for 53.3% of the datasets); maxi-1185 mum number of boosting iterations (nrounds), with values 50, 100 and1186 150; and learning rate (η), with values 0.3 and 0.4.1187 63. xgbLinear is the extreme gradient boosting with booster=gblinear1188 and linear regression as objective function (xgboost package). Its hy-1189 perparameters are the L2 (square loss) regularization term on weights1190 (λ, with values 0, 0.1 and 0.0001), bias (α, with values 0 and 0.1), and1191 number of boosting iterations (nrounds, tuned as xgbTree).1192 XVII. Neural networks (NET)1193 64. mlpWD is the classical multi-layer perceptron with one hidden layer1194 and weight decay (named mlpWeightDecay in the caret model list).1195 It uses the mlp function in the RSNNS package, with learnFunc =1196 67 BackpropWeightDecay. The tunable hyperparameters are the size1197 of the hidden layers (5 odd values between 1 and 5) and the weight1198 decay (values 0, 0.1, 0.042, 0.01778 and 0.007498).1199 65. mlpWDml is the same network with three hidden layers (RSNNS pack-1200 age, named mlpWeightDecayML in the caret model list), tuning four1201 hyperparameters: the sizes of the three hidden layers (layer1,layer21202 and layer3, tuned with values 1, 3 and 5 each one) and the weight1203 decay (same values as mlpWD).1204 66. avNNet is the model averaged neural network (caret package). A1205 committee of 5 (argument repeats) multi-layer perceptron neural net-1206 works of the same size trained using different random seeds, being av-1207 eraged to give an output [81]. The boolean argument linout is set to1208 have linear output neurons for regression, and MaxNWts is adjusted to1209 allow the number of weights required by the dataset inputs. The hy-1210 perparameters are the network size, tuned with 7 odd values between1211 1 and 13, and the weight decay (with values 0, 0.1 and 0.0001).1212 67. rbf is the radial basis function network (RSNNS package) which does1213 a linear combination of basis functions, each centered around a pro-1214 totype [56]. The information is locally codified (opposed to globally1215 in the MLP), the training should be faster and the network is more1216 interpretable, although the output might be undefined if a test pattern1217 does not activate any prototype. The only hyperparameter is the size1218 of the hidden layer (10 odd values from 1 to 19).1219 68. grnn is the generalized regression neural network [61], a special type1220 of RBF network implemented by the Matlab neural network toolbox.1221 68 After a clustering of the training set, the nodes of the hidden layer store1222 the cluster centers, although the Matlab implementation uses so many1223 clusters as training patterns. The output for a test pattern is a weighted1224 sum of the Gaussian functions centered in the cluster centers, scaled1225 by the cluster populations. During training, whenever a pattern is1226 assigned to a cluster, the weight of the Gaussian function corresponding1227 to that cluster is updated using the desired output. The Gaussian1228 spread is the only hyperparameter (13 values between 0.01 and 2):1229 large (resp. small) values lead to smooth (resp. close) approximations.1230 69. elm is the extreme learning machine [63] implemented by the elmNN1231 package. The only hyperparameters are the number of hidden neurons1232 (nhid), with 40 odd values between 1 and 79 (the last value was se-1233 lected in 11.7% of the datasets), and the activation function (actfun),1234 with 4 values: sin,radbas,purelin and tansig, selected with similar1235 frequencies.1236 70. kelm is the ELM neural network with Gaussian kernel [63] using the1237 publicly available Matlab code12. The hyperparameters are regulariza-1238 tion Cand kernel spread σ, tuned with values {2i}14 −5and {2i}8 −16, with1239 20 and 25 values, respectively.1240 71. pcaNNet is a multi-layer perceptron neural network with one hidden1241 layer trained on the PCA-mapped training patterns, implemented by1242 the caret and nnet packages. The principal components which account1243 for more than 95% of the data variance are used for training. Each test1244 pattern is mapped to the principal component space and the trained1245 12http://www.ntu.edu.sg/home/egbhuang/elm kernel.html (visited March 29, 2017). 69 pcaNNet model gives an output. The tunable hyperparameters are the1246 size of the hidden layer, with 7 values between 1 and 13, and the1247 weight decay of the network, with values 0, 0.1 and 0.0001.1248 72. bdk is the supervised bi-directional Kohonen network, implemented1249 in the kohonen package [66]. The bdk combines Kohonen maps and1250 counterpropagation networks using two maps, for inputs and output1251 respectively. In each iteration, the direct (resp. inverse) pass updates1252 only the weights of the input (resp. output) map, using a weighted sim-1253 ilarity measurement (Euclidean distance for regression) which involves1254 both maps, leading to a bi-directional updating. The test output is the1255 weight of the winner node of the output map. The hyperparameters1256 are the sizes of both maps (xdim and ydim, with 3 values from 3 to1257 17) and the initial weight (xweight) given to the input map during the1258 distance calculation for the output map, and to the output map for1259 updating the input map, tuned with values 0.5, 0.75 and 0.9).1260 XVIII. Deep learning (DL)1261 73. dlkeras is the deep learning neural network implemented by the Keras1262 module [67] of the Python programming language, with three hidden1263 layers tuned with 50 and 75 neurons for each layer (nh1,nh2 and nh3,1264 with 8 combinations). The deep learning methods [82, 83] are very1265 popular, specially for image classification, and they are included in this1266 comparison for regression tasks.1267 74. dnn is the deep belief network implemented in R by the DeepNet pack-1268 age [68]. It uses three hidden layers, tuning their number of neurons1269 using 3 values for each layer (27 combinations). The weights are ini-1270 70 tialized using stacked autoencoder (SAE), which in our experiments1271 gave better results than deep belief network (DBN). Hidden and out-1272 put neurons have hyperbolic tangent and linear activation functions,1273 respectively.1274 XIX. Support vector regression (SVR)1275 75. svr is the ε-support vector regression with Gaussian, accessed via the1276 C++ interface of the LibSVM library [9]. We tuned the regularization1277 hyperparameter Cand the kernel spread γwith values {2i}14 −5and1278 {2i}8 −16, with 20 and 25 values, respectively.1279 76. svmRad is another implementation of SVR (named svmRadial in the1280 caret model list) with Gaussian kernel, which uses the (ksvm function1281 in the kernlab package [70] for regression (argument type=eps-svr).1282 This implementation also uses LibSVM, and it tunes the regularization1283 hyperparameter C, with 20 values in the set {2i}15 −4, and the kernel1284 spread σ. Although we specify 25 values for σ, the getModelInfo1285 function only lists 6 values in the set {2−i}7 5.1286 77. rvmRad is the relevance vector machine [71] with Gaussian kernel1287 (kernlab package), named rvmRadial in the caret model list. The1288 RVM has the same functional form as the SVM, but it uses a Bayesian1289 learning framework which reduces the number of basis functions, com-1290 pared to the SVM, while keeping an accurate prediction. This regres-1291 sion model avoids the tunable regularization hyperparameter (C) of1292 the SVM, but it uses a method similar to expectation-maximization1293 which, unlike SMO, may fall in local minima. The value of the Gaus-1294 sian spread σis estimated by the getModelInfo function, which only1295 71 lists one value.1296 References1297 [1] M. Fern´andez-Delgado, E. Cernadas, S. Barro, D. Amorim, Do we need1298 hundreds of classifiers to solve real classification problems?, J. Mach.1299 Learn. Res. 15 (2014) 3133–3181.1300 [2] D. H. Wolpert, The lack of a priori distinctions between learning algo-1301 rithms, Neural Computation 9 (1996) 1341–1390.1302 [3] K. Bache, M. Lichman, UCI machine learning repository (2013).1303 URL http://archive.ics.uci.edu/ml1304 [4] G. E. P. Box, D. R. Cox, An analysis of transformations, J. Royal Stat.1305 Soc. Series B (Methodological) 26 (2) (1964) 211–252.1306 [5] R Team, R: A language and environment for statistical computing, Vi-1307 enna, Austria, ISBN 3-900051-07-0 (2008).1308 URL https://www.R-project.org1309 [6] M. Kuhn, Caret: classification and regression training, R package1310 (2016).1311 URL http://topepo.github.io/caret/train-models-by-tag.html1312 [7] Theano Team, Theano: A Python framework for fast computation of1313 mathematical expressions, arXiv e-prints.1314 [8] Python Software Foundation, Python Language (2017).1315 URL https://www.python.org1316 72 [67] F. Chollet, Keras: The Python Deep Learning library (2015).1444 URL https://keras.io1445 [68] X. Ron, Deepnet package (2015).1446 URL https://cran.r-project.org/package=deepnet1447 [69] P. Buehlmann, B. Yu, Boosting with the L2 loss: regression and classi-1448 fication, J. Am. Stat. Assoc. 98 (2003) 324–339.1449 [70] A. Karatzoglou, Kernlab package (2015).1450 URL https://cran.r-project.org/package=kernlab1451 [71] M. Tipping, Sparse Bayesian learning and the relevance vector machine,1452 J. Mach. Learn. Res. 1 (2001) 211–244.1453 [72] T. Hothorn, F. Leish, A. Zeileis, K. Hornik, The design and analysis1454 of benchmark experiments, J. Comput. and Graph. Stat. 13 (3) (2005)1455 675–699.1456 [73] S. Garc´ıa, A. Fern´andez, A. Ben´ıtez, F. Herrera, Statistical comparisons1457 by means of non-parametric tests: a case study on genetic based machine1458 learning, in: Proc. II Congreso Espa˜nol de Inform´atica (CEDI 2007),1459 2007, pp. 95–104.1460 [74] M. Hollander, D. Wolfe, Nonparametric Statistical Methods, John Wiley1461 & Sons, 1973.1462 [75] J. Demˇsar, Statistical comparisons of classifiers over multiple data sets,1463 J. Mach. Learn. Res. 7 (2006) 1–30.1464 79 [76] C. Dunnett, A multiple comparison procedure for comparing several1465 treatments with a control, J. Am. Stat. Assoc. 50 (1955) 1096–1121.1466 [77] J. D. Gibbons, S. Chakraborti, Nonparametric Statistical Inference,1467 CRC Press, 2011.1468 [78] T. Pohlert, The pairwise multiple comparison of mean ranks package1469 (PMCMR), r package (2014).1470 URL http://CRAN.R-project.org/package=PMCMR1471 [79] T. Colton, Statistical in medicine, Little Brown and Co., New York, NJ,1472 1974.1473 [80] D. MacKay, Bayesian interpolation, Neural Computation 4 (1992) 415–1474 447.1475 [81] B. Ripley, Pattern recognition and neural networks, Cambridge Univ.1476 Press, 1996.1477 [82] G. E. Hinton, S. Osindero, Y.-W. Teh, A fast learning algorithm for1478 deep belief nets, Neural Comput. 18 (7) (2006) 1527–1554.1479 [83] W. Liu, Z. Wang, X. Liu, N. Zeng, Y. Liu, F. Alsaadi, A survey of deep1480 neural network architectures and their applications, Neurocomputing1481 234 (2017) 11–26.1482 80