Inference on Treatment Effects

Inference on treatment effects often involves testing omnibus hypotheses and contrasts. These can be formulated using the general linear hypothesis:

\[ H_0: K'\beta = 0 \]

where \(K'\) is a contrast matrix. When the variance-covariance parameters in \(G\) and \(R\) are known, the estimate of \(\beta\) is:

\[ \hat{\beta} = (X^TV^{-1}X)^{-1}X^TV^{-1}y \]

And its variance is:

\[ C = (X^TV^{-1}X)^{-1} \]

The sampling distribution of \(K'\hat{\beta}\) is:

\[ K'\hat{\beta} \sim N(0, K'CK) \]

However, in practical situations, the matrices \(G\) and \(R\) are unknown and must be estimated using methods like Maximum Likelihood (ML) or Restricted ML (REML). The estimate of \(\beta\) is obtained by plugging in the estimated covariance matrices \(\hat{V}\), where:

\[ \hat{V} = Z\hat{G}Z^T + \hat{R} \]

The resulting estimate of \(\beta\) is:

\[ \hat{\beta} = (X^T\hat{V}^{-1}X)^{-1}X^T\hat{V}^{-1}y \]

And its estimated variance is:

\[ \hat{C} = (X^T\hat{V}^{-1}X)^{-1} \]

When testing the null hypothesis \(H_0: K'\beta = 0\), an approximate F-statistic is used. The F-statistic is given by:

\[ F = \frac{(K'\hat{\beta})' [K'\hat{C}K']^{-1} (K'\hat{\beta})}{\text{rank}(K)} \]

\(F\) follows an approximate F-distribution \(F(v_1, v_2)\) under \(H_0\), where \(v_1 = \text{rank}(K) \geq 1\) represents the numerator degrees of freedom (df), \(v_2\) is the denominator df.

When \(\text{rank}(K) = 1\), the F-statistic simplifies to the square of the t-statistic:

\[ F = t^2 \]where \(t = \frac{K'\hat{\beta}}{\sqrt{K'\hat{C}K}}\) with \(v_2\) df.

In balanced designs, where data is analyzed using a variance components model—commonly applied in experimental animal research—\(v_2\)​ can be precisely determined through degrees of freedom decomposition, as applied in analysis of variance (ANOVA):

\[ v_2 = n - 1 - [\text{rank}(X) - 1] - [\text{rank}(Z) - 1] \]

However, for general cases, such as unbalanced designs or models with correlated random intercept and slope effects, \(v_2\) must be approximated using methods like Satterthwaite’s approximation (Fai and Cornelius, 1996) or Kenward-Roger’s method (Kenward and Roger, 1997), as implemented in the lmerTest package.

Power Calculation Under the Alternative Hypothesis

Under the alternative hypothesis \(H_A: K'\beta \neq 0\), the F-statistic follows a non-central distribution \(F(v_1, v_2, \phi)\), where \(\phi\) is the non-centrality parameter that measures the departure from the null hypothesis \(H_0\). The non-centrality parameter \(\phi\) is given by:

\[ \phi = (K'\hat{\beta})' [K'\hat{C}K']^{-1} (K'\hat{\beta}) \]

Once the distribution of the F-statistic under \(H_A\) is known, the power of the test can be calculated as the conditional probability of rejecting \(H_0\) when \(H_A\) is true:

\[ \text{Power} = P(\text{reject } H_0: F > F_{\text{crit}} \mid H_A) \]

Where: \(F_{\text{crit}}\) is the critical value of the F-statistic used to reject \(H_0\), determined such that \(P(F > F_{\text{crit}} \mid H_0) = \alpha\), where \(\alpha\) is the type I error rate.

The determination of the degrees of freedom \(v_1\) and \(v_2\), as well as the non-centrality parameter \(\phi\), are critical steps for power calculation. Notably, \(\phi\) resembles the numerator of the F-statistic but with population parameters for \(\beta\), \(G\), and \(R\) replacing their sample estimates \(\hat{\beta}\), \(\hat{G}\), and \(\hat{R}\).

Generally, power analysis requires specifying the following components:

• The Design to be evaluated, which determines the matrices \(X\) (fixed effects) and \(Z\) (random effects).
• The Treatment Effects, which determine \(\beta\) (the fixed effect coefficients).
• The Variance Components, which determine \(G\) (the covariance matrix of the random effects) and \(R\) (the covariance matrix of the residuals).

A key aspect of conducting a valid power analysis is obtaining reasonable estimates for the magnitude of the parameters that will be used in the model. This includes:

• Treatment Effects (\(\beta\)): The size of the treatment effect(s) you expect to detect. This can be obtained from previous studies, pilot experiments, or subject-matter expertise.
• Variance Components (\(G\) and \(R\)): These represent the variability in the random effects and residual errors. Variance estimates are often obtained from:
  • Pilot studies or preliminary data, which can provide initial estimates of variability in random effects (e.g., subject-to-subject variability or group-level variability).
  • Literature on similar experiments, where variance components have been reported.
  • Subject-matter expertise, where researchers provide estimates based on their knowledge of the system being studied.

Performing a power analysis with unrealistic parameter magnitudes can lead to incorrect conclusions, either overestimating the likelihood of detecting a treatment effect or requiring an unnecessarily large sample size.

Step1: Creating a Design Object

Design objects in pwr4exp can be created using functions that generate several standard experimental designs available in the package. These functions include:

• designCRD(treatments, label, replicates, formula, beta, sigma2): For completely randomized design (CRD).
• designRCBD(treatments, label, blocks, formula, beta, VarCov, sigma2, ...): For randomized complete block design (RCBD).
• designLSD(treatments, label, squares, reuse, formula, beta, VarCov, sigma2, ...): For Latin square design (LSD).
• designCOD(treatments, label, squares, formula, beta, VarCov, sigma2, ...): For crossover design (COD), which is a special case of LSD with time periods and individuals as blocks. Period blocks are reused when replicating squares.
• designSPD(trt.main, trt.sub, label, replicates, formula, beta, VarCov, sigma2, ...): For split-plot design (SPD).

Step2: Calculating Power or Sample Size

Once the design object is created, calculating power or sample size is straightforward. Power for omnibus tests, including main effects and their interactions (if specified in the model during the creation of the design object), can be calculated using:

• pwr.anova(design, alpha = 0.05, ...)

The required inputs include:

1. design, the object created in Step 1;
2. alpha, indicating the Type I error rate, with a default value of 0.05.

For specific contrasts, power can be calculated using:

• pwr.contrast(design, alpha = 0.05, spec, method, ...)

The syntax of the emmeans package is inherited to specify contrasts of interest. The required inputs are:

1. design;
2. alpha;
3. spec, an argument inherited from emmeans, which specifies the names of the factors over which the contrasts are performed.
4. method, another argument inherited from emmeans, which specifies the method of contrasts (e.g., “pairwise”, “trt.vs.ctrl”, “poly”).

The minimal sample size needed to achieve a target power can be determined using:

• find_sample_size(design.quote, alpha = 0.05, target.power = 0.8, n_init = 2, n_max = 99)

This function calculates the minimum sample size necessary by incrementally checking integers from n_init to n_max. The required inputs are:

1. design.quote: A quoted design object with an unknown and unevaluated replication argument, to be evaluated with varying values;
2. alpha;
3. target.power: A single value specifying the target power for all effects, or a vector specifying individual target power levels for each effect, with a default value of 0.05;
4. n_init: The initial number of replications for the iterative process, with a default value of 2;
5. n_max: The maximum number of replications for the iterative process, with a default value of 99.

Currently, sample size determination is available only for omnibus tests and not for specific contrasts in pwr4exp.

Example 1. Completely Randomized Design

In this example, we will create a CRD with one treatment factor. The design parameters are as follows:

1. Treatments: 1 treatment factor with 4 levels;
2. Replicates: 8 experimental units per treatment.
3. Mean and effect size: The mean of control is 35, and the effects of other three treatments are -5, +2, and +3.
4. Error variance: The variance of response variable is 15.

Create the CRD

crd <- designCRD(
  treatments = 4,
  replicates = 8,
  beta = c(35, -5, 2, 3),
  sigma2 = 15
)

Power for the omnibus test

The power of the omnibus test (i.e., F-test) can be calculated using the pwr.anova function. Under the type I error rate of 0.05, the power for testing an overall difference among treatments is 0.95467.

pwr.anova(design = crd)
#> Power analysis of Completely Randomized Design
#>     NumDF DenDF non-centrality alpha   power
#> trt     3    28         20.267  0.05 0.95467

Power for specific contrasts

To assess the power for specific contrasts, use the pwr.contrast function. For example, to calculate the power for detecting differences between treatments and the control:

pwr.contrast(design = crd, specs =  ~ trt, method = "trt.vs.ctrl")
#>      contrast estimate df non-centrality alpha     power
#> 1 trt2 - trt1       -5 28       6.666667  0.05 0.7028739
#> 2 trt3 - trt1        2 28       1.066667  0.05 0.1694975
#> 3 trt4 - trt1        3 28       2.400000  0.05 0.3216803

To calculate the power for detecting linear or polynomial trends across the treatment levels:

pwr.contrast(design = crd, specs =  ~ trt, method = "poly")
#>    contrast estimate df non-centrality alpha     power
#> 1    linear       16 28       6.826667  0.05 0.7130735
#> 2 quadratic        6 28       4.800000  0.05 0.5617849
#> 3     cubic      -18 28       8.640000  0.05 0.8098383

Example 2. Randomized Complete Block Design

In this example, we will create an RCBD with two treatment factors. The design parameters are as follows:

1. Treatments: A 2x2 factorial design.
2. Replicates: 8 blocks.
3. Mean and effect size: The mean of the control (A1B1) is 35. The effect of A2 alone is an increase of 5 units, and the effect of B2 alone is an increase of 3 units. The interaction between A2 and B2 introduces an additional effect of -2 units, meaning the combined effect of A2 and B2 is 2 units below the sum of their individual effects. The corresponding cell means are:\[ \begin{array}{c|c|c} & B1 & B2 \\ \hline A1 & 35 & 38 \\ A2 & 40 & 41 \\ \end{array} \]
4. Variance among blocks: 11.
5. Error variance: 4. The total variance of the response variable (15) is decomposed into variance between blocks (11) and variance within blocks (4).

Create the RCBD

rcbd <- designRCBD(
  treatments = c(2, 2),
  blocks = 8,
  beta = c(35, 5, 3, -2),
  VarCov = 11,
  sigma2 = 4
)

Power for the omnibus test

pwr.anova(design = rcbd)
#> Power Analysis of Randomized Complete Block Design
#>           NumDF DenDF non-centrality alpha   power
#> facA          1    21             32  0.05 0.99969
#> facB          1    21              8  0.05 0.76950
#> facA:facB     1    21              2  0.05 0.27138

Power for specific contrasts

The power for detecting the effect of facA, both overall and within each level of facB, can be assessed as follows:

# across all levels of facB
pwr.contrast(design = rcbd, specs = ~ "facA", method = "pairwise")
#> NOTE: Results may be misleading due to involvement in interactions
#>        contrast estimate df non-centrality alpha    power
#> 1 facA1 - facA2       -4 21             32  0.05 0.999691

# at each level of facB
pwr.contrast(design = rcbd, specs = ~ facA|facB, method = "pairwise")
#>        contrast facB estimate df non-centrality alpha     power
#> 1 facA1 - facA2    1       -5 21             25  0.05 0.9974502
#> 2 facA1 - facA2    2       -3 21              9  0.05 0.8160596

Sample Size Determination

To determine the number of blocks required to achieve the target power, the find_sample_size function can be used. First, we create a quoted design object where blocks = n remains unevaluated:

rcbd_quote <- quote(
  designRCBD(
  treatments = c(2, 2),
  blocks = n,
  beta = c(35, 5, 3, -2),
  VarCov = 11,
  sigma2 = 4
  )
)

The optimal sample size for the target power within the range of n_init and n_max can be determined as follows:

find_sample_size(design.quote = rcbd_quote, n_init = 2, n_max = 99)
#>           alpha     power best_n
#> facA       0.05 0.8212779      3
#> facB       0.05 0.8207219      9
#> facA:facB  0.05 0.8115100     33

Example 3. Latin Square Design

In this example, we extend the design from Example 2 by introducing another blocking factor, thus creating an LSD. The treatment structure and effect sizes remain the same as in Example 2. The design controls for two sources of variability (row and column blocks) while evaluating the treatment effects. In the LSD, the total variance (15) is further decomposed into three components:

• Variance between row blocks (11),
• Variance between column blocks (2), and
• Residual error variance (2).

In this example, we will customize the labels for the two factors as follows: “temp” with levels “T1” and “T2”, and “dosage” with levels “D1” and “D2”.

Create the LSD

lsd <- designLSD(
  treatments = c(2, 2),
  label = list(temp = c("T1", "T2"), dosage = c("D1", "D2")),
  squares = 4,
  reuse = "both",
  beta = c(35, 5, 3, -2),
  VarCov = list(11, 2),
  sigma2 = 2
)

Power for the omnibus test

pwr.anova(design = lsd)
#> Power Analysis of Latin Square Design
#>             NumDF DenDF non-centrality alpha   power
#> temp            1    33            128  0.05 1.00000
#> dosage          1    33             32  0.05 0.99979
#> temp:dosage     1    33              8  0.05 0.78387

Power for specific contrasts

The power for detecting the effect of dosage, both overall and within each level of temp, can be assessed as follows:

# the effect of dosage across all levels of temp
pwr.contrast(design = lsd, specs = ~ "dosage", method = "pairwise")
#> NOTE: Results may be misleading due to involvement in interactions
#>   contrast estimate df non-centrality alpha     power
#> 1  D1 - D2       -2 33             32  0.05 0.9997892

# the effect of dosage at each level of temp
pwr.contrast(design = lsd, specs = ~ dosage|temp, method = "pairwise")
#>   contrast temp estimate df non-centrality alpha     power
#> 1  D1 - D2   T1       -3 33             36  0.05 0.9999429
#> 2  D1 - D2   T2       -1 33              4  0.05 0.4927485

Sample Size Determination

To determine the number of squares required to achieve the target power, the find_sample_size function can be used. First, we create a quoted design object where squares = n remains unevaluated:

lsd_quote <- quote(
  designLSD(
  treatments = c(2, 2),
  squares = n,
  reuse = "both",
  beta = c(35, 5, 3, -2),
  VarCov = list(11, 2),
  sigma2 = 2
  )
)

The optimal sample size for the target power within the range of n_init and n_max can be determined as follows:

find_sample_size(design.quote = lsd_quote, n_init = 2, n_max = 99)
#>           alpha     power best_n
#> facA       0.05 1.0000000      2
#> facB       0.05 0.9618851      2
#> facA:facB  0.05 0.8705999      5

Example 4: Split-plot Design

In this example, we will create a SPD with two treatment factors, one at each level. The design parameters are as follows:

1. Treatments: One main plot factor having 2 levels, and another factor with 3 levels at the sub-plot level.

2. Label: The two treatment factors are labeled as “Main” with levels “Main1” and “Main2,” and “Sub” with levels “Sub1” and “Sub2,” corresponding to the main plot and subplot, respectively.

3. Replicates: There are 10 main plots, with 5 main plots per “Main” treatment level, resulting in 10 experimental units (blocks) for each “Sub” treatment level.

4. The mean of the control (Main1 with Sub1) is 20. The effect of Main2 alone is an increase of 2 units. The effects of Sub2 and Sub3 alone are increases of 2 and 4 units, respectively. The interaction between Main2 and Sub2 is zero. The interaction between Main2 and Sub3 introduces an additional effect of 2 units, meaning the combined effect of Main2 and Sub3 is 2 units above the sum of their individual effects. The corresponding cell means are:\[ \begin{array}{c|c|c} & Sub1 & Sub2 & Sub3 \\ \hline Main1 & 20 & 22 & 24\\ Main2 & 22 & 24 & 28 \\ \end{array} \]

5. The total variance (15) is assumed to decompose into 4 for whole-plot error and 11 for subplot error.

Create the SPD

spd <- designSPD(
  trt.main = 2,
  trt.sub = 3, 
  replicates = 10, 
  label = list(Main = c("Main1", "Main2"), Sub = c("Sub1", "Sub2", "Sub3")),
  beta = c(20, 2, 2, 4, 0, 2),
  VarCov = list(4),
  sigma2 = 11
)

Power for the omnibus test

pwr.anova(spd)
#> Power Analysis of Split Plot Design
#>          NumDF DenDF non-centrality alpha   power
#> Main         1    18         4.6377  0.05 0.53114
#> Sub          2    36        23.0303  0.05 0.98924
#> Main:Sub     2    36         1.2121  0.05 0.14311

Power for specific contrasts

The power for detecting the effect of subplot treatment under each level of main plot treatment, can be assessed as follows:

pwr.contrast(design = spd, specs = ~ Sub|Main, method = "trt.vs.ctrl")
#>      contrast  Main estimate df non-centrality alpha     power
#> 1 Sub2 - Sub1 Main1        2 36       1.818182  0.05 0.2592167
#> 2 Sub3 - Sub1 Main1        4 36       7.272727  0.05 0.7467531
#> 3 Sub2 - Sub1 Main2        2 36       1.818182  0.05 0.2592167
#> 4 Sub3 - Sub1 Main2        6 36      16.363636  0.05 0.9758744

Example 5: A customized design

In this example, we will create a COD arranged at the “subplot” level of a split-plot layout.

Sixteen subjects, split evenly between two breeds (8 subjects per breed), are enrolled in a 4x4 crossover design. Within each breed, the subjects are further divided into two squares (for a total of four squares), with each square consisting of four subjects. Each subject follows one of four treatment sequences and receives treatments across four periods. This hybrid SPD-COD combines elements of a SPD and a COD to evaluate treatment effects across different breeds. The Latin Square structure ensures that within each square, treatments are balanced across periods, while the split-plot element accounts for breed as a whole-plot factor. The overall structure is as follows:

• Treatment structure

  • Whole-plot level: Breed (2 breeds, with 8 subjects in each breed)
  • Sub-plot level: A 2x2 factorial design with two treatment factors, facA and facB.

• Experimental units

  • Whole-plot factor: Subjects (randomly selected from two breeds)
  • Sub-plot factor: Measurements taken from subjects at each period.

• Variance decomposition The overall variance (15) is decomposed into three main components: subject variance (7), period variance (4), and random error (4).

Construct the design data frame

To create the design object for this hybrid SPD-COD, we first need to construct a data frame containing all the independent variables, structured to resemble the actual experimental data. This data frame can be generated using the internal function df.cod, which builds a data frame for a COD with the specified treatment structure and number of replications. Once the base data frame is generated, a column for the “Breed” variable is manually added. It is important to note that no randomization occurs at this stage—the data frame simply represents the data structure and is not the actual randomized experimental data. Its purpose is to create the design layout.

df_spd_cod <- pwr4exp:::df.cod(
  treatments = c(2, 2),
  squares = 4
)
## Create main plot factor, i.e., breed
df_spd_cod$Breed <- rep(c("1", "2"), each = 32)

## Check data structure
head(df_spd_cod, n = 4); tail(df_spd_cod, n = 4)
#>   subject period trt facA facB square Breed
#> 1       1      1   2    2    1      1     1
#> 2       2      1   3    1    2      1     1
#> 3       3      1   4    2    2      1     1
#> 4       4      1   1    1    1      1     1
#>    subject period trt facA facB square Breed
#> 61      13      4   1    1    1      4     2
#> 62      14      4   2    2    1      4     2
#> 63      15      4   3    1    2      4     2
#> 64      16      4   4    2    2      4     2

Specify the model formula

Next, an appropriate model formula must be specified to fit the experimental design. The formula captures the interactions between breed, the two treatment factors (facA and facB), and the random effects for subjects and periods.

formula <- y ~ Breed*facA*facB + (1|subject) + (1|period)

Fixed effects

We then specify the fixed effects for the model, where beta contains the baseline intercept and the effects of each factor and their interactions.

beta = c(
  `(Intercept)` = 35,        # Baseline (mean of Breed1_A1_B1)
  Breed2 = -5,           # Effect of the second breed alone
  facA2 = -5,               # Effect of A2 alone
  facB2 = 1,                # Effect of B2 alone
  `Breed2:facA2` = 1,         # Interaction between Breed2 and A2
  `Breed2:facB2` = 0,         # Interaction between Breed2 and B2
  `facA2:facB2` = 2,             # Interaction between A2 and B2
  `Breed2:facA2:facB2` = 1       # Three-way interaction between Breed2, A2, and B2
)

Create the design object

After constructing the data frame, model formula, and fixed effects, the design object can be created using the designCustom function. The variance-covariance structure and error variance are also provided to complete the design.

SPD_COD <- designCustom(
  design.df = df_spd_cod,
  formula = formula,
  beta = beta,
  VarCov = list(7, 4),
  sigma2 = 4,
  design.name = "hybrid SPD COD"
)

Power for the omnibus test

pwr.anova(SPD_COD)
#> Power Analysis of hybrid SPD COD
#>                 NumDF DenDF non-centrality alpha   power
#> Breed               1    14          9.031  0.05 0.79790
#> facA                1    39         42.250  0.05 0.99999
#> facB                1    39         20.250  0.05 0.99238
#> Breed:facA          1    39          2.250  0.05 0.30997
#> Breed:facB          1    39          0.250  0.05 0.07768
#> facA:facB           1    39          6.250  0.05 0.68372
#> Breed:facA:facB     1    39          0.250  0.05 0.07768

Power for specific contrasts

We can also calculate the power for specific contrasts. For example, to assess the power for detecting differences between facA1 and facA2 at each level of facB within each breed:

pwr.contrast(SPD_COD, ~facA|facB|Breed, "pairwise")
#>        contrast facB Breed estimate df non-centrality alpha     power
#> 1 facA1 - facA2    1     1        5 39             25  0.05 0.9982139
#> 2 facA1 - facA2    2     1        3 39              9  0.05 0.8328312
#> 3 facA1 - facA2    1     2        4 39             16  0.05 0.9737940
#> 4 facA1 - facA2    2     2        1 39              1  0.05 0.1641134