Rの{future.apply}や{furrr}をつかってmultisessionで並列処理しようとしたときに、"object 'xxx' not found"（xxxはグローバル環境にあるオブジェクト名）とエラーがでてちょっと困ったので暫定的な解決策のメモ

要旨

• Rの{future.apply}や{furrr}をつかってmultisessionで並列処理しようとしたときに、object 'xxx' not found（xxxはグローバル環境にあるオブジェクト名）とエラーがでたときの暫定的な解決策について。
• {future.apply}ならfuture.globals引数に参照できなかったグローバル環境変数を名前をベクトルで与えるか、名前付きリスト形式でグローバル環境変数を与えるといい。
• {furrr}なら.options引数のfurrr_options()のglobals引数に参照できなかったグローバル環境変数を名前をベクトルで与えるか、名前付きリスト形式でグローバル環境変数を与えるといい。
• いずれの場合も、グローバル環境のオブジェクト名や名前付きリストをリストアップするのが面倒なら、ls(.GlobalEnv)やas.list(.GlobalEnv)（ls(globalenv())やas.list(globalenv()でも可）でグローバル環境変数名や名前付きリストを取得して引数に渡してもよい。

Rの{future.apply}や{furrr}をつかってmultisessionで並列処理しようとしたときに、object 'xxx' not found（xxxはグローバル環境にあるオブジェクト名）とエラーがでてちょっと困ったので暫定的な解決策のメモ。

{future.apply}や{furrr}はともに、最初にだまっていればグローバル環境変数を取り込んでくれているはずなのでこんなことはいらないはずだけど、条件を踏むとエラーが再現される。

単純な例を次にあげる。

full_fomula <- formula(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species)
full_model <- lm(full_fomula, data = iris)

として、irisデータの回帰モデルをfull_modelとしてオブジェクトとして持っておく、これの決定係数はsummary(full_model)$r.squaredとして参照できるが、ちょっとへんてこなことをするとsummary(lm(as.formula(full_model[["call"]]), data = iris))$r.squaredとして求める場面を考える。

as.fomula(full_model[["call"]])というのはfull_model[["call"]]に格納されているオブジェクトの式であるlm(formula = full_fomula, data = iris)を呼び出して、as.fomula()関数をつかって式部分だけ（full_fomula）を取り出しているという、ここだけならただただ回りくどいことをしているだけ。

たとえばこうすると次のように、データをブートストラップして決定係数を複数もとめて回るという処理が簡単にできる。

sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris[sample(1:nrow(iris), replace = TRUE), ]))$r.squared
})

##  [1] 0.8822899 0.8935367 0.8807054 0.8778784 0.8748053 0.8722960 0.8768117
##  [8] 0.8533426 0.8951152 0.8880692

もっと単純にブートストラップしないでただループする場合は次のような例となる。今回の例は一回グローバル環境変数をオブジェクト（今回の場合、full_fomula）を含む式の再評価を使うことがミソとなっている。

sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
})

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

これはうまくいく。同様に{purrr}のmap()をつかった場合もうまくいく。

# ベクトル形式で答えが欲しいのでmap_dblを使う
purrr::map_dbl(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
})

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

同じようなことを{future}のplan()でmultisessionを指定して、並列処理できるようにして、{future.apply}を使ってコードを書いてみる。

library(future)
plan("multisession")
library(future.apply)
future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
})

## Error in eval(x$formula): object 'full_fomula' not found

とエラーがでる。これは{furrr}を使った場合でも同様のエラーがでる。

library(furrr)
future_map_dbl(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
})

## Error:
## ℹ In index: 1.
## Caused by error:
## ! object 'full_fomula' not found

これはplan()がsequentialの場合はエラーが出ない。

plan("sequential")
future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
})

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

つまり、“multisession”で環境をコピーするときになにかエラーが起きている様子。

furrr_options()のhelpには次のように書いてあり、

Global variables

globals controls how globals are identified, similar to the globals argument of future::future(). Since all function calls use the same set of globals, furrr gathers globals upfront (once), which is more efficient than if it was done for each future independently.

とあり、グローバル環境変数はコピーされているはずなのだがよくわからない。なんか自分が勘違いしてる？

関数内にグローバル環境変数を用いることが悪いわけじゃなく、次のような例の場合はうまくいく。

plan("multisession")
future_sapply(1:10, function(x){
  summary(lm(full_model, data = iris))$r.squared
})

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model), data = iris))$r.squared
})

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

future_sapply(1:10, function(x){
  summary(lm(full_fomula, data = iris))$r.squared
})

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

ないと言われたら与えればいいので、任意のグローバル環境変数を与えて解決することとする。{future.apply}ならfuture.globals引数に、{furrr}なら.options引数のfurrr_options()のglobals引数に参照できなかったグローバル環境変数を名前をベクトルで与えるか、名前付きリスト形式でグローバル環境変数を与えるといい。

ここで単純にオブジェクトがないよ、と先の例で言われたfull_fomulaだけをオブジェクトとして与えると

plan("multisession")
future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
}, future.globals = "full_fomula")

## Error in as.formula(full_model[["call"]]): object 'full_model' not found

このようにエラーがでる。

そんなわけで使うオブジェクトを渡してやる必要があり、

plan("multisession")
future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
}, future.globals = c("full_fomula", "full_model"))

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

となるといい。

これは名前付きのリストでもいいので、

plan("multisession")
future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
}, future.globals = list("full_fomula" = full_fomula, "full_model" = full_model))

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

とかける。

オブジェクト2つくらいならそんなに面倒くさくないが、必要なオブジェクトが増えたりしたときや、コードを書き換えたときにここの部分も書き換えなければいけないのはちょっと手間に感じる。

そこで、ls(.GlobalEnv)やas.list(.GlobalEnv)でグローバル環境変数名や名前付きリストを取得して引数に渡すと楽になる（ls(globalenv())やas.list(globalenv()でも可）（必要のないオブジェクトもコピーしてしまうので無駄なものもあるが）。

plan("multisession")
future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
}, future.globals = ls(.GlobalEnv))

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

future_sapply(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
}, future.globals = as.list(globalenv()))

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

これは{furrr}では例えば次のようになる。

future_map_dbl(1:10, function(x){
  summary(lm(as.formula(full_model[["call"]]), data = iris))$r.squared
}, .options = furrr_options(globals = ls(.GlobalEnv)))

##  [1] 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123 0.8673123
##  [8] 0.8673123 0.8673123 0.8673123

たぶん、{globals}に起因する問題な気がするが、ちょっとまだ良くわからない。

手元にLinuxがないのでmulticoreでもこのような事態になるかは分からないが、仕組み上はならないんじゃないかと思っているが、どうなんだろう。

オブジェクトが入れ子になって再評価するのが引っ掛かっている気がする（eval()でも再現する、というか最初に気づいたのはeval()を使っていたときに気づいた）が、詳細はわからない。

上記の例だとlm()にfull_fomulraを入れて回せばいいんだけど、もう少し決定係数について汎用性を持たせて、glmでも、lmでも{MASS}のpolrでもとりあえず{performance}のr2()をつかって決定係数を求めて、ブートストラップして95％区間を求めたいと思った場合、オブジェクト側がglmなのか、lmなのかとりあえずどうでもいいからデータを入れ替えて再評価して回すというのを考えると次のようになる。

bootstrap_r2 <- function(model, alpha = 0.05, B = 1000){
  rsq <- sapply(1:B, function(x){
    data <- model[["model"]]
    model_call <- model
    model_call[["call"]][["data"]] <- data[sample(1:nrow(data), replace = TRUE), ]
    unname(performance::r2(eval(model_call[["call"]]))[[1]])
  })
  
  CI_lower <- quantile(rsq, probs = (alpha / 2))
  CI_upper <- quantile(rsq, probs = (1 - alpha / 2))
  
  c(CI_low = CI_lower, CI_upper = CI_upper)
}

1000回くらいのブートストラップだが、lmモデルくらいならサクッと計算できるが、polrモデルだと途端に重たくなる。

そこで次のように{future.apply}を使って問題解決したいと思うと上に書いたようにグローバル環境変数を読み込まなければいけない。また、polrをMASS::polrで読み込んでいないでlibrary(MASS)で読み込んでいる場合にはこれも行方不明になるので{future.apply}ならfuture.packagesに文字列でパッケージ名を指定してやる必要がある。{future}であれば、furrr_options()のpackages引数に指定してやる必要がある。これもわざわざ指定してやるのが手間なのでserch()で読み込んでいるパッケージリストを手に入れ引数に与えてやる。

bootstrap_r2_future <- function(model, alpha = 0.05, B = 1000){
  oplan <- future::plan("multisession")
  on.exit(future::plan(oplan), add = TRUE)
  
  packages <- search()
  packages <- packages[grep("package:", packages)] |> 
    gsub("package:", "", x = _)
  
  rsq <- future.apply::future_sapply(1:B, function(x){
    data <- model[["model"]]
    model_call <- model
    model_call[["call"]][["data"]] <- data[sample(1:nrow(data), replace = TRUE), ]
    unname(performance::r2(eval(model_call[["call"]]))[[1]])
  }, future.globals = ls(globalenv()), future.seed = TRUE, future.packages = packages)
  
  CI_lower <- quantile(rsq, probs = (alpha / 2))
  CI_upper <- quantile(rsq, probs = (1 - alpha / 2))
  
  c(CI_low = CI_lower, CI_upper = CI_upper)
}

これで分散処理できる。

適当にHairEyeColorでpolrモデルを作ってみる。

d <- data.frame(HairEyeColor)
d_name <- colnames(d[1:3])
d <- lapply(d_name, function(x){
  rep(d[[x]], d$Freq) 
}) |> cbind.data.frame()
colnames(d) <- d_name
library(MASS)
polr_fomula <- formula(Hair ~ Eye + Sex)
model_polr <- polr(data = d, formula = polr_fomula)

これで先のbootstrap_r2_future()で処理すると95％区間が求められる。

bootstrap_r2_future(model_polr)

##    CI_low.2.5% CI_upper.97.5% 
##      0.1466325      0.2642870

力技の解決策なので、もっとスマートな解決策があればいいと思うが……