INTRODUCTION
Plant growth processes are influenced by soil, agricultural practices and climatic factors leading to significant variations of final yields. In rainfed production systems such as those that include winter wheat, climate variability is responsible for as much as 80% of the variability in the agricultural production (Hoogenboom, Reference Hoogenboom2000). Understanding how growth and development of rainfed wheat respond to the climate and especially to drought variability could provide early gross estimates of final yield production, which are of great importance for the agricultural policy and market planning, especially for countries in arid and semi-arid environments (Hoogenboom, Reference Hoogenboom2000; Simane et al., Reference Simane, Peacock and Struik1993). However, seasonal forecasts provide much less valuable information related to global trends and therefore assuming production only from a yield forecast may be quite simplistic.
Several techniques have been used for the prediction of pre-harvest crop yields including visual field estimates, multiple frame-based sample surveys, analog-year approaches, remote sensing, process-based simulation crop models and regression models (Becker-Reshef et al., Reference Becker-Reshef, Vermote, Lindeman and Justice2010). Regression models are among the most commonly used techniques for yield prediction due to the simpler structure and the use of fewer and more accessible environmental data. These models have the advantage that they capture both weather and management aspects of the yield variation, but their use may be limited on large spatial scales due to the inserted error, which is attributed to spatial variability of climate conditions (Kaufmann and Snell, Reference Kaufmann and Snell1997). Regression models that use meteorological parameters (mainly precipitation, radiation and temperature), soil conditions and agricultural practices as explanatory variables, have shown adequate capability to predict wheat yields (Lobell and Burke, Reference Lobell and Burke2010; Olesen et al., Reference Olesen, Bocher and Jensen2000; Wassenaar et al., Reference Wassenaar, Lagacherie, Legros and Rounsevell1999; You et al., Reference You, Rosegrant, Wood and Sun2009), while others focus particularly on the relationship between drought (in terms of water stress) and final yields (Hlavinka et al., Reference Hlavinka, Trnka, Semeradova, Dubrovsky, Zalud and Mozny2009; Mavromatis, Reference Mavromatis2007; Richter and Semenov, Reference Richter and Semenov2005; Stephens et al., Reference Stephens, Walker and Lyons1994; Yamoah et al., Reference Yamoah, Walters, Shapiro, Francis and Hayes2000).
The response of winter wheat yields to drought and water stress have been thoroughly investigated by several studies, which have indicated that the vegetative-reproductive stage is more sensitive than the seedling and the tillering stages, not only in arid and semi-arid areas but even in temperate and humid climate conditions (Austin et al., Reference Austin, Morgan, Ford and Blackwell1980; Blum, Reference Blum1998; El Hafid et al., Reference El Hafid, Smith, Karrou and Samir1998; Hlavinka et al., Reference Hlavinka, Trnka, Semeradova, Dubrovsky, Zalud and Mozny2009; Kimurto et al., Reference Kimurto, Kinyua and Njoroge2003; Simane et al., Reference Simane, Peacock and Struik1993). Combined drought effects during the seedling stage and the reproductive stage accounted for the highest yield reduction (El Hafid et al., Reference El Hafid, Smith, Karrou and Samir1998; Kimurto et al., Reference Kimurto, Kinyua and Njoroge2003). The occurrence of drought after sowing or during the seedling stage and the early growth stages lead to partial failure of seed germination and preservation of the young plants (Kimurto et al., Reference Kimurto, Kinyua and Njoroge2003). Drought occurrence at a later growth stage of the crop (i.e. in the reproductive stage) leads to lower respiration and photosynthetic rates, kernel shrivelling and lower ability to confront foliar diseases, which tend to spread and intensify towards and after flowering (Blum, Reference Blum1998). According to Kimurto et al. (Reference Kimurto, Kinyua and Njoroge2003), drought and consequently water stress during the reproductive stage was associated with lower number of seeds per head, increased number of sterile florets per head, reduced number of reproductive tillers, smaller ears, reduced number of spikelets per head and reduced kernel weights.
The objective of this study was to develop regression models for the prediction of grain yields of non-irrigated winter wheat in a semi-arid Mediterranean environment based on drought variability. Twenty-five years of yield data from four soils in northern Greece with different fertility–productivity degree and climate data covering the wheat-growing season were used for this purpose. The statistical resampling technique using subsets of the data (bootstrapping) was used to estimate the coefficients of the models, to assess the uncertainty of the selected explanatory variables and to validate the models. A further step was also performed to assess the positive or negative signal of yield gap (YG) between two growing seasons.
MATERIALS AND METHODS
Experimental design, yields and soil properties
Wheat (Triticum aestivum L.), cv. Yecora, was continuously grown for 25 years (1980–2004) without irrigation in four different soil-textured fields: a sandy loam (SL), a clay (C), a clay loam (CL) and a sandy clay loam (SCL) (Table 1) of the Aristotle University Farm of Thessaloniki in northern Greece (40o32΄Ν, 22ο59΄Ε) (Lithourgidis et al., Reference Lithourgidis, Damalas and Gagianas2006). Conventional tillage practices were used before sowing in each growing season. Sowing took place about mid-November in rows (spaced 16 cm apart) at a seeding rate of 150 kg ha−1. Nitrogen (N) at 120 kg ha−1 and P2O5 at 60 kg ha−1 as ammonium sulfo-phosphate (20–10–0) were incorporated into the soil before sowing each growing season. Weed control was achieved with appropriate herbicides registered for weed control in wheat. Each experimental field covered a 4-ha area. Wheat was harvested after mid-June and grain yield was adjusted to 13% grain moisture. Grain yield was determined by harvesting the total area of each experimental field and expressed as Mg ha−1. Straw was baled and removed after harvest and crop residues were incorporated into the soil.
Soil properties (0–30 cm) of the four experimental fields (Lithourgidis et al., Reference Lithourgidis, Damalas and Gagianas2006).
The summary of yields statistics over the 25-year period of the study for each soil is given in Table 2 (Lithourgidis et al., Reference Lithourgidis, Damalas and Gagianas2006). Differences among the mean yields of the four soils were observed, while the observed range (min–max: 0.85–4.72) covers the respective observed range on a national scale in Greece (Bakker et al., Reference Bakker, Govers, Ewert, Rounsevell and Jones2005). Taking into account that the four soils were under the same climatic conditions and were subjected to the same agricultural practices, the mean yield differences over an adequately long period of time (e.g. 25 years in this case) can be attributed to the intrinsic fertility–productivity of each soil. Considering the properties of each soil from Table 1 and the mean yields of Table 2, the soils that present low intrinsic fertility–productivity (i.e. the coarser SL and the finer soil C) showed higher yield variability according to the CV values (Table 2) than the soils with high intrinsic fertility–productivity (i.e. the SCL and the CL soil).
Summary of yield statistics (Mg ha−1) of the 25-year period (1980–2004) for each soil type (Lithourgidis et al., Reference Lithourgidis, Damalas and Gagianas2006).
*Different letters indicate statistically significant differences at p = 0.05.
SD: standard deviation; CV: coefficient of variation.
Agrometeorological parameters
Daily data of the incident solar radiation Rs (MJ m−2), relative sunshine hours n/N, mean relative humidity RH (%), precipitation P (mm), wind speed at two meters above the soil surface U2 (m sec−1), maximum temperature T max and minimum temperature T min (°C) were provided by the Hellenic National Meteorological Service for the period 1980–2004. The reference evapotranspiration ETo (mm day−1) was calculated using FAO-56 Penman–Monteith method (Allen et al., Reference Allen, Pereira, Raes and Smith1998). The mean monthly values of the agro-meteorological parameters are given in Table 3. A primary analysis of the monthly values of the agro-meteorological parameters and their combinations during the wheat-growing season from November to May was carried out to define the most important parameters, which are involved in final crop yield. June covers the last days of the maturity stage and for this reason it was not included in the calculations due to its insignificant effect on crop growth.
Mean monthly values of the agro-meteorological parameters for the period 1980–2004 in the study area.
T: temperature; P: precipitation; n/N: relative sunshine hours; RH: relative humidity; U2: wind speed; Rs: incident solar radiation; ETo: reference evapotranspiration.
The parameter of precipitation versus reference evapotranspiration (P/ETo) ratio was found to be the most crucial parameter, indicating the strong effects of water availability in final wheat yields under the specific climatic conditions. For this reason, this parameter was selected to be introduced in the models. This ratio has been used as a basic parameter in drought indices (Dubrovsky et al., Reference Dubrovsky, Svoboda, Trnka, Hayes, Wilhite, Zalud and Hlavinka2008; Tsakiris and Pangalou, Reference Tsakiris, Pangalou, Iglesias, Garrote, Cancelliere, Cubillo and Wilhite2009) and indirectly expresses the soil water regime for crop growth including all the meteorological parameters through ETo. Moreover, this ratio was adopted by UNEP (1992) for the classification of climatic environments according to drought classes. The annual value of P/ETo is 0.46 in the region (using the last row data of Table 3). It ranges between 0.2 and 0.5 and corresponds to a semi-arid climatic environment according to UNEP (1992).
For the elimination of normality departures, the parameter of P/ETo before its inclusion in the prediction models was transformed using the following equation:
$\begin{equation}
x_j = \ell n [ {1 + ( {P_j /ET_{oj} } )} ],\end{equation}
$Statistics of the transformed monthly precipitation versus reference evapotranspiration ratios for the wheat-growing season (period 1980–2004).
*Different letters indicate statistically significant differences at p = 0.05.
SD: standard deviation; CV: coefficient of variation.
Model development
The transformed ratios of P/ETo and the mean yields of each soil were used as independent variables in multiple regression analysis for the development of yield prediction models. The general form of the multiple regression models is the following:
$\begin{equation}
y_i\, {=}\, b_o\, {+}\, \sum\limits_{j\, {=}\, 1}^k {b_j x_{j\,i} }\, {+}\, b_{k\, {+}\, 1} \cdot \bar y\, {+}\, e_i ,\quad {\rm where}\quad x_{j\,i}\, {=}\, \ln [ {1\, {+}\, ( {P/ET_o } )_{j\,i} }]\quad {\rm and}\quad k\, {=}\, 7,\end{equation}
$A preliminary multiple regression analysis of equation (2) on the full dataset was performed using the least squares to assess the statistical significance (p-value from t-test) of each explanatory variable (Table 5). The analysis was performed with full variables selection and with stepwise backward selection to exclude the less significant explanatory variables.
Average values and standard deviations of the regression coefficients estimated by the bootstrapping procedure and results of the statistical tests from the validation procedure.
*t-test results using the full dataset.
†Results of the bootstrapping procedure (N = 10,000).
‡The values outside parenthesis are the statistical tests of the validation procedure using the observed mean yield as explanatory variable, while inside parenthesis using the modelled mean yield obtained by the clay function of equation (9).
Statistical techniques of resampling using subsets of the initial data such as jackknife, cross-validation and bootstrapping are common techniques to assess the uncertainty of the selected explanatory variables and to validate models by using subsets of the initial data. In jackknife and cross-validation, the resampling is based on pre-designed selection schemes, while bootstrapping is based on the generation of a large number of new datasets from the initial one by randomly sampling data with replacement (Efron and Tibshirani, Reference Efron and Tibshirani1993). In our study, the bootstrapping technique was selected and applied in the same way as described in a previous similar study by Prost et al. (Reference Prost, Makowski and Jeuffroy2008) using the ‘glm’ function of the R statistical software (Venables and Ripley, Reference Venables and Ripley2002). The procedure was performed for two cases of equation (2), with full variables (model I) and with reduced variables (model II) according to the statistical significance (p-value from t-test) of each explanatory variable. All observations (N = 100) were used to generate 10,000 bootstrap samples of N size from the initial dataset in both cases. The procedure produced 10,000 estimations of each coefficient of equation (2), while their average values were used in the final form of the two models.
Statistical tests were used to validate the prediction accuracy of the final models using the full dataset. These tests were: the correlation coefficient (R), the mean bias error (MBE), the root mean square error (RMSE), the mean absolute error (MAE), the variance of the distribution of differences that expresses the variability of (C–O) distribution from MBE (s2d) and the index of agreement (d), which are given by the following equations:
$\begin{equation}
R = \frac{{\sum\limits_{i = 1}^N {( {C_i - \bar C} )( {O_i - \bar O} )} }}{{\sqrt {\sum\limits_{i = 1}^N {( {C_i - \bar C} )^2 } } \cdot \sqrt {\sum\limits_{i = 1}^N {( {O_i - \bar O} )^2 } } }},\end{equation}
$
$\begin{equation}
{\rm MBE} = \frac{1}{N}\sum\limits_{i = 1}^N {( {C_i - O_i } )} ,\end{equation}
$
$\begin{equation}
{\rm RMSE} = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {( {C_i - O_i } )^2 } } ,\end{equation}
$
$\begin{equation}
{\rm MAE} = \frac{1}{N}\sum\limits_{i = 1}^N {| {C_i - O_i } |} ,\end{equation}
$
$\begin{equation}
s_d^2 = \frac{1}{{( {N - 1} )}}\sum\limits_{i = 1}^N {( {C_i - O_i - {\rm MBE}} )^2 } ,\end{equation}
$
$\begin{equation}
d = 1 - \frac{{\sum\limits_{i = 1}^N {( {C_i - O_i } )^2 } }}{{\sum\limits_{i = 1}^N {( {| {C_i - \bar O_i } | + | {O_i - \bar O_i } |} )^2 } }},
\end{equation}
$RESULTS AND DISCUSSION
Evaluation of the models
The general model of equation (2) was used in the bootstrapping procedure using full variables selection (model I) and reduced variables selection by excluding the effects of January and May (model II) according to the statistical significance of the explanatory variables estimated by the p-value (Table 5). The observed mean yields (Table 2) were used in the term y. The average values of the regression coefficients and their standard deviation according to the bootstrapping procedure are given in Table 5. The comparisons between observed and predicted yields with both models are given in Figures 1a and b. The statistical tests from the validation procedure are also given in Table 5 (the values outside parenthesis), where both models indicated adequate and similar prediction accuracy. The exclusion of January and May in the second model had very little effect on the prediction accuracy. The residuals of both models did not indicate any serial autocorrelation (Box–Pierce test) and they were normally distributed (Shapiro–Wilk test) at 95.0% confidence level. These results are in accordance with findings of previous studies on water stress and drought responses of wheat cultivars (Austin et al., Reference Austin, Morgan, Ford and Blackwell1980; Blum, Reference Blum1998; El Hafid et al., Reference El Hafid, Smith, Karrou and Samir1998; Hlavinka et al., Reference Hlavinka, Trnka, Semeradova, Dubrovsky, Zalud and Mozny2009; Kimurto et al., Reference Kimurto, Kinyua and Njoroge2003; Simane et al., Reference Simane, Peacock and Struik1993). Moreover, drought effects at the mid-end of the ripening stage of wheat (mid-end of May) were found insignificant even though the lowest values of P/ETo were observed during this stage.
Observed versus predicted yields for (a) the model I (full variables) using the observed y, (b) the model II (reduced variables) using the observed y, (c) the model I using the y(clay) and (d) the model II using the y(clay).
For detailed predictions at farm level, the term of mean yield y can be derived by the farmers’ experience or better by the use of the farmers’ yield diary, if this exists. A more generalized approach is to describe this term as a function of soil properties, which needs a long-term experiment with huge amount of laboratory measurements and a respective amount of sampling sites. Here, an attempt is made to describe y as a function of specific soil properties using the data of this study for the site-specific environmental conditions. Taking into account the above considerations, a multiple variable analysis (linear correlations) was performed on the parameters of Table 1. The results showed that soil clay content presented the higher number of statistically significant non-zero correlations at 95.0% confidence level with the other parameters and especially (a) with the cation exchange capacity, which indicates the ability of soils to absorb and retain nutrients and (b) with those which are related with the soil hydraulic properties (water-holding capacity, permanent wilting point and available water). These results set this parameter as the most representative in describing other soil properties and consequently the term y. Using the mean yields (Mg ha−1) of the four soils (Table 2), y was expressed as a function of the clay content (g 100 g−1) (Figure 2):
$\begin{eqnarray}
&& \bar y({\rm clay}) = \exp \left[ {7.595 + \frac{{ - 37.59}}{{{\rm clay}}} - 1.486 \cdot \ell n\left( {{\rm clay}} \right)} \right] \nonumber \\
&& \quad ({\rm R} = 0.99)\quad (11.3 < {\rm clay} < 59.2).\end{eqnarray}
$Mean yield y as a function of soil clay content and comparison with the yield data of the four soils (mean values and standard error).
Yield gap
The YG between two growing seasons can be predicted using the general form of equation (2) by abstracting the respective modelled yields:
$\begin{equation}
{\rm YG}_{(i,i - z)} = y_i - y_{i - z} \approx \sum\limits_{j = 1}^k {b_j x_{j\,i} } - \sum\limits_{j = 1}^k {b_j x_{j\,i - z} } ,\end{equation}
$Equation (10) is free from the term y and estimates the mean YG between two growing seasons for the selected soils, which were used for the calibration of equation (2). The YG is of primary importance in operational policy and market planning, where yield predictions need to be performed for huge agricultural lands.
The coefficients of model I were incorporated in the second part of equation (10), which was used to predict all the cases of YG between the years of the period 1980–2004 (e.g. for 2004, YG was calculated between 2004 and 2003, 2004 and 2002,. . .etc. and the same procedure was repeated for 2003, 2002. . .). The predicted YG signals were compared with the mean observed (yi − yi − z) signals of the four soils from the respective years (Figure 3). Figure 3 was divided in four parts to optimize the signals presentation, where YG function presented 84% predictive accuracy to identify the correct signal (+ or −) and adequate performance to predict the mean YG of the four soils (Figure 3).
Predicted YG versus the mean observed YG (mean values and standard error) of the four soils between all the years of the period 1980–2004.
CONCLUSIONS
The results of this study indicated that regression models that use the monthly ratio of precipitation versus reference evapotranspiration and the term of intrinsic mean yield y, which indirectly describes the effects of soil properties, can adequately predict grain yields of non-irrigated winter wheat in a semi-arid Mediterranean environment. The results of the models identified adequately the effects of monthly drought variability on the final yields. Firstly, the most crucial period was found to be the vegetative-reproductive stage (i.e. February to April), and secondly the seedling stage and the early growth stage (i.e. November to December). Modelling the term of intrinsic mean yield y under the site-specific conditions as a function of specific soil properties revealed that soil clay content was the most representative parameter in describing other soil properties and consequently the term y of the soils. The basic conclusion of this attempt was that yield is usually non-linearly correlated with most of the soil properties, indicating that future attempts for yield estimations using multiple regression models need to describe the explanatory variables related to soil properties as complex non-linear functions. With the proposed models, YG predictions between two growing seasons of the selected soils were adequately accurate and the prediction of the mean YG (free of the term mean yield) showed adequate performance.
Acknowledgements
The authors are thankful to the administration of the Hellenic National Meteorological Service for its contribution to this study by providing the meteorological data and to the biometrician Andrea Benazzo from Ferrara University (Department of Biology and Evolution) in Italy for his contribution to improve the statistical analysis.
