Long-term fertilizer field trials : comparison of three mathematical response models

Accession to the European Union caused a drop of nearly 60 per cent from 1994 to 1995 in prices of wheat, barley and oats in Finland. The economic use of fertilizer therefore decreased accordingly. To calculate the effect of the price changes on the economic optima, the physical production function must be known. Three physical production functions, the quadratic, the linear response and plateau (LRP) and the exponential function were estimated for this purpose. The models differed little in respect of the R >dj value (0.82-0.90) but the calculated optimum varied, depending on the production function. Data on a long-term field trial (21 years) were analysed. The field trial was established in 1973 to demonstrate the effect of mineral fertilizer in crop production. The crops grown in the trial were barley, wheat and oats. Different varieties were included in the models.

ntroduction Fertilizer field trials yield a vast amount of data that cannot be analysed without the use of a proper mathematical model. Fertilizer recommendations and political decisions on fertilizer use depend directly on the parameters estimated from the data. However, in agronomy there is no standard system for choosing between models to represent the relationship between fertilizer input and yield formation (Cerrato and Blackmer 1990).
All the models chosen, the quadratic model, the exponential model usually known as the Mitscherlich model, and a modified form of the plateau model, are commonly used in crop response analyses (Bock and Sikora 1990, Cerrato and Blackmer 1990, Frank et al. 1990, Paris 1992, Sumelius 1993. Even though the models give comparable R 2 values, they may give different optimal fertilizer rates. Optima based on the quadratic function or on the exponential function have been criticised for giving excessively high optimal fertilizer rates. A model based on a linear response and plateau (LRP) function gives lower optima (Paris 1992). Frank et al. (1990) found that the costs of wrong specification of optimal input amounts are substantial. They claimed that neither polynomial nor plateau growth should be assumed a priori on agronomic grounds. Paris (1992) pointed out that most people interpret the von Liebig hypothesis as a linear relation to limiting nutrient, and that this interpretation is not consistent with von Liebig's concept. According to Paris, the von Liebig hypothesis claims that there is a direct, but not necessarily linear, relation to the limiting nutrient. Incorrect specification of optimal input has not only economic but also environmental impacts. Excessive application rates lead to nutrient losses to the environment by leaching or volatilization. Miettinen (1993) studied the effect of fertilizer policy and leaching and found that imposition of input quotas caused farmers the lowest additional costs in barley production. The loss of profits was greater with other measures and the decrease in nitrogen leaching was smaller. Evaluation of these input quotas (or nutrient quotas) again requires the use of mathematical models.
We here compare three models (quadratic, LRP, exponential) in an effort to establish their benefits and drawbacks and to calculate optima according to them. Some of the potential benefits and drawbacks are assumed in advance (Table 1). The use of these models in practice and the theory behind them are also discussed. To test the models we used data on a long-term field trial (21 years) in which compound NPK fertilizers were used. We found that nitrogen explained most of the yield response; fertilization with phosphorus did not increase the yield, though at the lowest levels the yield was determined by both nitrogen and phosphorus (Yli-Halla 1991).

Material and methods
The fertilizer field trial was located on the same plot, with five fertilizer levels being applied each year. The trial was established in 1973 and the data used are from the period 1973-1993. New varieties replaced old ones but the cultivation techniques remained the same throughout the period. The field trial was conducted at Kemira Agro's research farm at Kotkaniemi, 40 km northeast of Helsinki. The change in varieties was gradual. Barley (29 varieties), wheat (21 varieties) and oats (16 varieties) were cultivated. The five fertilizer levels are referred to as A, B, C, D and E in increasing order of fertilizer application, with level A as zero. The primary purpose of the trial was demonstrative, and the treatments were not randomized.
The amounts of fertilizer applied in the trial were reduced from 1986 onwards ( Table 2), except at level A, at which it was zero throughout. The NPK fertilizer applied was the compound fertilizer most commonly used in Finland. The composition and amounts of this fertilizer are listed in Table 2.
The nutrients N, P and K are considered as linear combinations and N is assumed to be the main limiting nutrient. Partial correlation anal-Vol. 6(1997): 151-160. Year 1973-19791980-19851986-19881989  yses also show a higher value for the N fertilizer than for the compound fertilizer as a whole (Table 3). The N fertilizer applied was chosen as the independent variable in the model, but note that P and K were also applied as macronutrients in the field trial. Differences in the composition of the compound fertilizer had little effect on output (Fig. 1). Annual effects were analysed with F tests. Inclusion of annual effects allowed for annual physical production functions. The effect of residual N on levels B, C, D and E was assumed to be nil from year to year. The effect of residual P was assumed to be close to nil at level C. At level A, yield decreased from 1973 onwards.
The functions are quadratic, LRP and exponential. The LRP and exponential functions were chosen because they reflect a level at which no further increase in yield response is achieved. The LRP function showed a "breakpoint" according to von Liebig's Law of the Minimum. The exponential function increases asymptotically to a maximal yield level.
The indices i and j denote variety and year variationsrespectively, y stands for the yield response in kg/ha and x for the N fertilizer applied in kg/ha. a, b and c are parameters and e are error terms that are assumed to be normally distributed and with zero mean. This model, in which all coefficients vary both among varieties and in time, is easy to calculate and is the one most commonly used in fertilizer trial analyses.
The indices i and j denote variety and year variations respectively, y stands for the yield re- 153 AGRICULTURAL AND FOOD SCIENCE IN FINLAND sponse in kg/ha and x for N fertilizer applied in kg/ha. The index b on x stands for the level of fertilizer that achieves the breakpoint in the yield response to additional fertilizer nitrogen, y is the yield response level achieved in a given year for a specific variety. The function imposes a maximal yield level at which further application of a specific input will not give any further increase in the yield.
The indices i and j denote variety and year variations, respectively, y stands for the yield response in kg/ha and x for the N fertilizer applied in kg/ha. A factor, 0.01, was included in the model because of a problem with decimals when printing the results, and this factor was thereafter considered in the results. The function gives an asymptotic maximal yield level. The growth to the maximal level diminishes exponentially, which means that the parameter c. is negative.
The quadratic function was computed with IMSL subroutines; RGLM for fitting the multivariate linear regression model and RSTAT for computing and printing statistics. These subroutines were retrieved from the MATLAB program. The quadratic function could thereafter be estimated simply by changing the dummy variable effects. For estimation of the LRP and the exponential functions we used the subroutine DUNLSF. This IMSL command solves a nonlinear least squares problem using a modified Levenberg-Marquardt algorithm and a finite-difference Jacobian. The Levenberg-Marquardt method is a modification of the Gauss-Newton algorithm for solving non-linear least squares problems (Appendix 1). From one current point another point is calculated by the trust region approach (Dennis and Schnabel 1983). This procedure is repeated until stopping criteria are satisfied. Each separate estimation was made by a different FORTRAN program. All effects were analysed with F tests.
For estimation of the LRP model, a program loop was used to restrict the movement of values of the parameter x b| , with values outside the domain of the field trial  giving a penalty to the minimization of the non-linear least squares. The first restricting penalty loop (algorithm x b = (200-x b ) 2 ) we tried did not perform effectively enough on all data. The second restricting penalty loop, x h = xb (Fig. 2), was more effective and x. finally varied between -4 bij J and 4, which meant optimal solutions between 0 and 200 kg of N per hectare. These penalty loops served to bring the routines to a clear stop. The first algorithm probably did not work well because the sample range was too small to reflect the decreasing response phase in all cases. The barley data were nevertheless calculated with the first algorithm, and so barley reflects a higher plateau response than wheat or oats (Fig. 3).
We assumed that all species and varieties in the trial would respond differently to fertilizer, and that all species and varieties in the trial would give different yields without fertilizer. We also assumed that every year would give different yield responses. F tests were performed to test our assumptions.
The R a 2 d value will give the estimated function response to the distribution pattern, is the R a 2 d value corrected for the number of degrees of freedom and is therefore better suited for comparing different models with different numbers ofparameters. The number of parameters is also a criterion in which the number of parameters as well as the flexibility of the function should be related to the information they give. Preference should be given to function estimates that are easy both to explain and to calculate.

Economic analysis
Making marginal cost equal to marginal revenue gives us the maximal economic return. This occurs, according to the production theory, when the marginal physical product (MPP) is equal to the price ratio, w/p, where w is the unit price of fertilizer and p the unit selling product price. Due to the short growing season in Finland, annual variations in the corresponding yield responses are high. Knowledge of the physical production function is therefore valuable. Another factor adding to the interest of the physical production function interesting is the change in price relations in Finland caused by accession to the European Union in 1995. The ensuing switch from an agricultural policy favouring price support to one favoring more direct support meant a drop of about 60 per cent in the price of cereals. The maximal yield is the point at which the slope of the production function equals zero. The fertilizer input for maximal economic output, when costs are included, is lower than that for maximal yield. The maximal economic return calculated with the LRP function is approximated to equal zero, or breakpoint. The first-order derivatives of the exponential function show four solutions, depending on the signs selected for the parameters. We chose the solution that is most likely to occur on the basis of diminishing returns (Table 4). The optimal economic output is calculated from the mean LS estimate, which gives a lower result than the average of the sum of the optima calculated with the included effects.

Results
The values estimated for the quadratic function were checked with F tests. All the effects were significant except for the variations in the response of the oats variety to N. The second degree variety variations were not significant either. Estimation of the quadratic function succeeded without major drawbacks but resulted in many optima that were explicit (outside the sample range).
The parameters estimated (Table 5) are the values that can be fitted to the models presented. Note that the intercepts (parameter a, LRP and quadratic model) for oats and barley are higher than those for wheat. Wheat also responds less to N fertilizer than does oats or barley (Fig.  3). The restricting penalty algorithm used for barley did not allow us to include as many effects as did the algorithm used for wheat and oats. This is probably why the breakpoint for barley is higher than that for oats and wheat. The parameters a for the exponential model show (Table 5) the asymptotic yield level. The standard error in the parameter of the asymptotic level for wheat is remarkably high, possibly implying that the annual and variety variations do not explain the entire yield variation for wheat at high fertilizer levels.
A comparison of the distribution pattern and the LRP model estimated for oats shows that a reduction in the yield response for oats at high fertilizer levels cannot be reflected by the LRP model. Oats shows a clear decrease in yield response at high fertilizer rates, but the LRP model does not have a decreasing phase at which MPPcO. The residuals now cross the zero line twice and the optimum should therefore be found between the crossings. The yearly variation in the breakpoint at the x-axis was not significant for the LRP model, but the resulting yield plateau or N response had significant variations. There were too few replications for each variety in the data on the barley experiments. The LRP and the exponential results for barley were therefore calculated without considering variations due to variety.
Estimation of the exponential function succeeded well. The R 2 values were higher than those for the quadratic function. The exponential function does not have a phase for a decreasing yield response and so is it less suitable for yield response analysis than the quadratic function. Only a few of the optima obtained were implicit. Moreover, the sample range was too Vol. 6(1997): 151-160. small when the exponential function was used on the data for wheat. The number of parameters varies according to the number of effects included (Table 6). When only R * was compared, the quadratic function performed better than the exponential or LRP function on oats and barley. The exponential function performed best on wheat, but there was a problem with the sample range, which may explain the high error term for the coefficient a for wheat. The LRP function performed only slightly worse than the others. R 2 , , however, is not a adj very good criterion for selecting a model, as shown by, among others, Kvålseth (1985), McGuirk and Driscoll (1995). It would be more appropriate to test the models against each other. This would be an interesting topic for a follow-up study. We did not do so here because our main interest was in finding out which effects to include, and also because we lacked reliable methods for testing the models against each other.
It is striking that, according to the quadratic and exponential models, a decrease in optimal fertilizer use was a result oflower prices in 1995 (Table 7). According to the LRP model, however, changes in price do not affect economic optimal fertilizer use. This is because the MPP is zero after the breakpoint, which tells us that any additional fertilizer input gives no further increase in yield. Optimal fertilizer use according to the LRP model is, however, still lower than according to the quadratic or exponential model. Newer varieties show an increase in optima.

Discussion
Recommendations based on an "incorrect" model have an impact on both profitability and on the environment. Preference has been given to functions that enable us to relate the parameters achieved to biological and physical processes. Biological theories and economic theories are not, however, always comparable. In agronomy, there is no satisfactory relation between growth inputs and yield formation. Attempts have been made to construct deterministic models, but as van Keulen and Stol (1991) remarked: "These models cannot be used in practice. They are an aid to structure thinking about the system." Stochastic models are more often used, but they may lack any idea with which the results can be explained.
The optimization of fertilizer use depends on the aim ofa particular optimum. Possible optima are: a) biological (includes quality and quantity) b) economic (includes micro-and macroeconomic optima) c) environmental.
The biological optimum depends on the quantity measures and will give the biologically maximal yield. The biologically optimal input of fertilizers may be different if quality is taken into account. Economic optima are calculated to give the maximal economic return for the individual farmer at the microlevel. The optimal macroinput of fertilizer is the optimal input for a certain number of farmers with different production options (functions). The results of these experimental data cannot be used as such at the macrolevel. Here we have concentrated on deriving the economic and biological quantity optima for this experiment, and considered only the variable fertilizer costs.
The quadratic function has been criticized for giving excessively high estimates (Ackello- Ogutu et al. 1985, Paris 1992. We, however, found that it was easy to both calculate and use for optima. The quadratic function has a diminishing phase that was accurate if comparisons are made with the distribution pattern of the data on oats. The high estimates are partly due to the inflexibility of the quadratic function but partly to the size of the domain. The lack of the diminishing phase in the exponential function is a disadvantage. A diminishing phase would, however, require inclusion of a fourth parameter in the model. The advantage of a four-parameter exponential function over the quadratic function would be that the former could include a steeper increasing phase and a flatter diminishing phase. The disadvantage of a fourth parameter is the decrease in the degree of freedom. The quadratic function is symmetric in the increasing and diminishing phases. The distribution pattern appears not to be symmetric. The diminishing phase is interesting only ifit affects the earlier (increasing and rational) phases; otherwise it has no significance for the optima. A non-linearfunction with three parameters that could include a steeper increasing phase and a flatter diminishing phase would probably perform well. These non-linear functions are available but they do not necessarily have any relevance to biological theories. These Vol. 6(1997): 151-160.
functions might be an interesting topic for a follow-up study.
The LRP model gives optimal input levels that are closer to the recommendations used today than the optimal input levels given by the exponential or quadratic models. According to the exponential function, the optimal use of fertilizer is highly dependent on the unit price relation between input and output. This is not so according to the LRP function. The LRP function is, however, quite cumbersome to calculate and therefore probably not that useful.
If the error term were brought further into the calculation of the optimal solutions in the form of an economic risk, the optimizations obtained would be more accurate. This was done for the quadratic function by Carmer et al. (1991). It is, however, a very time-consuming task and the results should be considered in the perspective of variability in actual growing conditions. Optimising under varying conditions can give us higher optima because of the price relation and attitudes to risk. The average curve is still, however, the expectation curve.
In conclusion, we stress that none of these three models, quadratic, LRP or exponential, performed satisfactorily in all respects. The quadratic model performed best for oats because of the diminishing phase. The LRP performed best in view of the recommendations, and the exponential function performed best in respect of the steep increasing phase. As a topic for further research new models should be tried out or new methods introduced that include other factors affecting economic output.