Toward a standardized statistical methodology comparing optimum nitrogen rates among management practices: A bootstrapping approach

There are a range of approaches to compare differences between or among optimum nitrogen (N) fertilizer rates resulting from different management practices; however, this goal lacks statistical standardization. To provide the statistical rigor needed to give clear recommendations for greater or less N need based on specific management practices, we propose a bootstrapping approach that resamples residuals with replacement. While bootstrapping is not new to data processing in agronomic fields, we provide an example of how to conduct residual‐resampled bootstrapping with nonlinear regression to identify differences in response curves, optimum N rates, and maximum yields using the FertBoot package in R. Our example dataset provides clear evidence of the value of the bootstrapping approach, as it can aid in determining significant differences between even relatively small differences in optimum N rate. We encourage adoption of this approach as a way to accurately evaluate differences in optimum fertilizer levels between or among treatments to better inform future agronomic decision making.

Without the statistical ability to determine if there are significant differences between or among optimum N rates or maximum yield as determined from the nonlinear models, we cannot properly assess the validity of the calculated fertilizer replacement value associated with an organic input or management practice. This highlights a need for a standardized statistical approach that allows a meaningful comparison between or among optimum N rates from different treatments.
The difference in optimal N rate is an important variable of interest, but traditional statistical approaches do not allow researchers to identify if this difference is statistically true. An approach using bootstrapping residuals with nonlinear regression offers that ability. Bootstrapping is a statistical method used to address the problem of finding a sampling distribution of a variable when the probability distribution is unknown (Efron, 1979). To estimate the sampling distribution of a variable, the bootstrap approach uses independent resampling of the residuals with replacement (resampling usually upwards of 1,000 or more times) to produce new datasets. Each bootstrapping run develops a new response curve for each treatment by calculating a new Y value for each X value by adding the original predictive value to a randomly selected residual value. In an example with nonlinear regression modeling, we obtain model coefficients, plus the optimum N rate and maximum yield (the plateau yield) from the new dataset. A randomization test can then be applied to the bootstrap sample to determine if there were differences between or among optimum N rates and maximum yield (Davidson & Hinkley, 1997).
Bootstrapping is not new to agronomic research. For example, in agronomic-focused journals, we see it used with metaanalyses (Wortman, 2016;Wortman et al., 2017;Zhang et al., 2020), field-crops research (King & Blesh, 2018;Wang et al., 2013;Zheng et al., 2009), and ecological modeling (Heikkinen & Mäkipää, 2010), in conjunction with artificial neural networks (Zeng et al., 2016), as cross-validation of variables (Corstanje et al., 2009), and in spectroscopy (Yang et al., 2020), and hyperspectral methods (Kawamura et al., 2013;Li et al., 2014) in order to determine confidence intervals or compare regression results. Bootstrapping has also been used to determined confidence intervals for optimum fertilizer rates (Hernandez & Mulla, 2008;Qin et al., 2018). The bootstrapping technique proposed here can be used in any study where researchers want to identify how differences in management (e.g., tillage, cover cropping, crop rotation, N fertilizer source) affect optimum N rates. This approach is specifically valuable to ensuring that the fertilizer replacement value of input or management changes are real. While the bootstrap-

Core Ideas
• A standardized approach to determining statistical differences in optimum N rates is needed. • A bootstrapping approach that resamples residuals with replacement is a useful approach. • The FertBoot package in R is an available tool to enable this analysis.
ping approach allows for comparison of nonlinear regression equations, it does rely on the following assumptions: (a) the residuals of the original model can mimic the variability well enough and (b) the dependence within the same block can be ignored. This approach has some limitations: it is computationally intensive and some researchers may feel uncomfortable with the slight variations in the results (confidence intervals and P values) due to randomness in the analyses.
Here, we provide an example of how to conduct residualresampled bootstrapping with nonlinear regression to identify differences in response curves, optimum N rates, and maximum yields from a green manure cover crop study using the FertBoot package in R (Ma & Francis, 2021).

Experimental databases
An unpublished N response dataset was used for this analysis that compares the N replacement value of red clover (Trifolium incarnatum L.) as a green manure cover crop. The study evaluates red clover in a corn-corn (CC) rotation. In this CC study, corn was planted on 25 May 2017 and red clover was interseeded on 26 June 2017. In 2018, corn was planted on 24 May and red clover was terminated on 1 June with 2,4dichlorophenoxyacetic acid or dicamba with glyphosate. In 2018, urea coated with a urease inhibitor (Agrotain [Koch Agronomic Services]) was surface applied at eight different rates (0,45,90,135,180,225,270, and 315 kg-N ha -1 ). Red clover was not interseeded again in 2018. The experimental design was a randomized, complete block, split-plot design, with two whole-plot treatments (with and without red clover) and N rates as the split-plot factor. The study was conducted at the University of Wisconsin Arlington Agricultural Research Station on Plano silt loam (fine-silty, mixed, superactive, mesic Typic Argiudolls) soil. Raw yield data from each study are reported in Table 1 (corn yields reported at 15.5% moisture).
T A B L E 1 Nitrogen fertilizer rates and corn grain yield for each replicate (R) following no cover crop (None) or red clover (RC) in 2019 the corn-corn rotation study. Red clover was interseeded in 2018 and terminated before corn planting in 2019

Statistical approach
In our analysis, we compared optimum N fertilizer rates from quadratic-plateau models. Preliminary analysis comparing quadratic, linear-plateau, and quadratic-plateau curves (NLIN in SAS version 9.4 [SAS Institute, Cary, NC]) and confirmed with easynls package in R (version 4.0.3; R Core Team, 2020)] was used to identify that quadratic-plateau response curves (Eq. 1) had the lowest RMSE among models for these datasets. This was expected as quadratic-plateau models have been previously shown to be the model of best fit for corn N response analysis (Cerrato & Blackmer, 1990). All data for each cover crop treatment (eight N rates and four replicates) were used in the analysis. Quadratic-plateau parameter estimates and standard errors and 95% confidence intervals of the quadratic-plateau model are produced with the easynls package. Results of this analysis are also used to determine the optimum N rate and the maximum yield. For a quadratic plateau model, where x ij and y ij denote N rate and yield of the jth observation with treatment i, , j = 1, ..., n i , the model is where x m,i = −b i /2c i is the optimal N rate and y m,i = a i + b i x m,i + c i x 2 m,i is the maximum yield response to N rate for the ith treatment. We used the log-likelihood ratio test to determine if quadratic-plateau models (produced from easynls) were different between the two cover crop treatments. Models were considered different if at least one coefficient estimate was significantly different (P < .05) between models (Bates & Watts, 1988).
The FertBoot package in R (Ma & Francis, 2021) was used to determine significant differences in model coefficients, optimum N, and maximum yield between cover crop treatments. Within FertBoot, quadratic plateau models were fitted, and residuals of the model were resampled (with replacement) 1,000 times to produce a population estimate dataset of a, b, and c coefficients, and of optimum N rates and maximum Y values. If an optimum N value greater than the largest N rate used in the study was determined, that value was replaced with the largest N rate used in the study, which was then used to recalculate the maximum yield. The parameter estimates of the bootstrap model are calculated from the average of estimates among bootstrapped samples. The bootstrapping approach will produce different coefficients for the plateau curves compared with the original model but ultimately produces a more robust measurement of the variability. The ggplot2 package in R was used to plot the bootstrapped models and their 95% confidence intervals with the original data points. Density curves were produced using the bootstrapped output of optimum N fertilizer values with ggplot (geom_density) in R to visualize differences in optimal N rates. Confidence intervals (95%) of the mean of each response variable (model coefficients, optimum N rate, and maximum yield) were constructed from the 2.5 and 97.5 percentiles (Davison & Hinkley, 1997). Differences between optimum N rates and maximum yields were determined using a two-sample bootstrap test (using 1,000 replicates) at P < .05 (Boos, 2003).

RESULTS AND DISCUSSION
The dataset for determining the N replacement value of red clover in a CC rotation produced quadratic-plateau curves with large 95% confidence intervals (Table 2) of the coefficients for both the none and red clover treatments. However, results from the log-likelihood ratio test indicate that the two models are significantly different (P = .0048) based on the original dataset. The optimum N rates differed between the two treatments by 20.2 kg ha -1 using the original curves, and by 16.9 kg ha -1 using the bootstrapped model ( Figure 1, Table 2). The two-sample bootstrap test determined that the difference of 16.9 kg ha -1 was significantly different than zero (P > .001), and the density plots revealed a clear separation of peaks ( Figure 2). Evidence of statistical differences allows us T A B L E 2 Parameter estimates from quadratic-plateau response curves determined from the original analysis and bootstrapped analysis for no cover crop (None, R 2 = .68 from original analysis) and red clover (RC) cover cropped (RC, R 2 = .72 from original analysis) treatments in 2018 in the corn-corn study to discuss the results as being "real" rather than simply speculating based on the 95% intervals of the model coefficients and relative size of the difference. We can then appropriately assess the agronomic and economic relevance of this significant difference. From an outreach perspective, as we build datasets comparing the yield or N effect of red clover cover crops, we can decide to include a difference (e.g., 16 kg ha -1 ) into the average response, or if we do not determine a difference in optimum N rate, we can include a difference of zero. Interestingly, this approach determined a significant difference between the maximum yields of 10.5 Mg ha -1 (none) and 10.0 Mg ha -1 (red clover), providing the ability to assess the impact of the cover crop on crop production (in this case a yield reduction).

CONCLUSIONS
There does not currently appear to be a universally accepted or used approach to determine statistically different optimum N rates between or among quadratic-plateau models in agronomic research. The large variation in fertilizer response data in field studies, and the generally limited amount of individual observations that can be produced, can make some methods for comparisons between treatments of interest statistically inappropriate or cumbersome using traditional methods of evaluation. Here we proposed the approach of bootstrapping residuals with replacement as a functional way to determine statistical differences in optimum N fertilizer rates and maximum yields resulting from different management practices. This bootstrapping approach is statistically sound and can be the primary means of analysis for studies comparing optimum N rates with non-linear regression. To facilitate ease of exploration of the method proposed here, the R package FertBoot has been published to the CRAN repository. Our example analysis demonstrated that the bootstrapped quadratic-plateau models produced similar coefficients to the original model and also produced large 95% confidence intervals. However, while alternative methods may rely too heavily on the 95% confidence interval or standard error terms to determine if differences between optimum N fertilizer rates are significantly different, the two-sample bootstrap test following bootstrapping residuals of the fertilizer response curve models determined significant differences between the optimum N rates of comparison. Along with use of this statistical approach, we also promote an increase in the number of field F I G U R E 1 Fertilizer response curve for corn following no cover crop (None) or red clover (RC) in the 2018 corn-corn study. The colored bands are the bias corrected (BCa) 95% confidence intervals determined by bootstrapping residuals. Results from field experiments were fit to quadratic-plateau models, which were subjected to the bootstrapping residuals procedure.
F I G U R E 2 Density curves showing the bootstrapped outputs for optimal N fertilizer rates in the no clover (None) and clover treatment (RC) from the 2018 corn-corn study.
replicates and N rates used in similar studies to those analyzed here. Furthermore, we envision this approach could be applied to many research questions beyond those considered here, such as in developing regional fertilizer recommendations based on soil, management, or a combination thereof.

D A T A AVA I L A B I L I T Y S T A T E M E N T
Supplementary materials related to running FertBoot can be found at https://ruarklab.soils.wisc.edu/publications/ supplementary-materials/.

A C K N O W L E D G M E N T S
The authors would like to thank Lorna Yen for their article "An introduction to the bootstrapping method" (https://towar dsdatascience.com/an-introduction-to-the-bootstrap-method-58bcb51b4d60) that served as a source of inspiration for this paper. Funding for this research was provided by the Wisconsin Fertilizer Research Council, Project 280-15 and 282-16.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.