Comparing a Random Forest Based Prediction of Winter Wheat Yield to Historical Yield Potential

Predicting wheat yield is crucial due to the importance of wheat across the world. When modeling yield, the difference between potential and actual yield consistently changes because of advances in technology. Considering historical yield potential would help determine spatiotemporal trends in agricultural development. Comparing current and historical yields in Denmark is possible because yield potential has been documented throughout history. However, the current national winter wheat yield map solely uses soil properties within the model. The aim of this study was to generate a new Danish winter wheat yield map and compare the results to historical yield potential. Utilizing random forest with soil, climate, and topography variables, a winter wheat yield map was generated from 876 field trials carried out from 1992 to 2018. The random forest model performed better than the model based only on soil. The updated national yield map was then compared to yield potential maps from 1688 and 1844. While historical time periods are characterized by numerous low yield potential areas and few highly productive areas, current yield is evenly distributed between low and high yields. Advances in technology and farm practices have exceeded historical yield predictions, mainly due to the use of fertilizer, irrigation, and drainage. Thus, modeling yield projections could be unreliable in the future as technology progresses.


Introduction
Developing accurate national yield models is challenging due to variations in growth conditions across large areas and changes through time. A range of models have been created to simulate crop growth because of the importance of predicting yield accurately (e.g., CERES-Wheat: [1]; WOFOST: [2]; Sirius: [3]). Correctly modeling yield is necessary for an array of crops, but a large number of studies have focused on wheat. Wheat is a staple food globally, corresponding to approximately 30% of total world cereal production, and provides a large portion of both daily calories and protein for people across the world [4]. Wheat is grown in nearly all parts of Europe [5], spanning a wide range of soil types and climate zones, meaning the models are required to be diverse and robust.
The importance of wheat has driven research to determine what affects yield potential [6][7][8][9][10]. Many decisions made by farmers alter yields (crop variety, planting date, irrigation, etc.); however, many factors affecting yield potential are independent from the farmer's decisions. Three major components affecting yield are: soil, climate, and topography. Soil properties, such as structure and clay percentage, either hinder or promote growth [11]. Fluctuations in temperature and precipitation have been found to be key reasons for uncertainty in yield predictions [12]. Topography can impact yield by affecting water flow, energy supply, and accumulation of soil and nutrients in an area [13,14]. Y opt = 4.78 + 6.67 × (% Clay + % Silt) + 0.11 × R H − 0.21 × (% Clay + % Silt) × R H (1) where Y opt is predicted yield in t ha −1 and R H is (% Clay + % Silt)/% Organic Carbon. The model includes clay, silt, and organic carbon of the topsoil (0-30 cm) only. This study has two main objectives: (1) to update the current national winter wheat yield map of Denmark, and (2) to compare current yield to historical yield potential. The first objective uses an array of soil, climate, and topography variables with random forest to generate a more accurate yield map. The updated yield map is compared to the existing yield map to ensure the addition of climate and topography have increased yield prediction accuracy. The second objective will use historical yield potential maps from 1688 and 1844 for comparisons with the newly developed yield map. Comparing yield through time will help determine how advances in technology have driven agriculture in Denmark.

Study Area
Denmark is a Northern European country lying between 54.3-57.5°N and 8.0-15.2°E with a temperate climate. The mean annual temperature is 9°C with summer days averaging to 17°C and winter days averaging to 1°C. Mean annual precipitation is 770 mm, with the west receiving more rain than the east [27]. Multiple glaciations during the Weichselian geological period have resulted in highly variable soil distribution and formation [28]. Soils are mainly loamy Luvisols in eastern and central regions and sandy Podzols in western regions [29]. Moraine landscapes dominate the east and west, with the west being glaciofluvial flood plains and older eroded moraine landscapes. Raised seabeds dominate the north. The country covers 43,000 km 2 of area, and due to many glacial advances, topography is relatively flat with a mean elevation of 32 m and a maximum elevation of 171 m.

Winter Wheat Fields
Winter wheat (Triticum aestivum L., spp. vulgare) yields were acquired from the Danish National Field Trials managed by the Danish Agricultural Advisory Service [30]. The data from 876 fields from 1992 to 2018 were used in this study ( Figure 1). Bornholm was excluded due to the lack of fields on the island. Each field follows an experimental design that uses 4 or 5 randomized blocks (30 m 2 plot sizes) with every block receiving a different level of mineral nitrogen fertilizer. The maximum yield (in kg ha −1 ) per field was used in this study. The average mineral nitrogen fertilizer applied was 202 kg N ha −1 . Fields covered all climate regions and soil types in Denmark, as different locations were selected each year. An average of 32 fields was selected each year (range between 14 and 60). Yield increased gradually from 1992 to 2018 (Figure 2A). For every year, average yield increases by 70 kg ha −1 (linear regression: p-value < 0.001). Thus, the trend caused by differences between years was removed ( Figure 2B). This increase could be attributed to many factors (advances in farming technology, better management techniques, higher annual temperature, etc.). The trend was removed by subtracting a correction factor from each yield value. The correction factor is calculated by subtracting the global average from the year's average, resulting in different correction factors for each year.

Yield Modeling
Random forest (RF) is used to model winter wheat yields from field trials. This method is a supervised decision tree algorithm that bags un-pruned trees by using randomly selected covariates at each split [31]. The RF method grows a large number of trees from a set of independent observations. Randomness is used to select a variable at each split in the tree, and this variable is from a random subset of all predictor variables. The RF model only generates the mean value for each tree, but using quantile regression forest [32], the distribution of values predicted for each tree can be utilized. The difference in values from quantiles (0.1 and 0.9) is used as the prediction interval for the model. Both methods are implemented from the ranger package in R [33].
The RF model was generated 100 times on different training sets to determine the model performance on different input data. Winter wheat fields are divided into training and test datasets. The training dataset represents 75% of the original data (n = 657) and was used to build the model; the test dataset represents the remaining 25% of the original data (n = 219) and was used to assess model performance. The final winter wheat yield map was generated using all of the data. Variable importance was calculated by averaging variable importance for all 100 model runs. The method used for calculating variable importance within the RF model was the permutation method.
The RF model was compared to the existing winter wheat yield map of Denmark [26]. The previous method used the same source to attain winter wheat yield values as the RF model. Values from the winter wheat yield map based only on soil data are extracted at each field location to determine whether the RF model improved yield predictions. To ensure that the correction factor for a year did not create a bias towards the RF model, uncorrected and corrected values were compared to both the RF and soil-only models. Model performance was assessed using: mean squared error (MSE), root mean squared error (RMSE), Lin's concordance correlation coefficient (CCC), and coefficient of determination (R 2 ).

Predictor Variables
Covariates are needed to build the RF model to find the association between the environment and yield. Covariates are divided into three groups (Table 1): soil (12 covariates), topography (12 covariates), and climate (4 covariates). Covariates were chosen based on expert knowledge and previous studies indicating the predictor variables affecting crop growth. All covariates were generated or resampled to have a 30.4-m spatial resolution. Danish variables created for digital soil mapping typically have a 30.4-m resolution because the digital elevation model (DEM) was derived from 1.6-m resolution national LiDAR data and resampled to 30.4-m resolution [34].
Climate variables are based on temperature, solar radiation, and precipitation. Daily mean temperature in°C at 2 m above the ground and global solar radiation in MJ/m 2 were downloaded from Aarhus University's meteorological database [41]. Downloaded climate data were from 1985 to 2013 at a resolution of 40 km. Kriging was performed on 40-km grids to generate maps of: annual number of growing days above 10°C, average annual temperature, and global solar radiation during the growing season (between 1 April and 31 October). The solar radiation map was combined with a DEM solar radiation map, generated from the Area Solar Radiation tool in ArcMap, to create a topographically corrected solar radiation map. The precipitation map was derived from average total rainfall during the growing seasons between 1961 and 1990 [42].  1 Topography variables are derivatives of the digital elevation model [34] and were generated in SAGA. 2 Raw data from AGRO Climate Database [41]; kriging was used on points to obtain the map.

Historical Comparison
Historical yield potential values are used to assess temporal changes of production through Denmark. The 1688 ( [23]) and 1844 ( [43]) land registers have been mapped using previously digitized data [44]. The proper data corrections were made to ensure both time periods aligned spatially. Hard grain ha −1 values from 1688 and 1844 were assigned to each parish from 1820. Land in the registers was evaluated based on yield potential using a unit of measure called barrels of hard grain, but different definitions were used to quantify hard grain in 1688 than in 1844. Thus, values were normalized from 0 to 1 to determine relative differences. In both definitions, hard grain was based on field measurements and evaluations which could be expressed in terms of production.
The three time periods (1688, 1844, and current) are compared to determine the spatial trend in agricultural production. The winter wheat yield values were averaged for each parish and normalized. The first two time periods assess yield potential, and the current values determine actual yield. While the comparison is based on relative values, the distribution of values for each time period will indicate how national farming has progressed. The comparisons will illuminate areas that have always been viewed as highly productive and areas that have recently been developed due to advances in farming.

Random Forest Predictor Variables
The predictor variable importance indicates that using climate is crucial when predicting winter wheat yield ( Table 2). All four climate variables were within the top 10 predictor variables and significantly correlated with yield. Organic carbon at a depth of 0-30 cm was significant and negatively correlated with yield, but had low importance. Organic carbon was one of the three properties used in the soil-only model, but neither organic carbon variable (depths 0-30 and 100-200 cm) was within the top 15 predictors. Clay at a depth of 0-30 cm was used in the soil-only model and was an important predictor for yield; however, clay at a depth of 100-200 cm was also crucial in the RF model, but was not utilized in the soil-only model. Winter wheat roots have been shown to reach 2 m depth [45], indicating subsoil properties are necessary for predicting winter wheat yield. This shows that the soil-only model might be too simple to accurately predict yield values. Table 2. Top 15 predictors ordered by importance in the random forest model using permutation as the method to determine variable importance. Correlation is indicated by a + for positive correlation, -for negative correlation, and C for categorical variable. A bold + or -indicates whether predictors are significantly correlated with winter wheat yield (linear regression: p-value < 0.05). -Topography variables did not have a large influence on yield predictions. Elevation was the only variable within the top 15 variables. While some of the topography variables were significantly correlated with yield, the overall importance of these variables was low. Contrary to our results, studies have found that topographic position has a crucial effect on wheat yield. While topography explained a lower amount of variation in yield than soil properties, Kravchenko and Bullock [46] indicated elevation was the most influential topographic variable on yield. However, due to the Danish landscape being relatively flat, the addition of topography variables does not play a large role in yield predictions. With soil being an important factor for predicting yield, and topography being a large component in modeling soil, topography is indirectly associated with yield. All of the individual soil properties used in this study have incorporated topography when being modeled and mapped. Thus, adding topography variables alone is likely redundant and seems to have only made the RF model more complicated and computationally expensive.

Rank Predictor Variables Importance Correlation
The addition of climate variables in modeling yield is an important aspect many researchers are focusing on, due to shifts in climate patterns. In varying climate simulations, yield was predicted to decrease with an increase in temperature variability, while changes in yield were mediated by soil type when precipitation was varied [47]. In yield forecasting, models incorporating seasonal temperature and precipitation were able to reliably predict yield variability for 25-38% of the globally harvested area before harvesting [48]. Kristensen et al. [49] stated winter wheat yield in Denmark will likely decrease with the projected climate change and Webber et al. [50] concluded the same results but for all of Europe. Yield has been predicted to decrease by 550 kg ha −1 from 1985 to 2040 with varying degrees of yield reduction depending on the location and soil type. Using the same field trial data source, Ørum et al. [51] found many aspects of climate were important for predicting winter wheat yield in Denmark. Average temperature and precipitation, during certain months on varying soil types, impact plant growth significantly [52]. Besides climate variables, atmospheric carbon dioxide (CO 2 ) should be included in modeling winter wheat yield because the increase in CO 2 and the resulting interaction with environmental variables are expected to increase yield [53]. Broberg et al. [54] elevated CO 2 levels from 372 to 605 ppm and found winter wheat yield increased by 26%. While the seasonal climate during the growing period does not change drastically in Denmark, climate variables and CO 2 need to be included in future yield models.

Model Comparison
The RF model used for predicting winter wheat yield performed better than the previous winter wheat yield model based only on soil properties in all model performance evaluations (Table 3). The correction factor applied to the winter wheat values did not bias the RF model, since the soil-only model performed better when comparing to corrected values to uncorrected values. The RF model's MSE was half the one of the soil-only model, and the RF model's RMSE was 726 kg ha −1 lower than that of the soil-only model. Lin's CCC constitutes a measurement for both precision and accuracy [55], and a valuable way to measure agreement between two methods measuring the same variable. The coefficient was ranked poor for the soil-only model (0.21) and nearly ranked fair to good (0.38) for the RF model using the Fleiss interpretation [56]. The coefficient of determination (R 2 ) is 0.06 higher for the RF model than for the soil-only model. This indicates the RF model explains 6% more variation than the previous model. Looking at model performances, while the R 2 did not increase drastically, the incorporation of climate and topography are beneficial in modeling winter wheat yield, as indicated by every other performance measurement.
Distributions of the predicted yields for both models do not account for extreme values, as seen in the raw data ( Figure 3). The RF model did not perform well when calculating extreme values, but the predictions are centered around the observed data, whereas the soil-only model consistently predicts low yields. The ability to predict extreme values could be improved by incorporating farm management data. Integrating information on irrigation and the previous crops planted have shown to be significant factors that impact yield [51]. However, generating national maps of farm management would be time-consuming and expensive. The placement of field trials could be another reason why the model does not have a higher R 2 value and topography does not have a large contribution. The range of slope for agricultural land in Denmark is between 0 and 34°, while the range of slope for the field trials is between 0 and 12°. The fields change every year but the placement is not completely random, resulting in a biased sampling scheme that does not consider extreme soil and topography locations. Other studies have compared the use of ML algorithms to linear regressions in crop yield studies. Linear regression has out-performed neural networks when using only one input, but neural networks consistently out-performed linear regression when multiple inputs were used [57]. Random forest has been found to produce smaller RMSE values than linear regression when dealing with regional and global scale crop yield predictions [58]. In a study calculating solar radiation from both statistical models and neural networks, the neural networks provided the best estimates on the training data but performed poorly in a test site with different weather conditions. In contrast, the statistical models achieved adequate results in both the training and test sites [59]. The ability to find complex relationships between multiple inputs is the benefit of ML algorithms; however, these algorithms do fall short when applied to datasets with few inputs or when new data arises in test datasets.
The final map from the RF model along with prediction interval (difference in quantiles 0.1 and 0.9) is shown in Figure 4A,B. The pattern indicates that the best yield for winter wheat is in the east with loamy soils and the worst in the west with fluvial sands. This result is consistent with the soil fertility map from 1844 [60] and the existing yield map ( Figure 4C). Areas with the highest and lowest yields also displayed the lowest prediction intervals. This indicates that these areas are typically stable and do not fluctuate as much as other areas. The stability in yield prediction is probably due to soil quality either being good or poor, resulting in consistent yield no matter the climate. The largest prediction intervals are associated with areas that have intermediate yield values. The large prediction interval likely indicates that the yearly fluctuations are greatest in these areas, since the soil quality is average and not a reliable buffer. Thus, climate would have a larger influence on yield.
The two yield maps were normalized between 0 and 1 and subtracted from each other ( Figure 4D). The maps were normalized because the soil-only model predicts negative values in certain areas. These negative values appear for locations where the organic carbon content is below 1% and much smaller than the clay and silt contents. Thus, the soil-only yield map needed to be normalized. The RF model predicted higher values in the east and lower values in the west compared to the previous model. Areas with average yield values are roughly in the same location in both models. The main differences are in the extreme values where the yield is either the highest or lowest. As seen in Figure 4C, the previous model generates a smoother map because only three soil properties were used in the model. The soil-only model originated from work that was focusing on the effects soil organic matter has on crop yield when optimal nitrogen was allocated [26], creating a model that did not incorporate environmental factors. The RF model has sharper boundaries because of variables that take into account glacial boundaries. The lowest values in the RF model from the middle of Jutland are indicative of the boundary between the young moraine and outwash plain landscape types [24]. The fine soil particles (clay and silt) are typically above 20% in the young moraine, resulting in a larger supply of plant available water, and below 10% in the outwash plain, resulting in a smaller supply of plant available water. These differences in soil genesis and soil texture lead to large yield differences in Jutland. The incorporation of variables explaining climate and geology into the model are crucial for depicting the actual pattern seen, rather than the simplified version only including soil.

Historical Comparisons
The three time periods are mapped by parish in Figure 5. Spatial patterns are similar for the first two time periods (1688 and 1844). Comparing the two historical periods, areas in the east had higher yield potentials in 1844 than 1688. Both time periods were similar in the west, and no clear pattern is seen in the north regarding which time period had a higher yield potential. The differences are similar between current yield and both historical yield potentials. The differences indicate that the changes between 1844 and 1688 are minor compared to the drastic changes in current yield. However, areas that have been the least productive are constantly the least productive. Although advances in technology have made the western part of Denmark more productive, the soils in the west are still far less productive than the fertile soils in the east.
Observing the overall trend in parish production, current national yield is more evenly distributed than in previous time periods ( Figure 6). In 1688 and 1844, the data were skewed towards a large portion of parishes having low yield potential and a few having high yield potential. The slight increase in yield potential from 1688 and 1844 could be attributed to the large land reform between the two time periods. During this land reform in the late 18th century, the land that was once owned by the king, church, and manors, and farmed by peasant farmers, was redistributed to the people to allow farmers to cultivate their own land. The increase in agriculture productivity has been associated with the peasant emancipation and the reforms [61]. As the years have progressed, Denmark has become more productive. Current yield is evenly spread across parishes from low to high, with most of the parishes having average yields. The differences in historical yield and current yield could be attributed to new cultivars, technology spreading, and climate differences. Selective breeding has increased yield through the years with no yield plateau being discovered [62]. The introduction of machinery for ploughing, irrigation, and harvesting, doubled the agricultural capital in Denmark [63]. With drainage pipes and marling being introduced in the mid-19th century, land improved for agriculture, and production quantities grew by 75% [63]. Currently, about 50% of agricultural land is drained in Denmark [64], resulting in more area to be harvested and higher yields compared to historical yield potential. From 1930 to 1970, tile drains were mainly being placed in wetland areas, increasing the amount of arable land in parishes [65,66]. One qualifier to be the lowest class of soil fertility in the historical registers was waterlogged soil [23]. Draining waterlogged soils would cause these soils that once were given a low hard grain value to be given a much higher hard grain value. The loamy soils in southeastern Denmark are now highly productive areas due to the introduction of tile drains. The increase in CO 2 levels through the years can contribute to the increase in yield through time. The CO 2 level during the 18th century was around 280 ppm and current day is around 400 ppm [67]. The change in temperature is also a cause for differences between historical and current production. The temperature during 1875 to 1900 was around 2°C lower compared to average annual temperature data in Denmark from 2010 to 2017 [68]. Advances in crop genetics, the introduction of new technologies and management practices after the 1844 register, and higher ambient CO 2 levels and temperatures, are major reasons why the current yield values are larger than historical values.

Conclusions
The winter wheat yield map of Denmark generated in this study incorporates soil, climate, and topography to produce a more accurate yield prediction than previous efforts. Using a machine learning algorithm and adding climate data were major aspects that helped increase the variation explained in yield predictions compared to the previous statistical model. The previous effort, only using three topsoil properties to model winter wheat yield, was too simple to correctly predict yield. The random forest model did not predict extreme values well, but the values were centered around the mean. Even though extreme values will not be predicted, the new winter wheat yield map could be applicable for future yield predictions due to the centered values being more reliable predictions. The ease with which machine learning allows additional data also means the map can be updated more readily when new data are collected.
Comparing the current yield to historical yield potentials from the 17th and 19th centuries provided a way to determine the spatiotemporal trends in Danish agriculture. The differences in yield potential between the historical data are minor when comparing the current productivity to historical periods. The advances in agriculture have proven to out-perform what was once thought of as a field's yield potential. The agriculture potential will keep rising as technological breakthroughs continue to occur. Since predicting how technology and new cultivars will change yield potential is challenging, future work should use climate projections when creating yield maps to make the predictions viable in future years. The current yield map, while better than the previous soil-only model, will likely be obsolete in future years because of changes in climate. Thus, the trend in yield changes between time periods should be evaluated so that future yield predictions incorporate temporal trends, generating predictions that last longer into the future.