mlr3を使ったパイプライン処理による前処理等の自動化

# １．はじめに

Rによる機械学習をする場合のフレームワークであるmlr3については、あまり紹介した投稿がないので作成します。
パイプライン処理とは、機械学習によるデータの学習前に、特徴量に施す前処理を組み合わせて処理を行う自動化処理のことで、Tidymodelではrecipeパッケージに相当します。
過去の投稿は次のとおり。
1. [mlr3を使ってみる](https://qiita.com/kenkenvw/items/c2c27d917a95b4a8b8a5)
2. [mlr3MBOを用いたベイズ最適化について（Bayesian Optimization）](https://qiita.com/kenkenvw/items/409d9f9a24d99f2f1cde)
3. [mlr3MBOを用いたベイズ最適化によるハイパーパラメータチューニング](https://qiita.com/kenkenvw/items/82efea5010cd374b58e9)
4. [mlr3を使った機械学習](https://qiita.com/kenkenvw/items/3d80d5d66c4baf2b6ed6)
5. [MLR３【R】とscikit-learn【Python】の機械学習フレームワークの比較](https://qiita.com/kenkenvw/items/bc2e78a7270ce09c4bc4)

**環境**
```
Ubuntu 22.04.5 LTS
R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
```

# ２．PipeOp: パイプライン演算子
パイプラインを実行するmlr3pipelinesの基本クラスはPipeOpで、「pipeline operator 」の略であり、関数は短縮形であるpo()が用いられます。mlr3ではpo()による処理を逐次的または並列的に組み合わせて、パイプライン処理として実行することができます。

```r
library(mlr3)
library(mlr3pipelines)
library(mlr3learners)
library(mlr3viz)
```

実行できる処理の一覧
```r
po()

#> <DictionaryPipeOp> with 72 stored values
#> Keys: adas, blsmote, boxcox, branch, chunk, classbalancing,
#>   classifavg, classweights, colapply, collapsefactors,
#>   colroles, copy, datefeatures, encode, encodeimpact,
#>   encodelmer, featureunion, filter, fixfactors, histbin,
#>   ica, imputeconstant, imputehist, imputelearner,
#>   imputemean, imputemedian, imputemode, imputeoor,
#>   imputesample, kernelpca, learner, learner_cv,
#>   learner_pi_cvplus, learner_quantiles, missind,
#>   modelmatrix, multiplicityexply, multiplicityimply, mutate,
#>   nearmiss, nmf, nop, ovrsplit, ovrunite, pca, proxy,
#>   quantilebin, randomprojection, randomresponse, regravg,
#>   removeconstants, renamecolumns, replicate, rowapply,
#>   scale, scalemaxabs, scalerange, select, smote, smotenc,
#>   spatialsign, subsample, targetinvert, targetmutate,
#>   targettrafoscalerange, textvectorizer, threshold, tomek,
#>   tunethreshold, unbranch, vtreat, yeojohnson
```
mlr3pipelinesに含まれるPipeOpsの最新リストとドキュメントはhttps://mlr-org.com/pipeops.html

## ２．１ DATA
アメリカのSaratoga Housing Dataを用います。
```r
library(mosaicData)
data("SaratogaHouses", package = "mosaicData")
saratoga_data = SaratogaHouses
saratoga_data$fireplaces=as.factor(saratoga_data$fireplaces)
```

まずタスクを作成します。
```r
task_house = TaskRegr$new(id="saratoga_houses",backend = saratoga_data,
                          target = "price")
task_house

#> <TaskRegr:saratoga_houses> (1728 x 16)
#> * Target: price
#> * Properties: -
#> * Features (15):
#>   - fct (7): centralAir, fireplaces, fuel, heating,
#>     newConstruction, sewer, waterfront
#>   - int (6): age, bedrooms, landValue, livingArea,
#>     pctCollege, rooms
#>   - dbl (2): bathrooms, lotSize
```

## ２．２　Graphオブジェクト
パイプラインは、Graph オブジェクトによって定義されます。Graphは エッジ を持ちます。
最初の前処理事例として、カテゴリカル特徴量であるfireplacesは5つの水準を持ちますが、3と4が極端に少ない水準数であり、特徴量としてはノイズ要因になるので、この場合は3水準にまとめたいと思います。
```r
table(saratoga_data$fireplaces)

#>  0   1   2   3   4 
#> 740 942  42   2   2 
```
po("collapsefactors")のGraphを用いて、自動的に0.01％以下の閾値で低水準化できます。
graphオブジェクトもLearnerと同様trainメソッドでタスクに適用します（この場合はタスクをリストで引き渡します）。
```r
ps = po("collapsefactors", no_collapse_above_prevalence = 0.01) 
result = ps$train(list(task_house))　
table(result$output$data(,"fireplaces"))

#> fireplaces
#>   0   1   2 
#> 740 942  46 
```
少数水準は2のクラスにまとめられ、42から46になりました。

次に、カテゴリカル特徴量であるファクター列には全てワンホットエンコーディングを適用します。
po("encode")を用います。自動的に適用列はファクター列になります。
```r
ps = po("encode") 
result = ps$train(list(task_house))
as.data.table(result$output$data())[1:5,8:11]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/624a7a53-9046-b12d-b6cd-53c10e3eab6a.png)

次に数値列にはスケール化を適用します。
po("scale")を用いますが、対象を数値列に適用します。
affect_columns()を使用しますが、留意点は関数しか受け付けないため、select_type()という関数を用います。
列のtypeは次のとおり指定します。
{'logical','integer','numeric','character','factor','ordered','POSIXct'}
```r
ps = po("scale", affect_columns= selector_type(c('integer','numeric')))
result = ps$train(list(task_house))
as.data.table(result$output$data())[1:5,1:5]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/d47bd457-96b7-5124-c67a-4b0ed1e1362d.png)

## ２．３　Graphの連結（シーケンシャルパイプライン）
Graphを構築するパイプラインは、%>>%-演算子（”二重矢印 “と読む）を使用してPipeOps同士を接続できます。
前処理として、最初から最後までデータに逐次的に適用されるシーケンシャルパイプラインを作成します。
```r
graph = po("scale", affect_columns=selector_type(c('integer','numeric'))) %>>% po("collapsefactors",no_collapse_above_prevalence = 0.01) %>>% po("encode")   
graph$plot(horizontal = TRUE)   
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/ec7dcead-40e5-4b4c-3729-9eaea20f5a3a.png)
適用する順番を間違えないように注意してください。上記の場合、スケール化を最後に持ってくると、ワンホット化した数字もスケール化されてしまいます。

データが意図どおり変換されているか確認します。
```r
result = graph$train(task_house) #連結後のgraphにはlistでなくそのままtaskを渡す
result[[1]]$data() #[[1]]でdata()メソッドを適用する階層を１段階飛ばしている。
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/0efa092e-5e3e-4fd0-8417-913b4d955282.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/51263de3-ebca-01cc-a607-d959638592ff.png)


## ２．４　ターゲットの変換
ターゲットのprice列は対数変換する必要がありますが、これもパイプライン処理に組み込みます。
ppl("targettrafo")を用いてターゲット列の変換をします。
targetmutateにターゲットに適用する関数を指定しますが、必ず元のスケールに戻す関数も指定します。

ppl()は、目的に合わせて既にpo()を複数組み合わせ作成されているgraphオブジェクトの集合です。
次のppl()がmlr3には用意されています。

> + ppl(”bagging”, graph)： mlr3pipelinesでは、バギングとは、異なるデータサンプルに対してグラフを複数回実行し、その結果を平均化する
> + ppl(”branch”, graphs)： PipeOpBranch を使用して、与えられたグラフから異なるパスブランチを作成する
> + ppl(”greplicate”, graph, n)： グラフ（単一のPipeOpでもよい）をn回複製するグラフを作成する。このパイプラインは、各PipeOpに接尾辞を追加することで、IDの衝突を回避する
> + ppl(”ovr”,graph)： 多クラス分類タスクを複数のバイナリ分類タスクに変換するためのOne-versus-rest分類。これらのタスクは、与えられたグラフによって評価されます。グラフは、学習器（または予測を出す学習器を含むパイプライン）である必要があります。バイナリ・タスクで行われた予測は、元のタスクに必要な多クラス予測に結合される
> + ppl(”robustify”)： 任意のタスクを任意のLearnerと互換性を持たせるための一般的な前処理ステップを実行します。
> + ppl(”stacking”, base_learners, super_learner)： スタッキングは、1つ以上のモデル（base_learners）からの予測を、後続のモデル（super_learner）の特徴として使用するプロセスです。
> + ppl(”targettrafo”, graph)： タスクの予測対象を変換するグラフを作成し、（targetmutate.trafo ハイパーパラメータに渡された関数を使用して）トレーニング中に適用された変換が、（targetmutate.inverter ハイパーパラメータに渡された関数を使用して）結果の予測で反転されることを保証する。


特にppl(”robustify”)は、機械学習で一般的に用いられる前処理がワンセットになっている優れものです。次の順番で前処理が逐次的に進みます。

>　ppl("robustify")とは、欠測値の推定とカテゴリカル値の特徴量エンコーディング等のシンプルで再利用可能なパイプラインを提供します。このパイプラインには、以下のPipeOpsが含まれています。
> 1. po("removeconstants") - 定数特徴量が削除されます。 
> 2. po("colapply") - 文字および順序特徴量がカテゴリとしてエンコードされ、日付/時刻特徴が数値としてエンコードされます。
> 3. po("imputehist") - 欠測値の数値の特徴量がヒストグラムサンプリングによって補完されます。
> 4. po("imputesample") - 欠測値の論理値を持つ特徴量が経験分布からのサンプリングによって補完されます。
> 5. po("missind") - 欠損補完したデータのインジケータとして、新しく論理変数を持つ列が追加されます。
> 6. po("imputeoor") - カテゴリ特徴量の欠落値が、新しいレベルでエンコードされます。 
> 7. po("fixfactors") - 予測およびトレーニング中に同じレベルが存在するように、カテゴリ特徴量のレベルを修正します (これは、空の要素レベルを削除することを含む場合があります)。
> 8. po("imputesample") - 前のステップでレベルを削除することで導入されたカテゴリ特徴量の欠損値を、経験分布からサンプリングして補完します。
> 9. po("collapsefactors") - max_cardinality引数（デフォルトは1000）で制御される水準レベル数未満になるまで、カテゴリ特徴量の水準レベルを（トレーニングデータで最もレアな因子から）折りたたみます。
> 10. po("encode") - カテゴリ特徴量がワンホットエンコーディングされます。
> 11. po("removeconstants") - 前のステップで作成された可能性のある定数特徴量が削除されます。
> ppl("robustify") は、オプションで task と learner を指定できます。
> 例えば、学習器が "missings "プロパティを持つ場合や、そもそもタスクに欠損値がない場合は、欠損値を代入しません。デフォルトでは、タスクと学習器が提供されない場合、防除的に設定されます：すべての欠損値を補完し、すべての特徴量タイプを数値に変換します。

Learnerには、線型回帰を指定しました。（LearnerもGraphと同じR6クラスなので、同じように扱えます）

```r
lrn_lm = lrn("regr.lm")
log_lm = ppl("targettrafo",
  graph = lrn_lm,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1))
```

これを、前述の前処理用graphと連結します。
```r
graph_lrn = graph %>>% log_lm 
graph_lrn$plot()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/f51fbe94-e9ba-da67-1239-af64a83c05f0.png)
このように、レゴブロックみたいに必要な処理をパイプラインとして連結することができます（自由度が高い）。
以上で、前処理と学習器が全て一体となったパイプライン処理が完成しました。

## ２．５　パイプラインによる学習と評価

データを訓練データとテストデータに４：１の割合で分割し、学習した訓練データを用いて、学習には用いられていないテストデータで評価します。
まず学習と評価に入る前に作成したgraphオブジェクトをLearnerオブジェクトに変換します。
（このままでは評価に必要なpredictメソッドが使えないため、Learnerクラスのメソッドの継承が必要です）
```r
splits = partition(task_house, ratio = 0.8) #データの分割
graph_lrn = as_learner(graph_lrn) #Learnerオブジェクトに変換
graph_lrn$train(task_house,row_ids = splits$train)
result_lm = graph_lrn$predict(task_house,row_ids = splits$test) #テストデータを適用
```
結果を図示化します。
```r
plot(result_lm)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/293176e8-7da8-c225-946c-430d4269aa85.png)
RMSEでテストデータを評価します。scoreメソッドに評価関数を指定します。
```r
result_lm$score(msr("regr.rmse"))

#> regr.rmse 
#>  80133.61 
```

## ２．６　ノンシーケンシャルパイプラインによるアンサンブルモデルの構築
mlr3では、逐次的なパイプライン処理だけでなく並列的な処理もパイプライン処理として構築可能です。
　上記の線型回帰モデル(lm)の他、決定木モデル(rpart)、サポートベクターマシーン(SVM)、k近傍法(knn)の機械学習モデルとこれらを混合したスタッキングアンサンブルモデルを構築します。
　スタッキングアンサンブル手法は、予測性能を大幅に向上させることができる、Kaggle等でよく使われるアンサンブル技法です。スタッキングの基本的な考え方は、複数のモデルからの予測を、後続のモデルの特徴として使用し、そのモデルでこれらの予測を組み合わせるというものです。ここでは後続モデルとして、ランダムフォレストモデルを用いた混合モデルとそれぞれの単モデルと比較してみます。
 + k近傍法(knn)
```r
lrn_knn = lrn("regr.kknn")
log_knn = ppl("targettrafo",
  graph = lrn_knn,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1))
graph_knn = graph %>>% log_knn
```

+ サポートベクターマシーン(SVM)
```r
lrn_svm = lrn("regr.svm")
log_svm = ppl("targettrafo",
  graph = lrn_svm,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1))
graph_svm = graph %>>% log_svm
```

+ 決定木モデル(rpart)
```r
lrn_rpt = lrn("regr.rpart")
log_rpt = ppl("targettrafo",
  graph = lrn_rpt,
  targetmutate.trafo = function(x) log(x+1),
  targetmutate.inverter = function(x) list(response = exp(x$response)-1))
graph_rpt = graph %>>% log_rpt
```
 
+ スタッキングアンサンブルモデル
まずpo("learner_cv"）で、それぞれアンサンブルモデルの学習器を定義付けします。
```r
po_lm_cv = po("learner_cv",
  learner = graph_lrn,
  resampling.folds = 2, id = "lm_cv"
)

po_knn_cv = po("learner_cv",
  learner = graph_knn,
  resampling.folds = 2, id = "knn_cv"
)

po_svm_cv = po("learner_cv",
  learner = graph_svm,
  resampling.folds = 2, id = "svm_cv"
)

po_rpt_cv = po("learner_cv",
  learner = graph_rpt,
  resampling.folds = 2, id = "rpt_cv"
)
```

次にgunion()で並列処理を行い、po("featureunion") で結果を結合します。そして、po("learner", lrn("regr.ranger")) でそれの予測モデルの結果を、ランダムフォレストモデルでアンサンブルします。
```r
vortex_lrn = gunion(list(
  po_lm_cv,
  po_knn_cv,
  po_svm_cv,
  po_rpt_cv
  )) %>>% po("featureunion") %>>% po("learner", lrn("regr.ranger")) 
vortex_lrn$plot()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/d7ac132e-b613-dea5-4884-cb0416bd607a.png)
gunion()で並列処理後、po("featureunion") で結果を結合していることがわかります。
また各並列処理の各過程（例えばlm_cv)においても、前述の前処理と学習器が全て一体となったパイプライン処理がそれぞれ実行されています。いわばパイプライン処理のパイプライン処理となっています。

## ２．7　アンサンブルモデルとの比較
mlr3では簡単にモデル間の比較が可能なベンチマーク関数を使います。比較したい項目をベクトル化して渡します。
```r
rsmp_cv3 = rsmp("cv", folds = 3) #3分割交差検証法
design = benchmark_grid(task_house, c(vortex_lrn,graph_lrn,graph_knn,graph_rpt,graph_svm), rsmp_cv3)
bmr = benchmark(design)
```

結果の図示化はautoplot()一発で可能です。autoplot()には、指標となる目的関数だけを指定します。
```r
autoplot(bmr,measure = msr("regr.rmse"))
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/bf08f6e3-bbad-8dcd-f247-7a8e55ce6bf6.png)
左側から、アンサンブルモデル、線型回帰モデル、k近傍法モデル、決定木モデル、SVMモデルとなっています。
```r
bmr$aggregate(measures =msr("regr.rmse"))[,.(learner_id,regr.rmse)]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161102/b86b1063-aedb-7a67-e996-db68a285d05a.png)
目的関数はRMSEなので値が低いほど良く、アンサンブルモデルが一番良い性能を示しました。
また、スタッキングモデルには、あらかじめppl(”stacking”, base_learners, super_learner)が用意されているので、簡単に実装できます。

# ３　まとめ
　このようにmlr3ではTidymodelと比べてもパイプライン処理の自由度が高く、簡単に自動化ができます。
また、mlr3は完全にオブジェクト指向の考え方に従っており、従来のRに慣れた人には使いにくいかもしれませんが、機械学習のフレームワークとして、Tidymodelの他に、このMlr3も十分に選択肢に入ってくると思っています。
Pythonに近い感覚でオブジェクト指向のmlr3は使えますし、Rなのでtidy文法も混ぜられます。個人的にはRの文法はごった煮感が強く、変態的な言語だとは思っていますが、そこが好きな言語になっています。

# ４　Enjoy