Soil Science Society of America Journal 66:1134-1142 (2002)
© 2002 Soil Science Society of America
DIVISION S-1—SOIL PHYSICS
Accounting for Soil Spatial Autocorrelation in the Design of Experimental Trials
M. Fagrouda and
M. Van Meirvenne*,b
a Ecole Nationale d'Agriculture de Meknès, BP S40, 50000 Meknès, Morocco
b Dep. of Soil Management and Soil Care, Ghent Univ., Coupure 653, 9000 Gent, Belgium
* Corresponding author (marc.vanmeirvenne{at}rug.ac.be)
 |
ABSTRACT
|
|---|
Soil heterogeneity complicates the design and analysis of field experiments. Block designs were developed for this purpose. However, the analysis of experimental results supposes that the residuals from the treatment are spatially independent and that within block variation is random. Experience indicates that this rarely is the case in field experiments, because of the strong spatial autocorrelation of soil properties. This paper applies geostatistical tools, such as variogram analysis and conditional stochastic simulation, to investigate the optimal experimental plot size and shape and to decide which experimental design is to be preferred. The methodology is illustrated using a case study of water-use efficiency under semiarid conditions in Morocco. It was found that under these conditions an experiment with 16 treatments would use best a plot size of 4 by 8 m oriented north-south, configured according to an incomplete block design with 8 plots per block oriented in two rows in the east-west direction.
Abbreviations: AWC, available-water capacity • ccdf, conditional cumulative distribution function • CV, coefficient of variation • Db, bulk density • Ei, average efficiency • FCw, field capacity • GM, gravimetric moisture • NSR, nugget/sill ratio • OC, organic C • PWPw, permanent-wilting point • SGSIM, Sequential Gaussian Simulation (SGSIM)
|
INTRODUCTION
|
|---|
SOIL SPATIAL VARIATION is an important source of external variation that affects the design of field experiments on, for example, water-use efficiency. Consequently, intensive research (e.g., Cochran and Cox, 1957) has been conducted to develop improved methods for determining treatment effects when soil variability affects plot yield.
Traditionally, well-designed experiments are based on the concepts of replication, blocking, and randomization. Replication implies that a given number of treatments are applied under identical conditions to different experimental plots. The aim of using replications is to estimate the experimental error. However, since in water-use efficiency trials the underlying experimental material is soil, it will rarely be homogenous. Therefore experimental plots are combined into blocks and blocks are replicated. To overcome the effect of spatial dependence of soil properties, plots are randomized within blocks. Randomization of plots within blocks is expected to equalize error over all treatment differences thus allowing the use of statistical methods of analysis, such as analysis of variance and regression analysis. However, the underlying hypothesis is that spatial variability is random; i.e., spatially unstructured, within blocks. However, in the presence of a significant spatial correlation over small distances, the assumption of independence between plots is violated. In such a situation, a field researcher may be faced with contradictory results: clear difference in crop yields between experimental plots but no significant treatment effect (see e.g., Bhatti et al., 1991; Van Es and Van Es, 1992).
Consequently, there has been much interest in methods to analyze experiments that take into account spatial correlation (e.g., Wilkinson et al., 1983; Zimmerman and Harville, 1991; Zhang et al., 1994; Brownie and Gumperz, 1997). These methods adjust the observations by values measured on neighboring plots. Smith (1938) was the first to model empirically variations in soil fertility to compare experimental designs. A simple method for comparing different experimental designs was by placing different experimental plans in a uniformity trial and then calculating the residual variation for each of these plans (Claustriaux and Rousseau, 1974). Patterson and Hunter (1983) used two different models to examine the effect of block sizes in variety trials with incomplete blocks but they did not take plot size or shape into account.
An alternative approach to the quantification of heterogeneity is by applying geostatistics which is used mostly to analyze the structure of the spatial variability of regionalized variables such as soil properties and to optimize their interpolation. Only a few studies (including Ersboll, 1996; Van Es and Van Es, 1993) have related geostatistics to the design of field experiments.
The practical purpose of this paper is to present a methodology to account for soil spatial autocorrelation in the design of experiments on water-use efficiency under the semiarid conditions of Morocco. Conditional stochastic simulations was used to generate a number of equiprobable realizations of the most important soil property, available water capacity (AWC). These realizations were then used to analyze the experimental setup in terms of (i) plot size and plot shape, and (ii) the design to be preferred.
|
MATERIALS AND METHODS
|
|---|
Study Area and Soil Sampling
The study site is located within the experimental farm of the National School of Agriculture in Meknes, Morocco, (central coordinates: 33°52' N long., 5°33' W lat., altitude: 625 m above sea level) situated on the Sais plateau (Fig. 1)
. It has been used previously to conduct experiments on yield response (Corbeels et al., 1998). This area has a semiarid climate with a temperate winter. Rains occur mainly between October and May, with an average monthly maximum of 80 mm in December. The average annual rainfall is 526 mm, but extreme fluctuations occur between years.
In August 1995, we sampled the topsoil (0–40 cm) at 100 locations within a 150 by 150 m field according to the pair random sampling method (Van Meirvenne and Hofman, 1991). Therefore, fifty samples were located randomly within the field and 50 additional samples were each located randomly within a circle of 5 m around one of the 50 previously located points. This configuration allows assessing the spatial variability at both short and longer distances in respect to the dimension of the experimental field. The radius of 5 m was chosen intuitively to ensure a sufficient number of samples intervals within this distance.
Disturbed soil samples were analyzed for the three main particle-size fractions (clay, silt, and sand) by the usual pipette-sieve method, CaCO3, organic C (OC), and the gravimetric moisture (GM) content at -0.033 and -1.5 MPa using a pressure membrane apparatus. These two GM contents were taken to represent field capacity (FCw) and the permanent-wilting point (PWPw) on a weight basis, respectively. The difference between FCw and PWPw was multiplied with the soil bulk density (Db) to obtain the volumetric AWC (Cassel and Nielsen, 1986). The latter was further multiplied by 1000 to yield the AWC within 1-m soil depth, expressed as milimeters, assuming that Db remained unchanged within the soil profile.
Since the soil was too hard and too dry to take undisturbed soil samples, Db could not be observed directly. Therefore we used published Db and soil data from within the studied region to build a pedotransfer function. It was found that Db could best be predicted from the percentages of silt (S), CaCO3 (C), and OC by:
 | [1] |
Variations in the AWC were found to be responsible for yield variations in the study area (Fagroud et al., 1997); therefore, we used it as the target variable to analyze the design of experimental plots.
Geostatistics
The sample variogram 
(h) is commonly used to estimate the structure of the spatial variation of a regionalized variable Z. It can be obtained from:
 | [2] |
where z(xi) and z(xi+h) are the observed values of Z at locations xi and xi+h respectively, h is a separation vector and n(h) is the number of paired comparisons at that lag (Journel and Huijbergts, 1978).
A theoretical model is fitted to the experimental values of the variogram, usually using three parameters: the range, the sill, and the nugget effect. The range is the separation distance beyond which two observations are independent of each other. The sill is the variogram value corresponding to the range. The discontinuity at the origin is called the nugget effect and arises from a combination of random errors and sources of variation at distances smaller than the shortest sampling interval (Goovaerts, 1998).
Based on the variogram model, conditional stochastic simulation was used to generate a number of sets of values, called realizations, which aim to reproduce each the sample histogram and the variogram model. Such a set of L realizations can be synthesized according to the objectives of the study. Several conditional stochastic simulation algorithms are available, (Deutsch and Journel, 1998) and we used Sequential Gaussian Simulation (SGSIM) that proceeds as follows (Fig. 2)
:
- Data are transformed according to y(xi) =
(z(xi)), with z(xi) the original data, (·) a transformation function, and y(xi) the normal scores having a standard normal (or Gaussian) histogram (Goovaerts, 1997).
- The sample variogram of the normal scores is computed and modeled.
- The L simulations are performed in normal score space as follows:
- define a random path visiting each of the m unsampled locations to be simulated only once;
- at each unsampled location x0, estimate the parameters (mean and variance) of the Gaussian conditional cumulative distribution function (ccdf) by simple kriging using the normal score variogram model and the mean value of the normal scores. The conditioning information consists of n neighboring data of both original normal score data y(xi) and values y(l)(x0) simulated at previously estimated locations, l being the realization number (l = 1,..,L);
- randomly draw a simulated value y(l)(x0) from the ccdf, and add it to the data set;
- proceed to the next location along the random path, and repeat the two previous steps;
- loop until all m locations are simulated;
- proceed with the next simulation by repeating the previous steps, until all L realizations are available.
- The results are finally back-transformed to the original variable space by applying the inverse of the normal score transform (·) to the simulated y-values.
Each of these L simulations is a realization of the unknown spatial distribution of Z. Differences between them provide a measure of the spatial uncertainty about Z which was used to optimize the design of our experiment.
|
RESULTS
|
|---|
Exploratory Data Analysis
Some descriptive statistics of the variables of interest are summarized in Table 1. As can be observed, the study site mainly consists of calcareous (on average 6.5% CaCO3) soils with a clayey texture and a low OC content (on average 1.5%). Moderate variations (coefficient of variations [CVs] between 9 and 30%) were observed for all these variables, with the exception of CaCO3 (CV of 110%).
View this table:
[in this window]
[in a new window]
|
Table 1. Some statistical parameters of the studied variables of 100 topsoil samples (CV = coefficient of variation).
|
|
Because of the clayey texture, the average PWPw and FCw values differ only by 10%. Consequently the average AWC within 1 m is low (129 mm), it exceeds 150 mm only at 10% of the locations and has a minimum of only 70 mm. Because of the extreme dryness in summer, soils dry out well beyond PWPw and thus require quite some precipitation (or irrigation) amounts before they can support plant growth. But too large water additions quickly result in drainage losses because of the small AWC of the top 1 m. Although both PWPw and FCw display moderate variations (CVs of 6.9 and 5.5%, respectively), their combination into AWC within 1 m results in a CV similar to soil texture (15%). Hence the importance of an accurate mapping of water holding characteristics to improve the management efficiency of these soils (El Jaafari et al., 1993; Corbeels et al., 1998).
Variogram and Simulations
A normal score transformation was applied to the AWC data and the omni-directional experimental variogram of normal scores exhibits a clear spatial structure (Fig. 3)
. This variogram was modeled by a spherical model (McBratney and Webster, 1986) with a range of 80 m and a nugget effect which represents 
42% of the total variance indicating an important proportion of microscale and random variation. Since AWC was calculated from the combination of several variables (including a pedotransfer function), this relatively large nugget effect could be explained partly by the accumulated uncertainty related to the analytical determinations and the pedotransfer function.
Five hundred realizations of the spatial distribution of AWC were generated over a 100 by 100 m subpart of the study area using a grid spacing of 1 by 1 m. A map of the mean values of these 500 realizations, which represents the estimates of the expected value at every location, is shown in Fig. 4
(top left). Available water capacity tends to decrease towards the lower part of the area (small y-coordinates), which could indicate some overall nonstationarity This was, however, not taken into account since we used a small neighborhood, assuming local stationarity.
View larger version (77K):
[in this window]
[in a new window]
|
Fig. 4. (a) Mean of 500 realizations of available water capacity (AWC), (b) means of aggregating AWC in plots of 4 x 8 m, (c) residual mean map from incomplete blocks of eight treatments and (d) corresponding residual variogram.
|
|
Plot Size and Shape
Plot Configurations
We considered 24 different combinations of plot sizes and plot shapes with properties presented in Table 2. Plot sizes were chosen to avoid very small plots (<8 m2) or large ratios of length/width since these were considered unrealistic in practice. Figure 5
illustrates some of these configurations for a 20 by 20 m subpart of the study area. Each of the 500 realizations created by SGSIM at a grid spacing of 1 by 1 m were aggregated into these 24 different plot configurations. Goovaerts (1999) provides some arguments why aggregating simulated values to block estimates is to be preferred over block kriging.
View this table:
[in this window]
[in a new window]
|
Table 2. Plot sizes and plot shapes considered and number of plots which could be created within the 100 by 100 m study area (see also Fig. 5).
|
|
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 5. Some of the 24 plot configurations considered (see Table 2) for a 20 by 20 m subpart of the simulated areas (the top plot represents the original grid spacing (1 by 1 m) of the simulated maps before aggregation).
|
|
Evaluation Criterion
To evaluate the efficiency of the aggregation into plots, the average AWC of every plot was computed for each of the 500 realizations and for each of the 24 plot configurations of Table 2 using:
 | [3] |
where p and q are the width and length of the plots considered, respectively (i.e., p = 2, 4, 5, 8, and 10 m; q = 2, 4, 5, 8, and 10 m), i and j are column and row number of the plot configuration and k, l represent the column and row number of the simulated map. To compare the degree of efficiency, the average AWC values were used to compute experimental, and model theoretical, variograms of every plot configuration for a given realization. The nugget/sill ratio (NSR) was used as a criterion indicating the extent to which the experimental errors between plots were randomly distributed in space (Bhatti et al., 1991; Ersboll, 1996). This calculation was repeated for each of the 500 realizations, yielding an empirical distribution of NSR values for every plot configuration. Therefore the average NSR and its confidence intervals were considered. The smaller the average NSR, the less random errors remain between plots; thus, the less the plot configuration meets the condition that experimental errors can be considered to be random. So the largest NSR represents the configuration which best meets the underlying hypothesis of classical ANOVA techniques.
Determination of the Most Efficient Plot Size and Plot Shape
Figure 6
shows the average NSR versus plot width, pooled for all lengths. The maximum NSR was reached for plot widths of 4 m, but the 95% confidence intervals were large (from 36.5 to 55.1%) and they overlapped with the intervals of 5 and 8 m. The smallest NSR was found for a width of 10 m with narrow intervals (between 29.2 and 30.3%).
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 6. Average Nugget/Sill ratio as a function of plot width for AWC within 1 m. Bars indicate the empirical confidence intervals with a 95% confidence level.
|
|
Figure 7
illustrates, for two fixed widths (4 and 10 m), the effect of plot shape on the NSR. It appeared that for a width of 10 m, the NSR showed only small fluctuations between the considered plot lengths. In contrast, a large variation of NSR was observed for a plot width of 4 m with a clear maximum for a plot length of 8 m. Therefore, it was concluded that plots of 4 x 8 m could be considered to be most efficient in terms of obtaining the largest spatially unstructured between plot variance. The corresponding map of the mean AWC per plot of 4 by 8 m is given in Fig. 4 (top right) where it can be observed that the orientation of the overall variability of AWC is more or less parallel to the longest direction of the plots.
Experimental Design
Complete versus Incomplete Block Designs
Planning a new experiment involves deciding on the experimental design. Two elements are to be considered in the design of an experiment: (i) the nature of the blocks (complete or incomplete) and (ii) the number of plots per block and their configuration within the blocks (orientation and shape). Blocks of experimental plots are called complete when the number of plots is equal to the number of treatments (e.g., different levels of irrigation amounts). Since this involves extensive layouts when the number of treatments to be investigated is large and when blocks are being replicated several times, incomplete block designs were created. In incomplete blocks, the number of experimental plots per block is smaller than the number of treatments, but each treatment occurs equally frequently over all blocks (called "balanced incomplete blocks arrangements" by Cochran and Cox, 1957).
To compare various designs, we examined the following two statistics:
- The parameters of the variogram model of the plot residuals as summarized in the NSR. Since the ranges of the variogram models greatly differ, this ratio was standardized as NSR/Range.
- The average efficiency, Ei, of a given incomplete block design i compared with the corresponding complete block design with the same configuration:
 | [4] |
with RMSi as the residual mean square of the ith incomplete block design and RMSc the residual mean square for a complete block design.
Block Configuration and Comparison Criterions
We examined the layout of 7 experimental designs to be applied in a field experiment with 16 treatments, as shown in Table 3. The designs with incomplete blocks and the lattice square were constructed using plans taken from Cochran and Cox (1957). For each design, various combinations of plot configurations within blocks were considered as shown in Fig. 8
for eight treatments (as in the incomplete block designs). Blocks were oriented either in the east-west (horizontal) or in the north–south direction. For the designs with 8, 10, 12, and 16 plots per block all four configurations were compared; whereas for the designs with 4, 5, and 6 plots per blocks only two configurations (I and III) were considered because Configurations II and IV (with double lines or double columns) would become too small for practical implementations. An overview of the various configurations is presented in Table 4 together with the block sizes, given that the dimensions of every experimental plot was 4 by 8 m as determined in the previous chapter. In total 20 configurations belonging to seven different designs were evaluated. It is important to note that the relative superiority of Configurations III and IV compared with I and II is because of the overall spatial tendency in AWC (see Fig. 4 top left). Thus, blocks directed east-west would contain less spatial variability, so one could expect a larger part of the total variation to be attributed to differences between blocks.
View larger version (10K):
[in this window]
[in a new window]
|
Fig. 8. Configurations for one block with eight treatment plots as in an incomplete block design. Each block consist of one column (type I), one row (type III), two columns (Type II), or two rows (Type IV).
|
|
Each of the considered designs was overlaid on each of the 500 simulated maps of AWC. Each time, 50 randomizations of treatment numbers within the blocks were generated using a pseudo random number generator. So for each design of Table 4, 25000 different configurations were investigated by an analysis of variance and the residuals were computed by subtracting block effects from plot values (we eliminated only the block effects since no treatment effect was applied). Again, we calculated and modeled the variograms of the residuals for each of the 25000 configurations and retained their parameters using Variowin (Pannatier, 1996). As before, this allowed us to calculate the average, and its confidence intervals, of the two evaluation criteria.
The ratio NSR/Range indicates (Fig. 9)
that the optimal number of plots per block should be neither too high (NSR small) nor too small (significant spatial correlation) since it represents a compromise between two parameter requirements indicating a weak spatial dependence of the residual errors: a large NSR and a short range. The best result was obtained for an incomplete block design with eight plots per block, yielding an average NSR to the range ratio of 2.35, and the more one moves from this optimal number of plots, the more this statistics decreases, reaching on average 0.13 and 0.26 for 4 and 16 plots respectively.
The average Ei (Eq. [4]) for the incomplete block designs with eight plots compared with a complete block design of the simulations are shown in Fig. 10
. All incomplete block designs resulted in an average Ei of more than 100, suggesting their superiority over a complete block design. Although no significant differences were found between the different incomplete block designs, it appears that Design 9 (Configuration III) gives better results than the other configurations.
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 10. Average relative efficiency (Ei) of the incomplete block designs compared with the corresponding complete block design for a design with eight parcels per block. Design number refers to Table 4. Bars indicate the empirical 95% confidence intervals.
|
|
Therefore, we concluded that the best experimental design to investigate the effect of the AWC in the top 1 m of our experimental field consists of an incomplete block design with eight plots per block configured according to Type III (Fig. 8). For illustrative purposes, the corresponding residual map is presented on Fig. 4c. An almost complete random pattern was obtained. To quantify this impression, the experimental variogram of the residuals was calculated (Fig. 4d). Its model approached a pure nugget effect (85% NSR) indicating that the residuals of Fig. 4c behave almost purely randomly in space. So the hypothesis underlying an ANOVA analysis was approached closely.
|
CONCLUSIONS
|
|---|
The design of experimental trials involves two important decisions: the dimension and shape of experimental plots and the design of plots into blocks.
To investigate the most suitable plot dimension, SGSIM was used to generate 500 realizations of AWC to which were applied 24 different plot configurations. For each of these combinations plot average values of AWC the NSR of the experimental variogram was used as a statistic to select the most optimal plot size, which was in our case study 4 by 8 m, oriented north-south (parallel to the orientation of major variation of AWC over the study area). This plot dimension resulted in the largest average NSR, 52.7%.
To identify the preferred experimental design, the 500 simulated maps of AWC were overlaid by 20 configurations belonging to seven different designs, considering 16 treatments. Again variogram analysis was used to identify the design which accounted for most of the spatial autocorrelation in AWC. It was found that the experimental design to adopt, for a 16 treatments trial in the studied field, should be an incomplete block design with eight plots per block.
Because of the overall tendency of AWC, showing a stronger variation in the north-south direction, it was found that the optimal configuration was the one with two lines of plots per block oriented in the east-west direction (Type IV).
|
ACKNOWLEDGMENTS
|
|---|
The first author thanks the Belgian Development Cooperation (DGIS) for partly funding this research through a grant.
Received for publication July 25, 2001.
|
REFERENCES
|
|---|
- Bhatti, A.U., D.J. Mulla, F.E. Koehler, and A.H. Gurmani. 1991. Identifying and removing spatial correlation from yield experiments. Soil Sci. Soc. Am. J. 55:1523–1528.[Abstract/Free Full Text]
- Brownie, C., and M.L. Gumperz. 1997. Validity of spatial analyses for large field trials. J. Agric. Biol. Environ. Stat. 2:1–23.
- Cassel, D.K., and D.R. Nielsen. 1986. Field capacity and available water capacity. p. 901–926. In A. Klute (ed.) Methods of soil analysis: Part 1. 2nd ed. Agron. Monogr. 9. ASA and SSSA, Madison, WI.
- Claustriaux, J.J., and G. Rousseau. 1974. Analyse statistique des résultats d'essais à blanc. (In French.) Biometrie Praximetrie. 14:57–79.
- Cochran, W.G., and C. Cox. 1957. Experimental designs, 2nd ed. John Wiley & Sons, New York.
- Corbeels, M., G. Hofman, and O.Van Cleemput. 1998. Analysis of soil water use by wheat grown on a cracking clay soil in a semi arid Mediterranean environment: Weather and nitrogen effects. Agric. Water Manage. 38:147–167.
- Deutsch, C.V., and A.G. Journel. 1998. GSLIB: Geostatistical software library and user's guide, 2nd ed. Oxford University Press, New York.
- El Jaafari, S., R. Paul, Ph. Lepoivre, J. Semal, and E. Laitat. 1993. Résistance à la sécheresse et réponse à l'acide abscissique: analyse d'une approche synthétique. (In French.) Agricultures. 2:256–183.
- Ersboll, A.K. 1996. Spatial experimental design. p. 112–125. In J.M.C. Ocerin et al. (ed.) III HaRMA meeting; Cordoba: Design of experiments and statistical education in agriculture. European Community, Brussels.
- Fagroud, M., M. Van Meirvenne, and M. Corbeels. 1997. Variabilité spatiale de la texture et conséquence pour la récolte de tournesol dans un champ à Meknès, Maroc. (In French.) Pédologie-Themata. 2:61–63.
- Goovaerts, P. 1997. Geostatistics for natural resources evaluation. Oxford University Press, New York.
- Goovaerts, P. 1998. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Biol. Fertil. Soils. 27:615–334.
- Goovaerts, P. 1999. Geostatistical tools for deriving block-averaged values of environmental attributes. J. Geographical Information Sci. 5:88–96.
- Journel, A.G., and C.J. Huijbergts. 1978. Mining geostatistics. Academic Press, New York.
- McBratney, A.B., and R. Webster. 1986. Choosing functions for semi variograms of soil properties and fitting them to sampling estimates. J. Soil Sci. 37:617–639.
- Pannatier, Y. 1996. VARIOWIN: Software for spatial data analysis in 2D. Springer-Verlag, New York.
- Patterson, H.D., and E.A. Hunter. 1983. The efficiency of incomplete block designs in National list and recommended list cereal variety trials. J. Agric. Sci. 101:427–433.
- Smith, F.H. 1938. An empirical law of describing heterogeneity in the yield of agricultural crops. J. Agric. Sci. 28:1–23.
- Van Es, H.M., and C.L. Van Es. 1992. Soil spatial correlation and its implications for field experimental design. p. 275–286. In Proceeding of Pedometrics-92: "Developments in spatial statistics for soil sciences". International Agricultural Centre, Wageningen, the Netherlands.
- Van Es, H.M., and C.L. Van Es. 1993. Spatial nature of randomisation and its effect on the outcome of field experiments. Agron. J. 85:420–428.[Abstract/Free Full Text]
- Van Meirvenne, M., and G. Hofman. 1991. Sampling strategy for quantitative soil mapping. Pedologie. 41:263–275.
- Wilkinson, G.N., S.R. Eckert, T.W. Hancock, and O. Mayo. 1983. Nearest neighbour analysis of field experiments. J. Royal Statistical Society, Ser. B. 45:151–211.
- Zimmerman, D.L., and D.A. Harville. 1991. A random field approch to the analysis of field-plot experiments and other spatial experiments. Biometrics 47:223–239.
- Zhang, R., A.W. Warrick, and D.E. Myers. 1994. Soil heterogeneity and plot shape effect. Geoderma 62:183–197.
TOP
ABSTRACT
INTRODUCTION
