琥珀色呑んだくれ備忘録

メモとか備忘録とか

既存のGithubのレポジトリをRStudioでcloneする

RStudioの操作だけで完結したい人向け。タイトル通りのことしたくて戸惑ったので備忘録。

(1) File > New Project... > Version Control > Git
f:id:kato-satoshi-0:20180126003129p:plain

(2) 既存のレポジトリのアドレスをコピー
f:id:kato-satoshi-0:20180126003414p:plain

(3) コピーしたアドレスを(1)に張り付け
f:id:kato-satoshi-0:20180126003449p:plain

(4) clone完了
f:id:kato-satoshi-0:20180126003609p:plain

あとは、ローカルの変更を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で喋ってきたので、そちらのスライドも参照して欲しい。

www.slideshare.net
大雑把な手順は、下記のとおり:

  1. データのスケーリング
  2. 予測モデル(random forest)の構築
  3. ルール抽出
  4. 各事例への推薦の計算

なお、日本語の解説は、接点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))

f:id:kato-satoshi-0:20180122185241p:plain
この方法は総当たり探索的なアプローチなので、変数(ルール)の数とツリーの数が増えると計算量が爆発的に増える。したがって変数選択とツリーの絞り込みのほか、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)

f:id:kato-satoshi-0:20180122191903p:plain
変数を減らしたので予測精度が下がっているが、「予測ミス(擬陽性)を起こすサンプルは本質的には改善されて良いのでは」という立場をとって許容範囲と思うことにして目をつぶることにする。

feature tweaking

ここからが提案手法の処理。手順は、random forest予測モデルを学習したのち、(1)ルールを抽出する+変更ルールを算出する(2)個別のサンプルの予測と推薦案を算出する(3)標準化の復元、の1+3ステップとなる。

# (0)random forest予測モデルを学習する
forest.rf <- randomForest(X.train, true.y, ntree=100)

# (1)ルールを抽出する+変更ルールを算出する
es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE)

# (2)個別のサンプルの予測と推薦案を算出する
tweaked <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = TRUE)

# (3)標準化の復元
dt <- descale.tweakedFeature(tweaked, X.test)
(1)ルールを抽出する+変更ルールを算出する

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]]

f:id:kato-satoshi-0:20180122194147p:plain
xgboostにも同様にルールセットを取り出すxgb.model.dt.tree()があるので、そのうち対応したい。

抽出されたルールセットと変更ルール(候補案)のセットは対応しているため、抽出と同時に計算してしまう。元論文のアルゴリズムは各サンプルの評価時にこれを毎回計算する構造になっており、素直に実装してしまうと計算量が大変になるので、今から実装しようと思う人は注意。
epsiron は変更の強さを指定するハイパーパラメタ。原理上、どうやっても予測を覆せないサンプルがあり、この量は変更可能な提案が可能なサンプル数に影響する。

es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3)

f:id:kato-satoshi-0:20180122200255p:plain
ntreeはルールセットを抽出する弱学習器の数で、省略するとすべてのツリーを使用する。デフォルトでは先頭からntreeで指定したツリーを解析対象にするが、データ不均衡に対処するために訓練セットをアンダーサンプリングして学習した幾つかのモデルをrandomForest::combine(...)で結合したい場合がある。この場合は、先頭から順番に解析してしまうと困るので、resampleオプションをTRUEに指定することでランダムに弱学習器を評価することができる。

es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE)

f:id:kato-satoshi-0:20180122200705p:plain

(2)個別のサンプルの予測と推薦案を算出する

予測対象のデータセットが label.from と予測されたとき、label.to と予測されるための改善案を計算する。学習データと同じスケールで標準されていないとマトモな推薦が出来ないので注意。

tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam",
                 .dopar = FALSE)

f:id:kato-satoshi-0:20180122202613p:plain
各サンプルについてルールを全部評価するので結構な計算時間がかかるため、並列化オプション.doparはデフォルトではTRUEになっている。

tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam",
                 .dopar = TRUE)
(3)標準化の復元

元のサンプルの操作方向とその量が知りたいのでスケーリングを解除する。元のサンプルの特徴量(original)、予測を変えるための変更後の特徴量(suggest)、予測の変更に必要な努力量=両者の差分(diff)のlistを返す。
>|r|
dt <- descale.tweakedFeature(tweaked, X.test)
dt %>% str(1)
|

f:id:kato-satoshi-0:20180122204026p:plain

個別のサンプルへの推薦パターンの可視化

個別の推薦結果を可視化する。操作量を相対的に比較したいため、標準化されたままのものを可視化している。サンプルによって違う推薦パターンが見えるので面白い。

# plot suggestion for each instance
which(tweaked$predict == "spam")
plot.suggest(tweaked, 4)
plot.suggest(tweaked, 11)

f:id:kato-satoshi-0:20180122203751p:plain
f:id:kato-satoshi-0:20180122203745p:plain
操作量の大きい順に、0でない変数だけ可視化したければ、それぞれ .ordered = TRUEと、.nonzero.only = TRUEを指定する。

plot.suggest(tweaked, 15)
plot.suggest(tweaked, 15, .ordered = TRUE, .nonzero.only = TRUE)

f:id:kato-satoshi-0:20180122203928p:plain
f:id:kato-satoshi-0:20180122203934p:plain

変数重要度の可視化

論文中ではユーザーによる納得感(アンケート調査?)をもとに評価しているが、それだと機械的には評価できないので、「予測サンプル集団において操作回数が多かった変数は重要なんじゃないか」という考えと、「操作された変数の操作の方向性が揃っている変数は、その意味合いを解釈できそう」と考えて、これらを各変数の重要度と思うことにして、それぞれ可視化する。

# Plot population importances
pp <- plot.tweakedPopulation(tweaked, type = "absoluteSum")
pp <- plot.tweakedPopulation(tweaked, type = "direction")
pp <- plot.tweakedPopulation(tweaked, type = "frequency")

type = "frequency"は、サンプル集団において各変数が操作をされた集計値を可視化する。
f:id:kato-satoshi-0:20180122205329p:plain
type = "direction"は、各サンプルの推薦で各変数がどのような操作をされたかを集計して可視化する。相対的な重要性を比較したいので、標準化された状態の推薦量を使用する。
f:id:kato-satoshi-0:20180122210838p:plain
変更量の総和を可視化する"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で選択された係数\beta^*の推定量が非ゼロの変数X_iと推定量ゼロの変数X_jの組を考え、\beta^*のもとで全ての組(i,j)について\beta_i=0\beta_j\not=0とおいて目的関数を解きなおす。\beta_j\not=0と再推定されたら、「X_jX_iの代替となる変数である」と提示され、やっぱり\beta_j=0なら「X_jX_iの置き換えにはならない」とする。
また、代替となる変数同士の優劣は、目的関数の変化量を指標スコアとして相対的に順位づけされる。

で、ざっと論文を読んだ後、著者のPythonコードを見たらRでも簡単にできそうだったので、シンプルなOLS + soft-thresholdな罰則項という最小限だけをRで実装した(github)。

www.slideshare.net
ローカルに持ってきてDESCRIPTIONを用意してビルドすればパッケージになる(はず)。二分グラフを縦にプロットする方法で地味に手こずったが、この論文の可視化以外にもどこかで役に立つだろう。

面白いなと思った点は、目的関数をタスクに合わせてフレキシブルに設計できるところ。実際に、論文では損失関数をloglossに差し替えて二値分類タスクに適用・実験している。正則化項のほうを弄れば、例えばadaptive Lassoなどへの拡張もできるのではないかと思う。

また、各変数の相関行列と代替性の指標スコアとを比較すると、必ずしも大きさの順番がが一致しないのも興味深い。相関の高い変数が上位に来るのはわかりやすい一方、罰則の強さ(\lambda)の選び方によっては、割と相関の高くない変数でも代替候補として挙げられてくる。これはもしかしたら、前処理としてマルチコを除く作業をこの手法によってモデル構築後の二次工程に回せるかもしれないし、データによっては発見的な議論につながるかもしれない。
\lambdaの選び方について、スライドの例ではcross-validationで求めた値を使った。しかし、変数が増えていくと見なければいけない情報が多すぎるので、強めの罰則をかけてエッセンシャルな変数だけに絞って分析するほうがシンプルかな、とは思った。この辺りは好みもあるかもしれないが、元の論文中では推奨値に関する議論は特に言及がなかったと思う。

代替候補の変数は必ずオリジナルの変数よりも目的関数の評価が下がることもあり、おそらく予測はオリジナルの\beta^*のモデルでやったほうがよい。そういう理由でprediction()は実装しなかった。この方法で予測したいとか2値分類について同様の分析がしたい方は本家のPython版を利用されたい。

推定値がズレる縮小推定の問題や、サンプルサイズよりも変数の数が大きいときに見過ごすといった、Lassoの癖というか性質に対処するアプローチの一例については、以前紹介したことがあるのでそちらも参照されたい。

www.slideshare.net

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 plot(PDP)を眺めなさい、と書いてある。ある変数の動きに注目して予測値がどう変化するか見てみよう、という分析だ。

 これを発展(?)させたforestFloorは、決定木の各ノードでの予測値の変化量に着目して、random forestにおける変数が動いたとき予測値がどちら向きにどのくらい動くか定量化しようというもの。 

www.slideshare.net

スライド中のコードは、ほとんどマニュアルの丸写しなので割愛する。

 forestFloorでは、ある決定木でのある事例の予測値(もしくは、あるクラスへの所属確率)= ブートストラップ平均 + ノード1での変化量 + ...  + 終端ノードの変化量、と解釈する。この解釈の下で、ブートストラップによって当該の観察事例が出てこない木や、観察事例が終端ノードに到達するまでに注目する変数を使わなかった、というケースは除外して、注目する変数による予測値の変化量をすべての木について足し合わせることで、ある観察事例に対するある変数の貢献度(feature contribution)を計算する。

 このアイデアが面白いのは、回帰(→予測値)であろうが分類(→あるクラスへの所属確率)であろうが、弱学習器に決定木を使っているアンサンブル手法であれば貢献度が計算できる点だと思う。 弱点は、論文を眺めた限りでは、他の変数の影響を排除できているのかよく判らない点と、「それで結局予測はどうなるんだ」と訊かれそうな点と、原理上アンサンブルツリー縛りがある点。

 とくに後者2つで、一般性(や有用性)という観点から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で関数名の変更があったようだ。

Partial Dependence Plots - mlr tutorial

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)
)


反省点その1、これほどかっこいいフレーズなのに当日聴衆には全く刺さらなかった(というか完全に滑った)。プレゼンスキルの改善がである。その2、圧縮センシングのお得感とか、スパース推定後の解釈性とか、もうちょっと使い道に触れた導入にすべきだった。最後の"rigorous lasso"では、計算時間は少なからず増えているわけで、実際の現場感としてここまでの精度が必要か?とか、いずれ費用対効果やロバスト化の議論はしてみたい。

random forestを予測以外の目的で使う

数年ほど前には最強と言われて一世を風靡していたrandom forestだが、予測以外にも使い道が提案されている。Rのパッケージから紹介したい。

 予測全体の把握と仮説ルールの抽出

決定木分析が便利な理由の一つは「どういうルールでその予測が成り立っているのか」を極めて簡単に可視化出来る点。inTreesパッケージは予測ルールを集計し、適当に枝狩りして全体を要約することで、アンサンブルモデルにおいても決定木と同じような情報を可視化してくれる。

www.slideshare.net

このパッケージではもう一つ、すべての木から取り出した1つずつの枝をトランザクションとみなしてアソシエーション分析する機能を提供する。メジャーな関連ルールを取り出す方に力点が置かれているが、変数同士の稀な関係が意外なアウトカムを導いてる方が、(研究上は)見つかったときに面白いので、欲を言えばLift値も評価してくれると良かった。

ノンパラメトリックな欠損値の補完

randomForestパッケージには、rfImputeという欠損値を補完する関数が用意されている。ところが、rfImputeは同じルールにより同じ予測をされる別の観察データを使って補完するため、目的変数が欠損しているデータには適用できない。

既知のデータにおいて欠測部分を推定するにはこれで構わないのだが、予測をしたいデータの前処理に使おうとするとこれでは困る。

missForestでは、補完対象の変数を目的変数とみなすことで、random forestにより欠損値を直接予測する。

www.slideshare.net

もともと予測したい変数も補完プロセスに含めることができるので半教師学習と言える。その場合、イテレーションの回数と過学習の関係が気になるが、元の論文にはその言及はなかったように思う。