Spower {Spower} | R Documentation |
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,
lastSpower = NULL,
select = NULL,
control = list()
)
## S3 method for class 'Spower'
print(x, ...)
Arguments
... |
expression to use in the simulation that returns a Internally the first expression is passed to either For expected power computations the arguments to this expression can themselves
be specified as a function to reflect the prior uncertainty. For instance, if
|
power |
power level to use. If set to |
sig.level |
alpha level to use. If set to If the return of the supplied experiment is a
|
interval |
search interval to use when |
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
|
integer |
a logical value indicating whether the search iterations use integers or doubles. If missing, automatically set to |
parallel |
for parallel computing for slower simulation experiments
(see |
cl |
see |
packages |
see |
ncores |
see |
predCI |
predicting confidence interval level
(see |
predCI.tol |
predicting confidence interval consistency tolerance
for stochastic root solver convergence (see |
verbose |
logical; should information be printed to the console? |
check.interval |
logical; check the interval range validity
(see |
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 |
lastSpower |
a previously returned Note that if the object was not stored use |
select |
a character vector indicating which elements to
extract from the provided stimulation experiment function. By default, all elements
from the provided function will be used, however if the provided function contains
information not relevant to the power computations (e.g., parameter estimates,
standard errors, etc) then these should be ignored. To extract the complete
results post-analysis use |
control |
a list of control parameters to pass to
|
x |
object of class |
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 this results in observed/retrospective power (not recommended), while if only sample size is included as the observed quantity, but the effect sizes are treated as unknown, then this results in post-hoc power (Cohen, 1988)
- 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
(argumentbeta_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 rphilip.chalmers@gmail.com
See Also
update
, SpowerCurve
,
getLastSpower
, is.CI_within
,
is.outside_CI
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 (not run)
# p_t.test(n = 50, d = .5) |> Spower(replications=30000)
# alternatively, increase precision from previous object.
# Here we add 20000 more replications on top of the previous 10000
p_t.test(n = 50, d = .5) |>
Spower(replications=20000, lastSpower=out) -> out2
out2$REPLICATIONS # total of 30000 replications for estimate
# 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)
# Similiar to above for precision improvements, however letting
# the root solver continue searching from an early search history.
# Usually a good idea to increase the maxiter and lower the predCI.tol
p_t.test(n = NA, d = .5) |>
Spower(power=.8, interval=c(2,500), lastSpower=out,
maxiter=200, predCI.tol=.008) #starts at last iteration in "out"
# 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
###############################
# Using CIs instead of p-values
###############################
# CI test returning TRUE if psi0 is outside the 95% CI
ci_ind.t.test <- function(n, d, psi0=0, conf.level=.95){
g1 <- rnorm(n)
g2 <- rnorm(n, mean=d)
CI <- t.test(g2, g1, var.equal=TRUE,conf.level=conf.level)$conf.int
is.outside_CI(psi0, CI)
}
# returns logical
ci_ind.t.test(n=100, d=.2)
ci_ind.t.test(n=100, d=.2)
# simulated prospective power
ci_ind.t.test(n=100, d=.2) |> Spower()
# compare to pwr package
pwr::pwr.t.test(n=100, d=.2)
############################
# Equivalence test power using CIs
#
# H0: population d is outside interval [LB, UB] (not tolerably equivalent)
# H1: population d is within interval [LB, UB] (tolerably equivalent)
# CI test returning TRUE if CI is within tolerable equivalence range (tol)
ci_equiv.t.test <- function(n, d, tol, conf.level=.95){
g1 <- rnorm(n)
g2 <- rnorm(n, mean=d)
CI <- t.test(g2, g1, var.equal=TRUE,conf.level=conf.level)$conf.int
is.CI_within(CI, tol)
}
# evaluate if CI is within tolerable interval (tol)
ci_equiv.t.test(n=1000, d=.2, tol=c(.1, .3))
# simulated prospective power
ci_equiv.t.test(n=1000, d=.2, tol=c(.1, .3)) |> Spower()
# higher power with larger N (more precision) or wider tol interval
ci_equiv.t.test(n=2000, d=.2, tol=c(.1, .3)) |> Spower()
ci_equiv.t.test(n=1000, d=.2, tol=c(.1, .5)) |> Spower()
####
# superiority test (one-tailed)
# H0: population d is less than LB (not superior)
# H1: population d is greater than LB (superior)
# set upper bound to Inf as it's not relevant, and reduce conf.level
# to reflect one-tailed test
ci_equiv.t.test(n=1000, d=.2, tol=c(.1, Inf), conf.level=.90) |>
Spower()
# higher LB means greater requirement for defining superiority (less power)
ci_equiv.t.test(n=1000, d=.2, tol=c(.15, Inf), conf.level=.90) |>
Spower()