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.

# Load packages

```
library(DeclareDesign)
library(tidyverse)
```

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

# Summarize simulations

```
summary_df <-
simulations %>%
group_by(N, estimand_label) %>%
summarize(power = mean(p.value <= 0.05),
mean_estimate = mean(estimate))
```

# Power plot

```
ggplot(summary_df, aes(N, power, group = estimand_label, color = estimand_label)) +
geom_point() +
geom_line() +
geom_hline(yintercept = 0.8) +
theme_bw()
```