# Automatic generation of generalised regular factorial designs

Andre Kobilinsky, Hervé Monod, R.A. Bailey

► To cite this version:

Andre Kobilinsky, Hervé Monod, R.A. Bailey. Automatic generation of generalised regular factorial designs. Computational Statistics and Data Analysis, 2017, 113, pp.311-329. 10.1016/j.csda.2016.09.003. hal-01608689

HAL Id: hal-01608689

https://hal.science/hal-01608689

Submitted on 25 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Automatic generation of generalised regular factorial designs☆

André Kobilinskyª, Hervé Monodª,*, R.A. Baileyᵇ,ᶜ

## ARTICLE INFO

**Article history:**

Received 29 January 2016

Received in revised form 5 September 2016

Accepted 6 September 2016

Available online xxxx

**Keywords:**

Backtrack algorithm

Design key

Hierarchy constraint

Ineligible factorial term

Model specification

Multi-stratum design

## ABSTRACT

The R package *planor* enables the user to search for, and construct, factorial designs satisfying given conditions. The user specifies the factors and their numbers of levels, the factorial terms which are assumed to be non-zero, and the subset of those which are to be estimated. Both block and treatment factors can be allowed for, and they may have either fixed or random effects, as well as hierarchy relationships. The designs are generalised regular designs, which means that each one is constructed by using a design key and that the underlying theory is that of finite abelian groups. The main theoretical results and algorithms on which *planor* is based are developed and illustrated, with the emphasis on mathematical rather than programming details. Sections 3–5 are dedicated to the elementary case, when the numbers of levels of all factors are powers of the same prime. The ineligible factorial terms associated with users’ specifications are defined and it is shown how they can be used to search for a design key by a backtrack algorithm. Then the results are extended to the case when different primes are involved, by making use of the Sylow decomposition of finite abelian groups. The proposed approach provides a unified framework for a wide range of factorial designs.

© 2016 The Author(s). Published by Elsevier B.V. This is an open access article under the

CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

## 1. Introduction

Factorial designs that may include several block and treatment factors date back to the pioneering work of Fisher and Yates at Rothamsted Experimental Station (Yates, 1933, 1937; Fisher, 1942), followed by Finney (1945) and Bose (1947). Since then, the construction of fractional designs has been a constantly active field of research in the theory of design of experiments. It has also been widely applied in many different application areas, including food research, biology, industry, and – more recently – computer experiments.

The designs we are interested in are known today as regular factorial designs. Their construction, which is based on algebra and group theory, gives a large class of orthogonal factorial designs, including as special cases block and row-column designs (Yates, 1937) as well as fractional designs (Finney, 1945). Regular fractional designs became a standard method of construction very early. Their main principles in standard cases have been explained in numerous text books including classics (Kempthorne, 1952; Cochran and Cox, 1957) and more recent ones (Ryan, 2007; Bailey, 2008; Morris, 2011; Cheng, 2014).

<footnote>ª MalAGE, INRA, Université Paris-Saclay, 78350 Jouy-en-Josas, France</footnote>

<footnote>ᵇ School of Mathematics and Statistics, University of St Andrews, St Andrews KY16 9SS, UK</footnote>

<footnote>ᶜ School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E1 4NS, UK</footnote>

<footnote>☆ This paper also includes Supplementary material providing proofs (Appendix A) (and additional information on the algorithms (Appendices B, C).</footnote>

<footnote>* Correspondence to: Hervé Monod, MalAGE bât. 210, INRA, Domaine de Vilvert, 78350 Jouy-en-Josas, France. Fax: +33 1 34 65 22 17.</footnote>

<footnote>E-mail address: herve.monod@inra.fr (H. Monod).</footnote>
Of course, there are situations where non-regular factorial designs are needed. Examples include incomplete-block designs where the block size does not divide the number of treatments, and orthogonal arrays, such as those obtained from Hadamard matrices (Hedayat et al., 1999) and by computer intensive algorithms (Kuhfeld and Tobias, 2005) or those given by Bose and Bush (1952) and Addelman and Kempthorne (1961), whose number of experimental units does not divide the number of treatments. See Grömping and Hu (2014) and Grömping and Bailey (2016) for further discussion of the definition of regularity. However, regular designs still have wide applicability, and these are what the software planor deals with.

An important notion in fractional designs is resolution (Box and Hunter, 1961a,b). If R is a positive integer, a fraction of resolution R allows the estimation of all factorial effects up to interactions of order strictly smaller than R/2, assuming that all interactions of order strictly larger than R/2 are zero. More discriminating criteria such as minimum aberration (Fries and Hunter, 1980) or maximum estimation capacity (Cheng and Mukerjee, 1998) have been developed. The construction of optimal designs with respect to these criteria is still an active field of research, which includes analytical as well as algorithmic issues. However, resolution and aberration imply that the model of interest be symmetric with respect to all factors.

Most work in this area thus deals with problems which are highly symmetric with respect to the factors and models of interest. However, there is also a need to find generic and user-friendly methods of construction adapted to much more flexible problem specifications, allowing for unconstrained numbers of levels and flexible model assumptions.

The construction of a regular fraction involves two steps: first, finding defining relationships or generators of the fraction ensuring that factorial effects of interest will be estimable; second, generating the actual design. Despite its apparent simplicity in standard cases, the first step still represents a major challenge in general situations. It has been applied and programmed mainly in the case of symmetric designs, which require all factors to have the same number of levels or, at least, numbers of levels that are powers of the same prime.

Algorithms were developed and studied in the 1970s and 1980s (Patterson, 1965, 1976; Bailey et al., 1977; Franklin and Bailey, 1977; Patterson and Bailey, 1978; Bailey, 1985; Franklin, 1985) and some of them implemented in statistical packages (SAS, 2010; Payne, 2012; Grömping, 2014). Lewis (1982) tabulated generators for asymmetric factorial designs with resolution 3. This paper extends these approaches to a generalised class of regular designs. We follow the theoretical framework of Kobilinsky (1985) and Kobilinsky and Monod (1991, 1995), although we occasionally modify a term or notation for simplicity. We also make intensive use of pseudofactors (Monod and Bailey, 1992). The generalised class includes symmetric as well as asymmetric designs, with estimation planned in one or more strata. For example, our results apply to the construction of generalised split-plot or criss-cross designs as well as the more usual fractional designs. The method allows the user to define the model and specify what should be estimated. The approach and results have been implemented in the R package planor (Kobilinsky et al., 2012; Monod et al., 2012) based on the initial APL version by Kobilinsky (2005).

## 2. Overview of the design search issue

### 2.1. Examples

We start with three examples to illustrate the diversity of situations we want to consider, all based on real applications from our consulting experience.

We follow standard practice in using the same notation A for a factor and for its main effect, and the notation A.B for the interaction between factors A and B. Hierarchy relationships are denoted as in Bailey (2008, Chapter 10). In particular:

* B ≪ A means that factor B is or must be nested in factor A, in other words B must be finer than or equivalent to A and, reciprocally, A must be coarser than or equivalent to B, so that each level of B occurs with a single level of A;

* A ∧ C denotes the product or infimum of A and C, which is the factor whose levels are the combinations of levels of A and C.

Note the following two properties:

* A ∧ C ≪ A and A ∧ C ≪ C;

* if A ≪ B and A ≪ C, then A ≪ B ∧ C.

**Example 1.** There are four treatment factors $F_1, F_2, F_3, F_4$ with 6, 4, 3, 4 levels respectively. A complete factorial design would require 288 experimental units but we assume this is much larger than possible, so that a smaller fractional design is looked for.

The experimenter intends to analyse the data using a model that consists of the general mean, the four main effects of $F_1, F_2, F_3, F_4$ and the interaction $F_1.F_3$, all with fixed effects, with good reasons to consider the other interactions as negligible. The factorial terms he or she wants to estimate are the four main effects only. The factorial terms in the model and those that must be estimated are listed in the following sets $\mathcal{M}$ and $\mathcal{E}$:

$$ \mathcal{M} = \{\mu, F_1, F_2, F_3, F_4, F_1.F_3\} \quad \text{and} \quad \mathcal{E} = \{F_1, F_2, F_3, F_4\} = \mathcal{M} \setminus \{\mu, F_1.F_3\}, $$

where $\mu$ denotes the general mean.
**Example 2.** A row-and-column design has to be constructed with two columns (factor C), three rows (factor R) and two units in each of the six cells defined by a row and a column. There are three treatment factors: two 2-level factors D, E and one 3-level factor A. The experimenter wants to estimate the interactions D.A, E.A, considering a model that includes the general mean, row, column and cell effects as well as all main effects of treatment factors and all interactions between two treatment factors. The sets $\mathcal{M}$ and $\mathcal{E}$ now are:

$$ \underline{\mathcal{M}} = \{\mu, C, R, C.R, D, E, A, D.A, E.A, D.E\} \text{ and } \underline{\mathcal{E}} = \{D.A, E.A\}. $$

An additional constraint is that, for practical reasons, the factor A must remain constant on each row. We call this a hierarchy constraint, imposing that factor R be nested in factor A, and denote it by $R \prec A$.

**Example 3.** The experimental units consist of four blocks, each containing two subblocks of four units. This structure can be described by three factors P, Q, U with four, two and four levels respectively. The levels of P define the blocks, the levels of the infimum $P \wedge Q$ define the subblocks, and the levels of the infimum $P \wedge Q \wedge U$ determine the units. In addition, there are four treatment factors A, B, C, D with two levels. There is again a hierarchy constraint: we assume that the levels of A cannot be varied between the four units of a given subblock, which is denoted by $P \wedge Q \prec A$.

The experimenter is interested in the main effects and two-factor interactions of A, B, C, D. Thus the set $\mathcal{E}$ is

$$ \underline{\mathcal{E}} = \{A, B, C, D, A.B, A.C, A.D, B.C, B.D, C.D\}. $$

In addition to $\mathcal{E}$, the model must contain the block and subblock effects and the general mean. However, the hierarchy constraint means that the main effect of A is necessarily confounded with subblock effects so that it cannot be estimated in a model with fixed subblock effects.

If a proper randomisation is performed, block and subblock effects can instead be considered as centred random effects and the analysis can be decomposed into three levels of variability, or strata (Bailey, 2008), associated with P, $P \wedge Q$ and $P \wedge Q \wedge U$ respectively. The blocks stratum consists of all contrasts between blocks. The between-subblocks stratum consists of all contrasts between subblocks which are orthogonal to blocks. The bottom stratum, which includes residual variability only, consists of all contrasts orthogonal to subblocks.

The objective now is to estimate all effects of interest except the main effect A in the bottom stratum, and to estimate the main effect A in the between-subblocks stratum. These requirements can be described by using two pairs of model-estimate lists, one for each stratum where an effect will be estimated. The first pair deals with the bottom stratum by pretending that P and $P \wedge Q$ have fixed effects and omitting A from the effects to be estimated. Thus

$$ \underline{\mathcal{M}}_1 = \{\mu, P, Q, P.Q\} \cup \underline{\mathcal{E}} \text{ and } \underline{\mathcal{E}}_1 = \underline{\mathcal{E}} \setminus \{A\}. $$

The second pair deals with the between-subblocks stratum by pretending that block effects, but not subblock effects, are fixed, and declaring that the only effect to be estimated in this stratum is the main effect A. Thus

$$ \underline{\mathcal{M}}_2 = \{\mu, P\} \cup \underline{\mathcal{E}} \text{ and } \underline{\mathcal{E}}_2 = \{A\}. $$

## 2.2. Factorial terms and model specifications

Let $F_1, \dots, F_h$ denote all the genuine block and treatment factors involved in the experiment, that is, those which have a direct meaning for the experimenter. Denote by T the set of all n combinations of levels of $F_1, \dots, F_h$, with $n = n_1 \cdots n_h$ where $n_i$ is the number of levels of $F_i$.

The column vector $\tau$ of block and treatment effects belongs to the vector space $\mathbb{R}^n$. In the standard analysis of variance (ANOVA), this vector space is decomposed into mutually orthogonal subspaces $W_i$ associated with the $2^h$ subsets of factors $\{F_i : i \in I\}$, for $I \subseteq \{1, \dots, h\}$. These subspaces are given by the recurrence relation

$$ W_i = \nu_i \cap \left( \bigoplus_{J \in I} W_J \right)^{\perp}, $$

where $\nu_0$ is the subspace spanned by the all-one vector and $\nu_i$ is the subspace spanned by the indicator vectors of the level-combinations of all factors in $\{F_i : i \in I\}$ (Bailey, 2008). According to this decomposition, the effects can be decomposed into factorial effects, as given by the equation

$$ \tau = \sum_{I \subseteq \{1, \dots, h\}} S_I \tau, $$

where $S_I \tau$ is the orthogonal projection of $\tau$ onto $W_i$. By convention, the subsets $\{F_i : i \in I\}$ are called factorial terms and denoted by $\mathcal{F}(I)$. In examples, as in Section 2.1, we follow a more usual practice for writing factorial terms, so that $F_1 F_2 = \mathcal{F}(\{1, 2\})$ and $\mu = \mathcal{F}(\emptyset)$. The order of a factorial term $\mathcal{F}(I)$ is given by the cardinality of the subset I. Factorial terms of order 1 are called main effects, and factorial terms of order 2 or more are called interactions.
When constructing a fractional design it is necessary to assume that some factorial terms $\mathcal{F}(I)$ are negligible, that is, that $S_I \tau$ is zero. In the examples above, the model set $\mathcal{M}$ contains the non-negligible effects and its subset $\mathcal{E}$ contains the non-negligible effects that the experimenter wants to estimate. In multi-stratum experiments, such as Example 3, it is necessary to consider several such pairs of model and estimate sets, one pair for each stratum in which any fixed effect should be estimated.

## 2.3. Ingredients of the search

In this paper, generating a factorial design means specifying the combination of levels of the factors that must be allocated to each experimental unit. The situations that we consider generalise the examples of Section 2.1, using the generic factorial decomposition of Section 2.2. Their components consist of

(a) the list of genuine factors $F_1, \dots, F_h$, together with their numbers of levels and any hierarchy constraints;

(b) one or more joint model and estimate specifications ($\mathcal{M}, \mathcal{E}$), where $\mathcal{M}$ contains the factorial terms in the model and $\mathcal{E}$ contains the terms to estimate ($\mathcal{E} \subseteq \mathcal{M}$);

(c) the size of the experiment, i.e. the number $N$ of experimental units.

Note that we allow for two ways to handle the case when a factor $F_j$ at $n_i \times n_j$ levels is nested in a factor $F_i$ at $n_i$ levels: either $F_j$ is declared as a $n_j$-level factor and its actual levels are the levels of $F_i \wedge F_j$, or $F_j$ is declared as a $n_1 n_2 \dots n_{j-1}$-level factor and the hierarchy constraint $F_j \preceq F_i$ must be specified. We advise using the first option when $F_j$ is a natural refinement of $F_i$ and the second option when the hierarchy is due to experimental constraints rather than innate relationships between factors. Thus, in Example 3, we used the first option for the hierarchy between blocks and subblocks, with factors $P$ and $Q$ in the roles of $F_i$ and $F_j$ respectively. In contrast, we used the second option for the hierarchy between $A$ and the block factors (with factors $A$ and $P \wedge Q$ in the roles of $F_i$ and $F_j$ respectively), because these factors have no relationship except for experimental constraints.

In regular factorial designs, the number of units $N$ is constrained by the other specifications. For a solution to exist, $N$ must be a multiple of $p_1^{\rho_1} \cdots p_l^{\rho_l}$, where $p_1, \dots, p_l$ are the prime numbers that divide $n$ and $\rho_1, \dots, \rho_l$ are lower bounds on the exponents which depend on the model and estimate specifications. The algorithm described in this paper assumes that $N$ is given by the user (or by a higher-level algorithm). So it is up to the user to propose for $N$ a value of the form $N = q \cdot p_1^{\rho_1} \cdots p_l^{\rho_l}$, where $q$ is coprime to $p_1, \dots, p_l$ and $\rho_1 \le r_1, \dots, \rho_l \le r_l$. In practice, the user may proceed by trial and error by testing different values of $N$, provided the computing time is not too long. In Examples 2 and 3, the number of units is imposed by the problem, with $N = 2^2 \cdot 3 = 12$ and $N = 2^5 = 32$, respectively. In Example 1, the complete factorial design has size $2^5 \cdot 3^2 = 288$ but we look for smaller design sizes. We know that $\rho_1 \le 5$ and $\rho_2 \le 2$ and so we can proceed by trial and error to find them.

# 3. Elementary regular factorial designs

In Sections 3–5, the number of units and the numbers of levels of all factors are powers of the same prime $p$. Thus $N = p^r$ and $n = p^s$ for some scalars $r$ and $s$. This occurs in Example 3 with $p=2$, but not in Examples 1 and 2, which both involve primes 2 and 3. The designs are constructed using elementary abelian groups of exponent $p$, and so they are called *elementary designs* when they have to be distinguished from those considered in Section 6 and later.

The cyclic group of order $p$, denoted $C_p$, is identified with the integers modulo $p$ under addition. The experimental units are identified with the elements of a product group $U \cong (C_p)^r$, and the combinations of factor levels with the elements of a product group $T \cong (C_p)^s$. A design $d$ is a function from $U$ to $T$, allocating combination $d(u)$ to unit $u$. From now on, elements of both $U$ and $T$ are regarded as column vectors and $v^\top$ denotes the transpose of a vector $v$.

## 3.1. Pseudofactors

The canonical projections $V_1, \dots, V_r$ from $U$ onto the cyclic group $C_p$ are called the *unit pseudofactors*. If $u = (u_1, \dots, u_r)^\top$ is a unit in $U$ then $V_i(u) = u_i$. The canonical projections $A_1, \dots, A_s$ from $T$ onto the cyclic group $C_p$ are called the *treatment pseudofactors*, even though some of the genuine factors involved may be block factors. If $t = (t_1, \dots, t_s)^\top$ is in $T$, then $A_j(t) = t_j$.

Both kinds of pseudofactor must be viewed as technical intermediates in the design construction. Each genuine factor considered in Section 2 is the product of one or more treatment pseudofactors. If $F_i$ is the product of a family $(A_j)_{j \in J}$ of pseudofactors, this family is said to be a *decomposition of $F_i$ into pseudofactors*. For instance, if $F_i = A_1 \wedge A_2$, then $F_i$ is said to be decomposed into two pseudofactors $A_1, A_2$, which means that $F_i(t) = (A_1(t), A_2(t))^\top$ for every $t \in T$. The set of pseudofactors in the decomposition of $F_i$ is denoted by $\mathcal{P}(F_i)$. Alternatively, we use $\mathcal{P}(i)$ to denote the indices of the pseudofactors in $\mathcal{P}(F_i)$.

The decomposition of factors into pseudofactors induces a decomposition of factorial terms into pseudofactorial terms. Thus, the factorial term $\prod_{i \in I} F_i$ is decomposed into the pseudofactorial terms $\prod_{j \in J} A_j$ such that $J \cap \mathcal{P}(i) \neq \emptyset$ if and only if $i \in I$. We use the notation $\mathcal{P}(\cdot)$ to denote the set of pseudofactorial terms that decompose a given factorial term.
**Example 3 (Continued).** There are 32 = 2⁵ units and thus five unit pseudofactors $V_1, \dots, V_5$ at two levels. The nine treatment pseudofactors ($A_1, \dots, A_9$) will rather be denoted by $P_1, P_2, Q, U_1, U_2, A, B, C$ and $D$ to keep the correspondence with the genuine factors more explicit. With this notation, we have $\mathcal{P}(P) = \{P_1, P_2\}$, $\mathcal{P}(Q) = \{Q\}$, $\mathcal{P}(U) = \{U_1, U_2\}$, $\mathcal{P}(A) = \{A\}, \dots$ and $\mathcal{P}(D) = \{D\}$.

For the decomposition of factorial terms, we give the following examples, which include two main effects and three interactions:

$$
\begin{align*}
\tilde{\mathcal{P}}(A) &= \{A\}, \\
\tilde{\mathcal{P}}(P) &= \{P_1, P_2, P_1.P_2\}, \\
\tilde{\mathcal{P}}(A.B) &= \{A.B\}, \\
\tilde{\mathcal{P}}(P.Q) &= \{P_1.Q, P_2.Q, P_1.P_2.Q\}, \quad \text{and} \\
\tilde{\mathcal{P}}(P.U) &= \{P_1.U_1, P_1.U_2, P_1.U_1.U_2, P_2.U_1, P_2.U_2, P_2.U_1.U_2, P_1.P_2.U_1, P_1.P_2.U_2, P_1.P_2.U_1.U_2\}.
\end{align*}
$$

3.2. Characters

Our methods are based on the theory of duals of abelian groups, which can be found in Ledermann (1977). The group homomorphisms from $T$ into $C_p$ are all the linear combinations $A = a_1A_1 + \cdots + a_sA_s$ of the treatment pseudofactors, with $a_j \in C_p$. These homomorphisms are called the characters of $T$ and they make up a group $T^*$ called the dual of $T$. Each character can be represented by its vector of coefficients $\mathbf{a} = (a_1, \dots, a_s)^\top$ in the product group $(C_p)^s$, and the group $T^*$ can be consequently identified with this product group. The characters and dual of $U$ are defined and represented similarly. The elements of $T^*$ will be called treatment characters and the elements of $U^*$ unit characters.

Each character $A$ of $T^*$ is associated with a pseudofactorial effect $e_\tau(A)$. This means a precise linear combination of treatment effects in $\mathbb{R}^n$ or $\mathbb{C}^n$, either the general mean of $\tau$ if $A = 0$ or a contrast if $A \neq 0$ (see Kobilinsky, 1985 or Pistone and Rogantin, 2008 for more details). The important point is that each pseudofactorial effect belongs to a unique pseudofactorial term in the ANOVA decomposition of the treatment effects, and this term is easy to identify by the non-zero coefficients of $A$. For example, if there is only one non-zero coefficient $a_j$, then $A = a_jA_j$ and $e_\tau(A)$ belongs to the main effect of pseudofactor $A_j$. If $a_j \neq 0$ and $a_k \neq 0$ then the effect $e_\tau(a_jA_j + a_kA_k)$ belongs to the interaction between $A_j$ and $A_k$, and so on. When $p = 2$, each pseudofactorial term has one degree of freedom, so each includes a single character and a single pseudofactorial effect. In the general case, a pseudofactorial term of order $q$ includes $(p-1)^q$ characters and the same number of pseudofactorial effects.

Throughout this paper, additive notation is used for characters and their associated pseudofactorial effects. Multiplicative notation, in which $A_1^{a_1} \cdots A_s^{a_s}$ is used instead of $a_1A_1 + \cdots + a_sA_s$, is more common. Then $A_j^{\delta_j}$ belongs to the main effect of pseudofactor $A_j$, while $A_j^{\delta_j} A_k^{\delta_k}$ belongs to the interaction $A_j.A_k$, and so on. This has the disadvantage that $A_1A_2$ might be interpreted as an interaction or as one of the characters whose effect is part of that interaction.

3.3. Link between factorial terms and characters

Consider now the genuine factors $F_i$, for $i = 1, \dots, h$. Let $E_i$ be the subset of $T^*$ consisting of all characters $\sum a_jA_j$ involving pseudofactors in $\mathcal{P}(F_i)$ only. Let $\tilde{E}_i = E_i \setminus \{0\}$, so that $\tilde{E}_i$ consists of all the non-zero elements of $E_i$. By extension, put $E_\emptyset = \{0\}$ and, for each main effect or interaction $\mathcal{F}(I)$, put

$$
E_I = \bigoplus_{i \in I} E_i \quad \text{and} \quad \tilde{E}_I = \bigoplus_{i \in I} \tilde{E}_i,
$$

where $E \oplus E' = \{A+B: A \in E, B \in E'\}$. The first set $E_I$ includes all the characters associated with pseudofactors coming from the decomposition of the factors $F_i$, for $i \in I$. Those among them which have at least one non-zero coefficient for each factor make up the set $\tilde{E}_I$, which is therefore the subset of $T^*$ associated with the factorial term $\mathcal{F}(I)$.

**Example 3 (Continued).** For brevity, we give only one example based on the interaction P.Q. The corresponding set of characters $E_I$ is

$$
\{0, P_1, P_2, P_1 + P_2\} \oplus \{0, Q\} = \{0, P_1, P_2, P_1 + P_2, Q, P_1 + Q, P_2 + Q, P_1 + P_2 + Q\},
$$

while $\tilde{E}_I$ is restricted to

$$
\{P_1, P_2, P_1 + P_2\} \oplus \{Q\} = \{P_1 + Q, P_2 + Q, P_1 + P_2 + Q\}.
$$

3.4. Elementary regular design and its key matrix

In a regular factorial design, the treatment pseudofactors are algebraically derived from the unit ones.
**Definition 3.1** (Regular Design). If $U \cong (C_p)^r$ and $T \cong (C_p)^s$, a design $d$ is $C_p$-regular if there are coefficients $\phi_{ij}$ and $\alpha_j$ in $C_p$ such that, for $j = 1, \dots, s,$

$$A_j \circ d = \phi_{1j}V_1 + \dots + \phi_{ij}V_i + \dots + \phi_{rj}V_r + \alpha_j, \quad (1)$$

where $V_1, \dots, V_r$ are the unit pseudofactors and $A_1, \dots, A_s$ the treatment pseudofactors.

In a $C_p$-regular factorial design, the combination $t$ allocated to unit $u$ satisfies $t = \Phi^\top u + t_0$, where

$$\Phi = \begin{pmatrix} \phi_{11} & \cdots & \phi_{1s} \\ \vdots & & \vdots \\ \phi_{r1} & \cdots & \phi_{rs} \end{pmatrix}$$

and $t_0 = (\alpha_1, \dots, \alpha_s)^\top$. The definition below follows the definition of the design key $K$ given by Patterson (1976) and Bailey et al. (1977). As Cheng (2014) does, we have adapted their definition so that the subscripts on the entries in $\Phi$ appear in the usual order; note that Kobilinsky and Monod (1991) use its transpose.

**Definition 3.2** (Key Matrix). Let $d$ be a regular design as specified in Definition 3.1. The matrix $\Phi$ is called the key matrix of $d$.

Define the mapping $\psi: U \to T$ by $\psi(\mathbf{u}) = d(\mathbf{u}) - t_0 = \Phi^\top \mathbf{u}$. Then $\psi$ is a group homomorphism, in the sense that if $\mathbf{u}$ and $\mathbf{u}'$ are in $U$ then $\psi(\mathbf{u} + \mathbf{u}') = \psi(\mathbf{u}) + \psi(\mathbf{u}')$. The dual of the homomorphism $\psi$, denoted by $\varphi$, is the homomorphism from $T^*$ to $U^*$ sending a character $A$ of $T$ to the character $\varphi(A)$ of $U$ defined by

$$(\varphi(A))(\mathbf{u}) = A(\psi(\mathbf{u})) = A(\Phi^\top \mathbf{u}).$$

It is clear from (1) that if $A = a_1A_1 + \dots + a_sA_s$ then $\varphi(A) = B = b_1V_1 + \dots + b_rV_r$, where

$$b_i = \sum_j \phi_{ij} a_j \quad \text{for } i = 1, \dots, r.$$

Representing the characters $A$ and $B$ by their column vectors $\mathbf{a}$ and $\mathbf{b}$ gives $\mathbf{b} = \Phi \mathbf{a}$. The $j$th column $(\phi_{1j}, \dots, \phi_{rj})^\top$ of $\Phi$ is the unit character $\varphi(A_j)$ in $U^*$ which is the image of the treatment character $A_j$ in $T^*$. We denote it by $\tilde{A}_j$.

### 3.5. Confounding

For a vector $\mathbf{c}$ in $\mathbb{R}^n$ indexed by the genuine combinations, let $\mathbf{c}^{(d)}$ denote the vector in $\mathbb{R}^N$ indexed by the units and defined by $(\mathbf{c}^{(d)})_{\mathbf{u}} = (\mathbf{c})_{d(\mathbf{u})}$. Two treatment effects $\mathbf{c}_1^\top \tau$ and $\mathbf{c}_2^\top \tau$ are said to be confounded with each other in design $d$ if there is a constant $\gamma$ such that $\mathbf{c}_1^{(d)} = \gamma \mathbf{c}_2^{(d)}$. In this case, it is impossible to estimate the effects $\mathbf{c}_1^\top \tau$ and $\mathbf{c}_2^\top \tau$ separately. Treatment effects $\mathbf{c}_1^\top \tau$ and $\mathbf{c}_2^\top \tau$ are said to be orthogonal to each other if $\mathbf{c}_1^{(d)} \top \mathbf{c}_2^{(d)} = 0$.

The basic statistical properties of regular factorial designs are given by the following proposition (see e.g. Kobilinsky and Monod, 1995). Here $\text{Ker}(\varphi)$ denotes the kernel of $\varphi$, which is $\{A \in T^* : \varphi(A) = 0\}$.

**Proposition 3.1.** Let $A$ and $B$ denote two characters in $T^*$. A regular design with key matrix $\Phi$ and corresponding homomorphism $\varphi$ satisfies the following four properties:

(i) The pseudofactorial effect $e_\tau(A)$ is confounded with the general mean $e_\tau(0)$ if and only if $A \in \text{Ker}(\varphi)$;

(ii) The pseudofactorial effects $e_\tau(A)$ and $e_\tau(B)$ are confounded with each other if and only if $A - B \in \text{Ker}(\varphi)$;

(iii) The sets of mutually confounded pseudofactorial effects are given by the cosets of the subgroup $\text{Ker}(\varphi)$;

(iv) The pseudofactorial effects $e_\tau(A)$ and $e_\tau(B)$ are orthogonal to each other if $A$ and $B$ are in different cosets of $\text{Ker}(\varphi)$.

If $A \in T^*$, the effect $e_\tau(A)$, is estimable if and only if it is not confounded with any other non-zero effect. Proposition 3.1 shows that this occurs if and only if all the other characters $B$ in the same coset of $\text{Ker}(\varphi)$ as $A$ are assumed to have no effect on the expectation of the response to be measured ($e_\tau(B) = 0$).

## 4. Conditions on the design key matrix

### 4.1. Ineligible characters and factorial terms

#### 4.1.1. Ineligible characters due to model specifications

For the experimenter, the model and the effects to estimate consist of factorial terms defined on the genuine factors, as illustrated by the examples in Section 2. As shown in Section 3.3, each such factorial term $\mathcal{F}(I)$ is associated with a set $\tilde{E}_I$ of characters in $T^*$. Thus, for each pair of model-estimate sets $(\mathcal{M}, \mathcal{E})$ of factorial terms, there is an associated pair $(\mathcal{M}, \mathcal{E})$ of character sets. The model set $\mathcal{M}$ is the union of the sets $\tilde{E}_I$, for all factorial terms $\mathcal{F}(I)$ in $\mathcal{M}$, while the estimate set $\mathcal{E}$ is the union of the sets $\tilde{E}_I$ for all terms $\mathcal{F}(I)$ in $\mathcal{E}$. If $E$ and $E'$ are any two subsets of $T^*$, we write $E - E' = \{A - B : A \in E, B \in E'\}$.
<table><caption>**Table 1**
Construction of the ineligible set $\underline{I}$ by symmetric differences (Example 4).</caption><thead><tr><th rowspan="2">Estimate set <i>&epsilon;</i></th><th colspan="6">Model set <i>M</i></th></tr><tr><th>&mu;</th><th>A</th><th>B</th><th>C</th><th>A.B</th><th>B.C</th></tr></thead><tbody><tr><td>A</td><td>A</td><td>&mu;</td><td>A.B</td><td>A.C</td><td>B</td><td>A.B.C</td></tr><tr><td>B</td><td>B</td><td>A.B</td><td>&mu;</td><td>B.C</td><td>A</td><td>C</td></tr><tr><td>C</td><td>C</td><td>A.C</td><td>B.C</td><td>&mu;</td><td>A.B.C</td><td>B</td></tr><tr><td>A.B</td><td>A.B</td><td>B</td><td>A</td><td>A.B.C</td><td>&mu;</td><td>A.C</td></tr></tbody></table>

**Definition 4.1** (Ineligible Characters). Let $(\mathcal{M}, \mathcal{E})$ be model-estimate sets of factorial terms and $(\mathcal{M}, \mathcal{E})$ be the associated model-estimate sets of characters. Put

$$ \underline{I} = \{A - B : A \in \mathcal{E}, B \in \mathcal{M}, A \neq B\} = (\mathcal{E} - \mathcal{M}) \setminus \{0\}. $$

Then $\underline{I}$ is called the set of ineligible characters with respect to $(\mathcal{M}, \mathcal{E})$ or, equivalently, to $(\mathcal{M}, \mathcal{E})$.

**Proposition 4.1.** Consider a regular design with key matrix $\Phi$ and corresponding homomorphism $\varphi$, and the model consisting of the factorial terms in $\mathcal{M}$, considered as fixed effects. All factorial terms in $\mathcal{E}$ are estimable if and only if $\text{Ker}(\varphi) \cap \underline{I} = \emptyset$, where $\underline{I}$ is the set of ineligible characters with respect to $(\mathcal{M}, \mathcal{E})$.

All proofs are in the Supplementary material (see Appendix A).

### 4.1.2. Ineligible factorial terms

Not every set of characters of $T$ is a union of the subsets $\tilde{E}_i$ associated with factorial terms $\mathcal{F}(I)$. However, Proposition 4.2 shows that, for any given model-estimate pair $(\mathcal{M}, \mathcal{E})$, there is a set $\underline{I}$ of factorial terms such that

$$ \underline{I} = \bigcup_{\mathcal{F}(K) \in \underline{I}} \tilde{E}_K. \qquad (2) $$

Thus it makes sense to say that $\mathcal{F}(K)$ is an ineligible factorial term if $\mathcal{F}(K) \in \underline{I}$. Moreover, $\underline{I}$ can be calculated explicitly, as the union of the sets $\underline{I}(I, J) \setminus \{\mathcal{F}(\emptyset)\}$ for $\mathcal{F}(I) \in \mathcal{E}$ and $\mathcal{F}(J) \in \mathcal{M}$, where $\underline{I}(I, J) = \{\mathcal{F}(K) : I \Delta J \subseteq K \subseteq (I \Delta J) \cup L\}$ and $L = \{i \in I \cap J : n_i > 2\}$, so that $\underline{I}(I, J) = \{\mathcal{F}(K) : I \Delta J \subseteq K \subseteq I \cup J\}$ if $p \ge 3$. Here $I \Delta J$ denotes the symmetric difference between $I$ and $J$, which is $(I \setminus J) \cup (J \setminus I)$. Proposition 4.3 shows that, in three important special cases,

$$ \underline{I} = \{\mathcal{F}(I \Delta J) : \mathcal{F}(I) \in \mathcal{E}, \mathcal{F}(J) \in \mathcal{M}, I \neq J\}. \qquad (3) $$

**Proposition 4.2.** Let $I$ and $J$ be subsets of $\{1, \dots, h\}$. Then

$$ \tilde{E}_I - \tilde{E}_J = \bigcup_{\mathcal{F}(K) \in \underline{I}(I, J)} \tilde{E}_K. $$

The model $\mathcal{M}$ is usually complete in the sense that, when it contains an effect, it also contains all effects marginal to it. For instance, if it contains the interaction $F \cdot G$, it also contains the main effects $F$ and $G$ and the general mean $\mu$.

**Proposition 4.3.** If the model $\mathcal{M}$ is complete, or if the estimate-set $\mathcal{E}$ is complete, or if all genuine factors have two levels, then the set $\underline{I}$ of ineligible factorial terms is given by Eq. (3).

Thus the first step to determine the set $\underline{I}$ of ineligible characters is to determine the set $\underline{I}$ of ineligible factorial terms. If any of the conditions in Proposition 4.3 is satisfied, this is done by identifying the sets $K = I \Delta J$ for $\mathcal{F}(I)$ in $\mathcal{E}$ and $\mathcal{F}(J)$ in $\mathcal{M}$. Otherwise, the sets $\underline{I}(I, J)$ of ineligible factorial terms must be used, for $\mathcal{F}(I)$ in $\mathcal{E}$ and $\mathcal{F}(J)$ in $\mathcal{M}$. This step can be performed before the decomposition into pseudofactors. In a second step, the ineligible characters can be deduced from Eq. (2).

**Example 4.** For sake of brevity, we give an example simpler than Example 3. Suppose that there are three factors $A, B, C$, that

$$ \mathcal{M} = \{\mu, A, B, C, A.B, B.C\} \quad \text{and} \quad \mathcal{E} = \{A, B, C, A.B\}. $$

Table 1 gives for each model term $\mathcal{F}(I)$ and each estimate term $\mathcal{F}(J)$ the associated ineligible factorial term $\mathcal{F}(I \Delta J)$. Therefore $\underline{I}$ includes all non-mean terms in this table, that is $\{A, B, C, A.B, A.C, B.C, A.B.C\}$.

**Example 5.** Suppose that there are two factors $A, B$ and we choose $\mathcal{M} = \{\mu, B, A.B\}$ and $\mathcal{E} = \{B\}$. Because it does not contain $A$, the model $\mathcal{M}$ is not complete, and, because it does not contain $\mu$, $\mathcal{E}$ is not complete. Proposition 4.2 implies that $\underline{I} = \{A, B\}$ if factor $B$ has two levels, whereas $\underline{I} = \{A, B, A.B\}$ if factor $B$ has three or more levels. Here are some possibilities.
(a) If factors A and B have two levels, then $\mathcal{M} = \{0, B, A + B\}$, $\mathcal{E} = \{B\}$ and $\mathcal{I} = \{A, B\}$. So the only ineligible factorial terms are the main effects of A and B. We could construct a design by confounding the character $A + B$.

(b) If factor A has four levels and factor B has two levels, then $\mathcal{M} = \{0, B, A_1 + B, A_2 + B, A_1 + A_2 + B\}$, $\mathcal{E} = \{B\}$ and $\mathcal{I} = \{B, A_1, A_2, A_1 + A_2\}$. We could construct a design confounding any one of $A_1 + B, A_2 + B, A_1 + A_2 + B$.

(c) If factors A and B have three levels, then $\mathcal{M} = \{0, B, 2B, A + B, A + 2B, 2A + B, 2A + 2B\}$, $\mathcal{E} = \{B, 2B\}$ and $\mathcal{I} = \{A, 2A, B, 2B, A + B, A + 2B, 2A + B, 2A + 2B\}$. So the ineligible factorial terms are the main effects of A and B and the interaction A.B. If we confound either part of the A.B interaction then B is confounded with the other part and so cannot be estimated.

(d) If factor A has two levels and factor B has four levels, then $\mathcal{M} = \{0, B_1, B_2, B_1 + B_2, A + B_1, A + B_2, A + B_1 + B_2\}$, $\mathcal{E} = \{B_1, B_2, B_1 + B_2\}$ and $\mathcal{I} = \{A, B_1, B_2, B_1 + B_2, A + B_1, A + B_2, A + B_1 + B_2\}$. Again, no non-zero character is eligible for confounding.

As this example shows, if the model and the estimate-set are both incomplete, the ineligible factorial terms may depend on the numbers of levels of the factors.

## 4.2. Hierarchy constraints

In practice, besides constraints of ineligibility, it may be necessary to satisfy hierarchy constraints between factors, such as those shown in Examples 2 and 3. It is assumed that all constraints are of the form $F_{i_1} \wedge \dots \wedge F_{i_k} \prec F_{i_0}$ or can be deduced from such constraints. This assumption is satisfied in most practical cases. Recall that $\mathcal{P}(i)$ denotes the pseudofactors $A_j$ that decompose the genuine factor $F_i$, or more precisely the set of indices of these pseudofactors.

**Proposition 4.4.** Let $F_{i_0}, F_{i_1}, \dots, F_{i_k}$ be $k+1$ genuine factors. Let $J = \mathcal{P}(i_1) \cup \dots \cup \mathcal{P}(i_k)$, which is the set of pseudofactors that decompose $F_{i_1}, \dots, F_{i_k}$. For a regular design with design key matrix $\Phi$, the following three conditions are equivalent:

(i) $F_{i_1} \wedge \dots \wedge F_{i_k} \prec F_{i_0}$;

(ii) $\bigwedge_{j \in J} A_j \prec A_{j_0}$, for all $j_0 \in \mathcal{P}(i_0)$;

(iii) if $j_0 \in \mathcal{P}(i_0)$ then the column $\tilde{A}_{j_0}$ of $\Phi$ is a linear combination of the columns $\tilde{A}_j$, for $j \in J$.

**Proposition 4.4** shows that each hierarchy constraint between genuine factors generates one or more hierarchy constraints ($A_{j_1} \wedge \dots \wedge A_{j_l} \prec A_{j_0}$) between pseudofactors, each of which must be satisfied during the design search. In the sequel, it will be assumed that the pseudofactors are ordered so that, for any such constraint, $j_0$ is greater than $j_1, \dots, j_l$. The set of all coarser pseudofactors $A_{j_0}$ involved in such constraints, denoted by $\mathcal{H}_+$, results from the decomposition of the coarser factors in the original constraints. The set of all pseudofactors in the finer part of the constraint ($A_{j_1} \wedge \dots \wedge A_{j_l} \prec A_{j_0}$) will be denoted by $\mathcal{H}_{<z}$, where, depending on the context, $z$ may refer to the coarser pseudofactor $A_{j_0}$ or to its index $j_0$.

**Example 3 (Continued).** Recall that there is a unique hierarchy constraint $P \wedge Q \prec A$. It follows that $\mathcal{H}_+ = \{A\}$ and $\mathcal{H}_{<A} = \{P, Q\}$.

## 4.3. Extension to multi-stratum experiments

Multi-stratum experiments provide many examples of mixed models, where some factorial terms are assumed to be random and centred rather than fixed. For example, the block factors in Examples 2 and 3 are assumed to be random while the treatment factors are assumed fixed. If the proper randomisation is applied, the block effects can indeed be considered random and centred.

When the block factors form a poset block structure on the experimental units, each stratum is defined by a suitable subset of the block factors. In Example 3 these subsets are $\{P\}, \{P, Q\}$ and $\{P, Q, U\}$. In Example 2 they are $\{R\}, \{C\}, \{R, C\}$ and $\{R, C, Z\}$, where Z is an implicit two-level factor for the units within each cell.

Let J be the set of indices for such a stratum subset. Then the stratum can also be identified to the subset $\mathcal{G}(J)$ of random block factorial terms $\mathcal{F}(K)$, for $K \subseteq J$ and $K \not\subseteq J'$ for all stratum subsets $J', J' \subseteq J$. In Example 3 these subsets of factorial terms are $\{P\}, \{Q, P.Q\}$ and $\{U, P.U, Q.U, P.Q.U\}$. In Example 2 they are $\{R\}, \{C\}, \{R.C\}$ and $\{Z, R.Z, C.Z, R.C.Z\}$. Note that the factorial terms which include the units factors U or Z correspond to residual error terms and are not considered below.

If the treatment factorial term $\mathcal{F}(I)$ is subject to the hierarchy constraint $\bigwedge_{j \in J} F_j \preceq \bigwedge_{i \in I} F_i$ then $\mathcal{F}(I)$ should, if possible, be estimated in the stratum defined by J. To do this, we define a model-set $\mathcal{M}_J$ and estimate-set $\mathcal{E}_J$ associated with this stratum, as follows. The model-set $\mathcal{M}_J$ is made by removing from $\mathcal{M}$ all the block factorial terms in $\mathcal{G}(J')$ for any stratum subset $J'$ with $J' \supseteq J$. The estimate-set $\mathcal{E}_J$ contains all treatment factorial terms $\mathcal{F}(I)$ in $\mathcal{E}$ which are subject to the hierarchy constraint $\bigwedge_{j \in J} F_j \preceq \bigwedge_{i \in I} F_i$ but are not subject to any constraint $\bigwedge_{j \in K} F_j \preceq \bigwedge_{i \in I} F_i$ for which K is a subset of the (indices of the) block factors defining a stratum and $K \subseteq J$.

For each stratum, the set $\mathcal{I}_J$ is derived from $\mathcal{E}_J$ and $\mathcal{M}_J$ by Eq. (3) if $\mathcal{M}_J$ is complete or if $\mathcal{E}_J$ is complete or if all factors involved in $\mathcal{M}_J$ have two levels, and otherwise is derived from $\mathcal{E}_J$ and $\mathcal{M}_J$ by using Proposition 4.2. The union $\mathcal{I}$ of these sets $\mathcal{I}_J$ is the full set of ineligible terms, from which the full set $\mathcal{I}$ of ineligible characters must be deduced.
**Example 3 (Continued).** For simplicity, we use subscripts 1 and 2 instead of {$P, Q, U$} and {$P, Q$}. We have

$$ 
\begin{align*}
\mathcal{I}_1 &= \{A, B, C, D, A.B, A.C, A.D, B.C, B.D, C.D, \\
&\qquad A.B.C, A.B.D, A.C.D, B.C.D, A.B.C.D, \\
&\qquad B.P, C.P, D.P, A.B.P, A.C.P, A.D.P, B.C.P, B.D.P, C.D.P, \\
&\qquad B.Q, C.Q, D.Q, A.B.Q, A.C.Q, A.D.Q, B.C.Q, B.D.Q, C.D.Q, \\
&\qquad B.P.Q, C.P.Q, D.P.Q, A.B.P.Q, A.C.P.Q, A.D.P.Q, B.C.P.Q, B.D.P.Q, C.D.P.Q\}
\end{align*}
$$

and

$$ \mathcal{I}_2 = \{A, B, C, D, A.B, A.C, A.D, A.B.C, A.B.D, A.C.D, A.P\}. $$

The set of model-based ineligible factorial terms is given by $\mathcal{I}_{1,2} = \mathcal{I}_1 \cup \mathcal{I}_2 = \mathcal{I}_1 \cup \{A, A.P\}$. Note that the interactions A.Q and A.P.Q are absent from $\mathcal{I}_{1,2}$, whereas they would be included among the ineligible terms if the latter were based on $\mathcal{M}$ and $\mathcal{E}$.

### 4.4. Ineligibility due to combinatorial requirements

In some situations, it is required that all the combinations of levels of some factors are in the experiment, whatever the model and estimate constraints are. Then all the terms including the corresponding pseudofactors have to be included in the ineligible set. For instance, in a row-and-column design, all combinations of levels of the row and column factors must be present. If the rows are defined by a factor $R$, the columns by a factor $C$, then $\mathcal{I}$ must include the factorial terms $R, C$ and $R.C$ whatever the model specifications, and $\mathcal{I}$ will include the corresponding characters.

**Example 3 (Continued).** The experimental units consist of all levels combinations of $P, Q$ and $U$. To ensure that the design has all these combinations, the set of ineligible terms must become $\mathcal{I} = \mathcal{I}_{1,2} \cup \{P, Q, P.Q, U, P.U, Q.U, P.Q.U\}$.

## 5. Search for key matrices of elementary designs

### 5.1. Main steps

Section 4 showed that the first step in the search for a key matrix solution is the determination of the set $\mathcal{I}$ of ineligible factorial terms. Further technical details are given in Section 5.2. The next step is to deduce from $\mathcal{I}$ a reduced set $\mathcal{R}$ of representative ineligible treatment characters, as described in Section 5.3. The third step is the search for one or more key matrices $\Phi$ satisfying condition (iii) of Proposition 4.4 (hierarchies) and Eq. (5) of Definition 5.1 (ineligibility). This can be done by the backtrack algorithm described in Section 5.4.

### 5.2. Determining the ineligible factorial terms

When Eq. (3) holds, the elements of $\mathcal{I}$ can be identified by representing each main effect or interaction $\mathcal{F}(I)$ by a vector of dimension $h$ over $C_2$ with $i$th coordinate equal to 1 if $i \in \mathcal{I}$, to 0 otherwise. With this representation, if $\mathbf{x}$ represents $\mathcal{F}(I)$ and $\mathbf{z}$ represents $\mathcal{F}(J)$, then the vector representing $\mathcal{F}(I \Delta J)$ is simply $\mathbf{x} + \mathbf{z}$, whose $a$th coordinate is 1 if $x_a \neq z_a$ and is 0 otherwise.

**Example 4 (Continued).** The associated vectors $\mathbf{x}, \mathbf{z}$ and $\mathbf{x} + \mathbf{z}$ for $\mathbf{x} \neq \mathbf{z}$ are

* for $\mathcal{M}(\mathbf{x})$: $\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 0 \\ 0 \end{pmatrix}$;

* for $\mathcal{E}(\mathbf{z})$: $\begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 0 \\ 0 \end{pmatrix}$;

* for $\mathcal{I}(\mathbf{x} + \mathbf{z})$: $\begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 0 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix}$.

### 5.3. Reduced set of ineligible characters

Once the full set $\mathcal{I}$ of ineligible factorial terms has been determined, the ineligible characters can be deduced as the union of the sets $\tilde{E}_I$, for $I \in \mathcal{I}$. If $\mathcal{I}$ is the set of ineligible characters, the homomorphism $\varphi$ must satisfy

$$ \varphi(A) \neq 0 \quad \text{for every character } A \in \mathcal{I}. \qquad (4) $$

However, if $A$ and $B$ in $\mathcal{I}$ are such that $A$ is an integer multiple of $B$, that is $A = kB$, then the inequality $\varphi(A) \neq 0$ clearly implies $\varphi(B) \neq 0$. In the search for $\varphi$, the inequality (4) has therefore to be checked only for an adequately chosen subset $\mathcal{R}$ of $\mathcal{I}$ called a reduced ineligible set.
**Definition 5.1** (Reduced Set). A reduced ineligible set $\mathcal{R}$ is any subset of $\mathcal{I}$ such that condition (4) is equivalent to the apparently weaker condition

$$ \varphi(A) \neq 0 \quad \text{for every } A \in \mathcal{R}. \qquad (5) $$

We now indicate how such a reduced ineligible set can be selected. Let $\langle A \rangle$ be the cyclic subgroup generated by A. The relation

$$ \langle A \rangle = \langle B \rangle \iff \exists \delta, \delta' \in \mathbb{N} \text{ such that } A = \delta B \text{ and } B = \delta' A \qquad (6) $$

defines an equivalence relation on $\mathcal{I}$. Clearly it is enough to check inequality (4) for only one representative in each equivalence class.

When $T^*$ and $U^*$ are elementary abelian p-groups, all non-zero characters have order p, the non-zero equivalence classes contain $p-1$ characters, and the representatives are often chosen as those whose first non-zero coordinate is 1. In Section 5.4, the representatives are chosen like this if there is only one non-zero coordinate, but otherwise chosen to have last non-zero coordinate $-1$.

**Example 3** (Continued). In this example, $p = 2$ so that each equivalence class has a single element: hence the set of ineligible characters cannot be reduced.

**Example 6.** If factors A and B both have three levels and $\mathcal{I}$ contains their two-factor interaction then a possible reduced set $\mathcal{R}$ contains the characters $A - B (= A + 2B)$ and $2A - B (= 2A + 2B)$ for this term but neither $A + B$ nor $2A + B$.

### 5.4. Elementary backtrack search

Searching for the key matrix $\Phi$ is equivalent to searching for its columns $\tilde{A}_1, \dots, \tilde{A}_s$ among the set $U^*$ of unit characters. For the homomorphism $\varphi$ to satisfy (5), these characters must satisfy:

$$ \text{for every } a_1A_1 + \dots + a_sA_s \in \mathcal{R}, \quad a_1\tilde{A}_1 + \dots + a_s\tilde{A}_s \neq 0. \qquad (7) $$

In addition, because of the hierarchy constraints, some columns must be linear combinations of other ones of smaller index, as shown in Proposition 4.4. More precisely, for the indices $j$ in the subset $\mathcal{H}_+$, the columns $\tilde{A}_j$ must satisfy

$$ \tilde{A}_j = \sum_{k \in \mathcal{H}_{<j}} a_k \tilde{A}_k \qquad (8) $$

for some values $a_k$. Note that all indices $k$ in $\mathcal{H}_{<j}$ are strictly smaller than $j$.

The columns $\tilde{A}_j$ can be found successively by the backtrack search presented in Algorithm 1 below. The process may either end at the first success or continue until there is no more admissible $\tilde{A}_1$ to select.

Algorithm 1 involves non-trivial calculations only when determining the set $a\tilde{A}_j$ in Step 1. A unit character is considered admissible for column $j$ if it satisfies the inequalities (7) involving it and the previous columns $\tilde{A}_1, \dots, \tilde{A}_{j-1}$. Let $\mathcal{R}_j$ be the subset of characters $a_1A_1 + \dots + a_jA_j$ in $\mathcal{R}$ having $a_j$ as the last non-zero coordinate. The inequalities in (7) to consider when searching for $\tilde{A}_j$ are all those involving a character in $\mathcal{R}_j$. They can be written:

$$ \text{for every } a_1A_1 + \dots + a_{j-1}A_{j-1} - A_j \text{ in } \mathcal{R}_j, \quad \tilde{A}_j \neq a_1\tilde{A}_1 + \dots + a_{j-1}\tilde{A}_{j-1}, \qquad (9) $$

which includes $\tilde{A}_j \neq 0$ if $A_j \in \mathcal{R}$. The admissible characters satisfying (9) are looked for among $U^*$ unless $j$ belongs to $\mathcal{H}_+$. In that case, the hierarchy constraints (8) allow the search to be restricted to the subgroup of $U^*$ generated by the columns in $\mathcal{H}_{<j}$.

Optimising the backtrack algorithm is a complex task which is not the subject of this paper. In Appendix B of the Supplementary material (see Appendix A), however, we describe a few tricks implemented in planor (Kobilinsky, 2005) to make the backtrack search run faster.

### 5.5. End of the search

It is sometimes necessary to go beyond a success to find other solutions, for instance if the whole set of solutions is searched for, or if the backtrack column search is part of a backtrack process among Sylow components as will be described in Section 7.

When the number of factors involved increases, the time taken by the search may become very long, especially if there is no solution or only a small number of solutions. So it can be necessary to stop a search which takes too much time. In that case, which must clearly be distinguished from a true failure, it is useful to know the index of the last column reached as this indicates the kind of experimental design obtainable in a reasonable time.
**Algorithm 1 Backtrack search**

jprev ← 0 and j ← 1

while j > 0 do

{Step 1: determine or update the admissible set}

if jprev < j then {forward case}

determine the set aA<sub>j</sub> of currently admissible unit characters for column j of Φ

else {backward case}

delete the current character in column j from aA<sub>j</sub>

end if

{Step 2: determine and make the next move}

if aA<sub>j</sub> is empty then

j ← j − 1 {no solution for column j, so move backward}

else

set column j to the first element in aA<sub>j</sub>

if j < s then

j ← j + 1 {solution for column j, so move forward}

else {j = s so that all columns have been found}

save the current key matrix in the solution set

either stop or continue to find more solutions

end if

end if

jprev ← j

end while

**Example 3 (Continued).** For completeness, we give the first solution for Φ found by the planor algorithm:

$$
\Phi = \begin{array}{c}
\begin{array}{ccccccccc}
\tilde{p}_1 & \tilde{p}_2 & \tilde{q} & \tilde{U}_1 & \tilde{U}_2 & \tilde{A} & \tilde{B} & \tilde{C} & \tilde{D} \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 1
\end{array} \\
\begin{array}{c}
v_1 \\
v_2 \\
v_3 \\
v_4 \\
v_5
\end{array}
\end{array}
$$

In its current implementation, *planor* finds 9216 solutions, including solutions obtained from others by permutation of *B*, *C*, *D*, or *P*₁ and *P*₂, or *U*₁ and *U*₂.

It took planor 0.073 s to find the first solution and 2.193 s to complete an exhaustive search, using a laptop computer (7.7
Gio RAM, Intel Core i7-5600U CPU@2.60 GHz x 4 processor; Linux Ubuntu operating system; R version 3.3.0; planor version
1.0.1). See Appendix C.3 of the Supplementary material for more details (see Appendix A).

**6. Generalised regular factorial designs**

This section extends the previous methods to cover designs involving more than one prime. The underlying theory is
given by Bailey (1977, 1985), Kobilinsky (1985), Kobilinsky and Monod (1995) and Pistone and Rogantin (2008).

6.1. Pseudofactors

Let $p_1, \dots, p_l$ denote distinct prime numbers. The experimental units are now identified with the elements of a product group $U = (C_{p_1})^{r_1} \times \dots \times (C_{p_l})^{r_l}$, and the genuine combinations with the elements of a product group $T = (C_{p_1})^{s_1} \times \dots \times (C_{p_l})^{s_l}$. Note that the theory extends to prime powers but, for simplicity, we avoid this level of generality here. The number $N$ of units and the number $n$ of genuine combinations factorise into $N = N_1 \times \dots \times N_l$ and $n = n_1 \times \dots \times n_l$, where $N_k = p_k^{r_k}$ and $n_k = p_k^{s_k}$.

As in Section 3, the treatment pseudofactors associated with this decomposition of *T* are a finer decomposition of the
genuine factors *F*<sub>*i*</sub>.

Let $r = \sum_k r_k$ and $s = \sum_k s_k$. In what follows, we denote the *r* unit pseudofactors by $V_i$, for $i = 1, \dots, r$ and the *s* treatment pseudofactors by $A_j$, for $j = 1, \dots, s$. We denote by $\pi(V_i)$ and $\pi(A_j)$ the numbers of levels of pseudofactors $V_i$ and $A_j$. Occasionally, to stress the structure introduced by the different primes, we denote the unit pseudofactors at $p_k$ levels by $V_{[k,i]}$ and the treatment pseudofactors at $p_k$ levels by $A_{[k,j]}$, for $k = 1, \dots, l$, for $i = 1, \dots, r_k$ and for $j = 1, \dots, s_k$. In the examples, though, we use $p_k$ rather than $k$ as the first index between square brackets, because it improves clarity. Unless specified otherwise, the pseudofactors are ordered in the natural lexicographic order induced by their double index, which puts together all the pseudofactors associated with the same prime.
To define the pseudofactors properly, all levels are embedded into the same cyclic group $C_M$ where $M = \prod_{k=1}^l p_k$. Unit and treatment pseudofactors are considered as mappings into that cyclic group $C_M$. That is, if $\mathbf{u} = (u_1, \dots, u_r)^\top$ is a unit in $U$ and $\mathbf{t} = (t_1, \dots, t_s)^\top$ a genuine combination in $T$, then

$$V_i(\mathbf{u}) = \frac{M}{\pi(V_i)} u_i \quad \text{and} \quad A_j(\mathbf{t}) = \frac{M}{\pi(A_j)} t_j. \qquad (10)$$

**Example 1 (Continued).** The two primes involved in this example are 2 and 3, so that all genuine factors are decomposed into pseudofactors at two or three levels. We have $M = 6$, $s_2 = 5$ and $s_3 = 2$, and the pseudofactors are $A_{[2,1]}$, $A_{[2,2]}$, $A_{[2,3]}$, $A_{[2,4]}$, $A_{[2,5]}$, $A_{[3,1]}$ and $A_{[3,2]}$, where the first index denotes the prime and the second index runs from 1 to $s_2$ or $s_3$. The association between factors and pseudofactors is given by

$$F_1 = A_{[2,1]} \wedge A_{[3,1]}; \quad F_2 = A_{[2,2]} \wedge A_{[2,3]}; \quad F_3 = A_{[3,2]} \quad \text{and} \quad F_4 = A_{[2,4]} \wedge A_{[2,5]}.$$

An equivalent notation keeps the lexicographic order of the pseudofactors but identifies them with a single index. In the $\mathcal{P}(\cdot)$ notation, the association reads

$$\mathcal{P}(F_1) = \{A_1, A_6\}; \quad \mathcal{P}(F_2) = \{A_2, A_3\}; \quad \mathcal{P}(F_3) = \{A_7\} \quad \text{and} \quad \mathcal{P}(F_4) = \{A_4, A_5\}.$$

Consider, for example, the combination identified with $\mathbf{t} = (1, 0, 1, 0, 2, 1)^\top \in (C_2)^5 \times (C_3)^2$. Following Eq. (10), we have $A_{[2,1]}(\mathbf{t}) = 3$, $A_{[2,2]}(\mathbf{t}) = 0$, $A_{[2,3]}(\mathbf{t}) = 3$, $A_{[2,4]}(\mathbf{t}) = 3$, $A_{[2,5]}(\mathbf{t}) = 0$, $A_{[3,1]}(\mathbf{t}) = 4$ and $A_{[3,2]}(\mathbf{t}) = 2$.

An option involving prime powers is to decompose the factors as little as possible. Then the pseudofactors are $A'_{[2,1]}$, $A'_{[2,2]}$, $A'_{[2,3]}$, $A'_{[3,1]}$ and $A'_{[3,2]}$, at respectively 2, 4, 4, 3 and 3 levels, with $M = 12$, so that

$$F_1 = A'_{[2,1]} \wedge A'_{[3,1]}; \quad F_2 = A'_{[2,2]}; \quad F_3 = A'_{[3,2]} \quad \text{and} \quad F_4 = A'_{[2,3]}.$$

This is not implemented in planor but discussed briefly in Section 8.

**Example 2 (Continued).** The two primes involved are 2 and 3 once again. Using the double index notation, we have:

$$C = A_{[2,1]}; \quad R = A_{[3,1]}; \quad D = A_{[2,2]}; \quad E = A_{[2,3]} \quad \text{and} \quad A = A_{[3,2]}.$$

However, we will rather use $C$, $D$, $E$, $R$ and $A$ to denote the pseudofactors, which are confounded with factors in this example.

## 6.2. Characters

The characters from $T$ into $C_M$ are all the linear combinations $A = a_1A_1 + \dots + a_sA_s$ of the treatment pseudofactors, with $a_j \in C_{\pi(A_j)}$. They belong to the group $T^*$ which is the dual of $T$, so that $T^* \cong (C_{p_1})^{s_1} \times \dots \times (C_{p_l})^{s_l}$. The characters and dual of $U$ are defined and represented in the same way.

As in Section 3, each character $A$ of $T^*$ is associated with a pseudofactorial effect, denoted by $e_\tau(A)$, which belongs to a unique pseudofactorial term in the ANOVA decomposition of the genuine effects. This term is identified by the non-zero coefficients of the character $A$. We use the same definitions and interpretations as before for the character subsets $E_i$, $\tilde{E}_i$, $E_l$ and $\tilde{E}_l$.

## 6.3. Generalised regular factorial designs and key matrices

**Definition 3.1** can be generalised to the more general setting of the present section (Kobilinsky and Monod, 1995; Pistone and Rogantin, 2008). The design $d$ should satisfy an appropriate generalisation of Eq. (1). If $\mathbf{u} \in U$ then $\pi(A_j)A_j(d(\mathbf{u})) = 0 \pmod{M}$, by (10), for $j = 1, \dots, s$. When $\mathbf{u} = \mathbf{0}$ then $A_j(d(\mathbf{u})) = \alpha_j$, and so $M$ must divide $\pi(A_j)\alpha_j$. Now consider the unit $\mathbf{u}$ defined by $u_i = 1$ and $u_k = 0$ if $k \neq i$. Then $A_j(d(\mathbf{u})) = \phi_{ij}M/\pi(V_i) + \alpha_j$, from (10). Hence $\pi(A_j)\phi_{ij}M/\pi(V_i) = 0 \pmod{M}$ and so $\pi(V_i)$ divides $\pi(A_j)\phi_{ij}$. If $\pi(V_i) \neq \pi(A_j)$ then $\pi(V_i)$ divides $\phi_{ij}$ and so $\phi_{ij}M/\pi(V_i) = 0 \pmod{M}$.

**Definition 6.1 (Regular Design).** A factorial design $d$ with $U$ and $T$ as sets of units and genuine combinations respectively is called regular if there are coefficients $\phi_{ij}$ and $\alpha_j$ in $C_M$ such that, for $j = 1, \dots, s$,

$$A_j \circ d = \phi_{1j}V_1 + \dots + \phi_{ij}V_i + \dots + \phi_{rj}V_r + \alpha_j, \qquad (11)$$

where $V_1, \dots, V_r$ are the unit pseudofactors, $A_1, \dots, A_s$ are the treatment pseudofactors, and the following two conditions are satisfied:

$$M \text{ divides } \pi(A_j)\alpha_j \quad \text{for all } j = 1, \dots, s, \qquad (12)$$

$$\pi(V_i) \text{ divides } \pi(A_j)\phi_{ij} \quad \text{for all } i = 1, \dots, r \text{ and } j = 1, \dots, s. \qquad (13)$$
Fix j, and put $p = \pi(A_j)$. If $\mathbf{t} = d(\mathbf{u})$ and (11)–(13) hold then

$$t_j = \frac{p}{M} A_j(t) = \sum_{i=1}^{r} \frac{p}{M} \phi_{ij} \frac{M}{\pi(V_i)} + \frac{p}{M} \alpha_j = \sum_{i}' \phi_{ij} + \beta_j \pmod{p},$$

where $\beta_j = p\alpha_j/M$ (mod $p$) and the summation in $\sum'$ is restricted to those $i$ for which $\pi(V_i) = p$, and $\phi_{ij}$ is interpreted modulo $p$ if $\pi(V_i) = p$.

**Proposition 6.1.** A factorial design $d$ is regular if and only if the combination $\mathbf{t} = (t_1, \dots, t_s)^\top$ allocated to unit $\mathbf{u} = (u_1, \dots, u_r)^\top$ satisfies $\mathbf{t} = \Phi^\top \mathbf{u} + \mathbf{t}_0$, where the calculation of $t_j$ is performed modulo $\pi(A_j)$,

$$\Phi = \begin{pmatrix} \phi_{11} & \cdots & \phi_{1s} \\ \vdots & \ddots & \vdots \\ \phi_{r1} & \cdots & \phi_{rs} \end{pmatrix},$$

$\mathbf{t}_0 = (\beta_1, \dots, \beta_s)^\top$, $\phi_{ij} = 0$ if $\pi(V_i) \neq \pi(A_j)$, and, for $j = 1, \dots, s$, $\beta_j \in C_{\pi(A_j)}$ and $\phi_{ij} \in C_{\pi(A_j)}$. In particular, $\Phi$ is block diagonal:

$$\Phi = \mathrm{diag}(\check{\Phi}_1, \dots, \check{\Phi}_l), \qquad (14)$$

where the block $\check{\Phi}_k$ corresponds to the prime $p_k$.

**Definition 6.2 (Key Matrix).** The matrix $\Phi$ is called the key matrix of $d$.

Concerning the characters, we have the same relationships as in Section 3.4. The mapping $\psi: U \to T$ defined by $\mathbf{u} \mapsto \mathbf{t} = \Phi^\top \mathbf{u}$ is a group homomorphism from $U$ into $T$. If the dual of the homomorphism $\psi$ is denoted by $\varphi$ and if $\varphi(A) = B$ (with $A \in T^*$ and $B \in U^*$), then we have $\mathbf{b} = \Phi \mathbf{a}$, where $\mathbf{a}$ and $\mathbf{b}$ are the vectors of coefficients of $A$ and $B$.

## 6.4. Decomposition into Sylow subgroups

For $k = 1, \dots, l$, the elements of $U$ of order $p_k$ or 1 form a subgroup $\tilde{U}_k$ isomorphic to $(C_{p_k})^{r_k}$. These subgroups are called the Sylow subgroups of $U$. By the fundamental theorem of abelian groups, $U$ is the direct sum $\tilde{U}_1 \oplus \cdots \oplus \tilde{U}_l$. If $r_k = 0$ then $\tilde{U}_k = \{0\}$. Similarly, $T$ has Sylow subgroups $\tilde{T}_k$ isomorphic to $(C_{p_k})^{s_k}$ for $k = 1, \dots, l$, and $T = \tilde{T}_1 \oplus \cdots \oplus \tilde{T}_l$.

Likewise, the duals are direct sums of their Sylow subgroups, which are the duals of those of $U$ and $T$: that is, $U^* = \tilde{U}_1^* \oplus \cdots \oplus \tilde{U}_l^*$ and $T^* = \tilde{T}_1^* \oplus \cdots \oplus \tilde{T}_l^*$. For instance, the character $A$ in $T^*$ associated with the element $(\tilde{A}_1, \dots, \tilde{A}_l)$ of $\tilde{T}_1 \oplus \cdots \oplus \tilde{T}_l$ is the mapping $(\tilde{t}_1, \dots, \tilde{t}_l) \mapsto \tilde{A}_1(\tilde{t}_1) + \cdots + \tilde{A}_l(\tilde{t}_l)$, provided all the characters take their values in the common cyclic group $C_M$. We can write $A = (\tilde{A}_1, \dots, \tilde{A}_l)$ or $A = \tilde{A}_1 + \cdots + \tilde{A}_l$.

The decomposition into Sylow subgroups corresponds to the block diagonal decomposition of $\Phi$ in Eq. (14), because $\check{\Phi}_k$ is the restriction of $\Phi$ to $\tilde{U}_k$ and $\tilde{T}_k$.

**Definition 6.3 (Primary Components).** The character $\tilde{A}_k$ is called the $p_k$-primary component of the character $A$.

**Definition 6.4 (Sylow Components).** The diagonal blocks $\check{\Phi}_1, \dots, \check{\Phi}_l$ of the matrix $\Phi$ associated with the distinct primes $p_1, \dots, p_l$ are called the Sylow components of $\Phi$.

For the homomorphism $\varphi$, we have $\varphi(A) = \varphi(\tilde{A}_1) + \cdots + \varphi(\tilde{A}_l)$, and Proposition 6.1 implies that

$$\varphi(A) = \check{\varphi}_1(\tilde{A}_1) + \cdots + \check{\varphi}_l(\tilde{A}_l), \qquad (15)$$

where each $\check{\varphi}_k$ is the homomorphism from $\tilde{T}_k^*$ to $\tilde{U}_k^*$ associated with the matrix $\check{\Phi}_k$. The search for a matrix $\Phi$ meeting the requirements can thus be decomposed into the search for its Sylow components. The Sylow component $\check{\Phi}_k$ does not appear explicitly in $\Phi$ if there is no treatment pseudofactor with $p_k$ levels ($r_k \neq 0$ and $s_k = 0$) or if no unit pseudofactor has $p_k$ levels ($r_k = 0$ and $s_k \neq 0$). In the latter case, the regular designs associated with $\Phi$ give a constant value to genuine factors having $p_k$ levels, and this is usually prohibited unless the design is part of a larger one.

Section 7.5.1 shows that the Sylow components can often be searched for independently. Section 7.5.3 provides a backtrack search method otherwise.

# 7. Search for key matrices of generalised regular designs

## 7.1. Main steps

The search for a key matrix follows the same main steps as those in Section 5.1. The first step is to determine the set $\mathcal{I}$ of ineligible factorial terms. Let $I$ and $J$ be different subsets of $\{1, \dots, h\}$. If $i \in I \cap J$ and $n_i$ is not prime then $\mathcal{P}(i)$ contains
<caption>**Table 2**</caption>

Terms in the ineligible set (Example 2), obtained by Proposition 4.3 (redundant terms and the mean are crossed out).

<table><thead><tr><th rowspan="2">Set <i>&#x211D;</i></th><th colspan="10">Model <i>&#x211D;</i></th></tr><tr><th>&mu;</th><th>C</th><th>R</th><th>C.R</th><th>A</th><th>D</th><th>E</th><th>D.A</th><th>E.A</th><th>D.E</th></tr></thead><tbody><tr><td>D.A</td><td>D.A</td><td>C.D.A</td><td>D.R.A</td><td>C.D.R.A</td><td>D</td><td>A</td><td>D.E.A</td><td>&#x2213;</td><td>D.E</td><td>E.A</td></tr><tr><td>E.A</td><td>E.A</td><td>C.E.A</td><td>E.R.A</td><td>C.E.R.A</td><td>E</td><td>D.E.A</td><td>A</td><td>D.E</td><td>&#x2213;</td><td>D.E.A</td></tr></tbody></table>

at least two indices, and so $\tilde{E}_I - \tilde{E}_J$ contains some characters with at least one non-zero coefficient $a_j$ for some $j$ in $\mathcal{P}(i)$, and also some characters for which $a_j = 0$ for all $j$ in $\mathcal{P}(i)$. As in Proposition 4.2,

$$ \tilde{E}_I - \tilde{E}_J = \bigcup_{\mathcal{F}(K) \in \mathcal{I}(I,J)} \tilde{E}_K. $$

Hence Proposition 4.3 remains true in this more general setting, and so $\mathcal{I}$ can be determined as in Section 4.1. Concerning the second step, a reduced ineligible set of characters $\mathcal{R}$ can be deduced from $\mathcal{I}$ by the following steps.

(a) Deduce from $\mathcal{I}$ a set $\mathcal{I}_p$ of ineligible pseudofactorial terms when decomposing the factors into pseudofactors. This set can be reduced as explained in Section 7.3.

(b) Deduce from $\mathcal{I}_p$ the equivalence classes to be considered for $\mathcal{R}$.

(c) Select one representative character in each class and eliminate representatives having a proper multiple in $\mathcal{I}$.

In Examples 1 and 2, the models are complete and so the first main step is based on Proposition 4.3.

**Example 1 (Continued).** Proceeding as in Example 4 (Section 4.1.2), we find the ineligible set

$$ \mathcal{I} = \{F_1, F_2, F_3, F_4, F_1.F_2, F_1.F_3, F_1.F_4, F_2.F_3, F_2.F_4, F_3.F_4, F_1.F_2.F_3, F_1.F_3.F_4\}. $$

**Example 2 (Continued).** The set $\mathcal{M}$ is complete, so Proposition 4.3 shows that the terms in the ineligible set specified by $\mathcal{M}$ and $\mathcal{E}$ are those given in Table 2. Furthermore, the design must contain all the level combinations of factors C (columns) and R (rows), so the factorial terms C, R and C.R are also ineligible (see Section 4.4). Therefore

$$ \mathcal{I} = \{C, D, E, R, A, C.R, D.E, D.A, E.A, D.E.A, C.D.A, C.E.A, D.R.A, E.R.A, C.D.R.A, C.E.R.A\}. \quad (16) $$

## 7.2. Main principle of the reduction

Section 5.3 explained how the set $\mathcal{I}$ of ineligible characters can be reduced for elementary regular designs. The basic principle is that, if characters A and B are ineligible and there is a $\delta \in \mathbb{N}$ such that $A = \delta B$, then it is sufficient to check that $\varphi(A) \neq 0$. Thus the character B can be omitted from the set of ineligible characters. It follows that, to be parsimonious, a reduced ineligible set $\mathcal{R}$ must include only one representative per equivalence class defined by (6).

The same principle applies to generalised regular designs, but it is possible to go further. Indeed the relation

$$ \langle A \rangle \subseteq \langle B \rangle \iff \exists \delta \in \mathbb{N} \text{ such that } A = \delta B \quad (17) $$

defines a partial order on equivalence classes and it is clear that (4) has to be checked only for representatives of minimal classes. The class of A is minimal if and only if $\langle A \rangle$ contains no proper subgroup $\langle B \rangle$ with $B \in \mathcal{I}$. A reduced ineligible set $\mathcal{R}$ can thus be obtained by picking one representative in each equivalence class and avoiding representatives having a proper multiple in $\mathcal{I}$.

## 7.3. Reduction of the set of ineligible pseudofactorial terms

Suppose that the interaction $F_1.F_2$ is ineligible and that $F_1 = A_1 \wedge A_2, F_2 = B_1 \wedge B_2$. The character set $\tilde{E}_{(1,2)}$ of the interaction $F_1.F_2$ is the union of the character sets of the nine pseudofactor interactions in $\tilde{\mathcal{P}}(F_1.F_2)$, where

$$ \tilde{\mathcal{P}}(F_1.F_2) = \{A_1.B_1, A_2.B_1, A_1.A_2.B_1, A_1.B_2, A_2.B_2, A_1.A_2.B_2, A_1.B_1.B_2, A_2.B_1.B_2, A_1.A_2.B_1.B_2\}. $$

We call a term such as $A_1$ or $A_1.B_1$ a pseudofactorial term. Note that a pseudofactorial term is ineligible if and only if it is part of an ineligible factorial term.

In some cases, only a few of the pseudofactor interactions such as the nine above need to be considered when determining $\mathcal{R}$. This property is related to the Sylow decomposition of T and was not relevant for the elementary regular designs. It motivates the construction of a reduced set $\mathcal{I}_p$ of ineligible pseudofactorial terms, intermediate between the set $\mathcal{I}$ of ineligible factorial terms and the reduced set $\mathcal{R}$ of ineligible characters.
In the pseudofactor notation $A_{[k,j]}$ introduced in Section 6.1, each pseudofactorial term can be expressed as a product $\prod_{k \in K} \prod_{j \in J_k} A_{[k,j]}$, where the index $k$ varies over a subset $K$ of $\{1, \dots, l\}$ and, for each $k, J_k$ is a non-empty subset of $\{1, \dots, s_k\}$. The associated characters in $T^*$ are the linear combinations $\sum_{k \in K} \sum_{j \in J_k} a_{[k,j]} A_{[k,j]}$ for which every coefficient $a_{[k,j]}$ is non-zero. The set of these characters is the sum $\bigoplus_{k \in K} E_k$ of the pseudofactorial sets $E_k$ of characters associated with the pseudofactorial terms $\prod_{j \in J_k} A_{[k,j]}$.

**Definition 7.1 (Support).** The support of the character $(\tilde{A}_1, \dots, \tilde{A}_l)$ is the set $\{k : 1 \le k \le l \text{ and } \tilde{A}_k \ne 0\}$. The support of the pseudofactorial term $\prod_{k \in K} \prod_{j \in J_k} A_{[k,j]}$ is the set $K$.

**Proposition 7.1.** Let $\prod_{k \in K} \prod_{j \in J_k} A_{[k,j]}$ be an ineligible pseudofactorial term. If $L$ is any proper subset of $K$ then the pseudofactorial term $\prod_{k \in L} \prod_{j \in J_k} A_{[k,j]}$ is different from $\mu$ and different from $\prod_{k \in K} \prod_{j \in J_k} A_{[k,j]}$. If there is any such subset $L$ such that $\prod_{k \in L} \prod_{j \in J_k} A_{[k,j]}$ is ineligible, then the ineligible set of characters can be reduced by removing all characters associated with $\prod_{k \in K} \prod_{j \in J_k} A_{[k,j]}$.

**Corollary 7.1.** Consider any two pseudofactors that decompose the same factor $F_i$ and have different prime numbers of levels. Any pseudofactorial term that includes both pseudofactors can be omitted from $\mathcal{I}_p$.

A more thorough elimination can proceed according to Algorithm 2 below, where `iptI` is the initial set $\mathcal{I}_p$ deduced directly from $\mathcal{I}_p$, `iptR` denotes the reduced set under construction, and `iptq`, `iptK` are temporary subsets of `iptI`.

<caption>**Algorithm 2 Reduction of ineligible pseudofactorial terms**</caption>

```text
iptI ← complete set of ineligible pseudofactorial terms
iptR ← ∅
for q = 1, ..., l-1 do
    iptq ← subset of elements in iptI with support of size q
    iptR ← iptR ∪ iptq
    iptI ← iptI \ iptq
    for each pseudofactorial term pft in iptq do
        determine the support L of pft
        iptK ← subset of elements in iptI whose restriction to the support L equals pft
        iptI ← iptI \ iptK
end for
end for
return iptR
```

**Example 1 (Continued).** The factorial terms in $\mathcal{I}$ are expanded as functions of the pseudofactors. In this process,

* $F_1$ gives $A_{[2,1]}$, $A_{[3,1]}$, $(A_{[2,1]} \cdot A_{[3,1]})$,

* $F_2$ gives $A_{[2,2]}$, $A_{[2,3]}$, $A_{[2,2]} \cdot A_{[2,3]}$,

* $F_1 \cdot F_2$ gives $A_{[2,1]} \cdot A_{[2,2]}$, $A_{[3,1]} \cdot A_{[2,2]}$, $(A_{[2,1]} \cdot A_{[3,1]} \cdot A_{[2,2]})$, $A_{[2,1]} \cdot A_{[2,3]}$, $A_{[3,1]} \cdot A_{[2,3]}$, $(A_{[2,1]} \cdot A_{[3,1]} \cdot A_{[2,3]})$, $A_{[2,1]} \cdot A_{[2,2]} \cdot A_{[2,3]}$, $A_{[3,1]} \cdot A_{[2,2]} \cdot A_{[2,3]}$, $(A_{[2,1]} \cdot A_{[3,1]} \cdot A_{[2,2]} \cdot A_{[2,3]})$,

and so on, where the terms between parentheses involve different primes for the same genuine factor and so may be omitted immediately, by Corollary 7.1.

Then the algorithm starts with the ineligible pseudofactorial terms with only one non-zero primary component (support size $q = 1$). Here only 22 of the 31 pseudofactorial terms with support $\{2\}$ are ineligible:

$$ \begin{aligned} & A_{[2,1]}, A_{[2,2]}, A_{[2,3]}, A_{[2,4]}, A_{[2,5]}, A_{[2,1]}, A_{[2,2]}, A_{[2,1]}, A_{[2,3]}, A_{[2,1]}, A_{[2,4]}, \\ & A_{[2,1]}, A_{[2,5]}, A_{[2,2]}, A_{[2,3]}, A_{[2,2]}, A_{[2,4]}, A_{[2,2]}, A_{[2,5]}, A_{[2,3]}, A_{[2,4]}, A_{[2,3]}, A_{[2,5]}, \\ & A_{[2,4]}, A_{[2,5]}, A_{[2,1]}, A_{[2,2]}, A_{[2,3]}, A_{[2,1]}, A_{[2,4]}, A_{[2,5]}, A_{[2,2]}, A_{[2,3]}, A_{[2,4]}, \\ & A_{[2,2]}, A_{[2,3]}, A_{[2,5]}, A_{[2,2]}, A_{[2,4]}, A_{[2,5]}, A_{[2,3]}, A_{[2,4]}, A_{[2,5]}, \\ & \text{and } A_{[2,2]}, A_{[2,3]}, A_{[2,4]}, A_{[2,5]}. \end{aligned} $$

However, all three of the pseudofactorial terms with support $\{3\}$, namely $A_{[3,1]}$, $A_{[3,2]}$ and $A_{[3,1]} \cdot A_{[3,2]}$, are ineligible. When considered as subsets of pseudofactors, all pseudofactorial terms with support $\{2, 3\}$ include one or more of these pseudofactorial terms with support $\{3\}$. Proposition 7.1 shows that they can be eliminated.

Since this first reduced set of ineligible elements includes only elements with one non-zero primary component, the same is true of any reduced ineligible set deduced from it. Section 7.5.1 shows that this very often occurs in practice and that it allows us to make the search separately for each prime. But it is not always true, as shown by Example 2.
<table><caption>**Table 3**
Representatives of equivalence classes in the ineligible set $\mathcal{I}_p$ in Example 2.</caption><thead><tr><th>Support</th><th>Representative characters</th></tr></thead><tbody><tr><td>{2}</td><td>C, D, E, D + E</td></tr><tr><td>{3}</td><td>R, A</td></tr><tr><td>{2, 3}</td><td>C + D + R + 2A, C + D + 2R + 2A, C + E + R + 2A, C + E + 2R + 2A</td></tr></tbody></table>

**Example 2** (Continued). The pseudofactorial terms are confounded with the factorial ones in this example so they are given in (16). Those with support of size one may include C, D or E for $p_1 = 2$, or R or A for $p_2 = 3$, which yields D, E, C, D.E, R and A. Among the other pseudofactorial terms, we can eliminate D.A, D.R.A, D.E.A, E.A, E.R.A and C.R, which have D, E, C or D.E as 2-primary component, and C.D.A and C.E.A, which have A as 3-primary component. The remaining terms are C.D.R.A and C.E.R.A. So the reduction of pseudofactorial terms leads to the set

$$ \mathcal{I}_p = \{C, D, E, D.E, R, A, C.D.R.A, C.E.R.A\}. $$

Section 3.3 shows that the subsets of characters associated with the factorial terms C.D.R.A and C.E.R.A are

$$ \{C + D + R + A, C + D + 2R + A, C + D + R + 2A, C + D + 2R + 2A\} $$

$$ \text{and } \{C + E + R + A, C + E + 2R + A, C + E + R + 2A, C + E + 2R + 2A\} $$

respectively. Picking one representative in each equivalence class gives the ten elements in Table 3, four of which have two non-zero primary components. Each equivalence class has one element if the coefficients of the three-level factors R and A are both 0, or two otherwise. If factors R and A are both involved, the representative selected is the one whose coefficient of A is $-1 = 2$ (mod 3).

### 7.4. Reduction of the set of ineligible characters

The first step in getting the reduced ineligible set is to select one representative in each equivalence class for the relation (6). This is easy if there is some canonical way of selecting unambiguously the representative in each class. The following proposition shows that a canonical representative can be formed by picking the canonical representative of each primary component.

**Proposition 7.2.** If $A = (\tilde{A}_1, \dots, \tilde{A}_l)$ then the equivalence class $\mathcal{A}$ of $A$ is the product of the equivalence classes $A_1, \dots, A_l$ of its primary components $\tilde{A}_1, \dots, \tilde{A}_l$: that is, $\mathcal{A} = A_1 \times \dots \times A_l$.

When looking for representatives of minimal classes, the following proposition is useful.

**Proposition 7.3.** Let $A = (\tilde{A}_1, \dots, \tilde{A}_l)$ and $B = (\tilde{B}_1, \dots, \tilde{B}_l)$ be the Sylow decompositions of elements $A$ and $B$ of $T^*$. Then $\langle A \rangle \subseteq \langle B \rangle$ if and only if $(\tilde{A}_k) \subseteq (\tilde{B}_k)$ for $k = 1, \dots, l$.

In other words, the class of $A$ is contained in the class of $B$ in the sense defined by (17) if and only if, for each $k \le l$, the class of $\tilde{A}_k$ is contained in that of $\tilde{B}_k$; that is, there exists an integer $\delta_k$ such that $\tilde{A}_k = \delta_k \tilde{B}_k$.

Given an ineligible pseudofactorial term $\prod_{k \in K} \prod_{j \in J_k} A_{[k,j]}$, we seek a set of representatives of the equivalence classes in the corresponding set $E$ of characters, by which we mean a set containing exactly one element in each equivalence class. As in Section 7.3, $E = \bigoplus_{k \in K} E_k$, where $E_k$ is the set of characters associated with $\prod_{j \in J_k} A_{[k,j]}$.

**Proposition 7.4.** For each $k$ in $K$, let $\mathcal{C}_k$ be a set of representatives of the equivalence classes in $\tilde{E}_k$. Put $\mathcal{C} = \bigoplus_{k \in K} \mathcal{C}_k$. Then $\mathcal{C}$ is a set of representatives of the equivalence classes in $\bigoplus_{k \in K} \tilde{E}_k$.

Thus the search for $\varphi$ is reduced to the search for the primary homomorphisms $\tilde{\varphi}_k$ for each prime $p_k$. Denote by $\mathcal{I}_k$ the set of ineligible characters $A$ whose support is $\{k\}$. Then $\tilde{\varphi}_k$ must satisfy

$$ \tilde{\varphi}_k(A) \neq 0 \quad \text{for all } A \text{ in } \mathcal{I}_k, $$

for $k = 1, \dots, l$. Sometimes this necessary condition is also sufficient for condition (4) to be satisfied; sometimes it is not. We discuss the two cases in Section 7.5.

### 7.5. Dependencies between Sylow components of the key matrix

#### 7.5.1. A condition leading to independent searches for each Sylow component

For $k = 1, \dots, l$, denote by $\tilde{\mathcal{R}}_k$ the set of all non-zero $p_k$-primary components $\tilde{A}_k$ of the characters in $\mathcal{R}$.
**Definition 7.2.** A subset of characters in $T^*$ is thin if all of its elements have support of size one.

**Proposition 7.5.** If $\mathcal{R}$ is thin then condition (5) on $\varphi$ is equivalent to the conjunction of the $l$ conditions

$$ \check{\varphi}_k(\check{A}_k) \neq 0 \quad \text{for every } \check{A}_k \in \check{\mathcal{R}}_k, \qquad (18) $$

for $k = 1, \dots, l$.

Hence if $\mathcal{R}$ is thin then searching for $\varphi$ satisfying (5) is equivalent to searching separately (and independently) for the primary homomorphisms $\check{\varphi}_k$ satisfying (18). In practice it is easy to check directly whether $\mathcal{R}$ is thin. But a question naturally arises: is this often satisfied in practice? **Proposition 7.6** gives a positive answer by giving a mild condition under which $\mathcal{R}$ is thin. On the contrary, **Example 7** illustrates a practical situation in which it is not.

A subset $\mathcal{S}$ of a group is said to be closed under integer multiplication if

$$ A \in \mathcal{S} \implies \delta A \in \mathcal{S} \quad \text{for every integer } \delta. $$

It is easy to show that the subsets closed under integer multiplication are unions of subgroups.

In practice, a set like $\mathcal{M}$ defining the model is often a union of subgroups of $T^*$ and is thus closed under integer multiplication. This is always true when $\mathcal{M}$ is complete. As to the set $\mathcal{E}$ of effects to estimate, it is usual that, if it contains an interaction, it also contains all effects marginal to it except for the mean. For instance, if it contains A.B.C, it also contains the main effects A, B, C and the two-factor interactions A.B, A.C, B.C. Under these assumptions, $\mathcal{E} \cup \{0\}$ and $\mathcal{M}$ are both closed under integer multiplication, and so is the difference $\mathcal{E} - \mathcal{M}$, which is $\mathcal{I} \cup \{0\}$.

**Proposition 7.6.** If $\mathcal{I} \cup \{0\}$ is closed under integer multiplication then there exists a thin reduced ineligible set $\mathcal{R}$. In particular, this holds if $\mathcal{E} \cup \{0\}$ and $\mathcal{M}$ are both closed under integer multiplication.

This result can be generalised easily to ineligible sets of characters derived from one or more sets of the form $\mathcal{I}_j$ described in Section 4.3: if, for every subset $J$ defining a stratum, the subsets $\mathcal{E}_j \cup \{0\}$ and $\mathcal{M}_j$ are closed under integer multiplication, then the reduced set $\mathcal{R}$ is thin.

**Proposition 7.6** gives as a particular case the following classical result (Bailey, 1985): the design is of resolution $R$ if all its Sylow components are.

### 7.5.2. Counter-examples

There are, however, situations, like in some criss-cross experiments, when $\mathcal{E}$ includes an interaction but not the main effects of the corresponding factors and when $\mathcal{R}$ is not thin. **Examples 2 and 7** were constructed to illustrate this situation and its different consequences.

**Example 2 (Continued).** As shown in Section 7.3, the reduced set contains characters of support size 2 and so $\mathcal{R}$ is not thin. However, this example still allows for separate solutions of the Sylow components of $\mathcal{F}$, as we now explain.

The unit pseudofactors are $V_{[2,1]}$, $V_{[2,2]}$, and $V_{[3,1]}$. There is no loss of generality in putting $\check{\varphi}_2(C) = V_{[2,1]}$ and $\check{\varphi}_3(R) = V_{[3,1]}$. For example, one possibility for $\mathcal{F}$ is

$$ V_{[2,1]} \quad \left( \begin{array}{ccccc} \check{c} & \check{D} & \check{E} & \check{R} & \check{A} \\ 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 \end{array} \right). $$

It is clear that $\check{\varphi}_2(D)$ and $\check{\varphi}_2(E)$ must be two of $V_{[2,1]}$, $V_{[2,2]}$ and $V_{[2,1]} + V_{[2,2]}$. Since the character A is in $\mathcal{R}$, $\check{\varphi}_3(A)$ must be $V_{[3,1]}$ or $2V_{[3,1]}$, whatever $\check{\varphi}_2$ is. If $\check{\varphi}_3(A)$ is $V_{[3,1]}$ then $\varphi(C + D + R + 2A) = \varphi(C + D)$ and $\varphi(C + E + R + 2A) = \varphi(C + E)$, so $\check{\varphi}_2(D)$ and $\check{\varphi}_2(E)$ must both be different from $\check{\varphi}_2(C)$. If $\check{\varphi}_3(A)$ is $2V_{[3,1]}$ then exactly the same is true due to $\varphi(C + D + R + A)$ and $\varphi(C + E + R + A)$.

In **Example 2**, the solutions for $\check{\varphi}_2$ do not depend on the solution for $\check{\varphi}_3$ or vice versa. Here is a similar example with 36 units where the choice for the two-level factors depends on the choice previously made for the three-level ones.

**Example 7.** This is a small modification of **Example 2**. Now there are six units in each of the six cells, and the factor A is no longer constrained to be coarser than R. The sets $\mathcal{M}$, $\mathcal{E}$ and $\mathcal{I}$ are unchanged, and $\mathcal{M}$ is complete.

The unit pseudofactors are $V_{[2,1]}$, $V_{[2,2]}$, $V_{[3,1]}$ and $V_{[3,2]}$. There is no loss of generality in putting $\check{\varphi}_2(C) = V_{[2,1]}$ and $\check{\varphi}_3(R) = V_{[3,1]}$. Then $\check{\varphi}_2(D)$ and $\check{\varphi}_2(E)$ must be two of $V_{[2,1]}$, $V_{[2,2]}$ and $V_{[2,1]} + V_{[2,2]}$, while $\check{\varphi}_3(A)$ can be any non-zero combination of $V_{[3,1]}$ and $V_{[3,2]}$. If $\check{\varphi}_3(A)$ is $V_{[3,1]}$ then $\varphi(C + D + R + 2A) = \varphi(C + D)$ and $\varphi(C + E + R + 2A) = \varphi(C + E)$, so $\check{\varphi}_2(D)$ and $\check{\varphi}_2(E)$ must both be different from $\check{\varphi}_2(C)$. If $\check{\varphi}_3(A)$ is not a multiple of $V_{[3,1]}$ then there is no such constraint on $\check{\varphi}_2(D)$ and $\check{\varphi}_2(E)$.
In this case, the search cannot be made independently in the Sylow components. There are the following three
fundamentally different possibilities for $\Phi$.

$$
\begin{array}{c}
\begin{pmatrix}
\check{C} & \check{D} & \check{E} & \check{R} & \check{A} \\
1 & 0 & 1 & 0 & 0 \\
0 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 \\
0 & 0 & 0 & 0 & 0
\end{pmatrix}
\quad
\begin{array}{c}
\check{V}_{[2,1]} \\
\check{V}_{[2,2]} \\
\check{V}_{[3,1]} \\
\check{V}_{[3,2]}
\end{array}
\quad
\begin{pmatrix}
\check{C} & \check{D} & \check{E} & \check{R} & \check{A} \\
1 & 0 & 1 & 0 & 0 \\
0 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 1
\end{pmatrix}
\quad
\begin{array}{c}
\check{V}_{[2,1]} \\
\check{V}_{[2,2]} \\
\check{V}_{[3,1]} \\
\check{V}_{[3,2]}
\end{array}
\quad
\begin{pmatrix}
\check{C} & \check{D} & \check{E} & \check{R} & \check{A} \\
1 & 1 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 1
\end{pmatrix}
\end{array}
$$

More generally, counter-examples arise when, say, $A = \check{A}_1 + \check{A}_2$ belongs to $\mathcal{R}$ but neither $\check{A}_1$ nor $\check{A}_2$ does. This may happen if an interaction has to be estimated but not the associated main effects. In that case, it may well be necessary and sufficient that either $\check{\phi}_1(\check{A}_1) \neq 0$ or $\check{\phi}_2(\check{A}_2) \neq 0$. So there may be solutions with $\check{\phi}_1(\check{A}_1) \neq 0$ and solutions with $\check{\phi}_1(\check{A}_1) = 0$. Thus the equivalence in Proposition 7.5 is not satisfied.

7.5.3. Backtrack search in the non-independent case

When $\mathcal{R}$ is not thin, as in Example 7, the Sylow components of $\varphi$ and $\Phi$ can be found by backtrack search. Assume that
$\check{\varphi}_j$ has already been defined for $j = 1, \dots, k-1$. Let $\mathcal{R}_{[k]}$ be the subset of elements $A = (\check{A}_1, \dots, \check{A}_k, 0, \dots, 0)$ in $\mathcal{R}$ having
$\check{A}_k$ as last non-zero primary component. If $A \in \mathcal{R}_{[k]}$ then $\varphi(A) = (\check{\varphi}_1(\check{A}_1), \dots, \check{\varphi}_k(\check{A}_k), 0, \dots, 0)$ and the choice of $\check{\varphi}_k$
ensures that $\varphi(A) \neq 0$. If there is any index $j$ between $1$ and $k-1$ such that $\check{\varphi}_j(\check{A}_j) \neq 0$, then $\varphi(A) \neq 0$ whatever the choice
of $\check{\varphi}_k$. Such elements therefore need not be considered in the search for $\check{\varphi}_k$, which can proceed with reduced ineligible set
$\mathcal{R}_{\ll k\rr}$ consisting of those elements $A$ in $\mathcal{R}_{[k]}$ for which $\check{\varphi}_j(\check{A}_j) = 0$ for $j = 1, \dots, k-1$. If it succeeds and $k < l$, it goes on to
find $\check{\varphi}_{k+1}$. If it fails and $1 < k$, it goes back and tries to find another choice for $\check{\varphi}_{k-1}$. The search finally fails if it goes back to
$k = 1$ and fails to find another $\check{\varphi}_1$. It finally succeeds if it reaches $k = l$ and finds an admissible $\check{\varphi}_l$.

In any case, the elementary step in the search for $\varphi$ is the search for the primary homomorphisms $\check{\varphi}_k$ for each prime $p_k$.

8. Discussion

Quite apart from the computational aspects, this paper shows great unity between different types of factorial design: fractional or not; one prime or many; blocked, split-plot, row-column, criss-cross, and so on. The approach using one or more model-estimate pairs ($M$, $\tilde{E}$) gives a unified framework. The set of ineligible factorial terms is at the centre of this framework, since it synthesises all the constraints associated with the users' specifications. The other central component is the design key, which determines the combinatorial and statistical properties of the design. Indeed the design problem essentially consists of finding a design key adapted to the set of ineligible factorial terms.

A few remarks must be made from a statistical point of view. Of course, once an initial design has been generated, it then needs to be randomised. Since the two steps are quite independent, we only focused on the first one in this paper. Another point is that we made no distinction between the key matrices, provided they are solutions to the design specifications of Section 2. To cope with finer criteria such as minimum aberration or maximum estimation capacity (see e.g. Mukerjee and Wu, 2006), the approach developed here gives the possibility (up to computational constraints) to get all solutions and then select the best ones either according to such a criterion or by looking in detail at properties of the designs more relevant to the application at hand. An efficient alternative for a user of R is to use the FrF2 R package (Grömping, 2014), which makes better use of such considerations but is restricted to factors at two levels.

The framework could be even more general. For example, if a factor has four levels, it is possible to associate it with the cyclic group $C_4$ rather than using two pseudofactors with two levels each; similarly for other primes and other powers. Several authors have extended the theory to this more general setting, and it was implemented in the initial version of planor. However, the work of Voss (1988, 1993) suggests that there is no practical benefit from the more general framework.

The algorithmic approach presented here to generate designs is based on backtracking, which aims at a complete exploration of the possible solutions. The drawback is that the computational burden becomes too hard when the number of factors or the degree of fractionating becomes too high. So there is clearly a need to improve the speed of the algorithm. There are many directions to do so, but we want to stress two of them.

Cheng and Tsai (2013) have shown how templates may be used for the design key in certain situations. Such a template enables us to fix one or more columns in the matrix $\Phi$. For instance, in the matrix given for Example 3 in Appendix B of the Supplementary material (see Appendix A) we lose nothing by making the columns for $\check{A}$, $\check{B}$ and $\check{C}$ the same as those for $Q$, $U_1$ and $U_2$ respectively.

The search could also be accelerated by making use of symmetries between factors or pseudofactors, with respect to the
design specifications. To do so efficiently, it might be better to implement the search in a language like GAP (2016), which
is expressly designed to cut down searches in this way.

The R package planor is available on the CRAN (Monod et al., 2012). It deals with the whole class of generalised regular
factorial designs presented here. In addition to generating such designs, it can randomise them appropriately if the block
factors and their hierarchy relationships define an orthogonal block structure. Application to the main three examples of
the paper and to one higher dimensional one is presented in Appendix C of the Supplementary material (see Appendix A).
A more detailed presentation will be the subject of another paper.

Acknowledgements

This paper is dedicated to the first author, André Kobilinsky, who died before the present version was completed. André
conducted most of this research and turned it into the first version of the *planor* software as early as 1995. The R package
*planor* has been developed much more recently with the helpful support of Annie Bouvier (INRA, UR MalAGE). The second
and third authors are grateful to the Isaac Newton Institute for Mathematical Sciences, Cambridge, UK. A decisive step in
the writing was their joint work there during the 2011 INI programme on Design and Analysis of Experiments. Support was
also provided by the French Research Agency (ANR), project Escapade (ANR-12-AGRO-0003).

**Appendix A. Supplementary material**

Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.csda.2016.09.003.

**References**

Addelman, S., Kempthorne, O., 1961. Some main-effect plans and orthogonal arrays of strength two. Ann. Math. Statist. 32, 1167-1176.

Bailey, R.A., 1977. Patterns of confounding in factorial designs. Biometrika 64, 597–603.

Bailey, R.A., 1985. Factorial design and Abelian groups. Linear Algebra Appl. 70, 349–368.

Bailey, R.A., 2008. Design of Comparative Experiments. In: Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press, Cambridge.

Bailey, R.A., Gilchrist, F.H.L., Patterson, H.D., 1977. Identification of effects and confounding patterns in factorial designs. Biometrika 64, 347–354.

Bose, R.C., 1947. Mathematical theory of the symmetrical factorial design. Sankhyā 8, 107–166.

Bose, R.C., Bush, K.A., 1952. Orthogonal arrays of strength two and three. Ann. Math. Statist. 23, 508–524.

Box, G.E.P., Hunter, J.S., 1961a. The 2<sup>k-p</sup> fractional factorial designs. Part I. Technometrics 3, 311–351.

Box, G.E.P., Hunter, J.S., 1961b. The 2<sup>k-p</sup> fractional factorial designs. Part II. Technometrics 3, 449–458.

Cheng, C.-S., 2014. Theory of Factorial Design: Single- and Multi-Stratum Experiments. Chapman and Hall/CRC, Boca Raton.

Cheng, C.-S., Mukerjee, R., 1998. Regular fractional factorial designs with minimum aberration and maximum estimation capacity. Ann. Statist. 26, 2289–2300.

Cheng, C.-S., Tsai, P.-W., 2013. Templates for design key construction. Statist. Sinica 23, 1419–1436.

Cochran, W.G., Cox, G.M., 1957. Experimental Designs, second ed. John Wiley & Sons, New York.

Finney, D.J., 1945. The fractional replication of factorial experiments. Ann. Eugenics 12, 291–301.

Fisher, R.A., 1942. The theory of confounding in factorial experiments in relation to the theory of groups. Ann. Eugenics 11, 341–353.

Franklin, M.F., 1985. Selecting defining contrasts and confounded effects in $p^{n-m}$ factorial experiments. Technometrics 27, 165-172.

Franklin, M.F., Bailey, R.A., 1977. Selection of defining contrasts and confounded effects in two-level experiments. Appl. Stat. 26, 321-326.

Fries, A., Hunter, W.G., 1980. Minimum aberration 2<sup>k-p</sup> designs. Technometrics 22, 601–608.

Grömping, U., 2014. R package FrF2 for creating and analyzing fractional factorial 2-level designs. J. Stat. Softw. 56, 1–56.

Grömping, U., Bailey, R.A., 2016. In: Kunert, J., Müller, C.H., Atkinson, A.C. (Eds.), Regular fractions of factorial arrays. MODA 11—Advances in Model-Oriented Design and Analysis. Springer, pp. 143–151.

Grömping, U., Hu, X., 2014. Generalized resolution for orthogonal arrays. Ann. Statist. 42, 918–939.

Hedayat, A.S., Sloane, N.J.A., Stufken, J., 1999. Orthogonal Arrays. Springer-Verlag, New York.

Kempthorne, O., 1952. The Design and Analysis of Experiments. John Wiley & Sons, New York.

Kobilinsky, A., 1985. Confounding in relation to duality of finite Abelian groups. Linear Algebra Appl. 70, 321–347.

Kobilinsky, A., 2005. PLANOR: program for the automatic generation of regular experimental designs. Version 2.2 for Windows. Technical report. MIA Unit, INRA Jouy-en-Josas.

Kobilinsky, A., Monod, H., 2012. PLANOR: an R package for the automatic generation of regular fractional factorial designs. Version 1.0. Technical report. MIA Unit, INRA Jouy-en-Josas.

Kobilinsky, A., Monod, H., 1991. Experimental design generated by group morphisms: an introduction. Scand. J. Statist. 18, 119–134.

Kobilinsky, A., Monod, H., 1995. Juxtaposition of regular factorial designs and the complex linear model. Scand. J. Statist. 22, 223–254.

Kuhfeld, W.F., Tobias, R.D., 2005. Large factorial designs for product engineering and marketing research applications. Technometrics 47, 132-141.

Ledermann, W., 1977. Introduction to Group Characters. Cambridge University Press, Cambridge.

Lewis, S.M., 1982. Generators for asymmetrical factorial experiments. J. Statist. Plann. Inference 6, 59–64.

Monod, H., Bailey, R.A., 1992. Pseudofactors: normal use to improve design and facilitate analysis. Appl. Stat. 41, 317–336.

Monod, H., Bouvier, A., Kobilinsky, A., 2012. A quick guide to PLANOR, an R package for the automatic generation of regular factorial designs. Technical report. MIA Unit, INRA Jouy-en-Josas. https://cran.r-project.org/web/packages/planor.

Morris, M.D., 2011. Design of Experiments. Chapman & Hall/CRC Press, Boca Raton.

Mukerjee, R., Wu, C.F., 2006. A Modern Theory of Factorial Designs. In: Springer Series in Statistics, Springer, New York.

Patterson, H.D., 1965. The factorial combination of treatments in rotation experiments. J. Agric. Sci. 65, 171–182.

Patterson, H.D., 1976. Generation of factorial designs. J. R. Stat. Soc. Ser. B Stat. Methodol. 38, 175–179.

Patterson, H.D., Bailey, R.A., 1978. Design keys for factorial experiments. Appl. Stat. 27, 335–343.

Payne, R., 2012. A Guide to Anova and Design in GenStat©, fifteenth ed. VSN International, Hemel Hempstead.

Pistone, G., Rogantin, M.-P., 2008. Indicator function and complex coding for mixed fractional factorial designs. J. Statist. Plann. Inference 138, 787-802.

Ryan, T.P., 2007. Modern Experimental Design. In: Wiley Series in Probability and Statistics, John Wiley & Sons, Hoboken.

SAS Institute Inc., 2010. SAS/QC©9.2: User's Guide, second ed. SAS Institute Inc., Cary, NC.

The GAP Group. 2016. GAP — Groups, Algorithms, and Programming, Version 4.8.4. http://www.gap-system.org.

Voss, D.T., 1988. Single-generator generalized cyclic factorial designs as pseudofactor designs. Ann. Statist. 16, 1723-1726.

Voss, D.T., 1993. Comparing the classes of single replicate generalized cyclic designs with and without pseudofactors. J. Statist. Plann. Inference 35, 113-120.

Yates, F., 1933. The principles of orthogonality and confounding in replicated experiments. J. Agric. Sci. 23, 108–145.

Yates, F., 1937. The Design and Analysis of Factorial Experiments. Technical Communication 35. Harpenden: Imperial Bureau of Soil Science.