通常のクロスバリデーションでは乱数によって結果がかなり違うため、そこで今回は一つ抜き交差検証（Leave-One-Out Cross-Validation, LOOCV）の安定性のシミュレーションを行うことにした。
LOOCVの計算の方法はこちら。
また、結果にどのくらいの安定性があるか、計算のたびに異なる乱数を与え、LOOCVのLassoの計算を100回反復した。1) 独立変数側で選ばれた変数の回数、2)選ばれた時の推定値の平均とその標準偏差を計算した。
安定性のシミュレーション
# 必要なライブラリの読み込み
library(haven)
library(tidyr)
library(nestedcv)
library(doParallel)
# データの用意
auto <- haven::read_dta("http://www.stata-press.com/data/r9/auto.dta")
auto <- auto %>% drop_na() # 全カラムに対してNAがない行を抽出
n <- nrow(auto)
price <- auto$price
X <- auto[, c("mpg", "rep78", "headroom", "trunk", "weight", "length",
"turn", "displacement", "gear_ratio", "foreign")]
X$foreign <- as.integer(X$foreign)
X <- as.matrix(X)
# 乱数シードの設定
set.seed(123)
repeat_times <- 100 # 実際の実行時には100回に設定
seeds <- sample(1:10000, repeat_times) # 指定された回数の乱数シードを生成
# 結果を保存するリストを初期化
coefficients_list <- list()
# 並列処理
num_cores <- detectCores() - 2
# 推定を実行
for (i in 1:repeat_times) {
set.seed(seeds[i])
fit2 <- nestcv.glmnet(
y = price,
x = X,
family = "gaussian",
alphaSet = 1,
outer_method = "LOOCV", # LOOCVを指定
n_outer_folds = nrow(X), # データポイント数に一致させる
n_inner_folds = 10, # 内側のCVのfold数
cv.cores = num_cores # 並列処理のためのコア数
)
coefficients_list[[i]] <- fit2$final_coef
}
# 変数のリスト（全ての可能な変数名）
fixed_row_names <- c("(Intercept)", "foreign", "headroom", "rep78",
"displacement", "weight", "turn", "gear_ratio",
"length", "trunk", "mpg")
# 各リスト要素を処理し、欠けている行をNAで埋める関数
process_coefficients <- function(coef_list, fixed_row_names) {
df <- as.data.frame(coef_list)
df_fixed <- data.frame(matrix(NA, nrow = length(fixed_row_names), ncol = ncol(df)))
rownames(df_fixed) <- fixed_row_names
colnames(df_fixed) <- colnames(df)
df_fixed[rownames(df), ] <- df
return(df_fixed)
}
# coefficients_listを処理
processed_list <- lapply(coefficients_list, process_coefficients, fixed_row_names = fixed_row_names)
# リストの各要素を結合して1つのデータフレームにする
final_df <- do.call(cbind, lapply(processed_list, function(x) x["coef"]))
# 変数が選ばれた回数の計算
selected_counts <- apply(final_df, 1, function(row) sum(!is.na(row)))
# 選ばれた時の平均値とその標準偏差の計算
coef_stats <- t(apply(final_df, 1, function(row) {
valid_values <- row[!is.na(row)]
mean_value <- mean(valid_values)
sd_value <- sd(valid_values)
c(mean = mean_value, sd = sd_value)
}))
# 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係数の計算
jaccard_indices <- combn(1:repeat_times, 2, function(index) {
set1 <- rownames(processed_list[[index[1]]])[!is.na(processed_list[[index[1]]]$coef)]
set2 <- rownames(processed_list[[index[2]]])[!is.na(processed_list[[index[2]]]$coef)]
jaccard_index(set1, set2)
})
mean_jaccard_index <- mean(jaccard_indices, na.rm = TRUE)
# 結果をデータフレームにまとめる
results_df <- data.frame(
Variable = rownames(final_df),
SelectedCount = selected_counts,
Mean = coef_stats[, "mean"],
SD = coef_stats[, "sd"]
)
# Jaccard係数の結果を表示
print(paste("Mean Jaccard Index:", mean_jaccard_index))
# 選ばれた回数と統計量の結果を表示
print(results_df)
Variable SelectedCount Mean SD (Intercept) (Intercept) 100 850.690282 3144.3525824 foreign foreign 100 2981.745823 181.9827829 headroom headroom 100 -462.979344 72.1114951 rep78 rep78 100 132.656012 20.0136089 displacement displacement 100 13.215741 0.6015858 weight weight 100 2.379150 0.7957266 turn turn 80 -61.469561 29.2796136 gear_ratio gear_ratio 62 -140.013670 110.0054683 length length 32 -39.067285 17.9979634 trunk trunk 26 31.002440 17.1044275 mpg mpg 24 -9.062924 5.2083239
Mean Jaccard Index: 0.798967598058507
