Simulate power and bias of ANCOVA vs change score with the same measurement errors of pre an post and some correlation

# In obs studies, ANCOVA has more power than the change score, but bias is variable. In RCT, ANCOVA increases efficiency

simulate_data <- function(n = 100) {

  # Simulate baseline scores and add measurement error

  baseline <- rnorm(n, mean = 50, sd = 10) + rnorm(n, mean = 0, sd = 2)

  

  # Introduce a treatment effect with some noise, ensuring a correlation of r=0.5

  treatment_effect <- 5

  noise <- rnorm(n, 0, sqrt(10^2 - (0.5 * 10)^2))

  post <- baseline * 0.5 + treatment_effect + noise + rnorm(n, mean = 0, sd = 2)

  

  group <- ifelse(runif(n) > 0.5, "treatment", "control")

  post[group == "control"] <- post[group == "control"] - treatment_effect

  data <- data.frame(id = 1:n, group = group, baseline = baseline, post = post)

  return(data)

}

n_simulations <- 1000

alpha <- 0.05

power_ANCOVA <- 0

power_change_score <- 0

# To store estimated effects for both methods across simulations

estimates_ANCOVA <- numeric(n_simulations)

estimates_change_score <- numeric(n_simulations)

for (i in 1:n_simulations) {

  # Simulate data

  data <- simulate_data()

  # ANCOVA

  fit <- lm(post ~ baseline + group, data = data)

  estimates_ANCOVA[i] <- coef(fit)["grouptreatment"]

  

  # Change Score

  change_score <- data$post - data$baseline

  fit_change <- lm(change_score ~ data$group)

  estimates_change_score[i] <- coef(fit_change)["data$grouptreatment"]

  

  # Power calculations

  if (coef(summary(fit))["grouptreatment", "Pr(>|t|)"] < alpha) {

    power_ANCOVA <- power_ANCOVA + 1

  }

  if (coef(summary(fit_change))["data$grouptreatment", "Pr(>|t|)"] < alpha) {

    power_change_score <- power_change_score + 1

  }

}

# Bias calculations

bias_ANCOVA <- mean(estimates_ANCOVA) - 5

bias_change_score <- mean(estimates_change_score) - 5

power_ANCOVA <- power_ANCOVA / n_simulations

power_change_score <- power_change_score / n_simulations

print(paste("Power (ANCOVA):", power_ANCOVA))

print(paste("Power (Change Score):", power_change_score))

print(paste("Bias (ANCOVA):", bias_ANCOVA))

print(paste("Bias (Change Score):", bias_change_score))