Movatterモバイル変換

Title:	Power Analysis for Research Experiments
Version:	1.0.1
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-11-11 11:20:28 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-11-11 11:50:02 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

AR(1) Correlation Structure

Description

Re-exportsnlme::corAR1 from thenlme package.

Usage

corAR1(value, form, fixed)

Arguments

Details

In the originalnlme::corAR1 function, a covariatet must be an integer class.Inpwr4exp,t can also be a factor class, which is then converted to an integer internally for chronological order.The class oft in the model formula, if present, is not converted. For example, a time covariate can be fitted as a factorin the model formula, whereas it is converted to an integer in the correlation formula.

See Also

corClassesSeenlme::corAR1 for original documentation.

ARMA(p,q) Correlation Structure

Description

Re-exportsnlme::corARMA from thenlme package.

Usage

corARMA(value, form, p, q, fixed)

Arguments

Details

In the originalnlme::corARMA function, a covariatet must be an integer class.Inpwr4exp,t can also be a factor class, which is then converted to an integer internally for chronological order.The class oft in the model formula, if present, is not converted. For example, a time covariate can be fitted as a factorin the model formula, whereas it is converted to an integer in the correlation formula.

See Also

corClassesSeenlme::corARMA for original documentation.

Continuous AR(1) Correlation Structure

Description

Re-exportsnlme::corCAR1 from thenlme package.

Usage

corCAR1(value, form, fixed)

Arguments

See Also

corClassesSeenlme::corCAR1 for original documentation.

Correlation Structure Classes

Description

Standard classes of correlation structures (corStruct) available from thenlme packageand re-exported bypwr4exp for convenience when specifying correlation structures inmkdesign.

All arguments are identical to the correspondingnlme functions. For more details on the original implementations,seenlme::corClasses.

Note: In the originalnlme::corAR1,nlme::corARMA, andnlme::corSymm functions,the covariatet in the correlation formula~ t or~ t | g must be an integer class. Inpwr4exp,the covariate can also be a factor class, which is then converted to an integer internally for sorting purposes.The class of the covariate variable in the model formula, if present, will not be converted. For example, a time covariatecan be fitted as a factor in the model formula, whereas it is converted to an integer in the correlation formula temporarily for matrix sorting.

Details

Available standard classes:

corAR1

autoregressive process of order 1.

corARMA

autoregressive moving average process, with arbitrary ordersfor the autoregressive and moving average components.

corCAR1

continuous AR(1)

corCompSymm

compound symmetry structure corresponding to a constant correlation.

corExp

exponential spatial correlation.

corGaus

Gaussian spatial correlation.

corLin

linear spatial correlation.

corRatio

Rational quadratics spatial correlation.

corSpher

spherical spatial correlation.

corSymm

general correlation matrix, with no additional structure.

References

Pinheiro, J.C., and Bates, D.M. (2000) "Mixed-Effects Models in S and S-PLUS", Springer.

Compound Symmetry Correlation Structure

Description

Re-exportsnlme::corCompSymm from thenlme package.

Usage

corCompSymm(value, form, fixed)

Arguments

See Also

corClassesSeenlme::corCompSymm for original documentation.

Exponential Correlation Structure

Description

Re-exportsnlme::corExp from thenlme package.

Usage

corExp(value, form, nugget, metric, fixed)

Arguments

See Also

corClassesSeenlme::corExp for original documentation.

Gaussian Correlation Structure

Description

Re-exportsnlme::corGaus from thenlme package.

Usage

corGaus(value, form, nugget, metric, fixed)

Arguments

See Also

corClassesSeenlme::corGaus for original documentation.

Linear Correlation Structure

Description

Re-exportsnlme::corLin from thenlme package.

Usage

corLin(value, form, nugget, metric, fixed)

Arguments

See Also

corClassesSeenlme::corLin for original documentation.

Rational Quadratics Correlation Structure

Description

Re-exportsnlme::corRatio from thenlme package.

Usage

corRatio(value, form, nugget, metric, fixed)

Arguments

See Also

corClassesSeenlme::corRatio for original documentation.

Spherical Correlation Structure

Description

Re-exportsnlme::corSpher from thenlme package.

Usage

corSpher(value, form, nugget, metric, fixed)

Arguments

See Also

corClassesSeenlme::corSpher for original documentation.

General Correlation Structure

Description

Re-exportsnlme::corSymm from thenlme package.

Usage

corSymm(value, form, fixed)

Arguments

See Also

corClassesSeenlme::corSymm for original documentation.

Creation of Standard Experimental Designs

Description

These functions facilitate the creation of standard experimental designscommonly used in agricultural studies for power analysis. Unlikemkdesignwhich requires a pre-existing data frame, these functions allow users todirectly specify key design characteristics to generate experimental layouts.Quantitative parameters describing fixed and random effects remain consistentwith those inmkdesign.

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", "none"),  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 thelabel argument is provided, the formula is automatically updated with thespecified treatment factor names.

designRCBD

Randomized Complete Block Design.Experimental units are grouped into blocks, with each treatment appearingexactly 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. Iflabel is provided, thefixed effect parts of the formula are automatically updated with the specifiednames. 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. Iflabel is provided,the fixed effect parts of the formula are automatically updated with the specifiednames. The names of row ("row") and column ("col") block factors are not changeable.

designCOD

Crossover Design, which is a special case of LSDwith time periods and individuals as blocks. Period blocks are reused whenreplicating squares.The default formula is~ trt + (1|subject) + (1|period) for one factorand~ facA*facB + (1|subject) + (1|period) for two factors. Iflabelis provided, the fixed effect parts of the formula are automatically updatedwith the specified names. Note that "subject" and "period" are the names forthe two blocking factors and cannot be changed.

designSPD

Split Plot Design.The default formula includes the main effects of all treatment factors atboth the main and sub-plot levels, their interactions, and the random effectsof main plots:~ . + (1|mainplot). Iflabel is provided, the fixedeffect 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 thesubplot 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 matricesfor fixed and random effects, variance-covariance matrices for random effectsand residuals, etc.

• Internally calculated higher-level parameters (deParam), including variance-covariancematrix 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 objectcrd <- 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 testpwr.anova(crd)## power of contrastpwr.contrast(crd, which = "trt", contrast = "pairwise") # pairwise comparisonspwr.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", "none"))

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 individualterms in the model, which is why the argument is calledterm and nota name likeform, 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, correlation = NULL)

Arguments

Value

A list containing:

• corStruct: An initialized correlation structure.

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

• R: A block-diagonal residual variance-covariance structure, not yet scaledby 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 withvariance-covariance structures for random effects (G-side) and residuals (R-side).

Usage

mkStruct(formula, data, correlation)

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 internallycalculated 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 packagessuch asagricolae,AlgDesign,crossdes, orFrF2.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 iftemplate = TRUE. Templates serve as guides for specifying inputs:

  • Template forbeta: Represents the sequence of model coefficients.

  • Template formeans: 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 variableat 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 formeans andbeta are identical.See the "Examples" for illustrative scenarios.

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

• correlation: Various residual correlation structures can be specified following the instructions fromcorClasses constructors.

  Note: In the originalcorAR1,corARMA, andcorSymm functions,the covariatet in the correlation formula~ t or~ t | g must be an integer class. Inpwr4exp,the covariate can also be a factor class, which is then converted to an integer internally for sorting purposes.The class of the covariate variable in the model formula, if present, will not be converted. For example, a time covariatecan be fitted as a factor in the model formula, whereas it is converted to an integer in the correlation formula temporarily for matrix sorting.

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-covariancematrix 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 variablesdf1 <- 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 xdf1$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 requiredmkdesign(~ 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 zmkdesign(~ 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 measuresn_subject <- 6n_trt <- 3n_hour <- 8trt <- 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)## Templatestemp <- 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.6design <- 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 requireddf2$hour <- as.integer(df2$hour)mkdesign(formula = ~ trt * hour, data = df2)$fixeff$means## Polynomial terms of time in the modelmkdesign(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 individualterms in the model, which is why the argument is calledterm and nota name likeform, indicating a formula.

Power of omnibus tests

Description

Calculates the statistical power for testing the overall effects of treatmentfactors 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), denominatordegrees of freedom (DenDF), type I error rate (sig.level), andpower.

See Also

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

Examples

# generate an RCBDrcbd <- 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 testpwr.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 eachby 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(byalternative).When multipleby 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 conductedpwr.contrast(rcbd, which = "facA") # Marginal effect of facApwr.contrast(rcbd, which = "facA", by = "facB") # Conditional effect of facA within levels of facB# Custom contrast vector, identical to pairwise comparisonpwr.contrast(rcbd, which = "facA", contrast = c(1, -1))pwr.contrast(rcbd, which = "facA", by = "facB", contrast = c(1, -1))# A single factor combining two factorspwr.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-modelformula.

Usage

subbars(term)

Arguments

Value

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

Note

This function is called recursively on individualterms in the model, which is why the argument is calledterm and nota name likeform, indicating a formula.