Spatial variability of crop responses to agronomic inputs in on-farm precision experimentation

Within-field variability of crop yield levels has been extensively investigated, but the spatial variability of crop yield responses to agronomic treatments is less understood. On-farm precision experimentation (OFPE) can be a valuable tool for the estimation of in-field variation of optimal input rates and thus improve agronomic decisions. Therefore, the objectives of this study were to investigate the spatial variability of optimal input rates in OFPE and the potential economic benefit of site-specific input management. Mixed geographically weighted regression (GWR) models were used to estimate local yield response functions. The methodology was applied to investigate the spatial variability in corn response to nitrogen and seed rates in four cornfields in Illinois, USA. The results showed that spatial heterogeneity of model parameters was significant in all four fields evaluated. On average, the RMSE of the fitted yield decreased from 1.2 Mg ha−1 in the non-spatial global model to 0.7 Mg ha−1 in the GWR model, and the r-squared increased from 10 to 68%. The average potential gain of using optimized uniform rates of seed and nitrogen was US$ 65.00 ha−1, while the added potential gain of the site-specific application was US$ 58.00 ha−1. The combination of OFPE and GWR proved to be an effective tool for testing precision agriculture’s central hypothesis of whether optimal input application rates display adequate spatial variability to justify the costs of the variable rate technology itself. The reported results encourage more research on response-based input management recommendations instead of the still widespread focus on yield-based algorithms.


Introduction
Site-specific technologies, including yield monitoring, remote sensing imaging, and variable rate input application, have become increasingly available to farmers in recent decades. Additionally, many farmers have access to software tools to process this information and use it to guide site-specific management decisions. This decision-making process is based primarily on the knowledge about yield response to crop input management in agronomic trials, which involve changing management practices and subsequent monitoring of the effects on the system output (Pringle et al. 2004a, b).
The potential for economic gains from the implementation of site-specific crop management technologies depends on the assumption that yield response to an agronomic input can be described as a function of managed input strategies, field characteristics, and weather. Because field characteristics are spatially dependent, how yields respond to managed inputs changes over space the shape of the traditional "yield curve" plot with input application rate on the horizontal axis and yield on the vertical axis will vary among sites within a field (Bullock and Bullock 1994). However, most variable rate recommendations are calculated using the same methods, and calibrations developed for uniform field management. These are based on agronomic trials that are usually planned to be representative of the conditions in the geographical region of interest to the researchers aiming to infer the management insights coming from small-scale plots to larger regions. To make the recommendations more accessible to users, these models are usually very simplified versions of the true yield response function ). In the case of nitrogen fertilizer, this simplification has led to the widespread use of yield levels as a proxy for estimating the optimal input application rates, even though academic research provides only weak evidence of any correlation between yield levels and economically optimal input application rates (Bachmaier and Gandorfer 2009;Rodriguez et al. 2019;Scharf et al. 2006).
To improve site-specific management requires information about how crop yield responds to varying treatments and how those responses vary over space (Bullock and Bullock 1994). Variable-rate technology can be used to systematically control input levels in highly mechanized, large-scale production systems (Piepho et al. 2011). In addition, these operations make possible the running of large-scale, on-farm precision experiments (OFPE), that generate large amounts of site-specific response data, which can be used to understand the spatial variation of optimal input application rates . Moreover, because OFPE data are gathered in the same fields for which management recommendations are desired, field and site-specific yield response functions can be estimated. Thus, the ultimate purpose of OFPE is to develop site-specific input application (Piepho et al. 2011). OFPEs also allow testing the fundamental hypothesis of precision agriculture, which is that the rate at which inputs are applied can be profitably varied within fields to match site-specific requirements (Bachmaier and Gandorfer 2009;Lark and Wheeler 2003). Of course, estimating site-specific yield response requires spatial data analysis (Bullock and Lowenberg-DeBoer 2007;Hurley et al. 2005).
Most of the examples of OFPE consider only the effect of a single factor on crop yield (Kindred et al. 2017;Piepho et al. 2011;Pringle et al. 2004a, b). The statistical analysis of this type of trial often involves using geostatistical interpolation methods to estimate the effects of all tested treatment levels in all points of a regular grid. With that, yield estimates for each treatment level are obtained not only for the treatment tested on that point but for all other treatments as well. Then local response function parameters need to be estimated using the interpolated values for each point in the interpolated grid. The confidence of each estimate will depend on the distance of neighboring experimental units with the same treatment level, and for that, systematic designs are preferable over randomized designs (Pringle et al. 2004a, b).
The main limitation of these geostatistical methods is that they rely on interpolating yield maps using only one level of the treatment. If the treatment has five levels, for example, only 20% of the data is used at every interpolation. In a factorial design with four levels of the first treatment and five levels of the second, totaling twenty combinations, only 5% of the observations would be used in each one of the twenty interpolations. Due to the suboptimal use of the neighboring plots in the estimation of the yield response function parameters, the geostatistical interpolation is not well suited for the analysis of continuous variables and factorial designs with many levels (Pringle et al. 2004a, b).
Most spatial inference methods, including most generalized least square and spatial econometrics methods, were developed to focus on the implications of spatial dependence (Anselin 2010). However, spatial heterogeneity has been overlooked, even though it is crucially important to model spatial data appropriately (Geniaux and Martinetti 2018;Murakami and Griffith 2019).
To model the spatial heterogeneity as local yield responses, it is necessary to have some local variable as an input, either as part of the design or some covariable that was measured in the field (Bachmaier and Gandorfer 2009;Thöle et al. 2013). Besides the need for additional data, the main limitation of this approach is that the results will be dependent on the correlation between the covariables and the yield response. The hypotheses that can be tested with such a model are restricted to whether the yield response function has any interaction with the spatial covariables, but it is not possible to test whether there is significant spatial variability in the parameters of the yield response function itself. Since there will always remain variables that can not be observed, drawing conclusions based only on the known variables may lead to wrong generalizations.
In recent years, spatially varying coefficient (SVC) models have attracted considerable attention in various fields of applied sciences. However, their use in agronomic research is still limited (Cai et al. 2014;Murakami et al. 2018;Trevisan et al. 2019a). These methods have been proposed to investigate the spatial variability or no stationarity of coefficient estimates in regression models. In other words, SVC methods are equivalent to the direct estimation of the site-specific yield production function parameters. The estimates, along with associated inference diagnostics, can be mapped to the original measurement locations or a new set of locations. Among the SVC models, the geographically weighted regression (GWR) has been one of the methods commonly applied. In the GWR, the parameters of the models vary in space and can be mapped and interpreted as a spatial variable (Brunsdon et al. 1996;Fotheringham 1997).
The GWR method allows the specification of more complex and less restrictive models, in which all parameters are estimated for each location. The neighboring points are included in the estimation even if they are from different treatment levels by applying weighted least squares estimation to neighboring subsamples. Weights are estimated via a distance-decay kernel, similar to the weights given by the semivariogram to neighbors in kriging interpolation.
Therefore, this study proposes to integrate the recent advances in precision agriculture and methods of spatial analysis to develop a new methodology for understanding the spatial variability of crop responses, and thereby test the viability and profitability of precision agriculture. The novelty of the proposed methodology is the use of GWR as a data analysis technique applied to OFPE, without the need for spatial covariates. The hypothesis being tested is that yield response function with spatially varying coefficients is adequate relative to the alternative of using a yield response function with only global parameters. The specific objectives of this study were to investigate the spatial variability of optimal input application rates, and the potential economic benefit of joint use of OFPE and sitespecific management.

OFPE design
The fields used to generate the datasets for this study come from Data-Intensive Farm-Management (DIFM) project on-farm field trials. DIFM is a multi-university project supported by the United States Department of Agriculture-National Institute of Food and Agriculture (USDA-NIFA) to use precision agriculture technology to generate original, high-quality, full-field, on-farm trial data at low cost . Four Illinois fields representing typical U.S. Corn Belt maize production systems were used for the trials (Table 1), with Fields 1 and 2 hosting trials in 2017, and Fields 3 and 4 in 2018.
The OFPEs were designed to generate data for localized estimation of crop response to guide site-specific crop management (Piepho et al. 2011). Each experiment had two managed input factors, seeding rate (S) and nitrogen fertilizer application rate (N), with at least four levels of each factor. The allocation of the rates in the experimental units followed a completely randomized factorial design replicated over the entire field.
All treatment levels were implemented in the field using variable rate planters and fertilizer applicators. Each experimental unit ("plot") was given a dimension to fit the swath width of the machinery available (Table 1). The ranges of variation for the tested rates were chosen according to the producer's experience and expectations for the field, typically allowing a 20% variation in each direction around the status quo rate (Table 2). Other farming practices were kept constant throughout the field and were conducted by the farmers by standard protocols for the region.

Data analysis
Treatment applications were monitored using the feedback sensors in the variable rate applicators, generating the as-applied data. The as-applied data was filtered to remove extreme values, defined as values distant from the average more than three times the interquartile range. Yield data were collected during harvest using combine yield monitoring systems. Yield data were filtered to remove global and spatial outliers, taking into account the boundary of the plots (Trevisan et al. 2019b). The data from the headlands and borders of the field were also discarded. Each plot was subdivided into ten to twenty smaller regular polygons, depending on the plot's original length (Fig. 1). The area of observational units was constant within each field and ranged from 37 to 87 m 2 among fields because of differing machinery dimensions. The average value of all points within a polygon was used as the observational unit in the final dataset, which in each field consisted of approximately 3000 observations (Table 1).
All steps and models were individually applied to each field. Letting I represent the exact number of observations in a field, with observations are indexed by i = 1, …, I. Let the (longitude, latitude) geographic coordinates of the centroid of an observational unit (a "point") be denoted (l i , t i ). Based on weighted least squares, the GWR method modifies the estimates of the regression coefficients as a function of the geographic position (l i , t i ). In matrix notation, the regression coefficients were estimated by Eq. 1: Table 2 Treatment rates and yield of the four Illinois cornfields used in the OFPE trials *Base N nitrogen rate applied at a uniform rate, at or before planting. Total N total nitrogen rate consisting of the base N added to the experimental rates applied at side-dressing. Yield grain yield estimated by the yield monitor system; count represents the number of raw sensor observations, followed by the field average yield and the standard deviation within parenthesis where Y is the (I × 1) vector whose components are the Yield i , and X i is the (I × 1) vector of the values of the J explanatory variable at the location (l i , t i ). X is the (I ×J ) matrix in which X i is the ith column, for i = 1, …, I. In this notation, the GWR model incorporates local variation in the estimated coefficients b(l i , t i ) using a geographic weighting matrix G(l i , t i ). This diagonal matrix incorporates the distance d ij between predictors X i at (l i , t i ) and dependent observations Y j at the location (l j , t j ). As d ij increases, the explanatory variables are expected to decrease in influence over the response variable (Plant 2012).
One particular case of GWR is the moving window regression, which applies ordinary least squares estimation to neighboring subsamples by specifying a kernel weighting scheme in which all weights are equal to one within the kernel and zero otherwise (Lloyd 2009). Using distance-decay weighting provides added flexibility to local regression modeling, allowing more data to have local influence, and tends to yield more smoothly varying coefficient surfaces. This decay is captured by elements g ij of a matrix G, for example, as would be defined by the Gaussian spatial kernel function g ij in Eq. 2: Function g ij calculates the weights given to the regression points at each location, with a bandwidth h adjusted by cross-validation and selected based on the corrected Akaike information criterion (AICc) ). In addition to the Gaussian, the exponential and the bi-square spatial kernels were also compared. The spatial kernel can use a fixed distance for the bandwidth h, or it can be adapted at each location to include a minimum number of neighbors. Although the observations in the final dataset were regularly spaced in the field, an adaptive distance bandwidth was chosen to improve the results at the borders of the experiment, where a fixed kernel would include fewer observations. An example of the weights at one location is presented in Fig. 2. Fig. 2 Illustration of the weights given to neighbor points for the parameter estimates at the center of the field using a Gaussian kernel. NR stands for nitrogen rate and represents the input variation according to the trial design Literature often suggests that non-linear functions such as the quadratic-plateau function, are more suited to describe yield response in corn than are their linear-in-variables counterparts (Pahlmann et al. 2016;Scharf et al. 2005). Although nonlinear-in-variables models may be necessary for describing the response to a wide range of rates, a linearin-variables function can be an adequate approximation over a sufficiently small subset of the function's domain. Since, as was the case in the experiments reported here, the ranges of rates tested in whole field OFPEs are restricted to reduce potential profit losses caused by experimental rates being far different from economically suboptimal rates, linear-in-variables models may be deemed best after the trade-off between goodness-of-fit and research expenses is considered. The polynomial function in Eq. 3 presents the full model considered to describe yield response for each field: In Eq. 3, Y i is the observed yield, S i is the seed rate, and N i is the nitrogen fertilizer application rate at a location i ∈ {1, …, I}. The βs in Eq. 3 can assume values according to one of the three scenarios: The standard GWR considers that all model parameters vary spatially, meaning that case (c) is assumed for every k. But in the research reported here, a mixed GWR model was applied, meaning that case (b) could hold for some k and case (c) for others-that is, that some terms could be local and others global. By limiting the number of coefficients that can vary over space, the mixed GWR approach can sidestep multicollinearity problems, which can lead to artificial spatial patterns of the coefficients (Geniaux and Martinetti 2018). This could be a problem, for example, when including coefficients of both the linear and the quadratic terms, for example, β 1(i) and β 4(i) in Eq. 3, as spatially varying. It has also been demonstrated that setting some parameters of a yield response function to a common value for all subareas within a field and year can improve both the goodness of fit and interpretability of a model (Bachmaier and Gandorfer 2012;Pahlmann et al. 2016). Based on these recommendations and to reduce the number of scenarios tested, model selection was performed setting all interaction and higher-order terms (i.e., terms 3, …, 8) as global variables.
The mixed GWR was also used to test whether a local parameter was needed for the intercept, S i , and N i terms in Eq. 3. Model selection was based on the AICc criterion . The change in the root mean squared error (RMSE) when a term was added as a local or a global variable was used to represent a term's importance. The significance of the coefficient at each location was also tested using the p-values from pseudot-tests, with the alpha adjusted by the Fotheringham-Byrne method (Lu et al. 2014).
The "MGWRSAR" R package (Geniaux and Martinetti 2018) was used for model calibration and the term inclusion. The "GWmodel" R package (Gollini et al. 2015) was used to provide additional diagnostics and to test parameter significance at every (3) location. Both packages were used because MGWRSAR computed the results much faster, which was especially important for testing multiple values of kernel type, bandwidth, and parameters in the model, while GWmodel provided more detailed output. Table 3 summarizes the equations used to compare optimal rates under the various scenarios. The assumptions described above about coefficients implied location-specific yield response functions f i (N, S), where units for yield, N, and S were ( Mg ha −1 ), ( kg ha −1 ), and ( kseed ha −1 ), where kseed denotes thousands of seeds. The status quo seed rate S sq and nitrogen fertilizer application rate N sq were those that the producers said they would have applied had they not taken part in the on-farm precision experiments.

Comparison of optimal rate scenarios
The following prices were used to calculate the economically optimal rates: corn price of p = US$ 160.00 Mg −1 (US$ 4.00 bu −1 ), corn seed price of w S = US$ 3.50 kseeds −1 , and nitrogen price of w N = US$ 0.88 kg −1 (US$ 0.40 lbs −1 ). These values imply approximate price ratios w S /p = 0.022 Mg kseeds −1 (0.88 bu kseeds −1 ), and w N /p = 0.0055 Mg kg −1 (0.10 bu lbs −1 ). The quantitative results reported below were calculated by replacing the expected yield response function f i (N, S) in Table 3 with the predicted expected yield response function in Eq. 3, and solving the maximization problems using standard first-order conditions from multivariate calculus, obtaining the values of the arguments S and N that maximize (argmax) the expected net revenue function (Bullock and Bullock 1994). In these scenarios it is assumed that all equipment needed is already available at no extra cost. For simplicity, costs associated with information management and data processing are also not considered in the economic analysis. The Table 3 Summary of the equations used in the comparison of optimal rates scenarios Equation Description The expected yield response function Status quo expected yield also referred to as reference yield The expected net revenue function Optimal nitrogen fertilizer and seeding plan under uniform management Expected yield under optimal uniform management r u * i = py u i − w N N u − w S S u Expected net revenue under optimal uniform management The optimal site-specific management plan Expected yield under optimal site-specific management

OFPE implementation
Variable-rate applicators generally apply inputs with some degree of error, rarely exactly matching to the target rate. Nevertheless, monitors can record very accurately the "asapplied" and "as-planted" rates, which are the actual rates at which the inputs "went into the ground". Therefore, the as-planted and as-applied maps presented highly accurate implementations of the trial designs (Fig. 3a, b). Historically, the main reason for using the targeted rates or grouping the as-applied rates into a few levels was to analyze variance using treatments as categorical variables. Since the proposed methodology treats treatments as continuous variables in all steps, such categorization was not necessary. The yield maps of the four fields reveal that most yield variation could be attributed to factors other than the treatment rates since there are no obvious patterns of treatmentyield relations (Fig. 3c). Field 1 presents a wide range of yield values, but not a clear spatial structure of the variability. Yields in Field 2 were the most spatially homogeneous among the four, with a yield range of only 2.0 Mg ha −1 . Spatial variability is more evident in Fields 3 and 4, with Field 3 on average producing the lowest yields and Field 4 the highest.
The factorial combinations of the experiment's several nitrogen and seed rates resulted in many (N, S) treatment levels in each field's trial design, making a visual interpretation of the results difficult. The quality of the yield data is currently a major concern in OFPE and also affect how the results look in the map. Errors in yield data can come from the different lengths of time needed for combines' internal threshing mechanisms to process the grains, which can reduce the contrast between the yields at different input rates (Lark et al. 1997;Lark and Wheeler 2003). To improve yield data quality, the convoluted yield monitoring process was taken into account in trial design, with plot lengths sufficiently long to allow the harvester to pass through any treatment plot for at least 30 s. Therefore, the effects of the yield convolution in the treatment transition zones were considered to be small, since most of the points included in each location were outside of the transition zones. To maintain methodological simplicity, the points in the transitions were kept in the data analyzed, and no explicit effect was included in the model. The importance of these effects should be explored in future research.
The dispersion of yield values for the same input rate is much higher in OFPE than what is commonly observed in small plot research (Fig. 4). The noise in the data comes from many sources, including the errors associated with the yield monitoring process Fig. 4 Scatterplots of the relationship between a yield and as-planted seed rates, and b yield and as-applied nitrogen rates, in the four on-farm precision experiments and the spatial variability of yield that would be observed even if inputs were applied at a uniform rate throughout the field. The only scenario in which an overall response can be observed was for the nitrogen rate in Field 1. The distribution of points also supports the decision of using linear regression to represent the yield responses.

GWR results
The optimal bandwidth converged to about the same number of neighbors independently of the kernel type and the field. This optimal bandwidth determined that between 200 and 300 parameters were necessary for the Gaussian and exponential kernels, and the number of parameters in the bi-square kernel was between 1000 and 1500. There are indications in the literature that the cross-validation tends to lead to an overfitted model with a too-small bandwidth, which can lead to instability and overestimation of the spatial variability in the fitted parameters (Murakami et al. 2018).
This seemed to be the case for the bi-square kernel. The Gaussian kernel presented lower AICc than the exponential in all fields, and it was therefore selected for subsequent steps.
In each field, the optimal size of the adaptive bandwidth parameter was found to allow the inclusion of about twenty neighbors. With the Gaussian kernel, these twenty neighbors captured about 67% of the total weights, and another eighty neighboring observations were needed to capture 99% of the total weights (Fig. 2). The convergence to the same bandwidth size even though the support area of the observation units varied from 37 to 87 m 2 , further demonstrates the cross-validation used to optimize the bandwidth is taking into account the variance-bias trade-off more than the spatial variability of the inputs Murakami et al. 2018). A discussion of alternative methods of parameter estimation is beyond the scope of the present article, but these nuances of the GWR method may significantly impact results and should be further explored.
A summary of the comparisons between the fitted GWR model and the alternative model using ordinary least squares (OLS) is presented in Table 4. The AICc and BIC values show that the GWR model fits the data significantly better than does the OLS model. The largest difference was observed in Field 3, where the RMSE decreased from 1.88 to 0.81 Mg ha −1 . The relative importance of each variable, as measured by the relative decrease in the RMSE, is one of the keys to understanding why there is such a significant difference between the two models (Fig. 5). The spatially varying intercept ( 0 ) explained 25% to 50% of observed yield variability; this percentage is related to the spatial structure of the field variability. For example, Field 3 had the most extensive range of the predicted reference yield and was also the field for which 0 explained most of the variability. Comparing the results in Table 4 and Fig. 5, it becomes evident that the explanatory power of the GWR model is also more related to the spatial structure of the yield variability than to the magnitude of the crop response to the treatments.
The global response caused a reduction in the RMSE greater than the local terms only in explaining the variability of the responses for seed rate ( 1g > 1l ) in Field 4 and nitrogen rate ( 2g > 2l ) in Field 1. For all other field and input combinations, the local terms were more prominent in explaining yield than the global terms. The spatial variability in response to S and N were consistent across all fields, with 1l explaining on average 3% and 2l explaining on average 5% of the total variability. The unexplained variations, represented by the error term, are caused mainly by the short distance variance caused by the noise in the yield monitoring system, and by other sources of variation such as crop damage by machinery and yield losses caused by pests and diseases.
The outcome of accounting for most of the yield variability by the variation in factors other than controllable inputs have been reported in other studies (Kindred et al. 2017). The limited range of variation in the treatment rates, in comparison to what is used in small plots, also plays a significant role in the explanatory power of the model. While in small plot trials, the nitrogen rates may vary from zero to twice the usual rates (0 to 200% of the status quo rate), in the OFPE reported here, the variation was of only 20% (80% to 120% of the status quo rate). The small range of variation in N and S rates, coupled with the broader importance of spatial effects, the large number of observations, and the spatial correlation between observations make the analysis and interpretation of OFPE results challenging. Initially, input ranges were chosen in an attempt to minimize field implementation costs and encourage farmer's participation. Trial design scenarios with more extensive input rate ranges are planned in future studies.
Although the 1 and 2 coefficients explained only about 10% of yield variation, these are the most important model parameters from a crop management perspective and they will be the focus of the discussion (Fig. 6). The values of 1 ranged from about −150 to 150 kg kseeds −1 . If coefficients β k = 0 for k = 3, 4, …, 8, then β 1 represents the change in the expected yield in kg ha −1 when the seed rate is increased by 1.0 kseed ha −1 (the marginal expected product of kseed), and 2 represents the marginal expected product of N. Since the grain equivalent cost of increasing the seed rate by 1.0 kseeds ha −1 is about 22 kg ha −1 (buying one kseed costs w S = $3.50, and selling 22 kg (i.e., 0.022 Mg) of grain provides 160 × 0.022 ≈ $3.50 in revenue). Thus, the "break-even" marginal product of kseed is = be 1 22 kg ha −1 . It is also possible to compare the distributions between fields, evidencing for example that Field 3 has a wider variation than Field 2. Similar results can be seen for the 2 values, with an overall range within -40 to 40 kg kg −1 . The grain-equivalent cost of increasing the nitrogen rate by 1.0 kg ha −1 is w N /p ≈ $5.50 kg ha −1 , which, if β 3 = … = β 8 = 0, is the break-even marginal product of nitrogen, called be 2 . The spatial variation of 1 and 2 can also be interpreted as proxies or latent variables that account for the combined effects that crop development, soil parameters, weather factors, and their interactions exert on the crop response to each input. The use of proxies of yield response is a common practice in most of the agronomic recommendations, since measuring all the factors influencing yields is either prohibitively expensive or virtually impossible. To be effective, the proxy variable has to be well correlated with the crop response, which may not be always the case for variables such as electrical conductivity or vegetation indices. The advantage of this type of OFPE and the GWR analysis is that the crop itself is used as a "sensor" to estimate the spatial variabilities of the crop response and the optimal rates. The use of soil and weather characteristics may still be necessary for modeling yield response under different growing conditions, such as in other fields and years. The main difference is that the focus is no longer yield prediction; instead, the target variable could be the yield response, optimal rates, or simply the 1 and 2 . The distribution of these parameters is the main difference between the yield response functions in each location, and are the main drivers of the differences in optimal input rates. Therefore the statistical significance of these variations is strong evidence of the within-field spatial variability of crop responses to agronomic inputs. The maps in Fig. 7 show the regions where that was the case. The significance is reflected in the raw p-values, and in the p-values adjusted by the Fotheringham-Byrne procedure (Lu et al. 2014). The percentage of locations with significantly different yield responses to (N, S) varies among fields, with Field 3 showing the most substantial variation. With the conservative protection given by the Fotheringham-Byrne procedure, less than 10% of the locations remain as significant in Fields 1 and 2. In general, no more than 50% of the locations are significantly different from the average. If temporal variability can be ignored, these maps could be used as management zones, clearly showing which parts would benefit the most from site-specific management. Considering temporal variability is possible, but beyond the scope of the research reported here.

Comparison of crop management scenarios
The predicted yields at the status quo rates showed smooth variations across the field (Fig. 8). These are the best estimations of the reference yield, and the variation is mainly due to differences in soil parameters which affect yield both directly, and indirectly through interactions among soil characteristics and weather. The estimated optimal seed rates ranged from less than 70 kseeds ha −1 to the maximum of the experimental seed rates, slightly more than 90 kseeds ha −1 (Fig. 9a). The number of fields evaluated was not sufficient to conclude how weather influenced crop response variability in each year. By examining the cumulative frequency distributions of optimal seed rates (Fig. 9a), it is possible to conclude that the optimal rate was lower than the status quo rate in over 85% of Field 2, but in less than 15% of Field 3. The range of variation in optimal nitrogen rates was about 100 kg ha −1 for Field 3 and only 40 kg ha −1 for Field 1 (Fig. 9b). However, the optimal nitrogen rate was higher than the status quo rate in the totality of Field 1, and the estimates of optimum rates were limited by the maximum rates used in the experiment, thus restricting the perceived variability.
For the four fields' geographic regions, the Maximum Return to Nitrogen MRTN program's recommended rates, which are the "best practices" promoted by several midwestern universities, ranged between 180 and 220 kg ha −1 (Sawyer et al. 2006). For Fields 2 and 4, the higher end of this range contained this study's estimated economically optimal N rates. However, the empirical results predict that Fields 1 and 3 could have benefitted from even higher rates. Note that Field 3 had the lowest yields (Table 2) and the highest optimal nitrogen rates among the four fields, providing additional evidence against yield-based N management recommendation algorithms.
The patterns of the spatial variability for the optimal rates of seed and nitrogen (Fig. 10a, b) are consistent with the significance of the local coefficients presented in  Fig. 7. In some regions, there is an indication of a weak negative correlation between both inputs, with low optimal seed rates being more frequently associated with higher optimal nitrogen rates. There were also many locations where the optimal rates were either high or low for both inputs. It is possible that in the regions with a negative correlation, the main driver for the optimal nitrogen rates was the soil availability, while in the region with positive correlation, the main driver was the plant demand . For reasons outlined above, understanding the reasons behind the crop response is important, and a key limitation of GWR analysis is that it cannot provide such understanding. To draw robust conclusions about the causes of yield response, more fields and weather scenarios should be investigated, and formal statistical methods need to be applied. This is left for future research.
The resulting range of the spatial variability of optimal input rates was shorter than the range of the reference yield variability. This again raises questions about the efficacy of using yield-based management zones to establish prescriptions for variable rate input rate management; for the four fields examined here, yield-based algorithms would recommend similar management strategies for areas having similar yields but different yield responses to managed inputs (Rodriguez et al. 2019). The results also provide Fig. 10 Spatial distribution of predicted optimal rates for a seed (kseeds ha −1 ), and b nitrogen (kg ha −1 ) and predicted corn yield (Mg ha −1 ) if the optimized rates of seed and nitrogen had been used in the fields understanding about why other strategies, such as the nitrogen-rich strip and the ramped calibration methods used to calibrate sensor-based nitrogen applications, often fail to provide reliable nitrogen requirement predictions (Colaço and Bramley 2019;Roberts et al. 2011). Since they are allocated with no previous knowledge of crop response, and the range of the spatial variability is limited, recommendations based on the nitrogenrich strip will be highly variable, depending on the field to which it was applied. Figure 11 illustrates the impact of using site-specifically optimized input rates. The difference between the reference yield and the optimized yield reflects the potential for improving yields by applying site-specifically according to the optimal rate of each portion of the field. There are no distinguishable differences in yield for the different scenarios in Field 2, indicating the status quo rates were already close to the optimal uniform rates for the whole field. The yield improvements were also observed in the high yielding areas of Field 1. The more compact the distributions are, the smaller is the yield variability. Comparing the cumulative frequency distribution of the reference yield and the optimized yield for Field 1, it is also possible to note that the variability of the optimized yields is higher than the variability in the reference yield. This shows that the objective of variable rate application is not to make the yield more homogenous, but rather to accept the spatial variability and take advantage of it by maximizing profits at every location within the field.
The spatial variability of the optimized yields ( Fig. 10c) closely follows the variations in the reference yield (Fig. 8). This is related to the observation that most of the yield variability cannot be explained by controllable inputs (Kindred et al. 2017). In some areas, such as in Field 4, it is possible to observe some "hot-spots" with high yields, associated with high optimal rates of the inputs (Fig. 10). These are the locations where the local estimates of crop response had the higher values and could be due to a combination of factors that make those areas highly responsive, or due to random chance and the noise in the measurements. Statistical methods such as the bootstrap may be used to estimate the confidence intervals of these estimations and refine the local predictions by excluding some observations from the analysis (Harris et al. 2017).

Economic results
The estimation of any opportunity index of the potential economic benefits of spatially adjusting crop inputs requires the spatial resolution used to characterize the variability Fig. 11 Cumulative frequency distribution of yields comparing the observed yield, the expected yield if the status quo rates had been applied and, the expected yield if the optimized rates of seed and nitrogen had been applied in each field to be compatible with the resolution at which the recommendations can be applied (Leroux and Tisseyre 2018). The economic opportunity presented here implicitly considers the machinery sizes as given in the analysis because the spatial resolution of the data was defined by the ability of the available machinery to apply the assigned treatments. Table 5 reports that, except for Field 4, adjusting the uniform rates of both inputs would have an economic impact substantially greater than the loss of revenue due to the suboptimal trial rates. Considering the eight possible combinations of fields and inputs, in five of them, profits from the status quo rate were within US$ 11.00 from the maximized profit. The spatial homogeneity of Fields 1 and 2 resulted in small differences in the local yield response to the agronomic inputs. For Fields 3 and 4, taking into account the spatial variability of optimal rates would have generated a positive economic impact greater than adjusting the optimal uniform rate for both seed and nitrogen rates. The average increase in net revenue of using the site-specific optimal seed and nitrogen rates would be US$ 17.00 ha −1 and US$ 36.00 ha −1 , respectively. These results are similar to the overall average profit of US$ 30.00 ha −1 reported with the use of crop sensors to drive variable rate nitrogen application (Colaço and Bramley 2018).
Although the number of fields used in this work was adequate for demonstrating the method and making initial considerations, the extrapolation of these results will require more data. Nevertheless, the tendency of OFPE resulting in higher economic benefits on fields with higher yield variability, even though yield levels and optimal input rates are not necessarily correlated, is a promising result. This tendency means that fields with higher yield variability should be prioritized for the implementation of whole field OFPE, while in fields with less variability, the experiments may be replicated only in a small part of the field.
The main question that remains is the extent to which the results from 1 year of data can be used to guide the decisions of the next year, especially considering that some input strategies, such as seed rate, cannot be changed in-season. Answering that question requires knowledge of the temporal stability in the spatial variability of crop responses to each input, which is not extensively explored in the scientific literature. The majority of the related studies focus only on the temporal stability of the spatial variability in yields, which interacted with the weather of each growing season (Maestrini and Basso 2018). Many studies and reviews have focused on within-and among-field variations in optimal nitrogen rates, and have concluded that most of the variability is related to Table 5 Summary of net revenue losses due to yield loss in suboptimal treatment rates, and potential economic benefits if an optimized uniform or variable rate had been used instead of the farmer's status quo rates *All values are in US$ ha −1 at a corn price of US$ 160.00 Mg −1 (US$ 4.00 bu −1 ) SR seed rate,NR nitrogen rate, SRNR seed and nitrogen rate, VR-UR difference in profit of applying the optimal variable rate compared to a uniform optimal rate  dynamic variables, such as precipitation and temperature (Puntel et al. 2019); and their interaction with field-specific variables, such as previous crop, tillage practice, soil drainage class, and N form and timing Tao et al. 2018).
In terms of the model described in Eq. 3, the results available in the literature suggest that at least the global s will be season-specific, and therefore may need to be adjusted for every season. Depending on the degree of temporal variability in 1 and 2 , there could be no need to repeat the whole field experiment every year. Instead, a small subset of plots in representative areas of the field could be used to adjust the other parameters of the equation. Even more promising, the season-specific parameters may be adjusted by the use of other tools such as crop modeling systems and pre-plant soil nitrate tests. Optimizing the trial designs to account for the temporal variability should also be explored in future research. The relative importance of temporal and spatial variability in the overall crop response variability is likely to be different for each field, crop and, input considered. The methods presented here have the potential to be applied in many other scenarios to improve management decisions by accounting for these sources of variability.

Conclusions
The main contribution of this work is to demonstrate an alternative method to test and characterize the spatial variability of crop response to inputs. The combination of OFPE and GWR proved to be an effective methodology to test precision agriculture central hypothesis' of whether there is significant within-field variability in optimal rates. It also allows changing the focus from yield-based to response-based variable-rate prescriptions for crop input application. Future research on trial design and models with spatially varying coefficients for OFPE is advised.
Incorporating spatial heterogeneity of yield responses into model parameters improved model performance in all four fields evaluated. On average, the RMSE of the fitted yield decreased from 1.2 Mg ha −1 in the non-spatial model to 0.7 Mg ha −1 in the GWR model, and the r-squared increased from 10 to 68%. In 10% to 50% of the observations, the coefficients of the local parameters were found to be significantly different from the average, providing further evidence of the need for increased knowledge about local yield response functions.
In Fields 1 and 2 the greatest benefits of OFPE would come from optimizing the field's uniform rate, while in Fields 3 and 4 the highest revenue increases would account for the spatial variability in crop response, and implementing the site-specifically optimal variable rates would result in the highest increase in revenue. The average potential gain of using optimized uniform rates of seed and nitrogen was US$ 65.00 ha −1 , while the added potential gain of using variable rate application was US$ 58.00 ha −1 . permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.