R/eval_design.R

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

#'@title 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, \code{eval_design}
#'assumes a signal-to-noise ratio of 2, but this can be changed with the
#'\code{effectsize} or \code{anticoef} parameters.
#'
#'@param design The experimental design. Internally, \code{eval_design} rescales each numeric column
#'to the range [-1, 1], so you do not need to do this scaling manually.
#'@param 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.
#'@param alpha Default `0.05`. The specified type-I error.
#'@param blocking Default `NULL`. If `TRUE`, \code{eval_design} will look at the rownames (or blocking columns) to determine blocking structure. Default FALSE.
#'@param anticoef The anticipated coefficients for calculating the power. If missing, coefficients
#'will be automatically generated based on the \code{effectsize} argument.
#'@param 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 \code{eval_design}
#' normalizes the input data), assuming that the root mean square error is 1. If you do not specify \code{anticoef},
#' the anticipated coefficients will be half of \code{effectsize}. If you do specify \code{anticoef}, \code{effectsize} will be ignored.
#'@param 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.
#'@param contrasts Default \code{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`.
#'@param detailedoutput If `TRUE``, return additional information about evaluation in results. Default FALSE.
#'@param 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`.
#'@param 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`.
#'@param 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.
#'@param high_resolution_candidate_set Default `NA`. If you have continuous numeric terms and disallowed combinations, the closed-form I-optimality value
#' cannot be calculated and must be approximated by numeric integration. This requires sampling the allowed space densely, but most candidate sets will provide
#' a sparse sampling of allowable points. To work around this, skpr will generate a convex hull of the numeric terms for each unique combination of categorical
#' factors to generate a dense sampling of the space and cache that value internally, but this is a slow calculation and does not support non-convex candidate sets.
#' To speed up moment matrix calculation,  pass a higher resolution version of your candidate set here with the disallowed combinations already applied.
#' If you generated your design externally from skpr, there are disallowed combinations in your design, and need correct I-optimalituy values, you must pass your candidate set here.
#'@param moment_sample_density Default `20`. The density of points to sample when calculating the moment matrix to compute I-optimality. Only
#'required if the design was generated outside of skpr and there are disallowed combinations. It is much faster and potentially more accurate to
#' provide your own candidate set with higher resolution continuous factors to `high_resolution_candidate_set`.`
#'@param ... Additional arguments.
#'@return 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 \code{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.
#'
#'If the design is has no blocking or restrictions on randomization, the model assumed is:
#'
#'\eqn{y = X \beta + \epsilon}.
#'
#'If the design is a split-plot design, the model is as follows:
#'
#'\ifelse{html}{\eqn{y = X \beta + Z b}\out{<sub>i</sub>} + \eqn{\epsilon}\out{<sub>ij</sub>}}{\eqn{y = X \beta + Z b_{i} + \epsilon_{ij}}},
#'
#'Here, \eqn{y} is the vector of experimental responses, \eqn{X} is the model matrix, \eqn{\beta} is
#'the vector of model coefficients, \eqn{Z_{i}} are the blocking indicator,
#'\eqn{b_{i}} is the random variable associated with the \eqn{i}th block, and \eqn{\epsilon}
#'is a random variable normally distributed with zero mean and unit variance (root-mean-square error is 1.0).
#'
#'\code{eval_design} calculates both parameter power as well as effect power, defined as follows:
#'
#'1) Parameter power is the probability of rejecting the hypothesis \eqn{H_0 : \beta_i = 0}, where \eqn{\beta_i} is a single parameter
#'in the model
#'2) Effect power is the probability of rejecting the hypothesis \eqn{H_0 : \beta_{1} = \beta_{2} = ... = \beta_{n} 0} for all \eqn{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 \code{conservative = TRUE}, \code{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.
#'@export
#'@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)
eval_design = function(
    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,
    high_resolution_candidate_set = NA,
    moment_sample_density = 20,
    ...
) {
    if (missing(design)) {
        stop("skpr: No design detected in arguments.")
    }
    if (missing(model) || (is.numeric(model) && missing(alpha))) {
        if (is.numeric(model) && missing(alpha)) {
            alpha = model
        }
        if (is.null(attr(design, "generating.model"))) {
            stop(
                "skpr: No model detected in arguments or in design attributes."
            )
        } else {
            model = attr(design, "generating.model")
        }
    }
    user_specified_varianceratio = TRUE
    if (is.null(varianceratios)) {
        user_specified_varianceratio = FALSE
        if (!is.null(attr(design, "varianceratios"))) {
            varianceratios = attr(design, "varianceratios")
        } else {
            varianceratios = 1
        }
    }
    if (!is.null(attr(design, "splitcolumns"))) {
        if (
            varianceratios[length(varianceratios)] != 1 &&
                !user_specified_varianceratio
        ) {
            warning(
                "Lowest level of varianceratios cannot be set to anything other than 1 (value of ",
                varianceratios[length(varianceratios)],
                " was set during design generation). Setting run-to-run variance to 1."
            )
            varianceratios[length(varianceratios)] = 1
        }
    }
    input_design = design
    candidate_set = attr(input_design, "candidate_set")

    if(!is.null(candidate_set)) {
      normalized_candidate_set = normalize_design(candidate_set)
    }
    args = list(...)
    if ("RunMatrix" %in% names(args)) {
        stop("skpr: RunMatrix argument deprecated. Use `design` instead.")
    }

    if (is.null(blocking)) {
        blocking = FALSE
        if (!is.null(attr(design, "blocking"))) {
            blocking = attr(design, "blocking")
        }
        if (!is.null(attr(design, "splitplot"))) {
            blocking = blocking || attr(design, "splitplot")
        }
    }

    #detect pre-set contrasts
    presetcontrasts = list()
    for (x in names(design)[
        lapply(design, class) %in% c("character", "factor")
    ]) {
        if (!is.null(attr(design[[x]], "contrasts"))) {
            presetcontrasts[[x]] = attr(design[[x]], "contrasts")
        }
    }
    # reorder levels for the conservative calculation (if not a balanced design)
    if (conservative) {
        for (x in names(design)[
            lapply(design, class) %in% c("character", "factor")
        ]) {
            number_levels = table(design[[x]])
            if (length(unique(number_levels)) != 1) {
                if (identical(contrasts, contr.sum)) {
                    number_levels = table(design[[x]])
                    order_levels = names(number_levels)[order(number_levels)]
                    design[[x]] = factor(
                        design[[x]],
                        levels = rev(order_levels)
                    )
                } else if (identical(contrasts, contr.treatment)) {
                    number_levels = table(design[[x]])
                    order_levels = names(number_levels)[order(number_levels)]
                    design[[x]] = factor(design[[x]], levels = order_levels)
                }
            }
        }
    }
    if (is.null(advancedoptions$aliaspower)) {
        aliaspower = 2
    } else {
        if (!is.numeric(advancedoptions$aliaspower)) {
            stop("skpr: advancedoptions$aliaspower must be a positive integer")
        }
        aliaspower = advancedoptions$aliaspower
    }
    nointercept = attr(
        stats::terms.formula(model, data = design),
        "intercept"
    ) ==
        0
    #covert tibbles
    run_matrix_processed = as.data.frame(design)
    #Detect externally generated blocking columns and convert to rownames
    run_matrix_processed = convert_blockcolumn_rownames(
        run_matrix_processed,
        blocking,
        varianceratios
    )
    varianceratios = attr(run_matrix_processed, "tempvar")
    attr(run_matrix_processed, "tempvar") = NULL
    zlist = attr(run_matrix_processed, "z.matrix.list")

    #Remove skpr-generated REML blocking columns if present
    run_matrix_processed = remove_skpr_blockcols(run_matrix_processed)

    #----- Convert dots in formula to terms -----#
    model = convert_model_dots(run_matrix_processed, model)

    #----- Rearrange formula terms by order -----#
    model = rearrange_formula_by_order(model, data = run_matrix_processed)
    if (nointercept) {
        model = update.formula(model, ~ -1 + .)
    }

    #---- Reduce run matrix to terms in model ---#
    run_matrix_processed = reduceRunMatrix(run_matrix_processed, model)

    #---Develop contrast lists for model matrix---#
    #Variables used later: contrastslist, contrastslist_cormat
    contrastslist = list()
    contrastslist_cormat = list()
    for (x in names(run_matrix_processed[
        lapply(run_matrix_processed, class) %in% c("character", "factor")
    ])) {
        if (!(x %in% names(presetcontrasts))) {
            contrastslist[[x]] = contrasts
        } else {
            contrastslist[[x]] = presetcontrasts[[x]]
        }
        contrastslist_cormat[[x]] = contr.simplex
    }
    if (length(contrastslist) < 1) {
        contrastslist = NULL
        contrastslist_cormat = NULL
    }

    #------Normalize/Center numeric columns ------#
    run_matrix_processed = normalize_design(run_matrix_processed)

    #-Generate Model Matrix & Anticipated Coefficients-#
    #Variables used later: anticoef
    attr(run_matrix_processed, "model.matrix") = model.matrix(
        model,
        run_matrix_processed,
        contrasts.arg = contrastslist
    )

    if (!missing(anticoef) && !missing(effectsize)) {
        warning(
            "User defined anticipated coefficients (anticoef) detected; ignoring effectsize argument."
        )
    }
    if (missing(anticoef)) {
        default_coef = gen_anticoef(run_matrix_processed, model, nointercept)
        anticoef = anticoef_from_delta(default_coef, effectsize, "gaussian")
        if (nointercept) {
            anticoef = anticoef[-1]
        }
    }
    if (length(anticoef) != dim(attr(run_matrix_processed, "model.matrix"))[2]) {
        stop("skpr: Wrong number of anticipated coefficients")
    }

    #-----Generate V inverse matrix-----X
    #Variables used later: V, vinv, degrees_of_freedom, parameter_names
    if (blocking) {
        V = convert_rownames_to_covariance(
            run_matrix_processed,
            varianceratios,
            user_specified_varianceratio
        )
        # varianceratios = attr(V,"tempvar")
        attr(run_matrix_processed, "tempvar") = NULL
        vinv = solve(V)
        #Below code detects the split-plot columns, and calculates the adjusted degrees of freedom for each term
        degrees_of_freedom = calculate_degrees_of_freedom(
            run_matrix_processed,
            nointercept,
            model,
            contrasts
        )
        if (!nointercept) {
            extract_intnames_formula = ~1
            parameter_names = list()
            parameter_names[["(Intercept)"]] = "(Intercept)"
        } else {
            extract_intnames_formula = ~ -1
            parameter_names = list()
        }
        for (i in (length(parameter_names) + 1):length(degrees_of_freedom)) {
            currentterm = names(degrees_of_freedom)[i]
            extract_intnames_formula = update.formula(
                extract_intnames_formula,
                as.formula(paste0("~. + ", currentterm))
            )
            newcolnames = suppressWarnings(colnames(model.matrix(
                extract_intnames_formula,
                data = design,
                contrasts.arg = contrastslist
            )))
            parameter_names[[currentterm]] = newcolnames[
                !(newcolnames %in% unlist(parameter_names))
            ]
        }
    } else {
        vinv = NULL
        degrees_of_freedom = NULL
        parameter_names = NULL
    }
    factornames = attr(terms(model), "term.labels")
    levelvector = calculate_level_vector(
        run_matrix_processed,
        model,
        nointercept
    )
    effectresults = effectpower(
        run_matrix_processed,
        levelvector,
        anticoef,
        alpha,
        vinv = vinv,
        degrees = degrees_of_freedom
    )
    parameterresults = parameterpower(
        run_matrix_processed,
        levelvector,
        anticoef,
        alpha,
        vinv = vinv,
        degrees = degrees_of_freedom,
        parameter_names = parameter_names
    )

    typevector = c(
        rep("effect.power", length(effectresults)),
        rep("parameter.power", length(parameterresults))
    )
    if (!nointercept) {
        effectnamevector = c("(Intercept)", factornames)
    } else {
        effectnamevector = factornames
    }
    parameternamevector = colnames(attr(run_matrix_processed, "model.matrix"))
    namevector = c(effectnamevector, parameternamevector)
    powervector = c(effectresults, parameterresults)

    results = data.frame(
        parameter = namevector,
        type = typevector,
        power = powervector
    )

    if (length(namevector) != length(typevector)) {
        warning("Number of names does not equal number of power calculations")
    }

    attr(results, "model.matrix") = attr(run_matrix_processed, "model.matrix")
    attr(results, "anticoef") = anticoef

    modelmatrix_cor = model.matrix(
        model,
        run_matrix_processed,
        contrasts.arg = contrastslist_cormat
    )
    if (ncol(modelmatrix_cor) > 2) {
        if (!blocking) {
            V = diag(nrow(modelmatrix_cor))
        }
        if (!nointercept) {
            correlation.matrix = abs(cov2cor(solve(
                t(modelmatrix_cor) %*% solve(V) %*% modelmatrix_cor
            ))[-1, -1])
            colnames(correlation.matrix) = colnames(modelmatrix_cor)[-1]
            rownames(correlation.matrix) = colnames(modelmatrix_cor)[-1]
        } else {
            correlation.matrix = abs(cov2cor(solve(
                t(modelmatrix_cor) %*% solve(V) %*% modelmatrix_cor
            )))
            colnames(correlation.matrix) = colnames(modelmatrix_cor)
            rownames(correlation.matrix) = colnames(modelmatrix_cor)
        }
        attr(results, "correlation.matrix") = round(correlation.matrix, 8)
        tryCatch(
            {
                if (ncol(attr(run_matrix_processed, "model.matrix")) > 2) {
                    amodel = aliasmodel(model, aliaspower)
                    if (amodel != model) {
                        aliasmatrix = suppressWarnings({
                            model.matrix(
                                aliasmodel(model, aliaspower),
                                design,
                                contrasts.arg = contrastslist_cormat
                            )[, -1]
                        })
                        A = solve(t(modelmatrix_cor) %*% modelmatrix_cor) %*%
                            t(modelmatrix_cor) %*%
                            aliasmatrix
                        attr(results, "alias.matrix") = A
                        attr(results, "trA") = sum(diag(t(A) %*% A))
                    } else {
                        attr(
                            results,
                            "alias.matrix"
                        ) = "No alias matrix calculated: full model specified"
                        attr(
                            results,
                            "trA"
                        ) = "No alias trace calculated: full model specified"
                    }
                }
            },
            error = function(e) {
            }
        )
    }
    attr(results, "generating.model") = model
    attr(results, "runmatrix") = run_matrix_processed
    attr(results, "model.matrix") = modelmatrix_cor
    attr(results, "blocking") = blocking
    attr(results, "varianceratios") = varianceratios
    attr(results, "alpha") = alpha
    attr(results, "contrastslist") = contrastslist

    factors = colnames(modelmatrix_cor)
    classvector = sapply(lapply(run_matrix_processed, unique), class) ==
        "factor"

    moment_matrix = get_moment_matrix(
      design = input_design,
      candidate_set_normalized = normalized_candidate_set,
      factors = factors,
      classvector = classvector,
      model = model,
      moment_sample_density = moment_sample_density,
      high_resolution_candidate_set = high_resolution_candidate_set
    )
    if(all(!is.na(moment_matrix))) {
      attr(results, "moments.matrix") = moment_matrix
    }
    attr(results, "A") = AOptimality(modelmatrix_cor)

    if (!blocking) {
        attr(results, "variance.matrix") = diag(nrow(modelmatrix_cor)) *
            varianceratios
        if(all(!is.na(moment_matrix))) {
          attr(results, "I") = IOptimality(
              modelmatrix_cor,
              momentsMatrix = moment_matrix,
              blockedVar = diag(nrow(modelmatrix_cor))
          )
        }
        attr(results, "D") = 100 * DOptimalityLog(modelmatrix_cor)
        attr(results, "T") = sum(diag(t(modelmatrix_cor) %*% modelmatrix_cor))
        attr(results, "E") = min(unlist(eigen(
            t(modelmatrix_cor) %*% modelmatrix_cor
        )["values"]))
    } else {
        attr(results, "z.matrix.list") = zlist
        attr(results, "variance.matrix") = V
        if(all(!is.na(moment_matrix))) {
          attr(results, "I") = IOptimality(
              modelmatrix_cor,
              momentsMatrix = moment_matrix, #Fix here maybe, need blocked moment matrix?
              blockedVar = V
          )
        }
        deffic = DOptimalityBlocked(modelmatrix_cor, blockedVar = V)
        if (!is.infinite(deffic)) {
            attr(results, "D") = 100 *
                DOptimalityBlocked(modelmatrix_cor, blockedVar = V)^(1 /
                    ncol(modelmatrix_cor)) /
                nrow(modelmatrix_cor)
        } else {
            attr(results, "D") = 100 *
                DOptimalityBlockedLog(modelmatrix_cor, blockedVar = V)^(1 /
                    ncol(modelmatrix_cor)) /
                nrow(modelmatrix_cor)
        }
    }
    if (detailedoutput) {
        if (nrow(results) != length(anticoef)) {
            results$anticoef = c(
                rep(NA, nrow(results) - length(anticoef)),
                anticoef
            )
        } else {
            results$anticoef = anticoef
        }
        results$alpha = alpha
        results$trials = nrow(run_matrix_processed)
    }

    #For conservative coefficients, look for lowest power results from non-conservative calculation and set them to one
    #and the rest to zero. (If equally low results, apply 1 -1 pattern to lowest)

    if (conservative) {
        #at this point, since we are going to specify anticoef, do not use the effectsize argument
        #in the subsequent call. Do replicate the magnitudes from the original anticoef
        conservative_anticoef = calc_conservative_anticoef(results, effectsize)
        results = eval_design(
            design = design,
            model = model,
            alpha = alpha,
            blocking = blocking,
            anticoef = conservative_anticoef,
            detailedoutput = detailedoutput,
            varianceratios = varianceratios,
            contrasts = contrasts,
            conservative = FALSE,
            reorder_factors = reorder_factors
        )
    }

    if (reorder_factors) {
        if (detailedoutput) {
            results$reordered_factors = NA
        }
        results_temp = results
        sum_effect_power = sum(
            with(results_temp, power[type == "effect.power"]),
            na.rm = TRUE
        )
        for (i in seq_len(ncol(design))) {
            if (is.factor(design[, i])) {
                number_factors = length(levels(design[, i]))
                temp_design = design
                cycle = list()
                cycle[[1]] = levels(temp_design[, i])
                order_cycle = seq_len(number_factors)
                for (j in seq_len(number_factors)[-1]) {
                    order_cycle = c(order_cycle[-1], order_cycle[1])
                    cycle[[j]] = levels(temp_design[, i])[order_cycle]
                }
                for (j in seq_len(number_factors)) {
                    temp_design[, i] = factor(temp_design[, i], cycle[[j]])
                    adjusted_sum_power = eval_design(
                        design = temp_design,
                        model = model,
                        alpha = alpha,
                        blocking = blocking,
                        detailedoutput = detailedoutput,
                        varianceratios = varianceratios,
                        contrasts = contrasts,
                        conservative = conservative,
                        reorder_factors = FALSE
                    )
                    new_sum_power = sum(with(
                        adjusted_sum_power,
                        power[type == "effect.power"]
                    ))
                    if (new_sum_power < sum_effect_power) {
                        design[, i] = factor(design[, i], cycle[[j]])
                        results_temp = adjusted_sum_power
                        sum_effect_power = new_sum_power
                        if (detailedoutput) {
                            results_temp$reordered_factors[
                                results_temp$parameter == names(temp_design)[i]
                            ] = paste0(cycle[[j]], collapse = " ")
                        }
                    }
                }
            }
        }
        results = results_temp
    }
    if (!inherits(results, "skpr_eval_output")) {
        class(results) = c("skpr_eval_output", class(results))
    }
    if (any(is.na(results$power))) {
        warning(
            "skpr: NA indicates not enough degrees of freedom to estimate power for those terms."
        )
    }
    recommend_analysis_method = function(blocking) {
    }
    #Add recommended analysis method
    contrast_string = deparse(substitute(contrasts))
    attr(results, "contrast_string") = sprintf("`%s`", contrast_string)
    if (!blocking) {
        attr(results, "parameter_analysis_method_string") = "`lm(...)`"
        attr(
            results,
            "effect_analysis_method_string"
        ) = r"{`car::Anova(fit, type = "III")`}"
    } else {
        attr(
            results,
            "parameter_analysis_method_string"
        ) = "`lmerTest::lmer(...)`"
        attr(
            results,
            "effect_analysis_method_string"
        ) = r"{`car::Anova(fit, type = "III")`}"
    }
    return(results)
}