Spower: Simulation-based Power Analyses

In Spower: Power Analyses using Monte Carlo Simulations

Simulation-based Power Analyses

Description

General purpose function that provides power-focused estimates for a priori, prospective/post-hoc, compromise, sensitivity, and criterion power analysis. Function provides a general wrapper to the SimDesign package's runSimulation and SimSolve functions. As such, parallel processing is automatically supported, along with progress bars, confidence/prediction intervals for the results estimates, safety checks, and more.

Usage

Spower(
  ...,
  power = NA,
  sig.level = 0.05,
  interval,
  beta_alpha,
  replications = 10000,
  integer,
  parallel = FALSE,
  cl = NULL,
  packages = NULL,
  ncores = parallelly::availableCores(omit = 1L),
  predCI = 0.95,
  predCI.tol = 0.01,
  verbose = TRUE,
  check.interval = FALSE,
  maxiter = 150,
  wait.time = NULL,
  control = list()
)

## S3 method for class 'Spower'
print(x, ...)

Arguments

...	expression to use in the simulation that returns a numeric vector containing only p-value information, where the first p-value in this vector is treated as the focus for all analyses other than prospective/post-hoc power, or a similarly structure logical vector when utilizing confidence intervals (CIs). Internally the first expression is passed to either SimSolve or runSimulation depending on which element (including the power and sig.level arguments) is set to NA. For instance, Spower(p_t.test(n=50, d=.5)) will perform a prospective/post-hoc power evaluation since power = NA by default, while Spower(p_t.test(n=NA, d=.5), power = .80) will perform an a priori power analysis to solve the missing n argument. For expected power computations the arguments to this expression can themselves be specified as a function to reflect the prior uncertainty. For instance, if d_prior <- function() rnorm(1, mean=.5, sd=1/8) then Spower(p_t.test(n=50, d=d_prior()) will compute the expected power over the prior sampling distribution for d
power	power level to use. If set to NA then the empirical power will be estimated given the fixed ... inputs (e.g., for prospective/post-hoc power analysis)
sig.level	alpha level to use. If set to NA then the empirical alpha will be estimated given the fixed conditions input (e.g., for criterion power analysis)
interval	search interval to use when SimSolve is required. Note that for compromise analyses, where the sig.level is set to NA, if not set explicitly then the interval will default to c(0,1)
beta_alpha	(optional) ratio to use in compromise analyses corresponding to the Type II errors (beta) over the Type I error (alpha). Ratios greater than 1 indicate that Type I errors are worse than Type II, while ratios less than one the opposite. A ratio equal to 1 gives an equal trade-off between Type I and Type II errors
replications	number of replications to use when runSimulation is required
integer	a logical value indicating whether the search iterations use integers or doubles. If missing, automatically set to FALSE if interval contains non-integer numbers or the range is less than 5, as well as when sig.level = NA, though in general this should be set explicitly
parallel	for parallel computing for slower simulation experiments (see runSimulation for details).
cl	see runSimulation
packages	see runSimulation
ncores	see runSimulation
predCI	predicting confidence interval level (see SimSolve)
predCI.tol	predicting confidence interval consistency tolerance for stochastic root solver convergence (see SimSolve). Default converges when the power rate CI is consistently within .01/2 of the target power
verbose	logical; should information be printed to the console?
check.interval	logical; check the interval range validity (see SimSolve). Disabled by default
maxiter	maximum number of stochastic root-solving iterations
wait.time	(optional) argument to indicate the time to wait (specified in minutes if supplied as a numeric vector). See SimSolve for details and See timeFormater for further specifications
control	a list of control parameters to pass to runSimulation or SimSolve
x	object of class 'Spower'

Details

Five types of power analysis flavors can be performed with Spower, which are triggered based on which supplied input is set to missing (NA):

A Priori

Solve for a missing sample size component (e.g., n) to achieve a specific target power rate

Prospective (and Post-hoc)

Estimate the power rate given a set of fixed conditions. If estimates of effect sizes and other empirical characteristics (e.g., observed sample size) are supplied instead this results in post-hoc/observed/retrospective power (not recommended)

Sensitivity

Solve a missing effect size value as a function of the other supplied constant components

Criterion

Solve the error rate (argument sig.level) as a function of the other supplied constant components

Compromise

Solve a Type I/Type II error trade-off ratio as a function of the other supplied constant components and the target ratio q = \beta/\alpha (argument beta_alpha)

Prospective and compromise analyses utilize the runSimulation function, while the remaining three approaches utilize the stochastic root solving methods in the function SimSolve. See the example below for a demonstration with an independent samples t-test analysis.

Value

an invisible tibble/data.frame-type object of class 'Spower' containing the power results from the simulation experiment

Author(s)

Phil Chalmers [e-mail]

See Also

update, SpowerCurve, getLastSpower

Examples

############################
# Independent samples t-test
############################

# Internally defined p_t.test function
args(p_t.test)    # missing arguments required for Spower()
# help(p_t.test)  # additional information

# p_* functions generate data and return single p-value
p_t.test(n=50, d=.5)
p_t.test(n=50, d=.5)

# test that it works
Spower(p_t.test(n = 50, d = .5), replications=10)

# also behaves naturally with a pipe
p_t.test(n = 50, d = .5) |> Spower(replications=10)



# Estimate power given fixed inputs (prospective power analysis)
out <- Spower(p_t.test(n = 50, d = .5))
summary(out)   # extra information
as.data.frame(out)  # coerced to data.frame

# increase precision
p_t.test(n = 50, d = .5) |> Spower(replications=30000)

# previous analysis not stored to object, but can be retrieved
out <- getLastSpower()
out   # as though it were stored from Spower()

# Same as above, but executed with multiple cores (not run)
p_t.test(n = 50, d = .5) |>
   Spower(replications=30000, parallel=TRUE, ncores=2)

# Solve N to get .80 power (a priori power analysis)
p_t.test(n = NA, d = .5) |>
  Spower(power=.8, interval=c(2,500)) -> out
summary(out)  # extra information
plot(out)
plot(out, type = 'history')

# total sample size required
ceiling(out$n) * 2

# same as above, but in parallel with 2 cores
out.par <- p_t.test(n = NA, d = .5) |>
  Spower(power=.8, interval=c(2,500), parallel=TRUE, ncores=2)
summary(out.par)

# similar information from pwr package
(pwr <- pwr::pwr.t.test(d=.5, power=.80))
ceiling(pwr$n) * 2

# If greater precision is required and the user has a specific amount of time
# they are willing to wait (e.g., 5 minutes) then wait.time can be used. Below
# estimates root after searching for 1 minute, and run in parallel
#  with 2 cores (not run)
p_t.test(n = NA, d = .5) |>
  Spower(power=.8, interval=c(2,500), wait.time='1', parallel=TRUE, ncores=2)

# Solve d to get .80 power (sensitivity power analysis)
p_t.test(n = 50, d = NA) |> Spower(power=.8, interval=c(.1, 2))
pwr::pwr.t.test(n=50, power=.80) # compare

# Solve alpha that would give power of .80 (criterion power analysis)
#    interval not required (set to interval = c(0, 1))
p_t.test(n = 50, d = .5) |> Spower(power=.80, sig.level=NA)

# Solve beta/alpha ratio to specific error trade-off constant
#   (compromise power analysis)
out <- p_t.test(n = 50, d = .5) |> Spower(beta_alpha = 2)
with(out, (1-power)/sig.level)   # solved ratio

# update beta_alpha criteria without re-simulating
(out2 <- update(out, beta_alpha=4))
with(out2, (1-power)/sig.level)   # solved ratio

##############
# Power Curves
##############

# SpowerCurve() has similar input, though requires varying argument
p_t.test(d=.5) |> SpowerCurve(n=c(30, 60, 90))

# solve n given power and plot
p_t.test(n=NA, d=.5) |> SpowerCurve(power=c(.2, .5, .8), interval=c(2,500))

# multiple varying components
p_t.test() |> SpowerCurve(n=c(30,60,90), d=c(.2, .5, .8))

################
# Expected Power
################

# Expected power computed by including effect size uncertainty.
# For instance, belief is that the true d is somewhere around d ~ N(.5, 1/8)
dprior <- function(x, mean=.5, sd=1/8) dnorm(x, mean=mean, sd=sd)
curve(dprior, -1, 2, main=expression(d %~% N(0.5, 1/8)),
      xlab='d', ylab='density')

# For Spower, define prior sampler for specific parameter(s)
d_prior <- function() rnorm(1, mean=.5, sd=1/8)
d_prior(); d_prior(); d_prior()

# Replace d constant with d_prior to compute expected power
p_t.test(n = 50, d = d_prior()) |> Spower()

# A priori power analysis using expected power
p_t.test(n = NA, d = d_prior()) |>
  Spower(power=.8, interval=c(2,500))
pwr::pwr.t.test(d=.5, power=.80) # expected power result higher than fixed d


###############
# Customization
###############

#   Make edits to the function for customization
if(interactive()){
    p_my_t.test <- edit(p_t.test)
    args(p_my_t.test)
    body(p_my_t.test)
}

# Alternatively, define a custom function (potentially based on the template)
p_my_t.test <- function(n, d, var.equal=FALSE, n2_n1=1, df=10){

    # Welch power analysis with asymmetric distributions
    # group2 as large as group1 by default

    # degree of skewness controlled via chi-squared distribution's df
    group1 <- rchisq(n, df=df)
    group1 <-  (group1 - df) / sqrt(2*df)   # Adjusted mean to 0, sd = 1
    group2 <- rnorm(n*n2_n1, mean=d)
    dat <- data.frame(group = factor(rep(c('G1', 'G2'),
                                     times = c(n, n*n2_n1))),
    				  DV = c(group1, group2))
    obj <- t.test(DV ~ group, dat, var.equal=var.equal)
    p <- obj$p.value
    p
}

# Solve N to get .80 power (a priori power analysis), using defaults
p_my_t.test(n = NA, d = .5, n2_n1=2) |>
  Spower(power=.8, interval=c(2,500)) -> out

# total sample size required
with(out, ceiling(n) + ceiling(n * 2))

# Solve N to get .80 power (a priori power analysis), assuming
#   equal variances, group2 2x as large as group1, large skewness
p_my_t.test(n = NA, d=.5, var.equal=TRUE, n2_n1=2, df=3) |>
  Spower(power=.8, interval=c(30,100)) -> out2

# total sample size required
with(out2, ceiling(n) + ceiling(n * 2))

# prospective power, can be used to extract the adjacent information
p_my_t.test(n = 100, d = .5) |> Spower() -> post

Spower documentation built on April 4, 2025, 5:11 a.m.