Substitutes for the non-existent square lattice designs for 36 varieties

Square lattice designs are often used in trials of new varieties of various agricultural crops. However, there are no square lattice designs for 36 varieties in blocks of size six for four or more replicates. Here we use three different approaches to construct designs for up to eight replicates. All the designs perform well in terms of giving a low average variance of variety contrasts. Supplementary materials are available online.


Introduction
In variety-testing programmes, later-stage trials can involve multiple replications of up to 100 varieties: see Patterson, Williams and Hunter (1978). Denote the number of varieties by v. Even at a well-run testing centre, variation across the experimental area makes it desirable to group the plots (experimental units) into homogeneous blocks, usually too small to contain all the varieties. As R. A. Fisher wrote in a letter in 1938, ". . . on any given field agricultural operations, at least for centuries, have followed one of two directions", so that variability among the plots is well captured by blocking in one or both of these directions, with no need for more complicated spatial correlations: see Fisher et al. (1990, p. 270). Thus, on land which has been farmed for centuries, or where plots cannot be conveniently arranged to allow blocking in two directions (rows and columns), it is reasonable to assume the following model for the yield Y ω on plot ω: Here V (ω) denotes the variety planted on ω and B(ω) denotes the block containing ω. The variety constants τ i are the unknown parameters of interest, and the block constants β j are unknown nuisance parameters. The quantities ε ω are independent identically distributed random variables with zero mean and common (unknown) variance σ 2 . For management reasons, it is often convenient if the blocks can themselves be grouped into replicates, in such a way that each variety occurs exactly once in each replicate. Such a block design is called resolvable. Let r be the number of replicates. Yates (1936Yates ( , 1937 introduced square lattice designs for this purpose. In these, v = n 2 for some positive integer n, and each replicate consists of n blocks of n plots. The design is constructed by first listing the varieties in an abstract n × n square array S. The rows of S form the blocks of the first replicate, and the columns of S form the blocks of the second replicate. If r > 2 then r −2 mutually orthogonal Latin squares L 1 , . . . , L r−2 of order n are needed. This is not possible unless r ≤ n+1: see Street and Street (1987, Chapter 6). For replicate i, where i > 2, superimpose Latin square L i−2 on the array S: the n positions where any given letter of L i−2 occurs give the set of varieties in one block. Orthogonality implies that each block has one variety in common with each block in each other replicate. Thus these designs belong to the class of affine resolvable designs defined by Bose (1942), and this construction is a special case of that given by Bailey, Monod and Morgan (1995). Moreover, all pairwise variety concurrences are in {0, 1}, where the concurrence of varieties i and j is the number of blocks in which varieties i and j both occur. If r = n + 1 then all pairwise concurrences are equal to 1 and so the design is balanced.
Equireplicate incomplete-block designs are typically assessed using the A-criterion: see Shah and Sinha (1989). Denote by Λ the v × v concurrence matrix : its (i, j)-entry is equal to the concurrence of varieties i and j, which is r when i = j. The scaled information matrix is I − (rk) −1 Λ, where I is the identity matrix and k is the block size. The constant vectors are in the null space of this matrix. The eigenvalues for the other eigenvectors, counting multiplicities, are the canonical efficiency factors. Denote their harmonic mean as A. (John and Williams (1995) call this E, but many authors, including several cited in Section 5, use E to denote the smallest canonical efficiency factor.) Under model (1), the average variance of the estimator of a difference τ i − τ j between two distinct varieties is 2σ 2 /(rA). If the variance in an ideal design with the same number of plots but no need for blocking is σ 2 0 , then this average variance would be 2σ 2 0 /r. Hence A ≤ 1, and a design maximizing A, for given values of v, r and k, is called A-optimal. Cheng and Bailey (1991) showed that if r ≤ n + 1 then square lattice designs are A-optimal (even over non-resolvable designs).
The class of square lattice designs also has some practical advantages. Adding or removing a replicate gives another square lattice design, which permits last-minute changes in the planning stage. It also means that if a whole replicate is lost (for example, if heavy rain during harvest flattens the plants in the last replicate) then the remaining design is A-optimal for its size.
If n ∈ {2, 3, 4, 5, 7, 8, 9} then there is a complete set of n − 1 mutually orthogonal Latin squares of order n: see Street and Street (1987, Chapter 6). These give square lattice designs for n 2 varieties in rn blocks of size n for r ∈ {2, . . . , n + 1}. However, there is not even a pair of mutually orthogonal Latin squares of order 6, so square lattice designs for 36 varieties are available for two or three replicates only. This gap in the catalogue of good resolvable block designs is pointed out in many books: for example, Cochran and Cox (1957); John and Williams (1995). Patterson and Williams (1976a) used computer search to find an efficient resolvable design for 36 varieties in four replicates of blocks of size six. All pairwise variety concurrences are in {0, 1, 2}. It has A = 0.836, which compares well with the unachievable upper bound of 0.840 for the non-existent square lattice design.
In this paper we present three new methods of contructing efficient resolvable block designs for 36 varieties in 6r blocks of size six, for r ∈ {4, . . . , 8}. These methods are in Sections 3-5. In each case, Supplementary Material gives the design as a plain text file, which can easily be imported into a spreadsheet or statistical software.
The concurrence graph of an incomplete-block design has a vertex for each variety. The number of edges between vertices i and j is equal to the concurrence of varieties i and j. Although the methods in Sections 3-5 are very different, their designs for eight replicates all have the same concurrence graph, which we describe in Section 2. The final sections compare the new designs and discuss further work.

The Sylvester graph
The Sylvester graph Σ is a graph on 36 vertices with valency 5. See Brouwer, Cohen and Neumaier (1989) and R. F. Bailey (2019). Here we give enough information about it to show how the designs in Section 3 are constructed, using the approach in Cameron and van Lint (1991, Chapter 6).
Consider the complete graph K 6 . It has a set A of six vertices, labelled 1 to 6. There is an edge between every pair of distinct vertices. Let B be this set of edges. A 1-factor is a partition of A into three edges (subsets of size two). For example, one 1-factor consists of the pairs {1, 2}, {3, 6} and {4, 5}. For brevity, we write this in the slightly non-standard way 12|36|45 in Table 1. Let C be the set of 1-factors. A 1-factorization is a set of five elements of C with the property that each edge is contained in just one of them. For example, d 1 in Table 1 is a 1-factorization; here we use the symbol || to separate the five 1-factors in d 1 . Let D be the set of 1-factorizations.
It was shown by Sylvester (1844) that (a) |A| = |D| = 6 and |B| = |C| = 15; (b) any two elements of D share exactly one element of C.
(Sylvester used the terms duads, synthemes and synthematic totals for edges, 1-factors and 1-factorizations respectively.) Table 1 shows the six 1-factorizations, labelled as d 1 to d 6 . Each of these can be considered as a schedule for a tournament involving six teams which takes place over five weekends so that each pair of teams meets exactly once. Now we construct the Sylvester graph Σ as follows. The vertex set consists of the cells of the 6 × 6 square array S with rows labelled by the elements of A and columns labelled by the elements of D. Given distinct d i and d j in D, the unique 1-factor they have in common defines six edges in Σ, each joining a vertex in column d i to a vertex in column d j . For example, d 3 and d 4 have the one-factor 12|34|56 in common, so we put edges between Thus each vertex in Σ is joined to five other vertices, one in each other row and one in each other column. Figure 1(b) shows the five edges at the vertex a = (3, d 3 ). We shall call this set of six vertices the starfish centred at a.
It can be shown that the graph Σ has no triangles or quadrilaterals. One consequence of this is that, given any vertex, the vertices at distances one and two from it in the graph are precisely all the other vertices in different rows and different columns.
Denote by Adj(Σ) the adjacency matrix of the graph Σ. We call a block design for 36 varieties in 48 blocks of size six a Sylvester design if there is a permutation of the varieties that takes the concurrence matrix for the design to 7I + J + Adj(Σ), where J is the all-1 matrix. This means that the concurrences are 2 on each edge of Σ and 1 for every other pair of varieties. In particular, a Sylvester design is a regular-graph design, as defined by John and Mitchell (1977), where the graph is the Sylvester graph Σ.
Two block designs (for the same set of varieties) are isomorphic if one can be converted into the other by a permutation of varieties and a permutation of blocks. An isomorphism from a block design to itself is called an automorphism.
If two block designs are isomorphic then their canonical efficiency factors are the same and their automorphism groups have the same order, but nei- ther converse need be true. In particular, all Sylvester designs have the same canonical efficiency factors, and hence the same value of the A-criterion, but they are not all isomorphic. Of the three that we construct in this paper, no two are isomorphic, as we discuss in Section 6.
3 New designs constructed from the Sylvester graph

The new designs
Figure 1(a) shows that if we choose two different vertices in the same column then their starfish will not overlap. Thus each column gives what we call a galaxy of six starfish, which together include every vertex just once. In other words, we can think of a galaxy as a Latin square of order 6. Figure 2 shows the galaxy of starfish centred on vertices in column d 3 . Just as with a square lattice design, we can identify the varieties with the 36 vertices and use this Latin square to construct a single replicate of six blocks of size six.
However, unlike in a square lattice design, the Latin squares defined by different columns are not orthogonal to each other. If two vertices are joined by an edge then they both occur in the two starfish which they define. Thus if we use galaxies of starfish from two or more columns then some variety concurrences will be bigger than one. On the other hand, a consequence of the lack of triangles and quadrilaterals is that if two or more galaxies are used as replicates then there is no other way that two varieties can concur in more than one block.
We therefore propose the following resolvable designs. The design Γ r consists of the galaxies of starfish from r columns, where 0 ≤ r ≤ 6; Γ 0 (a design with no blocks) and Γ 1 (a disconnected design) are used in the following constructions, but are not themselves suitable designs. For 1 ≤ r ≤ 7, the design Γ R r consists of Γ r−1 together with another replicate whose blocks are the rows of S, while the design Γ C r consists of Γ r−1 together with another replicate whose blocks are the columns of S. The design Γ C 7 was used by Bailey, Cameron and Nilson (2018). For 2 ≤ r ≤ 8, the design Γ RC r consists of Γ R r−1 together with another replicate whose blocks are the columns of S. In particular, Γ RC 2 is the square lattice design whose blocks are the rows and columns.
The automorphisms of Σ consist of the symmetric group S 6 acting simultaneously on rows and columns of the array, as well as a further involution transposing it. It follows that, for a design consisting of m galaxies (possibly with rows, and possibly with columns), it does not matter which m galaxies we choose.
When r = 2 then Γ RC 2 , Γ R 2 and Γ C 2 are square lattice designs, and hence A-optimal, but Γ 2 is not. When r = 3 then Γ RC 3 is a square lattice design, and hence A-optimal, but none of the others is. When r ≥ 4 then none of the designs is a square lattice design, so we need to calculate the canonical efficiency factors, and hence A.
As discussed in Section 6.1, for each value of r the designs Γ R r and Γ C r have the same canonical efficiency factors, so we do not include Γ R r in further comparisons.
Another useful consequence of the lack of triangles and quadrilaterals in Σ is that the relations 'same row', 'same column', 'joined in the graph' and 'other' form a 4-class association scheme on the set of vertices. It follows that Γ 6 , Γ C 7 and Γ RC 8 are partially balanced with respect to this association scheme, and so their canonical efficiency factors can be calculated using the methods in Bailey (2004). Table 2 shows the results. In fact, Γ 6 and Γ RC 8 are also partially balanced with respect to the 3-class association scheme obtained by merging the classes 'same row' and 'same column'. Moreover, Γ RC 8 is a Sylvester design. For all the other designs, we calculated A as an exact rational number by using the DESIGN package (Soicher, 2019) in GAP (The GAP Group, 2019). The method used for such exact calculation of block design efficiency measures is described in Appendix B of Soicher (2013a). These results were verified by using GAP to find the exact Moore-Penrose inverse of the scaled information matrix, calculate its trace, divide this by 35, and then invert this as an exact rational number. Table 3 shows the results to four decimal places. This shows that, apart from the square lattice designs Γ RC 2 and Γ C 2 , the design Γ RC r always beats Γ C r and Γ r . Moreover, Γ RC 4 does very slightly better than the design found by Patterson and Williams (1976a).
The final column of Table 3 shows the value of A for a square lattice  Section 3 Section 4 Section 5 square design. This exists only for r = 2 and r = 3, when Γ RC r is an example. For 4 ≤ r ≤ 7 there is no square lattice design, so the value shown gives an unachievable upper bound; in every case, A for Γ RC r comes very close to this. The computer search algorithm used by Patterson and Williams (1976a) to obtain the efficient design for 36 varieties has been extensively developed, both in the range of design types that can be constructed and in the algorithmic approach. Significant improvements in computer speed have also facilitated search procedures. CycDesigN Version 6.0 (VSNI, 2016) is a computer package for the generation of optimal or near-optimal experimental designs, as measured by the A-criterion. The package has been written in Visual C++ and uses simulated annealing in the design search process. CycDesigN can be used to construct efficient resolvable block designs for 36 varieties in blocks of size six with a range of values for r. Hence running CycDesigN separately for r = 3 through 8 gives designs Θ r with the results in Table 3. For r = 3, . . . , 7, the design Θ r has pairwise variety concurrences in {0, 1, 2}, while Θ 8 has concurrences in {1, 2}. In fact, Θ 8 is a Sylvester design. For r = 4, the improvement from Patterson and Williams (1976a), namely A = 0.836 to that in Table 3 (A = 0.839) is representative of overall developments in computer speed and search methods throughout the years. Because Θ r is not constructed by simply omitting a replicate from Θ r+1 , we do not show all these designs here. Figure 4 shows Θ 8 . Plain-text versions

New designs constructed from semi-Latin squares
Let ∆ be a resolvable block design for 36 varieties in 6r blocks of size 6. Its dual design ∆ ′ is obtained by interchanging the roles of blocks and varieties, so it has 6r varieties in 36 blocks of size r. If the varieties of ∆ are identified with cells of the 6 × 6 square array S, then the blocks of ∆ ′ also form a 6 × 6 square array. When each variety in ∆ ′ occurs exactly once in each row and once in each column then ∆ ′ is called a (6 × 6)/r semi-Latin square: see Yates (1935) and Preece and Freeman (1983). The term orthogonal multiarray is also used: see Brickell (1984). One way of constructing such a semi-Latin square is to superpose r Latin squares with disjoint alphabets. Not all semi-Latin squares arise in this way, but the resolvability of ∆ forces the r replicates to be a collection of r Latin squares when ∆ ′ is a semi-Latin square.
For 0 ≤ r ≤ 6, denote by ∆ r the design for 36 varieties in 6r blocks of size six given by the Latin squares L 1 , . . . , L r ; ∆ 0 is a design with no blocks, but, for 0 < r ≤ 6, ∆ r is the dual of the design called X r by Soicher (2013a). Table 2 of Soicher (2013a) gives A ′ for X 2 , . . . , X 6 : from this, A can be calculated from Equation (2).
As in Section 3, we can add to ∆ r−1 another replicate whose blocks are the rows of S, to obtain ∆ R r , or we can add another replicate whose blocks are the columns of S, to obtain ∆ C r . Adding both of these extra replicates to ∆ r−2 gives ∆ RC r . Unlike the situation in Section 3, choosing a different r-subset of {L 1 , . . . , L 6 } may give a design with a value of A different from that for ∆ r , but computation shows that, in every case, the highest value of the A-criterion arises from taking {L 1 , . . . , L r } as our r-subset. Table 3 shows the values of the A-criterion for ∆ r , ∆ C r and ∆ RC r , calculated exactly using the DESIGN package and rounded to four decimal places. We found that, for each r = 2, . . . , 7, the canonical efficiency factors of ∆ R r are the same as those of ∆ C r , and so their A-values are the same. Note that, just as in Section 3, ∆ RC r always beats ∆ R r , ∆ C r and ∆ r , apart from the fact that ∆ RC 2 , ∆ R 2 and ∆ C 2 are all square lattice designs. Note also that, to four decimal places of the A-criterion, ∆ RC r is always at least as good as Γ RC r , and sometimes better. Figure 5 shows the design ∆ RC 8 , with the replicates defined in the order columns, rows, L 1 , . . . , L 6 . The varieties are numbered 1 to 6 in row 1, then 7 to 12 in row 2, and so on. A plain-text version of this design is available in the Supplementary Material. For a design with r replicates, use the first r of the given replicates.
Although the design X 6 of Soicher (2013a) was not constructed using the Sylvester graph, it turns out that ∆ RC 8 is a Sylvester design, and so has the same canonical efficiency factors as Γ RC Now the square lattice design Γ R 2 is isomorphic to Γ C 2 , and we found, using the DESIGN package, that Γ R 7 is isomorphic to Γ C 7 . For 3 ≤ r ≤ 6, it turns out that the designs Γ R r and Γ C r are not isomorphic, but they have the same canonical efficiency factors.
We found that, for r = 2, 3, 5, 7, ∆ R r is isomorphic to ∆ C r . For r = 4, 6, the designs ∆ R r and ∆ C r are not isomorphic, but they have the same canonical efficiency factors. We do not yet know a theoretical reason for this.  the first method started from the Sylvester graph, neither of the others did; indeed, the second method used numerical optimization. CycDesigN also calculates upper bounds for the A-criterion. For r ≤ 7 these are same as those shown in the final column of Table 3. For r = 8 it gives an upper bound of 0.854931, compared to the A-criterion for all the Sylvester designs, which is equal to 0.854929 to six decimal places. The proximity of these values gives further support to the conjecture that the Sylvester designs are A-optimal.
For a given value of r, which design should be used in practice? Table 3 shows that the A-values for Γ RC r , Θ r and ∆ RC r are extremely close, but Θ r is always at least as good as the other two. If the user is concerned about the possible loss of one replicate, then Table 4 shows that Θ r is still at least as good as the other two when r ≤ 7 but Γ RC 8 might be preferred to Θ 8 . If the 36 varieties are replaced by 36 treatments consisting of all combinations of two factors with six levels each, and r ≤ 6, then levels of these two factors can be identified with the blocks of any two replicates which together form a square lattice design. For the designs in Sections 3 and 5, these are either rows and columns, or one of rows and columns combined with any other replicate. Table 3 shows that the second possibility gives higher values for the A-criterion.
Bibliographic note An extended abstract for this paper is in .
crete Mathematics) is gratefully acknowledged.