simulate_linear_model.r

# sample size
n <- 500

# simulate data for each predictor 
# (note predictor b is categorical)
predictor_a <- rnorm(n, 0, 1) 
predictor_b <- rbinom(n, 1, .5)

# population-level beta values
b_predictor_a <- -0.20
b_predictor_b <- 0
b_int <- 0.20

# simulate outcome variable
y <- b_predictor_a * predictor_a + b_predictor_b * predictor_b + 
  (predictor_a * predictor_b * b_int) + rnorm(n, 0, 1)

# collate into data frame
data <- data.frame(
  y = y, 
  predictor_a = predictor_a, 
  predictor_b = predictor_b
)

# fit linear model
summary(lm(y ~ predictor_a * predictor_b, data = data))