Introduction

In the context of contemporary agriculture, enhancing the efficiency of nitrogen fertilization stands as a crucial objective. Incorrect fertilizer decisions can be costly in terms of farm economics and external costs due to environmental pollution (Meyer-Aurich et al., 2010b; Sutton et al., 2011). Thus, a careful application of nitrogen is necessary for achieving the best possible crop yield while simultaneously mitigating environmental concerns and effectively utilizing resources (Paul & Ford, 2002). However, limited understanding of the input–output relationship results in inaccuracies when fitting crop yield response functions. This inaccuracy can introduce biases in determining economically optimal fertilizer application rates, as well as optimal rates from agronomic, ecological, or social perspectives (Dhakal & Lange, 2021).

The complexity of this endeavor is magnified when exploring within-field site-specific nitrogen fertilization (Anselin et al., 2004; Meyer-Aurich et al., 2010a). Unlike uniform nitrogen application, site-specific nitrogen fertilization acknowledges the inherent variability within fields, taking into account factors such as soil properties, topography, and crop variability within a single agricultural area (Maestrini & Basso, 2018). Determining the optimum nitrogen rate under such circumstances requires adaptable and nuanced approaches that can capture and accommodate the spatial variations present across the field (Basso et al., 2016).

The selection of a suitable mathematical model that accurately characterizes the intricate interplay between nitrogen application and crop yield is fundamental to making informed decisions about nitrogen management. However, the process of selecting the optimal functional form is not a one-size-fits-all endeavor; it heavily relies on the specific characteristics of the crop, soil, and environmental conditions prevalent at the particular site (Dhakal & Lange, 2021). This introduces uncertainties with respect to model selection that influences, for example, net revenues and nitrogen balances (Henke et al., 2007). The choice of an inappropriate model carries implications for both the chances of correctly identifying optimal fertilizer rate and the marginal costs associated with a reduction of fertilizer rates. Unlike continuous models, where the marginal opportunity cost of nitrogen fertilizer is zero at the economically optimal nitrogen rate, a linear plateau function exhibits a constant slope or yield response within its linear range. Consequently, optimal fertilizer rates exhibit substantial variations across different functions. Meyer-Aurich and Karatay (2022), for instance, identified a range of optimal fertilizer rates of up to 86 kg N ha−1, highlighting the sensitivity to the choice of functions.

Studies have been conducted in order to address the uncertainty with regards to the estimation of economic optimum nitrogen rates (EONR). Jaynes (2011) compared different models using a Monte Carlo approach to compute confidence intervals and found considerable statistical uncertainty in the true EONR based on any of the yield response functions and even greater uncertainty when comparison are made across functions. Hernandez and Mulla (2008) proposed the use of profile-likelihood and bootstrap methodologies to calculate confidence intervals as standard procedures of describing uncertainty with regards to EONR estimation. Bachmaier and Gandorfer (2012) estimated likelihood confidence intervals of EONR for the quadratic and linear-plateau function in the context of site-specific nitrogen aplication to capture the EONR uncertainity. They found extremely wide confidence intervals with little pratical use. Nigon et al. (2019) computed the EONR and its confidence intervals for the quadratic-plateau function with applications that consider variable costs and/or externalities. Even though different authors argue for different functional forms, there is often no statistical proof for the one or other function (Cerrato & Blackmer, 1990; Henke et al., 2007). Nor does agronomic reasoning provide robust evidence about the true functional form. In a recent study, Miguez and Poffenbarger (2022) explored the use of model averaging with information criteria weights, advocating for its potential to yield less biased estimates of the true agronomic optimum nitrogen rate (AONR) when compared to the traditional practice of selecting a single model. This approach considers the uncertainty associated with model selection, eliminating the need for a definitive choice of the most appropriate model, a decision that often rests on the farmer or modeler. While the concept of model averaging holds inherent appeal, the impact of this approach on the confidence intervals of Economic Optimum Nitrogen Rate (EONR) has not yet been thoroughly investigated, nor has there been a determination of the optimal number and specifications of models to include.

This study contributes to the ongoing discourse on nitrogen yield response modeling, with a specific focus on EONR. We consider commonly employed models such as linear-plateau, quadratic-plateau, Mitscherlich, and quadratic functions, as highlighted in previous literature (Lyons et al., 2019). Our investigation extends beyond the choice of individual functional forms to explore the performance of model averaging against the selected models. Performance measures, including bias, mean squared errors, and variance, are employed for a comprehensive evaluation.

Moreover, our study integrates experimental data to undertake a thorough exploration of the model averaging approach, particularly in the context of site-specific nitrogen fertilization. We aim to shed light on how model averaging impacts confidence intervals for optimal nitrogen application rates, providing valuable insights for practical applications. By addressing these aspects, our research adds depth to the understanding of nitrogen yield response modeling and the implications of employing a model averaging strategy.

Methodology

Data and experimental design

The data (Fig. 1) was collected in the Sieblerfeld test field (5 ha), located in Upper Bavaria, Germany. This field exhibits two distinct yield zones. The high-yield zone features sandy loam soil with a available field capacity of 160 mm in the rooted soil zone, while the low-yield zone has loamy sand soil with a available field capacity of 100 mm in the rooted soil zone (0–90 cm) (Liebler, 2003). The trial employed a randomized complete block design with four blocks in each yield zone. To explore the relationship between yield and fertilizer rate, 11 plots from each block in both zones were randomly chosen. These plots received varying N rates (0, 80, 100, 120, 140, 160, 180, 200, 220, 240, and 260 kg/ha) and yield measurements were conducted using a plot combine. Plots were sized at 12.5 m2 (10 m length, 1.25 m width). Fertilizer rates were split, if single N rates were higher than 60 kg N/ha. Above 140 kg N/ha the fertilizer rates were split in four application rates to avoid nutrient losses. The source of N was calcium ammonium nitrate (27%N). The experiment was established in 2002 with the crop winter wheat (Triticum aestivum L.). More details on the experiment can also be found in Bachmaier and Gandorfer (2012) and Bachmaier and Gandorfer (2009). Table 1 presents the analysis of variance (ANOVA) results of the experimental data.

Fig. 1
figure 1

Yield response data from high and low yield zones

Full size image
Table 1 Analysis of variance (ANOVA) results of the experimental data
Full size table

Block effect elimination

To ensure the robustness and reliability of our results, it was essential to account for the block effect (Table 1). To address this, we calculated the overall yield mean, representing the average yield across all experimental units in the study. Additionally, the average yield for each of the four blocks separately was computed. To eliminate the block effect from the original yield data, we subtracted the overall yield mean (\(OYM\)) from the block-specific yield means (\(BM\)). This subtraction yielded the block effect for each block, representing the deviation of a block’s yield from the overall mean. By subtracting this effect from the original yield data, we effectively removed the block-specific variability while retaining the treatment-specific information. Mathematically, this block effect (\(BE\)) elimination can be expressed as follows for each block \(k\):

$$BE_{k} ~ = ~BM_{k} - OYM$$
(1)

where \(BE_k\) represents the effect for block \(k\), \(BM_k\) is the average yield for block \(k\), and \(OYM\) is the overall yield mean. The corrected yield (CY) for each experimental unit in block was calculated as follows:

$${\text{CY}}_{ik} = {\text{OY}}_{ik} - {\text{}}BE_{k}$$
(2)

where \({\text{O}\text{Y}}_{ik}\) is the original yield for each experimental unit in block \(k\).

Bootstrap resampling for EONR calculation

Since the data used to assign values to the parameters was collected in two distinct yield zones, the null hypotesis \({H}_{0}: {EONR}_{high}-{EONR}_{low}=0\), againts the precision agriculture hypotesis \({H}_{1}: {EONR}_{high}-{EONR}_{low}\ne 0\) was tested. To this end, EONR was calculated by utilizing bootstrap resampling techniques. This methodology involved resampling the data within each treatment group independently, maintaining the treatment-specific structure of the data. We performed this resampling procedure iteratively to create multiple bootstrap samples (10,000 samples) for each treatment group. The sampling frame was stratified by treatment within the yield zones. By repeatedly resampling and calculating EONR for these simulated datasets, we aimed to account for the variability and uncertainty inherent in the original data. Subsequently, the bootstrap samples were compared using t-statistics to assess the significance of the differences between the EONR values derived from the high and low yield zones. This approach was used to determine whether the observed differences were likely due to true differences in the EONR values or if they could have occurred by chance.

The t statistics to assess the significance of the differences between the EONR values derived from the high and low yield zones using bootstrap techniques showed that the EONR from the high yield is higher than the EONR from the low yield zone (\({mean}_{high}=218.00\) kg/ha, \({mean}_{low}=181.11\) kg/ha, \({mean}_{high} -{mean}_{low}= 36.89\), \({p}_{value} <0.001\), confidence interval [36.54–37.26]). This provides a reason to analyse the data from the different yield zones separetely. Thus, data was simulated using only estimates from the high yield zone and covered a range of nitrogen rates, starting from 0 and going up to 260 kg/ha, with increments of 20 kg/ha, spanning the range of typical nitrogen rates used in agricultural settings.

Data generation using simulation

The data generating process was framed using a standard approach for planning and reporting simulation studies, the ADEMP structure: Aims (A), Data-generating mechanisms (D), Estimands (E), Methods (M), Performance measures (P) (Morris et al., 2019). The aim is to understand how the average model performs using data generated by producing parametric draws from a known model. The method is to compare the average model with the selected functional equation, namely, Mitscherlich, Linear-Plateau, Quadratic, and Quadratic-Plateau, using the EONR as our estimand (the parameter of interest). The performance measures are the bias, variance and mean squared error.

To generate synthetic datasets parametric bootstrap approach was employed (Genest & Rémillard, 2008), reflecting different scenarios of crop yield response. Thus, the parameters (coefficients) were estimated using the selected functional equation as the underlying model.

Four simulation procedures were performed to generate yield data using the selected functional equation (Mitscherlich, Linear-Plateau, Quadratic, and Quadratic-Plateau as the underlying model). In each simulation, each nitrogen level had six repeatitions. The data generation process was repeated for each model 10,000 times to create 10,000 replicates for each model. To simulate real-world variability in the data, we introduced random noise in the data by using a residual standard deviation estimated using the real data. Simulations in which the models did not converge or that resulted in EONRs with little practical application (i.e., EONR > 350 kg/ha) were discarded. To avoid convergence problems, especially while estimating the underlying true model, the true model’s parameter values (i.e. parameters estimates from the real data used to simulate the data) were used as starting values. This was done to provide the best possible approximation of the true model. As such, for data generated by the Mitscherlich function, 68 simulations (out of 10, 000) were eliminated due to convergence problems. For data generated by the Quadratic and Quadratic-plateau functions models did not converge in 92 and 87 simulations, respectively. The analysis done in the data simulated using linear plateau parameters presented the higest number of simulations in which the models did not converge (175).

Calculating economic optimum nitrogen rate (EONR)

Economic Optimum Nitrogen Rate (EONR) is the nitrogen rate at which the net return is maximized. It is a fundamental concept in agricultural economics that aims to balance crop yield potential with the cost of nitrogen fertilizer application. This approach plays a pivotal role in optimizing agricultural practices to achieve maximum yield while considering economic efficiency (Sieling et al., 2023). At its core, the EONR calculation is rooted in the principles of optimizing input-output relationships in agriculture. By identifying the nitrogen application rate that yields the highest net return, farmers can make informed decisions that enhance both productivity and profitability. Notably, in addition to agronomic factors (e.g. yield intercept, slope, plateau level), economic considerations are paramount for the EONR calculation (Nigon et al., 2019). The price of the crop (\(p\)) directly affects the economic value of yield increase, as it influences the market value of the produce. In addition to crop prices, the cost of nitrogen fertilizer (\(r\)) is a crucial economic factor. This cost directly affects the trade-off between yield enhancement and input expenses.

The EONRs estimation was based on each of the selected models for each simulated dataset using nonlinear regression techniques, such as the the crop yield \({y}_{i}=f\left(N\right)+{\epsilon }_{1}\). All functions considered in this paper (except the linear-plateau) have a smooth first and second derivative over the domain N ≥ 0 and hence the EONR can be obtained by standard differential calculus, involving setting the first derivative of the net return function p*f(N) – rN with respect to N to zero and solving for N.

The linear-plateau function can be denoted as follows:

$$f\left( N \right) = \min \left( {\beta _{0} + \beta _{1} N,~\,\mu } \right)$$
(3)

where \(\mu\) is the plateau (this represents the maximum attainable yield without additional nitrogen input). Parameter \(\beta _{0}\) represents the intercept, while parameter \(\beta _{1}\) signifies the slope of the line until it reaches the breakpoint, which is the point where the linear response intersects with the plateau line.  The breakpoint \(\left({N}_{s}\right)\) for linear-plateau can be calculated as:

$${N}_{s} = \frac{\mu - \beta _{0} }{\beta _{1} }$$
(4)

It is important to recognize the linear-plateau model’s distinct characteristics. The presence of a break point introduces a non-smooth first derivative, necessitating careful treatment during optimization processes. In certain scenarios, the EONR can be zero. Specifically, when the price of nitrogen (\(r\)) exceeds the product of crop price (p) and the slope of the linear plateau (\({\beta }_{1}\)) (\(r>p{\beta }_{1}\)), EONR becomes zero. This insight carries implications for sustainable nitrogen management, particularly in contexts where excessive nitrogen application is a concern. In all scenarios considered in our simulations, however, \(r<p{\beta }_{1}\), so that EONR > 0. In such condition, the EONR is equivalent to the breakpoint. Hence, the EONR for the linear-plateau was estimated by Eq. (4), setting EONRLP = Ns.

The quadratic and quadratic-plateau functions are, respectively, specified as:

$$f\left( N \right) = \beta _{0} + \beta _{1} N + \beta _{2} N^{2}$$
(5)
$$f\left( N \right) = \left\{ {\begin{array}{*{20}c} {\beta _{0} + \beta _{1} N + \beta _{2} N^{2} ,~~if~N < ~N^{*} } \\ {\beta _{0} + \beta _{1} N^{*} + \beta _{2} \left( {N^{*} } \right)^{2} ,~~if~N \ge ~N^{*} } \\ \end{array} } \right.$$
(6)

where \(N\) is the nitrogen rate, \(N^{*} = - \frac{{\beta _{1} }}{{2\beta _{2} }},\,\beta _{0}\) is the yield at the origin (i.e the yield achieved when no nitrogen is applied), \({\beta _{1} }\) is the linear regression coefficient for N, and \({\beta _{2} }\) is the quadratic regression coefficient. The EONR is the rate where the parabola for the expected yield have the gradient \(\frac{r}{p}\) (Bachmaier & Gandorfer, 2009). Thus, we need to find the nitrogen rate such that:

$$f^{\prime}\left( N \right)~ = ~\frac{r}{p}$$
(7)

Since \(f^{\prime}\left( N \right) = 0 + \beta _{1} + 2\beta _{2} N\) the first-order condition that that must be satisfied to ensure the profit-maximazation nitrogen rate is:

$$\beta _{1} + 2\beta _{2} N = \frac{r}{p}$$
(8)

After solving N, the EONR is calculated using the equation below:

$$EONR_{Q} = \frac{{\frac{r}{p} - \beta _{1} }}{{2\beta _{2} }}$$
(9)

where \(r\) is the price of nitrogen, \(p\) is the price of crop. In this case, the term \(\left( {\frac{r}{p}} \right)\) represents the relative cost of nitrogen compared to the crop price, which as estimated at 0.086. Thus, this formulation reflects the fact that, as the relative cost of nitrogen increases, the economic incentive to apply more nitrogen decreases. Higher relative nitrogen costs would lead to a lower optimal nitrogen application rate. Finally, the Mitscherlich function is expressed as:

$$f\left( N \right) = \beta _{0} + (\beta _{1} - \beta _{0} )e^{{ - \beta _{2} N}}$$
(10)

Since \(f^{\prime}\left( N \right) = (\beta _{1} - \beta _{0} )e^{{ - \beta _{2} N}} \left( { - \beta _{2} } \right)\), the the first-order condition is:

$$(\beta _{1} - \beta _{0} )e^{{ - \beta _{2} N}} \left( { - \beta _{2} } \right) = ~\frac{r}{p}$$
(11)

Then,

$$e^{{ - \beta _{2} N}} = ~\frac{{ - r}}{{p(\beta _{1} - \beta _{0} )\beta _{2} }}$$
(12)

After solving for N, the corresponding EONR was calculating using Eq. 13:

$$EONR_{M} = - \frac{{\ln \left[ {\frac{{ - r}}{{p(\beta _{1} - \beta _{0} )\beta _{2} }}} \right]}}{{\beta _{2} }}$$
(13)

Model-averaging approach

The commonly used method of choosing the “best” model based on the information criteria, might not always be obvious (Miguez & Poffenbarger, 2022; Piepho & Williams, 2021), which has a large impact on the estimated EONR, resulting in under- or over-estimation of EONR (Lyons et al., 2019). In such cases, model averaging might be a good option for obtaining predictions of EONR (Buckland et al., 1997; Burnham & Anderson, 2002). Therefore, we study the potential of model-averaging approach in estimating EONR and the expected yield as accurately as possible. This was done considering two scenarios. The first scenario includes the the true model from the set of models considered in the averaging process, while the second was done excluding the true model. To estimate the model-averaged EONR for a given set of candidate models, we used the Eq. 14 (Burnham et al., 2011):

$$\widehat{{\bar{\theta }}} = \mathop \sum \limits_{{i = 1}}^{M} w_{i} \hat{\theta }_{i}$$
(14)

where \({{\hat{\theta }}}_{i}\) is the estimated EONR obtained from model i and \(w_{i}\) is the weight associated with that model (Buckland et al., 1997; Fletcher, 2018). The weight associated with model i is calculated as:

$$w_{i} = \frac{{\exp \left( { - \frac{1}{2}\Delta AIC_{i} } \right)}}{{\mathop \sum \nolimits_{{i = 1}}^{M} \exp \left( { - \frac{1}{2}\Delta AIC_{i} } \right)}}$$
(15)

where \(\varDelta {AIC}_{i}\equiv {AIC}_{i}-\text{m}\text{i}\text{n}\left({AIC}_{i}\right)\). The Akaike Information Criterion is given by \(\text{A}\text{I}\text{C} = -2 \text{l}\text{o}\text{g} \widehat{L}+ 2p\), where \(\widehat{L}\) is the the maximized value of the likelihood function of the model under consideration and \(p\) is the number of parameters of the model.

Mean squared error and bias estimation

We obtained the misspecification bias (Eq. 16) that arises in estimating EONR under model i by calculating the absolute difference between the the sample mean of the estimates of EONR (\(\widehat{{\uptheta }}\)) for model i and the true EONR (\({\uptheta })\), obtained from estimating one of the candidate models from the real data (Eq. 14). For example, for the data generated by the quadratic function, we took the EONR determined for the quadratic function using the real data as the true EONR. We also assume that the EONR estimated by the corresponding underlying function will be the closest possible to the true EONR (Miguez & Poffenbarger, 2022):

$${\text{Bias}}\left( {\hat{\theta }_{i} } \right) = \left( {\frac{1}{n}\mathop \sum \limits_{{j = 1}}^{n} \hat{\theta }_{{i\left( j \right)}} } \right){\text{~}} - \theta$$
(16)

where \(n\) is the number of simulations, \({\hat{\theta }_{{i\left( j \right)}} }\) represents the estimated EONR for model i in the \(j\)th simulation run and \(\theta\) is the true EONR. The mean squared error was partitioned into variance and squared bias as follows (Molina & Rao, 2010):

$${\text{MSE}}\left( {\hat{\theta }_{i} } \right) = {\text{Var}}\left( {\hat{\theta }_{i} } \right) + {\text{Bias}}\left( {\hat{\theta }_{i} } \right)^{2}$$
(17)

where \({\hat{\theta }}\) represents the estimated EONR and \({\text{Var}}{\hat{\theta }}\) is the variance for the individual functions, given by:

$${\text{Var}}\left( {\hat{\theta }_{i} } \right) = \frac{1}{{n - 1}}\mathop \sum \limits_{{j = 1}}^{n} \left( {\hat{\theta }_{{i\left( j \right)}} - \left( {\frac{1}{n}\mathop \sum \limits_{{j = 1}}^{n} \hat{\theta }_{{i\left( j \right)}} } \right){\text{~}}} \right)^{2}$$
(18)

As outlined by Buckland et al. (1997), in the case of averaged model, the bias arising from any single model becomes the contribution of the model selection uncertainty. Hence, they proposed the variance estimator:

$${\text{Var}}\left( {\hat{\theta }} \right)_{{average~model}} = \left\{ {\mathop \sum \limits_{{1 = 1}}^{M} w_{i} \sqrt {{\text{Var}}\left( {\hat{\theta }_{i} } \right) + \left( {\hat{\theta }_{i} - \widehat{{\bar{\theta }}}} \right)^{2} } } \right\}^{2}$$
(19)

In the equation above, a variance component for model selection uncertainty \({\left({\widehat{{\uptheta }}}_{i}-\widehat{{\bar{\theta }}}\right)}^{2}\)is incorporated, where \(\widehat{{\bar{\theta }}}\) is the the estimated EONR averaged over all models (Burnham & Anderson, 2002). We included this variance estimator in the simulation in order to assess whether it provides a valid astimator of the MSE. In each simulation run, the variance \(\text{V}\text{a}\text{r}\left({\widehat{{\uptheta }}}_{i}\right)\) in (19) was calculated by squaring the asymptotic standard error values obtained from the R package nlraa: nonlinear regression for agricultural applications (Miguez, 2021).

Real-world experimental data analysis

Real-world experimental data was used also used to estimate EONR for each model using actual observations. The selected models (Mitscherlich, Linear-Plateau, Quadratic, Quadratic-Plateau) were fitted to the experimental data, including the average model. The EONRs were calculated for each fitted model using the procedures previously described. The analysis took into account the block effect. This was done by using Eqs. (1) and (2). Confidence intervals for the EONR estimates were determined through the utilization of bootstrapping methods. The bootstrapping technique was employed to compute estimates for the nonlinear models by resampling the residuals. The boot package in the R programming environment facilitated this process (Canty, 2002). A total of 10,000 bootstrap replicates were generated, serving to effectively capture the variability inherent in the estimate. Subsequently, the boot.ci function, a part of the boot package, was utilized to calculate the confidence intervals based on the bias-corrected and accelerated (BCa) bootstrapped confidence intervals. This comprehensive approach ensured robust and reliable estimation of the confidence intervals for the EONR, offering a thorough understanding of the potential range of the estimate (Hernandez & Mulla, 2008). Standard error of the EONRs were also computed through bootstrapping procedures.

Results and discussion

Simulation results

The results from our simulations show that the linear-plateau funtion is, in general, negatively biased compared to the other models (Table 2). This is in line with the results found by Miguez and Poffenbarger (2022). The bias for the linear-plateau function ranged from − 65.0 to − 0.4 kg N ha−1 (Table 2). The quadratic-plateau function presented smaller bias, when fitted with data generated by the quadratic function. Likewise, the quadratic function presented less bias for data generated by the quadratic-plateau function. Overall, the EONR estimated by the corresponding underlying function was the closest possible to the true EONR, as expected.

Nevertheless, the Mitscherlich function presented slightly higher bias when compared to the average model for data generated using the Mitscherlich function parameters estimates, mainly when the true model is included in the analysis. Overall, the average model produced small biases compared to the individual functions (except when compared to the bias estimates of the data generating function) and tend to have smaller variance and mean squared errors. When the true model is included in the averaging process, the mean of the variances calculated based on Eq. (19) were slightly higher but not too distant to the MSE calculated using Eq. (17) and presented in Table 2. These mean variances are 175.25, 52.19, 598.6 and 273.34 when the Linear-Plateau, Mitscherlich, Quadratic, Quadratic-plateau are the true models, respectively. However, when the set of fitted models does not include the true model, the variances are 853, 183, 937 and 582 when the Linear-Plateau, Mitscherlich, Quadratic, Quadratic-plateau are the true models, respectively.

Small biases were reported in different studies on the performance of model averaging procedures (Lukacs et al., 2010; Wheeler & Bailer, 2007, 2009). Thus, it seems to be a promising solution to address the uncertainity with respect to EONR estimation, especially because EONR estimation is sensitive to model choice (Buckland et al., 1997). The quadratic and quadratic-plateau were, on average, less biased compared to Mitscherlich and linear-plateau. Previous research has shown that the quadratic-plateau is often a superior model to describe yield response and EONR estimation in wheat and corn production (Cerrato & Blackmer, 1990; Ortuzar-Iragorri et al., 2010; Scharf et al., 2005). However, in the presences of yield penalties with higher N rates, quadratic yield response functions have been found to be superior (Meyer-Aurich et al., 2010b).

The uncertainties in yield predictions across EONRs are also central. Quantifying MSE and variance provides valuable insights for agronomists, highlighting potential variations in crop responses to nitrogen fertilization. Understanding these uncertainties is crucial for informed decision-making in optimizing nitrogen application rates. Additionally, the exploration of rebound effects contributes to a comprehensive perspective, offering practical guidance for sustainable and efficient nitrogen management practices in real-world agronomic contexts.

Table 2 EONR (kg N ha−1), bias, variance and mean squared errors (MSE) for the different models
Full size table

For purposes of comparison, we present the results of EONR, bias, variance and mean squared, when the set of fitted models does not include the true model. The exclusion of the true model increases the biases, variances and mean squared errors for the average model. However, compared to other model, the average model presents reduced bias. Notably, exception occurs when the true underlying model is quadratic and quadratic plateau, in which both the quadratic and quadratic palteu are less biased than the averaged model, similar to the results obtained when the true model in included. The average model consistently exhibits lower bias compared to the linear-plateau and Mitscherlich models. This implies that model averaging can still be useful even when the true model is not in the set of fitted models. It can be used to reduce model selection bias effects on coefficient estimates (Burnham & Anderson, 2002) but it is required that the set of fitted models includes at least one that approximates well the true model (Buckland et al., 1997). Hence, if there is no a firm foundation of prior knowledge regarding the most appropriate model, it may be essential to use a rich set of potential models (Piepho & Williams, 2021). Luckily, a promising set of yield response models have alredy been documented in the literature on crop yield response to nutrient application (Dhakal & Lange, 2021; Lyons et al., 2019).

Results from the experimental data

The confidence intervals for each yield function are shown in Fig. 2 and listed in Table 3. The results reveal that estimates of EONR for the different response functions differ substantialy but there is considerable overlap among the confidence intervals. The overlaping, however, does not imply non-superiority of one model to another in estimating the EONR (Payton et al., 2003). The linear-plateau presented the lowest estimates for EONR. This was also reported by Jaynes (2011). Likewise, some model presented wider confidence intervals compared to others. For example, despite having similar EONR, the quadratic and quadratic-plateau functions had different confidence interval widths. As expected, the quadratic-plateau presented a wider confidence interval. Although narrower confidence band widths are a desirable characteristic as they reflect greater precision in EONR estimates for each response function, these widths do not indicate the accuracy of the various EONR estimates in determining the actual EONR value (Lyons et al., 2019). This is because they do not take model uncertainty into account. The confidence interval widths for the Mitscherlich function where the widest (159.1 kg/ha and 106.6 kg/ha for the high and low yield zones, respectively). This is consistent with the high standard error estimates (Table 3) and the simulation results in which the Mitscherlich function presented the highest mean squared error and variance (Table 2). This result supports the idea that choosing the “best” model based on model ranking with information criteria may be misleading (Lukacs et al., 2010; Yang, 2006) since the Mitscherlich function presented the smallest AIC and and the highest weight (Table 3). In fact, one of the advantages of a model averaged estimator is that it has a better precision and reduced bias compared to the estimator from the selected best model (Burnham & Anderson, 2002).

Given the range of variations in the calculated EONR and the associated uncertainties determined from the response functions (Henke et al., 2007), it is challenging to definitively identify which function provides the most accurate estimation of the true EONR. Particularly noteworthy is that the confidence interval of the averaged model encompasses all other functions’ estimates and exhibits an overlap with each function’s interval. Consequently, owing to its relatively reduced bias (Table 2), selecting the averaged model could be a more cautious approach compared to opting for any of the individual models like linear-plateau, Mitscherlich, quadratic, or quadratic plateau. In the context of precision farming, where accurate and site-specific recommendations are paramount, the integration of model averaging introduces a layer of complexity. The bootstrap confidence interval, with a width of 78 kg/ha in high yield zone raises legitimate concerns about the precision of nitrogen rate determinations (Jaynes, 2011). The practical utility of model averaging hinges on the assumption that one of the competing models closely approximates the true underlying model (Piepho & Williams, 2021). However, in real-world scenarios characterized by diverse and dynamic field conditions, this assumption may not hold uniformly.

Despite these challenges, the process of model averaging remains a valid and valuable procedure to select the best approximating model from a set of competing models (Buckland et al., 1997). It allows for a nuanced consideration of uncertainties and performance metrics, providing a comprehensive view of potential outcomes. In our simulation study, where the underlying model was both known and incorporated into the analysis, the favorable performance exhibited by the averaged model underscores the importance of accurate model specification. As precision agriculture continues to evolve, future applications of model averaging should carefully weigh the trade-offs between model complexity, practical feasibility, and the dynamic nature of field conditions.

Moreover, agricultural practices have multifaceted impacts that extend beyond economic considerations, encompassing energy usage and environmental sustainability. While our study primarily focuses on the economic dimensions, we recognize the broader relevance of energy and environmental factors in the context of fertilizer application. Thus, future research could explore the intricate interplay between economic efficiency, energy consumption, and environmental implications in agricultural systems.

Fig. 2
figure 2

Yield response functions, EONRs and bootstrap confidence intervals

Full size image
Table 3 Computed EONR (kg N ha−1), AIC, model weights, standard errors (SE) and 95% confidence bootstrap intervals widths for yield response functions
Full size table

Conclusion

In this study, a comparative analysis of biases, Economic Optimum Nitrogen Rate (EONR) estimations, variance, mean squared erros and confidence intervals was conducted to gain insights into the performance disparities among distinct functional forms, as well as to explore the potential advantages of employing model averaging to optimize nitrogen application in crop cultivation. Overall, the simulation results indicate that notable biases exist when contrasting the various yield functions with the averaged model, particularly in the case of the linear-plateau model. Nonetheless, upon scrutinizing the empirical data, it becomes evident that the calculated confidence intervals for the averaged model tend to overlap with the estimated ranges of all functions. This finding suggests that the averaged model might be less biased and well-suited for the determination of EONR. The practical utility of model averaging can be realized if models that approximate the true model are incorporated into the analysis. This has implications for further research on optimal fertilizer rates in agriculture. Instead of arguing for the one or other functional form researchers are encouraged to follow the model averaging approach to prevent creating biases on optimal fertilizer levels and marginal fertilizer costs. This may help to better identify the economic potentials of site-specific fertilizer management strategies. The unbiased estimation of marginal fertilizer costs is becoming more and more important as the reduction of fertilizer supply and surplus is needed to contribute to reduced emissions of reactive N. The model averaging approach can help to better communicate these costs among stakeholders.