::: {.content-visible when-format="html"}

::: callout-note


- The PDF version of this paper can be found on arxiv: <https://arxiv.org/abs/2311.09705>
- The paper uses [edibble]{.pkg} R package version 1.1.0 which can be installed as:

```{r, eval = FALSE}
devtools::install_version("edibble", version = "1.1.0", repos = "http://cran.rstudio.org")
```


:::

:::




# Introduction

The critical role of data collection is well captured in the expression "garbage in, garbage out" -- in other words, if the collected data is rubbish, then no analysis, however complex it may be, can make something out of it. Experiments offer the highest degree of control in the method of data collection and are of a critical importance for validating or investigating hypotheses across numerous fields (e.g., agriculture, biology, chemistry, business, marketing, engineering).  A proper design of experiments is critical to ensure that the desired inference is valid and efficient based on the experimental data. Conducting experiments is usually expensive; thus improper designs are at best inefficient use of resources and, at worst, a complete waste of resources. Given such high stakes, it is clear that improving the practice of experimental design can translate to large gains by increasing the inherent value of the experimental data.

However, more holistically, experiments involve multiple steps that extend beyond generating an experimental design table or layout. An appropriate experimental design cannot be generated without context and subject matter expertise to identify the appropriate experimental factors [@hahnExperimentalDesignComplex1984; @steinbergExperimentalDesignReview1984]. Furthermore, the experimental design is useless if it is not carried out as intended or if there are issues with data entry. Generating an experimental design should not be seen as an isolated activity, but rather viewed as an interconnected activity within the life of the whole experiment. Even a fit-for-purpose analysis of experimental data is not an independent activity void of the knowledge of the experimental design. Naturally, we can develop systems for constructing experimental designs that support an integrated approach for the entire lifecycle of an experiment.

<!---

A recent survey on teaching of experimental designs by @smuckerProfilesTeachingExperimental2023 suggest that the use of software in the teaching of experimental design is ubiquitous.

--->

Many experimental design software programs do not offer an integrated approach. Rather, the software tightly integrates the algorithmic aspect into a given experimental structure, thus losing the flexibility to specify other structures [@tanakaCurrentStateProspects2022a]. Unfortunately, this tight integration causes friction when considering alternative algorithms and requires users to rebuild the specification from scratch for another software solution. In addition, because the processes to specify the design are generally not modular, users must process all aspects of the experimental design simultaneously, burdening their cognitive load. As an alternative framework to specify experimental designs, @tanakaUnifiedLanguageExperimental2023 presented a computational framework, called "the grammar of experimental designs", which allows for intermediate constructs of experimental designs. This framework adopts a process-based approach to support the gathering of required information for the entire experiment in a modular manner. 

In @tanakaUnifiedLanguageExperimental2023, an experimental design is treated as a mutable object that progressively builds to the final design by a series of functions that modify a targetted element of the experimental design. This system is applicable to a broad range of experimental designs with fixed levels for each experimental factor. The system can be understood as an analogy to grammar in linguistics; a vast array of experimental designs (sentences) is expressed by combining key functions (words) through the shared understanding of the rules (grammar). In this system, an intermediate construct of an experimental design is represented as a pair of directed acyclic graphs (DAGs), where the nodes of one DAG represent the experimental factors, whereas the nodes of the other DAG represent the levels of those factors. Users explicitly define the experimental factors, their roles (such as treatment, unit or record), and relationships between the nodes. These purposeful manipulations of the DAGs are invoked by a small collection of object methods. The names of these methods (e.g., "set units" and "assign treatments") are semantically aligned to raise conscious awareness of the user (and the reader) to the elements of the experimental design.  The [edibble]{.pkg} package in the [R]{.proglang} language [@rstats] was used to demonstrate the utility of the framework in @tanakaUnifiedLanguageExperimental2023, however, little explanation was given to the usage of the package. 



This paper extensively demonstrates the main features of the [edibble]{.pkg} package. @sec-usage provides an overview of its usage, @sec-define describes various methods for defining the experimental structure. @sec-examples presents additional examples based on real experiments. @sec-contrast contrasts the existing systems. Finally, we conclude the paper with a discussion in @sec-conclude.


```{r}
#| echo: false
options(pillar.print_min = 6)
options(pillar.print_max = 10)
```

# Usage overview {#sec-usage}

The [edibble]{.pkg} package is available on the Comprehensive R Archive Network ([https://cran.r-project.org/](https://cran.r-project.org/)) with the latest developments available on GitHub at [https://github.com/emitanaka/edibble](https://github.com/emitanaka/edibble). The code presented in this paper is based on version `r packageVersion("edibble")` of the [edibble]{.pkg} package.

The ultimate aim of the [edibble]{.pkg} package is to produce an experimental design table. The name reflects this aim, with  [edibble]{.pkg} standing for **e**xperimental **d**esign t**ibble**,  where [tibble]{.class} is a special type of [data.frame]{.class} from [tibble]{.pkg} package [@tibble]. The key functions to achieve this aim are illustrated in @fig-workflow. An example shown next demonstrates a quick overview of its usage. 

![The general workflow for using [edibble]{.pkg}. The design starts with either no data or with data. We then set the scene by defining the experimental factors with its roles and assign links between the factors and levels. The information is stored as an [edbl_graph]{.class} object enclosed in an [edbl_design]{.class} object that also encloses the [Provenance]{.class} object that peforms low-level manipulations of the graphs. Once the full experimental structure is defined, then the table can be produced and exported. The record factors can be validated or simulated. The usage details of the functions are presented in @sec-define.](images/concrete-syntax){#fig-workflow}

Consider Example 1, derived from the example in @baileyDesignComparativeExperiments2008 for a quick demonstration of the [edibble]{.pkg} package. The main aim of this experiment is to determine the best feed type for calves to gain weight.

::: callout-note 

## Example 1: Calf feeding

There were 8 pens, with 10 calves in each pen. The experimenter was interested in comparing the effects of the four types of feed on the calves' weights. Each calf was individually weighed. The four feeds were combinations of two types of hay, which were put into the pen, and two types of anti-scour treatments, which were administered individually to each calf. 

:::

To specify this experiment in [edibble]{.pkg}, we composed it using a series of functions as shown below. The pipe (`%>%`) function, imported from [magrittr]{.pkg} [@magrittr], allows for a series of nested operations and can be substituted with the native pipe (`|>`) available from [R]{.proglang} version 4.1.0 onwards. This style of coding would be familiar to users of [tidyverse]{.pkg} [@wickhamWelcomeTidyverse2019].

```{r}
library(edibble)
calf <- design("Calf feeding") %>% # <1>
  set_units(pen = 8, # <2>
            calf = nested_in(pen, 10)) %>% # <2> 
  set_rcrds(weight = calf) %>% # <3>
  # set_trts(feed = 4) %>% # <4>
  set_trts(hay = 2, # <5>
           antiscour = 2) %>%  # <5>
  allot_table(hay ~ pen, # <6>
              antiscour ~ calf) # <6>
```

1. We begin the design specification by initiating the [edbl_design]{.class} object using [design]{.fct}. An optional title is provided, which is encoded and persistently displayed at various outputs (e.g., print or export of the object).
2. As described in the example, we had 8 pens with 10 calves in each pen. This part of the code specifies the unit factors "pen" and "calf". The right hand side shows the number of levels as a single numerical value. The function [nested_in]{.fct} is a helper function that encodes the calf is nested in pen.
3. As each calf is individually weighed, we set the records such that the weight is recorded for each calf (`weight = calf`). Here, the right hand side has to be a unit factor that has been previously defined. 
4. We are told that we are comparing four types of feed; hence, we set the treatments to `feed = 4`, but we realise later that this is not the entire description. We comment this line out and instead write our understanding in the next lines. You can, of course, choose to delete this line, but it can be useful to keep a record of your process.  
5. Here, we write our new understanding that the feed is composed of two types of hay and two types of anti-scour treatments.
6. In the final step, [allot_table]{.fct} is a short hand for [allot_trts]{.fct}, [assign_trts]{.fct}, and [serve_table]{.fct}; as such, multiple processes occur in this call. The relationship between the units and treatments is specified: hay type is alloted and assigned randomly to the pen, and anti-scour treatment type is alloted and assigned randomly to the calf within the pen. The design object is then translated into a table.

Printing the object above shows the table below, in which each experimental factor is translated into a column. 

```{r}
calf
```

The resultant object also contains the [edbl_design]{.class} object (hereafter referred plainly as the design object). The user may also plot the internal graphs using the [plot]{.fct} function. By default, this shows the factor graph as an interactive plot by internally using the [visNetwork]{.pkg} package [@visNetwork]. The level graph can be seen by adding the argument as `which = "levels"` in the [plot]{.fct} function. The static plots of the factor and level graphs are shown in @fig-graphs. The level graph can contain many nodes since every level of each experimental factor, except the record factors, is shown. 

```{r}
#| warning: false
#| message: false
#| eval: false
plot(calf) # factor graph is default
plot(calf, which = "levels") # level graph
```

![(A) Factor graph for the calf experiment. The shape and color of the nodes correspond to the role of the factor (red diamond = treatment, blue circle = unit, yellow three-dimensional cylinder = record). The labels correspond to the name of the factor and arrows show the relationship of the factors. Here we have the hay treatment factor alloted to pen and antiscour treatment factor alloted to calf. The weight is recorded for each calf. (B) A level graph for the calf experiment. The colors represent different experimental factors. The dropdown menu can be used to select a particular node or group of nodes that belong to the same experimental factor. Selecting a node highlights all connecting nodes, which can be useful for verifying the links between nodes.](images/graphs){#fig-graphs}







In experiments, there would be a response of interest. A response is often not required to generate an experimental layout. It is generally assumed that the smallest unit in the experimental design table is the observational unit (although not always so). In [edibble]{.pkg}, users can optionally specify their intention of what to record (e.g., responses) on which particular unit factor, as was done above, where the intention to capture the `weight` of each `calf` is specified in the function [set_rcrds]{.fct}. We can further add the expectations of the values for the `weight` factor. For example, below, we use [expect_rcrds]{.fct} to encode a data validation that the weight should be a numeric value between 0 and 10. 

```{r}
calfr <- calf %>% 
  expect_rcrds(weight > 0, weight < 10)
```

The benefit of encoding the data validation rule is that it is now interoperable. Two functionalities in [edibble]{.pkg} take advantage of this encoding. The first is when exporting the design through [export_design]{.fct}. This outputs an Excel file with the column "weight" in the data sheet for calves. The cells in this column are empty, with the intention that the data collector enters the data in this Excel sheet. These cells have the data validation embedded such that if the entry is not a numeric value between 0 and 10, it will result in an error in the data entry. See more details in @sec-expect.



The other functionality in [edibble]{.pkg} that takes advantage of the data validation encoding is when simulating data for record factor(s). There are two approaches to simulating data in [edibble]{.pkg}: one that requires users to extensively describe the simulation process ([simulate_process]{.fct} and then [simulate_rcrds]{.fct}, and the other that automatically writes the simulation scheme for you while ensuring to generate valid values: [autofill_rcrds]{.fct}. When users describe the simulation data on their own, it is possible to violate the valid values. In this case, the values are censored (by default to missing values) in [simulate_rcrds]{.fct}. More details on the simulation capabilities are presented in @sec-simulate.


# Defining structure {#sec-define}



An _experimental structure_ always consists of a _unit structure_; in other words, an experiment cannot be defined without units. If the experiment is comparative, then it must consist of a _treatment structure_ in addition to the mapping between the unit and treatment factors. This mapping can be defined on two levels: a high-level mapping between factors (usually for humans to understand the broad picture), and a low-level mapping between factor levels (usually assigned algorithmically). 

There are three roles for the experimental factors that can be encoded in [edibble]{.pkg}: treatment, unit, and record [see more details in @tanakaUnifiedLanguageExperimental2023]. The functions, [set_trts]{.fct}, [set_units]{.fct}, and [set_rcrds]{.fct} are used to encode the treatment, unit, and record factors, respectively, into the design object. The calls to these functions are generally associative, that is, the order does not matter nor does the number of calls to it matter (e.g., you do not need to define all units in one call and instead call on [set_units]{.fct} repeatedly). The exception is when the new factor is dependent on a previously defined factor; in this case, the dependent factor must be defined later.

Once the factors are defined, the relationship between them and their levels must be defined. The high-level grammar restricts the users in the type of relationships that can be defined (e.g., you cannot assign treatments to record factors), although low-level tweaks to the grammar (meant for developers) can bypass these restrictions (not presented in this paper). 

The main functions in [edibble]{.pkg} are explained in detail In the following subsections. In this section, we assume that the experimental aim is to find the best wheat variety from a wheat field trial.

## Initialisation

```{r}
library(edibble)
```

A new design constructed using [edibble]{.pkg} must start by initialising the design object. An optional title of the design may be provided as input. This information persists as metadata in the object and is displayed in various places (e.g., print output and exported files). 

### With data 

If you have existing data, you can use it as a base to build the design object. Below, we convert the data `gilmour.serpentine` in [agridat]{.pkg} [@agridat] into an [edbl_table]{.class} object.

```{r}
as_edibble(agridat::gilmour.serpentine, 
           .title = "Wheat Experiment", 
           .units = c(col, row, rep), .trts = gen) 
```

The main benefit of converting an existing experimental data into an [edbl_table]{.class} format is that other functionalities offered by [edibble]{.pkg} can be used. 

For the remainder of the paper, we assume experiments have no prior data that is used directly in the experimental design (as is commonly the case when you are conducting an experiment from scratch). 

### With no data

When you have no data, you start by simply initialising the design object.

```{r}
#| results: "hide"
design("Wheat field trial")
```

At this point, there is nothing particularly interesting. The design object requires the user to define the experimental factor(s) as described next.



## Units

At minimum, the design requires units to be defined via [set_units]{.fct}. In the code below, we initialise a new design object and then set a unit called "site" with 4 levels. The left hand side (LHS) and the right hand side (RHS) of the function input correspond to the factor name and the corresponding value, respectively. Here, the value is a single integer that denotes the number of levels of the factor. Note that the LHS can be any arbitrary (preferably syntactically valid) name. Selecting a name that succinctly describes the factor is recommended. Acronyms should be avoided where reasonable. We assign this design object to the variable called `demo`.

```{r}
demo <- design("Demo for defining units") %>% 
  set_units(site = 4)
```

At this point, the design is in a _graph form_. The print of this object shows a prettified tree that displays the title of the experiment, the factors, and their corresponding number of levels. Notice the root in this tree output corresponds to the title given in the object initialisation.

```{r}
demo
```



To obtain the _design table_, you must call on [serve_table]{.fct} to signal that you wish the object to be transformed into the _tabular form_. The transformation for `demo` is shown below, where the output is a type of [tibble]{.class} with one column (the "site" factor), four rows (corresponding to the four levels in the site), and the entries corresponding to the actual levels of the factor (name derived as "site1", "site2", "site3", and "site4" here). The first line of the print output is decorated with the title of the design object, which acts as a persistent reminder of the initial input. The row just under the header shows the role of the factor denoted by the upper case letter (here, U = unit) with the number of levels in that factor displayed via [pillar]{.pkg} [@pillar] with methods encoded for a custom vector class using the [vctrs]{.pkg} package [@vctrs]. If the number of levels exceed a thousand, then the number is shown with an SI prefix rounded to the closest digit corresponding to the SI prefix form (e.g., 1000 is shown as 1k and 1800 is shown as ~2k). The row that follows shows the class of the factor (e.g., character or numeric).

```{r}
serve_table(demo)
```



If particular names are desired for the levels, then the RHS value can be replaced with a vector like below where the levels are named "Narrabri", "Horsham", "Parkes" and "Roseworthy".

```{r}
design("Character vector input demo") %>% 
  set_units(site = c("Narrabri", "Horsham", "Parkes", "Roseworthy")) %>% 
  serve_table()
```

The RHS value in theory be any vector. Below the input is a numeric vector, and the corresponding output will be a [data.frame]{.class} with a numeric column. 

```{r}
design("Numeric vector input demo") %>% 
  set_units(site = c(1, 2, 3, 4)) %>% 
  serve_table()
```

In the instance that you do want to enter a single level with a numeric value, this can be specified using [lvls]{.fct} on the RHS. 

```{r}
design("Single numeric level demo") %>% 
  set_units(site = lvls(4)) %>% 
  serve_table()
```





### Multiple units

We can add more unit factors to this study. Suppose that we have 72 plots. We append another call to [set_units]{.fct} to encode this information. 

```{r}
demo2 <- demo %>% 
  set_units(plot = 72)
```

However, we did not defined the relationship between `site` and `plot`; so it fails to convert to the tabular form.

```{r}
#| error: true
serve_table(demo2)
```

The relationship between unit factors can be defined concurrently when defining the unit factors using helper functions. One of these helper functions is demonstrated next.

### Nested units {#sec-nested}

Given that we have a wheat trial, we imagine that the site corresponds to the locations, and each location would have its own plots. The experimenter tells you that each site contains 18 plots. This nesting structure can be defined by using the helper function [nested_in]{.fct}. With this relationship specified, the graph can be reconciled into a tabular format, as shown below.

```{r}
demo %>% 
  set_units(plot = nested_in(site, 18)) %>% 
  serve_table()
```

In the above situation, the relationship between unit factors have to be apriori known, but there are situations in which the relationship may become cognizant only after defining the unit factors. In these situations, users can define the relationships using the functions [allot_units]{.fct} and [assign_units]{.fct} to add the edges between the relevant unit nodes in the factor and level graphs, respectively (see @sec-allot and @sec-assign for more details of these functions). 



```{r}
demo2 %>% 
  allot_units(site ~ plot) %>% 
  assign_units(order = "systematic-fastest") %>% 
  serve_table()
```

The code above specifies the nested relationship of `plot` to `site`, with the assignment of levels performed systematically. The systematic allocation of `site` levels to `plot` is done so that the `site` levels vary the fastest, which is not the same systematic ordering as before. If the same result as before is desirable, users can define `order = "systematic-slowest"`, which offers a systematic assignment where the same levels are close together. 


### Crossed units {#sec-crossed}

Crop field trials are often laid out in rectangular arrays. The experimenter confirms this by alerting to us that each site has plots laid out in a rectangular array with 6 rows and 3 columns. We can define crossing structures using [crossed_by]{.fct}.

```{r}
design("Crossed experiment") %>% 
  set_units(row = 6,
            col = 3,
            plot = crossed_by(row, col)) %>% 
  serve_table()
```

The above table does not contain information on the site. For this, we need to combine the nesting and crossing structures, as shown next.

### Complex unit structures

Now, suppose that there are four sites (Narrabri, Horsham, Parkes, and Roseworthy), and the 18 plots at each site are laid out in a rectangular array of 3 rows and 6 columns. We begin by specifying the site (the highest hierarchy in this structure). The dimensions of the rows and columns are specified for each site (3 rows and 6 columns). The plot is a result of crossing the row and column within each site. 

```{r}
#| echo: false
options(pillar.print_min = 20)
```


```{r}
complex <- design("Complex structure") %>% 
  set_units(site = c("Narrabri", "Horsham", "Parkes", "Roseworthy"),
            col = nested_in(site, 6),
            row = nested_in(site, 3),
            plot = nested_in(site, crossed_by(row, col))) 

serve_table(complex)
```

You may realise that the labels for the rows do not start with "row1" for Horsham. The default output displays distinct labels for the unit levels that are actually distinct. This safeguards for instances where the relationship between factors is lost, and the analyst will have to guess what units may be nested or crossed. However, nested labels may still be desirable. You can select the factors to show the nested labels by naming these factors as arguments for the `label_nested` in [serve_table]{.fct} (below shows the nesting labels for `row` and `col` -- notice `plot` still shows the distinct labels).

```{r}
serve_table(complex, label_nested = c(row, col))
```

You later find that the dimensions of Narrabri and Roseworthy are larger. The experimenter tells you that there are in fact 9 columns available, and therefore 27 plots at Narrabri and Roseworthy. The number of columns can be modified according to each site, as below, where `col` is defined to have 9 levels at Narrabri and Roseworthy but 6 levels elsewhere.

```{r}
complexd <- design("Complex structure with different dimensions") %>% 
  set_units(site = c("Narrabri", "Horsham", "Parkes", "Roseworthy"),
             col = nested_in(site, 
                      c("Narrabri", "Roseworthy") ~ 9,
                                                . ~ 6),
             row = nested_in(site, 3),
            plot = nested_in(site, crossed_by(row, col))) 

complextab <- serve_table(complexd, label_nested = everything())
table(complextab$site)
```

You can see above that there are indeed nine additional plots at Narrabri and Roseworthy. The argument for `label_nested` supports [tidyselect]{.pkg} [@tidyselect] approach for selecting factors. 

## Treatments {#sec-trt}

Defining treatment factors is only necessary when designing a comparative experiment. The treatment factors can be set similar to the unit factors using [set_trts]{.fct}. Below, we define an experiment with three treatment factors: variety (a or b), fertilizer (A or B), and amount of fertilizer (0.5, 1, or 2 t/ha).

```{r}
factrt <- design("Factorial treatment") %>% 
  set_trts(variety = c("a", "b"),
           fertilizer = c("A", "B"),
           amount = c(0.5, 1, 2)) 
```

The links between treatment factors need not be explicitly defined. It is automatically assumed that treatment factors are crossed (i.e., the resulting treatment is the combination of all treatment factors) with the full set of treatments shown via [trts_table]{.fct}. For the above experiment, there are a total of 12 treatments with the levels given below.

```{r}
trts_table(factrt)
```

The `factrt` cannot be served as an [edbl_table]{.class} object, since there are no units defined in this experiment and how these treatments are administered to the units.


### Conditional treatments {#sec-conditioned}

In some experiments, certain treatment factors are dependent on another treatment factor. A common example is when the dose or amount of a treatment factor is also a treatment factor. In the field trial example, we can have a case in which we administer no fertilizer to a plot. In this case, there is no point crossing with different `amount`s; in fact, the amount of no fertilizer should always be 0. We can specify this conditional treatment structure by describing this relationship using the helper function, [conditioned_on]{.fct}, as below. The "." in the LHS is a shorthand to mean all levels, except for those specified previously. 

```{r}
factrtc <- design("Factorial treatment with control") %>% 
  set_trts(variety = c("a", "b"),
           fertilizer = c("none", "A", "B"),
           amount = conditioned_on(fertilizer,
                                    "none" ~ 0,
                                         . ~ c(0.5, 1, 2)))
```

We can see below that the variety is crossed with other factors, as expected, but the amount is conditional on the fertilizer.

```{r}
trts_table(factrtc)
```


## Links

In [edibble]{.pkg}, each experimental factor is encoded as a node in the factor graph along with its levels as nodes in the level graph (@fig-graphs). The edges (or links) can only be specified after the nodes are created. The links define the relationship between the experimental factors and the direction determining the hierarchy with the nodes [see @tanakaUnifiedLanguageExperimental2023 for more information]. Often, these links are implicitly understood and not explicitly encoded, thus making it difficult to utilise the information downstream. By encoding the links, we can derive information and validate processes downstream. 

Users specify these links using functions that are semantically aligned with thinking in the construction of an experimental design. There are three high-level approaches to defining these links as summarised in @tbl-links and explained in more detail next.


| Approach | Functions | Modifies | Purpose |
| --- | --- | --- | ------ | 
| Within role group | [nested_in]{.fct}, [crossed_by]{.fct}, [conditioned_on]{.fct} | Both factor and level graphs  | Links between the nodes of the same role only. | 
| Allotment | [allot_trts]{.fct}, [allot_units]{.fct}, [set_rcrds]{.fct}, [set_rcrds_of]{.fct} | Factor graph only | Capture high-level links that are typically apriori known by the user.  |  
| Assignment | [assign_trts]{.fct}, [assign_units]{.fct} | Level graph only | Determine links between nodes, often algorithmically. | 

: The table lists the three main approaches to specify links between experimental factors or their levels, the functions relevant for each approach, which graph the functions modifies, and the general purpose of the functions.   {#tbl-links}

### Within role group {#sec-within}

The helper functions, [nested_in]{.fct} and [crossed_by]{.fct}, are demonstrated in @sec-nested and @sec-crossed, to construct nested and crossed units, respectively. The helper function, [conditioned_on]{.fct}, demonstrated in @sec-conditioned, constructs a conditional treatment structure. These helper functions concurrently draw links between the relevant nodes in both factor and level graphs. These links would be apriori known to the user and these helper functions are just semantically designed to make it easier for the user to specify the links between nodes. These helper functions only construct links between nodes belonging to the same role (i.e., unit or treatment). 

### Allotment {#sec-allot}

Links specified using an allotment approach designate high-level links between factors. In other words, this approach only draws edges between nodes in the factor graph, and almost always, these edges are intentionally formed by the user. The purpose of this approach is to capture a user's high-level intention or knowledge.

For demonstration, we leverage the previously defined unit (`complexd`) and treatment structures (`factrtc`). These structures can be combined to obtain the combined design object as below. 

```{r}
complexd + factrtc
```

The above design object does not describe the links between the treatments and units. The function [allot_trts]{.fct} ascribes the links between treatments to units in the factor graph.

```{r}
alloted1 <- (complexd + factrtc) %>% 
  allot_trts(    fertilizer ~ row,
             amount:variety ~ plot)
```




### Assignment {#sec-assign}

The [assign_trts]{.fct} (often algorithmically) draw links between the treatment and unit nodes in the level graph (conditioned on the existing links in the factor graph). An overview of the algorithm in [assign_trts]{.fct} is shown in @fig-assign-alg. 

![An overview of the treatment assignment algorithm in [assign_trts]{.fct}. For the given experimental structure, the algorithm iterates over the defined allotments. Each treatment to unit allotment then generates an order of treatment levels assigned to units based on the ordering algorithm selected.](images/assignment-algorithm){#fig-assign-alg}

There are five in-built assignment algorithms: "systematic-fastest" (synonym for "systematic"), "systematic-random-fastest" (synonym for "systematic-random"), "systematic-slowest", "systematic-random-slowest", and "random". The variation in systematic assignment results in repeated ordering with respect to the unit order, without regard to any unit structure. When the number of units is not divisible by the total number of treatments, the earlier treatment levels would have an extra replicate.
The "systematic-random-fastest" and "systematic-random-slowest" are systematic variants that ensure equal chances for all treatment levels to obtain an extra replicate by randomising the order of treatment levels before the systematic allocation of treatment to units proceeds. The "fastest" and "slowest" variants determine if treatment levels are fast or slow in varying across order of the unit (slow varying meaning that the same treatment levels will be closer together in unit order, whereas fast varying means the same treatment levels are spread out in unit order). An example, with three treatments allotted to ten units, illustrating the in-built assignment algorithms is shown in @fig-actual-assign. 

![The left shows the treatment structure (3 treatment levels) and unit structure (2 blocks with 5 units nested in each, respresented by the rounded rectangle and numbers within it representing the unit order) that collectively make up the experimental structure. The treatment is alloted to unit factor. The actual assignment depends on the ordering algorithm selected with six variations shown on the right.](images/actual-assignment){#fig-actual-assign}

Building on the previously defined structure and allotment, we define an algorithm to assign links between unit and treatment levels using the function [assign_trts]{.fct}. Below, we use a systematic ordering for the first allotment (fertilizer to row) then a random ordering for the second allotment (interaction of amount and variety to plot). An optional seed number is provided to ensure the generated design could be reproduced. The generated design is shown in @fig-design-output (A). 

```{r}
design1 <- alloted1 %>% 
  assign_trts(order = c("systematic", "random"),
              seed = 2023) %>% 
  serve_table(label_nested = c(row, col))
```

```{r design2}
#| echo: false
design2 <- (complexd + factrtc) %>% 
  allot_table(fertilizer ~ row,
                  amount ~ plot,
                 variety ~ plot, 
              order = c("systematic", "random", "random"),
              label_nested = c(row, col),
              seed = 2023)
```

```{r design3}
#| echo: false
design3 <- alloted1 %>% 
  assign_trts(order = c("systematic", "random"),
              seed = 2023, 
              constrain = list(row = "site", plot = "row")) %>% 
  serve_table(label_nested = c(row, col))
```


```{r, echo = FALSE}
#| echo: false
#| message: false
#| warning: false
#| label: fig-design-output
#| fig-cap: "(A)-(C) shows the treatment allocation to plots from the three generated designs (`design1`, `design2`, `design3`). A plot is represented as a tile laid out in a rectangular array for each site. The color of the tile shows the type of fertilizer that the plot was allocated. The opacity of the color represents the amount of fertilizer allocated to the plot (the amount is also written in each tile in grey). The variety (a or b) allocated to the plot is written on each tile. (D) shows the number of plots that were alloted particular level of the treatment for each of the three designs. "
#| fig-height: 8
#| fig-width: 8
library(ggplot2)
library(dplyr)
library(patchwork)
g <- design1 %>% 
  tibble::as_tibble() %>% 
  ggplot(aes(col, row, fill = fertilizer)) +
  geom_tile(color = "black", aes(alpha = amount)) +
  geom_text(aes(label = variety), size = 6, nudge_y = 0.2) +
  geom_text(aes(label = amount), size = 3.5, nudge_y = -0.3, color = "darkgrey") +
  facet_wrap(~site, scales = "free_x") +
  theme_void(base_size = 12) +
  scale_fill_manual(values = c("#1b85b8", "#c3cb71", "white")) +
  scale_alpha_continuous(breaks = c(0, 0.5, 1, 2)) +
  ggtitle("design1") +
  ggh4x::force_panelsizes(cols = c(0.6, 0.9)) +
  theme(axis.text = element_text(), strip.text = element_text(size = 12, face = "bold"),
        axis.text.x = element_text(angle = 90),
        plot.title = element_text(family = "mono", margin = margin(b = 10)),
        #legend.margin = margin(t = unit(100, "lines")),
        legend.title = element_text(face = "bold", margin = margin(l = unit(10, "lines"))))

gcount <- bind_rows(mutate(design1, design = "design1"),
          mutate(design2, design = "design2"),
          mutate(design3, design = "design3")) %>% 
  count(fertilizer, amount, variety, design) %>% 
  ggplot(aes(variety, n)) +
  geom_col(aes(group = interaction(fertilizer, amount)), fill = "white", position = "dodge", color = "black") +
  geom_col(aes(fill = fertilizer, alpha = amount, group = interaction(fertilizer, amount)), position = "dodge", color = "black") + 
  facet_wrap(~design) +
  theme_bw(base_size = 12) +
  scale_fill_manual(values = c("#1b85b8", "#c3cb71", "white")) +
  scale_alpha_continuous(breaks = c(0, 0.5, 1, 2)) +
  theme(strip.text = element_text(family = "mono"),
        panel.grid.major.x = element_blank()) +
  labs(y = "Number of plots") +
  scale_y_continuous(expand = expansion(add = c(0, 0.5))) +
  guides(fill = "none", alpha = "none")

(g + (g %+% tibble::as_tibble(design2) + ggtitle("design2"))) / ((g %+% tibble::as_tibble(design3) + ggtitle("design3")) + gcount) + plot_layout(guides = "collect", heights = c(1, 1)) + 
  plot_annotation(tag_levels = "A",
                  tag_prefix = "(",
                  tag_suffix = ")") & 
  theme(legend.position = 'bottom', 
        legend.box.margin = margin(t = unit(0, "lines")),
        plot.margin = margin(b = unit(0, "lines")))
```


While allotment (high-level allocation) and assignment (actual allocation) are distinguished in the system to provide flexibility to the user for defining these processes separately, it is likely that many users would concurrently define these processes. The [allot_table]{.fct} function offers a shorthand that combines the call to [allot_trts]{.fct}, [assign_trts]{.fct}, and [serve_table]{.fct} into one call. 

To illustrate the difference when treatment interaction is alloted to a unit (like the second allotment in `allotment1`), below, we have a different allotment where the amount of fertilizer and variety are allotted to plot in a separate allotment. A separate allotment can be assigned using different algorithms and is considered independent of other allotments (unless the treatment factor is conditional on another treatment factor). The resulting assignment is shown in @fig-design-output (B) for the design below. 

```{r design2}
```

The assignment algorithms in the system use the default constraint, which takes the nesting structure defined in the unit structure (i.e. row is nested in site and plot is crossed by row and column and nested in site). This constraint is used to define the nature of "random" assignment. For example, in the code below, we relax this constraint such that the `plot` factor is constrained within a `row` (default was `row`, `col` and `site`), which in turn is contained within the `site`. This difference in constraints results in a different path in the algorithm (as shown in the overview in @fig-assign-alg), which generates the design shown in @fig-design-output (C). 

```{r design3}
```

The above three different designs (`design1`, `design2` and `design3`) share the same unit and treatment structure, but the allotment and/or assignment algorithm differed.  One result of this is that the treatment replications, as shown in @fig-design-output (D), differ across the generated designs with the most ideal distribution seen in `design3` (if all fertilizer and amount combinations are of equal interest and fertilizer allocation is restricted to the row; arguably, it is better to remove the latter constraint, if practically feasible, so the units with the control treatment can be assigned for other treatment levels to obtain a more even distribution).  The difference in `design1` and `design2` is that the amount and variety were allocated as an interaction in the former but independently in the latter. The latter process does not ensure near-equal replication of the treatment levels, so it is not surprising that in @fig-design-output (D), `design2` has  the least uniform treatment distribution. 


Finding or creating the most appropriate assignment algorithm is one of the challenging tasks in the whole workflow. The default algorithm is unlikely to be optimal for the given structure, and the user is encouraged to modify this step to suit their own design. @sec-composition presents an example of this. 



## Records {#sec-rcrds}

The values of a record factor is unknown prior to the execution of the experiment; thus, the method used to define it fundamentally differs from other factors. The record factor is defined using [set_rcrds]{.fct} with the arguments in the form of `record = unit`, where `record` corresponds to user-specified name of the record, and `unit` corresponds to an existing unit factor. The LHS of the input argument is always a new factor, as it was the case for [set_units]{.fct} and [set_trts]{.fct}.  Multiple record factors can be specified as separate arguments within the same function call as shown below.

```{r}

record1 <- design1 %>% 
  set_rcrds(biomass = plot,
              yield = plot, 
           rainfall = site)
```

The above specification can be awkwardly lengthy when we expect multivariate data (i.e., multiple responses of the same unit). In this instance, it may be easier to specify the responses as follows.

```{r}
#| results: hide
design1 %>% 
  set_rcrds_of(plot = c("biomass", "yield"),
               site = "rainfall")
```

In the above specification, the LHS corresponds to an existing unit factor, while the RHS is a character vector where each element corresponds to the name of the record factor. Notice that in this specification, the RHS elements are quoted as they do not yet exist in the design. By contrast, the RHS specification in [set_rcrds]{.fct} is unquoted because these factors already exist. Because the intention of the LHS specification is to point to an existing unit factor, it differs from the convention of other functions that prefix with `set_`. As a signal to the different input structure, the function appends a suffix `_of`, which is supposed to act as a signal to the user that this is different from other `set_` functions. 

Record factors generally do not play a role in the assignment of treatments to units; however, they play a critical role in the analysis stage. The specification of these record factors allows for the user's intention to be explicit. The added benefits also include the encoding of the expected and simulated values, as described in detail next.

### Expected values {#sec-expect}

If the record factors are defined in the design object, like in `record1`, then the user may specify the expected values of the records. For example, below we specify that `biomass` is greater than or equal to zero, yield is between 0 and 10, and rainfall is a factor with two levels: high or low. 

```{r}

expect1 <- record1 %>% 
  expect_rcrds(biomass >= 0,
               yield > 0, yield < 10,
               factor(rainfall, levels = c("high", "low")))
```

Once the expected values are encoded, they can be used in various functions. For example, suppose we export this design as below.

```{r}
#| message: false
#| warning: false
export_design(expect1, file = "mydesign.xlsx", overwrite = TRUE)
```

The exported design automatically embeds the data validation in the corresponding cells, as shown in @fig-export.

![The exported design in the Excel sheet encodes the data validation defined in the design object such that the unexpected values cannot be entered in the corresponding entry. The data entry is designed to remove redundancy of duplicate entries (e.g. the rainfall status recorded for site, but yield and biomass on the plot) by creating separate data sheets for the corresponding unit factor.](images/export-design){#fig-export}

### Simulated values {#sec-simulate}

Another benefit of encoding the expected values of the record factors is that you can do a lazy simulation using [autofill_rcrds]{.fct}. This function randomly chooses a simulation scheme (including the variables that influence the record), while ensuring that it keeps to the expected values. 

```{r}
set.seed(2023)
sim1auto <- autofill_rcrds(expect1) 
```

The aforementioned lazy simulation is designed for quick diagnostics of planned analytical methods and not for any serious simulation study. For proper simulation schemes, users can freely enter their own schemes using the function [simulate_process]{.fct}. In this function, the user defines a series of functions, where each function returns either 1) a vector of the same size as the data or 2) a [data.frame]{.class} of the same row dimension as the data. For 1) the name of the function must correspond to the name of the record factor, and for 2) the name must start with a dot and the column names of the returning object must match the names of the record factors.

In the example below, we define two simulation processes called "yield" and ".multi". The first simulation process is a function that has arguments that define the variety main effect and combined effect of fertilizer and its amount; the return object is a numeric vector where random noise is added on top of the fixed treatment effects. The second simulation process is a function that simulates yield and biomass from a multivariate normal distribution.

```{r}
process1 <- expect1 %>% 
  simulate_process(
    yield = function(v = c(2, 4),
                     af = c(1, 2.5, 3, 0.5, 1, 5, 0)) {
      veffects <- setNames(v, c("a", "b"))
      afeffects <- setNames(af, c("A:0.5", "A:1", "A:2",
                                  "B:0.5", "B:1", "B:2", "none:0"))
      # must return vector of the same size at the data
      veffects[variety] + afeffects[paste0(fertilizer, ":", amount)] + rnorm(n())
    },
    
    .multi = function() {
      Sigma <- matrix(c(2, 0.9, 0.9, 1), nrow = 2)
      yveffects <- c("a" = 3, "b" = 9)
      bveffects <- c("a" = 1, "b" = 3)
      res <- cbind(yveffects[variety], bveffects[variety])
      res <- res + mvtnorm::rmvnorm(n(), mean = c(0, 0), sigma = Sigma)
      res <- as.data.frame(res)
      colnames(res) <- c("yield", "biomass")
      # must return a data.frame where the column names match the
      # record factor names and the number of rows match the 
      # row dimension of data
      res
    }
  )
```

After defining the simulation process, the actual simulation of the record is performed through a call to [simulate_rcrds]{.fct}. Below, we simulate the actual records by naming the processes we want to use along with any parameters of the function. Optionally, users can define values to censor the simulated values if they lie outside the expected values.

```{r}
sim1default <- simulate_rcrds(process1, yield = with_params())
sim1var <- simulate_rcrds(process1, 
                          yield = with_params(v = c(-1, 5), 
                                              .censor = c(0, 10)))
sim1multi <- simulate_rcrds(process1, .multi = with_params(), .seed = 1)
sim1censor <- simulate_rcrds(process1, 
                .multi = with_params(.censor = list(yield = c(0, 10),
                                                    biomass = 0)),
                .seed = 1)
```

The resulting simulation values from the above calls are shown in @fig-simulation. We can see in @fig-simulation (A) that the `sim1auto` shows a huge difference in yield with respect to fertilizer type (reflecting that the automated simulation process must have included fertilizer type in the data generation process); `sim1default` shows some difference in yield across variety (the actual difference is 2 in the simulation process); and `sim1var` shows a variation in the latter process where the difference in the main variety effects are now more pronounced (actual difference is 6). @fig-simulation (B) shows the simulation results from `sim1multi` and `sim1censor`; the difference in the latter is that the values are censored to 0 or 10 in the yield and 0 in biomass (the extra points in yield at 0 and 10, and biomass at 0 are easily seen in the plot).

```{r}
#| echo: false
#| label: fig-simulation
#| warning: false
#| fig-height: 5
#| fig-cap: "(A) show the simulated yields from `sim1auto`, `sim1default` and `sim1var` colored by fertilizer type with shape showing the variety, and size reflecting the amount of fertilizer. (B) show the scatterplot of biomass and yield from `sim1multi` and `sim1censor`. The only difference in the two simulated values are the censored values that are shown as red points in `sim1censor`."
library(dplyr)
g1 <- bind_rows(mutate(sim1auto, data = "sim1auto"),
          mutate(sim1default, data = "sim1default"),
          mutate(sim1var, data = "sim1var")) %>% 
  mutate(variety = factor(variety),
         amount = factor(amount),
         fertilizer = factor(fertilizer, c("A", "B", "none"))) %>% 
  ggplot(aes(data, yield)) +
  ggbeeswarm::geom_quasirandom(aes(color = fertilizer, size = amount, shape = variety,
                                group = variety),
                            alpha = 0.5, dodge.width=0.3) +
  scale_color_manual(values = c("#1b85b8", "#c3cb71", "grey")) +
  theme_classic(base_size = 12) +
  labs(x = "") +
  theme(axis.text.x = element_text(family = "mono"))

ind <- which(is.na(sim1multi$yield) | is.na(sim1multi$biomass))

g2 <- bind_rows(mutate(sim1multi, data = "sim1multi"),
          mutate(sim1censor, data = "sim1censor")) %>% 
  mutate(amount = factor(amount),
         fertilizer = factor(fertilizer, c("A", "B", "none")),
         data = factor(data, levels = c("sim1multi", "sim1censor"))) %>% 
  ggplot(aes(yield, biomass)) +
  geom_point(aes(shape = variety, size = amount), color = "red", show.legend = FALSE, 
             data = ~filter(., data == "sim1censor") %>% 
               slice(ind)) + 
  geom_point(aes(color = fertilizer, size = amount, shape = variety), alpha = 0.5) + 
  facet_wrap(~data) +
  scale_color_manual(values = c("#1b85b8", "#c3cb71", "grey")) +
  theme_classic(base_size = 12) +
  theme(strip.text = element_text(family = "mono"))

g1 / g2 + plot_layout(guides = "collect") +
  plot_annotation(tag_levels = "A",
                  tag_prefix = "(",
                  tag_suffix = ")")
```

# Worked examples {#sec-examples}

In addition to the example presented in @sec-usage, we demonstrate two examples from real experiments: one in ecology and one in linguistics.

## Common garden design

@cochraneClimateWarmingDelays2015 studied the temperature and moisture impact on seedling emergence in a common garden. Below is a digest of the experimental design description.

::: callout-note

## Example 2: Common garden design 

Twelve shelters were constructed on site to manipulate water availability and assess the effects on seed germination, seedling emergence, and seedling growth and survival. Two raised garden beds were placed below each shelter. Water was collected from the roof of each shelter and passively irrigated to the garden beds in each of the three water treatments: reduced rainfall R (80% of the total), natural rainfall N (100% of the total), and increased rainfall I (140% of the total). Each shelter contained one garden bed fitted with 24 warming chambers (W), and one garden bed without warming chambers (C).

Four non-sprouting, obligate seeding species (*Banksia media*, *Banksia coccinea*, *Banksia baxteri*, and *Banksia quercifolia*) were chosen for the study. Six discrete populations of each species were selected to represent a rainfall cline (High, Medium or Low). The term 'population' is used to describe plants originating from a particular geographic and climatic location.

The experimental units were arranged within each garden bed in a balanced array of three columns and eight rows, with seeds from all 24 populations sown in each garden bed. Each row contained seeds from a single species. Each column in the array contained seeds from eight populations corresponding to a particular level of rainfall. Within these constraints, the populations were divided into two sets of 12, so that the same 12 populations always occurred together in a set of four rows.


The allocation of seeds to the experimental units was fully randomised, subject to the constraints of the design.


:::

Based on the above description, my understanding of the unit structure was follows.

```{r}
garden_units1 <- design("Common garden experiment") %>% 
  set_units(shelter = 12,
            bed = nested_in(shelter, 2),
            row = nested_in(bed, 8),
            col = nested_in(bed, 3),
            plot = nested_in(bed, crossed_by(row, col)))
```

Initially, I did not understand the full treatment structure and wrote only partial treatment factors. From a further dissection of the description, I inferred that rainfall was the other treatment condition. At this point it was unclear to me how their so-called "population" is actually mapped to the rainfall and so I ignore it for the moment.

```{r}
garden_trts1 <- set_trts(water = c("R", "N", "I"),
                         chamber = c("W", "C"), 
                         species = c("media", "coccinea",
                                     "baxteri", "quercifolia"),
                         rainfall = c("High", "Medium", "Low"))
```





Based on the description, I inferred the following mapping between the treatment factors to unit factors. The default treatment assignment uses the nesting structures to define the constraint and allocates using the "random" order. A random generation of this design is shown in @fig-garden-design1.

```{r}
garden_design1 <- (garden_units1 + garden_trts1) %>% 
  allot_table(water ~ shelter,
              chamber ~ bed,
              species ~ row,
              rainfall ~ col,
              label_nested = c(row, col), 
              seed = 2023)
```

```{r fig-garden-design1}
#| echo: false
#| fig-height: 10
#| fig-width: 8
#| fig-cap: "The generated experimental design in `garden design1`."
garden_design1 %>% 
  mutate(bed2 = attr(bed, "label-nested")) %>% 
  mutate(rainfall = factor(rainfall, levels = c("High", "Medium", "Low"))) %>% 
  ggplot(aes(row, col)) +
  geom_tile(aes(fill = rainfall)) +
  geom_point(aes(color = species), size = 3) +
  scale_fill_manual(values = c("#088DA5", "#6abac9", "#cde8ed")) +
  scale_color_manual(values = c("#E69F00", "#F0E442", "#D55E00", "#CC79A7")) +
  ggnewscale::new_scale_fill() +
  geom_rect(aes(fill = chamber), color = "black",
            xmin = 0.5, xmax = 8.5, ymin = 3.5, ymax = 4.8,
            data = ~slice(., 1, .by = bed)) + 
  geom_text(aes(label = bed),
            x = 4.5, y = 4.2, color = "black",
            data = ~slice(., 1, .by = bed)) +
  scale_fill_manual(values = c("#f9e5e5", "#e06666")) +
  ggnewscale::new_scale_fill() +
  geom_rect(xmin = 8.5, xmax = 9.1, ymin = 0.5, ymax = 4.8,
            aes(fill = water), color = "black",
            data = ~filter(., 
                           bed2 == "bed2",
                           row == "row8",
                           col == "col2")) +
  geom_text(aes(label = shelter),
            x = 8.8, y = 2.6, angle = 270, size = 3.5,
            data = ~filter(., 
                           bed2 == "bed2",
                           row == "row8",
                           col == "col2")) +
  scale_fill_manual(values = c("#2986cc", "#94c2e5", "#e9f2f9")) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  facet_grid(shelter ~ bed2) +
  theme_void(base_size = 12) +
  theme(strip.background.x = element_blank(), 
        strip.text = element_blank(),
        legend.position = "bottom",
        legend.direction = "vertical",
        legend.title = element_text(family = "mono", face = "bold", margin = margin(r = 20)),
        axis.text = element_text(color = "black"),
        axis.text.y = element_text(margin = margin(r = unit(1, "lines"))),
        axis.text.x = element_text(margin = margin(t = unit(1, "lines"))),
        panel.background = element_rect(color = 'black'),
        panel.spacing.y = unit(1.4, "lines"),
        plot.margin = margin(t = unit(14, "lines"), b = unit(1, "lines"), r = unit(30, "lines"))) +
  coord_cartesian(clip = 'off') 
```

The design above does not necessarily ensure that "the same 12 populations always occurred together in a set of four rows". This split by four rows sounded like the bed was blocked into two parts across the rows. Below, we create a block factor with two levels in each bed, then the block is assigned to rows in a systematic manner. 

```{r}
garden_units2 <- garden_units1 %>% 
  set_units(block = nested_in(bed, 2)) %>% 
  allot_units(block ~ row) %>% 
  assign_units(order = "systematic-slowest")
```

Table A1 in the supplementary material in @cochraneClimateWarmingDelays2015 provides a list of the population IDs. Upon inspection of these IDs, it became clear that they were encoded by species, rainfall cline, and replicate number (e.g., "Bm H1", "Bm H2", and "Bm L1" are population of species _Banksia media_ from site with mean annual precipitation of 574, 557, and 301, respectively). There were two replicates for each species and the rainfall cline.

```{r}
garden_trts2 <- garden_trts1 %>% 
  set_trts(rep = 1:2)
```

Although not specified in the description, the replicate number was likely allotted to the blocks. The population is then determined from the species, rainfall, and replication number assigned. This assignment ensures that "the same 12 populations always occurred together in a set of four rows". The resulting design, shown in @fig-garden-design2, now aligns with all of the text description.

```{r}
garden_design2 <- (garden_units2 + garden_trts2) %>%
  allot_table(water ~ shelter,
              chamber ~ bed,
              species ~ row,
              rainfall ~ col,
              rep ~ block, # The missing allotment in garden_design1
              label_nested = c(col, row),
              seed = 2023) %>% 
  dplyr::mutate(population = paste0("B", substr(species, 1, 1), " ", 
                                    substr(rainfall, 1, 1), rep))
```


```{r fig-garden-design2}
#| echo: false
#| fig-height: 10
#| fig-width: 8
#| fig-cap: "The generated experimental design in `garden design2`. The text label in the tile correpond to the population ID in the supplementary material in @cochraneClimateWarmingDelays2015."
garden_design2 %>% 
  mutate(bed2 = attr(bed, "label-nested"))  %>% 
  ggplot(aes(row, col)) +
  geom_tile(aes(fill = chamber)) +
  geom_text(aes(label = population), size = 2) +
  geom_text(aes(label = bed),
             x = 4.6, y = 4, color = "black",
            data = ~filter(., row == "row4", col == "col3")) +
  scale_fill_manual(values = c("#f9e5e5", "#e06666")) +
  ggnewscale::new_scale_fill() +
  geom_rect(xmin = 8.6, xmax = 9.2, ymin = -0.2, ymax = 4.2,
            aes(fill = water), color = "black",
            data = ~filter(., 
                           bed2 == "bed2",
                           row == "row8",
                           col == "col2")) +
  geom_text(aes(label = shelter),
            x = 8.9, angle = 270,
            data = ~filter(., 
                           bed2 == "bed2",
                           row == "row8",
                           col == "col2")) +
  scale_fill_manual(values = c("#2986cc", "#94c2e5", "#e9f2f9")) +
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  facet_grid(shelter ~ bed2) +
  theme_void(base_size = 12) +
  theme(strip.background.x = element_blank(), 
        strip.text = element_blank(),
        legend.position = "bottom",
        legend.title = element_text(family = "mono", face = "bold", margin = margin(l = 20)),
        axis.text = element_text(color = "black"),
        axis.text.y = element_text(margin = margin(r = unit(1, "lines"))),
        axis.text.x = element_text(margin = margin(t = unit(1, "lines"))),
        #strip.background.y = element_rect(color = "black"),
        panel.background = element_rect(color = 'black'),
        panel.spacing.y = unit(1.4, "lines"),
        #strip.text.y = element_text(angle = 270, margin = margin(2, 2, 2, 2)),
        plot.margin = margin(t = unit(14, "lines"), b = unit(1, "lines"), r = unit(30, "lines"))) +
        #axis.text.x = element_text(angle = 30, hjust = 1, vjust = 0)) +
  #expand_limits(y = c(1, 4.8))  +
  coord_cartesian(clip = 'off') 
```



## Evaluation of Japanese compositions  {#sec-composition}

@williamsExperimentalDesignPractice2021a describe a study on the evaluation of Japanese language compositions written by Japanese beginners by assessors of different backgrounds. The background of this experiment is given in @imakiAnatomyCompositionHow2014, with a reduced summary provided below.

::: callout-note

## Example 3: Evaluation of Japanese composition


Ten compositions (C1–C10) were conditionally drawn from the Taiyaku Database of the National Institute of the Japanese Language on a topic about "special event in my country". The compositions were selected from writers who were Japanese beginners (studied Japanese for less than 300 hours), with a final selection consisting of three compositions each from Singapore, Brazil, Finland, and a single composition from Cambodia. The selected compositions have an average of 600 characters. 

The study investigated the ratings from people of three different backgrounds: native Japanese-language teachers (NT), non-native Japanese language teachers (NNT), and native Japanese who are not language teachers (NG). While it was intended to involve 10 raters from each category, availability resulted in 10, 8, and 11 raters, respectively. 

Each rater was asked to evaluate the 10 compositions in four aspects (accuracy, structure and form, context, and richness) plus an overall score. Each rater was asked to rate the composition in a particular order (but different order across raters). The raters were given the option of taking breaks during the assessment process.

Possible carry-over effect was accounted for by using a set of three different $10 \times 10$ Williams Latin squares design [@williamsExperimentalDesignsBalanced1949a] where the columns represent the raters and the rows represents the order of composition presented. As there were 29 raters, one column was not used.  

:::

From the description of the experiment, I inferred the following structure.

```{r}
composition_str <- design("Japanese composition") %>% 
  set_units(background = c("NT", "NNT", "NG"),
            rater = nested_in(background, 
                              "NT" ~ 10,
                              "NNT" ~ 8,
                              "NG" ~ 11),
            order = 1:10,
            assess = crossed_by(rater, order)) %>% 
  set_trts(composition = paste0("C", 1:10)) %>% 
  set_rcrds_of(assess = c("accuracy", "structure", "context", 
                          "richness", "overall")) %>% 
  allot_trts(composition ~ assess)
```

The algorithms in [edibble]{.pkg} does not contain an algorithm to generate a Williams Latin square design [@williamsExperimentalDesignsBalanced1949a], thus we require a custom treatment assignment order algorithm. The [crossdes]{.pkg} package [@crossdes] has the function [williams]{.fct} to generate the Williams Latin squares design; we leverage this function in the algorithm. 

In [edibble]{.pkg}, a custom treatment assignment order algorithm is implemented as an S3 generic function, [order_trts]{.fct}, within [assign_trts]{.fct}. A new algorithm is created by writing an S3 method as a function with an argument that takes in the name of the algorithm, treatments table, units table, any constraint on the unit that the treatment is applied, and other objects (such as the provenance object that contains the entire graph structure). Suppose we name our new algorithm as "williams" then we can write the S3 method [order_trts.williams]{.fct} as follows. 

```{r}
order_trts.williams <- function(name, trts_table, units_table, constrain, ...) {
  # the experimental unit must be a unit crossed by two unit factors
  stopifnot(length(constrain)==2)
  ntrts <- nrow(trts_table)
  
  # assumes unit is ordered correctly
  units <- lapply(units_table[as.character(constrain)], unique) 
  nunits <- sapply(units, length)
  
  # smaller unit is assumed to be row, larger unit is the column
  unit_row <- as.character(constrain[order(nunits)==1])
  unit_col <- as.character(constrain[order(nunits)==2])
  
  # the row unit must be exactly equal to the number of treatments
  stopifnot(sort(nunits)[1] == ntrts)
  
  # generate the William's Latin square design
  williams_design <- t(crossdes::williams(ntrts))
  
  # number of times to replicate William's design
  nwilliams <- ceiling(sort(nunits)[2]/ntrts)
  
  # randomise treatment for each replicate
  out <- do.call(rbind, lapply(1:nwilliams, function(i) {
    data.frame(..trt.. = sample(1:ntrts)[as.vector(williams_design)],
               ..row.. = as.vector(row(williams_design)),
               ..col.. = as.vector(col(williams_design)) + (i - 1) * ntrts)
  })) %>% 
    subset(..col.. <= sort(nunits)[2])
  
  # translate row and column to their unit ids
  out[[unit_row]] <- units[[unit_row]][out$..row..]
  out[[unit_col]] <- units[[unit_col]][out$..col..]
  
  # combine and return the allocated treatment id vector 
  # ensuring the treatment order matches the units_table
  units_table$..id.. <- 1:nrow(units_table)
  out <- merge(units_table, out)
  out[order(out$..id..), "..trt..", drop = TRUE]
}
```

The above algorithm is used by setting the argument `order = "williams"`. The above function contains checks to ensure the experimental structure is as expected before proceeding to the allocation. This check makes it clear what the expected experimental structure is for the selected ordering algorithm. The resulting design is illustrated in @fig-composition.

```{r}
composition_design <- composition_str %>% 
  assign_trts(order = "williams", seed = 2023) %>% 
  serve_table()
```


```{r fig-composition}
#| echo: false
#| fig-height: 6
#| fig-cap: "An experimental design for the evaluation of Japanese compositions. The plot shows the order in which the ten composition was evaluated by each rater with the rater grouped by their background. "
composition_design %>% 
  mutate(order = factor(order),
         rater = factor(rater)) %>% 
  ggplot(aes(order, rater)) +
  geom_tile(aes(fill = composition), color = "black") +
  geom_text(aes(label = composition)) +
  theme_void(base_size = 12) +
  theme(axis.text = element_text(color = 'black')) +
  facet_wrap(~background, scales = "free_y", ncol = 1) + 
  ggh4x::force_panelsizes(rows = c(11, 8, 10)) +
  guides(fill = "none") +
  scale_fill_manual(values = c("#77AADD", "#99DDFF", "#44BB99", "#BBCC33", "#AAAA00", "#EEDD88", "#EE8866", "#FFAABB", "#DDDDDD", "white")) +
  theme(axis.title.x = element_text(color = "black"))
```



## Menu Designs 

While selecting a design from a set of "named" or "recipe" designs are not encouraged, there are cases when it is convenient to refer to these. In [edibble]{.pkg},  functions that prefix with `menu_` generate the full recipe code in their basic terms. For example, a completely randomised design (CRD) can be generated as:

```{r}
menu_crd(n = 10, t = 2, seed = 1)
```

If no argument is supplied, then the code automatically generates a random number:

```{r}
menu_crd()
```

The above output is a code print out designed so that users can copy and paste the output and easily modify the numbers or factor names to their context. If you would like to actually run the code from the menu output, then the [takeout]{.fct} function allows you to do this easily: 

```{r}
takeout(menu_crd(n = 4, t = 2))
```

If no argument is supplied for [takeout]{.fct}, a random recipe design is selected from the menu. The full menu can be seen from [scan_menu]{.fct}.

```{r}
#| message: true
scan_menu()
```




# Comparisons with other packages {#sec-contrast}

There are currently `r n_distinct(ctv:::.get_pkgs_from_ctv_or_repos("ExperimentalDesign", repos = "https://cran.rstudio.org")[[1]])` R-packages in the CRAN Task View of Design of Experiments and Analysis of Experimental Data (as of `r Sys.Date()`). The list is a comprehensive, although not exhaustive,  list of R-packages that can construct experimental designs, and some included in the list have an analytical focus or serve as a repository of experimental data [@tanakaCurrentStateProspects2022a]. 

The packages used to construct experimental designs can be broadly categorised into three approaches: specialist, recipe or general purpose. The first type is exclusive to a domain or has a limited experimental structure. The second often makes use of named experimental designs, e.g., the [agricolae]{.pkg} package [@agricolae] specifies a Latin square design and a split-plot design by the function [design.lsd]{.fct} and [design.split]{.fct}, respectively. Alternatively, a package can contain a catalogue of pre-optimised designs, e.g., [FrF2]{.pkg} [@FrF2]. The general purpose approach can construct a more flexible range of experimental designs using either a randomisation or optimisation methods. E.g. [AlgDesign]{.pkg} package [@AlgDesign] includes a [optFederov]{.fct} function finds an optimal design using Federov's exchange algorithm under certain optimal criterion. 



A specialist or recipe approach clearly does not easily allow users to deviate away from a named experimental design. A general purpose approach seems to offer the most flexibility but it often assumes a *mise en place*; in other words, the ingredients are prepared before "cooking".  A prime example of this is optimal or model-based design, such as in [AlgDesign]{.pkg}, which requires the initial data structure and specification of the algorithm type and/or model. There are often helper functions to obtain the required initial data structure (such as [gen.factorial]{.fct} in [AlgDesign]{.pkg}).

All of the above mentioned packages have an emphasis on the algorithm to generate the experimental design. On the other hand, [edibble]{.pkg} package serves as a holistic approach in which segments of the process are modularised. The algorithms, catalogues, or recipes can be integrated into the system as shown in @sec-composition; in this sense, the system is largely complementary to existing systems and offers fundamentally different functionalities. A prime example is the downstream usage of record factors to enter data validation rules (see @sec-rcrds), which are not offered in the other experimental design packages. 


# Discussion {#sec-conclude}



This paper showcases the capabilities of the [edibble]{.pkg} package, an implementation of the grammar of experimental designs system by @tanakaUnifiedLanguageExperimental2023 in the [R]{.proglang} language. While it is conceivable to implement the system in other programming languages, the concrete implementation as an [R]{.proglang} package should aid in a faster adoption by the scientific community to make their experiments transparent.  

The demonstrated downstream benefits of [edibble]{.pkg} include the ability to specify the record (including responses) of unit factors ([set_rcrds]{.fct}), and the expected values these records ([expect_rcrds]{.fct}) can take, which in turn can be encoded as data validation rules in the exported data ([export_design]{.fct}). These functionalities are also leveraged in the simulation of record factors ([simulate_rcrds]{.fct}), including in a lazy, automated fashion ([autofill_rcrds]{.fct}). These functionalities aid in better planning and management, and consequently, in the workflow for constructing the experimental design. 


Software packages for experimental designs often focus on the algorithmic components of the design. In contrast, [edibble]{.pkg} package reframes the specification of the experimental design by conceiving experimental designs as a mutable object that can be progressively built. This reframing allows separation in the process of specifying the experimental structure and the assignment algorithm. The benefit of this approach is that  users can easily consider various assignment algorithms without respecifying their structure. This approach also aligns with the common mantra that a good experimental design considers the structure first rather than just selecting a recipe design from the menu  [@baileyDesignComparativeExperiments2008]. 

@meadStatisticalPrinciplesDesign2012 attest Procrustean designs that compromise the experimental objective are still common. The [edibble]{.pkg} package takes a holistic approach that aims to encourage users to think of the experimental structure. By engaging user cognition, this can potentially translate into a better practice of experimental design. 

The [edibble]{.pkg} package can also aid in the teaching of experimental design. @smuckerProfilesTeachingExperimental2023 found from their survey of experimental design courses that nearly all cover randomisation, replication, blocking design principles, and multi-factor experiments, and the use of software for both instruction and assessment is ubiquitous. By adhering closely to the principles of [tidyverse]{.pkg} [@wickhamWelcomeTidyverse2019], which has widespread adoption in modern day practice and teaching [@cetinkaya-rundelEducatorPerspectiveTidyverse2022], [edibble]{.pkg} can integrate more easily in the workflow of everyday users or learners by leveraging their knowledge in a familiar system.


The [edibble]{.pkg} package is fairly flexible in constructing various fixed experimental structures, and it is fairly easy to wrap existing experimental design algorithms into the system, as shown in @sec-composition. The ability to incorporate other algorithms into the system provides an opportunity for users to select the most appropriate algorithm for their structure. However, there is no guarantee that the generated experimental design will be appropriate or optimal. Afterall, a tool is only as effective as the skillfulness of the hands that wield it. The [edibble]{.pkg} package aims to cognitively engage the user with the hope that possible broader concerns in the experimental design are detected prior to the execution of the experiment, but ultimately, the user should ensure to perform their own diagnostic checks that the experimental design output is appropriate for their own objectives.



# Computational details {.unnumbered}

The paper was written using [Quarto]{.proglang} version `r quarto::quarto_version()` with [knitr]{.pkg} version `r packageVersion("knitr")`, Pandoc version `r rmarkdown::pandoc_version()` and [R]{.proglang} `r gsub("R ", "", version$version.string)`. Some of the figures were drawn using [ggplot2]{.pkg} version `r packageVersion("ggplot2")`. The code uses [edibble]{.pkg} version `r packageVersion("edibble")`.


# Acknowledgements {.unnumbered}

I thank Francis Hui and many others that supported the idea and provided feedback on the system. 

# References {.unnumbered}




