skpr (version 0.57.0)

eval_design: Calculate Power of an Experimental Design

Description

Evaluates the power of an experimental design, for normal response variables, given the design's run matrix and the statistical model to be fit to the data. Returns a data frame of parameter and effect powers. Designs can consist of both continuous and categorical factors. By default, eval_design assumes a signal-to-noise ratio of 2, but this can be changed with the effectsize or anticoef parameters.

Usage

eval_design(design, model, alpha, blocking = FALSE, anticoef = NULL,
  effectsize = 2, varianceratios = 1, contrasts = contr.sum,
  conservative = FALSE, detailedoutput = FALSE, ...)

Arguments

design

The experimental design. Internally, eval_design rescales each numeric column to the range [-1, 1], so you do not need to do this scaling manually.

model

The model used in evaluating the design. It can be a subset of the model used to generate the design, or include higher order effects not in the original design generation. It cannot include factors that are not present in the experimental design.

alpha

The specified type-I error.

blocking

Default FALSE. If TRUE, eval_design will look at the rownames to determine blocking structure.

anticoef

The anticipated coefficients for calculating the power. If missing, coefficients will be automatically generated based on the effectsize argument.

effectsize

The signal-to-noise ratio. Default 2. For continuous factors, this specifies the difference in response between the highest and lowest levels of the factor (which are -1 and +1 after eval_design normalizes the input data), assuming that the root mean square error is 1. If you do not specify anticoef, the anticipated coefficients will be half of effectsize. If you do specify anticoef, effectsize will be ignored.

varianceratios

Default 1. The ratio of the whole plot variance to the run-to-run variance. For designs with more than one subplot this ratio can be a vector specifying the variance ratio for each subplot. Otherwise, it will use a single value for all strata.

contrasts

Default contr.sum. The function to use to encode the categorical factors in the model matrix. If the user has specified their own contrasts for a categorical factor using the contrasts function, those will be used. Otherwise, skpr will use contr.sum.

conservative

Specifies whether default method for generating anticipated coefficents should be conservative or not. TRUE will give the most conservative estimate of power by setting all but one level in each categorical factor's anticipated coefficients to zero. Default FALSE.

detailedoutput

If TRUE, return additional information about evaluation in results. Default FALSE.

...

Additional arguments.

Value

A data frame with the parameters of the model, the type of power analysis, and the power. Several design diagnostics are stored as attributes of the data frame. In particular, the modelmatrix attribute contains the model matrix that was used for power evaluation. This is especially useful if you want to specify the anticipated coefficients to use for power evaluation. The model matrix provides the order of the model coefficients, as well as the encoding used for categorical factors.

Details

This function evaluates the power of experimental designs.

Power is calculated under a linear regression framework: you intend to fit a linear model to the data, of the form

\(y = X \beta + \epsilon\), (plus blocking terms, if applicable)

where \(y\) is the vector of experimental responses, \(X\) is the model matrix, \(\beta\) is the vector of model coefficients, and \(\epsilon\) is the statistical noise. eval_design assumes that \(\epsilon\) is normally distributed with zero mean and unit variance (root-mean-square error is 1), and calculates both parameter power and effect power. Parameter power is the probability of rejecting the hypothesis \(\beta_i = 0\), where \(\beta_i\) is a single parameter in the model, while effect power is the probability of rejecting the hypothesis \(\beta_{i} = 0\), where \(\beta_{i}\) is the set of all parameters associated with the effect in question. The two powers are equivalent for continuous factors and two-level categorical factors, but they can be different for categorical factors with three or more levels.

When using conservative = TRUE, eval_design first evaluates the power with default coefficients. Then, for each multi-level categorical, it sets all coefficients to zero except the level that produced the lowest power, and then re-evaluates the power with this modified set of anticipated coefficients.

Examples

# NOT RUN {
#Generating a simple 2x3 factorial to feed into our optimal design generation
#of an 11-run design.
factorial = expand.grid(A = c(1, -1), B = c(1, -1), C = c(1, -1))

optdesign = gen_design(candidateset = factorial,
                      model= ~A + B + C, trials = 11, optimality = "D", repeats = 100)

#Now evaluating that design (with default anticipated coefficients and a effectsize of 2):
eval_design(design = optdesign, model= ~A + B + C, alpha = 0.2)

#Evaluating a subset of the design (which changes the power due to a different number of
#degrees of freedom)
eval_design(design = optdesign, model= ~A + C, alpha = 0.2)

#Halving the signal-to-noise ratio by setting a different effectsize (default is 2):
eval_design(design = optdesign, model= ~A + B + C, alpha = 0.2, effectsize = 1)

#With 3+ level categorical factors, the choice of anticipated coefficients directly changes the
#final power calculation. For the most conservative power calculation, that involves
#setting all anticipated coefficients in a factor to zero except for one. We can specify this
#option with the "conservative" argument.

factorialcoffee = expand.grid(cost = c(1, 2),
                              type = as.factor(c("Kona", "Colombian", "Ethiopian", "Sumatra")),
                              size = as.factor(c("Short", "Grande", "Venti")))

designcoffee = gen_design(factorialcoffee,
                         ~cost + size + type, trials = 29, optimality = "D", repeats = 100)

#Evaluate the design, with default anticipated coefficients (conservative is FALSE by default).
#(Setting detailedoutput = TRUE provides information on the anticipated
#coefficients that were used:)
eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05, detailedoutput = TRUE)

#Evaluate the design, with conservative anticipated coefficients:
eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05, detailedoutput = TRUE,
            conservative = TRUE)

#which is the same as the following, but now explicitly entering the coefficients:
eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05,
            anticoef = c(1, 1, 1, 0, 0, 1, 0), detailedoutput = TRUE)


#If the defaults do not suit you, enter the anticipated coefficients in manually.
eval_design(designcoffee,
           model = ~cost + size + type, alpha = 0.05, anticoef = c(1, 1, 0, 0, 1, 0, 1))

#You can also evaluate the design with higher order effects, even if they were not
#used in design generation:
eval_design(designcoffee, model = ~cost + size + type + cost * type, alpha = 0.05)

#Split plot designs can also be evaluated by setting the blocking parameter as TRUE.

#Generating split plot design
splitfactorialcoffee = expand.grid(caffeine = c(1, -1),
                                  cost = c(1, 2),
                                  type = as.factor(c("Kona", "Colombian", "Ethiopian", "Sumatra")),
                                  size = as.factor(c("Short", "Grande", "Venti")))

coffeeblockdesign = gen_design(splitfactorialcoffee, ~caffeine, trials = 12)
coffeefinaldesign = gen_design(splitfactorialcoffee,
                              model = ~caffeine + cost + size + type, trials = 36,
                              splitplotdesign = coffeeblockdesign, splitplotsizes = 3)

#Evaluating design
eval_design(coffeefinaldesign, ~cost + size + type + caffeine, 0.2, blocking = TRUE)

#We can also evaluate the design with a custom ratio between the whole plot error to
#the run-to-run error.
eval_design(coffeefinaldesign, ~caffeine + cost + size + type + caffeine, 0.2, blocking = TRUE,
            varianceratios = 2)

#If the design was generated outside of skpr and thus the row names do not have the
#blocking structure encoded already, the user can add these manually. For a 12-run
#design with 4 blocks, the rownames will be as follows:

manualrownames = paste(c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4), rep(c(1, 2, 3), 4), sep = ".")

#If we wanted to add this blocking structure to the design coffeeblockdesign, we would
#simply set the rownames to this character vector:

rownames(coffeeblockdesign) = manualrownames

#Deeper levels of blocking can be specified with additional periods.
# }