Increasing Sample Size Doesn't Fix Bias

increasing_sample_size_bias.R

# Title: Increasing Sample Size Doesn't Fix Bias

# Description: Demonstrating that increasing sample size doesn't fix bias. 
# Note: this code was written with the help of Gemini

# Setup ----

#.. Libraries ----

library(tidyverse) # ol faithful
library(broom) # tidying output from model 

#... Functions ----

# Simulating Data 

sim_data <- function(n = 250, # sample size 
                     beta_trt = 1.5, # treatment effect
                     # Parameters for Z1 
                     z1_mean = 5, z1_sd = 2, 
                     # Parameters for Z2
                     z2_size = 1, z2_prob = 0.5, 
                     # Confounder - Effect on X
                     z1_on_x = 0.05, z2_on_x = 0.2,
                     # Confounder - Effect on Y
                     z1_on_y = 0.5, z2_on_y = 0.3){
  
  # Creating the Dataframe
  
  df <- data.frame(
    z1 = rnorm(n = n, mean = z1_mean, sd = z1_sd), 
    z2 = rbinom(n = n, size = z2_size, prob = z2_prob)
  ) %>% 
    dplyr::mutate(
      prob = plogis(z1_on_x*z1 + z2_on_x*z2), 
      x = rbinom(n = n, size = 1, prob = prob), 
      y = beta_trt*x + z1_on_y*z1 + z2_on_y*z2 + rnorm(n = n, mean = 0, sd = 1)
    )
  
  # Return the dataframe
  
  return(df)
}

# Simulation Setup ----

set.seed(456) # setting seed
sample_sizes <- c(100, 500, 1000, 10000) # various sample sizes
true_beta <- 1.5 # true effect 
iterations <- 1000 # number of times to repeat it 

# Simulating Data & Fitting Model ----

# Running the 1,000 simulations per N
results_long <- map_df(sample_sizes, function(n_val) {
  
  replicate(iterations, {
    data <- sim_data(n = n_val, beta_trt = true_beta)
    
    # Updated formulas to include the .x suffix
    m_biased  <- glm(y ~ x, data = data) %>% broom::tidy() %>% filter(term == "x") %>% select(estimate)
    m_correct <- glm(y ~ x + z1 + z2, data = data) %>% broom::tidy() %>% filter(term == "x") %>% select(estimate)
    
    c(Biased = m_biased, Correct = m_correct)
  }, simplify = FALSE) %>%
    bind_rows() %>%
    mutate(N = n_val)
})


# Calculating Bias ----

# Calculating bias for both the correct and incorrect models. 
# These formulas are based on the Morris et al. (2019) paper. 

# Summary with MCSE calculation
bias_summary <- results_long %>%
  mutate(
    bias_for_one = Biased.estimate - true_beta,
    bias_correct_model = Correct.estimate - true_beta,
    squared = (bias_for_one - mean(bias_for_one))^2,
    squared_correct_model = (bias_correct_model - mean(bias_correct_model))^2
  ) %>% 
  group_by(N) %>%
  summarise(
    # Incorrect Model 
    bias = mean(bias_for_one), 
    sqrt(sum(squared))*(1 / (1000*999)),
    # Correct Model 
    bias_cm = mean(bias_correct_model), 
    sqrt(sum(squared_correct_model))*(1 / (1000*999))
  ) 

print(as.data.frame(bias_summary))