![Partial Dependency Plot (PDP) using randomForest
library(randomForest)
library(forestFloor)
#grow a forest, remember to include inbag
multireg.rfo=randomForest::randomForest(X,Y,
keep.inbag=TRUE,
ntree=1000, sampsize=500,
replace=TRUE, importance=TRUE)
names.X <- c("X2","X3","X1","X5","X6","X4")
# randomForest::partialPlot()
par(mfrow=c(2, 3))
for (i in seq_along(names.X)) {
partialPlot(multireg.rfo, X, names.X[i], xlab=names.X[i],
main = names.X[i], ylim=c(-4,10))
}
par(mfrow=c(1, 1))
注目する変数がある値のとき、
「平均的な予測値」を可視化する](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-8-320.jpg)
![# (続き)
library(rgl)
rgl::plot3d(ff$X[,2],ff$X[,3],apply(ff$FCmatrix[,2:3],1,sum),
#add some colour gradients to ease visualization
#box.outliers squese all observations in a 2 std.dev box
#univariately for a vector or matrix and normalize to [0;1]
col=fcol(ff,2,orderByImportance=FALSE))
#add grid convolution/interpolation
#make grid with current function
grid23 = convolute_grid(ff,Xi=2:3,userArgs.kknn= alist(k=25,kernel="gaus"),grid=50,zoom=1.2)
#apply grid on 3d-plot
rgl::persp3d(unique(grid23[,2]),unique(grid23[,3]),grid23[,1],alpha=0.3,
col=c("black","grey"),add=TRUE)
Feature Contribution
Rglパッケージに渡すことで変数間の相互作用も観察できる
X2の大きさでグラデーションを作ってプロット](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-19-320.jpg)
![Irisデータで3クラス分類
## Not run: example of plot_simplex3
library(randomForest)
library(forestFloor)
library(utils)
data(iris)
X = iris[,!names(iris) %in% "Species"]
Y = iris[,"Species"]
as.numeric(Y)
rf.test42 = randomForest(X,Y,
keep.forest=TRUE,
replace=FALSE,
keep.inbag=TRUE,
samp=15,
ntree=100)
ff.test42 = forestFloor(rf.test42,X,
calc_np=FALSE,binary_reg=FALSE)
plot(ff.test42,plot_GOF=TRUE,cex=.7,
colLists=list(c("#FF0000A5"),
c("#00FF0050"),
c("#0000FF35")))
show3d(ff.test42,1:2,3:4,plot_GOF=TRUE)
各変数の各クラスへの所属確率への貢献度](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-23-320.jpg)
![所属確率の可視化（3-class限定）
# つづき
pars = plot_simplex3(ff.test42,Xi=c(1:3),restore_par=FALSE,zoom.fit=NULL,
var.col=NULL,fig.cols=2,fig.rows=1,fig3d=FALSE,includeTotal=TRUE,auto.alpha=.4, set_pars=TRUE)
pars = plot_simplex3(ff.test42,Xi=0,restore_par=FALSE,zoom.fit=NULL,
var.col=alist(alpha=.3,cols=1:4),fig3d=FALSE,includeTotal=TRUE, auto.alpha=.8,set_pars=FALSE)
for (I in ff.test42$imp_ind[1:4]) {
#plotting partial OOB-CV separation(including interactions effects)
#coloured by true class
pars = plot_simplex3(ff.test42,Xi=I,restore_par=FALSE,zoom.fit=NULL,
var.col=NULL,fig.cols=4,fig.rows=2,fig3d=TRUE,includeTotal=FALSE,label.col=1:3,
auto.alpha=.3,set_pars = (I==ff.test42$imp_ind[1]))
#coloured by varaible value
pars = plot_simplex3(ff.test42,Xi=I,restore_par=FALSE,zoom.fit=TRUE,
var.col=alist(order=FALSE,alpha=.8),fig3d=FALSE,includeTotal=(I==4), auto.alpha=.3,set_pars=FALSE)](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-25-320.jpg)
![余談： グラデーション作成関数 fcol()
library(forestFloor)
X = data.frame(matrix(rnorm(1000), nrow=1000, ncol=4))
X[] = lapply(X,jitter,amount = 1.5)
plot(X)](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-26-320.jpg)
![余談： グラデーション作成関数 fcol()
library(forestFloor)
X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4))
X[] = lapply(X,jitter,amount = 1.5)
plot(X,col=fcol(X,1))
１番目の変数の勾配でカラー生成](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-27-320.jpg)
![余談： グラデーション作成関数 fcol()
library(forestFloor)
X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4))
X[] = lapply(X,jitter,amount = 1.5)
plot(X,col=fcol(X,2))
２番目の変数の勾配でカラー生成](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-28-320.jpg)
![余談： グラデーション作成関数 fcol()
library(forestFloor)
X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4))
X[] = lapply(X,jitter,amount = 1.5)
plot(X,col=fcol(X,1:2))
１＋２番目の変数の線形な勾配でカラー生成](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-29-320.jpg)
![余談： グラデーション作成関数 fcol()
library(forestFloor)
X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4))
X[] = lapply(X,jitter,amount = 1.5)
plot(X,col=fcol(X,1:3))
１～３番目の変数のPCA勾配でカラー生成](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-30-320.jpg)
![余談： グラデーション作成関数 fcol()
library(forestFloor)
X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4))
X[] = lapply(X,jitter,amount = 1.5)
plot(X,col=fcol(X,1:4))
１～４番目の変数のPCA勾配でカラー生成](https://image.slidesharecdn.com/sensitivityanalysisusingforestfloor-160730022428/85/forestFloor-randomForest-31-320.jpg)
forestFloorパッケージを使ったrandomForestの感度分析
- 1. ランダムフォレストで感度分析＋ “sensitivity analysis using forestFloor” forestFloor packageを紹介します (S. H. Welling, et al., ArXiv e-prints, June 2016.) 第55回R勉強会＠東京（#TokyoR） feature contribution forestFloorを使ってみる
- 2. 自己紹介 専門 • データ分析／シミュレーション／統計モデリング／数理モデリング • 生態学／環境科学 → 臨床検査／医療情報学／疫学 1cm
- 3. やりたいこと X1 X2 X3 X4 X5 X6 X7 X8 ： モデル C1 C2 C3 C4 C5 ：? （例） 変数X6 が増加したとき、 分類クラスC4への所属確率は増えるのか❓減るのか❓ X6 は、見ておく and/or 操作する価値がある変数か？ 注目した変数が予測にどう影響するか知りたい
- 4. ランダムフォレスト 学習データのランダムサブセットで構築した様々な決定木の集合（＝森）の 予測結果 を統合する 分類 → 多数決 回帰 → 平均 特定の説明変数への依存が少ないため、クエリデータの説明変数が 欠損していても良い出力を与える ALL DATA Random subset Random subset Random subset …
- 6. 変数重要度を求めたら、その次は… • In data mining applications the input predictor variables are seldom equally relevant. Often only a few of them have substantial influence on the response; the vast majority are irrelevant and could just as well have not been included. It is often useful to learn the relative importance or contribution of each input variable in predicting the response. • After the most relevant variables have been identified, the next step is to attempt to understand the nature of the dependence of the approximation f(X) on their joint values. in, Hastie, Tibshirani, Friedman (2008), ESLII(2nd) pp367- 変数が予測にどう影響するかを知りたい ⇒感度分析
- 7. 多変量回帰の例 デモ用データ 𝑌 = 2 𝑋1 + 2 sin 𝑝𝑖 𝑋2 + 3 𝑋3 + 𝑋2 2 + 𝜀 +0(X4 + 𝑋5 + 𝑋6) #simulate data obs=1500 vars = 6 X = data.frame(replicate(vars,runif(obs)))*2-1 Y = with(X, X1*2 + 2*sin(X2*pi) + 3* (X3+X2)^2 ) #X4, X5,X6 are noises Yerror = 1 * rnorm(obs) var(Y)/var(Y+Yerror) Y= Y+Yerror
- 8. Partial Dependency Plot (PDP) using randomForest library(randomForest) library(forestFloor) #grow a forest, remember to include inbag multireg.rfo=randomForest::randomForest(X,Y, keep.inbag=TRUE, ntree=1000, sampsize=500, replace=TRUE, importance=TRUE) names.X <- c("X2","X3","X1","X5","X6","X4") # randomForest::partialPlot() par(mfrow=c(2, 3)) for (i in seq_along(names.X)) { partialPlot(multireg.rfo, X, names.X[i], xlab=names.X[i], main = names.X[i], ylim=c(-4,10)) } par(mfrow=c(1, 1)) 注目する変数がある値のとき、 「平均的な予測値」を可視化する
- 9. PDP focuses on marginal averages ① 𝑥i がある値のとき、残りの変数の 違いによる y の値を平均する ② すべての 𝑥i のついて y の平均を 算出し、線でつなぐ
- 10. やりたいこと Partial dependence • 注目する変数 vs 「平均的な予測値」をプロット • 変数同士に相互作用がある場合、うまくいかない ※ PDPで十分、というデータも当然ある
- 11. • ランダムフォレストの感度分析に便利な指標を提供 • feature contributionによる『効果』の可視化 • 回帰だけでなく、2クラス分類・多クラス分類もｏｋ • （ Cross-Validationを使うための拡張方法も併せて提案 ） forestFloor: Visualizes Random Forests with Feature Contributions https://cran.r-project.org/web/packages/forestFloor/index.html
- 12. Feature Contribution Random forestにおける弱学習器の各ノードに着目することで、 ある説明変数が、ある観察値にどう寄与したか？を可視化する • 回帰なら予測値 • 分類ならクラスへの所属確率 観察値そのもの（PDP)ではなく、移動方向を可視化しているところがポイント library(randomForest) #grow a forest, remember to include inbag multireg.rfo=randomForest::randomForest(X,Y, keep.inbag=TRUE, ntree=1000, sampsize=500, replace=TRUE, importance=TRUE) names.X <- c("X2","X3","X1","X5","X6","X4") # disaggregation using ICEbox::ice() require(forestFloor) #compute ff multireg.ff <- forestFloor(multireg.rfo,X) #print forestFloor print(multireg.ff) #plot partial functions of most important variables first Col <- fcol(multireg.ff,1) plot(multireg.ff,col=Col,orderByImportance=TRUE) X2の値でカラー付け
- 15. 回帰木 ⇒ random forest回帰 回帰木 ある観察事例を、終端ノード（葉）によって、 ひとつの予測値にマップする random forest回帰 ある観察事例について、 各各弱学習器でマップされた予測値を平均する F1(Xi) → yi1 F2(Xi) → yi2 F3(Xi) → yi3 F4(Xi) → yi4 F5(Xi) → yi5 : ある観察事例の予測値 ｙ𝑖 = 1 𝑁 𝑚=1 𝑀 𝑦𝑖𝑚
- 16. 回帰木を解釈すると.. Xi 変数： x2 貢献= -7.4 ある葉（N5)に落ちた事例 Xi の観察値 y^ 5 ＝ 0.14 -7.4 +4.0 ＝-3.0 変数： x1 貢献= +4 事例 Xi の予測値＝各ノードでの増分 (increment)の総和で表現される -7.4 Xi -7.4 +4
- 17. ある観察事例に対する feature contribution Xi 変数： x2 貢献= -7.4 ある葉（N5)に落ちた事例 Xi の観察値に対する 変数x2の貢献は、 L = -7.4 事例 Xi の観察値に対して、ある変数による増分 (increment)を考える -7.4Xi
- 18. ある観察事例に対する feature contribution (random forest回帰) すべてのモデルについて、 事例 Xi の観察値に対する、変数 xl の貢献を足しあげる 𝐹𝑖𝑙 = 𝑗=1 𝑛 𝑡𝑟𝑒𝑒 𝑘⊆𝐻 𝑖𝑗𝑙 𝐿 𝑖𝑗𝑘 𝑛 𝑡𝑟𝑒𝑒 あるモデルにおける ある変数の貢献 Bootstrapの時点で事例 Xi が除外されたもの、変数 xl が登場しないものは計算されない
- 19. # (続き) library(rgl) rgl::plot3d(ff$X[,2],ff$X[,3],apply(ff$FCmatrix[,2:3],1,sum), #add some colour gradients to ease visualization #box.outliers squese all observations in a 2 std.dev box #univariately for a vector or matrix and normalize to [0;1] col=fcol(ff,2,orderByImportance=FALSE)) #add grid convolution/interpolation #make grid with current function grid23 = convolute_grid(ff,Xi=2:3,userArgs.kknn= alist(k=25,kernel="gaus"),grid=50,zoom=1.2) #apply grid on 3d-plot rgl::persp3d(unique(grid23[,2]),unique(grid23[,3]),grid23[,1],alpha=0.3, col=c("black","grey"),add=TRUE) Feature Contribution Rglパッケージに渡すことで変数間の相互作用も観察できる X2の大きさでグラデーションを作ってプロット
- 20. ある事例 Xi がクラスA に所属する確率は P(Xi) ある事例 Xi がクラスB に所属する確率は 1 - P(Xi) -P=0 2-class分類木の場合 P=1 sklearn.randomForestClassifier などは、弱学習器における確率値出力をしてくれるらし い sklearn.randomForestClassifier would rather pass on the probabilistic vote from terminals nodes and 弱学習器 = あるクラスへの所属確率をマップする関数、と考える クラスA に所属する確率
- 21. ある事例 Xi がクラス m に所属する確率 Pm(Xi) -Pm=0 multi-class分類木の場合 Pm=1 弱学習器 = あるクラスへの所属確率をマップする関数、と考える
- 22. 3-class分類の決定木 L L は、ノード5における特徴変数のContribution ⇒ 事例Xi がクラスm に所属する確率の増分
- 23. Irisデータで3クラス分類 ## Not run: example of plot_simplex3 library(randomForest) library(forestFloor) library(utils) data(iris) X = iris[,!names(iris) %in% "Species"] Y = iris[,"Species"] as.numeric(Y) rf.test42 = randomForest(X,Y, keep.forest=TRUE, replace=FALSE, keep.inbag=TRUE, samp=15, ntree=100) ff.test42 = forestFloor(rf.test42,X, calc_np=FALSE,binary_reg=FALSE) plot(ff.test42,plot_GOF=TRUE,cex=.7, colLists=list(c("#FF0000A5"), c("#00FF0050"), c("#0000FF35"))) show3d(ff.test42,1:2,3:4,plot_GOF=TRUE) 各変数の各クラスへの所属確率への貢献度
- 25. 所属確率の可視化（3-class限定） # つづき pars = plot_simplex3(ff.test42,Xi=c(1:3),restore_par=FALSE,zoom.fit=NULL, var.col=NULL,fig.cols=2,fig.rows=1,fig3d=FALSE,includeTotal=TRUE,auto.alpha=.4, set_pars=TRUE) pars = plot_simplex3(ff.test42,Xi=0,restore_par=FALSE,zoom.fit=NULL, var.col=alist(alpha=.3,cols=1:4),fig3d=FALSE,includeTotal=TRUE, auto.alpha=.8,set_pars=FALSE) for (I in ff.test42$imp_ind[1:4]) { #plotting partial OOB-CV separation(including interactions effects) #coloured by true class pars = plot_simplex3(ff.test42,Xi=I,restore_par=FALSE,zoom.fit=NULL, var.col=NULL,fig.cols=4,fig.rows=2,fig3d=TRUE,includeTotal=FALSE,label.col=1:3, auto.alpha=.3,set_pars = (I==ff.test42$imp_ind[1])) #coloured by varaible value pars = plot_simplex3(ff.test42,Xi=I,restore_par=FALSE,zoom.fit=TRUE, var.col=alist(order=FALSE,alpha=.8),fig3d=FALSE,includeTotal=(I==4), auto.alpha=.3,set_pars=FALSE)
- 26. 余談： グラデーション作成関数 fcol() library(forestFloor) X = data.frame(matrix(rnorm(1000), nrow=1000, ncol=4)) X[] = lapply(X,jitter,amount = 1.5) plot(X)
- 27. 余談： グラデーション作成関数 fcol() library(forestFloor) X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4)) X[] = lapply(X,jitter,amount = 1.5) plot(X,col=fcol(X,1)) １番目の変数の勾配でカラー生成
- 28. 余談： グラデーション作成関数 fcol() library(forestFloor) X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4)) X[] = lapply(X,jitter,amount = 1.5) plot(X,col=fcol(X,2)) ２番目の変数の勾配でカラー生成
- 29. 余談： グラデーション作成関数 fcol() library(forestFloor) X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4)) X[] = lapply(X,jitter,amount = 1.5) plot(X,col=fcol(X,1:2)) １＋２番目の変数の線形な勾配でカラー生成
- 30. 余談： グラデーション作成関数 fcol() library(forestFloor) X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4)) X[] = lapply(X,jitter,amount = 1.5) plot(X,col=fcol(X,1:3)) １～３番目の変数のPCA勾配でカラー生成
- 31. 余談： グラデーション作成関数 fcol() library(forestFloor) X = data.frame(matrix(rnorm(1000),nrow=1000,ncol=4)) X[] = lapply(X,jitter,amount = 1.5) plot(X,col=fcol(X,1:4)) １～４番目の変数のPCA勾配でカラー生成
- 32. まとめ • ランダムフォレストの感度分析に便利な指標を提供 • feature contributionによる『効果』の可視化 • グラデーション関数の提供 • 回帰だけでなく、2クラス分類・多クラス分類もｏｋ • （ Cross-Validationを使うための拡張方法も併せて提案 ） • もとのPDPは決定木ベースでなくてもいける
- 33. 参考文献 • randomForest • Breiman, L. (2001). Random forests. Machine learning, 45(1), 5–32. • Hastie, Tibshirani, Friedman (2008), Partial dependency plot. IN: The Elements of Statistical Learning.(2nd) pp367- • http://statweb.stanford.edu/~tibs/ElemStatLearn/ • forestFloor • CRAN • https://cran.r-project.org/web/packages/forestFloor/index.html • Official Site • http://forestfloor.dk/ • Welling et al. (2016). Forest Floor Visualizations of Random Forests. ArXiv e-prints, June 2016.“ • http://arxiv.org/abs/1605.09196 • Palczewska et al (2014). Interpreting random forest classification models using a feature contribution method. • http://arxiv.org/abs/1312.1121 • ICEbox • CRAN • https://cran.r-project.org/web/packages/ICEbox/index.html • Goldstein et al, (2015). Peeking Inside the Black Box: Visualizing Statistical Learning With Plots of Individual Conditional Expectation. Journal of Computational and Graphical Statistics, 24(1): 44-65 • http://arxiv.org/abs/1309.6392