Title:	Power Analysis for Research Experiments
Version:	1.0.0
Description:	Provides tools for calculating statistical power for experiments analyzed using linear mixed models. It supports standard designs, including randomized block, split-plot, and Latin Square designs, while offering flexibility to accommodate a variety of other complex study designs.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	https://github.com/an-ethz/pwr4exp, https://an-ethz.github.io/pwr4exp/
BugReports:	https://github.com/an-ethz/pwr4exp/issues
Depends:	R (≥ 2.10), stats
Suggests:	agricolae, AlgDesign, crossdes, FrF2, knitr, rmarkdown
VignetteBuilder:	knitr
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.2
Imports:	emmeans, MASS, Matrix, methods, nlme, numDeriv
NeedsCompilation:	no
Packaged:	2025-03-17 10:24:20 UTC; wangkai
Author:	Kai Wang [aut, cre, cph], Mutian Niu [aut, cph]
Maintainer:	Kai Wang <kai.wang@usys.ethz.ch>
Repository:	CRAN
Date/Publication:	2025-03-17 10:40:03 UTC

Computes power of t-test for one-dimensional contrast matrices

Description

Computes power of t-test for one-dimensional contrast matrices

Usage

contrast1D(
  object,
  L,
  method = c("Satterthwaite"),
  sig.level = 0.05,
  alternative = c("two.sided", "one.sided"),
  strict = TRUE
)

Arguments

Value

A data frame with columns for effect size, degrees of freedom, significance level, power, and test type.

Computes power of F-test for multi-dimensional contrast matrices

Description

Computes power of F-test for multi-dimensional contrast matrices

Usage

contrastMD(object, L, sig.level = 0.05, eps = sqrt(.Machine$double.eps))

Arguments

Value

A data frame

Creation of Standard Experimental Designs

Description

These functions facilitate the creation of standard experimental designs commonly used in agricultural studies for power analysis. Unlike mkdesign which requires a pre-existing data frame, these functions allow users to directly specify key design characteristics to generate experimental layouts. Quantitative parameters describing fixed and random effects remain consistent with those in mkdesign.

Usage

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

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

designLSD(
  treatments,
  label,
  squares = 1,
  reuse = c("row", "col", "both"),
  formula,
  beta = NULL,
  means = NULL,
  vcomp,
  sigma2,
  template = FALSE,
  REML = TRUE
)

designCOD(
  treatments,
  label,
  squares = 1,
  formula,
  beta = NULL,
  means = NULL,
  vcomp,
  sigma2,
  template = FALSE,
  REML = TRUE
)

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

Arguments

Details

Each function creates a standard design as described below:

designCRD

Completely Randomized Design. By default, the model formula is ~ trt for one factor and ~ facA*facB for two factors, unless explicitly specified. If the label argument is provided, the formula is automatically updated with the specified treatment factor names.

designRCBD

Randomized Complete Block Design. Experimental units are grouped into blocks, with each treatment appearing exactly once per block (i.e., no replicates per treatment within a block). The default model formula is ~ trt + (1|block) for one factor and ~ facA*facB + (1|block) for two factors. If label is provided, the fixed effect parts of the formula are automatically updated with the specified names. The block factor is named "block" and not changeable.

designLSD

Latin Square Design. The default formula is ~ trt + (1|row) + (1|col) for one factor and ~ facA*facB + (1|row) + (1|col) for two factors. If label is provided, the fixed effect parts of the formula are automatically updated with the specified names. The names of row ("row") and column ("col") block factors are not changeable.

designCOD

Crossover Design, which is a special case of LSD with time periods and individuals as blocks. Period blocks are reused when replicating squares. The default formula is ~ trt + (1|subject) + (1|period) for one factor and ~ facA*facB + (1|subject) + (1|period) for two factors. If label is provided, the fixed effect parts of the formula are automatically updated with the specified names. Note that "subject" and "period" are the names for the two blocking factors and cannot be changed.

designSPD

Split Plot Design. The default formula includes the main effects of all treatment factors at both the main and sub-plot levels, their interactions, and the random effects of main plots: ~ . + (1|mainplot). If label is provided, the fixed effect parts of the formula are automatically updated with the specified names. The experimental unit at the main plot level (i.e., the block factor at the subplot level) is always named as "mainplot".

Value

A list object containing all essential components for power calculation. This includes:

• Structural components (deStruct): including the data frame, design matrices for fixed and random effects, variance-covariance matrices for random effects and residuals, etc.

• Internally calculated higher-level parameters (deParam), including variance-covariance matrix of beta coefficients (vcov_beta), variance-covariance matrix of variance parameters (vcov_varpar), gradient matrices (Jac_list), etc.

See Also

mkdesign, pwr.anova, pwr.contrast

Examples

# Evaluate the power of a CRD with one treatment factor
## Create a design object
crd <- designCRD(
  treatments = 4, # 4 levels of one treatment factor
  replicates = 12, # 12 units per level, 48 units totally
  means = c(30, 28, 33, 35), # means of the 4 levels
  sigma2 = 10 # error variance
)

## power of omnibus test
pwr.anova(crd)

## power of contrast
pwr.contrast(crd, which = "trt", contrast = "pairwise") # pairwise comparisons
pwr.contrast(crd, which = "trt", contrast = "poly") # polynomial contrasts

# More examples are available in `vignette("pwr4exp")`
# and on https://an-ethz.github.io/pwr4exp/

Create a data frame for Crossover design

Description

Create a data frame for Crossover design

Usage

df.cod(treatments, label, squares)

Arguments

Value

a data.frame representing the data structure of the design

Create a data frame of completely randomized design

Description

Create a data frame of completely randomized design

Usage

df.crd(treatments, label, replicates)

Arguments

Value

a data.frame representing the data structure of the design

Create a data frame for Latin square design

Description

Create a data frame for Latin square design

Usage

df.lsd(treatments, label, squares = 1, reuse = c("row", "col", "both"))

Arguments

Value

a data.frame representing the data structure of the design

Create a data frame of randomized complete block design

Description

Create a data frame of randomized complete block design

Usage

df.rcbd(treatments, label, blocks)

Arguments

Value

a data.frame representing the data structure of the design

Create data frame for split-plot design

Description

Create data frame for split-plot design

Usage

df.spd(trt.main, trt.sub, label, replicates)

Arguments

Value

a data.frame representing the data structure of the design

Expand terms with '||' notation into separate '|' terms

Description

Expand terms with '||' notation into separate '|' terms

Usage

expandDoubleVerts(term)

Arguments

Value

the modified term

Convert variables to factors as necessary

Description

Convert variables to factors as necessary

Usage

factorize(x, frloc, char.only = FALSE)

Arguments

Value

a copy of the data frame with factors converted

Determine random-effects expressions from a formula

Description

Determine random-effects expressions from a formula

Usage

findbars(term)

Arguments

Value

pairs of expressions that were separated by vertical bars

Note

This function is called recursively on individual terms in the model, which is why the argument is called term and not a name like form, indicating a formula.

An exemplary dataset of a 4x4 crossover design with 2 squares

Description

Milk yield records from 8 cows over 4 different periods in a 4x4 crossover design. The design includes 2 Latin squares, each consisting of 4 cows and 4 periods.

Usage

milk

Format

A data frame with 32 rows and 4 variables:

Cow

Factor: Cow index (8 levels)

Period

Factor: Period index (4 levels)

Treatment

Factor: Treatment index (4 levels)

MilkYield

Numeric: milk yield recordings (in kg)

Source

Simulated data for package demonstration purposes.

Residual Variance-Covariance Matrices

Description

Residual Variance-Covariance Matrices

Usage

mkRTrms(data, corcall = NULL)

Arguments

Value

A list containing:

• corcall

• corframe: A processed data frame with indexed grouping variables and ordering for correlation structures.

• R: A block-diagonal residual variance-covariance structure, not yet scaled by the residual variance

Design Matrices and Variance Components for Random Effects

Description

Adapted from lme4, this function constructs the design matrix (Z), variance-covariance matrix (G), etc.

Usage

mkReTrms(
  bars,
  fr,
  drop.unused.levels = TRUE,
  reorder.terms = FALSE,
  reorder.vars = FALSE
)

Arguments

Value

A list with the following components:

• Zt: Transposed random-effects design matrix.

• G: Variance-covariance matrix for random effects.

• Gind: Index mapping variance-covariance parameters to their positions in G.

• G_temp: List of individual variance component matrices for each random effect.

• flist: List of grouping factors used in random effects.

• cnms: Column names of the random-effects design matrix.

• Ztlist: List of per-term transposed design matrices for random effects.

• nl: Number of levels for each grouping factor.

Building Design Matrices and Covariance Structures for Linear Mixed Models

Description

Constructs design matrices for fixed and random effects, along with variance-covariance structures for random effects (G-side) and residuals (R-side).

Usage

mkStruct(formula, data, corcall)

Arguments

Value

A list containing:

• data: Processed data frame with NA values omitted.

• fxTrms: Fixed-effects design structure, including model frame and design matrix.

• reTrms: Random-effects structure (if applicable), including grouping factors, design matrices, and variance-covariance matrix.

• rTrms: Residual structure (R-side variance-covariance components).

• formula: Expanded model formula.

Create a Design Object for Power Calculation

Description

Generate a design object for power analysis by specifying a model formula and data frame. This object is not a true experimental design as created by design generation procedures, where randomization and unit allocation are performed. Instead, it serves as an object containing all necessary information for power analysis, including design matrices, assumed values of model effects, and other internally calculated parameters.

Usage

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

Arguments

Details

• data: A long-format data frame is required, as typically used in R for fitting linear models. This data frame can be created manually or with the help of design creation packages such as agricolae, AlgDesign, crossdes, or FrF2. It should include all independent variables specified in the model (e.g., treatments, blocks, subjects). All the irrelevant variables not specified in the model are ignored.

• template: Templates are automatically generated when only the formula and data are supplied, or explicitly if template = TRUE. Templates serve as guides for specifying inputs:

  • Template for beta: Represents the sequence of model coefficients.

  • Template for means: Specifies the order of means (for categorical variables) and/or regression coefficients (for continuous variables), depending on the scenario:

    • Categorical variables without interactions: Requires marginal means for each level of the categorical variable(s).

    • Interactions among categorical variables: Requires conditional (cell) means for all level combinations.

    • Numerical variables without interactions: Requires regression coefficients. The intercept must also be included if there are no categorical variables in the model.

    • Interactions among numerical variables: Requires regression coefficients for both main effects and interaction terms. The intercept must also be included if there are no categorical variables in the model.

    • Categorical-by-numerical interactions: Requires regression coefficients for the numerical variable at each level of the categorical variable, as well as marginal means for the levels of the categorical variable.

    Note: For models containing only numerical variables, the inputs for means and beta are identical. See the "Examples" for illustrative scenarios.

  • Template for vcomp: Represents a variance-covariance matrix, where integers indicate the order of variance components in the input vector.

• correlation: Various residual correlation structures can be specified following the instructions from nlme::corClasses.

  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, factor is also supported.

Value

A list object containing all essential components for power calculation. This includes:

• Structural components (deStruct): including design matrices for fixed and random effects, variance-covariance matrices for random effects and residuals, etc.

• Internally calculated higher-level parameters (deParam), including variance-covariance matrix of beta coefficients (vcov_beta), variance-covariance matrix of variance parameters (vcov_varpar), gradient matrices (Jac_list), etc.

Examples

# Using templates for specifying "means"

# Create an example data frame with four categorical variables (factors)
# and two numerical variables
df1 <- expand.grid(
  fA = factor(1:2),
  fB = factor(1:2),
  fC = factor(1:3),
  fD = factor(1:3),
  subject = factor(1:10)
)
df1$x <- rnorm(nrow(df1))  # Numerical variable x
df1$z <- rnorm(nrow(df1))  # Numerical variable z

## Categorical variables without interactions
# Means of each level of fA and fB are required in sequence.
mkdesign(~ fA + fB, df1)$fixeff$means

## Interactions among categorical variables
# Cell means for all combinations of levels of fA and fB are required.
mkdesign(~ fA * fB, df1)$fixeff$means

## Numerical variables without and with interactions, identical to beta.
# Without interactions: Regression coefficients are required
mkdesign(~ x + z, df1)$fixeff$means

# With interactions: Coefficients for main effects and interaction terms are required.
mkdesign(~ x * z, df1)$fixeff$means

## Categorical-by-numerical interactions
# Marginal means for each level of fA, and regression coefficients for x
# at each level of fA are required.
mkdesign(~ fA * x, df1)$fixeff$means

## Three factors with interactions:
# Cell means for all 12 combinations of the levels of fA, fB, and fC are required.
mkdesign(~ fA * fB * fC, df1)

# A design that mixes the above-mentioned scenarios:
# - Interactions among three categorical variables (fA, fB, fC)
# - A categorical-by-numerical interaction (fD * x)
# - Main effects for another numerical variable (z)
# The required inputs are:
# - Cell means for all combinations of levels of fA, fB, and fC
# - Means for each level of fD
# - Regression coefficients for x at each level of fD
# - Regression coefficients for z
mkdesign(~ fA * fB * fC + fD * x + z, df1)$fixeff$means

# Using templates for specifying "vcomp"

# Assume df1 represents an RCBD with "subject" as a random blocking factor.
## Variance of the random effect "subject" (intercept) is required.
mkdesign(~ fA * fB * fC * fD + (1 | subject), df1)$varcov

# Demonstration of templates for more complex random effects
## Note: This example is a demo and statistically incorrect for this data
## (no replicates under subject*fA). It only illustrates variance-covariance templates.
## Inputs required:
## - Variance of the random intercept (1st)
## - Covariance between the intercept and "fA2" (2nd)
## - Variance of "fA2" (3rd)
mkdesign(~ fA * fB * fC * fD + (1 + fA | subject), df1)$varcov

# Power analysis for repeated measures

## Create a data frame for a CRD with repeated measures
n_subject <- 6
n_trt <- 3
n_hour <- 8
trt <- c("CON", "TRT1", "TRT2")
df2 <- data.frame(
  subject = as.factor(rep(seq_len(n_trt * n_subject), each = n_hour)), # Subject as a factor
  hour = as.factor(rep(seq_len(n_hour), n_subject * n_trt)),           # Hour as a factor
  trt = rep(trt, each = n_subject * n_hour)                           # Treatment as a factor
)

## Templates
temp <- mkdesign(formula = ~ trt * hour, data = df2)
temp$fixeff$means  # Fixed effects means template

## Create a design object
# Assume repeated measures within a subject follow an AR1 process with a correlation of 0.6
design <- mkdesign(
  formula = ~ trt * hour,
  data = df2,
  means = c(1, 2.50, 3.50,
            1, 3.50, 4.54,
            1, 3.98, 5.80,
            1, 4.03, 5.40,
            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)
)

pwr.anova(design)  # Perform power analysis

## When time is treated as a numeric variable
# Means of treatments and regression coefficients for hour at each treatment level are required
df2$hour <- as.integer(df2$hour)
mkdesign(formula = ~ trt * hour, data = df2)$fixeff$means

## Polynomial terms of time in the model
mkdesign(formula = ~ trt + hour + I(hour^2) + trt:hour + trt:I(hour^2), data = df2)$fixeff$means

Omit terms separated by vertical bars in a formula

Description

Remove the random-effects terms from a mixed-effects formula, thereby producing the fixed-effects formula.

Usage

nobars(term)

Arguments

Value

the fixed-effects part of the formula

Note

This function is called recursively on individual terms in the model, which is why the argument is called term and not a name like form, indicating a formula.

Power of omnibus tests

Description

Calculates the statistical power for testing the overall effects of treatment factors and their interactions, i.e., power of F-test.

Usage

pwr.anova(object, sig.level = 0.05, type = c("III", "II", "I", "3", "2", "1"))

Arguments

Value

a data frame with numerator degrees of freedom (NumDF), denominator degrees of freedom (DenDF), type I error rate (sig.level), and power.

See Also

mkdesign, designCRD, designRCBD, designLSD, designCOD, designSPD, pwr.summary and pwr.contrast

Examples

# generate an RCBD
rcbd <- designRCBD(
  treatments = c(2, 2),
  label = list(facA = c("1", "2"), facB = c("1", "2")),
  blocks = 12,
  formula = ~ facA*facB + (1|block),
  means = c(32, 35, 30, 37),
  vcomp = 4,
  sigma2 = 6
)
# power of omnibus test
pwr.anova(rcbd)

Power of contrasts

Description

Computes the statistical power of t-tests for comparisons among means.

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

Value

For each by condition, returns a data frame containing the contrast value (effect), degrees of freedom (df), type I error rate (sig.level), power, and the test direction (by alternative). When multiple by conditions are present, the results are returned as a list.

Examples

rcbd <- designRCBD(
  treatments = c(2, 2),
  label = list(facA = c("1", "2"), facB = c("1", "2")),
  blocks = 12,
  formula = ~ facA*facB + (1|block),
  means = c(32, 35, 30, 37),
  vcomp = 4,
  sigma2 = 6
)

# If contrast is not specified, pairwise comparisons are conducted
pwr.contrast(rcbd, which = "facA") # Marginal effect of facA
pwr.contrast(rcbd, which = "facA", by = "facB") # Conditional effect of facA within levels of facB

# Custom contrast vector, identical to pairwise comparison
pwr.contrast(rcbd, which = "facA", contrast = c(1, -1))
pwr.contrast(rcbd, which = "facA", by = "facB", contrast = c(1, -1))

# A single factor combining two factors
pwr.contrast(
  rcbd,
  which = "facA*facB",
  contrast = list(
    A1B1vs.A2B1 = c(1, -1, 0, 0),
    A1B1vs.A2B2 = c(1, 0, 0, -1)
  )
)

Power for model coefficients

Description

Computes the statistical power for testing (t-test) model coefficients.

Usage

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

Arguments

Value

a data frame containing model coefficients, degrees of freedom (df), type I error rate (sig.level), power, and the test direction (alternative).

Examples

rcbd <- designRCBD(
  treatments = c(2, 2),
  label = list(facA = c("1", "2"), facB = c("1", "2")),
  blocks = 12,
  formula = ~ facA*facB + (1|block),
  means = c(32, 35, 30, 37),
  vcomp = 4,
  sigma2 = 6
)
pwr.summary(rcbd)

Substitute Bars

Description

Substitute the '+' function for the '|' and '||' function in a mixed-model formula.

Usage

subbars(term)

Arguments

Value

the formula with all | and || operators replaced by +

Note

This function is called recursively on individual terms in the model, which is why the argument is called term and not a name like form, indicating a formula.