クロスバリデーションの基本的な方法ホールドアウト法やり方と変数選択と推定の安定性のシミュレーションを行った。 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