PDPとICEを利用したEDA（探索的データ解析）｜R

2019.12.14

この記事は約6分で読めます。

この記事では、PDPとICEを利用したEDA（探索的データ解析）について解説します。

EDAのアプローチとして、予測力の高いモデルを使って特徴量重要度を出すことが有効な場合が多々あります。

以下は、Bostonデータに対するランダムフォレストモデルboston.rfをもとに特徴量重要度を出力するRコード例です。

#ランダムフォレストモデルの作成
#データの準備
xy <- Boston
colnames(xy)[ncol(xy)] <- "y"

#データを学習用とテスト用に分割
n <- nrow(xy)
set.seed(111, sample.kind = "Rounding")
test.id <- sample(n, round(n / 4))
test <- xy[test.id, ]
train <- xy[-test.id, ]

#ランダムフォレストモデルの学習
library(randomForest)
set.seed(111)
boston.rf <- randomForest(y ~ .,
                          data = train,
                          importance = TRUE)

#テストデータに対する予測
pred <- predict(boston.rf, newdata = test)
plot(test$y, pred, main = boston.rf$call)
curve(identity, add = TRUE)
rms <- function(act, pred) {
  sqrt(mean((act - pred) ^ 2))
}
cat(" RMSE =", rms(test$y, pred))

#特徴量重要度の出力
boston.imp <-
  sort(importance(boston.rf, type = 2, scale = FALSE)[, 1],
       decreasing = TRUE)
barplot(boston.imp, names.arg = rownames(boston.imp))

特徴量重要度は、予測に対する寄与の大きさを把握するのに参考になります。

しかし、特徴量がとる値によって予測値がどのくらい変化するのか、また特徴量と予測値の相関関係（正の相関または負の相関なのか）といった寄与方向を上図から読み取ることはできません。

そこで、考えられたのがPDP（Partial Dependence Plot）とICE（Indivisual Conditional Expectation）です。

PDPとICEは特徴量の値の変化による予測値の変動幅、寄与方向の参考になります。

歴史的にはPDPが先に登場したそうですが、説明を簡単にするためにICEから説明します。

ICEでは、着目している特徴量以外は観測値のまま固定し、着目している特徴量の値を可能な範囲で動かしながら、用意した予測モデルに基づいて予測値を計算します。

上記処理を観測の数だけ繰り返し、着目している特徴量を横軸に、予測値を縦軸にとってプロットしたものがICEです。

ランダムフォレストモデルboston.rfで特徴量重要度が最も高かった特徴量lstatのICEを描きます。

ICEは、pdpパッケージの関数partialとplotPartialを利用することで簡単に描けます。

#ICE
library(pdp)
ice <- partial(object = boston.rf,
               pred.var = c("lstat"),
               ice = T)
plotPartial(ice)

ICEからはlstatが低いほどyhatは高くなっていることが読み取れます。

図の中央の赤線はICEの平均をとったものです。

この赤線がPDPです。

PDPのみを描くRコードは次のとおりです。

#PDP
pdp <- partial(object = boston.rf,
               pred.var = c("lstat"))
plot(pdp, type = "l", col = "red")

すべての特徴量のPDPをまとめて出力したい場合は、次のRコードを実行すればOKです。

#PDPをまとめて出力
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 5))
for (i in 1:(ncol(xy) - 1)) {
  xname <- colnames(xy)[i]
  pdp <- partial(object = boston.rf,
                 pred.var = c(xname))
  plot(pdp, type = "l", col = "red")
  abline(h = mean(train$y))
}
par(oldpar)