在R中使用具有bootstrap推断和插补估计的CMAverse的不同结果

disho6za  于 2023-04-18  发布在  Bootstrap
关注(0)|答案(2)|浏览(350)

在R中使用CMAverse,当我运行相同的模型两次时,我似乎不能得到相同的结果。示例代码:

# 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).

如何让R在重新运行相同的模型时给予相同的结果?

uidvcgyl

uidvcgyl1#

首先运行这个“set.seed(123)”,然后运行模型,你会得到相同的结果

lyr7nygr

lyr7nygr2#

我没有运行你的整个脚本,因为CMAverse不适用于我运行的R版本,但我猜你没有在每次运行代码时都运行set.seed()?如果你希望每次调用rnorm()时都有相同的值,你必须每次运行set.seed()。例如:

set.seed(1)
rnorm(10)
# [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684  0.4874291
# [8]  0.7383247  0.5757814 -0.3053884
rnorm(10)
# [1]  1.51178117  0.38984324 -0.62124058 -2.21469989  1.12493092 -0.04493361 -0.01619026
# [8]  0.94383621  0.82122120  0.59390132

将产生具有不同值的向量,而:

set.seed(1)
rnorm(10)
# [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684  0.4874291
# [8]  0.7383247  0.5757814 -0.3053884
set.seed(1)
rnorm(10)
# [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684  0.4874291
# [8]  0.7383247  0.5757814 -0.3053884

我会一遍又一遍地重复相同的向量。如果这不能解决你的问题,让我知道,我会进一步探讨你的问题。

相关问题