These functions are used to create design objects for the further evaluation of statistical power.
Usage
designCRD(treatments, label, replicates, formula, beta, sigma2)
designRCBD(treatments, label, blocks, formula, beta, VarCov, sigma2, ...)
designLSD(
treatments,
label,
squares = 1,
reuse = c("row", "col", "both"),
formula,
beta,
VarCov,
sigma2,
...
)
designCOD(treatments, label, squares = 1, formula, beta, VarCov, sigma2, ...)
designSPD(
trt.main,
trt.sub,
label,
replicates,
formula,
beta,
VarCov,
sigma2,
...
)
designCustom(design.df, formula, beta, VarCov, sigma2, design.name, ...)
Arguments
- treatments
An integer-valued vector specifying the treatment structure, in which the length of the vector indicates the number of treatment factors, and each value represents the number of levels for each factor. A maximum of two factors is allowed, and they are arranged in a factorial design. For instance,
treatments = n
specifies one treatment factor with n levels, andtreatments=c(2,3)
creates a "2x3" factorial design of two treatment factors with 2 and 3 levels, respectively.- label
Optional. A list of character vectors specifying the names of treatment factors and factor levels. Each vector in the list represents a treatment factor, where the name of the vector specifies the name of the factor, and the values in the vector are the labels for that factor's levels. If not provided, factors and levels for one and two treatment factors are labeled as
list(trt = c("1", "2", ...))
andlist(facA = c("1", "2", ...), facB = c("1", "2", ...))
, respectively.- replicates
The number of experimental units per treatment in a completely randomized design or the number of experimental units (main plots) per treatment of main plot factors.
- formula
A model formula for testing treatment effects in post-experimental data analysis. Use the syntax of
lm
for fixed effects and lmer for random effects. The response variable is always denoted asy
. By default, all interaction terms between treatment factors are included in the formula.- beta
A numeric vector of expected model coefficients, representing the effect sizes. The first element represents the intercept term, corresponding to the mean of the reference level for categorical variables. Subsequent elements correspond to the effect sizes of the independent variables in the order they appear in the model matrix. For categorical variables, each coefficient represents the difference between a non-reference level and the reference level (intercept), as
contr.treatment
contrast coding is used for constructing the model matrix. Ensure thatbeta
aligns with the columns of the model matrix, including any dummy variables created for categorical predictors.- sigma2
error variance.
- blocks
The number of blocks.
- VarCov
Variance-covariance components of random effects. For multiple random effect groups, supply the variance (for a single random effect term) or variance-covariance matrix (for two or more random effect terms) of each group in a list, following the order in the model formula.
- ...
Additional arguments passed to the
anova
function inlmerTest
. The type of ANOVA table (default is Type III) and the method for computing denominator degrees of freedom (default is Satterthwaite's method) can be modified. For balanced designs, the choice of sum of squares (SS) and degrees of freedom (df) does not affect the results.- squares
The number of replicated squares. By default, 1, i.e., no replicated squares.
- reuse
A character string specifying how to replicate squares when there are multiple squares. Options are: "row" for reusing row blocks, "col" for reusing column blocks, or "both" for reusing both row and column blocks to replicate a single square.
- trt.main
An integer-valued vector specifying the treatment structure at main plot level for a split plot design, similar to
treatments
.- trt.sub
An integer-valued vector specifying the treatment structure at sub plot level for a split plot design, similar to
treatments
.- design.df
Required input for creating a customized design. A data frame with all independent variables of the design as columns, representing the actual data structure (long format data frame) without response variables.
- design.name
Optional input for creating a customized design. A character.
Value
a list with the design name, data structure (data frame), model formula, and a pseudo model object with the expected fixed and random effects.
Details
Each function creates a specific design as described below:
designCRD
Completely Randomized Design. By default, the model formula is
y ~ trt
for one factor andy ~ facA*facB
for two factors, unless explicitly specified. If thelabel
argument is provided, the formula is automatically updated with the specified treatment factor names.designRCBD
Randomized Complete Block Design. The default model formula is
y ~ trt + (1|block)
for one factor andy ~ facA*facB + (1|block)
for two factors. Iflabel
is provided, the fixed effect parts of the formula are automatically updated with the specified names. The label of block factor ("block") in the formula is not changeable.designLSD
Latin Square Design. The default formula is
y ~ trt + (1|row) + (1|col)
for one factor andy ~ facA*facB + (1|row) + (1|col)
for two factors. Iflabel
is provided, the fixed effect parts of the formula are automatically updated with the specified names. The labels 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
y ~ trt + (1|subject) + (1|period)
for one factor andy ~ facA*facB + (1|subject) + (1|period)
for two factors. Iflabel
is provided, the fixed effect parts of the formula are automatically updated with the specified names. Note that "subject" and "period" are the labels 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:
y ~ . + (1|mainplot)
. Iflabel
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 denoted as "mainplot".designCustom
Customized Design.
Examples
# Example 1: 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
# mean of level1, and the means of other levels minus level1, respectively
beta = c(30, -2, 3, 5),
sigma2 = 10 # error variance
)
## power of omnibus test
pwr.anova(crd)
#> Power analysis of Completely Randomized Design
#> NumDF DenDF non-centrality alpha power
#> trt 3 44 34.8 0.05 0.999
## power of contrast
pwr.contrast(crd, specs = "trt", method = "pairwise") # pairwise comparisons
#> contrast estimate df non-centrality alpha power
#> 1 trt1 - trt2 2 44 2.4 0.05 0.3285822
#> 2 trt1 - trt3 -3 44 5.4 0.05 0.6228308
#> 3 trt1 - trt4 -5 44 15.0 0.05 0.9661638
#> 4 trt2 - trt3 -5 44 15.0 0.05 0.9661638
#> 5 trt2 - trt4 -7 44 29.4 0.05 0.9995822
#> 6 trt3 - trt4 -2 44 2.4 0.05 0.3285822
pwr.contrast(crd, specs = "trt", method = "poly") # polynomial contrasts
#> contrast estimate df non-centrality alpha power
#> 1 linear 20 44 24.0 0.05 0.9976700
#> 2 quadratic 4 44 4.8 0.05 0.5726028
#> 3 cubic -10 44 6.0 0.05 0.6685119
# Example 2: Evaluate the power of an RCBD with 2 x 2 factorial treatments
# Treatment factors are A (A1 vs. A2) and B (B1 vs. B2).
# To illustrate how to provide `beta`, treatment means are presented:
# B1 B2
# A1 20 24
# A2 17 22
#
# From these means, we calculate:
# 1. the mean of reference level (A1B1): 20
# 2. the effect of A2 alone: Effect_A2 = A2B1 - A1B1 = 17 - 20 = -3
# 3. the effect of B2 alone: Effect_A2 = A1B2 - A1B1 = 24 - 20 = 4
# 4. the interaction effect of A2 and B2:
# Interaction_A2B2 = A2B2 - A2B1 - A1B2 + A1B1 = 22 - 17 - 24 + 20 = 1, representing
# the additional effect of combining A2B2 compared to what would be expected
# from the sum of individual effects of A2 and B2.
# The `beta` vector is constructed as:
# beta = c(mean_A1B1, Effect_A2, Effect_B2, Interaction_A2B2)
# beta = c(20, -3, 4, 1)
## Create a design object
rcbd <- designRCBD(
# 2x2 factorial design
treatments = c(2, 2),
# Specify treatment names
label = list(A = c("A1", "A2"), B = c("B1", "B2")),
# 12 blocks, totaling 48 experimental units
blocks = 12,
# Mean of the reference level and effect sizes as calculated above
beta = c(20, -3, 4, 1),
# Variance of block effects (between-block variance)
VarCov = 30,
# Error variance (within-block variance)
sigma2 = 20
)
## power of omnibus test
pwr.anova(rcbd)
#> Power Analysis of Randomized Complete Block Design
#> NumDF DenDF non-centrality alpha power
#> A 1 33 3.75 0.05 0.46821
#> B 1 33 12.15 0.05 0.92258
#> A:B 1 33 0.15 0.05 0.06636
## power of B2 vs. B1 at each level of A
pwr.contrast(rcbd, specs = ~B|A, method = "pairwise")
#> contrast A estimate df non-centrality alpha power
#> 1 B1 - B2 A1 -4 33 4.8 0.05 0.5663008
#> 2 B1 - B2 A2 -5 33 7.5 0.05 0.7574784
# More examples are available in the package vignette("pwr4exp")
# and on the package website: https://an-ethz.github.io/pwr4exp/