Agronomy Journal 93:531-539 (2001)
© 2001 American Society of Agronomy
PRODUCTION PAPERS
Statistical Methods for Predicting Responses to Applied Nitrogen and Calculating Optimal Nitrogen Rates
David Makowskia,
Daniel Wallacha and
Jean-Marc Meynardb
a INRA, B.P. 27, 31326 Castanet-Tolosan Cedex, France
b Laboratoire d'agronomie, INRA INA P-G, BP 01 78850 Thiverval-Grignon, France
Corresponding author (makowski{at}toulouse.inra.fr)
Received for publication February 15, 2000.
 |
ABSTRACT
|
|---|
Models of response to applied N can be useful for deriving improved N dose recommendations. Here we show how response model parameters can be estimated and how model predictions and model N dose recommendations can be evaluated. For parameter estimation, we use a statistical approach based on random parameter models. Two methods for evaluating models are applied. The first method is to calculate mean squared error of prediction (MSEP) by cross validation, and the second is to perform nonparametric regressions to evaluate the profitability of calculated optimal N rates. The proposed methods are used with a data set consisting of 37 winter wheat (Triticum aestivum L.) experiments. Different functions taking into account end-of-winter mineral soil N are evaluated. The results show that the different functions all have similar MSEP values for predictions of yield and grain protein content and lead to N recommendations of similar profitability. However, there are substantial differences in MSEP for residual mineral N at harvest. One of these models is then compared with a model that does not include any site-year characteristic and with a model that does not have random parameters. We find that using the model without a site-year characteristic leads to predictions that are less accurate and optimal N rates that are less profitable by F17 to F105 ha-1. Another result is that the gross margin obtained with the optimal N rates calculated using the model without random parameters is lower by F438 to F550 ha-1.
Abbreviations: E, expectation F, French franc LP, linear-plus-plateau MSEP, mean squared error of prediction PL, plateau-plus-linear PQ, plateau-plus-quadratic Q, quadratic QP, quadratic-plus-plateau
 |
INTRODUCTION
|
|---|
MODELS predicting the responses of winter wheat to applied N can be useful for making better N dose recommendations to farmers. Numerous functions have been proposed for describing yield response to applied N (Anderson and Nelson, 1975; Cerrato and Blackmer, 1990; Bullock and Bullock, 1994; Colwell, 1994), but few have been proposed for describing the responses of grain protein content (Murray and Nunn, 1987; Fowler et al., 1989) and residual mineral soil N at harvest (Jauregui and Paris, 1985). In a recent study, Makowski et al. (1999) defined several sets of functions that relate yield, grain protein content, and residual mineral N at harvest to N applied to winter wheat. The authors showed that their functions give satisfactory fits to the data when fitted site-year by site-year. The purpose of this paper is to show how the functions of Makowski et al. (1999) can be used for predicting the responses to applied N and calculating optimal N rates as a function of site-year characteristics. Solving this problem supposes the definition of appropriate statistical methods for estimating model parameters and evaluating model quality using standard N fertilizer trials. Such methods are presented in this paper and are applied to a set of winter wheat experiments in the Paris Basin of France. The proposed methods are quite general and can be used with different response functions.
Parameter estimation in N response models is generally based on ordinary least squares (e.g., Mombiela et al., 1981; Sain and Jauregui, 1993). This method is easily implemented. However, the underlying assumption that model errors are independent is unrealistic for data like N response data where several measurements (corresponding to different N doses) are made for the same site-year. In this paper, we propose an alternative method for parameter estimation. Our approach consists of using a statistical model called a random parameter model, which is a specific type of mixed-effects model. A random parameter model assumes that the form of the modeled response is common to all site-years but that parameter values vary between site-years. The parameters are thus treated as having some random variation around their mean values. This approach is consistent with the data that are usually observed, and it automatically introduces correlations among errors for the same site-year. This method was applied by Wallach (1995a)(1995b) to a simple quadratic (Q) model for yield as a function of applied fertilizer.
We show in this paper how random parameter models can be used with the more complex models of Makowski et al. (1999). Then we show how these models can be used for predicting the response to applied N and calculating optimal N rates. The possibility of calculating optimal N rates that depend explicitly on prices, yield, grain protein content, residual mineral N, and field characteristics is one of the major attractions of the models of Makowski et al. (1999). This is not possible with the balance-sheet method currently used in France and other countries to derive N fertilizer recommendations (Stanford, 1973; Rémy and Hébert, 1977; Neeteson, 1990; Meynard et al., 1997).
Another problem treated in this paper is that of model evaluation. Many different response models can be developed by using different functions. To select a particular model for practical use, it is necessary to define criteria for evaluating the different possible models. In numerous studies, response models are evaluated by calculating R2 values (Anderson and Nelson, 1975; Cerrato and Blackmer, 1990; Sain and Jauregui, 1993). Another approach is to study graphically the distribution of the model residuals (Cerrato and Blackmer, 1990; Bullock and Bullock, 1994). These approaches evaluate the fit of the proposed functions to past data. Often, however, the goal of developing the models is either to predict crop response to applied N for new site-years or to propose optimal N fertilizer rates. Specific tests to evaluate model quality in each of these cases have been proposed. For prediction, the MSEP is an appropriate criterion. It can be estimated using cross validation (Efron, 1983). If a model is developed to calculate N rates that are optimal for a particular objective function, for example gross margin, it is the value of the objective function when the model-based recommendations are applied that is of interest. Antoniadou and Wallach (2000) propose this criterion and show how it can be estimated. In this paper, we use MSEP and the method of Antoniadou and Wallach (2000) for evaluating our models.
These methods are applied in a case study that allows us to compare the different functions proposed by Makowski et al. (1999) to study the interest of using end-of-winter mineral soil N as an explanatory variable and to study the interest of using random parameter models for modeling responses to applied N.
 |
METHODOLOGICAL FRAMEWORK
|
|---|
Random Parameter Model
We present here a random parameter model based on one set of functions proposed by Makowski et al. (1999). It is convenient to present this model in two stages. At the first stage, the responses of N uptake, yield, grain protein content, and residual mineral N to applied N are modeled for a particular site-year. Note that N uptake is not of direct interest but is used as an intermediate variable. At the second stage, the variation among site-years is ascribed to random variation in the parameters.
First Stage
N uptake for a particular site-year is related to applied N by a linear-plus-hyperbola function:
 | (1) |
 | (2) |
where u is the measured N uptake value for a particular site-year, x the applied N dose, and
u denotes the error of the model for N uptake. This function has six parameters: C0 (N recovery coefficient), U0 (N uptake when no fertilizer is applied), B (N uptake per yield unit), YMAX (maximal yield), T (parameter that determines the asymptotic level of N uptake), and XMAX (value of applied N at which the hyperbolic part of the N uptake response begins). To ensure continuity, we require that XMAX = (BYMAX - U0)/C0, and so the number of parameters to specify in Eq. [1] and [2] is reduced to five. It is assumed that
u is independent for different N doses and has identical normal distributions,
u
N(0,
2u).
Yield is related to N uptake by a linear-plus-plateau (LP) function. It is assumed that the plateau begins at the N dose XMAX. Combining this with Eq. [1] and [2] leads to the following function for yield vs. applied N:
 | (3) |
 | (4) |
where y is the measured yield value for a particular site-year, A is a parameter, and
y denotes the error of the model for yield,
y
N
.
Grain protein content is related to the ratio (N uptake/yield) by a linear function. Combining this linear function with Eq. [1] through [4] leads to the function for grain protein content vs. applied N:
 | (5) |
 | (6) |
where p is the measured grain protein content for a particular site-year, P1 and P2 are parameters, and
p denotes the error of the model for grain protein content,
p
N
.
Finally, residual mineral N at harvest is related to applied N by a plateau-plus-linear (PL) function:
 | (7) |
 | (8) |
where r is the measured residual mineral N at harvest for a particular site-year and RMIN and R are parameters. RMIN is the minimal value of r, and R is the increase in residual mineral N per unit increase in applied N.
r denotes the error of the model for residual mineral N at harvest. Equations [7] and [8] assume that the increase in residual mineral N begins at XMAX where maximal yield is reached (Machet and Mary, 1990; Chaney, 1990; Makowski et al., 1999). Once again it is assumed that
r is independent and has identical normal distributions,
r
N
. Finally, we assume that the model errors for different responses are independent.
Second Stage
At this stage, the variability of the parameters between site-years is modeled. Let
denote the vector of the 10 parameters of the model,
= [U0, C0, T, YMAX, B, A, P1, P2, R, RMIN]t, where t is the transpose matrix operator. We assume that
has a normal distribution
 | (9) |
where µ is the vector of the expected value of
, µ = [E(U0), E(C0), E(T), E(YMAX), E(B), E(A), E(P1), E(P2), E(RMIN), E(R)]t, and
is a (10 x 10) variancecovariance matrix. The 10 diagonal elements of
are the variances of the 10 elements of
, indicating how widely they vary from site-year to site-year. The off-diagonal elements of
are the covariances between the elements of
, indicating the extent to which they tend to vary together. Some diagonal elements can be zero if the parameters are not random, and some off-diagonal elements can be zero if the parameters are not correlated.
It is possible to extend the model by modeling how µ depends on site-year characteristics. This is of interest when it is thought that part of the variability of
between site-years can be explained by the variability of some site-year characteristics. For example, some of the variability between site-years in the parameter U0 (N uptake in the absence of applied N) can generally be explained by differences in amounts of end-of-winter mineral soil N (
) (Stanford, 1973). This dependence between U0 and
can be modeled as
 | (10) |
where
1 and
2 are two fixed parameters. Field characteristics other than mineral soil N could also be easily introduced. Introducing site-year characteristics in the model increases the number of parameters in µ to be estimated. For example, if Eq. [10] is used, then µ is a function of 11 parameters because now E(U0) is not given by a single parameter but rather depends on
1 and
2.
Estimation
The method proposed by Lindstrom and Bates (1990) and Davidian and Giltinan (1995) for nonlinear mixed-effects models can be used to obtain estimators of
(vector of parameters on which µ depends),
,
2u,
2y,
2p, and
2r. This method is based on a linearization of the nonlinear model. The estimators are derived by a generalization of maximum likelihood techniques for the linear mixed-effects model to the nonlinear case. This method is now quite easily applied, for example, by using the NLME function of the SPLUS software (S-PLUS4, 1997). Estimates of
and
will be denoted
and
, and the corresponding estimate of µ will be denoted [µ(
)].
Model Predictions
Let
(x,
;
) be the value of yield calculated from Eq. [3], [4], and [10] with
y = E(
y) = 0 and with a particular value of
drawn at random from the distribution N[µ(
),
]. As the predictor of yield for a site-year with characteristic
and for applied N rate x, we use an approximation to the expectation over
of
(x,
;
), denoted
(x,
;
;
). Specifically,
 | (11) |
where
k is a parameter vector drawn at random from the distribution N[µ(
),
], and Q is the total number of generated parameter vectors.
The same approach is used to predict grain protein content and residual mineral N at harvest.
Evaluation of Model Predictions
A natural criterion for evaluating accuracy of model predictions is the MSEP (Efron, 1983; Wallach and Goffinet, 1987). The smaller the value of MSEP, the smaller the prediction error. Efron (1983) compared cross validation and bootstrap methods for estimating MSEP. He found the bootstrap somewhat better, but the computing time required is large, so in this paper we use cross validation. For the yield predictions described above, the cross validation estimator of MSEP, denoted M
EPcvy, is
 | (12) |
where xij is the jth applied N rate for site-year i, yi(xij) is the corresponding measured yield value, ni is the number of N treatments in the ith site-year, I is the total number of site-years, and
-i and
-i are cross validation estimators of
and
. These estimators are obtained in the same way as
and
but after deleting the ith site-year from the data set. That is, for each site-year, we estimate the parameters from the other site-years. Then we use the resulting model to predict yield for each of the observed N doses.
The same approach can be followed to estimate MSEP for grain protein content or residual mineral N at harvest.
Calculation of Optimal Decision Rules
The random parameter model can be used to calculate optimal decision rules. For example, suppose that we have a model that predicts yield as a function of end-of-winter mineral soil N. This model can be used to estimate the expected value of the gross margin by
 | (13) |
where
(x,
) is the estimated value of the gross margin obtained for applied N dose x and site-year characteristic
, q is the grain price per unit of yield, and c is the price of one unit of N fertilizer. Maximizing
(x,
) for a particular value of
gives us an optimal N rate, XOPT(
). Maximizing
(x,
) for each
gives us a calculated optimal decision rule that specifies optimal N rates as a function of
. In the case of a model that contains no site-year characteristics, the optimal decision rule is simply reduced to a fixed N rate, the same for all site years.
With our model, it is not possible to obtain an analytical expression for XOPT(
), but it can be quite easily calculated numerically. To do so, one calculates
(x,
;
;
), using Eq. [11], for a series of x values. Then XOPT(
) is the value of x that maximizes
(x,
).
Evaluation of Decision Rules
An evaluation of the accuracy of model predictions does not provide direct information on the quality of calculated decision rules. Evaluating the quality of decision rules involves two problems. The first problem is to define an appropriate criterion for judging decision rules. The second problem is to estimate this criterion.
A natural criterion for evaluating the quality of a decision rule is the expectation of the objective function if the decision rule were implemented. This criterion is more directly related to the value of the calculated N doses than is MSEP. One way to estimate this criterion is to model response to applied N and then use that model to evaluate the results of a decision rule (Bullock and Bullock, 1994). This approach has two disadvantages. First, the accuracy of the criterion estimator depends on the quality of the parametric model used for modeling the response to applied N. Second, there is a risk of positive bias in the estimator of the criterion when both the decision rule and the estimator of the objective function are based on the same model.
Antoniadou and Wallach (2000) showed how the expected objective function value that results from applying a particular decision rule can be estimated using nonparametric regression. This approach does not involve a parametric model. The estimator of the objective function is estimated directly from data. In an extensive simulation study, the authors found that this method is at least as good and often much better than the method based on a parametric model. The bias in the nonparametric estimator is always very small. Briefly, the method involves four steps:
1. Calculate the value of the objective function for each applied N value for which data exist. For example, if the objective function is defined by Eq. [13], we calculate mi(xij) = q yi(xij) - c xij, i = 1, ..., I, j = 1, ..., ni, where mi(xij) is the gross margin value obtained for the jth level of applied N on the ith site-year using the corresponding measured yield value.
2. Derive the differences, XDIFFij = xij - XOPT(
i), i = 1, ..., I, j = 1, ..., ni, between the applied N values xij and the calculated optimal N rate for the ith site-year XOPT(
i).
3. Define a translated objective function m*i
= mi
.
4. Fit a nonparametric regression to the data to estimate E
. The resulting value, Ê
, is the estimate of the expected gross margin value when the calculated optimal decision rule XOPT(
i) is applied. Such a nonparametric regression can be performed, for instance, with the SPLUS function LOESS (S-PLUS4, 1997).
These four steps are used in the following sections to estimate the expected gross margin values that result from applying various optimal decision rules.
 |
MATERIALS AND METHODS
|
|---|
Data
The data used in this study were from 37 winter wheat experiments carried out between 1990 and 1996 on commercial farms in the Paris Basin of France (see Makowski et al. (1999) for a complete description). Each experiment was from a different site. Two soil types were represented, a loam soil and a chalky soil. Common winter wheat varieties were used. Each experiment consisted of five to eight different N fertilizer rates for a total of 224 N treatments. Nitrogen fertilizer was applied in two applications during the growing season. For each N treatment, N uptake (amount of N in grain, straw, and root at harvest), yield (adjusted to 150 g kg-1 grain moisture content), grain protein content, and residual mineral N in the soil at harvest
at 0- to 90-cm depth were measured. Range values for these measurements were N uptake, 73.9 to 426 kg ha-1; yield, 3.44 to 11.54 Mg ha-1; grain protein content, 6.15 to 15.9%; and residual mineral N at harvest, 11 to 211 kg ha-1, respectively. The variances of these four variables were N uptake, 5995.23 kg2 ha-2; yield, 2.90 Mg2 ha-2; grain protein content, 3.61%2; and residual mineral N at harvest, 771.36 kg2 ha-2, respectively. In each experiment, end-of-winter mineral soil N
in the 0- to 90-cm layer was measured during tillering (Feb.) before the N application. Values of end-of-winter mineral soil N were in the range 40 to 180 kg ha-1.
Models
Eight different models (Table 1) were studied. The first six are random parameter models based on the response functions proposed by Makowski et al. (1999) and taking into account end-of-winter mineral soil N as a site-year characteristic. The model LP-PL-RAND is defined by Eq. [1] through [10]. The five other models, namely QP-PL-RAND (QP, quadratic-plus-plateau), Q-PL-RAND, LP-PQ-RAND (PQ, plateau-plus-quadratic), QP-PQ-RAND, and Q-PQ-RAND, differ from LP-PL-RAND by the function describing the response of yield to N uptake, the function describing the response of residual mineral N to applied N, or both. In all six models, the linear-plus-hyperbola function defined by Eq. [1] and [2] was used for describing the response of N uptake to applied N, and the characteristic
was related to the parameter U0 as shown in Eq. [10].
In QP-PL-RAND and QP-PQ-RAND, a QP function was used for yield vs. N uptake. Combining this function with the linear-plus-hyperbola function for N uptake vs. applied N (Eq. [1] and [2]) gives the following function for yield vs. applied N:
 | (14) |
 | (15) |
and the following function for grain protein content vs. applied N:
 | (16) |
 | (17) |
In Q-PL-RAND and Q-PQ-RAND, a Q function was used for yield vs. N uptake. Combining this function with the linear-plus-hyperbola function for N uptake vs. applied N gives the following function for yield vs. applied N:
 | (18) |
 | (19) |
and the following function for grain protein content vs. applied N:
 | (20) |
 | (21) |
In the models LP-PQ-RAND, QP-PQ-RAND, and Q-PQ-RAND, the function used for residual mineral N at harvest vs. applied N was a PQ function defined by
 | (22) |
 | (23) |
As will be seen below, several of these functions are more or less equally good based on our evaluations. To test specific modeling options, we chose one set of functions giving good results, namely the set of functions used in LP-PL-RAND. This set of functions was used as a basis for the model LP-PL-RAND-0 (model with random parameters but without characteristic
) and LP-PL-FIXED (model without random parameters but with characteristic
). In LP-PL-FIXED, the model errors were supposed independent with different variances for the different types of response.
The estimation of all 55 elements of
was not numerically feasible for the random parameter models defined in Table 1. We reduced the number of elements in two steps. First we used the likelihood ratio test provided by the NLME routine (S-PLUS4, 1997) to identify elements of
that have variances not significantly different from 0 at a 5% level. This led us to treat T as a fixed rather than random parameter. The other nine variances were all significantly different from 0. In the second step, we drastically reduced the number of covariances to be estimated. Our approach here was to assume a covariance not equal to 0 only if there was some convincing argument in favor of that. This led us to estimate eight covariances (cov), namely cov(YMAX,B), cov(U0,C0), cov(U0,RMIN), cov(U0,R), cov(C0,RMIN), cov(C0,R), cov(RMIN,R), and cov(P1,P2). The first seven covariances were assumed not equal to 0 based on agronomic considerations: YMAX and B tend to be correlated because the climate has an effect on these two parameters (Boiffin et al., 1981). Soil properties and root development can influence both N uptake and residual mineral N at harvest, and consequently can induce correlations between U0, C0, RMIN, and R (Meynard et al., 1981; Wibawa, 1992). The correlation between P1 and P2 was found to be very high (-0.87), and so these two parameters were considered correlated.
and
were calculated for the different models using the NLME function of the SPLUS software (S-PLUS4, 1997). The parameters of the fixed parameter model, LP-PL-FIXED, were calculated using least squares for nonlinear models.
Optimal Decision Rules
The eight models were used to calculate optimal decision rules for three different objective functions. The first objective function, m1(x), is
 | (24) |
where 750 is the approximate current grain price in France (F Mg-1), and 2.5 is the fertilizer cost (F kg-1). The second expression for gross margin is
 | (25) |
where g[p(x)] is a grain price function that depends on grain protein content, p(x), as
 | (26) |
 | (27) |
 | (28) |
The function g[p(x)] is based on actual practice by a French co-operative, which reduces the price paid to farmers for grain with protein content <13%. The third expression for gross margin is
 | (29) |
where h[r(x)] is a penalization function that depends on the residual mineral N, r(x), as
 | (30) |
 | (31) |
Currently no such penalty exists in France, but the cost of F250 ha-1 seems reasonable because that is approximately the cost of planting a cover crop during winter to avoid excessive leaching of NO-3. The threshold of 35 kg ha-1 was chosen because that is approximately the level of residual mineral N that induces a water NO-3 concentration higher than the permissible NO-3 concentration of 50 mg L-1 for drinking water in Europe (EEC, 1980).
Each model in Table 1 was used to calculate optimal N rates for the gross margins in Eq. [24], Eq. [25], and Eq. [26] for each value of
in the data set. For the random parameter models, the expected values of the gross margins for each x value were estimated as the average over 2000 generated values of
. Then optimal N rates were calculated for each value of
from the estimated expected values of gross margin by testing x values in the range of 0 to 350 kg ha-1. The optimal N rates calculated with LP-PL-RAND-0 were the same for all site-years because this model does not include any site-year characteristics.
Evaluation
Each model in Table 1 was evaluated according to six criteria. Three of these criteria were MSEP values calculated by cross validation to evaluate the accuracy of predictions of yield
, grain protein content
, and residual mineral N
. For the random parameter models, predictions were estimated as in Eq. [11] with Q = 2000. For the model LP-PL-FIXED, predictions were derived by simply fixing the parameters at their estimated values.
The other three criteria were the expected values of the three versions of gross margin that would result from applying the optimal N rates calculated with the models. These criteria were estimated using a local regression as proposed by Antoniadou and Wallach (2000). These nonparametric regressions were performed for each model in Table 1 with a span fixed to 1 (i.e., 100% of the measurements were used to perform each regression) and with a tricube weight function (S-PLUS4, 1997). F-tests and graphical studies of residuals showed that a smaller span value was not justified.
 |
RESULTS AND DISCUSSION
|
|---|
Mean Squared Error of Prediction
The values of
and
are quite similar for all six different random parameter models that take into account end-of-winter mineral soil N (Table 2). This result indicates that the accuracy of predicting yield and grain protein content is not very dependent on the type of function used for describing the responses to applied N. On the other hand, the values of
obtained with the same six models are quite different. The highest value of
(worst prediction) is obtained with Q-PQ-RAND; the value obtained with QP-PQ-RAND is a bit lower; and the values of
obtained with LP-PL-RAND, QP-PL-RAND, Q-PL-RAND, and LP-PQ-RAND are much lower and more or less similar (Table 2). This may be due to the assumption that the N rate at which the residual mineral N begins to increase is equal to the N rate that maximizes yield, which may not be reasonable for the models QP-PQ-RAND and Q-PQ-RAND (Makowski et al., 1999).
Table 2 shows that the accuracy of predicting yield, grain protein content, and residual mineral N at harvest is improved by including end-of-winter mineral soil N to explain the variability in U0. This is shown by the fact that the values of
,
, and
obtained with LP-PL-RAND are 15, 20, and 19% lower, respectively, than the values obtained with LP-PL-RAND-0.
Finally, Table 2 shows that the accuracy of model predictions is similar for the random parameter model LP-PL-RAND and the least squares model LP-PL-FIXED. The MSEP values for yield and grain protein content obtained with LP-PL-RAND are slightly higher than those obtained with LP-PL-FIXED (by 0.7 and 4.7%, respectively). The MSEP value for residual mineral N at harvest obtained with LP-PL-FIXED is 6.5% higher than the value obtained with LP-PL-RAND. These marginal differences may be a result of the fact that our data are well balanced, having similar numbers of measurements for each site-year (Davidian and Giltinan, 1995). The advantage of a random parameter model would presumably be greater for highly unbalanced data.
Optimal N Rates
Table 2 gives the average, minimal, and maximal values of the optimal N rates calculated for the 37 site-years of the data set (and therefore 37 values of end-of-winter mineral soil N) using the eight different models of Table 2. The range in optimal N rate with end-of-winter mineral soil N is large (e.g., 90225 kg ha-1 for gross margin m1 using LP-PL-RAND). Figure 1 shows that optimal N rate calculated with model LP-PL-RAND decreases roughly linearly as a function of end-of-winter mineral soil N. A similar result was obtained with the other models that take this characteristic into account (not shown). This relationship is logical. Increased end-of-winter mineral soil N means that the soil is providing more N, decreasing the requirement for applied N. A linear relationship between optimal fertilizer rates and soil test was also found by Mombiela et al. (1981).

View larger version (18K):
[in this window]
[in a new window]
|
Fig. 1. Optimal N rates as a function of end-of-winter mineral soil N calculated with model LP-PL-RAND for three versions of gross margin (m1, m2, and m3). Grain price is fixed, and high levels of residual mineral N at harvest are not penalized in m1. In m2, grain price depends on grain protein content, and in m3, high levels of residual mineral N at harvest are penalized
|
|
Table 2 shows large differences in optimal N rates between the six random parameter models including end-of-winter mineral soil N. For instance, the average value of the optimal N rates calculated for the gross margin m1 varies from 158 to 186 kg ha-1 among these models. Important differences between these models are observed for gross margins m2 and m3 as well. These results show that the optimal N rate depends strongly on the functions used for describing the responses to applied N. Such discrepancies between response functions have been already reported in previous studies (Cerrato and Blackmer, 1990; Bullock and Bullock, 1994).
With all the random parameter models, the optimal N rates calculated for gross margin m2 are higher than those obtained for gross margin m1 while the optimal N rates for gross margin m3 are lower (Table 2 and Fig. 1). This is logical. Gross margin m2 favors high protein content, and thus high applied N, whereas m3 favors low residual N, and thus low applied N.
Table 2 also shows that the optimal N rates calculated with model LP-PL-FIXED vary only marginally with gross margin type and are much lower than the optimal N rates calculated using LP-PL-RAND. The average values of the optimal N rates calculated using LP-PL-FIXED for m1, m2, m3 are only equal to 89, 90, and 89 kg ha-1, respectively, whereas these rates are equal to 186, 213, and 178 kg ha-1, respectively, when LP-PL-RAND is used (Table 2). Thus, the type of statistical model used for modeling response to applied N can have a strong influence on the values of the calculated optimal N rates. This result is consistent with the result of Babcock (1992), who showed analytically for a simple LP yield response function that a random plateau (i.e., random maximal yield) induces an increase in the calculated optimal N rates if the price of N fertilizer is low relative to the price of the crop. Similar results are obtained here with a more complex model including several random parameters.
The drastic differences observed between the optimal N rates calculated with LP-PL-RAND and LP-PL-FIXED are due to the influence of the random parameters on the shape of the predicted gross margin response curves as illustrated in Fig. 2 for a site-year with end-of-winter mineral soil N equal to 80 kg ha-1. Figure 2a shows that gross margin m1 vs. applied N has a discontinuous derivative when model LP-PL-FIXED is used but not when LP-PL-RAND is used. This result is due to the fact that, with LP-PL-RAND, the gross margin response curve is calculated by averaging 2000 LP response curves (Eq. [11]), resulting in a smooth average response curve. Similar results were obtained analytically by Berck and Helfand (1990) for a simple response function. Differences between LP-PL-RAND and LP-PL-FIXED are important for other versions of gross margin as well (Fig. 2b and 2c).

View larger version (21K):
[in this window]
[in a new window]
|
Fig. 2. Responses of gross margins (a) m1, (b) m2, and (c) m3 to applied N predicted for a site-year with end-of-winter mineral soil N equal to 80 kg ha-1 using models LP-PL-FIXED and LP-PL-RAND. Grain price is fixed and high levels of residual mineral N at harvest are not penalized in m1. In m2, grain price depends on grain protein content, and in m3, high levels of residual mineral N at harvest are penalized
|
|
Profitability of Optimal N Rates
Table 2 gives the nonparametric estimates of the gross margin values obtained when the calculated optimal N rates are applied. The results discussed above show important differences between optimal N rate values. Here we compare the monetary results induced by these different N rates to determine which model is better for deriving N fertilizer recommendations. It is seen that the gross margin values obtained with N recommendations based on the six random parameter models taking into account end-of-winter mineral soil N are quite similar. The differences between the minimal and maximal gross margin values among the six models are only equal to F39, F17, and F44 ha-1 for m1, m2, and m3, respectively (Table 2). Thus, the optimal N rates calculated using the different functions, though quite different, give very similar values of profitability.
Compared with the values obtained by applying the optimal N rates calculated with LP-PL-RAND, the gross margin values obtained by applying the optimal N rates of LP-PL-RAND-0 are lower by F102, F17, and F105 ha-1 for m1, m2, and m3 (Table 2). Thus, the optimal N rates calculated by using the model taking into account end-of-winter mineral soil N are more profitable. However, it should be noted that the monetary gain is relatively small compared with gross margin values. Moreover, this gain will be reduced by the cost of obtaining measurements of end-of-winter mineral soil N, which is about F90 for one field. As field sizes in the Paris Basin are usually in the range 5 to 10 ha, this cost is in the range of F9 to F18 ha-1.
The small values of optimal N calculated using LP-PL-FIXED lead to very low gross margin values (Table 2). Thus, compared with the values obtained with the optimal N rates calculated using LP-PL-RAND, the gross margin values obtained by applying the optimal rates calculated using LP-PL-FIXED are lower by F477, F550, and F438 ha-1 for m1, m2, and m3, respectively. These results show that the decision rules calculated using a random parameter model can be more profitable.
 |
SUMMARY AND CONCLUSIONS
|
|---|
This paper shows how the functions of Makowski et al. (1999) that describe yield, grain protein content, and residual mineral N responses to applied N can be used for predicting responses to applied N and calculating optimal N rates. This is done by treating these functions as random parameter models. The proposed method employs standard statistical software and can be applied with response functions other than those considered in this paper.
Two types of model evaluation are used here. The first estimates the predictive accuracy of a model in terms of MSEP. The second type of evaluation measures the quality of the optimal N rates calculated using the model. The criterion used in the second type of evaluation is the value of the objective function that would be achieved if the calculated optimal N rates were applied. Both MSEP and the value of the calculated optimal strategy are estimated directly from the data without making any assumption that the tested model, or some other parametric model, is correct. Thus, they provide objective tools for comparing different models and different modeling strategies. These criteria and the associated estimators have been taken from the literature and adapted to this study.
The argument in favor of random parameter models is that they provide a realistic description of the structure of the data. In the case study treated here, it is further shown that the profitability of optimal N rates calculated using a random parameter model can be substantially higher (by F438F550 ha-1 with the set of functions considered in the application) than the profitability of optimal N rates derived from a simpler statistical model. On the other hand, the use of a random parameter model does not improve the accuracy of model predictions in our example.
An advantage of the functions of Makowski et al. (1999) is that site-year characteristics can conveniently be introduced into the model. The model parameters are meaningful, which allows one to use agronomic knowledge to decide how to introduce site-year characteristics. We showed in this paper how the model parameter representing N uptake when no fertilizer is applied can be related to end-of-winter mineral soil N measurements. This site-year characteristic was chosen because it is often available and seems likely to be one of the most important characteristics to introduce. In the case study, a model taking into account end-of-winter mineral soil N was compared with a model based on the same functional forms but without a site-year characteristic. It was found that the MSEP values of the model without a site-year characteristic are about 20% higher and that the gross margin values obtained by applying the optimal N rates calculated with this model are F17 to F105 ha-1 lower.
Other field characteristics such as the nature of the buried crop residues, soil organic matter, or soil texture could be easily introduced in our model. Including all of these characteristics may improve model performance. However, such improvement would not be systematic because including new site-year characteristics means increasing the number of model parameters to be estimated, which may degrade model performance.
In the case study, the different functional forms proposed by Makowski et al. (1999) were compared with end-of-winter mineral soil N in the model. Four of the models, namely LP-PL-RAND, QP-PL-RAND, Q-PL-RAND, and LP-PQ-RAND, all gave similar results. Based on these results, there does not seem to be any convincing reason, in terms of predictive accuracy or gross margin value, to choose one over another.
Numerical results are, of course, specific to the data set considered. However, the general conclusions (advantage of the random parameter models, usefulness of end-of-winter mineral soil N, choice of response functions) may serve as indicators for other studies.
 |
REFERENCES
|
|---|
- Anderson, R.L., and L.A. Nelson. 1975. A family of models involving intersecting straight lines and concomitant experimental designs useful in evaluating response to fertilizer nutrients. Biometrics 31: 303318.[Web of Science]
- Antoniadou, T., and D. Wallach. 2000. Evaluating decision rules for nitrogen fertilization. Biometrics 56:194199.
- Babcock, B.A. 1992. The effects of uncertainty on optimal nitrogen applications. Rev. Agric. Econ. 14:271280.
- Berck, P., and G. Helfand. 1990. Reconciling the von Liebig and differentiable crop production functions. Am. J. Agric. Econ. 72: 985996.
- Boiffin, J., J. Caneill, J-M. Meynard, and M. Sebillotte. 1981. Elaboration du rendement et fertilisation azotée du blé d'hiver en Champagne crayeuse: I. Protocole et méthode d'étude d'un problème technique régional. Agronomie (Paris) 1:549558.
- Bullock, D.G., and D.S. Bullock. 1994. Quadratic and quadratic-plus-plateau models for predicting optimal nitrogen rate of corn: A comparison. Agron. J. 86:191195.[Abstract/Free Full Text]
- Cerrato, M.E., and J-M. Blackmer. 1990. Comparison of models for describing corn yield response to nitrogen fertilizer. Agron. J. 82:138143.[Abstract/Free Full Text]
- Chaney, K. 1990. Effect of nitrogen fertilizer rate on soil nitrate nitrogen content after harvesting winter wheat. J. Agric. Sci. (Cambridge) 75:533537.
- Colwell, J.D. 1994. Estimating fertilizer requirements: A quantitative approach. CAB Int., Wallingford, UK.
- Davidian, M., and D.M. Giltinan. 1995. Nonlinear models for repeated measurement data. Chapman and Hall, London, UK.
- Efron, B. 1983. Estimating the error rate of a prediction rule: Improvement on cross-validation. J. Am. Stat. Assoc. 78:316331.[Web of Science]
- [EEC] European Economic Community. 1980. Off. J. EEC 23:1129.
- Fowler, D.B., J. Bryfon, and R.J. Baker. 1989. Nitrogen fertilization of no-till winter wheat and rye: II. Influence on grain protein. Agron. J. 81:7277.[Abstract/Free Full Text]
- Jauregui, M.A., and Q. Paris. 1985. Spline response functions for direct and carry-over effects involving single nutrient. Soil Sci. Soc. Am. J. 49:140145.[Abstract/Free Full Text]
- Lindstrom, M.J., and D.M. Bates. 1990. Nonlinear mixed effects models for repeated measures data. Biometrics 46:673687.[Web of Science][Medline]
- Machet, J-M., and B. Mary. 1990. Effet de différentes successions culturales sur les risques de pertes en nitrate en région de grande culture. p. 395405. In R. Calvet (ed.) Nitratesagricultureeau. Institut National de la Recherche Agronomique (INRA), Paris, France.
- Makowski, D., D. Wallach, and J-M. Meynard. 1999. Models of yield, grain protein, and residual mineral nitrogen responses to applied nitrogen for winter wheat. Agron. J. 91:377385.[Abstract/Free Full Text]
- Meynard, J-M., J. Boiffin, J. Caneill, and M. Sebillotte. 1981. Elaboration du rendement et fertilisation azotée du blé d'hiver en Champagne crayeuse: II. Types de réponse à la fumure azotée et application de la méthode du bilan prévisionnel. Agronomie (Paris) 1:795806.
- Meynard, J-M, E. Justes, J-M. Machet, and S. Recous. 1997. Fertilisation azotée des cultures annuelles de plein champ. p.183199. In G. Lemaire and B. Nicolardot (ed.) Gestion de l'azote dans les agrosystèmes. Institut National de la Recherche Agronomique (INRA), Paris, France.
- Mombiela, F., J.J. Nicholaides III, and L.A. Nelson. 1981. A method to determine the appropriate mathematical form for incorporating soil test levels in fertilizer response models for recommendation purposes. Agron. J. 73:937941.[Abstract/Free Full Text]
- Murray, A.W.A, and P.A. Nunn. 1987. A nonlinear function to describe the response of % nitrogen in grain to applied nitrogen fertilizer. Aspects of Appl. Biol. 15:219225.
- Neeteson, J.J. 1990. Development of N fertilizer recommendations for arable crops in the Netherlands in relation to nitrate leaching. Fert. Res. 26:291298.
- Rémy, J.C., and J. Hébert. 1977. Le devenir des engrais azoté dans le sol. C. R. Acad. Agric. Fr. 63:700710.
- Sain, G.E., and M.A. Jauregui. 1993. Deriving fertilizer recommendations with a flexible functional form. Agron. J. 85:934937.[Abstract/Free Full Text]
- S-PLUS4. 1997. Guide to statistics. MathSoft, Seattle, WA.
- Stanford, G. 1973. Rationale for optimum nitrogen fertilization in corn production. J. Environ. Qual. 2:159166.[Abstract/Free Full Text]
- Wallach, D. 1995a. Regional optimization of fertilization using a hierarchical linear model. Biometrics 51:338346.
- Wallach, D. 1995b. Empirical Bayes optimal fertilizer decisions. J. Appl. Stat. 22:507516.
- Wallach, D., and B. Goffinet. 1987. Mean squared error of prediction in models for studying ecological and agronomic systems. Biometrics. 43:561573.[Web of Science]
- Wibawa, G. 1992. Approche, par enquête et expérimentation, de l'effet de l'état structural du sol sur la nutrition azotée et l'élaboration du rendement de l'orge de brasserie. Ph.D. diss. Institut National Agronomique Paris-Grignon (INA P-G), Paris.