pwr4exp is an R package for performing statistical power analysis for experiments analyzed using linear mixed models (LMM). It supports power calculations for omnibus F-tests (ANOVA-like tests) and specific contrasts (t-tests) for Gaussian response variables, employing the Satterthwaite approximation for degrees of freedom.
Key features:
mkdesign
) or specialized functions (e.g.,
designCRD
, designRCBD
, designLSD
,
designCOD
, designSPD
) for standard
designs.The workflow involves two main steps:
Creating a design object: Defining the experimental structure and specifying parameters (fixed effects, (co-)variances).
Calculating power: Computing the power of F-tests or t-tests for the effects of interest.
A design object in pwr4exp is a list
containing design matrices, fixed effects, variance-covariance
components, and other higher-level parameters created by the
design-generating functions. It is not an experimental design in the
traditional sense, where randomization and unit allocation are
performed, but rather a container for all the information necessary to
conduct power analysis.
Two approaches are provided to define a design:
mkdesign
: A design object can be created by
providing a model formula and a data frame that outlines the
experimental layout.designCRD
,
designRCBD
, designLSD
, designCOD
,
or designSPD
generate design objects for common
experimental setups. These functions allow specification of design
characteristics (e.g., number of treatments, blocks, or replicates) and
automatically handle the creation of the data frame.mkdesign
functionThe 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).
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)
Randomized complete block design (RCBD)
Latin square design (LSD)
Crossover design
This is a special case of LSD where time periods and individuals act as blocks. Period blocks are reused when replicating squares.
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
andtrt.sub
: Define main-plot and subplot treatments for SPD, following the same rule fortreatment
.
label
: (optional) A list assigning labels to factor levels. If not specified, default names are assigned. For one treatment factor, the default islist(trt = c("1", "2", ...))
. For two factors, the default islist(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 toTRUE
, templates forbeta
,means
, andvcomp
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.
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
beta
TemplateThe
beta
template displays the order of model coefficients in the fitted model.
This is particularly useful when specifying expected model
coefficients directly as effect size measures and when the model
includes complex interactions. For example, the expected values of the
model coefficients for ~ fA*fB + x
should be provided in
the following sequence:
means
TemplateFor treatment factors, it may be more more convenient to provide
means
instead of beta
. The means
template represents either marginal or cell means for factors, depending
on the presence of interactions. For numerical variables included in the
model, regression coefficients remain necessary alongside alongside the
means for factors.
Factors without interactions
If no interactions are present, marginal means for each factor level are required.
For example, for the model ~ fA + fB
, means should be
reported in the order specified by the template:
Interactions among factors
When factors interact, conditional (cell) means must be provided for all possible combinations of their levels.
For example, for a model with three interacting factors
(fA
, fB
, fC
), cell means are
required for all 12 level combinations in the following order:
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
Numerical variables
For numerical variables, regression coefficients are required. If the model only includes numerical variables, the intercept must also be included. In this case, the values in means are identical to beta:
When numerical variables interact, regression coefficients for both main effects and interaction terms must be included.
Factor-by-numerical interactions
For interactions between factor and numerical variables, both Marginal means for the factor, and Regression coefficients for the numerical variable at each level of the factor must be provided.
For instance, in the model ~ fA * x
, the means for
levels of fA
and the regression coefficients for
x
at each level of fA
must be provided:
Combining Multiple Situations
The rules outlined for factor and numerical variables apply when multiple types of variables and interactions appear in the same model.
For example, consider a model including interactions among three
factors (fA
, fB
, fC
),
factor-by-numerical interaction (fD * x
), and main effects
for another numerical variable (z
).
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
The required elements and their order in means are:
1. Means for each level of fD
(positions 1–3)
2. Regression coefficient for z
(position 4)
3. Regression coefficients for x
at each level of
fD
(positions 5–7)
4. Cell means for all combinations of levels of fA
,
fB
, and fC
(positions 8–19)
vcomp
TemplateThe
vcomp
template represents variance-covariance matrices, where integer values indicate the order of variance components in the input vector.
For a single random effect, a template is unnecessary since it corresponds to a single variance value. For multiple random effects, the template specifies the sequence of variance and covariance components, ensuring proper alignment when specifying variance components in the model.
For example, consider a model that includes both a random intercept
and a random slope for fA
by subject
:
mkdesign(~ fA * fB * fC * fD + (1 + fA | subject), df1)$varcov
#> $subject
#> (Intercept) fA2
#> (Intercept) 1 2
#> fA2 2 3
The template specifies the following required inputs:
Variance of the random intercept
Covariance between the random intercept and slop
Variance of the random slop (fA2
)
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.
Once a design object has been created, calculating statistical power is straightforward.
The statistical power for F-tests (ANOVA-like tests)
can be computed using the pwr.anova
function.
Usage
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
).
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.
In this example, we create a CRD with one treatment factor.
Design parameters
Step 1: Creating the CRD
Step 2: Power for omnibus test
To compute the power for the overall F-test (ANOVA):
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
In this example, we create a Randomized Complete Block Design (RCBD) with two treatment factors in a 2×2 factorial design.
Design Parameters
A
and B
).beta
values are as follows
(Optional):A1B1
): 35A2
alone: 5 unitsB2
alone: 3 unitsA2:B2
): -2 units (i.e., the combined
effect of A2
and B2
is 2 units below the
additive effects).Step 1: Creating the RCBD
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
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
:
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.
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:
Treatments:
Main plot factor with 2 levels
Subplot factor with 3 levels
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.
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} \]
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.
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.
pwr4exp
is developed based on LMM theory. The general
form of an LMM can be expressed as:
\[ y = X\beta + Zu + \varepsilon \]
where: \(y\) represents the observations of the response variable, \(\beta\) represents the fixed effect coefficients, \(u\) denotes the random effects, where \(u \sim N_q(0, G)\), \(\varepsilon\) represents the random errors, where \(\varepsilon \sim N_n(0, R)\), \(X_{(n \times p)}\) and \(Z_{(n \times q)}\) are the design matrices for the fixed and random effects, respectively.
It is assumed that \(u\) and \(\varepsilon\) are independent, and the marginal distribution of \(y\) follows a normal distribution \(y \sim N_n(X\beta, V)\), where:
\[ V = ZGZ^T + R \]
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).
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:
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:
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.
Fai, A. H., & Cornelius, P. L. (1996). Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments. Journal of Statistical Computation and Simulation, 54(4), 363–378.
Giesbrecht, F. G., & Burns, J. C. (1985). Two-stage analysis based on a mixed model: Large-sample asymptotic theory and small-sample simulation results. Biometrics, 41(2), 477–486.
Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2(6), 110–114.
Stroup, W. W. (2002). Power analysis based on spatial effects mixed models: A tool for comparing design and analysis strategies in the presence of spatial variability. Journal of Agricultural, Biological, and Environmental Statistics, 7(4), 491–511.
Stroup, W. W. (2012). Generalized linear mixed models: Modern concepts, methods, and applications (1st ed.). CRC Press.