tuneRF関数の挙動の検証
RのrandomForestパッケージのtuneRF関数が遅いというお話があったので，調べてみた．
tuneRF関数は，ランダムフォレストを用いて予測モデルを構築する際に使用する特徴量の個数を貪欲的な方法により求める．
tuneRF関数だけでなく，Rで機械学習のアルゴリズムをチューニングする方法については，id:TJOさんの記事「Rで機械学習するならチューニングもグリッドサーチ関数orオプションでお手軽に」が詳しいので，是非参照してほしい．
検証に使用するデータと環境
UCI Machine Learning RepositoryのBank Marketingデータセットを使用する．このデータセットは，ポルトガルの銀行で電話でダイレクトマーケティングを実施した際に収集したデータを用いて，予測モデルのアルゴリズムに投入できるように特徴量を構築したもの．
まずは，データを取得して解凍する．
$ wget http://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank.zip $ unzip bank.zip
同じディレクトリに"bank-full.csv", "bank-names.txt", "bank.csv"という3つのファイルが解凍される．"bank-full.csv"がすべてのレコードを含んでおり，ここではこのデータを使用する．
データをRに読み込んで，行頭やサイズ，各列のデータ型を確認する．
> bank <- read.csv("bank-full.csv", sep=";") > head(bank, 3) age job marital education default balance housing loan contact day 1 58 management married tertiary no 2143 yes no unknown 5 2 44 technician single secondary no 29 yes no unknown 5 3 33 entrepreneur married secondary no 2 yes yes unknown 5 month duration campaign pdays previous poutcome y 1 may 261 1 -1 0 unknown no 2 may 151 1 -1 0 unknown no 3 may 76 1 -1 0 unknown no > dim(bank) [1] 45211 17 > sapply(bank, class) age job marital education default balance housing loan "integer" "factor" "factor" "factor" "factor" "integer" "factor" "factor" contact day month duration campaign pdays previous poutcome "factor" "integer" "factor" "integer" "integer" "integer" "integer" "factor" y "factor"
行数が45,211，列数が17であることを確認できる．最後の列(y)は，ダイレクトマーケティングに顧客が反応したかどうかを表すカテゴリ変数であり，この項目を目的変数として使用する(Yes: 反応した，No: 反応していない)．残りの16個の列を説明変数と使用する．上記の結果を見ると，16個のうち9項目がカテゴリ変数であることを確認できる．
なお，以下で検証に使用する環境は次の通り．
まずは単純にtuneRF関数を実行してみる
まずは特別な設定を行わずにtuneRF関数を実行してみよう．tuneRF関数の第1引数には説明変数，第2引数には目的変数を指定する．また，doBest引数をTRUEに指定すると，評価が最も良いモデルを返すようになる．
> set.seed(123) > system.time(tuneRF(bank %>% select(-y), bank$y, doBest=TRUE)) mtry = 4 OOB error = 9.46% Searching left ... mtry = 2 OOB error = 9.78% -0.03295162 0.05 Searching right ... mtry = 8 OOB error = 9.48% -0.001402197 0.05 ユーザ システム 経過 77.131 1.725 81.135
上記を見ると，おおよそ1分20秒程度実行に要したことを確認できる．この結果，特徴量の個数が
- 4個のときに，Out-of-Bag誤差(OOB error)は9.46%
- 2個のときに，Out-of-Bag誤差は9.78%
- 8個のときに，Out-of-Bag誤差は9.48%
となり，特徴量の個数が4個のときにOut-of-Bag誤差が最少となり，この個数に設定するのが良さそうであることがわかる*1．
構築する決定木の個数を増やしてみる
ntreeTry引数はデフォルトでは50となっており，50個の決定木を構築することがわかる．500個の決定木を構築するように指定してみよう．
> set.seed(123) > system.time(tuneRF(bank %>% select(-y), bank$y, ntreeTry=500, doBest=TRUE)) mtry = 4 OOB error = 9.1% Searching left ... mtry = 2 OOB error = 9.53% -0.04689018 0.05 Searching right ... mtry = 8 OOB error = 9.18% -0.008503401 0.05 ユーザ システム 経過 237.656 2.943 244.725
構築する決定木の個数が10倍になったので，実行時間も約10倍になるのかと思うところであるが，3倍程度に抑えられていることを確認できる．
tuneRF関数の引数と仕様
tuneRF関数にはいくつかの引数を指定できる．
> args(tuneRF) function (x, y, mtryStart = if (is.factor(y)) floor(sqrt(ncol(x))) else floor(ncol(x)/3), ntreeTry = 50, stepFactor = 2, improve = 0.05, trace = TRUE, plot = TRUE, doBest = FALSE, ...) NULL
| 引数 | 説明 | データ型 |
| x | 説明変数 | data.frame |
| y | 目的変数 | factor or numeric |
| mtryStart | 選択する特徴量の個数の初期値 | integer |
| ntreeTry | 構築する決定木の個数 | integer |
| stepFactor | 特徴量の個数を減少または増加させるファクター | integer |
| improve | モデルの評価指標の改善率 | numeric |
| trace | チューニングの過程をトレースして画面に表示するかどうかを指定するフラグ | logical |
| doBest | 評価が最も良いモデルを返すかどうかを指定 | logical |
また，tuneRF関数自体もそれほど長くないので，中身を読んでみる．
tuneRF <- function(x, y, mtryStart=if(is.factor(y)) floor(sqrt(ncol(x))) else floor(ncol(x)/3), ntreeTry=50, stepFactor=2, improve=0.05, trace=TRUE, plot=TRUE, doBest=FALSE, ...) { if (improve < 0) stop ("improve must be non-negative.") # yが因子かどうかを判定 classRF <- is.factor(y) errorOld <- if (classRF) { # yが因子の場合はクラス分類を実行し，誤分類率により評価 randomForest(x, y, mtry=mtryStart, ntree=ntreeTry, keep.forest=FALSE, ...)$err.rate[ntreeTry,1] } else { # yが因子でない場合は回帰を実行し，平均二乗平均誤差により評価 randomForest(x, y, mtry=mtryStart, ntree=ntreeTry, keep.forest=FALSE, ...)$mse[ntreeTry] } if (errorOld < 0) stop("Initial setting gave 0 error and no room for improvement.") # traceがTRUEに指定された場合は，チューニング過程を画面に出力 if (trace) { cat("mtry =", mtryStart, " OOB error =", if (classRF) paste(100*round(errorOld, 4), "%", sep="") else errorOld, "\n") } oobError <- list() oobError[[1]] <- errorOld names(oobError)[1] <- mtryStart # 以下，特徴量の個数を増減させて予測モデルを構築し結果を評価(left: 個数を削減，right: 個数を増加) for (direction in c("left", "right")) { if (trace) cat("Searching", direction, "...\n") Improve <- 1.1*improve mtryBest <- mtryStart mtryCur <- mtryStart while (Improve >= improve) { mtryOld <- mtryCur mtryCur <- if (direction == "left") { # 個数の削減は，現在の個数をstepFactor引数で割ったものを新しい個数とする max(1, ceiling(mtryCur / stepFactor)) } else { # 個数の増加は，現在の個数にstepFactor引数を掛けたものを新しい個数とする min(ncol(x), floor(mtryCur * stepFactor)) } if (mtryCur == mtryOld) break errorCur <- if (classRF) { # yが因子の場合はクラス分類を実行し，誤分類率により評価 randomForest(x, y, mtry=mtryCur, ntree=ntreeTry, keep.forest=FALSE, ...)$err.rate[ntreeTry,"OOB"] } else { # yが因子でない場合は回帰を実行し，平均二乗平均誤差により評価 randomForest(x, y, mtry=mtryCur, ntree=ntreeTry, keep.forest=FALSE, ...)$mse[ntreeTry] } if (trace) { cat("mtry =",mtryCur, "\tOOB error =", if (classRF) paste(100*round(errorCur, 4), "%", sep="") else errorCur, "\n") } oobError[[as.character(mtryCur)]] <- errorCur # 予測精度(誤差)の改善率の算出 Improve <- 1 - errorCur/errorOld cat(Improve, improve, "\n") if (Improve > improve) { errorOld <- errorCur mtryBest <- mtryCur } } } mtry <- sort(as.numeric(names(oobError))) res <- unlist(oobError[as.character(mtry)]) res <- cbind(mtry=mtry, OOBError=res) # plotがTRUEに指定された場合は，横軸を特徴量の個数，縦軸を評価値とする折れ線グラフをプロット if (plot) { plot(res, xlab=expression(m[try]), ylab="OOB Error", type="o", log="x", xaxt="n") axis(1, at=res[,"mtry"]) } # doBestがTRUEに指定された場合は，評価が最も良いモデルを返り値とする if (doBest) res <- randomForest(x, y, mtry=res[which.min(res[,2]), 1], ...) res }
以上のtuneRF関数の実装を読むと，特徴量の個数と評価結果を保持するノードを木構造で表現していくことがわかる．
そのため，この処理の高速化の方針として，
- 木の左側の計算と右側の計算を並列化する
- 個々のrandomForest関数で決定木を構築する処理を並列化する
の2つが考えられる．これ以外の方法は，randomForest関数，およびその内部のC++やFortranのコードを改変することになる．
なお，前者は，木が深い時に特に有効であり，後者は木の深さに関係なく，ある程度の個数の決定木を構築するならば有効な手段だと思われる．
*1:ランダムフォレストはOut-of-Bag誤差が求まるのでクロスバリデーション不要とも言われるが，実際はクロスバリデーションなども組み合わせて，もう少し精緻に検証することが望ましいと思われる．
