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 and Schaalje (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 et al. 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:

Sources of Variation	Symbol
Treatments	\(\alpha,\tau\)
Block	\(\beta\)
Interaction	\(\gamma\)
Row block	\(\eta\)
Column block	\(\delta\)

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:

library(gexp)
crd <- gexp(mu = 15,
            err = matrix(rep(0, 6),
                         nrow = 6),
            r = 3,
            fe = list(alpha = c(1, -2)))
summary(crd)
#> Database 
#>   r  X1 Y1
#> 1 1 x11 16
#> 2 2 x11 16
#> 3 3 x11 16
#> 4 1 x12 13
#> 5 2 x12 13
#> 6 3 x12 13

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:

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)
#> Database 
#>   r  X1 Y1 Y2
#> 1 1 x11 16  8
#> 2 2 x11 16  8
#> 3 3 x11 16  8
#> 4 1 x12 13  9
#> 5 2 x12 13  9
#> 6 3 x12 13  9

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:

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)
#> Database 
#>   r Block  X1 Y1
#> 1 1     1 x11 18
#> 2 1     2 x11 20
#> 3 1     3 x11 22
#> 4 1     1 x12 15
#> 5 1     2 x12 17
#> 6 1     3 x12 19

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:

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)
#> Database 
#>   Row Column  X1 Y1
#> 1   1      1 x11 24
#> 2   2      1 x12 22
#> 3   3      1 x13 28
#> 4   1      2 x12 22
#> 5   2      2 x13 28
#> 6   3      2 x11 27
#> 7   1      3 x13 28
#> 8   2      3 x11 27
#> 9   3      3 x12 25

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:

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)
#> Database 
#>    r  X1  X2 Y1
#> 1  1 x11 x21 20
#> 2  2 x11 x21 20
#> 3  1 x12 x21 15
#> 4  2 x12 x21 15
#> 5  1 x11 x22 16
#> 6  2 x11 x22 16
#> 7  1 x12 x22  7
#> 8  2 x12 x22  7
#> 9  1 x11 x23 18
#> 10 2 x11 x23 18
#> 11 1 x12 x23 15
#> 12 2 x12 x23 15

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:

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)
#> Database 
#>    r  X1  X2 Y1
#> 1  1 x11 x21 20
#> 2  2 x11 x21 20
#> 3  3 x11 x21 20
#> 4  1 x12 x21 15
#> 5  2 x12 x21 15
#> 6  3 x12 x21 15
#> 7  1 x11 x22 16
#> 8  2 x11 x22 16
#> 9  3 x11 x22 16
#> 10 1 x12 x22  7
#> 11 2 x12 x22  7
#> 12 3 x12 x22  7

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:

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)
#> Database 
#>    r dose Y1
#> 1  1    0  2
#> 2  2    0  2
#> 3  3    0  2
#> 4  1    5 17
#> 5  2    5 17
#> 6  3    5 17
#> 7  1   10 32
#> 8  2   10 32
#> 9  3   10 32
#> 10 1   15 47
#> 11 2   15 47
#> 12 3   15 47

The plot below shows the simulated data.

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:

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)
#> Database 
#>    r dose  Y1
#> 1  1    0   2
#> 2  2    0   2
#> 3  3    0   2
#> 4  1    5 117
#> 5  2    5 117
#> 6  3    5 117
#> 7  1   10 432
#> 8  2   10 432
#> 9  3   10 432
#> 10 1   15 947
#> 11 2   15 947
#> 12 3   15 947

The plot below shows the simulated data.

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:

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)
#> Database 
#>    r  N P  Y1
#> 1  1  0 2   5
#> 2  2  0 2   5
#> 3  1  5 2  80
#> 4  2  5 2  80
#> 5  1 10 2 155
#> 6  2 10 2 155
#> 7  1 15 2 230
#> 8  2 15 2 230
#> 9  1  0 4   9
#> 10 2  0 4   9
#> 11 1  5 4 144
#> 12 2  5 4 144
#> 13 1 10 4 279
#> 14 2 10 4 279
#> 15 1 15 4 414
#> 16 2 15 4 414
#> 17 1  0 6  13
#> 18 2  0 6  13
#> 19 1  5 6 208
#> 20 2  5 6 208
#> 21 1 10 6 403
#> 22 2 10 6 403
#> 23 1 15 6 598
#> 24 2 15 6 598
#> 25 1  0 8  17
#> 26 2  0 8  17
#> 27 1  5 8 272
#> 28 2  5 8 272
#> 29 1 10 8 527
#> 30 2 10 8 527
#> 31 1 15 8 782
#> 32 2 15 8 782

The plot below shows the simulated data.

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.

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)
#> Database 
#>    r  N P     Y1
#> 1  1  0 2     17
#> 2  2  0 2     17
#> 3  1  5 2    972
#> 4  2  5 2    972
#> 5  1 10 2   3727
#> 6  2 10 2   3727
#> 7  1 15 2   8282
#> 8  2 15 2   8282
#> 9  1  0 4     57
#> 10 2  0 4     57
#> 11 1  5 4   3542
#> 12 2  5 4   3542
#> 13 1 10 4  13727
#> 14 2 10 4  13727
#> 15 1 15 4  30612
#> 16 2 15 4  30612
#> 17 1  0 6    121
#> 18 2  0 6    121
#> 19 1  5 6   7776
#> 20 2  5 6   7776
#> 21 1 10 6  30231
#> 22 2 10 6  30231
#> 23 1 15 6  67486
#> 24 2 15 6  67486
#> 25 1  0 8    209
#> 26 2  0 8    209
#> 27 1  5 8  13674
#> 28 2  5 8  13674
#> 29 1 10 8  53239
#> 30 2 10 8  53239
#> 31 1 15 8 118904
#> 32 2 15 8 118904

The plot below shows the simulated data.

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:

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:

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

Plot:

plot(rcbd_designs)

RCBD with repetition

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

With the gexp package:

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

Plot:

plot(rcbd_designc)

LSD

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

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

Plot:

plot(lsd_design)

Factorial - CRD

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

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

The layout plot:

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:

split_design <- gexp(r = 3,
                     type = 'SPE')
summary(split_design)
#> Database 
#>    r  X1  X2   Y1
#> 1  1 x11 x21 30.6
#> 2  2 x11 x21 28.1
#> 3  3 x11 x21 28.6
#> 4  1 x12 x21 27.0
#> 5  2 x12 x21 31.4
#> 6  3 x12 x21 25.5
#> 7  1 x13 x21 28.3
#> 8  2 x13 x21 30.8
#> 9  3 x13 x21 28.1
#> 10 1 x11 x22 31.5
#> 11 2 x11 x22 26.6
#> 12 3 x11 x22 29.1
#> 13 1 x12 x22 29.0
#> 14 2 x12 x22 32.5
#> 15 3 x12 x22 27.5
#> 16 1 x13 x22 29.6
#> 17 2 x13 x22 31.0
#> 18 3 x13 x22 27.6

The layout plot:

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:

crd_plot <- gexp(r = 3,
                 fe = list(f1 = rep(1, 2)))
plot(crd_plot,
     dynamic = TRUE)

References

Aquino, Luiz Henrique de. 1992. “Técnica Experimental Com Animais i.” Unpublished manuscript.

Ferreira, Daniel Furtado. 2008. Estatística Multivariada. 1st ed. Edited by UFLA.

Naes, T., A. H. Aastveit, and N. S. Sahni. 2007. “Analysis of Split-Plot Designs: An Overview and Comparison of Methods.” Qual. Reliab. Engng. Int. 23: 801–20.

Rencher, Alvin C., and Bruce G. Schaalje. 2007. Linear Models in Statistics. 2nd ed. Edited by John Wiley & Sons.

gexp: package overview

Simulating a Completely Randomized Design (CRD)

Univariate case

Multivariate Case

Simulating a Randomized Completely Block Design (RCBD)

Simulating a Latin Square Design (LSD)

Factorial Experiments

CRD

Split-plot experiments

CRD

Simulating a CRD with a quantitative factor

Non-Orthogonal Contrast

Linear Effect

Quadratic effect

Factorial CRD with two quantitative factors

Linear response surface

Quadratic response surface

Experimental planning

Non-interactive layouts

CRD

RCBD without repetition

RCBD with repetition

LSD

Factorial - CRD

Split-plot - CRD

Interactive

CRD

References