Power-curves are useful when researchers are interested in the range of power effects when varying different population and sample properties. For instance, researchers may be interested in how power improves as the magnitude of the population mean difference for an independent samples \(t\)-test increases, or how quickly power improves when increasing the sample size (\(N\)). The following simulation assumes that the population variances are equal.
Define the functionslibrary(SimDesign)
#SimFunctions(comments = FALSE)
### Define design conditions. Here we want a finer range of values for plotting
Design <- createDesign(mean_diff = seq(0, 1, by = .1),
sample_size = c(10, 20, 30))
#--------------------------------------------------------------------
Generate <- function(condition, fixed_objects = NULL) {
Attach(condition)
dat1 <- rnorm(sample_size)
dat2 <- rnorm(sample_size, mean=mean_diff)
dat <- list(dat1=dat1, dat2=dat2)
dat
}
Analyse <- function(condition, dat, fixed_objects = NULL) {
ret <- c(p = t.test(dat$dat1, dat$dat2,
var.equal=TRUE)$p.value)
ret
}
Summarise <- function(condition, results, fixed_objects = NULL) {
ret <- EDR(results, alpha = .05)
ret
}
#--------------------------------------------------------------------
### Run the simulation (using all available cores)
res <- runSimulation(Design, replications=10000, verbose=FALSE, parallel=TRUE,
generate=Generate, analyse=Analyse, summarise=Summarise)# summary of simulation object
summary(res)$sessionInfo
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.3.3 (2024-02-29)
os Ubuntu 20.04.6 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Etc/UTC
date 2025-04-17
pandoc 2.9.1.1 @ /usr/bin/ (via rmarkdown)
quarto 1.7.25 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
audio 0.1-11 2023-08-18 [1] CRAN (R 4.3.1)
beepr 2.0 2024-07-06 [1] CRAN (R 4.3.3)
brio 1.1.5 2024-04-24 [1] CRAN (R 4.3.3)
cli 3.6.4 2025-02-13 [1] CRAN (R 4.3.3)
codetools 0.2-20 2024-03-31 [4] CRAN (R 4.3.3)
digest 0.6.37 2024-08-19 [1] CRAN (R 4.3.3)
dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.3.1)
evaluate 1.0.3 2025-01-10 [1] CRAN (R 4.3.3)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.3.3)
future 1.40.0 2025-04-10 [1] CRAN (R 4.3.3)
future.apply 1.11.3 2024-10-27 [1] CRAN (R 4.3.3)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.1)
globals 0.17.0 2025-04-16 [1] CRAN (R 4.3.3)
glue 1.8.0 2024-09-30 [1] CRAN (R 4.3.3)
htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.3.3)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.3.3)
knitr 1.50 2025-03-16 [1] CRAN (R 4.3.3)
lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.1)
listenv 0.9.1 2024-01-29 [1] CRAN (R 4.3.2)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.1)
parallelly 1.43.0 2025-03-24 [1] CRAN (R 4.3.3)
pbapply 1.7-2 2023-06-27 [1] CRAN (R 4.3.1)
pillar 1.10.2 2025-04-05 [1] CRAN (R 4.3.3)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.1)
progressr 0.15.1 2024-11-22 [1] CRAN (R 4.3.3)
R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.2)
R.oo 1.27.0 2024-11-01 [1] CRAN (R 4.3.3)
R.utils 2.13.0 2025-02-24 [1] CRAN (R 4.3.3)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.3)
rlang 1.1.6 2025-04-11 [1] CRAN (R 4.3.3)
rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.3.3)
sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.3.3)
SimDesign * 2.19.3 2025-04-16 [1] local
testthat 3.2.3 2025-01-13 [1] CRAN (R 4.3.3)
tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.1)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.3)
vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.1)
xfun 0.52 2025-04-02 [1] CRAN (R 4.3.3)
yaml 2.3.10 2024-07-26 [1] CRAN (R 4.3.3)
[1] /home/phil/R/x86_64-pc-linux-gnu-library/4.3
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────
$save_info
named character(0)
$seeds
[1] 1140350788 312928385 866248189 1909893419 554504146 884616499
[7] 803234389 1158971242 934673902 1632225031 1867003471 651436069
[13] 1147691737 57500133 1587828514 2070518142 2119708588 1406053153
[19] 1263737763 960850250 1378461094 1215204756 127650324 100235820
[25] 997765463 2051398777 1881620140 303569271 1358398662 695822375
[31] 1773292330 1971672344 647971342
$ncores
[1] 47
$number_of_conditions
[1] 33
$date_completed
[1] Thu Apr 17 21:30:12 2025
$total_elapsed_time
[1] 15.00s
$Design.ID
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33# print results
print(res)# A tibble: 33 × 7
mean_diff sample_size p REPLICATIONS SIM_TIME SEED COMPLETED
<dbl> <dbl> <dbl> <dbl> <chr> <int> <chr>
1 0 10 0.0495 10000 0.56s 1140350788 Thu Apr 17 21:…
2 0.1 10 0.0541 10000 0.37s 312928385 Thu Apr 17 21:…
3 0.2 10 0.0698 10000 0.34s 866248189 Thu Apr 17 21:…
4 0.3 10 0.1002 10000 0.46s 1909893419 Thu Apr 17 21:…
5 0.4 10 0.1365 10000 0.61s 554504146 Thu Apr 17 21:…
6 0.5 10 0.1847 10000 0.41s 884616499 Thu Apr 17 21:…
7 0.6 10 0.2445 10000 0.52s 803234389 Thu Apr 17 21:…
8 0.7 10 0.3149 10000 0.45s 1158971242 Thu Apr 17 21:…
9 0.8 10 0.3864 10000 0.57s 934673902 Thu Apr 17 21:…
10 0.9 10 0.4686 10000 0.51s 1632225031 Thu Apr 17 21:…
# ℹ 23 more rowsA power curve is created by placing the detection rates on a \(y\)-axis and including different factors on the x-axis and other aesthetics (shading, colours, etc). Here we use ggplot2 to construct suitable power-curves.
library(ggplot2)
ggplot(res, aes(mean_diff, p, colour=factor(sample_size))) + geom_point() + geom_line() +
xlab('Mean Difference') + ylab('Detection Rate') + ggtitle('Empirical Power Curves') +
scale_color_discrete('N')Compare these results to analytic formulas, we can see that the two are essentially coinciding with good accuracy.
library(pwr)
pwr.t.test(d=0.5, n=30, sig.level=0.05, type="two.sample", alternative="two.sided")
Two-sample t test power calculation
n = 30
d = 0.5
sig.level = 0.05
power = 0.478
alternative = two.sided
NOTE: n is number in *each* groupsubset(res, mean_diff == .5 & sample_size == 30)# A tibble: 1 × 7
mean_diff sample_size p REPLICATIONS SIM_TIME SEED COMPLETED
<dbl> <dbl> <dbl> <dbl> <chr> <int> <chr>
1 0.5 30 0.4827 10000 0.47s 303569271 Thu Apr 17 21:30…pwr.t.test(d=0.5, n=10, sig.level=0.05, type="two.sample", alternative="two.sided")
Two-sample t test power calculation
n = 10
d = 0.5
sig.level = 0.05
power = 0.185
alternative = two.sided
NOTE: n is number in *each* groupsubset(res, mean_diff == .5 & sample_size == 10)# A tibble: 1 × 7
mean_diff sample_size p REPLICATIONS SIM_TIME SEED COMPLETED
<dbl> <dbl> <dbl> <dbl> <chr> <int> <chr>
1 0.5 10 0.1847 10000 0.41s 884616499 Thu Apr 17 21:29…In situations where analytic power estimation is not available, and therefore simulations are used as a suitable stand-in, it’s important to evaluate the precision of the observed power estimates from the simulation study. As such, one reasonable (and extremely general approach) is to apply bootstrapping methodology to the within replication; for simulated power analyses this involves bootstrapping the \(R\) replicated \(p\)-values to obtain suitable non-parametric confidence intervals.
The following implements the empirical (i.e., basic) bootstrap CI method for the simulation above, and instead plots the point-wise 95% CIs estimates for each condition in the Design input. Notice that the replications were reduced from 10000 to 1000 to demonstrate the potential uncertainty when using too few replications when estimating power via simulation.
### Run the simulation (using all available cores)
res <- runSimulation(Design, replications=1000, verbose=FALSE, parallel=TRUE,
generate=Generate, analyse=Analyse, summarise=Summarise,
boot_method = 'basic')
res# A tibble: 33 × 9
mean_diff sample_size p BOOT_p_2.5 BOOT_p_97.5 REPLICATIONS SIM_TIME
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 0 10 0.052 0.038 0.065 1000 0.78s
2 0.1 10 0.064 0.048 0.079 1000 0.68s
3 0.2 10 0.075 0.058 0.091 1000 0.69s
4 0.3 10 0.09 0.072 0.10702 1000 0.66s
5 0.4 10 0.129 0.109 0.15 1000 0.73s
6 0.5 10 0.206 0.181 0.229 1000 0.69s
7 0.6 10 0.252 0.226 0.27802 1000 0.70s
8 0.7 10 0.299 0.27097 0.327 1000 0.69s
9 0.8 10 0.395 0.367 0.424 1000 0.70s
10 0.9 10 0.489 0.458 0.52002 1000 0.66s
# ℹ 23 more rows
# ℹ 2 more variables: SEED <int>, COMPLETED <chr>library(ggplot2)
ggplot(res, aes(mean_diff, p, colour=factor(sample_size))) + geom_point() + geom_line() +
geom_ribbon(aes(ymin=BOOT_p_2.5, ymax=BOOT_p_97.5, fill=factor(sample_size)), alpha=.3, linetype=0) +
xlab('Mean Difference') + ylab('Detection Rate') + ggtitle('Empirical Power Curves with 95% CIs') +
scale_fill_discrete('N') + scale_colour_discrete('N')Finally, when increasing the replications back to 10000 it become visually clear how precise the original point-estimates presented in the previous section were.
### Run the simulation (using all available cores)
res <- runSimulation(Design, replications=10000, verbose=FALSE, parallel=TRUE,
generate=Generate, analyse=Analyse, summarise=Summarise,
boot_method = 'basic')
res# A tibble: 33 × 9
mean_diff sample_size p BOOT_p_2.5 BOOT_p_97.5 REPLICATIONS SIM_TIME
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 0 10 0.0528 0.0484 0.057102 10000 6.43s
2 0.1 10 0.053 0.0486 0.0574 10000 6.35s
3 0.2 10 0.0737 0.068498 0.078902 10000 6.29s
4 0.3 10 0.0962 0.090497 0.10170 10000 6.36s
5 0.4 10 0.1328 0.1261 0.1394 10000 6.15s
6 0.5 10 0.1898 0.18250 0.1977 10000 6.07s
7 0.6 10 0.2539 0.24520 0.26240 10000 6.10s
8 0.7 10 0.3249 0.3159 0.33450 10000 6.36s
9 0.8 10 0.3959 0.38550 0.4059 10000 6.11s
10 0.9 10 0.475 0.46530 0.4848 10000 6.45s
# ℹ 23 more rows
# ℹ 2 more variables: SEED <int>, COMPLETED <chr>library(ggplot2)
ggplot(res, aes(mean_diff, p, colour=factor(sample_size))) + geom_point() + geom_line() +
geom_ribbon(aes(ymin=BOOT_p_2.5, ymax=BOOT_p_97.5, fill=factor(sample_size)), alpha=.3, linetype=0) +
xlab('Mean Difference') + ylab('Detection Rate') + ggtitle('Empirical Power Curves with 95% CIs') +
scale_fill_discrete('N') + scale_colour_discrete('N')Power-curves can also be useful for obtaining estimates to satisfy particular power rates when sample size planning. Note that this is a rather inefficient approach to solving sample sizes, though it has been very popular (see ?SimDesign::SimSolve for a stochastic root-solving alternative).
In the following, a \(d=.5\) effect size is assumed in the population, where the value of \(N\) associated with a target power rate of \(.80\) is desired. To obtain an \(\hat{N}\) that is roughly associated with this target power a range of \(N\) inputs are first evaluated, the adjacent points closest to .80 are isolated, and finally the \(\hat{N}\) is solved using linear interpolations to obtain an approximate target sample size.
Design <- createDesign(mean_diff = .5, # equivalent to Cohen's "medium" d effect size
sample_size = seq(10, 100, by=5))
res <- runSimulation(Design, replications=10000, verbose=FALSE, parallel=TRUE,
generate=Generate, analyse=Analyse, summarise=Summarise)
# target N somewhere between 60 and 65
res# A tibble: 19 × 7
mean_diff sample_size p REPLICATIONS SIM_TIME SEED COMPLETED
<dbl> <dbl> <dbl> <dbl> <chr> <int> <chr>
1 0.5 10 0.1836 10000 0.85s 540368552 Thu Apr 17 21:…
2 0.5 15 0.2715 10000 0.51s 895502048 Thu Apr 17 21:…
3 0.5 20 0.3334 10000 0.49s 1200947712 Thu Apr 17 21:…
4 0.5 25 0.4101 10000 0.47s 1390946546 Thu Apr 17 21:…
5 0.5 30 0.4776 10000 0.36s 297095334 Thu Apr 17 21:…
6 0.5 35 0.5336 10000 0.39s 26408929 Thu Apr 17 21:…
7 0.5 40 0.6075 10000 0.59s 953293844 Thu Apr 17 21:…
8 0.5 45 0.6491 10000 0.47s 896131004 Thu Apr 17 21:…
9 0.5 50 0.6958 10000 0.42s 1166181780 Thu Apr 17 21:…
10 0.5 55 0.7486 10000 0.50s 205322161 Thu Apr 17 21:…
11 0.5 60 0.7749 10000 0.50s 1570516689 Thu Apr 17 21:…
12 0.5 65 0.8085 10000 0.36s 1818143070 Thu Apr 17 21:…
13 0.5 70 0.8371 10000 0.55s 1103983110 Thu Apr 17 21:…
14 0.5 75 0.8573 10000 0.40s 627018390 Thu Apr 17 21:…
15 0.5 80 0.8755 10000 0.42s 1596415590 Thu Apr 17 21:…
16 0.5 85 0.8951 10000 0.49s 1728962223 Thu Apr 17 21:…
17 0.5 90 0.9171 10000 0.42s 1920077889 Thu Apr 17 21:…
18 0.5 95 0.927 10000 0.52s 1388819422 Thu Apr 17 21:…
19 0.5 100 0.9442 10000 0.46s 828181090 Thu Apr 17 21:…# solve using linear interpolated power
pmin <- min(which(res$p > .80))
pick <- res[c(pmin-1, pmin),]
pick# A tibble: 2 × 7
mean_diff sample_size p REPLICATIONS SIM_TIME SEED COMPLETED
<dbl> <dbl> <dbl> <dbl> <chr> <int> <chr>
1 0.5 60 0.7749 10000 0.50s 1570516689 Thu Apr 17 21:3…
2 0.5 65 0.8085 10000 0.36s 1818143070 Thu Apr 17 21:3…After the empirical power curve has been built, and it’s clear the target power is within the range specified when building the power curve, interpolations can be used to solve for \(\hat{N}\). In this case a linear interpolation is used between the closest points, though alternatives are of course possible (e.g., fitting higher-order polynomials, splines, or other non-linear models).
mod <- lm(p ~ sample_size, pick)
coef(mod)(Intercept) sample_size
0.37170 0.00672 fn <- function(N, power = .8)
predict(mod, newdata = data.frame(sample_size=N)) - power
root <- uniroot(fn, interval = c(10, 100))
root # sample size estimate from the root solver$root
[1] 63.7
$f.root
1
-1.11e-16
$iter
[1] 2
$init.it
[1] NA
$estim.prec
[1] 6.1e-05# build power curve and include interpolation estimates
ggplot(res, aes(sample_size, p)) + geom_point() + geom_line() +
xlab('N') + ylab('Detection Rate') +
geom_abline(intercept = .80, slope=0, col='red', linetype='dashed') +
geom_vline(xintercept=root$root, col='blue', linetype='dashed') +
annotate("text", x=80, y=.6, label = paste0('N estimate = ', round(root$root, 2))) +
ggtitle('Power curve with linear interpolations')Finally, in this case one can compare the \(\hat{N}\) with the \(N\) associated with a power of .80 given the well-known non-central \(t\)-distribution using the pwr package. Had the associated non-central distribution been unavailable further Monte Carlo simulations around the estimate would be required (e.g., using runSimulation() given the integer values for sample_size given 63 and 64).
# true N from root-solving the deterministic power function
pwr.t.test(d=0.5, sig.level=0.05, power=.8, type="two.sample", alternative="two.sided")
Two-sample t test power calculation
n = 63.8
d = 0.5
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* groupRetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4