Spower: Simulation-based Power Analyses
In Spower: Power Analyses using Monte Carlo Simulations
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.
The function SpowerBatch, on the other hand, can be used to
run Spower across
different simulation combinations, returning a list of results instead.
Can also be used as a pre-computing step before using
SpowerCurve, and shares the same syntax specification (see
SpowerCurve for further examples).
SpowerCurve draws power curves that either a) estimate the power given a
set of varying conditions or b) solves a set of root conditions
given fixed values of power. Confidence/prediction intervals are
included in the output to reflect the estimate uncertainties, though note
that fewer replications/iterations are used compared to
Spower as the goal is visualization of competing
variable inputs rather than precision of a given input.
Usage
Spower(
...,
power = NA,
sig.level = 0.05,
interval,
beta_alpha,
sig.direction = "below",
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, ...)
## S3 method for class 'Spower'
as.data.frame(x, ...)
SpowerBatch(
...,
interval = NULL,
power = NA,
sig.level = 0.05,
beta_alpha = NULL,
sig.direction = "below",
replications = 10000,
integer,
parallel = FALSE,
cl = NULL,
ncores = parallelly::availableCores(omit = 1L),
predCI = 0.95,
predCI.tol = 0.01,
verbose = TRUE,
check.interval = FALSE,
maxiter = 150,
wait.time = NULL,
select = NULL,
control = list()
)
## S3 method for class 'SpowerBatch'
print(x, ...)
## S3 method for class 'SpowerBatch'
as.data.frame(x, ...)
SpowerCurve(
...,
interval = NULL,
power = NA,
sig.level = 0.05,
sig.direction = "below",
replications = 2500,
integer,
plotCI = TRUE,
plotly = TRUE,
parallel = FALSE,
cl = NULL,
ncores = parallelly::availableCores(omit = 1L),
predCI = 0.95,
predCI.tol = 0.01,
verbose = TRUE,
check.interval = FALSE,
maxiter = 50,
wait.time = NULL,
select = NULL,
batch = NULL,
control = list()
)
Arguments
...
expression to use in the simulation that returns a numeric
vector containing either the p-value (under the null hypothesis), the
probability of the alternative hypothesis in the Bayesian setting,
where the first numeric value
in this vector is treated as the focus for all analyses other than prospective/post-hoc power.
This corresponds to the alpha value used to flag samples
as 'significant' when evaluating the null hypothesis
(via p-values; P(D|H_0)),
where any returned p-value less that sig.level indicates significance.
However, if sig.direction = 'above' then only values
above sig.level are flagged as significant, which is useful
in Bayesian posterior probability contexts that focus on the alternative
hypothesis, P(H_1|D).
Alternatively, a logical vector can be returned (e.g., when using confidence intervals (CIs) or
evaluating regions of practical equivalence (ROPEs)), where the average of these
TRUE/FALSE vector corresponds to the empirical power.
For SpowerCurve and SpowerBatch,
first expression input must be identical to ... in
Spower, while the remaining named inputs must match the arguments
to this expression to indicate which variables should be modified in the
resulting power curves. Providing NA values is also supported to
solve the missing component. Note that only the first three named
arguments in SpowerCurve will be plotted using
the x-y, colour, and facet wrap aesthetics, respectively. However,
if necessary the data can be extracted for further visualizations via
ggplot_build to provide more customized control
power
power level to use. If set to NA (default)
then the empirical power will be estimated given the fixed ... inputs
(e.g., for prospective/post-hoc power analysis). For
SpowerCurve and SpowerBatch this can be a vector
sig.level
alpha level to use (default is .05). If set to NA then the value will
be estimated given the fixed conditions input
(e.g., for criterion power analysis). Only used when the value returned
from the experiment is a numeric (e.g., a p-value, or a
posterior probability; see sig.direction).
If the return of the supplied experiment is a
logical then this argument will be entirely ignored. As such,
arguments such as conf.level should be included
in the simulation experiment definition itself
to indicate the explicit inferential
criteria, and so that this argument can be manipulated should the need arise.
interval
required search interval to use when SimSolve is called
to perform stochastic root solving.
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 q = \beta/\alpha = 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
sig.direction
a character vector that is either 'below'
(default) or 'above' to indicate which direction relative to
sig.level is considered significant. This is useful, for instance,
when forming cutoffs for Bayesian
posterior probabilities organized to show support
for the hypothesis of interest (P(H_1|D)). As an example,
setting sig.level = .95 with sig.direction = 'above'
flags a sample as 'significant' whenever the
posterior probability is greater than .95.
replications
number of replications to use when
runSimulation is required. Default is 10000,
though set to 2500 for SpowerCurve
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
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.
Default is 150, though set to 50 for SpowerCurve
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
lastSpower
a previously returned Spower object to be updated.
Use this if you want to continue where an estimate left off but wish to increase the
precision (e.g., by adding more replications, or by letting the stochastic root solver
continue searching).
Note that if the object was not stored use getLastSpower
to obtain the last estimated power object
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 SimResults to allow manual
summarizing of the stored results (applicable only with prospective/post-hoc power)
control
a list of control parameters to pass to
runSimulation or SimSolve
x
object of class 'Spower'. If SpowerBatch were used
the this will be a list
plotCI
logical; include confidence/prediction intervals in plots?
plotly
logical; draw the graphic into the interactive plotly
interface? If FALSE the ggplot2 object will be returned instead
batch
if SpowerBatch were previously used to perform the computations
then this information can be provided to this batch argument to avoid
recomputing
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 (argument beta_alpha)
To understand how the package is structured, the first expression in the ... argument,
which contains the simulation experiment definition for a single sample,
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 the simulation experiment arguments can
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
Value
an invisible tibble/data.frame-type object of
class 'Spower' containing the power results from the
simulation experiment
a ggplot2 object automatically rendered with
plotly for interactivity
Author(s)
Phil Chalmers rphilip.chalmers@gmail.com
See Also
update, SpowerCurve,
getLastSpower, is.CI_within,
is.outside_CI
Spower, SpowerBatch
Examples
############################
# Independent samples t-test
############################
# Internally defined p_t.test function
args(p_t.test) # missing arguments required
# 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()
##############################################
# SpowerBatch() examples
##############################################
## Not run:
# estimate power given varying sample sizes
p_t.test(d=0.2) |>
SpowerBatch(n=c(30, 90, 270, 550), replications=1000) -> nbatch
nbatch
# can be stacked to view the output as data.frame
as.data.frame(nbatch)
# plot with SpowerCurve()
SpowerCurve(batch=nbatch)
# equivalent, but re-runs the computations
p_t.test(d=0.2) |> SpowerCurve(n=c(30, 90, 270, 550), replications=1000)
# estimate power given varying sample sizes and effect size
p_t.test() |> SpowerBatch(n=c(30, 90, 270, 550),
d=c(.2, .5, .8), replications=1000) -> ndbatch
ndbatch
# plot with SpowerCurve()
SpowerCurve(batch=ndbatch)
## End(Not run)
##############################################
# SpowerCurve() examples
##############################################
# estimate power given varying sample sizes
gg <- p_t.test(d=0.2) |> SpowerCurve(n=c(30, 90, 270, 550))
# Output is a ggplot2 (rendered with plotly by default); hence, can be modified
library(ggplot2)
gg + geom_text(aes(label=power), size=5, colour='red', nudge_y=.05) +
ylab(expression(1-beta)) + theme_grey()
# Increase precision by using 10000 replications. Parallel computations
# generally recommended in this case to save time
p_t.test(d=0.2) |> SpowerCurve(n=c(30, 90, 270, 550), replications=10000)
# estimate sample sizes given varying power
p_t.test(n=NA, d=0.2) |>
SpowerCurve(power=c(.2, .4, .6, .8), interval=c(10, 1000))
# get information from last printed graphic instead of saving
gg <- last_plot()
gg + coord_flip() # flip coordinates to put power on y-axis
# estimate power varying d
p_t.test(n=50) |> SpowerCurve(d=seq(.1, 1, by=.2))
# estimate d varying power
p_t.test(n=50, d=NA) |>
SpowerCurve(power=c(.2, .4, .6, .8), interval=c(.01, 1))
#####
# vary two inputs instead of one (second input uses colour aesthetic)
p_t.test() |> SpowerCurve(n=c(30, 90, 270, 550),
d=c(.2, .5, .8))
# extract data for alternative presentations
build <- ggplot_build(last_plot())
build
df <- build$plot$data
head(df)
ggplot(df, aes(n, power, linetype=d)) + geom_line()
# vary three arguments (third uses facet_wrap ... any more than that and
# you're on your own!)
p_t.test() |> SpowerCurve(n=c(30, 90, 270, 550),
d=c(.2, .5, .8),
var.equal=c(FALSE, TRUE))
########################################
# If objects were precomputed using SpowerBatch() then
# these can be plotted instead
p_t.test(d=0.2) |>
SpowerBatch(n=c(30, 90, 270, 550), replications=1000) -> nbatch
nbatch
as.data.frame(nbatch)
# plot the results, but avoid further computations
SpowerCurve(batch=nbatch)
Spower documentation built on Sept. 9, 2025, 5:46 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| 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.
The function SpowerBatch, on the other hand, can be used to
run Spower across
different simulation combinations, returning a list of results instead.
Can also be used as a pre-computing step before using
SpowerCurve, and shares the same syntax specification (see
SpowerCurve for further examples).
SpowerCurve draws power curves that either a) estimate the power given a
set of varying conditions or b) solves a set of root conditions
given fixed values of power. Confidence/prediction intervals are
included in the output to reflect the estimate uncertainties, though note
that fewer replications/iterations are used compared to
Spower as the goal is visualization of competing
variable inputs rather than precision of a given input.
Usage
Spower(
...,
power = NA,
sig.level = 0.05,
interval,
beta_alpha,
sig.direction = "below",
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, ...)
## S3 method for class 'Spower'
as.data.frame(x, ...)
SpowerBatch(
...,
interval = NULL,
power = NA,
sig.level = 0.05,
beta_alpha = NULL,
sig.direction = "below",
replications = 10000,
integer,
parallel = FALSE,
cl = NULL,
ncores = parallelly::availableCores(omit = 1L),
predCI = 0.95,
predCI.tol = 0.01,
verbose = TRUE,
check.interval = FALSE,
maxiter = 150,
wait.time = NULL,
select = NULL,
control = list()
)
## S3 method for class 'SpowerBatch'
print(x, ...)
## S3 method for class 'SpowerBatch'
as.data.frame(x, ...)
SpowerCurve(
...,
interval = NULL,
power = NA,
sig.level = 0.05,
sig.direction = "below",
replications = 2500,
integer,
plotCI = TRUE,
plotly = TRUE,
parallel = FALSE,
cl = NULL,
ncores = parallelly::availableCores(omit = 1L),
predCI = 0.95,
predCI.tol = 0.01,
verbose = TRUE,
check.interval = FALSE,
maxiter = 50,
wait.time = NULL,
select = NULL,
batch = NULL,
control = list()
)
Arguments
... |
expression to use in the simulation that returns a Alternatively, a For |
power |
power level to use. If set to |
sig.level |
alpha level to use (default is If the return of the supplied experiment is a
|
interval |
required 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 |
sig.direction |
a character vector that is either |
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.
Default is 150, though set to 50 for |
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 |
plotCI |
logical; include confidence/prediction intervals in plots? |
plotly |
logical; draw the graphic into the interactive |
batch |
if |
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)
To understand how the package is structured, the first expression in the ... argument,
which contains the simulation experiment definition for a single sample,
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 the simulation experiment arguments can
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
Value
an invisible tibble/data.frame-type object of
class 'Spower' containing the power results from the
simulation experiment
a ggplot2 object automatically rendered with
plotly for interactivity
Author(s)
Phil Chalmers rphilip.chalmers@gmail.com
See Also
update, SpowerCurve,
getLastSpower, is.CI_within,
is.outside_CI
Spower, SpowerBatch
Examples
############################
# Independent samples t-test
############################
# Internally defined p_t.test function
args(p_t.test) # missing arguments required
# 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()
##############################################
# SpowerBatch() examples
##############################################
## Not run:
# estimate power given varying sample sizes
p_t.test(d=0.2) |>
SpowerBatch(n=c(30, 90, 270, 550), replications=1000) -> nbatch
nbatch
# can be stacked to view the output as data.frame
as.data.frame(nbatch)
# plot with SpowerCurve()
SpowerCurve(batch=nbatch)
# equivalent, but re-runs the computations
p_t.test(d=0.2) |> SpowerCurve(n=c(30, 90, 270, 550), replications=1000)
# estimate power given varying sample sizes and effect size
p_t.test() |> SpowerBatch(n=c(30, 90, 270, 550),
d=c(.2, .5, .8), replications=1000) -> ndbatch
ndbatch
# plot with SpowerCurve()
SpowerCurve(batch=ndbatch)
## End(Not run)
##############################################
# SpowerCurve() examples
##############################################
# estimate power given varying sample sizes
gg <- p_t.test(d=0.2) |> SpowerCurve(n=c(30, 90, 270, 550))
# Output is a ggplot2 (rendered with plotly by default); hence, can be modified
library(ggplot2)
gg + geom_text(aes(label=power), size=5, colour='red', nudge_y=.05) +
ylab(expression(1-beta)) + theme_grey()
# Increase precision by using 10000 replications. Parallel computations
# generally recommended in this case to save time
p_t.test(d=0.2) |> SpowerCurve(n=c(30, 90, 270, 550), replications=10000)
# estimate sample sizes given varying power
p_t.test(n=NA, d=0.2) |>
SpowerCurve(power=c(.2, .4, .6, .8), interval=c(10, 1000))
# get information from last printed graphic instead of saving
gg <- last_plot()
gg + coord_flip() # flip coordinates to put power on y-axis
# estimate power varying d
p_t.test(n=50) |> SpowerCurve(d=seq(.1, 1, by=.2))
# estimate d varying power
p_t.test(n=50, d=NA) |>
SpowerCurve(power=c(.2, .4, .6, .8), interval=c(.01, 1))
#####
# vary two inputs instead of one (second input uses colour aesthetic)
p_t.test() |> SpowerCurve(n=c(30, 90, 270, 550),
d=c(.2, .5, .8))
# extract data for alternative presentations
build <- ggplot_build(last_plot())
build
df <- build$plot$data
head(df)
ggplot(df, aes(n, power, linetype=d)) + geom_line()
# vary three arguments (third uses facet_wrap ... any more than that and
# you're on your own!)
p_t.test() |> SpowerCurve(n=c(30, 90, 270, 550),
d=c(.2, .5, .8),
var.equal=c(FALSE, TRUE))
########################################
# If objects were precomputed using SpowerBatch() then
# these can be plotted instead
p_t.test(d=0.2) |>
SpowerBatch(n=c(30, 90, 270, 550), replications=1000) -> nbatch
nbatch
as.data.frame(nbatch)
# plot the results, but avoid further computations
SpowerCurve(batch=nbatch)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.