r-dopar中代码的过程优化

Question

我正在尝试优化我的代码以多次运行glms，我想利用并行化，无论是使用foreach还是其他更有效的方式。

正如你所看到的；CCD_ 2循环需要大约800秒来运行270000 glms；而CCD_ 3和CCD_。

谢谢你的帮助。

Jinesh

library(data.table)
library(parallel)
library(doParallel)
library(foreach)
scen_bin <- expand.grid(n = c(10, 20, 30), rate1 = c(0.1, 0.2, 0.3),
rate2 = c(0.5, 0.6, 0.9))
rep <- 10000
scen_sims <- rbindlist(replicate(rep, scen_bin, simplify = FALSE),
idcol = TRUE)
scen_sims[, `:=`(glm, list(c(1L, 2L)))]
for (i in 1:270000) {
set(scen_sims, i, 8L, list(glm(formula = c(rbinom(scen_sims$drug[i], 1L, scen_sims$Treatment_Rates[i]),
rbinom(scen_sims$control[i], 1L, scen_sims$Comparator_Rates[i])) ~ factor(c(rep("Trt",
scen_sims$drug[i]), rep("Cont", scen_sims$control[i]))), family = "binomial")))
}
split_scen_sims <- split(scen_sims, seq(1, 270000, length.out = 1000))

jh <- foreach(x = 1:1000, .packages = c("data.table")) %dopar% {
jh <- split_scen_sims[[x]]
for (i in 1:270000) {
set(jh, i, 8L, list(glm(formula = c(rbinom(jh$n[i], 1L, jh$rate1[i]), rbinom(jh$n[i],
1L, jh$rate1[i])) ~ factor(c(rep("Trt", jh$n[i]), rep("Cont", jh$n[i]))),
family = "binomial")))
}
return(jh)
}

Answer 1

首先要注意的是，在循环中使用提取函数$会使其性能不佳。最好是1(生成一个函数，然后2(使用常规的data.table调用。

fx_make_glm = function(drug, treat_rate, control, Comparator_Rates){
glm(formula = c(rbinom(drug, 1L, treat_rate),
rbinom(control, 1L, Comparator_Rates)) ~
factor(c(rep("Trt", drug), rep("Cont", control))), 
family = "binomial")
}

这将大大简化其余部分——我将使用Map，它将循环通过感兴趣的变量的每个元素：

scen_sims[, glm := list(Map(fx_make_glm, n, rate1, n, rate2))]

不幸的是，这仍然没有提供理想的性能：(

Unit: seconds
expr  min   lq mean median   uq  max neval
OP_loop 3.01 3.21 3.21   3.22 3.26 3.36     5
map_call 2.64 2.89 2.90   2.92 2.96 3.08     5

我选择的并行包是future.apply——只需将future_放在*apply系列前面，就可以进行并行评估：

library(future.apply)
plan(multiprocess)
system.time({
scen_sims[, glm := list(future_Map(fx_make_glm, n, rate1, n, rate2))]
})
user  system elapsed 
1.22    0.13    3.22 
## truncated the microbenchmark call
Unit: seconds
expr  min   lq mean median   uq  max neval
OP_loop 2.93 2.98 3.08   3.00 3.18 3.32     5
map_call 2.65 2.70 2.94   2.89 3.18 3.25     5
future_map_call 2.84 3.24 3.37   3.43 3.49 3.85     5

我使用的是2核/4线程的Windows。如果我在Linux上，我会尝试plan(multicore)，看看分叉进程是否更有效率。

数据生成：

library(data.table)
## generate data
scen_bin <- expand.grid(n = c(10, 20, 30), rate1 = c(0.1, 0.2, 0.3),
rate2 = c(0.5, 0.6, 0.9))
rep <- 50L
scen_sims <- rbindlist(replicate(rep, scen_bin, simplify = FALSE),
idcol = TRUE)
scen_sims[, `:=`(glm, list(c(1L, 2L)))]