Introduction
Optimisation of fertiliser use is increasingly becoming critical due to the growing need to balance agricultural production with demand for nutritious food and global environmental concerns (Giller 2020; Jones et al., Reference Jones, Cross, Withers, DeLuca, Robinson, Quilliam, Harris, Chadwick and Edwards-Jones2013; Palm et al., Reference Palm, Neill, Lefebvre and Tully2017). This is particularly important in sub-Saharan Africa (SSA), where fertiliser application rates have been historically low, nutrient depletion rates are high, and deficiencies of macronutrients and micronutrients are now widespread (Kihara et al., Reference Kihara, Sileshi, Nziguheba, Kinyua, Zingore and Sommer2017, Reference Kihara, Bolo, Kinyua, Rurinda and Piikki2020). Fertiliser application rates often involve blanket recommendations of nitrogen (N), phosphorus (P) and to a lesser extent of potassium (K) for a whole region (Ichami et al., Reference Ichami, Shepherd, Sila, Stoorvogel and Hoffland2018; Snapp et al., Reference Snapp, Blackie and Donovan2003). This gradually exhausted soil nutrient reserves through nutrient mining, leading to imbalances and deficiencies of micronutrients and consequently poor crop response to NPK fertiliser (Ichami et al., Reference Ichami, Shepherd, Sila, Stoorvogel and Hoffland2018; Kihara et al., Reference Kihara, Sileshi, Nziguheba, Kinyua, Zingore and Sommer2017, Reference Kihara, Bolo, Kinyua, Rurinda and Piikki2020; Tamene et al., Reference Tamene, Sileshi, Ndengu, Mponela, Kihara, Sila and Tondoh2019).
Governments, international agricultural research institutions, and donors continue to invest in fertiliser trials across SSA (Kaizzi et al., Reference Kaizzi, Mohammed, Nouri, Wortmann and Sones2016; Snapp et al., Reference Snapp, Blackie and Donovan2003). Recent efforts in a number of African countries have also led to the development of economically optimum nutrient rates (EORs) and fertiliser optimisation tools (Kaizzi et al., Reference Kaizzi, Mohammed, Nouri, Wortmann and Sones2016). All these efforts recognise that site-specific nutrient management needs to be guided by robust nutrient dose-response models. However, practitioners sometimes apply response functions that are not suitable for the data, resulting in EORs that are at odds with existing N and P recommendations. For example, the EOR was recently reported to be 6 kg N ha−1 against the recommended rate of 75 kg N ha−1 for a maize crop with an asymptotic yield of 3.7 t ha−1 in the Kenyan Rift Valley (Supplementary Table S1). In the same area, EOR was estimated at 0 and 8 kg P ha−1 for a maize crop with an asymptotic yield of 4.7 and 6.1 t ha−1, respectively (Supplementary Table S1). Similarly, the EOR was reported to be 27 kg N ha−1 and 0 kg P ha−1 for a maize crop with an asymptotic yield of 4.7 and 3.4 t ha−1, respectively, in the Central zone of Tanzania (Supplementary Table S1). Close examination of these reports reveals that very low EOR values are artefacts of errors in parameter estimation (Supplementary Table S1). Note that a 5 t ha−1 maize grain harvest requires the uptake of approximately 100 kg N ha−1, 24 kg P ha−1 and 85 kg K ha−1 (Nalivata et al., Reference Nalivata, Kibunja, Mutegi, Tetteh, Tarfa, Dicko, Ouattara, Cyamweshi, Nouri, Bayu, Wortmann, Wortmann and Sones2016). Against that background, all of the EORs in Supplementary Table S1 may be too low to replace what is removed through crop harvest and therefore speed up nutrient mining.
Derivation of EORs depends on accurate estimation of the optimum nutrient rate and the maximum yield response, which often vary with the cultivar used, soil type and climate. However, the following principles hold good in most situations (Soffe, Reference Soffe2011): (1) some yield (Y 0 ) is always obtained where no external nutrient has been applied, depending on the level of indigenous nutrients in the soil; (2) there is a finite maximum agronomic yield (Y max ) which is obtained from a certain level of nutrient input; and (3) Y max is reached at a point where the agronomic use efficiency of the nutrient is so low that further application would not result in yield increase. Therefore, the EOR is found below the nutrient rate that gives Y max depending on the ratio of crop value to nutrient cost (Soffe, Reference Soffe2011). The EORs can also vary widely across a set of response functions applied to the same dataset. The uncertainty associated with the estimation of EORs is, however, often overlooked (Hernandez and Mulla, Reference Hernandez and Mulla2008; Morris et al., Reference Morris, Murrell, Beegle, Camberato, Ferguson, Grove, Ketterings, Kyveryga, Laboski, McGrath, Meisinger, Melkonian, Moebius-Clune, Nafziger, Osmond, Sawyer, Scharf, Smith, Spargo, van Es and Yang2018). The assumptions behind the different response functions and the disagreements among their predictions have also received little attention (Cerrato and Blackmer, Reference Cerrato and Blackmer1990).
Crop response can significantly vary from site to site or season to season within a site (Albarenque et al., Reference Albarenque, Basso, Caviglia and Melchiori2016; Xu et al., Reference Xu, He, Pampolino, Qiu, Zhao and Zhou2019), and striking examples are summarised in Figure 1. Therefore, researchers face tremendous challenges in selecting a mathematical function sufficiently general to avoid the need to develop separate models for each cropping season or soil type. Most response functions require fitting nonlinear regression models whose solution space can be discerned via trial and error (Archontoulis and Miguez, Reference Archontoulis and Miguez2015). As such, it is often difficult to know whether or not the estimated parameters represent the underlying reality (ibid.; Bachmaier and Gandorfer, Reference Bachmaier and Gandorfer2012). This is rarely appreciated by practitioners, and sometimes arbitrarily selected models have been applied with over-confidence for development of EORs. Therefore, the objectives of this paper were to (1) point out common problems in the development of dose-response models and (2) provide corrective measures to facilitate future trial design and data analysis.
Figure 1. Variability in maize grain yield with site and season in Uganda and (data from Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Musinguzi, Ebanyat, Hyuha and Wortmann2012) and Tanzania (data from Senkoro et al., Reference Senkoro, Marandu, Ley and Wortmann2017). The red lines in each figure represent yields that have not yet reached Y max . The blue and black lines represent yields that are declining or remaining constant after reaching their maximum, respectively.
Methods
In 2016–2018, the author of this article facilitated a series of training workshops on the analysis of data from trials for optimising fertiliser recommendations in a number of projects across southern, eastern and western Africa. During these workshops, a number of critical problems were identified when analysing nutrient-response data. Motivated by this first-hand experience, the author examined a number of legacy data on fertiliser response and undertook a critical review of the literature on the subjects published between 2000 and 2020 across SSA. The data used include maize yield response to N in Niger (from Maman et al., Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018), Malawi and Rwanda (from Wortmann et al., Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018), Tanzania (from Senkoro et al., Reference Senkoro, Marandu, Ley and Wortmann2017) and Uganda (from Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Musinguzi, Ebanyat, Hyuha and Wortmann2012); wheat yield response to N in Tanzania and Rwanda (from Cyamweshi et al. Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018); barley yield response to N (from Agegnehu et al., Reference Agegnehu, Nelson and Bird2016) and P (from Abdulkadir et al. 2017) in Ethiopia; and rice yield response to N, P and Zn in Nigeria (from Daudu et al., Reference Daudu, Ugbaje, Oyinlola, Tarfa, Yakubu, Amapu and Wortmann2018) and Uganda (from Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Odama and Wortmann2014). In order to facilitate an informed choice, a brief overview of available response functions is provided below, highlighting their underlying assumptions, strengths and weaknesses. The performances of the different functions were also compared with re-analysis of published datasets on maize, wheat, rice and barley yield responses to N and P from African countries. The datasets were chosen because they represent the typical data used in developing nutrient dose-response models in SSA. Wherever specific examples are cited in this article, they were by no means intended to criticise the authors of the work but only to highlight the issues.
In addition, data on wheat yield response to N on two sites (Betzendorf and Wolfsdorf in Germany) were extracted from Gandorfer and Rajsic (Reference Gandorfer and Rajsic2008). These datasets were chosen because they contain excellent data to the estimate parameters of most response functions (Figure 2) and demonstrate some of the problems encountered in model fitting. The performances of the different models were also compared in terms of the predicted maximum agronomic yield (Y max ) and the nutrient rate (X max ) that gives Y max . Throughout this analysis, statistical inferences were based on the 95% confidence limits (CLs). The width of the 95% CLs was used as a measure of uncertainties around estimated values. Parameters were deemed unreliable when their CLs cover both negative and positive values.
Figure 2. Wheat grain yield (Mg ha−1) response to increasing rates of applied N at the Wolfsdorf (a) and Betzendorf (b) sites, (c) the estimated maximum agronomic yield (Y max ) and (d) the N rates that give the maximum yield (X max ). The dotted vertical grey lines represent the hypothesised boundaries of segments A, B and C. The vertical green lines represent economically optimum rates (EORs) of 201 kg N ha−1 at Wolfsdorf and 109 kg N ha−1 at Betzendorf, calculated using the quadratic model by Gandorfer and Rajsic (Reference Gandorfer and Rajsic2008).
Common Problems and Corrective Measures
Choice of mathematical functions
The development of nutrient dose-response models involves choosing a mathematical function that provides the most accurate prediction of crop response (e.g., yield Y
i
in t ha−1) using a range of nutrient rates (X
i
in kg ha−1). Typically, dose responses show the following patterns: (1) initially Y
i
increases from Y
0
(control yield) at an increasing rate with increase in X
i
; (2) when an inflection point is reached, Y
i
increases at a decreasing rate until Y
max
is reached; and (3) after Y
max
, Y
i
either remains constant or decreases towards 0 if the nutrient rate is increased further. Decreases in Y
i
after Y
max
may occur due to toxicities, nutrient imbalances, salinity, lodging, increased susceptibility to disease, or greater respiration (Agegnehu et al., Reference Agegnehu, Angaw and Agajie2012; Gill et al., Reference Gill, Kanwal and Aziz2004; Khan et al., Reference Khan, Mulvaney and Ellsworth2014; Kindred et al., Reference Kindred, Milne, Webster, Marchant and Sylvester-Bradley2014). For example, high N rates often cause lodging (Agegnehu et al., Reference Agegnehu, Angaw and Agajie2012), whereas high rates of P induce zinc (Zn) deficiency, resulting in yield reduction in maize (Gill et al., Reference Gill, Kanwal and Aziz2004). High potassium (K) rates, when applied as KCl, depress crop yields due to its high salt index (Khan et al., Reference Khan, Mulvaney and Ellsworth2014). The presence of excessive K in the soil can also result in magnesium (Mg) and calcium (Ca) deficiencies. For example, use of high amounts of K fertilisers in sugarcane induced Mg deficiency in South Africa (Rhodes et al., Reference Rhodes, Miles and Hughes2018). Accordingly, three segments may be observed in the response domain, hereafter referred to as A, B and C (Figures 2a, 2b). Segment A represents the region where a steep linear response occurs with incremental nutrient addition due to the higher agronomic use efficiency (AUE) of applied nutrients (Soffe, 2003). AUE is normally calculated as
${{{Y_N} - {Y_O}} \over {{N_f}}}$
where Y
N
is the grain yield (in kg ha−1) in the plot without nutrient limitation, Y
O
is the yield (in kg ha−1) in plots where the nutrient in question was omitted and N
f
is the nutrient input applied (in kg ha−1) (Ladha et al., Reference Ladha, Pathak, Krupnik, Six and van Kessel2005). Segment B represents the region where the rate of response decreases and yield slowly reaches Y
max
as AUE stagnates. Segment C represents the region where AUE decreases rapidly and yields either decline or remain the same with further addition of nutrients. No single function can model these varying response patterns adequately. Therefore, a variety of response functions have been developed by different workers. The functions commonly used in SSA can be classified into four broad categories based on the shape of the assumed curve: (1) linear, quadratic and square root functions; (2) response surface functions; (3) asymptotic functions; and (4) plateau functions. Here, the underlying assumptions, strengths and limitations of the various functions will be discussed to aid informed choice of functions by practitioners.
Linear, quadratic and square root functions
The simple linear function is given as:
$$Y = {b_0} + {b_1}X + \varepsilon $$
where b 0 is the intercept and b 1 is the slope of the line, X is the nutrient rate, and ε is the random error term. In studies in SSA, Cyamweshi et al. (Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018), Daudu et al. (Reference Daudu, Ugbaje, Oyinlola, Tarfa, Yakubu, Amapu and Wortmann2018), Kaizzi et al. (Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Odama and Wortmann2014), Maman et al. (Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018), Senkoro et al. (Reference Senkoro, Tetteh, Kibunja, Ndungu-Magiroi, Quansah, Marandu, Ley, Mwangi and Wortmann2018), Sereme et al. (Reference Serme, Ouattara, Bandaogo and Wortmann2018) and Wortmann et al., (Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018) applied this function to derive nutrient-response models for various crops. The main advantage of the linear function is that it is very easy to fit. Its main limitation is that it does not have Y max contrary to biological reality, and as such, EORs cannot be calculated. Another limitation of this function is that a spurious linear fit can be found when the nutrient rate that achieves Y max was not included in the trial design. Therefore, this function should not be used at all in deriving EORs.
The quadratic and square root functions are simple extensions of equation 1, the difference being only in the last term in equation 2a. The quadratic function assumes a smooth linear increase in response followed by a decrease in response with increasing nutrient rates.
$$Y = {b_0} + {b_1}X - {b_2}{X^2} + \varepsilon $$
Equation 2a implies that there is a critical nutrient level (X
max
) that gives Y
max
beyond which yields decrease; X
max
is calculated as
${X_{max}} = \displaystyle{{ - {b_1}} \over {2{b_2}}}$
. A few authors (e.g., Agegnehu et al., Reference Agegnehu, Nelson, Bird and van Beek2015; Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Odama and Wortmann2014; Wortmann et al., Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018) have used equation 2a in SSA. The main strength of the quadratic function is that it can realistically represent situations where inadequate, adequate or excessive levels of nutrients have been used (Webb, Reference Webb2009). The other advantage of the quadratic function is that it provides an easy way to calculate the EOR (e.g., Bachmaier and Gandorfer, Reference Bachmaier and Gandorfer2012; Gandorfer and Rajsic, Reference Gandorfer and Rajsic2008; Webb, Reference Webb2009). Its limitation is that it makes spurious predictions at the extremes of the nutrient rates, and, as such, it may not be used for extrapolation beyond the observed data. The extrapolation beyond the observed data depends on the crop and the nutrient. For example, maize can support excessive N application rates while potatoes or wheat do not. It has also been shown to consistently over-estimate both X
max
and the yield depression at high nutrient rates (Cerrato and Blackmer, Reference Cerrato and Blackmer1990).
Equation 2a can be reformulated as a square-root function as follows:
$$Y = {b_0} + {b_1}X - {b_2}^{1/2} + \varepsilon $$
where b
0
, b
1
and ε are defined as in equation 1 and b2 is the quadratic or square root coefficient representing yield depression, and X
max
is calculated as
${X_{max}} = \displaystyle{{{b_2}^2} \over {4{b_1}^2}}$
.
The square root function is a purely empirical function with no theoretical support. Its advantage is the ease to fit it just like the quadratic function. Its disadvantage is that X and X 1/2 are co-linear. Therefore, b 1 and b 2 remain non-significant due to co-linearity.
Response surface functions
When two or more nutrients are considered simultaneously, the quadratic function can be extended to response surface models. For example, where a factorial combination of N and P has been tested, equation 2a can be extended as follows (Webb, Reference Webb2009):
$$Y = {b_0} + {b_1}N + {b_2}{\rm{P}} + {b_3}{N^2} + {b_4}{P^2} + {b_5}NP$$
where b 0 is the intercept, b 1 and b 2 are the linear coefficients, b 3 and b 4 are the quadratic coefficients, and b 5 is the interaction coefficient. Similarly, where a factorial combination of three nutrients (e.g., N, P and K) has been tested, equation 3a can be extended as follows:
$$Y = {b_0} + {b_1}N + {b_2}P + {b_3}K + {b_4}{N^2} + {b_5}{P^2} + {b_6}{K^2} + {b_7}NP + {b_8}NK + {b_9}PK + {b_{10}}NPK$$
where b 0 is the intercept, b 1 to b 3 are the linear coefficients, b 4 to b 6 are the quadratic coefficients, and b 7 to b 10 are the interaction coefficients. A response surface function similar to equation 4 was used by Akinnifesi et al. (Reference Akinnifesi, Makumba, Sileshi, Ajayi and Mweta2007) to determine N and P rates and total seasonal rainfall that optimises maize yield response in Malawi. The advantage of these functions is that they can provide insights into interactive effects between nutrients. Their limitation is that they require skills and specialised software to fit and create the desired graphics.
Asymptotic family of functions
This family of functions includes the simple asymptotic and Mitscherlich functions, which assume that dose responses follow Mitscherlich’s law of diminishing returns.
The simple asymptotic function is used with the assumption that response to applied nutrients diminishes with increasing application rates (Patterson, Reference Patterson1956). This has been the most frequently used function in recent publications from SSA (e.g., Cyamweshi et al., Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018; Daudu et al., Reference Daudu, Ugbaje, Oyinlola, Tarfa, Yakubu, Amapu and Wortmann2018; Essel et al., Reference Essel, Abaidoo, Opoku and Ewusi-Mensah2020; Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Odama and Wortmann2014, Reference Kaizzi, Cyamweshi, Kibunja, Senkoro, Nkonde, Maria and Wortmann2018; Maman et al., Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018; Senkoro et al., Reference Senkoro, Tetteh, Kibunja, Ndungu-Magiroi, Quansah, Marandu, Ley, Mwangi and Wortmann2018; Wortmann et al., Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018). This function has been erroneously reported as a ‘curvilinear to plateau’ in most of these publications and as ‘asymptotic quadratic-plus-plateau’ in Essel et al. (Reference Essel, Abaidoo, Opoku and Ewusi-Mensah2020), although it is neither a plateau nor a quadratic response function. This function is given as:
$$Y = a - b{c^X}$$
where a is the asymptotic yield, b is the amplitude (yield increase due to nutrient application), and c is the curvature coefficient (Patterson, Reference Patterson1956); a, b and c must all be positive. In the literature, a large number of spurious values of a and b have been reported. For example, in the database used in Tesfahunegn and Wortmann (Reference Tesfahunegn and Wortmann2017), 1.9%, 6.1% and 25% of the b values for response to N, P and K are ≤ 0 (see the extreme values on the left-hand side of Supplementary Figure S1). Over 31% of c values for N and 65% of c values for P and K are significantly biased downward (c < 0.90). Similarly, in Senkoro et al. (Reference Senkoro, Marandu, Ley and Wortmann2017), 42% and 58% of the b values for response to P and K respectively are ≤ 0. Likewise, in Ndungu-Magiroi et al. (Reference Ndungu-Magiroi, Wortmann, Kibunja, Senkoro, Mwangi, Wamae, Kifuko-Koech and Msakyi2017), 42% of the b values for response to K are ≤ 0. These are biased parameter estimates arising as artefacts of fitting the asymptotic function to data for which it is not suited. A downward bias in c results in a significant downward bias in b and the resultant EOR (see Supplementary Table S1).
The Mitscherlich function has been formulated variously by different workers (Ferreira et al., Reference Ferreira, Zocchi and Baron2017; Harmsen, Reference Harmsen2000; Sorensen, Reference Sorensen1983). The following is the notation commonly reported in the literature:
$$Y = a\left( {1 - exp\left( { - cX} \right)} \right)$$
where a is the asymptotic yield and c is a parameter which controls the steepness of the relationship between X and Y. This model predicts zero yield at zero nutrient input – that is, for X = 0, Y 0 = 0 contrary to the first principle of response (Soffe, Reference Soffe2011) and empirical data (Figures 3 and 4). It also tends to underestimate Y max compared with the other formulations (Table 1).
Figure 3. Maize grain yield (Mg ha−1) response to application rates of N in Niger (from Maman et al., Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018); Tanzania (from Senkoro et al., Reference Senkoro, Marandu, Ley and Wortmann2017); Uganda (from Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Musinguzi, Ebanyat, Hyuha and Wortmann2012); Malawi, Zambia and Kenya (from Wortmann et al., Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018); and Mozambique. The red lines in each figure represent yields that have not yet reached Y max . The blue and black lines represent yields that are declining or remaining constant after reaching their maximum, respectively.
Figure 4. Wheat grain yield (t ha−1) response to N and P application rates in Tanzania and Rwanda (from Cyamweshi et al., Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018), and barley response (e–h) to application rates of N on Holetta and Robge sites (from Agegnehu et al., Reference Agegnehu, Nelson and Bird2016), Welmera, and on two soil types in Jimma (unpublished data) in Ethiopia, and response of rice to application rates of N in Nigeria (from Daudu et al., Reference Daudu, Ugbaje, Oyinlola, Tarfa, Yakubu, Amapu and Wortmann2018) and N in Uganda (from Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Odama and Wortmann2014). The red lines in each figure represent yields that have not yet reached their maximum (asymptote or plateau). The blue and black lines represent yields that are declining or remaining constant after reaching their maximum, respectively.
Table 1. Differences in parameter estimates with different formulations of the asymptotic family of functions. Data on wheat response to N at the Wolfsdorf site in Germany were used to illustrate the differences
Values in red represent parameters estimated with large uncertainty (i.e., 95% CLs covering negative and positive values).
AIC = Akaike information criterion.
A modification of this formulation as follows relaxes the constraint that Y = 0 for all X = 0 by adding the parameter b to account for the indigenous nutrient available in the soil (Gomes, Reference Gomes1953):
$$Y = a\left( {1 - exp\left( { - c\left( {X + b} \right)} \right)} \right)$$
where a is the asymptotic yield, b represents the inherent soil nutrient (in kg ha−1), and c is defined as in equation 6a. Thus, b is the estimated amount of the available amount of the nutrient in the soil at the start of the experiment. Equation 6b is called the Mitscherlich-Baule function (Harmsen, Reference Harmsen2000).
Another modification of Mitscherlich function is given as follows (Sorensen, Reference Sorensen1983):
$$Y = a\left( {1 - b*exp\left( { - cX} \right)} \right)$$
where a and c are defined as in 6a and
$b = 1 - \left( {{{{Y_0}} \over a}} \right)$
. When multiplied by 100, b represents the percentage deficiency of the nutrient in question (Sorensen, Reference Sorensen1983).
A third modification of the Mitscherlich function is given as follows (Dobermann et al., Reference Dobermann, Cruz and Cassman1996).
$$Y = b + a\left( {1 - exp\left( { - cX} \right)} \right)$$
where a is the maximum yield increase due to applied nutrient (in t ha−1), b is the yield (in t ha−1) under no nutrient input (i.e., Y 0 ), and c is a constant related to the use efficiency of soil and fertiliser nutrients. Note that a in equation 6d is the same as b in equation 5.
When dose responses to factorial combinations of two nutrients (e.g., N and P) are modelled, equation 6b can be modified as follows (Harmsen, Reference Harmsen2000):
$$Y = a\left( {1 - exp\left( { - {c_N}\left( {{b_N} + N} \right)} \right)} \right)\left( {1 - {\rm{exp}}\left( { - {c_P}\left( {{b_P} + P} \right)} \right)} \right)$$
where a represent the asymptotic yield, b N and b P represent the indigenous soil available N and P, and c N and c P quantify the increase in yield per unit of N and P applied. Van der Velde et al. (Reference van Der Velde, Folberth, Balkovič, Ciais, Fritz, Janssens, Obersteiner, See, Skalský, Xiong and Peñuelas2014) fitted this function to maize yield response to N and P across 741 locations in Africa.
The main advantage of asymptotic functions is that they have a solid theoretical foundation (Paine et al., Reference Paine, Marthews, Vogt, Purves, Rees, Hector and Turnbull2012), and their parameters can be associated with biologically meaningful processes. Second, they can be used for extrapolation outside the range of observed data (Archontoulis and Miguez, Reference Archontoulis and Miguez2015) because their predictions tend to be more robust than those of the linear and quadratic functions. Their main limitation is that their parameters are estimated through an approximate procedure that starts with best-guess initial values, adjusting parameters through an iterative process. The algorithms often fail to converge with some data, depending on the choice of initial values and the optimisation method (ibid.). Even when the algorithm converges, it could fail to arrive at an optimal solution, and some parameters may still be biased (Tables 1 and 2; Supplementary Figure S2). As such all parameters are treated as approximate values.
Table 2. Effect of trial design on the uncertainty around parameter estimates of the asymptotic function (equation 5)
† The X i represents the sample size, which is N rates in this case.
‡ DF = degrees of freedom.
# AIC = Akaike information criterion; NE = not estimable due to fewer degrees of freedom. Values in red represent parameters estimated with large uncertainty (i.e., 95% CLs covering negative and positive values).
Plateau functions
Plateau functions follow the von Liebig’s law of the minimum where dose response is assumed to increase linearly until the nutrient reaches a critical dose and a ‘plateau’ is reached after this point (Ferreira et al., Reference Ferreira, Zocchi and Baron2017). This pattern may be described by either a linear-plateau or quadratic-plateau function. The linear-plateau function implies a region of linear response followed by a plateau (Anderson and Nelson, Reference Anderson and Nelson1975) as follows:
$$Y = \left\{ {\matrix{ {{b_0} + {b_1}X + \varepsilon \;if\;X \lt {X_{max}}} \cr {{Y_{max}} + \varepsilon \;if\;X \ge {X_{max}}} \cr } } \right.\;$$
where b 0 , b 1 and ε are defined as in equations 1 and 2, Y max is the plateau yield, and X max is the ‘join point’ (i.e., the critical point after which increase in nutrient rates can no longer increase yields). Zheng et al. (Reference Zheng, Mmari, Nishigaki, Kilasara and Funakawa2018) applied this function to determine the optimal soil inorganic N availability to maize in Tanzania.
The quadratic-plateau model implies a region of quadratic response followed by a plateau (Anderson and Nelson, Reference Anderson and Nelson1975; Bullock and Bullock, Reference Bullock and Bullock1994) when the curve reaches its maximum point.
$$Y = {\rm{ }}\left\{ {\matrix{ {{b_0} + {b_1}X + {b_2}{X^2} + \varepsilon \;if\;X \lt {X_{max}}} \cr {{Y_{max}} + \varepsilon \;if\;X \ge {X_{max}}} \cr } } \right.$$
where b
0
, b
1
, b
2
and ε are defined as in equation 2a and Y
max
and X
max
are defined as in equation 7. The join point (X
max
) is related to b
1
and b
2
as
${X_t} = \displaystyle{{ - {b_1}} \over {2{b_2}}}$
just like the critical nutrient level in equation 2a. Similarly, Y
max
is related to b
0
, b
1
and b
2
as
${Y_{max}} = {b_0} - {{b_1^2} \over {4{b_2}}}$
. By substituting Y
max
in equation 8a, the quadratic-plateau function can also be reparametrised as follows:
$$Y = \left\{ {\matrix{ {{b_0} + {b_1}X + {b_2}{X^2}\;if\;X \lt {X_{max}}} \cr {{b_0} - {{b_1^2} \over {4{b_2}}}\;if\;X \ge {X_{max}}} \cr } } \right.$$
The main advantage of plateau functions is that they have a solid theoretical foundation grounded in von Liebig’s law. Their main limitation is that they are more difficult to solve because their parameters have to be estimated using nonlinear optimisation. Another limitation is that the linear-plateau functions result in lower optima compared with the other functions (Alivelu et al., Reference Alivelu, Srivastava, Rao, Singh, Selvakumari and Raju2003; Figure 2). The quadratic plateau function failed to converge in all datasets analysed in Supplementary Table S2.
Parameter estimation and model evaluation
All model parameters including Y max and X max are random variates, and as such, they are subject to various sources of error. Y max , X max can differ widely depending on the trial design and the choice of the response function (Figures 2c, 2d). Failure of the algorithms to converge was one of the common problems encountered during the training workshops. This has been reported widely for the asymptotic function in Cyamweshi et al. (Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018), Daudu et al. (Reference Daudu, Ugbaje, Oyinlola, Tarfa, Yakubu, Amapu and Wortmann2018), Kaizzi et al. (Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Musinguzi, Ebanyat, Hyuha and Wortmann2012), Maman et al. (Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018), Senkoro et al. (Reference Senkoro, Marandu, Ley and Wortmann2017), and Wortmann et al. (Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018). Non-convergence was often caused by inadequate data, incorrect initial parameter values and choice of the minimisation algorithm (Archontoulis and Miguez, Reference Archontoulis and Miguez2015).
Data are said to be inadequate if information is available only for a limited interval of the response domain and/or the sample size is small (e.g., X i ≤ 5) relative to 3–4 parameters to be estimated in most equations. This is particularly true for nonlinear functions (equations 5–8), whose parameters are computed using asymptotic formulae. The software commonly used in recent publications uses the Levenberg-Marquardt-Nash (LMN) algorithm. Unless the starting values are very good, the LMN algorithm takes large, uncontrolled steps and fails to converge. When the algorithm for the asymptotic function failed to converge, some researchers (e.g., Cyamweshi et al., Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018; Daudu et al., Reference Daudu, Ugbaje, Oyinlola, Tarfa, Yakubu, Amapu and Wortmann2018; Kaizzi et al., Reference Kaizzi, Byalebeka, Semalulu, Alou, Zimwanguyizza, Nansamba, Odama and Wortmann2014; Maman et al., Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018; Senkoro et al., Reference Senkoro, Tetteh, Kibunja, Ndungu-Magiroi, Quansah, Marandu, Ley, Mwangi and Wortmann2018; Sereme et al., 2017; Wortmann et al., Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018) have fitted linear models. This is inappropriate because a linear model does not yield Y max and X max ; hence, no EOR can be determined. A more appropriate option is to use another model suited to the data at hand from the set of models given in equations 2–8.
Although rigorous evaluation is an important element of model development, the R2 was the only metric used to judge the performance of dose-response functions in many publications from SSA. Even when they have the same R2 values, models can widely differ in their estimates of Y max , X max , and EORs (Tables 1 and 2; Figure 2). The R2 is inappropriate for assessing the performance of nonlinear models (Bachmaier and Gandorfer, Reference Bachmaier and Gandorfer2012; Cerrato and Blackmer, Reference Cerrato and Blackmer1990; Spiess and Neumeyer, Reference Spiess and Neumeyer2010) or for comparisons among linear, quadratic and nonlinear models. In nonlinear regression, the residual variance and explained variance do not add up to total variance. As such, the R2 does not necessarily fall between 0 and 1; it can take a negative value in nonlinear models (see Table 1 and Supplementary Table S2). That is why the R2 is called ‘pseudo R2’. It can also remain close to 1.0 even when the model is poor and thus unable to distinguish between bad and good models (Tables 1 and 2). This property has also been demonstrated by Spiess and Neumeyer (Reference Spiess and Neumeyer2010), who compared the R2, the bias-corrected Akaike information criterion (AIC) and Bayesian information criterion (BIC) using a predetermined log-logistic model with known parameters (i.e., the correct model) and nine other sigmoid models differing in their numbers of parameters. After performing thousands of simulations, Spiess and Neumeyer (Reference Spiess and Neumeyer2010) demonstrated that the R2 results in choice of the correct model only in 28–43% of the time.
The most important steps in evaluation of models involve the examination of the omnibus test, the significance of model coefficients, and the goodness-of-fit statistics. If the omnibus test is not significant, one or more of the parameters are likely to be nonsignificant, and the whole model may be useless for prediction. The significance of model coefficients should be judged by their 95% CLs. If the CLs cover negative and positive values, the estimate is unreliable (see parameters in red in Tables 1 and 2). When selecting the best model from a cohort of models, the sample corrected Akaike information criterion (AIC c ) is the most appropriate metric (Spiess and Neumeyer, Reference Spiess and Neumeyer2010). However, for trial designs with X i ≤ 5, the AICc cannot be estimated because the degrees of freedom are fewer (see Table 2; Supplementary Figure S2) relative to the parameters (p) to be estimated and the denominator (N-p−1) in equation 9 becomes zero. As a result, the AIC c is undefined.
$$AI{C_c} = AIC + \;{{{\rm{2p}}\left( {{\rm{p + 1}}} \right)} \over {N - p - 1}}$$
where N = sample size (X i ) and p = number of parameters.
Trial design issues
During the series of training workshops facilitated by the author and data from the literature, dose-response models could not be fitted in many instances due to the limitations of the trial design. The common problems are briefly summarised below.
Inadequate control of interactions between nutrients
Interactions between nutrients (Aulakh and Malhi, Reference Aulakh and Malhi2005; Rhodes et al., Reference Rhodes, Miles and Hughes2018; Rietra et al., Reference Rietra, Heinen, Dimkpa and Bindraban2017) and co-limitation are common in nature (ibid.; Tamene et al., Reference Tamene, Sileshi, Ndengu, Mponela, Kihara, Sila and Tondoh2019; Weil, Reference Weil, Weil and Brady2017). Interaction among plant nutrients can yield antagonistic or synergistic outcomes that influence nutrient use efficiency (Rietra et al., Reference Rietra, Heinen, Dimkpa and Bindraban2017). For example, significant interactions exist between applied N, P, K, sulphur and Zn, and these interactions are often synergistic (Aulakh and Malhi, Reference Aulakh and Malhi2005). Ca, Mg and K show antagonistic interactions such that high levels of one or more of these nutrients in the soil can result in decreased uptake of the other (Rhodes et al., Reference Rhodes, Miles and Hughes2018; Rietra et al., Reference Rietra, Heinen, Dimkpa and Bindraban2017). Therefore, dose response of K can only be determined using factorial designs that have carefully balanced Ca, Mg and K taking into account their concentrations in the exchangeable complex. In most studies reviewed, truly factorial trial designs were not applied; many studies varied one macronutrient (mostly N or P) holding the other nutrient constant. Some studies held only N and P rates constant when determining the dose response for K (e.g., Ndungu-Magiroi et al., Reference Ndungu-Magiroi, Wortmann, Kibunja, Senkoro, Mwangi, Wamae, Kifuko-Koech and Msakyi2017; Senkoro et al., Reference Senkoro, Marandu, Ley and Wortmann2017). Due to nutrient interactions and the confounding effects of indigenous soil supply, sometimes unusual patterns of response may be observed such as those on the Selian site in Figures 1b, 4a, and 4b. The negative parameter values for K response and decline in sole maize yield with increase in K rate on the low potential site in Kenya (e.g., Figure 2 in Ndungu-Magiroi et al., Reference Ndungu-Magiroi, Wortmann, Kibunja, Senkoro, Mwangi, Wamae, Kifuko-Koech and Msakyi2017) are also indicative of such problems.
Inadequate definition of the response domain
The response domain is the X i range in which the function is valid – that is, where segment A–C can be observed (Figure 2a, 2b). However, in many trials, nutrient application was limited within segment A–B – that is, the linear part of the response function (Figures 3–4). In some studies, high nutrient rates were deliberately excluded from the trial designs as these rates were considered uneconomical a priori. For example, the maximum N rate for maize was 75 kg N ha−1 in the FURP trials in Kenya (KARI, 1994). This represents a severe case of truncation of the response domain. In the most recent fertiliser optimisation trials for maize, wheat, barley and rice across SSA (Cyamweshi et al., Reference Cyamweshi, Nabahungu, Senkoro, Kibunja, Mukuralinda, Kaizzi, Mvuyekure, Kayumba, Ndungu-Magiroi, Koech and Wortmann2018; Maman et al., Reference Maman, Traore, Garba, Dicko, Gonda and Wortmann2018; Wortmann et al., Reference Wortmann, Senkoro, Cyamweshi, Kibunja, Nkonde, Munthali, Nalivata, Nabahungu and Kaizzi2018), the maximum N rate was set at 120 kg N ha−1 although these cereals positively respond to up to 200 kg N ha−1 (Figures 3 and 4). This means that N rates were limited to segment A and Y max was not achieved on many sites (represented by red lines in Figures 3 and 4). If the trial design does not cover segments A–C, model parameters cannot be determined correctly (see below).
On inherently fertile sites such as Dareda and Selian in Tanzania (e.g., Figures 1b, 3b), Musanze (Figures 3d, 4c) and Gahunga (Figure 4d) in Rwanda, and Eldoret in Kenya (Figure 3g), initial yield increases due to external nutrient inputs may not be rapid because the high level of indigenous soil nutrients has already raised the curve beyond the inflection point. On such sites, nutrient rates need to be increased further to reach Y max and allow accurate estimation of the EOR.
Inadequate number and spacing of nutrient levels
The number and spacing of X i in the trial design can affect the quality of the final model to predict responses accurately. In most trials, four or five equally spaced nutrient levels were used (Figures 1, 3–4). When X i < 6, the algorithm may fail to converge, model parameters may be biased or can have large uncertainty due to the inadequate degrees of freedom. When the design involves wide spacing of nutrient levels, especially in segment A (e.g., 0 and 50 kg N ha−1 in Figure 1a), it may not be possible to capture the inflection points between segments A and B and hence large uncertainty (very wide confidence limits) in parameter estimations and predictions may occur (see the Abongomola and Kwera sites Supplementary Figure 2).
To demonstrate this point, the wheat yield data from Wolfsdorf and Betzendorf were re-analysed assuming five scenarios. Scenario 1 represents the complete data where the response domain covers segments A–C. Scenario 2 was created by excluding the highest rate (200 kg N ha−1), and Scenario 3 by excluding the two high rates (160 and 200 kg N a−1) to represent a design where only segments A–B were covered. This is similar to most designs in recent publications (Figures 3 and 4). In addition, Scenario 4 was created by removing 40 and 160 kg N ha−1 to represent a trial design with wider spacing between levels similar to the one in Figure 1a. I fitted the commonly used asymptotic function (equation 5) to analyse these data and a dataset from Uganda collected using the same design as Figure 1a (Supplementary Figure S2).
These analyses revealed striking differences in parameter estimates between scenarios (Table 2). With removal of the higher N rates, a and b were biased upwards, with increasing uncertainty (wider 95% CLs) of the parameter from Scenario 1 to Scenario 3 (Table 2). Scenario 4 on the other hand led to downward bias of a and b (Table 2). The small sample size (X i = 4) and inadequate spacings of the X i (e.g., in the data from Uganda) were responsible for the poor parameter estimates and large uncertainties in predictions of yield (Supplementary Figure S2a). These results highlight the fact that trial designs with truncation of X i and wide spacing can produce biased estimates of model parameters of the asymptotic functions and subsequent EORs. Therefore, I recommend a minimum of six carefully spaced nutrient levels with closer spacing in segments A–B (where responses are expected to change rapidly), and wider spacing in segments C (where responses change slowly). Although one may not know the limits of segments A, B and C a priori, it is possible to define these segments based on earlier reports or data (e.g., Figures 1, 3 and 4).
Inadequate consideration of soil type in trial design
Soil type, which varies widely across SSA (Jones et al., Reference Jones, Breuning-Madsen, Brossard, Dampha, Deckers, Dewitte, Hallett, Jones, Kilasara, Le Roux, Micheli, Montanarella, Spaargaren, Tahar, Thiombiano, van Ranst, Yemefack and Zougmore2013), plays a key role in crop yield response (Elias et al., Reference Elias, Okoth and Smaling2019; Kihara et al., Reference Kihara, Sileshi, Nziguheba, Kinyua, Zingore and Sommer2017; Sileshi et al., Reference Sileshi, Akinnifesi, Debusho, Beedy, Ajayi and Mong’omba2010), nutrient use efficiency (Sileshi et al., Reference Sileshi, Jama, Vanlauwe, Negassa, Harawa, Kiwia and Kimani2019), and the profitability of fertiliser use (Kihara et al., Reference Kihara, Bolo, Kinyua, Rurinda and Piikki2020). Figure 1 highlights differences in response that may be associated with soil type. For example, response to applied N is still increasing at Kapchorwa and Dareda located on Nitisols and Phaeozems, respectively, soils which are ranked excellent for maize. On the other hand, yields have levelled off at Bulindi located on Ferralsols, Kawanda, Ngetta and Tororo located on Plinthosols, which are ranked as marginal for maize cultivation (see page 45 in the harmonised soil atlas of Africa; Jones et al., Reference Jones, Breuning-Madsen, Brossard, Dampha, Deckers, Dewitte, Hallett, Jones, Kilasara, Le Roux, Micheli, Montanarella, Spaargaren, Tahar, Thiombiano, van Ranst, Yemefack and Zougmore2013). In the literature reviewed, most of the nutrient dose-response models were developed without sufficient information on the soil type and the inherent soil nutrient supply. The size of Y 0 and Y max depends on the indigenous nutrient available in the soil, which in turn depends on the soil type. For example, Andosols, Cambisols, Chernozems, Fluvisols, Luvisols, Nitisols, Phaeozems, and Vertisols are considered ‘excellent’ soils for maize production (Jones et al., Reference Jones, Cross, Withers, DeLuca, Robinson, Quilliam, Harris, Chadwick and Edwards-Jones2013), achieving high values of Y 0 in SSA. On the other hand, Ferralsols, Acrisols, Alisols, Arenosols, Leptosols, Lixisols and Plinthosols, which are considered ‘marginal’ or ‘poor’ (Jones et al., Reference Jones, Cross, Withers, DeLuca, Robinson, Quilliam, Harris, Chadwick and Edwards-Jones2013) often achieve very low values of Y 0 . It must also be noted a huge variation in the content of total and available N, P and K within a soil reference group. For example, the availability of P may vary across and within the same soil group (Batjes, Reference Batjes2011). In addition, the availability of mineral N in the soil during the crop cycle may depend on the cropping history or the land use. For example, with repeated cultivation yield declines more rapidly on Acrisols and Ferralsols than on Cambisols; the rate of decline depending on soil cover (Sileshi et al., Reference Sileshi, Akinnifesi, Debusho, Beedy, Ajayi and Mong’omba2010). Therefore, interpreting fertiliser treatments without taking into account sites which have different soil types, cropping history and properties may generate conclusive results (Elias et al., Reference Elias, Okoth and Smaling2019). In future, nutrient-response trials should be designed taking into account the predominant soil types, agroecology, land use and cropping history in an area. It is strongly recommended that Reference Soil Groups (IUSS Working Group WRB, 2014) and cropping history be carefully considered when establishing nutrient-response trials.
Knowledge of the gap between the nutrients provided by the inherent soil supply and the crop demand is also crucial when establishing fertiliser trials. If the site is inherently fertile, the response to nutrient inputs will be small. Thus, at low nutrient input levels, responses may be too small and the EOR may not be correctly determined. For example, the control yields were already very high (>3 t ha−1) on the Selian site in Tanzania, which was located on Chernozems – one of the most productive soil in SSA (Figure 1b). Thus, the maximum rate of 120 kg N ha−1 applied to maize was probably too low to determine Y max as the inherent soil nutrient supply was probably sufficient to sustain yields >3 t ha−1. On such sites, unless the trial design included higher nutrient rates (e.g., ≥200 kg N ha−1), it may not be possible to establish the underlying crop response curve (Figure 1b; Supplementary Table S2).
Conclusions
The review and analyses have demonstrated that the choice of inappropriate response functions and inadequate trial designs can create uncertainty around model parameters and EORs. If the objective is to determine response to a single nutrient, the modified Mitscherlich function should be used in preference to the other functions as it provides additional information on the proportional deficiency of the nutrient in question and the expected yield in the no-input control. In any case, the choice of any given functions should be justified using the AIC c criterion. In the past, fertiliser trials have focused on the response of one nutrient by holding the other nutrient constant. The use of factorial designs is recommended to address interactions between two or more nutrients. When dose responses to factorial combinations of two nutrients are modelled, the recommended functions are equations 3 and 6e. Where response to three or more nutrients is to be modelled, the recommended function is equation 4. It is also recommended that a minimum of six carefully spaced nutrient levels be used for proper estimation of EORs. In addition, Reference Soil Groups and cropping history should be carefully considered when locating field trials.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S0014479721000193
Acknowledgements
I sincerely thank scientists who participated in the training workshops on fertiliser optimisation and shared their challenges in analysis of trials. Many thanks also to the reviewers whose comments helped to improve the manuscript substantially.
Funding Support
This work was supported, in whole or in part, by the Bill & Melinda Gates Foundation [INV-005460]. Under the grant conditions of the Foundation, a Creative Commons Attribution 4.0 Generic License has already been assigned to the Author Accepted Manuscript version that might arise from this submission.
Compliance with Ethical Standards
This study received no funding from any research grant. The author has no conflict of interest.
Open access
