MB4 GLM simulation

mb4_sim.R

library(tidyverse) 
n_sim <- 100
sims <- expand_grid(n_total = seq(50,500,25), 
                    i = 1:n_sim) %>%
  mutate(idx = 1:n()) %>%
  split(.$idx) %>%
  map_df(function (df) {
    cntl_sim <- tibble(choice = c(rbinom(n = df$n_total/2, size = 1, p = .68), 
                                  rbinom(n = df$n_total/2, size = 1, p = .5)), 
                       condition = c(rep("social", df$n_total/2), 
                                     rep("nonsocial", df$n_total/2)))
    nocntl_sim <- tibble(choice = rbinom(n = df$n_total/2, size = 1, p = .68))
    
    cntl <- summary(glm(choice ~ condition, family = "binomial", data = cntl_sim))
    nocntl <- summary(glm(choice ~ 1, family = "binomial", data = nocntl_sim))
    
    df$p_cntl <- cntl$coefficients["conditionsocial", "Pr(>|z|)"]
    df$p_nocntl <- nocntl$coefficients["(Intercept)", "Pr(>|z|)"]
    return(df)
  })
  
sims %>%
  group_by(n_total) %>%
  summarise(cntl_power = mean(p_cntl < .05), 
            nocntl_power = mean(p_nocntl < .05)) %>%
  pivot_longer(names_to = "condition", values_to = "power", contains("power")) %>%
  mutate(condition = ifelse(condition == "cntl_power", "With Social Control", "No Social Control")) %>%
  ggplot(aes(x = n_total, y = power, col = condition)) + 
  geom_line() + 
  geom_hline(yintercept = .8, lty = 2) + 
  xlab("N Total") + ylab("Power") +
  theme(legend.position = "bottom")