モンテカルロ交差検証(Monte Carlo cross-validation, MCCV)、もしくは反復ランダムサブサンプリング検証(Repeated random sub-sampling validation)と呼ばれる方法は、訓練データと検証データをランダムに分割し、ホールドアウト法と同じく訓練データに基づいて検証データモデル評価を行うという手順を複数回(例えば100回)行う。ランダム分割の数が無限大に近づくにつれて、p-個抜き交差検証の結果に近づくとされている。反復交差検証とよく似た手順で行われるが、反復交差検証はk-fold交差検証を繰り返すのに比べ、モンテカルロ交差検証はホールドアウト法を繰り返すところに違いがある。
モンテカルロ交差検証の実施
# 必要なライブラリの読み込み
library(haven)
library(tidyr)
library(caret)
library(glmnet)
library(dplyr)
# データの準備
set.seed(123)
auto_data <- haven::read_dta("http://www.stata-press.com/data/r9/auto.dta")
auto_data <- auto_data %>% drop_na() # 全カラムに対してNAがない行を抽出
auto_data$foreign <- as.integer(auto_data$foreign)
price_data <- auto_data$price
X_data <- auto_data[, c("mpg", "rep78", "headroom", "trunk", "weight", "length",
"turn", "displacement", "gear_ratio", "foreign")]
X_data <- as.matrix(X_data)
# MCCVの設定
train_control <- trainControl(method = "LGOCV",
p = 0.8, # トレーニングセットの割合
number = 100, # 繰り返し回数
savePredictions = "final")
# Lasso回帰モデルのトレーニング
lasso_model <- train(x = X_data, y = price_data,
method = "glmnet",
trControl = train_control,
tuneLength = 10) # チューニングするλの数
# 結果の表示
print(lasso_model)
# 最良のλを取得
best_lambda <- lasso_model$bestTune$lambda
# 最終モデルの係数を取得
final_coefs <- coef(lasso_model$finalModel, s = best_lambda)
print(final_coefs)
結果。
RMSE was used to select the optimal model using the smallest value.
The final values used for the model were alpha = 0.1 and lambda = 111.2479.
11 x 1 sparse Matrix of class "dgCMatrix"
s1
(Intercept) 5491.703265
mpg -30.009245
rep78 219.079762
headroom -592.963110
trunk 31.700221
weight 2.529418
length -24.302025
turn -93.892082
displacement 14.249337
gear_ratio -401.315437
foreign 2989.835976
変数選択と推定量の安定性のシミュレーション
異なる乱数シードを設定し、50回反復計算をして、1) 独立変数側で選ばれた変数の回数、2)選ばれた時の推定値の平均とその標準偏差を計算した。
# 必要なライブラリの読み込み
library(haven)
library(tidyr)
library(caret)
library(glmnet)
library(dplyr)
library(doParallel)
# データの準備
set.seed(123)
auto_data <- haven::read_dta("http://www.stata-press.com/data/r9/auto.dta")
auto_data <- auto_data %>% drop_na() # 全カラムに対してNAがない行を抽出
auto_data$foreign <- as.integer(auto_data$foreign)
price_data <- auto_data$price
X_data <- auto_data[, c("mpg", "rep78", "headroom", "trunk", "weight", "length",
"turn", "displacement", "gear_ratio", "foreign")]
X_data <- as.matrix(X_data)
# 結果を保存するリストを初期化
coef_list <- list()
selected_vars_list <- list()
n_repeats <- 50
# 並列処理
num_cores <- detectCores() - 2
cl <- makeCluster(num_cores)
registerDoParallel(cl)
# 50回の反復計算
for (i in 1:n_repeats) {
set.seed(i + 1000) # 異なる乱数シードを設定
# MCCVの設定
train_control <- trainControl(method = "LGOCV",
p = 0.8, # トレーニングセットの割合
number = 1, # 繰り返し回数
savePredictions = "final",
allowParallel = TRUE)
# Lasso回帰モデルのトレーニング
lasso_model <- train(x = X_data, y = price_data,
method = "glmnet",
trControl = train_control,
tuneLength = 10) # チューニングするλの数
# 最良のλを取得
best_lambda <- lasso_model$bestTune$lambda
# 最終モデルの係数を取得
final_coefs <- coef(lasso_model$finalModel, s = best_lambda)
coef_list[[i]] <- as.vector(final_coefs)
names(coef_list[[i]]) <- rownames(final_coefs) # 係数に名前をつける
# 選ばれた変数を保存
final_coefs_matrix <- as.matrix(coef(lasso_model$finalModel, s = best_lambda))
selected_vars <- rownames(final_coefs_matrix)[final_coefs_matrix != 0]
selected_vars_list[[i]] <- selected_vars
}
# 並列処理の停止
stopCluster(cl)
# 変数名を取得
variables <- rownames(final_coefs_matrix)
# 結果1: 各変数が選ばれた回数の集計
selected_counts <- sapply(variables, function(var) {
sum(sapply(coef_list, function(coefs) !is.na(coefs[var]) && coefs[var] != 0))
})
# 結果2: 各変数の推定値の平均と標準偏差の計算
coef_stats <- sapply(variables, function(var) {
coefs <- unlist(lapply(coef_list, function(coefs) if (!is.na(coefs[var]) && coefs[var] != 0) coefs[var] else NA))
coefs <- coefs[!is.na(coefs)] # NAを除外
c(mean = mean(coefs), sd = sd(coefs))
})
# Jaccard係数の計算
jaccard_index <- function(set1, set2) {
intersect_len <- length(intersect(set1, set2))
union_len <- length(union(set1, set2))
if (union_len == 0) return(NA)
return(intersect_len / union_len)
}
jaccard_indices <- combn(1:length(selected_vars_list), 2, function(index) {
set1 <- selected_vars_list[[index[1]]]
set2 <- selected_vars_list[[index[2]]]
jaccard_index(set1, set2)
})
mean_jaccard_index <- mean(jaccard_indices, na.rm = TRUE)
# データフレームにまとめる
selected_counts_df <- data.frame(Variable = variables, SelectedCount = selected_counts)
coef_stats_df <- data.frame(Mean = coef_stats["mean",], SD = coef_stats["sd",])
# 結果の表示
print(selected_counts_df)
print(coef_stats_df)
cat("平均Jaccard係数: ", mean_jaccard_index, "\n")
選ばれた変数の回数
Variable SelectedCount (Intercept) (Intercept) 50 mpg mpg 36 rep78 rep78 50 headroom headroom 50 trunk trunk 32 weight weight 50 length length 28 turn turn 39 displacement displacement 50 gear_ratio gear_ratio 40 foreign foreign 50
Jaccard係数
平均Jaccard係数: 0.8055363
推定値の平均と標準偏差
Mean SD (Intercept) 4080.293882 4276.907768 mpg -24.196628 13.124276 rep78 172.274918 57.859475 headroom -472.151445 161.348569 trunk 42.094461 28.146215 weight 2.718052 1.326874 length -51.858017 27.296527 turn -78.406699 40.667062 displacement 12.104590 1.320900 gear_ratio -292.793462 102.845914 foreign 2818.097140 535.273449
