RでK-fold Cross Validationを実装してみた

以下の記事を参考にK-fold クロスバリデーションを実装してみました。解く問題はkNN法のハイパーパラメータのkを決定する問題です。

Cross-Validation for Predictive Analytics Using R - MilanoR

作業概要

• irisデータセット（n = 150）を使用。
• irisデータセットを80%：20%で訓練データ（n = 120）とテストデータ（n = 30）に分ける。
• 訓練データを更に80%：20%で訓練データとバリデーションデータに分ける。（5-foldのクロスバリデーション）
• 訓練データ（n = 120）で最適なkを決定する際には各foldのaccuracyの平均を使用する。
• 最後に訓練データとテストデータを合わせて（irisデータ全体：n = 150）でもう一度クロスバリデーションを実施し、モデルのaccuracyの推定値を計算する

クロスバリデーションのやり方については、以下の記事を参考にしました。

Chapter 29 Cross validation | Introduction to Data Science

データの分割方法イメージは以下の通りです。

実装コード

今回作成したコードは以下の通りです。

library(tidyverse)
library(caret)

set.seed(2019)

### Choosing best kNN's k
# split iris data into train data and test data
data_size <- nrow(iris)
train_rate <- 0.8
train_i <- sample(1:data_size, size  = data_size * train_rate, replace = FALSE)
train_data <- iris[train_i, ]
test_data <- iris[-train_i, ]

# get index for spliting train data into each fold
fold_i <- sample(rep(1:n_fold, length.out = nrow(train_data)))

# parameters for cross validation
kmax <- 20
k_range <- 1:kmax
n_fold <- 5

# initialise matrix which holds validation result(accuracy)
accuracy <- matrix(NA, nrow = n_fold, ncol = kmax)

### cross validation 
for (i in 1:n_fold) {
  validate_i <- which(fold_i == i)
  train_temp <- train_data[-validate_i, ]
  validate_temp <- train_data[validate_i, ]
  # build model
  models <- apply(t(k_range), 2, function(k) knn3(Species ~ ., train_temp, k = k))
  preds <- mapply(function(model) predict(model, newdata = validate_temp, type = "class"), models)
  accuracy[i, ] <- colMeans(validate_temp$Species == preds)
}

### plot cross validation accuracy
av <- colMeans(accuracy)
df <- data.frame(cbind(k = 1:kmax, av))
ggplot(data=df, aes(x = k)) +
  geom_line(aes(y = av) , color = "red") +
  labs(x = "k-fold's k", y = "validation accuracy") +
  coord_cartesian(ylim = c(0, 1))
  
### choose best k
kbest <- which(av == max(av))
if (length(kbest) != 1) {
  # choose lower k
  kbest <- kbest[1]
}
kbest
max(av)

### estimate test accuracy
# get fold index for full dataset (iris)
fold_i <- sample(rep(1:n_fold, length.out = nrow(iris)))
# initialise accuracy vecotr
accuracy <- rep(NA, length.out = n_fold)

# cross validation of full dataset (iris)
for (i in 1:n_fold) {
  # get index of fold i
  test_i <- which(fold_i == i)
  # split iris data into train data and test data
  train_temp <- iris[-test_i, ]
  test_temp <- iris[test_i, ]
  # build model and predict
  model <- knn3(Species ~ ., train_temp, k = kbest)
  pred <- predict(model, newdata = test_temp, type = "class")
  # calculate accuracy 
  accuracy[i] <- mean(test_temp$Species == pred)
}

# calculate test accuracy (estimate)
test_accuracy <- mean(accuracy)
test_accuracy

実装のポイント

クロスバリデーションのメイン部分は以下です。

### cross validation 
for (i in 1:n_fold) {
  validate_i <- which(fold_i == i)
  train_temp <- train_data[-validate_i, ]
  validate_temp <- train_data[validate_i, ]
  # build model
  models <- apply(t(k_range), 2, function(k) knn3(Species ~ ., train_temp, k = k))
  preds <- mapply(function(model) predict(model, newdata = validate_temp, type = "class"), models)
  accuracy[i, ] <- colMeans(validate_temp$Species == preds)
}

ポイント以下です。

• 事前に設定したインデックス(fold_i)を使用して、バリデーションデータに重複がないように訓練データとバリデーションデータを作成する。

実行結果概要

クロスバリデーションを行い、accuracyの平均をとった結果が以下のグラフです。

### plot cross validation accuracy
av <- colMeans(accuracy)
df <- data.frame(cbind(k = 1:kmax, av))
ggplot(data=df, aes(x = k)) +
  geom_line(aes(y = av) , color = "red") +
  labs(x = "k-fold's k", y = "validation accuracy") +
  coord_cartesian(ylim = c(0, 1))

大きな違いは出ませんでしたが、この中で最もaccuracyが高いkを最適なkとします。accuracyが同じ場合はkが小さいほうを採用します。
なお、最適なkの決定方法としてよく使われる"one-standard error"ルールは簡略化のために今回使用しておりません。

"one-standard error"ルールについては以下を参照ください。
Cross-Validation for Predictive Analytics Using R - MilanoR

### choose best k
kbest <- which(av == max(av))
if (length(kbest) != 1) {
  # choose lower k
  kbest <- kbest[1]
}
kbest
max(av)

> kbest
[1] 6
> max(av)
[1] 0.9583333

計算の結果、最適なkは6となりました。

クロスバリデーションの実装としては以上ですが、追加で、モデルのaccuracyの推定値を計算しました。
訓練データにテストデータをあわせたデータ全体でもう一度クロスバリデーションを行い、その平均accuracyをモデルのaccuracyの推定値とします。

# cross validation of full dataset (iris)
for (i in 1:n_fold) {
  # get index of fold i
  test_i <- which(fold_i == i)
  # split iris data into train data and test data
  train_temp <- iris[-test_i, ]
  test_temp <- iris[test_i, ]
  # build model and predict
  model <- knn3(Species ~ ., train_temp, k = kbest)
  pred <- predict(model, newdata = test_temp, type = "class")
  # calculate accuracy 
  accuracy[i] <- mean(test_temp$Species == pred)
}

# calculate test accuracy (estimate)
test_accuracy <- mean(accuracy)
test_accuracy

> test_accuracy
[1] 0.96

なお、max(av)はモデルaccuracyの推定値としては使用できません。max(av)は最適なkを使用するためのバリデーションデータでの精度です。最終的なモデルaccuracyの推定値は、最初に取っておいたテストデータを含めた全量でクロスバリデーションを行った際のaccuracyを使用します。

以上