\documentclass[a4paper]{article}
\usepackage{Sweave}

\usepackage[english]{babel}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{url}
\usepackage{a4wide}
\usepackage{color}

\newcommand{\PLANOR}{\texttt{PLANOR }}

\begin{document}

\title{A quick guide to PLANOR, an R package for the automatic
  generation of regular factorial designs} 

\SweaveOpts{engine=R}
%\VignetteIndexEntry{A quick guide to PLANOR}
\author{Monod H., Bouvier A., Kobilinsky A.\\[2mm] 
  INRA, UR 341, Unit\'e MIA-Jouy \\ 
  Math\'ematiques
et Informatique Appliqu\'ees \\ F78352 Jouy en Josas Cedex\\ France}
\date{\today}
\maketitle
\tableofcontents
\medskip



\section{Introduction}

The \PLANOR R package generates regular factorial designs for a wide
and flexible range of user specifications. The main motivation for
its creation was to generate full and fractional factorial
designs, but \PLANOR can also be used to construct and randomize
complete block designs, Latin squares, split-plot designs,
\emph{etc}. The main limitation is that \PLANOR generates
orthogonal designs only, which excludes most incomplete block designs
and diverges from an optimal design approach. 

\PLANOR is based on algebraic methods of construction and more
specifically on the key matrix method
\cite{bailey85},\cite{franklin85},\cite{patterson78}, described in detail
in \cite{kobiasu97}, \cite{kobimonod91}, \cite{kobimonod94}, and 
illustrated in \cite{cliquet94}. This method produces so-called {\em
  regular} designs in which factorial effects are either estimable
independently or completely confounded. The \PLANOR R package
originates from the \PLANOR software which was written in the APL
language by Andr\'e Kobilinsky. The initial \PLANOR manual
\cite{kobi05} has been adapted to the \PLANOR R package \cite{kobi11}
and gives more details on the theory than this short guide.

To generate a design with \PLANOR, the user provides information on
the design factors, on the anova model to be used when analysing the
results, and on the design size. He or she then asks \PLANOR to search
for one or more designs meeting the requirements. One \PLANOR function
gives a design solution directly. The design can be randomized at once
according to a block structure formula given by the user or it can be
randomized later by a specific function. Alternatively, the solutions, if
any, can be obtained as a list of design key matrices. Several
specific functions then allow to investigate the solutions' properties
and to print and store the resulting designs. \PLANOR can manage
factors with different numbers of levels. It can take into account
hierarchical relationships among factors. It is also possible to
control the confounding of treatments effects with block effects, like
in split-plot or criss-cross experiments.

This vignette presents the basic usage of \PLANOR. A more
comprehensive presentation is under preparation, as well as additional
package functions. More details are also available through the help
functions of the \PLANOR package. For an introduction to the design of
experiments, many textbooks are available in the statistics
literature. For the \PLANOR user, we particularly recommend
\cite{bai08} and \cite{kobiasu97} (in French).

Please note that \PLANOR is still under construction. We advise to
check that the designs obtained by \PLANOR behave as expected before
using them for a real experiment, by inspecting them and conducting
analysis on simulated data, for example.

\section{Construction and randomization of orthogonal designs~: one-step examples}
In \PLANOR, the experiment requirements are specified in three parts~:
\emph{(i)} the factors~; \emph{(ii)} the model and (optionally) the
subset of factorial effects to estimate~; \emph{(iii)} the design
size. All this information can be provided to the
\texttt{regular.design} or to the \texttt{planor.designkey} functions.
The function \texttt{regular.design} is simpler to use because it
integrates all the steps included in \PLANOR for design construction.
Indeed it gives the factorial design directly as a ready-to-use
dataframe, if it finds a solution to the user's specifications. We
illustrate this function by the construction of a few well known
classes of orthogonal designs.

\subsection{Full factorial design}
Suppose we want to construct a full factorial design for three factors
$A$, $B$, $C$ at 2, 2, 3 levels respectively. Then the R code is:
<<full1, keep.source=TRUE>>=
library("planor")
ABCfull <- regular.design(factors=list(A=1:2,B=1:2,C=1:3), 
                          model=~A*B*C, 
                          nunits=2*2*3, 
                          randomize=~UNITS)
print(ABCfull)
@ 

The first argument describes the factors. To get a full factorial
design, the model formula includes all interactions and the size of
the experiment is the product of the numbers of levels of all
factors. We assume that the user wants a completely randomized design,
so the randomization formula is limited to the \texttt{UNITS}
level. The result is a dataframe. Based here on simulated data, the
analysis gives:
<<full2, keep.source=TRUE>>=
set.seed(123)
ABCfull$Y <- runif(2*2*3)
ABCfull.aov <- aov(Y~A*B*C, data=ABCfull)
summary(ABCfull.aov)
@ 

\subsection{Complete block design}
From a combinatorial point of view, the complete block design is a
special case of the full factorial design. From a statistical point of
view, the main differences are that treatment and block effects are
usually assumed to be additive (no interaction) and that the
randomization takes into account the blocks.

Suppose there are 5 treatments and 4 blocks. The R code becomes:
<<cbd1, keep.source=TRUE>>=
CBD <- regular.design(factors=list(
               Block=1:4, Treatment=c("T1","T2","T3","T4", "T5")), 
                      model=~Block+Treatment, 
                      nunits=4*5, 
                      randomize=~Block/UNITS)
print(CBD)
@ 

Based on simulated data again, the analysis gives:
<<cbd2, keep.source=TRUE>>=
CBD$Y <- runif(20)
CBD.aov <- aov(Y~Block+Treatment, data=CBD)
summary(CBD.aov)
@ 

\subsection{Latin square}
The Latin square involves three factors at $n$ levels in $n^2$ units,
with all three factors being pairwise orthogonal. From a combinatorial
point of view, the Latin square is a fractional factorial design. From
a statistical point of view, it is often used to study a treatment
factor of interest and to control two block factors that are
completely crossed.

Suppose that in a sensory experiment, there are 4 products to compare
and 4 judges. Each judge tastes the four products during four
consecutive periods. The design must ensure that the main effects of
the product, judge and period factors are estimable and orthogonal,
assuming the interactions are negligible. The appropriate
randomization consists in permuting the judge and period
labels at random, independently. 

The R code is:
<<LS1, keep.source=TRUE>>=
LS <- regular.design(factors=list(
  Judge=c("J1","J2","J3","J4"), Period=1:4, Product=c("P1","P2","P3","P4")), 
                       model= ~Judge + Period + Product, 
                       nunits=4*4, 
                       randomize=~Judge+Period)
print(LS)
@ 

The analysis gives:
<<LS2, keep.source=TRUE>>=
LS$Y <- runif(16)
LS.aov <- aov(Y~Judge + Period + Product, data=LS)
summary(LS.aov)
@ 

\subsection{Fractional factorial design of resolution 5}
For the last example in this section, we consider the type of design
for which \PLANOR was originally conceived: a highly fractionated
factorial design. Such designs allow to cope with a large number of
factors in (much) fewer units than required by a full factorial
design. They have been used for real experiments for a long time
\cite{boxhun78} and, more recently, for computer experiments
(\emph{e.g.} \cite{coumon11}, \cite{lurtou09}).

A fractional design of resolution~5 is generated below for 10~factors
at 4~levels, assuming that $2^{10}=1024$ units are available instead
of $4^{10}=1~048~576$ for a full factorial design. A fraction of
resolution~5 guarantees that all terms can be estimated from a model
with main effects and two-factor interactions. It can be generated as
follows~:
<<FFD, keep.source=TRUE>>=
FFD <- regular.design(factors=LETTERS[1:10], nlevels=4, 
                      resolution=5, 
                      nunits=2^10)
print(dim(FFD))
print(FFD[1:5,])
@ 
Note that in \PLANOR syntax, it is equivalent, but shorter, to
specify \texttt{resolution = 5} rather than 
\begin{verbatim}
model = (A+B+C+D+E+F+G+H+I+J)^2.
\end{verbatim}

\section{Construction of regular factorial designs through the search for a
  design key}
We now adopt a more progressive way to construct the design. For this
reason, we focus on the \texttt{planor.designkey} function rather than
\texttt{regular.design}. In that case, design construction involves
two main steps~:
\begin{enumerate}
\item the search for key matrices (function \texttt{planor.designkey});
\item then the design generation and randomization (function \texttt{planor.design}).
\end{enumerate}
We start by a short technical subsection. It cannot go
into details, but we hope it helps to make a link with other
approaches to the construction of regular factorial designs.

\subsection{A very short technical point}
A key matrix of base $p$ in \PLANOR is a matrix of integers modulo
$p$, where $p$ is a prime. It encodes the information required to
construct a regular factorial design for factors at $p$
levels. 

Consider for example a design for $4$ factors $A$, $B$, $C$, $D$ at
$p=2$ levels in $2^3=8$ units, whereas a full factorial design would
require $2^4=16$ units. It is possible to construct a design which
allows to estimate the main effects of the factors assuming the three-
and four-factor interactions are negligible. The solution is explained
in many books on factorial designs (\emph{e.g.} \cite{boxhun78})~:
\begin{itemize}
\item assimilate the factors' levels to $0, 1 \bmod 2$;
\item make a full factorial design on $A$, $B$, $C$;
\item add the level of $D$ on each unit by the equation $D=A+B+C \bmod
  2$, called the \emph{defining relationship} of the design.
\end{itemize}
Then it can be shown that the interaction $A.B.C.D$ is confounded with
the general mean, the main effect $A$ is confounded with the
interaction $B.C.D$, etc.

In \PLANOR, this construction is encoded in the following key matrix
of base $2$~:
\begin{equation*}
  K= 
  \left(
    \begin{array}{cccc}
      1&0&0&1\\
      0&1&0&1\\
      0&0&1&1
    \end{array}
  \right)
\end{equation*}
The rows of $K$ are associated with three factors $U_1$, $U_2$, $U_3$
which are called the units factors. The idea is that the set of units
can be identified to the full factorial design on these units
factors. The columns of $K$ are associated with the treatment factors
$A$, $B$, $C$, $D$. Here the first column of $K$ means that in the
design, we must have $A=U_1$ (modulo 2). The second, third and fourth
columns mean $B=U_2$, $C=U_3$ and $D=U_1+U_2+U_3$, respectively. It
follows that the defining relationship $D=A+B+C \bmod 2$ will be
satisfied.

The core algorithm in \PLANOR basically constructs $K$ by searching
for its columns successively, using a backtrack algorithm. However,
there is also much pre-processing to turn the factors and model
specifications into appropriate constraints on the columns of $K$. In
particular, all factors are automatically decomposed into
pseudofactors which all have a prime number of levels, and the whole
problem is decomposed according to the different prime numbers
involved.

A detailed presentation of the methodology implemented in \PLANOR is
under preparation \cite{kobbai12}. See also the references given in
the introduction or \cite{pisrog08} for the extension of regular
factorial designs to the case when different primes are involved.

\subsection{Fractional designs with 2-level factors}
\subsubsection{Search for a design key}
Consider an experiment to study four treatment factors $A$, $B$, $C$,
$D$ at two levels, using two blocks of size four. A full factorial
design on the treatment factors would require 16 units. Only eight are
available so that a fractional design must be used. In addition, some
treatment effects are necessarily confounded with the block effect.

At first, we may look for a design adapted to the model that includes
the main effects of the block and treatment factors, as well as the
interactions between pairs of treatment factors~:
<<ex1, keep.source=TRUE>>=
ex1Key <- planor.designkey(factors=c("block","A","B","C","D"),nlevels=rep(2,5), 
                           model=~block+(A+B+C+D)^2, 
                           nunits=2^3)
@ 
It turns out that \PLANOR fails to find a solution. There is indeed
no solution to the problem.

For the second try, we keep the same model but relax the implicit
constraint to estimate all factorial terms in the model. This is done
by using the \texttt{estimate} argument of the
\texttt{planor.designkey} function. This argument is optional~: by
default, it is considered that all terms in the \texttt{model} formula
must be estimated. In contrast, we only require below that the main
effects of the treatment factors be estimable. It follows that we now
allow for designs in which two-factor interactions are mutually
confounded.
<<ex1a, keep.source=TRUE>>=
ex1Key <- planor.designkey(factors=c("block","A","B","C","D"),nlevels=rep(2,5), 
                           model=~block+(A+B+C+D)^2, 
                           estimate=~A+B+C+D, 
                           nunits=2^3)
@

During the search, the backtrack algorithm looks successively for new
columns to add to the key matrix. Succinct information is given to
check the algorithm progress (default argument
\texttt{verbose=TRUE}). The search stops as soon as all columns of the
key matrix have been found (default argument \texttt{max.sol=1}).

An alternative to using \texttt{planor.designkey} directly is to
provide the information on the experiment step by step with
the functions \texttt{planor.factors} and \texttt{planor.model}. The
idea is to store the results of these functions in R objects and use
them as arguments to \texttt{planor.designkey}. This may be
convenient, for example, when one wants to explore several possible
models and design sizes with the same set of factors.
<<ex1a.bis, keep.source=TRUE>>=
ex1Fac <- planor.factors(factors=c("block","A","B","C","D"), nlevels=rep(2,5),
                         block=~block)
ex1Mod <- planor.model( model=~block+(A+B+C+D)^2, estimate=~A+B+C+D )
ex1Key <- planor.designkey(factors=ex1Fac, model=ex1Mod, nunits=2^3)
@ 

\subsubsection{Design-key properties}
In both cases, the key matrix solution is stored in the object
\texttt{ex1Key}. Its detailed properties can be obtained by two
different functions. The \texttt{summary} function prints the key
matrix and the defining relationships associated with this key
matrix. More detailed information on the aliasing between factorial
effects is given by the function \texttt{alias}.

Note that we have used the optional \texttt{block} argument of
\texttt{planor.factors} (also available in \texttt{planor.designkey}).
It specifies the factors that should be considered as block (or
nuisance) factors. In \PLANOR, the distinction between
\emph{treatment} and \emph{block} factors is taken into account when
studying confounding and aliasing properties.

<<ex1b, keep.source=TRUE, results=verbatim>>=
summary(ex1Key, show="dtb")
alias(ex1Key)
@

\subsubsection{Design generation}
Last but not least, a design can be generated by the function
\texttt{planor.design}. The design itself is the object in slot
\texttt{design} of the more complex object generated by
\texttt{planor.design}. An option allows the design to be randomized,
according to a block structure formula that the user must specify
(option \texttt{randomize}).
<<ex1c, keep.source=TRUE, results=verbatim>>=
ex1Des <- planor.design(ex1Key)
print(getDesign(ex1Des))
ex1Rand <- planor.design(ex1Key, randomize=~block/UNITS)
print(getDesign(ex1Rand))
@

\subsection{Fractional designs with 3-level factors}
We keep the same example but with 3-level factors and a few
more options. The results are not shown for sake of brevity.

% results=hide => les resultats n' apparaitront pas
<<ex2, results=hide, keep.source=TRUE>>=
# ***************** EXAMPLE 2 *****************
# Four 3-level treatment factors and one 3-level block factor
# Model: block+(A+B+C+D)^2  -   Estimate: A+B+C+D
# N = 3^3 = 27 units
#
ex2Key <- planor.designkey(factors=c(LETTERS[1:4],"block"),
                           nlevels=rep(3,5),  
                           block=~block,
                           model=~block+(A+B+C+D)^2,
                           estimate=~A+B+C+D,
                           nunits=3^3, base=~A+B+C, max.sol=2)
summary(ex2Key)
summary(ex2Key)
ex2Des <- planor.design(ex2Key[2])
@

Two optional arguments of \texttt{planor.designkey} have been used,
first to specify that $A$, $B$ and $C$ should be used as basic
factors, and second to ask for two solutions whereas the default is
one. Both solutions are examined by \texttt{summary} and the second
one, say, is chosen by the user to generate a factorial design. When
\emph{basic} factors are specified, they are identified to the units
factors $U_i$ \cite{kobi05}. As a consequence, all combinations of the
basic factors are guaranteed to be included in the design. When
relevant, using basic factors is recommended because it can speed up
the search.

The following lines also work; they illustrate that the basic factors
need not be part of the model but they must have been declared in
\texttt{planor.factors}.
<<ex2, results=hide, keep.source=TRUE>>=
ex2Fac <- planor.factors(factors=c(LETTERS[1:4], "block", "BASE"),
                         nlevels=rep(3,6) )
ex2Mod <- planor.model(model=~block+(A+B+C+D)^2, 
                       estimate=~A+B+C+D )
ex2Key <- planor.designkey(factors=ex2Fac, 
                           model=ex2Mod, 
                           nunits=3^3,
                           base=~A+B+BASE, 
                           max.sol=2)
@

\subsection{Asymmetric fractional factorial designs}
A regular fractional factorial design is called mixed or asymmetric when the numbers
of levels of the factors involve several different prime numbers. The
asymmetric designs constructed in \PLANOR consist of the cross products of
designs based on each prime. This does not allow for a great flexibility
in terms of confounding, but it enlarges the scope of situations
that can be addressed. 

% cf Planor Manual, Exemple 2 page 12
<<mix, keep.source=TRUE, results=verbatim>>=
# Four treatment factors at 6, 6, 4, 2 levels and one 6-level block factor
# Model: block+(A+B+C+D)^2 ; Estimate: A+B+C+D\n")
# N = 144 = 2^4 x 3^2 experimental units
mixKey <- planor.designkey(factors=c( LETTERS[1:4], "block"), 
                           nlevels=c(6,6,4,2,6), 
                           block=~block,
                           model=~block+(A+B+C+D)^2, 
                           estimate=~A+B+C+D, 
                           nunits=144,
                       base=~A+B+D, max.sol=2)
summary(mixKey)
mixPlan <- planor.design(key=mixKey, select=c(1,1), randomize=~block/UNITS)
print(getDesign(mixPlan)[1:25,])
@

The algorithm starts by decomposing the factors into pseudofactors
that all have a prime number of levels. Then it performs a similar
decomposition of the \texttt{model} and \texttt{estimate} terms. After
these initial steps, separate key-matrix searches are performed, one for
each prime involved in the problem. The prime decompositions are
automatic and transparent to the user. The recomposition when
generating a design is transparent too. In contrast, most
information on the search process and on the fraction properties are
given according to the prime decompositions.

\subsection{Split-plot designs}
In a split-plot experiment, there are two treatment factors
\texttt{variety} and \texttt{fert}, say, at $m$ and $n$ levels
respectively. The block structure consists of $r$ blocks each
containing $m$ sub-blocks of size $n$ and the factor \texttt{variety}
is constrained to be constant within sub-blocks. 

In an orthogonal split-plot design, each variety occupies one
sub-block of each block, and each sub-block contains the $n$ distinct
levels of factor \texttt{fert}. In \PLANOR, this design can be
constructed by defining the block structure as a cross between a
\texttt{block} and a \texttt{subblock} factor. The \texttt{hierarchy}
argument is used to specify that \texttt{variety} must be constant
within the combinations of \texttt{block} and \texttt{subblock}. Two
model-estimate pairs are given to the \texttt{listofmodels}
argument. First, the main effect of \texttt{fert} and the interaction
between \texttt{fert} and \texttt{variety} must be estimable when
blocks and sub-blocks are included in the model. Second, the main
effect of \texttt{variety} must be estimable between sub-blocks, that
is, when blocks but not sub-blocks are included in the model. The
command below calculates the design key of a split-plot design with
$r=2$, $n=2$, $m=2$.
<<splitplota, keep.source=TRUE, results=verbatim>>=
splitKey <- planor.designkey(factors=list(block=1:2,
                               subblock=1:2,
                               variety=LETTERS[1:2],
                               fert=c("organic","mineral")), 
                             block=~block+subblock,
                             hierarchy=list(~variety/(block*subblock)), 
                             listofmodels=
                             list(c( ~block*subblock+variety*fert, ~fert+fert:variety),
                                  c( ~block+variety,               ~variety)), 
                             nunits=2*2*2, 
                             base=~block+subblock)
summary(splitKey)
alias(splitKey)
print(getDesign(planor.design(splitKey, randomize=~block/subblock/UNITS)))
@ 

An alternative command to get the split-plot is given below. The main
difference is that the \texttt{subblock} factor now takes $rm$ levels
and is considered as nested in \texttt{block} rather than crossed with
it.
<<splitplotb, keep.source=TRUE, results=verbatim>>=
splitKey <- planor.designkey(factors=list(block=1:2,
                               subblock=1:4,
                               variety=LETTERS[1:2],
                               fert=c("organic","mineral")), 
                             block=~block+subblock,
                             hierarchy=list( ~block/subblock, ~variety/subblock), 
                             listofmodels=
                             list(c( ~subblock+variety*fert, ~fert+fert:variety),
                                  c( ~block+variety,               ~variety)), 
                             nunits=2*2*2, 
                             base=~subblock)
print(getDesign(planor.design(splitKey, randomize=~block/subblock/UNITS)))
@ 

\subsection{Fractional designs with nested factors and a complex
  block structure}
We now consider an experiment with concrete and more complex
specifications. This example stems from an experiment to study the
cleaning of surfaces by a robot, see \cite{kobi05}, example~3 on
page~3. There are five treatment factors at 2 levels. The block
structure consists of four plates with 2 rows and 4 columns per plate,
resulting in 32 experimental units. In addition, the design must cope
with experimental constraints between treatment and block factors. The
treatment factors concentration (\texttt{conc}) and temperature
(\texttt{Tact}) must remain constant within a plate. The treatment
factors denoted by \texttt{nsoil} and \texttt{qsoil} must remain
constant within each column of each plate. Only treatment factor
rugosity (\texttt{Rug}) can be modified freely between experimental
units.

To begin with, we show how to specify user-defined factor levels, by
providing a list to the \texttt{factors} argument of
\texttt{planor.factors}. Then, experimental constraints are
specified through the \texttt{hierarchy} argument of
\texttt{planor.factors}.
% echo=FALSE => les commandes R n' apparaitront pas
<<ex3a, keep.source=TRUE, results=verbatim>>=
# ************ ROBOT1A EXAMPLE *************
# Block structure: 4 plates / (2 rows x 4 columns)
# Treatments: 4 2-level factors
# Hierarchy 1: conc constant in plate
# Hierarchy 2: Tact constant in plate
# Hierarchy 3: nsoil constant in plate x column
# Hierarchy 4: qsoil constant in plate x column
# N = 32 units
#
robotFac <- planor.factors( factors=list(
                            conc=c(1,3),
                            Tact=c(15,30),
                            nsoil=c("curd","Saint-Paulin"),
                            qsoil=c("0.01g","0.10g"),
                            Rug=c(0.25,0.73),
                            plate=1:4,
                            row=1:2,
                            col=1:4),
                           hierarchy=list(~conc/plate,
                           ~Tact/plate,
                           ~nsoil/(plate*col),
                           ~qsoil/(plate*col)))
@

This example requires several model-estimate combinations. The main
model-estimate pair contains all the treatment factorial effects but
no block effect. It guarantees that all treatment combinations will be
present in the design, since all treatment factorial effects are
required to be estimable in the model with no block effect. The second
model-estimate pair (\texttt{listofmodels} argument) ensures that the
\texttt{Rug} factor is orthogonal to block factors.
%Note that all factors in the design must
%appear in the model-estimate formulae~; it is not sufficient to declare
%them as basic factors.

<<ex3b, keep.source=TRUE, results=verbatim>>=
robotMod <- planor.model( model=~nsoil*qsoil*Rug*conc*Tact,
                          listofmodels=list(c(~plate+row+col+Rug, ~Rug)) )
@

The \texttt{base} option of the \texttt{planor.designkey} function is
used here to impose that experimental units be associated with the
combinations of the block factors.
<<ex3c, results=verbatim>>=
robotKey <- planor.designkey(factors=robotFac, model=robotMod,
                             nunits=32, base=~plate+row+col)
summary(robotKey[1])
robotDes <- planor.design(robotKey[1], randomize=~plate/(row*col))
print(getDesign(robotDes))
@


%---------------------------------------------------------------------------
\section*{Acknowledgements}
%---------------------------------------------------------------------------
This vignette was typed using the Sweave package (Leisch, 2002\cite{sweave02}).

%---------------------------------------------------------------------------
\bibliographystyle{smfplain}      %\bibliographystyle{smfalpha ou acm}
\bibliography{planor}
%---------------------------------------------------------------------------

%---------------------------------------------------------------------------
\end{document}
%---------------------------------------------------------------------------
