DeclareDesign code from Alex Coppock to simulate SEs from pilot studies

pilot-se.R

library(tidyverse)
library(DeclareDesign)

# Original parameters
tau <- 1  # treatment effect
se <- tau / (qnorm(0.95) + qnorm(0.80))  # standard error for 80% power
n_planned <- 500  # sample size per condition in planned full study
# calculate required standard deviation to yield 80% power given the above
sigma <- se * sqrt(2 * n_planned / 2)

# Declare the design
design <- 
  declare_model(N = 2*n_planned,
                potential_outcomes(Y ~ rnorm(N, tau*Z, sigma))) +
  declare_assignment(Z = complete_ra(N)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z)

# We're after power for a one-tailed test
diagnosands <- 
  declare_diagnosands(
    p_one_tailed_power = mean(pt(statistic, df, lower.tail = FALSE) <= 0.05))

# This function simulates and gives bootstrapped SEs for simulation uncertainty
diagnosis <- diagnose_design(design,
                             diagnosands = diagnosands)

# compare over sample sizes
n_pilot_values <- c(10, 30, 60, 90, 150)
pilots <- redesign(design, n_planned = c(10, 30, 60, 90, 150))
simulations <- simulate_designs(pilots)

# plot the standard errors
ggplot(simulations, aes(std.error)) +
  geom_histogram() +
  facet_wrap(vars(n_planned))