A modified delta yield approach for estimation of economic optimal nitrogen rate ( EONR ) for wheat ( Triticum aestivum L . ) and barley ( Hordeum vulgare L . )

A Swedish field trial database was mined for information on economic optimal nitrogen rate (EONR). A total of 100 trials with wheat (Triticum aestivum L.) and 47 trials with barley (Hordeum vulgare L.) were used to parameterise prediction models for EONR based on yields in plots with no nitrogen (N) fertilisation, intended to reflect N mineralisation, and yields in plots with a high N rate, aimed as a proxy variable for yield potential, i.e. a modification of the delta yield (∆Y) approach. When the prediction models were applied to new sites and years, predictions had mean absolute error (MAE) = 10 kg N ha-1 for wheat and = 9 kg N ha-1 for barley. Performing modified ∆Y experiments can complement current N fertilisation trials with many rates, in order to improve the spatial representation of EONR estimations. Moreover, ∆Y experiments can potentially be used for in-season EONR estimation, in which case the accuracy of the EONR predictions depends also on the uncertainty in yield predictions made at the time of supplementary fertilisation.


Introduction
The EONR concept Crops respond to nitrogen (N) application by increasing yield up to a limit at which other factors limit growth.The economic optimal N rate (EONR) is defined as the fertilisation rate where the slope of the response curve equals the price ratio of the fertiliser N and the grain produced.This EONR definition is based on a simple economic model that is commonly used as the basis for N rate recommendations (for further discussion on the economics of the EONR concept, see Kyveryga et al. 2007).Different types of univariate functions can be fitted to data from field trials.Rajsic and Weersink (2008) reviewed use of ex ante N rates (i.e.in-season N recommendations) and ex post N rates (post-harvest determination of maximum economic rate of nitrogen) in Canadian maize (Zea mays L.) production and found that the ex ante N rates were generally higher than the ex post N rates.They compared different response model types, but found no evidence that the discrepancy between ex ante and ex post N rates was due to the type of model used to describe the N response of maize.However, earlier studies report a systematic difference in EONR depending on model type (e.g.Bullock and Bullock 1994).Whether it is a good choice or not, second-grade polynomial is the model type used by the Swedish Board of Agriculture for deriving N recommendations for cereals (Albertsson et al. 2016).

Environmental relevance of EONR
The EONR has been demonstrated to be a threshold for increased environmental risks.Lord and Mitchell (1998) found that fertilisation rates above EONR increased N leaching, and Delin and Stenberg (2014) confirmed this in Swedish experiments.Both studies found that the leaching is not linearly related to N fertilisation rate.Instead, nitrate leaching shows a weak trend in relation to N rate for N rates <EONR and a strong accelerating trend in relation to N rate for N rates >EONR.Thus, it is crucial from both an economic and an environmental perspective not to exceed EONR.
The ∆Y concept for EONR prediction Kachanoski et al. (1996) demonstrated that EONR is poorly correlated with the yield level in maize (cf.Raun et al. 2011) and consequently they developed the ∆Y concept, where ∆Y is defined as the difference between yield in plots receiving the EONR (Yeonr) and yield in plots with no N application (Y0), see Equation 1.
The ∆Y has been demonstrated to be a much better indicator of EONR than yield goal (Lory andScharf 2003, Fixen 2006) or Yeonr alone (Nel and Bloem 2006).What is interesting to note is that relationships parameterised with field trial data from the USA and relationships parameterised with South African field trial data are remarkably similar (see comparison by Nel and Bloem 2006).The theory behind the ∆Y concept for EONR estimation is described in detail in Kachanoski (2009).

Yield predictions for in-season EONR predictions
The ∆Y concept is directly applicable to post-harvest EONR predictions, but a significant drawback of the approach is that the yield in experimental plots is not known at the time of the final N application, so the EONR prediction models cannot be applied for in-season predictions of current-year EONR without including a yield prediction step.Using yields predicted at the time of fertilisation instead of measured yields inevitably means less accurate EONR predictions, but this is still an interesting approach.There are promising methods for in-season prediction of wheat (Triticum aestivum L.) yield, based on proximal crop canopy (multi-or hyperspectral) reflectance sensors (e.g.Solie et al. 2012, Øvergaard et al. 2013, Engström and Piikki 2016, Montesinos-López et al. 2017), proximal crop canopy fluorescence sensors (Zecha et al. 2017) and satellite-borne hyperspectral imaging (Wang et al. 2014).Previous work has indicated that predictions made at later growth stages are generally more accurate than earlier predictions and that hyperspectral measurements are generally more accurate than predictions based on vegetation indices calculated from a few wavelength bands.

Aim
The aim of the present study was to mine a field trial database comprising cereal yield data for a period of 13 years at different locations in Sweden.Specific objectives were: 1. To examine the magnitude of the variation in EONR, soil N supply (indicated by the yield in plots with no N application) and yield potential (indicated by the yield in plots with a high N rate) for wheat and barley (Hordeum vulgare L.).
2. To evaluate the performance of a modified ∆Y approach when applying parameterised models for new sites and years (i.e.sites and years not included in the parameterisation dataset).
3. To simulate the accuracy of the EONR predictions in relation to increasing input data uncertainty.The results from this exercise were intended to hint at prediction accuracy when using sensor-based yield estimates at the time of fertilisation, instead of observed yield values at harvest. 4. To develop and visualise prediction models for EONR in wheat and barley, using a modified ∆Y approach.

Field trial data
A subset of data from the Swedish field trial database (Field Research Unit, Swedish University of Agricultural Sciences) consisting of 100 winter wheat trials and 47 barley trials was used.The wheat trials were conducted in 1999-2011 and the barley trials in 2002-2011.The selection criteria for barley trials were: highest N rate at least 160 kg N ha -1 , EONR at most 160 kg N ha -1 , at least four N levels and that the trial site was georeferenced.The selection criteria for wheat trials were: highest N rate at least 180 kg N ha -1 , EONR at most 300 kg N ha -1 , at least four N levels, cereal pre-crop and that the trial site was georeferenced.No restrictions were imposed on cultivar or N strategy (number of N fertilisation occasions; only the total amount of N was considered).For sources of agronomic details on the individual experiments, see Appendix A.

EONR calculation
A second-grade polynomial (Equation 2) was fitted for each trial and EONR was determined as the N rate where the derivative of the polynomial equalled the price ratio of the fertiliser N and the grain produced.A price ratio of 10 was used as a case example, but the derived model can easily be parameterised for other price ratios.
Grain yield = a + b × N rate + c × (N rate) 2 Equation 2Modification of the ∆Y concept In order to apply a parameterised relationship between EONR and ∆Y to new sites and years, it would be tempting to carry out simple two-rate N experiments, one N rate being 0 and the other being EONR.This is obviously not possible, since the EONR is not known at the time when the trial plots are being fertilised.In the present study, we therefore modified the ∆Y concept and used yield at one high N rate instead of Yeonr (Equation 3).The chosen rates were 200 kg N ha -1 for barley and 300 kg N ha -1 for wheat.The yields at the predetermined high N rates were denoted Yhigh or Y200 and Y300 for barley and wheat, respectively.
∆ Y = Yhigh -Y0 Equation 3When an experiment with two N rates has been carried out, it is a waste of information not to use the two yields as separate predictor variables, instead of merging them into one (the ∆Y).Their effect on EONR may not be of the same magnitude and that information would be lost if the relationship were to be based on ∆Y.We therefore parameterised a bivariate linear function for EONR prediction (Equation 4).A model with more parameters is more flexible and has the potential to better parameterise the variation in the response variable.On the other hand, using more parameters carries an increased risk of overfitting, which would mean less robust models (Hastie et al. 2009).
For comparison, univariate linear prediction models were also parameterised, one for each predictor (Equations 5-6): EONR = a + b × Y0 Equation 5EONR = a + b × Yhigh Equation 6k-fold cross-evaluation For a prediction model to be of any value, it has to make reliable predictions when applied to new sites and years not present in the parameterisation dataset.Therefore a k-fold cross-evaluation strategy was applied.A fold (subset of data) consisted of all trials performed in one year.One such fold was withheld from model parameterisation and the resulting model was then applied to the withheld fold.Additionally, any trials conducted in the neighbourhood (<1 km) of a trial in the withheld fold were omitted from parameterisation.The procedure was repeated k times, until EONR was predicted for all trials.This evaluation design was applied because it has been demonstrated that the exclusion of nearby sites is crucial for the identification of robust models (Piikki et al. 2016).

Evaluation measures
The mean absolute error (MAE; Equation 7) and the modelling efficiency (ME; Equation 8; Nash and Sutcliffe 1970) were used to quantify the modelling accuracy in absolute (MAE) and relative (ME) terms.The MAE indicates the magnitude of errors, but does not reveal whether there is any general over-or under-estimation of EONR.The ME is defined in the half-open interval (-∞, 1).ME = 0 indicates that using the model is not better than using the mean of the parameterisation dataset, while ME = 1 means that the model predicted all EONR values correctly.The letter o in Equations 7 and 8 denotes observed values (i.e.EONR determined in field trials) and p denotes predicted EONR.
Simulation of accuracy loss caused by yield data uncertainty Random errors (noise) were added to the predictors for simulation of model accuracy in relation to yield data uncertainty.The errors were randomly sampled from Gaussian distributions with a mean of zero and a standard deviation of 400, 800, 1200, 1600, 2000 and 2400 kg ha -1 .
MAE values were determined by the k-fold cross-evaluation strategy described above.In this case, however, the models were parameterised with non-noisy yield data and predictions were made using the yield data with added noise as input.The sampling of errors and the k-fold cross-evaluation were repeated 100 times for each error level.The MAE values averaged for the 100 simulation runs were then plotted against the standard deviations of the added errors.For reference, the range of reported errors in earlier studies was added to the plot; e.g.Øvergaard et al. (2013) reported a root mean squared error of prediction (RMSEP) of 1560 kg ha -1 for wheat yield predictions based on hyperspectral measurements in field trials in Norway, while Engström and Piikki (2016) reported an MAE of 1100-1600 kg ha -1 for predictions based on a handheld hyperspectral sensor in Sweden.The spectral measurements by Øvergaard et al. (2013) were conducted during anthesis (growth stage [GS] 50-60; Zadoks et al. 1974).
That is somewhat later than normal for supplementary N fertilisation of wheat in Sweden.The measurements by Engström and Piikki (2016) were carried out between GS 37 and GS 63.In Swedish wheat production, the last N application is most often carried out between GS 37 and GS 45, but sometimes as late as GS 65 (mid-anthesis).
In barley, the second fertilisation is often carried out between GS 32 and GS 37 and, if the crop requirement is judged to be large, an additional third application is made during GS 37-49.

Parameterisation of final models
Final models were parameterised by all trial data.These models can be expected to produce EONR predictions of the accuracy indicated by the cross-evaluation when applied to new sites (in Sweden) and years.

Software
All data analyses were carried out using the statistical software R (R Core Team 2017), the geographical information system (GIS) software ArcMap 10.4 (Esri Inc., Redlands, CA, USA) and the spreadsheet software Excel 2013 (Microsoft, Redmond, WA, USA).

Quantifying variation
The total range of variation in EONR in the dataset (both temporal and spatial variation) was 245 kg N ha -1 for wheat and 112 kg N ha -1 for barley (Table 1).The variation in EONR between locations in the same year and the variation in mean EONR between years can be seen in Figure 1.For wheat, the magnitude of the spatial variation seemed somewhat larger than the temporal variation, although in some years (1999 and 2003) the average EONR was clearly different from in other years.A similar trend was not found for barley, where EONR also differed considerably between years.
There was also considerable variation in the yield with no N fertilisation (Y0) and in the yield with no assumed N limitation (Y300 and Y200).The total variation in yield with no fertilisation was 4623 kg ha -1 (wheat) and 4248 kg ha -1 (barley), disregarding one extremely low yield observation of 401 kg ha -1 , while with high fertilisation rates it was 9840 kg ha -1 (wheat) and 8945 kg ha -1 (barley).
Table 1.Descriptive statistics on the economic optimal nitrogen rate (EONR), the yield at 0 kg N ha -1 (Y0) and the yield at 300 kg N ha -1 (Y300) (wheat) and at 200 kg N ha -1 (Y200) (barley Finding the best predictor(s) For wheat, models based on Y0 and Y300 resulted in EONR predictions with the smallest MAE (10 kg N ha -1 ) and the largest ME (0.80) (Fig. 2).Using both predictors was a considerable improvement over using only one of the predictor variables.For barley too, using both predictors (Y0 and Y200) was superior to using either of the predictor variables alone (Fig. 2).Slight under-prediction of high values and over-prediction of low values was observed.This is a very common phenomenon in prediction modelling and happens when the model explains most, but not all, of the variation in the parameterisation dataset (some generalisation is made).Extremes are then predicted as less extreme.A model that is parameterised to capture almost all variation in the training dataset (including the noise) is likely to be very complex and less robust (less accurate when applied to new datasets).

Simulated effects of input data uncertainty
The simulated effects of uncertainty in input data (intended to mimic the uncertainty in yield predictions from e.g.proximal sensing) are presented in Figure 3.The prediction error increased rapidly with increased uncertainty in input data and the models with two predictors were more sensitive than those with one predictor.As can be inferred from Figure 3a, with the uncertainty of yield predictions reported in the literature, the MAE of EONR in wheat is 25-35 kg N ha -1 when the best model (Equation 3) is used.

Models for EONR prediction
The final prediction models (parameterised using all trial data, no trials withheld) for EONR prediction are visualised, and their equations are given, in Figure 4.As expected, the EONR decreased with increasing Y0 and increased with increasing Yhigh for both crops.The same colour scale is used for each row of plots and it is evident that the models based solely on Y0 do not accurately capture the variation in EONR.The Yhigh predictors seem to explain somewhat more of the site × year EONR variation than the Y0 predictors, but Y0 and Yhigh together were superior as a basis for EONR prediction.Fig. 4. Univariate and bivariate linear regression models for prediction of economic optimal nitrogen rate (EONR).The models were parameterised by all trial data (no trials left out).The filled circles show the EONR of the individual trials and the colour scale shows the prediction surface.The circle size and the colour scale of the surface are the same within each row of plots but different between the rows (i.e.comparable between predictor sets of the same crop, but not comparable between crops).Y0 = yield at 0 kg N ha -1 , Y200 = yield at 200 kg N ha -1 and Y300 = yield at 300 kg N ha -1 SD of added noise (kg ha -1 ) SD of added noise (kg ha -1 )

Discussion
Using both predictors (Y0 and Yhigh) proved superior to using only one, but the models based on two yield variables were more sensitive to uncertainty in input data, such as can be expected when using in-season yield predictions instead of post-harvest yield determinations.For wheat, an error level of 25-35 kg N ha -1 was estimated for the simulated case of bivariate EONR modelling using sensor-based yield predictions made just before the last N application as input data.Whether this is an acceptable error level or not depends on which alternative methods and tools there are to decide on an appropriate N rate.
Another source of prediction uncertainty (relevant in any decision support for supplementary N rates that takes economics into account) is that it might be difficult to know the price of the produce at the time of fertilisation, and thus the price ratio on which the EONR depends.
The spatial variation in EONR is large, from field scale to within large regions.The soil N supply by mineralisation depends on soil type (Delin and Lindén 2002), weather conditions and the previous crop, and the variation within fields can be >100 kg N ha -1 (Wetterlind et al. 2008).Moreover, crop yield has been demonstrated to vary considerably, by several tonnes per hectare, within fields (Algerbo et al. 2003, Blackmore 2003).It is evident from the data presented here (Fig. 1) and data presented by Stenberg et al. (2009) that the spatial variation at national scale is also considerable.In Sweden, a number of N fertilisation field trials are conducted across the country each year to support the Swedish Board of Agriculture in their work of deriving recommendations on N rates to cereals and other crops.These experiments usually have a randomised block design (with 4-7 N levels and four blocks).This means that EONR is determined with high precision at the relatively few sites ( ~10 per year) where the trials are conducted, but these sites are not necessarily a good representation of a region for which recommendations are derived.A simpler form of trial for EONR prediction (such as the modified ∆Y trials examined in the present study) could be a valuable complement.What is lost in precision (how well EONR is predicted) in the simpler trials should be more than gained in better coverage of the spatial variation.If combined with yield predictions, it would be possible also to estimate EONR at the time of supplementary fertilisation (at the cost of some loss in prediction accuracy, as discussed above).There are currently commercial companies in Sweden that conduct measurements with handheld hyperspectral crop canopy sensors in Y0 plots in farmers' fields.Such a service could gain from adding high-N plots and using the modified ∆Y models for EONR prediction.Nitrogen experiments are thereby another tool that a farmer can use to optimise N fertilisation to small-grain crops.As none of the tools available is perfect, the best option is to combine them, for example by splitting the N application between three occasions for the best possibility to adjust it to the needs of the current crop, using absolute N rate decided from two-rate experiments in which yield is estimated with a proximal sensor and using satellite-based vegetation index maps (Söderström et al. 2017) to vary the N rate spatially within fields.
The k-fold evaluation strategy, where the parameterisation and the evaluation data are always both temporally and spatially independent of each other, allowed us to draw conclusions on expected accuracy for untested sites and during coming years.This does not mean that the models developed here should be regarded as valid for all times; rather, they can be regarded as 'living models' that should be updated yearly using newly performed manyrate N fertilisation trials, in order to keep them up-to-date regarding cultivars and agricultural practices.
We found that linear functions explained the variation in EONR well (the bivariate case), whereas Kachanoski et al. (1996) and Nel and Bloem (2006) used non-linear ∆Y models.One reason for this difference may be that we used yield at a predetermined high N rate as the high N rate predictor, while the previous studies used yield at EONR.

Conclusions
We developed a modified ∆Y modelling approach to predict EONR for wheat and barley based on two-rate N trials.The Y0 treatment (0 kg N ha -1 ) is intended to give an indication of the soil N supply and the Yhigh treatment (200 kg N a -1 for barley and 300 kg N ha -1 for wheat) is used as a proxy for the yield potential.The models were parameterised with data from a total of 147 N fertilisation trials conducted in wheat and barley at different locations in Sweden over a period of 13 years.The model evaluation was designed to inform about prediction accuracy when the method is applied for new sites and years.
We observed considerable variation in EONR between sites and years; the EONR for wheat differed by up to 245 kg N ha -1 and that for barley by up to 112 kg ha -1 between individual experiments.The great variation means that recommendations for N fertilisation based on average values can be rather erroneous for a certain site (farm/field/ management zone/point location) and year.It also means that trials at only a few sites to support national advisory services would be very sensitive to where the trials are located.The method presented here -of performing simple field trials with only two N levels -could be a valuable complement to the many-rate N fertilisation trials currently conducted, in order to better cover the spatial variation across the country in crop N response and to provide better advice for farmers based on local conditions.
Provided that accurate input data are available, EONR could be predicted with MAE of 10 kg N ha -1 for wheat and 9 kg N ha -1 for barley, using a bivariate linear model.The prediction error increases with increasing uncertainty level in the input yield data.The models presented here can be used to extend the current plan for field trials to achieve better geographical coverage.EONR predictions can also be made at the time of fertilisation, but depend on estimation of yield (e.g. by use of proximal optical crop canopy sensors).
A drawback of the method, as with all N fertilisation trials, is that EONR determination is based on grain yield and the yield is not known until harvest.In order to minimise the environmental risk of nitrate leaching and to improve the profitability in crop production, it is necessary to adjust the N fertilisation rate to both local and current conditions.This requires the development of accurate yield prediction methods.

Fig. 1 .
Fig. 1.Spatial and temporal variation in yield in (a-b) trial plots without fertilisation (Y0) and (c-d) trial plots receiving an N rate of 300 or 200 kg N ha -1 ).Plots (e-f) show the economic optimum nitrogen rate (EONR).Each circle represents one field trial and the solid lines represent the time series of yearly averages.

Fig. 2 .
Fig. 2. Predicted values of economic optimal nitrogen rate (EONR) plotted against EONR determined in field trials for wheat (a-c) and barley (d-f).The predictions are from the k-fold cross-evaluation of models.The models were linear models with different predictor sets (given in brackets).The dashed line is the 1:1 line and the solid line shows a linear regression model between the predicted EONR and the EONR determined from experimental data.Y0 = yield at 0 kg N ha -1 , Y200 = yield at 200 kg N ha -1 and Y300 = yield at 300 kg N ha -1 .MAE = mean absolute error; ME = Nash-Sutcliffe modelling efficiency

Fig. 3 .
Fig. 3. Mean absolute error (MAE) of predictions of economic optimal nitrogen rate (EONR) in relation to the magnitude of random errors added to the predictor variables.The grey area in diagram (a) indicates the range of reported error levels for yield prediction based on proximal sensing (Øvergaard et al. 2013, Engström and Piikki 2016).Y0 = yield at 0 kg N ha -1 , Y200 = yield at 200 kg N ha -1 , and Y300 = yield at 300 kg N ha -1