Abstract
Guidelines are needed to develop proper statistical analyses procedures and select appropriate models of covariance structures in response to expected temporal variation in long-term experiments. Cumulative yield, its temporal variance, and coefficient of variation were used in estimating and describing covariance structures in conventional and organic cropping systems of a long-term field experiment in a randomized complete block design. An 8-year database on 16 treatments (conventional and organic cropping systems, crop rotations, and tillage) was subjected to geostatistical, covariance structure, variance components, and repeated measures multivariate analyses using six covariance models under restricted maximum likelihood. Differential buildup of the cumulative effects due to crop rotations being repeated over time was demonstrated by decreasing structured and unstructured variances and increasing range estimates in the geostatistical analyses. The magnitude and direction of relationships between cumulative yield and its temporal variance, and coefficient of variation shaped the covariance structures of both cropping systems, crop rotations, and phases within crop rotations and resulted in significant deviations of organic management practices from their conventional counterparts. The unstructured covariance model was the best to fit most factor-variable combinations; it was the most flexible, but most costly in terms of computation time and number of estimated parameters.
1. Introduction
Long-term experiments (LTEs) that include fixed and random factors are valuable tools in understanding the effects of spatiotemporal variation and in developing guidelines for proper management practices [1, 2]. Moreover, LTEs can delineate risks and stability of cropping practices due to year-to-year variability in biotic and abiotic stresses [3]. A major advantage of an LTE over short-term experiment is that it provides insight into the causes of the changes in the slope of the responses, the causes of the inflection points, and the magnitude of the long-term change in crop yield and its variance, whereas a short-term experiment focuses only on the initial trajectories of variables under study [4]. In addition, there is more interest in examining not only treatment main effects but also, more importantly, the treatment × time interaction through which cumulative effects of treatments can be quantified and contrasted [5].
Crop yield in LTEs with repeated measures involving diverse cropping systems and crop rotations can be expressed as a function of large-scale deterministic structure and small-scale stochastic structure [6]. Such spatial structures can be confounded by strong and recurring temporal variation [7, 8] in which case farmers may not be able to fully and sustainably manage spatial variation [7]. Statistical analysis of LTEs is generally more complicated than that of short-term experiments and requires the proper use of multivariate statistical analyses procedures such as geostatistics [9], variance components, and repeated measures analyses [10, 11]. The serial correlation among measurements taken on the same plots in LTEs over many years can be accounted for by employing restricted maximum likelihood (REML) procedure [10–12]. In this case, repeated measures analysis provides models which are more parsimonious than those provided by other procedures such as Multivariate Analysis of Variance (MANOVA). Another advantage of repeated measures analysis is that it can accommodate several error terms and can handle incomplete designs and missing values [10].
Management factors, including crop rotations, phases within a crop rotation, and tillage and fertilizer applications, influence crop yields and their temporal variances; their effects are more pronounced under contrasting conventional and organic cropping systems [4]. Therefore, it is necessary to distinguish between the effects of management practices and those of spatio-temporal variation [2, 9]. Moreover, it is necessary to quantify the spatial variation in crop yield in cropping systems and how this variability is structured, how stable are the patterns of variability through time, and what are the underlying factors causing these patterns? Differences in management practices among cropping systems were found to contribute differently to the goal of achieving yield potential of conventional and organic cropping systems; additionally, cropping systems were found to influence the temporal variation and stability of crop yield [8, 13]. The objective of this study was to quantify covariance structures in conventional and organic cropping systems in a long-term experiment based on cumulative yield and its temporal variance and coefficient of variation in the presence of spatio-temporal variation.
2. Materials and Methods
2.1. Field Experiment
In a long-term split-plot experiment with four replicates in a randomized complete block design and comprising 16 treatments (Trt), that is, combinations of conventional and organic cropping systems, 2-year (corn-soybean) and 4-year (corn-soybean-wheat-alfalfa) crop rotations; all phases (Phase) of each crop rotation (i.e., corn-soybean, CS and soybean-corn, SC in 2-year crop rotations; alfalfa-corn-soybean-wheat, ACSW, corn-soybean-wheat-alfalfa, CSWA, soybean-wheat-alfalfa-corn, SWAC, and wheat-alfalfa-corn-soybean, WACS in 4-year crop rotations) were used in each of 8 years. In addition, conventional (CT) and strip tillage (ST); and fertility with nitrogen “+N” and without nitrogen “−N” fertilizer were used to estimate treatment effect on three secondary dependent variables derived from annual yield measurements. The dependent variables were cumulative yield, temporal yield variance, and CV% estimated for each factor during eight years in the presence of spatial and temporal variation. The first and second four years of the experiment comprised Cycle 1 and Cycle 2 of the crop rotations, respectively. Detailed information about the experiment is available elsewhere [13].
2.2. Data Collection
The procedures followed for data collection were reported earlier [8, 13] and will be presented briefly herein. Grain yield of corn, soybean, and wheat and dry matter produced by alfalfa were expressed in . Cumulative yield, yield variance, and coefficient of variation were estimated for each cropping system, crop rotation, and phase within crop rotation, over eight years, and used in subsequent geostatistical and multivariate statistical analyses.
2.3. Statistical Analyses
Descriptive statistics of cumulative yield, temporal yield variance, and CV% (i.e., least square means and minimum and maximum values) were calculated for each treatment “Trt” and phase within treatment “Phase(Trt).” These variables were used in (1) geostatistical analyses for the whole LTE and for each cycle separately and in developing two-dimensional map of spatial variation for each variable following interpolation using ordinary kriging [9, 14], (2) covariance data matrices based on cumulative yield, temporal variance, and coefficient of variation for each of two cycles in each cropping system which were used in extracting the first two factors in a principal components analysis and in testing the relationships between cropping systems, crop rotations, and phases within crop rotations among cropping systems. Matrix correlation statistics were calculated using Mantel test as performed by NTSYSpc [15] and verified by structural equation modelling [16, 17]. Statistical tests were then applied to determine whether the covariance structures implied by the structural equation models adequately fit the actual covariance structures of the data matrices based on a discrepancy function after performing 1000 permutations for each one of the matrix pairs; a nonsignificant () test result was considered as an adequate model fit, (3) variance components analysis using treatments as a fixed factor and Phase(Trt) as a random factor, where the variance components estimates associated with the random effects are independent of the fixed effects while at the same time taking into account their estimates, and (4) repeated measures multivariate analyses using first- (AR1) and second-order (AR2) autoregressive, first- (ANT1) and second-order (ANT2) antedependence, split-plot in time (SPT), and unstructured (UNS) covariance models under REML to approximate the relationship between observed values and errors associated with experimental units [10]. The choice of the most appropriate covariance structure for each variable was based on the largest Wald statistic value, and smallest deviance value, estimated as likelihood, after fitting with various within-subject covariance models [10, 18]. The Wald statistic has asymptotically a chi-squared distribution with one degree of freedom, provides an alternative test procedure to the likelihood ratio test, and, when divided by the degrees of freedom of its associated treatment sums of squares, is the F-statistic [10, 11]. Statistical analyses were performed using relevant modules in GenStat, version 10 [10], NTSYSpc [15] and STATISTICA, version 10. [17].
3. Results
Cumulative yield, considered as a function of large-scale deterministic structure and small-scale stochastic structure, varied among cropping systems, among treatments, and among phases of a crop rotation within treatments. Treatments within the conventional cropping system resulted in significantly larger cumulative yield and were associated with larger temporal yield variance and smaller CV% than those within the organic cropping system (Table 1). Minimum and maximum values for all variables, except minimum CV%, followed the same trend. The range (maximum-minimum) in treatment cumulative yield for the conventional cropping system (7.9) and the organic cropping system were almost similar in magnitude. However, the range in treatment temporal yield variance for the conventional cropping system (3.45) was much larger than the respective range for the organic cropping system (1.3), whereas the respective values for CV% were 3.1 and 17.3%. The ST-N and CT+N treatment combinations, regardless of crop rotation, produced minimum cumulative yield and temporal yield variance, respectively; the first was associated with the SC phase of 2-year crop rotation whereas, the second was independent of the 4-year crop rotation phases. The CT+N treatment combination, regardless of cropping system or crop rotation, produced maximum cumulative yield which was not always associated with maximum temporal yield variance or CV%; however, different treatment combinations among the conventional and the organic cropping systems resulted in maximum temporal yield variance and CV%. Notably, the SC and WACS crop rotation phases achieved maximum values in all variables.
3.1. Geostatistical Analyses
At the whole LTE level, a spherical model provided the best fit for cumulative yield with values increasing over time (from 0.35 to 0.90), whereas, exponential models provided the best fit for temporal yield variance ( from 0.25 to 0.57) and CV% ( from 0.23 to 0.49). This relatively improved model fit for cumulative yield was associated with increasing values of model parameters, except the proportion of sample variance that is explained by the spatially structured variance (i.e., ) which decreased over time from 0.89 to 0.64 (Table 2). The exponential model parameters for temporal yield variance and CV% displayed similar decrease (i.e., the nugget which is spatially unstructured variance and the structured variance), increase (i.e., range, m), and opposite trends over time. The exponential model fit ranged from 0.25 to 0.57 and from 0.23 to 0.49, for temporal yield variance and CV%, respectively.
Model parameters for cumulative yield in Cycle 1 and Cycle 2 were comparable; both temporal yield variance and CV% were fitted with exponential models, however, with smaller values for temporal yield variance in Cycle 1 (0.23) than in Cycle 2 (0.56). The respective values for CV% were 0.36 and 0.52 (Table 2). These changes in model fit were manifested over time by ~50% decrease in unstructured variance in both temporal yield variance and CV%, by 30 and 55% decrease in structured variance in temporal yield variance and CV%, respectively, and by 43% decrease in the proportion of sample variance that is explained by the spatially structured variance in temporal yield variance. These temporal changes were also accompanied by a 5.34-fold increase in the model range for temporal yield variance and 3.65-fold increase in the model range of CV%. Differences in most model parameters between both cycles, in addition to model fit (i.e., ), were larger in temporal yield variance and CV% than those observed in cumulative yield.
Basic statistics for all three variables and for each cycle are presented in Table 3, which provided quantitative information and assessment of cumulative yield, temporal yield variance, and CV%. The data indicated that there were major differences between both Cycles except for CV%, where data were generally comparable. The ±95% confidence intervals of the mean were numerically closer to the maximum rather than minimum values of each variable; however, there were a few exceptions, where the distribution is either positively or negatively skewed (data not presented). Basic statistics of cumulative yield suggested the presence of considerable variation (11.9 to 57.6 ) that may or may not have been associated with the same level of variation in temporal yield variance [2.0–13.5 ] or CV% (34.7%–111.4%). There were 22 statistically significant bivariate correlation coefficients out of the maximum 48 between cumulative yield, temporal yield variance, and CV% (Table 4). A minimum of two and a maximum of four significant correlation coefficients were found for each treatment across cropping systems and variables. Ten and 12 of these correlation coefficients were associated with the conventional cropping system and the organic cropping system, respectively; 12 and 10 with CT and ST, respectively, 12 and 10 with +N and –N, respectively, and 9 and 13 with 2-year and 4-year rotations, respectively. Eight positive correlation coefficients were found between cumulative yield and temporal yield variance mostly in 4-year rotations in both cropping systems regardless of tillage or fertilizer treatments. There were fewer significant correlations between cumulative yield and its CV%, three of which were negative and associated with the 2-year rotations, whereas the largest number of significant correlation coefficients was between temporal yield variance and CV% across several treatment combinations and in both cropping systems.
3.2. Covariance Structures
When the whole LTE was considered, the first and second factors in the principal components analyses accounted for 62.9% and 36.3% of total variation (Table 5). The respective values for the conventional system were 60.1% and 39.4%, and for the organic system were 52.2 and 47.1%. Factor 1 in the whole LTE and in both cropping system was dominated by positive loadings of all variables except CV1. Factor 2 displayed an almost mirror image of Factor 1 with positive loadings of TYV and CV% for the LTE while in the conventional system it had negative loadings of TYV and CV%, and in the organic system it had positive loadings of all variables. The supplementary variables in LTE, conventional and organic systems, followed the same trend in loading on Factor 1 and Factor 2, but with relatively smaller values. The communality values, which is the proportion of variance that each variable has in common with other variables, were large for TY, TYV, and CV, on the basis of their statistical relationships, whereas, they were smaller and comparable in values among Cycle 1 and Cycle 2 for each supplementary variable, however, with a few exceptions. These include a large difference between communalities of TYV in Cycle 1 (0.68) and Cycle 2 (0.37), and for CV% in Cycle 1 where it had a larger common variance in LTY (0.79) and conventional (0.83) and organic systems (0.77) than the respective values of 0.37, 0.04, and 0.20 in Cycle 2.
Conventional and organic cropping systems pairwise dissimilarity (i.e., covariance) matrix correlation statistics (Table 6) indicated that matrix sums of squares differed substantially in magnitude between the conventional and organic cropping systems for most pairwise comparisons. The -values ranged from negative (−0.31 for WACS phase) to positive (0.96 for CSWA phase), with ~50% of these values larger than 0.50 and 30% of which were found between tillage-fertility treatments of conventional and organic cropping systems. The probability of the normalized Mantel statistic being equal to, or smaller than the observed , ranged from 0.12 to 0.99, whereas the value of the same statistic based on 1000 permutations ranged from 0.175 to 0.999. According on these statistics, about 50% of the pairwise matrix comparisons were highly significant suggesting that substantial differences existed between their pairwise covariance structures.
3.3. Variance Components
Results of the REML-based mixed model analyses indicated that differences between treatments in cumulative yield, temporal yield variance, and CV%, as a fixed factor, were highly significant (Table 7). The -values suggest that variation between treatments in cumulative yield was much larger than the variation within treatments (i.e., error variance) as compared to those for temporal yield variance and CV%. Also, differences between phases of a crop rotation within treatments, as a random factor, were highly significant and accounted for 59.8%, 50.1% and 79.6% of total variances in cumulative yield, temporal yield variance, and CV%, respectively.
3.4. Repeated Measures Analyses
A summary of the extensive repeated measures analyses, which accounted for serial correlations among observations taken on the same plot over time, is presented in Table 8. These models were the best covariance models that fit the data; they were selected on the basis of their smallest deviance and largest Wald statistic. Three models (ANT1, AR1, and UNS, in increasing order of importance) were the best models selected from the initial 72 output models based on combinations of 2 systems × 2 factors × 3 variables × 6 covariance models. The ANT1 model was the best to fit cumulative yield for the Time × Phase(Trt) in the conventional cropping system, the AR1 model was the best to fit cumulative yield for Time × Trt in both the conventional cropping system and the organic cropping system, and the UNS model was the best to fit the remaining eight factor-variable combinations.
4. Discussion
The diversity-productivity-stability relationship is gaining importance in agricultural research and crop production [1]. Often, the lack of temporal stability of the spatial structure in field experiments may indicate that the factors controlling crop yield are dynamic [19]. Therefore, an LTE is necessary to allow both mechanisms and temporal dynamics to be identified and quantified [4]. The current analysis, which was focused at the time × treatment interaction, highlighted the importance of indirect effects not apparent in the short term, as well as responses in cumulative yield and its temporal variance and coefficient of variation that may change over the long term. The least squares means and associated statistics (Table 1) are reliable indicators of temporal variance as influenced by fixed and random factors in the experiment. In addition, the results point to the consistency of performance of the organic cropping system in relation to the conventional system whether at the minimum (67%), mean (68%), or maximum (67%) yield levels. Differences between cropping systems and between treatments are in agreement with earlier findings [8], whereas differences between phases within treatments emphasize the significance of multiple crop effects on all variables in this study. Similar conclusions can be derived from comparative assessments of temporal yield variance and CV%.
4.1. Geostatistical Analyses
The comprehensive geostatistical and covariance analyses of the current LTE prevented random factors (Table 2) from interfering with the desired inferences on the fixed factors [5]. Experimental error variances were reduced substantially due to the presence of all phases of each crop rotation every year over the duration of the experiment [8, 12], and all crops within each crop rotation were exposed to the same environmental fluctuations. The assumption that not only yield (and cumulative yield) but also its temporal variance and coefficient of variation are considered random variables that are functions of large-scale deterministic structure and small-scale stochastic structure [7] was validated in this study (Table 2). Interactions of time × treatments varied with the age of the LTE [5] as deduced from geostatistical model parameters for the whole LTE and for both cycles (Table 2) and from the magnitude of correlation coefficients between the variables under study (Table 4). Total variance of crop response in LTEs, with crop rotations in their design, can be separated into spatial (representing between-plots) and temporal (representing within-plots) variances [7]. Spatial variance provided information on the main effects and interactions of experimental factors, and temporal variance provided information on changes in the magnitude of the effect of their interactions with time. Obviously, there was differential build up of the cumulative effects as the crop rotations were repeated over time (Cycle 1 versus Cycle 2; Table 2) as quantified by the decreasing structured and unstructured variances and increasing range estimates in the geostatistical analyses, especially for temporal yield variance and CV% [7, 20].
Patterns of nonsignificant correlation coefficients between pairs of cumulative yield, temporal yield variance, and CV% (~50% of all pairwise correlations; Table 4) as compared with loadings and amount of variation explained in the principal components analyses (Table 5) suggest that treatments within the 2-year crop rotation may have shaped the relationships between cumulative yield (TY) and temporal yield variance (TYV) in both cropping systems. A similar conclusion can be based on the non-significant pairwise correlation coefficients between cumulative yield and CV% in both cropping systems, with treatments within 2-year crop rotations having more influence on shaping the relationship between TY and CV% in the organic cropping system; whereas the same effect can be attributed to treatments within 4-year crop rotation in the conventional cropping system.
4.2. Covariance Structures
Covariance structures based on principal components (Factors in Table 5) and matrix correlation (Table 4) analyses are indicators of the power of each variable (i.e., cumulative yield, temporal yield variance, and CV% in both cycles) in shaping the covariance matrix and in determining the magnitude of its deviation from its pair within cropping systems or within crop rotations. Based on PCA of these variable (and their association with supplementary variables) on Factor 1 (which explained more variance than Factor 2 in both conventional and organic systems, Table 5), it can be concluded that these cropping systems differed substantially in the structure of their covariances (Table 6) and consequently in the structure of their principal components. The covariance matrices of these cropping systems, quantified by their respective sums of squares (Table 6), were independent of each other (, ) and the difference between their covariances was associated positively (, ) and negatively (, ) with the covariance matrix of conventional and organic cropping systems, respectively. Additionally, the level of significance, with and without permutations, was significantly correlated with covariance matrix of the organic system (, ; , , resp.) but not with the conventional cropping system. Organic management practices over the duration of the experiment may have altered the magnitude and direction of relationships between measured variables and resulted in significant deviations from the conventional cropping system.
4.3. Variance Components Analyses
The estimates of variance components (Table 7) are in general agreement with the large differences between minimum and maximum values derived from the least mean squares analyses (Table 1), while the variation between phases within treatments accounted for >50% of variation in cumulative yield, temporal yield variance, and CV%. Estimates of variance components can be used to in the design of future long-term experiments [5, 11]; they provide the basis for estimating the optimum number of replications and samples and may identify which variables to monitor and statistically analyze. Based on values (Table 7), it is obvious that all three variables (CV%, cumulative yield, and temporal yield variance) contributed, in decreasing order, to the precision with which the analyses were performed [13].
4.4. Repeated Measures Analyses
Three covariance models were selected, on the basis of the deviance and Wald statistics [9] to fit factor-treatment combinations in the experiment (Table 8). The ANT1 covariance model is considered as a generalization of the multivariate analysis of variance that accommodates the types of covariance structures that typically occur with repeated measurements. According to ANT1 characteristics, variances in cumulative yield for the Time × Phase(Trt) in the conventional cropping system at different times are assumed unequal [11, 18]. The AR1 model was the best fit for cumulative yield based on Time × Trt in the conventional cropping system and the organic cropping system; therefore, it is assumed that all variances are equal, and presumably, this is due to the absence of Phase(Trt) effects. Finally, the UNS covariance model was the most flexible and was the best to fit temporal yield variance and CV% of most factor-variable combinations, but it was the most costly in terms of the number of estimated parameters [11, 12, 21]. Therefore, it is assumed that variances and covariances are considered arbitrary; the model allows for unequal variances and unequal covariances among the data. The UNS model was the best fit for temporal yield variance and CV% of both cropping system and of the cumulative yield for Time × Phase(Trt) in the organic cropping system.
5. Conclusions
A long-term experiment provided an opportunity to explore the effect of 16 combinations of management practices on the consistency of the spatial and temporal structures of cumulative crop yield and its temporal variance and coefficient of variation in conventional, and organic cropping systems. Geostatistical, covariance and variance components analyses resulted in mapping spatial variation, partitioning total variation into its components, and in selecting appropriate covariance models to describe the correlation between random effects in the experiment using restricted maximum likelihood in conjunction with repeated measures. Differential build up of cumulative treatment effects was observed as the crop rotations were repeated over time. The unstructured covariance model, with the smallest deviance and largest Wald statistic, was the best to fit most factor-variable combinations; it was the most flexible, but most costly in terms of computation time and the number of estimated parameters.
Abbreviations
| 2-year: | Two-year crop rotation |
| 4-year: | Four-year crop rotation |
| ACSW: | Alfalfa-corn-soybean-wheat |
| ANT1 and 2: | First- and second-order antedependence covariance model, respectively |
| AR1 and 2: | First - and second-order autoregressive covariance model, respectively |
| CNV: | Conventional cropping system |
| CS: | Corn-soybean |
| CSWA: | Corn-soybean-wheat-alfalfa |
| CT: | Conventional tillage |
| CV%: | Coefficient of variation |
| CumY: | Cumulative yield over time |
| Cycle 1, 2: | First- and second-cycle of 2-year or 4-year rotation, respectively |
| DMRT: | Duncan multiple range test |
| LTE: | Long-term experiment |
| MANOVA: | Multivariate analysis of variance |
| MSS: | Matrix sums of squares |
| +N: | With nitrogen fertilizer |
| −N: | No nitrogen fertilizer |
| ORG: | Organic cropping system |
| Phase(Trt): | Phase of a crop rotation within a treatment combination |
| REML: | Restricted error maximum likelihood |
| RM: | Repeated measurements |
| SC: | Soy-corn |
| SPT: | Split-plot in-time covariance model |
| ST: | Strip tillage |
| SWAC: | Soybean-wheat-alfalfa-corn |
| Trt.: | Treatment |
| TYV.: | Temporal yield variance |
| UNS: | Unstructured covariance model |
| WACS: | Wheat-alfalfa-corn-soybean. |
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Thanks are due to Jim Eklund, Steve Wagner, Charles Hennen, and Scott Larson for their assistance in the field work. This research was funded by USDA-ARS Project no. 3645-61600-001-00D, Morris, MN, USA. The use of trade, firm, or corporation names in this publication is for information and convenience of the reader. Such use does not constitute an official endorsement or approval by the United States Department of Agriculture or the Agricultural Research Service of any product or service to the exclusion of others that may be suitable. USDA is an equal provider and employer.