miceの並列処理で時間短縮

「欠測処理でRのmiceパッケージのmice()を使うとして、処理時間の問題皆さんどうしてますかね。」と言う質問に、「多重代入法はM個の代入を独立して実行するのが望ましいので、マルチコアを使える演算環境ならば、並列処理することができます…Flexible Imputation of Missing Data第2版の5.5節Parallel computationが参考になると思います。」と言うリプライがあって、参照先を見ると

doParallelパッケージのcomplete関数とibind関数を組み合わせて何とかする
parlMICE関数を使う
micemdパッケージのpar.mice関数を使う（追記：同様に使えるfuturemice関数の方が使いやすいそうです）

方法が示唆されていたのですが、一番手軽そうなparlMICEを使う方法の検索したら最初に出てきた記事の内容が古くなっていたので、例をつくってみました。

library(mvtnorm)
library(mice)

# 練習用のデータセットをつくる
mdata <- with(list(n = 100,
    mu = rep(0, 4),
    sigma = matrix(rep(c(1.0, rep(0.5, 4)), 4)[1:16], 4, 4),
    beta = matrix(c(1, -2, 3, -4), 4, 1)), {
    data <- rmvnorm(n, mean = mu, sigma = sigma)
    colnames(data) <- paste("x", 1:ncol(data), sep = "")
    y <- -1 + data %*% beta + rnorm(n, sd = 2)
    nhanes <- ampute(data, prop = 0.8)
    cbind(y, nhanes$amp)
})

# 並行処理のための準備
# n.core: 平行処理で利用するコア数
# n.imp.core: コアごとの生成データセット数
# methodやpredictorMatrixはmice関数と同様（blocksは上手く動かない？）
imp <- parlmice(mdata, n.imp.core = 20, cluster.seed = 223, n.core = 3)
# もしくは
# imp <- futuremice(mdata, n.imp.core = 20, cluster.seed = 223, n.core = 3)
imp$method

# 欠損値補定
fit <- with(imp, lm(y ~ x1 + x2 + x3 + x4))
pool(fit)