井出草平の研究ノート

Lasso ホールドアウト法と変数選択と推定の安定性のシミュレーション[R]

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

変数選択と推定量の安定性のシミュレーション

以下のエントリと同じ条件で反復計算をした。

ides.hatenablog.com

異なる乱数シードを設定し、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"

 \text{Mean Jaccard} = \frac{1}{n} \sum_{i=1}^n J(A_i, B_i)

推定値の平均と標準偏差

> 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