（修正あり）arulesの結果をdata.frameで探索する
R標準のデータフレーム＋dplyr等で、探索的にルール抽出⇔眺めるのにパッケージを作った。本家のarulesだとちょっとやりにくいなあと思っている人むけ。
使いかた
arules::inspect() の代わりに inspectDF() するだけ。
arules::DATAFRAME() でデータフレームとして出力できるとのこと。
抽出されたルールをdataframe化する関数としてはすでに arules::DATAFRAME() があると思いますが、新たに inspectDF を作られた理由ってありますでしょうか？
— データポエマー[,-5] (@bob3bob3) September 16, 2018
自作する前によく探すの、大事。
require(dplyr) require(arules) require(inspectDF) data(Groceries) params <- list(confidence=0.001, support=0.001, maxlen=7, minlen=2) glo.apriori <- apriori(Groceries, parameter = params, control = list(verbose=FALSE))
glo.inspectDF <- inspectDF(glo.apriori) glo.inspectDF %>% str #> 'data.frame': 40943 obs. of 8 variables: #> $ rule : chr "Rule 1" "Rule 2" "Rule 3" "Rule 4" ... #> $ LHS : chr "whole milk" "other vegetables" "other vegetables" "rolls/buns" ... #> $ RHS : chr "sparkling wine" "artif. sweetener" "bathroom cleaner" "nuts/prunes" ... #> $ n.lhs : num 1 1 1 1 1 1 1 1 1 1 ... #> $ support : num 0.00102 0.00102 0.00102 0.00102 0.00102 ... #> $ confidence: num 0.00398 0.00525 0.00525 0.00553 0.00583 ... #> $ lift : num 0.712 1.615 1.914 1.647 1.024 ... #> $ count : num 10 10 10 10 10 10 10 10 10 10 ...
取り出したルールは、グラフ描画用のユーティリティ関数で眺める。
set.seed(0) glo.inspectDF %>% arrange(support, confidence) %>% head(60) %>% plotRuleGraph()
インストール
GitHubから.
install.packages("devtools") # if you have not installed "devtools" package devtools::install_github("katokohaku/inspectDF")
ソースコードも同じ場所。
モチベーション
本家のsubset(), sort()とかでなく、通常のdata.frameの操作としてやりたい。
本家の inspect() がコンソールに全部出力するのがツライ
print(glo.apriori) #> set of 40943 rules rules.lhs %>% summary() #> rule LHS RHS n.lhs #> Length:8035 Length:8035 Length:8035 Min. :2.000 #> Class :character Class :character Class :character 1st Qu.:2.000 #> Mode :character Mode :character Mode :character Median :3.000 #> Mean :2.739 #> 3rd Qu.:3.000 #> Max. :5.000 #> support confidence lift count #> Min. :0.001017 Min. :0.01815 Min. : 0.4379 Min. : 10.00 #> 1st Qu.:0.001118 1st Qu.:0.14961 1st Qu.: 2.0673 1st Qu.: 11.00 #> Median :0.001322 Median :0.25175 Median : 2.6852 Median : 13.00 #> Mean :0.001740 Mean :0.30231 Mean : 2.9595 Mean : 17.11 #> 3rd Qu.:0.001830 3rd Qu.:0.41935 3rd Qu.: 3.5531 3rd Qu.: 18.00 #> Max. :0.022267 Max. :1.00000 Max. :22.1394 Max. :219.00
arulesは非常によくできたパッケージですが、高度にwrapされ過ぎていて、気軽に弄ろうと思うとツライ。本家のprint() も summary()も、探索的なルール分析にはちょっと不親切だと思いませんか。探索的なルール分析にはsubset(), sort()などが用意されているものの、表示はことごとく上記のそれです。ツライ。
また、本家の inspect() はprint()で返すため、コンソールにすべての出力が強制表示される。いつ終わるとも知れぬコンソール出力を眺めるの、ツライ。あと、標準では2番目の列名が無名なので、この後そのまま探索的に弄ろうとすると困る、というか矢印要る？ ツライ。
glo.inspect <- glo.apriori %>% head(10) %>% inspect() #> lhs rhs support confidence lift #> [1] {honey} => {whole milk} 0.001118454 0.733333333 2.870009 #> [2] {whole milk} => {honey} 0.001118454 0.004377238 2.870009 #> [3] {soap} => {whole milk} 0.001118454 0.423076923 1.655775 #> [4] {whole milk} => {soap} 0.001118454 0.004377238 1.655775 #> [5] {tidbits} => {soda} 0.001016777 0.434782609 2.493345 #> [6] {soda} => {tidbits} 0.001016777 0.005830904 2.493345 #> [7] {tidbits} => {rolls/buns} 0.001220132 0.521739130 2.836542 #> [8] {rolls/buns} => {tidbits} 0.001220132 0.006633499 2.836542 #> [9] {cocoa drinks} => {whole milk} 0.001321810 0.590909091 2.312611 #> [10] {whole milk} => {cocoa drinks} 0.001321810 0.005173100 2.312611 #> count #> [1] 11 #> [2] 11 #> [3] 11 #> [4] 11 #> [5] 10 #> [6] 10 #> [7] 12 #> [8] 12 #> [9] 13 #> [10] 13 glo.inspect %>% str #> 'data.frame': 10 obs. of 7 variables: #> $ lhs : Factor w/ 7 levels "{cocoa drinks}",..: 2 7 4 7 6 5 6 3 1 7 #> $ : Factor w/ 1 level "=>": 1 1 1 1 1 1 1 1 1 1 #> $ rhs : Factor w/ 7 levels "{cocoa drinks}",..: 7 2 7 4 5 6 3 6 7 1 #> $ support : num 0.00112 0.00112 0.00112 0.00112 0.00102 ... #> $ confidence: num 0.73333 0.00438 0.42308 0.00438 0.43478 ... #> $ lift : num 2.87 2.87 1.66 1.66 2.49 ... #> $ count : num 11 11 11 11 10 10 12 12 13 13
使い道
dplyrあたりと組み合わせて使うと、tidyに探索できるのではないか。
require(stringr) rules.lhs <- glo.inspectDF %>% filter(str_detect(LHS, pattern = "yogurt|sausage")) %>% arrange(confidence, lift) %>% filter(n.lhs > 1) rules.lhs %>% head() #> rule LHS RHS n support #> 1 Rule 97 whole milk,yogurt UHT-milk 2 0.001016777 #> 2 Rule 98 whole milk,yogurt red/blush wine 2 0.001016777 #> 3 Rule 99 whole milk,yogurt house keeping products 2 0.001016777 #> 4 Rule 100 whole milk,yogurt liver loaf 2 0.001016777 #> 5 Rule 7175 whole milk,yogurt chewing gum 2 0.001118454 #> 6 Rule 7176 whole milk,yogurt cling film/bags 2 0.001118454 #> confidence lift count #> 1 0.01814882 0.5425339 10 #> 2 0.01814882 0.9444108 10 #> 3 0.01814882 2.1767518 10 #> 4 0.01814882 3.5698730 10 #> 5 0.01996370 0.9485170 11 #> 6 0.01996370 1.7530626 11
プロットは hideaki さんのqiitaの記事をそのまま流用。エッジ上の円の面積はルールのsupport値、色の濃さはルールのconfidence値を表す。
adujust.support.sizeオプションでヒューリスティックなサイズの調整ができる。
require(stringr) rules.lhs <- glo.inspectDF %>% filter(str_detect(RHS, pattern = "yogurt")) %>% arrange(confidence, lift) %>% filter(n.lhs > 1) %>% head(15) par(mfrow = c(1,2)) set.seed(0) rules.lhs %>% plotRuleGraph(label = "default") set.seed(0) rules.lhs %>% plotRuleGraph(label = "adjusted rule size", adujust.support.size = 4000)
今後
{support, confidence, lift}と{size, color}の対応は可変にして、サイズも自動調整したい。
既存のGithubのレポジトリをRStudioでcloneする
RStudioの操作だけで完結したい人向け。タイトル通りのことしたくて戸惑ったので備忘録。
(1) File > New Project... > Version Control > Git
(2) 既存のレポジトリのアドレスをコピー
(3) コピーしたアドレスを(1)に張り付け
(4) clone完了
あとは、ローカルの変更をpushしたり、別の端末からpushされた変更をpullしたり。(1)からVersion Controlを選べばいいというのに気が付くまで、時間がかかった。
元ネタ？はこちら：
qiita.com
random forestを使った解釈性にもとづく推薦システム
[1706.06691] Interpretable Predictions of Tree-based Ensembles via Actionable Feature Tweakingという論文が提案されている。KDD2017でYahoo Researchの研究者によって発表されたもの、とのこと。Ensemble treesによって、あるサンプルがラベルAと予測されたとき、別のラベルBと予測されたかったら元のサンプルのどの特徴量にどのような操作を加えればよいか「改善案」を算出する。手元にあると便利そうだったのでRで実装したので紹介したい。
どんなことができるかの概要は、先日のTokyoRで喋ってきたので、そちらのスライドも参照して欲しい。
大雑把な手順は、下記のとおり：
- データのスケーリング
- 予測モデル（random forest）の構築
- ルール抽出
- 各事例への推薦の計算
なお、日本語の解説は、接点QBさんによる丁寧な記事がすでにあるので、そちらを参照されたい。
setten-qb.hatenablog.com
前準備
ライブラリとソースのロード。この記事で使用するコードはgithubに置いてあるので、適当にローカルにおいて読み込む。そのうち時間を作ってパッケージ化したい。
github.com
require(tidyverse) require(magrittr) require(randomForest) devtools::install_github("hoxo-m/pforeach") require(pforeach) source("./R/tweak_feature.R")
pforeachに隅々まで依存しているので、別途インストールが必須。
訓練セットとテストセットの準備
例によってirisを使いたいところだが、今回はspamデータセットを使ってspamと分類されるサンプルをnonspamに分類されるように推薦を行う。irisは多クラスだしルールも最小限なんで良いんだけど、アヤメの花びらを切って貼ってもカキツバタにはならないので今回はスキップ。実用上、変更できない説明変数が予測モデルに含まれる場合には、禁止リストをかませるとか適切な層別化するといった工夫が別途必要になる。
# data preparation ------------------------------------------------------------- set.seed(777) data(spam, package = "kernlab") dataset <- sample_frac(spam) n.test <- floor(NROW(dataset) *0.1) dataset.train <- chop(dataset, n.test) dataset.test <- tail(dataset, n.test) dim(dataset);dim(dataset.train);dim(dataset.test)
chop()はデータを後ろから指定したサイズだけ切り落とすための自作ユーティリティ。
random forestの学習とモデル選択
訓練セットで学習して重要度のプロットと弱学習器の数の評価を行う。この辺りはルーチンなので変わったことはしていないが、ある程度予測精度が担保されていることが前提なので、ここで予測精度がダメな場合は再検討。
# bilding randomForest ----------------------------------------- X <- dataset.train[, 1:(ncol(dataset.train)-1)] true.y <- dataset.train[, ncol(dataset.train)] forest.all <- randomForest(X, true.y, ntree=500) forest.all par(mfrow=c(1,2)) varImpPlot(forest.all) # to view varImp, x3 & x4 should be removed. plot(forest.all) par(mfrow=c(1,1))
この方法は総当たり探索的なアプローチなので、変数（ルール）の数とツリーの数が増えると計算量が爆発的に増える。したがって変数選択とツリーの絞り込みのほか、max_depthを減らす工夫が重要。
今回は上位12変数と100ツリーくらいあれば良さそうなので、適当に切り落とす。
# model shrinkage based on importance ------------------------- top.importance <- forest.all$importance %>% data.frame %>% tibble::rownames_to_column(var = "var") %>% arrange(desc(MeanDecreaseGini)) %>% head(12) dataset.train.fs <- dataset.train %>% select(top.importance$var) dataset.test.fs <- dataset.test %>% select(top.importance$var)
データのスケーリングと縮約済みモデルの学習
この手法はハイパーパラメータを一つだけにするため、各変数を標準化する必要がある。実用上、テストセットは訓練セットのスケーリングパラメータで標準化される必要があるのと、推薦後の改善は標準化前の尺度に復元する必要があり、毎回attributeを取ってくるのは面倒なので、それぞれrescale()とdescale()として、ユーティリティ関数を用意した。
# scaling feature-selected data --------------------------- X.train <- scale( dataset.train.fs ) X.test <- rescale( dataset.test.fs, scaled = X.train ) dataset.test.fs[1:6, 1:6] descale(X.test, scaled = X.train)[1:6, 1:6] descale(X.test, scaled = X.test)[1:6, 1:6]
ちゃんと復元されていることを確認して欲しい。スケーリングした訓練データでモデルの学習を行う。
forest.rf <- randomForest(X.train, true.y, ntree=100) forest.all forest.rf plot(forest.rf)
変数を減らしたので予測精度が下がっているが、「予測ミス（擬陽性）を起こすサンプルは本質的には改善されて良いのでは」という立場をとって許容範囲と思うことにして目をつぶることにする。
feature tweaking
ここからが提案手法の処理。手順は、random forest予測モデルを学習したのち、（１）ルールを抽出する＋変更ルールを算出する（２）個別のサンプルの予測と推薦案を算出する（３）標準化の復元、の1+3ステップとなる。
# （０）random forest予測モデルを学習する forest.rf <- randomForest(X.train, true.y, ntree=100) # （１）ルールを抽出する＋変更ルールを算出する es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE) # （２）個別のサンプルの予測と推薦案を算出する tweaked <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = TRUE) # （３）標準化の復元 dt <- descale.tweakedFeature(tweaked, X.test)
（１）ルールを抽出する＋変更ルールを算出する
random forestではgetTree()で弱学習器を取ってこれる。これをつかって各ラベルに分類されるルールセットは、ラベルごとにdata.frameのlistとして整理しなおして抽出する。
# test sampling (for demo) ep <- getRules.randomForest(forest.rf, k=2, label.to = NULL) ep %>% str(2) ep[["spam"]]$path[[1]]
xgboostにも同様にルールセットを取り出すxgb.model.dt.tree()があるので、そのうち対応したい。
抽出されたルールセットと変更ルール（候補案）のセットは対応しているため、抽出と同時に計算してしまう。元論文のアルゴリズムは各サンプルの評価時にこれを毎回計算する構造になっており、素直に実装してしまうと計算量が大変になるので、今から実装しようと思う人は注意。
epsiron は変更の強さを指定するハイパーパラメタ。原理上、どうやっても予測を覆せないサンプルがあり、この量は変更可能な提案が可能なサンプル数に影響する。
es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3)
ntreeはルールセットを抽出する弱学習器の数で、省略するとすべてのツリーを使用する。デフォルトでは先頭からntreeで指定したツリーを解析対象にするが、データ不均衡に対処するために訓練セットをアンダーサンプリングして学習した幾つかのモデルをrandomForest::combine(...)で結合したい場合がある。この場合は、先頭から順番に解析してしまうと困るので、resampleオプションをTRUEに指定することでランダムに弱学習器を評価することができる。
es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE)
（２）個別のサンプルの予測と推薦案を算出する
予測対象のデータセットが label.from と予測されたとき、label.to と予測されるための改善案を計算する。学習データと同じスケールで標準されていないとマトモな推薦が出来ないので注意。
tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = FALSE)
各サンプルについてルールを全部評価するので結構な計算時間がかかるため、並列化オプション.doparはデフォルトではTRUEになっている。
tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = TRUE)
（３）標準化の復元
元のサンプルの操作方向とその量が知りたいのでスケーリングを解除する。元のサンプルの特徴量(original)、予測を変えるための変更後の特徴量(suggest)、予測の変更に必要な努力量=両者の差分(diff)のlistを返す。
>|r|
dt <- descale.tweakedFeature(tweaked, X.test)
dt %>% str(1)
|
個別のサンプルへの推薦パターンの可視化
個別の推薦結果を可視化する。操作量を相対的に比較したいため、標準化されたままのものを可視化している。サンプルによって違う推薦パターンが見えるので面白い。
# plot suggestion for each instance which(tweaked$predict == "spam") plot.suggest(tweaked, 4) plot.suggest(tweaked, 11)
操作量の大きい順に、0でない変数だけ可視化したければ、それぞれ .ordered = TRUEと、.nonzero.only = TRUEを指定する。
plot.suggest(tweaked, 15) plot.suggest(tweaked, 15, .ordered = TRUE, .nonzero.only = TRUE)
変数重要度の可視化
論文中ではユーザーによる納得感（アンケート調査?）をもとに評価しているが、それだと機械的には評価できないので、「予測サンプル集団において操作回数が多かった変数は重要なんじゃないか」という考えと、「操作された変数の操作の方向性が揃っている変数は、その意味合いを解釈できそう」と考えて、これらを各変数の重要度と思うことにして、それぞれ可視化する。
# Plot population importances pp <- plot.tweakedPopulation(tweaked, type = "absoluteSum") pp <- plot.tweakedPopulation(tweaked, type = "direction") pp <- plot.tweakedPopulation(tweaked, type = "frequency")
type = "frequency"は、サンプル集団において各変数が操作をされた集計値を可視化する。
type = "direction"は、各サンプルの推薦で各変数がどのような操作をされたかを集計して可視化する。相対的な重要性を比較したいので、標準化された状態の推薦量を使用する。
変更量の総和を可視化する"absoluteSum"オプションも作ってみたが、これはあまり有意義と感じなかった。
Lassoの非ゼロ変数を代替する候補を可視化する
最近、モデルの解釈性（interpretability）について色々と研究が出ており、興味もあってぼんやりと追いかけている。
“Finding Alternate Features in Lassos”という論文が提案されている：
Satoshi Hara and Takanori Maehara (2016): "Finding Alternate Features in Lasso", in NIPS 2016 workshop on Interpretable Machine Learning for Complex Systems, December 10th, 2016. （arXiv, github）
著者の日本語ページにもコンパクトな解説が掲載されている。
Lassoには変数間の相関が高いとき、それらの変数の中から一つだけを選択したがる性質がある(「見落としの問題」と呼ぶそうだ）。この問題に対して、この論文では、LASSOによって選択された変数がもし使われなかったとしたら、かわりにどの変数が選択されるかと考える。著者らによると、代替となる変数を列挙提示することで、納得感や結果の解釈性が得られるのではないか、と期待している。顧客に選択肢を残しておくというのは良いんじゃないだろうか。
具体的な手続きはシンプルで、Lassoで選択された係数の推定量が非ゼロの変数
と推定量ゼロの変数
の組を考え、
のもとで全ての組
について
、
とおいて目的関数を解きなおす。
と再推定されたら、「
は
の代替となる変数である」と提示され、やっぱり
なら「
は
の置き換えにはならない」とする。
また、代替となる変数同士の優劣は、目的関数の変化量を指標スコアとして相対的に順位づけされる。
で、ざっと論文を読んだ後、著者のPythonコードを見たらRでも簡単にできそうだったので、シンプルなOLS + soft-thresholdな罰則項という最小限だけをＲで実装した(github)。
ローカルに持ってきてDESCRIPTIONを用意してビルドすればパッケージになる（はず）。二分グラフを縦にプロットする方法で地味に手こずったが、この論文の可視化以外にもどこかで役に立つだろう。
面白いなと思った点は、目的関数をタスクに合わせてフレキシブルに設計できるところ。実際に、論文では損失関数をloglossに差し替えて二値分類タスクに適用・実験している。正則化項のほうを弄れば、例えばadaptive Lassoなどへの拡張もできるのではないかと思う。
また、各変数の相関行列と代替性の指標スコアとを比較すると、必ずしも大きさの順番がが一致しないのも興味深い。相関の高い変数が上位に来るのはわかりやすい一方、罰則の強さ(）の選び方によっては、割と相関の高くない変数でも代替候補として挙げられてくる。これはもしかしたら、前処理としてマルチコを除く作業をこの手法によってモデル構築後の二次工程に回せるかもしれないし、データによっては発見的な議論につながるかもしれない。
の選び方について、スライドの例ではcross-validationで求めた値を使った。しかし、変数が増えていくと見なければいけない情報が多すぎるので、強めの罰則をかけてエッセンシャルな変数だけに絞って分析するほうがシンプルかな、とは思った。この辺りは好みもあるかもしれないが、元の論文中では推奨値に関する議論は特に言及がなかったと思う。
代替候補の変数は必ずオリジナルの変数よりも目的関数の評価が下がることもあり、おそらく予測はオリジナルののモデルでやったほうがよい。そういう理由でprediction()は実装しなかった。この方法で予測したいとか2値分類について同様の分析がしたい方は本家のPython版を利用されたい。
推定値がズレる縮小推定の問題や、サンプルサイズよりも変数の数が大きいときに見過ごすといった、Lassoの癖というか性質に対処するアプローチの一例については、以前紹介したことがあるのでそちらも参照されたい。
mxnetパッケージをUbuntuで使うためのセットアップ
Ubuntu14.04でRからmxnet使いたくなってセットアップ。ガイド見ながらやっていけば終わり...と思ったらinstall_depsでいくつか引っかかった。Windows/Macばりに甘やかされたい。
次からコピペでやりたいので自分用にメモしておく。
http://mxnet.readthedocs.io/en/latest/how_to/build.html#r-package-installationmxnet.readthedocs.io
# install mxnet [~]$ sudo apt-get install build-essential git libatlas-base-dev libopencv-dev [~]$ git clone --recursive https://github.com/dmlc/mxnet [~]$ cd mxnet [~/mxnet]$ make
今回、makeは特に失敗せず完了。makeがおわったら以下を実行する。インストールガイドどおりにコマンドラインからやっても良い。
require(devtools); require(methods); require(httr) setwd("mxnet/R-package/") set_config(use_proxy(url="プロキシ使ってれば指定",port=8080)) options(repos=c(CRAN='https://cran.rstudio.com')); install_deps(dependencies = TRUE)
install_depsが動作するディレクトリに注意。今回はプロキシ指定とかが(なぜか)上手くいかなかったので、コマンドラインじゃなくスクリプトであれこれ試行錯誤した。
いろいろ入れてる途中でimagerでfftw3.hがないよと言われて失敗したので追加。
[~]$ sudo apt-get install fftw3 fft3-dev install_depsが沈黙したら下記を実行。 >|sh| [~]$ cd mxnet [~/mxnet]$ make rpkg
あとはインストール。
[~]$ R CMD INSTALL mxnet_0.7.tar.gz
以上をシェルスクリプトにまとめれば終わり。
random forestの感度分析＋
教科書を読むと、random forestでモデルを学習させて、ある良い予測精度のものが得られたら、まずは各変数の重要度をみて、つぎにpartial dependency partial dependence plot（PDP）を眺めなさい、と書いてある。ある変数の動きに注目して予測値がどう変化するか見てみよう、という分析だ。
これを発展（？）させたforestFloorは、決定木の各ノードでの予測値の変化量に着目して、random forestにおける変数が動いたとき予測値がどちら向きにどのくらい動くか定量化しようというもの。
スライド中のコードは、ほとんどマニュアルの丸写しなので割愛する。
forestFloorでは、ある決定木でのある事例の予測値（もしくは、あるクラスへの所属確率）＝ ブートストラップ平均 ＋ ノード1での変化量 ＋ ... ＋ 終端ノードの変化量、と解釈する。この解釈の下で、ブートストラップによって当該の観察事例が出てこない木や、観察事例が終端ノードに到達するまでに注目する変数を使わなかった、というケースは除外して、注目する変数による予測値の変化量をすべての木について足し合わせることで、ある観察事例に対するある変数の貢献度(feature contribution)を計算する。
このアイデアが面白いのは、回帰(→予測値)であろうが分類(→あるクラスへの所属確率)であろうが、弱学習器に決定木を使っているアンサンブル手法であれば貢献度が計算できる点だと思う。 弱点は、論文を眺めた限りでは、他の変数の影響を排除できているのかよく判らない点と、「それで結局予測はどうなるんだ」と訊かれそうな点と、原理上アンサンブルツリー縛りがある点。
とくに後者２つで、一般性(や有用性)という観点からPDP＆ICEboxに譲ってしまうと思われる。Hastieたちの教科書では、partial dependenceは任意の教師あり学習に対して、ある変数についてのmarginal averageが計算できると説明されている。これに加えて、ICEBoxを使うと予測値がの動きを分解(disaggregate)でき、かつ、予測値が大きく変化する変数値のポイントを検出できる。一応、forestFloorでも似たようなことはできるし、論文でもICEboxについて言及してはいるのだが、若干歯切れの悪い流れであった。。
ちなみに、これらはmlRにて任意の教師あり学習に適用可能なwrapperとして実装されていて、下記のリンク先ではSVMをつかった謎のデータirisの3クラス分類での例が示されている。
Partial Prediction Plots - mlr tutorial
(追記)
上記のリンクがいろいろと切れてると思ったら、mlr 2.8→2.9で関数名の変更があったようだ。
LassoとOracle Property
「Oracle Property」。巷にこれほどかっこいいフレーズはそうそうないと思うので紹介を試みる。
さて、Lassoの罰則項はどういう作用を持っているか？ L1 normの縮小推定を解説する図はいろんな教科書にも載っているのだが、自分は最初見たときに実感を持てなかった。Rで実際に条件を変えながらOLSと比較することで、一助になれば幸いである。
www.slideshare.net
今回使用したコード。 foreachのブロックは後半で結構時間がかかるので、イテレーションを減らすか並列化されたい。
# comparison among lm, lasso, adaptive lasso, and rigorous lasso require(dplyr) require(foreach) #install.packages("glmnet") require(glmnet) #install.packages("msgps") require(msgps) #install.packages("hdm") require(hdm) # functions --------------------------------------------------------------- # generate data according to model: y ~ 0 + B_i*x_i + N(0,1) Data_with <- function(n,p,s, silent=TRUE){ # data with sample size(n), number of variables(p) and non-zero vector(s) if(! silent){ cat(paste0(n, " obs. of ", p, " variables.\n", "Last ", length(s)," variables are non-zero.")) } beta <- c(rep(0,p-length(s)),s) X <- matrix(rnorm(n*p), ncol=p) res <-list( X = X, y = 0 + X%*%beta + rnorm(n), beta = beta, non.zero.coeff = which(beta>1e-8), non.zero.value = s ) if(! silent){ str(res) } invisible(res) } # effect of shrinkage on estimates (lm vs lasso) ---------------------------------- test_data <- Data_with(n=10^3, p=10, s= c(1,10,100)) str(test_data,max.level = 2) attach(test_data) system.time( lm.reg <- lm(y~X) ) summary(lm.reg) system.time( lasso_cv <- cv.glmnet(x = X, y=y, alpha = 1) ) coef(lasso_cv, s="lambda.1se") test_data <- Data_with(n=10^3, p=10, s= c(100,100,100)) str(test_data,max.level = 2) attach(test_data) system.time( lm.reg <- lm(y~X) ) summary(lm.reg) system.time( lasso_cv <- cv.glmnet(x = X, y=y, alpha = 1) ) par(mfrow=c(1,2)) plot(lasso_cv) abline(v=log(lasso_cv$lambda.1se), lty=2, col="blue") # plot solution path (Lasso) ------------------------------------------------------ lasso.reg <- glmnet(x = X, y=y, alpha = 1) plot(lasso.reg, xvar="lambda", label=TRUE) abline(v=log(lasso_cv$lambda.1se), lty=2, col="blue") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) par(mfrow=c(1,1)) fit.coeff <- coef(lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x) coef(lasso_cv, s="lambda.1se") fit.coeff %>% str # using lasso with oracle property (adaptive lasso) ------------------------- system.time( alasso.reg <- msgps(X= X, y=c(y), penalty="alasso",gamma=1, lambda=0) ) par(mfrow=c(2,2)) plot(alasso.reg, stand.coef=FALSE, criterion="cp", main="Mallows' Cp") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) plot(alasso.reg, stand.coef=FALSE, criterion="aicc",main="bias-corrected AIC (AICc)") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) plot(alasso.reg, stand.coef=FALSE, criterion="gcv", main="generalized cross validation (GCV) ") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) plot(alasso.reg, stand.coef=FALSE, criterion="bic", main="BIC") abline(h= c(1,2,3), lty=2, col= c(1,2,3)) par(mfrow=c(1,1)) summary(alasso.reg) ss=300 # measure: estimate error vs sample size ------------- sample_sizes <- c(5:20*10, 400, 600, 800, 1000,2000,3000, 4000) itr <- 1:10 conditions <- expand.grid(ss=sample_sizes, itr=itr) %>% arrange(ss) sse.sampsi <- foreach(Ci=1:NROW(conditions), .combine = rbind) %do% { cond <- conditions[Ci,] this_data <- Data_with(n=cond$ss, p=10, s= c(10,10,10), silent=TRUE) beta <- this_data$beta fit.lm <- lm(this_data$y~this_data$X) sse.lm <- sum( (beta - coef(fit.lm)[-1])^2 ) fit.lasso_cv <- cv.glmnet(x = this_data$X, y=this_data$y, alpha = 1) coeff.lasso_cv <- coef(fit.lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x) sse.lasso_cv <- sum( (beta[coeff.lasso_cv$beta[-1]] - coeff.lasso_cv$Var[-1])^2 ) fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1, lambda=0) coeff.alasso <- coef(fit.alasso) %>% data.frame %>% select(BIC) %>% unlist sse.alasso <- sum( (this_data$beta - coeff.alasso[-1])^2 ) data.frame(sample.size=cond$ss, lm=sse.lm, lasso=sse.lasso_cv, alasso=sse.alasso) } sse.summ.sampsi <- sse.sampsi %>% group_by(sample.size) %>% summarize_each(funs(mean)) sse.summ.sampsi with(sse.summ.sampsi,{ plot (lm ~sample.size, type="l", lwd=2, col="green", ylim=c(0,0.35), ylab="mean SSE", log="x") lines(lasso~sample.size, type="l", lwd=2, col="red") legend("topright", legend = c("lm","lasso"), col = c("green","red"), lty = 1, lwd=2) }) with(sse.summ.sampsi,{ plot (lm ~sample.size, type="l", lwd=2, col="green", ylim=c(0,0.35), ylab="mean SSE", log="x") lines(lasso ~sample.size, type="l", lwd=2, col="red") lines(alasso~sample.size, type="l", lwd=2, col="blue") legend("topright", legend = c("lm","lasso","a-lasso"), col = c("green","red","blue"), lty = 1, lwd=2) }) # example(rlasso) --------------------------------------------------------- # use Lasso with rigorous penalty system.time( nonpost_lasso.reg <- rlasso( formula, data, post = FALSE, # FALSEで、罰則だけrigorousなLasso intercept = TRUE, # TRUEで切片項(β0）の罰則付き推定をしない penalty = list(homoscedastic = FALSE, # FALSEで誤差の等分散性を仮定しない X.dependent.lambda = FALSE, # TRUEでデータ依存的なλの選択をする lambda.start = NULL, c = 1.1, gamma = 0.1/log(n)), control = list(numIter =15, tol = 10^-5, threshold = NULL), ...) ) # use Post-Lasso procedure with rigorous penalty system.time( post_lasso.reg <- rlasso( formula, data, post = TRUE, # TRUEでPost-Lasso＋rigorousな罰則 intercept = TRUE, # TRUEで切片項(β0）の罰則付き推定をしない penalty = list(homoscedastic = FALSE, # FALSEで誤差の等分散性を仮定しない X.dependent.lambda = FALSE, # TRUEでデータ依存的なλの選択をする lambda.start = NULL, c = 1.1, gamma = 0.1/log(n)), control = list(numIter =15, tol = 10^-5, threshold = NULL), ...) ) # comparison results (set all=FALSE to summarize) summary(nonpost_lasso.reg, all=FALSE) summary(post_lasso.reg, all=FALSE) coef(post_lasso.reg) # measure: estimate error vs dimension size ------------- dim_sizes <- c(1:10*10, 1:5*200, 2:5*1000, 1e4, 2e4) itr <- 1:30 conditions <- expand.grid(ds=dim_sizes, itr=itr) %>% arrange(ds) sse.dim_sizes <- foreach(Ci=1:NROW(conditions), .combine = rbind) %do% { cond <- conditions[Ci,] this_data <- Data_with(n=1000, p=cond$ds, s= c(10,10,10), silent=TRUE) beta <- this_data$beta fit.lm <- lm(this_data$y~this_data$X) sse.lm <- sum( (beta - coef(fit.lm)[-1])^2 ) fit.lasso_cv <- cv.glmnet(x = this_data$X, y=this_data$y, alpha = 1) coeff.lasso_cv <- coef(fit.lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x) sse.lasso_cv <- sum( (beta[coeff.lasso_cv$beta[-1]] - coeff.lasso_cv$Var[-1])^2 ) fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1) coeff.alasso <- coef(fit.alasso) %>% data.frame %>% select(BIC) %>% unlist sse.alasso <- sum( (this_data$beta - coeff.alasso[-1])^2 ) plasso.reg <- rlasso(c(this_data$y) ~ this_data$X, post=TRUE, intercept=TRUE) coeff.plasso <-coef(plasso.reg) sse.plasso <-(this_data$beta - coeff.plasso)^2 %>% sum data.frame(dim.size=cond$ds, lm=sse.lm, lasso=sse.lasso_cv, alasso=sse.alasso, plasso=sse.plasso) } sse.summ.dim_sizes <- sse.dim_sizes %>% group_by(dim.size) %>% summarize_each((funs(mean))) sse.summ.dim_sizes with(sse.summ.dim_sizes,{ plot( lm ~dim.size, type="l", lwd=2, col="green", ylab="mean SSE", log="x") lines(lasso ~dim.size, type="l", lwd=2, col="red") lines(alasso~dim.size, type="l", lwd=2, col="blue") legend("topright", legend = c("lm","lasso","a-lasso"), col = c("green","red","blue"), lty = 1, lwd=2) abline(v=1000, lty=3) }) with(sse.summ.dim_sizes,{ plot(lasso~dim.size, type="l", lwd=2, col="red", ylim=c(0,0.1), ylab="mean SSE",log="x") lines(alasso~dim.size, type="l", lwd=2, col="blue") legend("topleft", legend = c("lasso","a-lasso"), col = c("red","blue"), lty = 1, lwd=2) abline(v=1000, lty=3) }) with(sse.summ.dim_sizes,{ plot(lasso~dim.size, type="l", lwd=2, col="red", ylim=c(0,0.1), ylab="mean SSE",log="x") lines(alasso~dim.size, type="l", lwd=2, col="blue") lines(plasso~dim.size, type="l", lwd=2, col="orange") legend("topleft", legend = c("lasso","a-lasso","post-lasso"), col = c("red","blue","orange"), lty = 1, lwd=2) abline(v=1000, lty=3) }) # measure: performance vs dimension size ------------- this_data <- Data_with(n=1000, p=2e4, s= rep(10,3),silent = FALSE) system.time( plasso.reg <- rlasso(c(this_data$y) ~ this_data$X, post=TRUE, intercept=TRUE) ) #coef(fit.lasso_cv) system.time( fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1) )
反省点その１、これほどかっこいいフレーズなのに当日聴衆には全く刺さらなかった（というか完全に滑った）。プレゼンスキルの改善がである。その２、圧縮センシングのお得感とか、スパース推定後の解釈性とか、もうちょっと使い道に触れた導入にすべきだった。最後の"rigorous lasso"では、計算時間は少なからず増えているわけで、実際の現場感としてここまでの精度が必要か？とか、いずれ費用対効果やロバスト化の議論はしてみたい。
