都道府県別の商業動態統計調査のデータの分析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)
#