1.動機

Rの機械学習のフレームワークであるMLR３は、マイナーな存在ではあるが、実は隠れた優れものである（と思っている）。それをPythonの定評のある優れたフレームワークのscikit-learnと比較してみたい。
データ・セットはkaggleの糖尿病データを用いて、それぞれ予測モデルを組み立ててみる。
注意：データ・セットや、それぞれの機械学習のアルゴリズムの説明は省きます。
以下からデータセットはダウンロードしてください。
kaggle-dataset

データはダウンロードしてinputフォルダに入れておく。

本稿の目的は、RとPythonの機械学習フレームワークを用いて、下記４モデルの性能を比較することとする。

1. ロジスティック回帰
2. サポートベクターマシーン（SVM）
3. lightGBM（決定木のアンサンブル勾配ブースティング）
4. スタッキングアンサンブルモデル（上記1-3をベースに、final_estimatorをロジスティック回帰とするスタッキングアンサンブル）

なお、ハイパーパラメータのチューニングはlightGBMのみとし、ベイズ最適化を用いる

目的変数はdiabetesの糖尿病が発病するかどうかの分類予測モデルの作成とする。
目的関数は分類性能を示すROC_AUCにした。

1.1　使用するライブラリ

1.1.1　R

library(tidyverse)
library(mlr3)
library(mlr3learners)
library(mlr3extralearners)
library(mlr3viz)
library(mlr3pipelines)
library(bbotk)
library(mlr3mbo)
library(mlr3tuning)

1.1.2　Python

import pandas as pd
import numpy as np

from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
import optuna

from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import roc_auc_score
from sklearn.metrics import RocCurveDisplay

from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
import lightgbm as lgb
from sklearn.ensemble import StackingClassifier

from sklearn.model_selection import StratifiedKFold
from sklearn.compose import ColumnTransformer

2.DATA

2.1 R

diabetes_df = read_csv('input/diabetes_prediction_dataset.csv') %>%
  mutate(diabetes = as.factor(diabetes)) %>%
  mutate_if(is.character, as.factor)

hypertension：高血圧とheart_disease：心血管疾患は、数値であるがファクター化する

diabetes_df$hypertension=as.factor(diabetes_df$hypertension)
diabetes_df$heart_disease = as.factor(diabetes_df$heart_disease)

2.1.1 Taskの作成

MLR3では分類や回帰モデルによって、まずデータのタスクを作成する必要がある。

tsk_db =TaskClassif$new(id="diabetes", backend = diabetes_df,target = "diabetes" )

圧倒的に陽性データが少ない不均衡データなので、訓練時に偏りがないよう層化属性を付加する。

tsk_db$set_col_roles("diabetes", c("target", "stratum"))
tsk_db

#> <TaskClassif:diabetes> (100000 x 9)
#> * Target: diabetes
#> * Properties: twoclass, strata
#> * Features (8):
#>   - dbl (4): HbA1c_level, age, blood_glucose_level, bmi
#>   - fct (4): gender, heart_disease, hypertension, smoking_history
#> * Strata: diabetes

タスクを作成するとデータの概要が簡単に見ることができる。数値データ列、ファクター列がそれぞれ４列ある２クラスの分類モデルであることがわかる。またターゲットには層化属性がついていることも。

2.2 Python

hypertension：高血圧とheart_disease：心血管疾患を文字列に変換。

diabetes_df = pd.read_csv('./input/diabetes_prediction_dataset.csv')
diabetes_df['hypertension'] = diabetes_df['hypertension'].astype(str)
diabetes_df['heart_disease'] = diabetes_df['heart_disease'].astype(str)
diabetes_df.info()

#> <class 'pandas.core.frame.DataFrame'>
#> RangeIndex: 100000 entries, 0 to 99999
#> Data columns (total 9 columns):
#>  #   Column               Non-Null Count   Dtype 
#> ---  ------               --------------   ----- 
#>  0   gender               100000 non-null  object 
#>  1   age                  100000 non-null  float64
#>  2   hypertension         100000 non-null  object 
#>  3   heart_disease        100000 non-null  object 
#>  4   smoking_history      100000 non-null  object 
#>  5   bmi                  100000 non-null  float64
#>  6   HbA1c_level          100000 non-null  float64
#>  7   blood_glucose_level  100000 non-null  int64 
#>  8   diabetes             100000 non-null  int64 
#> dtypes: float64(3), int64(2), object(4)
#> memory usage: 6.9+ MB

3.モデルの準備とハイパーパラメータ探索

3.1　R

3.1.1　モデルの準備

まず学習器オブジェクト(learner)を作成する。予測値はすべて確率を吐き出すように設定。

lg_lrn  = lrn("classif.logistic",
              predict_type = "prob")
svm_lrn = lrn("classif.svm",
               predict_type = "prob")

lightGBMのハイパーパラメータのチューニング設定。

lgm_lrn = lrn("classif.lightgbm",
             num_leaves = to_tune(10,200), 
             learning_rate = 0.06150108, 
             num_iterations = 539,
             bagging_fraction = to_tune(0.5,1),
             feature_fraction = to_tune(0.5,1),
             min_data_in_leaf = to_tune(1, 20),
             predict_type = "prob"
)

ハイパーパラメータのチューニングは、num_leaves、bagging_fraction 、feature_fraction、min_data_in_leaf の４パラメータとした。

3.1.2　前処理の設定

MLR3では前処理をパイプライン処理を用いて自動化できる。そのためにpipelineオブジェクトを作成してから、learner化（学習器オブジェクト）する。
カテゴリー列についてラベルエンコーディングするpipelineを作成した。

# ラベルエンコーディングを行うPipelineを作成 
pl =  po("colapply", affect_columns=selector_type("factor"), applicator = as.numeric)

#pipelineオブジェクトをlearner化する
lg_lrn_pl = as_learner(pl %>>% lg_lrn)  
svm_lrn_pl = as_learner(pl %>>% svm_lrn)
lgm_lrn_pl = as_learner(pl %>>% lgm_lrn)
lg_lrn_pl$id="logistic"
svm_lrn_pl$id="svm"
lgm_lrn_pl$id="lightGBM"

Rでは、ファクター化すると自動的に水準が振られるため、ラベルエンコーディングはas.numeric（数値化）するだけで良い。

3.1.3　ベイズ最適化

lightGBMのハイパーパラメータチューニングをする。

bayesopt_ego = mlr_loop_functions$get("bayesopt_ego")
surrogate = srlrn(lrn("regr.km", covtype = "matern5_2",
  optim.method = "BFGS", control = list(trace = FALSE)))
acq_function = acqf("ei")
acq_optimizer = acqo(opt("nloptr", algorithm = "NLOPT_GN_ORIG_DIRECT"),
  terminator = trm("stagnation", iters = 100, threshold = 1e-5))

上記のベイズ最適化のためのガウス過程回帰の設定は細かくできる。基本的には上記設定から変更しなくても良いと思う。
5分割交差検証法で、50回イテレーションでチューニングする。

rsmp_cv = rsmp("cv", folds = 5)
msr_auc = msr("classif.auc")
tuner = tnr("mbo",
  loop_function = bayesopt_ego ,
  surrogate = surrogate,
  acq_function = acq_function,
  acq_optimizer = acq_optimizer)
instance = tune(
  tuner,
  tsk_db,
  lgm_lrn_pl,
  rsmp_cv,
  msr_auc,
  50
)

私のPCで約10分ほどかかった。結果は次のとおり。

instance$result_learner_param_vals

#> $colapply.applicator
#> function (x, ...)  .Primitive("as.double")
#> 
#> $colapply.affect_columns
#> selector_type("factor")
#> 
#> $classif.lightgbm.num_iterations
#> [1] 539
#> 
#> $classif.lightgbm.verbose
#> [1] -1
#> 
#> $classif.lightgbm.convert_categorical
#> [1] TRUE
#> 
#> $classif.lightgbm.learning_rate
#> [1] 0.06150108
#> 
#> $classif.lightgbm.num_threads
#> [1] 1
#> 
#> $classif.lightgbm.num_leaves
#> [1] 10
#> 
#> $classif.lightgbm.min_data_in_leaf
#> [1] 10
#> 
#> $classif.lightgbm.bagging_fraction
#> [1] 0.7328989
#> 
#> $classif.lightgbm.feature_fraction
#> [1] 0.5323333

autoplot(instance)

MLR3では、チューニング結果を簡単に図示化できる。num_leavesは明らかに小さい方が良い。その他は、ばらつきが大きくなっている。それではベイズ最適化の結果を代入する。

lgm_lrn_pl$param_set$set_values(
  classif.lightgbm.num_leaves = 10,
  classif.lightgbm.min_data_in_leaf = 10,
  classif.lightgbm.bagging_fraction =  0.7328989,
  classif.lightgbm.feature_fraction =  0.5323333,
  classif.lightgbm.learning_rate = 0.06150108,
  classif.lightgbm.num_iterations = 539
)

3.1.4　スタッキングアンサンブルモデルの作成

パイプライン処理を使ってスタッキングモデルを組み立てる。

stack_pl = ppl("stacking",base_learners = c(svm_lrn_pl,lgm_lrn_pl,lg_lrn_pl),
  super_learner = lrn("classif.logistic",predict_type = "prob")
)
stack_pl$plot(horizontal = TRUE)

次にpipelineオブジェクトをlearner化（学習器）しておく。
（これはしなくてもベンチマーク時などには自動的にlearner化（学習器）に変換されるみたいだが、pipelineオブジェクトと区別を意識するために必要）

stack_pl=as_learner(stack_pl)
stack_pl$id="stack"

3.2　Python

3.2.2　モデルの準備

# ロジスティック回帰モデルを作成
model1 = LogisticRegression()

# サポートベクターマシーン（SVM)モデルを作成
model2 = SVC(probability=True)

3.2.3　前処理とベイズ最適化

前処理の準備

# キャラクター列の名前を取得
label_cols = diabetes_df.select_dtypes(include=['object']).columns.tolist()

# 使用するデータをセットする
data = diabetes_df.copy()

層化属性の設定と訓練データとテストデータの分割。

# トレーニングとテストでのデータ分割 ハイパーパラメータの探索
X_train, X_test, y_train, y_test = train_test_split(
                                        data.drop('diabetes', axis = 'columns'),
                                        data['diabetes'],
                                        train_size   = 0.9,
                                        stratify = data['diabetes'],
                                        random_state = 123,
                                        shuffle      = True )

次にラベルエンコーディングを実行する。pipline処理でラベルエンコーディングを適用しようとしたが、エラーがでた（原因はよくわからない）。

# LabelEncoderによる変換
label_encoder = LabelEncoder()
# データの変換
for col in  label_cols:
    X_train[col] = label_encoder.fit_transform(X_train[col])
    X_test[col] = label_encoder.fit_transform(X_test[col])

ベイズ最適化はoptunaを用いる。

def objective(trial):
    objective     = 'binary'
    metric        = 'binary_logloss'
    boosting_type = 'gbdt'
    learning_rate = 0.06150108
    n_estimators =  539
    num_leaves    = trial.suggest_int('num_leaves', 10, 200)
    min_data_in_leaf = trial.suggest_int('min_data_in_leaf', 10, 20)
    feature_fraction = trial.suggest_uniform('feature_fraction', 0.5, 1.0)
    bagging_fraction = trial.suggest_uniform('bagging_fraction', 0.5, 1.0)
    clf = lgb.LGBMClassifier(
        objective     = objective,
        metric        = metric  ,
        boosting_type = boosting_type,
        learning_rate = learning_rate,
        num_leaves    = num_leaves ,
        min_data_in_leaf = min_data_in_leaf,
        feature_fraction = feature_fraction,
        n_estimators = n_estimators,
        bagging_fraction =bagging_fraction
    )
    # クロスバリデーションの設定
    cv = StratifiedKFold(n_splits=5, shuffle=True)
    score = cross_val_score(clf, X_train, y_train, n_jobs=-1, cv=cv, scoring='roc_auc')
    auc = score.mean()
    return auc

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=50)

print("Best trial:")
trial = study.best_trial
print(f"  Value: {trial.value}")
print("  Params: ")
for key, value in trial.params.items():
    print(f"    {key}: {value}")

#> Execution Time: 245m 0.4s

何故か245分もかかった。Rの10分と雲泥の差がある。ほぼ同じことをしているはずなのに。

# ベストパラメータ
study.best_trial.params

#> {'num_leaves': 10,
#>  'min_data_in_leaf': 13,
#>  'feature_fraction': 0.6180631243210197,
#>  'bagging_fraction': 0.541423615086831}

パラメータの傾向はRの結果とそう大きくは変わらない。

チューニング結果を入れる。

model3 = lgb.LGBMClassifier(
    objective     = 'binary',
    metric        = 'binary_logloss' ,
    boosting_type = 'gbdt',
    num_leaves = 10,
    min_data_in_leaf = 13,
    feature_fraction = 0.6180631243210197,
    bagging_fraction = 0.541423615086831  ,
    learning_rate = 0.06150108,
    n_estimators =  539   
)

テストデータで検証する。

_ = model3.fit(X_train,y_train)
pred = model3.predict_proba(X_test)[:, 1]

aucを求める。

roc_auc = roc_auc_score(y_test,pred)
roc_auc

#> 0.9791605914496947

4.モデル性能の比較

モデル間の性能を交差検証法により比較する

4.1　R

ベンチマーク用のグリッドサーチ空間を作成する。benchmark_gridを使って(タスク、学習器{ここでは４つ}、リサンプリングの方法：rsmp_cv = rsmp("cv", folds = 5) ）を順番に指定する。

design = benchmark_grid(tsk_db,c(lg_lrn_pl,svm_lrn_pl,lgm_lrn_pl,stack_pl), rsmp_cv)

ベンチマークを実行する。

bmr = benchmark(design)

上記の実行は34分ほどかかった。

aggregateメソッドを用いて、目的関数を適用する（結果のデーターフレームの内、必要な3列だけを表示）

bmr$aggregate(msr("classif.auc"))[, .(task_id, learner_id, classif.auc)]

結果がデータフレームの一覧表で表示された。しかし、これは交差検証した平均の結果なので、さらに結果の分布を調べる。

autoplot(bmr, measure = msr("classif.auc"))

サポートベクターマシーンの性能は低く、しかもばらつきが大きい。一番性能が良いのは、lightGBMであるが、スタッキングアンサンブルモデルのばらつきがより小さいため、この程度の性能差であれば、より信頼性の高いスタッキングモデルの方が良いとも判断できる。

4.2　Python

層化指定した交差検証法を指定する（MLR3の場合はタスク段階で指定しているため、交差検証法を用いるときには気にしなくて良い）

skf = StratifiedKFold(n_splits=5, shuffle=True)

交差検証するときに、データ・セットはターゲットと説明変数群に切り分ける必要があるので、元データを壊さないようにコピーする。
（参照セマンティクスを使用しているため。MLR3も同じ）

# 使用するデータをセットする
data = diabetes_df.copy()
data_main = data.drop('diabetes', axis = 'columns')
data_target = data['diabetes']

前処理を適用する。

# データの変換
for col in  label_cols:
    data_main[col] = label_encoder.fit_transform(data_main[col])

model1（ロジスティック回帰）のスコアを求める。
MLR3のbenchmark関数と良く似ているcross_val_score関数に、学習器、データ、ターゲット、リサンプリング、目的関数をセットする

score1 = cross_val_score(model1, data_main, data_target, n_jobs=-1, cv=skf, scoring='roc_auc')
auc1 = score1.mean()
auc1  

#> 0.9527627643844422

model2（サポートベクターマシーン）のスコアを求める。

score2 = cross_val_score(model2, data_main, data_target, n_jobs=-1, cv=skf, scoring='roc_auc')
auc2 = score2.mean()
auc2  

#> 0.9610031147540983

model3（lightGBM）のスコアを求める。

score3 = cross_val_score(model3, data_main, data_target, n_jobs=-1, cv=skf, scoring='roc_auc')
auc3 = score3.mean()
auc3  

#> 0.9790715429122467

最後にスタッキングアンサンブルモデルのスコアを求める。
StackingClassifier関数を用いる。

# stackingによるアンサンブル学習
from sklearn.ensemble import StackingClassifier
estimators = [
    ("logstic", model1),
    ("svm", model2),
    ("lightGBM", model3),
]
clf = StackingClassifier(
    estimators=estimators, final_estimator=LogisticRegression())
clf

model4（スタッキングモデル）のスコアを求める。

score4 = cross_val_score(clf, data_main, data_target, n_jobs=-1, cv=skf, scoring='roc_auc')
auc4 = score4.mean()
auc4 

#> 0.9789842140790743

5.まとめと考察

RとPythonによる４モデルの性能比較結果をまとめる。

model	R	Python
Logistic	0.9613	0.9528
SVM	0.9303	0.9610
lightGBM	0.9794	0.9791
Stacking ensemble	0.9748	0.9790

LogisticはRがPythonよりも少し良くなった。
Logisticは、GLMで統計的なモデルであり、理論的には同じ値になるはずであるが、Rが0.01ほど良い結果となった。これはクロスバリデーション時のデータ割当の違いによるブレかもしれない。

次にRのSVMが性能低くなった。これについては、クロスバリデーションによるブレの範囲を超えているようにも思える。RとPythonではデフォルトのハイパーパラメータが異なるとかアルゴリズムが異なるとかの違いがあるのだろうか？

lightGBMは、ハイパーパラメータの設定は両者で若干異なってはいたものの、同じ結果となった。

Stacking ensembleだが、上記の結果が異なるので、当然に違ってくるのだが、その差は小さくなっている。スタッキングすると汎化性能が上がる場合があるので、うまく内部で処理ができているのではないか。

最後に、RのMLR3は、Pythonのscikit-learnと比較してもテーブルデータであれば、勝るとも劣らない使い勝手と性能があると思いました。

特に視覚化に関してはscikit-learnよりはるかに使い勝手が良い。同じことをscikit-learnですると、コードがさらに数行程度以上必要になり面倒くさいが、MLR3ではautoplot関数一発である。

Rの中でも機械学習フレームワークはtidymodelがメインで、MLR3はまだまだマイナーな存在ですが、MLR3をもっと推していきたい（結論）。