DeclareDesign simulation to see how power changes as we increase the number of arms (with corrections)

power_depends_on_arms

library(DeclareDesign)
library(tidyverse)
library(geomtextpath)

design_1 <-
  declare_model(N = 1000,
                U = rnorm(N),
                Y_Z_0 = U,
                Y_Z_1 = U + 0.2) +
  declare_assignment(Zk = complete_ra(N = N, num_arms = k)) +
  declare_measurement(Z = if_else(Zk == "T1", 0, 1),
                      Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Zk, term = TRUE)

designs_1 <- redesign(design_1, k = 2:10)
simulations_1 <- simulate_designs(designs_1, sims = 1000)


simulations_1 <-
  simulations_1 |>
  filter(term != "(Intercept)") |>
  group_by(k, sim_ID) |>
  mutate(p_bonferroni = p.adjust(p = p.value, method = "bonferroni"),
         p_holm = p.adjust(p = p.value, method = "holm"),
         p_fdr = p.adjust(p = p.value, method = "fdr"),
         "Naive p-values" = p.value <= 0.05,
         "Bonferroni correction" = p_bonferroni <= 0.05,
         "Holm correction" = p_holm <= 0.05,
         "FDR correction" = p_fdr <= 0.05,
         )

gg_df_1 <-
  simulations_1 |>
  filter(term == "ZkT2") |>
  group_by(k) |>
  summarise(tidy(lm_robust(
    cbind(`Naive p-values`,
          `Bonferroni correction`,
          `Holm correction`,
          `FDR correction`) ~ 1
  )))


g_1 <-
  ggplot(gg_df_1) +
  aes(k, estimate, color = outcome) +
  geom_ribbon(
    aes(ymin = conf.low, ymax = conf.high, fill = outcome),
    color = NA,
    alpha = 0.1
  ) +
  geom_textpath(aes(label = outcome), offset = unit(x = 2, units = "points")) +
  geom_line() +
  theme_minimal() +
  theme(legend.position = "none",plot.background = element_rect(fill = "white"),
        panel.grid.minor = element_blank()) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1),breaks = seq(0, 1, 0.2)) +
  scale_x_continuous(name = 'Total number of arms (k)', breaks = 2:10,
                     sec.axis = sec_axis(~1000/., name = 'Number of subjects in each arm',
                                         breaks = round(1000/2:10))) +
  labs(y = "Power for treatment 1 versus control",
       title = "How much does adding an extra treatment arm hurt my power?",
       subtitle = "Total N = 1000, all (k - 1) treatments have a 0.2SD effect")

ggsave("corrections_1.png", plot = g_1, width = 6.5, height = 5)



design_2 <-
  declare_model(N = 500*k,
                U = rnorm(N),
                Y_Z_0 = U,
                Y_Z_1 = U + 0.2) +
  declare_assignment(Zk = complete_ra(N = N, num_arms = k)) +
  declare_measurement(Z = if_else(Zk == "T1", 0, 1),
                      Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(Y ~ Zk, term = TRUE)

designs_2 <- redesign(design_2, k = 2:10)
simulations_2 <- simulate_designs(designs_2, sims = 1000)


simulations_2 <-
  simulations_2 |>
  filter(term != "(Intercept)") |>
  group_by(k, sim_ID) |>
  mutate(p_bonferroni = p.adjust(p = p.value, method = "bonferroni"),
         p_holm = p.adjust(p = p.value, method = "holm"),
         p_fdr = p.adjust(p = p.value, method = "fdr"),
         "Naive p-values" = p.value <= 0.05,
         "Bonferroni correction" = p_bonferroni <= 0.05,
         "Holm correction" = p_holm <= 0.05,
         "FDR correction" = p_fdr <= 0.05,
  )

gg_df_2 <-
  simulations_2 |>
  filter(term == "ZkT2") |>
  group_by(k) |>
  summarise(tidy(lm_robust(
    cbind(`Naive p-values`,
          `Bonferroni correction`,
          `Holm correction`,
          `FDR correction`) ~ 1
  )))



g_2 <-
  ggplot(gg_df_2) +
  aes(k, estimate, color = outcome) +
  geom_ribbon(
    aes(ymin = conf.low, ymax = conf.high, fill = outcome),
    color = NA,
    alpha = 0.1
  ) +
  geom_textpath(aes(label = outcome), offset = unit(x = 2, units = "points")) +
  geom_line() +
  theme_minimal() +
  theme(legend.position = "none",plot.background = element_rect(fill = "white"),
        panel.grid.minor = element_blank()) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1),breaks = seq(0, 1, 0.2)) +
  scale_x_continuous(name = 'Total number of arms (k)', breaks = 2:10,
                     sec.axis = sec_axis(~500*., name = 'Total number of subjects',
                                         breaks = 500*2:10)) +
  labs(y = "Power for treatment 1 versus control",
       title = "How much does adding an extra treatment arm hurt my power?",
       subtitle = "Total N = 500 * k, all (k - 1) treatments have a 0.2SD effect")

ggsave("corrections_2.png", plot = g_2, width = 6.5, height = 5)