This document shows one way to do simulations to do power for a continuous interaction (inspired by this tweet https://twitter.com/NathanKalmoe/status/1063839116742283264) Here we’re simulating an average treatment effect of 0.2 SDs and an interaction of 0.05SDs. That means people who score “3” on our covariate have an effect of 0.35 SDs and those who score “-3” have an effect of 0.05 SDs.We do the simulations in
DeclareDesign. The power for the interaction term’s not great at small sample sizes, even though this is a pretty hefty interaction effect.
Declare your design
design <- declare_population( N = N, pid_7 = sample(-3:3, size = N, replace = TRUE), noise = rnorm(N) ) + declare_potential_outcomes(Y ~ 0.2 * Z + 0.2 * pid_7 + 0.05 * pid_7 * Z + noise) + declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0), interaction = 0.05) + declare_assignment(prob = 0.5) + declare_reveal(Y, Z) + declare_estimator( Y ~ Z + pid_7 + Z * pid_7, model = lm_robust, term = c("Z", "Z:pid_7"), estimand = c("ATE", "interaction") )
Simulate at varying sample sizes
designs <- redesign(design, N = c(500, 1000, 3000, 5000)) simulations <- simulate_design(designs, sims = 1000)
summary_df <- simulations %>% group_by(N, estimand_label) %>% summarize(power = mean(p.value <= 0.05), mean_estimate = mean(estimate))
ggplot(summary_df, aes(N, power, group = estimand_label, color = estimand_label)) + geom_point() + geom_line() + geom_hline(yintercept = 0.8) + theme_bw()