R中CMAverse用Bootstrap推断与插补估计重复运行结果不一致的解决办法

在R语言中使用CMAverse包，采用imputation（插补）估计方法与bootstrap（自助法）推断进行因果中介分析时，重复运行同一模型却得到不同结果。示例代码如下：

# Load packages and set seed.
library(CMAverse)
library(tidyverse)
library(magrittr)
library(janitor)
set.seed(1)

# Simulate data containing a continuous baseline confounder C1, a binary baseline confounder C2, a binary exposure A, a binary mediator M and a binary outcome Y. 
expit <- function(x) exp(x)/(1+exp(x))
n <- 10000
C1 <- rnorm(n, mean = 1, sd = 0.1)
C2 <- rbinom(n, 1, 0.6)
A <- rbinom(n, 1, expit(0.2 + 0.5*C1 + 0.1*C2))
M <- rbinom(n, 1, expit(1 + 2*A + 1.5*C1 + 0.8*C2))
Y <- rbinom(n, 1, expit(-3 - 0.4*A - 1.2*M + 0.5*A*M + 0.3*C1 - 0.6*C2))
data <- data.frame(A, M, Y, C1, C2)

# Run causal mediation analysis.
model <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
                               mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
                               mreg = list("logistic"), yreg = "logistic",
                               astar = 0, a = 1, mval = list(1), yval=1,
                               estimation = "imputation", inference = "bootstrap", nboot = 10)

# Get the summary of the model.
summary <- model %>% summary() 
summary$summarydf %>% clean_names() %>% mutate_at(vars(estimate, x95_percent_cil, x95_percent_ciu), ~format(round(., digits=2), nsmall=2, trim=TRUE)) %>% mutate(estimate=paste0(estimate, " (", x95_percent_cil, "-", x95_percent_ciu, ")")) %>% select(-x95_percent_cil, -x95_percent_ciu) %>% mutate(p_val=format(round(p_val, digits=3), nsmall=3))
 # The pure natural direct effect is 1.42 (1.27-1.87) and the pure natural indirect effect is 0.92 (0.89-0.98).

# Rerun the causal mediation analysis.
model <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
               mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
               mreg = list("logistic"), yreg = "logistic",
               astar = 0, a = 1, mval = list(1), yval=1,
               estimation = "imputation", inference = "bootstrap", nboot = 10)

# Get the summary of the model.
summary <- model %>% summary() 
summary$summarydf %>% clean_names() %>% mutate_at(vars(estimate, x95_percent_cil, x95_percent_ciu), ~format(round(., digits=2), nsmall=2, trim=TRUE)) %>% mutate(estimate=paste0(estimate, " (", x95_percent_cil, "-", x95_percent_ciu, ")")) %>% select(-x95_percent_cil, -x95_percent_ciu) %>% mutate(p_val=format(round(p_val, digits=3), nsmall=3))
 # The pure natural direct effect is 1.43 (1.08-1.75) and the pure natural indirect effect is 0.91 (0.85-0.99).

问题根源在于bootstrap自助法和imputation插补过程都会生成随机数，仅在脚本开头设置一次set.seed()无法控制后续重复运行时的随机过程。要实现重复运行结果一致，有两种可行方法：

• 方法一：每次调用cmest()前显式设置相同随机种子
  在两次运行cmest()的代码前都添加相同的set.seed()语句，确保每次bootstrap抽样和插补的随机逻辑完全一致。修改后的关键代码片段如下：

  # 第一次运行模型前设置种子
  set.seed(1)
  model <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
                               mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
                               mreg = list("logistic"), yreg = "logistic",
                               astar = 0, a = 1, mval = list(1), yval=1,
                               estimation = "imputation", inference = "bootstrap", nboot = 10)
  
  # 第二次运行模型前再次设置相同种子
  set.seed(1)
  model <- cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
               mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
               mreg = list("logistic"), yreg = "logistic",
               astar = 0, a = 1, mval = list(1), yval=1,
               estimation = "imputation", inference = "bootstrap", nboot = 10)
  

• 方法二：使用withr::with_seed()锁定单次函数调用的随机环境
  若不想全局重复设置种子，可以用withr包的with_seed()函数，将每次cmest()调用包裹起来，限定随机种子仅对该次调用生效，避免影响其他代码的随机过程：

  library(withr)
  
  # 第一次运行
  model <- with_seed(1, cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
                               mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
                               mreg = list("logistic"), yreg = "logistic",
                               astar = 0, a = 1, mval = list(1), yval=1,
                               estimation = "imputation", inference = "bootstrap", nboot = 10))
  
  # 第二次运行
  model <- with_seed(1, cmest(data = data, model = "rb", outcome = "Y", exposure = "A",
               mediator = "M", basec = c("C1", "C2"), EMint = TRUE,
               mreg = list("logistic"), yreg = "logistic",
               astar = 0, a = 1, mval = list(1), yval=1,
               estimation = "imputation", inference = "bootstrap", nboot = 10))
  

修改后，两次运行cmest()得到的因果中介分析结果将完全一致。

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