Rで何かをしたり、読書をするブログ

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

UCI Machine Learning Repository の Chess (King-Rook vs. King-Pawn) のデータの分析3 - ランダムフォレストによる分類。素晴らしい性能

www.crosshyou.info

の続きです。前回はglmnetパッケージを使い、Elastic-Netで分類してみました。ROCのAUCが0.996と非常に素晴らしい値でした。正直、Elastic-Netでここまで高精度の分類器が作れたので、もういいかな、と思いましたが、今回は、tidymodelsパッケージを使い、rangerエンジンでランダムフォレストモデルを試してみます。

まず、トレーニング用のデータとテスト用のデータにわけます。

モデルを作ります。rand_forest()で、set_engine()でrangerを指定します。特徴量の重要度も見たいので、importance = "impurity"を加えておきます。mtryとmin_nはあとでチューニングできるように、ここでは、tune()としておきます。

レシピを作ります。

step_dummy()でファクター型の変数を0/1のダミー変数にして、step_zv()でゼロ分散の特徴量を除外するというレシピです。

ワークフローは、モデルとレシピを合体させただけです。

fold_cv()関数でクロスバリデーションの各foldを作ります。

チューニング・グリッドを作ります。grid_regular関数を使いました。

tune_grid()関数でチューニングを実行します。結構、時間かかりました。私のLaptop(プロセッサ:Intel(R) Core(TM) i5-1035G1 CPU @ 1.00GHz (1.19 GHz) / 実装 RAM:8.00 GB (7.71 GB 使用可能))だと、5分かかりました。

最適なパラメータを確認します。

mtry = 20, min=n = 2 が最適パラメータでした。この最適なパラメータで、最終ワークフローを作成し、train_dataで学習します。

vipパッケージのvip関数で、重要な変数を調べてみます。

X21, X10, X33が重要度が大きいですね。

前回のElastic-Netのモデルでは、どの変数が重要だったか、みてみましょう。

X21, X10, X33とランダムフォレストと同じ順番で重要となっていますね。

テスト用データで予測しましょう。

ROCのAUCを計算します。roc_auc()関数です。

おお!1.000です!完璧ですね。

混合行列をみてみます。

いや~!凄いですね。7個しか間違えてないです。ランダムフォレスト、素晴らしいですね。Sensitivityとか、Specificityを計算してみようと思います。

カスタム関数を作って、一括計算しました。

どれも素晴らしい性能を表していますね。正解率は、99.3%、F1-Scoreが99.3%という精度です。

今回は以上です。

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# トレーニング用のデータとテスト用のデータに分ける
train_data <- df[train, ]
test_data <- df[-train, ]
#
# モデルの作成
rf_model <- rand_forest(
  mtry = tune(),
  min_n = tune(),
  trees = 500
) |>
  set_engine("ranger", importance = "impurity") |>
  set_mode("classification")
#
# レシピの作成
rf_recipe <- recipe(target ~ ., data = train_data) |>
  step_dummy(all_nominal_predictors()) |>
  step_zv(all_predictors())
#
# ワークフローの作成
rf_wf <- workflow() |>
  add_model(rf_model) |>
  add_recipe(rf_recipe)
#
# クロスバリデーションの設定
set.seed(123)
folds <- vfold_cv(train_data, v = 10, strata = target)
#
# チューニング・グリッドを作成
rf_grid <- grid_regular(
  mtry(range = c(2, 20)),
  min_n(range = c(2, 20)),
  levels = 5
)
#
# チューニングの実行
set.seed(123)
rf_tuned <- tune_grid(
  rf_wf,
  resamples = folds,
  grid = rf_grid,
  metrics = metric_set(roc_auc, accuracy)
)
#
# 最適なパラメータの確認
best_params <- select_best(rf_tuned, metric = "roc_auc")
best_params
#
# 最終ワークフローと最終学習モデル
final_rf_wf <- finalize_workflow(rf_wf, best_params)
final_rf_fit <- fit(final_rf_wf, data = train_data)
#
# 重要な変数の表示
library(vip)
final_rf_fit |>
  extract_fit_parsnip() |>
  vip(num_features = 20)
#
# Elastic-Netの変数の重要度
coef_abs <- coef(c0.8, s = "lambda.1se") |> abs() |> round(4)
coef_abs[rev(order(coef_abs[ , 1])), ]
#
# テスト用データで予測
rf_pred <- predict(final_rf_fit, test_data, type = "prob") |>
  bind_cols(predict(final_rf_fit, test_data)) |>
  bind_cols(test_data |> select(target))
rf_pred
#
# ROCのAUC
roc_auc(rf_pred, truth = target, .pred_nowin)
#
# 混合行列
conf_mat(rf_pred, truth = target, estimate = .pred_class)
#
# Sensitivityなどを計算する
conf_mat_metrics <- function(TN, FN, FP, TP) {
  
  Accuracy <- (TN + TP) / (TN + FN + FP + TP)
  Sensitivity <- TP / (TP + FN)
  Specificity <- TN / (TN + FP)
  Precision <- TP / (TP + FP)
  F1_Score <- 2 * (TP / (TP + FP)) * (TP / (TP + FN)) / 
    *1 + (TP / (TP + FN)))
  
  tibble(
    metric = c("Accuracy", "Sensitivity", "Specitificty",
               "Precision", "F1_Score"),
    value = c(Accuracy, Sensitivity, Specificity, Precision, F1_Score)
  )
  
}
conf_mat_metrics(449, 1, 6, 503)
#

 

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Beautiful landscape of 300 years old natural green grass field and higher mountains and close up of red Oleander flowers and trees. Photo です。)

 

*1:TP / (TP + FP

UCI Machine Learning Repository の Chess (King-Rook vs. King-Pawn) のデータの分析2 - glmnetのElastic-Netで分類してみる。ROCのAUCが0.996と高性能でした。

www.crosshyou.info

の続きです。前回はデータをRに読み込ませて前処理をしました。

今回は実際に分類器を作ります。glmnetパッケージでElastic-Netでやってみようと思います。

glmnetパッケージを読み込みます。

glmnetは特徴量はマトリックス型、targetはベクトル型にする必要があるので、dfをxとyにわけます。

トレーニング用のインデックスを作成します。

alpha = 0.2, alpha = 0.5, alpha = 0.8 の3つのalpha水準でモデルを作りました。

この3つの中で、どれが一番よかったか調べます。

alpha = 0.8 の cvm が0.182で一番小さいので、c0.8を使って予測しましょう。

その前にこのモデルは、どんな係数か確認します。

alpha = 0.8 なので、LASSOよりのElastic-Netですが、あんまり特徴量の削減はされていないですね。

predict()関数で予測します。

混合行列をtable()関数で作ります。

お!結構あたっていますね。TN(True Negative), FN(False Negative), FP(False Negative), TP(True Positive)をあてはめます。

正解率(Accuracy)を計算します。

97.2%です。たいしたもんですね。

Sensitivityを計算します。

96.8%です。

Specitifityを計算します。

97.6%です。

Precisionを計算します。

97.8%です。

F-1 Scoreを計算します。

97.3%です。

どうでしょうか?結構いい分類結果だと思います。

glmnetパッケージには、assess.glmnet()関数、confusion.glmnet()関数がありますので、それも使ってみます。

$classが0.02815433でした。1 - 0.02815433 = 0.9718457 で先ほどのAccuracyと一致します。

$aucが0.9959663です。これはほぼ完璧レベルですね。

confusion.glmnet()関数も実行してみます。

confusion.glmnet()関数は、混合行列を作り、Accuracy(Percent Correct)を計算する関数でした。

roc.glmnet()関数を使ってROC曲線も描けます。

今回は以上です。AUCが0.996ですから、とても素晴らしいROC曲線ですね。

次回は

 

www.crosshyou.info

です。

はじめから読むには、

www.crosshyou.info

です。

今回のコートは以下になります。

#
# まず、tidymodelsを使わないElastic-Net
#
library(glmnet)
#
# はじめに、dfをx, y にわける
x <- df |>
  select(-target) |> 
  makeX()
#
y <- df$target
#
# トレーニング用のインデックス
set.seed(1919)
train <- sample(1:length(y), 0.7 * length(y), replace = FALSE) 
#
# alpha = 0.2
set.seed(1)
c0.2 <- cv.glmnet(x[train, ], y[train], family = "binomial",
                  alpha = 0.2)
#
# alpha = 0.5
set.seed(1)
c0.5 <- cv.glmnet(x[train, ], y[train], family = "binomial",
                  alpha = 0.5)
#
# alpha = 0.8
set.seed(1)
c0.8 <- cv.glmnet(x[train, ], y[train], family = "binomial",
                  alpha = 0.8)
#
# c0.2, c0.5, c0.8の中でどれが一番か調べます
tibble(
  cvm = c0.2$cvm,
  lambda = c0.2$lambda
) |> 
  filter(lambda == c0.2$lambda.1se)
#
tibble(
  cvm = c0.5$cvm,
  lambda = c0.5$lambda
) |> 
  filter(lambda == c0.5$lambda.1se)
#
tibble(
  cvm = c0.8$cvm,
  lambda = c0.8$lambda
) |> 
  filter(lambda == c0.8$lambda.1se)
#
# c0.8が一番でした。
# c0.8(lambda.1se)の係数を確認します。
coef(c0.8, s = "lambda.1se")
#
# c0.8(lambda.1se)で予測してみます。
c0.8_pred <- predict(c0.8, newx = x[-train, ], s = "lambda.1se",
                     type = "class")
head(c0.8_pred)
#
# 実際のyとの比較(混合行列)
table(c0.8_pred, y[-train])
#
# TN, FN, FP, TPを当てはめる
TN <- 444; FN <- 16; FP <- 11; TP <- 488
#
# Accuracy(正解率): 
# 全体のうち、正しく予測できた割合
(TN + TP) / (TN + FN + FP + TP)
#
# Sensitivity(Recall, True Positive Rate):
# 実際にwonのものを正しくwonとできた割合
TP / (TP + FN)
#
# Specificity(True Negative Rate):
# 実際にnowinのものを正しくnowinとできた割合
TN / (TN + FP)
#
# Precision(適合率):
# wonと予測したもののうち、実際にwonだった割合
TP / (TP + FP)
#
# F-1 Score:
# PrecisionとRecallの調和平均
# 予測が当たっているか(Precision), 見逃していないか(Recall)の両方を同時に評価
2 * (TP / (TP + FP)) * (TP / (TP + FN)) / 
  *1 + (TP / (TP + FN)))
#
# assess.glmnet()
assess.glmnet(c0.8, newx = x[-train, ], newy = y[-train],
              family = "binomial", s = "lambda.1se")
#
# confusion.glmnet()
confusion.glmnet(c0.8, newx = x[-train, ], newy = y[-train],
                 family = "binomial", s = "lambda.1se")
#

#
# roc.glmnet()でROC曲線
roc <- roc.glmnet(c0.8, newx = x[-train, ], newy = y[-train], 
                  family = "binomial", s = "lambda.min")
#
plot(roc$FPR, roc$TPR, type = "l", 
     xlab = "False Positive Rate", 
     ylab = "True Positive Rate")
abline(0, 1, lty = 2, col = "red")

 

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Natural landscape of rice field, close up of purple violet flowers, there is a white cloud in the high blue sky, photo です。)

 

*1:TP / (TP + FP

UCI Machine Learning Repository の Chess (King-Rook vs. King-Pawn) のデータの分析1 - データの読み込みと前処理

今回からしばらくは、UCI Machine Learning Repository の Chess () のデータを使って分類問題を練習してみたいと思います。

Chess (King-Rook vs. King-Pawn) - UCI Machine Learning Repository

Shapiro, A. (1983). Chess (King-Rook vs. King-Pawn) [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C5DK5C.

"King-Rook vs. King-Pawn on A7", KEKPA7と省略されるそうです。私はチェスは詳しくないので、よくわかりません。

Copilotに聞いてみました。

####ここからCopilotの回答####

KRKPA7(King–Rook vs. King–Pawn on A7)とは、
チェスの特定のエンドゲーム局面を指す専門用語で、白がキング+ルーク、黒がキング+a7のポーンという構成の終盤を意味します。


📝 KRKPA7 のポイント

● どんな局面?

  • 黒のポーンが a7 にあり、あと1手で昇格(クイーン化) できる非常に危険な状況。
  • 白は キング+ルーク でこのポーンの昇格を阻止しつつ、勝ち切れるかどうかがテーマ。

● なぜ有名なのか?

  • このエンドゲームは、機械学習の研究データセットとして非常に有名で、UCI Machine Learning Repository に「KRKPA7」として登録されています。
  • 3196局面が収録され、各局面が
    • White-can-win(白が勝てる)
    • White-cannot-win(白が勝てない)
      の2クラスに分類されています。

● 研究での利用

  • 1980年代から、知識表現・機械学習・概念学習の研究で頻繁に使われてきた古典的データセットです。

🎯 まとめ

KRKPA7 は、
「白:キング+ルーク vs 黒:キング+a7ポーン」
という特定のチェス終盤を扱った、歴史的にも重要な研究用データセットの略称です。

チェスの実戦でも理論的にも興味深いテーマで、
「ポーンが昇格寸前のとき、ルーク側は勝てるのか?」
という問題を体系的に分析するために作られました。


####ここまでCopilotの回答####

ということだそうです。

Rで分析してみようと思います。tidyverseとtidymodelsパッケージを読み込んでおきます。

データを読み込みます。

データはこういうテキストファイルでした。

glimpse()関数でデータを確認します。

X37が目的変数です。TRUE/FALSEのロジカル型の変数が多いですね。

X37をtargetと名前を変えて、ファクター型にしました。頻度を確認します。

勝ったゲームが1669で、負けたゲームが1527でした。勝率は52.2%です。

TRUE/FALSEのロジカル型の変数は、1/0の数値型に、文字列型の変数はファクター型なおします。

summary()関数でデータの確認をしてみましょう。

すべての変数が0/1の数値型、またはファクター型に変わっています。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

今回のコードは以下になります。

#
# tidyverseの読み込み
library(tidyverse)
#
# tidymodelsの読み込み
library(tidymodels)
tidymodels_prefer()
#
# データの読み込み
df <- read_delim("kr-vs-kp.txt", col_names = FALSE, delim = ",")
#
# データの確認
glimpse(df)
#
# X37をtargetに名称を変更し、ファクター型にする
df <- df |> 
  rename(target = X37) |> 
  mutate(target = as.factor(target))
#
#
# targetの頻度
df |> 
  count(target) |> 
  mutate(prop = n / sum(n))
#
# logicalはnumericに、characterはfactorに
df <- df |> 
  mutate(across(where(is.logical), as.numeric)) |> 
  mutate(across(where(is.character), as.factor))
#
# データの確認
summary(df)
#

 

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Landscape of natural green grass field, there are a few small ponds, close up of blue Begonia flowers and yellow Begonia flowers, photo. です。)

 

都道府県別の商業動態統計調査のデータの分析9 - Elastic-Net(ペナルティ付き線形モデル)をパワーアップ

www.crosshyou.info

の続きです。前回は線形モデル、ペナルティ付き線形モデル、決定木モデル、ランダムフォレストモデル、サポートベクター回帰(SVR)モデル、k-NN回帰モデル、ニューラルネットワークモデルの7つのモデルの性能評価をしました。その結果、ランダムフォレストモデルが一番の性能でした。

今回は、ペナルティ付き線形モデルの性能をもう少し上げることを考えましょう。

前回までの線形モデルは、target = 切片 + beta1 x 特徴量1 + beta2 x 特徴量2 + ... betaP x 特徴量P のように、特徴量どうしの相互作用や特徴量の2乗項、3乗項を考慮していないものでした。

なので、相互作用や2乗項、3乗項を加えてみて、性能が上がるかを確認してみます。はじめのどのような相互作用をいれたらいいかを確認したいので、線形モデルのモデルのサマリーをみてみます。

v1, v5, v6が統計的に有意な変数ですね。なので、この3つの相互作用と全てのの特徴量の2乗項を加えてみます。

レシピを新しくしたので、ワークフローも新しくします。

チューニング・グリッドも新しくしました。pennaltyの刻みを多くしました。

このチューニング・グリッドでチューニングします。

最適なパラメータを確認します。

この最適なパラメータで最終ワークフローを作ります。

この最終ワークフローで学習します。

predict()関数で予測します。

この新しい予測を従来の結果のデータフレームに追加します。

それでは、結果をみてみましょう。まずは、RMSEからです。

glmnet2が新しいペナルティ付き線形モデルの予測です。2.33で、ランダムフォレストモデル(ranger)を上回り、2位になりました。

MAEを見てみます。

MAEは見事に一位になりました。

R-squaredはどうでしょうか?

R-squaredは4位でした。それでもパワーアップ前のglmnetよりはよくなっています。

そして、総合結果の発表です。

glmnet2が2位になりました。パワーアップの効果がでましたね。

最後に、新しいペナルティ付き線形モデルの係数を確認します。

extract_fit_engine()関数でモデルを取り出し、coef()関数で係数を確認しました。

target = 切片 + beta1 * v1 + beta2 * v3 + beta3 * v6 + beta4 * v1:v5 + beta6 * v3:v6というモデルでした。v1とv5, v3とv6の相互作用がありました。単純な線形モデルに相互作用や2乗項を加えて、Elastic-Net(ペナルティ付き線形モデル)のパワーアップができました。

今回は以上です。 

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# 線形モデルのモデルを見る
extract_fit_engine(lm_fit) |> 
  summary()
#
# v1, v5, v6との相互作用と各特徴量の2乗項を加えたレシピ
glmnet_rec2 <- recipe(target ~ ., data = train_data) |> 
  # 相互作用
  step_interact(
    terms = ~ v1:v2 + v1:v3 + v1:v4 + v1:v5 + v1:v6 + v1:v7 + v1:v8 + v5:v2 +
      v5:v3 + v5:v4 + v5:v6 + v5:v7 + v5:v8 + v6:v2 + v6:v3 + v6:v4 + v6:v7 +
      v6:v8
  ) |> 
  # 各特徴量の2乗項
  step_poly(v1, v2, v3, v4, v5, v6, v7, v8, degree = 2) %>%
  # Elastic Net 用に標準化
  step_normalize(all_predictors())
#
# 新しいワークフロー
glmnet_wf2 <- workflow() |> 
  add_model(glmnet_mod) |> 
  add_recipe(glmnet_rec2)
#
# 新しいチューニング・グリッド
glmnet_grid2 <- grid_regular(
  penalty(range = c(-4, 0)),
  mixture(range = c(0.2, 0.8)),
  levels = c(40, 7)
)
#
# 新しいパラメータのチューニング
set.seed(123)
glmnet_tuned2 <- tune_grid(
  glmnet_wf2,
  grid = glmnet_grid2,
  resamples = folds,
  control = control_resamples(save_pred = TRUE)
)
#
# 新しいベストパラメータ
glmnet_params2 <- select_best(glmnet_tuned2, 
                              metric = "rmse")
glmnet_params2
#
# 新しい最終ワークフロー
glmnet_wf_f2 <- finalize_workflow(glmnet_wf2, glmnet_params2)
#
# 新しいfit
glmnet_fit2 <- fit(glmnet_wf_f2, train_data)
#
# 新しい予測
glmnet_pred2 <- test_data |> 
  select(target) |> 
  bind_cols(predict(glmnet_fit2, test_data)) |> 
  mutate(model = "glmnet2")
#
# 新しい予測も追加
pred_results <- pred_results |> 
  bind_rows(glmnet_pred2)
#
# 新しい結果発表
# RMSE
results_rmse <- pred_results |> 
  group_by(model) |> 
  rmse(truth = target, estimate = .pred) |> 
  arrange(.estimate) |> 
  mutate(score = scale(.estimate) |> as.numeric() * -1)
results_rmse
#
# MAE
results_mae <- pred_results |> 
  group_by(model) |> 
  mae(truth = target, estimate = .pred) |> 
  arrange(.estimate) |> 
  mutate(score = scale(.estimate) |> as.numeric() * -1)
results_mae
#
# R-Square
results_rsq <- pred_results |> 
  group_by(model) |> 
  rsq(truth = target, estimate = .pred) |> 
  arrange(desc(.estimate)) |> 
  mutate(score = scale(.estimate) |> as.numeric())
results_rsq
#
# 13-4. 総合評価
results_rmse |> 
  bind_rows(results_mae, results_rsq) |> 
  group_by(model) |> 
  summarize(total_score = mean(score)) |> 
  arrange(desc(total_score))
#
# 新しいペナルティ付き線形モデルの係数を確認する
# Step 1. エンジンを取り出す
glmnet_engine <- extract_fit_engine(glmnet_fit2)
#
# Step 2. coef()で係数を表示
coef(glmnet_engine, s = glmnet_params2$penalty)
#

(冒頭の画像は、Bing Image Creator で検索しました。プロンプトは Landscape of natural great river, close up of white Holly flower trees, Photo です。)

 

読書記録 - 「楽毅 (二)」 宮城谷 昌光 著 (新潮文庫)

楽毅は軍を率いても、自分一人での戦闘でも、行政でも超一級品の能力ですが、それでも使える君主、中山王が愚かなために中山国は滅んでいきます。この話は現代の日本の企業でも社長がダメなために、その下の重役、幹部、社員がいくら頑張っても会社は傾いてしまう、というようなことと同じだなと思いました。

中山国が滅んでいく、という大きな流れの中で、楽毅自身もその流れには抗うことができないということがわかっていて、それでも力の限り抗い、趙に戦いを挑む姿が、私の心に響きました。楽毅の一巻目の冒頭に「人がみごとに生きることは、むずかしいものだな。」と楽毅がつぶやいた言葉がこの二巻目でも大きなテーマのように思いました。

 

都道府県別の商業動態統計調査のデータの分析8 - 7つのモデルの予測結果比較

www.crosshyou.info

の続きです。前回は、tune_grid()関数でハイパーパラメータのチューニングをしました。

今回はチューニング結果を見るところからはじめます。collect_metrics()関数を使います。

線形モデル(lm)をみてみます。

rmseは2.30です。R-squaredは0.311です。

次は、ペナルティ付き線形モデル(glmnet)です。

rmseの最小は、2.22でした。

次は決定木モデル(rpart)です。

2.48が最小です。

次はランダムフォレストモデル(ranger)です。

2.08です。

次は、サポートベクター回帰(SVR)モデルです。(kernlab)

2.17が最小でした。

k-NN回帰モデル(kknn)はどうでしょうか?

2.18が最小でした。

最後はニューラルネットワークモデル(nnet)です。

次は、select_best()関数で最適なパラメータを取り出します。

この最適なパラメータで、最終ワークフローを作成します。finalize_workflow()関数を使います。

最終ワークフローが作成できたら、fit()関数で、train_dataをつかって学習します。

これでモデルが学習できました。あともう少しです。predict()関数でtest_dataの予測をします。

予測結果をみていきましょう。まず全ての予測結果を統合します。

RMSEから見ていきます。

決定木モデル(rpart)が一番小さいRMSEでした。このRMSEの値でZスコアを計算しています。このZスコアは後で総合評価を決定するときに使います。2番目はランダムフォレストモデル(ranger)で、3番目はペナルティ付き線形モデル(glmnet)でした。

MAEはどうでしょうか?

RMSEでは一番成績がよかった決定木モデル(rpart)がMAEでは一番悪い結果です。1番MAEが小さいのはランダムフォレストモデル(ranger)です。

R-squaredの結果を見てみます。

ランダムフォレストモデル(ranger)が一番の成績です。

それでは、RMSE、MAE、R-squaredの3つの総合をみてみます。scoreの平均値で決定します。

ランダムフォレストモデル(ranger)が一番良いモデルでした。2番目は線形モデル(lm)で、3番目はペナルティ付き線形モデル(glmnet)でした。線形モデル、あなどれませんね。

今回は以上です。

次回は

www.crosshyou.info

です。

 

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。


# 7. 性能確認
# 7-1. 線形モデル(lm)
collect_metrics(lm_tuned)
#
# 7-2. ペナルティ付き線形モデル(glmnet)
collect_metrics(glmnet_tuned) |> 
  filter(.metric == "rmse") |> 
  arrange(mean) |> 
  relocate(mean)
#
# 7-3. 決定木モデル(rpart)
collect_metrics(rpart_tuned) |> 
  filter(.metric == "rmse") |> 
  arrange(mean) |> 
  relocate(mean)
#
# 7-4. ランダムフォレストモデル(ranger)
collect_metrics(ranger_tuned) |> 
  filter(.metric == "rmse") |> 
  arrange(mean) |> 
  relocate(mean)
#
# 7-5. サポートベクター回帰(SVR)モデル(kernlab)
collect_metrics(kernlab_tuned) |> 
  filter(.metric == "rmse") |> 
  arrange(mean) |> 
  relocate(mean)
#
# 7-6. k-NN回帰モデル(kknn)
collect_metrics(kknn_tuned) |> 
  filter(.metric == "rmse") |> 
  arrange(mean) |> 
  relocate(mean)
#
# 7-7. ニューラルネットワークモデル(nnet)
collect_metrics(nnet_tuned) |> 
  filter(.metric == "rmse") |> 
  arrange(mean) |> 
  relocate(mean)
#
#
# 8ベストパラメータ
# 8-1. 線形モデル(lm)
# 線形モデルはパラメータは無し

# 8-2. ペナルティ付き線形モデル(glmnet)
glmnet_params <- glmnet_tuned |> 
  select_best(metric = "rmse")
glmnet_params
#
# 9-3. 決定木モデル(rpart)
rpart_params <- rpart_tuned |> 
  select_best(metric = "rmse")
rpart_params
#
# 9-4. ランダムフォレストモデル(ranger)
ranger_params <- ranger_tuned |> 
  select_best(metric = "rmse")
ranger_params
#
# 9-5. サポートベクター回帰(SVR)モデル(kernlab)
kernlab_params <- kernlab_tuned |> 
  select_best(metric = "rmse")
kernlab_params
#
# 9-6. k-NN回帰モデル(kknn)
kknn_params <- kknn_tuned |> 
  select_best(metric = "rmse")
kknn_params
#
# 9-7. ニューラルネットワークモデル(nnet)
nnet_params <- nnet_tuned |> 
  select_best(metric = "rmse")
nnet_params
#
# 10. 最終ワークフロー
# 10-1. 線形モデル(lm)
lm_wf_f <- lm_wf #線形モデルはワークフローと同じ
#
# 10-2. ペナルティ付き線形モデル(glmnet)
glmnet_wf_f <- finalize_workflow(glmnet_wf, glmnet_params)
#
# 10-3. 決定木モデル(rpart)
rpart_wf_f <- finalize_workflow(rpart_wf, rpart_params)
#
# 10-4. ランダムフォレストモデル
ranger_wf_f <- finalize_workflow(ranger_wf, ranger_params)
#
# 10-5. サポートベクター回帰(SVR)回帰モデル(kernlab)
kernlab_wf_f <- finalize_workflow(kernlab_wf, kernlab_params)
#
# 10-6. k-NN回帰モデル(kknn)
kknn_wf_f <- finalize_workflow(kknn_wf, kknn_params)
#
# 10-7. ニューラルネットワークモデル(nnet)
nnet_wf_f <- finalize_workflow(nnet_wf, nnet_params)
#
#
# 11. 最終トレーニング
# 11-1. 線形モデル(lm)
lm_fit <- fit(lm_wf_f, train_data)
#
# 11-2. ペナルティ付き線形モデル(glmnet)
glmnet_fit <- fit(glmnet_wf_f, train_data)
#
# 11-3. 決定木モデル(rpart)
rpart_fit <- fit(rpart_wf_f, train_data)
#
# 11-4. ランダムフォレストモデル(ranger)
ranger_fit <- fit(ranger_wf_f, train_data)
#
# 11-5. サポートベクター回帰(SVR)モデル(kernlab)
kernlab_fit <- fit(kernlab_wf_f, train_data)
#
# 11-6. k-NN回帰モデル(kknn)
kknn_fit <- fit(kknn_wf_f, train_data)
#
# 11-7. ニューラルネットワークモデル(nnet)
nnet_fit <- fit(nnet_wf_f, train_data)
#
#
# 12. テスト用データで予測
# 12-1. 線形モデル(lm)
lm_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(lm_fit, test_data)) |> 
  mutate(model = "lm")
#
# 12-2. ペナルティ付き線形モデル
glmnet_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(glmnet_fit, test_data)) |> 
  mutate(model = "glmnet")
#
# 12-3. 決定木モデル(rpart)
rpart_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(rpart_fit, test_data)) |> 
  mutate(model = "rpart")
#
# 12-4. ランダムフォレストモデル(ranger)
ranger_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(ranger_fit, test_data)) |> 
  mutate(model = "ranger")
#
# 12-5. サポートベクター回帰(SVR)モデル(kernlab)
kernlab_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(kernlab_fit, test_data)) |> 
  mutate(model = "kernlab")
#
# 12-6. k-NN回帰モデル(kknn)
kknn_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(kknn_fit, test_data)) |> 
  mutate(model = "kknn")
#
# 12-7. ニューラルネットワークモデル(nnet)
nnet_pred <- test_data |> 
  select(target) |> 
  bind_cols(predict(nnet_fit, test_data)) |> 
  mutate(model = "nnet")
#
#
# 12-9. 全ての予想を統合
pred_results <- lm_pred |> 
  bind_rows(glmnet_pred, rpart_pred, ranger_pred, kernlab_pred, kknn_pred,
            nnet_pred)
#
# 13. 結果発表
# 13-1. RMSE
results_rmse <- pred_results |> 
  group_by(model) |> 
  rmse(truth = target, estimate = .pred) |> 
  arrange(.estimate) |> 
  mutate(score = scale(.estimate) |> as.numeric() * -1)
results_rmse
#
# 13-2. MAE
results_mae <- pred_results |> 
  group_by(model) |> 
  mae(truth = target, estimate = .pred) |> 
  arrange(.estimate) |> 
  mutate(score = scale(.estimate) |> as.numeric() * -1)
results_mae
#
# 13-3 R-Square
results_rsq <- pred_results |> 
  group_by(model) |> 
  rsq(truth = target, estimate = .pred) |> 
  arrange(desc(.estimate)) |> 
  mutate(score = scale(.estimate) |> as.numeric())
results_rsq
#
# 13-4. 総合評価
results_rmse |> 
  bind_rows(results_mae, results_rsq) |> 
  group_by(model) |> 
  summarize(total_score = mean(score)) |> 
  arrange(desc(total_score))
#

 

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Landscape of natural spring season higher land, close up of Hyacinthus orientalis, photo です。)

 

都道府県別の商業動態統計調査のデータの分析7 - tidymodelによるチューニング - 勾配ブースティング回帰モデルは向いてなかった

www.crosshyou.info

の続きです。前回はチューニング・グリッドを作成しましたから、今回は tune_grid()関数で実際にチューニングをしていきます。

線形モデルは、fit_resample()関数を使います。

futureパッケージとtuneパッケージを読み込み、plan()関数で並列処理をしています。チューニングの計算時間を少しでも短くするためです。

ペナルティ付き線形モデル(glmnet)以降は、tune_grid()関数を使います。

警告メッセージが表示されていますが、気にしなくていいです。Fold06のデータで、推定値が定数になってしまった、という意味ですが、他のFoldでは推定できているので大丈夫です。

勾配ブースティング回帰モデル(xgboost)でたくさん警告メッセージが出てしまいました。Copilotに相談したら、レシピとチューニング・グリッドがダメだったようです。

レシピは

xgboost_rec <- recipe(target ~ ., data = train_data) |>
  step_zv(all_predictors()) |>        # ゼロ分散の列を削除
  step_dummy(all_nominal_predictors()) |>  # カテゴリ変数をダミー化
  step_normalize(all_numeric_predictors())  # 数値を標準化

チューニンググリッドは、

xgboost_grid <- grid_space_filling(
  mtry(range = c(1L, 8L)),
  min_n(range = c(2L, 20L)),
  tree_depth(range = c(3L, 8L)),
  learn_rate(range = c(-1.5, -0.5)),   # 0.03〜0.3
  loss_reduction(range = c(0, 5)),     # gamma >= 0 が安定
  sample_size(range = c(0, 1)),        # 0〜1 の実数を自動生成
  size = 20
)

このように変えたらいいらしいです。

変えてやってみます。

変えてもだめでした。。。

さらにCopilotに相談したら、train_dataのglimpse()を見せて、ということだったので、見せたら、

ということでした。

そして、ここでは割愛しますが、さらに調整したコードを走らせたのですが、結果は改善しませんでした。

そして、最終的なCopilotの結論は、

でした。

ということなので、勾配ブースティング回帰モデルはあきらめます。

Copilotによる詳しい説明は、以下のようでした。

🔍 なぜこのデータで xgboost が木を作れないのか
① 特徴量が 8 列しかない(低次元)
xgboost は本来:
- 数百〜数千の特徴量
- 非線形な関係
- カテゴリ変数の大量の dummy
- ノイズの多いデータ
で威力を発揮します。
しかし今回のデータは:
- 数値 8 列
- 112 行
- 分布も比較的素直
→ 木が分割する余地がほとんどない

② v3, v4, v7, v8 に “0 が大量に含まれる”
例:
- v4: 0 が多い
- v8: 0 が多い
- v3, v7: 負の値が多く偏りが強い
fold によっては:
- ほぼ 0 しかない
- ほぼ負の値しかない
- 分散が極端に小さい
→ xgboost は「分割しても意味がない」と判断
→ 定数予測

③ 112 行 × 5-fold CV → fold が小さすぎる
1 fold あたり 22 行。
この規模だと:
- 目的変数の分布が fold ごとに大きく変わる
- 特徴量の分散も fold ごとに変わる
- → 分割が成立しない fold が大量に発生
- → 全部同じ値を返す

④ このデータは “ほぼ線形” の可能性が高い
あなたのデータを見る限り:
- v1〜v8 が target と線形関係に見える
- 変数間の相関も高そう
- 非線形性が弱い
xgboost は非線形を掘りに行くモデルなので、
線形に近いデータでは 木が成長しない → 定数予測 が起きます。

🔥 つまり、xgboost が悪いのではなく、
このデータが xgboost に向いていない

今回のデータは、勾配ブースティング回帰モデルのような高性能モデルは必要ない、ということですね。

今回は以上です。

次回は、

www.crosshyou.info

です。

 

はじめから読むには、

www.crosshyou.info

です。

今回のコードは以下になります。

#
# 6. 各モデルのチューニング
# 6-1. 線形モデル(lm)
library(future)
library(tune)
plan(multisession, workers = parallel::detectCores() - 1) # 並列処理スタート

set.seed(123)
lm_tuned <- fit_resamples(
  lm_wf,
  resamples = folds,
  control = control_resamples(save_pred = TRUE)
)
#
# 6-2. ペナルティ付き線形モデル(glmnet)
set.seed(123)
glmnet_tuned <- tune_grid(
  glmnet_wf,
  grid = glmnet_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)
# 6-3. 決定木モデル(rpart)
set.seed(123)
rpart_tuned <- tune_grid(
  rpart_wf,
  grid = rpart_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)
#
# 6-4. ランダムフォレストモデル(ranger)
set.seed(123)
ranger_tuned <- tune_grid(
  ranger_wf,
  grid = ranger_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)
# 6-5. サポートベクター回帰(SVR)モデル(kernlab)
library(kernlab)
set.seed(123)
kernlab_tuned <- tune_grid(
  kernlab_wf,
  grid = kernlab_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)
#
# 6-6. k-NN回帰モデル(kknn)
set.seed(123)
kknn_tuned <- tune_grid(
  kknn_wf,
  grid = kknn_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)
#
# 6-7. ニューラルネットワークモデル(nnet)
set.seed(123)
nnet_tuned <- tune_grid(
  nnet_wf,
  grid = nnet_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)
#
# 6-8. 勾配ブースティング回帰モデル(xgboost)
set.seed(123)
xgboost_tuned <- tune_grid(
  xgboost_wf,
  grid = xgboost_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)

plan(sequential) # 並列処理エンド
#
# 勾配ブースティング回帰モデルのレシピを再定義
xgboost_rec <- recipe(target ~ ., data = train_data) |>
  step_zv(all_predictors()) |>        # ゼロ分散の列を削除
  step_dummy(all_nominal_predictors()) |>  # カテゴリ変数をダミー化
  step_normalize(all_numeric_predictors())  # 数値を標準化
#
# 勾配ブースティング回帰モデルのワークフローを再定義
xgboost_wf <- workflow() |> 
  add_model(xgboost_mod) |> 
  add_recipe(xgboost_rec)
#
# 勾配ブースティング回帰モデルのチューニンググリッドを再定義
xgboost_grid <- grid_space_filling(
  mtry(range = c(1L, 8L)),
  min_n(range = c(2L, 20L)),
  tree_depth(range = c(3L, 8L)),
  learn_rate(range = c(-1.5, -0.5)),   # 0.03〜0.3
  loss_reduction(range = c(0, 5)),     # gamma >= 0 が安定
  sample_size(range = c(0, 1)),        # 0〜1 の実数を自動生成
  size = 20
)
#
# 勾配ブースティング回帰モデルの再チューニング
plan(multisession, workers = parallel::detectCores() - 1) # 並列処理スタート

set.seed(123)
xgboost_tuned <- tune_grid(
  xgboost_wf,
  grid = xgboost_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)

plan(sequential) # 並列処理エンド
#
# さらに再定義
xgboost_rec <- recipe(target ~ ., data = train_data) %>%
  step_normalize(all_numeric_predictors())

xgboost_wf <- workflow() |> 
  add_model(xgboost_mod) |> 
  add_recipe(xgboost_rec)

folds <- vfold_cv(train_data, v = 5, repeats = 3)

xgboost_grid <- grid_space_filling(
  tree_depth(range = c(3L, 6L)),
  learn_rate(range = c(-1.3, -0.7)),  # 0.05〜0.2
  min_n(range = c(2L, 10L)),
  loss_reduction(range = c(0, 2)),
  sample_size(range = c(0, 1)),
  mtry(range = c(1L, 8L)),
  size = 20
)

plan(multisession, workers = parallel::detectCores() - 1) # 並列処理スタート

set.seed(123)
xgboost_tuned <- tune_grid(
  xgboost_wf,
  grid = xgboost_grid,
  resamples = folds,
  control = control_grid(save_pred = TRUE)
)

plan(sequential) # 並列処理エンド
#

 

 

(冒頭の画像は、Bing Image Creator で生成しました。プロンプトは、Natural scene landscape, close up rainbow color flowers, green grass fields, one tall tree, blue sky, photo です。)