Decomposing global crop yield variability

Our dataset includes yield and area time series extracted from the FAOSTAT database [19] for four crops species (maize, rice, soybean and wheat) in all FAOSTAT world regions excluding Melanesia, Micronesia, Polynesia and the Caribbean. We thus consider 17 regions for maize, rice and soybean and 18 for wheat (adding Northern Europe). For most regions the time series extend from 1961 to 2012, except in Central Asia (1992–2012).

We define a yield residual ${{r}_{it}}$ at year t in a region i as

$$\begin{eqnarray}&&{{r}_{it}}={{Y}_{it}}-{{\mu }_{it}}\;,\end{eqnarray} \tag{ 1 }$$

where ${{Y}_{it}}$ is the observed yield (t/ha) and ${{\mu }_{it}}$ is the expected yield value (t/ha) at year t for a given crop in the ith region. Values of ${{\mu }_{it}}$ are estimated using linear, quadratic or cubic regression models. These models are fitted by ordinary least squares, using the R function lm for each yield time series for each crop species at each geographical unit. The model showing the lowest Akaike Criteria (AIC) is selected to detrend yield time series for each crop and each region. We check for absence of autocorrelation in residuals.

For simplicity, we use the word yield 'anomaly' instead of yield 'residual' to characterize detrended yield time series.

We estimate regional yield interannual variability ${{V}_{i}}$ for a given crop from the variance of regional yield anomalies as

$$\begin{eqnarray}&&{{V}_{i}}=\frac{1}{T}\sum \limits_{t=1}^{T} {{\left( {{r}_{it}}-\overline{{{r}_{i}}} \right)}^{2}}=\frac{1}{T}\sum \limits_{t=1}^{T} r_{it}^{2},\end{eqnarray} \tag{ 2 }$$

where $\overline{{{r}_{i}}}$ is time-average of anomalies (equal to zero by construction). Note that we also compute a standard deviation for each region ${{s}_{i}}=\sqrt{\frac{1}{T}\mathop \sum \nolimits_{t=1}^{T} r_{it}^{2}}$ (figure 1.)

Zoom In Zoom Out Reset image size

We define below a mathematical expression (equation (5)) connecting regional yield anomalies (equation (1)) to global yield anomalies. We calculate the variance of global yield in a second step (equation (6)). Two limit cases are presented at the end of the section.

Let R_t be the global anomaly (i.e., the difference between global yield and global expected yield value at the world level at time t) for a given crop, ${{R}_{t}}={{\bar{Y}}_{t}}-{{\bar{\mu }}_{t}}$, with ${{\bar{Y}}_{t}}=\mathop \sum \nolimits_{i=1}^{N} {{\omega }_{it}}{{Y}_{it}}$, ${{\bar{\mu }}_{t}}=\mathop \sum \nolimits_{i=1}^{N} {{\omega }_{it}}{{\mu }_{it}}$, ${{\omega }_{it}}=\frac{{{S}_{it}}}{\mathop \sum \nolimits_{i=1}^{N} {{S}_{it}}}$ the proportion of the world cultivated area allocated to the ith region at year t, and N the total number of regions where the considered crop is grown.

The squared value of ${{R}_{t}}$ is defined by:

$$\begin{eqnarray*}&&\begin{array}{l} R_{t}^{2}={{\left( {{{\bar{Y}}}_{t}}-{{{\bar{\mu }}}_{t}} \right)}^{2}} \\ ={{\left( \sum \limits_{i=1}^{N} {{\omega }_{it}}{{Y}_{it}}-\sum \limits_{i=1}^{N} {{\omega }_{it}}{{\mu }_{it}} \right)}^{2}}={{\left( \sum \limits_{i=1}^{N} {{\omega }_{i,t}}\left( {{Y}_{it}}-{{\mu }_{it}} \right) \right)}^{2}}\;. \\ \end{array}\end{eqnarray*}$$

As

$$\begin{eqnarray*}&&\begin{array}{} {{\left( \mathop \sum \nolimits_{i=1}^{N} {{\omega }_{it}}\left( {{Y}_{it}}-{{\mu }_{it}} \right) \right)}^{2}}={\mkern 1mu} \mathop \sum \nolimits_{i=1}^{N} {{\left( {{\omega }_{it}}\left( {{Y}_{it}}-{{\mu }_{it}} \right) \right)}^{2}}\qquad \qquad \quad \ \ \; \\ \qquad \qquad \qquad \qquad \quad \ \ \ {\mkern 1mu} +{\mkern 1mu} {\rm 2}\sum j\gt i {{\omega }_{it}}{{\omega }_{jt}}\left( {{Y}_{it}}-{{\mu }_{it}} \right)\left( {{Y}_{jt}}-{{\mu }_{jt}} \right), \\ \end{array}\end{eqnarray*}$$

we have

$$\begin{eqnarray}&&R_{t}^{2}=\sum \limits_{i=1}^{N} {{\omega }^{2}}_{it}{{r}_{it}}^{2}+2\mathop \sum \limits_{j\gt i} {{\omega }_{it}}{{\omega }_{jt}}{{c}_{ijt}}\end{eqnarray} \tag{ 3 }$$

with $r_{it}^{2}={{\left( {{Y}_{it}}-{{\mu }_{it}} \right)}^{2}}$ (squared regional yield anomaly in t²/ha²), ${{c}_{ijt}}$, the product of regional yield anomalies in regions i and j, i.e., ${{c}_{ijt}}=\left( {{Y}_{it}}-{{\mu }_{it}} \right)\left( {{Y}_{jt}}-{{\mu }_{jt}} \right)$ in t²/ha².

We introduce ${{\delta }_{it}}$ as the difference between squared yield anomaly in the ith region $r_{it}^{2}$ and the mean of squared regional yield anomalies over all regions at each time step (i.e., the interregional yield variance ${{v}_{t}}=\frac{1}{N}\mathop \sum \nolimits_{i=1}^{N} r_{it}^{2}$), further noted ${{\delta }_{it}}=r_{it}^{2}-{{v}_{t}}$ (t²/ha²). This difference can be negative (or positive) if the squared residual in a given region is higher (smaller) than the average over all regions.

We then get from equation (3)

$$\begin{eqnarray}\begin{array}{lllllllllllllll} R_{t}^{2} & = & \sum \limits_{i=1}^{N} {{\omega }^{2}}_{it}\left( {{\delta }_{it}}+{{v}_{t}} \right)+2\mathop \sum \limits_{j\gt i} {{\omega }_{it}}{{\omega }_{jt}}{{c}_{ijt}} \\ {} & = & {{v}_{t}}\sum \limits_{i=1}^{N} {{\omega }^{2}}_{it}+\sum \limits_{i=1}^{N} {{\omega }^{2}}_{it}{{\delta }_{it}}+2\mathop \sum \limits_{j\gt i} {{\omega }_{it}}{{\omega }_{jt}}{{c}_{ijt}}. \\ \end{array}\end{eqnarray} \tag{ 4 }$$

As the Simpson diversity index can be defined by ${{S}_{t}}=\mathop \sum \nolimits_{i=1}^{N} {{\omega }^{2}}_{it}$, we finally get from (4)

$$\begin{eqnarray}&&R_{t}^{2}={{v}_{t}}{{S}_{t}}+\sum \limits_{i=1}^{N} {{\omega }^{2}}_{it}{{\delta }_{it}}+2\mathop \sum \limits_{j\gt i} {{\omega }_{it}}{{\omega }_{jt}}{{c}_{ijt}}.\end{eqnarray} \tag{ 5 }$$

Equation (5) decomposes squared global yield anomalies at a given year as a sum of three terms (in t²/ha²): ${{\Gamma }_{1t}}={{v}_{t}}{{S}_{t}}$, ${{\Gamma }_{2t}}=\mathop \sum \nolimits_{i=1}^{N} {{\omega }^{2}}_{it}{{\delta }_{it}}$, and ${{\Gamma }_{3t}}={\rm 2}\sum j\gt i {{\omega }_{it}}{{\omega }_{jt}}{{c}_{ijt}}$. The yield variance at the global scale is then defined by the time average of squared global yield anomalies $V=\frac{1}{T}\mathop \sum \nolimits_{t=1}^{T} R_{t}^{2}$ expressed as

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} V=\frac{1}{T}\sum \limits_{t=1}^{T} \left( {{v}_{t}}{{S}_{t}}+\sum \limits_{i=1}^{N} {{\omega }^{2}}_{it}{{\delta }_{it}}+2\mathop \sum \limits_{j\gt i} {{\omega }_{it}}{{\omega }_{jt}}{{c}_{ijt}} \right) \\ \quad =\overline{{{\Gamma }_{1\bullet }}}+\overline{{{\Gamma }_{2\bullet }}}+\overline{{{\Gamma }_{3\bullet }}}\;. \\ \end{array}\end{eqnarray} \tag{ 6 }$$

Equation (6) relates global yield variance V to three additive components $\overline{{{\Gamma }_{1\bullet }}},\;\overline{{{\Gamma }_{2\bullet }},\;}\overline{{{\Gamma }_{3\bullet }}}$. Each of these components depends on yield anomalies and on the proportions of cultivated areas and are all expressed in t²/ha². Note that although none of the individual components actually corresponds to a variance, their sum does. The three components can be estimated independently using regional yield anomalies and proportion of cultivated areas derived from FAOSTAT data. Importantly we also compute values of ${{\Gamma }_{2it}}={{\omega }^{2}}_{it}{{\delta }_{it}}$, i = 1,..., N, i.e., the regional contributions to component 2, for each region, but only present area proportions for the most important regions in figure 3.

There are a few limit cases of equation (5) that can be explored. In particular we investigate two extreme situations in terms of land allocation. In the first limit case, all cultivated areas are allocated in one region only (region k):

Limit case 1:

$$\begin{eqnarray*}&&\begin{array}{} {{\omega }_{kt}}=1,\;{{S}_{t}}=1,\;{{\Gamma }_{1t}}={{v}_{t}},\;{{\Gamma }_{2t}}={{\delta }_{kt}}=r_{kt}^{2}-{{v}_{t}},\;{{\Gamma }_{3}}=0, \\ R_{t}^{2}=r_{kt}^{2}\ {\rm and}\ V=\frac{1}{T}\mathop \sum \nolimits_{t=1}^{T} r_{kt}^{2}\;. \\ \end{array}\end{eqnarray*}$$

In this case, the global yield variance is simply equal to the yield variance in the region where all the cultivated area is allocated (region k).

In the second limit case, all cultivated areas are evenly distributed among the N world regions:

Limit case 2:

$$\begin{eqnarray*}&&\begin{array}{lllllllllllllll} {{\omega }_{it}}=\frac{1}{N},\;i\;=\;1,\ldots ,\;{{S}_{t}}=\frac{1}{N},\;{{\Gamma }_{1t}}=\frac{{{v}_{tt}}}{N},\;{{\Gamma }_{3t}}={\rm 2}\frac{1}{{{N}^{2}}}\sum j\gt i {{c}_{ijt}}, \\ {{\Gamma }_{2t}}=\mathop \sum \nolimits_{i=1}^{N} {{\omega }^{2}}_{it}{{\delta }_{it}}=\mathop \sum \nolimits_{i=1}^{N} \frac{1}{{{N}^{2}}}\left( r_{it}^{2}-{{v}_{t}} \right)=\frac{1}{{{N}^{2}}}\mathop \sum \nolimits_{i=1}^{N} r_{it}^{2}-\frac{1}{{{N}^{2}}}\mathop \sum \nolimits_{i=1}^{N} {{v}_{t}}. \\ \end{array}\end{eqnarray*}$$

As

$$\begin{eqnarray*}&&\begin{array}{lllllllllllllll} \;{{v}_{t}}=\frac{1}{N}\mathop \sum \nolimits_{i=1}^{N} r_{it}^{2}, \\ \end{array}\end{eqnarray*}$$

we get

$$\begin{eqnarray*}&&\begin{array}{} {{\Gamma }_{2t}}=\frac{1}{{{N}^{2}}}\mathop \sum \nolimits_{i=1}^{N} r_{it}^{2}-\frac{1}{{{N}^{2}}}\mathop \sum \nolimits_{i=1}^{N} {{v}_{t}} \\ \ =\;\frac{1}{{{N}^{2}}}N{{v}_{t}}-\frac{1}{{{N}^{2}}}N{{v}_{t}}=0 \\ \end{array}\end{eqnarray*}$$

In this second limit case, we thus get $R_{t}^{2}=\frac{{{v}_{t}}}{N}+2\frac{1}{{{N}^{2}}}\sum j\gt i {{c}_{ijt}}$ and $V=\frac{1}{T}\mathop \sum \nolimits_{t=1}^{T} \frac{{{v}_{t}}}{N}+2\frac{1}{{{N}^{2}}}\sum j\gt i \mathop \sum \nolimits_{t=1}^{T} {{c}_{ijt}}$. The value of $V$thus depends on the mean of squared regional yield anomalies and on co-variation between yield anomalies in distinct regions. If the covariance between ${{r}_{it}}$ and ${{r}_{jt}}$ is negligible across time, then $V=\frac{1}{T}\mathop \sum \nolimits_{t=1}^{T} \frac{{{v}_{t}}}{N}$.

Regional yield anomalies and proportions of cultivated areas are calculated from FAOSTAT data for the period 1961–2012. This enables us to calculate regional yield variance ${{V}_{i}}$ (equation (2)) for each crop, and global yield variance V (equation (6)) with its three components, for each crop over the entire time period (see result section) and over shorter time windows (see SI). Note that a comparison of the results obtained for different time periods shows that the respective values of each components are fairly stable in time.

Figure 1 shows the time distribution of regional yield anomalies (equation (1)) for maize, rice, soybean, and wheat (panels a to d). These distributions differ among crop species; differences between minimum to maximum yield anomalies range from 2 (t/ha) for soybean to 6 (t/ha) for rice. Maize and rice regional yield anomaly distributions are asymmetrical towards negative values in the topmost region (i.e., North America and Oceania respectively). Soybean yield anomaly distributions are more similar among regions and wheat asymmetry is linked to one region only with highest yield positive anomaly (i.e., Northern Europe).

Standard deviations of yield anomalies calculated per region per crop over the totality of the time period are also displayed (figures 1(e)–(h)). We find large inter-regional differences especially for rice (yield anomaly standard deviation ranging from 0.073–1.094 t/ha) but also for maize (0.047–0.638 t/ha). Standard deviation confidence intervals (calculated from 1000 bootstrap samples) indicate that regions with highest standard deviation are significantly different from the average standard deviation calculated over all regions (dotted line in figures 1(e)–(h)). Groups of regions that are significantly above average are clearly distinct from the ones that are below average with a third group of regions not being significantly different from average (this is one region for maize, four for rice and soybean and six for wheat). Although some regions do not significantly differ in terms of standard deviation, at least two groups of regions with low and high standard deviations can be identified for each crop.

Figure 2 illustrates the evolution of cultivated area proportions for regions totalizing at least 75% of global cultivated areas over the study period. The Simpson diversity index which was first developed to characterize species diversity [20] is here used to reflect the evolution of the spatial repartition of cultivated areas among world regions (see inserts in figure 2). The index indicates heterogeneity for high values (i.e., a concentration of production in a few regions) and homogeneity at low values (i.e., a more even distribution of hectares among regions). Highest values of the Simpson index are obtained for soybean and to a lesser extent, rice. It indicates that soybean production is the most spatially concentrated followed by rice with maize and wheat exhibiting comparable levels of heterogeneity in particular at the end the time period. While maize and rice Simpson indices remained somewhat constant during the study period, wheat and most notably soybean Simpson indices show a decreasing trend with some differences. Wheat index shows a monotonous decreasing trend followed by a stabilizing period. On the other hand, the spatial spread of soybean production appears to be marked by three consecutive periods.

Zoom In Zoom Out Reset image size

These patterns can be explained by the dynamics of area proportions (figures 2(a)–(d)). A close look at the Simpson index for soybean reveals a short increase until 1968 (reflecting an increase in the concentration of production in North America) followed by an overall decrease of the index during 1968–1993. The trend observed in 1968–1993 is due to the accession of South America to significant proportions of total cultivated areas (i.e., from less than 2% in 1960 to almost 30% in 1993, see figure 2(c)). After 1993, proportions of soybean areas allocated to Eastern Asia and to North America still follow a decreasing trend, and the soybean production becomes more and more concentrated in South America (about 45% of the total soybean area in 2010). But, although the replacement of Eastern Asia and North America by South American areas initiated a concentration of production, soybean production is less concentrated in 2012 than it was in 1961 when only two regions totalized above 90% of total areas (i.e., Eastern Asia and North America for about one fourth of actual surfaces).

American maize production occupies more than 20% of global maize areas over the totality of the time period. Towards the end of the period, Eastern Asia maize areas (China among three other countries) also reach about 20% of world maize areas and, together with Eastern Africa areas (ranging from Ethiopia to Madagascar) replacing South and Central America and Eastern Europe areas (figure 2(a)). The actual repartition of land among regions is fairly spread out with a total of seven regions producing 84% of global maize production on 78% of total areas, on average over the study period. Note that maize Simpson index has the lowest values of all crop species (i.e., values ranging between 0.11 and 0.13), suggesting that maize production is less spatially aggregated than the other crops.

Rice is predominantly produced in Asia; in Eastern Asia (about 40% of total production on average over the study period), Southern Asia (∼29%) and Southeastern Asia (∼23%). Rice production is spatially concentrated as revealed by the relatively high Simpson index (about 0.3 over the time period); on average, the three above-mentioned regions totalized above 90% of the production on above 90% of global cultivated areas (figure 2(b)). The most important region in terms of cultivated areas and production is Southern Asia covering a stable 40% of global areas throughout the period. A notable development is the increase -in particular after the 1990 s- in total and in proportion of South Eastern Asia rice areas (especially Thailand and Indonesia, among 11 countries). Eastern Asia areas remained somewhat constant but decreased in proportion (figure 2(b)).

Five regions totalize about 76% of global wheat areas over the study period. Wheat Simpson index decreases from about 0.2 to about 0.13 over the time period. Wheat production is thus more spread out than rice and soybean and about equivalent to maize. Wheat prominent feature is the collapse of Eastern Europe areas in 1992. Area proportions show an increase followed by a decrease in both North America and Eastern Asia. Southern Asia (India and eight other countries) is the only major producing region with a sharp increase of its cultivated area during the study period, representing about 23% of global surfaces at the end of the time period.

The order and values of global YIV components differ among crop species (figures 3(a)–(d)). In particular, maize YIV decomposition is very different from that of rice, soybean and wheat. In the following we focus on maize and rice because of their contrasted variance decomposition.

Zoom In Zoom Out Reset image size

Component 1,$\overline{{{\Gamma }_{1\bullet }}}$ (equation (6)), is bigger for rice (0.042) than maize (0.014). This is because of higher Simpson index values for rice compared to maize (i.e., rice is more than twice as spatially concentrated as maize). Also, average squared anomalies over all producing regions, ${{v}_{t}}$ is slightly bigger on average for rice than maize (mainly because of the 2008 rice yield anomaly in Oceania).

Component 2 describes how regional yield anomalies compare to the inter-regional average. This component is positive and small for maize and strongly negative for rice (figures 3(a)–(b)). Component 2 is negative (positive) when regions with large proportions of cultivated area show small (high) squared regional yield anomalies. For rice, the three major producing regions are located in Asia and all show a negative regional component ${{\bar{\Gamma }}_{2i}}$ (figure 3(b)) meaning that these three regions have smaller squared anomalies than average. For maize, North America is both the largest and the most variable producing region (figure 1(e)). Consequently, the value of the regional component ${{\bar{\Gamma }}_{2i}}$ (equation (5)) calculated for North America is positive and much higher that the regional components calculated for the other maize producing regions (figure 3(a)).

Component 3 $\overline{{{\Gamma }_{3\bullet }}}$, equal to the weighted sum of all yield regional covariances, is negative for maize (i.e., on average negative yield anomalies tend to be compensated by positive anomalies in other maize producing regions) and positive and small for rice.

Soybean and wheat yield variance components are similar (figures 3(c), (d)). Global variances V are lower than those of maize and rice. Component 2 is negative for both crops indicating that the majority of the production is located in regions with below-average squared anomalies, in particular North America for soybean (the most important region on average although outpaced by South America post 2000, see figure 2) and Southern Asia for wheat. Soybean regional yield anomalies tend to be negatively correlated thus decreasing global yield variance. Wheat component 3 values are very close to zero. Note that the relative importance of each component within the decomposition is unchanged when varying the time window of the classification (figure S1).

Figure 4 shows values of global variance and of its components in case of a uniform repartition of cultivated areas among regions. These results correspond to the second limit case described in section 2.2. In this case, the same proportion of cultivated area is allocated to each region, Simpson index is equal to its minimum value (i.e., 1/N), and component 2 is equal to zero. This limit case is theoretical as proportions of cultivated area in all regions cannot be equal in reality. Our underlying hypothesis is that a more even distribution of land will decrease global yield variance by distributing risk over a large number of regions. We use this limit case to determine how much global yield variance could be reduced or increased by using an even repartition of cultivated areas at the world scale.

Zoom In Zoom Out Reset image size

Results show that an even distribution of cultivated land decreases maize and soybean global variance to roughly half and a third of its original value respectively. It also decreases, but to a lesser extent, global wheat variance. On the contrary, rice global yield variance is increased in case of a uniform land use distribution (figure 4). These contrasted results are due to the fact that an allocation of 1/Nth of the total areas to all N regions has two effects on the decomposition. First, component 2 is automatically brought to zero (see methods section, limit case 2). Second, component 1 automatically decreases, because the Simpson index is at its lowest value. Component 1 is then proportional to inter-regional yield variance $\left( {{\Gamma }_{1,t}}=\frac{{{\nu }_{t}}}{N} \right)$.

Component 2 was initially highly negative for rice, and this negative value is set to zero when using of a uniform distribution of cultivated area. The use of a null value for component 2 instead of a negative one is only partially compensated by a decrease of component 1, and the global yield variance for rice is thus higher compared to its original value (figure 4(b)). For maize, component 2 was originally highly positive, and this positive value is set to zero when using of a uniform distribution of cultivated area. For this crop, the use of a null value for component 2 instead of a positive one leads to a strong decrease of the global maize yield variance (figure 4(a)). Although component 2 was initially negative for soybean and wheat, an even repartition of cultivated lands leads to a smaller global yield variance for these two crops. This is due to the fact that the reduction of component 1 is able to compensate the increase of component 2 for these two crops, resulting in an overall variance decrease (figures 4(c)–(d)).