IntroductionTop

Soybean (Glycine max (L.) Merrill) is one of the most economically significant crops worldwide (Kulcheski et al., 2016), and multiple linear regression models are constantly being developed to partially explain yield variations in that crop. When modeling soybean yield, variables concerning agricultural meteorology (Penalba et al., 2007; Tao et al., 2008), agriculture (Zheng et al., 2009), management (Lobell et al., 2005), vegetation index (Mercante et al., 2010) and soil parameters (Garcia-Paredes et al., 2000) are often utilized. Models differ according to the nature of the explanatory variables used in the modeling process, for which a review on model categories in use can be found in Vera-Diaz et al. (2008).

Considering that determining the values of certain variables is often a burdensome and arduous task, in some cases analyses are carried out on small samples[1]. This fact may call into question the inferences being made, since traditional inference methods are asymptotic and standard errors and confidence intervals may be biased in small samples, as explained by Hao & Naiman (2010). Adopting more parsimonious models and determining influential points, are procedures that can also provide misleading results when working with small sample sets. Kamo et al. (2013) explain that the Akaike information criterion - AIC (Akaike, 1973) used for model selection presents a bias that cannot be ignored, especially with small samples, given that it is derived from asymptotic properties. Regarding the diagnostic measures of overall influences one problem is related to its cutoff points. According to Martin & Roberts (2010) they are based on large sample theory and therefore may not be suitable for small samples.

An alternative to traditional inference methods is the use of the bootstrap, a simulation method developed by Efron (1979) which uses resampling with replacement of the sample data set to perform statistical inferences such as hypothesis testing and determination of confidence intervals (Dubreuil et al., 2014). The bootstrap method has applications in regression analysis (Rahman, 2014), model selection (Al-Marshadi, 2011) and definition of global influence diagnostics (Beyaztas & Alin, 2013).

By comparing the results obtained from bootstrap methods with results of asymptotic methods, Chaves-Neto & Faria (2015) conclude that bootstrap performed well in samples of all sizes and was higher than the asymptotic method in small samples. Although bootstrap is a well-known technique and is frequently employed in agricultural studies – as seen in works by Sabaghnia et al. (2010), García-Gallego et al. (2015), Losada et al. (2015) and Sutton et al. (2016) – the development of statistical and computing models has led to the study of new techniques based on the bootstrap method.

The objective of this work was to utilize bootstrap methods to select explanatory variables, investigate the existence of influential points through diagnostic analysis, and obtain confidence intervals for the parameters of a multiple linear regression model for soybean yield considering physical and chemical soil properties as explanatory variables.

Material and methodsTop

ResultsTop

Descriptive statistics of the explanatory variables indicated homogeneous behavior of the variables, with no multicollinearity found. The multiple linear regression model of soybean yield, estimated through OLS considering all explanatory variables (Eq. [3]), showed an adjusted coefficient of determination (R²_Adj) of 0.41 and root mean square error (RMSE) of 0.33.

It could be observed that estimates for those parameters associated with SRP₁, SRP₃, Des₁ and Des₂ variables showed negative signs, indicating that an increase in the value of these variables implies a reduction in soybean yield (Eq. [3]). The parameters estimation associated to SRP₂ and Des₃ variables from the Eq. [3] showed different signals from the expected scenario, since it indicates a direct relation from such variables towards soybean productivity (Eq. [3]). The positive estimate signal, from the associated parameter of the variable K, indicates that, while maintaining other variables constant, an increase in one unit in the K variable produces an increase in soybean productivity, at a rate of 1.208 t/ha (Eq. [3]). The bootstrap intervals were determined with reliability of 95% for the parameters of the multiple linear regression model using the techniques of bootstrap percentile by Efron and BC bootstrap (Table 1).

It was observed that the vast majority of the confidence intervals, determined by the bootstrap technique, contained zero indicating, that with exception of the variable P, the other explanatory variables may not be individually significant. In search for a more appropriate multiple linear regression model it was applied the model selection method using bootstrap considering 1000 resamples (Table 2). It was observed that, of the 1000 models for which bootstrap resamples had been adjusted, by applying the backward selection method with statistical Akaike – AIC to each of them, the result showed that in 91% of the models the predictor variable P was selected, indicating that phosphorus is an important soil attribute for soybean yield prediction. Furthermore, it was observed that in 100% of models in which phosphorus had been selected, its estimated parameter was negative, which ensures that when other variables are held constant an increase in phosphorus level implies reduction of soybean yield.

Other variables selected for most models were Des₂ with a selection percentage of 87%, Ca with 81% and SRP₁ with 79%. Analyzing the signs of the estimated parameters associated with these variables in the models in which they were selected it is highlighted that in 94% of models in which the Ca variable was selected the sign of its estimated parameter was positive, suggesting the increase in value of this variable contributes for increasing soybean yield. For those estimated parameters associated with SRP₁ and Des₂ variables, in 98% of models in which they were selected the signals were negative. It is clear that some variables may not be useful to explain soybean yield behavior. For example, among the 1000 models obtained, the Des₁ variable was selected in only 500 and additionally for 180 of those the estimated parameter sign was positive and for 320 of those the sign was negative, thus, this set of oscillations is a guarantee that this variable is not significant, therefore, can be deleted without causing damage to the modeling. A similar case occurs with SRP₃ variable, as well as being selected in only 460 models, the appropriate sign of its estimated parameter cannot be identified considering that 230 models had a positive sign and 230 had a negative sign. As per the parameters estimates associated to SRP₂ and Des₃ variables, it showed opposite signals from the expected scenario (Eq. [3]), it is desirable to verify the importance of such variables for modeling purposes. Although the positive signals from the associated estimated parameters to such variables appear in a great part of the models (80% and 78%, respectively), the selection percentages were not very elevated (61% and 42%, respectively) and, therefore, there were evidences that they were not significant and could be removed from modeling (Table 2). In view of these observations four models were set to be analyzed, namely M₈₁, M₇₉, M₇₅ and M₇₁ (Table 3), and each of them was determined according to a number of explanatory variables selected in accordance with how many times they had been selected in the bootstrap models.

Regressors present in the M₈₁ model can explain only 37% of the soybean yield variation, a result lower than that obtained when considering the model containing all the explanatory variables (R²_Adj = 0.41). The M₇₅ (R²_Adj = 0.42) and M₇₁ (R²_Adj = 0.49) models provided a greater degree of explanation between the explanatory variables and soybean yield than the full model, while the M₇₉ (R²_Adj = 0.41) model provided an equivalent level of explanation, however, these models had a higher RMSE compared to the complete model (RMSE = 0.33) and that difference is most evident in the M₇₉ model (RMSE = 0.39). As the M₇₁ model explained 49% of the soybean yield variation and RMSE of this model (RMSE = 0.34) is close to RMSE of the complete model (RMSE = 0.33) the M₇₁ model was chosen as best adjusted model to soybean yield and analysis was performed using JaB to investigate the existence of influential points.

It is noteworthy to mention no points were detected as influential when value 1 is established as the cutoff point (Fig. 2). The same is true when considering the criteria that detects point i as influent if D_i is higher than the median of the distribution F of Snedecor with freedom degrees of p = 6 and n - p = 24 once the cutoff point is 2.50 to these, thus they were also not detected as influential points. Considering 4/n ≈ 0.13 as cutoff point, the points 15, 23 and 29 were detected as influential indicating these points can change the estimation of the parameters in the regression model, so it is important to investigate the model behavior without the use of these points. It should be emphasized that only point 23 was detected as being influential through analysis using Cook’s distance (D_i) with JaB methodology.

Figure 2. Determination of influential points using Cook’s distance (Di) with JaB methodology. Points that cross the dashed line are considered influent according to Cook’s distance (Di), while points outside the bootstrap interval are considered influent from jackknife-after-bootstrap (JaB) analysis.

JaB graphswere created to help identify influential points, they give a visual interpretation of how a particular point affects the bootstrap distribution for the estimation of parameters in M₇₁ (Fig. 3). Observing the graphs in Fig. 3, it was noted that points 10, 15, 23 and 29 were detected as influent.

Figure 3. Jackknife-after-bootstrap (JaB) plots for the parameter estimates associated with the explicative variables SRP1: soil penetration resistances from 0 to 0.1 m (a), Ca: calcium (b), Mg: magnesium (c), P: phosphorus (d), Mn: manganese (e) and Des2: soil densities from 0.1 to 0.2 m (f) of the model containing the variables selected in at least 71% of the bootstrap models (M71). ui↑ {s}: Values of the relative jackknife influence function in growing order. Pk: k-th percentile of the bootstrap distribution, k = {5,10,16,50,84,90,95}. The vertical dashed lines that pass through values -2 and 2 represent the cutoff interval. The points allocated outside this interval are shown in the upper part of each graph and their relative jackknife influence values are highlighted with vertical dashed lines. The red lines represent the variations in these percentiles and the intersections with the blue dotted line representative of point i represent the percentiles obtained when the bootstrap distribution is formed with the bootstrap replicates calculated from the bootstrap samples that do not contain point i.

Two new models were adjusted to the variables P, Des₂, Ca, SRP₁, Mn, Mg as to measure the effect of influential points in modeling. The M_{71-{15,23,29}} model was adjusted to the data set without points (15, 23, 29), as these were detected as influential by traditional Cook distance method with cutoff point of 4/n. The M_{71-{10,15,23,29}} was also adjusted to the data set without points (10, 15, 23, 29) for these were considered as influential by analysis using JaB (Table 4).

The M_{71-{15,23,29}} model adjusted to data set without sample elements 15, 23 and 29, which were identified as influential by the traditional method, is more explicative than M₇₁ model attained from the complete set of points for after removal of these points the percentage of soybean yield variation that can be explained by the regressors increased from 49% to 63% (Table 4). When considering M_{71-{10,15,23,29}} model prepared without points 10, 15, 23 and 29 it was observed that the adjusted coefficient of determination (0.65) was higher than the adjusted coefficient of determination obtained from M_{71-{15,23,29}} model, resulting in a more explanatory model. When comparing RMSE of M_{71-{15,23,29}} and M_{71-{10,15,23,29}} models (Table 4) it was emphasized that the identification of point 10 as influential and its withdrawal from the data set as a result corroborated reduction of this statistic, thus resulting in a more accurate model for making predictions. In view of these results it was decided to choose the M_{71-{10,15,23,29}} model as the best model suited to soybean yield and determine bootstrap confidence intervals for the parameters associated with the explanatory variables (Table 5).

By comparing the confidence intervals of parameters of the M_{71-{10,15,23,29}} model with the respective intervals obtained from the multiple linear regression model generated with all the explanatory variables and all sampling points (Table 1) it can be observed that regardless of the bootstrap method used the confidence intervals of parameters for the M_{71-{10,15,23,29}} model had lower amplitude, indicating estimates of this model was more accurate.

DiscussionTop

The average soybean productivity in the monitored area (4.305 t/ha) is considered high compared with other regions, according to data from CONAB (2015) in the agricultural year 2013/2014 average productivity in Brazil was 2.854 t/ha and in Paraná was 2.950 t/ha.

The negative sign of estimates for parameters associated with SRP₁, SRP₃, Des₁ and Des₂ variables (Eq. [3]) are expected once soil density (Des) shows a direct relationship with SRP (Busscher et al., 1997), and as SRP has great influence on plant growth, root growth and crop yields vary inversely proportional to its value (Freddi et al., 2006). The parameters estimation associated to SRP₂ and Des₃ variables show opposite signals from the expected scenario; however, since it is verified that multicollinearity was non-existent, it is also prudent to investigate the significance of such variables. The positive estimate signal from the associated parameter to K variable is expected, once and in accordance with Pettigrew (2008), potassium is one of major nutrients considered essential for crop growth and yield development.

The comparison of confidence intervals can be done in terms of their amplitudes according to Paes (1998), to whom a high amplitude interval indicates reduced accuracy of estimation as compared with a range of lower amplitude, thus comparing the two techniques bootstrap confidence intervals (Table 1) it is clear the intervals obtained by Efron percentile technique showed lower amplitude and therefore is the most accurate. Given the fact that zero is present in most of the confidence intervals (Table 1) it is prudent to investigate whether there are irrelevant variables and/or influential points in the data set for they cause an increase in the parameter variance (Rao, 1971; Meloun & Militký, 2001) and as a result confidence intervals tend to have a greater range and loss of accuracy.

The fact that predictor variable P is selected in a large share of models (Table 2) and that its signs of the estimated parameters are negative in all of them can be explain by the high phosphorus values ​​found (on average 12 mg/dm³) which according to Popp et al. (2002), may indirectly decrease yields due to micronutrients imbalance. The high percentage (94%) of times when the sign of the estimated parameter associated with variable Ca is positive is also expected, because calcium deficiency is among the main factors that inhibit root growth as reported by Oliveira et al. (2009), especially in latosols. Such a deficiency would have the plant vulnerable to biotic, biological and nutritional stresses and consequently would lead to reduced productivity (Dourado Neto et al., 2014). As SRP₁ and Des₂ variables are used to assess the state of soil compaction, their effect on soybean yield is the opposite, for plants exhibit alterations in depth, branch and distribution of roots in response to soil compaction (Rosolem et al., 2002), which undermines the efficient use of nutrients and water and limits crop yield (Alakukku & Elomen, 1995).

The model selection method using bootstrap is effective in determining the significant variables resulting in a more parsimonious model. Although the model determined by this method (M₇₁) has been the same selected by the conventional method using Akaike, the application of this methodology serve to attest the model selected by the Akaike criterion is not super-parameterized, which can occur when the amount of samples is small.

Analyzing Fig. 3a, the graph of the bootstrap distribution of the parameter estimates associated with the variable SRP₁, it is seen that points 15 and 10 are detected as influential. Point 15 has a negative influence (-3.7) and its removal reduces bootstrap distribution amplitude, a fact that occurs mainly due to a shift in the initial percentiles if one considers the empirical distribution formed with 3000 replicates, P₅ = -0.373, P₁₀ = -0.330, P₁₆ = -0.302 and considers the empirical distribution formed only by bootstrap replicates with bootstrap samples not containing point 15 (1124 samples), P₅ = -0.336, P₁₀ = -0.295, P₁₆ = -0.270. The influence of point 10 is positive (2.6). It is observed that when considering the bootstrap distribution formed with those bootstrap samples that do not contain point 10 (1039 samples) the values considered percentile decrease, causing distribution displacement and reduction of its range from 0.865 to 0.727.

JaB graph in Fig. 3b for Ca variable indicated point 23 as negative influence (-2.8) and point 15 as positive influence (3.5), thus withdrawal of these points also causes changes in the empirical distribution of bootstrap estimates. After disregarding bootstrap replicates obtained from bootstrap samples which had point 23, the initial percentiles increased and the distribution range went from 1.214 to 0.999; and after disregarding those replicas obtained from samples containing point 15, there was a reduction in values of final percentiles, which also reduces the amplitude of the empirical distribution. Analyses of other graphs (Figs. 3c through 3f) are similar and indicate point 15 has a negative influence on bootstrap distributions of the parameters associated with the Mg, Mn and Des₂ variables; it also indicates point 23 has a positive influence on bootstrap distributions of the parameters associated with Mg and Mn variables and a negative influence on bootstrap distribution of the parameter associated with variable P. Point 29 has a positive influence on bootstrap distribution of the parameter associated with variable P and point 10 has a positive influence on bootstrap distribution of the parameter associated with variable Des₂.

Comparing all the points that were detected as influential in JaB graphs (Fig. 3), it is clear the sampling member 15 stands out due to its influence on most bootstrap distributions of estimated parameters, only the distribution of the parameter associated with the variable P is not influenced by excluding this element. The sample elements 23 and 10 also stood out as influential on various confidence intervals and the sample element 29 is the least influential of the four. By taking sample element 29 out of the empirical distribution of the bootstrap replicates obtained from bootstrap samples only the bootstrap distribution of parameter estimates associated with the variable P is influenced, seeing its range reduced from 0.163 to 0.112.

Regarding diagnostic analysis it is seen that the influential points determination method using JaB methodology together with Cook’s distance (Fig. 2) do not identify some points clearly highlighted as influential by traditional methods and JaB graphics. Thus, JaB graphics prove to be a great alternative to identify influential points, as well as identifying the influential points with greater accuracy compared to traditional analysis they also provide information on bootstrap distributions of parameter estimates, making it possible to see what happens to confidence intervals when the influential samples are excluded.

Regarding the significance of the explanatory variables in M_{71-{10,15,23,29}} model it is observed that only the parameter associated with Mn variable showed confidence intervals containing zero in both bootstrap techniques (Table 5) giving signs it may be irrelevant. This suspicion can be ruled out once there is evidence that the parameter signal associated with this variable is negative. It is necessary to simply notice that the zero has appeared at the high end of intervals and to also notice that in the variable selection method (Table 2) variable Mn is selected in 75% of models (750 models) and has a negative sign in 94% of them (705 models). By comparing techniques of confidence interval determination used it is observed they presented a similar behavior, though Efron percentile method stands out for providing intervals of lower amplitude. Cunha & Colosimo (2003) also highlights Efron percentile method to determine confidence intervals for regression models with measurement errors, since according to the authors this method is evidenced by its greater simplicity with equal performance compared to the others.

It is noteworthy that the bootstrap methods are fundamental to obtaining a more descriptive and more accurate model, as aside from the model M_{71-{10,15,23,29}} to furnish a higher percentage of explanation of the soybean productivity (65%) than the initial model in Eq. [3] (41%), it furnishes a lower RMSE, being more accurate. It is important to highlight the explanatory influence of the model M_{71-{10,15,23,29}} to be under satisfactory terms taking into account to be built only by physical and chemical features of the soil. The soybean productivity percentage variation not covered by such model (35%) is due to variables not considered, for example, the agricultural-meteorological, since climate has a significant impact upon the growth and development of crops (Hoogenboom, 2000). It is valuable to note the non-inclusion of the agricultural-meteorological variables in this essay to be due to the fact of the limitation on the spatial representation of the results to be obtained from the collected data in weather stations. (Junges & Fontana, 2011).

The results showed that the bootstrap methods enabled us to select the physical and chemical soil properties, which were significant in the construction of the soybean yield regression model, construct the confidence intervals of the parameters and identify the points that had great influence on the estimated parameters.