|
A Comparison of Yield Prediction Approaches Using Long-Term Multi-Crop Site-Specific Data
Joaquin Casanova1,*, Garett Heineck1, David Huggins1
Published in Journal of the ASABE 67(3): 601-615 (doi: 10.13031/ja.15216). 2024 American Society of Agricultural and Biological Engineers.
1 USDA ARS, Pullman, Washington, USA.
* Correspondence: joaquin.casanova@usda.gov
The authors have paid for open access for this article. This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License https://creative commons.org/licenses/by-nc-nd/4.0/
Submitted for review on 1 June 2022 as manuscript number NRES 15216; approved for publication as a Research Article and as part of the Artificial Intelligence Applied to Agricultural and Food Systems Collection by Associate Editor Dr. Jonathan Aguilar and Community Editor Dr. Yiannis Ampatzidis of the Natural Resources & Environmental Systems Community of ASABE on 18 January 2024.
Mention of company or trade names is for description only and does not imply endorsement by the USDA. The USDA is an equal opportunity provider and employer.
Highlights
Abstract. Growers in the inland Pacific Northwest face numerous challenges in managing cropping systems. Climate variability, soil degradation, and topography all lead to significant spatial and temporal variability in yield. Often, yield modeling approaches such as deep learning can be “black boxes” or suffer from parameter uncertainty and instability, as with process-based crop models. To explore an alternative, we examined nearly two decades of crop rotation data from the R.J. Cook Agronomy Farm, along with soil properties, topography, weather, and multispectral data. We then tested two modeling approaches to estimate yield: linear modeling (LM) and Bayesian hierarchical modeling (BHM). We found BHM with spatial and temporal random effects performed best in predicting relative yield, both using soil variables as predictors or remotely sensed data. Since the BHM approach handles missing data, offers the possibility of farmer knowledge to be incorporated into prior probabilities, and gives uncertainty, this methodology lends itself well to decision support tools and on-farm study design.
Growers throughout the inland Pacific Northwest face unique challenges that are linked to the remarkable precipitation gradient and variation in soils found throughout the region (Jennings et al., 1990). Declining soil health, highly variable topography (Yang et al., 1998), and pressures of climate change (Kaur et al., 2017) will necessitate changes in management and agroecological shifts (Karimi et al., 2017) and increased variability in crop yields. Clearly, it would be useful to have easily accessible decision support tools that incorporate management decisions such as crop rotations, fertilization regimes, simple climatic variables, and topography, that also provide estimates of uncertainty. Such a tool would enable prescription mapping for fertilizer and testing of different rotations and allow for the optimization of yield under constraints of risk tolerance for any point along the precipitation gradient. Modeling will be a useful tool in sustainably intensifying agriculture in this region. In addition, knowledge of uncertainty in productivity across a field could allow researchers to improve their sampling strategies in on-farm research, where data and access are not always available (Lacoste et al., 2021).
There is a long history of crop modeling efforts that can predict crop development and yield over large spatial scales. Process-based models, such as CropSyst (Stöckle et al., 2003), DSSAT (Jones et al., 2003), APSIM (Keating et al., 2003), WOFOST (Van Diepen et al., 1989), and AquaCrop (Steduto et al., 2009), mechanistically simulate the growth of the crop and soil processes, often requiring multiple input files of meteorological data, soil properties, and historical crop management. These models can be applied to large spatial scales and incorporate remotely-sensed data (Kaur et al., 2022). These models simulate crop growth using a suite of parameters determining phenological development and biomass assimilation according to thermal time, radiation, and abiotic stressors. Mechanistic models are widely used and well tested, but the large number of parameters means there is often a steep learning curve, and the parameters themselves are unstable and have high uncertainty. This has been shown through Monte Carlo simulation and evaluation of different parameterization schemes (Lamsal et al., 2018; Kawakita et al., 2020).
Machine learning provides nonlinear and empirical modeling methods for yield estimation (Van Klompenburg et al., 2020). There are many techniques within this field (Duda and Hart, 2006). Random forest regression aggregates many regression trees, which split the data into progressively smaller groups based on predictor variables until the final groupings, where different regressions are applied to generate the predicted values. This method has been trained and tested using Australia-wide field trials (Newman and Furbank, 2021). Support vector machines find an optimal hypersurface between data points of different classes; this can be adapted to regression and has been used in rice yield modeling (Ghandi et al., 2016). Artificial neural networks take multidimensional data points and apply nonlinear functions to each dimension at “nodes,” then weight and sum the result before passing to the next layer of nodes until the final output node. Deep learning approaches, which are artificial neural networks with a large number of layers and nodes, have been used for yield prediction on a large and diverse maize dataset (Khaki and Wang, 2019). These techniques have been compared for wheat biomass estimation in Wang et al. (2016). A major criticism of machine learning or deep learning approaches is their black-box nature (Castelvecchi, 2016); there have also been reported difficulties in understanding the impact of individual inputs on the model (Khaki and Wang, 2019).
Linear models (LM) have been used as a purely statistical way to predict yield, as in Lobell and Burke (2010), where climatic variables are used as predictors for yield across a region and perform well compared to DSSAT. Similarly, Peng et al. (2016) modeled forage yield in South Korea. An extension of simple linear models is to incorporate the non-independence of the data through fixed and random effects, which have a hierarchical structure. This random error structure allows the estimation of error variance between and within different categories in the data. For instance, regions can be analyzed for yield variability without predictor variables (Blasch et al., 2020). Random effects can also incorporate spatial and temporal relationships (McCulloch and Searle, 2004). Often, yield measured at one location in a field is autocorrelated with nearby spatial locations and with measurements from preceding and following years. When using modern mixed model software, there are readily available functions that allow for the decay of correlation over distance in both space and time (Pinheiro et al., 2017). Spatial random effects processes can be incorporated, as in Chandra (2013), where district-level yields in Uttar Pradesh, India, were gap filled with a linear model with spatial random effects following a conditional autoregressive model (CAR(1)).
In the Bayesian view of LM, parameter estimation (e.g., slopes and intercepts) is based on the Bayes theorem, where posterior distributions for the parameters are found by updating prior distributions by the observations (Broemeling, 2017). Bayesian hierarchical models (BHM) extend this by considering multiple levels of distributions among observations and parameters (Kass and Steffey, 1989). More salient to agroecology, a spatiotemporal process can be made implicit within this hierarchy, with its own parameters to estimate (Wikle et al., 1998; Wikle et al., 2019). For example, this modeling framework has been used for modeling invasive species spread (Hooten et al., 2007) and frost risk, with the outcomes indicating that BHM has resulted in higher prediction accuracy over LM (Crimp et al., 2015).
In a basic spatio-temporal BHM, the measured variable z is modeled as a known form of distribution conditional on a latent variable y and parameters ?1. Next in the hierarchy, y is a variable conditional on parameters ?2. These parameters describe the distributions of any random effects as well as the coefficients of a linear predictor. Finally, the parameters ?1 and ?2follow a known distribution, p.
(1)
The latent field y could be a linear combination of predictor variables, x, plus a random effect that could be spatially or temporally correlated. In this way, the variable of interest is modeled hierarchically. In terms of linear equations, this looks like:
(2)
where
x = covariates
e = independent and identically distributed (iid) measurement error
? = spatiotemporally-correlated process noise
?1 and ?2comprise a vector of parameters describing the prior distributions for these noise terms and the fixed effects ß. The specification of these priors for ?1 and ?2gives this model the additional benefit of incorporating experiential knowledge derived from farmers. As an example, knowing qualitatively that low-elevation portions of the field have greater variability could be incorporated by increasing the uncertainty on the prior for the slope of the elevation covariate. BHM also permits gap-filling: by modeling missing data with BHM, the uncertainty associated with these data can be incorporated into the top-level estimates in the hierarchy.
To approximate the full distributions in a BHM, the Integrated Nested Laplace Approximation (INLA) uses analytical functions from the exponential family for the variables using the Laplace approximation (Martins et al., 2012). The approximations in INLA are designed to speed computations for large datasets as an alternative to Markov-chain Monte Carlo. Examples include crop modeling, crop insurance pricing in Brazil (Ozaki et al., 2008), and the effects of landscape position on biomass crop yield (Theleman et al., 2010).
The objectives of this paper are threefold. First, fit and compare the results from two classes of yield models satisfying the requirements of simplicity and explicability (LM and BHM) with two sorts of data (soil chemistry or remote sensing). Second, examine for the statistical models the impact of different spatiotemporal random effect structures/processes. Finally, we also evaluate the effects of BHM gap filling compared to more conventional interpolation. The goal is to develop a simple yield modeling framework that works from commonly available data, incorporates spatiotemporal effects, and provides uncertainty estimates.
Materials and Methods
This research uses data from years of crop rotation data at R. J. Cook Agronomy Farm (Huggins, 2015) as well as imagery courtesy of the U.S. Geological Survey from the Landsat 5 and 8 missions. The first part of this section describes the long-term research at this site, and the second describes the data processing and analysis done in R 4.1.2 (R Core Team, 2013).
Field Experiments at Cook Agronomy Farm
The R.J. Cook Agronomy Farm is a 37-ha USDA Long-Term Agroecosystem Research (LTAR) site, with varied topography, including Thatuna, Palouse, and Naff soil types (NRCS, 2013), located near Pullman, Washington. 369 georeferenced points on a non-aligned grid were established in 1999 (fig. 1). These sample points were used for harvest sampling every growing season, and a subset of 177 were chosen for soil sampling (153 cm deep by 4.1cm in diameter cores) in 1999, 2008, and 2015. A series of crop rotation experiments have been conducted there under no-till management. Crop species grown include winter wheat, spring wheat, canola, barley, lentils, winter pea, and alfalfa. It is important to note the complexity of these rotational species in modeling, as some crops are taken for grain yield (wheat and peas) and some for forage (alfalfa). These differences make conceptualizing “yield” difficult in a modeling context. More information can be found in Huggins and Uberuaga (2010) and Huggins (2015), but a summary of the spatially extensive, long-term dataset is given here. Table 1 gives a summary of variables from the Cook Agronomy Farm dataset considered in this study.
 |
| Figure 1. R.J. Cook Agronomy Farm, East field (46°47'N, 117°5'W), Northeast of Pullman, WA. This field represents a 37-ha watershed used for long-term agroecological research. Black dots are the 369 yield sampling locations (yearly), a portion of which are also soil sampling locations (1999, 2008, 2015). Contours are 2m elevation derived from the DEM. |
Harvest samples were taken at each site in a 2 m2 area and separated into grain and residue portions. Air-dried and oven-dried masses were measured, and finely-milled samples of grain and residue were analyzed for carbon and nitrogen content.
Soil samples were taken with a Giddings probe and divided into 10 cm depth increments for analysis. Profiles were described by layer type and morphological characteristics, and bulk density was measured. Subsamples were dried at 55°C and dry combusted for total carbon and nitrogen measurement with a CNS Leco analyzer. The pH of dried, sieved soil samples was measured with a 1:1 soil and water slurry with a Denver Instrument model 250 pH ISE conductivity benchtop meter and an Accumet #13-62-631 saturated KCl-filled glass electrode (McFarland, 2016).
A Trimble RTK 4400 was used to collect high-accuracy elevation data. The results were processed in the System for Automated Geoscientific Analyses (SAGA; Conrad et al., 2015), using local polynomial interpolation to produce digital elevation models at a range of raster sizes. From this, other topographic variables were derived, including elevation, aspect, transformed aspect (TRASP), slope, analytical hillshade, and topographic wetness index.
| Table 1. Variables comprising the Cook Agronomy Farm dataset. | ||
| Type | Variables | |
| Crop | Crop (crop type), GrainYield (dry grain per area), GYZ (Grain Yield, Z-normalized) | |
| Management | Nf, Kf, Pf (applied nitrogen, phosphorous, and potassium per area) | |
| Topography | Aspect, TRASP (transformed aspect), Face, Elevation, Slope, AnalyticalHillshade, TopographicalWetnessIndex, ConvergenceIndex | |
| Meteorology | Tave (season average daily temperature), Precip (total precip in calendar year before harvest) | |
| Soil | C (profile mass carbon per area), H (profile mass hydrogen ions per area), N (profile mass nitrogen per area) | |
| Remote sensing | Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), Plant Senescence Response Index (PSRI), Green Red Vegetation Index (GRVI), Green Normalized Difference Vegetation Index (GNDVI), Normalized Difference Tillage Index (NDTI) | |
Data Analysis
The process of data aggregation, processing, modeling, and model evaluation is shown in figure 2.
Preparation
Several datasets were used for this analysis: meteorological data from the PULLMAN 2 NW weather station including daily precipitation, snowfall, minimum and maximum temperature; soil core bulk density, total carbon, nitrogen, and pH by depth from samples taken in 1998/1999, 2008, and 2015; management data included planting and harvest dates and fertilizer applications by crop and harvest year; high-resolution topographical data (10 m2) including elevation, transformed aspect, slope, and topographic wetness index; dry grain yield per area; and Landsat 5 or 8 imagery for cloud-free passes during the month preceding harvest. Landsat surface reflectance imagery had a spatial resolution of 30 m per pixel, and the yield sampling points each had associated Landsat pixel values. Notably, there was no useful Landsat data for 2012; the only available data was Landsat 7, for which data was corrupted by a scan line correction error. Common Landsat measures of the Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), Plant Senescence Response Index (PSRI), Green Red Vegetation Index (GRVI), Green Normalized Difference Vegetation Index (GNDVI), and Normalized Difference Tillage Index (NDTI) were chosen as initial predictors before variable selection. Indices were taken as the value 30 days before harvest, which generally corresponds to anthesis and full vegetation cover and is a representative value. Others have found remote sensing indices around this time period to be most correlated with yield (Tuvendorj et al., 2019; Huang et al., 2014). Ongoing work with the same dataset is exploring the within-season changes in multispectral indices. However, it produces far more predictor variables than is practical for the methodology explored in this paper.
The first step was to convert soil core composition by depth to totals for the entire soil column. Soil profile total carbon and total nitrogen were obtained by multiplying each layer C and N concentration by layer thickness and bulk density, then summing over the entire profile to 153 cm to provide a single value of carbon grams per cm2 and nitrogen grams per cm2. Because pH is log-scale, it was first converted into mols H+ per liter, then grams of H+ per cm3, then multiplied by layer thickness and layer porosity (derived from bulk density). Those layer-wise values were summed over the column to get a single value of grams H+ per cm2 as an indicator of acidity. Next, because it varied slowly over time, soil core data were interpolated using inverse distance weighting (IDW) using the R package gstat (Pebesma, 2004) over time and space to match the yield sampling locations. A more rigorous approach, in the Bayesian context, is to add an additional model to the hierarchy to impute the missing covariates (Gómez-Rubio et al., 2019), which is the subject of a separate analysis.
Since the Cook Agronomy Farm dataset only contains harvest values, not in-season biomass, daily weather values were summarized to get yearly values. Maximum and minimum daily temperatures were averaged, first by day, then over the calendar year prior to harvest since this would broadly correspond to growing degree days. Precipitation and snowfall (scaled by 13 to equivalent water) were likewise summed over the calendar year before harvest, as in this dryland cropping system, most of the precipitation falls in a short span of time. Management values (Nf, Pf, and Kf for nitrogen, phosphorous, and potassium grams per square meter) were summed over the harvest year by crop type. As with remote sensing variables, including within-season timing of management variables produces far more variables than is practical for the methodology described here and is the subject of a separate analysis.
Several of the measured or derived topographical variables were used: elevation, slope, transformed aspect, topographical wetness index, and analytical hillshade. Of the topographical variables, the transformed aspect was converted into a categorical variable (Face): values over 0.5 were labeled North, and otherwise were labeled South.
Regardless of crop type, yield data typically follows a normal distribution and was scaled to z-scores; all other variables were similarly scaled, but over the entire dataset. This helped in model computation, as yields of different crops were very different, and variables had many different ranges. At the end of this process, the data is in a single array, containing data for every georeferenced point during each harvest year. There are large gaps in the soil data compared with the yield data, since the soil data was collected approximately every 8 years and at a subset of points. This data was gap-filled using Inverse Distance Weighting (IDW). There were also some gaps in the yield data due to missing samples or rotations where yield was not taken (i.e., alfalfa, fallow), and these were not filled. Consequently, we created two datasets: one including soil variables (C, N, and H) and another using remote sensing indices to determine variable importance (table 2). Otherwise, both datasets contain the same topographic, management, and climate information.
Analysis
Two modeling strategies were evaluated: LM with R package lme4 (Bates et al., 2009) and BHM with R package INLA (Martins et al., 2012). For LM, we considered four random effects structures: a normally distributed error, a sample location random effect, a year random effect, or both. For BHM, we considered a normally distributed error, an order one spatial conditional autoregressive model (Besag), an order one temporal random walk (RW(1)), or both Besag and RW(1) together. Models with space/time interaction could not be stably computed because there was only at most one data point per space-time point.
 |
| Figure 2. Flow chart of data processing and modeling. Diagram describes how the models were parameterized for soil and remote sensing variables. IDW = Inverse Distance Weighting; BHM = Bayesian Hierarchical Model; LM = Linear Model. |
| Table 2. Model linear equations after variable selection. C, H, N = soil column carbon, hydrogen ion, and nitrogen; Kf, Nf, Pf = elemental applied potassium, nitrogen, and phosphorous per area. | |
| Model | Variables |
| Soil variables | C, H, N, Crop, Elevation, Face, Kf, Nf, Pf, Precip, Slope, Tave, Crop:Elevation, Crop:Face, Crop:H, Crop:N, Crop:Slope |
| Remote sensing variables | Crop, Elevation, Face, Kf, Pf, Precip, GNDVI, GRVI, NDTI, NDVI, PSRI, SAVI, Crop:Elevation, Crop:GNDVI, Crop:GRVI, Crop:NDTI, Crop:NDVI, Crop:PSRI |
The first step was to select variables, which was done via a two-stage process for the LM including soil information and the LM including remotely sensed data. First, a full model was constructed, including crop interactions with topography and other variables. Crop type interaction with weather or management led to rank deficiency for certain crop types, so these terms were not included. Then, this model was reduced by an exhaustive search for the combination optimizing Aikake Information Criterion (AIC), using dredge from the MuMIn package in R (Barton, 2009). This combination of fixed effects (detailed in the Results section) was used in the LM and BHM.
To test model performance, we looked at several metrics. First, we examined root mean square error (RMSE), mean absolute difference (MAD) of predicted normalized yield, and Deviance Information Criterion (DIC). There is agreement that random cross-validation is optimistic and block cross-validation is pessimistic, so we examined both (Oliveira et al., 2021). Two cross-validations were tested: a five-fold random cross-validation, and a block cross-validation, using 6 contiguous space-time blocks (dividing the dataset into two contiguous space-blocks and three contiguous time-blocks). In addition, we looked at plots of predicted versus observed for the cross-validation and spatial plots of predictions to assess any spatial patterns or general tendencies. To show the uncertainty estimation of the BHM, we plotted the standard deviation of the probability distribution of predicted relative yield over space and time.
 |
| Figure 3. Pairs plot of soil chemistry, weather, and yield, color-coded by crop type. Spring wheat = SW; Winter wheat = WW; Spring canola = SC; Winter canola = WC; Spring barley = SB; Spring pea = SP; Winter barley = WB; Winter pea = WP; Winter triticale = WT; Winter lentil = WL; Garbanzo Beans = GB. |
Finally, to demonstrate the missing data capabilities of BHM, we performed an additional test on the soil variable dataset, where most data was missing. BHM has unique gap-filling properties not present in IDW or other nonprobabilistic methods. First, linear models were developed for C, N, and H using the non-gap-filled dataset and fit using a BHM, with explanatory values of crop type, topographic wetness index, crop/topographic wetness index interaction, applied nitrogen, and annual precipitation. Then, the marginal distributions of these BHMs were sampled n = 50 times to get data to fill the gaps. These sets of gap-filled data were then fit to n BHMs, separately. Then the n models were averaged in a Bayesian sense. Finally, we examined the estimated fixed effects of C, N, and H for both the BHMs fit with either IDW-filled or BHM-filled data to see the effect of the gap-filling method.
Results and Discussion
Important Variables
Within each data class, some variables conformed to a Gaussian distribution, such as normalized yield (GYZ), elevation, and slope; however, some did not, such as transformed aspect (TRASP). This is an important feature to note due to the impact that the non-normal data has on LM. Initial examination of figure 3 shows GYZ has a relationship with soil chemistry variables C, N, and H (P < 0.001). The strong correlation between C and N is to be expected due to organic matter in the soil contributing to the total amount of nitrogen each year. However, the correlation between C and H (pH) was crop-dependent; for example, winter canola (WC) had a strong correlation (r = 0.72, P < 0.001), whereas spring canola (SC) did not (r = 0.1, P < 0.05). This crop dependent variability adds to the complexity of yield predictions, especially when considering the confounding effect of variables such as annual precipitation, which vary independently from crop type. There is also collinearity present, particularly between C and N. For remote sensing variables (fig. 4), strong relationships exist for all indices considered, which is expected because all are indications of biomass or senescence, and would therefore correlate with yield (Tuvendorj et al., 2019). Collinearity also exists between remote sense variables because these indices have similar formulations (e.g., SAVI and NDVI). Finally, relationships between GYZ and topography are significant, with low correlation coefficients indicating that generally, the relationship between GYZ was crop dependent (fig. 5). For example, elevation negatively influences spring peas (SP) (r = -0.297, P < 0.001) but positively impacts winter pea yields (WP) (r = 0.307, P < 0.001) (fig. 5). Not shown are average temperature, total precipitation, and fertilizer regime, for which crop-type interaction was not considered due to rank deficiency. However, these variables all have a significant predictive power for normalized yield. Variable selection using dredge (an exhaustive search of combinations of variables, choosing the one that minimized the Aikake Information Criterion) yielded the equations presented in table 2. Most variables considered in the soil variables model were good predictors, with crop-type interaction terms. Similarly, the remote sensing indices and their crop-type interactions were good predictors. However, for the remote sensing model, applied N and average temperature do not appear after selection to minimize AIC. This may be because the remote sensing indices were able to explain more of the variation in yield than soil variables.
 |
| Figure 4. Pairs plot of remote sensing indices and yield. |
Model Performance
Model fit for the two types of models (LM vs. BHM) and explanatory variables (soil vs. remote sensing), cross-validation method (random vs. block), and random effects structure was tested using the measures RMSE, MAD, and DIC. The random-effects structures assigned for LM vs. BHM were slightly different, but both still contained the same spatial (field location ID) and temporal (year) information. The goodness of fit measures for each model, validation method, and random effects structure are given in table 3. Of the three measures of fit, RMSE and MAD were highly correlated (r = 0.99, P <0.001). MAD serves to indicate model error magnitude, while DIC gives complexity-penalized fit. A direct comparison of DIC for LM vs. DIC for BHM may not be valid since the latter has actual estimates of the likelihood and posterior marginals for the parameters, so the calculations are different. Regardless, DIC is a valid metric because the distribution of GYZ is approximately normal.
 |
| Figure 5. Pairs plot of normalized yield and topographic variables. |
Overall, the best performing model was the remote sensing BHM with Besag spatial effects and random walk temporal effects, and the worst performing was the LM using soil variables and no random effects. Evaluation of model performance through random and block cross-validation (table 3) reveals several points. First, as expected, random cross-validation gives more favorable DIC metrics for all models compared with block cross-validation for both LM and BHM. Random cross-validation only consistently gives better MAD scores for LM. Taken as a whole, this shows the importance of considering spatio-temporal effects in measuring model performance. Second, incorporating spatio-temporal effects in either the LM or BHM improves model performance, with the best model fit occurring with the addition of both space and time. Third, considering MAD scores, the BHM outperforms the equivalent LM. For example, the remote sensing LM with space and time random effects performs worse than the corresponding BHM. This is likely due to the correlation structures in the random effects of the latter, absent in the former, and the BHM’s ability to handle non-Gaussian distributions. Finally, the BHM models have more similar performance for random and block cross-validation than the LM. This means that whether the training data is contiguous in space-time or not, the BHM is performing about the same, whereas the LM performs differently depending on the continuity of the data. This implies the BHM is better able to predict new spatial or temporal blocks than the equivalent LM.
We examined the model performance further for selected models in figure 6 by plotting predicted vs. observed normalized yield for a few different models. Ideally, there would be a 1:1 relationship between predicted and observed values, with little scatter in the data points. Without spatio-temporal effects, the slope is quite low, with less variation present in the predicted values. This makes sense as the random effects structures add to the total variation in the predicted values. Models that considered spatio-temporal effects had improved slope estimates and R2 values. This was true for both LM and BHM. Looking at spatio-temporal trends by comparing figure 7 (observations by year over space) to predictions in figure 8 (BHM predictions with remote sensing variables and no spatio-temporal effects) and figure 9 (BHM predictions with remote sensing variables and spatio-temporal effects) shows that qualitatively, the spatio-temporal random effects provide smoother estimates over space and time. A striking example can be seen by comparing observed vs. predicted values for the years 2004 and 2008 in figures 7-9. The central geo-referenced points of the field in these years have contrasting normalized observed yield values (fig. 7), which are missed by the BHM with no spatio-temporal effects (fig. 8) but are easily observable in the BHM with spatio-temporal effects (fig. 9).
| Table 3. Model performance metrics. | |||||
| Model | Cross Validation | Random Effects | RMSE | MAD | DIC |
| Soil LM | Random | None | 0.8915 | 0.7009 | 13257.9644 |
| Soil LM | Block | None | 0.9438 | 0.7467 | 13752.4554 |
| Soil BHM | Random | None | 0.8748 | 0.6879 | 13301.8879 |
| Soil BHM | Block | None | 0.8705 | 0.6844 | 13796.2852 |
| Soil LM | Random | Point ID | 0.8476 | 0.6606 | 12549.2811 |
| Soil LM | Block | Point ID | 0.9116 | 0.7153 | 13010.5923 |
| Soil BHM | Random | Point ID (Besag) | 0.7961 | 0.6212 | 12813.9137 |
| Soil BHM | Block | Point ID (Besag) | 0.7926 | 0.6181 | 13275.2466 |
| Soil LM | Random | Year | 0.7875 | 0.612 | 11671.722 |
| Soil LM | Block | Year | 0.8329 | 0.65 | 12118.1259 |
| Soil BHM | Random | Year (RW(1)) | 0.7716 | 0.6039 | 12047.5914 |
| Soil BHM | Block | Year (RW(1)) | 0.7685 | 0.5973 | 12498.6251 |
| Soil LM | Random | Point ID + Year | 0.7274 | 0.5533 | 11105.9181 |
| Soil LM | Block | PointID + Year | 0.776 | 0.5908 | 11512.9177 |
| Soil BHM | Random | PointID + Year (Besag + RW(1)) | 0.674 | 0.5178 | 11216.9102 |
| Soil BHM | Block | PointID + Year (Besag + RW(1)) | 0.672 | 0.5122 | 11622.0256 |
| RS LM | Random | None | 0.8488 | 0.6694 | 12067.9404 |
| RS LM | Block | None | 0.9152 | 0.7254 | 12506.0768 |
| RS BHM | Random | None | 0.8425 | 0.6658 | 12928.949 |
| RS BHM | Block | None | 0.8376 | 0.6624 | 13398.3579 |
| RS LM | Random | Point ID | 0.7935 | 0.6146 | 11295.4495 |
| RS LM | Block | Point ID | 0.8671 | 0.6758 | 11694.2401 |
| RS BHM | Random | PointID (Besag) | 0.7426 | 0.5712 | 12177.8698 |
| RS BHM | Block | PointID (Besag) | 0.7388 | 0.5739 | 12601.9564 |
| RS LM | Random | Year | 0.7879 | 0.6185 | 11103.8602 |
| RS LM | Block | Year | 0.8407 | 0.6639 | 11523.726 |
| RS BHM | Random | Year (RW(1)) | 0.7691 | 0.6045 | 12028.7832 |
| RS BHM | Block | Year (RW(1)) | 0.7659 | 0.6003 | 12476.8801 |
| RS LM | Random | Point ID + Year | 0.7131 | 0.545 | 10394.8457 |
| RS LM | Block | Point ID + Year | 0.7618 | 0.5845 | 10775.7241 |
| RS BHM | Random | Point ID + Year (Besag + RW(1)) | 0.6486 | 0.4991 | 10880.1533 |
| RS BHM | Block | Point ID + Year (Besag + RW(1)) | 0.6469 | 0.4954 | 11272.3022 |
Figure 10 shows the standard deviation over space for the model shown in figure 9, demonstrating that areas of the field are more variable over time. Strips of high uncertainty are often associated with crops in the rotation that have fewer data for model fitting, but some areas, like the central-east area, have higher uncertainty than their neighbors regardless.
Missing Data Handling
Figure 11 shows the results of the test of missing data capabilities. The posterior marginals of C, N, and H are shown for data filling techniques. There is considerable spread in the estimates using Bayesian techniques because it properly incorporates the uncertainty in the missing data. Fitting a BHM using the IDW-filled data gives artificially low uncertainty in the estimates, which will translate to artificially low estimates of uncertainty in yield estimates. This shows the advantage of the BHM approach in the presence of missing covariates (Gómez-Rubio et al., 2019; Shelton et al., 2012).
Model Utility and Applications
The BHM with spatio-temporal effects was shown to not only have higher predictive power than an LM but also produce more realistic yield predictions compared to a BHM without spatio-temporal effects. The selected BHM model can utilize a wide variety of explanatory variables that are commonly available to researchers and farmers alike (satellite, topographic, and soil pH and nitrogen) to make relatively accurate yield predictions for a plethora of crop species.
Agricultural researchers will find this modeling technique useful, especially for on-farm research, where plots may be several acres in size and results can be aggregated over a large area. A BHM has been successfully used for yield predictions across a farmer's network in Iowa; however, within the field, spatial variability and climatic factors were not considered (Laurent et al., 2019). The ability to a priori map yield variability in a field will allow researchers to observe potential variation across a large spatial scale. This may improve research plot placement under limited replication, for instance, for different crop treatments (e.g., SP vs. WP), by reducing the possibility of confounding treatment effects with the inherent variation. Furthermore, the BHM has other benefits over the LM in that it should be possible to incorporate farmer knowledge into field maps to add confidence to predictions. There have been several reports demonstrating that farmers’ experiential knowledge, gained over decades or generations, is both highly predictive of in-field yield variation and is correlated with remote sense data such as NDVI (Oliver et al., 2010; Heijting et al., 2011). Future research could use this experiential data as an explanatory input to improve BHM. A grower would find this information directly useful in management decisions; where variance is high, more inputs (e.g., nitrogen, irrigation) can still risk low yields and thus may not be worth their investment. Or they might be more risk-tolerant and assess that the potential gains due to increased fertilizer, for example, are worth the potential for low yields regardless. Though strictly speaking, the models developed here were descriptive, not predictive, they are useful in several ways. First, getting in-season pre-market estimates of yield is useful for grower financial planning. Second, a dominant correlate of the next year’s yield is the previous year’s yield. The model described here can be combined with another level in the hierarchy that predicts the next year’s yield using the previous year’s yield. Finally, for policymakers and researchers, regional-level maps of yield with uncertainty can be used to guide recommendations to growers for rotational changes, management strategies, and marketing.
 |
| Figure 6. Predicted vs. observed normalized yield for selected models. Numbers in parentheses indicate (R2, slope) of the points. BHM = Bayesian hierarchical model; LM = linear model; RE = random effect; RS = remote sensing; RW(1) = one degree random walk. |
Conclusions
Two classes of yield prediction models were compared (LM and BHM), with different error structures and sources of predictive data. For a long-term, multi-crop study, the BHM using spatio-temporal random errors and RS information performs best by RMSE, MAD, and DIC. This study highlights the importance of considering spatio-temporal effects in long-term crop models as well as the superiority of the BHM approach for this application. Such a model provides a predictive framework that can function even if covariates are missing and, through the specification of priors, can incorporate farmer knowledge. The BHM also provides distributions for predictions, which is useful for policymakers and growers alike who need to assess risk in concrete terms. Knowing the uncertainty of the yield allows one to assess if additional expense on fertilizer or prescription agriculture to increase yield is worth the risk of potential low yields due to random effects. Further study is needed to evaluate the power of multivariate models (e.g., grain yield and grain protein together) and grouped spatial and temporal effects, where the spatial error itself is correlated in time. Within-season variation of remote sensing indices and weather could provide more predictive accuracy, and this is another area of active research. Additionally, the models described here can be extended to make predictions about future years’ yields, making them a more valuable tool for in- or pre-season management decisions.
 |
| Figure 7. Observations of normalized yield over all study years. |
 |
| Figure 8. Predicted normalized yield from remote sensing BHM with no spatio-temporal effects. |
 |
| Figure 9. Predictions of normalized yield from best performing BHM (includes spatio-temporal effects). |
 |
| Figure 10. Uncertainty, as measured by standard deviation of normalized yield, for BHM predictions with RS information and spatial/temporal error structures. |
 |
| Figure 11. Posterior marginals for soil properties’ fixed effects using IDW filling (dashed) and BHM filling (solid). |
Acknowledgments
This research was a contribution from the Long-Term Agroecosystem Research (LTAR) network. LTAR is supported by the United States Department of Agriculture. This research used resources provided by the SCINet project of the USDA Agricultural Research Service, ARS project number 0500-00093-001-00-D.
References
Barton, K. (2009). MuMIn: Multi-model inference. R package version 1.0.0. Retrieved from http://r-forge.r-project.org/projects/mumin/
Bates, D., Maechler, M., Bolker, B., Walker, S., Christensen, R. H., Singmann, H.,... Grothendieck, G. (2009). Package ‘lme4’. Retrieved from http://lme4.r-forge.r-project.org
Blasch, G., Li, Z., & Taylor, J. A. (2020). Multi-temporal yield pattern analysis method for deriving yield zones in crop production systems. Precis. Agric., 21, 1263-1290. https://doi.org/10.1007/s11119-020-09719-1
Broemeling, L. D. (2017). Bayesian analysis of linear models. CRC Press. https://doi.org/10.1201/9781315138145
Castelvecchi, D. (2016). Can we open the black box of AI? Nature News, 538(7623). https://doi.org/10.1038/538020a
Chandra, H. (2013). Exploring spatial dependence in area-level random effect model for disagreggate-level crop yield estimation. J. Appl. Stat., 40(4), 823-842. https://doi.org/10.1080/02664763.2012.756858
Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L.,... Böhner, J. (2015). System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev., 8(7), 1991-2007. https://doi.org/10.5194/gmd-8-1991-2015
Crimp, S., Bakar, K. S., Kokic, P., Jin, H., Nicholls, N., & Howden, M. (2015). Bayesian space-time model to analyse frost risk for agriculture in Southeast Australia. Int. J. Climatol., 35(8), 2092-2108. https://doi.org/10.1002/joc.4109
Duda, R. O., & Hart, P. E. (2006). Pattern classification. John Wiley & Sons.
Gandhi, N., Armstrong, L. J., Petkar, O., & Tripathy, A. K. (2016). Rice crop yield prediction in India using support vector machines. Proc. 2016 13th Int. Joint Conf. on Computer Science and Software Engineering (JCSSE) (pp. 1-5). IEEE. https://doi.org/10.1109/JCSSE.2016.7748856
Gómez-Rubio, V., Cameletti, M., & Blangiardo, M. (2019). Missing data analysis and imputation via latent Gaussian Markov random fields. arXiv preprint arXiv:1912.10981. https://doi.org/10.48550/arXiv.1912.10981
Heijting, S., de Bruin, S., & Bregt, A. K. (2011). The arable farmer as the assessor of within-field soil variation. Precis. Agric., 12(4), 488-507. https://doi.org/10.1007/s11119-010-9197-y
Hooten, M. B., Wikle, C. K., Dorazio, R. M., & Royle, J. A. (2007). Hierarchical spatiotemporal matrix models for characterizing invasions. Biometrics, 63(2), 558-567. https://doi.org/10.1111/j.1541-0420.2006.00725.x
Huang, J., Wang, H., Dai, Q., & Han, D. (2014). Analysis of NDVI data for crop identification and yield estimation. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 7(11), 4374-4384. https://doi.org/10.1109/JSTARS.2014.2334332
Huggins, D. R. (2015). The Cook Agronomy Farm LTAR: Knowledge Intensive Precision Agro-ecology. Proc. American Geophysical Union, Fall Meeting 2015, (pp. H53B-1665).
Huggins, D. R., & Uberuaga, D. P. (2010). Field heterogeneity of soil organic carbon and relationships to soil properties and terrain attributes. In Climate friendly farming: Center for Sustaining Agriculture and Natural Resources Research Report 2010.
Jennings, M. D., Miller, B. C., Bezdicek, D. F., & Granatstein, D. (1990). Sustainability of dryland cropping in the Palouse: An historical view. J. Soil Water Conserv., 45(1), 75-80.
Jones, J. W., Hoogenboom, G., Porter, C. H., Boote, K. J., Batchelor, W. D., Hunt, L. A.,... Ritchie, J. T. (2003). The DSSAT cropping system model. Eur. J. Agron., 18(3-4), 235-265. https://doi.org/10.1016/S1161-0301(02)00107-7
Karimi, T., Stöckle, C. O., Higgins, S. S., Nelson, R. L., & Huggins, D. (2017). Projected dryland cropping system shifts in the Pacific Northwest in response to climate change. Front. Ecol. Evol., 5. https://doi.org/10.3389/fevo.2017.00020
Kass, R. E., & Steffey, D. (1989). Approximate Bayesian inference in conditionally independent hierarchical models (parametric empirical Bayes models). J. Am. Stat. Assoc., 84(407), 717-726. https://doi.org/10.1080/01621459.1989.10478825
Kaur, H., Huggins, D. R., Carlson, B., Stockle, C., & Nelson, R. (2022). Dryland fallow vs flex-cropping decisions in inland Pacific Northwest of USA. Agric. Syst., 201, 103432. https://doi.org/10.1016/j.agsy.2022.103432
Kaur, H., Huggins, D. R., Rupp, R. A., Abatzoglou, J. T., Stockle, C. O., & Reganold, J. P. (2017). Agro-ecological class stability decreases in response to climate change projections for the Pacific Northwest, USA. Front. Ecol. Evol., 5.
Kawakita, S., Takahashi, H., & Moriya, K. (2020). Prediction and parameter uncertainty for winter wheat phenology models depend on model and parameterization method differences. Agric. For. Meteorol., 290, 107998. https://doi.org/10.1016/j.agrformet.2020.107998
Keating, B. A., Carberry, P. S., Hammer, G. L., Probert, M. E., Robertson, M. J., Holzworth, D.,... Smith, C. J. (2003). An overview of APSIM, a model designed for farming systems simulation. Eur. J. Agron., 18(3), 267-288. https://doi.org/10.1016/S1161-0301(02)00108-9
Khaki, S., & Wang, L. (2019). Crop yield prediction using deep neural networks. Front. Plant Sci., 10. https://doi.org/10.3389/fpls.2019.00621
Lacoste, M., Cook, S., McNee, M., Gale, D., Ingram, J., Bellon-Maurel, V.,... Hall, A. (2022). On-Farm experimentation to transform global agriculture. Nature Food, 3(1), 11-18. https://doi.org/10.1038/s43016-021-00424-4
Lamsal, A., Welch, S. M., White, J. W., Thorp, K. R., & Bello, N. M. (2018). Estimating parametric phenotypes that determine anthesis date in Zea mays: Challenges in combining ecophysiological models with genetics. PLoS One, 13(4), e0195841. https://doi.org/10.1371/journal.pone.0195841
Laurent, A., Kyveryga, P., Makowski, D., & Miguez, F. (2019). A framework for visualization and analysis of agronomic field trials from on-farm research networks. Agron. J., 111(6), 2712-2723. https://doi.org/10.2134/agronj2019.02.0135
Liaw, A., & Wiener, M. (2002). Classification and regression by randomForest. R News, 2(3), 18-22.
Lindgren, F., Lindström, J., & Rue, H. (2010). An explicit link between Gaussian fields and Gaussian Markov random fields; The SPDE approach. Preprints Math. Sci., 3. Retrieved from https://lup.lub.lu.se/search/files/5312973/1581115.pdf
Lobell, D. B., & Burke, M. B. (2010). On the use of statistical models to predict crop yield responses to climate change. Agric. For. Meteorol., 150(11), 1443-1452. https://doi.org/10.1016/j.agrformet.2010.07.008
Martins, T. G., Simpson, D., Lindgren, F., & Rue, H. (2012). Bayesian computing with INLA: New features. Comput. Stat. Data Anal., 67. https://doi.org/10.1016/j.csda.2013.04.014
McCulloch, C. E., & Searle, S. R. (2005). Generalized linear mixed models. https://doi.org/10.1002/0470011815.b2a10021
McFarland, C. R. (2016). Liming no-till soils and determining lime requirement in the Palouse region. MS thesis. Washington State University.
Morari, F., Zanella, V., Gobbo, S., Bindi, M., Sartori, L., Pasqui, M.,... Ferrise, R. (2021). Coupling proximal sensing, seasonal forecasts and crop modelling to optimize nitrogen variable rate application in durum wheat. Precis. Agric., 22(1), 75-98. https://doi.org/10.1007/s11119-020-09730-6
Newman, S. J., & Furbank, R. T. (2021). Explainable machine learning models of major crop traits from satellite-monitored continent-wide field trial data. Nat. Plants, 7, 1354-1363. https://doi.org/10.1038/s41477-021-01001-0
NRCS. (2013). Whitman county, WA soil survey. Retrieved from http://websoilsurvey.sc.egov.usda.gov/App/HomePage.htm
Oliveira, M., Torgo, L., & Santos Costa, V. (2021). Evaluation procedures for forecasting with spatiotemporal data. Mathematics, 9(6), 691. https://doi.org/10.3390/math9060691
Oliver, Y. M., Robertson, M. J., & Wong, M. T. (2010). Integrating farmer knowledge, precision agriculture tools, and crop simulation modelling to evaluate management options for poor-performing patches in cropping fields. Eur. J. Agron., 32(1), 40-50. https://doi.org/10.1016/j.eja.2009.05.002
Ozaki, V. A., Ghosh, S. K., Goodwin, B. K., & Shirota, R. (2008). Spatio-temporal modeling of agricultural yield data with an application to pricing crop insurance contracts. Am. J. Agric. Econ., 90(4), 951-961. https://doi.org/10.1111/j.1467-8276.2008.01153.x
Pebesma, E. J. (2004). Multivariable geostatistics in S: The gstat package. Comput. Geosci., 30(7), 683-691. https://doi.org/10.1016/j.cageo.2004.03.012
Peng, J. L., Kim, M. J., Kim, B. W., & Sung, K. I. (2016). A yield estimation model of forage rye based on climate data by Locations in South Korea using General Linear Model. J. Korean Soc. Grassl. Forage Sci., 36(3), 205-214. https://doi.org/10.5333/KGFS.2016.36.3.205
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., Heisterkamp, S., Van Willigen, B., & Maintainer, R. (2017). Package ‘nlme’. Linear and nonlinear mixed effects models. R package version 3.1.
R. Core Team. (2013). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Retrieved from http://www.R-project.org/
Shelton, A. O., Dick, E. J., Pearson, D. E., Ralston, S., & Mangel, M. (2012). Estimating species composition and quantifying uncertainty in multispecies fisheries: Hierarchical Bayesian models for stratified sampling protocols with missing data. Can. J. Fish. Aquat. Sci., 69(2), 231-246. https://doi.org/10.1139/f2011-152
Steduto, P., Hsiao, T. C., Raes, D., & Fereres, E. (2009). AquaCrop — The FAO crop model to simulate yield response to water: I. Concepts and underlying principles. Agron. J., 101(3), 426-437. https://doi.org/10.2134/agronj2008.0139s
Stöckle, C. O., Donatelli, M., & Nelson, R. (2003). CropSyst, a cropping systems simulation model. Eur. J. Agron., 18(3-4), 289-307. https://doi.org/10.1016/S1161-0301(02)00109-0
Stöckle, C. O., Higgins, S., Nelson, R., Abatzoglou, J., Huggins, D., Pan, W.,... Brooks, E. (2018). Evaluating opportunities for an increased role of winter crops as adaptation to climate change in dryland cropping systems of the U.S. Inland Pacific Northwest. Clim. Change, 146, 247-261. https://doi.org/10.1007/s10584-017-1950-z
Thelemann, R., Johnson, G., Sheaffer, C., Banerjee, S., Cai, H., & Wyse, D. (2010). The effect of landscape position on biomass crop yield. Agron. J., 102(2), 513-522. https://doi.org/10.2134/agronj2009.0058
Tuvdendorj, B., Wu, B., Zeng, H., Batdelger, G., & Nanzad, L. (2019). Determination of appropriate remote sensing indices for spring wheat yield estimation in Mongolia. Remote Sens., 11(21), 2568. https://doi.org/10.3390/rs11212568
van Diepen, C. A., Wolf, J., van Keulen, H., & Rappoldt, C. (1989). WOFOST: A simulation model of crop production. Soil Use Manag., 5(1), 16-24. https://doi.org/10.1111/j.1475-2743.1989.tb00755.x
van Klompenburg, T., Kassahuna, A., & Catal, C. (2020). Crop yield prediction using machine learning: A systematic literature review. Comput. Electron. Agric., 177. https://doi.org/10.1016/j.compag.2020.105709
Wang, L., Zhou, X., Zhu, X., Dong, Z., & Guo, W. (2016). Estimation of biomass in wheat using random forest regression algorithm and remote sensing data. Crop J., 4(3), 212-219. https://doi.org/10.1016/j.cj.2016.01.008
Wikle, C. K., & Hooten, M. B. (2010). A general science-based framework for dynamical spatio-temporal models. TEST, 19, 417-451. https://doi.org/10.1007/s11749-010-0209-z
Wikle, C. K., Berliner, L. M., & Cressie, N. (1998). Hierarchical Bayesian space-time models. Environ. Ecol. Stat., 5(2), 117-154. https://doi.org/10.1023/A:1009662704779
Wikle, C. K., Zammit-Mangion, A., & Cressie, N. (2019). Spatio-temporal statistics with R (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781351769723
Yang, C., Peterson, C. L., Shropshire, G. J., & Otawa, T. (1998). Spatial variability of field topography and wheat yield in the Palouse region of the Pacific Northwest. Trans. ASABE, 41(1), 17-27. https://doi.org/10.13031/2013.17147
Zammit-Mangion, A. (2022). Package IDE: Integro-difference equation spatio-temporal models.