Power and Simulation

Lisa DeBruine, Cristian Mesquida, Daniël Lakens

https://scienceverse.org/talks/2026-power-sim/
https://github.com/scienceverse/talks

Abstract

Power analyses help us to plan studies and interpret findings. The statistical power of a specific test depends on the effect size, sample size, and alpha criterion. In this talk, I will cover common mistakes with conducting and reporting power analyses, and present simulation-based methods for conducting power analyses for any inferential test. I will cover fixed effects and mixed effects models, with additional materials at https://debruine.github.io/data-sim-workshops/.

Power

• false positive rate: the probability a test concludes there is an effect when there is no effect; Type I Error Rate
• alpha : the false positive rate we accept for a test
• false negative rate: the probability a test concludes there is no effect when there is one; Type II Error Rate
• beta: the false negative rate we accept for a test
• true positive rate: the probability a test concludes there is an effect when there is one
• power: the true positive rate (1-beta), given a test, effect, sample size, and alpha

Types of Power Analysis

• A Priori: calculates the sample size required to achieved a pre-specified power for a pre-specified effect size

• Sensitivity (Effect Size): calculates the smallest effect size required to detect a pre-specified sample size with a pre-specified power

• Sensitivity (Power): calculates power given a pre-specified effect size and sample size

• Post Hoc: calculates the power associated with the observed effect size and the sample size actually obtained

Reporting

• power type: a priori, power/sample size sensitivity, post hoc
• power: value between alpha and 100%
• sample size: number of subjects/items
• effect size: numeric value
• effect size metric: e.g., Cohen’s d, r, \(\eta_{p}^{2}\), raw
• statistical test: full specification
• alpha: 0-1, often assumed to be 0.05
• software: e.g., G*Power, Superpower, simulation!

How would you write up this power analysis?

pwr::pwr.t.test(          
  power = 0.80,           
  d = 0.5,                
  sig.level = 0.01,       
  type = "two.sample",    
  alternative = "greater" 
)

     Two-sample t test power calculation 

              n = 81.65515
              d = 0.5
      sig.level = 0.01
          power = 0.8
    alternative = greater

NOTE: n is number in *each* group

Model Reporting Example

We conducted an a priori analysis to determine the sample size with 80% power to detect the smallest effect size of interest (SESOI), which was a Cohen’s d of 0.5 determined by a cost-benefit analysis of the intervention. We used the pwr package (Champely, 2020) for an independent samples t-test, assuming equal variance and using a one-tailed prediction, with a critical alpha of 0.01. This determined that we would need 82 subjects in each group.

The Reality

A sample size of 22 participants was predetermined to yield a power of 80% and to detect an effect of comparable size to previous studies using similar independent and dependent measures (CITE).

The Reality

A G*Power (Version 3.1; Faul et al., 2009) analysis using pilot data from 20 participants for a Bonferroni-corrected (for two comparisons) one-sample t test (d = 0.45, α = .05, p = .95) determined that we should recruit 67 participants.

The Reality

We decided to obtain a sample size that would allow us to detect a within-subjects difference of at least a medium effect size (Cohen’s d > 0.5) with 80% power, using a two-tailed t-test. A power analysis using G*Power (Version 3.1.9.3; Faul et al., 2007) indicated that this required a sample size of 34.

We Could Do Better

Only 11% of 1295 power analyses from 1937 Psychological Science papers published between 2014 and 2026 have complete reporting.

Component	Percent Missing1
power type	0%
power	15%
sample size	19%
effect size	32%
effect size metric	28%
statistical test	22%
alpha	67%
software	71%

1. Preliminary analysis using Metacheck, from Mesquida, Lakens & DeBruine (in prep)

Power Simulation

How do you calculate power for a design where we think we will only be able to recruit twice as many younger people as older people, and expect the older subjects to have a lower score with more variance?

Simulate Data Sets

sim_data_func <- function(n, rep = 100) {
  faux::sim_design(
    between = list(age_group = c("young", "old")),
    n = c(2*n, n),
    mu = c(100, 90),
    sd = c(20, 25),
    dv = "score",
    long = TRUE,
    rep = rep
  )
}

sim_data <- sim_data_func(n = 30)

rep	age_group	n	mean	sd
1	young	60	101.0	17.4
1	old	30	89.4	26.6
2	young	60	101.6	16.9
2	old	30	83.0	26.5

Analysis Function

analysis_func <- function(data) {
  tt <- t.test(score ~ age_group, data = data, 
               alternative = "greater", 
               var.equal = FALSE)
  broom::tidy(tt)[, 1:8]
}

Test on the first data set

analysis_func(sim_data$data[[1]])

estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high
11.6	101.0	89.4	2.2	0.018	41.9	2.6	Inf

Power Sensitivity

sim_results <- sim_data |>
  mutate(analysis = purrr::map(data, analysis_func)) |>
  select(-data) |>
  unnest(analysis)

power <- mean(sim_results$p.value < .05)

A Priori

Iterate the process above over a sequence of possible Ns.

A Priori

Decrease the range and increase the density of sampling and number of replications to get a more accurate estimate.

Write Up

We conducted an a priori analysis to determine the sample size with 80% power to detect the smallest effect size of interest (SESOI), which was a 10-point raw difference in score between younger and older subjects, where younger subjects had an SD of 20 and older subject had an SD of 25. We used a custom data simulation script in R (see shared code) with 250 replications for a Welch two-sample t-test using a one-tailed prediction, with a critical alpha of 0.05. This determined that we would need approximately 51 subjects in the older group and 102 subjects in the younger group.

Mixed Effects

xkcd #2682

Mixed Effects

xkcd #2682

Faux

Shiny app or scriptable

Extensive Documentation

Simulate a Crossed Design

• 100 subjects
• 20 items
• two conditions: control and experimental
• dependent variable = reaction time in ms
• SESOI = 10 ms
• variability between subjects and between items in average RT and the effect of condition

Random factors

Start with a small number of subjects to more easily view the resulting table.

subj_n <- 2 # number of subjects

lmem_dat <- add_random(subj = subj_n)

subj
subj1
subj2

Crossed random factors

Cross subjects with a small number of items.

item_n <- 2 # number of items

lmem_dat <- add_random(subj = subj_n) |>
  add_random(item = item_n)

subj	item
subj1	item1
subj1	item2
subj2	item1
subj2	item2

Nested random factors

You can also nest items in faux (but this is not our design).

lmem_dat <- add_random(subj = subj_n) |>
  add_random(item = item_n, .nested_in = "subj")

subj	item
subj1	item1
subj1	item2
subj2	item3
subj2	item4

Fixed factors (between)

You can set fixed factors that vary between subject or items (but this is not our design).

lmem_dat <- add_random(subj = subj_n) |>
  add_random(item = item_n) |>
  add_between(condition = c("control", "experimental"), .by = "item")

subj	item	condition
subj1	item1	control
subj1	item2	experimental
subj2	item1	control
subj2	item2	experimental

Fixed factors (within)

Here, we set a factor that varies within subjects and items.

lmem_dat <- add_random(subj = subj_n) |>
  add_random(item = item_n) |>
  add_within(condition = c("control", "experimental"))

subj	item	condition
subj1	item1	control
subj1	item1	experimental
subj1	item2	control
subj1	item2	experimental
subj2	item1	control
subj2	item1	experimental
subj2	item2	control
subj2	item2	experimental

Contrast Coding

My name	Other names	R	Julia	faux
Treatment	Treatment (2), Dummy (1, 4, 6), Simple (5)	contr.treatment	DummyCoding	contr_code_treatment
Anova	Deviation (2), Contrast (1), Simple (4)	contr.treatment - 1/k	HypothesisCoding	contr_code_anova
Sum	Sum (1, 2, 6), Effects (3), Deviation (4, 5), Unweighted Effects (7)	contr.sum	EffectsCoding	contr_code_sum
Difference	Contrast (3), Forward/Backward (4), Repeated (5)	MASS::contr.sdif	SeqDiffCoding	contr_code_difference
Helmert	Reverse Helmert (1, 4), Difference (5), Contrast (7)	contr.helmert / (column_i+1)	HelmertCoding	contr_code_helmert
Polynomial	Polynomial (5), Orthogonal Polynomial (4), Trend (3)	contr.poly	HypothesisCoding	contr_code_poly

Faux Explainer | RePsychLing Explainer

Codings

lmem_dat <- add_random(subj = subj_n) |>
  add_random(item = item_n) |>
  add_within(condition = c("control", "experimental")) |>
  add_contrast("condition", "treatment", colnames = "treat") |>
  add_contrast("condition", "anova", colnames = "anova") |>
  add_contrast("condition", "sum", colnames = "sum")

subj	item	condition	treat	anova	sum
subj1	item1	control	0	-0.5	1
subj1	item1	experimental	1	0.5	-1
subj1	item2	control	0	-0.5	1
subj1	item2	experimental	1	0.5	-1
subj2	item1	control	0	-0.5	1
subj2	item1	experimental	1	0.5	-1
subj2	item2	control	0	-0.5	1
subj2	item2	experimental	1	0.5	-1

Design

# define this to make room later

design <- add_random(subj = subj_n) |>
  add_random(item = item_n) |>
  add_within(condition = c("control", "experimental")) |>
  add_contrast("condition", "treatment", colnames = "cond")

subj	item	condition	cond
subj1	item1	control	0
subj1	item1	experimental	1
subj1	item2	control	0
subj1	item2	experimental	1
subj2	item1	control	0
subj2	item1	experimental	1
subj2	item2	control	0
subj2	item2	experimental	1

Fixed Effects

intercept <- 600 # model intercept (mean for control condition)
cond_eff  <-  10 # condition effect size

lmem_dat <- design |>
  mutate(dv = intercept + 
           cond * cond_eff)

subj	item	condition	cond	dv
subj1	item1	control	0	600
subj1	item1	experimental	1	610
subj1	item2	control	0	600
subj1	item2	experimental	1	610
subj2	item1	control	0	600
subj2	item1	experimental	1	610
subj2	item2	control	0	600
subj2	item2	experimental	1	610

Error Term

error_sd <- 20 # SD of trial-level error (residuals)

lmem_dat <- design |>
  add_ranef(err = error_sd) |>
  mutate(dv = intercept + err +
           cond * cond_eff)

subj	item	condition	cond	err	dv
subj1	item1	control	0	-6.70	593.30
subj1	item1	experimental	1	-28.93	581.07
subj1	item2	control	0	33.56	633.56
subj1	item2	experimental	1	-25.67	584.33
subj2	item1	control	0	-10.48	589.52
subj2	item1	experimental	1	-13.40	596.60
subj2	item2	control	0	18.31	618.31
subj2	item2	experimental	1	-32.30	577.70

Random Intercepts

subj_sd      <-   8 # SD of subject-level intercepts
item_sd      <-   4 # SD of item-level intercepts

lmem_dat <- design |>
  add_ranef(err = error_sd) |>
  add_ranef(.by = "subj", subj_i = subj_sd) |>
  add_ranef(.by = "item", item_i = item_sd) |>
  mutate(dv = intercept + subj_i + item_i + err +
           cond * cond_eff)

subj	item	condition	cond	err	subj_i	item_i	dv
subj1	item1	control	0	19.46	-6.76	-3.19	609.51
subj1	item1	experimental	1	-27.82	-6.76	-3.19	572.23
subj1	item2	control	0	-42.91	-6.76	0.06	550.38
subj1	item2	experimental	1	-7.91	-6.76	0.06	595.38
subj2	item1	control	0	9.38	-5.79	-3.19	600.41
subj2	item1	experimental	1	25.35	-5.79	-3.19	626.37
subj2	item2	control	0	-25.40	-5.79	0.06	568.86
subj2	item2	experimental	1	-7.71	-5.79	0.06	596.56

Random Slopes

subj_cond_sd <-   5 # SD of subject-level condition effect size
item_cond_sd <-  15 # SD of item-level condition effect size
subj_cors    <-   0.5 # correlation between subject intercept and slope
item_cors    <-  -0.5 # correlation between item intercept and slope

lmem_dat <- design |>
  add_ranef(err = error_sd) |>
  add_ranef(.by = "subj", 
            subj_i = subj_sd, 
            subj_cond = subj_cond_sd, 
            .cors = subj_cors) |>
  add_ranef(.by = "item", 
            item_i = item_sd, 
            item_cond = item_cond_sd, 
            .cors = item_cors) |>
  mutate(dv = intercept + subj_i + item_i + err +
           cond * (cond_eff + subj_cond + item_cond))

Random Slopes

subj	item	condition	cond	err	subj_i	subj_cond	item_i	item_cond	dv
subj1	item1	control	0	24.08	21.20	1.30	1.70	-9.83	646.98
subj1	item1	experimental	1	27.70	21.20	1.30	1.70	-9.83	652.09
subj1	item2	control	0	13.80	21.20	1.30	-0.90	21.20	634.10
subj1	item2	experimental	1	-5.75	21.20	1.30	-0.90	21.20	647.06
subj2	item1	control	0	-9.94	16.84	4.07	1.70	-9.83	608.60
subj2	item1	experimental	1	-28.29	16.84	4.07	1.70	-9.83	594.49
subj2	item2	control	0	-15.51	16.84	4.07	-0.90	21.20	600.44
subj2	item2	experimental	1	-28.15	16.84	4.07	-0.90	21.20	623.07

Parameters

subj_n       <- 100   # number of subjects
item_n       <-  20   # number of items
intercept    <- 600   # model intercept (mean for control condition)
cond_eff     <-  10   # condition effect size
error_sd     <-  20   # SD of trial-level error (residuals)
subj_sd      <-   8   # SD of subject-level intercepts
item_sd      <-   4   # SD of item-level intercepts
subj_cond_sd <-   5   # SD of subject-level condition effect size
item_cond_sd <-  15   # SD of item-level condition effect size
subj_cors    <-   0.5 # correlation between subject intercept and slope
item_cors    <-  -0.5 # correlation between item intercept and slope

Full Data Simulation Code

lmem_dat <- add_random(subj = subj_n) |>
  add_random(item = item_n) |>
  add_within(condition = c("control", "experimental")) |>
  add_contrast("condition", "treatment", colnames = "cond") |>
  add_ranef(err = error_sd) |>
  add_ranef(.by = "subj", 
            subj_i = subj_sd, 
            subj_cond = subj_cond_sd, 
            .cors = subj_cors) |>
  add_ranef(.by = "item", 
            item_i = item_sd, 
            item_cond = item_cond_sd, 
            .cors = item_cors) |>
  mutate(dv = intercept + subj_i + item_i + err +
           cond * (cond_eff + subj_cond + item_cond))

LMEM Analysis

lmer(dv ~ cond + 
       (1 + cond | subj) + 
       (1 + cond | item),
     data = lmem_dat) |> summary()

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: dv ~ cond + (1 + cond | subj) + (1 + cond | item)
   Data: lmem_dat

REML criterion at convergence: 35619.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.3660 -0.6485  0.0197  0.6725  3.4046 

Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 subj     (Intercept)  64.79    8.049         
          cond         20.07    4.480   0.44  
 item     (Intercept)  10.78    3.283         
          cond        196.50   14.018   -0.28 
 Residual             395.43   19.885         
Number of obs: 4000, groups:  subj, 100; item, 20

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)  599.565      1.177  53.557 509.521   <2e-16 ***
cond           3.999      3.228  19.740   1.239     0.23    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
cond -0.178

Power simulation

1. Wrap dataset creation and analysis in a function that returns the values you care about as a data frame
2. Iterate this function and save the output as a data frame
3. Summarise the output (e.g., power, range of effect sizes)

Function Outline

sim_func <- function(iteration = 0, subj_n = 100, item_n = 20, cond_eff = 0) {
  # variables not in function arguments

  # simulate the data
  lmem_dat <- ...
  
  # run the analysis
  mod <- lmer(...)
  
  # return a table of fixed effects only
  broom.mixed::tidy(mod, effects = "fixed") |>
    mutate(iteration = iteration)
}

Full Function

sim_func <- function(iteration = 0, subj_n = 100, item_n = 20, cond_eff = 0) {
  # variables not in function arguments
  intercept    <- 600
  error_sd     <-  20
  subj_sd      <-   8
  item_sd      <-   4
  subj_cond_sd <-   5
  item_cond_sd <-  15
  subj_cors    <-   0.5
  item_cors    <-  -0.5

  # simulate the data
  lmem_dat <- add_random(subj = subj_n) |>
    add_random(item = item_n) |>
    add_within(condition = c("control", "experimental")) |>
    add_contrast("condition", "treatment", colnames = "cond") |>
    add_ranef(err = error_sd) |>
    add_ranef(.by = "subj", 
              subj_i = subj_sd, 
              subj_cond = subj_cond_sd, 
              .cors = subj_cors) |>
    add_ranef(.by = "item", 
              item_i = item_sd, 
              item_cond = item_cond_sd, 
              .cors = item_cors) |>
    mutate(dv = intercept + subj_i + item_i + err +
             cond * (cond_eff + subj_cond + item_cond))
  
  # run the analysis
  mod <- lmer(dv ~ cond + 
       (1 + cond | subj) + 
       (1 + cond | item),
     data = lmem_dat)
  
  # return a table of fixed effects only
  broom.mixed::tidy(mod, effects = "fixed") |>
    mutate(iteration = iteration)
}

Test the Function

sim_func()

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	599.62	1.16	515.29	45.88	0.00
fixed	cond	2.21	2.86	0.77	19.51	0.45

# again to make sure iterations differ
sim_func()

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	601.10	1.20	501.51	43.59	0.00
fixed	cond	-8.00	3.05	-2.62	20.01	0.02

# change effect size
sim_func(cond_eff = 10)

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	600.49	1.24	483.81	40.21	0.00
fixed	cond	10.98	3.18	3.46	19.65	0.00

Iterate

set.seed(8675309)

sim_results <- map_df(1:10, sim_func, cond_eff = 10)

effect	term	estimate	std.error	statistic	df	p.value	iteration
fixed	(Intercept)	598.97	0.99	604.31	67.08	0.00	1
fixed	cond	11.69	3.18	3.67	20.22	0.00	1
fixed	(Intercept)	602.87	1.09	552.87	57.20	0.00	2
fixed	cond	7.64	2.85	2.68	19.88	0.01	2
fixed	(Intercept)	600.08	1.43	419.96	33.76	0.00	3
fixed	cond	6.06	4.42	1.37	19.22	0.19	3
fixed	(Intercept)	600.82	1.33	453.44	43.60	0.00	4
fixed	cond	8.67	4.36	1.99	19.76	0.06	4
fixed	(Intercept)	598.51	1.35	444.71	41.22	0.00	5
fixed	cond	13.69	3.31	4.14	20.13	0.00	5

Summarise

alpha <- 0.05 # critical alpha

sim_results |>
  summarise(power = mean(p.value < alpha), .by = term)

term	power
(Intercept)	1.0
cond	0.7

Change Parameters

Try increasing the SD for subject and item intercepts and slopes by 20%.

sim_func2 <- function(iteration = 0, subj_n = 100, item_n = 20, cond_eff = 0) {
  # variables not in function arguments
  intercept    <- 600
  error_sd     <-  20
  subj_sd      <-   8 * 1.2
  item_sd      <-   4 * 1.2
  subj_cond_sd <-   5 * 1.2
  item_cond_sd <-  15 * 1.2
  subj_cors    <-   0.5
  item_cors    <-  -0.5

  # continue the function as before...
}

term	power
(Intercept)	1.0
cond	0.5

Write Up

We conducted a sensitivity power analysis using custom simulation code and the R packages faux (DeBruine, 2026) and lme4 (Bates, Maechler, Bolker & Walker, 2015) to determine the power. With 100 subjects and 20 items, an effect of condition of 10 milliseconds, using treatment coding (control = 0, experimental = 1), the expected parameters in Table 1, and a critical alpha of 0.05, simulation using 10 replications resulted in a power of 70%.

Table 1: Model parameters for simulation were guessed at random

Model Parameter	Value
SD of trial-level error (residuals)	20
SD of subject-level intercepts	8
SD of item-level intercepts	4
SD of subject-level condition effect	5
SD of item-level condition effect	15
Correlation between subject-level intercept and condition effect	0.5
Correlation between item-level intercept and condition effect	-0.5

Further Resources

PsyPag Simulation Summer School

Data Simulation Workshops

Thank You!

https://scienceverse.org/talks/2026-power-sim/