Title:	Power Analyses using Monte Carlo Simulations
Version:	0.3.1
Description:	Provides a general purpose simulation-based power analysis API for routine and customized simulation experimental designs. The package focuses exclusively on Monte Carlo simulation experiment variants of (expected) prospective power analyses, criterion analyses, compromise analyses, sensitivity analyses, and a priori/post-hoc analyses. The default simulation experiment functions defined within the package provide stochastic variants of the power analysis subroutines in G*Power 3.1 (Faul, Erdfelder, Buchner, and Lang, 2009) <doi:10.3758/brm.41.4.1149>, along with various other parametric and non-parametric power analysis applications (e.g., mediation analyses). Additional functions for building empirical power curves, reanalyzing simulation information, and for increasing the precision of the resulting power estimates are also included, each of which utilize similar API structures.
Depends:	SimDesign (≥ 2.20.0), R (≥ 4.1.0), stats
Imports:	cocor, car, polycor, parallelly, methods, ggplot2, plotly, lavaan, EnvStats
Suggests:	knitr, rmarkdown, bookdown, VGAM, copula, pwr
VignetteBuilder:	knitr
License:	GPL (≥ 3)
Encoding:	UTF-8
Repository:	CRAN
Maintainer:	Phil Chalmers <rphilip.chalmers@gmail.com>
URL:	https://github.com/philchalmers/Spower
BugReports:	https://github.com/philchalmers/Spower/issues?state=open
RoxygenNote:	7.3.2
NeedsCompilation:	no
Packaged:	2025-07-16 19:57:01 UTC; phil
Author:	Phil Chalmers [aut, cre]
Date/Publication:	2025-07-17 09:10:02 UTC

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

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)

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()

Draw power curve from simulation functions

Description

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

SpowerCurve(
  ...,
  interval = NULL,
  power = NA,
  sig.level = 0.05,
  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,
  control = list()
)

Arguments

Value

a ggplot2 object automatically rendered with plotly for interactivity

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

Spower

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))

Get previously evaluated Spower execution

Description

If the result of Spower was not stored into an object this function will retrieve the last evaluation.

Usage

getLastSpower()

Value

the last object returned from Spower

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

Spower

Evaluate whether a confidence interval is within a tolerable interval

Description

Return TRUE if an estimated confidence interval falls within a tolerable interval range. Typically used for equivalence, superiority, or non-inferiority testing.

Usage

is.CI_within(CI, interval)

Arguments

Value

logical

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

is.outside_CI, Spower

Examples

CI <- c(.2, .4)
LU <- c(.1, .3)
is.CI_within(CI, LU)        # not within tolerable interval
is.CI_within(CI, c(0, .5))  # is within wider interval

# complement indicates if CI is outside interval
!is.CI_within(CI, LU)

#####
# for superiority test
is.CI_within(CI, c(.1, Inf))  # CI is within tolerable interval

# for inferiority test
is.CI_within(CI, c(-Inf, .3))  # CI is not within tolerable interval

Evaluate whether parameter is outside a given confidence interval

Description

Returns TRUE if parameter reflecting a null hypothesis falls outside a given confidence interval. This is an alternative approach to writing an experiment that returns a p-value.

Usage

is.outside_CI(P0, CI)

Arguments

Value

logical

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

is.CI_within, Spower

Examples

p0 <- .3
CI <- c(.2, .4)
is.outside_CI(p0, CI)

# complement indicates if p0 is within CI
!is.outside_CI(p0, CI)

p-value from comparing two or more correlations simulation

Description

Function utilizes cocor to perform correlation comparison for independent, overlapping, and non-overlapping designs.

Usage

p_2r(
  n,
  r.ab1,
  r.ab2,
  r.ac1,
  r.ac2,
  r.bc1,
  r.bc2,
  r.ad1,
  r.ad2,
  r.bd1,
  r.bd2,
  r.cd1,
  r.cd2,
  n2_n1 = 1,
  two.tailed = TRUE,
  type = c("independent", "overlap", "nonoverlap"),
  test = "fisher1925",
  gen_fun = gen_2r,
  ...
)

gen_2r(n, R, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

# independent (same x-y pairing across groups)
p_2r(100, r.ab1=.5, r.ab2=.6)



   # estimate empirical power
   p_2r(n=100, r.ab1=.5, r.ab2=.6) |> Spower()

   # estimate n required to reach 80% power
   p_2r(n=NA, r.ab1=.5, r.ab2=.6) |>
        Spower(power=.80, interval=c(100, 5000))



# overlap (same y, different xs)
p_2r(100, r.ab1=.5, r.ab2=.7,
          r.ac1=.3, r.ac2=.3,
          r.bc1=.2, r.bc2=.2, type = 'overlap')

# nonoverlap (different ys, different xs)
p_2r(100, r.ab1=.5, r.ab2=.6,
          r.ac1=.3, r.ac2=.3,
          r.bc1=.2, r.bc2=.2,
          r.ad1=.2, r.ad2=.2,
          r.bd1=.4, r.bd2=.4,
          r.cd1=.2, r.cd2=.2,
          type = 'nonoverlap')

p-value from one-way ANOVA simulation

Description

Generates continuous multi-sample data to be analyzed by a one-way ANOVA, and return a p-value. Uses the function oneway.test to perform the analyses. The data and associated test assume that the conditional observations are normally distributed and have have equal variance by default, however these may be modified.

Usage

p_anova.test(
  n,
  k,
  f,
  n.ratios = rep(1, k),
  two.tailed = TRUE,
  var.equal = TRUE,
  means = NULL,
  sds = NULL,
  gen_fun = gen_anova.test,
  ...
)

gen_anova.test(n, k, f, n.ratios = rep(1, k), means = NULL, sds = NULL, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_anova.test

Examples

# n=50 in 3 groups, "medium" effect size
p_anova.test(50, k=3, f=.25)

# explicit means/sds
p_anova.test(50, 3, means=c(0,0,1), sds=c(1,2,1))


  # compare simulated results to pwr package
  pwr::pwr.anova.test(f=0.28, k=4, n=20)
  p_anova.test(n=20, k=4, f=.28) |> Spower()

p-value from chi-squared test simulation

Description

Generates multinomial data suitable for analysis with chisq.test.

Usage

p_chisq.test(
  n,
  w,
  df,
  correct = TRUE,
  P0 = NULL,
  P = NULL,
  gen_fun = gen_chisq.test,
  ...
)

gen_chisq.test(n, P, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_chisq.test

Examples

# effect size w + df
p_chisq.test(100, w=.2, df=3)

# vector of explicit probabilities (goodness of fit test)
p_chisq.test(100, P0 = c(.25, .25, .25, .25),
                   P = c(.6, .2, .1, .1))

# matrix of explicit probabilities (two-dimensional test of independence)
p_chisq.test(100, P0 = matrix(c(.25, .25, .25, .25), 2, 2),
                   P = matrix(c(.6, .2, .1, .1),2,2))


    # compare simulated results to pwr package

    P0 <- c(1/3, 1/3, 1/3)
    P <- c(.5, .25, .25)
    w <- pwr::ES.w1(P0, P)
    df <- 3-1
    pwr::pwr.chisq.test(w=w, df=df, N=100, sig.level=0.05)

    # slightly less power when evaluated empirically
    p_chisq.test(n=100, w=w, df=df) |> Spower(replications=100000)
    p_chisq.test(n=100, P0=P0, P=P) |> Spower(replications=100000)

    # slightly differ (latter more conservative due to finite sampling behaviour)
    pwr::pwr.chisq.test(w=w, df=df, power=.8, sig.level=0.05)
    p_chisq.test(n=NA, w=w, df=df) |>
           Spower(power=.80, interval=c(50, 200))
    p_chisq.test(n=NA, w=w, df=df, correct=FALSE) |>
           Spower(power=.80, interval=c(50, 200))

    # Spower slightly more conservative even with larger N
    pwr::pwr.chisq.test(w=.1, df=df, power=.95, sig.level=0.05)
    p_chisq.test(n=NA, w=.1, df=df) |>
           Spower(power=.95, interval=c(1000, 2000))
    p_chisq.test(n=NA, w=.1, df=df, correct=FALSE) |>
           Spower(power=.95, interval=c(1000, 2000))

p-value from (generalized) linear regression model simulations with fixed predictors

Description

p-values associated with (generalized) linear regression model. Requires a pre-specified design matrix (X).

Usage

p_glm(
  formula,
  X,
  betas,
  test,
  sigma = NULL,
  family = gaussian(),
  gen_fun = gen_glm,
  ...
)

gen_glm(formula, X, betas, sigma = NULL, family = gaussian(), ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

p_lm.R2

Examples

X <- data.frame(G = factor(rep(c('control', 'treatment'), each=50)),
                C = sample(50:100, 100, replace=TRUE))
head(X)

# ANCOVA setup
p_glm(y ~ G + C, test="Gtreatment = 0",
  X=X, betas=c(10, .3, 1), sigma=1)

# ANCOVA setup with logistic regression
p_glm(y ~ G + C, test="Gtreatment = 0",
  X=X, betas=c(-2, .5, .01), family=binomial())

# ANCOVA setup with poisson regression
p_glm(y ~ G + C, test="Gtreatment = 0",
  X=X, betas=c(-2, .5, .01), family=poisson())



# test whether two slopes differ given different samples.
#   To do this setup data as an MLR where a binary variable S
#   is used to reflect the second sample, and the interaction
#   effect evaluates the magnitude of the slope difference
gen_twogroup <- function(n, dbeta, sdx1, sdx2, sigma, n2_n1 = 1, ...){
  X1 <- rnorm(n, sd=sdx1)
  X2 <- rnorm(n*n2_n1, sd=sdx2)
  X <- c(X1, X2)
  N <- length(X)
  S <- c(rep(0, n), rep(1, N-n))
  y <- dbeta * X*S + rnorm(N, sd=sigma)
  dat <- data.frame(y, X, S)
  dat
}

# prospective power using test that interaction effect is equal to 0
p_glm(formula=y~X*S, test="X:S = 0",
	  n=100, sdx1=1, sdx2=2, dbeta=0.2,
	  sigma=0.5, gen_fun=gen_twogroup) |> Spower(replications=1000)

p-value from Kruskal-Wallis Rank Sum Test simulation

Description

Simulates data given two or more parent distributions and returns a p-value using kruskal.test. Default generates data from Gaussian distributions, however this can be modified.

Usage

p_kruskal.test(
  n,
  k,
  means,
  n.ratios = rep(1, k),
  gen_fun = gen_kruskal.test,
  ...
)

gen_kruskal.test(n, k, n.ratios, means, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

# three group test where data generate from Gaussian distributions
p_kruskal.test(n=30, k=3, means=c(0, .5, .6))

# generate data from chi-squared distributions with different variances
gen_chisq <- function(n, k, n.ratios, means, dfs, ...){
  dat <- vector('list', k)
  ns <- n * n.ratios
  for(g in 1:k)
 	dat[[g]] <- rchisq(ns[g], df=dfs[g]) - dfs[g] + means[g]
  dat
}

p_kruskal.test(n=30, k=3, means=c(0, 1, 2),
   gen_fun=gen_chisq, dfs=c(10, 15, 20))


  # empirical power estimate
  p_kruskal.test(n=30, k=3, means=c(0, .5, .6)) |> Spower()
  p_kruskal.test(n=30, k=3, means=c(0, 1, 2), gen_fun=gen_chisq,
         dfs = c(10, 15, 20)) |> Spower()

p-value from Kolmogorov-Smirnov one- or two-sample simulation

Description

Generates one or two sets of continuous data group-level data and returns a p-value under the null that the groups were drawn from the same distribution (two sample) or from a theoretically known distribution (one sample).

Usage

p_ks.test(n, p1, p2, n2_n1 = 1, two.tailed = TRUE, parent = NULL, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_t.test

Examples

# two-sample test from two Gaussian distributions with different locations
p1 <- function(n) rnorm(n)
p2 <- function(n) rnorm(n, mean=-.5)
p_ks.test(n=100, p1, p2)

# one-sample data from chi-squared distribution tested
#   against a standard normal distribution
pc <- function(n, df=15) (rchisq(n, df=df) - df) / sqrt(2*df)
p_ks.test(n=100, p1=pc, parent=pnorm, mean=0, sd=1)


  # empirical power estimates
  p_ks.test(n=100, p1, p2) |> Spower()
  p_ks.test(n=100, p1=pc, parent=pnorm, mean=0, sd=1) |> Spower()

p-value from global linear regression model simulation

Description

p-values associated with linear regression model using fixed/random independent variables. Focus is on the omnibus behavior of the R^2 statistic.

Usage

p_lm.R2(n, R2, k, R2_0 = 0, k.R2_0 = 0, R2.resid = 1 - R2, fixed = TRUE, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

p_glm

Examples

# 5 fixed IVs, R^2 = .1, sample size of 95
p_lm.R2(n=95, R2=.1, k=5)

# random model
p_lm.R2(n=95, R2=.1, k=5, fixed=FALSE)

p-value from Mauchly's Test of Sphericity simulation

Description

Perform simulation experiment for Mauchly's Test of Sphericity using the function mauchlys.test, returning a p-value. Assumes the data are from a multivariate normal distribution, however this can be modified.

Usage

p_mauchly.test(n, sigma, gen_fun = gen_mauchly.test, ...)

gen_mauchly.test(n, sigma, ...)

mauchlys.test(X)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

sigma <- diag(c(1,2,1))
sigma

p_mauchly.test(100, sigma=sigma)

# Null is true
sigma.H0 <- diag(3)
p_mauchly.test(100, sigma=sigma.H0)



    # empirical power estimate
    p_mauchly.test(100, sigma=sigma) |> Spower()

    # empirical Type I error estimate
    p_mauchly.test(100, sigma=sigma.H0) |> Spower()

p-value from McNemar test simulation

Description

Generates two-dimensional sample data for McNemar test and return a p-value. Uses mcnemar.test.

Usage

p_mcnemar.test(
  n,
  prop,
  two.tailed = TRUE,
  correct = TRUE,
  gen_fun = gen_mcnemar.test,
  ...
)

gen_mcnemar.test(n, prop, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_mcnemar.test

Examples

# from ?mcnemar.test
Performance <- matrix(c(794, 86, 150, 570),
		   nrow = 2,
		   dimnames = list("1st Survey" = c("Approve", "Disapprove"),
		               "2nd Survey" = c("Approve", "Disapprove")))
(prop <- prop.table(Performance))

# one sample + test and resulting p-value
p_mcnemar.test(n=sum(Performance), prop=prop)



# post-hoc power (not recommended)
Spower(p_mcnemar.test(n=sum(Performance), prop=prop))

p-value from three-variable mediation analysis simulation

Description

Simple 3-variable mediation analysis simulation to test the hypothesis that X -> Y is mediated by the relationship X -> M -> Y. Currently, M and Y are assumed to be continuous variables with Gaussian errors, while X may be continuous or dichotomous.

Usage

p_mediation(
  n,
  a,
  b,
  cprime,
  dichotomous.X = FALSE,
  two.tailed = TRUE,
  method = "wald",
  sd.X = 1,
  sd.Y = 1,
  sd.M = 1,
  gen_fun = gen_mediation,
  ...
)

gen_mediation(
  n,
  a,
  b,
  cprime,
  dichotomous.X = FALSE,
  sd.X = 1,
  sd.Y = 1,
  sd.M = 1,
  ...
)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_mediation

Examples

# joint test H0: a*b = 0
p_mediation(50, a=sqrt(.35), b=sqrt(.35), cprime=.39)
p_mediation(50, a=sqrt(.35), b=sqrt(.35), cprime=.39, dichotomous.X=TRUE)



  # power to detect mediation
  p_mediation(n=50, a=sqrt(.35), b=sqrt(.35), cprime=.39) |>
    Spower(parallel=TRUE, replications=1000)

  # sample size estimate for .95 power
  p_mediation(n=NA, a=sqrt(.35), b=sqrt(.35), cprime=.39) |>
    Spower(power=.95, interval=c(50, 200), parallel=TRUE)

p-value from proportion test simulation

Description

Generates single and multi-sample data for proportion tests and return a p-value. Uses binom.test for one-sample applications and prop.test otherwise.

Usage

p_prop.test(
  n,
  h,
  prop = NULL,
  pi = 0.5,
  n.ratios = rep(1, length(prop)),
  two.tailed = TRUE,
  correct = TRUE,
  exact = FALSE,
  gen_fun = gen_prop.test,
  ...
)

gen_prop.test(
  n,
  h,
  prop = NULL,
  pi = 0.5,
  n.ratios = rep(1, length(prop)),
  ...
)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_prop.test

Examples

# one sample, 50 observations, tested against pi = .5 by default
p_prop.test(50, prop=.65)

# specified using h and pi
h <- pwr::ES.h(.65, .4)
p_prop.test(50, h=h, pi=.4)
p_prop.test(50, h=-h, pi=.65)

# two-sample test
p_prop.test(50, prop=c(.5, .65))

# two-sample test, unequal ns
p_prop.test(50, prop=c(.5, .65), n.ratios = c(1,2))

# three-sample test, group2 twice as large as others
p_prop.test(50, prop=c(.5, .65, .7), n.ratios=c(1,2,1))

# Fisher exact test
p_prop.test(50, prop=matrix(c(.5, .65, .7, .5), 2, 2))


    # compare simulated results to pwr package

    # one-sample tests
    (h <- pwr::ES.h(0.5, 0.4))
    pwr::pwr.p.test(h=h, n=60)

    # uses binom.test (need to specify null location as this matters!)
    Spower(p_prop.test(n=60, h=h, pi=.4))
    Spower(p_prop.test(n=60, prop=.5, pi=.4))

    # compare with switched null
    Spower(p_prop.test(n=60, h=h, pi=.5))
    Spower(p_prop.test(n=60, prop=.4, pi=.5))

    # two-sample test, one-tailed
    (h <- pwr::ES.h(0.67, 0.5))
    pwr::pwr.2p.test(h=h, n=80, alternative="greater")
    p_prop.test(n=80, prop=c(.67, .5), two.tailed=FALSE,
      correct=FALSE) |> Spower()

    # same as above, but with continuity correction (default)
    p_prop.test(n=80, prop=c(.67, .5), two.tailed=FALSE) |>
      Spower()

    # three-sample joint test, equal n's
    p_prop.test(n=50, prop=c(.6,.4,.7)) |> Spower()

p-value from correlation simulation

Description

Generates correlated X-Y data and returns a p-value to assess the null of no correlation in the population. The X-Y data are generated assuming a bivariate normal distribution.

Usage

p_r(n, r, rho = 0, method = "pearson", two.tailed = TRUE, gen_fun = gen_r, ...)

gen_r(n, r, ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_r

Examples

# 50 observations, .5 correlation
p_r(50, r=.5)
p_r(50, r=.5, method = 'spearman')

# test against constant other than rho = .6
p_r(50, .5, rho=.60)



    # compare simulated results to pwr package

    pwr::pwr.r.test(r=0.3, n=50)
    p_r(n=50, r=0.3) |> Spower()

    pwr::pwr.r.test(r=0.3, power=0.80)
    p_r(n=NA, r=0.3) |> Spower(power=.80, interval=c(10, 200))

    pwr::pwr.r.test(r=0.1, power=0.80)
    p_r(n=NA, r=0.1) |> Spower(power=.80, interval=c(200, 1000))

p-value from tetrachoric/polychoric or polyserial

Description

Generates correlated X-Y data and returns a p-value to assess the null of no correlation in the population. The X-Y data are generated assuming a multivariate normal distribution and subsequently discretized for one or both of the variables.

Usage

p_r.cat(
  n,
  r,
  tauX,
  rho = 0,
  tauY = NULL,
  ML = TRUE,
  two.tailed = TRUE,
  score = FALSE,
  gen_fun = gen_r,
  ...
)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_r

Examples

# 100 observations, .5 correlation, tetrachoric estimate
p_r.cat(100, r=.5, tauX=0, tauY=1)

# Wald test
p_r.cat(100, r=.5, tauX=0, tauY=1, score=FALSE)

# polyserial estimate (Y continuous)
p_r.cat(50, r=.5, tauX=0)

p-value from Scale Test simulation

Description

Simulates data given one or two parent distributions and returns a p-value testing that the scale of the type distributions are the same. Default implementation uses Gaussian distributions, however the distribution function may be modified to reflect other populations of interest. Uses ansari.test or mood.test for the analysis.

Usage

p_scale(
  n,
  scale,
  n2_n1 = 1,
  two.tailed = TRUE,
  exact = NULL,
  test = "Ansari",
  parent = function(n, ...) rnorm(n),
  ...
)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

# n=30 per group,
#  Distributions Gaussian with sd=1 for first group and sd=2 for second
p_scale(30, scale=2)
p_scale(30, scale=2, test='Mood')

# compare chi-squared distributions
parent <- function(n, df, ...) rchisq(n, df=df) - df
p_scale(30, scale=2, parent=parent, df=3)


  # empirical power of the experiments
  p_scale(30, scale=2) |> Spower()
  p_scale(30, scale=2, test='Mood') |> Spower()

  p_scale(30, scale=2, parent=parent, df=3) |> Spower()
  p_scale(30, scale=2, test='Mood', parent=parent, df=3) |> Spower()

p-value from Shapiro-Wilk Normality Test simulation

Description

Generates univariate distributional data and returns a p-value to assess the null that the population follows a Gaussian distribution shape. Uses shapiro.test.

Usage

p_shapiro.test(dist)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

# 50 observations drawn from normal distribution (null is true)
p_shapiro.test(rnorm(50))

# 50 observations from slightly skewed chi-squared distribution (power)
p_shapiro.test(rchisq(50, df=100))


    # empirical Type I error rate estimate
    p_shapiro.test(rnorm(50)) |> Spower()

    # power
    p_shapiro.test(rchisq(50, df=100)) |> Spower()

p-value from independent/paired samples t-test simulation

Description

Generates one or two sets of continuous data group-level data according to Cohen's effect size 'd', and returns a p-value. The data and associated t-test assume that the conditional observations are normally distributed and have have equal variance by default, however these may be modified.

Usage

p_t.test(
  n,
  d,
  mu = 0,
  r = NULL,
  type = c("two.sample", "one.sample", "paired"),
  n2_n1 = 1,
  two.tailed = TRUE,
  var.equal = TRUE,
  means = NULL,
  sds = NULL,
  gen_fun = gen_t.test,
  ...
)

gen_t.test(
  n,
  d,
  n2_n1 = 1,
  r = NULL,
  type = c("two.sample", "one.sample", "paired"),
  means = NULL,
  sds = NULL,
  ...
)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_t.test

Examples

# sample size of 50 per group, "medium" effect size
p_t.test(n=50, d=0.5)

# point-biserial correlation effect size
p_t.test(n=50, r=.3)

# second group 2x as large as the first group
p_t.test(n=50, d=0.5, n2_n1 = 2)

# specify mean/SDs explicitly
p_t.test(n=50, means = c(0,1), sds = c(2,2))

# paired and one-sample tests
p_t.test(n=50, d=0.5, type = 'paired')
p_t.test(n=50, d=0.5, type = 'one.sample')


  # compare simulated results to pwr package

  pwr::pwr.t.test(d=0.2, n=60, sig.level=0.10,
             type="one.sample", alternative="two.sided")
  p_t.test(n=60, d=0.2, type = 'one.sample', two.tailed=TRUE) |>
         Spower(sig.level=.10)

  pwr::pwr.t.test(d=0.3, power=0.80, type="two.sample",
                  alternative="greater")
  p_t.test(n=NA, d=0.3, type='two.sample', two.tailed=FALSE) |>
         Spower(power=0.80, interval=c(10,200))




###### Custom data generation function

# Generate data such that:
#   - group 1 is from a negatively distribution (reversed X2(10)),
#   - group 2 is from a positively skewed distribution (X2(5))
#   - groups have equal variance, but differ by d = 0.5

args(gen_t.test)   ## can use these arguments as a basis, though must include ...

# arguments df1 and df2 added; unused arguments caught within ...
my.gen_fun <- function(n, d, df1, df2, ...){
 	 group1 <- -1 * rchisq(n, df=df1)
	     group2 <- rchisq(n, df=df2)
	     # scale groups first given moments of the chi-square distribution,
	     #   then add std mean difference
	     group1 <- ((group1 + df1) / sqrt(2*df1))
	     group2 <- ((group2 - df2) / sqrt(2*df2)) + d
	     dat <- data.frame(DV=c(group1, group2),
	        			   group=gl(2, n, labels=c('G1', 'G2')))
	     dat
}

# check the sample data properties
df <- my.gen_fun(n=10000, d=.5, df1=10, df2=5)
with(df, tapply(DV, group, mean))
with(df, tapply(DV, group, sd))

library(ggplot2)
ggplot(df, aes(group, DV, fill=group)) + geom_violin()

p_t.test(n=100, d=0.5, gen_fun=my.gen_fun, df1=10, df2=5)



  # power given Gaussian distributions
  p_t.test(n=100, d=0.5) |> Spower(replications=30000)

  # estimate power given the customized data generating function
  p_t.test(n=100, d=0.5, gen_fun=my.gen_fun, df1=10, df2=5) |>
    Spower(replications=30000)

  # evaluate Type I error rate to see if liberal/conservative given
  # assumption violations (should be close to alpha/sig.level)
  p_t.test(n=100, d=0, gen_fun=my.gen_fun, df1=10, df2=5) |>
    Spower(replications=30000)

p-value from variance test simulation

Description

Generates one or or more sets of continuous data group-level data to perform a variance test, and return a p-value. When two-samples are investigated the var.test function will be used, otherwise functions from the EnvStats package will be used.

Usage

p_var.test(
  n,
  vars,
  n.ratios = rep(1, length(vars)),
  sigma2 = 1,
  two.tailed = TRUE,
  test = "Levene",
  correct = TRUE,
  gen_fun = gen_var.test,
  ...
)

gen_var.test(n, vars, n.ratios = rep(1, length(vars)), ...)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

See Also

gen_var.test

Examples

# one sample
p_var.test(100, vars=10, sigma2=9)

# three sample
p_var.test(100, vars=c(10, 9, 11))
p_var.test(100, vars=c(10, 9, 11), test = 'Fligner')
p_var.test(100, vars=c(10, 9, 11), test = 'Bartlett')


  # power to detect three-group variance differences
  p_var.test(n=100, vars=c(10,9,11)) |> Spower()

  # sample size per group to achieve 80% power
  p_var.test(n=NA, vars=c(10,9,11)) |>
         Spower(power=.80, interval=c(100, 1000))

p-value from Wilcox test simulation

Description

Simulates data given one or two parent distributions and returns a p-value. Can also be used for power analyses related to sign tests.

Usage

p_wilcox.test(
  n,
  d,
  n2_n1 = 1,
  mu = 0,
  type = c("two.sample", "one.sample", "paired"),
  exact = NULL,
  correct = TRUE,
  two.tailed = TRUE,
  parent1 = function(n, d) rnorm(n, d, 1),
  parent2 = function(n, d) rnorm(n, 0, 1)
)

Arguments

Value

a single p-value

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

# with normal distributions defaults d is standardized
p_wilcox.test(100, .5)
p_wilcox.test(100, .5, type = 'paired')
p_wilcox.test(100, .5, type = 'one.sample')

# using chi-squared distributions (standardizing to 0-1)
p_wilcox.test(100, .5, type = 'one.sample',
   parent1 = function(n, d) rchisq(n, df=10) - 10 + d)
p_wilcox.test(100, .5,
   parent1 = function(n, d) (rchisq(n, df=10) - 10)/sqrt(20) + d,
   parent2 = function(n, d) (rchisq(n, df=10) - 10)/sqrt(20))

Update compromise or prospective/post-hoc power analysis without re-simulating

Description

When a power or compromise analysis was performed in Spower this function can be used to update the compromise or power criteria without the need for re-simulating the experiment. For compromise analyses a beta_alpha criteria must be supplied, while for prospective/post-hoc power analyses the sig.level must be supplied.

Usage

## S3 method for class 'Spower'
update(object, sig.level = 0.05, beta_alpha = NULL, predCI = 0.95, ...)

Arguments

Value

object of class Spower with updated information

Author(s)

Phil Chalmers rphilip.chalmers@gmail.com

Examples

########
## Prospective power analysis update

# Estimate power using sig.level = .05 (default)
out <- p_t.test(n = 50, d = .5) |> Spower()

# update power estimate given sig.level=.01 and .20
update(out, sig.level=.01)
update(out, sig.level=.20)


########
## Compromise analysis update

# Solve beta/alpha ratio to specific error trade-off constant
out <- p_t.test(n = 50, d = .5) |> Spower(beta_alpha = 2)

# update beta_alpha criteria without re-simulating
update(out, beta_alpha=4)

# also works if compromise not initially run but prospective/post-hoc power was
out <- p_t.test(n = 50, d = .5) |> Spower()
update(out, beta_alpha=4)