發表文章

目前顯示的是 9月, 2023的文章

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 simu...

A Visual Diagnostic Tool for Causal Inference

 # Simulating the data with a heterogeneous treatment effect. Fitting a misspecified model without interaction. Generating residuals. Plotting residuals vs fits colored/faceted by treatment. Refitting a correct model with interaction. Re-plotting residuals vs fits. This shows how the diagnostic plot can detect model misspecification. library(tidyverse) ## Simulate data set.seed(8)   n <- 1000 x <- rnorm(n) t <- rbinom(n, 1, exp(x) / (1 + exp(x))) y1 <- 0.5 * x + rnorm(n) y0 <- - 0.5 * x + rnorm(n) y_obs <- t * y1 + (1 - t) * y0 d <- tibble(   y = y_obs,   p = predict(lm(t ~ x)),    t = as.factor(t),   f = predict(lm(y_obs ~ t + x)),   r = resid(lm(y_obs ~ t + x)),   f_correct = predict(lm(y_obs ~ t*x)),   r_correct = resid(lm(y_obs ~ t*x)) ) ## Plots (p_naive <- ggplot(d, aes(x = f, y = r)) +   geom_point(alpha = 0.25) +   geom_smooth(     color = "black",     method = "loess",  ...