Dose-response models to guide site-specific nutrient management and lessons for fertiliser trial design in sub-Saharan Africa

Summary

Optimisation of fertiliser use and site-specific nutrient management are increasingly becoming critical because of the growing need to balance agricultural productivity with the growing demand for food and environmental concerns. Trials to determine responses of crops to fertilisers have been widely conducted in sub-Saharan Africa (SSA) with increasing emphasis on the development of economically optimum rates (EORs). Computation of EORs depends on accurate estimation of both the optimum nutrient rate and the agronomic maximum yield response; however, estimation of nutrient-response parameters and EORs is beset by a number of problems. Therefore, the objectives of this paper were to (1) point out common problems in the development and use of nutrient dose-response models and (2) provide corrective measures to facilitate future trial design and data analysis. This review outlines the underlying assumptions, strengths and limitations of the various response functions in order to facilitate informed choices by practitioners. Using specific examples, it also shows that (1) the commonly used trial designs do not allow examination of interactions between two or more nutrients and (2) trial designs with ≤5 nutrient levels and wide spacing between the levels result in large uncertainty in dose-response parameters. The key recommendations emerging from the review are as follows: (1) factorial designs and response surface models should be used more widely to address interactions between nutrients; (2) a minimum of six carefully spaced nutrient levels should be used to correctly estimate dose-response parameters; and (3) when locating field trials, Reference Soil Groups and cropping history should be carefully considered to produce site-specific EORs.

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., 2013; Palm et al., 2017). 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., 2017, 2020). 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., 2018; Snapp et al., 2003). 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., 2018; Kihara et al., 2017, 2020; Tamene et al., 2019).

Governments, international agricultural research institutions, and donors continue to invest in fertiliser trials across SSA (Kaizzi et al., 2016; Snapp et al., 2003). 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., 2016). 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⁻¹ against the recommended rate of 75 kg N ha⁻¹ for a maize crop with an asymptotic yield of 3.7 t ha⁻¹ in the Kenyan Rift Valley (Supplementary Table S1). In the same area, EOR was estimated at 0 and 8 kg P ha⁻¹ for a maize crop with an asymptotic yield of 4.7 and 6.1 t ha⁻¹, respectively (Supplementary Table S1). Similarly, the EOR was reported to be 27 kg N ha⁻¹ and 0 kg P ha⁻¹ for a maize crop with an asymptotic yield of 4.7 and 3.4 t ha⁻¹, 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⁻¹ maize grain harvest requires the uptake of approximately 100 kg N ha⁻¹, 24 kg P ha⁻¹ and 85 kg K ha⁻¹ (Nalivata et al., 2016). 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, 2011): (1) some yield (Y ₀ ) 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, 2011). 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, 2008; Morris et al., 2018). The assumptions behind the different response functions and the disagreements among their predictions have also received little attention (Cerrato and Blackmer, 1990).

Crop response can significantly vary from site to site or season to season within a site (Albarenque et al., 2016; Xu et al., 2019), 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, 2015). As such, it is often difficult to know whether or not the estimated parameters represent the underlying reality (ibid.; Bachmaier and Gandorfer, 2012). 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., 2012) and Tanzania (data from Senkoro et al., 2017). 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., 2018), Malawi and Rwanda (from Wortmann et al., 2018), Tanzania (from Senkoro et al., 2017) and Uganda (from Kaizzi et al., 2012); wheat yield response to N in Tanzania and Rwanda (from Cyamweshi et al. 2018); barley yield response to N (from Agegnehu et al., 2016) and P (from Abdulkadir et al. 2017) in Ethiopia; and rice yield response to N, P and Zn in Nigeria (from Daudu et al., 2018) and Uganda (from Kaizzi et al., 2014). 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 (2008). 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⁻¹) 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⁻¹ at Wolfsdorf and 109 kg N ha⁻¹ at Betzendorf, calculated using the quadratic model by Gandorfer and Rajsic (2008).

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⁻¹) using a range of nutrient rates (X _i in kg ha⁻¹). Typically, dose responses show the following patterns: (1) initially Y _i increases from Y ₀ (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., 2012; Gill et al., 2004; Khan et al., 2014; Kindred et al., 2014). For example, high N rates often cause lodging (Agegnehu et al., 2012), whereas high rates of P induce zinc (Zn) deficiency, resulting in yield reduction in maize (Gill et al., 2004). High potassium (K) rates, when applied as KCl, depress crop yields due to its high salt index (Khan et al., 2014). 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., 2018). 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⁻¹) in the plot without nutrient limitation, Y _O is the yield (in kg ha⁻¹) in plots where the nutrient in question was omitted and N _f is the nutrient input applied (in kg ha⁻¹) (Ladha et al., 2005). 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:

(1) $$Y = {b_0} + {b_1}X + \varepsilon $$

where b ₀ is the intercept and b ₁ is the slope of the line, X is the nutrient rate, and ε is the random error term. In studies in SSA, Cyamweshi et al. (2018), Daudu et al. (2018), Kaizzi et al. (2014), Maman et al. (2018), Senkoro et al. (2018), Sereme et al. (2018) and Wortmann et al., (2018) 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.

(2a) $$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., 2015; Kaizzi et al., 2014; Wortmann et al., 2018) 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, 2009). The other advantage of the quadratic function is that it provides an easy way to calculate the EOR (e.g., Bachmaier and Gandorfer, 2012; Gandorfer and Rajsic, 2008; Webb, 2009). 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, 1990).

Equation 2a can be reformulated as a square-root function as follows:

(2b) $$Y = {b_0} + {b_1}X - {b_2}^{1/2} + \varepsilon $$

where b ₀ , b ₁ 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 ₁ and b ₂ 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, 2009):

(3) $$Y = {b_0} + {b_1}N + {b_2}{\rm{P}} + {b_3}{N^2} + {b_4}{P^2} + {b_5}NP$$

where b ₀ is the intercept, b ₁ and b ₂ are the linear coefficients, b ₃ and b ₄ are the quadratic coefficients, and b ₅ 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:

(4) $$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 ₀ is the intercept, b ₁ to b ₃ are the linear coefficients, b ₄ to b ₆ are the quadratic coefficients, and b ₇ to b ₁₀ are the interaction coefficients. A response surface function similar to equation 4 was used by Akinnifesi et al. (2007) 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, 1956). This has been the most frequently used function in recent publications from SSA (e.g., Cyamweshi et al., 2018; Daudu et al., 2018; Essel et al., 2020; Kaizzi et al., 2014, 2018; Maman et al., 2018; Senkoro et al., 2018; Wortmann et al., 2018). 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. (2020), although it is neither a plateau nor a quadratic response function. This function is given as:

(5) $$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, 1956); 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 (2017), 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. (2017), 42% and 58% of the b values for response to P and K respectively are ≤ 0. Likewise, in Ndungu-Magiroi et al. (2017), 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., 2017; Harmsen, 2000; Sorensen, 1983). The following is the notation commonly reported in the literature:

(6a) $$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 contrary to the first principle of response (Soffe, 2011) 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⁻¹) response to application rates of N in Niger (from Maman et al., 2018); Tanzania (from Senkoro et al., 2017); Uganda (from Kaizzi et al., 2012); Malawi, Zambia and Kenya (from Wortmann et al., 2018); 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⁻¹) response to N and P application rates in Tanzania and Rwanda (from Cyamweshi et al., 2018), and barley response (e–h) to application rates of N on Holetta and Robge sites (from Agegnehu et al., 2016), 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., 2018) and N in Uganda (from Kaizzi et al., 2014). 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, 1953):

(6b) $$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⁻¹), 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, 2000).

Another modification of Mitscherlich function is given as follows (Sorensen, 1983):

(6c) $$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, 1983).

A third modification of the Mitscherlich function is given as follows (Dobermann et al., 1996).

(6d) $$Y = b + a\left( {1 - exp\left( { - cX} \right)} \right)$$

where a is the maximum yield increase due to applied nutrient (in t ha⁻¹), b is the yield (in t ha⁻¹) under no nutrient input (i.e., Y ₀ ), 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, 2000):

(6e) $$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. (2014) 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., 2012), 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, 2015) 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., 2017). 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, 1975) as follows:

(7) $$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 ₀ , b ₁ 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. (2018) 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, 1975; Bullock and Bullock, 1994) when the curve reaches its maximum point.

(8a) $$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 ₀ , b ₁ , b ₂ 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 ₁ and b ₂ 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 ₀ , b ₁ and b ₂ 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:

(8b) $$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., 2003; 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. (2018), Daudu et al. (2018), Kaizzi et al. (2012), Maman et al. (2018), Senkoro et al. (2017), and Wortmann et al. (2018). Non-convergence was often caused by inadequate data, incorrect initial parameter values and choice of the minimisation algorithm (Archontoulis and Miguez, 2015).

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., 2018; Daudu et al., 2018; Kaizzi et al., 2014; Maman et al., 2018; Senkoro et al., 2018; Sereme et al., 2017; Wortmann et al., 2018) 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 R² was the only metric used to judge the performance of dose-response functions in many publications from SSA. Even when they have the same R² values, models can widely differ in their estimates of Y _max , X _max , and EORs (Tables 1 and 2; Figure 2). The R² is inappropriate for assessing the performance of nonlinear models (Bachmaier and Gandorfer, 2012; Cerrato and Blackmer, 1990; Spiess and Neumeyer, 2010) 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 R² 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 R² is called ‘pseudo R²’. 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 (2010), who compared the R², 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 (2010) demonstrated that the R² 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, 2010). 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.

(9) $$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, 2005; Rhodes et al., 2018; Rietra et al., 2017) and co-limitation are common in nature (ibid.; Tamene et al., 2019; Weil, 2017). Interaction among plant nutrients can yield antagonistic or synergistic outcomes that influence nutrient use efficiency (Rietra et al., 2017). For example, significant interactions exist between applied N, P, K, sulphur and Zn, and these interactions are often synergistic (Aulakh and Malhi, 2005). 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., 2018; Rietra et al., 2017). 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., 2017; Senkoro et al., 2017). 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., 2017) 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⁻¹ 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., 2018; Maman et al., 2018; Wortmann et al., 2018), the maximum N rate was set at 120 kg N ha⁻¹ although these cereals positively respond to up to 200 kg N ha⁻¹ (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⁻¹ 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⁻¹), and Scenario 3 by excluding the two high rates (160 and 200 kg N  a⁻¹) 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⁻¹ 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., 2013), plays a key role in crop yield response (Elias et al., 2019; Kihara et al., 2017; Sileshi et al., 2010), nutrient use efficiency (Sileshi et al., 2019), and the profitability of fertiliser use (Kihara et al., 2020). 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., 2013). 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 ₀ 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., 2013), achieving high values of Y ₀ in SSA. On the other hand, Ferralsols, Acrisols, Alisols, Arenosols, Leptosols, Lixisols and Plinthosols, which are considered ‘marginal’ or ‘poor’ (Jones et al., 2013) often achieve very low values of Y ₀ . 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, 2011). 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., 2010). 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., 2019). 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⁻¹) 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⁻¹ 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⁻¹. On such sites, unless the trial design included higher nutrient rates (e.g., ≥200 kg N ha⁻¹), 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.