Tidymodelsでデータの前処理内容を確認する

カテゴリ: r

tidymodelsはRにおける統計モデリングや機械学習を便利にするためのフレームワークです。 tidymodelsを利用すると

  • パイプ演算子による処理の流れが明瞭なモデリング
  • パッケージごとに異なる学習・予測インターフェースの統一

といったメリットを享受でき、徐々にはやってきている印象です。

tidymodelsはパッケージでもあります。 tidyverseのように複数のパッケージの集合体で、特徴量エンジニアリング担当のrecipesパッケージ、学習・予測担当のparsnipパッケージなどが属します。

もっと詳しい概要や、実際の分析の流れはUryuさんのスライドを参考にしてください。

tidymodelsで覚えるRでのモデル構築と運用

ここではrecipesパッケージで前処理を行った後、前処理の内容を調査する方法を紹介します。たとえば、

  • 中心化した時の平均値ってどんなもんだったんだろ?
  • PCAかけたけど、固有ベクトルってどんな感じだったのかなあ?

といったことを確認する方法を紹介します。

例として、ダイヤモンドの価格のデータセット(ggplot2::diamonds)を使ってみましょう。簡単のため、数値の列だけを扱います。

前処理までのざっとした流れは以下のとおり。実際に則すようにデータ分割から始めていますが、本質的なところは、recipes::recipe関数から始まる前処理方法の定義と、続くrecipes::prep関数による前処理定義の学習です。

重要: 後で前処理の各ステップの検証をしやすくするため、step_*関数のid引数に、stepごとに固有の名前を指定してください。検証不要なステップのidは省略可能です。

# パイプ演算子の読み込み
`%>%` <- magrittr::`%>%`

# シード固定
set.seed(71)

# データ分割
split <- ggplot2::diamonds %>%
  dplyr::select(where(is.numeric)) %>%
  rsample::initial_split(p = .9)

# 訓練データの取り出し
tr <- rsample::training(split)

# 前処理方法を定義
rec <- recipes::recipe(tr, price ~ .) %>%
  recipes::step_center(-price, id = 'center') %>%
  recipes::step_scale(-price, id = 'scale') %>%
  recipes::step_pca(-price, id = 'pca') %>%
  recipes::step_log(price)

# 前処理方法を学習
prep <- recipes::prep(rec)

前処理方法の学習が完了したら、prepオブジェクトの中身を除いてみましょう。大きなリストなので、ここでは名前だけを表示してみます。

names(prep)
#> [1] "var_info"       "term_info"      "steps"          "template"      
#> [5] "retained"       "tr_info"        "orig_lvls"      "last_term_info"

"steps"といういかにも怪しげな要素がありますね。ではprep$stepsの1番目の要素の構造をstr関数で覗いてみましょう。

str(prep$steps[1L])
#> List of 1
#>  $ :List of 7
#>   ..$ terms  :List of 1
#>   .. ..$ : language ~-price
#>   .. .. ..- attr(*, ".Environment")=<environment: 0x562acefa2100> 
#>   .. ..- attr(*, "class")= chr [1:2] "quosures" "list"
#>   ..$ role   : logi NA
#>   ..$ trained: logi TRUE
#>   ..$ means  : Named num [1:6] 0.799 61.752 57.456 5.733 5.736 ...
#>   .. ..- attr(*, "names")= chr [1:6] "carat" "depth" "table" "x" ...
#>   ..$ na_rm  : logi TRUE
#>   ..$ skip   : logi FALSE
#>   ..$ id     : chr "center"
#>   ..- attr(*, "class")= chr [1:2] "step_center" "step"

どうやらprep$stepsはリストで、1番目の内容を見ると、idという要素があり、"center"という文字列が入っています。前処理を定義した時、最初のステップにrecipes::step_centerを引数id = "center"と共に使ったことと対応していますね。また、meansという要素もあり、ここには中心化に使った平均値が入っているようです。

prep$steps[[1L]]$means
#>      carat      depth      table          x          y          z 
#>  0.7985464 61.7516324 57.4559705  5.7325532  5.7360282  3.5391435

念のため、トレーニングデータであるtrの各列の平均値を見てみましょう。

colMeans(tr)
#>        carat        depth        table        price            x            y 
#>    0.7985464   61.7516324   57.4559705 3937.0699734    5.7325532    5.7360282 
#>            z 
#>    3.5391435

priceの平均値の有無こそ違いますが、他は一致しています。一安心ですね。 priceの有無は中心化では予測対象なので対象外としましたが、colMeans関数は全変数を対象としていることに由来します。

さて、もう少し汎用的に検証できるよう、inspect_step関数を定義してみます。これは第一引数に学習済みの定義、第二引数に任意の前処理ステップのidをとります。

inspect_step <- function(recipe, step_id) {
  recipe$steps[
    purrr::map_chr(recipe$steps, "id") == step_id
  ][[1L]]
}

実際にPCAの結果を取り出してstr関数で構造を見てみます。 idに"pca"が入っていることを確認できます。

str(inspect_step(prep, "pca"))
#> List of 10
#>  $ terms    :List of 1
#>   ..$ : language ~-price
#>   .. ..- attr(*, ".Environment")=<environment: 0x562acefa2250> 
#>   ..- attr(*, "class")= chr [1:2] "quosures" "list"
#>  $ role     : chr "predictor"
#>  $ trained  : logi TRUE
#>  $ num_comp : num 5
#>  $ threshold: logi NA
#>  $ options  : list()
#>  $ res      :List of 4
#>   ..$ sdev    : num [1:6] 1.988 1.133 0.827 0.222 0.162 ...
#>   ..$ rotation: num [1:6, 1:6] 0.494469 -0.000991 0.120718 0.500085 0.493796 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:6] "carat" "depth" "table" "x" ...
#>   .. .. ..$ : chr [1:6] "PC1" "PC2" "PC3" "PC4" ...
#>   ..$ center  : logi FALSE
#>   ..$ scale   : logi FALSE
#>   ..- attr(*, "class")= chr "prcomp"
#>  $ prefix   : chr "PC"
#>  $ skip     : logi FALSE
#>  $ id       : chr "pca"
#>  - attr(*, "class")= chr [1:2] "step_pca" "step"

せっかくなので各主成分の固有ベクトルを取り出してみましょう。

inspect_step(prep, "pca")$res$rotation
#>                 PC1          PC2          PC3         PC4          PC5
#> carat  0.4944690257 -0.045142961  0.027733335 -0.69872947  0.503950252
#> depth -0.0009910619 -0.734418543 -0.670622228  0.04053755  0.045526726
#> table  0.1207179046  0.669542455 -0.732751827  0.01440279  0.001342399
#> x      0.5000854569 -0.008609363  0.070109159 -0.06174404 -0.429135337
#> y      0.4937960687 -0.009923418  0.087127512  0.70974748  0.494690171
#> z      0.4969975516 -0.100656952 -0.008059259  0.04870922 -0.561323575
#>                PC6
#> carat  0.102492101
#> depth -0.084724689
#> table  0.002367648
#> x     -0.746294618
#> y     -0.003989857
#> z      0.652180892

現状、recipesパッケージではこの手の痒いところに手が届く関数が整備されていないようです。とはいえ、オブジェクトの内部構造はよく練られているので、困ったらstr関数を使って検証してみると良いでしょう。これはモデル解釈の時にも同じことが言えます。

tidymodelsでもxgboostを解釈したい

Enjoy!