Overview

The time period of 1965–2010 over which yield trends were evaluated is the most relevant for use in estimating historical rates of yield gain to inform projections of future food security because it represents a contemporary period when public and private sectors in many developed and developing countries have invested heavily in application of modern genetics, agronomy and supporting basic sciences and information technologies to develop and gain adoption of technological advances in crop production.

Evaluated mathematical models were as follows: linear (L), quadratic with upper plateau (QP), linear piecewise (PW), linear with upper (LUP) or lower plateau (LLP), and compound-rate exponential (EXP) (Fig. 3). Except for the United States and Africa, where yield trends were analysed at the subnational and regional levels, respectively, the present study focuses on average yield at the national level because this is the spatial scale at which most projections of food security are made. Criteria for the selection of the best-fit trend are based on the assumption that the most appropriate statistical model is one that can describe the observed time trend in yield with the lowest error and minimal bias in distribution of residuals, and which uses the smallest possible number of estimated parameters while meeting the assumptions of normality, independence and homogenous variance in regression analysis37 (Fig. 4).

Analysis of trends in cropland harvested area

Long-term (1965–2011) data in crop-harvested area were retrieved from FAOSTAT and used to analyse changes in cropland area of staple crops (cereal, oil, sugar, fibre, pulses, tuber and root crops) and three major cereal crops (rice, wheat and maize; FAOSTAT Database–Agricultural Production ( http://faostat.fao.org/). A three-phase linear model was fitted to the observed trends using the PROC NONLIN procedure in SAS Version 9.1 (SAS Institute, 9.1 Foundation for 64-bit Microsoft® Windows®, SAS Institute, Cary, NC):

where y is crop-harvested area (Mha), x is year and x 1, 2 are the breakpoint years. The fitted trilinear models have coefficient of determination (r2) of 0.97 and 0.93 (root mean square error (r.m.s.e.)=9.5 and 6.9 Mha) for area of staple crops and three major cereal crops, respectively. A total of 47 years of yield data were used in the regression analysis (1965–2011 time period). All estimated parameters were significant (Student’s t-test, P<0.01), except for parameter c in the fitted trend for cropland area of the three major cereal crops (Student’s t-test, P=0.20). An initial period of increase in harvested area occurred until 1980, followed by a period of little or no increase in harvested area of staple and major cereal crops, respectively, which lasted until early 2000s (Fig. 1). This period was followed by an unprecedented rate of expansion in harvested cropland area during the last decade, as indicated by the statistically significant higher rates of increase in harvested area of staple and major cereal crops during the last decade (parameter d in (equation 1)) compared with earlier rates of increase during the first two decades of the green revolution (parameters b and c in (equation 1); Student’s t-test; P<0.01).

Change in crop-harvested area during the last 10 years (2002–2011) was calculated for staple crops and, separately, for four crops (rice, wheat, maize and soybean) that accounted, altogether, for 83% of observed increase in staple crop area:

where net and relative changes are the absolute (Mha) and relative (%) differences in harvested area between two periods (Mha), and HA 2002–2003 and HA 2010–2011 are the calculated 2-year average harvested area for the 2002–2003 and 2010–2011 intervals, respectively. A net change of 85 Mha was observed for staple crops between 2002–2003 and 2010–2011 intervals, indicating that an 8% increase in global cropland area has occurred in only 10 years (Supplementary Table S1). This change was mostly accounted for by increased harvested area of maize and rice in Africa (+9 Mha), soybean in South America (+15 Mha) and maize, rice and soybean in South and Southeastern Asia (+16 Mha). In addition, remarkable was the little relative change of harvested crop areas (negative in some cases) observed for North and South Africa, West Asia, Europe and North America (Supplementary Table S1).

Published projections in grain yields of cereal crop yields

Previously published projections on yield trajectories, based on exponential compound rates, were retrieved from the literature12,13,14,15,16,17,18 (Supplementary Table S2). In many of these projections, countries were grouped into categories according to their geographic location or degree of economic development, and a crop-specific exponential yield gain rate was assumed for each category. Other studies simplified this approach by assuming a single, worldwide exponential yield gain rate for each crop. A common feature of previous studies is the lack of recognition of a biophysical limit on crop yields determined by solar radiation, temperature and water supply from both rainfall and irrigation. To illustrate the discrepancy between reported projections and observed historical trends, reported projected trajectories for the US average maize yield were plotted and compared against the projected yields based on the historical (1965–2011) yield trend (Fig. 2).

Data on grain yield data for cereal crops

Average annual yield data from 1965 to 2010 were retrieved for rice, wheat and maize in selected countries and regions, resulting in a total of 36 (crop–country or –region) cases that include a wide range of production environments and yield levels (FAOSTAT Database–Agricultural Production ( http://faostat.fao.org/; National Agricultural Statistics Service—Crops US state and county databases ( http://www.nass.usda.gov/index.asp; Supplementary Table S3). Available wheat and maize yield data in the USA also includes 2011 National Agricultural Statistics Service—Crops US state and county databases ( http://www.nass.usda.gov/index.asp). Hence, a total of 46 years of yield data were used for the yield trend analysis, except for the United States where a total of 47 years of yield data were used. In the case of the United States, trends were analysed for major producing regions: rice in California and south-central, wheat in the southern, central and northern Great Plains, and maize in the eastern and western Corn Belt. These regions accounted for 77, 47 and 85% of total US production of rice, wheat and maize, respectively. Separate analyses were performed for rainfed and irrigated maize yields in the western US Corn Belt, where both irrigated and rainfed production are important (54 and 46% of total maize production in western US Corn Belt, respectively). The analysis of maize yield trends in Africa was conducted for three regions: East, Central and West Africa. The average yield for a given region in the United States or Africa was calculated as the sum of total crop production in all the states (the United States) or countries (Africa) within the region divided by total harvested area.

Statistical analysis of yield trends

Data on average annual grain yield were plotted against year for each crop–region case. Six statistical models, extensively used in the literature to describe yield trajectories, were evaluated for their performance of fitting trends in grain yields since the onset of the green revolution (Fig. 3):

where y is grain yield (kg ha−1), x is year, x i is the initial year of the time series (1965 in the present study) and y 0 is the yield plateau level (kg ha−1).

SAS Version 9.1 programmes and procedures were used for all statistical analyses. The six models were fitted to observed yields using the PROC MIXED, PROC REG and PROC NONLIN procedures. Estimates of model coefficients (and associated confidence intervals), coefficient of determination (r2) and the r.m.s.e. were calculated for each crop–region model combination (Supplementary Tables S4 and S7). In the case of nonlinear models, r2 was calculated as the square correlation between predicted and observed values38. Violations of assumptions of regression analysis were identified by performing Shapiro–Wilks test for normality, Levene’s test for variance homogeneity and the Durbin–Watson’s (D) test for serial correlation39 (Supplementary Table S8). Residuals were plotted against year and the significance of linear- and quadratic models fitted to these plots was evaluated to detect potential biases along the time series (Figs 6, 7, 8). Finally, the Box–Cox method was used to identify suitable power transformations that can improve models performance40 (Supplementary Table S9).

The following criteria were used to identify the best-fit model for each crop–region yield trend: highly significant (Student’s t-test; P<0.01) model parameters, highest r2 and lowest r.m.s.e. compared with other models, and independent, normally distributed yield residuals with homogenous variance and unbiased distribution when residuals were plotted against year (Fig. 4). An additional criteria was that none of the applied data transformations improve model fit by >5% based on comparison of calculated r2 among models based on transformed and untransformed data37 (Fig. 4). In the process of selecting the best-fit model, the PW and LUP or LLP models were mutually exclusive, that is, the LUP or LLP models were discarded if estimated parameters of the PW model were all significant and no bias was detected in the distribution of yield residuals. In 22 out of the 36 cases, one statistical model clearly outperformed the others by a difference in r.m.s.e. ≥5% and this was the chosen best-fit model (Supplementary Tables S4 and S7). When differences in r.m.s.e. among two or more models were <5% (14 out of 36 cases), the two models with lowest r.m.s.e. were selected as best-fit models, unless there was a justification to discard one of them, for example, failure to meet assumptions of regression analysis (EXP in Canada) or when the model with higher r.m.s.e. also had the larger number of parameters (QP for wheat in the Netherlands and maize in Italy).

Selected best-fit models exhibited highly significant coefficients (Student’s t-test; P<0.01) and higher goodness of fit compared with other fitted models (higher r2 and smaller r.m.s.e.; Supplementary Tables S4 and S7). Out of 45 selected best-fit models, only 1 and 4 did not meet assumptions of variance homogeneity (LLP model for rice in Vietnam) and normality (L and EXP models for maize in the United States), respectively (Levene and Shapiro–Wilks tests; P<0.01) (Supplementary Table S8). In contrast, ten selected best-fit models exhibited significant positive serial autocorrelation: rice in Bangladesh, China, Indonesia, Republic of Korea, Philippines, Vietnam, wheat in China and maize in central and west Africa and Italy (D-test; P<0.01; Supplementary Table S8). Residual plots of the best-fit models did not exhibit any obvious trends over time, except for four cases with severe positive serial autocorrelation: rice and wheat in China and rice in Indonesia and Philippines (0.45<D<0.80), which are possibly related to vigorous public sector varietal improvement and agronomic research programmes that promote rapid adoption of improved rice varieties and associated fertilizer and pest management practices and/or the frequency of the yield survey (Figs 6, 7, 8 and Supplementary Table S8). Serial correlation affects the estimated variances but does not affect the value of the estimated model coefficients, which results in the estimators looking more accurate than they actually are38. To test how serial correlation may have affected the significance of the estimators, we assumed the estimated variance to be 60% of the true variance, which is expected for a first-order serial correlation coefficient (AR(1)) of 0.5 (ref. 39), similar to the estimated AR(1) of 0.6 found for the above four crop–country cases. Results indicate that the parameter estimates of the best-fit models were still highly significant despite serial correlation (Student’s t-test; P<0.01). EXP models exhibited a remarkable bias on their residuals plots in 75% of the 36 cases (with significant positive autocorrelation in 58% of the cases) and did not meet assumptions of regression analysis in 31% of the cases (Figs 6, 7, 8 and Supplementary Table S8).

Other evidence of the robustness of selected best-fit models was that their fit, based on the original untransformed data, was not improved after data transformation using the Box–Cox method40, except for rice in Republic of Korea and maize in east Africa (Supplementary Table S9). In both cases, model fitness increased slightly after applying a reciprocal transformation to the yield data (r2=0.90 versus 0.84 (Republic of Korea) and 0.57 versus 0.49 (East Africa), with and without data transformation, respectively). None of the non-selected models fit to the transformed data outperformed the fit of the selected models shown in Table 1. And although the presence of abnormally low- or high-yield years during the last years of the time series can potentially hinter the identification of yield plateaus, this is not a large concern in the present study because all cases that exhibited an upper yield plateau have ≥10 years of yield data after the breakpoint year, except for maize in France with 7 years. Perhaps, more important, the 99% confidence interval of the breakpoint year (x 0 ) was within the 1965–2010 interval in all cases. In addition, visual inspection of the crop–country trends with statistically significant yield plateaus did not have unusually high or low yields around the breakpoint year. Therefore, it is unlikely that the upper yield plateaus identified in the present study are the consequence of a few yield outliers at the end of the time series. The presence of autocorrelation and changes in variance over time can also hinder the identification of yield plateaus29. However, out of the 21 selected LUP or LLP models in the present study, variance was not homogenous in only one case (rice in Vietnam) and autocorrelation was significant in six cases (rice in China, Republic of Korea and Vietnam, and maize in central and west Africa, and Italy) (Supplementary Table S8). And even in the six cases with significant autocorrelation or non-homogenous variance, visual inspection of the residuals clearly indicated that LUP or LLP models outperform linear models (Figs 6, 7, 8). Although the yield residual versus time relationship of the linear model exhibited a strong trend, there was no detectable pattern in the yield residuals of the LUP or LLP models, with the only exception of maize in West Africa where there was no detectable trend in both linear and LLP models.

Yield trends from cases in which there was evidence of yield plateaus were re-analysed to determine the number of years after the breakpoint year (x 0 ) that were needed to identify a yield plateau (hereafter called x n ). The LUP model was successively fitted to yield trends in which data-years after x 0 were removed and then added one by one until estimated x 0 became statistically significant (Student’s t-test; P<0.01). Values of x n varied across cases, depending on the magnitude of year-to-year variation in yields along the entire trend line (Supplementary Table S6). In fact, there was a negative relationship between x n and r2 of the fitted LUP model to each crop–country case (y=20−14.9 x; r2=0.63; Student’s t-test; P<0.01). A conservative value of x n =7 years can be taken as the minimum number of years needed, after the breakpoint year, to detect a statistically significant yield plateau in production systems that exhibited high r2 (>0.85) such as rice in China and California, wheat in northwestern Europe and maize in Italy and France. In production systems with larger variability in yield along the yield trends, a greater number of years is required, ranging from 9 years in Japan (r2=0.63) to 17 years in East Africa (r2=0.49). The total number of years with 'flat' yields required to detect a statistically significant yield plateau, including the breakpoint year, was calculated as x n +1 and this is the value reported in the Results section.

Calculation of absolute and relative yield gain rates

Absolute yield gain rates (kg ha−1 per year) were calculated from the first derivative of the best-fit model, for each crop–region case, at three points in the time interval: 1970, 1990 and 2010 (Supplementary Table S5). For those cases in which there were two best-fit models, the model with lowest r.m.s.e. was chosen to describe the yield trend. Despite the slightly better performance of the EXP model, the L model was chosen to describe trends in maize yield in the United States because fit of the EXP and L models was very similar (difference in r.m.s.e. ≤2%), there was no evidence of changes in the absolute yield gain rate between 1965 to 2011 as tested using piecewise with one, two or three breakpoint years and the EXP model underperformed compared with other models in all the other cases. Moreover, using an exponential model would imply predicted average maize yields by 2030 that are 35 and 37% higher than the 2010 yield for maize in the eastern and (irrigated) western US Corn Belt, which will require average (2010–2030) yield gain rates of 180 and 220 kg ha−1 per year. These rates are 50 and 70% higher than the observed yield gain rates during the 1965–2010 interval. This is unlikely to occur given the slow down in yield gain rates observed in other intensive cropping systems and lack of increase observed for irrigated maize in contest winners and farmers’ fields in the western US Corn Belt3,19. In fact, although no yield plateau was detected for irrigated maize in the western US Corn Belt for the 1965–2011 interval, the linear regression for the last 10-year period (2002–2011) has a slope indistinguishable from zero (Student’s t-test; P=0.26, n=10).

Relative rates of yield gain (% per year) were calculated for 1970, 1990 and 2010 as the ratio between the absolute yield gain rate and the yield-trend-predicted yield in the particular year, expressed as percentage (Supplementary Table S5). Exponential compound rates (% per year) were also derived for each crop–region case from the parameter b in (equation 9), expressed as percentage. The objective was to highlight the limitations of using compound rates of yield gain to predict future yield trajectories: out of a total of 29 cases that have compound rate >1% per year, there was evidence of yield plateaus or decreasing yield gain rates in 41% of the cases, decreasing relative rates of yield gain during the entire time series in 66% of the cases and decrease in the relative rate of yield gain in the time interval between 1990 and 2010 in all cases.