クロスバリデーションの基本的な方法ホールドアウト法やり方と変数選択と推定の安定性のシミュレーションを行った。 n-fold クロスバリデーションがよくつかわれているため、ホールドアウト法はあまり使われないので実践向きではない。
今回は訓練データを70%、テストデータを30%と一番よく選ばれる比率で行った。
# 必要なライブラリの読み込み
library(haven)
library(dplyr)
library(glmnet)
library(caret)
# データの読み込みと前処理
auto <- haven::read_dta("http://www.stata-press.com/data/r9/auto.dta")
auto <- auto %>% drop_na() # 全カラムに対してNAがない行を抽出
auto$foreign <- as.factor(auto$foreign)
data <- auto[, c("price", "mpg", "rep78", "headroom", "trunk", "weight", "length", "turn", "displacement", "gear_ratio", "foreign")]
# 訓練データとテストデータに分割（70%:30%）
set.seed(123) # 再現性のためにシードを設定
train_index <- createDataPartition(auto$price, p = 0.7, list = FALSE)
train_data <- data[train_index, ]
test_data <- data[-train_index, ]
# モデルの訓練
x_train <- model.matrix(price ~ ., data = train_data)[, -1]
y_train <- train_data$price
lasso_model <- cv.glmnet(x_train, y_train, alpha = 1)
# 最適なlambdaを取得
best_lambda <- lasso_model$lambda.min
# モデルの係数を取得
coef_lasso <- coef(lasso_model, s = best_lambda)
# テストデータでの予測
x_test <- model.matrix(price ~ ., data = test_data)[, -1]
y_test <- test_data$price
# 予測を行う
predictions <- predict(lasso_model, s = best_lambda, newx = x_test)
# モデルの性能評価
mse <- mean((predictions - y_test)^2)
rmse <- sqrt(mse)
mae <- mean(abs(predictions - y_test))
r2 <- cor(predictions, y_test)^2
# 結果の表示
cat("最適なlambda: ", best_lambda, "\n")
cat("MSE: ", mse, "\n")
cat("RMSE: ", rmse, "\n")
cat("MAE: ", mae, "\n")
cat("R-squared: ", r2, "\n")
# 係数を表示
print("Lasso回帰の係数:")
print(coef_lasso)
係数
s1 (Intercept) -3222.784180 mpg . rep78 337.196361 headroom -347.170167 trunk . weight 1.820159 length . turn . displacement 14.260552 gear_ratio . foreign1 2583.081529
統計量類
最適なlambda: 140.1621 MSE: 5518769 RMSE: 2349.206 MAE: 1813.332 R-squared: 0.337947
変数選択と推定量の安定性のシミュレーション
以下のエントリと同じ条件で反復計算をした。
異なる乱数シードを設定し、50回反復計算をして、1) 独立変数側で選ばれた変数の回数、2)選ばれた時の推定値の平均とその標準偏差を計算した。
library(haven)
library(dplyr)
library(glmnet)
library(caret)
library(doParallel)
# データの読み込みと前処理
auto <- haven::read_dta("http://www.stata-press.com/data/r9/auto.dta")
auto <- auto %>% drop_na() # 全カラムに対してNAがない行を抽出
auto$foreign <- as.factor(auto$foreign)
data <- auto[, c("price", "mpg", "rep78", "headroom", "trunk", "weight", "length", "turn", "displacement", "gear_ratio", "foreign")]
# 50回のシミュレーション結果を格納するリスト
selected_variables <- list()
coefficients <- list()
# 並列処理の設定
num_cores <- detectCores() - 2 # システム全体のコア数から2つを引いた数を使用
cl <- makeCluster(num_cores)
registerDoParallel(cl)
# 50回のシミュレーション
for (i in 1:50) {
set.seed(i) # 異なる乱数シードを設定
# 訓練データとテストデータに分割（70%:30%）
train_index <- createDataPartition(data$price, p = 0.7, list = FALSE)
train_data <- data[train_index, ]
test_data <- data[-train_index, ]
# モデルの訓練
x_train <- model.matrix(price ~ ., data = train_data)[, -1]
y_train <- train_data$price
lasso_model <- cv.glmnet(x_train, y_train, alpha = 1)
# 最適なlambdaを取得
best_lambda <- lasso_model$lambda.min
# モデルの係数を取得
coef_lasso <- as.matrix(coef(lasso_model, s = best_lambda))
# 選ばれた変数とその係数を保存
non_zero_coefs <- coef_lasso[coef_lasso != 0, , drop = FALSE]
selected_variables[[i]] <- rownames(non_zero_coefs)
coefficients[[i]] <- non_zero_coefs # 切片も含む
}
# 並行処理の停止
stopCluster(cl)
# Jaccard係数の計算
calculate_jaccard <- function(set1, set2) {
intersect_len <- length(intersect(set1, set2))
union_len <- length(union(set1, set2))
return(intersect_len / union_len)
}
# すべての組み合わせに対するJaccard係数を計算
jaccard_matrix <- matrix(0, nrow = 50, ncol = 50)
for (i in 1:50) {
for (j in 1:50) {
jaccard_matrix[i, j] <- calculate_jaccard(selected_variables[[i]], selected_variables[[j]])
}
}
# 平均Jaccard係数を計算
mean_jaccard <- mean(jaccard_matrix[upper.tri(jaccard_matrix, diag = FALSE)])
print(paste("平均Jaccard係数:", mean_jaccard))
# 選ばれた変数の回数をカウント
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)) {
return(coef[rownames(coef) == var, 1])
} else {
return(NA)
}
}))
coefs <- coefs[!is.na(coefs)]
mean_coef <- mean(coefs, na.rm = TRUE)
sd_coef <- sd(coefs, na.rm = TRUE)
data.frame(variable = var, mean = mean_coef, sd = sd_coef)
})
variable_stats <- do.call(rbind, variable_stats)
variable_counts <- as.data.frame(variable_counts)
colnames(variable_counts) <- c("Variable", "Count")
print(variable_counts)
variable_stats <- as.data.frame(variable_stats)
print(variable_stats)
選ばれた変数の回数
> print(variable_counts)
Variable Count
1 (Intercept) 50
2 displacement 49
3 foreign1 50
4 gear_ratio 22
5 headroom 49
6 length 13
7 mpg 21
8 rep78 31
9 trunk 13
10 turn 27
11 weight 50
Jaccard係数
今回はJaccard係数も計算した。
Jaccard係数（Jaccard index）とは、2つの集合の類似度を測るための統計的指標の一つである。具体的には、2つの集合がどれだけ共通しているかを表す指標であり、次のように計算される。
"平均Jaccard係数: 0.695805872132403"
推定値の平均と標準偏差
> print(variable_stats)
variable mean sd
1 (Intercept) 3360.632722 5965.710162
2 displacement 12.196286 5.103783
3 foreign1 2874.097212 629.855815
4 gear_ratio -572.481917 647.989326
5 headroom -419.618083 234.511503
6 length -81.376900 33.351025
7 mpg -45.288926 36.645942
8 rep78 203.456938 259.791928
9 trunk 55.313944 63.119901
10 turn -110.594330 111.530918
11 weight 2.535542 1.807896
