Rのcaretパッケージで機械学習ー回帰問題

Introduction

Rの機械学習って色々なパッケージに分かれていてややこしい．そんなとき200種類以上の機械学習手法を簡単な記述で扱える便利なパッケージcaretを見つけたので基本的な使い方についてまとめる．

データの準備

今回はボストンデータを使用．
chasとradはカテゴリー変数であるため，factor型に変換しておく

library(MASS)
data = Boston

caretパッケージのインストールと呼び出し

#install.packages("caret")
library(caret)

並列化処理

モデルの学習を高速化するため，複数のコアで同時に処理をしてもらうための命令文．caretを使用するときはこちらも呪文として唱えておくとよい．

#install.packages("doParallel")
library(doParallel)

# 全コア数の半分を使用
cores = detectCores() / 2
cl <- makePSOCKcluster(cores)
registerDoParallel(cl)

データの分割

訓練データとテストデータの分割にはcaretパッケージ内のcreateDataPartitionを使用．pで分割割合を指定．今回は8 : 2に分割．

set.seed(11)

# データ分割用行番号作成
partition = createDataPartition(data[,1], p = 0.8, list = FALSE)

data_train = data[partition,]
data_test = data[-partition,]

X_train = data_train[,colnames(data) != 'medv']
y_train = data_train['medv']

X_test = data_test[,colnames(data) != 'medv']
y_test = data_test['medv']

データの前処理

標準化

データの標準化（平均0，標準偏差１）にはpreProcessを使用．

# 標準化モデルの構築
scale_model = preProcess(X_train, method = c('center', 'scale'))

# 学習データの標準化
X_train_scale = predict(scale_model, newdata = X_train)
# テストデータの正規化
X_test_scale = predict(scale_model, newdata = X_test)

ダミー変数化

カテゴリー変数のダミー変数化（One-Hot Encoding）にはdummyVarsを使用．

# ダミー変数化モデルの構築
dummy_model = dummyVars(~., data = X_train, fullRank = TRUE)

# 学習データのダミー変数化
X_train_scaledummy = data.frame(predict(dummy_model, newdata = X_train_scale))
# テストデータのダミー変数化
X_test_scaledummy = data.frame(predict(dummy_model, newdata = X_test_scale))

データの一部

モデルの構築

caretパッケージ内で使用できるアルゴリズム一覧はこのサイトで確認できる．
caretパッケージ内アルゴリズム一覧
またmodelLookupでアルゴリズムの詳細が確認できる．
今回はランダムフォレストでモデルを構築してみる．

modelLookup('rf')

Out:

学習

学習はtrainで簡単に行える．
最低限以下の3つの引数が必要．

• formula：モデルの構造を記述したオブジェクト
• data：目的変数と説明変数が含まれたデータ
• method：アルゴリズムの種類

上記に加え今回はパラメータのチューニングをクロスバリデーションで行うため，以下の2つの引数も追加．

• tuneLength：各パラメータのチューニング回数
• trControl：パラメータチューニングの方法

引数をあらかじめ設定

formula = medv ~ .  # 目的変数medv, 説明変数すべて
method = 'rf'  # ランダムフォレスト
tL = ncol(X_train_scaledummy) - 1  # 今回は説明変数はそれほど多くないため，説明変数の数-1とする 
tC = trainControl(method = 'cv', number = 5)  # 5分割のクロスバリデーション

最適なハイパーパラメータを持つモデルが生成される．

set.seed(11)
model_rf <- train(formula, data = cbind(y_train, X_train_scaledummy),
                  method = method, tuneLength = tL, trControl = tC)

学習結果

学習の結果，9個の特徴量を含むモデルのRMSEが最も小さく，最良のモデルといえる．

model_rf

Out:

Random Forest 

406 samples
 20 predictor

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 325, 324, 325, 326, 324 
Resampling results across tuning parameters:

  mtry  RMSE      Rsquared   MAE     
   2    4.288801  0.8206868  2.898890
   3    3.771919  0.8551542  2.530535
   4    3.536236  0.8685549  2.393931
   5    3.407725  0.8756650  2.290953
   6    3.359574  0.8769825  2.256613
   7    3.299778  0.8802924  2.224651
   8    3.311583  0.8787008  2.231915
   9    3.277041  0.8808178  2.220307
  10    3.279955  0.8797193  2.233822
  11    3.340146  0.8748618  2.245376
  12    3.339566  0.8748642  2.267656
  13    3.333614  0.8750673  2.259978
  14    3.312243  0.8765812  2.244778
  15    3.326672  0.8755117  2.249137
  16    3.380582  0.8708817  2.269478
  17    3.368999  0.8716857  2.284721
  18    3.392529  0.8699610  2.282623
  19    3.392954  0.8693971  2.272084
  20    3.397931  0.8691550  2.279761

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was mtry = 9.

特徴量の重要度はvarImpで見ることができる．
今回はrmからageまでを含むモデルが選択されたことになる．

plot(varImp(model_rf))

Out:

モデルが格納されたオブジェクトをplotに入れれば，パラメータ毎の精度がグラフで確認できる．

plot(model_rf)

Out:

テストデータで検証

最後にテストデータの精度を確認する.
predictで予測値を算出し，postResampleで精度を確認．

pred = predict(model_rf, X_test_scaledummy)
postResample(y_test, pred)

Out:
     RMSE  Rsquared       MAE 
3.0827929 0.8790809 2.1232139

Memo

パラメータチューニングの評価指標のデフォルトはRMSEになっているが，別の評価指標を使用したい場合はtrainの引数にmetricを追加する．
以下のコードは評価指標にMAEを使用した例．

set.seed(11)
model_rf = train(formula, data = cbind(y_train, X_train_scaledummy),
                  method = method, tuneLength = tL, trControl = tC,
                  metric = 'MAE')

Conclusion

今回使用したcaretパッケージ内の関数．

• createDataPartition：訓練データとテストデータの分割
• preProcess：データの標準化
• dummyVars：データのダミー変数化
• modelLookup：機械学習アルゴリズムの詳細を確認
• trainControl：パラメータチューニングの方法
• train：モデルの学習
• varImp：特徴量重要度の表示
• postResample：精度の確認

今回はランダムフォレストを回帰問題に適用してみた．
次回は分類問題への適用や，手法の比較をまとめてみようと思う．

Code

# データ準備
library(MASS)
data = Boston

data[,'chas'] = as.factor(data[,'chas'])
data[,'rad'] = as.factor(data[,'rad'])

# パッケージインストール
#install.packages("caret")
library(caret)

# 並列化処理
#install.packages("doParallel")
library(doParallel)

cores = detectCores() / 2
cl <- makePSOCKcluster(cores)
registerDoParallel(cl)

# データの分割
set.seed(11)
partition = createDataPartition(data[,1], p = 0.8, list = FALSE)

data_train = data[partition,]
data_test = data[-partition,]

X_train = data_train[,colnames(data) != 'medv']
y_train = data_train['medv']

X_test = data_test[,colnames(data) != 'medv']
y_test = data_test['medv']

# データの前処理
## 標準化
scale_model = preProcess(X_train, method = c('center', 'scale'))

X_train_scale = predict(scale_model, newdata = X_train)
X_test_scale = predict(scale_model, newdata = X_test)

## ダミー変数化
dummy_model = dummyVars(~., data = X_train, fullRank = TRUE)

X_train_scaledummy = data.frame(predict(dummy_model, newdata = X_train_scale))
X_test_scaledummy = data.frame(predict(dummy_model, newdata = X_test_scale))

# モデルの構築
## 学習
formula = medv ~ .  # 目的変数medv, 説明変数すべて
method = 'rf'  # ランダムフォレスト
tL = ncol(X_train_scaledummy) - 1  # 今回は説明変数はそれほど多くないため，説明変数の数-1とする 
tC = trainControl(method = 'cv', number = 5)  # 5分割のクロスバリデーション

set.seed(11)
model_rf = train(formula, data = cbind(y_train, X_train_scaledummy),
                  method = method, tuneLength = tL, trControl = tC)

## 学習結果
model_rf

plot(varImp(model_rf))

plot(model_rf)

# テストデータの精度を確認
pred = predict(model_rf, X_test_scaledummy)
postResample(y_test, pred)