并行化中介分析遇foreach/dopar错误：循环内对象无法找到

需对约70万组变量组合开展中介分析，单循环代码在子集测试中运行正常，但组合量过大导致运行效率极低。尝试用foreach+doParallel并行化时抛出错误：Error in { : task 1 failed - "object 'model1' not found"，但单独指定i=1运行循环内代码可正常执行，推测worker进程无法访问model1变量。现寻求该错误的解决方法，同时欢迎其他并行化方案（如bbapply、分块处理等）建议。技术环境：MacOS Big Sur，R版本4.1.1。

单循环代码 #

library(mediation)
library(broom)

# Generate random data frames
df1 <- data.frame(x = rnorm(10), y = rnorm(10))
df2 <- data.frame(a = rnorm(10), b = rnorm(10))
df3 <- data.frame(q = rnorm(10), r = rnorm(10))

#specify column names to generate the desired combinations 
cols1 <- names(df1)
cols2 <- names(df2)
cols3 <- names(df3)

#generate the combinations
combinations <- expand.grid(cols1, cols2, cols3)

# Initialize a data frame to store the results
results <- data.frame(independent = character(),
                      mediator = character(),
                      dependent = character(),
                      ab = numeric(),
                      ac = numeric(),
                      bc = numeric(),
                      p.value = numeric(),
                      stringsAsFactors = FALSE)

# Loop through each combination and perform a nonparametric causal mediation analysis
for (i in 1:nrow(combinations)) {
  independent <- as.character(combinations[i, 1])
  mediator <- as.character(combinations[i, 2])
  dependent <- as.character(combinations[i, 3])
  
  
  # Combine the independent and mediator variables into a single data frame
  m.data <- data.frame(df1[independent], df2[mediator], df3[dependent])
 
  
  independent <- names(m.data)[1]
  mediator <- names(m.data)[2]
  dependent <- names(m.data)[3]
  
  #  # Fit a model for the independent variable and mediator
   
  model1 <- paste(mediator,dependent,sep = " ~ ") 
  model2 <- paste(dependent,"~", independent,"+",mediator,sep=' ') 
  
  model.M <- lm(model1, data=m.data) 
  model.Y <- lm(model2, data=m.data)
  
  fit <- mediate(model.M, model.Y, treat=dependent, mediator=mediator, 
                 boot=TRUE, sims=500)
  
  tidy_fit <- tidy(fit)
  # Extract the estimates of the total effect (ac) = fit$d0, the direct effect (ab),
  # and the indirect effect (bc) and their p-values
  ab <- tidy_fit$estimate[3]
  ac <- tidy_fit$estimate[1]
  bc <- tidy_fit$estimate[2]
  p.value <-  tidy_fit$p.value[1]
  
  # Add the results to the data frame
  results <- rbind(results, data.frame(independent = independent,
                                       mediator = mediator,
                                       dependent = dependent,
                                       ab = ab,
                                       ac = ac,
                                       bc = bc,
                                       p.value = p.value,
                                       stringsAsFactors = FALSE))
}

print(results)

并行化尝试代码 #

library(doParallel)
library(foreach)
cores <- 4

# Initialize a cluster
cl <- makeCluster(cores)


# Register the cluster
registerDoParallel(cl)
results <- foreach(i = 1:nrow(combinations), .combine = rbind) %dopar% {
  independent <- as.character(combinations[i, 1])
  mediator <- as.character(combinations[i, 2])
  dependent <- as.character(combinations[i, 3])
  
  independent_col <- df1[, independent]
  mediator_col <- df2[, mediator]
  dependent_col <- df3[, dependent]
  
  # Combine the independent and mediator variables into a single data frame
  m.data <- data.frame(df1[independent], df2[mediator], df3[dependent])
  # colnames(m.data) <- c("independent", "mediator", "dependent")
  
  independent <- names(m.data)[1]
  mediator <- names(m.data)[2]
  dependent <- names(m.data)[3]
  
  model1 <- paste(mediator,dependent,sep = " ~ ") 
  model2 <- paste(dependent,"~", independent,"+",mediator,sep=' ') 
  
  model.M <- lm(model1, data=m.data) 
  model.Y <- lm(model2, data=m.data)
  
  fit <- mediation::mediate(model.M, model.Y, treat = dependent, mediator = mediator,
                            boot = TRUE, sims = 500)
  
  tidy_fit <- tidy(fit)
  # Extract the estimates of the total effect (ac) = fit$d0, the direct effect (ab),
  # and the indirect effect (bc) and their p-values
  ab <- tidy_fit$estimate[3] #ADE
  ac <- tidy_fit$estimate[1] #ACME
  bc <- tidy_fit$estimate[2] #TOTAL EFFECT
  p.value <-  tidy_fit$p.value[1]
  
  # Add the results to the data frame
  data.frame(independent = independent,
             mediator = mediator,
             dependent = dependent,
             ab = ab,
             ac = ac,
             bc = bc,
             p.value = p.value,
             stringsAsFactors = FALSE)
}

stopCluster(cl)

解决方案 #

一、修复foreach并行化的错误 #

错误核心是worker进程无法获取主环境的依赖对象和包，修改后的代码如下：

library(doParallel)
library(foreach)
library(mediation)
library(broom)

cores <- 4
cl <- makeCluster(cores)
# 导出主环境的依赖数据到worker
clusterExport(cl, c("df1", "df2", "df3", "combinations"))
# 在每个worker上加载必要的包
clusterEvalQ(cl, {
  library(mediation)
  library(broom)
})

registerDoParallel(cl)

results <- foreach(i = 1:nrow(combinations), .combine = rbind) %dopar% {
  independent <- as.character(combinations[i, 1])
  mediator <- as.character(combinations[i, 2])
  dependent <- as.character(combinations[i, 3])
  
  m.data <- data.frame(df1[independent], df2[mediator], df3[dependent])
  
  # 直接从m.data列名构建公式，避免变量重复赋值的作用域混淆
  model1 <- paste(names(m.data)[2], names(m.data)[3], sep = " ~ ") 
  model2 <- paste(names(m.data)[3], "~", names(m.data)[1], "+", names(m.data)[2], sep=' ') 
  
  model.M <- lm(model1, data=m.data) 
  model.Y <- lm(model2, data=m.data)
  
  fit <- mediation::mediate(model.M, model.Y, treat = names(m.data)[1], mediator = names(m.data)[2],
                            boot = TRUE, sims = 500)
  
  tidy_fit <- tidy(fit)
  data.frame(independent = names(m.data)[1],
             mediator = names(m.data)[2],
             dependent = names(m.data)[3],
             ab = tidy_fit$estimate[3],
             ac = tidy_fit$estimate[1],
             bc = tidy_fit$estimate[2],
             p.value = tidy_fit$p.value[1],
             stringsAsFactors = FALSE)
}

stopCluster(cl)

关键修改点：

• 用clusterExport()将df1/df2/df3/combinations导出到worker进程
• 用clusterEvalQ()在worker上加载mediation和broom包（worker不会继承主环境的包加载状态）
• 直接从m.data列名构建公式，避免重复赋值变量导致的作用域混淆

二、其他并行化方案建议 #

1. furrr（基于future框架，语法简洁）

furrr是purrr的并行版，无需手动管理集群：

library(furrr)
library(mediation)
library(broom)

plan(multisession, workers = 4)
comb_list <- split(combinations, 1:nrow(combinations))

results <- future_map_dfr(comb_list, function(row) {
  independent <- as.character(row[1])
  mediator <- as.character(row[2])
  dependent <- as.character(row[3])
  
  m.data <- data.frame(df1[independent], df2[mediator], df3[dependent])
  
  model1 <- paste(names(m.data)[2], names(m.data)[3], sep = " ~ ") 
  model2 <- paste(names(m.data)[3], "~", names(m.data)[1], "+", names(m.data)[2], sep=' ') 
  
  model.M <- lm(model1, data=m.data) 
  model.Y <- lm(model2, data=m.data)
  
  fit <- mediation::mediate(model.M, model.Y, treat = names(m.data)[1], mediator = names(m.data)[2],
                            boot = TRUE, sims = 500)
  
  tidy_fit <- tidy(fit)
  data.frame(independent = names(m.data)[1],
             mediator = names(m.data)[2],
             dependent = names(m.data)[3],
             ab = tidy_fit$estimate[3],
             ac = tidy_fit$estimate[1],
             bc = tidy_fit$estimate[2],
             p.value = tidy_fit$p.value[1],
             stringsAsFactors = FALSE)
})

plan(sequential)

2. 分块处理+批量保存

适合内存有限或担心进程崩溃的场景，将70万组拆分成小块逐块处理并保存中间结果：

library(mediation)
library(broom)

chunk_size <- 7000
chunks <- split(combinations, ceiling(1:nrow(combinations)/chunk_size))

for (k in 1:length(chunks)) {
  chunk <- chunks[[k]]
  chunk_results <- data.frame()
  
  for (i in 1:nrow(chunk)) {
    independent <- as.character(chunk[i, 1])
    mediator <- as.character(chunk[i, 2])
    dependent <- as.character(chunk[i, 3])
    
    m.data <- data.frame(df1[independent], df2[mediator], df3[dependent])
    
    model1 <- paste(names(m.data)[2], names(m.data)[3], sep = " ~ ") 
    model2 <- paste(names(m.data)[3], "~", names(m.data)[1], "+", names(m.data)[2], sep=' ') 
    
    model.M <- lm(model1, data=m.data) 
    model.Y <- lm(model2, data=m.data)
    
    fit <- mediation::mediate(model.M, model.Y, treat = names(m.data)[1], mediator = names(m.data)[2],
                              boot = TRUE, sims = 500)
    
    tidy_fit <- tidy(fit)
    chunk_results <- rbind(chunk_results, data.frame(
      independent = names(m.data)[1],
      mediator = names(m.data)[2],
      dependent = names(m.data)[3],
      ab = tidy_fit$estimate[3],
      ac = tidy_fit$estimate[1],
      bc = tidy_fit$estimate[2],
      p.value = tidy_fit$p.value[1],
      stringsAsFactors = FALSE
    ))
  }
  
  saveRDS(chunk_results, paste0("mediation_chunk_", k, ".rds"))
}

all_results <- do.call(rbind, lapply(1:length(chunks), function(k) readRDS(paste0("mediation_chunk_", k, ".rds"))))

3. parallel::mclapply（Mac/Linux专属，无需手动集群）

基于fork机制，无需显式导出变量，代码更简洁：

library(parallel)
library(mediation)
library(broom)

cores <- 4
comb_list <- split(combinations, 1:nrow(combinations))

results_list <- mclapply(comb_list, function(row) {
  independent <- as.character(row[1])
  mediator <- as.character(row[2])
  dependent <- as.character(row[3])
  
  m.data <- data.frame(df1[independent], df2[mediator], df3[dependent])
  
  model1 <- paste(names(m.data)[2], names(m.data)[3], sep = " ~ ") 
  model2 <- paste(names(m.data)[3], "~", names(m.data)[1], "+", names(m.data)[2], sep=' ') 
  
  model.M <- lm(model1, data=m.data) 
  model.Y <- lm(model2, data=m.data)
  
  fit <- mediation::mediate(model.M, model.Y, treat = names(m.data)[1], mediator = names(m.data)[2],
                            boot = TRUE, sims = 500)
  
  tidy_fit <- tidy(fit)
  data.frame(independent = names(m.data)[1],
             mediator = names(m.data)[2],
             dependent = names(m.data)[3],
             ab = tidy_fit$estimate[3],
             ac = tidy_fit$estimate[1],
             bc = tidy_fit$estimate[2],
             p.value = tidy_fit$p.value[1],
             stringsAsFactors = FALSE)
}, mc.cores = cores)

results <- do.call(rbind, results_list)

额外优化建议 #

• 降低sims参数：如果精度要求不高，可将sims=500改为sims=200，大幅缩短单组分析时间
• 预合并数据：将df1/df2/df3合并为一个大表，避免循环中重复创建m.data
• 使用data.table：替换data.frame，提升行绑定和数据操作效率

内容的提问来源于stack exchange，提问作者freutopia