Application of machine learning to agricultural soil data

UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Centro de Investigación en Tecnoloxías da Información (CiTIUS) Tesis doctoral APPLICATION OF MACHINE LEARNING TO AGRICULTURAL SOIL DATA Presentada por: Manisha Sanjay Sirsat Dirigida por: Eva Cernadas García Manuel Fernández Delgado Santiago de Compostela, Julio de 2017 Dr. Eva Cernadas García, Profesora Titular de Universidad del Área de Ciencias de la Computación e Inteligencia Artificial de la Universidad de Santiago de Compostela Dr. Manuel Fernández Delgado, Profesor Titular de Universidad del Área de Ciencias de la Computación e Inteligencia Artificial de la Universidad de Santiago de Compostela HACEN CONSTAR: Que la memoria titulada APPLICATION OF MACHINE LEARNING TO AGRICULTURAL SOIL DATA ha sido realizada por Manisha Sanjay Sirsat bajo nuestra dirección en el Centro Singular de Investigación en Tecnoloxías da Información de la Universidad de Santiago de Compostela, y constituye la Tesis que presenta para obtar al título de Doctora. Santiago de Compostela, Julio de 2017 Eva Cernadas García Directora de la tesis Manuel Fernández Delgado Co-director de la tesis Manisha Sanjay Sirsat Autora de la tesis Education is the manifestation of the perfection already in man. Swami Vivekananda Acknowledgments I would like to express my gratitude and respect for my thesis advisors, Dr. Eva Cernadas García and Dr. Manuel Fernández Delgado for their guidance and support which contributed greatly to the timely completion of this project. I also extend my gratitude to Centro de Investigación en Tecnoloxías da Información (CiTIUS), University of Santiago de Compostela. My gratitude also goes to Dr. Razaullah Khan and Dr. Balaji Aglave for their advise and motivation. I deeply thank to my family Sakharam, Antika and Sanjay for their heart-warming kindness, timely encouragement, and endless patience. Finally, I am thankful to Erasmus Mundus Euphrates programme [project number 2013-2540/001-001-EMA2] for giving me wonderful opportunity to pursue my research work internationally. Santiago de Compostela, July 2017 Resumen El sector agrícola constituye uno de los principales pilares de la economía en la India, representando en torno al 47% del empleo (datos de 2013–2014) y el 11% del producto interior bruto de este país, que dedica el 60% de su territorio a labores agrícolas y ganaderas. El funcionamiento de este sector económico está condicionado por factores climáticos, así como por las prácticas agrícolas y las características del suelo. De hecho, la sobre-explotación y el uso excesivo, o incorrecto, de fertilizantes afecta a un porcentaje significativo del suelo agrícola indio, dos tercios del cual se puede calificar como degradado. La medición de parámetros del suelo tales como los índices de fertilidad de las diferentes localidades para distintos nutrientes, tipo de suelo, acidez y niveles de nutrientes, entre otros, resulta fundamental de cara a planificar el uso de fertilizantes y la selección de cultivos, pero requiere un considerable esfuerzo por parte de personal especializado en análisis químicos. Esta tesis evalúa la capacidad de los métodos de aprendizaje automático (machine learning) para la predicción de estos parámetros a partir de datos químicos del suelo publicados por los gobiernos de la India y del estado de Maharashtra. Estos datos permiten la predicción de los siguientes 12 parámetros del suelo de interés para la agricultura: 1-6. Índices locales de fertilidad para 6 nutrientes del suelo: carbono orgánico (OC), pentóxido de fósforo (P2O5), óxido de potasio (K2O), hierro (Fe), manganeso (Mn) y cinc (Zn). Para un nutriente y para cada terreno agrícola de la localidad (que se corresponde con un patrón para las tareas de aprendizaje automático), se cuantifica el nivel del nutriente, que es un dato de tipo numérico, en bajo, medio o alto, usando para ello los niveles definidos por el gobierno de la India. El índice local de fertilidad se calcula promediando de forma ponderada el número de terrenos de la localidad con niveles bajos, medios y altos (Nl,NmyNh) con pesos 1, 2 y 3 respectivamente, mediante la fórmula (Nl+2Nm+3Nh)/Nt, sendo Ntel número total de terrenos en la localidad. Este índice viii rápidas de la colección. Es el caso de extraTrees, cubist, random forest y generalized boosted regression model (gbm), situadas en las 10 primeras posiciones en el ranking de correlación y en posiciones medias del ranking de tiempos (25, 48, 51 y 37 respectivamente). Esta colección de regresores se ha empleado para predecir los valores numéricos de los parámetros del suelo descritos anteriormente. En este caso, al tratarse de métodos de regresión, se han suprimido el cultivo recomendado y tipo de suelo, ya que al poseer salidas discretas son intrínsecamente problemas de clasificación. Sin embargo, entre los problemas de regresión se han incluido los índices locales de fertilidad para K2OyZn, que no se consideraron en el capítulo 2 porque los datos disponibles sólo pertenecen a uno de los niveles de cuantificación (clases) definidos por el gobierno de la India. La metodología experimental y medidas de calidad son las mismas que en el capítulo 3. A diferencia de los experimentos con los datos de la UCI, en los datos de suelo el regresor extraTrees, que ya obtenía buenos resultados de correlación y velocidad, alcanza los mejores resultados en 7 de 10 problemas, y en los restantes funciona casi tan bien como el mejor regresor, superando el 90% de la correlación máxima. De hecho, el ranking de Friedman de los regresores sobre todos los problemas sitúa a extraTrees en primera posición con un valor de 2.3, lo cual significa que en promedio este regresor figura entre las posiciones 2 y 3 de la colección de 76 regresores. Sin embargo, en términos absolutos la mejor correlación promedio sobre todos los problemas, obtenida por extraTrees, es relativamente reducida (0.68), lo cual refleja la complejidad del problema de predicción abordado. Los 5 mejores regresores en términos de correlación pertenecen a la familia random forest: extraTrees, regularized random forest, random forest, quantile regression forest y Boruta. La svr y dos regresores de la familia gradient boosted machines, generalized boosted regression model (gbm) y gradient boosting with regression trees (bstTree) ocupan las siguientes posiciones. Los resultados son peores para elm-kernel, que se sitúa en las posiciones 10 y 14 en términos de RMSE y correlación respectivamente, aunque su correlación media (0.62) no es muy inferior a la obtenida por extraTrees. Considerando conjuntamente la correlación y la velocidad de los regresores, extraTrees proporciona el mejor compromiso entre ambos factores, ya que alcanza el mínimo ranking de Friedman de correlaciones, y el segundo menor ranking de Friedman de tiempos entre los 20 regresores más precisos. En un análisis individualizado para cada uno de los 10 problemas de regresión, se aprecia una cierta correspondencia con los resultados de las técnicas de clasificación (capítulo 2). En la clasificación de índices locales de fertilidad de OC yFe, donde los mejores κ alcanzan ix 90.65% y 69.35% respectivamente, las mejores correlaciones obtenidas en la regresión son superiores a 0.8. También en el índice local de Zn la correlación supera 0.8, pero este índice no se considera en la clasificación, por pertenecer todos los patrones disponibles a la misma clase. En los índices de fertilidad de P2O5yMn, y en el pH, donde κ alcanza 85.54%, 64.8% y 47.32% respectivamente, las correlaciones son intermedias (0.776, 0.758 y 0.694 respectivamente). El índice local de K2Otambién es intermedio (0.631), pero este índice tampoco se considera en la clasificación, por la misma razón que el índice de Zn. Y, finalmente, para los nutrientes N2O,P2O5yK2O, donde en la clasificación κ alcanza valores de 33.6%, 35.08% y 31.85% respectivamente, en la regresión se alcanzan las correlaciones más bajas (0.519, 0.487 y 0.517, respectivamente). Estos resultados confirman la dificultad de los problemas analizados, ya que las mejores correlaciones no alcanzan 0.9, y con los tres nutrientes N2O,P2O5y K2Ono superan 0.6. Se puede decir que la predicción es aceptable sólo en los índices locales de fertilidad de OC,Fe yZn, donde la correlación supera 0.8. En el resto de los problemas, la incertidumbre en la predicción es más elevada. Sin embargo, si no resulta imprescindible una aproximación precisa de los parámetros del suelo, las técnicas empleadas proporcionan una predicción aproximada para los valores numéricos mediante técnicas de regresión, y para los valores cuantificados mediante técnicas de clasificación. Contents 1 Introduction 1 1.1 Soil fertility index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Fertilizer nutrients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Soil reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Cropping cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Soil texture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.6 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.7 Outline of the work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Application of classification methods to agricultural soil data 9 2.1 Village-wise OC,P2O5,Mn and Fe fertility indices . . . . . . . . . . . . . . 9 2.2 Soil nutrients N2O,P2O5and K2O....................... 10 2.3 Soil pH ..................................... 12 2.4 Crop selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Soil type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.6 Classification methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.7 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.8 Global discussion of the results . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.9 Classification of village-wise OC,P2O5,Mn and Fe fertility indices . . . . . 26 2.10 Classification of soil nutrients N2O,P2O5and K2O.............. 27 2.11 Classification of soil pH ............................ 27 2.12 Classification of crop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.13 Classification of soil type . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.14 Comparison among regions . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 xii Contents 3 Application of regression methods to general datasets 31 3.1 Regression methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4 Application of regression methods to agricultural soil data 65 4.1 Regression problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Prediction of OC,P2O5,K2O,Fe,Mn and Zn village-wise fertility indices . . 68 4.3 Prediction of soil nutrients N2O,P2O5and K2O. . . . . . . . . . . . . . . . 74 4.4 Prediction of soil pH .............................. 76 4.5 Global discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5 Conclusions 85 Bibliography 89 List of Figures 101 List of Tables 103 CHAPTER 1 INTRODUCTION India is second largest country in the farm production and seventh largest in agricultural export [30]. Agriculture is the pillar of the Indian economy, offering 47% employment of total population in 2013-14 and having significant contribution in global food basket. India is the fastest growing exporter of agricultural products over a decade to more than 100 countries, mainly in the Middle East, Southeast Asia, countries of the South Asian Association for Regional Cooperation (SAARC), the European Union and the United States [128]. According to data of year 2011, India devotes 60.5% of its land1to agriculture, distributed among arable land (52.8%), land for permanent crops (4.2%) and pastures (3.5%). Share of agriculture and related activities was 11.3% of the Gross State Domestic product (GSDP) in 2013-14. The 11th five-year economic plan acknowledges the need of proper soil management in agriculture. Excessive and miscalculated use of fertilizer focused on increase production has led to soil degradation, and today nearly 66.67% of India’s agricultural land can be categorized as either degraded or sick [89]. In India, each state and union territory is responsible for the set-up of the soil testing facilities and maintaining the state soil database. Some states, for e.g. Gujarat initiated the “Soil Health Cards Programme” in 2006 to recommend fertilizers, crop rotations and to record the data on a national network which can be used to survey different soil types2. The farmers are assisted by several non-governmental organization and community groups, who are responsible for soil sample testing [89]. In year 2013-14, the cultivation areas of major crops were 15 and 57 millions of hector in Kharif and Rabi seasons, respectively [24]. However, agriculture in India is conditioned by the poor fertility of the soil, 1https://www.cia.gov/library/publications/the-world-factbook/geos/in.html 2http://indiagovernance.gov.in/news.php?id=204 2 Chapter 1. Introduction which depends on the levels of its nutrients. The physical, chemical and biological properties of the soil are useful to evaluate its fertility, to design a cultivation plan and to predict the crop productivity. The information technologies, and specifically Machine Learning (ML), offer new possibilities in the field of agriculture and may help in data evaluation for decision making. Maharashtra ranks second largest state in terms of population and geographical area in India and it is located at 15o38” to 22o01” North and 72o39” to 80o44” East where, total 64.14% of the people are employed in agriculture and allied activities. Agricultural calendar of Maharashtra is governed by monsoon, being the 60% of the farming rain-fed. This state is divided in 9 agro climatic zones, composed by 39% of shallow soils and 42.4% of degraded land. Maharashtra has 5 main regions; Vidarbha, Konkan, Marathwada, Paschim Maharashtra, and North Maharashtra, although only the three latter will be considered in the current study. Marathwada, is one of the most prominent agricultural regions in India, located at 19o 52’ 59.88” North and 75o19’ 59.88” East. Figure 1.1: Geographical representation of 6 regions of Maharashtra (India). Marathwada, Paschim Maharashtra and North Maharashtra are study areas, highlighted by red borders. The soil of Marathwada is made of basalt rock with scarlet, blackish and yellowish colors and semi-dry plateau, good in iron level, moisture retentive, and poor in nitrogen as well as organic matter, including variable climatic condition. However, temperature is mostly humid throughout the year, maximum and minimum temperature oscillates from 27◦C to 40◦C and 14◦C to 27◦C. The classification of soil according to its physical and chemical properties is 1.1. Soil fertility index 3 useful to maintain and enhance its productivity, to avoid soil degradation problems and to overcome environmental damage. The major challenge is to increase crop yield for solving global food security problem. However, soil quality and crop yield are negatively affected by changing trends of temperature and rainfall, insufficient water and light, agriculture practices and absence of nutrients. It is important to develop an effective nutrient management by means of an adequate soil analysis and a proper application of fertilizers. Hence the relevance of a research effort to classify soil parameters such as the fertility indices for several nutrients: organic carbon (OC), phosphorus pentoxide (P2O5), manganese (Mn) and iron (Fe), among others; soil pH, soil type, preferred crop and levels of nutrients as the nitrous oxide (N2O), phosphorus pentoxide (P2O5) and potassium oxide (K2O), which are relevant for fertilizer recommendation. The interest of predicting the levels of these magnitudes with ML techniques is to avoid the need to chemically measure these magnitudes, thus reducing the cost of the analysis and saving time of specialized technicians. The current study tries to enhance the accuracy of soil problem interpretation for Indian agriculture, although similar studies would benefit other nations around the globe. 1.1 Soil fertility index In India, land-man ratio is quickly decreasing, so there is a need increase agricultural production without harm to environment and sustainability. The per capita land in India has decreased from 0.48 to 0.41 ha./person since year 1951 to year 2001, and the prediction is to decrease even more, until 0.10 ha/person by the year 2025. Besides, urbanization and industrialization led to destruction of forest, and to a reduction in the cultivation land and in its quality. As a result, agricultural production is being affected unfavorably. This situation appeals for planning of soil fertility by supplying essential nutrients to the crop in sufficient amount and at right time for its best growth. Therefore, fertilizers are a great significance input for targeting high crop production. Moreover, imbalances in soil quality leads to crop health and lower/higher crop yield [87]. The Indian soil has generally low or medium levels of soil nutrients and organic matters. Consequently, an efficient soil nutrient management is a major concern for maintaining soil fertility and it mainly depends on optimum fertilizer rates as per nutrient demand of crop. The crop needs optimum quantity of various nutrients for its vegetative growth and ultimate yielding. Declining status of soil fertility and mismanagement of soil nutrients may be factors for food crises for the world’s population [44]. Generally, 4 Chapter 1. Introduction Indian soil fertility data are summarized for block and district level. These data are useful for decision making of application of suitable amount fertilizers, policy of fertilizer distribution and consumption in the view of changes in fertility levels. One of the objectives of the work developed in the current PhD. Thesis is to predict village-wise soil fertility indices, which can be used in preparing a village-wise soil fertility index map. This map would allow to compare levels of soil fertility among villages, and to make fertilizer recommendations. 1.2 Fertilizer nutrients Maharashtra is one of the major fertilizer-consuming states in India with annually around 1.8 to 2.0 million tonnes in terms of N2O,P2O5and K2Osoil nutrients3. Distorted levels of these nutrients might increase fertilizer demand and this can adversely affect soil fertility. The usage of these fertilizers in Maharashtra has rapidly increased from 56.85, 29.74, 16.14 4to 65.93, 32.13, 17.64 5(kg/ha) since 2012-13 to 2013-14. The soil Nis crucial nutrient fertilizer for food production moreover, responsible parameters for it’s availability in soil are its type, texture, soil pH, climate, etc. Generally, soils contain 0.02 to 0.44 % total N however, climate of study area is tropical so that soils are poor in OC, and consequently low in N. The phosphorus (P) plays a key role in substances are used as building blocks for genes and chromosomes, in plant root growth and biochemical processes that involve energy transfer. Generally, the Pcontent in agricultural crop ranges from 0.1 to 0.5%. The status of Pin Indian soils range between 134.6 to 310.4 kg/ha. There is significant relation between OC and Pdue to creation of sustainable soil environment with soil organic matter [121]. However, deficiency of Pleads to: breakdown of plant cell membranes, which reduces energy transfer; lowering of root ratio; poor seed/fruit setting; decreased disease resistance, and reduced tillering in cereals. 1.3 Soil reaction The soil reaction is another name for the pHwhich measures the acidity or alkalinity in the soil. The soil pH influences the solubility of nutrients, affecting the activity of microorganisms responsible for breaking down organic matter and most chemical transformations in the 3http://www.moef.nic.in/soer/state/SoE%20report%20of%20Maharashtra.pdf 4http://fert.nic.in/sites/default/files/Indian%20Fertilizer%20SCENARIO-2014.pdf 5http://fert.nic.in/sites/default/files/Indian%20Fertilizer%20SCENARIO-2014_0.pdf 1.4. Cropping cycle 5 soil, affects how plants grow, and modifies the availability of several plant nutrients. The pH change over time is influenced by factors including chemical composition of the soil, weathering and current agricultural practices, and it also fluctuates through the year. When soil acidity changes, the solubility of metal ions also changes, and the plant growth is really affected by the varying concentration of these metals in solution rather than by the acidity itself. The desirable pH range for optimum plant growth varies among crops. While some crops grow best between 6.0 and 7.0, others grow well under slightly acidic conditions. Soil properties that influence the need for and response to lime vary by region. Soils become acidic when basic elements such as calcium, magnesium, sodium and potassium held by soil colloids are replaced by hydrogen ions. Soils formed under conditions of high annual rainfall are more acidic than soils formed under more arid conditions. In India, most Southeastern soils are inherently more acidic than soils of the Midwest and far West. Soils formed under low rainfall conditions tend to be basic with soil pH readings around 7.0. Intensive farming over a number of years with nitrogen (N) fertilizers or manures can result in soil acidification. In the wheat-growing regions of Kansas and Oklahoma (USA), for example, which have soil pH of 5.0 and below, aluminum toxicity in wheat and good response to liming have been documented in recent years. A knowledge of the soil and the crop is important in managing soil pH for the best crop performance. The pH levels generally find in soil are listed in the figure 2.1 of the subsection 2.3. Significance of soil pH classification is to increase the soil quality and crop production. 1.4 Cropping cycle In Maharashtra, uncertain climatic situations, namely rainfall, affects on agricultural productivity that tends to certain type of farming only. This state is highly dependent on monsoon cycle for large crop yields. The major crop produced by Maharashtra is cotton and soybeans. Moreover, farmers generally use a single cropping due to drought conditions, climate change, soil fertility status and unawareness about soil nutrient status. However, crop rotation is an effective technique to control of weeds, pests, diseases, and more economical utilization of soil fertility. An accurate crop classification allows to predict the optimal crop based on available soil nutrients, whereas significance of optimal cropping system is to avoid environmental damage, consequently it improves soil quality, reduces the build-up of pests, spreads 12 Chapter 2. Application of classification methods to agricultural soil data are used to recommend the suitable amounts of fertilizers in the following way. Let C(X) be the predicted class for nutrient X, defined as C(X)= 1, 2 or 3 for classes low, medium and high, respectively, being the nutrient X=N2O,X=P2O5or X=K2O. Let also R(Y)be the amount, in kg/ha, of fertilizer Y, being Y=uria, Y=super phosphate or Y=muriate of potash the recommended fertilizers to correct the level of nutrients N2O,P2O5or K2O, respectively. The Mahatma-Phule Agricultural University [71] recommends an amount R(X,Y)of fertilizer Yto correct the level of nutrient Xcalculated using C(X)and a pre-defined reference amount F(Y)of fertilizer Y, according to an expression proposed by the technicians of the State Government of Maharashtra (see footnote 1): R(X,Y) = 6−C(X) 4F(Y)(2.1) Therefore, R(1,Y) = 5F(Y)/4, R(2,Y) = F(Y)and R(3,Y) = 3F(Y)/4 for C(X)=1, 2 and 3, respectively, so the recommended amount is 125%, 100% or 75% of F(Y)when the level of nutrient Xis low, medium and high, respectively. This expression allows to calculate the recommended amount of fertilizer Yusing the predicted class C(X)for nutrient X, given the reference amount F(Y). 2.3 Soil pH The pH is the scale of soil acidity or alkalinity, which affects the crop yield and all the soil parameters, because soil acidity is one of its major degradation problems. Specifically, the region of Marathwada has slightly alkaline soil (high pH), which leads to nutrient deficiency, low OC and high CaCO3levels, limits the crop growth and reduces the crop yield. The pH classification uses the inputs listed in the line labeled “pH” of Table 2.2. Although usually nine pH levels are considered (Figure 2.1), but for our purposes it is enough to discriminate between the four middle classes: slightly acidic (labeled SA), neutral (N), slightly alkaline (SAL) and moderately alkaline (MAL). The classification of pH into levels is useful to decide suitable crops and pesticides, and to evaluate microbial activity, nutrient levels and soil corrosion. 2.4 Crop selection The growth of a given crop needs a balanced supply of important nutrients, whose levels define the best crop for a given soil. The cropping cycle determines the way in which the soil 2.5. Soil type 13 Figure 2.1: Degree of acidity and alkalinity of soil [16], with nine classes (up) and with the four classes considered in this paper (down). parameters are enhanced: a good cycle is important for optimum yield and improvement of soil quality, erosion and moisture, organic matter and carbon storage [59]. There are several studies about crop classification [127], crop yield forecasting [81], mapping for crop rotation [93] and land cover classification on remotely sensed data using time-series analysis techniques [113]. The RF, linear discriminant analysis (LDA) and SVM were used to classify crops (walnut, table grape, almond and European plum) on four feature sets [96]. The SVM was also used for plant discrimination using satellite land images [45]. Our data includes four crops, which require neutral pH (in the range 6.5-7.5): bajra(R), cotton(I), cotton(R) and soybean(R), where R and I mean rainfed and irrigated, respectively. The Figure 2.3 (left panel) shows a geographical plot of the crop data, which are distributed across a wide area including the Aurangabad, Jalna, Bid, Osmanabad and Latur districts of the Marathwada region, inside of the state of Maharashtra. In the current paper, crop classification predicts which crop is suitable for the next stage of the cropping cycle using the inputs listed in the line labeled “Crop” of Table 2.2. 2.5 Soil type The soil classification according to its type allows to select the best soil for a particular crop. Soil has been classified [126] using LR, ANN, SVM, KNN, RF and DT in 5 types: 1) coarse loamy, mixed, mesic, lithic xerorthents; 2) fine, mixed, mesic, typic talcixerepts; 3) fine loamy, carbonatic, mesic, typic calcixerepts; 4) fine loamy, mixed, mesic, typic haploxerepts; and 5) 14 Chapter 2. Application of classification methods to agricultural soil data Figure 2.2: Soil map of Maharashtra [118]. The geographical study areas are highlighted with outlines: red (Marathwada), blue (North Maharashtra) and pink (Paschim Maharashtra). fine, mixed, mesic, typic haploxerepts. The study used 217 patterns collected from Kurdistan Province, North-West Iran, achieving accuracy and Cohen kappa [131] of 71% and 69% using DT and ANN, respectively. Figure 2.2 shows1the different types of soil in the state of Maharashtra (label 6 locates the soil type of Marathwada). In our study, three soil classes are considered according to its texture. Light soil has large proportion of sand, low parameter levels and ability to hold water. Sandy, peaty and chalky soils are subtypes of light soil. Medium (loam) soil contains silt, clay and humus (decayed matter) and it is suitable for several crops, being the predominant type in Marathwada. Heavy soil contains more moisture 1http://eusoils.jrc.ec.europa.eu/esdb_archive/EuDASM/Asia/images/maps/download/ in3010_3so.jpg (visited March, 29, 2017). 2.6. Classification methods 15 Figure 2.3: Maps with the geographical location of crop (left panel) and soil type (right panel) data over several districts in the Marathwada region (red outline in the map of Figure 2.2). Each point locates a different village, from where several patterns are recorded. and sticky lump, due to the high proportion of silt (slightly larger particles of rock) or clay (small particles of rock). The classification of soil type uses the inputs listed in the last line of Table 2.2, although the available data only contains patterns of classes light and medium. A geographical plot of the soil type data is shown in the Figure 2.3 (right panel). Locations are widely distributed among the Jalna and Bid districts of the Marathwada region. 2.6 Classification methods In order to develop the classification for the ten problems described in the previous section, we used 20 classifiers selected among the ones which provided the best performance in the comprehensive comparison [28]. These classifiers are implemented in the Java programming language using the Weka data mining software [47], in the R statistical computing language [103], in the C++ programming language and in the Matlab platform [72]. Henceforth, the suffix of the classifier name shows the implementation used: _w, _r, _c and _m mean Weka, R, C++ and Matlab, respectively. The classifiers, grouped by families, are described in the following enumeration alongside with their tunable metaparameters (the symbol # means ‘number of’, e.g. #inputs ‘number of inputs’). 16 Chapter 2. Application of classification methods to agricultural soil data I. Decision trees 1. j48_w is the Java implementation of the C4.5 decision tree [100] provided by the Weka data mining software. This binary decision tree is composed of internal and leaf nodes, where each one does a binary splitting of a given input. Each lead node gives as output a class label. The C4.5 method selects the input with the highest normalized information gain (also called Kullback-Leibler divergence, which is a measure of the non-symmetric difference between two class probability distributions) to be splitted for a new node of the tree. The method then proceeds by recurrence on the two subsets in which the node splits the training set, and the new nodes (internal or leafs) are added as children of the current node. When all the training patterns in the subset belong to the same class, a leaf node (whose output is that class) is created. When no input provides information gain, or when a training pattern of an unseen class appears, an upper level node is created for that class. For a test pattern, the decision tree starts from the root node and travels down the tree, according to the splittings in each nodes and the input values, until it reaches a leaf node, giving the class associated to this node as output. The only hyperparameter of the J48 tree is the pruning confidence threshold Cwith values 10, 20 and 30. 2. rt_w is the random tree provided by Weka. Each node of the tree splits ⌊log2(#inputs)+ 1⌋randomly chosen inputs at each node of the tree, without any pruning. 3. rpt_r is the recursive partitioning method [9] provided by the rpart R package. This method creates a decision tree by binary splitting the inputs and partitioning the training patterns in subsets in a recursive way. 4. dj48_w is the decorate ensemble of J48 tree base classifiers with high diversity [76], provided by Weka. Decorate is an acronym for diverse ensemble creation by oppositional relabeling of artificial training examples. This method iteratively generates an ensemble by adding trained base classifiers (in our case, J48 trees). The first one is trained on the original training set, but the following ones use a fraction of patterns which are artificially generated. These labels are chosen to be maximally different from the actual ensemble output, in order to increase its diversity. Any new base classifier which reduces the ensemble accuracy is rejected. The iterative classifier addition continues until a specified ensemble size is reached. In our case, this size is a hyperparameter tuned with values 10, 15 and 30. 2.6. Classification methods 17 II. Rule-based classifiers 5. dtnb_w is the hybrid decision table-naive Bayes classifier provided by Weka [46]. A decision table (DT) is a lookup table which associates observed frequencies of selected inputs (which are expectedly highly discriminative) to class probability estimates. These inputs are selected using forward search using as objective the maximization of cross-validation performance. The DTNB uses two disjoint subsets of inputs: one for DT and other for naive Bayes (NB). At each search iteration, DTNB decides whether split the input in one or the other subset. The decision is oriented to maximize the area-under-curve of DTNB in a cross-validation. Initially, all the inputs are modeled by DT, but during search the selected inputs are modeled using NB, and the remaining ones by DT. The output class for a test pattern is the one with the highest weighted sum of DT and NB probabilities for that class and pattern (the last one divided by the prior class probability). The hyperparameter is the type of cross-validation, tuned with 4 values: leave-one-out, 4-fold, 5-fold and 10-fold. 6. jrp_w: repeated incremental pruning to produce error reduction (RIPPER), provided by Weka [21]. The RIPPER classifier starts from an empty rule set, and it iteratively: 1) creates a new rule by adding conditions to it until the rule is 100% accurate; 2) tries all the values of the inputs and selects the condition with the highest information gain; and 3) prunes the rule set using a specific metric. Once the rule is created, all the training patterns covered by the rule, both positives and negatives, are discarded. These three steps are repeated until the description length of the rule set reaches a maximum size, until the whole training set has been discarded, or until the error rate overcomes 50%. Finally, an optimized rule set is created in the following way: for each rule, two variants are generated: one from an empty rule using steps 1-2 above, and the other by adding antecedents to the rule. Then, one of the three rules (the original one and the two variants), the one with the smallest description length, is selected to be included in the optimized rule set. Besides, rules that increase the description length of the rule set are removed. The hyperparameters are the number of folds for reduced error pruning (4 values from 3 to 10) and the number of optimization runs (2,3 and 4). III. Bagging ensembles 7. bg_r is the bagging ensemble [7] of decision tree base classifiers, provided by the ipred (bagging ensemble) and rpart (base classifiers) packages. Bagging is an acronym for 18 Chapter 2. Application of classification methods to agricultural soil data bootstrap aggregating, a method to build an ensemble of base learners (classifiers or regressors) based on training several base classifiers (in our case, rpart decision trees) on different random bootstrap samples of the training set. The sampling is with replacement, which means that training patterns can be repeated. A new test pattern is classified using a voting scheme over all the base classifiers. The bagging method, proposed initially for decision tree classifiers, aims to reduce the variance of the ensemble in order to avoid overfitting (i.e., a good learning of the training set which does not mean good learning of unseen test patterns) and to improve the stability of the single classifiers. 8. bgf_r is a bagging ensemble of flexible discriminant analysis (FDA) base classifiers which uses the bagFDA and fda functions provided by the caret and mda packages, respectively [48]. The FDA classifier is a generalization of LDA for non-linear regression using optimal class scorings. In our case, we use the multi-variate adaptive resonance splines (MARS), provided by the earth package, as basis functions for FDA. The hyperparameters are the polynomial degree (1 and 2) of the polynomial splines and the maximum number of terms to keep in the pruned model (10 values from 2 to 11). 9. bgp_w is a bagging ensemble of pruned partial C4.5 decision trees (PART) provided by Weka [32]. The PART trees combine C4.5 and RIPPER (see jrp_w classifier above) trees. Both create a starting set of rules, which are refined by dropping rules (C4.5) or adjust them (RIPPER) by dropping the tail of a rule based on empirical error on a separate training set until a stop criterion (minimum description length heuristic) is met. Both are two-step methods which do global optimization. On the contrary, the PART creates a rule, removes the training patterns covered by the rule, and continues creating new rules recursively until the whole training set is processed. A partial tree is a tree with branches to undefined subtrees. The training process integrates the tree creation and pruning to find a stable (which can not be simplified) subtree. Initially, a split is selected and the training set is divided in two subsets. Both subsets are expanded, first the one with the smallest average entropy, which probably result in a smaller subtree and a more general rule. The process continues recursively until a subset is expanded into a leaf. When an internal node in the tree has only leafs, the method searches whether the node can be replaced by a leaf. The only hyperparameter is the bag size P, given as the percentage of the training set, tuned with values 25%, 50%, 75% and 100%. 2.6. Classification methods 19 IV. Boosting ensembles 10. ab_r is an Adaboost.M1 [33] ensemble of classification trees implemented by the boosting function in the adabag package [2]. Adaboost is an acronym of adaptive boosting, a training method for classifier ensembles where each base classifier (in our case, classification trees) is biased to learn better those training patterns which have been misclassified by the previous classifiers of the ensemble. This is done by weighting each training pattern with the error on that sample. At each iteration, the weight of the base classifier with the smallest weighted error (sum of weights of the misclassified training patterns) is updated using that error. The weights decay are normalized again in each iteration. The test output is decided by a weighted voting among the base classifiers. V. Nearest neighbors 11. knn_r is the K-nearest neighbor [110] classifier, implemented by the knn function of the class R package. The output class for a test pattern is decided by voting over its K nearest training patterns. The number Kof neighbors is the only hyperparameter, tuned with 13 values from 1 to 37. VI. Neural networks 12. elm_m is the extreme learning machine (ELM). The ELM is a single hidden layer feed-forward network which assigns random values to the input weights, calculates the hidden neuron outputs H using the activation function and calculates the output weight matrix B as the product of the Moore-Penrose pseudo-inverse of the H matrix multiplied by the desired output [53]. The test output is the product of H (dependent on the test pattern) times B. This network does not require iterative training, it does not fall in local minima, it avoids tunable hyperparameters as learning rate or momentum, and its training is much simpler than classical neural networks and support vector machines. We used the publicly available Matlab code2, and the tunable hyperparameters are the transfer function (six functions: sinus, signum, hard limit, triangular basis, radial basis and sigmoid) and the number of neurons in the hidden layer (20 values between 3 and 200). 13. gelm_m is the ELM with Gaussian kernel tuning the regularization hyperparameter C with 20 values in the set {2i}14 −5and the kernel spread with 25 values in the set {2i}8 −16. 2http://extreme-learning-machines.org (visited March, 29, 2017). 20 Chapter 2. Application of classification methods to agricultural soil data 14. mlp_m is the classical multi-layer perceptron neural network implemented by the train function in the Matlab Neural Network Toolbox. The network has only one hidden layer, whose number of neurons is a hyperparameter tuned with 11 values from 3 to 30. 15. rbf_m is the radial basis function (RBF) neural network [95], implemented by the newrb Matlab function. The RBF is a two-layer network whose hidden layer is composed by radial basis neurons, usually with Gaussian activation, which are iteratively added during training. Starting from an empty layer, the training set is classified and the pattern with the highest root mean square error is selected as weight vector for a new hidden neuron, being the center of the neuron Gaussian activation. Then, the output weights between the hidden and output layers are calculated to minimize the error. The process is repeated, adding new hidden neurons until the error falls below a goal, or the maximum number of hidden neurons (which by default is the number of training patterns) is reached. The only hyperparameter is the spread of the Gaussian activations for the hidden neurons, tuned with 29 values from 0.1 to 70. 16. pnn_m is the probabilistic neural network [123], implemented by the newpnn Matlab function. This network has only one hidden layer whose neurons are radial basis functions centered in the training patterns. The weights of the output layer are the class labels for the training patterns. The only hyperparameter is the Gaussian spread, tuned with 10 values between 0.01 and 10. VII. Support vector machines 17. svm_c is the support vector classifier [15] with Gaussian kernel implemented by the LIBSVM library3in C++, tuning the regularization hyperparameter Cwith 20 values in the set {2i}14 −5, and the Gaussian spread with 25 values in the set {2i}8 −16. VIII. Random forests 18. rf_r is the random forest ensemble [8] of random tree base classifiers implemented by the randomForest package in R. Each random tree is trained using feature bagging, which randomly selects an input subset with √#inputs items at each candidate split. If some input is very useful to predict the output (i.e., the input is very related to the 3http://www.csie.ntu.edu.tw/∼cjlin/libsvm (visited March, 29, 2017). 2.7. Experimental setup 21 output), that input will be selected in many base trees, and they might be highly correlated, which is undesirable to achieve a diverse ensemble of trees. In order to avoid this, and to correct the overfitting of individual trees, they are trained in different bootstrap samples of the training set, and voting is used to give a class output for a test pattern. The base classifiers should not be strongly correlated among them, so that voting is expected to reduce its variance without increasing its bias. 19. rf_w is an alternative implementation of random forest provided by Weka. We tune the number of trees in the forest with values 100, 250, 500 and 750. 20. rtf_w is the rotation forest ensemble [112] of J48 base classifiers, implemented in Weka. The rotation forest simultaneously promotes classifier accuracy and ensemble diversity. The input set is randomly split into a number of subsets. For each input subset, PCA is applied on a bootstrap sample of the training set. Thus, each input subset represents a different rotation of the original input space. Each J48 base tree is trained on the whole collection of training patterns with all the principal components in order to promote accuracy. On the other hand, each classifier uses a different input subset of the principal components, in order to promote diversity. For a test pattern, each base tree gives a probability to assign that pattern to each class, so the ensemble assigns the test pattern to the class with the highest total probability summed over all the trees. The hyperparameters are the percentage of patterns to be removed in the bootstrap sample (5 values from 10% to 75%) and the number of iterations (5 values from 10 to 50). 2.7 Experimental setup In order to apply the classifiers described in subsection 2.6 to the datasets described in subsections 2.1-2.3, we used as quality measure to evaluate the classification performance the Cohen kappa [131], henceforth denoted by κ and measured in %. The Cohen κ evaluates the classification accuracy discarding the probability of classifier success by chance, being defined as: κ =100 pa−pe s−pe,pa= C ∑ i=1 nii pe=1 s C ∑ i=1 C ∑ j=1 nij! C ∑ k=1 nki!;s= C ∑ i=1 C ∑ j=1 nij 28 Chapter 2. Application of classification methods to agricultural soil data 0.75% where the diagonal value is 0.37%; N and SAL, being more probable to classify N patterns as SAL (10.45%) than the opposite (6.66%); SAL and MAL, with more MAL patterns (7.78%) classified as SAL than MAL patterns correctly classified (2.67%). The remaining non-diagonal values are below 1%. The SE and PP are only acceptable (above 68%) for the most populated classes N and SAL, while the less populated classes SA and MAL exhibit worse results (SE about 20-30%). pH SA N SAL MAL SE(%) PP(%) SA 0.37 0.64 0.27 0 29.2 26.9 N 0.75 26.92 10.45 0.27 70.1 74.6 SAL 0.21 6.66 39.66 1.44 82.7 68.2 MAL 0.05 1.87 7.78 2.67 21.6 61.0 Table 2.7: Confusion matrix (in %) of the best classifier (rf_r, with κ = 47.32% and accuracy= 69.63%) for soil pH classification. Labels SA, N, SAL and MAL mean slightly acidic, neutral, slightly alkaline and moderately alkaline respectively. 2.12 Classification of crop The rf_r achieves the best κ (88.13%, with accuracy above 90%) also for crop classification (Table 2.3), followed by rf_w (87.64%), bgp_w (86.37%) and rtf_w (86.34%). The ab_r and svm_c, which achieved good results in the previous problems, also work well. The confusion matrix of the rf_r for this problem (Table 2.8) has almost all the values outside the diagonal under 1%. Only the percentage (3.34%) of bajra(R) patterns classified as cotton(R) is higher than the corresponding diagonal term (2.89%), which reduces the SE of class bajra(R) to 45.1%. The reason is that the usual crop rotation in Marathwada is between both classes, so their patterns are very similar. Besides, cotton(R) is more important and populated than bajra(R), whose patterns tend to be classified as cotton(R). All the remaining classes have SE above 90% and PP above 80%. 2.13 Classification of soil type The best soil type classifier is dtnb_w (hybrid classifier of decision table and naive Bayes), with κ =97.82%, followed by the group of methods which were the best in the previous classification problems: rf_w and rtf_w (96.80%), ab_r (96.65%), rf_r (96.37%) and svm_c 2.14. Comparison among regions 29 Crop Bajra(R) Cotton(I) Cotton(R) Soybean(R) SE(%) PP(%) Bajra(R) 2.89 0.04 3.34 0.14 45.1 79.8 Cotton(I) 0 10.14 0.35 0.76 90.1 95.7 Cotton(R) 0.74 0.07 23.14 0.81 93.5 85.3 Soybean(R) 0 0.34 0.32 56.90 97.1 97.1 Table 2.8: Confusion matrix (in %) of rf_r for crop classification ( κ =88.13%, accuracy=93.09%). Labels R and I mean rainfed and irrigated respectively. (95.8%). The confusion matrix of the two best classifiers (Table 2.9) shows similarly low non-diagonal values for classes L and M, with SE and PP values above 97%. dtnb_w rf_w Soil L M SE(%) PP(%) L M SE(%) PP(%) L27.96 0.47 98.33 99.42 27.73 0.71 97.50 97.17 M 0.41 71.15 98.54 99.34 0.59 70.97 97.91 99.01 Table 2.9: The confusion matrix (in %) of the two best classifiers (dtnb_w and rf_w, κ = 97.82% and 96.80% respectively) for soil type classification. Labels L and M mean light and medium respectively. 2.14 Comparison among regions As well as the data from the region of Marathwada, we also have data available from regions Paschim-Maharashtra and North-Maharashtra, in the same state of Maharashtra. An interesting issue is how valid are the classifiers, trained using data from one region, to test data from different regions. In other words, what is the representativity and generalization ability of the trained classifiers with respect to regions? We developed experiments training and tuning the metaparameters of the classifiers with patterns from one region, and then testing the trained and tuned classifier with data from different regions. Specifically, we have data from Marathwada, North Maharashtra and Paschim Maharashtra (henceforth labeled as M, NM and PM, respectively) highlighted in Figure 2.2, for village-wise OC-F, P2O5-F, Mn-F, Fe-F and pH problems. The corresponding experiments for N2O,P2O5and K2O, crop and soil type classification could not be developed due to the lack of data for regions NM and PM. Table 2.10 reports the best classifier for each classification problem and combination of training and test regions, and the κ that it achieves (e.g. M-PM in the leftmost column means training and test using data from Marathwada and Paschim-Maharashtra, respectively). No 30 Chapter 2. Application of classification methods to agricultural soil data OC-F P2O5-F Mn-F Fe-F pH Regions Best κ Best κ Best κ Best κ Best κ M-NM rpt_r 17.21 bgp_r 15 bg_r 65.98 rpt_r 66.78 rtf_w 23.97 M-PM elm_m 9.66 rf_r 100 bg_r 66.62 rpt_r 70.52 rtf_w 66.22 NM-M rtf_w 32.60 rf_r 34.10 rt_w 45.70 svm_c 42.50 — — NM-PM pnn_m 14.99 bg_r 100 bg_r 100 j48_w 100 — — PM-M rt_w 12.73 bgf_r 34.57 bgf_r 45.15 elm_m 44.33 knn_r 69.61 PM-NM pnn_m 11.13 bg_r 100 bg_r 100 svm_c 100 rf_r 48.23 Table 2.10: Values of κ (in %) achieved by the best classifier training and testing with patterns of different regions (first column, see text for region labels). The symbol ‘—’ means that data are not available. classifier achieves good results for OC-F, no matter the combination of regions used for training and testing, and the best result is κ =32.6% with rtf_w. The results for P2O5-F are very good ( κ =100%) for NM-PM and PM-NM (both with bg_r) and for M-PM (with rf_r), but much worse (between 15% and 35%) for the remaining combinations. This is surprising, because M-PM works well, but PM-M works much worse, suggesting that data from region M are somehow representative for data from PM, but the opposite is not true. For Mn-F and Fe-F problems, both NM-PM and PM-NM work well (100% with bg_r for Mn-F, and with j48_w and svm_c for Fe-F), so the data seem to be valid between NM and PM. The performance of cases M-NM and M-PM are about 66-71%, so the data from M are somehow valid for the other two regions, but the inverse combinations (NM-M and PM-M) are about 42-46%, so NM and PM data are not so valid for region M. As a conclusion: 1) the data for OC-F classification are not valid among regions; 2) the data for P2O5-F, Mn-F and Fe-F classifications are valid only between North Maharashtra and Paschim Maharashtra, and from Maharashtra to Paschim-Maharashtra; 3) the data for P2O5-F classification are compatible from Marathwada to Paschim Maharashtra, but not the inverse; and 4) the data for Mn-F and Fe-F classifications are slightly compatible (about 66-70%) from Marathwada to the other regions, but not the inverse. CHAPTER 3 APPLICATION OF REGRESSION METHODS TO GENERAL DATASETS Previously to apply regression methods on the different datasets of the soil problem, we developed an experimental comparative of the most popular regression methods in order to know which are the best methods for non-specific datasets. For this work we used the benchmark datasets of the University of California at Irving (UCI) machine learning repository. The current chapter presents and discusses this comparison. 3.1 Regression methods We have applied a wide collection of 76 regressors which belong to several families. The majority of them (68 regressors) are selected from the caret model list1and implemented in R. Instead of using the interface provided by caret (train function), we run the regressors directly using the corresponding R packages (see the detailed list below), in order to have control on the execution of each single model, and to avoid the execution of some regressors not included in the caret list. In fact, we also included other four popular methods implemented in other platforms: deep learning neural network, using the module dlkeras in Python (named dlkeras); support vector regression, using the LibSVM library in C++ (named svr); generalized regression neural network and extreme learning machine with Gaussian kernels in Matlab (named grnn and elm-kernel respectively). Some regression models on the caret list gave errors (list 1http://topepo.github.io/caret/train-models-by-tag.html (visited March, 29, 2017). 32 Chapter 3. Application of regression methods to general datasets them). The model operation is optimized by tuning the set of hyperparameters specified in the caret regressor list. The hyperparameter values used are specified by the caret package (getModelInfo function), and they are different for each dataset. For each regressor, we use the list of tunable hyperparameters specified by the caret package in the previous link, and a number of values used for tuning them. These values are listed in the regressor description below. Note that for some regressors (e.g. gaussprRadial) and datasets, the caret function getModelInfo(...) returns a value list with less items than the number specified in values.txt, and even sometimes just one value is used, so although the caret website specifies that hyperparameter as tunable, in the practice it uses only one value, so it is not tuned at all. The regressors in other languages use pre-specified values equal for all datasets. The 76 regressors are described in the following list, grouped by families. I. Linear regression 1. lm is the linear model provided by the stats package [14]. Collinear inputs exhibit undefined regression coefficients (as returned by lm), so they are discarded for lm and many other regressors in the list. 2. rlm implements robust linear model (MASS package), fitted using iteratively re-weighted least squares with maximum likelihood type estimation, which is robust to outliers in the output although not in inputs [54]. The only hyperparameter is the psi functions: Huber, which provides a convex optimization problem; Hampel and Tukey bisquare, both with local minima. II. Generalized linear regression 3. glm is the generalized linear model provided by the stats package [25], which combines a probability distribution (e.g. Gaussian, binomial, Poisson, etc.), a linear predictor and the link function corresponding to the distribution (which may be non-Gaussian), which relates the output mean and the inputs. They are used to model positive values, categorical or ordinal data. 4. penalized is the penalized linear regression (penalized package), which fits generalized linear models with a combination of L1 and L2 penalties. The L1 penalty, also called least absolute shrinkage and selection operator (lasso), penalizes the sum of absolute values of the coefficients, thus reducing the coefficients of inputs which are not relevant 3.1. Regression methods 33 similarly to input selection. The L2 penalty (also called ridge) penalizes the sum of squared coefficients, reducing the consequences of input collinearity. The regression is regularized by weighting both penalties [41], whose weights (given by hyperparameters lambda1 and lambda2, with 5 and 4 values respectively) are tuned. 5. glmnet is the lasso and elastic-net regularized GLM provided by the glm package [120]. The glmnet model uses penalized maximum likelihood to fit a GLM for the lasso and elastic-net non-convex penalties. The mixing percentage alpha is tuned with 5 values, including alpha=1 (resp. alpha<1), which corresponds to lasso penalty, and alpha<1 for elastic-net penalty (alpha=0 corresponds to ridge regression penalty). The glmnet function already tries a number of values (100 by default) for the regularization parameter lambda, so it was not tuned despite of being included among the glm hyperparameters in the caret list. 6. glmStepAIC is the generalized linear model with stepwise feature selection [111] using the Akaike information criterion (stepAIC function in the MASS package). III. Least squares 7. nnls is the non-negative least squares regression (nnls package), it solves for xthe optimization problem minx|Ax −b|subject to x≤0 using the Lawson-Hanson NNLS method [67]. 8. krlsRadial is the radial basis function kernel regularized least squares regression (KRLS package), which uses Gaussian radial basis functions to learn the best fitting function minimizing the squared loss of a Tikhonov regularization problem. The krls method learns a closed form function so interpretable as ordinary regression models. The only hyperparameter is the kernel spread (10 values). The krls determines the trade-off between model fit and complexity (lambda parameter) by minimizing the sum of squared leave-one-out errors, so the caret getModelInfo function does not provide values for it, despite being listed as a tunable parameter. IV. Partial least squares 9. spls is the sparse partial least squares regression (spls package), which introduces sparse linear combinations of the inputs in the dimensionality reduction of PLS in order to 34 Chapter 3. Application of regression methods to general datasets avoid lack of consistency of PLS with high dimensional patterns [20]. The hyperparameters are the number of latent components (K) and the threshold eta (3 and 7 values, respectively). 10. simpls fits a PLS regression model with the simpls method (plsr function in the pls package, with method=simpls). The PLS method projects the inputs and the output to a new space and searches the direction in the input space which explains the maximum output variance, being particularly useful when there are more inputs than patterns and inputs are collinear. The simpls method [115] directly calculates the PLS factors as linear combinations of the inputs maximizing a covariance criterion with orthogonality and normalization constraints. The only hyperparameter is the number of components used by the simpls model (10 values). 11. kernelpls is the PLS regression with method=kernelpls [116] in the same function and package, using the same hyperparameter setting. 12. enpls.fs is an ensemble of sparse partial least squares regressors provided by the enpls package [134]. The number of components (1 value) is specified by the caret function getModelInfo for each dataset, but the threshold argument, specified as a hyperparameter by the caret list, is missing in the enpls.fit function. 13. plsRglm is the partial least squares generalized linear model (plsRglm package) with modele=pls-glm-family [4]. The hyperparameters are the number of extracted components and (4 values) and the input significance level (alpha.pvals.expli, 5 values). V. Least absolute shrinkage and selection operator (lasso) 14. lasso does lasso regression, using the enet function in the elasticnet package similarly to ridge, but using the lambda parameter equal to zero to obtain the lasso solution. 15. relaxo develops relaxed lasso (relaxo package), which generalizes the lasso shrinkage method for linear regression [73]. This method overcomes the trade-off between speed and convergence, specially for sparse high-dimensional patterns, in the L2-loss function of the regular lasso, providing sparser solutions with better prediction error. The relaxation and penalty hyperparameters phi and lambda are tuned using 7 and 3 values respectively. 3.1. Regression methods 35 VI. Ridge (or Tikhonov) regression 16. ridge (elasticnet package), which uses the LARS-EN algorithm to compute the elastic net [136] regression model. Elasticnet provides a model for regularization and input selection, grouping together the inputs which are strongly correlated, being more useful when the number of inputs is higher than the number of patterns, as opposed to lasso models (see below). The only hyperparameter is the quadratic penalty (or regularization) parameter (10 non-zero values). 17. spikeslab implements the spike and slab regression [56] uses the spikeslab package to compute weighted generalized ridge regression estimators using Bayesian spike and lab models. This model combines filtering for dimensionality reduction, model averaging using Bayesian model averaging, and variable selection using the gnet estimator. 18. foba develops ridge regression with forward, backward and sparse input selection [135] (foba package). We use the adaptive forward-backward greedy version of foba (default value ’foba’ for the ’type’ argument of the foba function), which does a backward step when the ridge penalized risk increases in less than the parameter nu (0.5 by default) multiplied by the ridge penalized risk reduction in the previous forward step. The hyperparameters are regularization for ridge regression and the number of selected inputs (sparsity) for the prediction (10 and 2 values, respectively). VII. Neural networks 19. mlpWeightDecay is the multi-layer perceptron with one hidden layer and weight decay (mlp function in the RSNNS package, with learnFunc=BackpropWeightDecay). The size of the hidden layers and weight decay are tunable hyperparameters (5 values each one). 20. mlpWeightDecayML is the same network with three hidden layers, tuning their sizes (3 values each one) and the weight decay (5 values). 21. avNNet is the model averaged neural network provided by the caret package. A committee of 5 neural networks [110] of the same size is trained using different random seeds, being averaged to give an output. The hyperparameters are the network size and the weight decay (7 and 3 values, respectively). 36 Chapter 3. Application of regression methods to general datasets 22. rbf is the radial basis function network (RSNNS package) which does a linear combination of basis functions centered around a prototype [1]. The information is locally codified (opposed to globally in the MLP), the training should be faster and the network is more interpretable, although the output might be undefined if a test pattern does not activate any prototype. The only hyperparameter is the size of the hidden layer (10 values). 23. grnn is the generalized regression neural network [124] implemented by the Matlab neural network toolbox. The GRNN is a special type of RBF network: after a clustering of the training set, the nodes of the hidden layer store the cluster centers (the Matlab implementation uses so many clusters as training patterns). The output for a test pattern is a weighted sum of the Gaussian functions centered in the cluster centers, scaled by the cluster populations. During training, whenever a pattern is assigned to that cluster the weight of the Gaussian function corresponding to that cluster is updated using the desired output. The Gaussian spread is a hyperparameter (14 values): large (resp. small) values lead to smooth (resp. close) approximations. 24. elm is the extreme learning machine, which has been already introduced in the section 2.6. For the regression problem we use the elmNN package [53]. The only hyperparameters are the number of hidden neurons (20 values) and the activation function (sinus, radial basis, linear and hyperbolic tangent). 25. elm-kernel is the ELM neural network but with Gaussian kernel [53] using the publicly available Matlab code2. The hyperparameters are regularization Cand kernel spread with values 2−5..214 and 2−16..28(20 and 25 values, respectively). 26. pcaNNet is a multi-layer perceptron neural network with one hidden layer trained on the PCA-mapped training patterns, using the caret and nnet packages. The principal components which account for more than 95% of the data variance are used for training. With a test pattern, it is mapped to the principal component space and the trained pcaNNet model gives an output. Tunable hyperparameters are the size of the hidden layer and the weight decay of the network (7 and 3 values, respectively). 27. bdk is the supervised bi-directional Kohonen network, using the kohonen package [75]. The bdk combines Kohonen maps and counterpropagation networks, using two maps, 2http://www.extreme-learning-machines.org (visited March, 29, 2017). 3.1. Regression methods 37 for inputs and output respectively. In each iteration, the direct (resp. inverse) pass updates only the weights of the input (resp. output) map, using a weighted similarity measurement (Euclidean distance for regression) which involves both maps, leading to a bi-directional updating. The test output is the weight of the winner node of the output map. The hyperparameters are the sizes of both maps (3 values each one) and the initial weight given to the input map in the distance calculation for the output map, and vice versa (2 values). VIII. Deep learning neural networks 28. dlkeras is the deep learning neural network using the Keras module [18] in Python, with three hidden layers tuned with 50 and 75 neurons for each layer (27 combinations). The deep learning methods [51, 69] are very popular, specially for image classification, and we included them in this comparison for regression tasks. 29. dnn implements a deep belief network in R using the DeepNet package using only one hidden layer and tuning the number of neurons with values from 5 to 60 with step 2. We tried with several hidden layers but the results are worse. The weights are initialized using stacked autoencoder (SAE), which gave better results than deep belief network (DBN). IX. Support vector machines 30. svr is the support vector machine for regression, with Gaussian kernel using the LibSVM library [15] with the C++ interface. We tuned the regularization hyparameter C and the kernel spread γ with values 2−5..214 and 2−16..28. 31. svmRadial is another implementation of SVR with Gaussian kernel (ksvm function in the kernlab package) for regression (type=eps-svr), using SMO to solve the quadratic SVM problem and tuning the same parameters as svr (5 and 4 values respectively). 32. rvmRadial is the relevance vector machine [129] with Gaussian kernel (same package). The RVM has the same functional form as the SVM, but using a Bayesian learning framework which reduces the number of basis functions, compared to the SVM, while keeping an accurate prediction. It also avoids tunable (e.g. regularization) hyperparameters of the SVM, but uses a method similar to Expectation-Maximization which 44 Chapter 3. Application of regression methods to general datasets number of principal components (1-3) and the threshold for retaining the input scores (2 values). It is provided by the superpc package. XVII. Generalized additive models 65. gam is the generalized additive model [133] using splines with the mgcv package. The GAM is a GLM whose linear predictor is a sum of smooth functions (penalized regression splines) of the covariates. The estimation of the spline parameters uses by default the generalized cross validation criterion. The hyperparameter is the select parameter with values true and false: in the former case, an extra penalty term is added to each function penalizing its wiggliness (waving). 66. gamboost is the boosted generalized additive model [12] using the mboost package. This method is a gradient boosting ensemble which minimizes (computing its negative gradient) a weighted sum of the loss function evaluated at the training patterns. The base regressors are component-wise models (P-splines with a B-spline base, by default). The only hyperparameter is the number of initial boosting iterations (mstop), with 10 values. XVIII. Gaussian processes 67. gaussprLinear implements Gaussian process regression with linear (vanilladot) kernel using the kernlab package (function gausspr). 68. gaussprRadial uses the same function with Gaussian (rbfdot) kernel, with spread automatically calculated (option kpar=1). 69. gaussprPoly is the same method with polynomial (polydot) kernel, tuning the kernel hyperparameters degree and scale (3 values each one). XIX. Quantile regression 70. rqlasso develops quantile regression with LASSO penalty, using the rq.lasso.fit in the rqPen package. The quantile regression models optimize the so-called quantile regression error, which uses the tilted absolute value instead the square as RMSE. This tilted function applies asymmetric weights to positive/negative errors, computing conditional quantiles of the predictive distribution. The qrlasso method fits a quantile regression model with the LASSO penalty [80], tuning the lambda hyperparameter (10 values). 3.1. Regression methods 45 71. rqnc performs non-convex penalized quantile regression, with the rq.nc.fit function in the previous package. This regressor performs penalized quantile regression using local linear approximation [137] to maximize the penalized likelihood for non-convex penalties. The two hyperparameters are lambdas (10 values) and 2 penalty types: MCP (minimax concave penalty) and SCAD (smoothly clipped absolute deviation). 72. qrnn is the quantile regression neural network [13], with the qrnn package. A QRNN is a neural network where the transfer function is a ramp, being the error function the quantile regression error. The hyperparameters are number of hidden neurons and the penalty for weight decay regularization (7 and 3 values respectively). XX. Other methods 73. lars is the least angle regression [26] using the Lars package. Lars is a model selection method which is less greedy than the typical forward selection methods. It starts with zero coefficients for all the inputs and finds the input imost correlated with the output, increasing step-by-step its coefficient until another input jhas high correlation with the current residual (error, or difference between desired and real output). Lars increases the coefficients of inputs iand jin the equiangular direction between inputs iand j until some other input kis so correlated with the residual as input j. Then, it proceeds in the equiangular direction among i,jand k, which is the “least angle direction”, and so on until all the coefficients are non-zero (i.e., all the inputs are in the model). The lasso and fraction options are specified for training and prediction respectively, and the fraction hyperparameter is tuned with 10 values. The number of terms is tuned with 10 values. 74. earth is the multivariate adaptive regression spline (MARS) using the earth package [34]. It uses a expansion of product spline functions to model non-linear data and interactions among inputs. The spline number and parameters are automatically determined from the data using recursive partitioning, and distinguishing between additive contributions of each input and interactions among them. The functions are added iteratively to reduce maximally the residual, until its change is too small or a number of iterations is reached. The maximum number of terms in the model is tuned with 15 values. 75. ppr performs the projection pursuit regression [36], in the stats package. The ppr models the output as a sum of averaging functions (mean, median, etc.) of linear combina- 46 Chapter 3. Application of regression methods to general datasets tions of the inputs. The coefficients are iteratively calculated to minimize a projection pursuit (fitting criterion, given by the fraction of unexplained variance which is explained by each function) until it falls below a predefined threshold. 76. sbc is subtractive clustering and fuzzy c-means rules [17], using the frbs package. This method uses SBC to get the cluster center of a fuzzy rule-based system for classification or regression. Initially, each training pattern is weighted by a potential function which decreases with its distances to the remaining centers. The center with the highest potential is selected as a cluster center, and the potential of the remaining centers are updated. The only hyperparameter is the neighborhood radius, tuned with 7 values usually between 0 and 1. The selection of new cluster center and potential updating is repeated until the potentials of the remaining patterns are below a pre-specified fraction of the potential of the first cluster center. Once all the centers are selected, they are optimized using fuzzy C-means. 3.2 Experimental setup We run the 76 previous regressors on a collection composed by 66 regression datasets selected from the UCI Machine Learning Repository [3]. For the current study, we selected 44 out of the 63 datasets labeled as “Regression” by the UCI repository 3, plus other two datasets (Breast cancer Wisconsin original and prognostic, included in the “Classification” UCI list), which are listed in Table 3.1. For space reasons, the table only lists the first two words of the original name of each dataset. The remaining 26 datasets were discarded due to diverse reasons. The “Amazon Access Samples” dataset does not explain what is the output. The “Educational Process Mining (EPM): A Learning Analytics Data Set” has time-series data with a very complex format. The “Challenger USA Space Shuttle O-Ring” has too few patterns (23) and inputs (3). The data are missing in the “Improved Spiral Test Using Digitized Graphics Tablet for Monitoring Parkinson’s Disease” dataset. The “NoisyOffice” dataset is difficult to read because it is composed by printed text images. We were unable to decide the output for the “Tennis Major Tournament Match Statistics” dataset. For large datasets, only the first 2000 patterns are used. The reason is that many regressors are not able to train and test with large datasets, and they give memory errors or simply they do not ever finish. 3http://archive.ics.uci.edu/ml (visited March, 29, 2017). 3.2. Experimental setup 47 Although we selected 46 UCI datasets, some of them generated several datasets, one for each data column which can be used as output for regression. This is the reason why some original datasets in column 1 give several datasets in column 2 (e.g., the Air quality dataset gives five datasets). Therefore, we achieved a list of 66 datasets. All the constant, repeated and collinear inputs, whose NA coefficients in a linear regression model using the lm(...) function in the R stats package, are removed from the dataset. For example, the “Blog feedback” dataset reduced its inputs from 280 to 13. On the contrary, those inputs with discrete values, they are replaced by dummy (also called indicator) inputs. For each discrete input with n values, it is replaced by n−1 dummy inputs with values zero or one: the first value of the original discrete input is codified as zero values for the n−1 dummy inputs; the second value is codified as 1 in the first dummy variable and zero in the remaining ones; the third value is codified as 1 in the second dummy input, and so on. Therefore, those datasets with discrete inputs increase the number of inputs, e.g. dataset “Student perf. mat.” increases its inputs from 32 to 77. In the Table 3.1 those #inputs column where appear two numbers (i.e. 8/23), the first one is the original #inputs, and the second one is the number of inputs used effectively by the regressors, after removing constant, repeated and collinear inputs, and after replacing discrete inputs by their corresponding dummy variables. Ten random partitions are generated for each dataset, using the 50% for training, 25% for validation (in hyperparameter tuning) and 25% for test. Each regressor is trained on the ten training partitions for each combination of its hyperparameter values, and tested on its corresponding validation partition [64]. The performance measure that we use is the root mean square error (RMSE), defined as: RMSE =s1 N N ∑ i=1 (yi−di)2(3.1) where Nis the number of test patterns, being yiand dithe regressor and the desired (right) outputs respectively. For each combination, the average RMSE over the ten partitions is calculated, and the one with the lowest RMSE is selected for testing. Finally, the regressor is trained, using this selected combination of its hyperparameter values, on the ten training partitions, and tested on the ten test partitions. The average RMSE over the ten test sets is used as the final quality measure. Some regressors which are specially sensitive to collinear inputs are trained, for each partition, using only those inputs which are not collinear. Although collinear inputs have been removed from the dataset in the initial preprocessing, sometimes it happens that for certain partition the training patterns have several collinear inputs, which are 48 Chapter 3. Application of regression methods to general datasets Original Datasets #Patterns #Inputs Original Datasets #Patterns #Inputs Air quality CO 1230 8 GPS trajectories 163 10 NMHC Greenhouse gas 2000 326 NO2 Housing 452 13 NOx Individual household 2000 6/5 O3 Insurance company 2000 85/439 Auto_MPG 398 8/23 Istanbul stock 536 8 Automobile 205 26/66 KEGG reaction 2000 27/25 Bike sharing day 731 13/30 KEGG relation 2000 22/17 hour 2000 14/42 Online news 2000 59/55 Blog feedback 2000 280/13 Online video video-trans. 2000 20/8 Breast Wisc. wbc 699 9/81 Parkinson speech 1040 26 B. W. diagnosis wdbc 569 30 Parkinsons tel. motor 2000 16 B. W. prognostic wpbc 198 33/32 total Buzz social buzz-twitter 2000 77 Physicochem. prop. 2000 9 Com. & crime com-crime 1994 122 Relative CT slices 2000 385/355 C. & C. unnorm. 2000 124/126 Servo 167 4/15 Computer hardware com-hd 209 7 Skillcraft master 2000 18 Concrete compressive cond-compress 1030 8 SML2010 2000 21/18 Concrete slump slump 103 9Student perf. mat 395 32/77 comp 7por 649 32/56 flow Tamilnadu electr. elec.hour. 2000 4/3 Condition based compress 2000 16/13 Twin gas sensor 2000 8 turbine UJIIndoorLoc floor 1999 528/373 Energy efficiency cool 768 8/7 lat heat long Fertility 100 9/16 UJIIndoorLoc-Mag curve-lat 2000 9 Forestfires 517 12/39 curve-long Gas sensor drift gas-drift 2000 129 line-lat 10 Gas sensor flow gas-flow 58 438/57 line-long Geogr. origin geo-lat 1059 116/72 wiki4HE 913 52/240 geo-long Wine quality red 1599 11 geo-music-lat 68 white 2000 11 geo-music-long Yacht hydro. 308 6 Table 3.1: Collection of 66 datasets from the UCI repository: original name in the UCI repository; datasets created from the original one; number of patterns and inputs, before and after preprocessing. not collinear considering the whole dataset. That is the reason why these inputs are discarded for these regressors, which otherwise give errors. This happens for regressors earth, evtree, foba, glm, icr, krlsRadial, lasso, lm, nnls, nodeHarvest, plsRglm, qrnn, ridge, SBC and spls. Obviously, discarding collinear inputs by partitions only happens for certain datasets. All the inputs and the output are pre-processed to have zero mean and standard deviation one. Some regressors gave errors for some datasets. In these cases, the programs are configured to give the mean of the desired outputs, which is zero. Most regressors have one, two, three or four 3.3. Results and discussion 49 tunable hyperparameters. As well as the RMSE, another performance measure that we also use is the correlation coefficient, defined as: Correlation =(N ∑ i=1 (yi−¯y))(N ∑ i=1 (di−¯ d)) σ (y) σ (d)(3.2) ¯y=1 N N ∑ i=1 yi,¯ d=1 N N ∑ i=1 di(3.3) σ (y) = v u u t 1 N n ∑ i=1 (yi−¯y)2!, σ (d) = v u u t 1 N n ∑ i=1 (di−¯ d)2!(3.4) where ¯yand σ (y)are the mean and standard deviation of y1,...,yN, respectively, and the same for d1,...,dN. Regressor #Datasets Datasets nodeHarvest 27 bike-day, UJ-lat, twin-gas-sensor, buzz-twitter, air-quality-NMHC, video-transcode, comcrime-unnorm, geo-lat, wine-red, air-quality-O3, skill-craft, KEGG-relation, blog-feedback, gas-drift, SML2010, UJ-MagLine-lat, wine-white, air-quality-NO2, com-crime, UJMagLine-long, air-quality-NOx, bike-hour, physico-protein, KEGG-reaction, park-totalUPDRS, park-motor-UPDRS and UJ-long qrnn 6 greenhouse-net, CT-slices, UJ-long , UJ-floor , UJ-lat, insurance-coil brnn 5 wiki4HE, insurance-coil, UJ-long , UJ-floor , UJ-lat randomGLM 3 gas-drift, insurance-coil and com-crime-unnorm rqlasso 3 wiki4HE , wbc , insurance-coil rqnc 2 wiki4HE, insurance-coil Table 3.2: List of regressors (and datasets) whose execution did not finish within 150 h. or required more than 128 GB RAM. 3.3 Results and discussion We run a collection of 76 regressors over 66 datasets, developing a total of 5016 experiments. These experiments were developed on the CiTIUS cluster, using Intel Xeon E5-2650L microprocessors with 64 GB RAM. As commented above, those regressors which gave errors for some partitions or combinations of the tunable hyperparameters are evaluated as if they give the mean output, which is zero because the desired output is preprocessed to have zero mean. Additionally, some regressors did not finish for certain datasets (see Table 3.2), either due to a huge memory consumption (more than 128 GB) or computation time (an upper limit 50 Chapter 3. Application of regression methods to general datasets 12 14 16 18 20 22 24 26 28 elm−kernel svr extraTrees rf RRF bstTree cubist gbm penalized bagEarth avNNet M5 brnn qrf cforest xgbTree earth gaussprPoly krlsRadial blackboost Friedman rank of the MSE (20 best models) 12 14 16 18 20 22 24 26 28 30 elm−kernel svr cubist extraTrees rf RRF bstTree gbm penalized avNNet qrf bagEarth brnn xgbTree M5 earth cforest dlkeras gaussprPoly krlsRadial Friedman rank of the correlation (20 best models) Figure 3.1: Friedman rank of RMSE (upper panel) and correlation (lower panel) for the 20 best regressors. of 150 hours or 6.25 days was set). These combinations of regressor-dataset which did not finish are 46 of 5016, which represents 0.92% of the total experiments. The nodeHarvest had a specially bad behavior, because it did not finished for 27 datasets, which represents 40.9% of the datasets, and randomGLM overcame the memory limit for 3 datasets. The remaining regressors did not finish within the 150 hours deadline. In order to develop a comparison of the different regressors over the collection of datasets, we used the Friedman rank [37] of the RMSE for the whole regressor collection. This rank 3.3. Results and discussion 51 RMSE rank Correlation rank Order Regressor Rank Regressor Rank Avg. correl. p-value 1 elm-kernel 11.4 elm-kernel 10.5 0.79198 – 2 svr 11.8 svr 11.7 0.79432 1 3 extraTrees 13.4 cubist 12.1 0.78843 1 4 rf 13.4 extraTrees 13.7 0.79071 1 5 RRF 14.3 rf 14.6 0.78632 1 6 bstTree 15.3 RRF 15.0 0.78611 1 7 cubist 15.5 bstTree 15.6 0.77836 1 8 gbm 17.8 gbm 17.4 0.77304 1 9 penalized 19.0 penalized 18.6 0.77478 0.99999 10 bagEarth 20.6 avNNet 20.4 0.76830 0.99001 11 avNNet 22.1 qrf 20.7 0.77053 0.98697 12 M5 22.6 bagEarth 21.9 0.75025 0.83450 13 brnn 22.7 brnn 23.7 0.70599 0.56961 14 qrf 23.7 xgbTree 25.1 0.75767 0.20598 15 earth 25.2 earth 26.4 0.74607 0.05424 16 xgbTree 25.3 M5 26.7 0.75375 0.08772 17 cforest 25.4 cforest 27.8 0.73645 0.02051 18 gaussprPoly 27.4 dlkeras 28.4 0.72999 0.01010 19 krlsRadial 27.9 gaussprPoly 29.5 0.69821 0.00290 20 blackboost 29.1 krlsRadial 29.9 0.65848 0.00129 21 treebag 30.5 ppr 30.4 0.72542 0.00084 22 bag 32.6 treebag 31.1 0.74875 0.00035 23 grnn 32.7 blackboost 31.7 0.68776 0.00011 24 kknn 33.1 rvmRadial 32.2 0.72535 0.00008 25 lars 33.3 kknn 33.3 0.73088 0.00001 26 foba 33.5 svmRadial 33.4 0.73238 0.00001 27 ppr 33.8 foba 33.8 0.68832 0.00001 28 spikeslab 33.9 grnn 33.9 0.72847 0 29 dlkeras 34.1 lars 34.2 0.71300 0 30 rvmRadial 34.2 gaussprRadial 34.7 0.72400 0 31 glmboost 35.2 simpls 36.1 0.70851 0 32 gaussprRadial 35.5 bag 36.3 0.69490 0 33 svmRadial 36.3 glmboost 36.8 0.67369 0 34 simpls 36.8 spikeslab 36.9 0.67195 0 35 kernelpls 38.2 ridge 37.4 0.66220 0 36 BstLm 38.7 kernelpls 37.7 0.70214 0 37 spls 39.6 lasso 39.1 0.70074 0 38 icr 39.7 plsRglm 39.2 0.67618 0 Table 3.3: Friedman rank of the RMSE (left) and of the correlation (right), average correlation and p-value of the Posthoc Friedman Nemenyi test comparing the best regressor to the remaining ones. Continued in Table 3.4. evaluates the position in which each regressor is, in average over all the datasets, when the measure is sorted by decreasing performance (in our case, by increasing RMSE). This rank is calculated as follows: let rij =0, for i=1,...,nand j=1,...,m, being nand mthe 52 Chapter 3. Application of regression methods to general datasets RMSE rank Correlation rank Order Regressor Rank Regressor Rank Avg. correl. 39 plsRglm 39.7 bayesglm 39.4 0.67650 40 ctree2 40.3 rqlasso 39.6 0.68176 41 evtree 41.2 glm 39.8 0.66703 42 SBC 41.2 spls 39.9 0.68100 43 rqlasso 41.6 BstLm 40.1 0.67797 44 rpart 41.6 SBC 40.8 0.66650 45 lasso 41.7 gaussprLinear 40.8 0.66061 46 ridge 42.0 lm 40.8 0.66703 47 glmStepAIC 42.7 glmStepAIC 41.1 0.66918 48 bayesglm 43.0 icr 41.6 0.67397 49 glm 43.2 rlm 41.8 0.66732 50 randomGLM 43.3 rqnc 42.0 0.68974 51 relaxo 43.4 rpart 42.2 0.70422 52 lm 44.2 gam 42.5 0.65881 53 gaussprLinear 44.3 evtree 42.5 0.67557 54 nodeHarvest 44.3 randomGLM 42.8 0.51610 55 rqnc 44.4 ctree2 43.6 0.68371 56 elm 44.7 relaxo 45.3 0.53902 57 rlm 44.8 nodeHarvest 46.5 0.39834 58 gamboost 45.1 elm 47.9 0.66668 59 gam 45.3 gamboost 48.2 0.34728 60 rbf 48.1 xgbLinear 48.9 0.54951 61 bstSm 48.3 pcaNNet 49.1 0.63487 62 xgbLinear 49.3 qrnn 50.1 0.34677 63 nnls 50.6 nnls 50.4 0.59542 64 qrnn 51.0 mlpWeightDecay 50.6 0.64055 65 enpls.fs 53.1 bstSm 51.3 0.26817 66 partDSA 56.3 rbf 51.9 0.55881 67 pcaNNet 57.6 enpls.fs 55.1 0.44590 68 mlpWeightDecay 57.8 partDSA 58.6 0.55021 69 pcr 59.0 bdk 61.1 0.48553 70 mlpWeightDecayML 61.2 superpc 61.5 0.39352 71 Boruta 61.4 pcr 64.0 0.39539 72 dnn 61.6 mlpWeightDecayML 64.0 0.46150 73 superpc 62.1 dnn 65.9 0.36021 74 bdk 62.5 Boruta 67.7 0.01491 75 bartMachine 72.2 bartMachine 73.8 -0.03007 76 glmnet 75.8 glmnet 74.7 -0.04054 Table 3.4: Continuation of Table 3.3. number of regressors and datasets, respectively. For each dataset j, the RMSE values of all the regressors are sorted increasingly. For each regressor ilet be rij =k, where kis the position of the regressor ifor dataset j. The Friedman rank of regressor iis the average of rij for j=1,...,m, i.e., the average position of regressor iover all sorting of RMSE for the 3.3. Results and discussion 53 different datasets. For example, a rank of 5 means that the regressor is, in average over all the datasets, the 5th best regressor. The Figure 3.1 plots the Friedman rank for the RMSE and correlation (in this case, the correlation must be sorted decreasingly) for the best 20 regressors. The extreme learning machine (elm-kernel) achieves the best rank, followed by the support vector regression (svr) with LibSVM and Gaussian kernel, extraTrees (a forest of extremely randomized regression trees) and the random forest (rf). Both the RMSE and correlation ranks are very similar, specially in the first positions. The Tables 3.3 and 3.4 report the complete Friedman ranks for RMSE and correlation for the 76 regressors. The column “Avg. correl.” reports the average correlation of each regressor over all the datasets. We do not report the average RMSE of each regressor over all the datasets because this measure has different ranges for each dataset, depending on the data complexity, so we can not average them directly. However, the correlation coefficient ranges between -1 and +1 for any dataset, so it can be averaged. RMSE Correlation Regressor #times best Regressor #times best Regressor #times best Regressor #times best penalized 17 SBC 2 penalized 16 SBC 2 svr 10 brnn 2 svr 11 brnn 2 extraTrees 9 qrnn 2 extraTrees 9 bstTree 2 cubist 6 elm-kernel 1 cubist 6 qrnn 2 avNNet 4 grnn 1 avNNet 5 elm-kernel 1 gbm 4 rf 1 qrf 4 qrf 4 bstTree 1 glm 3 bagEarth 2 gbm 3 Table 3.5: List of the regressors which achieve the best RMSE (left part) or correlation (right part) for some dataset. Following the suggestions in [22], we also developed a post-hoc Friedman-Nemenyi test, using the PMCMR package [99] provided by the R language, in order to evaluate the statistical significance of the differences in the correlation coefficients of the different regressors. The last column of Table 3.3 reports the p-values of the comparisons between the best regressor (elm-kernel) and the remaining ones. The p-values of the 9 best regressors are almost 1, and it decreases very fast in such a way that for cforest (position 17) the p-value is already under 0.05 (the threshold value). Therefore, the differences between elm-kernel and regressors under position 17 are statistically significant. Given the meaning of the Friedman rank, it is also interesting to see, for each regressor, how many datasets it achieves the best result. The Table 3.5 lists the regressors which achieve 60 Chapter 3. Application of regression methods to general datasets 0.66 0.68 0.7 0.72 0.74 0.76 0.78 svr elm−kernel extraTrees cubist bstTree penalized M5 bagEarth dlkeras simpls brnn lasso gaussprPoly foba rqlasso icr lm SBC gam krlsRadial Average correlation Figure 3.6: Average correlation (sorted decreasingly) of the best regressors of each family. Since the svmRadial in the kernlab package seems to use the same LibSVM implementation as svr, perhaps the cause of the difference between them is the reduced collection of hyperparameter values for svmRadial, which are generated by the getModelInfo function of the caret package. 3) All the random forest regressors work well (within the first 17 positions), except Boruta, which suggests that the feature selection is not working properly. 4) Among boosting regressors, the gbm (position 8, just one below bstTree) and xgbTree are competitive. 5) The R implementation of the deep learning (dnn) is not competitive to the Keras module in Python. The discussion about elapsed time spent by each regressor leads us to develop a Friedman rank of these times over all the datasets. We can not simply average the times of each regressor over all the datasets due to the large difference between the times spent for different datasets, whose number of patterns and inputs are very different. The measurements include training and test times for test stage, divided by the number (10) of test trials. We do not include the time spent for hyperparameter tuning, because it is obviously conditioned by the number of tunable parameters and the number of values tried for each hyperparameter. The Table 3.8 reports the Friedman rank for the times over all the datasets. The first conclusion is that, unfortunately, the fastest regressors are also worst ones: dnn, Boruta, pcr, glm and others. 3.3. Results and discussion 61 Family Regressor Pos. Regressor Pos. Regressor Pos. Regressor Pos. Neural networks avNNet 10 grnn 27 elm 58 pcaNNet 61 mlpwd 64 rbf 66 bdk 69 mlpwdML 72 Support vector machines rvmRadial 24 svmRadial 26 Prototype models kknn 25 Random forests rf 5 RRF 6 qrf 11 cforest 17 Boruta 74 Boosting ensembles gbm 8 xgbTree 14 blackboost 23 glmboost 33 BstLm 43 randomGLM 54 xgbLinear 60 bstSm 65 Generalized linear regression glm 42 glmStepAIC 46 glmnet 76 Bagging ensembles treebag 22 bag 32 Bayesian models bayesglm 37 bartMachine 75 Regression trees rpart 51 evtree 52 ctree2 55 nodeHarvest 57 partDSA 68 Deep learning dnn 72 Gaussian processes gaussprRadial 30 gaussprLinear 47 Least squares kernelpls 35 plsRglm 38 spls 41 enpls.fs 67 Ridge spikeslab 34 ridge 36 Partial least squares nnls 63 Lasso relaxo 56 Quantile regression rqnc 49 qrnn 62 Other methods lars 29 earth 16 ppr 21 Linear regression lm 45 rlm 50 Principal component analysis superpc 70 pcr 71 Generalized additive models gamboost 59 Table 3.7: Positions of the regressors, grouped by families, in the Friedman rank (data extracted from Tables 3.3-3.4), in the same order as Figure 3.6. The acronym mlpwd means mlpWeightDecay. Conversely, the two slowest regressors are the best ones: elm-kernel and svr. For illustrative purposes, the average times spent by elm-kernel and svr are 1772 seconds and 5705 seconds respectively, while the fastest regressor (dnn) spent only 3.21 seconds. Actually, the elmkernel is about 5 times faster than svr. Among the 20 regressors with the highest correlations, the M5 exhibits the best time rank, being in the position 16 of the correlation rank (average time 4.13 seconds) and in the position 6 of the time rank. Besides, regressor earth is in the positions 15 and 17 of the correlation and time ranks, respectively (average time 6.0 seconds). They are the only regressors which are among the 20 first of both ranks. Since there is no clear trade-off between correlation and time, we can discuss the times of the 20 best regressors according to correlation. The Figure 3.7 plots the time rank against the correlation rank for these regressors. The two best regressors elm-kernel and svr are very high on the left, which means low correlation rank (good performance) and high time rank (low speed). Perhaps the best trade-off between correlation and time might be provided by 62 Chapter 3. Application of regression methods to general datasets Pos. Regressor Rank Pos. Regressor Rank 1 dnn 2.55 39 gaussprLinear 36.09 2 Boruta 5.53 40 rqnc 38.82 3 pcr 7.89 41 icr 39.14 4 glm 12.11 42 avNNet 40.29 5 ridge 12.52 43 plsRglm 40.48 6 M5 14.02 44 treebag 42.30 7 glmnet 16.98 45 penalized 43.24 8 superpc 17.00 46 gamboost 44.79 9 kknn 19.55 47 glmStepAIC 45.17 10 rvmRadial 20.70 48 cubist 45.64 11 foba 20.91 49 brnn 46.79 12 gam 21.29 50 partDSA 47.12 13 kernelpls 21.91 51 rf 48.32 14 rpart 23.20 52 bag 49.00 15 rqlasso 23.73 53 qrnn 49.79 16 lm 23.74 54 bagEarth 50.17 17 earth 24.00 55 BstLm 51.11 18 lars 24.05 56 enpls.fs 51.17 19 elm 24.55 57 xgbLinear 52.61 20 simpls 24.92 58 randomGLM 52.91 21 spls 25.55 59 spikeslab 53.08 22 glmboost 26.09 60 mlpWeightDecay 53.20 23 ctree2 26.11 61 rbf 53.27 24 nnls 26.33 62 bartMachine 54.52 25 extraTrees 26.89 63 RRF 54.68 26 gaussprPoly 28.00 64 mlpWeightDecayML 54.74 27 rlm 28.23 65 qrf 56.00 28 bdk 28.80 66 krlsRadial 57.35 29 relaxo 30.32 67 evtree 57.64 30 ppr 30.58 68 bstTree 57.79 31 svmRadial 31.35 69 cforest 58.62 32 gaussprRadial 31.62 70 xgbTree 65.73 33 lasso 31.95 71 SBC 67.73 34 bayesglm 34.05 72 nodeHarvest 67.83 35 bstSm 35.38 73 grnn 69.58 36 pcaNNet 35.47 74 elm-kernel 72.97 37 gbm 35.56 75 dlkeras 73.52 38 blackboost 35.71 76 svr 73.74 Table 3.8: Friedman rank of the elapsed time spent by the regressors over all the datasets. extraTrees, whose positions are 5 and 25 in the correlation and time ranks respectively. In fact, extraTrees spent an average time of 7.17 seg., while the fastest regressor (dnn) spends 3.21 s. The cubist regressor provides an even better correlation rank (about 12), but at the cost of a much higher time rank (about 35) compared to extraTrees. There is a group of three regressors (rf, RRF and bstTree) with similar correlation rank (about 15) and medium-high 3.3. Results and discussion 63 10 12 14 16 18 20 22 24 26 28 30 10 20 30 40 50 60 70 80 elm−kernel svr cubist extraTrees rf RRF bstTree gbm penalized avNNet qrf bagEarth brnn xgbTree M5 earth cforest dlkeras gaussprPoly krlsRadial Correlation Friedman rank Time Friedman rank Figure 3.7: Friedman rank of the time (vertical axis) against the Friedman rank of the correlation for the 20 best regressors in Table 3.3. time rank (about 50-60). Another group, composed by gbm, penalized, avNNet, qrf, bagEarth and brnn) is in the middle of both ranks (correlation about 17-24, time about 35-60). Finally, there are other two groups: the first one, composed by M5, earth and gaussprPoly, has low time and correlation (lower right corner); the second group, composed by xgbTree, cforest, krlsRadial and dlkeras, is placed in the upper right corner, with high time and low correlation. Considering globally the results of the current comparative, the best results are provided by the neural networks and support vector machines, specifically the extreme learning machine (position 1) and the support vector regression (position 2), both with Gaussian kernels. Cubist and extraTrees also achieve average correlations above 0.78, and can be considered also as very competitive. The extraTrees provides the best trade-off between correlation and elapsed time, being the third fastest regressor among the 20 more accurate ones. Other regressor achieve average correlations about 0.77, such as bstTree, a gradient boosting ensemble of regression trees and penalized linear regression, or about 0.75, such as the bagging ensemble of multivariate adaptive regression splines (MARS) and the regression tree model (M5). None ot the two implementations of deep learning neural networks achieves good results, but the Keras module of Python provides an acceptable performance (about 0.72). Finally, other re- 64 Chapter 3. Application of regression methods to general datasets gressors as the Bayesian regularized neural network (brnn), Gaussian processes, least squares and partial least squares, ridge and lasso regression, quantile methods, linear regression, principal component analysis and generalized additive models achieve correlations between 0.70 and 0.65, and can be considered not so competitive. CHAPTER 4 APPLICATION OF REGRESSION METHODS TO AGRICULTURAL SOIL DATA After the comparison of regression methods over the UCI machine learning repository datasets developed in chapter 3, in the current chapter we will apply them to the prediction of soil parameters described in chapter 2. An automatic prediction of soil parameters should be accurate enough to be used instead of a direct measurement of these parameters. Most of the literature about prediction of soil parameters uses the concept of pedotransfer function (PTF): a function or algorithm which uses soil measurements to predict or estimate certain soil parameters whose measurement is time-consuming or expensive [5]. The PTF can be formulated using data mining, exploration and machine learning regression methods. Examples of PTF were developed by [40] and [66] for prediction of soil total nitrogen using global soil data and water retention of soil, respectively. The latter paper demonstrated that SVM outperformed ANN for this task. Multiple-linear regression and ANN were used [77] to predict soil hydraulic parameters such as: field capacity, permanent wilting point, available water capacity, saturated hydraulic conductivity and soil water content. The inputs were basic soil properties such as sand, silt, clay, bulk density and number of pores of various diameters. The multiple-linear regression method outperformed ANN, although the difference was not statistically significant. The performance of PTF is limited by the accuracy of the prediction methods, and by spatial and temporal variations of soil parameters. Several studies [58] focused on PTFs for the determination of water retention, and the saturated and unsaturated hydraulic conductivity in order to solve the groundwater problem, properties which are expensive and difficult to measure. 66 Chapter 4. Application of regression methods to agricultural soil data Extended nonlinear regression, MLR and ANN were used [79] to estimate water-retention PTFs on Australian soil datasets. The ROSETTA software1is a user friendly tool [117] to access five PTFs to estimate hydraulic properties. Likewise, the soil inference systems (SINFERS) allows to select the PTFs with the minimum variance using knowledge rules [92]. Consequently, research about PTFs is being demanded because they are used in all branches of soil science for describing any mathematical relationships among soil properties, and they allow to estimate missing soil parameters [91]. After reviewing the literature about PTFs, our objective is to use regression techniques as PTFs that predict village-wise soil fertility indices for OC,K2O,P2O5,Fe,Mn, and Zn, as well as soil pH and nutrients N2O,P2O5and K2O, using data from the Indian region of Marathwada. Prediction of soil and crop type can not be included into this regression framework because the outputs are discrete, so they can only be treated as classification problems. K2O(kg/ha) Zn (PPM) Low <0.6 108 Medium 0.6-1.2 108-280 High >1.2 280 Table 4.1: Intervals defined by the Indian Government [23] for the calculation of village-wise soil fertility indices of K2Oand Zn nutrients [84, 63]. 4.1 Regression problems Good farm practice targets to maintain the several soil parameters that are responsible to optimize the yields of crops in Eco-friendly ways. There is need of sustainable land management practices for maintaining yield potential of agricultural crops. The soils of Marathwada are intensively cultivated for crop production by introducing novel practices. However, application of heavy doses of chemicals fertilizers deteriorates the soil health. A major factor for soil productivity is fertility, which primarily deals with ability of soil to supply nutrients to plants. Fertility of agricultural soil is depleting due to intensive cultivation practices and inadequate use of chemical fertilizers. To solve soil problems, there is a need of knowledge about soil physical and chemical status. The village-wise soil fertility indices for OC,K2O, P2O5,Fe,Mn and Zn are not only helpful to choose correct fertilizer doses, but also to know about inherent excess and deficiency in them. The predicted values can be used to balance soil 1https://www.ars.usda.gov/pacific-west-area/riverside-ca/us-salinity-laboratory/docs/rosetta-model 4.1. Regression problems 67 nutrients up to critical level in soil. Inspired by the previous concept, we apply the collection of 76 regression techniques described in chapter 3, using the experimental setup described in section 3.2, to predict village-wise soil fertility indices of major (OC,P2O5,K2O) and micro (Fe,Mn,Zn) nutrients, denoted as OC-F, K2O-F, P2O5-F, Fe-F, Mn-F and Zn-F, which were already considered in chapter 2 for classification. We also include two additional village-wise soil fertility indices (for K2Oand Zn), which are missing in chapter 2 because the available patterns belong just to one of the classes. For the eight previous regression problems, there are 372 patterns with 11 inputs: EC,OC,N2O,P2O5,K2O,SO4,Cu,Fe,Mn,Zn and Boron. In order to calculate the K2Oand Zn village-wise soil fertility indices (see chapter 2), we used the formula defined in [105] and the limits listed in Table 4.1. Our study does not consider soil N2Oand Cu village-wise soil fertility indices because the available data only include patterns of one soil nutrient level. Fertilizer aCrop bNutrient 3.31 Bajra(R) 0.38 Uria 3.38 4.11 N2O 1.65 0.068 13.1 Cotton(R), (I) 0.75 Super-phosphate 6.83 2.84 P2O5 8.57 0.18 6.86 Soybean(R) 0.68 Murate of potash 6.17 4.46 K2O 3.96 0.13 Table 4.2: Values of parameters aand bfor the prediction of different soil nutrients and crops. In order to maintain the soil quality and to get an optimal crop yield, it is necessary to apply a suitable type and amount of fertilizer which avoid excess of deficiency of soil nutrients N2O,P2O5and K2O. The prediction of the levels of these nutrients is therefore useful to calculate the right fertilizer dose, without the need of direct measurements in soil testing laboratory. This prediction will be also useful for developing field applications which recommend amounts of specific fertilizers based on the predicted nutrients levels and the target crop. Our available data contain the levels of the three nutrients for years 2011 to 2015. In our study, we use the formula developed by the Mahatma Phule Agriculture University (India) to recommend the amount of nutrient fertilizer adequate for a particular crop. The amount Aof fertilizer to apply is calculated as A=aE −bP, being Ethe expected crop production (an input data) and Pthe predicted nutrient level. The values Aand Pare measured in kg/ha, while Eis measured in quintal/ha. Depending on the nutrient and crop, the parameters aand 68 Chapter 4. Application of regression methods to agricultural soil data btake different values (see Table 4.2) and the fertilizer is different. The values of aand bin Table 4.2 are set for black soil, being different for other soil types. For nutrients N2O,P2O5 and K2O, the corresponding fertilizer is different: uria (which contains 46% of N2O), super phosphate (46% of P2O5) and muriate of potash (60-62% of K2O), respectively. The following subsections discuss the results achieved by the collection of regressors for the village-wise soil fertility indices, for the soil nutrients and for soil pH. The last section discusses globally the results. 4.2 Prediction of OC,P2O5,K2O,Fe,Mn and Zn village-wise fertility indices The Table 4.3 reports the RMSE for the OC village-wise soil fertility index, sorted by increasing values. The extremely randomized regression trees (extraTrees) achieves the lowest RMSE (0.564), although regularized random forests (RRF), standard random forest (rf), random forest with feature selection (Boruta), gradient boosting of regressor trees (bstTree) and gradient boosted machine (gbm) also achieve values about 0.57. The remaining regressors are about 0.6, including elm-kernel and svr, which are the two best regressors in chapter 3. On the opposite side of the table, the worst regressor is icr, with RMSE=1.17. These values can be considered high, because the output is normalized to have zero mean and standard deviation one, so the RMSE should be low (e.g. under 0.5). The Figure 4.1 (upper panel) shows the correlations achieved by the 20 best regressors for soil OC-F problem. The four best regressors are the same as in Table 4.3, with changes in their positions: Boruta, RRF, rf and extraTrees achieved high correlation above 0.83 for the soil OC-F problem, while bstTree and gbm are about 0.82. However cforest, elm-kernel and kknn and 6 more regressors achieved correlations below 0.77. The lower panel of Figure 4.1 presents the scatter plot of the best regressor, Boruta, with RMSE=0.573 and correlation 0.8355. The results for P2O5village-wise soil fertility index are reported by the Table 4.4. Again, the extraTrees achieves the best RMSE alongside with RRF, gbm, quantile random forest (qrf), Boruta, svr and rf (about 0.63). The elm-kernel works even worse (0.683), and penalized, which also achieves good results in chapter 3, is about 0.75. The gaussprPoly is the worst regressor, with RMSE = 8.41, alongside with ridge, gam, glm, lasso and lm, with RMSE about 1.93, which also correspond to poor performances. The regressors from glmnet to gaussprPoly got the lowest performances with RMSE values above 1. 4.2. Prediction of OC , P2O5 , K2O , Fe , Mn and Zn village-wise fertility indices 69 Regressor RMSE Regressor RMSE Regressor RMSE Regressor RMSE extraTrees 0.564 treebag 0.641 lars 0.709 earth 0.772 RRF 0.57 bstSm 0.645 rqlasso 0.71 gam 0.772 rf 0.571 brnn 0.65 glmStepAIC 0.711 glm 0.772 Boruta 0.573 grnn 0.666 rqnc 0.718 lasso 0.772 bstTree 0.573 penalized 0.667 spls 0.722 enpls.fs 0.772 gbm 0.578 blackboost 0.671 foba 0.728 lm 0.772 bartMachine 0.584 rbf 0.677 dlkeras 0.729 randomGLM 0.776 qrf 0.592 ppr 0.677 xgbLinear 0.732 evtree 0.795 cubist 0.599 xgbTree 0.684 elm 0.735 superpc 0.811 nodeHarvest 0.601 simpls 0.689 plsRglm 0.748 mlpWeightDecay 0.836 svr 0.606 kernelpls 0.689 rlm 0.75 BstLm 0.838 krlsRadial 0.62 avNNet 0.689 gaussprPoly 0.759 nnls 0.845 svmRadial 0.623 SBC 0.689 bdk 0.759 pcaNNet 0.848 gamboost 0.626 relaxo 0.691 ridge 0.765 mlpWeightDecayML 0.975 elm-kernel 0.629 bagEarth 0.697 M5 0.766 dnn 1 cforest 0.632 qrnn 0.699 ctree2 0.769 pcr 1 rvmRadial 0.632 bag 0.7 bayesglm 0.769 glmnet 1.01 gaussprRadial 0.633 glmboost 0.702 gaussprLinear 0.771 partDSA 1.01 kknn 0.633 spikeslab 0.706 rpart 0.772 icr 1.17 Table 4.3: RMSE for the prediction of OC village-wise soil fertility index. Regressor RMSE. Regressor RMSE. Regressor RMSE. Regressor RMSE. extraTrees 0.631 gaussprRadial 0.7 randomGLM 0.814 dnn 1.01 RRF 0.634 brnn 0.704 nnls 0.814 glmboost 1.01 gbm 0.634 nodeHarvest 0.707 rqlasso 0.825 partDSA 1.01 qrf 0.635 blackboost 0.723 rqnc 0.837 dlkeras 1.08 Boruta 0.635 M5 0.741 foba 0.859 enpls.fs 1.21 svr 0.636 bag 0.744 gamboost 0.864 earth 1.27 rf 0.637 SBC 0.745 evtree 0.865 bstSm 1.33 cubist 0.643 xgbTree 0.746 BstLm 0.88 spikeslab 1.58 bstTree 0.647 rbf 0.751 bdk 0.893 rlm 1.61 bartMachine 0.654 penalized 0.752 bagEarth 0.899 xgbLinear 1.7 svmRadial 0.659 ppr 0.756 mlpWeightDecay 0.899 glmStepAIC 1.73 krlsRadial 0.672 elm 0.779 pcaNNet 0.904 ridge 1.76 elm-kernel 0.683 relaxo 0.782 icr 0.929 bayesglm 1.86 cforest 0.683 spls 0.786 superpc 0.939 gaussprLinear 1.91 rvmRadial 0.684 simpls 0.786 mlpWeightDecayML 0.952 gam 1.93 avNNet 0.686 kernelpls 0.786 pcr 0.958 glm 1.93 kknn 0.687 qrnn 0.797 lars 0.964 lasso 1.93 grnn 0.691 rpart 0.806 plsRglm 0.99 lm 1.93 treebag 0.691 ctree2 0.81 glmnet 1 gaussprPoly 8.41 Table 4.4: RMSE for the prediction of P2O5village-wise soil fertility index. The Table 4.5 reports the RMSE for the prediction of K2Ovillage-wise soil fertility index. The best RMSE (0.768 using extraTrees) in this table is higher than the two previous indices, so this dataset is more difficult. Boruta achieves RMSE (0.78), while RRF, rf and qrf are about 76 Chapter 4. Application of regression methods to agricultural soil data 0.47 0.475 0.48 0.485 0.49 0.495 0.5 0.505 0.51 0.515 penalized cforest svr extraTrees rf RRF bartMachine qrnn elm−kernel treebag rbf nodeHarvest svmRadial qrf brnn superpc nnls relaxo avNNet kernelpls Correlation −2 −1 0 1 2 3 −2 −1 0 1 2 3 Desired outputs Real outputs Figure 4.3: Correlation (upper panel) and scatter plot of penalized (lower panel) for the prediction of N2Olevel in the state of Marathwada. and krlsRadial (0.829) and RRF (0.83). The highest RMSE (1.72) is achieved by independent component regression (icr). 4.4 Prediction of soil pH ExtraTrees achieves the best RMSE (0.711) for the prediction of soil pH (see Table 4.12), followed with minor differences by cubist, RRF and rf (about 0.74), Boruta and qrf (about 4.4. Prediction of soil pH 77 Regressor RMSE. Regressor RMSE. Regressor RMSE. Regressor RMSE. elm-kernel 0.818 pcaNNet 0.876 foba 0.93 ppr 1.09 gaussprRadial 0.818 M5 0.878 rqnc 0.931 gamboost 1.12 svr 0.819 mlpWeightDecay 0.883 SBC 0.936 bdk 1.21 extraTrees 0.827 relaxo 0.884 blackboost 0.943 qrnn 1.24 rf 0.829 rpart 0.889 gaussprPoly 0.949 rlm 1.24 krlsRadial 0.829 bagEarth 0.899 pcr 0.949 enpls.fs 1.29 RRF 0.831 brnn 0.899 spikeslab 0.952 glmStepAIC 1.29 kknn 0.835 earth 0.901 ctree2 0.953 spls 1.3 treebag 0.837 simpls 0.904 glmnet 0.954 plsRglm 1.3 qrf 0.841 kernelpls 0.904 mlpWeightDecayML 0.954 xgbLinear 1.36 nodeHarvest 0.842 BstLm 0.904 partDSA 0.955 ridge 1.36 svmRadial 0.843 grnn 0.907 dnn 0.956 bayesglm 1.44 bartMachine 0.847 elm 0.907 evtree 0.962 gaussprLinear 1.47 Boruta 0.848 glmboost 0.908 cubist 0.968 gam 1.48 penalized 0.857 rbf 0.911 avNNet 0.97 glm 1.48 rvmRadial 0.859 nnls 0.913 lars 0.972 lasso 1.48 bstTree 0.865 bag 0.916 superpc 1 lm 1.48 cforest 0.867 xgbTree 0.918 dlkeras 1.02 bstSm 1.56 gbm 0.872 rqlasso 0.925 randomGLM 1.05 icr 1.72 Table 4.11: RMSE values for the prediction of K2Onutrient level. Regressor RMSE Regressor RMSE Regressor RMSE Regressor RMSE extraTrees 0.711 svmRadial 0.812 gamboost 0.894 superpc 1.03 cubist 0.746 M5 0.815 bstSm 0.903 bdk 1.06 RRF 0.747 rvmRadial 0.817 relaxo 0.904 BstLm 1.13 rf 0.748 evtree 0.823 simpls 0.908 glmboost 1.15 Boruta 0.749 earth 0.829 kernelpls 0.908 randomGLM 1.16 qrf 0.751 grnn 0.831 pcaNNet 0.913 enpls.fs 1.28 bartMachine 0.767 brnn 0.848 rqlasso 0.94 spikeslab 1.43 nodeHarvest 0.774 bag 0.85 spls 0.941 ridge 1.45 bstTree 0.776 bagEarth 0.859 rqnc 0.953 xgbLinear 1.55 gbm 0.78 SBC 0.861 mlpWeightDecay 0.953 rlm 1.55 xgbTree 0.782 rpart 0.862 lars 0.956 gaussprLinear 1.55 kknn 0.782 ppr 0.864 gaussprPoly 0.959 gam 1.56 treebag 0.783 rbf 0.868 pcr 0.976 glm 1.56 cforest 0.785 blackboost 0.871 mlpWeightDecayML 0.982 lasso 1.56 krlsRadial 0.797 qrnn 0.874 dnn 0.986 bayesglm 1.56 avNNet 0.798 elm 0.876 plsRglm 0.986 lm 1.56 gaussprRadial 0.801 ctree2 0.876 glmnet 0.99 nnls 1.57 elm-kernel 0.802 penalized 0.885 partDSA 0.99 glmStepAIC 1.59 svr 0.811 dlkeras 0.887 icr 0.991 foba 1.66 Table 4.12: RMSE values for the prediction of soil pH. 0.75). The remaining regressors are already above 0.77. The Figure 4.4 (upper panel) plots the low correlation values achieved by the 20 best regressors, started by extraTrees (0.6945, see its scatter plot on the right panel) and followed by a group about 0.65 (cubist, RRF, rf, 78 Chapter 4. Application of regression methods to agricultural soil data 0.58 0.6 0.62 0.64 0.66 0.68 extraTrees cubist RRF rf Boruta qrf bartMachine nodeHarvest bstTree gbm xgbTree cforest kknn treebag avNNet krlsRadial gaussprRadial elm−kernel svr svmRadial Correlation −4 −3 −2 −1 0 1 2 −4 −3 −2 −1 0 1 2 Desired outputs Real outputs Figure 4.4: Correlation of the twenty best regressors (upper panel) and scatter plot of extraTrees (lower panel), which achieves the best correlation for the prediction of soil pH. Boruta and qrf). The right panel of Figure 4.4 shows the scatter plot of extraTrees for the prediction of soil pH. 4.5 Global discussion Considering the results over all the soil datasets (see Tables 4.3-4.12), extraTrees achieved the best RMSE for 7 of 10 soil problems (OC,P2O5,K2O,Fe,Mn and Zn village-wise soil 4.5. Global discussion 79 fertility indices, and soil pH). Besides, penalized linear regression, random forest and elmkernel are the best regressors for soil N2O,P2O5and K2Onutrients, respectively. Perhaps the most relevant conclusion that we can draw from the results is that, besides being extraTrees the best regressors for more than a half of the datasets, other four regressors of the random forest family (rf, RRF, Boruta and qrf) are among the best five regressors for almost all the soil datasets. Therefore, this family can be considered the best for these datasets, confirming the good result of both random forests in the soil classification problems (see chapter 2). The svr and two gradient boosting ensembles (bstTree and gbm) are also among the ten best regressors for almost all the fertility indices, as well as for P2O5and pH. The M5 rule with nearest neighbors (cubist) is among the ten bests regressors for 6 datasets, being the second best for Zn village-wise soil fertility index and prediction of pH. OC-F P2O5-F K2O-F Fe-F Mn-F Regressor Corr. Regressor Corr. Regressor Corr. Regressor Corr. Regressor Corr. Boruta 0.835 extraTrees 0.776 extraTrees 0.631 extraTrees 0.816 extraTrees 0.758 RRF 0.834 RRF 0.775 Boruta 0.619 rf 0.812 RRF 0.741 rf 0.834 qrf 0.774 RRF 0.607 RRF 0.810 bstTree 0.739 extraTrees 0.832 gbm 0.774 rf 0.601 Boruta 0.807 rf 0.739 bstTree 0.821 Boruta 0.774 qrf 0.597 bstTree 0.789 Boruta 0.737 gbm 0.817 rf 0.772 elm-kernel 0.586 gbm 0.789 gbm 0.737 bartMachine 0.815 svr 0.772 bstTree 0.585 qrf 0.788 svr 0.734 nodeHarvest 0.815 cubist 0.770 gaussprRadial 0.581 cubist 0.775 cubist 0.730 qrf 0.811 bstTree 0.765 nodeHarvest 0.579 nodeHarvest 0.767 svmRadial 0.724 cubist 0.806 bartMachine 0.757 svr 0.578 bartMachine 0.762 qrf 0.723 Zn-F N2O P2O5K2O pH Regressor Corr. Regressor Corr. Regressor Corr. Regressor Corr. Regressor Corr. extraTrees 0.841 penalized 0.519 qrf 0.487 gaussprRadial 0.517 extraTrees 0.694 cubist 0.794 cforest 0.502 rf 0.476 svr 0.516 cubist 0.657 avNNet 0.790 svr 0.495 RRF 0.473 elm-kernel 0.514 RRF 0.655 svr 0.784 extraTrees 0.493 extraTrees 0.470 krlsRadial 0.511 rf 0.655 gbm 0.783 rf 0.489 gbm 0.464 extraTrees 0.500 Boruta 0.653 qrf 0.779 RRF 0.488 bartMachine 0.460 rf 0.495 qrf 0.651 rf 0.778 bartMachine 0.488 bstTree 0.458 RRF 0.492 bartMachine 0.630 Boruta 0.775 qrnn 0.486 treebag 0.446 kknn 0.486 nodeHarvest 0.626 RRF 0.774 elm-kernel 0.485 penalized 0.444 qrf 0.483 bstTree 0.618 krlsRadial 0.772 treebag 0.481 nodeHarvest 0.434 treebag 0.477 gbm 0.614 Table 4.13: Ten bests correlations for the prediction of each soil dataset. The Table 4.13 reports the ten best regressors according to correlation coefficient for the ten soil datasets. ExtraTrees achieves the best correlations for six of ten datasets (P2O5-F, K2O-F, Fe-F, Mn-F, Zn-F and pH), being the 4th or 5th in the remaining four datasets. The best regressors for the remaining datasets are Boruta, penalized, qrf and gaussprRadial for 80 Chapter 4. Application of regression methods to agricultural soil data RMSE rank Correlation rank Pos. Regressor Rank Regressor Rank Avg. correl. 1 extraTrees 1.7 extraTrees 2.3 0.68132 2 RRF 4.2 RRF 4.0 0.66519 3 rf 4.7 rf 4.3 0.66520 4 qrf 7.5 qrf 7.0 0.65677 5 svr 8.5 Boruta 8.4 0.64915 6 bstTree 9.2 svr 8.5 0.64234 7 gbm 9.5 gbm 9.7 0.64419 8 bartMachine 9.5 bstTree 9.7 0.64580 9 Boruta 9.8 bartMachine 9.9 0.64228 10 elm-kernel 12.9 nodeHarvest 12.9 0.62539 11 nodeHarvest 13.0 svmRadial 13.3 0.62436 12 cforest 14.1 krlsRadial 13.5 0.62748 13 svmRadial 14.1 cubist 13.5 0.62658 14 treebag 14.3 elm-kernel 14.0 0.62206 15 krlsRadial 14.7 treebag 15.6 0.61491 16 gaussprRadial 16.6 cforest 15.8 0.60960 17 cubist 17.3 gaussprRadial 16.4 0.61405 18 kknn 17.9 kknn 18.4 0.60867 19 rvmRadial 19.3 avNNet 19.3 0.59096 20 penalized 23.1 rvmRadial 19.4 0.59719 21 brnn 24.0 penalized 25.3 0.54830 22 avNNet 24.9 brnn 25.6 0.55449 23 grnn 25.3 grnn 27.1 0.56332 24 xgbTree 27.5 xgbTree 29.0 0.54622 25 bag 28.1 rbf 29.6 0.51867 26 rbf 28.1 ppr 30.5 0.52761 27 relaxo 31.3 bag 31.2 0.52418 28 kernelpls 31.8 qrnn 32.3 0.50870 29 blackboost 32.2 SBC 32.5 0.51816 30 simpls 32.8 kernelpls 34.3 0.48163 31 SBC 33.2 relaxo 34.3 0.48478 32 qrnn 33.5 blackboost 34.4 0.49932 33 ppr 35.7 pcaNNet 34.4 0.52631 34 M5 37.2 dlkeras 34.9 0.49958 35 bagEarth 37.4 simpls 35.3 0.48163 36 earth 38.5 mlpWeightDecay 36.2 0.49583 37 elm 39.3 bagEarth 37.5 0.48584 38 rpart 39.8 rpart 37.9 0.50542 Table 4.14: Friedman rank of the RMSE (left) and of the correlation (right), and average correlation. Continued in Table 4.15. OC-F, N2O,P2O5and K2O, respectively. The random forests achieve very good results: RRF is the 2nd or 3rd best for 7 datasets; rf is the 4th regressor in 4 datasets; Boruta is among the 5 bests for 6 datasets; and qrf is among the 10 bests for 9 of 10 datasets. The boosting regressors gbm and bstTree are among the 10 bests for 7 of 10 datasets; the svr, bartMachine and cubist 4.5. Global discussion 81 RMSE rank Correlation rank Order Regressor Rank Regressor Rank Avg. correl. 39 pcaNNet 41.4 M5 39.8 0.47795 40 rqlasso 42.3 earth 40.0 0.45912 41 evtree 42.4 elm 41.4 0.45311 42 mlpWeightDecay 42.6 ctree2 42.3 0.46354 43 ctree2 43.5 bdk 42.9 0.45099 44 rqnc 46.2 rqlasso 44.4 0.41581 45 lars 46.9 evtree 44.8 0.44705 46 superpc 47.2 rqnc 47.3 0.39901 47 nnls 47.6 spls 48.2 0.36216 48 gamboost 48.4 BstLm 49.1 0.35185 49 glmboost 48.6 gamboost 49.6 0.33209 50 BstLm 48.8 glmboost 50.3 0.33887 51 dlkeras 49.8 bstSm 50.4 0.31987 52 bstSm 50.3 lars 51.1 0.34588 53 spls 50.8 nnls 51.6 0.35760 54 pcr 51.4 superpc 51.9 0.32258 55 mlpWeightDecayML 52.1 spikeslab 54.2 0.29352 56 foba 53.3 foba 54.5 0.33074 57 bdk 53.3 plsRglm 55.2 0.31830 58 gaussprPoly 53.4 randomGLM 57.7 0.28005 59 dnn 54.4 glmStepAIC 58.2 0.25896 60 spikeslab 54.8 bayesglm 58.8 0.24860 61 plsRglm 58.4 ridge 59.0 0.25092 62 icr 59.3 gaussprPoly 59.4 0.26222 63 glmStepAIC 60.2 enpls.fs 59.5 0.27277 64 randomGLM 60.6 xgbLinear 59.7 0.22684 65 rlm 61.4 rlm 60.4 0.25463 66 xgbLinear 61.4 gaussprLinear 61.0 0.24405 67 enpls.fs 61.8 glm 61.2 0.24355 68 bayesglm 63.9 lm 62.2 0.24355 69 ridge 64.8 pcr 63.0 0.18370 70 gaussprLinear 65.6 lasso 63.2 0.24355 71 glm 66.4 icr 63.9 0.17313 72 lm 67.4 gam 64.2 0.24355 73 lasso 68.4 mlpWeightDecayML 65.0 0.18838 74 partDSA 69.0 dnn 71.9 0.06922 75 gam 69.4 glmnet 75.2 -0.03778 76 glmnet 76.0 partDSA 75.3 -0.04306 Table 4.15: Continuation of Table 4.14. for 6 of them; and nodeHarvest (5 datasets). Considering the best correlation values, they only overcome 0.8 for datasets OC-F, Fe-F and Zn-F, being between 0.6 and 0.8 in four datasets (P2O5-F, K2O-F, Mn-F and pH), and below 0.6 for N2O,P2O5and K2Onutrients. The low correlations for these three nutrients confirms the higher difficulty of these datasets, where the classification experiments (see chapter 2) already achieved Cohen κ lower than the village- 82 Chapter 4. Application of regression methods to agricultural soil data wise fertility indices, crop, soil type and pH. For some datasets the difference between the best correlation and the following values is relatively high: K2O-F, where extraTrees and Boruta achieve 0.631 and 0.619, respectively (difference 0.012), Zn-F, with difference 0.097 between extraTrees (0.891) and cubist (0.794); N2O, difference 0.017 between penalized and cforest; and pH, with difference 0.037 between extraTrees and cubist. Pos. Regressor p-value Pos. Regressor p-value Pos. Regressor p-value 1 extraTrees — 27 qrnn 0.0451546 53 randomGLM 0.00131494 2 rf 0.73373 28 ppr 0.0376353 54 lars 0.000768539 3 RRF 0.62318 29 bag 0.0376353 55 glmStepAIC 0.000768539 4 qrf 0.57075 30 blackboost 0.031209 56 superpc 0.000768539 5 svr 0.52052 31 pcaNNet 0.031209 57 spikeslab 0.000768539 6 Boruta 0.47268 32 mlpWeightDecay 0.0257481 58 lm 0.00058284 7 gbm 0.42735 33 evtree 0.0257481 59 enpls.fs 0.00058284 8 cubist 0.38467 34 rpart 0.0211339 60 rlm 0.00058284 9 bstTree 0.38467 35 bagEarth 0.0211339 61 ridge 0.00058284 10 elm-kernel 0.34470 36 earth 0.0172575 62 lasso 0.00058284 11 krlsRadial 0.34470 37 rbf 0.0172575 63 gam 0.00058284 12 gaussprRadial 0.30749 38 ctree2 0.0172575 64 gaussprLinear 0.00058284 13 svmRadial 0.30749 39 M5 0.0140193 65 bayesglm 0.00058284 14 bartMachine 0.30749 40 bdk 0.0140193 66 gaussprPoly 0.00058284 15 cforest 0.27304 41 relaxo 0.0091085 67 plsRglm 0.00058284 16 kknn 0.24132 42 elm 0.0091085 68 BstLm 0.00058284 17 nodeHarvest 0.24132 43 rqnc 0.00728456 69 glm 0.00058284 18 rvmRadial 0.18588 44 gamboost 0.00728456 70 xgbLinear 0.00058284 19 treebag 0.18588 45 kernelpls 0.00579536 71 glmnet 0.000182672 20 avNNet 0.16197 46 simpls 0.00579536 72 dnn 0.000182672 21 grnn 0.12122 47 bstSm 0.00458639 73 pcr 0.000182672 22 SBC 0.07566 48 nnls 0.00458639 74 icr 0.000182672 23 brnn 0.07566 49 rqlasso 0.00458639 75 partDSA 0.000182672 24 xgbTree 0.05390 50 foba 0.00282727 76 mlpWeightDecayML 0.000182672 25 penalized 0.05390 51 spls 0.00170625 26 dlkeras 0.05390 52 glmboost 0.00170625 Table 4.16: List of the p-values achieved by the Wilcoxon signed rank test comparing the correlations of extraTrees and the remaining 75 regressors over the 10 soil regression problems. The Tables 4.14 and 4.15 report the Friedman ranks for the RMSE and correlation, and the average correlation, for all the regressors over all the ten soil datasets. The extraTrees regressor achieves the first position, with a rank of 1.7 and 2.3 for RMSE and correlation respectively. This means that, in average over the ten datasets, extraTrees is between positions 1-2 for RMSE and 2-3 for correlation. However, the highest average correlation (0.681), achieved by extraTrees, reflects that the prediction is not very accurate, because an accurate prediction would require correlations about 0.9-0.95. Four regressors of the random forest family (RRF, rf, qrf, and Boruta) are placed in the first positions, alongside with the svr and 4.5. Global discussion 83 the two gradient boosting ensembles (bstTree and gbm). The elm-kernel is in position 10, with average correlation about 0.625, which is far from the best results. Other regressors with good results in some soil datasets (e.g. cubist) or in the UCI datasets (see chapter 3), e.g. penalized and avNNet, are in positions 10-20 on this ranking. The last positions of the ranking (Table 4.15) are for lm (linear regression), bstSm (gradient boosting with smoothing splines), glm (generalized linear models), icr (independent component regression), randomGLM (boosting ensemble of GLM), glmStepAIC (GLM with stepwise feature selection and the Akaike information criterion) and lasso (regression by least absolute shrinkage and selection operator). 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 80 extraTrees RRF rf qrf Boruta svr gbm bstTree bartMachine nodeHarvest svmRadial krlsRadial cubist elm−kernel treebag cforest gaussprRadial kknn avNNet rvmRadial Correlation Friedman rank Time Friedman rank Figure 4.5: Friedman rank of the time (vertical axis) against the Friedman rank of the correlation (horizontal axis) for the 20 best regressors over the ten soil data sets. The Table 4.16 reports the p-values for a Wilcoxon signed rank test [132] comparing the correlations achieved by extraTrees, which is the best regressor on the soil datasets according both to RMSE and correlation, to the correlations of the remaining regressors, sorted by decreasing order. The value in bold corresponds to the regressor (qrnn, position 27 of 76) from which the difference with respect to extraTrees is statistically significant for a 5%-confidence level (i.e., p<0.05). Since the difference extraTrees and the first regressors in the list is only statistically significant after position 27 (qrnn) of 76 regressors, it is clear that differences are not very high, in fact even the best regressors do not exhibit an excellent performance. 84 Chapter 4. Application of regression methods to agricultural soil data We also measured the elapsed time for each regressor and dataset. Since a simple averaging of times over datasets is not statistically acceptable, because times vary in different ranges for each dataset, we created a Friedman rank of the elapsed times for each regressor over all datasets. The Figure 4.5 plots the time against correlation (both in terms of Friedman rank) for the 20 best regressors in the correlation rank of Table 4.14. The figure shows that the best regressor (extraTrees) in terms of correlation, because it is placed on the left end of the plot (correlation rank about 2), is also the second fastest one because it is placed on the lower end (time rank 23), being only slower than kknn (time rank about 8), which however works much worse (correlation rank about 18). Among the other best regressors in chapter 3, the rf and RRF exhibit slightly lower correlation than extraTrees (they are on its right), but they are much slower (time rank above 60). BstTree and gbm are slightly slower (upper time rank) than extraTrees, but their correlation is much worse (correlation rank about 10). Finally, the elm-kernel and svr are much more slower than extraTrees (time rank above 55) with much lower correlation (ranks about 8 and 14, respectively). CHAPTER 5 CONCLUSIONS Agriculture is a major sector in the Indian economy, which is affected by changing trends in temperature and rainfall, insufficient water, agriculture practices and nutrient deficiencies. Adequate soil parameters and proper application of fertilizers may help to attenuate these problems. The current research supports the Indian Government to make decisions about improving soil quality and crop production. The soil quality depends on its type and pH, village-wise fertility indices of OC,P2O5,Mn and Fe, and on the selected crop. Thus, an automatic prediction of their values from measurements of N2O,P2O5,K2O,SO4and EC, among others, would reduce the cost of the chemical analysis and save time for specialized technicians. The prediction of levels for the soil nutrients N2O,P2O5and K2Owould also be very useful for the recommendation of suitable fertilizers. The work developed in this PhD. Thesis is oriented to use machine learning techniques to automatically predict these values for soils of the Indian state of Maharashtra. The results of this study might contribute to design agriculture strategies of the Indian Government to manage the soil fertility degradation, crop productivity and usage of fertilizers. Despite of being examples of regression problems, our first approach was to quantity the values of the magnitudes to be predicted into low, medium and high levels using thresholds defined by the Indian Government, transforming them into classification problems. We applied a wide and diverse collection of classifiers including decision trees, rule-based classifiers, bagging and boosting ensembles, random forests, neural networks, support vector machines and nearest neighbors classifiers. We achieved values of the Cohen κ about 97% and 90% for for soil classification and village-wise OC fertility index, respectively, above 85% for P2O5 92 Bibliography [36] J.H. Friedman and W. Stuetzle. Projection pursuit regression. Journal of the American Statistical Association, 76:817–823, 1981. [37] S. García, A. Fernández, A.D. Benítez, and F. Herrera. Statistical comparisons by means of non-parametric tests: A case study on genetic based machine learning. In Proceedings of the II Congreso Español de Informática (CEDI 2007). V Taller Nacional de Minería de Datos y Aprendizaje (TAMIDA), pages 95–104, 2007. [38] A. Gelman, A. Jakulin, M.G. Pittau, and Y.S. Su. A weakly informative default prior distribution for logistic and other regression models. The Annals of Applied Statistics, 2(4):1360–1383, 2009. [39] P. Geurts, D. Ernst, and L. Wehenkel. Extremely randomized trees. Machine Learning, 63(1):3–42, 2006. [40] M. J. Glendining, A. G. Dailey, D. S. Powlson, G. M. Richter, J.A. Catt, and A. P. Whitmore. Pedotransfer functions for estimating total soil nitrogen up to the global scale. European Journal of Soil Science, 62:13–22, 2011. [41] J.J. Goeman. L-1 penalized estimation in the cox proportional hazards model. Biometrical Journal, 52:70–84, 2010. [42] J.W. Groenigen, D. Huygens, P. Boeckx, Th.W. Kuyper, I.M. Lubbers, T. Rütting, and P.M. Groffman. The soil N cycle: new insights and key challenges. Soil, 1:235–256, 2015. [43] T. Grubinger, A. Zeileis, and K.P. Pfeiffer. Evtree: evolutionary learning of globally optimal classification and regression trees in R. Journal of Statistical Software, 61(1):1–29, 2014. [44] P. Gruhn, F. Goletti, and M. Yudelman. Integrated Nutrient Management, Soil Fertility, and Sustainable Agriculture: Current Issues and Future Challenges. International Food Policy Research Institute, 2000. [45] J.M. Guerrero, G. Pajares, M. Montalvo, J. Romeo, and M. Guijarro. Support vector machines for crop/weeds identification in maize fields. Expert Syst. Appl., 39:11149–11155, 2012. Bibliography 93 [46] M. Hall and E. Frank. Combining naive Bayes and decision tables. Proc. Artif. Intel. Soc. Conf., pages 318–319, 2008. [47] M. Hall, E. Frank, G.Ho., B. Pfahringer, P. Reutemann, and I.H. Witten. The Weka data mining software: an update. SIGKDD Explorations, 11:10–18, 2009. [48] T. Hastie, R. Tibshirani, and A. Buja. Flexible discriminant analysis by optimal scoring. Journal of the American Statistical Association, 89:1255–1270, 1993. [49] B. Heung, H. Chak Ho, J. Zhang, A. Knudby, C.E. Bulmer, and M.G. Schmidt. An overview and comparison of machine-learning techniques for classification purposes in digital soil mapping. Geoderma, 265:62–77, 2016. [50] M.G. Hill, P.G. Connolly, P. Reutemann, and D. Fletcher. The use of data mining to assist crop protection decisions on kiwifruit in New Zealand. Comput. Electron. Agric., 108:250–257, 2014. [51] Geoffrey E. Hinton, Simon Osindero, and Yee-Whye Teh. A fast learning algorithm for deep belief nets. Neural Comput., 18(7):1527–1554, 2006. [52] T. Hothorn, K. Hornik, and A. Zeileis. Unbiased recursive partitioning: A conditional inference framework. Journal of Computational and Graphical Statistics, 15(3):651–674, 2006. [53] G.-B. Huang, H. Zhou, X. Ding, and R. Zhang. Extreme learning machine for regression and multiclass classification. IEEE Trans. Systs., Man, and Cybern.-Part B: Cybern., 42(2):513–529, 2012. [54] P.J. Huber. Robust statistics. Wiley, 1981. [55] A. Hyvarinen and E. Oja. Independent component analysis: algorithms and applications. Neural networks, 13:411–430, 2000. [56] H. Ishwaran, J.S. Rao, and U.B. Kogalur. Spikeslab : prediction and variable selection using spike and slab regression. The R Journal, 2:68–73, 2010. [57] M.L. Jackson. Soil chemical analysis. Prentice Hall of India Pvt. Ltd., New Delhi, 1958. 94 Bibliography [58] S. K. Jain, V. P. Singh, and M. Th. V. Genuchten. Analysis of Soil Water Retention Data Using Artificial Neural Networks. Journal of Hydrologic Engineering, pages 415–420, 2004. [59] M.E. Jakubauskas, D.R. Legates, and J.H. Kastens. Crop identification using harmonic analysis of time-series AVHRR NDVI data. Comput. Electron. Agric., 37:127–139, 2002. [60] C.H. Jones. Activity of organic nitrogen as measured by the alkaline permanganate method. J. Ind. Eng. Chem., 24:438–441, 1912. [61] K. Hechenbichler and K.P. Schliep. Weighted k-nearest-neighbor techniques and ordinal classification. Technical report, Ludwig-Maximilians University Munich, 2004. [62] A. Kapelner and J. Bleich. bartMachine: machine learning with Bayesian additive regression trees. Journal of Statistical Software, 70(4):1–40, 2016. [63] J.C. Katyal and R.K. Rattan. Secondary and micronutrients: research gaps and future needs. Fertil. News, 48:9–20, 2003. [64] R. Kohavi. A study of cross-validation and bootstrap for accuracy estimation and model selection. International Joint Conference on Artificial Intelligence (IJCAI), 1995. [65] M.B. Kursa and W.R. Rudnicki. Feature selection with the Boruta package. Journal of Statistical Software, 36(11):1–13, 2010. [66] K. Lamorski, Y. Pachepsky, C. Slawin´ ski, and R. T. Walczak. Using Support Vector Machines to Develop Pedotransfer Functions for Water Retention of Soils in Poland. Soil Sci. Soc. Am. J., 72:1243–1247, 2008. [67] C.L. Lawson and R.J. Hanson. Solving least squares problems, volume 15 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), 1995. [68] Leggeit and D.P. Argyle. The DTPA-extractable iron, manganese, copper, and zinc from neutral and calcareous soils dried under different conditions. Soil Sci. Soc. Am. J., 47(3):518–522, 1983. Bibliography 95 [69] W. Liu, Z. Wang, X. Liu, N. Zeng, Y. Liu, and F.E. Alsaadi. A survey of deep neural network architectures and their applications. Neurocomputing, 234:11–26, 2017. [70] D.J.C. MacKay. Bayesian interpolation. Neural Computation, 4:415–447, 1992. [71] Mahatma Phule Agricultural University. Krishi Darshani. Parbhani, Maharashtra, India, 2016. Page 18 (in Marathi language). [72] Matlab. version 7.14 (R2012a). The MathWorks Inc., Natick, Massachusetts, 2012. [73] N. Meinshausen. Relaxed lasso. Computational Statistics and Data Analysis, pages 374–393, 2007. [74] N. Meinshausen. Node harvest. The Annals of Applied Statistics, 4(4):2049–2072, 2010. [75] W. Melssen, R. Wehrens, and L. Buydens. Supervised Kohonen networks for classification problems. Chemom. Intell. Lab. Syst., 83:99–113, 2006. [76] P. Melville and R.J. Mooney. Creating diversity in ensembles using artificial data. Inform. Fusion: special issue on diversity in multiclassifier systems, 6(1):99–111, 2004. [77] H. Merdun, O. Cinar, R. Meral, and M. Apan. Comparison of artificial neural network and regression pedotransfer functions for prediction of soil water retention and saturated hydraulic conductivity. Soil and Tillage Research, 90:108–116, 2006. [78] B.H. Mevik and H.R. Cederkvist. Mean squared error of prediction (msep) estimates for principal component regression (pcr) and partial least squares regression (plsr). Journal of Chemometrics, 18(9):422–429, 2004. [79] B. Minasny, A. B. McBratney, and K. L. Bristow. Comparison of different approaches to the development of pedotransfer functions for water-retention curves. Geoderma, 93:225–253, 1999. [80] I. Mizera and R. Koenker. Convex optimization in r. Journal of Statistical Software, 60(5):1–23, 2014. 96 Bibliography [81] M.S. Mkhabela, P. Bullock, S. Raj, S. Wang, and Y. Yang. Crop yield forecasting on the Canadian prairies using MODIS NDVI data. Agric. For. Meteorol., 151:385–393, 2011. [82] A.M. Molinaro, K.Lostritto, and M.J. van der Laan. Partdsa: deletion/substitution/addition algorithm for partitioning the covariate space in prediction. Bioinformatics, 26(10):1357–63, 2010. [83] A. Mucherino, P. Papajorgji, and P.M. Pardalos. A survey of data mining techniques applied to agriculture. Oper. Res., 9:121–140, 2009. [84] G.R. Muhr, N.P. Datta, S.N. Shankara, F. Dever, V.K. Lecy, and R.R. Donahue. Soil testing in India. U.S. Agency for International Development, Mission to India, 1965. [85] L.G. Naidu, V. Ramamuthy, G.S. Sidhu, and D. Sarkar. Emerging deficiency of potassium in soils and crops of india. Karnataka J. Agric. Sci., 24:12–19, 2011. [86] U.P. Narkhede and K.P. Adhiya. A study of clustering techniques for crop prediction - a survey. Am. Int. J. of Res. in Sci., Technol., Eng. and Math., 5(1):44–48, 2014. [87] National Research Council. Alternative agriculture. Technical report, National Academy of Sciences, Washington, DC, 1989. [88] V.P. Obade and R. Lal. Towards a standard technique for soil quality assessment. Geoderma, 265:96–102, 2016. [89] Planning Commission Government of India, editor. Eleventh Five Year Plan 2007-2012: Volume I: Inclusive Growth; Volume II: Social Sector Services; Volume III: Agriculture, Rural Development, Industry, Services, and Physical Infrastructure. Oxford University Press, New Delhi, 2008. [90] S.R. Olsen. Estimation of available phosphorus in soils by extraction with sodium bicarbonate. Circular series. U.S. Dept. of Agriculture, 1954. [91] Y. Pachepsky, K. Rajkai, and B. Tóth. Pedotransfer in soil physics: Trends and outlook - A review. AGROKÉMIA ÉS TALAJTAN, 64:339–360, 2015. [92] Y.A. Pachepsky and W.J. Rawls. Preface: Status of pedotransfer functions. Development fo Pedotransfer Functions in Soil Hydrology, 30:7–16, 2004. Bibliography 97 [93] S. Panigrahy and S.A. Sharma. Mapping of crop rotation using multidate Indian remote sensing satellite digital data. ISPRS J. of Photogrammetry and Remote Sens., 52:85–91, 1997. [94] X.E. Pantazi, D. Moshou, T. Alexandridis, R.L. Whetton, and A.M. Mouazen. Wheat yield prediction using machine learning and advanced sensing techniques. Comput. Electron. Agric., 121:57–65, 2016. [95] J. Park and I.W. Sandberg. Approximation and radial-basis-function networks. Neural Computation, 3:246–257, 1991. [96] M.A. Peña and A. Brenning. Assessing fruit-tree crop classification from Landsat-8 time series for the Maipo Valley, Chile. Remote Sens. Environ., 171:234–244, 2015. [97] A. Philibert, C. Loyce, and D. Makowski. Prediction of N2O emission from local information with random forest. Environmental Pollution, 177:156–163, 2013. [98] A. Philibert, C. Loyce, and D. Makowski. Predicting nitrous oxide emissions with a random-effects model. Environmental Modelling and Software, 61:12–18, 2014. [99] T. Pohlert. The pairwise multiple comparison of mean ranks package (PMCMR), 2014. R package. [100] R. Quinlan. C4.5: Programs for Machine Learning. Morgan Kaufmann Publishers, 1993. [101] R. Quinlan. Combining instance-based and model-based learning. In Proc. Intl. Conf. on Machine Learning, pages 236–243, 1993. [102] R.J. Quinlan. Learning with continuous classes. In 5th Australian Joint Conference on Artificial Intelligence, pages 343–348, 1992. [103] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2015. [104] P. Ramesh, N. R. Panwar, A. B. Singh, S. Ramana, S.K. Yadhav, R. Shrivastava, and A.S. Rao. Status of organic farming in india. Current Science, 98(9):1190–1194, 2010. 98 Bibliography [105] B. Rammoorthy and J.C. Bajaj. Available N,Pand Kstatus of Indian soils. Fertilizer news, 14(8):24–26, 1969. [106] M. Rashidi and M. Seilsepour. Modeling of soil total nitrogen based on soil organic carbon. ARPN J. of Agric. and Biol. Sci., 4(2):1–5, 2009. [107] D.W. Reeves. The role of soil organic matter in maintaining soil quality in continuous cropping systems. Soil and Tillage Research, 43:131–167, 1997. [108] Y. Ren, L. Zhang, and P. Suganthan. Ensemble classification and regression recent developments, applications and future directions. IEEE Computational intelligence magazine, pages 41–53, 2016. [109] L.A. Richards, L.E. Allison, L. Bernstein, C.A. Bower, J.W. Brown, M. Fireman, J.T. Hatcher, H.E. Hayward, G.A. Pearson, R.C. Reeve, A. Richards, and L.V. Wilcox. Diagnosis and improvement of saline and alkaline soils. Science, 12:800, 1954. [110] B.D. Ripley. Pattern Recognition and Neural Networks. Cambridge Univ. Press, 1996. [111] B.D. Ripley. Modern applied statistics with S. Springer, 2002. [112] J.J. Rodríguez and L.I. Kuncheva. Rotation forest: A new classifier ensemble method. IEEE Trans. on Pat. Anal. and Mach. Intel., 28(10):1619–1630, 2006. [113] J. Rogan, J. Franklin, D. Stow, J. Miller, C. Woodcock, and D. Roberts. Mapping land-cover modifications over large areas: A comparison of machine learning algorithms. Remote Sens. Environ., 112:2272–2283, 2008. [114] J.R. Romero, P.F. Roncallo, P.C. Akkiraju, I. Ponzoni, V.C. Echenique, and J.A. Carballido. Using classification algorithms for predicting durum wheat yield in the province of Buenos Aires. Comput. Electron. Agric., 96:173–179, 2013. [115] S. De Jong. SIMPLS: an alternative approach to partial least squares regression. Chemometrics and intelligent laboratory systems, 18:251–263, 1993. [116] S. De Jong. Comment on the PLS kernel algorithm. Journal of Chemometrics, 8:169–174, 1994. Bibliography 99 [117] M.G. Schaap, F. J. Leij, and M. Th. V. Genuchten. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology, 251:163–176, 2001. [118] J.L. Sehgal. Agro-ecological Regions of India. Technical bulletin (National Bureau of Soil Survey & Land Use Planning). Indian Council of Agricultural Research, 1990. [119] P.J. Sheela and K. Sivaranjani. A brief survey of classification techniques applied to soil fertility prediction. In Int. Conf. Eng. Trends in Sci. and Hum., pages 80–83, 2015. [120] N. Simon, J. Friedman, T. Hastie, and R. Tibshirani. Regularization paths for cox’s proportional hazards model via coordinate descent. Journal of statistical software, 39(5):1–13, 2011. [121] R. P. Singh and S. K. Mishra. Available macro nutrients (n, p, k, and s) in the soils of chiraigaon block of district varanasi (u.p) in relation to soil characteristics. Indian J.Sci.Res., 3(1):97–100, 2012. [122] L. Song, P. Langfelder, and S. Horvath. Random generalized linear model: a highly accurate and interpretable ensemble predictor. BMC Bioinformatics, 14(1):1–22, 2013. [123] D.F. Specht. Probabilistic neural networks. Neural Netw., 3:109–118, 1990. [124] D.F. Specht. A general regression neural network. IEEE Trans. on Neural Networks, 2:568–576, 1991. [125] B.V. Subbaiah and G.L. Asija. A rapid procedure for the estimation of available nitrogen in soil. Current Sci., 25:259–260, 1956. [126] R. Taghizadeh-Mehrjardi, K. Nabiollahi, B. Minasny, and J. Triantafilis. Comparing data mining classifiers to predict spatial distribution of USDA-family soil groups in Baneh region, Iran. Geoderma, 253–254:67–77, 2015. [127] K. Tatsumi, Y. Yamashiki, M.A. Canales Torres, and C.L. Ramos Taipe. Crop classification of upland fields using random forest of time-series Landsat 7 ETM+ data. Comput. Electron. Agric., 115:171–179, 2015. 100 Bibliography [128] The Agricultural and Processed Food Products Export Development Authority (APEDA). Indian Agro and Food Industry. Technical report, Incredible India, New Delhi, 2016. [129] M.E. Tipping. Sparse bayesian learning and the relevance vector machine. J. Mach. Learn. Res., 1:211–244, 2001. [130] M.-S. Turmel, A. Speratti, F. Baudron, N. Verhulst, and B. Govaerts. Crop residue management and soil health: A systems analysis. Agric. Systs, 134:6–16, 2015. [131] A.J. Viera and J.M. Garrett. Understanding interobserver agreement: the kappa statistic. Family Medicine, 37(5):360–363, 2005. [132] F. Wilcoxon. Individual comparisons by ranking methods. Biometrics Bulletin, 1(6):80–83, 1945. [133] S.N. Wood. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society, 1(73):3–36, 2011. [134] N. Xiao, D.S. Cao, M.Z. Li, and Q.S. Xu. Enpls: an R package for ensemble partial least squares regression. arXiv preprint, 2016. [135] T. Zhang. Adaptive forward-backward greedy algorithm for learning sparse representations. IEEE Trans. Inf. Theor., 57(7):4689–4708, 2011. [136] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society, 67:301–320, 2005. [137] H. Zou and R. Li. One-step sparse estimates in nonconcave penalized likelihood models. The Annals of Statistics, 36(4):1509–1533, 2008. List of Figures Fig. 1.1 Geographical representation of 6 regions of Maharashtra (India). Marathwada, Paschim Maharashtra and North Maharashtra are study areas, highlighted by red borders. . . . . . . . . . . . . . . . . . . . . . . . 2 Fig. 2.1 Degree of acidity and alkalinity of soil [16], with nine classes (up) and with the four classes considered in this paper (down). . . . . . . . . . . . . . . . 13 Fig. 2.2 Soil map of Maharashtra [118]. The geographical study areas are highlighted with outlines: red (Marathwada), blue (North Maharashtra) and pink (Paschim Maharashtra). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Fig. 2.3 Maps with the geographical location of crop (left panel) and soil type (right panel) data over several districts in the Marathwada region (red outline in the map of Figure 2.2). Each point locates a different village, from where several patterns are recorded. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Fig. 2.4 Intervals of κ (in %) of the different classifiers for each problem (the filled square shows the mean κ ). .......................... 25 Fig. 2.5 Intervals of percentages of the maximum κ for each problem, for all the problems and for each classifier. . . . . . . . . . . . . . . . . . . . . . . . 25 Fig. 3.1 Friedman rank of RMSE (upper panel) and correlation (lower panel) for the 20 best regressors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Fig. 3.2 Average RMSE (upper panel) and correlation (lower panel) over all the datasets for the 20 best regressors. . . . . . . . . . . . . . . . . . . . . . . 54 Fig. 3.3 Best correlation achieved by some regressor for each dataset. . . . . . . . . 56