David Camilo Corrales
email: dcorrales@unicauca.edu.co

Céline Schoving

Hélène Raynal

Philippe Debaeke

Etienne- Pascal Journet

Julie Constantin

A surrogate model based on feature selection techniques and regression learners to improve soybean yield prediction in southern France

Keywords: STICS, regression learners, filter, wrapper, embedded 1-Introduction

Empirical and process-based models are currently used to predict crop yield at field and regional levels. A mechanistic model named STICS (Multidisciplinary Simulator for Standard Crops) has been used to simulate soybean grain yield in several environments, including southern France. STICS simulates at a daily step the effects of climate, soil and management practices on plant growth, development and production. In spite of good performances to predict total aboveground biomass, poor results were obtained for final grain yield. In order to improve yield prediction, a surrogate model was developed from STICS dynamic simulations, feature selection techniques and regression learners. STICS was used to simulate functional variables at given growth stages and over selected phenological phases. The most representative variables were selected through feature selection techniques (filter, wrapper and embedded), and a subset of variables were used to train the regression learners Linear regression (LR), Support vector regression (SVR), Back propagation neural network (BPNN), Random forest (RF), Least Absolute Shrinkage and Selection Operator (LASSO) and M5 decision tree. The subset of variables selected by wrapper method combined with regression models SVR (R 2 = 0. 7102; subset of variables = 6) and LR (R 2 = 0. 6912; subset of variables = 14) provided the best results. SVR and LR models improved significantly the soybean yield predictions in southern France in comparison to STICS simulations (R 2 = 0.040).

regions being South-West and Center-East parts with an increasing contribution of organicallygrown production. One main objective in France is to achieve self-sufficiency at least in non-GMO soybean meals at 2025 horizon. Therefore, EU-27 members share a common objective:

"reducing markedly the dependency upon soybean imports by developing European production". Soybean crop requires few pesticides, no N-fertilizer and less irrigation than maize, results in low emissions of greenhouse gases, hence bringing environmental benefits. In addition, it could contribute as a summer crop to the diversification of winter cereal-based systems.

Grain yield in France slightly increased since the 80s (Terres [START_REF] Univia | Plantes riches en protéines[END_REF]. In 2019, yields were 2.61 t.ha -1 for France, 2.09 t.ha -1 for all Europe, but 3.19, 3.18, and 3.33 t.ha -1 for USA, Brazil and Argentina respectively (FAOSTAT, 2021). Climate change and its impacts on temperature, precipitation, and CO2 concentration, but also on water resources available for irrigation, will certainly impact the future production [START_REF] Porter | Making Sense of Methods and Measurements: Simple Linear Regression[END_REF][START_REF] Guilpart | Data-driven yield projections suggest large opportunities to improve Europe's soybean self-sufficiency under climate change[END_REF][START_REF] Kothari | First Soybean Multi-model Sensitivity Analysis to CO 2 , Temperature, Water, and Nitrogen[END_REF]. In addition, expanding soybean growing areas northward and introducing new cropping systems (e.g. double cropping with cereals, rainfed or irrigated soybean, etc.) will change the potential and attainable grain yields.

Therefore predicting soybean yield in various environments and a range of cropping systems will be necessary to evaluate the ability of France and European countries to achieve their objectives in terms of protein self-sufficiency by growing more soybean in cropland. Modeling can be efficient in yield analysis and investigation of the limiting factors due to easy manufacturing, testing, applying, understanding and interpretation of results [START_REF] Nehbandani | Estimating crop model parameters for simulating soybean production in Iran conditions[END_REF].

Yield prediction models are based on historical or future climate data for evaluating production potentials; also, yield prediction models assimilate remote sensing information when applied to in-season prediction. Nowadays, both statistical and mechanistic approaches are used in agricultural modelling, especially for yield prediction. Statistical approaches search and explore the relations between data to explain the variables of interest whereas mechanistic models are based on the description of biophysical processes. Dynamic crop models simulate daily growth and development in relation with environmental resources and agricultural inputs; they allow the testing of functional hypotheses and the identification of potential constraints to crop growth and yield [START_REF] Purcell | Simulating Soybean Yield Potential under Optimum Management[END_REF]. However, mechanistic and statistical approaches can be combined in order to improve the crop modeling predictions [START_REF] Casadebaig | SUNFLO, a model to simulate genotype-specific performance of the sunflower crop in contrasting environments[END_REF][START_REF] Casadebaig | A new approach to crop model calibration: Phenotyping plus post-processing[END_REF].

Statistical models from traditional Artificial Neural Networks (ANN) and Deep Learning (DL) have been used for soybean yield prediction. ANN models were proposed by [START_REF] Kaul | Artificial neural networks for corn and soybean yield prediction[END_REF] in order to predict Maryland soybean yield at state, regional, and local levels. ANN were developed using historical yield data . Field-specific rainfall data and Soil Rating for Plant Growth (SRPG) values were used for each location. The work developed in [START_REF] Maimaitijiang | Soybean yield prediction from UAV using multimodal data fusion and deep learning[END_REF] estimated the soybean grain yield through multispectral images (information type: canopy spectral, structure, thermal and texture features) and DLN in Columbia, Missouri. A Convolutional Neural Network (CNN) for soybean yield prediction in 15 states of CONUS (United States) is proposed in [START_REF] Sun | County-Level Soybean Yield Prediction Using Deep CNN-LSTM Model[END_REF]. The model was trained by crop growth and environment variables, which include weather data, MODIS Land Surface Temperature data, and MODIS Surface Reflectance data. In Cachoeira do Sul, Brazil, a Multi-Layer Perceptron (MLP) was used to adjust a predictive model for estimating the yield of soybean crop based on 9 vegetation indices [START_REF] Eugenio | Estimation of soybean yield from machine learning techniques and multispectral RPAS imagery[END_REF]. A soybean yield model was created by deep learning framework using CNN and recurrent neural networks [START_REF] Khaki | A CNN-RNN Framework for Crop Yield Prediction[END_REF]. Model was built based on environmental data and management practices from Corn Belt (including 13 states) in the United States. In southern Brazil different type of indices as Normalized Difference Vegetation Index (NVDI), Enhanced vegetation index (EVI), land surface temperature (LST) and precipitation were used to build a model using Long-Short Term Memory (LSTM), Neural Networks [START_REF] Schwalbert | Satellite-based soybean yield forecast: Integrating machine learning and weather data for improving crop yield prediction in southern Brazil[END_REF]. ANN was developed to evaluate the relative importance of predictor variables as vegetation indices (NDVI, red edge NDVI and simple ratio-SR) and elevation derived variables (slope, flow accumulation, aspect) for the prediction of soybean in Ontario, Canada [START_REF] Kross | Using Artificial Neural Networks and Remotely Sensed Data to Evaluate the Relative Importance of Variables for Prediction of Within-Field Corn and Soybean Yields[END_REF].

Traditional regression models/analysis were used with the same purpose. Authors in [START_REF] Stepanov | Predicting Soybean Yield at the Regional Scale Using Remote Sensing and Climatic Data[END_REF] used backward stepwise in order to build a regression model in Khabarovsk, Russia. Several data sources were used as Moderate Resolution Imaging Spectroradiometer (MODIS), arable land mask and meteorological stations. The NDVI was used to improve soybean yield predictions using the flexible Fourier transform model in United States [START_REF] Xu | Predicting Soybean Yield with NDVI Using a Flexible Fourier Transform Model[END_REF]. The objective was to explore the relationships between soybean yield and number of grains (NG) and thousand grains weight (TGW), generate equations to estimate yield in several countries in the years 2010 and 2019 [START_REF] Wei | Soybean Yield Estimation and Its Components: A Linear Regression Approach[END_REF]. A regression analysis was developed by [START_REF] Ma | Early Prediction of Soybean Yield from Canopy Reflectance Measurements[END_REF] in order to study the correlations between plant canopy reflectance and aboveground biomass for early prediction of soybean yield in Canada.

In addition, mechanistic models were calibrated in order to predict soybean yield. Authors in [START_REF] Robertson | Simulating growth and development of soybean in APSIM[END_REF] used Agricultural Production Systems Simulator (APSIM) with aim to simulate the soybean yield. The model was tested on an independent set of experiments, from northern Australia, with factors such as cultivars, sowing date, irrigation, soil type, plant population density row spacing varying. The research proposed by [START_REF] Jagtap | Adaptation and evaluation of the CROPGRO-soybean model to predict regional yield and production[END_REF] developed a procedure to simulate soybean yield and production by linking the CROPGROsoybean model with a regional resolution (about a 50 km grid cell) database of weather, soils, management, and varieties in the state of Georgia over the 1974-1995 time period. The CROPGRO-Soybean model was calibrated to estimate potential yields and yield gaps of soybean for 21 locations regions in India [START_REF] Bhatia | Analysis of potential yields and yield gaps of rainfed soybean in India using CROPGRO-Soybean model[END_REF]. Authors in [START_REF] Ovando | Evaluating accuracy of DSSAT model for soybean yield estimation using satellite weather data[END_REF] simulated soybean yield using a DSSAT model through weather data from Clouds and Earth's Radiant Energy System (CERES) and Tropic Rainfall Measurement Mission (TRMM) for 2006 -2016 in Oliveros, Argentina. In [START_REF] Battisti | Gauging the sources of uncertainty in soybean yield simulations using the MONICA model[END_REF], the model for Nitrogen and Carbon in Agroecosystems (MONICA) was used to simulate soybean grain yield for 14 sites in Southern Brazil. The models CSM-CROPGRO-Soybean and STICS (Multidisciplinary Simulator for Standard Crops) were used to simulate soybean yield responses under near (2041-2070) and distant (2071-2100) future climate scenarios in eastern Canada [START_REF] Jing | Modelling soybean yield responses to seeding date under projected climate change scenarios[END_REF].

The AgMip initiative started an inter-comparison of 10 soybean crop models at 5 locations in major global production areas with high quality observed data for calibration [START_REF] Kothari | First Soybean Multi-model Sensitivity Analysis to CO 2 , Temperature, Water, and Nitrogen[END_REF]. Among the tested models, STICS [START_REF] Brisson | Conceptual basis, formalisations and parameterization of the STICS crop model[END_REF], a widely used soil-plant crop model applied on a wide range of crops (305 papers in Web of Science), appeared as moderately performing with a prediction gap. Previous attempts to validate STICS on soybean concluded to good performance in aboveground biomass prediction but poor results on grain yield and protein concentration [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF]. In eastern Canada, [START_REF] Jégo | Calibration and performance evaluation of soybean and spring wheat cultivars using the STICS crop model in Eastern Canada[END_REF] obtained scattered results for biomass, LAI and yield, with Root Mean Square Error (RMSE) from 23 to 38 %. Heretofore, this model has been used less in soybean crops, but encouraging results were recently obtained through proper calibration on pea and faba bean, two other grain legumes [START_REF] Falconnier | Calibration and evaluation of the STICS soil-crop model for faba bean to explain variability in yield and N2 fixation[END_REF][START_REF] Falconnier | Contrasted response to climate change of winter and spring grain legumes in southwestern France[END_REF]. In addition, STICS describes with a good accuracy the dynamics of water and nitrogen in soils and a module was introduced for considering explicitly N2fixation in legumes.

Improving the STICS prediction for soybean could imply a thorough calibration of the model and/or a deep re-examination of the underlying biophysical processes. The experimental data required for such an improvement could be insufficient. The data-driven modeling approach has emerged as an alternative to model the biophysical system purely from the data available.

A data-driven model, also known as a meta-model or surrogate model, is a "model of the model" [START_REF] Cui | A recommendation system for meta-modeling: A metalearning based approach[END_REF]. A surrogate model is a statistical model trained from simulations/variables difficult to measure in field (e.g. leaf area index, aboveground crop biomass, N crop uptake, crop transpiration, etc.). The surrogate model can be deployed to replace or support the original biophysical simulation module to accurately approximate the simulation output. Figure 1 explains the interactions between surrogate and STICS model.

2-Materials and methods

Multidisciplinary Simulator for Standard Crops (STICS)

The STICS model simulates at a daily step the effects of climate, soil and management practices on plant growth, development and production (quantity and quality) and environmental impacts. The combination of these input variables is termed a USM (Unit of SiMulation). Each USM corresponds to one execution of the STICS model [START_REF] Brisson | STICS: a generic model for the simulation of crops and their water and nitrogen balances. I. Theory and parameterization applied to wheat and corn[END_REF]. STICS can be tuned to a single crop, two intercropped or several successive crop cycles. STICS has been evaluated over a large data set for 15 different crops and different conditions of soil and climate in France [START_REF] Coucheney | Accuracy, robustness and behavior of the STICS soil-crop model for plant, water and nitrogen outputs: Evaluation over a wide range of agro-environmental conditions in France[END_REF].

In order to calibrate STICS, the crop files contain species parameters, ecophysiological options (e.g. effect of photoperiod and/or cold requirements on crop phenology, potential radiation use efficiency) and cultivar specific parameters (e.g. flowering precocity, maximum number of grains per m 2 ). Crop temperature (calculated from weather variables) and photoperiod drive crop phenology. The model dynamically simulates (i) the development of the root system that takes up N and water according to root density over the whole soil profile and (ii) the establishment of the canopy that transpires water and intercepts light to produce the crop biomass [START_REF] Brisson | Conceptual basis, formalisations and parameterization of the STICS crop model[END_REF].

Study area and datasets

The data used in this work were collected by [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF] The eleven tested soybean varieties belonged to four maturity groups corresponding to different crop durations and potential yields: 000 (very early-maturing), 0, I and II (late-maturing). Three late-maturing varieties were tested in all experiments since 2010 (Ecudor, Santana, Isidor and Sarema). Detailed information on varieties and maturity groups are presented in Table A1 (Appendix A). Weather data were collected near to the experimental sites. Soil samples contain texture and physico-chemical analyzes. These data are essential to correctly initialize STICS with realistic values of soil moisture and mineral nitrogen (nitrate, ammonium). The water pressure of the soils was monitored in a micro-plot of Santana variety, at 30, 60 and 90 cm depth.

The dataset contains 227 simulation units (USM) created from combination of experimental sites, years and cropping practices (cultivar, water management and sowing date); We used the same training (105 USMs) and test (122 USMs) datasets as defined by [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF] For instance, preprocessed variables computed by phenology stages as Lai(n)_Iflo variable indicates the leaf area index Lai(n) at flowering (Iflo). Other variables are expressed over phenophases, for example Precip_cum_Iflo-Imat represents the cumulative precipitation (Precip) between two phenological stages: flowering (Iflo) and physiological maturity (Imat).

In addition, we define variables by phenophases and thresholds as Days_MinTemp_18_Ilev-Imat variable which indicates the number of days when minimum temperature was less than 18°C between emergence and physiological maturity stages. Table 2 presents a summary of preprocessed variables from STICS by stages or phases, descriptive values and thresholds.

Table A2 (Appendix A) lists all of the preprocessed variables.

STICS variables

Stages or phases

Descriptive value In order to create an interpretable soybean yield model with a minimum number of variables, we applied feature selection techniques to previous dataset based on three approaches: filter, embedded and wrapper.

Threshold Number of variables Nbgrmax - - - 1 Stlevdrp - - - 1 Stflodrp - - - 1 Stdrpmat - - - 1 Masec(n) Iflo, Idrp, Imat - - 3 Lai(n) Iflo, Idrp, Imat - - 3 Qnplante Iflo, Idrp, Imat - - 3 Qfix Iflo, Idrp, Imat - - 3 Zrac Iflo, Idrp, Imat - - 3 Jul Ilev-Imat, Iflo-Imat, Idrp-Imat, Iflo-Idrp Cum - 4 Raint Ilev-Imat, Iflo-Imat, Idrp-Imat, Iflo-Idrp Cum - 4 Etpp(n) Ilev-Imat, Iflo-Imat, Idrp-Imat, Iflo-Idrp Cum - 4 Precip Ilev-Imat, Iflo-Imat, Idrp-Imat, Iflo-Idrp Cum - 4 Ep Ilev-Imat, Iflo-Imat, Idrp-Imat,

Feature selection techniques

Feature selection is the process (automatic or manual) of selecting a subset of relevant variables which contribute most to learner [START_REF] Corrales | Feature selection for classification tasks: Expert knowledge or traditional methods?[END_REF]. Feature selection techniques can be grouped in three categories:

 Filter methods are based only on the intrinsic properties of the data [START_REF] Solorio-Fernández | A Supervised Filter Feature Selection method for mixed data based on Spectral Feature Selection and Information-theory redundancy analysis[END_REF]. Filter method computes an importance value between one independent variable and the dependent variable. Variables with highest importance values are selected based on user criteria. Filter methods are usually computationally less expensive than embedded and wrapper methods. We used classical feature selection methods based on Pearson [START_REF] Pearson | Notes on the History of Correlation[END_REF] and Spearman [START_REF] Spearman | The Proof and Measurement of Association Between Two Things[END_REF] coefficients and Information Gain [START_REF] Shannon | A mathematical theory of communication[END_REF]. In order to explain the statistical coefficients, independent variable is named X and dependent variable Y.

o Pearson coefficient measures the linear correlation between two variables. If both variables are linearly dependent, then their correlation coefficient is close to ±1. If the variables are uncorrelated, the correlation coefficient is 0 [START_REF] Pearson | Notes on the History of Correlation[END_REF]. When Pearson coefficient is used as filter method only positive values are considered following the equation ( 1):

𝑟 = | ∑ 𝑐𝑜𝑣(𝑥 𝑖 , 𝑦 𝑖 ) 𝑛 𝑖=1 ∑ 𝜎𝑥 𝑖 𝜎𝑦 𝑖 𝑛 𝑖=1 | (1)
Where 𝑥 𝑖 and 𝑦 𝑖 are the i th observations of independent and dependent variable respectively; 𝑐𝑜𝑣 corresponds to covariance and 𝜎 indicates the standard deviation of 𝑥 and 𝑦.

o Spearman coefficient through a monotonic function measures the correlation between two variables [START_REF] Spearman | The Proof and Measurement of Association Between Two Things[END_REF]. A monotonic function is defined as function which is either entirely increasing or decreasing. It is similar to Pearson coefficient except that it operates on the ranks of the data rather than the raw data [START_REF] Gauthier | Detecting trends using Spearman's rank correlation coefficient[END_REF]. The Spearman correlation rank is defined by equation ( 2):

𝜌 = 1 - 6 ∑ 𝑑 𝑖 2 𝑛 𝑖=1 𝑛(𝑛 2 -1) (2)
Where 𝑑 𝑖 is the difference between ranks for each 𝑥 𝑖 , 𝑦 𝑖 data pair and 𝑛 is the number of data pairs. o Entropy-based information gain discretizes the independent variable and subsequently the entropy is computed between 𝑥 and continuous 𝑦 variable [START_REF] Yang | A multi-filter enhanced genetic ensemble system for gene selection and sample classification of microarray data[END_REF] by equation ( 3):

𝐼𝑛𝑓𝑜𝐺𝑎𝑖𝑛 = 𝐻(𝑦) + 𝐻(𝑥) -𝐻(𝑦, 𝑥) (3)
Where 𝐻(𝑥) and 𝐻(𝑦) correspond to Shannon's Entropy for 𝑥 and 𝑦 variables.

𝐻(𝑦, 𝑥) is a joint Shannon's Entropy for a variable 𝑦 with a condition to 𝑥.

Detailed explanation of Shannon's Entropy is explained in [START_REF] Shannon | A mathematical theory of communication[END_REF].

 Embedded methods integrate the variables selection as part of training process into learner [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. We used as embedded methods the learners Random Forest [START_REF] Breiman | Random forests[END_REF], M5 decision tree [START_REF] Quinlan | Learning with continuous classes[END_REF] and Least Absolute Shrinkage and Selection Operator [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

 Wrapper methods selects a subset of variables according to performance criteria (regression tasks, measure of errors as mean absolute error, mean square error, root square mean error; classification tasks measures as accuracy, precision, overall, recall, etc.) of a learner [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. In this paper, we use the method Recursive Feature Elimination [START_REF] Guyon | Gene Selection for Cancer Classification using Support Vector Machines[END_REF].

Regression learners

In order to predict soybean yield in southern France crops, we used regression learners from different families of algorithms. They are explained briefly below.

 Linear regression (LR) explains the relationship between dependent variable and one or more independent variables by fitting a linear equation to observed data [START_REF] Porter | Making Sense of Methods and Measurements: Simple Linear Regression[END_REF]. Coefficients multiply the values of dependent values; the coefficient signs represent the direction of the relationship between a dependent variable and the independent variable.

 Support vector regression (SVR) is based on same principles as Support Vector

Machine [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. SVR determines a regression function in the feature space considering only data points within the decision boundary lines called support vectors.

In nonlinear data, a kernel function is used in order to transform the feature space into a linear hyperplane [START_REF] Brereton | Support Vector Machines for classification and regression[END_REF].

 Back propagation neural network (BPNN) calculates the gradient of the error function with respect to the weights of the neural network [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF]. The computed error is propagated in a backward manner from one layer to the other until the minimum Mean Squared Error (MSE) is attained and weights can be modified accordingly [START_REF] Deshwal | A Language Identification System using Hybrid Features and Back-Propagation Neural Network[END_REF].

In addition, we used each of Random forest, Least Absolute Shrinkage and Selection Operator and M5 decision tree as both an embedded method and a wrapper.

 Random forest (RF) builds several decision trees using a different bootstrap sample of data training set [START_REF] Breiman | Random forests[END_REF]. The decision trees are built using CART learner [START_REF] Breiman | Classification and regression trees[END_REF]. In regression tasks, RF final prediction is obtained by averaging the results of all the CART trees.

 Least Absolute Shrinkage and Selection Operator (LASSO) is a linear regression method which imposes a bound on the L1-norm of the regression coefficients, resulting in coefficient shrinkage [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. LASSO adds a L1 penalty equal to the absolute value of the magnitude of coefficients (Equation 1). Variables are discarded when the coefficients take values equal to zero. Larger penalties are expressed by coefficient values closer to zero. The objective function for finding the minimum is shown by Equation (4):

𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒𝛽 0 , 𝛽( 1 2𝑁 ∑ (𝑦 𝑖 -𝛽 0 -𝑥 𝑖 𝑇 𝛽) 2 𝑁 𝑖=1 + 𝜆 ∑ |𝛽 𝑗 | 𝑝 𝑗=1 ) ( 4
)
Where 𝑁 is the total number of observations, 𝜆 is a nonnegative regularization parameter corresponding to one value of Lambda, 𝑦 𝑖 is the dependent variable, 𝑝 is the number of independent variables 𝑥 𝑖 = (𝑥 𝑖1 , … , 𝑥 𝑖𝑝 ) 𝑇 , 𝛽 0 is the intercept, and 𝛽 𝑗 are the coefficients [START_REF] Shi | Prioritizing driving factors of household carbon emissions: An application of the LASSO model with survey data[END_REF].

 M5 is a conventional decision tree composed by different nodes such as root, intermediates and leaves [START_REF] Quinlan | Learning with continuous classes[END_REF]. Root and intermediate nodes are chosen over the dependent variable that maximizes the expected error reduction as a function of the standard deviation of output parameter [START_REF] Wang | Induction of model trees for predicting continuous classes[END_REF]. Leaves nodes predict the dependent variable through linear regression functions; they are fitted from data that follows one branch between root and leaf node.

3-Results

This section reports the subset of variables selected by feature selection methods and the evaluation of soybean yield models built from selected variables.

Feature selection

Filter methods

Filter methods were computed through R 'mlr' package [START_REF] Bischl | mlr: Machine Learning in R[END_REF]. With aim to create simple models with great explanatory predictive power with a minimum number of features, we defined two criteria to select the variables: (i) top-15 of variables with highest importance values for each filter method; (ii) variables selection by threshold based on importance values of top-15. Regarding to Pearson and Spearman coefficients, we selected the features with importance value greater than or equal to 0.6. Concerning to entropy-based information gain, we selected the features with importance value greater than or equal to 0.4. We defined these thresholds following the "conventional interpretation of the correlation coefficients" proposed in [START_REF] Schober | Correlation Coefficients: Appropriate Use and Interpretation[END_REF]. Values between 0.60 -0.79 are defined as "moderately correlated" and coefficient values between 0.70 -0.89 are interpreted as "strongly correlated". Table 3 shows Top-15 of soybean variables selected by filter methods (Appendix B contains the entire ranking). 5.

Pos

Embedded methods

We used as embedded methods the learners: Random forest, Least Absolute Shrinkage and Selection Operator and M5 decision tree. The variable selection process for each one is explained below.

Random Forest

Random Forest (RF) gathers a set of CART trees in order to obtain the soybean yield prediction by averaging the results of all of trees. We used R packages 'randomForest' to create RF model and 'randomForestExplainer' [START_REF] Liaw | Classification and regression by randomForest[END_REF] to design multi-way plot shown in Figure 3. Five hundred CART trees were built (ntree parameter) as result, 42 variables were sampled as candidates at each split. The multi-way plot focuses on three importance measures that derive from the structure of trees in the forest: (i) the mean_min_depth variable refers to the depth of first split on the variable to the top of the tree; (ii) the times_a_root variable measures the number of times a variable is set as top of a decision tree.

Least Absolute Shrinkage and Selection Operator

We used the R package 'caret' to build LASSO model [START_REF] Kuhn | Building predictive models in R using the caret package[END_REF]. Root Mean Square Error (RMSE) was used to select the optimal model using the smallest value. LASSO model was run with parameters fraction = 0.1 and lambda = 0.01 (Equation 1). The LASSO model set the regression coefficients of 21 variables to zero by imposing the L1 penalty. Table C1 (Appendix C) contains the regression coefficients of 66 variables calculated by LASSO.

M5 decision tree

M5 was built by R package 'RWeka' [START_REF] Hornik | Open-source machine learning: R meets Weka[END_REF]. The construction of M5 tree is based on recursive splitting of the standard deviation of Mafruit (dependent variable) that reach a node as a measure of the error at the node. The variable that maximizes the expected error reduction is selected for splitting at the node. The expected reduction of the error is obtained as a result of testing each variable at that node. To remove the problem of over fitting, M5 uses a method to prune back the over grown tree. Figure 4 shows the structure of the pruned M5 tree to predict the soybean yield (Mafruit) regarding the thirteen selected variables (Table 4). Left branches show conditions below the node threshold, for example, if the top node Lai(n)_Imat ≤ 1.103 and if Inn_avg_Iflo-Idrp ≤ 0.769, the linear model LM 1 must be selected. Thirty-seven nodes were created including nineteen linear models as decision rules. We defined 3 observations as minimum number at the leaf node. Appendix C contains the linear models generated by M5.

Wrapper methods

A greedy algorithm named Recursive Feature Elimination (RFE) was used as wrapper method [START_REF] Kohavi | Wrappers for feature subset selection[END_REF][START_REF] Guyon | Gene Selection for Cancer Classification using Support Vector Machines[END_REF]. RFE selects a subset of variables to improve a learner performance by removing the least important features. This process is repeated recursively (based on backward elimination) until the optimal number of features is obtained.

The recursive step plays a key role due the relative importance of each variable can change substantially over a different subset of variables during the backward elimination process (particularly for highly correlated variables) [START_REF] Granitto | Recursive feature elimination with random forest for PTR-MS analysis of agroindustrial products[END_REF]. We used the RFE version coded in R package 'caret' [START_REF] Kuhn | Building predictive models in R using the caret package[END_REF]. Six learners were selected (explained in section 2) to determine the performance of subsets of variables selected by RFE shown in RF andM5 (1421.25, 1268.27 and 1033.52 seconds). In contrast to Support Vector Regression which imposed considerable computational cost (5271.56 seconds) due to margin maximization to find the support vectors and nonlinear transformations of the feature space [START_REF] Yu | Classifying large data sets using SVMs with hierarchical clusters[END_REF].

Regression models

In order to examine the performance of subset of variables selected by feature selection methods, we used traditional statistical criteria to estimate the prediction accuracy of regression learners as Coefficient of determination (R 2 ), Mean Square Error (MSE), Root Mean Square Error (RMSE) and Mean Absolute Error (MAE). The metric R 2 was adopted to measure the correlation of the observed and predicted values.

Regression learners were trained with 105 observations and evaluated with 122 examples (explained in section 2). According to filter methods, Pearson, Spearman and Entropy coefficients selected 12, 14 and 5

variables respectively (Section 3). The regression learners LR, SVR and BPNN were trained with subset of variables selected by filter methods; SVR obtained the best results with subset of feature selected by Pearson coefficient (R 2 = 0.6389, MSE = 0.6539, RMSE = 0.8086 and MAE = 0.6267). Regression learners RF, LASSO and M5 decision tree act as ensemble approaches for feature selection due these learners are considered embedded methods as well.

In other words, embedded methods were trained with subset of variables selected by filter methods. Subsequently, embedded methods select a new subset of variables in the training process into learner. However, the ensemble of filter and embedded methods do not improve the performance obtained by SVR trained with subset of variables chosen by Pearson coefficient.

The number of variables selected by two embedded methods were higher respecting to others feature selection approaches. Random Forest and LASSO selected 87 and 66 variables respectively, whereas M5 decision Tree 13 variables (Section 3). Random Forest reached the best performance in the ensemble methods (R 2 = 0. 5020, MSE = 0. 8300, RMSE = 0. 9110 and MAE = 0.7258), however RF does not overcome the performance obtained by SVR and Pearson coefficient in the filter methods.

Concerning to wrapper method, Recursive Feature Elimination selected 6, 14 and 16 variables For real values of soybean yield equal to zero, RFE-SVR and RFE-LR predict negative soybean values or close to zero. RFE-SVR simulates Mafruit values equal to 0.066, -0.086 and -0.085 for 3 observations where Mafruit = 0 of the validation dataset (Béziers in 2011 with maturity group I in varieties Isidor, Santana and maturity group II in Ecudor variety). Similarly, LR-RFE predicts -0.227, -0.379 and -0.378 for same observations of validation dataset.

Comparative study

In order to demonstrate the performance of RFE-SVR and RFE-LR, we compared both regression models against STICS simulations developed in [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF]. The soybean yield model proposed by Schoving was calibrated with the same training set as presented in section 2. Table 6 shows the results of RFE-SVR, RFE-LR and STICS models evaluated from validation dataset (Section 2). [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF]. Validation dataset (Section 2) was used to evaluate the models. Statistical criteria used to estimate the performance of models: Coefficient of determination (R 2 )

Models

and Root Mean Square Error (RMSE).

The two regression models explained around 70 % of the grain yield variation and they achieved half of RMSE values obtained for STICS simulations developed in [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF]. Although total aboveground biomass was correctly simulated by [START_REF] Schoving | Analyse écophysiologique et modélisation dynamique des intéractions génotype x environnement x conduite de culture chez le soja[END_REF] (r² = 0.64), final grain yield of semi-indeterminate and indeterminate soybean cultivars was poorly represented (Table 6; r² = 0.04). This is probably because STICS uses the standard formalism of wheat and maize crops to simulate the final grain yield in soybean.

4-Discussion

Feature selection methods are relevant in order to reduce the computational complexity and improve the model generalization ability [START_REF] Maldonado | A wrapper method for feature selection using Support Vector Machines[END_REF]. High dimensional variables impose a high computational cost and a high cost of data acquisition. On the other hand, a low-dimensional representation reduces the risk of overfitting [START_REF] Famili | Data preprocessing and intelligent data analysis[END_REF][START_REF] Liu | FS_SFS: A novel feature selection method for support vector machines[END_REF]. The value added by applying feature selection techniques is to determinate a subset of available variables to build a good model, which is a combinatorial problem in the number of original variables [START_REF] Wolsey | Integer and combinatorial optimization[END_REF][START_REF] Guyon | An introduction to variable and feature selection[END_REF].

With this rationale, we used popular filter, embedded and wrapper methods in order to select the relevant variables to predict soybean yield.

Filter methods are independent of any learners and they are based on performance evaluation metrics calculated directly from the data. We used two correlation-based filters (Pearson and Spearman) followed by a measure based on information theory (Information Gain). Figure 7 presents Venus diagram of top-15 variables selected by filter methods. Correlation-based filters selected twelve same variables. The correlations found by Pearson and Spearman are equivalent to rank the common twelve variables; whereas Pearson inspects straight connections, Spearman evaluates monotonic connections (regardless of whether direct or not) [START_REF] Thirumalai | Analysing the concrete compressive strength using Pearson and Spearman[END_REF].

However, of the three filter methods share less features in common (eight variables) due to entropy-based information is designed to observe the amount of information gained between two discrete variables, and dataset used contains only numeric variables.

The principle of embedded methods (feature selection as part of the training process) is to reduce the computation time taken up for testing several subsets of variables which is done in wrapper methods [START_REF] Chandrashekar | A survey on feature selection methods[END_REF]. In this work, the number of variables selected by RF and LASSO were higher. LASSO selected 66 variables and Random Forest all of the initial variable set. Although M5 decision Tree selected less variables (13), the results

were worse than RF (R 2 = 0. 4010, MSE = 0. 9466, RMSE = 0. 9729 and MAE = 0. 7802).

Embedded methods showed an improvement when they were used as ensemble approaches of filter and wrapper techniques. (i) they evaluate variables iteratively with respect to performance learner. Therefore, variables selected by wrapper approach are more likely to suit the learner [START_REF] Kohavi | Wrappers for feature subset selection[END_REF](ii) wrapper approaches evaluate variables jointly and are effective in capturing intrinsic relationships such as interactions among multiple variables [START_REF] Freitas | Understanding the crucial role of attribute interaction in data mining[END_REF]. However computational cost is high, even for learners that exhibit a moderate complexity, the number of iterations required by the search of subset variables is high, especially as more complex search strategies are used [START_REF] Talavera | An evaluation of filter and wrapper methods for feature selection in categorical clustering[END_REF] [START_REF] Loyola-González | Black-Box vs. White-Box: Understanding Their Advantages and Weaknesses From a Practical Point of View[END_REF]. SVR builds a hyperplane or set of hyperplanes in a high-dimensional space, which are very hard to explain and to be understood by experts in practical applications [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF]. LR generates a linear equation to explain the correlation among variables in a language close to a human expert.

In this sense, RFE-LR is represented by the following linear equation:

Figure 1 .

 1 Figure 1. Steps of STICS and surrogate model to predict the soybean yield. The red arrows show the process of soybean yield simulations generated by STICS. (1) Parameters are calibrated and field observations are used to run STICS. (2) Soybean yield simulations 𝑌 1 , 𝑌 2 … , 𝑌 𝑛 generated by STICS. The blue arrows depict the steps to create the surrogate model. (3) The STICS simulations are used to select the relevant variables and train the surrogate model. (4) Soybean yield simulations 𝑌′ 1 , 𝑌′ 2 … , 𝑌′ 𝑛 generated by the surrogate model. In this sense, we proposed a surrogate model based on feature selection techniques and regression learners to predict soybean yield in southern France. The surrogate model is trained from the data produced by STICS simulations generated by (Schoving, 2020) (effects of climate, soil and management practices on dynamic variables of soybean crop functioning) to improve the prediction of final grain yield. This study progressed through three steps: a) Calculate crop variables at different phenological stages with STICS as evaluated by (Schoving, 2020) in southern France. b) Find the representative variables of soybean yield based on feature selection techniques. c) Build a regression model to predict soybean yield based on representative variables found by feature selection techniques.

  . Seventeen experimental sites were conducted during 2010-2018 from six regions in the south of France: Mauguio (2010), Béziers (2010 -2012), Mondonville (2010 -2014), Rivières (2010 -2014), En Crambade (2013 -2014) and Auzeville (2017 -2018) as shown in Figure 2.

Figure 2 .

 2 Figure 2. Locations in Southern France where experimental sites were conducted. Locations are depicted by red markers (Mauguio, Béziers, Mondonville, Rivières, En Crambade and Auzeville). This figure was created by Google Earth.

  Figure 3 presents multi-way plot for first 15 relevant variables. Lai(n)_Imat was the most used variable as top split criterion (148 times) followed by Masec(n)_Imat (67 times) and Ep_cum_Iflo-Imat (58 times). Variables as Qfix_Iflo and Precip_cum_Idrp-Imat were never used as top of decision trees and they have the longest distance (mean minimum depth of 6.728 and 6.398 respectively) to the top of decision trees considered less associated with the dependent variable Mafruit. Other variables as Qfix_Idrp, Lai(n)_Iflo and Precip_cum_Iflo-Idrp can be considered as intermediate nodes of the trees (times_a_root = 0) with mean minimum depth less than Precip_cum_Idrp-Imat and Qfix_Iflo (mean_min_depth = 4.856, 5.045 and 5.297 respectively).

Figure 3 .

 3 Figure 3. Multi-way plot between two measures of importance: mean_min_depth (x-axis) and times_a_root (y-axis). First 15 relevant variables are depicted. X-axis correspond to mean depth of first split on the variable, y-axis the number of trees in which the root is split on the variable.

Figure 4 .

 4 Figure 4. M5 pruned decision tree. Thirteen variables were selected and 19 decision rules (linear regressions) were created. At the leaf node, the number in parenthesis corresponds to observations for each linear regression.

Figure 5

 5 presents the time complexity of the variables subset selection by Recursive Feature Elimination and base learners. Wrapper methods were run on Windows 10 comprised of Intel Core i5-774HQ CPU 2.80GHz -16GB RAM based on sequential computing.

Figure 5 .

 5 Figure 5. Time complexity to select subset of variables by Recursive Feature Elimination (RFE) and base learners: Linear Regression (LR), Support Vector Regression (SVR), M5 decision tree, Random Forest (RF) and Backpropagation Neural Network (BPNN).

  for LR, SVR and BPNN. SVR achieved the best results compared to all of combination of feature selection approaches and learners (R 2 = 0. 7102, MSE = 0. 4170, RMSE = 0. 6458 and MAE = 0. 5230). Remain of learners that work as ensemble feature selection among wrapper and embedded methods (RF, LASSO and M5 decision tree), RFE proposes different subset of variables created in the backward elimination process. The embedded methods are tested with subset of variables proposed by wrapper method and they select a new subset of variables in the training process into learner. The best subset of variables are selected regarding to performance criteria of the embedded method. LASSO trained with 16 variables selected by RFE (Table 4) reached the best results (R 2 = 0. 6718, MSE = 0. 4760, RMSE = 0. 6899 and MAE = 0. 5747) compared to RF and M5 of the wrapper method. Besides, RFE improve the performance of LASSO compared to LASSO's version of filter and embedded. In summary, the wrapper methods RFE-SVR (R 2 = 0. 7102, MSE = 0. 4170, RMSE = 0. 6458 and MAE = 0. 5230) and RFE-LR (R 2 = 0. 6912, MSE = 0. 4807, RMSE = 0. 6933 and MAE = 0. 5469) achieved the best results from validation dataset. Figure 6 depicts the scatter plots of observed vs simulated soybean yield values by RFE-SVR and RFE-LR. a. Observed Mafruit values vs Simulated Mafruit values by RFE-SVR b. Observed Mafruit values vs Simulated Mafruit values by RFE-LR Figure 6. Scatter plot of observed (x-axis) vs simulated (y-axis) Mafruit values by RFE-SVR (14 variables) and RFE-LR (6 variables)

Figure 7 .

 7 Figure 7.Venus diagram of top-15 variables selected by filter methods: Pearson, Spearman and Entropy-based information Wrapper methods are widely recognized and considered a superior alternative for two reasons:

Table 2 .

 2 Preprocessed

	Iflo-Idrp	Cum	-	4

variables from Multidisciplinary Simulator for Standard Crops (STICS). STICS variables were preprocessed by crop phenology stages (Ilev: emergence, Idrp: grain filling onset, Iflo: flowering, Imat: physiological maturity), descriptive values (Cum: cumulative, Max: maximum, Avg: average) and thresholds (MinTemp, MaxTemp, Swfac and Inn).

Last column (Var) corresponds to number of variables processed from a STICS, phenology stages, descriptive values and thresholds. Table

A2

(Appendix A) lists all of the preprocessed variables.

Pearson coefficient Spearman coefficient Entropy-based information gain Variable Importance

				Variable	Importance	Variable	Importance
	1	Masec(n)_Imat	0.7923	Lai(n)_Imat	0.7699	Lai(n)_Imat	0.5929
	2	Lai(n)_Imat	0.7726	Masec(n)_Imat	0.7311	Masec(n)_Imat	0.5875
	3	Qnplante_Idrp	0.7710	Qnplante_Idrp	0.7137	Qnplante_Imat	0.5515
	4	Qnplante_Imat	0.7523	Masec(n)_Idrp	0.7003	Ep_cum_Idrp-Imat	0.4223
	5	Masec(n)_Idrp	0.7334	Ep_cum_Iflo-Imat	0.6886	Qfix_Imat	0.4072
	6	Ep_cum_Idrp-Imat	0.6819	Ep_cum_Idrp-Imat	0.6734	Inn_min_Iflo-Idrp	0.3799
	7	Qfix_Imat	0.6805	Qnplante_Imat	0.6714	Inn_avg_Iflo-Idrp	0.3652
	8	Ep_cum_Iflo-Imat	0.6650	Ep_cum_Ilev-Imat	0.6456	Ep_cum_Iflo-Imat	0.3486
	9			Lai(n)_Idrp	0.6439	Days_Swfac_0.6_Iflo-	0.3443
		Raint_cum_Idrp-Imat	0.6326			Imat
	10 Lai(n)_Idrp	0.6292	Raint_cum_Idrp-Imat	0.6288	Inn_avg_Ilev-Imat	0.3400
	11 AvgTemp_avg_Idrp-Imat	0.6157	Swfac_min_Idrp-Imat	0.6160	Inn_avg_Iflo-Imat	0.3380
	12 Swfac_avg_Iflo-Imat	0.6088	Swfac_min_Iflo-Imat	0.6136	Masec(n)_Idrp	0.3300
	13			AvgTemp_avg_Idrp-	0.6081	Qnplante_Idrp	0.3285
		MinTemp_avg_Idrp-Imat	0.5937	Imat
	14 Ep_cum_Ilev-Imat	0.5842	Swfac_min_Ilev-Imat	0.6009	Raint_cum_Idrp-Imat	0.3080
	15 Days_Swfac_0.6_Iflo-		MinTemp_avg_Idrp-	0.5977	Days_Swfac_0.6_Iflo-	0.3025
		Imat	0.5811	Imat		Idrp

Table 3 .

 3 Top-15 of variables selected by filter methods: Pearson, Spearman and entropy-based information gain. In this sense, 12 and 14 variables were selected by Pearson and Spearman respectively. The variables Lai(n)_Imat, Masec(n)_Imat, Qnplante_Imat, Ep_cum_Idrp-Imat and Qfix_Imat were selected by entropy-based information gain. The subdatasets with the selected variables by filter methods were used to train the regression learners presented in Section 2. Besides, RF, LASSO and M5 decision tree act as regression learners and feature selection techniques due these learners are considered embedded methods. In other words, RF, LASSO and M5 are trained with subset of variables selected by filter methods and subsequently the embedded methods/learners select a new subset of variables in the training process into learner. The results are presented in Table

Table 4

 4

	. The

Table 4 .

 4 Subset of variables selected by Recursive Feature Elimination (RFE) and base learners: Linear Regression (LR), Support Vector Regression (SVR), M5 decision tree, Random Forest (RF) and Backpropagation Neural Network (BPNN).

	Learner performance is based on Root Mean Square Error (RMSE).

Concerning time complexity, RFE is slower than filter and embedded methods, since RFE needs to evaluate performance criteria for each iteration besides the computational cost of the model training. In this sense, learners based on linear models as LR and LASSO obtained much less computational cost (59.42 and 137.53 seconds) than BPNN,

Table 5 .

 5 Table 5 presents the results of regression learners trained with subset of variables selected by feature selection methods. The underlined values correspond to statistical criteria obtained by best two learners and feature selection method. Results of regression learners trained with subset of variables selected by feature selection methods. Validation dataset

	Feature	Feature				Learners
	selection	selection	SC	LR	SVR	BPNN	RF	LASSO	M5
	approach	method
			R 2	0.6185 0.6389	0.1837	0.3487	0.6274	0.6193
		Pearson	MSE	0.9085 0.6539	2.1301	1.0325	0.7250	0.9211
		correlation	RMSE	0.9531 0.8086	1.4595	1.0161	0.8514	0.9597
			MAE	0.7227 0.6267	1.2375	0.8422	0.6552	0.7258
			No. Selected	12	12	12	12	12	12
			Variables
	Filter		R 2	0.5890 0.6394	0.1841	0.3463	0.6101	0.6276
		Spearman	MSE	1.0418 0.6948	1.4388	1.0818	0.8167	0.6776
		correlation	RMSE	1.0207 0.8335	1.1995	1.0401	0.9037	0.8231
			MAE	0.7889 0.6551	1.0158	0.8505	0.6862	0.6397
			No. Selected	14	14	14	14	14	14
			Variables
			R 2	0.3974 0.3916	0.2701	0.2622	0.3961	0.2502
		Entropy-based	MSE	0.9479 1.0085	4.0698	1.1859	0.9554	1.2713
		information	RMSE	0.9736 1.0042	2.0173	1.0890	0.9774	1.1275
		gain	MAE	0.7951 0.8215	1.7373	0.8815	0.7975	0.9212
			No. Selected	5	5	5	5	5	5
			Variables
			R 2	-	-	-	0.5020	0.1249	0.4010
	Embedded	Learners:	MSE	-	-	-	0.8300	6.1841	0.9466
		RF, LASSO	RMSE	-	-	-	0.9110	2.4867	0.9729
		and M5	MAE	-	-	-	0.7258	1.9206	0.7802
			No. Selected	-	-	-	42	66	87
			Variables
			R 2	0.6912 0.7102	0.1829	0.5020	0.6718	0.4010
	Wrapper	Recursive	MSE	0.4807 0.4170	1.4916	0.8300	0.4760	0.9466
		Feature	RMSE	0.6933 0.6458	1.2213	0.9110	0.6899	0.9729
		Elimination	MAE	0.5469 0.5230	0.9746	0.7258	0.5747	0.7802
		(RFE)	No. Selected	6	14	16	87	16	13
			Variables
	(Section 2) was used to evaluate the regression learners. Statistical criteria (SC) used to estimate the performance of regression
	learners: Coefficient of determination (R2), Mean Square Error (MSE), Root Mean Square Error (RMSE) and Mean Absolute
	Error (MAE). The underlined values correspond to statistical criteria obtained by best two learners and feature selection method.

Table 6 .

 6 Comparison of Support Vector Regression and Linear Regression (LR) trained with subset of variables selected by Recursive Feature Selection (RFS) vs soybean yield model proposed by

		Statistical criteria
		R 2	RMSE
	RFE-SVR	0.710	0.645
	RFE-LR	0.691	0.693
	(Schoving, 2020)	0.040	1.320

Table A1 .

 A1 Dataset description of seventeen experimental sites during 2010-2018 from six regions in the south of France. The dataset contains 227 simulation units (USM) created from combination of experimental sites, years and cropping practices (cultivar, water management and sowing date).

	Mafruit = -0.3095*MaxTemp_avg_Idrp-Imat + 0.0148*Qnplante_Idrp + 0.0057*Masec(n)_Imat + 0.0066*Ep_cum_Idrp-Imat + 0.0458 *Lai(n)_Imat + 0.0008*Raint_cum_Idrp-Imat + 8.5049 APPENDIX A	(9)
	RFE-LR model explains biologically the higher leaf area duration during grain filling through site group management the increase of grain yield (Mafruit) with aboveground biomass at maturity (Masec(n)_Imat), Experimental Year Maturity Variety Planting date Water	Soil	Observed variables
	precipitation amount (Raint_cum_Idrp-Imat), cumulative crop transpiration during grain filling Auzeville 000, I, II Blancas, Ecudor, Santana, Isidor, 80, 130 IRR, DRY	Clay, Loam	LAI, BNF, AGPN, AGB, GY, GNC, Oil
	(Ep_cum_Idrp-Imat), the mineral nitrogen accumulated by the plants at the onset of grain filling Sultana, RGT, Shouna
	000, I, II Blancas, Ecudor, Santana, Isidor, ES, (Qnplante_Idrp) and residual leaf area index at maturity (Lai(n)_Imat). All these variables 114, 155 IRR, DRY Pallador, Sultana, RGT, Shouna demonstrate that radiation, water and nitrogen resources are highly representative variables of Béziers I, II Isidor, Sumatra, Ecudor, Fukui, 55, 76, 112 IRR soybean grain yield. Further, the high temperatures may affect crop photosynthesis and grain S109554 filling (MaxTemp_avg_Idrp-Imat) which is in accordance with our knowledge of soybean 0, I, II Sarema, Isidor, Ecudor 67, 96, 132 IRR	Silty, Loam Loam Loam	Clay,	BNF, AGPN, AGB, GY, GNC, Oil, roots depth GY GY
	physiology (Grassini et al., 2021). 0, I, II Sarema, Isidor, Ecudor		76, 103, 131	IRR	Silt, Loam	GY, GNC
	En Crambade	I, II	Isidor, Santana, Ecudor		74, 115	IRR, DRY	Clay	LAI, BNF, AGPN, AGB, GY, GNC, Oil
		I, II	Isidor, Santana, Ecudor		73, 120	IRR, DRY	Clay	LAI, BNF, AGPN, AGB, GY, GNC, Oil
	Mauguio Mondonville 13 MinTemp_avg_Idrp-Imat I, II I, II 14 Ep_cum_Ilev-Imat 0, I, II 15 Days_Swfac_0.6_Iflo-Imat I, II 16 MinTemp_avg_Iflo-Imat I, II 17 Swfac_min_Idrp-Imat I, II 18 Days_MinTemp_18_Idrp-Imat Isidor, Sumatra, Ecudor, Fukui, 56 Days_Inn_0.6_Ilev-Imat 74, 98, 145 S109554 Isidor, Sumatra, Ecudor, Fukui, 61, 92, 138 0.5937 57 Inn_min_Iflo-Imat 0.5842 58 Inn_min_Iflo-Idrp S109554 Sarema, Isidor, Ecudor 80, 102, 124 0.5811 59 Inn_min_Idrp-Imat Isidor, Ecudor 76, 97, 124 0.5654 60 Precip_cum_Idrp-Imat Isidor, Santana, Ecudor 81, 147 0.5638 61 Inn_min_ Ilev-Imat Isidor, Santana, Ecudor 126 0.5618 62 AvgTemp_avg_Ilev-Imat Rivières I, II Isidor, Sumatra, Ecudor, Fukui, S109554 62, 99 19 Swfac_min_Iflo-Imat 0.5600 63 Jul_cum_Ilev-Imat 20 MaxTemp_avg_Idrp-Imat 0.5578 64 Masec(n)_Iflo	IRR IRR IRR IRR IRR, DRY 0.2903 0.2824 0.2684 0.2673 0.2453 0.2405 IRR, DRY 0.2320 IRR 0.1721 0.1698	Clay, Loam Silt, Loam Silt, Loam Silt, Loam Silt, Loam Silt, Loam Clay, Loam	GY GY GY GY, GNC LAI, BNF, AGPN, AGB, GY, GNC, Oil LAI, BNF, AGPN, AGB, GY, GNC, Oil GY
	0, I, II 21 AvgTemp_avg_Iflo-Imat	Sarema, Isidor, Ecudor 0.5551 65 Lai(n)_Iflo	70, 102, 131	IRR	0.1594	Clay, Loam	GY
	I, II 22 Swfac_min_Ilev-Imat	Isidor, Ecudor 0.5521	76, 108, 138 66 MaxTemp_max_Ilev-Imat	IRR	0.1560	Clay, Loam	GY, GNC
	I, II 23 Days_MinTemp_18_Iflo-Imat Isidor, Santana, Ecudor 0.5328 67 Inn_avg_Idrp-Imat 81, 126	IRR, DRY 0.1557	Clay, Loam	LAI, BNF, AGPN, AGB, GY, GNC, Oil
	I, II 24 Days_MaxTemp_28_Idrp-Imat	Isidor, Santana, Ecudor 68 0.5119	77, 126 Days_MaxTemp_28_Ilev-Imat	IRR, DRY 0.1442	Clay, Loam	LAI, BNF, AGPN, AGB, GY, GNC, Oil
	25 Ep_cum_Iflo-Idrp			0.5113	69 AvgTemp_avg_Iflo-Idrp	LAI: Leaf area index (m2.m-2), BNF: fixed nitrogen (kg.ha-1), AGPN: Total nitrogen (kg.ha-0.1030
	1), AGB: Biomass of aerial parts (t.ha-1), GY: Grain yield (t.ha-1), GNC: Nitrogen concentration in the grains (%), Oil: percentage of oil (%), roots depth: rooting depth (cm). Source: (Schoving, 26 Jul_cum_Idrp-Imat 0.5069 70 Zrac_Iflo 0.0934
	2020) . 27 Swfac_avg_Ilev-Imat		0.5054	71 Etpp(n)_cum_Idrp-Imat	0.0858
	28 Qfix_Idrp			0.4953	72 MaxTemp_avg_Iflo-Idrp	0.0856
	29 Swfac_avg_Idrp-Imat		0.4929	73 Qfix_Iflo		0.0821
	30 Swfac_avg_Iflo-Idrp		0.4909	74 Nbgrmax		0.0781
	31 Days_Swfac_0.6_ Ilev-Imat		0.4899	75 MaxTemp_max_Iflo-Imat	0.0773
	32 Days_Swfac_0.6_Iflo-Idrp		0.4834	76 Etpp(n)_cum_Ilev-Imat	0.0708
	33 Raint_cum_Iflo-Imat		0.4712	77 Etpp(n)_cum_Iflo-Imat	0.0701
	34 Days_Swfac_0.6_Idrp-Imat		0.4627	78 Days_MaxTemp_28_Iflo-Idrp	0.0643
	35 Zrac_Idrp			0.4592	79 MaxTemp_max_Iflo-Idrp	0.0618
	36 Zrac_Imat			0.4536	80 MaxTemp_avg_Ilev-Imat	0.0613
	37 Swfac_min_ Iflo-Idrp		0.4464	81 Etpp(n)_cum_Iflo-Idrp	0.0565
	38 MaxTemp_max_Idrp-Imat		0.4257	82 Raint_cum_Iflo-Idrp	0.0461
	39 Jul_cum_Iflo-Imat			0.3997	83 Jul_cum_Iflo-Idrp	0.0456
	40 Days_MinTemp_18_Ilev-Imat	0.3926	84 Stflodrp		0.0372
	41 Qnplante_Iflo			0.3866	85 Stlevdrp		0.0308
	42 Raint_cum_Ilev-Imat		0.3802	86 Stdrpmat		0.0074
	43 Precip_cum_Ilev-Imat		0.3714	87 Days_Inn_0.6_Idrp-Imat	0.0000
	44 Precip_cum_Iflo-Imat		0.3585

Table B1 .

 B1 Ranking of variables computed by Pearson correlation.

	Pos	Variable name	Var	Pos	Variable name	Var
			Importance			Importance
	1	Lai(n)_Imat	0.7699	45 Precip_cum_Iflo-Imat	0.3140
	2	Masec(n)_Imat	0.7311	46 Days_MinTemp_18_Iflo-Idrp	0.3018
	3	Qnplante_Idrp	0.7137	47 Inn_avg_Iflo-Idrp	0.2983
	4	Masec(n)_Idrp	0.7003	48 MinTemp_avg_Ilev-Imat	0.2851
	5	Ep_cum_Iflo-Imat	0.6886	49 Inn_min_Idrp-Imat	0.2817
	6	Ep_cum_Idrp-Imat	0.6734	50 Precip_cum_Ilev-Imat	0.2602
	7	Qnplante_Imat	0.6714	51 Days_Inn_0.6_Ilev-Imat	0.2599
	8	Ep_cum_Ilev-Imat	0.6456	52 Days_Inn_0.6_Iflo-Imat	0.2599
	9	Lai(n)_Idrp	0.6439	53 Days_Inn_0.6_Iflo-Idrp	0.2599
	10 Raint_cum_Idrp-Imat	0.6288	54 Inn_avg_Iflo-Imat	0.2570
	11 Swfac_min_Idrp-Imat	0.6160	55 Precip_cum_Iflo-Idrp	0.2407
	12 Swfac_min_Iflo-Imat	0.6136	56 Inn_min_Iflo-Idrp	0.2370
	13 AvgTemp_avg_Idrp-Imat	0.6081	57 Inn_avg_Ilev-Imat	0.2353
	14 Swfac_min_Ilev-Imat	0.6009	58 MinTemp_avg_Iflo-Idrp	0.2174
	15 MinTemp_avg_Idrp-Imat	0.5977	59 Masec(n)_Iflo	0.2164
	16 Qfix_Imat	0.5964	60 Inn_min_Iflo-Imat	0.2155
	17 MaxTemp_avg_Idrp-Imat	0.5937	61 Lai(n)_Iflo	0.2105
	18 Days_MinTemp_18_Idrp-Imat	0.5876	62 Inn_min_ Ilev-Imat	0.1991
	19 MinTemp_avg_Iflo-Imat	0.5758	63 Jul_cum_Ilev-Imat	0.1979
	20 Swfac_avg_Iflo-Imat	0.5655	64 Days_MaxTemp_28_Ilev-	0.1840
				Imat
	21 AvgTemp_avg_Iflo-Imat	0.5548	65 AvgTemp_avg_Ilev-Imat	0.1831
	22 Days_MinTemp_18_Iflo-Imat	0.5510	66 MaxTemp_max_Ilev-Imat	0.1801
	23 Ep_cum_Iflo-Idrp	0.5482	67 Inn_avg_Idrp-Imat	0.1667
	24 Jul_cum_Idrp-Imat	0.5380	68 Precip_cum_Idrp-Imat	0.1519
	25 Days_Swfac_0.6_Iflo-Imat	0.5361	69 AvgTemp_avg_Iflo-Idrp	0.1400
	26 Days_MaxTemp_28_Idrp-	0.5033	70 Qfix_Iflo	0.1263
		Imat
	27 Qfix_Idrp	0.4679	71 Nbgrmax	0.1160
	28 Jul_cum_Iflo-Imat	0.4667	72 Etpp(n)_cum_Idrp-Imat	0.1095
	29 MaxTemp_max_Idrp-Imat	0.4629	73 Raint_cum_Iflo-Idrp	0.1015
	30 Raint_cum_Iflo-Imat	0.4590	74 MaxTemp_max_Iflo-Imat	0.0949
	31 Days_Swfac_0.6_ Ilev-Imat	0.4546	75 Jul_cum_Iflo-Idrp	0.0849
	32 Swfac_avg_Ilev-Imat	0.4364	76 MaxTemp_avg_Ilev-Imat	0.0763
	33 Swfac_avg_Idrp-Imat	0.4319	77 Etpp(n)_cum_Ilev-Imat	0.0763
	34 Raint_cum_Ilev-Imat	0.4209	78 Etpp(n)_cum_Iflo-Idrp	0.0736
	35 Qnplante_Iflo	0.4206	79 Etpp(n)_cum_Iflo-Imat	0.0731
	36 Days_Swfac_0.6_Idrp-Imat	0.4148	80 Zrac_Iflo	0.0720
	37 Swfac_avg_Iflo-Idrp	0.4131	81 Days_MaxTemp_28_Iflo-Idrp	0.0538
	38 Swfac_min_ Iflo-Idrp	0.4105	82 Stflodrp	0.0465
	39 Days_Swfac_0.6_Iflo-Idrp	0.3993	83 MaxTemp_avg_Iflo-Idrp	0.0345
	40 Days_MinTemp_18_Ilev-Imat	0.3857	84 MaxTemp_max_Iflo-Idrp	0.0235
	41 Days_MaxTemp_28_Iflo-Imat	0.3552	85 Stlevdrp	0.0085
	42 Zrac_Idrp	0.3287	86 Stdrpmat	0.0085
	43 MaxTemp_avg_Iflo-Imat	0.3273	87 Days_Inn_0.6_Idrp-Imat	0.0000
	44 Zrac_Imat	0.3214

Table B2 .

 B2 Ranking of variables computed by Spearman correlation filter.

	Pos	Variable name	Var	Pos	Variable name	Var
			Importance			Importance
	1	Lai(n)_Imat	0.5929	45 Inn_min_Iflo-Imat	0.1745
	2	Masec(n)_Imat	0.5875	46 Inn_avg_Idrp-Imat	0.1744
	3	Qnplante_Imat	0.5515	47 MinTemp_avg_Iflo-Imat	0.1687
	4	Ep_cum_Idrp-Imat	0.4223	48 Inn_min_Idrp-Imat	0.1673
	5	Qfix_Imat	0.4072	49 Days_Swfac_0.6_ Ilev-Imat	0.1628
	6	Inn_min_Iflo-Idrp	0.3799	50 Days_MaxTemp_28_Ilev-	0.1595
				Imat
	7	Inn_avg_Iflo-Idrp	0.3652	51 Lai(n)_Iflo	0.1540
	8	Ep_cum_Iflo-Imat	0.3486	52 Precip_cum_Ilev-Imat	0.1213
	9	Days_Swfac_0.6_Iflo-Imat	0.3443	53 AvgTemp_avg_Ilev-Imat	0.1204
	10 Inn_avg_Ilev-Imat	0.3400	54 MaxTemp_avg_Ilev-Imat	0.1204
	11 Inn_avg_Iflo-Imat	0.3380	55 Etpp(n)_cum_Iflo-Imat	0.0993
	12 Masec(n)_Idrp	0.3300	56 Etpp(n)_cum_Ilev-Imat	0.0993
	13 Qnplante_Idrp	0.3285	57 Etpp(n)_cum_Iflo-Idrp	0.0993
	14 Raint_cum_Idrp-Imat	0.3080	58 Masec(n)_Iflo	0.0000
	15 Days_Swfac_0.6_Iflo-Idrp	0.3025	59 Qfix_Iflo	0.0000
	16 Precip_cum_Iflo-Imat	0.2988	60 Zrac_Iflo	0.0000
	17 Qnplante_Iflo	0.2895	61 Etpp(n)_cum_Idrp-Imat	0.0000
	18 Ep_cum_Ilev-Imat	0.2833	62 Precip_cum_Idrp-Imat	0.0000
	19 Swfac_avg_Iflo-Imat	0.2758	63 MaxTemp_avg_Iflo-Imat	0.0000
	20 Lai(n)_Idrp	0.2711	64 MaxTemp_max_Iflo-Imat	0.0000
	21 Swfac_min_Iflo-Imat	0.2476	65 Raint_cum_Ilev-Imat	0.0000
	22 Swfac_min_Ilev-Imat	0.2476	66 MinTemp_avg_Ilev-Imat	0.0000
	23 Zrac_Idrp	0.2457	67 Swfac_avg_Ilev-Imat	0.0000
	24 Zrac_Imat	0.2457	68 MaxTemp_max_Ilev-Imat	0.0000
	25 MinTemp_avg_Idrp-Imat	0.2374	69 Inn_min_ Ilev-Imat	0.0000
	26 Swfac_min_ Iflo-Idrp	0.2309	70 Jul_cum_Ilev-Imat	0.0000
	27 Swfac_min_Idrp-Imat	0.2289	71 Raint_cum_Iflo-Idrp	0.0000
	28 AvgTemp_avg_Iflo-Imat	0.2266	72 AvgTemp_avg_Iflo-Idrp	0.0000
	29 MaxTemp_avg_Idrp-Imat	0.2222	73 MinTemp_avg_Iflo-Idrp	0.0000
	30 Days_MinTemp_18_Ilev-Imat	0.2214	74 MaxTemp_avg_Iflo-Idrp	0.0000
	31 Jul_cum_Iflo-Imat	0.2191	75 MaxTemp_max_Iflo-Idrp	0.0000
	32 Days_MinTemp_18_Iflo-Imat	0.2063	76 Jul_cum_Iflo-Idrp	0.0000
	33 AvgTemp_avg_Idrp-Imat	0.2049	77 Days_Inn_0.6_Ilev-Imat	0.0000
	34 Days_MinTemp_18_Idrp-Imat	0.2035	78 Days_MaxTemp_28_Iflo-	0.0000
				Imat
	35 MaxTemp_max_Idrp-Imat	0.2000	79 Days_Inn_0.6_Iflo-Imat	0.0000
	36 Qfix_Idrp	0.1980	80 Days_Inn_0.6_Idrp-Imat	0.0000
	37 Ep_cum_Iflo-Idrp	0.1971	81 Days_MinTemp_18_Iflo-Idrp	0.0000
	38 Jul_cum_Idrp-Imat	0.1969	82 Days_MaxTemp_28_Iflo-Idrp	0.0000
	39 Days_MaxTemp_28_Idrp-Imat	0.1954	83 Days_Inn_0.6_Iflo-Idrp	0.0000
	40 Raint_cum_Iflo-Imat	0.1927	84 Nbgrmax	0.0000
	41 Swfac_avg_Idrp-Imat	0.1883	85 Stlevdrp	0.0000
	42 Days_Swfac_0.6_Idrp-Imat	0.1883	86 Stflodrp	0.0000
	43 Precip_cum_Iflo-Idrp	0.1872	87 Stdrpmat	0.0000
	44 Swfac_avg_Iflo-Idrp	0.1829

Table B3 .

 B3 Ranking of variables computed by filter of entropy-based information gain.

Acknowledgments

This research was supported by Horizon 2020 SusCrop-ERA-NET Cofound on Sustainable Crop Production: LegumeGap project (2019-2021) "Increasing productivity and sustainability of European plant protein production by closing the grain legume yield gap": https://www.suscrop.eu/projects-first-call/legumegap