---
title: "gexp: package overview"
author: "Allaman, I. B. <br> Faria, J. C."
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{gexp: package overview}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ../inst/biblio/bibliography.bib
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```
The `gexp` package builds on earlier functions developed by Faria, J. C. Its goal is to generate one or **p** response variables for a given experimental structure, with treatment effects specified by the user.

Packages such as **agricolae**, for example, focus mainly on generating experimental layouts. The `gexp` package goes further: it simulates responses under a specified design while letting the user control which treatments differ.

Supported designs include the completely randomized design (CRD), randomized complete block design (RCBD), Latin square design (LSD), factorial experiments under CRD, RCBD, and LSD, and split-plot arrangements under the same base designs [@Aquino:1992].

Function arguments follow a linear model in matrix form (@Ferreira:2008, @Rencher:2007):

$$Y = X\beta + E \hspace{5cm} (1)$$

where:

$Y$ is an $n \times p$ matrix with **n** observations on **p** response variables;

$X$ is an $n \times m$ model matrix, where **m** is the number of parameters affecting the response ($\mu$, $\alpha$, $\beta$, etc.);

$\beta$ is an $m \times p$ matrix of factor effects;

$E$ is an $n \times p$ matrix of errors, multivariate normal with mean **0** and covariance **$\Sigma$**;

For split-plot experiments, the package uses the following linear mixed model [@Naes:2007]:

$$Y = X\beta + Zu + E \hspace{5cm} (2)$$

where:

$Z$ is an $n \times q$ plot-error model matrix (**q** depends on the design: replication in CRD, block in RCBD, row:column in LSD);

$u$ is a $q \times p$ matrix of whole-plot errors, multivariate normal with mean **0** and covariance **$\Sigma_{plot}$**;

$E$ is an $n \times p$ matrix of subplot errors, multivariate normal with mean **0** and covariance **$\Sigma_{sub-plot}$**;

The remaining terms in equation (2) are defined as in equation (1).

Notation used in the package:

<center>
<table>
<tr class='head'>
<th>Sources of Variation</th>
<th>Symbol</th>
</tr>
<tr>
<td>Treatments</td>
<td>$\alpha,\tau$</td>
</tr>
<tr>
<td>Block</td>
<td>$\beta$</td>
</tr>
<tr>
<td>Interaction</td>
<td>$\gamma$</td>
</tr>
<tr>
<td>Row block</td>
<td>$\eta$</td>
</tr>
<tr class='bottom'>
<td>Column block</td>
<td>$\delta$</td>
</tr>
</table>
</center>

The **multivariate** case is illustrated with CRD only; other designs follow the same logic. 

# Simulating a Completely Randomized Design (CRD)
## Univariate case

Suppose we want to generate a **single** random variable under a CRD with a factor and two levels with three replicates. Then we can denote $Y_{ik}$ as the random variable observed in the k-th experimental unit (k = 1, 2, 3) that received the i-th level of factor X1 (i = 1,2). Also, assuming that this variable is under the effect of a constant ($\mu = 15$) and each level of the factor affects the response: $\alpha_1 = 1$ and $\alpha_2 = -2$ the response, respectively. The matrix layout is:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{12} \\
Y_{22} \\
Y_{13} \\
Y_{23}
\end{bmatrix}
&=&
\begin{array}{p{5cm}}
\begin{matrix}
\mu & \!\scriptsize{x11}\! & \scriptsize{x12}
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0} & 1\hphantom{0}  & 0 \\
1\hphantom{0} & 0\hphantom{0}  & 1 \\
1\hphantom{0} & 1\hphantom{0}  & 0 \\
1\hphantom{0} & 0\hphantom{0}  & 1 \\
1\hphantom{0} & 1\hphantom{0}  & 0 \\
1\hphantom{0} & 0\hphantom{0}  & 1
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\mu \\
\alpha_{1}  \\
\alpha_{2}  \\
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{12} \\
e_{22} \\
e_{13} \\
e_{23}
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{x11}\! & \scriptsize{x12}
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0} & 1\hphantom{0}  & 0 \\
1\hphantom{0} & 0\hphantom{0}  & 1 \\
1\hphantom{0} & 1\hphantom{0}  & 0 \\
1\hphantom{0} & 0\hphantom{0}  & 1 \\
1\hphantom{0} & 1\hphantom{0}  & 0 \\
1\hphantom{0} & 0\hphantom{0}  & 1
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
15 \\
1  \\
-2  \\
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{12} \\
e_{22} \\
e_{13} \\
e_{23}
\end{bmatrix} 
\\
&=&
\begin{bmatrix}
16 \\
13 \\
16 \\
13 \\
16 \\
13
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{12} \\
e_{22} \\
e_{13} \\
e_{23}
\end{bmatrix} 
\end{eqnarray}
$$ 

With zero mean and zero variance for the error, the simulated responses are:

$$
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{12} \\
Y_{22} \\
Y_{13} \\
Y_{23}
\end{bmatrix}
=
\begin{bmatrix}
16 \\
13 \\
16 \\
13 \\
16 \\
13
\end{bmatrix} 
$$

With the `gexp` package:

```{r}
library(gexp)
crd <- gexp(mu = 15,
            err = matrix(rep(0, 6),
                         nrow = 6),
            r = 3,
            fe = list(alpha = c(1, -2)))
summary(crd)
```

## Multivariate Case

Let us now assume that we want to generate **two** random variables by considering the same number of factors and the same number of repetitions of the previous case. Then we can denote $Y_{ikl}$ as the lth random variable (l = 1,2) observed in the k-th experimental unit (k = 1,2,3) that received the i-th factor level X1 (i = 1,2). Also, assuming that this variable is under the effect of a constant ($\mu_1 = 15$ and $\mu_2 = 6$) and each level of the factor affects the response: $\alpha_{11} = 1$ and $\alpha_{21} = -2$ the variable 1 and, $\alpha_{12} = 2$ and $\alpha_{22} = 3$ the variable 2, we will have the following matrix configuration:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{111} & Y_{112} \\
Y_{211} & Y_{212} \\
Y_{121} & Y_{122} \\
Y_{221} & Y_{222} \\
Y_{131} & Y_{132} \\
Y_{231} & Y_{232}
\end{bmatrix}
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{x11}\! & \scriptsize{x12}
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0}  & 1\hphantom{0} & 0 \\
1\hphantom{0}  & 0\hphantom{0} & 1 \\
1\hphantom{0}  & 1\hphantom{0} & 0 \\
1\hphantom{0}  & 0\hphantom{0} & 1 \\
1\hphantom{0}  & 1\hphantom{0} & 0 \\
1\hphantom{0}  & 0\hphantom{0} & 1
\end{array}\right]
\end{array}  
\cdot
\begin{bmatrix}
\mu_1       & \mu_2 \\
\alpha_{11} & \alpha_{12} \\
\alpha_{21} & \alpha_{22} \\
\end{bmatrix}
+
\begin{bmatrix}
e_{111} & e_{112} \\
e_{211} & e_{212} \\
e_{121} & e_{122} \\
e_{221} & e_{222} \\
e_{131} & e_{132} \\
e_{231} & e_{232}
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{x11}\! & \scriptsize{x12}
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0}  & 1\hphantom{0}  & 0 \\
1\hphantom{0}  & 0\hphantom{0}  & 1 \\
1\hphantom{0}  & 1\hphantom{0}  & 0 \\
1\hphantom{0}  & 0\hphantom{0}  & 1 \\
1\hphantom{0}  & 1\hphantom{0}  & 0 \\
1\hphantom{0}  & 0\hphantom{0}  & 1
\end{array}\right]
\end{array}   
\cdot
\begin{bmatrix}
15 & 6 \\
1  & 2 \\
-2 & 3 \\
\end{bmatrix}
+
\begin{bmatrix}
e_{111} & e_{112} \\
e_{211} & e_{212} \\
e_{121} & e_{122} \\
e_{221} & e_{222} \\
e_{131} & e_{132} \\
e_{231} & e_{232}
\end{bmatrix} 
\\
&=&
\begin{bmatrix}
16 & 8 \\
13 & 9 \\
16 & 8 \\
13 & 9 \\
16 & 8 \\
13 & 9
\end{bmatrix}
+
\begin{bmatrix}
e_{111} & e_{112} \\
e_{211} & e_{212} \\
e_{121} & e_{122} \\
e_{221} & e_{222} \\
e_{131} & e_{132} \\
e_{231} & e_{232}
\end{bmatrix}  
\end{eqnarray}
$$

With zero mean errors and a zero variance-covariance matrix, the simulated responses are:

$$
\begin{bmatrix}
Y_{111} & Y_{112} \\
Y_{211} & Y_{212} \\
Y_{121} & Y_{122} \\
Y_{221} & Y_{222} \\
Y_{131} & Y_{132} \\
Y_{231} & Y_{232}
\end{bmatrix}
=
\begin{bmatrix}
16 & 8 \\
13 & 9 \\
16 & 8 \\
13 & 9 \\
16 & 8 \\
13 & 9
\end{bmatrix} 
$$

With the `gexp` package:

```{r}
crd_mv <- gexp(mu = c(15, 6),
               err = matrix(c(0, 0,
                              0, 0,
                              0, 0,
                              0, 0,
                              0, 0,
                              0, 0),
                            ncol = 2),
               r = 3L,
               fe = list(alpha = matrix(c(1, -2,
                                          2, 3),
                                         ncol = 2)))
summary(crd_mv)
```

# Simulating a Randomized Completely Block Design (RCBD) 

Assuming we want to generate a **single** random variable under an RCBD with a factor and two levels being a block with three levels. Then we can denote $Y_ {ij}$ as the random variable observed in the j-th block (j = 1, 2, 3) that received the i-th level of factor X1 (i = 1, 2). Also, assuming that this variable is under the effect of a constant ($\mu = 15$) and each level of the factor affects the response: $\alpha_1 = 1$ and $\alpha_2 = -2$ the response, respectively, and each block effect is $\beta_1 = 2, \beta_2 = 4$ and $\beta_3 = 6$ respectively, we will have the following matrix configuration:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{12} \\
Y_{22} \\
Y_{13} \\
Y_{23}
\end{bmatrix}
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{x11}\!&\scriptsize{x12}\hphantom{0}&\scriptsize{b1}&\scriptsize{b2}&\scriptsize{b3}
\end{matrix}\\
\left[\begin{array}{c|cc|ccc} 
1 & 1\hphantom{0}& 0\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 0\\
1 & 0\hphantom{0}& 1\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 0\\
1 & 1\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 1\hphantom{0}& 0\\
1 & 0\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 1\hphantom{0}& 0\\
1 & 1\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 1\\
1 & 0\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 1\\
\end{array}\right]
\end{array}    
\cdot
\begin{bmatrix}
\mu \\
\alpha_1 \\
\alpha_2 \\
\beta_1  \\
\beta_2  \\
\beta_3  \\
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{12} \\
e_{22} \\
e_{13} \\
e_{23}
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{x11}\!&\scriptsize{x12}\hphantom{0}&\scriptsize{b1}&\scriptsize{b2}&\scriptsize{b3}
\end{matrix}\\
\left[\begin{array}{c|cc|ccc} 
1 & 1\hphantom{0}& 0\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 0\\
1 & 0\hphantom{0}& 1\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 0\\
1 & 1\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 1\hphantom{0}& 0\\
1 & 0\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 1\hphantom{0}& 0\\
1 & 1\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 1\\
1 & 0\hphantom{0}& 1\hphantom{0}& 0\hphantom{0}& 0\hphantom{0}& 1\\
\end{array}\right]
\end{array}     
\cdot
\begin{bmatrix}
15 \\
1  \\
-2  \\
2  \\
4  \\
6
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{12} \\
e_{22} \\
e_{13} \\
e_{23}
\end{bmatrix} 
\\
&=&
\begin{bmatrix}
18 \\
15 \\
20 \\
17 \\
22 \\
19
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{12} \\
e_{22} \\
e_{13} \\
e_{23}
\end{bmatrix} 
\end{eqnarray}
$$ 

With zero-mean errors:

$$
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{12} \\
Y_{22} \\
Y_{13} \\
Y_{23}
\end{bmatrix}
=
\begin{bmatrix}
18 \\
15 \\
20 \\
17 \\
22 \\
19
\end{bmatrix}
$$

With the `gexp` package:

```{r}
rcbd <- gexp(mu = 15,
             r = 1,
             err = matrix(rep(0, 6),
                          nrow = 6),
             fe = list(alpha = c(1, -2)),
             blke = c(2, 4, 6),
             design = 'RCBD')
summary(rcbd)
```

# Simulating a Latin Square Design (LSD)

Assuming we want to generate a **single** random variable under an LSD with a factor and three levels, with a blocking factor in the row sense with three levels and another blocking factor in the column sense with three levels. Then we can denote $Y_{rci}$ as the random variable observed in the c-th column (c = 1, 2, 3) and r-th row (r = 1, 2, 3) that received the i-th level of factor X1 (i = 1, 2, 3). Also, assuming this variable is under the effect of a constant ($\mu = 15$) and each level of the factor affects the response: $\alpha_1 = 1, \alpha_2 = -2$, and $\alpha_3 = 3$ for the response. Row effects are $\eta_1 = 2$, $\eta_2 = 3$, and $\eta_3 = 4$; column effects are $\delta_1 = 6$, $\delta_2 = 7$, and $\delta_3 = 8$. The matrix layout is:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{111} \\
Y_{122} \\
Y_{133} \\
Y_{212} \\
Y_{223} \\
Y_{231} \\
Y_{313} \\
Y_{321} \\
Y_{332} 
\end{bmatrix}
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{r1}\!&\scriptsize{r2}&\scriptsize{r3}&\scriptsize{c1}&\scriptsize{c2}&\scriptsize{c3}&\scriptsize{x11}&\scriptsize{x12}&\scriptsize{x13}
\end{matrix}\\
\left[\begin{array}{c|ccc|ccc|ccc}  
1&1&0&0\hphantom{0}&1&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\\
1&1&0&0\hphantom{0}&0&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\\
1&1&0&0\hphantom{0}&0&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\\
1&0&1&0\hphantom{0}&1&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\\
1&0&1&0\hphantom{0}&0&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\\
1&0&1&0\hphantom{0}&0&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\\
1&0&0&1\hphantom{0}&1&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\\
1&0&0&1\hphantom{0}&0&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\\
1&0&0&1\hphantom{0}&0&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\\
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\mu \\
\eta_1 \\
\eta_2 \\
\eta_3  \\
\delta_1  \\
\delta_2  \\
\delta_3 \\
\alpha_1 \\
\alpha_2 \\
\alpha_3 \\
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{122} \\
e_{133} \\
e_{212} \\
e_{223} \\
e_{231} \\
e_{313} \\
e_{321} \\
e_{332}
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix}
\mu & \!\scriptsize{r1}\!&\scriptsize{r2}&\scriptsize{r3}&\scriptsize{c1}&\scriptsize{c2}&\scriptsize{c3}&\scriptsize{x11}&\scriptsize{x12}&\scriptsize{x13}
\end{matrix}\\
\left[\begin{array}{c|ccc|ccc|ccc}  
1&1&0&0\hphantom{0}&1&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\\
1&1&0&0\hphantom{0}&0&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\\
1&1&0&0\hphantom{0}&0&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\\
1&0&1&0\hphantom{0}&1&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\\
1&0&1&0\hphantom{0}&0&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\\
1&0&1&0\hphantom{0}&0&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\\
1&0&0&1\hphantom{0}&1&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\\
1&0&0&1\hphantom{0}&0&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\\
1&0&0&1\hphantom{0}&0&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\\
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
15 \\
2 \\
3 \\
4  \\
6  \\
7  \\
8 \\
1 \\
-2 \\
3 \\
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{122} \\
e_{133} \\
e_{212} \\
e_{223} \\
e_{231} \\
e_{313} \\
e_{321} \\
e_{332}
\end{bmatrix}
\\
&=&
\begin{bmatrix}
24 \\
22 \\
28 \\
22  \\
28  \\
27  \\
28 \\
27 \\
25 \\
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{122} \\
e_{133} \\
e_{212} \\
e_{223} \\
e_{231} \\
e_{313} \\
e_{321} \\
e_{332}
\end{bmatrix}
\end{eqnarray}
$$ 

With zero-mean errors:

$$
\begin{bmatrix}
Y_{111} \\
Y_{122} \\
Y_{133} \\
Y_{212} \\
Y_{223} \\
Y_{231} \\
Y_{313} \\
Y_{321} \\
Y_{332} 
\end{bmatrix}
=
\begin{bmatrix}
24 \\
22 \\
28 \\
22  \\
28  \\
27  \\
28 \\
27 \\
25 \\
\end{bmatrix} 
$$

With the `gexp` package:

```{r}
lsd <- gexp(mu = 15,
            r = 1,
            err = matrix(rep(0, 9),
                         nrow = 9),
            fe = list(alpha = c(1, -2, 3)),
            rowe = c(2, 3, 4),
            cole = c(6, 7, 8),
            design = 'LSD')
summary(lsd)
```

# Factorial Experiments
## CRD

Suppose we want to generate a **single** random variable under a CRD in a factorial scheme of type 2 x 3 (two levels of factor X1 and three levels of factor X2) with 2 replicates. Then we can denote $Y_{ijk}$ as the random variable observed in the k-th experimental unit (k = 1, 2) that received the i-th level of factor X1 (i = 1,2) and j-th level of factor X2 (j = 1,2,3). Also, assuming this variable is under the effect of a constant ($\mu = 15$), and each level of factor X1 gives $\alpha_1 = 1$ and $\alpha_2 = -2$, each level of factor X2 are $\tau_1 = 1, \tau_2 = -1$ and $\tau_3 = 1$ and the effect of the interaction is $\gamma_{11} = 3,\gamma_{21} = 1,\gamma_{12} = 1,\gamma_{22} = -5,\gamma_{13} = 1$, and $\gamma_{23} = 1$, we obtain the following matrix layout:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{111} \\
Y_{211} \\
Y_{121} \\
Y_{221} \\
Y_{131} \\
Y_{231} \\
Y_{112} \\
Y_{212} \\
Y_{122} \\
Y_{222} \\
Y_{132} \\
Y_{232} 
\end{bmatrix}
&=&
\begin{array}{c}
\begin{matrix} 
\hspace{0.2cm}\mu&\scriptsize{x11}&\scriptsize{x12}\hspace{0.1cm}&\scriptsize{x21}&\scriptsize{x22}&\scriptsize{x23}&\scriptsize{x11x21}&\scriptsize{x12x21}&\scriptsize{x11x22}&\scriptsize{x12x22}&\scriptsize{x11x23}&\scriptsize{x12x23}\hspace{0.6cm}
\end{matrix}\\
\left[\begin{array}{c|cc|ccc|cccccc} 
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{000}&1\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&1\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&1\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 1\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&1\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{000}&1\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&1\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&1\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 1\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&1\hphantom{000}\\  
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\mu \\
\alpha_1\\
\alpha_2\\
\tau_1\\
\tau_2\\
\tau_3\\
\gamma_{11}\\
\gamma_{21}\\
\gamma_{12}\\
\gamma_{22}\\
\gamma_{13}\\
\gamma_{23}
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{211} \\
e_{121} \\
e_{221} \\
e_{131} \\
e_{231} \\
e_{112} \\
e_{212} \\
e_{122} \\
e_{222} \\
e_{132} \\
e_{232}  
\end{bmatrix}
\\
&=& 
\begin{array}{c}
\begin{matrix} 
\hspace{0.2cm}\mu&\scriptsize{x11}&\scriptsize{x12}\hspace{0.1cm}&\scriptsize{x21}&\scriptsize{x22}&\scriptsize{x23}&\scriptsize{x11x21}&\scriptsize{x12x21}&\scriptsize{x11x22}&\scriptsize{x12x22}&\scriptsize{x11x23}&\scriptsize{x12x23}\hspace{0.6cm}
\end{matrix}\\
\left[\begin{array}{c|cc|ccc|cccccc} 
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{000}&1\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&1\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&1\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 1\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&1\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{000}&1\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&1\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&1\hphantom{0000}& 0\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 1\hphantom{0000}&0\hphantom{000}\\  
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{000}&0\hphantom{0000}&0\hphantom{0000}&0\hphantom{0000}& 0\hphantom{0000}&1\hphantom{000}\\  
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
15 \\
1\\
-2\\
1\\
-1\\
1\\
3\\
1\\
1\\
-5\\
1\\
1
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{211} \\
e_{121} \\
e_{221} \\
e_{131} \\
e_{231} \\
e_{112} \\
e_{212} \\
e_{122} \\
e_{222} \\
e_{132} \\
e_{232}  
\end{bmatrix}
\\
&=&
\begin{bmatrix}
20 \\
15\\
16\\
7\\
18\\
15\\
20\\
15\\
16\\
7\\
18\\
15
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{211} \\
e_{121} \\
e_{221} \\
e_{131} \\
e_{231} \\
e_{112} \\
e_{212} \\
e_{122} \\
e_{222} \\
e_{132} \\
e_{232}  
\end{bmatrix}
\end{eqnarray}
$$

With zero-mean errors:

$$
\begin{bmatrix}
Y_{111} \\
Y_{211} \\
Y_{121} \\
Y_{221} \\
Y_{131} \\
Y_{231} \\
Y_{112} \\
Y_{212} \\
Y_{122} \\
Y_{222} \\
Y_{132} \\
Y_{232} 
\end{bmatrix}
=
\begin{bmatrix}
20 \\
15\\
16\\
7\\
18\\
15\\
20\\
15\\
16\\
7\\
18\\
15
\end{bmatrix}
$$

With the `gexp` package:

```{r}
FE_crd <- gexp(mu = 15,
               r = 2,
               err = matrix(rep(0, 12),
                            nrow = 12),
               fe = list(alpha = c(1, -2),
                         tau = c(1, -1, 1)),
               inte = c(3, 1, 1,
                        -5, 1, 1),
               type = 'FE')
summary(FE_crd)
```

# Split-plot experiments
## CRD

Suppose we want a single response under a CRD split-plot with two whole-plot levels of factor X1 and two subplot levels of factor X2, with three replicates. Let $Y_{ijk}$ denote the response in the subplot of experimental unit $k$ ($k = 1, 2, 3$) receiving level $i$ of X1 ($i = 1, 2$) and level $j$ of X2 ($j = 1, 2$). With overall mean $\mu = 15$, whole-plot effects $\alpha_1 = 1$ and $\alpha_2 = -2$, subplot effects $\tau_1 = 1$ and $\tau_2 = -1$, and interaction effects $\gamma_{11} = 3$, $\gamma_{21} = 1$, $\gamma_{12} = 1$, and $\gamma_{22} = -5$, the setup follows equation (2).


The matrix layout is therefore:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{111} \\
Y_{211} \\
Y_{121} \\
Y_{221} \\
Y_{112} \\
Y_{212} \\
Y_{122} \\
Y_{222} \\
Y_{113} \\
Y_{213} \\
Y_{123} \\
Y_{223} 
\end{bmatrix}
&=&
\begin{array}{c}
\begin{matrix} 
\mu&\scriptsize{x11}&\scriptsize{x12}&\scriptsize{x21}&\scriptsize{x22}&\scriptsize{x11x21}&\scriptsize{x12x21}&\scriptsize{x11x22}&\scriptsize{x12x22}
\end{matrix}\\
\left[\begin{array}{c|cc|cc|cccc}                                                  
1\hphantom{0}& 1\hphantom{0} & 0\hphantom{0} & 1\hphantom{0} & 0\hphantom{00} & 1\hphantom{000} & 0\hphantom{000} & 0\hphantom{000} & 0\hphantom{00}\\               
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{00}&0\hphantom{000}&1\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&1\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&0\hphantom{000}&1\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{00}&1\hphantom{000}&0\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{00}&0\hphantom{000}&1\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&1\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&0\hphantom{000}&1\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{00}&1\hphantom{000}&0\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{00}&0\hphantom{000}&1\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&1\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&0\hphantom{000}&1\hphantom{00}\\   
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\mu \\
\alpha_1\\
\alpha_2\\
\tau_1\\
\tau_2\\
\gamma_{11}\\
\gamma_{21}\\
\gamma_{12}\\
\gamma_{22}
\end{bmatrix}
+
\begin{array}{c}
\begin{matrix} 
\scriptsize{x11r1}&\scriptsize{x12r1}\hspace{0.1cm}&\scriptsize{x11r2}&\scriptsize{x12r2}&\scriptsize{x11r3}&\scriptsize{x12r3}\hspace{0.6cm}
\end{matrix}\\
\left[\begin{array}{cccccc}                                   
1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 1\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 1\hphantom{0} \\ 
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
u_{11} \\
u_{21} \\
u_{12} \\
u_{22} \\
u_{13} \\
u_{23} 
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{211} \\
e_{121} \\
e_{221} \\
e_{112} \\
e_{212} \\
e_{122} \\
e_{222} \\
e_{113} \\
e_{213} \\
e_{123} \\
e_{223}  
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix} 
\mu&\scriptsize{x11}&\scriptsize{x12}&\scriptsize{x21}&\scriptsize{x22}&\scriptsize{x11x21}&\scriptsize{x12x21}&\scriptsize{x11x22}&\scriptsize{x12x22}
\end{matrix}\\
\left[\begin{array}{c|cc|cc|cccc}                                                  
1\hphantom{0}& 1\hphantom{0} & 0\hphantom{0} & 1\hphantom{0} & 0\hphantom{00} & 1\hphantom{000} & 0\hphantom{000} & 0\hphantom{000} & 0\hphantom{00}\\               
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{00}&0\hphantom{000}&1\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&1\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&0\hphantom{000}&1\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{00}&1\hphantom{000}&0\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{00}&0\hphantom{000}&1\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&1\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&0\hphantom{000}&1\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{00}&1\hphantom{000}&0\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&1\hphantom{0}&0\hphantom{00}&0\hphantom{000}&1\hphantom{000}&0\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&1\hphantom{0}&0\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&1\hphantom{000}&0\hphantom{00}\\   
1\hphantom{0}&0\hphantom{0}&1\hphantom{0}&0\hphantom{0}&1\hphantom{00}&0\hphantom{000}&0\hphantom{000}&0\hphantom{000}&1\hphantom{00}\\   
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
15 \\
1\\
-2\\
1\\
-1\\
3\\
1\\
1\\
-5
\end{bmatrix}
+
\begin{array}{c}
\begin{matrix} 
\scriptsize{x11r1}&\scriptsize{x12r1}\hspace{0.1cm}&\scriptsize{x11r2}&\scriptsize{x12r2}&\scriptsize{x11r3}&\scriptsize{x12r3}\hspace{0.6cm}
\end{matrix}\\
\left[\begin{array}{cccccc}                                   
1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 1\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 1\hphantom{0} \\ 
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
u_{11} \\
u_{21} \\
u_{12} \\
u_{22} \\
u_{13} \\
u_{23} 
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{211} \\
e_{121} \\
e_{221} \\
e_{112} \\
e_{212} \\
e_{122} \\
e_{222} \\
e_{113} \\
e_{213} \\
e_{123} \\
e_{223}  
\end{bmatrix}
\\
&=&
\begin{bmatrix}
20 \\
15 \\
16 \\
7  \\
20 \\
15 \\
16 \\
7 \\
20 \\
15 \\
16 \\
7  
\end{bmatrix} 
+
\begin{array}{c}
\begin{matrix} 
\scriptsize{x11r1}&\scriptsize{x12r1}\hspace{0.1cm}&\scriptsize{x11r2}&\scriptsize{x12r2}&\scriptsize{x11r3}&\scriptsize{x12r3}\hspace{0.6cm}
\end{matrix}\\
\left[\begin{array}{cccccc}                                   
1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{000}& 0\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 1\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 1\hphantom{00}& 0\hphantom{0} \\
0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{000}& 0\hphantom{00}& 1\hphantom{0} \\ 
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
u_{11} \\
u_{21} \\
u_{12} \\
u_{22} \\
u_{13} \\
u_{23} 
\end{bmatrix}
+
\begin{bmatrix}
e_{111} \\
e_{211} \\
e_{121} \\
e_{221} \\
e_{112} \\
e_{212} \\
e_{122} \\
e_{222} \\
e_{113} \\
e_{213} \\
e_{123} \\
e_{223}  
\end{bmatrix} 
\end{eqnarray}
$$

With zero-mean whole-plot and subplot errors:

$$
\begin{bmatrix}
Y_{111} \\
Y_{211} \\
Y_{121} \\
Y_{221} \\
Y_{112} \\
Y_{212} \\
Y_{122} \\
Y_{222} \\
Y_{113} \\
Y_{213} \\
Y_{123} \\
Y_{223} 
\end{bmatrix} 
=
\begin{bmatrix}
20 \\
15 \\
16 \\
7  \\
20 \\
15 \\
16 \\
7 \\
20 \\
15 \\
16 \\
7  
\end{bmatrix}  
$$

With the `gexp` package:

```{r}
SPE_crd <- gexp(mu = 15,
                err = matrix(rep(0, 12),
                             nrow = 12),
                errp = matrix(rep(0, 6),
                              nrow = 6),
                r = 3,
                fe = list(alpha = c(1, -2),
                          tau = c(1, -1)),
                inte = c(3, 1, 1, -5),
                type = 'SPE')
summary(SPE_crd)
```
# Simulating a CRD with a quantitative factor
## Non-Orthogonal Contrast
### Linear Effect

Suppose we want to generate a **single** random variable under a CRD with one factor and four levels (0, 5, 10, 15) with three replicates. Then we can denote $Y_ {ik}$ as the random variable observed in the k-th experimental unit (k = 1, 2, 3) that received the i-th level of factor X1 (i = 1,2,3,4). In this case, it is only possible to adjust a regression up to the third degree. As the treatment is quantitative, instead of indicating the effects of treatments as deviations from the general mean, we will indicate the value we want for each regression coefficient. So, instead of having $\mu, \alpha_1, \alpha_2, \alpha_3$ and $\alpha_4$, we have $\beta_0$ (the intercept), $\beta_1$ (the linear coefficient), $\beta_2$ the quadratic coefficient and $\beta_3$ (the cubic coefficient). Since the effect we want is linear, then we can provide the following values for betas: $\beta_0 = 2, \beta_1 = 3, \beta_2 = 0$ and $\beta_3 = 0$. Therefore, we will have the following matrix configuration:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{31} \\
Y_{41} \\
Y_{12} \\
Y_{22} \\
Y_{32} \\
Y_{42} \\
Y_{13} \\
Y_{23} \\
Y_{33} \\
Y_{43}
\end{bmatrix}
&=&
\begin{array}{p{5cm}}
\begin{matrix}
\beta_0 & \!\scriptsize{x11}\! & \scriptsize{x12} & \scriptsize{x13} 
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\ 
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\beta_{0} \\
\beta_{1} \\
\beta_{2} \\
\beta_{3}
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{31} \\
e_{41} \\
e_{12} \\
e_{22} \\
e_{23} \\
e_{24} \\
e_{13} \\
e_{23} \\
e_{33} \\
e_{43}
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix}
\beta_0 & \!\scriptsize{x11}\! & \scriptsize{x12} & \scriptsize{x13} 
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\  
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
2 \\
3  \\
0  \\
0  \\
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{31} \\
e_{41} \\
e_{12} \\
e_{22} \\
e_{23} \\
e_{24} \\
e_{13} \\
e_{23} \\
e_{33} \\
e_{43} 
\end{bmatrix} 
\\
&=&
\begin{bmatrix}
2 \\
17 \\
32 \\
47 \\
2 \\
17 \\
32 \\
47 \\
2 \\
17 \\
32 \\
47
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{31} \\
e_{41} \\
e_{12} \\
e_{22} \\
e_{23} \\
e_{24} \\
e_{13} \\
e_{23} \\
e_{33} \\
e_{43}  
\end{bmatrix} 
\end{eqnarray}
$$ 

With zero mean and zero variance for the error, the simulated responses are:

$$
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{31} \\
Y_{41} \\
Y_{12} \\
Y_{22} \\
Y_{32} \\
Y_{42} \\
Y_{13} \\
Y_{23} \\
Y_{33} \\
Y_{43} 
\end{bmatrix}
=
\begin{bmatrix}
2 \\
17 \\
32 \\
47 \\
2 \\
17 \\
32 \\
47 \\
2 \\
17 \\
32 \\
47 
\end{bmatrix} 
$$

With the `gexp` package:

```{r}
level <- seq(0, 15, 5)
cont_crd <- matrix(c(level,
                     level^2,
                     level^3),
                   ncol = 3)

crd_l <- gexp(mu = 2,
              r = 3,
              err = matrix(rep(0, 12),
                           nrow = 12),
              fe = list(f1 = c(3, 0, 0)),
              fl = list(dose = level),
              contrasts = list(dose = cont_crd))

summary(crd_l)
```

The plot below shows the simulated data.

```{r echo = FALSE, results = 'hide'}
with(crd_l$dfm,
     plot(Y1 ~ dose,
          axes = FALSE,
          ylim = c(0, 50),
          xlab = 'dose'))
axis(1, pos = 0)
axis(2, pos = 0)
abline(h = 2, col = 'red')
abline(2, 3, col = 'red', lty = 2)
text(-0.5, 2, '2', col = 'red', xpd = TRUE)
```

### Quadratic effect

Using the same example, we now add only a quadratic effect ($\beta_2 = 4$). The matrix layout is:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{11} \\
Y_{21} \\
Y_{31} \\
Y_{41} \\
Y_{12} \\
Y_{22} \\
Y_{32} \\
Y_{42} \\
Y_{13} \\
Y_{23} \\
Y_{33} \\
Y_{43}
\end{bmatrix}
&=&
\begin{array}{p{5cm}}
\begin{matrix}
\beta_0 & \!\scriptsize{x11}\! & \scriptsize{x12} & \scriptsize{x13} 
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\ 
\end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\beta_{0} \\
\beta_{1} \\
\beta_{2} \\
\beta_{3}
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{31} \\
e_{41} \\
e_{12} \\
e_{22} \\
e_{23} \\
e_{24} \\
e_{13} \\
e_{23} \\
e_{33} \\
e_{43}
\end{bmatrix}
\\
&=&
\begin{array}{c}
\begin{matrix}
\beta_0 & \!\scriptsize{x11}\! & \scriptsize{x12} & \scriptsize{x13} 
\end{matrix}\\
\left[\begin{array}{c|cc}
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\
1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  \\
1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  \\
1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 \\
1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 \\  
\end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
2 \\
3  \\
4  \\
0  \\
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{31} \\
e_{41} \\
e_{12} \\
e_{22} \\
e_{23} \\
e_{24} \\
e_{13} \\
e_{23} \\
e_{33} \\
e_{43} 
\end{bmatrix} 
\\
&=&
\begin{bmatrix}
2 \\
117 \\
432 \\
947 \\
2 \\
117 \\
432 \\
947 \\
2 \\
117 \\
432 \\
947
\end{bmatrix}
+
\begin{bmatrix}
e_{11} \\
e_{21} \\
e_{31} \\
e_{41} \\
e_{12} \\
e_{22} \\
e_{23} \\
e_{24} \\
e_{13} \\
e_{23} \\
e_{33} \\
e_{43}  
\end{bmatrix} 
\end{eqnarray}
$$ 

With the `gexp` package:

```{r}
crd_q <- gexp(mu = 2,
              r = 3,
              err = matrix(rep(0, 12),
                           nrow = 12),
              fe = list(f1 = c(3, 4, 0)),
              fl = list(dose = level),
              contrasts = list(dose = cont_crd))

summary(crd_q)
```

The plot below shows the simulated data.

```{r echo = FALSE, results = 'hide'}
with(crd_q$dfm,
     plot(Y1 ~ dose,
          axes = FALSE,
          xlab = 'dose'))
axis(1, pos = 0)
axis(2, pos = 0)
curve(2+3*x+4*x^2, col = 'red', lty = 2, add = TRUE)
text(-0.5, 2, '2', col = 'red', xpd = TRUE)
``` 

# Factorial CRD with two quantitative factors
## Linear response surface

Suppose we want to generate a **single** random variable under a CRD with two quantitative factors with levels 0, 5, 10, 15 and 2, 4, 6, 8 respectively, being 2 replicates. Then we can denote $Y_{ijk}$ as the random variable observed in the k-th experimental unit (k = 1, 2) that received the i-th level of factor X1 (i = 1,2,3,4) and j-th level of factor X2 (j = 1,2,3,4). In these cases, we should note that there is only one intercept. Therefore, such an intercept must appear in the first factor to be declared. In the others, the value of the intercept must be declared as 0 (only to compute the length of the vector correctly). Therefore, consider the following values for the first factor: $\beta_0 = 1, \beta_1 = 3, \beta_2 = 0, \beta_3 = 0$. For the second factor we have: $\beta_4 = 2, \beta_5 = 0, \beta_6 = 0$. In matrix form:

$$
\begin{eqnarray}
\begin{bmatrix}
Y_{111} \\
Y_{211} \\
Y_{311} \\
Y_{411} \\
Y_{121} \\
Y_{221} \\
Y_{321} \\
Y_{421} \\
Y_{131} \\
Y_{231} \\
Y_{331} \\
Y_{431} \\
Y_{141} \\
Y_{241} \\
Y_{341} \\
Y_{441} \\
Y_{112} \\
Y_{212} \\
Y_{312} \\
Y_{412} \\
Y_{122} \\
Y_{222} \\
Y_{322} \\
Y_{422} \\
Y_{132} \\
Y_{232} \\
Y_{332} \\
Y_{432} \\
Y_{142} \\
Y_{242} \\
Y_{342} \\
Y_{442} \\
\end{bmatrix}
&=&
  \begin{array}{p{5cm}}
\begin{matrix}
\beta_0 & \!\scriptsize{x11}\! & \scriptsize{x12} & \scriptsize{x13} & \!\scriptsize{x21}\! & \scriptsize{x22} & \scriptsize{x23} 
\end{matrix}\\
\left[\begin{array}{c|cc}
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 8\hphantom{0}  & 8^2 & 8^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 8\hphantom{0}  & 8^2 & 8^3 \\ 
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 8\hphantom{0}  & 8^2 & 8^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 8\hphantom{0}  & 8^2 & 8^3 \\  
      \end{array}\right]
\end{array}
\cdot
\begin{bmatrix}
\beta_{0} \\
\beta_{1} \\
\beta_{2} \\
\beta_{3} \\
\beta_{4} \\
\beta_{5} \\
\beta_{6}
\end{bmatrix}
+
  \begin{bmatrix}
e_{111} \\
e_{211} \\
e_{311} \\
e_{411} \\
e_{121} \\
e_{221} \\
e_{321} \\
e_{421} \\
e_{131} \\
e_{231} \\
e_{331} \\
e_{431} \\
e_{141} \\
e_{241} \\
e_{341} \\
e_{441} \\
e_{112} \\
e_{212} \\
e_{312} \\
e_{412} \\
e_{122} \\
e_{222} \\
e_{322} \\
e_{422} \\
e_{132} \\
e_{232} \\
e_{332} \\
e_{432} \\
e_{142} \\
e_{242} \\
e_{342} \\
e_{442} \\ 
\end{bmatrix}
\\
&=&
  \begin{array}{c}
\begin{matrix}
\beta_0 & \!\scriptsize{x11}\! & \scriptsize{x12} & \scriptsize{x13} & \!\scriptsize{x21}\! & \scriptsize{x22} & \scriptsize{x23} 
\end{matrix}\\
\left[\begin{array}{c|cc}
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 8\hphantom{0}  & 8^2 & 8^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 8\hphantom{0}  & 8^2 & 8^3 \\ 
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 2\hphantom{0}  & 2^2 & 2^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 2\hphantom{0}  & 2^2 & 2^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 4\hphantom{0}  & 4^2 & 4^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 4\hphantom{0}  & 4^2 & 4^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 6\hphantom{0}  & 6^2 & 6^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 6\hphantom{0}  & 6^2 & 6^3 \\
      1\hphantom{0} & 0\hphantom{0}  & 0^2  & 0^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 5\hphantom{0}  & 5^2  & 5^3  & 8\hphantom{0}  & 8^2 & 8^3  \\
      1\hphantom{0} & 10\hphantom{0} & 10^2 & 10^3 & 8\hphantom{0}  & 8^2 & 8^3 \\
      1\hphantom{0} & 15\hphantom{0} & 15^2 & 15^3 & 8\hphantom{0}  & 8^2 & 8^3 \\   
      \end{array}\right]
\end{array} 
\cdot
\begin{bmatrix}
1 \\
3  \\
0  \\
0  \\
2  \\
0   \\
0
\end{bmatrix}
+
  \begin{bmatrix}
e_{111} \\
e_{211} \\
e_{311} \\
e_{411} \\
e_{121} \\
e_{221} \\
e_{321} \\
e_{421} \\
e_{131} \\
e_{231} \\
e_{331} \\
e_{431} \\
e_{141} \\
e_{241} \\
e_{341} \\
e_{441} \\
e_{112} \\
e_{212} \\
e_{312} \\
e_{412} \\
e_{122} \\
e_{222} \\
e_{322} \\
e_{422} \\
e_{132} \\
e_{232} \\
e_{332} \\
e_{432} \\
e_{142} \\
e_{242} \\
e_{342} \\
e_{442} \\  
\end{bmatrix} 
\\
&=&
  \begin{bmatrix}
 5 \\
20 \\ 
35 \\ 
50  \\
 9  \\
24  \\
39  \\
54  \\
13  \\
28  \\
43  \\
58  \\
17  \\
32  \\
47  \\
62  \\
 5  \\
20  \\
35  \\
50  \\
 9  \\
24  \\
39  \\
54  \\
13  \\
28  \\
43  \\
58  \\
17  \\
32  \\
47  \\
62  
\end{bmatrix}
+
  \begin{bmatrix}
e_{111} \\
e_{211} \\
e_{311} \\
e_{411} \\
e_{121} \\
e_{221} \\
e_{321} \\
e_{421} \\
e_{131} \\
e_{231} \\
e_{331} \\
e_{431} \\
e_{141} \\
e_{241} \\
e_{341} \\
e_{441} \\
e_{112} \\
e_{212} \\
e_{312} \\
e_{412} \\
e_{122} \\
e_{222} \\
e_{322} \\
e_{422} \\
e_{132} \\
e_{232} \\
e_{332} \\
e_{432} \\
e_{142} \\
e_{242} \\
e_{342} \\
e_{442}    
\end{bmatrix} 
\end{eqnarray}
$$ 

With the `gexp` package:

```{r}
level2 <- seq(2, 8, 2)
cont_crd2 <- matrix(c(level2, level2^2, level2^3),
                    ncol = 3)

crd_hb <- gexp(mu = 1,
               r = 2,
               err = matrix(rep(0, 32),
                            nrow = 32),
               fe = list(f1 = c(3, 0, 0),
                         f2 = c(2, 0, 0)),
               fl = list(N = level,
                         P = level2),
               inte = c(6, 0, 0, 0, 0, 0, 0, 0, 0),
               contrasts = list(N = cont_crd,
                                P = cont_crd2),
               type = 'FE')

summary(crd_hb)
```

The plot below shows the simulated data.

```{r echo = FALSE, results = 'hide'}
z <- matrix(crd_hb$dfm$Y1,
            ncol = 4,
            byrow = FALSE)[seq(1, 8, 2), ]
persp(level,
      level2,
      z,
      theta = 30,
      xlab = 'N',
      ylab = 'P',
      col = 'red')
``` 

## Quadratic response surface

Let's omit matrix algebra as it is analogous to the one presented above. Considering now $\beta_2 = 2$ and $\beta_5 = 3$, we have the following result using the `gexp` package.
 
```{r}
crd_hb2 <- gexp(mu = 1,
                r = 2,
                err = matrix(rep(0, 32),
                             nrow = 32),
                fe = list(f1 = c(3, 2, 0),
                          f2 = c(2, 3, 0)),
                fl = list(N = level,
                          P = level2),
                inte = c(2, 1, 0, 1, 8, 0, 0, 0, 0),
                contrasts = list(N = cont_crd,
                                 P = cont_crd2),
                type = 'FE')

summary(crd_hb2)
```

The plot below shows the simulated data.

```{r echo = FALSE, results = 'hide'}
z1 <- matrix(crd_hb2$dfm$Y1,
            ncol = 4,
            byrow = FALSE)[seq(1, 8, 2), ]
persp(level,
      level2,
      z1,
      theta = 30,
      xlab = 'N',
      ylab = 'P',
      col = 'red')
``` 

# Experimental planning
## Non-interactive layouts
### CRD

Suppose we want to plan a CRD with three treatments (default) and five replicates (default). When only the layout matters, treatment effects can be arbitrary.

With the `gexp` package:
 
```{r}
crd_design <- gexp()
plot(crd_design)
```

### RCBD without repetition

Consider an RCBD in which each block is a replicate, with three treatments (default) and three blocks (default).

With the `gexp` package:

```{r}
rcbd_designs <- gexp(r = 1,
                     design = 'RCBD')
```

Plot:

```{r fig.width = 7, fig.height = 7}
plot(rcbd_designs)
```

### RCBD with repetition

Using the same structure as above, we now use two replicates per block.

With the `gexp` package:

```{r}
rcbd_designc <- gexp(r = 2,
                     design = 'RCBD')
```

Plot:

```{r fig.width = 7, fig.height = 7}
plot(rcbd_designc)
```

### LSD

Consider an LSD with three treatments (default), three rows, and three columns.

```{r}
lsd_design <- gexp(r = 1,
                   design = 'LSD')
```

Plot:

```{r fig.width = 7, fig.height = 7}
plot(lsd_design)
```
 
### Factorial - CRD

Suppose we are interested in planning a 2x3 factorial (default) design with three replicates. 

```{r}
fat_design <- gexp(r = 3,
                   type = 'FE')
```

The layout plot:

```{r fig.width = 7, fig.height = 7}
plot(fat_design)
```

### Split-plot - CRD

To plan a CRD split-plot with two whole-plot treatments (default), three subplot treatments (default), and three replicates:
 
```{r}
split_design <- gexp(r = 3,
                     type = 'SPE')
summary(split_design)
```

The layout plot:

```{r fig.width = 7, fig.height = 7}
plot(split_design)
```

## Interactive
### CRD

If you have a sketch of the experimental area (PNG, JPG, or JPEG), `gexp` can overlay randomized treatment labels on the image. For example, with six tanks and two treatments in three replicates each: 

<video width="320" height="240" controls>
  <source src="../inst/video/crd.mp4" type="video/mp4">
</video>
 
```{r eval = FALSE}
crd_plot <- gexp(r = 3,
                 fe = list(f1 = rep(1, 2)))
plot(crd_plot,
     dynamic = TRUE)
```

# References