base-simpr-dd.R

library(simpr)
library(DeclareDesign)
library(tidyverse)

# base ------------------------------------------------------

sims <- 500
ns <- seq(100, 300, by = 50)
p.values <- matrix(NA, nrow = sims, ncol = length(ns))

for (j in seq_along(ns)) {
  n <- ns[j]
  for (i in 1:sims) {
    Z <- rbinom(n, 1, 0.5)
    U <- rnorm(n)
    Y <- 0.2 * Z + U
    p.values[i, j] <- summary(lm_robust(Y ~ Z))$coefficients[2, 4]
  }
}

power <- colMeans(p.values <= 0.05)
power

plot(ns, power)

# simpr -----------------------------------------------------

simpr_design <-
  blueprint(Z = ~ rbinom(n, 1, 0.5),
            U = ~ rnorm(n),
            Y = ~ 0.2 * Z + U)

simpr_simulations <-
  simpr_design %>%
  meta(n = seq(100, 300, by = 50)) %>%
  produce_sims(500) %>%
  fit(lm = ~ lm_robust(Y ~ Z, data = .)) %>%
  tidy_fits


simpr_simulations %>% 
  filter(term == "Z") %>% 
  group_by(n) %>% 
  summarize(power = mean(p.value < 0.05)) %>%
  ggplot(aes(n, power)) +
  geom_line()


# DeclareDesign ---------------------------------------------

dd_design <-
  declare_model(
    N = N,
    U = rnorm(N),
    potential_outcomes(Y ~  0.2 * Z + U)
  ) +
  declare_inquiry(ATE = mean(Y_Z_1 - Y_Z_0)) +
  declare_assignment(Z = simple_ra(N)) +
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Z, inquiry = "ATE")

dd_simulations <-
  dd_design %>%
  redesign(N = seq(100, 300, by = 50)) %>%
  simulate_designs(sims = 500)

dd_simulations %>% 
  group_by(N) %>% 
  summarize(power = mean(p.value < 0.05),
            rmse = sqrt(mean((estimate - estimand)^2))) %>%
  pivot_longer(c(power, rmse)) %>% 
  ggplot(aes(N, value)) +
  geom_line() +
  facet_wrap(~name, scales = "free_y")