Nothing
#' @title Sample Size and Effect Size Determination
#'
#' @description It searches by simulation
#' the sample size (given other factors,
#' such as effect sizes) or effect size
#' (given other factors, such as sample
#' size) with power to
#' detect an effect close to a target
#' value.
#'
#' @details
#' This is how to use [x_from_power()]:
#'
#' - Specify the model by [power4test()],
#' with `do_the_test = FALSE`, and set
#' the magnitude of the effect sizes
#' to the minimum levels to detect.
#'
#' - Add the test using [power4test()]
#' using `test_fun` and `test_args`
#' (see the help page of [power4test()]
#' for details). Run it on the
#' starting sample size or
#' effect size.
#'
#' - Call [x_from_power()] on the output
#' of [power4test()] returned from
#' the previous step. This
#' function will iteratively repeat
#' the analysis on either other sample
#' sizes, or other values for a
#' selected model parameter (the
#' effect sizes),
#' trying to achieve a goal (`goal`) for
#' a value of interest (`what`).
#'
#' If the `goal` is `"ci_hit"`, the
#' search will try to find a value (a sample
#' size, or a population value of
#' the selected model parameter) with
#' a power level close enough to the
#' target power, defined by having its
#' confidence interval for the power
#' including the target power.
#'
#' If the `goal` is `"close_enough"`,
#' then the search will try to find a
#' value of `x` with its level of
#' power (`"point"`), the upper bound
#' of the confidence interval for this
#' level of power (`"ub"`), or the
#' lower bound of the confidence interval
#' fro this level of power (`"lb"`)
#' "close enough" to the target level of
#' power, defined by having an absolute
#' difference less than the `tolerance`.
#'
#' If several values of `x` (sample
#' size or the population value of
#' a model parameter) have already been
#' examined by [power4test_by_n()] or
#' [power4test_by_es()], the output
#' of these two functions can also be
#' used as `object` by [x_from_power()].
#'
#' Usually, the default values of the
#' arguments should be sufficient.
#'
#' The results can be viewed using
#' [summary()], and the output has
#' a `plot` method ([plot.x_from_power()]) to
#' plot the relation between power and
#' values (of `x`) examined.
#'
#' A detailed illustration on how to
#' use this function for sample size can be found
#' from this page:
#'
#' <https://sfcheung.github.io/power4mome/articles/x_from_power_for_n.html>
#'
#' # Algorithms
#'
#' Two algorithms are currently available,
#' the simple (though inefficient)
#' bisection method, and a method that
#' makes use of the estimated crude power
#' curve.
#'
#' Unlike typical root-finding problems,
#' the prediction of the level of power
#' is stochastic. Moreover, the computational
#' cost is high when Monte Carlo or
#' bootstrap confidence intervals are
#' used to do a test because the estimation
#' of the power for one single value of
#' `x` can sometimes take one minute or
#' longer. Therefore, in addition to
#' the simple bisection method, a method,
#' named *power curve* method, was also
#' specifically developed for this
#' scenario.
#'
#' ## Bisection Method
#'
#' This method, `algorithm = "bisection"`,
#' basically starts with
#' an interval that probably encloses the
#' value of `x` that meets the goal,
#' and then successively narrows this
#' interval. The mid-point of this
#' interval is used as the estimate.
#' Though simple, there are cases in
#' which it can be slow. Nevertheless,
#' preliminary examination suggests that
#' this method is good enough for common
#' scenarios. Therefore, this method is
#' the default algorithm when `x` is
#' `n`.
#'
#' ## Power Curve Method
#'
#' This method, `algorithm = "power_curve"`,
#' starts with a crude
#' power curve based on a few points.
#' This tentative model is then used
#' to suggest the values to examine in
#' the next iteration. The form, not
#' just the parameters, of the
#' model can change across iterations,
#' as more and more data points are
#' available.
#'
#' This method can be used only with
#' the goal `"ci_hit"`.
#' This method is the default method
#' for `x = "es"` with `goal = "ci_hit"`
#' because the relation
#' between the power and the population
#' value of a parameter varies across
#' parameters, unlike the relation
#' between power and sample size. Therefore,
#' taking into account the working
#' power curve may help finding the
#' desired value of `x`.
#'
#' The technical internal workflow of
#' this method implemented in
#' [x_from_power()] can be found in
#' this page: <https://sfcheung.github.io/power4mome/articles/x_from_power_workflow.html>.
#'
#' @return
#' The function [x_from_power()]
#' returns an `x_from_power` object,
#' which is a list with the following
#' elements:
#'
#' - `power4test_trials`: The output of
#' [power4test_by_n()] for all sample
#' sizes examined, or of
#' [power4test_by_es()] for all
#' population values of the selected
#' parameter examined.
#'
#' - `rejection_rates`: The output of
#' [rejection_rates()].
#'
#' - `x_tried`: The sample sizes or
#' population values
#' examined.
#'
#' - `power_tried`: The estimated
#' rejection rates for all the values
#' examined.
#'
#' - `x_final`: The sample size or
#' population value in the
#' solution. `NA` if a solution not found.
#'
#' - `power_final`: The estimated power
#' of the value in the solution.
#' `NA` if a solution not found.
#'
#' - `i_final`: The position of the
#' solution in `power4test_trials`.
#' `NA` if a solution not found.
#'
#' - `ci_final`: The confidence interval
#' of the estimated power in the solution,
#' formed by normal approximation.
#' `NA` if a solution not found.
#'
#' - `ci_level`: The level of confidence
#' of `ci_final`.
#'
#' - `nrep_final`: The number of
#' replications (`nrep`) when estimating
#' the power in the solution.
#'
#' - `power_curve`: The output of
#' [power_curve()] when estimating the
#' power curve.
#'
#' - `target_power`: The requested
#' target power.
#'
#' - `power_tolerance`: The allowed
#' difference between the solution's
#' estimated power and the target
#' power. Determined by the number
#' of replications and the level of
#' confidence of the confidence intervals.
#'
#' - `x_estimated`: The value
#' (sample size or population value)
#' with the target power, estimated by
#' `power_curve`. This is used, when
#' solution not found, to determine the
#' range of the values to search when
#' calling the function again.
#'
#' - `start`: The time and date when
#' the process started.
#'
#' - `end`: The time and date when the
#' process ended.
#'
#' - `time_spent`: The time spent in
#' doing the search.
#'
#' - `args`: A named list of the arguments
#' of [x_from_power()] used in the search.
#'
#' - `call`: The call when this function
#' is called.
#'
#' @param object A `power4test` object,
#' which is the output of [power4test()].
#' Can also be a `power4test_by_n` object,
#' the output
#' of [power4test_by_n()], or
#' a `power4test_by_es` object, the
#' output of
#' [power4test_by_es()]. For these
#' two types of objects, the attempt
#' with power closest to the
#' `target_power` will be used as
#' `object`, and all other attempts in
#' them will be included in the estimation
#' of subsequent attempts and the final
#' output. Last, it can also be the
#' output of a previous call to
#' [x_from_power()], and the stored
#' trials will be retrieved.
#'
#' @param x For [x_from_power()],
#' `x` set the value to
#' be searched. Can be `"n"`, the sample
#' size, or `"es"`, the population value
#' of a parameter (set by `pop_es_name`).
#' For the `print` method of `x_from_power`
#' objects, this is the output of
#' [x_from_power()].
#'
#' @param pop_es_name The name of the
#' parameter. Required if `x` is `"es"`.
#' See the help page
#' of [ptable_pop()] on the names for
#' the argument `pop_es`.
#'
#' @param target_power The target power,
#' a value greater than 0 and less than
#' one.
#'
#' @param what The value for which is
#' searched: the estimate power (`"point"`),
#' the upper bound of the confidence
#' interval (`"ub"`), or the lower bound
#' of the confidence interval (`"lb"`).
#'
#' @param goal The goal of the search.
#' If `"ci_hit"`, then the goal is to
#' find a value of `x` with the
#' confidence interval of the estimated
#' power including the target power.
#' If `"close_enough"`, then the goal
#' is to find a value of `x` with the
#' value in `what` "close enough" to
#' the target power, defined by having
#' an absolute difference with the
#' target power less than `tolerance`.
#'
#' @param tolerance Used when the goal
#' is `"close_enough"`.
#'
#' @param ci_level The level of confidence
#' of the confidence intervals computed
#' for the estimated power. Default is
#' .95, denoting 95%.
#'
#' @param x_interval A vector of
#' two values, the minimum value
#' and the maximum values of `x`, in
#' the search for the values
#' (sample sizes or population values).
#' If `NULL`, default when `x = "es"`,
#' it will be determined internally.
#'
#' @param extendInt Whether `x_interval`
#' can be expanded when estimating the
#' the values to try. The value will
#' be passed to the argument of the
#' same name in [stats::uniroot()].
#' If `x` is `"n"`, then the default
#' value is `"upX"`.
#' That is, a value higher than
#' the maximum in `x_interval` is
#' allowed, if predicted by the tentative
#' model. Otherwise, the default value
#' is `"no"`. See the help page of
#' [stats::uniroot()] for further
#' information.
#'
#' @param progress Logical. Whether
#' the searching progress is reported.
#'
#' @param simulation_progress Logical.
#' Whether the progress in each call
#' to [power4test()], [power4test_by_n()],
#' or [power4test_by_es()]
#' is shown. To be passed to
#' the `progress` argument of these
#' functions.
#'
#' @param max_trials The maximum number
#' of trials in searching the value
#' with the target power. Rounded
#' up if not an integer.
#'
#' @param final_nrep The number of
#' replications in the final stage,
#' also the maximum number of replications
#' in each call to [power4test()],
#' [power4test_by_n()], or
#' [power4test_by_es()].
#'
#' @param final_R The number of
#' Monte Carlo simulation or
#' bootstrapping samples in the final
#' stage. The `R` in calling
#' [power4test()], [power4test_by_n()],
#' or [power4test_by_es()]
#' will be stepped up to this value
#' when approaching the target
#' power. Do not need to be very large
#' because the goal is to estimate
#' power by replications, not for high
#' precision in one single replication.
#'
#' @param seed If not `NULL`, [set.seed()]
#' will be used to make the process
#' reproducible. This is not always
#' possible if many stages of
#' parallel processing is involved.
#'
#' @param x_include_interval Logical.
#' Whether the minimum and maximum
#' values in `x_interval` are mandatory
#' to be included in the values
#' to be searched.
#'
#' @param power_curve_args A named
#' list of arguments to be passed
#' [power_curve()] when estimating
#' the relation between power and `x`
#' (sample size or effect size). Please
#' refer to [power_curve()] on available
#' arguments. There is one except:
#' `power_model` is mapped to
#' the `formula` argument of
#' [power_curve()].
#'
#' @param save_sim_all If `FALSE`,
#' the default, the data in each
#' `power4test` object for each
#' value of `x` is not saved,
#' to reduce the size of the output.
#' If set to `TRUE`, the size of the
#' output can be very large in size.
#'
#' @param algorithm The algorithm for
#' finding `x`. Can be `"power_curve"`
#' or `"bisection"`. The default algorithm
#' depends on `x`.
#'
#' @param control A named list of
#' additional
#' arguments to be passed to the
#' algorithm to be used. For advanced
#' users.
#'
#' @param check_es_interval If `TRUE`,
#' the default, and `x` is `"es"`,
#' a conservative probable
#' range of valid values for the selected
#' parameter will be determined, and it
#' will be used instead of `x_interval`.
#' If the range spans both positive and
#' negative values, only the interval
#' of the same sign as the population
#' value in `object` will be used.
#'
#'
#' @seealso [power4test()], [power4test_by_n()],
#' and [power4test_by_es()].
#'
#' @examples
#'
#' # Specify the population model
#'
#' mod <-
#' "
#' m ~ x
#' y ~ m + x
#' "
#'
#' # Specify the population values
#'
#' mod_es <-
#' "
#' m ~ x: m
#' y ~ m: l
#' y ~ x: n
#' "
#'
#' # Generate the datasets
#'
#' sim_only <- power4test(nrep = 5,
#' model = mod,
#' pop_es = mod_es,
#' n = 100,
#' do_the_test = FALSE,
#' iseed = 2345)
#'
#' # Do a test
#'
#' test_out <- power4test(object = sim_only,
#' test_fun = test_parameters,
#' test_args = list(pars = "m~x"))
#'
#' # Determine the sample size with a power of .80 (default)
#'
#' # In real analysis, to have more stable results:
#' # - Use a larger final_nrep (e.g., 400).
#'
#' # If the default values are OK, this call is sufficient:
#' # power_vs_n <- x_from_power(test_out,
#' # x = "n",
#' # seed = 4567)
#' power_vs_n <- x_from_power(test_out,
#' x = "n",
#' progress = TRUE,
#' target_power = .80,
#' final_nrep = 5,
#' max_trials = 1,
#' seed = 1234)
#' summary(power_vs_n)
#' plot(power_vs_n)
#'
#'
#' @importFrom graphics abline arrows par points text title
#'
#' @export
x_from_power <- function(object,
x,
pop_es_name = NULL,
target_power = .80,
what = c("point", "ub", "lb"),
goal = switch(what,
point = "ci_hit",
ub = "close_enough",
lb = "close_enough"),
ci_level = .95,
tolerance = .02,
x_interval = switch(x,
n = c(50, 2000),
es = NULL),
extendInt = NULL,
progress = TRUE,
simulation_progress = TRUE,
max_trials = 10,
final_nrep = 400,
final_R = 1000,
seed = NULL,
x_include_interval = FALSE,
check_es_interval = TRUE,
power_curve_args = list(power_model = NULL,
start = NULL,
lower_bound = NULL,
upper_bound = NULL,
nls_control = list(),
nls_args = list()),
save_sim_all = FALSE,
algorithm = NULL,
control = list()
) {
# Inputs
# - Target power
# Steps
# - Get the initial sample size and power.
# - Find the k sample sizes to explore.
# - Fit the curve.
# - See if the target power is in this range.
# - NO: Expand.
# - YES: Narrow.
# - Stop when convergence achieved.
# Outputs
# - Initial power4test object.
# - Final power4test object.
# - Final model by nls.
what <- match.arg(what)
goal <- match.arg(goal,
c("ci_hit", "close_enough"))
changed_to_point <- FALSE
if ((goal == "ci_hit") &&
(what != "point")) {
changed_to_point <- TRUE
what <- "point"
}
x <- match.arg(x,
choices = c("n", "es"))
if (!is.null(algorithm)) {
algorithm <- match.arg(algorithm,
c("bisection",
"power_curve"))
}
# ==== Set default algorithm ====
# what and goal take precedence
if (goal == "ci_hit") {
if (is.null(algorithm)) {
algorithm <- switch(x,
n = "bisection",
es = "power_curve")
}
}
changed_to_bisection <- FALSE
if (goal == "close_enough") {
# Only bisection is supported
if (isTRUE((algorithm != "bisection")) &&
(!is.null(algorithm))) {
# warning("Only bisection is supported when goal is 'close_enough'. ",
# "Switched automatically to bisection.")
changed_to_bisection <- TRUE
}
algorithm <- "bisection"
}
# what: The value to be examined.
# goal:
# - ci_hit: Only relevant for what == "point"
# - close_enough: Can be used for all what.
a <- abs(stats::qnorm((1 - ci_level) / 2))
power_tolerance_in_interval <- a * sqrt(target_power * (1 - target_power) / final_nrep)
power_tolerance_in_final <- a * sqrt(target_power * (1 - target_power) / final_nrep)
# ==== Sanity checks ====
if (target_power <= 0 || target_power >= 1) {
stop("'target power' (",
target_power,
") not between 0 and 1.")
}
if (x == "n") {
if (min(x_interval) < 0) {
stop("The minimum value of 'n_interval' cannot be negative")
}
}
if ((length(x_interval) != 2) &&
!is.null(x_interval)) {
stop("'x_interval' must be a vector with exactly two values.")
}
if ((x_interval[2] < x_interval[1]) &&
!is.null(x_interval)) {
stop("'x_interval' must be of the form c(minimum, maximum).")
}
if (is.null(extendInt)) {
if (x == "n") {
extendInt <- match.arg(extendInt,
choices = c("yes", "no", "upX", "downX"))
}
if (x == "es") {
extendInt <- match.arg(extendInt,
choices = c("yes", "no", "upX", "downX"))
}
}
max_trials <- ceiling(max_trials)
if (max_trials < 1) {
stop("'max_trials' must be at least 1 (after rounding, if necessary).")
}
if (!is.null(seed)) {
set.seed(seed)
}
# ==== Process the object ====
# Is it a compatible x_from_power object?
# - If yes and solution not found,
# - extract the stored by_x object
# - If yes and solution found,
# - The object is returned as is.
is_x_from_power <- FALSE
if (inherits(object, "x_from_power")) {
# ==== x_from_power object ====
# Throw an error if incompatible
# TODO:
# - Will run if goal and/or what changed
check_x_from_power_as_input(object,
x = x,
pop_es_name = pop_es_name,
final_nrep = final_nrep,
ci_level = ci_level)
# Check these for compatibility:
is_x_from_power <- TRUE
if (object$solution_found &&
(object$nrep_final == final_nrep) &&
(object$ci_level == ci_level) &&
(object$target_power == target_power)) {
if (progress) {
cat("\n--- Solution Already Found ---\n\n")
cat("Solution already found in object. It is returned as is.\n")
}
return(object)
}
object <- object$power4test_trials
}
# Is it a by_x object??
# - If yes and solution not found,
# - extract the power4test object closest to the solution
# - If yes and solution found,
# - This is possible because the by_x object
# may be generated directly by user,
# not extracted from an x_from_power object.
is_by_x <- FALSE
if (inherits(object, "power4test_by_n") ||
inherits(object, "power4test_by_es")) {
# ==== by_x object ====
is_by_x <- TRUE
object_by_org <- object
# Whether a solution exists will be checked later
i_org <- find_solution(
object_by_org,
target_power = target_power,
ci_level = ci_level,
what = what,
tol = tolerance,
goal = goal,
final_nrep = 1,
closest_ok = TRUE,
if_ties = "min")
object <- object_by_org[[i_org]]
} else {
object_by_org <- NA
}
# The object to be used below is always a power4test object
time_start <- Sys.time()
if (progress) {
cat("\n--- Setting ---\n\n")
cat("Algorithm: ", algorithm, "\n")
cat("Goal: ", goal, "\n")
tmp <- switch(what,
point = "(Estimated Power)",
ub = "(Upper bound of the confidence interval)",
lb = "(Lower bound of the confidence interval)")
cat("What: ", what, " ", tmp, "\n")
if (changed_to_bisection) {
catwrap(paste0("Note: Only bisection is supported when goal is 'close_enough'. ",
"Algorithm switched automatically to bisection."),
exdent = 0)
}
if (changed_to_point) {
cat("\n")
catwrap(paste0("Note: The goal 'ci_hit' only supports 'point'; ",
"'what' is changed to 'point'."),
exdent = 0)
}
if (progress) {
cat("\n--- Progress ---\n\n")
catwrap("- Set 'progress = FALSE' to suppress displaying the progress.",
exdent = 2)
if (simulation_progress) {
catwrap(paste0("- Set 'simulation progress = FALSE' ",
"to suppress displaying the progress ",
"in the simulation."),
exdent = 2)
}
}
}
# Initialize these flags because a solution
# may already be present before doing the search
ci_hit <- FALSE
solution_found <- FALSE
status <- NULL
technical <- NULL
# === Check Existing Solution ===
# If the input object is a by_x object,
# - Check if the solution is already in one of the
# power4test object.
# - If yes, skip the search and create the
# x_from_power object.
# ==== Solution already in by_x? ====
if (is_by_x) {
i_org_hit <- find_solution(
object_by_org,
ci_level = ci_level,
target_power = target_power,
what = what,
tol = tolerance,
goal = goal,
final_nrep = final_nrep,
closest_ok = FALSE,
if_ties = "min")
if (!is.na(i_org_hit) && !is.null(i_org_hit)) {
# ==== Solution in by_x. Skip the search ====
# Solution already in the input.
# DO not do the search
if (progress) {
cat("\n--- Solution Already Found ---\n\n")
cat("Solution already found in the object. Search will be skipped.\n")
}
# ==== Create the output ====
# Prepare the objects as if the search has completed
ci_hit <- ifelse(goal == "ci_hit",
TRUE,
NA)
solution_found <- TRUE
i2 <- i_org_hit
if (is_x_from_power || is_by_x) {
by_x_1 <- object_by_org
tmp1 <- rejection_rates(object_by_org,
level = ci_level,
all_columns = TRUE)
fit_1 <- power_curve(by_x_1,
formula = power_curve_args$power_model,
start = power_curve_args$start,
lower_bound = power_curve_args$lower_bound,
upper_bound = power_curve_args$upper_bound,
nls_control = power_curve_args$nls_control,
nls_args = power_curve_args$nls_args,
verbose = progress)
ci_hit <- ci_hit
x_tried <- get_x_tried(tmp1,
x = x)
x_out <- switch(x,
n = tmp1$n[i2],
es = tmp1$es[i2])
power_out <- tmp1$reject[i2]
nrep_out <- tmp1$nrep[i2]
ci_out <- unlist(tmp1[i2, c("reject_ci_lo", "reject_ci_hi")])
by_x_out <- by_x_1[[i2]]
}
}
}
# x_interval is the actual range of values
# to be searched, not necessarily the one
# set by users.
# ==== Fix x_interval ====
if (solution_found) {
x_interval <- range(x_tried)
} else {
if (check_es_interval) {
# ==== Find probable es interval ====
x_interval <- fix_es_interval(object = object,
x = x,
pop_es_name = pop_es_name,
x_interval = x_interval,
progress = progress)
}
}
x_max <- max(x_interval)
x_min <- min(x_interval)
# ==== Call the algorithm ====
if ((algorithm == "power_curve") && !solution_found) {
a_out <- do.call(alg_power_curve,
c(list(
object = object,
x = x,
pop_es_name = pop_es_name,
target_power = target_power,
x_max = x_max,
x_min = x_min,
progress = progress,
x_include_interval = x_include_interval,
x_interval = x_interval,
simulation_progress = simulation_progress,
save_sim_all = save_sim_all,
is_by_x = is_by_x,
object_by_org = object_by_org,
power_model = power_curve_args$power_model,
start = power_curve_args$start,
lower_bound = power_curve_args$lower_bound,
upper_bound = power_curve_args$upper_bound,
nls_control = power_curve_args$nls_control,
nls_args = power_curve_args$nls_args,
final_nrep = final_nrep,
final_R = final_R,
max_trials = max_trials,
ci_level = ci_level,
extendInt = extendInt,
power_tolerance_in_interval = power_tolerance_in_interval,
power_tolerance_in_final = power_tolerance_in_final),
control))
by_x_1 <- a_out$by_x_1
fit_1 <- a_out$fit_1
ci_hit <- a_out$ci_hit
x_tried <- a_out$x_tried
x_out <- a_out$x_out
power_out <- a_out$power_out
nrep_out <- a_out$nrep_out
ci_out <- a_out$ci_out
by_x_out <- a_out$by_x_out
i2 <- a_out$i2
solution_found <- a_out$solution_found
status <- a_out$status
technical <- list(iteration = a_out$iteration,
x_history = a_out$x_history,
reject_history = a_out$reject_history,
delta_tol = a_out$delta_tol,
last_k = a_out$last_k)
rm(a_out)
}
if ((algorithm == "bisection") && !solution_found) {
a_out <- do.call(alg_bisection,
c(list(
object = object,
x = x,
pop_es_name = pop_es_name,
target_power = target_power,
x_max = x_max,
x_min = x_min,
progress = progress,
x_include_interval = x_include_interval,
x_interval = x_interval,
simulation_progress = simulation_progress,
save_sim_all = save_sim_all,
is_by_x = is_by_x,
object_by_org = object_by_org,
final_nrep = final_nrep,
final_R = final_R,
extendInt = extendInt,
max_trials = max_trials,
R = attr(object, "args")$R,
digits = 3,
what = what,
goal = goal,
tol = tolerance
),
control))
by_x_1 <- a_out$by_x_1
fit_1 <- a_out$fit_1
ci_hit <- a_out$ci_hit
x_tried <- a_out$x_tried
x_out <- a_out$x_out
power_out <- a_out$power_out
nrep_out <- a_out$nrep_out
ci_out <- a_out$ci_out
by_x_out <- a_out$by_x_out
i2 <- a_out$i2
solution_found <- a_out$solution_found
status <- a_out$status
technical <- list(iteration = a_out$iteration,
x_history = a_out$x_history,
x_interval_history = a_out$x_interval_history,
f_interval_history = a_out$f_interval_history,
reject_history = a_out$reject_history,
f_history = a_out$f_history,
delta_tol = a_out$delta_tol,
last_k = a_out$last_k,
tol = a_out$tol)
rm(a_out)
}
# ** by_x_1 **
# The collection of all values tried and their results
# to be updated when new value is tried.
# Used after the end of the loop.
# ** fit_1 **
# The latest power curve
# To be updated whenever by_x_1 is updated.
# Used after the end of the loop.
if (progress) {
cat("\n\n--- Final Stage ---\n\n")
tmp <- format(Sys.time(), "%Y-%m-%d %X")
cat("- Start at", tmp, "\n")
cat("- Rejection Rates:\n\n")
tmp <- rejection_rates(by_x_1,
level = ci_level)
print(tmp)
cat("\n")
cat("- Estimated Power Curve:\n\n")
print(fit_1)
cat("\n")
}
# ==== Is solution found ====
# Is solution found?
# - The maximum number of replications reached.
if (goal == "ci_hit") {
# ==== Goal: ci_hit ====
if (isTRUE(ci_hit) && (nrep_out == final_nrep)) {
# ==== Solution found ====
# ** x_final, by_x_final, power_final, ci_final, nrep_final, i_final **
# The solution.
# ==== Create _final objects ====
x_final <- x_out
by_x_final <- by_x_out
power_final <- power_out
ci_final <- ci_out
nrep_final <- nrep_out
i_final <- i2
} else {
# ==== Solution not found ====
# No solution found.
# Set the NAs to denote this.
# ==== Set _final objects to NA====
x_final <- NA
by_x_final <- NA
power_final <- NA
ci_final <- NA
nrep_final <- NA
i_final <- NA
}
}
if (goal == "close_enough") {
# ==== Goal: close_enough ====
if (isTRUE(solution_found) && (nrep_out == final_nrep)) {
# ==== Solution found ====
# ** x_final, by_x_final, power_final, ci_final, nrep_final, i_final **
# The solution.
# ==== Create _final objects ====
x_final <- x_out
by_x_final <- by_x_out
power_final <- power_out
ci_final <- ci_out
nrep_final <- nrep_out
i_final <- i2
} else {
# ==== Solution not found ====
# No solution found.
# Set the NAs to denote this.
# ==== Set _final objects to NA====
x_final <- NA
by_x_final <- NA
power_final <- NA
ci_final <- NA
nrep_final <- NA
i_final <- NA
}
}
# ==== Compute extrapolated x? ====
# ** x_x **
# The estimated value based on power_curve.
# Used as a suggestion when no solution was found.
x_x <- NA
if (solution_found) {
# ==== Yes only if solution_found ====
x_tried <- switch(x,
n = as.numeric(names(by_x_1)),
es = sapply(by_x_1,
\(x) {attr(x, "pop_es_value")},
USE.NAMES = FALSE))
x_x <- estimate_x_range(power_x_fit = fit_1,
x = x,
target_power = target_power,
k = 1,
tolerance = 0,
power_min = .01,
power_max = .99,
interval = x_interval,
extendInt = "yes",
x_to_exclude = x_tried)
}
if (progress) {
if (solution_found && (nrep_out == final_nrep)) {
cat("\n")
x_out_str <- formatC(x_out,
digits = switch(x, n = 0, es = 4),
format = "f")
cat("- Final Value:", x_out_str, "\n\n")
cat("- Final Estimated Power:",
formatC(power_final, digits = 4, format = "f"), "\n")
cat("- Confidence Interval: [",
paste0(formatC(ci_final, digits = 4, format = "f"), collapse = "; "),
"]\n", sep = "")
cat("- CI Level: ",
formatC(ci_level*100, digits = 2, format = "f"), "%", "\n", sep = "")
} else {
cat("\n")
cat("- None of the value(s) examined",
"in the interval meet the goal.\n")
if (!is.na(x_x) &&
(goal == "ci_hit")) {
x_x_str <- formatC(x_x,
digits = switch(x, n = 0, es = 4),
format = "f")
cat("- The estimated required value is ",
x_x_str,
".\n", sep = "")
}
cat("\nTry changing the setup, such as:\n")
tmp <- switch(x,
n = paste0("[",
paste0(x_interval, collapse = ", "),
"]"),
es = paste0("[",
paste0(formatC(x_interval,
digits = 3,
format = "f"),
collapse = ", "),
"]"))
catwrap(paste0("- Changing 'x_interval' to a wider range to examine. ",
"(Current 'x_interval' is ",
tmp, ".)"),
exdent = 2)
catwrap(paste0("- Increasing the maximum number of trials by ",
"setting 'max_trials' to a larger number. ",
"(Current 'max_trials' is ",
max_trials,
".)"),
exdent = 2)
if (goal == "close_enough") {
catwrap(paste0("- Increasing the tolerance for the goal 'close_enough'. ",
"(Current 'tolerance' is ",
tolerance,
".)"),
exdent = 2)
}
}
}
# ==== Finalize the Output ====
my_call <- as.list(match.call())[-1]
args <- formals()
args <- utils::modifyList(args,
my_call)
args$object <- NULL
reject_1 <- rejection_rates(by_x_1,
level = ci_level)
time_end <- Sys.time()
out <- list(x = x,
pop_es_name = pop_es_name,
power4test_trials = by_x_1,
rejection_rates = reject_1,
x_tried = switch(x,
n = reject_1$n,
es = reject_1$es),
power_tried = reject_1$reject,
x_final = x_final,
power_final = power_final,
i_final = i_final,
ci_final = ci_final,
ci_level = ci_level,
nrep_final = nrep_final,
power_curve = fit_1,
target_power = target_power,
power_tolerance = power_tolerance_in_final,
x_estimated = x_x,
start = time_start,
end = time_end,
time_spent = difftime(time_end, time_start),
solution_found = solution_found,
what = what,
goal = goal,
args = args,
status = status,
technical = technical,
algorithm = algorithm,
call = match.call())
class(out) <- c("x_from_power", class(out))
return(out)
}
#' @rdname x_from_power
#'
#' @details
#'
#' The function [n_from_power()] is just
#' a wrapper of [x_from_power()], with
#' `x` set to `"n"`.
#'
#' @export
n_from_power <- function(object,
pop_es_name = NULL,
target_power = .80,
what = c("point", "ub", "lb"),
goal = switch(what,
point = "ci_hit",
ub = "close_enough",
lb = "close_enough"),
ci_level = .95,
tolerance = .02,
x_interval = c(50, 2000),
extendInt = NULL,
progress = TRUE,
simulation_progress = TRUE,
max_trials = 10,
final_nrep = 400,
final_R = 1000,
seed = NULL,
x_include_interval = FALSE,
check_es_interval = TRUE,
power_curve_args = list(power_model = NULL,
start = NULL,
lower_bound = NULL,
upper_bound = NULL,
nls_control = list(),
nls_args = list()),
save_sim_all = FALSE,
algorithm = NULL,
control = list()
) {
what <- match.arg(what)
goal <- match.arg(goal,
c("ci_hit", "close_enough"))
if ((goal == "ci_hit") &&
(what != "point")) {
what <- "point"
}
my_call <- match.call()
my_call$x <- "n"
my_call$what <- what
my_call$goal <- goal
my_call[[1]] <- quote(power4mome::x_from_power)
out <- eval(my_call,
envir = parent.frame())
out
}
#' @rdname x_from_power
#'
#' @details
#'
#' The function [n_region_from_power()] is just
#' a wrapper of [x_from_power()], with
#' `x` set to `"n"`, with two passes, one
#' with `what = "ub"` and one with
#' `what = "lb"`.
#'
#' @return
#' The function [n_region_from_power()]
#' returns a named list of two output of
#' [n_from_power()], of the class
#' `n_region_from_power`. The output
#' with `what = "ub"` is named `"below"`,
#' and the output with `what = "lb"` is
#' namd `"above"`.
#'
#' @export
n_region_from_power <- function(
object,
pop_es_name = NULL,
target_power = .80,
ci_level = .95,
tolerance = .02,
x_interval = c(50, 2000),
extendInt = NULL,
progress = TRUE,
simulation_progress = TRUE,
max_trials = 10,
final_nrep = 400,
final_R = 1000,
seed = NULL,
x_include_interval = FALSE,
check_es_interval = TRUE,
power_curve_args = list(power_model = NULL,
start = NULL,
lower_bound = NULL,
upper_bound = NULL,
nls_control = list(),
nls_args = list()),
save_sim_all = FALSE,
algorithm = NULL,
control = list()
) {
my_call <- match.call()
my_call$x <- "n"
my_call[[1]] <- quote(power4mome::x_from_power)
if (progress) {
tmp <- strwrap(paste0("Find the approximate region with power ",
"significantly below ",
target_power,
" ..."))
cat("\n")
cat("\n--- Phase 1 ---\n\n")
cat(tmp, sep = "\n")
cat("\n")
}
my_call$what <- "ub"
my_call$goal <- "close_enough"
out_ub <- eval(
my_call,
envir = parent.frame())
if (progress) {
tmp <- strwrap(paste0("Find the approximate region with power ",
"significantly above ",
target_power,
" ..."))
cat("\n")
cat("\n--- Phase 2 ---\n\n")
cat(tmp, sep = "\n")
cat("\n")
}
my_call$what <- "lb"
my_call$goal <- "close_enough"
tmp <- out_ub
tmp$what <- "lb"
tmp$goal <- "close_enough"
tmp$solution_found <- FALSE
tmp$call$what <- "lb"
tmp$call$goal <- "close_enough"
my_call2 <- my_call
my_call2$object <- tmp
out_lb <- eval(
my_call2,
envir = parent.frame())
out_lb$call <- my_call
out <- list(
below = out_ub,
above = out_lb
)
class(out) <- c("n_region_from_power", class(out))
attr(out, "call") <- match.call()
out
}
#' @rdname x_from_power
#'
#'
#' @param digits The number of digits
#' after the decimal when printing
#' the results.
#'
#' @param ... Optional arguments.
#' Not used for now.
#'
#' @details
#' The `print` method only prints
#' basic information. Call the
#' `summary` method of `x_from_power` objects
#' ([summary.x_from_power()]) and its
#' `print` method for detailed results
#'
#' @return
#' The `print`-method of `x_from_power`
#' objects returns the object `x`
#' invisibly.
#' It is called for its side effect.
#'
#' @export
print.x_from_power <- function(x,
digits = 3,
...) {
my_call <- x$call
cat("Call:\n")
print(my_call)
cat("\n")
solution_found <- !is.na(x$x_final)
predictor <- x$x
# cat("Predictor (x):",
# switch(predictor,
# n = "Sample Size",
# es = "Effect Size"),
# "\n")
# if (predictor == "es") {
# cat("Parameter Name (pop_es_name):",
# x$pop_es_name,
# "\n")
# }
goal <- x$goal
what <- x$what
algorithm <- x$algorithm
ci_level_str <- paste0(formatC(
x$ci_level * 100,
digits = 2,
format = "f"),
"%")
tmp1 <- c("Predictor(x):" =
switch(predictor,
n = "Sample Size",
es = "Effect Size"))
tmp1b <- c("Parameter:" =
switch(predictor,
n = "N/A",
es = x$pop_es_name))
tmp2 <- c("goal:" = goal)
tmp3 <- c("what:" = what)
tmp4 <- c("algorithm:" = algorithm)
tmp5 <- c("Level of confidence:" = ci_level_str)
tmp6 <- c("Target Power:" =
formatC(x$target_power, digits = digits, format = "f"))
tmp <- data.frame("Setting" = c(
tmp1,
tmp1b,
tmp2,
tmp3,
tmp4,
tmp5,
tmp6
))
print(tmp)
# cat("goal:", goal, "\n")
# cat("what:", what, "\n")
# cat("algorithm:", algorithm, "\n")
# cat("Level of Confidence (ci_level):",
# ci_level_str,
# "\n")
# cat("Target Power:",
# formatC(x$target_power, digits = digits, format = "f"),
# "\n")
if (solution_found) {
x_final_str <- formatC(x$x_final,
digits = switch(predictor,
n = 0,
es = digits),
format = "f")
cat("\n- Final Value of",
switch(x$x,
n = " Sample Size (n): ",
es = paste0("'", x$pop_es_name, "': ")),
x_final_str,
"\n\n",
sep = "")
ci_str <- paste0(
"[",
formatC(x$ci_final[1], digits = digits, format = "f"),
", ",
formatC(x$ci_final[2], digits = digits, format = "f"),
"]")
cat("- Final Estimated Power (CI): ",
formatC(x$power_final, digits = digits, format = "f"),
" ",
ci_str,
"\n",
sep = "")
# cat("- Confidence Interval of Power: [",
# formatC(x$ci_final[1], digits = digits, format = "f"),
# ", ",
# formatC(x$ci_final[2], digits = digits, format = "f"),
# "]\n",
# sep = "")
} else {
cat("\n- Solution not found.\n")
}
cat("\nCall `summary()` for detailed results.\n")
invisible(x)
}
#' @rdname x_from_power
#'
#' @return
#' The `print`-method of `x_from_power_region`
#' objects returns the object `x`
#' invisibly.
#' It is called for its side effect.
#'
#' @export
print.n_region_from_power <- function(
x,
digits = 3,
...) {
my_call <- attr(x, "call")
cat("Call:\n")
print(my_call)
cat("\n")
x_below <- x$below
x_above <- x$above
solution_found_below <- !is.na(x_below$x_final)
solution_found_above <- !is.na(x_above$x_final)
predictor <- x_below$x
goal <- x_below$goal
what <- x_below$what
algorithm <- x_below$algorithm
ci_level_str <- paste0(formatC(
x_below$ci_level * 100,
digits = 2,
format = "f"),
"%")
tmp1 <- c("Predictor(x)" = "Sample Size")
tmp2 <- c("Goal:" = "Power significantly below or above the target")
tmp4 <- c("algorithm:" = algorithm)
tmp5 <- c("Level of confidence:" = ci_level_str)
target_power_str <- formatC(x_below$target_power, digits = digits, format = "f")
tmp6 <- c("Target Power:" = target_power_str)
tmp <- data.frame("Setting" = c(
tmp1,
tmp2,
tmp4,
tmp5,
tmp6
))
print(tmp, right = FALSE)
if (solution_found_below) {
x_final_below_str <- formatC(x_below$x_final,
digits = switch(predictor,
n = 0,
es = digits),
format = "f")
ci_below_str <- paste0(
"[",
formatC(x_below$ci_final[1], digits = digits, format = "f"),
", ",
formatC(x_below$ci_final[2], digits = digits, format = "f"),
"]")
}
if (solution_found_above) {
x_final_above_str <- formatC(x_above$x_final,
digits = switch(predictor,
n = 0,
es = digits),
format = "f")
ci_above_str <- paste0(
"[",
formatC(x_above$ci_final[1], digits = digits, format = "f"),
", ",
formatC(x_above$ci_final[2], digits = digits, format = "f"),
"]")
}
solution_found <- solution_found_below &&
solution_found_above
if (solution_found_below ||
solution_found_above) {
cat("\nSolution: \n")
cat("\nApproximate region of sample sizes with power:\n")
if (solution_found) {
tmp <- paste0("- not significantly different from ",
target_power_str,
": ",
x_final_below_str,
" to ",
x_final_above_str)
cat(strwrap(tmp, exdent = 2), sep = "\n")
}
if (solution_found_below) {
tmp <- paste0("- significantly lower than ",
target_power_str,
": ",
x_final_below_str)
cat(strwrap(tmp, exdent = 2), sep = "\n")
}
if (solution_found_above) {
tmp <- paste0("- significantly higher than ",
target_power_str,
": ",
x_final_above_str)
cat(strwrap(tmp, exdent = 2), sep = "\n")
}
cat("\nConfidence intervals of the estimated power:\n")
if (solution_found_below) {
tmp <- paste0("- for the lower bound (",
x_final_below_str,
"): ",
ci_below_str)
cat(strwrap(tmp, exdent = 2), sep = "\n")
}
if (solution_found_above) {
tmp <- paste0("- for the upper bound (",
x_final_above_str,
"): ",
ci_above_str)
cat(strwrap(tmp, exdent = 2), sep = "\n")
}
if (!solution_found_below) {
cat("Solution not found for the lower region.")
}
if (!solution_found_above) {
cat("Solution not found for the upper region.")
}
} else {
cat("\n- Solution not found.\n")
}
cat("\nCall `summary()` for detailed results.\n")
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.