rpartの結果をtidyに扱う（変数重要度と分岐情報を取り出す）
『Rユーザのためのtidymodels[実践]入門』を読み進めている。
Rユーザのためのtidymodels[実践]入門 〜モダンな統計・機械学習モデリングの世界：書籍案内｜技術評論社
そのtidymodels関係で、決定木モデルを作成するrpartをtidyにするパッケージのお話。
kyphosisデータ
サンプルデータはrpartパッケージのkyphosisデータ。
# https://rdrr.io/cran/rpart/man/kyphosis.html data(kyphosis, package = "rpart") kyphosis %>% summary() # Kyphosis Age Number Start # absent :64 Min. : 1.00 Min. : 2.000 Min. : 1.00 # present:17 1st Qu.: 26.00 1st Qu.: 3.000 1st Qu.: 9.00 # Median : 87.00 Median : 4.000 Median :13.00 # Mean : 83.65 Mean : 4.049 Mean :11.49 # 3rd Qu.:130.00 3rd Qu.: 5.000 3rd Qu.:16.00 # Max. :206.00 Max. :10.000 Max. :18.00
broomパッケージ
まずはtidymodelsのbroomパッケージの関数について。
主に3つあり、前処理やモデルの実行結果をtidyにするtidy()関数、学習データに予測値と残差を付与するaugment()関数、モデル選択に利用するAICや決定係数を取得するglance()関数。
この記事ではtidy()関数を扱う。
挙動確認のため、試しにtidy()関数をlm()関数の結果に適用。偏回帰係数などがtidyに返ってくる。
lm_fit <- lm(Start ~ Age + Number, data = kyphosis) broom::tidy(lm_fit) # A tibble: 3 × 5 # term estimate std.error statistic p.value # <chr> <dbl> <dbl> <dbl> <dbl> # (Intercept) 16.3 1.54 10.6 7.88e-17 # Age 0.00427 0.00860 0.496 6.21e- 1 # Number -1.28 0.309 -4.15 8.55e- 5
一方、rpartの結果を渡すとエラーになってしまう。rpartクラスオブジェクトへのメソッドは無い、と。
fit_rpart <- rpart::rpart(Kyphosis ~ Age + Number + Start, data = kyphosis) broom::tidy(fit_rpart) # Error: No tidy method for objects of class rpart
broomstickパッケージ
そこで出てくるのがbroomstickパッケージ。broomstickとは箒の柄のこと。
broomstick.njtierney.com
2023年1月21日時点でCRANに登録されていないので、GitHubからいただいてくる。
このパッケージのtidy()関数を適用すると、変数重要度が返ってくる。
remotes::install_github("njtierney/broomstick") broomstick::tidy(fit_rpart) # A tibble: 3 × 2 # variable importance # <chr> <dbl> # Start 8.20 # Age 3.10 # Number 1.52
tidyrulesパッケージ
また、別のrpartに使えるパッケージでtidyrulesというものもある。 talegari.github.io
このパッケージのtidyRules()関数を適用すると、分岐情報などの実行結果が返ってくる。
tidyrules::tidyRules(fit_rpart) # A tibble: 5 × 6 # id LHS RHS support confidence lift # <int> <chr> <chr> <int> <dbl> <dbl> # 1 Start >= 8.5 & Start >= 14.5 absent 29 0.968 1.22 # 2 Start >= 8.5 & Start < 14.5 & Age < 55 absent 12 0.929 1.18 # 3 Start >= 8.5 & Start < 14.5 & Age >= 55 & Age >= 111 absent 14 0.812 1.03 # 4 Start >= 8.5 & Start < 14.5 & Age >= 55 & Age < 111 present 7 0.556 2.65 # 5 Start < 8.5 present 19 0.571 2.72
各変数の内容は以下の通り。
LHS: Rules.RHS: Predicted Class.support: Number of observation covered by the rule.confidence: Prediction accuracy for respective class. (laplace correction is implemented by default)lift: The result of dividing the rule's estimated accuracy by the relative frequency of the predicted class in the training set.
tidyrules/tidyrules_vignette.Rmd at master · talegari/tidyrules · GitHub
通常分岐情報を取得するには、以下のように内部変数にアクセスする必要があるところを楽にしてくれている。
rpart:::labels.rpart(fit_rpart, minlength = 0) # [1] "root" "Start>=8.5" "Start>=14.5" "Start< 14.5" "Age< 55" # [6] "Age>=55" "Age>=111" "Age< 111" "Start< 8.5" partykit::as.party(fit_rpart) %>% partykit:::.list.rules.party() %>% stringr::str_replace_all(pattern = "\\\"","'") %>% stringr::str_remove_all(pattern = ", 'NA'") %>% stringr::str_remove_all(pattern = "'NA',") %>% stringr::str_remove_all(pattern = "'NA'") %>% stringr::str_squish() # [1] "Start >= 8.5 & Start >= 14.5" # [2] "Start >= 8.5 & Start < 14.5 & Age < 55" # [3] "Start >= 8.5 & Start < 14.5 & Age >= 55 & Age >= 111" # [4] "Start >= 8.5 & Start < 14.5 & Age >= 55 & Age < 111" # [5] "Start < 8.5"
結びに
分岐情報取得の話は以前書いた。
rpartの決定木から分岐情報を取り出す - データ分析関連メモ（メモです）
スライドはこちら。
speakerdeck.com
以上
rsample::sliding_window()関数の引数
『Rユーザのためのtidymodels[実践]入門』を読み始めた。
Rユーザのためのtidymodels[実践]入門 〜モダンな統計・機械学習モデリングの世界：書籍案内｜技術評論社
第1章で時系列データを分析データと検証データに分割する関数rsample::sliding_window()関数の説明がある。
この関数の引数について自分向けに整理。
まずはデータの読み込み。データは参考書と同様。分割の仕方がわかり易いようにindexを追加している。
"S4248SM144NCEN"という販売量が入っている列名は、FRED（Federal Reserve Economic Data）のコードを意味しているのかもしれない。
Merchant Wholesalers, Except Manufacturers' Sales Branches and Offices: Nondurable Goods: Beer, Wine, and Distilled Alcoholic Beverages Sales (S4248SM144NCEN) | FRED | St. Louis Fed
data(drinks, package = "modeldata") drinks_annual <- drinks %>% dplyr::mutate(year = lubridate::year(date)) %>% dplyr::filter(dplyr::between(year, 1992, 1994)) %>% dplyr::mutate(index = dplyr::row_number(x = date)) drinks_annual %>% head() # date sales year index # <date> <dbl> <dbl> <int> # 1992-01-01 3459 1992 1 # 1992-02-01 3458 1992 2 # 1992-03-01 4002 1992 3 # 1992-04-01 4564 1992 4 # 1992-05-01 4221 1992 5 # 1992-06-01 4529 1992 6
次にsliding_window()の基本的な使い方について。
先ず、rsample::initial_time_split()関数をデータに適用して、rsplitクラスのオブジェクト（drinks_split_ts2）を取得する。
次に、rsplitクラスのオブジェクトにrsample::training()関数を適用して学習データ（train_ts_data2）を参照する。
drinks_split_ts2 <- drinks_annual %>% rsample::initial_time_split(prop = 0.68) train_ts_data2 <- drinks_split_ts2 %>% rsample::training()
この学習データ（train_ts_data2）にrsample::sliding_window()関数を適用すれば分析データと検証データに分割できるのだが、
この分割方法をrsample::sliding_window()関数の引数で指定することになる。
rsample::sliding_window()関数の引数は以下の通り（書籍から引用）。
lookbackが分析セット、assess系の2つは検証セットの操作をしている。
lookback……分析セットの件数（時系列の起点となる地震を含まない件数）
assess_start……検証セットの時系列上での起点
assess_stop……検証セットの時系列上での終点
引数のデフォルト値はそれぞれ0,1,1となっている。
Time-based Resampling — slide-resampling • rsample
sliding_window(
data,
...,
lookback = 0L,
assess_start = 1L,
assess_stop = 1L,
complete = TRUE,
step = 1L,
skip = 0L
)
ここから、引数に様々な値を与えて出力結果を確認していく。
デフォルトのまま実行すると、分析データは最初の1件から、検証データは分析データの次の1件で、1つずつずれていく。
fold_default <- rsample::sliding_window(data = train_ts_data2) rsample::analysis(fold_default$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-01-01 3459 1992 1 rsample::assessment(fold_default$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-02-01 3458 1992 2 rsample::analysis(fold_default$splits[[2]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-02-01 3458 1992 2 rsample::assessment(fold_default$splits[[2]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-03-01 4002 1992 3
lookback = 10とすると、分析データ11件（起点となる1件が足され、10 + 1で11件）、検証データはデフォルトのまま1件。
fold_look10 <- rsample::sliding_window(data = train_ts_data2, lookback = 10) rsample::analysis(fold_look10$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-01-01 3459 1992 1 # 1992-02-01 3458 1992 2 # 1992-03-01 4002 1992 3 # 1992-04-01 4564 1992 4 # 1992-05-01 4221 1992 5 # 1992-06-01 4529 1992 6 # 1992-07-01 4466 1992 7 # 1992-08-01 4137 1992 8 # 1992-09-01 4126 1992 9 # 1992-10-01 4259 1992 10 # 1992-11-01 4240 1992 11 rsample::assessment(fold_look10$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-12-01 4936 1992 12 rsample::analysis(fold_look10$splits[[2]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-02-01 3458 1992 2 # 1992-03-01 4002 1992 3 # 1992-04-01 4564 1992 4 # 1992-05-01 4221 1992 5 # 1992-06-01 4529 1992 6 # 1992-07-01 4466 1992 7 # 1992-08-01 4137 1992 8 # 1992-09-01 4126 1992 9 # 1992-10-01 4259 1992 10 # 1992-11-01 4240 1992 11 # 1992-12-01 4936 1992 12 rsample::assessment(fold_look10$splits[[2]]) # date sales year index # <date> <dbl> <dbl> <int> # 1993-01-01 3031 1993 13
assess_start = 5とすると、分析データはデフォルトと同様に1件、検証データは分析データの5件あとのindexが6であるレコード1件。
assess_stopはassess_start以上の数値指定しないとエラーになる。
fold_start5 <- rsample::sliding_window(data = train_ts_data2, assess_start = 5, assess_stop = 5) rsample::analysis(fold_start5$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-01-01 3459 1992 1 rsample::assessment(fold_start5$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-06-01 4529 1992 6
assess_start = 5, assess_stop = 10とすると、検証データは分析データから数えて5件目~10件目の6件を取る。
fold_start5_stop10 <- rsample::sliding_window(data = train_ts_data2, assess_start = 5, assess_stop = 10) rsample::analysis(fold_start5_stop10$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-01-01 3459 1992 1 rsample::assessment(fold_start5_stop10$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-06-01 4529 1992 6 # 1992-07-01 4466 1992 7 # 1992-08-01 4137 1992 8 # 1992-09-01 4126 1992 9 # 1992-10-01 4259 1992 10 # 1992-11-01 4240 1992 11
lookback = 10, assess_start = 5, assess_stop = 10とすると、分析データは11件、検証データはそこから5件目~10件目を取る。
fold_look10_start5 <- rsample::sliding_window(data = train_ts_data2, lookback = 10, assess_start = 5, assess_stop = 5) rsample::analysis(fold_look10_start5$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1992-01-01 3459 1992 1 # 1992-02-01 3458 1992 2 # 1992-03-01 4002 1992 3 # 1992-04-01 4564 1992 4 # 1992-05-01 4221 1992 5 # 1992-06-01 4529 1992 6 # 1992-07-01 4466 1992 7 # 1992-08-01 4137 1992 8 # 1992-09-01 4126 1992 9 # 1992-10-01 4259 1992 10 # 1992-11-01 4240 1992 11 rsample::assessment(fold_look10_start5$splits[[1]]) # date sales year index # <date> <dbl> <dbl> <int> # 1993-04-01 4377 1993 16
日付の途中を補完するtidyr::complete()関数
時系列データでよく使うので備忘録。
まずtidyr::complete()関数の使い方から。
組み合わせの欠損を埋めてくれる関数。
Complete a data frame with missing combinations of data — complete • tidyr
メモ：時系列とか連番のデータを補完するときはtidyrのcomplete()とfull_seq()が便利そう - Technically, technophobic.
サンプルデータを作成する。idが”111”と”222”のデータがあるが、それぞれyear_monthの途中が欠けている。
id“111”は3,5月、”222”は2,3月のレコードが無い。
df_1 <- tibble::tribble( ~id, ~year_month, ~cnt, "111", "2022-01-01", 1, "111", "2022-02-01", 2, "111", "2022-04-01", 4, "222", "2022-01-01", 1, "222", "2022-04-01", 4, "222", "2022-05-01", 5, ) |> dplyr::mutate(year_month = as.Date(year_month))
seq.Date()関数とtidyr::complete()関数を組み合わせると、補完することができる。
df_1 |> tidyr::complete(id, year_month = seq.Date(as.Date(min(year_month, na.rm = TRUE)), as.Date(max(year_month, na.rm = TRUE)), by = "months")) # A tibble: 10 × 3 # id year_month cnt # <chr> <date> <dbl> # 1 111 2022-01-01 1 # 2 111 2022-02-01 2 # 3 111 2022-03-01 NA # 4 111 2022-04-01 4 # 5 111 2022-05-01 NA # 6 222 2022-01-01 1 # 7 222 2022-02-01 NA # 8 222 2022-03-01 NA # 9 222 2022-04-01 4 # 10 222 2022-05-01 5
ここから気になって調べたところ。
seq.Date()関数の始まりと終わりの日付に、元のデータに存在しないものを指定しても、指定した範囲で補完してくれる。
例えば元のサンプルデータには2022年1月から2022年5月までしか存在しないが、seq.Date()関数に2021年12月から2022年6月まで指定すると次のように返ってくる。
df_1 |> tidyr::complete(id, year_month = seq.Date(as.Date("2021-12-01"), as.Date("2022-06-01"), by = "months")) # A tibble: 14 × 3 # id year_month cnt # <chr> <date> <dbl> # 1 111 2021-12-01 NA # 2 111 2022-01-01 1 # 3 111 2022-02-01 2 # 4 111 2022-03-01 NA # 5 111 2022-04-01 4 # 6 111 2022-05-01 NA # 7 111 2022-06-01 NA # 8 222 2021-12-01 NA # 9 222 2022-01-01 1 # 10 222 2022-02-01 NA # 11 222 2022-03-01 NA # 12 222 2022-04-01 4 # 13 222 2022-05-01 5 # 14 222 2022-06-01 NA
完全に範囲外の1900年の日付を指定すると、その範囲だけ組み合わせが発生する。
元々存在していた2022年のレコードは変わらず。
df_1 |> tidyr::complete(id, year_month = seq.Date(as.Date("1900-01-01"), as.Date("1900-02-01"), by = "months")) # A tibble: 10 × 3 # id year_month cnt # <chr> <date> <dbl> # 1 111 1900-01-01 NA # 2 111 1900-02-01 NA # 3 222 1900-01-01 NA # 4 222 1900-02-01 NA # 5 111 2022-01-01 1 # 6 111 2022-02-01 2 # 7 111 2022-04-01 4 # 8 222 2022-01-01 1 # 9 222 2022-04-01 4 # 10 222 2022-05-01 5
次に、idとは別に、idに対して紐づいている変数を持っている場合。
サンプルデータとして、各idにnameがあるデータを作成する。
df_2 <- tibble::tribble( ~id, ~name, ~year_month, ~cnt, "111", "xxx", "2022-01-01", 1, "111", "xxx", "2022-02-01", 2, "111", "xxx", "2022-04-01", 4, "222", "yyy", "2022-01-01", 1, "222", "yyy" ,"2022-04-01", 4, "222", "yyy", "2022-05-01", 5, ) |> dplyr::mutate(year_month = as.Date(year_month))
これをそのままtidyr::complete()に放り込むと、idとnameを含めて、すべての組み合わせが作成される。
df_2 |> tidyr::complete(id, name, year_month = seq.Date(as.Date(min(year_month, na.rm = TRUE)), as.Date(max(year_month, na.rm = TRUE)), by = "months")) # A tibble: 20 × 4 # id name year_month cnt # <chr> <chr> <date> <dbl> # 1 111 xxx 2022-01-01 1 # 2 111 xxx 2022-02-01 2 # 3 111 xxx 2022-03-01 NA # 4 111 xxx 2022-04-01 4 # 5 111 xxx 2022-05-01 NA # 6 111 yyy 2022-01-01 NA # 7 111 yyy 2022-02-01 NA # 8 111 yyy 2022-03-01 NA # 9 111 yyy 2022-04-01 NA # 10 111 yyy 2022-05-01 NA # 11 222 xxx 2022-01-01 NA # 12 222 xxx 2022-02-01 NA # 13 222 xxx 2022-03-01 NA # 14 222 xxx 2022-04-01 NA # 15 222 xxx 2022-05-01 NA # 16 222 yyy 2022-01-01 1 # 17 222 yyy 2022-02-01 NA # 18 222 yyy 2022-03-01 NA # 19 222 yyy 2022-04-01 4 # 20 222 yyy 2022-05-01 5
idとnameの組み合わせは動かしたくない場合は、tidyr::nesting()関数で囲って渡す。
df_2 |> tidyr::complete(tidyr::nesting(id, name), year_month = seq.Date(as.Date(min(year_month, na.rm = TRUE)), as.Date(max(year_month, na.rm = TRUE)), by = "months")) # A tibble: 10 × 4 # id name year_month cnt # <chr> <chr> <date> <dbl> # 1 111 xxx 2022-01-01 1 # 2 111 xxx 2022-02-01 2 # 3 111 xxx 2022-03-01 NA # 4 111 xxx 2022-04-01 4 # 5 111 xxx 2022-05-01 NA # 6 222 yyy 2022-01-01 1 # 7 222 yyy 2022-02-01 NA # 8 222 yyy 2022-03-01 NA # 9 222 yyy 2022-04-01 4 # 10 222 yyy 2022-05-01 5
rpartの決定木から分岐情報を取り出す
rpartで作成した決定木から分岐の情報を抽出する。
かわいいかわいいpalmerpenguinsをサンプルデータとして決定木を作成。
library(palmerpenguins) tree <- rpart::rpart( formula = species ~ ., data = penguins, method = "class") tree # n= 344 # # node), split, n, loss, yval, (yprob) # * denotes terminal node # # 1) root 344 192 Adelie (0.441860465 0.197674419 0.360465116) # 2) flipper_length_mm< 206.5 214 64 Adelie (0.700934579 0.294392523 0.004672897) # 4) bill_length_mm< 43.35 151 5 Adelie (0.966887417 0.033112583 0.000000000) * # 5) bill_length_mm>=43.35 63 5 Chinstrap (0.063492063 0.920634921 0.015873016) * # 3) flipper_length_mm>=206.5 130 7 Gentoo (0.015384615 0.038461538 0.946153846) # 6) island=Dream,Torgersen 7 2 Chinstrap (0.285714286 0.714285714 0.000000000) * # 7) island=Biscoe 123 0 Gentoo (0.000000000 0.000000000 1.000000000) *
rpart.plotはこんな感じ。
rpart.plot::rpart.plot(tree)
分岐情報は内部変数として持っているのでrpart:::labels.rpartでアクセスできる。
labels.rpart function - RDocumentation
「:::」はTriple Colon Operatorと呼ばれている内部変数にアクセスする演算子。
R: Double Colon and Triple Colon Operators
rpart:::labels.rpart(tree) # [1] "root" "flipper_length_mm< 206.5" # [3] "bill_length_mm< 43.35" "bill_length_mm>=43.35" # [5] "flipper_length_mm>=206.5" "island=bc" # [7] "island=a"
内部でabbreviateという関数が使用されており、stringsはa,b,cといった形に省略されてしまう。今回のサンプルでもislandの分岐基準が「"island=bc" 」「"island=a”」と省略されてしまっている。
abbreviate function - RDocumentation
minlengthを0とすると省略せずに返してくれる。
rpart:::labels.rpart(tree, minlength = 0) # [1] "root" "flipper_length_mm< 206.5" # [3] "bill_length_mm< 43.35" "bill_length_mm>=43.35" # [5] "flipper_length_mm>=206.5" "island=Dream,Torgersen" # [7] "island=Biscoe"
データフレームに色々まとめることもできる。
r - Extract variable labels from rpart decision tree - Stack Overflow
df_splits <- data.frame( splits = rpart:::labels.rpart(tree, minlength = 0), n = tree$frame$n, yval = tree$frame$yval ) df_splits # splits n yval var # 1 root 344 1 flipper_length_mm # 2 flipper_length_mm< 206.5 214 1 bill_length_mm # 3 bill_length_mm< 43.35 151 1 <leaf> # 4 bill_length_mm>=43.35 63 2 <leaf> # 5 flipper_length_mm>=206.5 130 3 island # 6 island=Dream,Torgersen 7 2 <leaf> # 7 island=Biscoe 123 3 <leaf>
ggplot2のgeom_boxplotについて備忘録
ggplot2の箱ひげ図を描く関数の備忘録。
geom_boxplot
penguinsデータで箱ひげ図を描く。
library(tidyverse) library(palmerpenguins) data(package = 'palmerpenguins') df <- penguins df %>% ggplot2::ggplot(ggplot2::aes(x = species, y = bill_length_mm)) + ggplot2::geom_boxplot()
A box and whiskers plot (in the style of Tukey) — geom_boxplot • ggplot2
平均の線を引く
箱の真ん中の線は平均値ではなく中央値。四分位で分けるのが箱ひげ図なので、箱の上が第3四分位、下が第1四分位、真ん中が第2四分位でつまり中央値になる。
それでも平均値を描画したい場合はstat_summaryを使用する。（stat_summaryあんまり調べてないので怪しい。）
df %>% ggplot2::ggplot(ggplot2::aes(x = species, y = bill_length_mm)) + ggplot2::geom_boxplot() + ggplot2::stat_summary(fun = mean, geom = "errorbar", aes(ymax = ..y.., ymin = ..y..), width = 0.75, size = 1, linetype = "dotted")
r - How do I plot the mean instead of the median with geom_boxplot? - Stack Overflow
ひげの先端を最大値最小値にする
ひげの先端が何を表しているのかは決まっている訳ではないらしい。
geom_boxplotはデフォルトでは箱の端からIQRの1.5倍内の最小最大がひげの先端になる。
ひげより先にある点は外れ値と解釈される。IQRとはInterquartile Rangeの略で、第3四分位数から第1四分位数を引いたもの。
四分位範囲 | 統計用語集 | 統計WEB
一方で、ひげの先端を最大値最小値として説明している場合もある。
4-2. 箱ひげ図の見方 | 統計学の時間 | 統計WEB
geom_boxplotのひげの先端を最大値最小値にするには、ひげの長さを調整するcoefという引数に値を設定する。coefに設定した数値がIQRに掛けられた結果がひげの長さになる。デフォルトは1.5。
適当に大きな数値を設定すれば最大値最小値が先端になるが、なんか違う気がする……。
df %>% ggplot2::ggplot(ggplot2::aes(x = species, y = bill_length_mm)) + ggplot2::geom_boxplot(coef = 1000)
箱の並び替え
並び替えはreorderを使う。箱ひげ図の真ん中の線である中央値で並び替える場合が多いはず。
df %>% dplyr::filter(!is.na(bill_length_mm)) %>% ggplot2::ggplot(ggplot2::aes(x = reorder(x = species, X = bill_length_mm, FUN = median), y = bill_length_mm)) + ggplot2::geom_boxplot()
ここではdplyr::filterでNAを除外しているが、除外しないとこうなる。
df %>% ggplot2::ggplot(ggplot2::aes(x = reorder(x = species, X = bill_length_mm, FUN = median), y = bill_length_mm)) + ggplot2::geom_boxplot()
並び替えられていない。これはreorderのmedianがNAを除外していないため、集計結果がNAになってしまって並び替えが機能していない。
先のようにdplyr::filterで除外するか、もしくはreorderが...（dot-dot-dot）でFUNに引数を渡せるのでna.rm=TRUEを設定する。
… optional: extra arguments supplied to FUN
reorder.default function - RDocumentation
df %>% ggplot2::ggplot(ggplot2::aes(x = reorder(x = species, X = bill_length_mm, FUN = median, na.rm = TRUE), y = bill_length_mm)) + ggplot2::geom_boxplot()
predict関数から学習データしか返ってこない（dot-dot-dotについて）
R言語のrpartにあるpredict関数を使っていた時の話。
penguinsデータをtrain233件、valid100件に分割。
# packages ---------------------------------------------------------------- library(palmerpenguins) data(package = 'palmerpenguins') # data preparation --------------------------------------------------------------- set.seed(2022) df <- penguins |> tidyr::drop_na(tidyselect::everything()) n_size <- nrow(df) split <- sample(n_size, size = 0.7 * n_size) train <- df[split,] valid <- df[-split,] train |> nrow() # [1] 233 valid |> nrow() # [1] 100
rpartで学習してvalidに適用したつもりだったが、なぜか233件のデータが返ってきた。validは100件のはず。
# modeling ---------------------------------------------------------------- tr <- rpart::rpart(formula = species ~ ., data = train, method = "class", parms = list(split = "gini")) pred_data <- predict(object = tr, data = valid, type = "prob") pred_data |> nrow() # [1] 233 # validの100件が返ってくるはずなのに233件 valid |> nrow() # [1] 100 # validの件数を確認するがやはり100件
RDocumentationに戻って確認すると、モデルを適用するデータを指定する引数は”newdata”だったが誤って”data”という引数名にしていたことが分かった。
newdata
data frame containing the values at which predictions are required.
predict.rpart function - RDocumentation
pred_data <- predict(object = tr, data = valid, # ここをdataではなくnewdataにしないといけなかった type = "prob")
疑問が2つ出てきて、1つ目は”newdata”に何も渡していないのに問題ないのかということ、2つ目は”data”という私が勝手に創り上げた引数に渡したオブジェクトはどこに行ったのかということ。
1つ目の”newdata”に関する疑問は、RDocumantationに書いてあった。”newdata”に何も与えていないと学習に使用した値が返ってくるとのこと。今回233件のデータが返ってきていたのはこれが起こっていたからだった。
If missing, the fitted values are returned.
predict.rpart function - RDocumentation
2つ目の”data”に関する疑問。これもRDocumentationに書いてあって「…」という引数が元凶だった。AdvancedRでは”dot-dot-dot”という名称で説明されている。
このdot-dot-dotは追加の引数をいくらでも受け取ることができる。主に他の関数で利用されることが想定される関数に、後々柔軟に引数を取れるよう設定される様子。一方でその柔軟性故に、引数名を誤った場合に警告やエラーなく受け取ってしまうことが注意点としてAdvancedRには挙げられている。今回は件数が明らかに違ったので気づくことができたが、関数によっては難しそう。AdvancedRではsum関数を例にとってどのような問題が起こるのか説明されている。
Functions can have a special argument ... (pronounced dot-dot-dot). With it, a function can take any number of additional arguments. 6 Functions | Advanced R
教訓として2つ。
- 引数のスペルを確認する。引数にdot-dot-dotがある場合特に注意。
- Documentちゃんと読む。
この記事の内容をスライドにしたのが以下のリンク先。
GoogleAnalytics4のBigQueryテーブルが作成される時間の分布
GoogleAnalytics4のBigQuery連携で日次転送を選択した場合、テーブルが作成される時間は一定ではない。
作成時間について公式ドキュメントでも見つけられなかったので、自分の環境で何時に作成されているのか分布を確認してみる。
191日分のテーブルしかないので参考まで。
メタ情報からテーブル作成時間を集計するクエリは以下。
#standardSQL WITH tbl AS ( SELECT EXTRACT(HOUR FROM TIMESTAMP_MILLIS(creation_time) AT TIME ZONE "Asia/Tokyo") AS creation_hour FROM `project_id.analytics_XXXXX.__TABLES__` ) SELECT creation_hour, COUNT(*) AS cnt FROM tbl GROUP BY creation_hour ORDER BY creation_hour
BigQueryで全テーブルのメタ情報を一括で取得する方法 | GMOアドパートナーズ TECH BLOG byGMO
集計した結果をDataPortalで可視化。横軸が時間。
7,8時台が多く、1~5時は自分の環境では今の所存在しない。
データサイズが大きいと作成時刻が遅くなるのかと思ったが、そうでもなさそう。
アクセスが少ないので影響が見てとれないのかもしれない。
縦軸が作成時刻で横軸はメタ情報にあるsize_bytes。
【2022年4月17日追記】
_TABLES_はパフォーマンスに問題があり、非推奨とのこと。（リンク先スライドの12ページ目）
データ管理に役立つメタデータに関する勉強会を社内外で開催しました - MonotaRO Tech Blog
