simulate_glm.R

# Show that interaction in glm() changes nature of main effect
# (only if a categorical predictor is dummy coded - not contrast coded)

# Returns the p-value associated with a predictor main effect, once
# with and once without interaction with a (non-predictive) categorical
# independent variable

simulate_glm <- function(N = 100, contrast_coding = FALSE) {
  iv1 <- rnorm(N) # related to DV
  iv2 <- rbinom(N, 1, 0.5) # not related to DV
  if (contrast_coding) {
    iv2 <- ifelse(iv2, -.5, .5) 
  }
  # construct a dependent variable that is correlated with IV1
  dv <- rbinom(N, 1, (iv1-min(iv1))/(max(iv1)-min(iv1)))

  m1 <- glm(dv ~ iv1, family = binomial)
  p1 <- summary(m1)$coefficients["iv1", "Pr(>|z|)"]

  m2 <- glm(dv ~ iv1 * iv2, family = binomial)
  p2 <- summary(m2)$coefficients["iv1", "Pr(>|z|)"]
  c(p1, p2)
}

sim1 <- replicate(1000, simulate_glm())
rowMeans(sim1 < .01)

sim2 <- replicate(1000, simulate_glm(contrast_coding = TRUE))
rowMeans(sim2 < .01)