Exploratory Model Selection for Spatially Designed Experiments – Some Examples

Exploratory model selection was used to find a response model that accounted for the spatial variability present in the experimental results from four examples of spatially designed field experiments. It was found that the class of differential gradients within incomplete blocks was useful for finding a response model that accounted for the spatial variability present in the first example. The class of orthogonal polynomial regressions of response on row and column position and interactions of the regressions was useful for discovering an appropriate response model for the data of examples two, three, and four. The results obtained from the selected response model were compared with standard textbook analyses. Considerable differences in residual mean squares, coefficients of variation, and F-values for treatment to residual mean squares were found. The increase in replication for the selected response model over the textbook response model is demonstrated. The increase can be many fold.


Introduction
Data from four different examples involving field experiments were examined to determine an appropriate response model that accounted for the variability present in the experimental results.An experimenter selects an experiment design plan that is thought appropriate for a forthcoming field experiment.Then the experiment is laid out in a field.It should be noted that the selected experiment design and the actual field layout determine the design for an experiment as far as spatial variation in an experiment is concerned.Also, events that occur during the conduct of an experiment need to be taken into account when analyzing the data.The direction of the spatial variation may not coincide with blocking pattern used for the experiment.Several types of events can occur during the course of an experiment that determines the pattern of variation.Since an experimenter is not blessed with the knowledge of the spatial variation patterns in an experiment, it is necessary to attempt to find a response model that accounts for the variability present in an experiment.In exploratory model selection, the data analyst sets up a class of plausible response models and then selects one that best accounts for the spatial variation present in the experiment.This needs to be done for each characteristic measured in an experiment as a different response model may be appropriate for each one.
The first example is an incomplete block experiment design arrangement with the incomplete blocks laid one below the other to form an 8-row by 15-column layout within each complete block (replicate).Examples two and three were laid out in a 15-row by 12-column arrangement.The fourth example was designed as a randomized complete block experiment design but laid out in an 8-row by 7-column arrangement.In the first example two classes of response models were examined for the character weight of grain.These were differential gradients within incomplete blocks and orthogonal polynomial regression of response on row and column position and interaction of row and column regressions.Since the latter did not account for the variation in grain weight, the differential regression method was used for the other six characters measured in this experiment.The class of orthogonal polynomial regressions for rows, columns, and interactions was appropriate for finding a response model for the data of examples two, three, and four.
A standard textbook model and analysis takes the blocking in the selected experiment design into account.To illustrate, for a Latin square experiment design, the model is row effect plus column effect plus treatment effect plus error (residual).For an incomplete block experiment design, the model is replicate (complete block) effect plus treatment effect plus incom-plete block within replicate effect plus error.The differential regression method for an incomplete block, say, considers the model to be replicate effect plus treatment effect plus incomplete block within replicate effect plus regression for trend within each incomplete block.That is, the trend regressions of response on position in the incomplete blocks are considered to vary from incomplete block to incomplete block.There is no basis to consider that the regressions are the same over an entire experiment or even over each replicate.Spatial variation is seldom that patterned in field experiments, especially those on farmers' fields.The orthogonal polynomial regression method replaces a row effect with polynomial regressions of response (dependent variable) on position (independent variable) of the response, i. e., functions of row effects.Fitting a row effect with r rows is equivalent to fitting an r − 1-th regression model, hence any comments about fitting high ordered polynomials are not well conceived in the context of explaining spatial variation.The procedure is to compute all r − 1 regressions and select a subset to take account of the variation present.Those not selected are relegated to the residual (error).

Incomplete block experiment design example
An incomplete block design experiment design example was obtained from Dr. Matthew Reynolds, International Center for Maize and Wheat Improvement (CIMMYT).There were v = 120 wheat genotypes arranged in b = 15 incomplete blocks of size k = 8 in each of r = 2 replicates.Seven different responses were obtained for each of the wheat genotypes.They were weight of grain (grainwt), rust infection index, grain weight per meter squared (grm2), maturity, anthesis, total green weight (TGW), and CTD X.Two types of response models for accounting for the spatial variation present in the experiment were examined.These were trend analyses using orthogonal polynomial regression coefficients for response on position in row and in column and interactions of these regressions (regression method) and differential orthogonal polynomial regressions on position within the incomplete blocks (differential gradient method).The first method was used only for weight of grain owing to the fact that the differential gradient method was so much more effective in accounting for the spatial variation present in the experiment.This means that there were spatial trends within each of the incomplete blocks and that the trends varied from incomplete block to incomplete block.
For the response weight of grain per plot (experimental unit), the standard textbook analysis for an incomplete block design (IBD) resulted in a residual (error) mean square of 52, 395 and a coefficient of variation of 10.9%.The residual mean square was associated with 91 = 240 − 1 (correction for mean) −1 (replicate) −119 (genotype or treatment) −28 (incomplete block within replicate) degrees of freedom.Since the experiment was laid out in a 15-row (block) by 8-column arrangement within each complete block (replicate), orthogonal polynomial regression coefficients of grain weight on position in row and in column were obtained.Also, the interactions of row and column regressions were examined.The Bozivich, Bancroft , and Hartley (1956) rule used by Federer, Crossa, and Franco (1998) to select those regressions to retain as blocking variables, was used to determine which regressions in the examples are to be selected.The rule is to keep the regressions whose F -value exceeds the 25% level and relegate the others to the residual effect.This resulted in the following response model where Ri is the i-th orthogonal polynomial regression coefficient of grain weight on row position and Cj is the j-th orthogonal polynomial regression coefficient of grain weight on column position.The asterisk denotes an interaction as used in SAS codes.The residual mean square for the above response model was 54,678 with 88 = 240 − 1 − 1 − 119 − 31 (one for each of the 31 regressions or interaction of regressions) degrees of freedom and a coefficient of variation of 11.2%.Thus the best regression model using the selection rule resulted in a residual mean that was larger than the textbook analysis for an incomplete block experiment design.The differential regression method was effective in accounting for more of the spatial variation than the above as shown below.Using the response model for weight of grain grainwt = mean + replicate + genotype + block (replicate) + Cj * block (replicate) + error, the following results (with the minimum listed in bold type) for j = 1, The above demonstrates that the response model that best explains the spatial variation present in an experiment must be determined for each characteristic measured.One size may not fit all!The same response model was obtained for weight of grain and for grain weight per meter squared as it should be since they are essentially the same.It also demonstrates that there is no ordering of polynomial regressions when it comes to explaining spatial variation in experiments.The so-called "hierarchical principle" discussed by Federer (2000) used in this context is misguided and inappropriate.
The data for the above example are available upon request.A SAS GLM code for the final models of the above example is described by Federer andWolfinger (1998, 2002).If the blocking variables in the model are considered to be random effects and inter-effect information is to be recovered (and it should be), this is easily accomplished using the SAS codes.

A 15-row by 12-column experiment design example
At another site, the 120 wheat genotypes in the preceding example were included in a 15-row by 12-column design along with two check genotypes each replicated 30 times.The 120 genotypes occurred once in the experiment.The experiment design was an augmented row-column experiment design.The polynomial regression method described above is appropriate for this type of layout.For the 60 responses from the two checks, not all row, column, and check effects have non-zero solutions.The rank of the design matrix is three less than required for solution of these effects, i. e., three of the effects of a variable are set equal to zero in order to obtain a fixed effects solution.It is possible to obtain an analysis of variance, ANOVA, using SAS PROC GLM but not least squares means.Note that the degrees of freedom for genotype in the Type I (nested analysis where confounding effects are eliminated for all effects above a given one) and Type III (confounding effects eliminated for all other variables) ANOVAs should have been 121 for genotype as there were 122 genotypes in the experiment.Owing to the fact that the rank was three less than needed for solutions for all effects, this shows up in the Type III ANOVA as well as in the last line of a Type I analysis.The coefficient of variation for this analysis is 12.8% and the F -value for genotype to residual mean squares is 1.26.
To obtain genotype means, some functions of row and column effects are required.Use is made of polynomial regression of responses on position, R1 to R12, and of column responses on position, C1 to C10, and interactions of these regression coefficients.Federer, Reynolds, and Crossa (2001) considered only interactions through quartic regressions.The response model given by them is This response model resulted in a residual mean square of 6,088 with 46 = 180 − 1 − 121 − 12 (12 regressions) degrees of freedom.The coefficient of variation is 8.2% and the F -statistic for genotype to residual mean squares is 3.34.
However, if the entire range of row regression by column regression interactions is examined, considerably more of the spatial variability in this experiment can be taken into account.Using the response model resulted in a residual mean square of 1,810 with 34 = 180 − 1 − 121 − 24 (24 regressions) degrees of freedom.Interactions of high degree polynomials are required to account for the spatial variation in this experiment.The coefficient of variation is 4.5% and the F -value for genotype over residual mean squares is 11.5, a considerable change over the previous two response models.The selection of the above response model resulted in 14, 820/1, 810 = 8.2 times more replication to obtain the same standard error of a mean than would have been obtained using the textbook rowcolumn-genotype response model.Federer, Reynolds, and Crossa (2001) should have considered higher than fourth degree polynomial regression interactions.The above model resulted in 6, 088/1810 = 3.4 times more replication than limiting interactions to fourth degree.A SAS PROC GLM and MIXED code for the above response models is given by Federer and Wolfinger (2002a).

A Second 15-row by 12-column experiment design example
The 120 wheat genotypes discussed in the above two examples were grown in a 15-row by 12-column experiment design at a third site.The weight of grain for this example has been given in Federer (1998).As in the preceding example the experiment design was such that not all row, column, and genotype effects have solutions as the rank is three less than required.
If two row or two column contrasts such as R13 and R14 or C10 and C11 are set equal to zero (eliminated from the model), then solutions for effects may be obtained.However, as pointed out above, the SAS PROC GLM code will do this automatically for an ANOVA.Doing this the residual mean square was 5,630.7 with 35 degrees of freedom and a coefficient of variation of 8.5%.The F -ratio of genotype to residual mean squares was 1.47.Limiting consideration to row-column interaction regressions to fourth degree polynomials as Federer, Crossa, and Franco (1998) did, the residual mean square was 3,449.1 with 44 degrees of freedom and a coefficient of variation of 6.7%.The F-statistic for genotypes was 2.44.Their response model was where Cj is the j-th polynomial regression coefficient of grain weight on column position and Ri is the i-th polynomial regression of grain weight on row position.If interactions of all column regressions by row regressions are screened by the rule used above, the resulting response model equation is Using this response model resulted in a residual mean square of 1,081.4 with 38 degrees of freedom and a coefficient of variation of 3.7%.The genotype F -ratio increased to 8.04.Thus, the preceding response model resulted in a residual mean square with 38 degrees of freedom and 5, 630.7/1, 081.4 = 5.2 times more replication than the textbook row-column-genotype response model and 3, 449.1/1,081.4 = 3.2 times more replication than the Federer, Crossa, and Franco (1998) response model.It is to be noted that the more patchy the spatial variation, the higher will be the degree of the polynomial regression interactions required to account for this.The above response model used fewer degrees of freedom, 20, for blocking variables than did the row-column-genotype response model, 23.

An 8-row by 7-columu experiment design
An experiment described in Federer and Schlottfeldt (1954) was designed as a randomized complete block experiment design (RCBD) with v = 7 treatments and r = 8 complete blocks.However, the experiment was laid out in an 8-row by 7-column arrangement, RCD.Owing to several sandy patches in the experimental area and to unfavorable moisture conditions, there was considerable spatial variation present in the experiment.The RCBD ANOVA resulted in a residual man square of 30,228.2 with 42 = 56 − 1 − 7 − 6 degrees of freedom, a coefficient of variation of 17.2%, and an F -value for treatments of 1.51.The RCD ANOVA produced a residual mean square of 7,351.8 with 36 degrees of freedom, a coefficient of variation of 8.5%, and an F -value for treatments of 2.71.The response model as used by Federer, Crossa, and Franco (1998) and limiting the investigations to interactions of fourth degree row and column regressions, produced a residual mean square of 4,204.5 with 33 degrees of freedom, a coefficient of variation of 6.4%, and an F -value for treatments of 6.36.Considering interactions of all polynomials resulted in the following response model equation: That is, nine interaction terms were added to the previous model.The resulting residual mean square is 1,320.4with 24 degrees of freedom.The coefficient of variation is 3.6% and the F -value for treatments is 20.49.As may be observed, tremendous differences exist between the analyses for the different response models.A standard textbook approach would use the RCBD analysis, and consideration of the spatial layout of the experiment would use the RCD analysis.Consideration of differential gradients would result in the response model The residual mean square for this model was 11,309.9 with 18 degrees of freedom, a coefficient of variation of 10.5%, and an F -value for treatments of 3.78.This model used 31 = 7 + 3(8) degrees of freedom for blocking variables, leaving only 18 for the residual.Furthermore, it was not as effective as the row-column-treatment model for controlling spatial variation.Using the next to last model above, was quite effective in accounting for the spatial variation in this experiment and effectively resulted in 30, 228.2/1, 320.4 = 22.9 times more replication than the RCBD analysis and 7, 351.8/1, 320.4 = 5.6 times more replication than the RCD analysis.A SAS PROC GLM and PROC MIXED code for the above models is given by Federer and Wolfinger (2002b).

Discussion
The above examples utilized a fixed effects approach to exploratory model selection.Federer and Wolfinger (2000) have presented two random effects procedures for model selection.A comparison of resulting models using the three procedures with above examples could be made.It is fairly certain that the resulting models would differ.Until the properties of the random effects selection procedures are known, this will not be done.It is known that the Bozivich, Bancroft, and Hartley (1956) procedure has only a small effect on the Type I error.It is possible that F at the 25% level is not optimal for reducing the effect on Type I errors.With present computing power, this could be investigated.
The model selection procedure utilized the fixed effects analyses.When obtaining treatment means, one should recover the information from the random blocking effects.As noted from remarks by several anonymous referees, they have a difficult time thinking of regression and gradient coefficients as random blocking effects.They appear to have no difficulty with considering row, column, and block effects as random, but they appear to fail to appreciate the fact that the regressions and gradients are merely functions of row and column effects that are random effects.Hence, in obtaining adjusted treatment means, all blocking effects should be considered as random effects and the information contained in them needs to be recovered in order to utilize efficient procedures.
Different models were obtained for each of the characters analyzed in the first example.This means that an experimenter should perform exploratory model selection for each characteristic being analyzed.Use of computer programs such as those given in the references, make this is relatively simple matter.
The fixed methods of regression and gradients used in this investigation have known degrees of freedom for the various parameters used in a response model.Several other procedures such as smoothing, Kriging, nearest neighbor, and autoregression have been proposed for spatial analyses (See Federer, Newton, and Altman, 1997, e.g.).The degrees of freedom for each of the parameters used in these methods are usually unknown.Hence, it is difficult to compare their ability to explain spatial variation in comparison to regression and gradient procedures.
Orthogonal polynomial regression was used in the above analyses.This involves using centered values for the independent (covariate) variable, position.It should be noted that interactions of non-centered covariates will not be the same as those from centered covariates.Also, instead of using orthogonal polynomial regressions, Fourier regression may be more appropriate in some situations, i.e., when the spatial variation is cyclical in nature.
Some statisticians, e.g., Gilmour (2000), appear to believe that the above exploratory model investigation is "post blocking that has gone too far."If the variation can be accounted for and if there are sufficient degrees of freedom, say 20-30, associated with the residual mean square, there should be no reason why a data analyst should not use procedures such as those described herein.Others believe that regression coefficients and gradients should always be considered as fixed effects.They appear to fail to realize that as used herein, they are functions of random variable effects and hence should be considered to be random effects.In field layouts, there are no valid reasons to consider that there will be a single regression and that all variation follows an orderly and systematical pattern.Even though an experimenter may try to select a uniform area in which to conduct the experiment, this is not always possible, e.g., conducting experiments on farmer's fields.
a residual mean square of 31,546 (listed in bold type) and a coefficient of variation of 8.5%.The ratio of this mean square to the IBD mean square is 31, 546/52, 395 = 0.602, which is 40% smaller than the standard textbook incomplete block analysis residual mean square.The addition of this term 63540 for the IBD response model to 1.99125 for a 24.4% reduction.Adding the terms C2*block (replicate) and C7*block (replicate) to the IBD model reduced the residual mean square to 1.55874 for a 40.9% reduction over that for the IBD residual mean square.The residual mean square for grain weight per meter squared, grm2, using the IBD response model was 802,741 and the coefficient of variation was 10.8%.Adding the term Cj*block (replicate) to this response model resulted in the following residual mean squares and coefficients of variation:Adding the two terms C1*block (replicate) and C7*block (replicate) to the IBD response model for grm2 resulted in a decrease in the residual mean square of 100(1 − 518, 515/802, 741) = 35.4%.The term C7*block (replicate) resulted in a decrease in the residual mean square of 100(1 − 518, 515/666, 921) = 22.3%.For the character anthesis, the IBD response model resulted in a residual mean square of 10.55307 and a coefficient of variation of 5.6%.When the term Cj*block (replicate) was added to the IBD response model, the following residual mean squares and coefficients of variation were obtained: all possible pairs added to the IBD response model resulted in selecting C1*block (replicate) plus C2*block (replicate) as the response model.Despite the small coefficient of variation, 5.6%, for the IBD model, the addition of the two terms C1*block (replicate) and C2*block (replicate) to the IBD response model resulted in a decrease in the residual mean square of 100(1 − 1.39405/10.55307)=86.8%.The addition of C2*block (replicate) resulted in a decrease of 100(1 − 1.39405/2.74642)=49.2% in the residual mean square.The IBD response model for the character maturity resulted in a residual mean square of 0.746855 and a coefficient of variation of 4.2%.Adding the term Cj*block (replicate) to the IBD response model resulted in the following residual mean squares and coefficients of variation: Cj*block (replicate) to the IBD response model resulted in adding C3*block (replicate) and C6*block (replicate) to the IBD model.Adding the term C2*block (replicate) to the IBD response model reduced the residual mean square to 0.642376 from 0.746855, or a reduction of 14.0%.Adding the two terms C3*block (replicate) and C6*block (replicate) further reduced the residual mean square to 0.607233, or a reduction of 18.7% over that obtained for the IBD response model.The IBD response model for the character CTD X resulted in a residual mean square of 0.146648 and a coefficient of variation of 13.8%.Adding Cj*block (replicate) to the IBD response model resulted in the following residual mean squares and coefficients of variation: Adding the pair C2*block (replicate) and C4*block (replicate) decreased the residual mean square by 100(1−0.068188/0.146648)=53.5% with the C4*block (replicate) term accounting for 100(1−0.068188/0.096373)=29.2% of the decrease.The response models resulting in minimum residual mean squares for each of the seven characteristics reported in this experiment are given below: where the first row is Cj, the second row is residual mean square, and the third row is the coefficient of variation.Adding the term C1*block (replicate) to the incomplete block response model was effective in accounting for a sizeable portion of the residual variation.These residual mean squares are associated with 61 = 91 − 30 (one for each of the 30 regressions in each of the 30 incomplete blocks) degrees of freedom.Since there are sufficient degrees of freedom to search further, a second term of Ch*block (replicate) for h not equal to j, is added to the above response model equation.All possible pairs were investigated.The residual mean square results with the associated coefficient of variation obtained are the minimum listed in bold type as given below on the next page.Adding the term C7*block (replicate) to the response model resulted in over the C1*block (replicate) term is sizeable, i.e., 31, 546/41, 604 = 0.758, or a 24% reduction.Using the textbook analysis for the incomplete block design (IBD) would require 52, 395/31, 546 = 1.66 times more replication than the differential gradient model with two terms in the response model.Also, even with allocating an additional 30 degrees of freedom to spatial variation, there are still 31 degrees of freedom associated with the residual mean square.For the response rust index, the IBD response model gave