eval_design: Calculate Power of an Experimental Design

In skpr: Design of Experiments Suite: Generate and Evaluate Optimal Designs

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 = NULL,
  alpha = 0.05,
  blocking = NULL,
  anticoef = NULL,
  effectsize = 2,
  varianceratios = NULL,
  contrasts = contr.sum,
  conservative = FALSE,
  reorder_factors = FALSE,
  detailedoutput = FALSE,
  advancedoptions = NULL,
  ...
)

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. If this is missing and the design was generated with skpr, the generating model will be used. 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	Default '0.05'. The specified type-I error.
blocking	Default 'NULL'. If 'TRUE', eval_design will look at the rownames (or blocking columns) to determine blocking structure. Default FALSE.
anticoef	The anticipated coefficients for calculating the power. If missing, coefficients will be automatically generated based on the effectsize argument.
effectsize	Default '2'. The signal-to-noise ratio. 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 'NULL'. The ratio of the whole plot variance to the run-to-run variance. If not specified during design generation, this will default to 1. For designs with more than one subplot this ratio can be a vector specifying the variance ratio for each subplot (comparing to the run-to-run variance). 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 (or multiple if they are equally low) level in each categorical factor's anticipated coefficients to zero. Default 'FALSE'.
reorder_factors	Default 'FALSE'. If 'TRUE', the levels will be reordered to generate the most conservative calculation of effect power. The function searches through all possible reference levels for a given factor and chooses the one that results in the lowest effect power. The reordering will be presenting in the output when 'detailedoutput = TRUE'.
detailedoutput	If 'TRUE“, return additional information about evaluation in results. Default FALSE.
advancedoptions	Default 'NULL'. A named list with parameters to specify additional attributes to calculate. Options: 'aliaspower' gives the degree at which the Alias matrix should be calculated.
...	Additional arguments.

Details

This function evaluates the power of experimental designs.

If the design is has no blocking or restrictions on randomization, the model assumed is:

y = X \beta + \epsilon.

If the design is a split-plot design, the model is as follows:

y = X \beta + Z b_i + \epsilon_ij,

Here, y is the vector of experimental responses, X is the model matrix, \beta is the vector of model coefficients, Z_{i} are the blocking indicator, b_{i} is the random variable associated with the ith block, and \epsilon is a random variable normally distributed with zero mean and unit variance (root-mean-square error is 1.0).

eval_design calculates both parameter power as well as effect power, defined as follows:

1) Parameter power is the probability of rejecting the hypothesis H_0 : \beta_i = 0, where \beta_i is a single parameter in the model 2) Effect power is the probability of rejecting the hypothesis H_0 : \beta_{1} = \beta_{2} = ... = \beta_{n} 0 for all n coefficients for a categorical factor.

The two power types are equivalent for continuous factors and two-level categorical factors, but they will differ for categorical factors with three or more levels.

For split-plot designs, the degrees of freedom are allocated to each term according to the algorithm given in "Mixed-Effects Models in S and S-PLUS" (Pinheiro and Bates, pp. 91).

When using conservative = TRUE, eval_design first evaluates the power with the default (or given) 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. If there are two or more equal power values in a multi-level categorical, two of the lowest equal terms are given opposite sign anticipated coefficients and the rest (for that categorical factor) are set to zero.

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.

Examples

#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)

#We do not have to input the model if it's the same as the model used
#During design generation. Here, we also use the default value for alpha (`0.05`)
eval_design(optdesign)

#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).
eval_design(designcoffee)

#Evaluate the design, with conservative anticipated coefficients:
eval_design(designcoffee, conservative = TRUE)

#which is the same as the following, but now explicitly entering the coefficients:
eval_design(designcoffee, anticoef = c(1, 1, 1, 0, 0, 1, 0))

#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)


#Generating and evaluating a 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, blocksizes = 3)

#Evaluating design (blocking is automatically detected)
eval_design(coffeefinaldesign, 0.2, blocking = TRUE)

#Manually turn blocking off to see completely randomized design power
eval_design(coffeefinaldesign, 0.2, blocking = FALSE)

#We can also evaluate the design with a custom ratio between the whole plot error to
#the run-to-run error.
eval_design(coffeefinaldesign, 0.2, 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, here is a vector indicated the blocks:

blockcolumn = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4)

#If we wanted to add this blocking structure to the design coffeeblockdesign, we would
#add a column with the format "Block1", "Block2", "Block3" ... and each one will be treated
#as a separate blocking layer.

coffeeblockdesign$Block1 = blockcolumn

#By default, skpr will throw out the blocking columns unless the user specifies `blocking = TRUE`.
eval_design(coffeeblockdesign, blocking=TRUE)

skpr documentation built on July 9, 2023, 7:23 p.m.