反復交差検証 Repeated Cross-Validationの安定性のシミュレーションを行った。 反復交差検証については以下のエントリを参照。
安定性のシミュレーションはLOOCVやnｰfolds CVの場合と同様に異なった乱数シードを使い、50回推定して、各変数が選ばれた回数と推定値の平均・標準偏差を計算した
データの読み込みと前処理
library(haven)
library(tidyr)
auto <- haven::read_dta("http://www.stata-press.com/data/r9/auto.dta")
auto <- auto %>% drop_na() # 全カラムに対してNAがない行を抽出
auto$foreign <- as.integer(auto$foreign)
並列計算の準備
doParallelパッケージを用いて1コアではなく複数コアで計算するようにして、計算時間を短縮させるようにした。
library(doParallel) # 並列処理の設定 num_cores <- detectCores() - 2 # システム全体のコア数から2つを引いた数を使用 cl <- makeCluster(num_cores) registerDoParallel(cl)
Lasso回帰の安定性を調べる
反復交差検証の反復回数はこちらで調べたように14回に設定した。
library(caret)
library(glmnet)
library(dplyr)
# 結果を格納するリスト
selected_variables <- list()
coefficients <- list()
# 50回のシミュレーション
for (i in 1:50) {
set.seed(i)
train_control <- trainControl(method = "repeatedcv",
number = 10, # 10-fold CV
repeats = 14, # 14回繰り返す
savePredictions = "final")
model <- train(price ~ mpg + rep78 + headroom + trunk + weight + length + turn + displacement + gear_ratio + foreign,
data = auto,
method = "glmnet", # Lasso回帰
trControl = train_control,
tuneLength = 10) # チューニングするλの数
# 最適なモデルの係数を取得
best_lambda <- model$bestTune$lambda
coef_best <- coef(model$finalModel, s = best_lambda)
# 選ばれた変数とその係数を保存
non_zero_coefs <- as.matrix(coef_best)
non_zero_coefs <- non_zero_coefs[non_zero_coefs[, 1] != 0, , drop = FALSE]
selected_variables[[i]] <- rownames(non_zero_coefs)
coefficients[[i]] <- non_zero_coefs
}
# 並行処理の停止
stopCluster(cl)
# 選ばれた変数の回数をカウント
variable_counts <- table(unlist(selected_variables))
# 選ばれた変数の推定値の平均と標準偏差を計算
variable_stats <- lapply(names(variable_counts), function(var) {
coefs <- unlist(lapply(coefficients, function(coef) if (var %in% rownames(coef)) coef[var, 1] else NA))
coefs <- coefs[!is.na(coefs)] # NAを除外
mean_coef <- mean(coefs)
sd_coef <- sd(coefs)
data.frame(variable = var, mean = mean_coef, sd = sd_coef)
})
variable_stats <- do.call(rbind, variable_stats)
変数が選ばれた回数
variable_counts_df <- as.data.frame(variable_counts)
colnames(variable_counts_df) <- c("Variable", "Count")
print(variable_counts_df)
Variable Count 1 (Intercept) 50 2 displacement 50 3 foreign 50 4 gear_ratio 32 5 headroom 50 6 length 12 7 mpg 16 8 rep78 50 9 trunk 12 10 turn 49 11 weight 50
推定値の平均と標準偏差
print(variable_stats)
variable mean sd 1 (Intercept) 349.73377 2330.9384052 2 displacement 13.57817 0.3616850 3 foreign 2839.84030 124.6804858 4 gear_ratio -203.81945 189.8846262 5 headroom -464.11233 54.2174071 6 length -12.02154 10.8404136 7 mpg -29.14755 8.0687846 8 rep78 159.76717 47.6848193 9 trunk 19.99706 10.3307869 10 turn -42.77568 22.0358426 11 weight 1.92157 0.2248898
