The mkdesign function

The function mkdesign is versatile for designs as long as the data frame is provided.

Usage

mkdesign(
  formula, data,
  beta = NULL, means = NULL,
  vcomp = NULL, sigma2 = NULL,
  correlation = NULL,
  template = FALSE,
  REML = TRUE
)

Arguments

• formula: A right-hand-side model formula that defines fixed and random effects. (Use lme4::lmer syntax for specifying random effects.)

• data: A data frame containing all independent variables required by the model, consistent with the design’s structure.

• means or beta: Expected fixed effects. Either a vector of model coefficients (beta) or expected means (means) must be provided. Templates for parameter ordering can be generated by setting template = TRUE.

• vcomp: Variance components for any random effects present in the model, provided as a numeric vector. A template for input ordering can be generated by setting template = TRUE.

• sigma2: The residual variance.

• correlation: (Optional) A correlation structure (nlme::corClasses) for repeated measures.

• template: When set to TRUE or when only the formula and data arguments are provided, templates for beta, means, and vcomp are generated to indicate the required input order.

• REML: (logical) Whether to use the REML or ML information matrix. The default is TRUE (REML).

Specific design functions

While mkdesign provides full flexibility, pwr4exp also offers specific functions for common experimental designs. These functions allow users to define a design based on treatment structure and replications without manually creating a data frame.

Completely randomized design (CRD)

designCRD(treatments, label, replicates, formula, beta, means, sigma2, template = FALSE)

Randomized complete block design (RCBD)

designRCBD(treatments, label, blocks, formula, beta, means, vcomp, sigma2, template = FALSE)

Latin square design (LSD)

designLSD(treatments, label, squares, reuse, formula, beta, means, vcomp, sigma2, template = FALSE)

Crossover design

This is a special case of LSD where time periods and individuals act as blocks. Period blocks are reused when replicating squares.

designCOD(treatments, label, squares, formula, beta, means, vcomp, sigma2, template = FALSE)

Split-plot design (SPD)

designSPD(trt.main, trt.sub, label, replicates, formula, beta, means, vcomp, sigma2, template = FALSE)

The inputs for these functions are similar to mkdesign, except that data is replaced by predefined design characteristics such as treatment levels and replications.

Key inputs

• treatments: A numeric vector specifying the number of levels for treatment factors. Multiple factors are arranged factorially. For example, treatments = c(3, 2) specifies two treatment factors one with 3 levels and another with 2 levels, forming a factorial treatment design.

• trt.main and trt.sub: Define main-plot and subplot treatments for SPD, following the same rule for treatment.

• label: (optional) A list assigning labels to factor levels. If not specified, default names are assigned. For one treatment factor, the default is list(trt = c("1", "2", ...)). For two factors, the default is list(facA = c("1", "2", ...), facB = c("1", "2", ...)), where “facA” and “facB” represent the two factors, and “1”, “2”, etc., represent the levels of each factor.

• replicates: Number of experimental units per treatment in CRD or number of main plots (i.e., the number of experimental units per main-plot treatment) in SPD.

• blocks: Number of blocks in RCBD.

• squares: Number of squares in LSD or COD.

• reuse: Specifies how squares are replicated in LSD. One of:

• "row": reuse row blocks
• "col": reuse column blocks
• "both": reuse both row and column blocks

• template: When set to TRUE, templates for beta, means, and vcomp are generated to indicate the required input order.

Each of these design-generating functions has a default formula based on the treatment structure (e.g., one factor or factorial factors). If formula is not specified, a default formula with main effects and all interactions (if applicable) is used internally. In RCBD, LSD, COD, and SPD designs, all block factors are fitted as random effects. The formula component in the output design object (a list) can be inspected.

Parameter templates

Templates help ensure the correct ordering of fixed effects (beta or means) and random effects (vcomp).

Below, an exemplary data frame is created to illustrate these templates. Note that this example does NOT represent a realistic design. The data frame includes:

• Four categorical variables (fA, fB, fC, fD),
• Two numerical variables (x, z).

 df1 <- expand.grid(
   fA = factor(1:2), # factor A with 2 levels
   fB = factor(1:2), # factor B with 2 levels
   fC = factor(1:3), # factor C with 3 levels
   fD = factor(1:3), # factor D with 3 levels
   subject = factor(1:10)  # 10 subjects
 )
 df1$x <- rnorm(nrow(df1))  # Numerical variable x
 df1$z <- rnorm(nrow(df1))  # Numerical variable z

mkdesign( ~ fA * fB + x, df1)$fixeff$beta
#> (Intercept)         fA2         fB2           x     fA2:fB2 
#>           1           2           3           4           5

mkdesign(~ fA + fB, df1, template = TRUE)$fixeff$means
#> fA1 fA2 fB1 fB2 
#>   1   2   3   4

mkdesign(~ fA * fB * fC, df1)$fixeff$means
#> fA1:fB1:fC1 fA2:fB1:fC1 fA1:fB2:fC1 fA2:fB2:fC1 fA1:fB1:fC2 fA2:fB1:fC2 
#>           1           2           3           4           5           6 
#> fA1:fB2:fC2 fA2:fB2:fC2 fA1:fB1:fC3 fA2:fB1:fC3 fA1:fB2:fC3 fA2:fB2:fC3 
#>           7           8           9          10          11          12

mkdesign(~ x + z, df1)$fixeff$means
#> (Intercept)           x           z 
#>           1           2           3

mkdesign(~ x * z, df1)$fixeff$means
#> (Intercept)           x           z         x:z 
#>           1           2           3           4

mkdesign(~ fA * x, df1)$fixeff$means
#>   fA1   fA2 fA1:x fA2:x 
#>     1     2     3     4

mkdesign(~ fA * fB * fC + fD * x + z, df1)$fixeff$means
#>         fD1         fD2         fD3           z       fD1:x       fD2:x 
#>           1           2           3           4           5           6 
#>       fD3:x fA1:fB1:fC1 fA2:fB1:fC1 fA1:fB2:fC1 fA2:fB2:fC1 fA1:fB1:fC2 
#>           7           8           9          10          11          12 
#> fA2:fB1:fC2 fA1:fB2:fC2 fA2:fB2:fC2 fA1:fB1:fC3 fA2:fB1:fC3 fA1:fB2:fC3 
#>          13          14          15          16          17          18 
#> fA2:fB2:fC3 
#>          19

mkdesign(~ fA * fB * fC * fD + (1 + fA | subject), df1)$varcov
#> $subject
#>             (Intercept) fA2
#> (Intercept)           1   2
#> fA2                   2   3

Correlation structures

Although specifying complex variance-covariance structures, such as unstructured nlme::corSymm, may seem impractical for power analysis, various correlation structures from the nlme package can be applies.

The available correlation structures, as outlined in nlme::corClasses, including

• corAR1 – First-order autoregressive

• corARMA – Autoregressive moving average

• corCAR1 – Continuous-time autoregressive

• corCompSymm – Compound symmetry

• corExp – Exponential spatial correlation

• corGaus – Gaussian spatial correlation

• corLin – Linear spatial correlation

• corSymm – Unstructured correlation

• corRatio – Ratio-based spatial correlation

• corSpher – Spherical spatial correlation

Note: In nlme::corAR1() and nlme::corARMA() when p=1 and q=0, the time variable must be an integer. However, in pwr4exp, this restriction has been released, factors are also supported.

Power Calculation for Omnibus F-Tests

The statistical power for F-tests (ANOVA-like tests) can be computed using the pwr.anova function.

Usage

pwr.anova(object, sig.level = 0.05, type = 3)

Arguments

• object: the design object created in the previous step.

• sig.level: significance level, default 0.05.

• type: the type of ANOVA table (default: 3).

Power Calculation for Specific Contrasts

To evaluate statistical power for t-tests on specific contrasts, use the pwr.contrast function.

Usage

pwr.contrast(object, which, by = NULL, 
  contrast = c("pairwise", "poly", "trt.vs.ctrl"),
  sig.level = 0.05, p.adj = FALSE, 
  alternative = c("two.sided", "one.sided"), strict = TRUE
)

Arguments

• object: The design object.

• which: The factor of interest. Multiple factors can be combined using : or * (e.g., "facA*facB" treats combined factor levels as a single factor).

• by: The variable to condition on.

• contrast: Contrast method, on of

  • "pairwise" - pairwise comparisons

  • "poly" - polynomial contrasts

  • "trt.vs.ctrl" - treatment vs. control comparison

  • Alternatively, a numeric vector for a single custom contrast, or a named list of contrast vectors. If a list is provided, each vector must match the number of factor levels in each by group. In multi-factor scenarios, factor levels are combined and treated as a single factor.

• sig.level: significance level (default: 0.05).

• p.adj: whether the sig.level should be adjusted using the Bonferroni method (default: FALSE).

• alternative: "two.sided" (default) or "one.sided".

• strict: Controls how power is calculated in two-sided tests (default: TRUE). If TRUE, includes the probability of rejection in the opposite direction of the true effect. If FALSE, power equals half the significance level when the true difference is zero.

Power Calculation for Model Coefficients

The pwr.summary function computes the statistical power for testing model coefficients (t-tests).

Usage

pwr.summary(object, sig.level = 0.05)

Arguments

• object: A design object
• sig.level: The significance level (default: 0.05)

Example 1. Completely Randomized Design

In this example, we create a CRD with one treatment factor.

Design parameters

1. Treatments: 1 treatment factor with 4 levels.
2. Replicates: 8 experimental units per treatment.
3. Mean and effect size: Expected means for four levels: 35, 30, 37, 38
4. Error variance: 15

Step 1: Creating the CRD

crd <- designCRD(
  treatments = 4,
  replicates = 8,
  means = c(35, 30, 37, 38),
  sigma2 = 15
)

Step 2: Power for omnibus test

To compute the power for the overall F-test (ANOVA):

pwr.anova(crd)
#> Power of type III F-test
#>     NumDF DenDF sig.level   power
#> trt     3    28      0.05 0.95467

Step 3: Power for specific contrasts

a) Pairwise Comparisons

To compute power for pairwise differences:

pwr.contrast(crd, which =  "trt", contrast = "pairwise")
#>             effect df sig.level      power alternative
#> trt1 - trt2      5 28      0.05 0.70287390   two.sided
#> trt1 - trt3     -2 28      0.05 0.16949749   two.sided
#> trt1 - trt4     -3 28      0.05 0.32168033   two.sided
#> trt2 - trt3     -7 28      0.05 0.93677955   two.sided
#> trt2 - trt4     -8 28      0.05 0.97860686   two.sided
#> trt3 - trt4     -1 28      0.05 0.07896844   two.sided

b) Polynomial Contrasts

To test for trends across treatment levels:

pwr.contrast(crd, which =  "trt", contrast = "poly")
#>           effect df sig.level     power alternative
#> linear        16 28      0.05 0.7130735   two.sided
#> quadratic      6 28      0.05 0.5617849   two.sided
#> cubic        -18 28      0.05 0.8098383   two.sided

c) Treatment vs. Control Comparison

The power for detecting differences between treatments and the control:

pwr.contrast(crd, which =  "trt", contrast = "trt.vs.ctrl")
#>             effect df sig.level     power alternative
#> trt2 - trt1     -5 28      0.05 0.7028739   two.sided
#> trt3 - trt1      2 28      0.05 0.1694975   two.sided
#> trt4 - trt1      3 28      0.05 0.3216803   two.sided

d) Custom Contrast: All Treatments vs. Control

In addition to pre-defined contrast method, custom contrast vectors can be specified. For example, to compare the average of all treatments against the control, use the contrast vector c(-1, 1/3, 1/3, 1/3). The power for this test can be computed as:

pwr.contrast(crd, which =  "trt", contrast = list(trts.vs.ctrl = c(-1, 1/3, 1/3, 1/3)))
#>                    effect df sig.level power alternative
#> trts.vs.ctrl 2.187139e-14 28      0.05  0.05   two.sided

Adjusting for Multiple Comparisons

In real-world analysis, P-values often need adjustment for multiple comparisons (MCP). While most post-hoc methods cannot be applied directly in power analysis, we can lower the significance level to approximate MCP adjustments.

Using a More Conservative Significance Level

pwr.contrast(crd, which =  "trt", contrast = "pairwise", sig.level = 0.01)
#>             effect df sig.level     power alternative
#> trt1 - trt2      5 28      0.01 0.4418907   two.sided
#> trt1 - trt3     -2 28      0.01 0.0546995   two.sided
#> trt1 - trt4     -3 28      0.01 0.1320866   two.sided
#> trt2 - trt3     -7 28      0.01 0.7946290   two.sided
#> trt2 - trt4     -8 28      0.01 0.9042775   two.sided
#> trt3 - trt4     -1 28      0.01 0.0194487   two.sided

Using the Bonferroni Correction

The pwr.contrast function includes an option for Bonferroni correction, though it may be overly conservative:

pwr.contrast(crd, which =  "trt", contrast = "pairwise", sig.level = 0.05, p.adj = TRUE)
#>             effect df   sig.level      power alternative
#> trt1 - trt2      5 28 0.008333333 0.41456682   two.sided
#> trt1 - trt3     -2 28 0.008333333 0.04782486   two.sided
#> trt1 - trt4     -3 28 0.008333333 0.11835238   two.sided
#> trt2 - trt3     -7 28 0.008333333 0.77333066   two.sided
#> trt2 - trt4     -8 28 0.008333333 0.89102508   two.sided
#> trt3 - trt4     -1 28 0.008333333 0.01655798   two.sided

Example 2. Randomized complete block design

In this example, we create a Randomized Complete Block Design (RCBD) with two treatment factors in a 2×2 factorial design.

Design Parameters

1. Treatments: A 2x2 factorial design (factors: A and B).
2. Replicates: 8 blocks.
3. Expected Means: \[ \begin{array}{c|c|c} & B1 & B2 \\ \hline A1 & 35 & 38 \\ A2 & 40 & 41 \\ \end{array} \] Corresponding beta values are as follows (Optional):

• Intercept (A1B1): 35
• Effect of A2 alone: 5 units
• Effect of B2 alone: 3 units
• Interaction (A2:B2): -2 units (i.e., the combined effect of A2 and B2 is 2 units below the additive effects).

4. Variance among blocks: 11.
5. Error variance: 4.
6. The total variance of the response variable (15) is decomposed into variance between blocks (11) and variance within blocks (4).

Step 1: Creating the RCBD

• Generate the templates to determine the correct parameter order. Provide beta or means, and vcomp according to the order below.

designRCBD(treatments = c(2, 2), blocks = 8, template = TRUE)
#> $fixeff
#> $fixeff$beta
#> (Intercept)       facA2       facB2 facA2:facB2 
#>           1           2           3           4 
#> 
#> $fixeff$means
#> facA1:facB1 facA2:facB1 facA1:facB2 facA2:facB2 
#>           1           2           3           4 
#> 
#> 
#> $varcov
#> $varcov$block
#>             (Intercept)
#> (Intercept)           1

• Create the design object

rcbd <- designRCBD(
  treatments = c(2, 2),
  blocks = 8,
  # beta = c(35, 5, 3, -2), # identical to means
  means = c(35, 40, 38, 41),
  vcomp = 11,
  sigma2 = 4
)

The default factor names and formula can be inspected in the design object:

unique(rcbd$deStruct$fxTrms$fixedfr)
#>   facA facB
#> 1    1    1
#> 2    2    1
#> 3    1    2
#> 4    2    2
rcbd$deStruct$formula
#> ~facA * facB + (1 | block)
#> <environment: 0x7fab420bb040>

Treatment factors and factor levels can be renamed and relabeled, and the model formula can be changed. For example, if interactions are removed, marginal means for both factors must be provided instead of cell means.

designRCBD(treatments = c(2, 2), 
           label = list(factorA = c("A1", "A2"), factorB = c("B1", "B2")), 
           blocks = 8, 
           formula = ~ factorA + factorB + (1|block), 
           template = TRUE)
#> $fixeff
#> $fixeff$beta
#> (Intercept)   factorAA2   factorBB2 
#>           1           2           3 
#> 
#> $fixeff$means
#> factorAA1 factorAA2 factorBB1 factorBB2 
#>         1         2         3         4 
#> 
#> 
#> $varcov
#> $varcov$block
#>             (Intercept)
#> (Intercept)           1

Step 2: Evaluate statistical power

To test overall effects of facA and facB

pwr.anova(rcbd)
#> Power of type III F-test
#>           NumDF DenDF sig.level   power
#> facA          1    21      0.05 0.99969
#> facB          1    21      0.05 0.76950
#> facA:facB     1    21      0.05 0.27138

To test differences between levels of facA, conditioned on facB:

pwr.contrast(rcbd, which = "facA", by = "facB")
#> $`facB = 1`
#>               effect df sig.level     power alternative
#> facA1 - facA2     -5 21      0.05 0.9974502   two.sided
#> 
#> $`facB = 2`
#>               effect df sig.level     power alternative
#> facA1 - facA2     -3 21      0.05 0.8160596   two.sided

Example 3. Latin square design

In this example, we extend Example 2 by introducing another blocking factor, transforming the design into an LSD.

• The treatment structure and effect sizes remain the same as in Example 2.

• The LSD controls for two sources of variability (row and column blocks) while evaluating treatment effects.

• The total variance (15) is further decomposed into:

  • Variance between row blocks: 11

  • Variance between column blocks: 2

  • Residual error variance: 2

Step 1: Creating the LSD

Generate the templates to determine the correct order of inputs.

designLSD(
  treatments = c(2, 2),
  squares = 4,
  reuse = "both",
  template = TRUE
)
#> $fixeff
#> $fixeff$beta
#> (Intercept)       facA2       facB2 facA2:facB2 
#>           1           2           3           4 
#> 
#> $fixeff$means
#> facA1:facB1 facA2:facB1 facA1:facB2 facA2:facB2 
#>           1           2           3           4 
#> 
#> 
#> $varcov
#> $varcov$row
#>             (Intercept)
#> (Intercept)           1
#> 
#> $varcov$col
#>             (Intercept)
#> (Intercept)           2

Either beta or means can be provided, as in Example 2. Note that although the variance component template sets the order for row and col variances in vcomp, these values don’t impact power as long as sigma2 is fixed because treatment effects are tested against residual error. However, in designs where the error terms for treatment factors include variance components beyond residual error (like the next example), variance components can affect power for main plot factors. It’s also recommended to use variance components as a guide to make an informed estimate of sigma2.

Create the LSD

lsd <- designLSD(
  treatments = c(2, 2),
  label = list(temp = c("T1", "T2"), dosage = c("D1", "D2")),
  squares = 4,
  reuse = "both",
  means = c(35, 40, 38, 41),
  vcomp = c(11, 2),
  sigma2 = 2
)

Once the design has been created, pwr.anova and pwr.contrast can be used to evaluate statistical power as demonstrated above.

Example 4: Split-plot Design

In this example, we create an SPD with two treatment factors, one at the main plot level and another at the subplot level. The design parameters are as follows:

1. Treatments:

   • Main plot factor with 2 levels

   • Subplot factor with 3 levels

2. Replicates:

   • Each main plot treatment has 5 plots, for a total of 10 plots.

   • This is a standard SPD, where the plots are the blocks at the subplot level and the subplot design follows an RCBD structure.

3. Expected cell means are: \[ \begin{array}{c|c|c} & trt.sub1 & trt.sub2 & trt.sub3 \\ \hline trt.main1 & 20 & 22 & 24\\ trt.main2 & 22 & 24 & 28 \\ \end{array} \]

4. Variance Components:

   • Main-plot error: 4

   • Subplot (residual) error: 11

   • Total variance: 15

Generate input templates

designSPD(
  trt.main = 2,
  trt.sub = 3, 
  replicates = 10, 
  template = TRUE
)
#> $fixeff
#> $fixeff$beta
#>        (Intercept)          trt.main2           trt.sub2           trt.sub3 
#>                  1                  2                  3                  4 
#> trt.main2:trt.sub2 trt.main2:trt.sub3 
#>                  5                  6 
#> 
#> $fixeff$means
#> trt.main1:trt.sub1 trt.main2:trt.sub1 trt.main1:trt.sub2 trt.main2:trt.sub2 
#>                  1                  2                  3                  4 
#> trt.main1:trt.sub3 trt.main2:trt.sub3 
#>                  5                  6 
#> 
#> 
#> $varcov
#> $varcov$mainplot
#>             (Intercept)
#> (Intercept)           1

Create the SPD

spd <- designSPD(
  trt.main = 2,
  trt.sub = 3, 
  replicates = 10, 
  means = c(20, 22, 22, 24, 24, 28),
  vcomp = 4,
  sigma2 = 11
)

One can verify that different values of main-plot error (vcomp) affect the power of the main plot factor, unlike in Example 3, where variance components did not influence power.

Example 5: Repeated measures

This example illustrates a repeated measures design, where three treatments (CON, TRT1, TRT2) are measured hourly over 8 hours. Within-subject correlations are modeled using an AR(1) structure with \(\rho = 0.6\) and \(\sigma^2 = 2\).

Design Details

1. Subjects: 6 per treatment group (total: 18 subjects).
2. Treatments: CON, TRT1, and TRT2.
3. Time Points: 8 hourly measurements.

4. Means:\[ \begin{array}{c|cccccccc} \textbf{Treatment} & \textbf{1} & \textbf{2} & \textbf{3} & \textbf{4} & \textbf{5} & \textbf{6} & \textbf{7} & \textbf{8} \\ \hline \text{CON} & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\ \text{TRT1} & 2.50 & 3.50 & 3.98 & 4.03 & 3.68 & 3.35 & 3.02 & 2.94 \\ \text{TRT2} & 3.50 & 4.54 & 5.80 & 5.84 & 5.49 & 4.71 & 4.08 & 3.78 \\ \end{array} \]

Step 1: Creating the design

Create a data frame for the design

n_subject = 6 # Subjects per treatment
n_trt = 3 # Number of treatments
n_hour = 8 # Number of repeated measures (time points)
trt = c("CON", "TRT1", "TRT2")

df.rep <- data.frame(
  subject = as.factor(rep(seq_len(n_trt*n_subject), each = n_hour)),
  hour = as.factor(rep(seq_len(n_hour), n_subject*n_trt)),
  trt = rep(trt, each = n_subject*n_hour)
)

Input templates

Either values of beta or means are required in the following order:

mkdesign(formula = ~ trt*hour, data = df.rep)
#> $fixeff
#> $fixeff$beta
#>   (Intercept)       trtTRT1       trtTRT2         hour2         hour3 
#>             1             2             3             4             5 
#>         hour4         hour5         hour6         hour7         hour8 
#>             6             7             8             9            10 
#> trtTRT1:hour2 trtTRT2:hour2 trtTRT1:hour3 trtTRT2:hour3 trtTRT1:hour4 
#>            11            12            13            14            15 
#> trtTRT2:hour4 trtTRT1:hour5 trtTRT2:hour5 trtTRT1:hour6 trtTRT2:hour6 
#>            16            17            18            19            20 
#> trtTRT1:hour7 trtTRT2:hour7 trtTRT1:hour8 trtTRT2:hour8 
#>            21            22            23            24 
#> 
#> $fixeff$means
#>  trtCON:hour1 trtTRT1:hour1 trtTRT2:hour1  trtCON:hour2 trtTRT1:hour2 
#>             1             2             3             4             5 
#> trtTRT2:hour2  trtCON:hour3 trtTRT1:hour3 trtTRT2:hour3  trtCON:hour4 
#>             6             7             8             9            10 
#> trtTRT1:hour4 trtTRT2:hour4  trtCON:hour5 trtTRT1:hour5 trtTRT2:hour5 
#>            11            12            13            14            15 
#>  trtCON:hour6 trtTRT1:hour6 trtTRT2:hour6  trtCON:hour7 trtTRT1:hour7 
#>            16            17            18            19            20 
#> trtTRT2:hour7  trtCON:hour8 trtTRT1:hour8 trtTRT2:hour8 
#>            21            22            23            24 
#> 
#> 
#> $varcov
#> NULL

Create the design

design.rep <- mkdesign(
formula = ~ trt*hour,
data = df.rep,
means =  c(1, 2.50, 3.5,
           1, 3.50, 4.54,
           1, 3.98, 5.80,
           1, 4.03, 5.4,
           1, 3.68, 5.49,
           1, 3.35, 4.71,
           1, 3.02, 4.08,
           1, 2.94, 3.78),
sigma2 = 2,
correlation = corAR1(value = 0.6, form = ~ hour|subject)
)

Step 2: Power calculation

Statistical power for main effects of treatment and time, and their interaction:

pwr.anova(design.rep)
#> Power of type III F-test
#>          NumDF  DenDF sig.level   power
#> trt          2 21.563      0.05 1.00000
#> hour         7 86.055      0.05 0.74687
#> trt:hour    14 86.055      0.05 0.38500

Statistical power for treatment differences at each time point:

pwr.contrast(design.rep, which = "trt", by = "hour", contrast = "trt.vs.ctrl", p.adj = TRUE)[1:2]
#> $`hour = 1`
#>                  effect       df sig.level     power alternative
#> trtTRT1 - trtCON    1.5 64.41176     0.025 0.3299823   two.sided
#> trtTRT2 - trtCON    2.5 64.41176     0.025 0.7765112   two.sided
#> 
#> $`hour = 2`
#>                  effect       df sig.level     power alternative
#> trtTRT1 - trtCON   2.50 64.41176     0.025 0.7765112   two.sided
#> trtTRT2 - trtCON   3.54 64.41176     0.025 0.9777118   two.sided

As shown in Example 5, any design can be created using the mkdesign function as long as the appropriate data frame is provided. However, the package currently does not support non-normal response variables.

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})^T [K\hat{C}K^T]^{-1} (K\hat{\beta})}{v_1} \]

\(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).

However, for more general cases, \(v_2\) must be approximated using methods.

The Satterthwaite approximation (Satterthwaite, 1946) for DF in t-tests can be calculated as outlined by Giesbrecht and Burns (1985):

\[ v_2 = \frac{2(k^T \hat{C} k)^2}{g^T A g} \]

where: \(g\) is the gradient of \(k^T C(\hat{\theta}) k\) with respect to \(\theta\), the variance-covariance parameters in \(V\), evaluated at \(\hat{\theta}\). Matrix \(A\) is the asymptotic variance-covariance matrix of \(\hat{\theta}\), obtained from the information matrix of ML or REML estimation of \(\hat{\theta}\) (Stroup, 2012).

In F-tests, \(v_2\) can be calculated by following the procedures described by Fai and Cornelius (1996). First, \(KC\hat{K}^T\) is decomposed to yield \(KC\hat{K}^T = P^T D P\), where \(P\) is an orthogonal matrix of eigen vectors, and \(D\) is a diagonal matrix of eigenvalues. Define \(k_m Ck_m^T\), where \(k_m\) is the \(m\)-th row of \(P K\), and let:

\[ v_m = \frac{2(D_m)^2}{g_m^T A g_m} \]

where: \(D_m\)is the \(m\)th diagonal element of \(D\), \(g_m\) is the gradient of \(k_m C k_m^T\) with respect to \(\theta\), evaluated at \(\hat{\theta}\).

Then let:

\[ E = \sum_{m=1}^{v_1} \frac{v_m}{v_m - 2} I(v_m > 2) \]

where \(I(v_m > 2)\) denotes the indicator function. The denominator DF \(v_2\) is calculated as:

\[ v_2 = \frac{2E}{E - v_1} \] The Satterthwaite approximation can be applied in power analysis by plugging in the assumed values of variance parameters (Stroup, 2002).

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})^T [K\hat{C}K^T]^{-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.

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.