causal-mediation

# Causal Mediation Analysis in R ## Overview Causal mediation analysis methods for decomposing total effects into direct and indirect effects. Covers traditional approaches, natural effect models, sensitivity analysis for unmeasured confounding, mediation with survival outcomes, and comprehensive causal mediation frameworks. ## Traditional Mediation (Baron-Kenny) ### Basic Approach ```r # Baron-Kenny steps for mediation # Total effect: Y = c*X + e # Mediator: M = a*X + e # Outcome with mediator: Y = c'*X + b*M + e # Indirect effect: a*b # Direct effect: c' # Step 1: Total effect (X -> Y) fit_total <- lm(outcome ~ treatment + covariates, data = df) c_total <- coef(fit_total)["treatment"] # Step 2: Effect on mediator (X -> M) fit_med <- lm(mediator ~ treatment + covariates, data = df) a <- coef(fit_med)["treatment"] # Step 3: Direct effect (X -> Y | M) fit_direct <- lm(outcome ~ treatment + mediator + covariates, data = df) c_prime <- coef(fit_direct)["treatment"] b <- coef(fit_direct)["mediator"] # Effects indirect_effect <- a * b direct_effect <- c_prime total_effect <- c_total proportion_mediated <- indirect_effect / total_effect cat("Total effect:", round(total_effect, 4), "\n") cat("Direct effect:", round(direct_effect, 4), "\n") cat("Indirect effect:", round(indirect_effect, 4), "\n") cat("Proportion mediated:", round(proportion_mediated * 100, 1), "%\n") ``` ### Sobel Test (Not Recommended for Small Samples) ```r # Sobel test for indirect effect sobel_test <- function(a, b, se_a, se_b) { se_ab <- sqrt(b^2 * se_a^2 + a^2 * se_b^2) z <- (a * b) / se_ab p <- 2 * pnorm(-abs(z)) return(list(indirect = a * b, se = se_ab, z = z, p = p)) } se_a <- summary(fit_med)$coefficients["treatment", "Std. Error"] se_b <- summary(fit_direct)$coefficients["mediator", "Std. Error"] sobel <- sobel_test(a, b, se_a, se_b) print(sobel) # Note: Bootstrap methods are preferred over Sobel test ``` ## Causal Mediation with mediation Package ### Basic Mediation Analysis ```r library(mediation) # Fit mediator model med_model <- lm(mediator ~ treatment + age + sex, data = df) # Fit outcome model out_model <- lm(outcome ~ treatment + mediator + age + sex, data = df) # Mediation analysis med_result <- mediate( med_model, out_model, treat = "treatment", mediator = "mediator", boot = TRUE, boot.ci.type = "bca", sims = 1000 ) summary(med_result) # Key outputs: # ACME: Average Causal Mediation Effect (indirect) # ADE: Average Direct Effect # Total Effect # Prop. Mediated # Plot effects plot(med_result) ``` ### Binary Outcomes ```r library(mediation) # Mediator model (continuous mediator) med_model <- lm(mediator ~ treatment + covariates, data = df) # Outcome model (binary) out_model <- glm(outcome ~ treatment + mediator + covariates, family = binomial(link = "logit"), data = df) # Mediation with simulation med_result <- mediate( med_model, out_model, treat = "treatment", mediator = "mediator", boot = TRUE, sims = 1000 ) summary(med_result) ``` ### Binary Mediator ```r library(mediation) # Binary mediator model med_model <- glm(mediator ~ treatment + covariates, family = binomial, data = df) # Outcome model out_model <- lm(outcome ~ treatment + mediator + covariates, data = df) # Mediation med_result <- mediate( med_model, out_model, treat = "treatment", mediator = "mediator", boot = TRUE, sims = 1000 ) summary(med_result) ``` ### Treatment-Mediator Interaction ```r library(mediation) # Outcome model with interaction out_model <- lm(outcome ~ treatment * mediator + covariates, data = df) # Mediation with interaction med_result <- mediate( med_model, out_model, treat = "treatment", mediator = "mediator", boot = TRUE, sims = 1000 ) summary(med_result) # Note: When interaction exists, effects may differ by treatment level # ACME(0), ACME(1): Indirect effect at control and treatment # ADE(0), ADE(1): Direct effect at control and treatment ``` ## Natural Effect Models with medflex ### Imputation-Based Approach ```r library(medflex) # Expand data for natural effect estimation expData <- neWeight( outcome ~ treatment + mediator + age + sex, data = df ) # Fit natural effect model neMod <- neModel( outcome ~ treatment0 + treatment1 + age + sex, expData = expData, se = "robust" ) # Effect decomposition neEffdecomp(neMod) # Natural Direct Effect (NDE): treatment1 = 0, treatment0 varies # Natural Indirect Effect (NIE): treatment0 = 0, treatment1 varies ``` ### Weighting-Based Approach ```r library(medflex) # Using inverse probability weighting expData <- neWeight( outcome ~ treatment + mediator + age + sex, data = df, weights = "estimated" ) neMod <- neModel( outcome ~ treatment0 + treatment1 + age + sex, expData = expData, se = "bootstrap", nBoot = 1000 ) summary(neMod) neEffdecomp(neMod) ``` ## Comprehensive Mediation with CMAverse ### Basic CMAverse Analysis ```r library(CMAverse) # Comprehensive mediation analysis cma_result <- cmest( data = df, outcome = "outcome", exposure = "treatment", mediator = "mediator", basec = c("age", "sex", "baseline"), EMint = TRUE, # Exposure-mediator interaction model = "rb", # Regression-based estimation = "imputation", inference = "bootstrap", nboot = 1000 ) summary(cma_result) # Decomposition: # CDE: Controlled Direct Effect # PNDE: Pure Natural Direct Effect # TNDE: Total Natural Direct Effect # PNIE: Pure Natural Indirect Effect # TNIE: Total Natural Indirect Effect # TE: Total Effect # PM: Proportion Mediated ``` ### Multiple Mediators ```r library(CMAverse) # Multiple mediators (parallel) cma_multi <- cmest( data = df, outcome = "outcome", exposure = "treatment", mediator = c("mediator1", "mediator2"), basec = c("age", "sex"), EMint = TRUE, model = "rb", estimation = "imputation", inference = "bootstrap", nboot = 500 ) summary(cma_multi) ``` ### Sequential Mediators ```r library(CMAverse) # Sequential/serial mediation # X -> M1 -> M2 -> Y cma_seq <- cmest( data = df, outcome = "outcome", exposure = "treatment", mediator = c("mediator1", "mediator2"), basec = c("age", "sex"), EMint = TRUE, model = "rb", mreg = list("linear", "linear"), # Mediator regressions yreg = "linear", # Outcome regression estimation = "imputation", inference = "bootstrap", nboot = 500 ) summary(cma_seq) ``` ## Mediation with Survival Outcomes ### Using mediation Package ```r library(mediation) library(survival) # Mediator model med_model <- lm(mediator ~ treatment + age + sex, data = df) # Survival outcome model surv_model <- coxph(Surv(time, event) ~ treatment + mediator + age + sex, data = df) # Mediation analysis (simulation-based) med_surv <- mediate( med_model, surv_model, treat = "treatment", mediator = "mediator", sims = 1000 ) summary(med_surv) ``` ### Survival Mediation with CMAverse ```r library(CMAverse) # Survival outcome cma_surv <- cmest( data = df, outcome = "event", event = "event", eventtime = "time", exposure = "treatment", mediator = "mediator", basec = c("age", "sex"), EMint = TRUE, model = "rb", yreg = "coxph", estimation = "paramfunc", inference = "bootstrap", nboot = 500 ) summary(cma_surv) ``` ## Sensitivity Analysis ### Sensitivity to Unmeasured Confounding ```r library(mediation) # Sensitivity analysis for sequential ignorability sens_result <- medsens( med_result, rho.by = 0.05, effect.type = "indirect", sims = 1000 ) summary(sens_result) plot(sens_result) # Key output: At what correlation (rho) between residuals # does the indirect effect become non-significant? ``` ### E-value for Mediation ```r library(EValue) # E-value for natural indirect effect # Convert to risk ratio scale if needed # For continuous outcome: nie_estimate <- med_result$d0 # NIE point estimate nie_se <- med_result$d0.se # NIE standard error nie_ci_lower <- med_result$d0.ci[1] # Standardized effect to approximate RR # (Simplified - actual conversion depends on outcome distribution) ``` ### Sensitivity in CMAverse ```r library(CMAverse) # Sensitivity analysis for unmeasured confounding cma_sens <- cmsens( object = cma_result, sens = "uc", # Unmeasured confounding eval_uc = list( ume_effect = seq(-0.5, 0.5, 0.1), # Effect on mediator uye_effect = seq(-0.5, 0.5, 0.1) # Effect on outcome ) ) summary(cma_sens) plot(cma_sens) ``` ## Moderated Mediation ### Conditional Process Analysis ```r library(mediation) # Outcome model with moderator interaction out_model_mod <- lm(outcome ~ treatment * moderator + mediator + covariates, data = df) # Mediation at different moderator levels med_low <- mediate( med_model, out_model_mod, treat = "treatment", mediator = "mediator", covariates = list(moderator = quantile(df$moderator, 0.25)), boot = TRUE, sims = 500 ) med_high <- mediate( med_model, out_model_mod, treat = "treatment", mediator = "mediator", covariates = list(moderator = quantile(df$moderator, 0.75)), boot = TRUE, sims = 500 ) # Compare summary(med_low) summary(med_high) # Test moderated mediation test.modmed(med_result, list(moderator = "low"), list(moderator = "high")) ``` ## Causal Diagrams for Mediation ```r library(dagitty) library(ggdag) # Simple mediation DAG dag_mediation <- dagitty(' dag { X -> M M -> Y X -> Y C -> X C -> M C -> Y } ') # Identify adjustment sets adjustmentSets(dag_mediation, exposure = "X", outcome = "Y") # Plot DAG ggdag(dag_mediation) + theme_dag() + labs(title = "Mediation DAG") # Check if mediation is identified # Need: No unmeasured confounding of X-Y, X-M, M-Y relationships ``` ## Reporting Mediation Results ```r # Create mediation summary table create_mediation_table <- function(med_result) { effects <- data.frame( Effect = c("ACME (Indirect)", "ADE (Direct)", "Total Effect", "Proportion Mediated"), Estimate = c(med_result$d0, med_result$z0, med_result$tau.coef, med_result$n0), CI_Lower = c(med_result$d0.ci[1], med_result$z0.ci[1], med_result$tau.ci[1], med_result$n0.ci[1]), CI_Upper = c(med_result$d0.ci[2], med_result$z0.ci[2], med_result$tau.ci[2], med_result$n0.ci[2]),