#### Abstract

Farmers' decision to adopt new management or production system depends on production risk. Grain yield data was used to assess production risk in a field experiment composed of two cropping systems (CNV and ORG), each with eight subsystems (two levels each of crop rotation (2-yr and 4-yr), tillage management (conventional, CT and strip, ST), and fertilizer input (fertilized, YF and non-fertilized, NF)). Statistical moments, cumulative yield (CY), temporal yield variance (TYV) and coefficient of variation (CV) were used to assess the risk associated with adopting combinations of new management practices in CNV and ORG. The mean-variance-skewness (M-V-S) statistics derived from yield data separated all 16 subsystems into three clusters. Both cropping systems and clustered subsystems differed as to their ability to maintain a constant yield over years, displayed different yield cumulative probabilities, exhibited significant and different M-V-S relationships, and differed as to the reliability of estimating TYV as a function of CY. Results indicated that differences in management among cropping systems and subsystems contributed differently to the goal of achieving yield potential as estimated by the cumulative density function, and that certain low-input management practices caused a positive shift in yield distribution, and may lower TYV and reduce production risk.

#### 1. Introduction

Production risk influences farmers’ decision to adopt a new management practice or a production system [1]; therefore, the sustainability of cropping systems is becoming increasingly important to farmers and researchers alike [2]. Although management systems with reduced chemical and energy inputs are of particular interest in assessing sustainability [3, 4], some of them can be less stable than others depending on the cropping system in question [2, 4]. Yield instability, whether caused by spatial and/or temporal variation, can be identified and measured based on performance of long-term experiments [5]. Farmers may not be able to fully manage spatial variation if temporal variation is a strong and recurring factor [6]. Temporal yield variation can be managed to some extent by making the right management decisions at the right time; however, it is equally important for famers, through investment and knowledge, to predict and control this variation [6, 7].

Quantifying treatment main effects in cropping systems experiments provides valuable information that can be augmented by examining the interaction between years and treatment [8] through which cumulative effects of treatments can be measured and compared; the latter can help unmask the real effects of fixed factors. Statistical analysis of long-term cropping systems experiments is generally more complicated than that of short-term experiments [4]; it requires a specific methodology of data processing and of biometric analysis. Principal components analysis (PCA), partial least squares (PLS) regression, and restricted (or residual) maximum likelihood (REML) are becoming increasingly used in analyzing and interpreting results of long-term experiments that include fixed and random factors, categorical and continuous variables, and spatial and temporal variation [4, 9, 10] The statistical analysis procedure of repeated measurements of annual crop yield during the course of a cropping systems experiment using REML is widely used [9, 10] to account for the serial correlation (i.e., variance-covariance) among measurements taken on the same experimental unit over time. Also, repeated measurement analysis provides models which are more parsimonious than the one implied by multivariate analysis of variance (MANOVA), accommodates several error terms, and can easily handle incomplete designs [10].

Crop yields and their temporal variances are influenced by management factors, especially crop rotations [4, 11]. However, cropping systems can have larger influence on the temporal variability and stability of crop yields [12]. Therefore, it is more appropriate to compare the performance of cropping systems and not just to quantify the effect of single factors in a cropping systems experiment [2, 4]. The objectives of this long-term study were to quantify multifactor and multivariate impact on the statistical determinants of the cumulative density function of total yield in conventional, organic cropping systems and their subsystems and model the impact of cumulative yield on temporal yield variance as influenced by management practices that can be easily implemented on the farm.

#### 2. Materials and Methods

##### 2.1. Field Experiment

A long-term cropping systems experiment was established in 2002 on a land area of about 3.0 ha as a split-plot randomized complete block design with four replications (i.e., blocks). The experimental site was uniformly cropped for one year prior to the start of the experiment. The predominant cropping system practiced by farmers in this part of the upper Midwest of the United States (45°41′N, 95°48′W at 370 meters above sea level) is based on conventional management of a 2 yr corn-soybean crop rotation using conventional tillage and external inputs (fertilizer, herbicides, etc.). Sixteen subsystems were formulated as combinations of two levels each of cropping systems (conventional and organic), crop rotation (2-yr of corn (*Zea mays* L.) and soybean (*Glycine max* (L.) Merr.), and 4-yr of corn-soybean-wheat (*Triticum aestivum* L.)-alfalfa/alfalfa (*Medicago sativa* L.)), tillage (conventional and strip tillage), and N fertilizer (fertilized and nonfertilized).

Fertilizer rates (inorganic N source for the conventional and animal manure N source for the organic cropping system) were determined for each crop on the basis of annual soil analyses and regional N recommendations. The Nitrogen Decision Aid software (http://www.ars.usda.gov/services/software/download.htm?softwareid=85) was used to determine annual N fertilizer rates for each crop and cropping system. Two full rotation cycles of the 4-yr crop rotation and four full rotation cycles of the 2-yr crop rotation were completed during the course of the experiment. All phases of each crop rotation were present each year; therefore, interaction with the environment was not a concern because all crops were exposed to the same annual environmental fluctuations and over the duration of the experiment. Details of the experimental design, crops, crop rotations, management practices, and inputs are presented elsewhere [4].

##### 2.2. Measurements

Grain yield in each of eight years (2002–2009) was measured from a central 15 m^{2} mechanically harvested strip per plot of corn, soybean, and wheat and adjusted to a moisture content of 15.5, 13.0, and 13.5%, respectively. Total dry matter yield of alfalfa was measured on two 0.5 m^{2} subplots per plot harvested three times per year, and adjusted to a moisture content of 15.0%. Measured yield data were expressed in Mg ha^{−1} for each subsystem and then used for further statistical analyses and modeling. Annual and total yields were used to calculate annual and cumulative yield means (years 1 to 8), medians, temporal yield variances, coefficients of variation, skewness, and kurtosis and to construct yield frequency distributions for each cropping system. The M-V-S estimates for subsystems were used to cluster them into three distinct groups (Clusters), based on results of a preliminary principal components analysis (see below).

##### 2.3. Statistical Analyses

Principal Components Analysis (PCA) was used to identify possible associations (i.e., loadings on the first PC) and quantify the amount of total variation accounted for by statistical descriptors, central moments, two cropping systems (Conventional and organic), their components (i.e., crop rotations, tillage, and fertility management options), and the resulting 16 subsystems. The correlation matrix between variables was used in the analyses to give variables equal weight after standardizing the variances to a value of 1. The results of PCA were used to classify the 16 subsystems into those with significant positive, negative, or no significant loadings on PC1. Two demarcation lines at the ±0.2 loadings separate these subsystems (Figure 1).

Repeated measurements analysis using REML with first-order autoregressive (AR1) models appropriate for equally spaced data (i.e., annual yield estimates) and an option to allow for heterogeneous variances over years, was used on annual yield data of all subsystems [9]. A linear mixed model of the form, was used to perform the analyses; where is a overall mean, is the random effect of block ; is the fixed main effect of treatment (i.e., subsystem) ; is the random plot error for the experimental unit in block receiving treatment ; is the fixed main effect of cycle ; represents the random effect that is associated with yield measurement taken in the th year of the study; is the fixed interaction effect that is associated with treatment and cycle ; represents the random interaction effect for treatment in the th year; is the random error that is associated with the yield measurement at year on the experimental unit in block receiving treatment . A cycle is defined as the factor representing the fixed effect of years.

##### 2.4. Modeling

Two modeling approaches were employed to evaluate and compare yield and temporal variance of both cropping systems and those of clustered subsystems; these were cumulative probabilities of obtaining a given yield level based on empirical cumulative density function (CDF) and on the normal distributions which were estimated using curve fitting. Empirical cumulative density functions were fitted to yield data, and the ±95% confidence ellipsis was constructed using a kernel density estimation procedure [9]. The central moments generated by curve fitting, in addition to statistical descriptors describing yield histograms, and the cumulative probability (labeled as *y2-y1*) of obtaining total yield within the category boundary (labeled as *x2-x1*) with the highest frequency, were used to compare and contrast the conventional and organic cropping systems and the clustered subsystems; Calibration PLS regression models were developed using 75% of cumulative yield data points to predict temporal yield variance as a measure of yield stability over time for each of the cropping systems and the clustered subsystems; a validation procedure was performed by deploying the calibration model to predict the remaining 25% data points [9]. Temporal yield variance was predicted (and validated) as a function of cumulative yield for each cropping system and for subsystems clustered in three groups. A model of the form,
was used to perform the analyses; where is a matrix of explanatory variables (i.e., cumulative yield data) for temporal yield variance given by the vector in (3), are -dimensional vectors called -loadings and is the residual matrix, and
where ( to 8 years, called the latent variables) and the are called the -loadings. In this model, the dependency among the -explanatory variables is broken up, and the relationship between and is transmitted through the latent variables .

Cumulative yields with significant regression coefficients were reported in the validation PLS models (4) to (8). The residual mean squares error (RMSE) and the coefficient of determination of the validation model () were reported for each PLS model as indicator of its significance and measure of its reliability, respectively. Relevant modules in GenStat Version 10 [9] were used in data processing, statistical analyses, and modeling.

#### 3. Results

##### 3.1. Total Variation

Principal components analysis captured a large portion () of total variation in one PC, and the validation model accounted for 0.74 of that variation (Figure 1). Differences between the conventional and organic cropping systems were the largest based on their loadings on PC1, followed, in decreasing order, by differences between fertilizer, tillage, and crop rotation treatments. Loadings of the sixteen subsystems on PC1 ranged from maximum positive (Conventional system, 2-yr rotation, conventional tillage, with N fertilizer; C2CTYF) to maximum negative (Organic system, 2-yr rotation, strip tillage, without N fertilizer; O2STNF). Most subsystems with the YF input loaded positively and most subsystems with the strip tillage nonfertilized- (STNF-) input loaded negatively on PC1. Loadings of statistical descriptors and central moments were clustered in two groups (Figure 1): mean yield and its descriptors, had positive loadings, while those describing the shape of yield distribution (i.e., skewness and kurtosis) as well as the coefficient of variation (CV) and RMSE had negative loadings on PC1. The eight subsystems within the ±0.2 bands, as well as crop rotations and the residual mean squares error (RMSE), explained the remaining portion () of total variation in PC2.

##### 3.2. Mean-Variance-Skewness (M-V-S) Analysis

The three-way relationship between S and the other two central moments (M-V) was positive and significant (r: S/M-V = 0.75; ; Figure 2(a)); however, the relationship was smaller in magnitude (0.45; ) when based on total yield and its variance and skewness estimates (Figure 2(b)). The clustering algorithm resulted in a complete separation between the three clusters in the first (Figure 2(a)) but not the second case (Figure 2(b)). Cluster 1 was composed of four conventional subsystems, Cluster 3 was composed of five organic subsystems, and the intermediate Cluster 2 was composed of a mixture of four conventional and three organic subsystems. The decreasing level of skewness (Figure 2(a)) was associated with increasing levels of both mean and variance suggesting that larger mean yields may tend to reduce skewness of yield distribution in spite of the increased variance. However, total yield seemed to have a much larger buffering capacity than mean yield thus resulting in a relatively lower three-way correlation (0.45; Figure 2(b)).

**(a)**

**(b)**

The magnitude and level of significance of bivariate M-V-S correlation coefficients varied between cropping systems and clusters within systems (Table 1). Significant positive -values (0.37 to 0.90) were found between mean and variance values across cropping systems and clusters, and significant and negative correlation coefficients were found between skewness and each of mean and variance; the former (−0.58 to −0.89) had larger magnitude and stronger level of significance than the latter (−0.33 to −0.74).

##### 3.3. Repeated Measurements Analysis

Results of repeated measurements analyses are summarized in Table 2. The -values (i.e., ratio between variance due to a fixed factor and that of the experimental error) and their level of significance illustrate the dynamic impact of years, subsystems, and their interactions. The steady increase in -values was largest for years followed, in decreasing order, by subsystems and their interaction. Largest changes in -value for years were observed at the beginning of the second rotation cycle (i.e., year 5) and at the last year of the experiment (year 8). This change coincided with similar changes, but smaller in magnitude, in -values for subsystems and their interaction with years. The magnitude of -value for the interaction between years and subsystems lagged behind, and it became of statistical significance at year 5. In comparison, -values for subsystems increased at a faster pace than those for years and became almost equal at the end of the experiment (year 8).

##### 3.4. Cumulative Density Function (CDF) and M-V-S of Total Yield

Most descriptive statistics and central moments differed significantly between conventional and organic cropping systems (Figures 3(a) and 3(b)). During the course of the experiment, conventional and organic cropping systems produced a total of 3746.18 and 2550.3 Mg of dry matter, respectively, with the respective means of 39.0 and 26.57 Mg ha^{−1}. Similarly, minimum and maximum yields of the conventional cropping system were much larger than those of the organic. The largest frequency in conventional cropping system was for the 35–40 Mg ha^{−1} yield category, followed by 30–35 Mg ha^{−1}; in the organic cropping system, it was 25–30 Mg ha^{−1}, followed by the 20–25 Mg ha^{−1} yield category. The distribution curve based on total yield of the organic cropping system was less skewed (0.098) than that of the conventional cropping system (0.27) which was positively skewed and with less kurtosis (−0.22 versus −0.79) although their respective medians were very close to their means. The larger mean yield of the conventional cropping system was associated with larger variance (68.4 Mg ha^{−1})^{2 }as compared to the organic cropping system (34.97); however, their long-term coefficient of variation values (21.19 and 22.26%, respectively) was almost equal.

**(a)**

**(b)**

Four subsystems, all belonging to the conventional cropping system, clustered together based on their M-V-S values to form Cluster 1 (Figure 4(a)). The fertilizer input (YF) was the factor common to all subsystems. The remaining factors (2-yr and 4-yr crop rotations, and conventional and strip tillage) had, more or less, the same impact on statistical descriptors and central moments of this cluster. The long-term mean-variance (M-V) values (45.0 and 37.15, resp.), but not skewness (), were the largest among all three clusters and were associated with the smallest coefficient of variation (CV = 13.53%). Yield distribution curve was slightly positively skewed, with the mean and median being almost equal. The largest frequency was found for the 45–50 Mg ha^{−1} yield category, followed, in decreasing order, by 40–45 and 35–40 Mg ha^{−1} yield categories. The second cluster (Figure 4(b)) was composed of seven subsystems, four of which belong to the conventional, and the remaining three belong to the organic cropping systems. The M-V-S values (32.0, 24.39, and 0.52, resp.) were intermediate, except for the skewness which was the largest among clusters. The nonfertilized (NF) and fertilized (YF) inputs were the common factors among the conventional and organic subsystems, respectively. The impact of rotation and tillage factors on statistical descriptors and central moments of this cluster was more or less the same. The largest frequency was found for the 25–30 Mg ha^{−1} yield category, followed, in decreasing order, by the 30–35 and 35–40 Mg ha^{−1} yield categories. The third cluster (Figure 4(c)) was composed of five organic subsystems and was characterized by the smallest M-V-S values (23.37, 21.0, and −0.12, resp.). This is the only cluster with a negative skewness similar to the organic cropping system. All subsystems, except O2STYF, share the nonfertilized (NF) input which may have contributed to the smallest yield and the largest coefficient of variation (CV = 19.62%) among clusters. Similar to Cluster 2, the impact of rotation and tillage factors on statistical descriptors and central moments of Cluster 3 was more or less the same. The largest frequency was found for the 20–25 Mg ha^{−1} followed by 25–30 Mg ha^{−1} yield category.

**(a)**

**(b)**

**(c)**

##### 3.5. PLS Validation Models

Temporal yield variance (TYV) of the conventional cropping system as a function of cumulative yield (CY) during eight years was predicted and validated by the following PLS model: The validation model accounted for 0.71 of total variation in temporal yield variance and included cumulative yield at four time points (i.e., years) during the course of the experiment. The first (CY1) had a negative, albeit very small regression coefficient ().

The remaining cumulative yields had positive impact on temporal yield variance. Except for the small (0.088) of CY6, the remaining two were relatively high. The largest positive impact on temporal yield variance was due to CY7 () which represents cumulative yield at the end of the experiment.

Temporal yield variance of the organic cropping system as a function of cumulative yield during the course of the experiment was predicted and validated by a simpler PLS model, as compared to temporal yield variance of the conventional cropping system: The validation model accounted for 0.45 of total variation in temporal yield variance and included cumulative yield at three time points during the course of the experiment. The first two (CY1 and CY2) had large and positive ( and 0.183, resp.) and the last (CY7) had negative () impact on TYV.

Temporal yield variance for all four subsystems (within the conventional cropping system) in Cluster 1 () was modeled as This model has a smaller and a slightly larger RMSE (1.56) as compared to the full model developed for the conventional cropping system. All seven cumulative yield measurements contributed to this model with a mixed positive (CY2, CY3, CY6, and CY7) and negative (CY1, CY4, and CY5) regression coefficients. A general trend of increasing positive impact over the years can be deduced from this model.

Only a small (), but significant portion of temporal yield variance of the heterogeneous subsystems in Cluster 2, was accounted for by the validation model This model is characterized by the largest RMSE (1.79) as compared to the remaining models. All regression coefficients were positive, and, except for CY3 (), they were small in magnitude. Finally, an intermediate portion of temporal yield variance () was accounted for by all cumulative yields of subsystems in Cluster 3 (all within the organic cropping system) A mixture of positive, large in magnitude (CY1, CY2, and CY5), and negative, but smaller in magnitude (CY4, CY6, and CY7) regression coefficients characterized this model and resulted in RMSE (1.55) smaller than the one (1.64) estimated for the full model of the organic cropping system.

#### 4. Discussion

Crop yield density estimation has been the subject of extensive empirical research, especially in risk analysis [1]. Numerous studies have been conducted [1, 13] arguing for the existence of nonnormality, skewness, and kurtosis in crop yields. On the other hand, there is little consensus on the impact of environmental variables on empirical distribution of crop yields [14, 15]. The current study was designed to allow each cropping system to come as close as possible to achieving its potential yield and to minimize, if not neutralize, any differential environmental effects on yield variation [6, 15]. Nevertheless, the statistical analyses of yield data and its secondary statistics from a long-term experiment present a number of challenges, including the possible temporal and spatial correlation, and the high probability of heterogeneity in the error variance across years [4].

Farmers are mostly risk averse; they are primarily interested in knowing the probability of getting low yields that would lead to estimating risk [1, 14]. A simple visual and quantitative measure of predicting low yield was provided in this study by the one-dimensional plot of all factors and central moments in the experiment (Figure 1). Negative loadings of a cropping system, subsystem, or input are indicative of lower yields. The negative loading of a central moment, when associated with a particular factor (i.e., system, subsystem, or input), is indicative of the magnitude of departure from potential yield or departure from its normal distribution. In addition, although maximum yields produced by the conventional and organic cropping systems were 48% and 50% higher than their respective means, their lower quartiles (LQs, as an additional measure of production risk) were only 18% and 16% smaller than their respective means (Figure 3); the LQ is expected to be produced by each cropping system at a cumulative probability of 25%.

The bivariate (Table 1) or three-way variance-covariance structures of M-V-S characteristic of each cropping system or clustered subsystems, although similar in sign (positive or negative), were different in magnitude and in their level of significance. These differences may suggest that similar management options among cropping systems contributed differently to the goal of achieving their yield potential, triggered changes in the bivariate and three-way M-V-S relationships, and influenced the yield distribution and its cumulative probability. The small skewness, kurtosis, and coefficient of variation values are a result of combining all crop yields for each cropping systems, or subsystem over time; however, statistical differences persisted due to differences between systems or between subsystems in response to management and input factors. Any input or management practice that reduces positive skewness, especially when associated with small coefficient of variation as a measure of static stability of yield over years, is important for the sustainability of any cropping system. Such input or management practice is expected to increase yield and to create a more favorable risk environment [16, 17]. In the graphical analyses (Figures 3 and 4), the intersection between the normal and empirical cumulative density function curves indicates the point at which positive skewness passes over to negative skewness, thus triggering a shift towards the negative side of the yield distribution curve, a change in yield distribution, and lower sustainability.

The design of this long-term experiment and the statistical analyses took into consideration the random temporal variation and prevented it from masking the real effects of fixed factors (i.e., subsystems, in this study) [6, 8]. Additionally, repeated measurements analysis allowed for heterogeneity across years and provided statistical evidence of the relative importance of years, subsystems, and their interaction (Table 2) [9, 18]. The interaction between years and subsystems became significant and persisted until the end of the experiment, when the magnitude of difference between subsystems became larger than the magnitude of differences between years. It can only be speculated that differences late in the experiment (e.g., year 6) are due to differences in cumulative (e.g., years 1–5) effects and possibly by the order in which they occurred [8, 18]. Posner et al. [17] noted that the effect of previous crops in a crop rotation would not be reflected fully in crop yields until the 3rd year in a 2-yr crop rotation and, as demonstrated by the current study, until the 5th year in a 4-yr crop rotation.

Posner et al. [17] questioned whether biologically diverse, low-input cropping systems can be as productive as conventional, high-input systems? Organically managed cropping systems have been reported to be less productive than conventional systems [3, 4]. Crop yields in the conventional cropping system may be more prone to temporal variation [4, 19] as indicated by its larger variance estimate (68.4; Figure 3(a)) as compared to the much lower variance estimated of the organic cropping system (34.97; Figure 3(b)). Large temporal variance values may result in larger but not stable total yields. In this study, yields of unfertilized controls still exceed 50% of the fertilized treatments after six years; nevertheless, average total yield of the organic cropping system was 70% of that of the conventional after eight years.

The reliability of using cumulative yield to predict temporal yield variance, although resulted in different model fit () and RMSE values for systems and subsystems (see (4)–(8)), can be attributed to eliminating the confounding effects of the differential magnitude of the means caused by differences among crops, years, management factors, and their interactions [16]. These estimates provided a powerful indicator of the cumulative influence of multiple factors on total crop yields during eight years. Additionally, the PLS regression models used to predict temporal yield variance, unlike ordinary multivariate linear regression models, took into consideration the variation in the predictors (i.e., cumulative yields), the predicted variable (i.e., temporal yield variance), and their joint variation [20]; therefore, PLS is considered more useful than ordinary multivariate linear regression in obtaining more parsimonious models for predicting yield variation. Depending on the cropping system or clustered subsystems, PLS models were effective in detecting the “variables” (i.e., cumulative yield over time) that explained up to 71% of total variation in temporal yield variance [9, 20].

This and a few other similar studies [12] demonstrated that cropping systems can influence the temporal variation and stability of crop yields; however, the multi level statistical analyses in this study provided reliable quantitative estimates of the deviation of modeled temporal variances from measured ones [13] and eliminated the confounding effects of the differential magnitude of the means caused by differences among crops, years, management factors, and their interactions. Furthermore, the multi level statistical analyses helped identify subsystems with high but less temporally stable yield (C2CTYF; Figure 2(a)), intermediate yield and temporal yield variance (e.g., O2CTYF and C4CTNF), and low but stable yield (O2STNF) within each of the conventional and organic cropping systems. Although inferences, based on the current analyses, will be restricted to the particular eight years during which the experiment was conducted, it will be possible to make general inferences when the variability due to sequences of environmental characteristics is quantified and accounted for in the statistical analyses [8]. Nevertheless, the data is invaluable in agronomic modeling studies and will be used in conjunction with long-term historic and simulated weather data [14] to predict crop yield and its temporal variation.

#### 5. Conclusions

Multivariate statistical analyses of annual and cumulative crop yields during the course of a long-term experiment identified stable and productive subsystems within each of conventional and organic cropping systems. The study provided evidences indicating that differences in management among cropping systems contributed differently to the goal of achieving the yield potential of each cropping system. Several statistical descriptors, including mean-variance-skewness of cumulative crop yields, were used to construct cumulative density functions and provided reliable indicators of yield stability in a long-term cropping systems experiment when subjected to repeated measurements, principal components, and partial least squares regression analyses. Results of the study may help farmers, agronomists, and crop consultants identify novel components of cropping systems that can be implemented on the farm with less external inputs and may result in reducing temporal variation of crop yields.