Comparative assessment of environmental variables and machine learning algorithms for maize yield prediction in the US Midwest

The study area includes 12 states in U.S. Midwest: Illinois, Indiana, Iowa, Kansas, Michigan, Minnesota, Missouri, Nebraska, North Dakota, Ohio, South Dakota, and Wisconsin. We obtained county-level maize yield statistics from the National Agricultural Statistics Service (NASS) of USDA for 2001 to 2016. For counties with both rainfed and irrigated maize yield, we used the reported overall yield value. Counties that had less than one percent cultivated maize areas were discarded.

We compiled an extensive set of variables extracted from gridded weather data, optical, thermal, and microwave remotely sensed observations, land surface data assimilation systems, drought indices, gridded soil maps, as well as crop progress reports (table 1). These include time-series of land surface temperature (LST), surface reflectance, and VIs including EVI, Green Chlorophyll Index (GCI) (Gitelson 2005), and Normalized Difference Water Index (NDWI) (Gao 1996), derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite images (Schaaf and Wang 2015, Wan et al 2015); air temperature, precipitation, and VPD from the Parameter-elevation Relationships on Independent Slopes Model (PRISM) climate dataset (Daly et al 2008, 2015); surface soil moisture from the European Space Agency (ESA) Climate Change Initiative (CCI) dataset (Gruber et al 2019); rootzone soil moisture, potential ET, and ET from the Global Land Data Assimilation System (GLDAS-v2.1) dataset (Rodell et al 2004); and ET from the Atmosphere-Land Exchange Inverse (ALEXI) remote sensing model (Anderson et al 2011b). We derived a water stress (WS) indicator defined as the ratio between ET and PET from the GLDAS dataset. We incorporated three drought indices including the Palmer Drought Severity Index (PDSI) (Abatzoglou et al 2014), Standardized Precipitation Index (SPI) derived from PRISM data (Guttman 1999), and the Soil Moisture Deficient Index (SMDI) derived from the CCI soil moisture time series (Narasimhan and Srinivasan 2005). A heat index, Extreme Degree Day (EDD), was computed as the cumulative temperature above 30 °C during summer months (June–August) using the PRISM temperature data (Lobell et al 2013, Jin et al 2016). Soil properties such as soil available water content (AWC), organic matter content (SOM) and nitrogen content were derived from gridded soil datasets SoilGrids250 m (Hengl et al 2017) and SoilGrid + (Ramcharan et al 2018). We included county-level statistics such as maize harvested area, irrigated fraction, crop progress, and previous 5 year yield mean extracted from the USDA NASS. The crop progress data at state-level were converted to the dates when a certain percentage of acreage reaches a specific growth stage. Detailed description of each variable and its relation to yield is provided in table S1 of supporting information(stacks.iop.org/ERL/15/064005/mmedia).

These datasets were masked to include all agricultural areas using MODIS land cover data (Friedl and Sulla-Menashe 2018) and spatially averaged for each county. The sequential variables were extracted from March to October. We did not use a specific maize mask since such information is not always available within the growing season in the US.

We examined six machine learning algorithms including least absolute shrinkage and selection operator (LASSO) regression, Support Vector Regression (SVR), Random Forest (RF), eXtreme Gradient Boost (XGBoost), as well as two deep neural networks—LSTM and CNN. LASSO is a linear regression model that performs L1 regularization and variable selection in favor of simple and sparse models. SVR is a variant of the Support Vector Machine (SVM) that optimizes hyperplanes to separate training samples and uses kernels to incorporate non-linear relationships. Random forest is an ensemble machine learning algorithm that creates regression trees using subsets of features and bootstrap samples. XGBoost is an algorithm that implements an optimized gradient boost technique over decision/regression trees (Chen and Guestrin 2016). Random forest and XGBoost are the most widely used machine learning algorithms in Geoscience and other areas (Anderson and Lucas 2018, Fernández-Delgado et al 2019, Hengl et al 2017; Shangguan et al 2017). LSTM network is a recurrent neural network (RNN) where a sequence of inputs is processed by nodes with directed connections (Lecun et al 2015). CNN is a specialized class of neural networks mostly used in image recognition. Figure S1 shows the LSTM and CNN architectures used in our analyses. Detailed descriptions of the models and parameterization can be found in table S2 and text S1.

In exploratory analysis, we assessed the relationship between individual environmental variables and yield. The temporal dynamics of both Pearson's ($r$) and Spearman's (${r_s}$) correlation coefficient was analyzed for sequential variables at an 8-day time interval as well as for non-sequential variables. The general trends of yield response to environmental variables were established using Multivariate Adaptive Regression (MARS) models with the year as a secondary predictor. The full dataset (2001–2016) was used for the exploratory analysis, while a separate assessment on the subset (2001–2006) involved only in the training sets showed similar results.

Based on the correlation analysis, we constructed three sets of environmental variables to build the feature space. A 'Minimum (Min) set' contains variables that are commonly used in previous yield prediction studies, including EVI, NDWI, daily average air temperature, average daily VPD, precipitation, daytime and nighttime LST, as well as the year and geographic coordinates to model spatiotemporal correlations. A 'Good set' extends the minimum set with variables that show a substantial correlation with yield (Section 3.1). These include ET, water stress, SMDI drought index, EDD, AWC, previous 5 year average yield, harvested and dented progress, etc. A third level—'Full set'—includes all the variables. In addition, we tested the effect of temporal aggregations over sequential variables at 8-day, 16-day, or monthly time steps. These resulted in nine experimental levels based on three levels of variable selection and three levels of temporal aggregation, with the total number of variables ranging between 58 and 891. Detailed information is provided in table S3.

A comparison of machine learning models and feature levels was performed based on the end-of-season yield estimation. Models were trained to predict yields from 2007 to 2016 for each county using data prior to the testing year. Model performance was evaluated using Root Mean Squared Error (RMSE), Mean Absolute Percentage Error (MAPE), Mean Absolute Error (MAE), correlation coefficient (r), and bias. The best-performing machine learning algorithm was further evaluated for within-season yield prediction skill, which was made twice a month from June to October using data up until the predicting date. The first monthly forecast was early in the month (around 2nd–8th day) and the second was late in the month (around 18th to 24th day). Features not available before the forecasting date were excluded. The seasonal county-level yield estimates were aggregated to national-level and compared to the forecasts of WASDE (details in text S2).

Linear correlation between input variables and maize yields provide important clues regarding the significance of predictors. The temporal evolution of the correlation for sequential variables are shown in figure 1(a). Overall, almost all variables exhibit statistically significant correlation with the yield during the growing season (p-value below 0.001) (figure S2). Based on their similarities in the temporal relationships with yield, variables are grouped into three categories (figure S3). The first group includes surface reflectance and VIs representing the amount of green foliage in the crop canopy. These variables possess strong correlations with yield during the whole growing season, especially in summer. The second group indicates crop heat stress, including air temperature, LST, and VPD. The last group is comprised of other variables that are associated with crop water stress, including ET, soil moisture, water stress index, and drought index SMDI. Besides, there are also high inter-correlations between variables within these groups (figure S4).

Zoom In Zoom Out Reset image size

Overall, VIs (GCI, NDWI, EVI) and the NIR band from MODIS show the strongest signal (figure 1(a)). The correlation between yield and GCI/NDWI/NIR is between 0.72–0.74 from mid-July to mid-August, attributed to their high sensitivity to canopy green leaf volume, which directly contributes to grain growth. The negative correlations between NIR/VIs and yield (positive for the visible and SWIR bands) in April and May are consistent with previous studies (Johnson 2016) and related to the increasing trend of yield from low (earlier planting dates) to high latitudes (later planting dates). Heat stress variables, especially daytime LST, maximum and mean VPD, and maximum daily air temperature, have strong negative correlations with yield in summer consistent with previous findings (Lobell et al 2013). The correlation between yield and daytime LST, maximum VPD, and maximum air temperature in the last week of July are −0.53, −0.47, and −0.41, respectively. The superior performance of LST over air temperature agrees with previous studies, as canopy temperatures incorporate the effect of soil water moisture and evaporative cooling, and are thus better indicators for heat stress and drought impact (Siebert et al 2014, Pede et al 2019). Among the variables related to water stress, ALEXI ET exhibits the highest correlation between mid-July and late-August (r ∼ 0.38–0.42). Among the drought indices, SMDI based on CCI microwave surface soil moisture explain more yield variation than SPI and PDSI (maximum r 0.30 vs. 0.28). But SMDI based on GLDAS root zone soil moisture shows only weak correlations with yield (maximum r = 0.22). Other variables, including nighttime LST, minimum VPD, potential ET, rainfall, and root zone soil moisture, show only weak correlations. Note that the correlation between rainfall and yield (maximum r = 0.26) is higher when rainfall is aggregated over a longer period (i.e. month) rather than 8-day (figure S5) since the effect of rainfall on yield is often cumulative.

Most non-sequential variables are highly correlated with yield (figure 1(b)). The previous five-year yield average has the highest correlation (r = 0.67) with yield, followed by the percentage of county area planted with maize and total harvested area. EDD negatively impacts yield with r = −0.46, stronger than the correlation with LST. The correlation between AWC and yield (r = 0.23) is also significant. Most crop progress variables show statistically significant correlations with yield. The dates of harvested and dented stage both exhibit high correlations with yield at r ∼ 0.24–0.28, indicating that longer growing season favors crop yield (figure 1). The positive correlation between planting/emergence/silking and yield is a manifest of the relationship between yield and latitude (Long et al 2017) since we did not find any significant correlation when fixing the effect of latitude (figure S6). A separate assessment using Spearman's correlation coefficient, which considers non-linear relationships, shows similar results for both sequential and non-sequential variables (figure S7).

Figure 1(c) shows the yield response curves of selected variables by fixing the year in 2008 (full result in figure S8). Most variables exhibit expected behaviors. The response curve of rainfall suggests that both drought conditions and excessive rainfall are harmful to yield (Lobell et al 2014, Li et al 2019a). High maximum temperature of July harms yield, especially when it exceeds 33.5 ºC. A linear reduction in yield is found with increasing extreme heat events (i.e. EDD). In addition, we also found increasing trends in the correlation between heat and water stress variables and yield since 2001 (figure S9), suggesting increasing sensitivity of yield to extreme drought and heat conditions (Hawkins et al 2013, Lobell et al 2014).

All machine learning models show reasonable accuracy in predicting end-of-season county-level maize yield (figures 2, S10). Overall, Lasso has the lowest accuracy for all feature levels (figure S11; p < 0.01), with average RMSE (2007–2016) between 16.7 and 24 bu ac⁻¹ and a large year-to-year fluctuation. SVR performs better than Lasso but not as good as the other models. XGBoost has the best overall performance (figure S11; p < 0.01) and the smallest year-to-year variation. It achieves the highest accuracy with the Good set and 16-day aggregation: RMSE = 14.8 bu ac⁻¹, MAE = 11.4 bu ac⁻¹, and MAPE = 9%. The superiority of XGBoost over RF is attributed to the boosting technique that iteratively reduces bias and variance. In additional, XGBoost also runs significantly faster than SVR and RF given the same computing resources (table S4). Interestingly, we did not find any significant difference between the performance of deep neural networks and RF and SVM (figure S11). LSTM only marginally outperforms RF and CNN (not statistically significant, figure S11), but it is not as good as XGBoost especially when a small number of features (monthly or Min set) was involved (figure 2). CNN and RNN are designed to capture specific patterns in structured data: CNN specializes in recognizing spatial locality in images and RNN in detecting sequential dependencies. However, the tabular style dataset used in this analysis has no spatial-related patterns for the convolutional kernels of CNN, although the sequential features might have favored LSTM, explaining its slightly better performance over CNN. This observation is consistent with previous studies (Wang et al 2017, Fernández-Delgado et al 2019, Zhang et al 2019).

Zoom In Zoom Out Reset image size

The yield prediction accuracy using the Good and Full sets is significantly higher than that of the Min set across all models, while the difference between the Good and Full sets is minimal (figures 2, S11). Additional information contained in the Good/Full set, such as previous yield records, ET, soil moisture, EDD, soil properties, and crop progress reduces the MAPE by more than 1% for RF, XGBoost, CNN, LSTM and 2%–5% for Lasso and SVR (figure 2(b)). The performance of non-neural-network algorithms does not always improve as sequential features are aggregated at smaller time steps (figures 2(b), S11). For example, RF achieves the lowest accuracy using 8-day aggregated observations. While frequent observations contain more information during rapid growth periods, aggregation over a longer time period reduces random errors and the frequency of missing data. On the contrary, the accuracy of CNN and LSTM increases with smaller temporal aggregation steps and more variables (figure 2(b)), indicating the ability of these deep neural networks to work on large feature spaces. Nevertheless, there is not a significant overall effect of temporal aggregation frequency (figure S11).

Feature importance was evaluated Using RF and XGBoost with the Full set and monthly aggregation. NDWI, GCI, and NIR reflectance in July and August were found to be the most important features for yield prediction (figure S12). Other important features include previous 5 year yield average, harvest area, percentage of irrigated area, potential ET, drought conditions (PDSI, SMDI), June and July average VPD, as well as soil properties (organic matter content, available water content). XGBoost also leverages the information from crop progress data such as the time of emergence and silking. The feature importance rankings are different each year, as the controlling environmental factors of yield vary (figures S13, S14).

Within-season yield forecasting was evaluated using XGBoost, 16-day temporal aggregation of features, and three variable selection levels (figure 3). Using only information from March to May, the early June prediction achieves a reasonable estimation accuracy with 23.9 bu ac⁻¹ RMSE (16% MAPE) based on the Good or Full sets, attributed to the strong correlations between yield and MODIS reflectance and ALEXI ET from March and May (figure 3). As more information becomes available in the summer (June–August), the prediction accuracy gradually improves (figure 3(a)). The most significant improvement occurs in July and early August when maize successively reaches the peak vegetative stages. Beyond late August, the prediction errors stabilize. Compared to the Min set, the Good and Full feature sets improve the prediction in all forecasting dates, especially in the early-season (inset in figure 3(a)). Again, there is not a significant difference in the performance between the Good and Full sets. In addition, we tested a remote sensing Evaporative Stress Index (ESI), defined as the climatological anomalies in the ratio between actual ET and potential ET (Anderson et al 2007, 2011b). Compared to other ET related variables, such as ALEXI ET and GLDAS water stress, ESI shows a higher correlation with yield early in the season (i.e. before June) (figure S15). However, adding ESI to the Good sets does not significantly improve the within-season yield forecasting (figure S15).

Zoom In Zoom Out Reset image size

The evolution of the forecasting accuracy during the growing season has significant spatiotemporal patterns (figures 3(b), S16,S17). The benefit of summer information is the most pronounced in southern states, such as Kansas, Missouri, southern Illinois, and Indiana (figures 3(b), S16). The prediction accuracy is higher in northern states, including North Dakota, Minnesota, Wisconsin, Iowa, and Michigan (figure 3(b)), where more than half of the counties could achieve RMSE below 10 bu ac⁻¹ by early September. While almost all the years see a decrease of RMSE from June to August, the decline is especially dramatic in 2012, when a summer drought significantly depressed yields (Grigg 2014) (figure S17). As a result, early-season forecasting in 2012 is inefficient with a significant overestimation of yield (figure S18). Although the yield forecast skill increases when late June and July data are provided, the late-season prediction accuracy of 2012 is still much lower compared to other years, especially for counties with significant irrigation presence (figure S18). This indicates that the data-driven models might not extrapolate well to unseen climatic conditions, and thus could lead to undermined performance.

The within-season yield forecasts produced by the Good set was further compared to that of the USDA WASDE estimates at national level. XGBoost significantly outperforms that of the WASDE in the early season (June and July) by 2 to 2.5 bu ac⁻¹ in RMSE, especially in extreme years such as 2012 (figure S19). XGBoost results are comparative to WASDE in August, and only slightly worse in September. Note that starting from August, the WASDE was primarily driven by surveys consisting of growers' assessment and field measurements, and thus became much more accurate.