こちらで利用したデータとコードを用いる。

データの準備

library(haven)
library(tidyr)
# データの読み込みと前処理
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)

必ず独立変数に入れる変数を指定する

ここではmpgを必ず独立変数に入れる変数に指定する。 mpgのペナルティをゼロに設定すると、必ず選ばれるようになる。

# ペナルティファクターの設定（）
penalty_factors <- rep(1, ncol(X))
penalty_factors[which(colnames(X) == "mpg")] <- 0

クロスバリデーションを行い最適なラムダを選択

library(glmnet)
set.seed(123)  # 再現性のためのシード設定
cv_model <- cv.glmnet(X, price, alpha = 1, penalty.factor = penalty_factors, nfolds = 10)

Lasso回帰

# 最適なλの取得
best_lambda <- cv_model$lambda.min

# 最適なλに対するモデルの係数を表示
final_model <- glmnet(X, price, alpha = 1, penalty.factor = penalty_factors)
coef(final_model, s = best_lambda)

確認

必ず独立変数に選ばれているか、50回反復して確認してみよう。先のエントリで50回の反復計算したところ、n-folds CVでmpgは8回選ばれていた。

# 変数名の取得
variables <- colnames(X)

# 乱数シードの設定
set.seed(1234)
seeds <- sample(1:10000, 50)

# 結果を保存するリストを初期化
coefficients_list <- list()

# 50回の推定を実行
for (i in 1:50) {
  set.seed(seeds[i])
  
  # ペナルティファクターの設定（mpgのペナルティをゼロに設定）
  penalty_factors <- rep(1, ncol(X))
  penalty_factors[which(colnames(X) == "mpg")] <- 0
  
  cv_fit <- cv.glmnet(X, price, alpha = 1, penalty.factor = penalty_factors, nfolds = 10)
  best_lambda <- cv_fit$lambda.min
  fit <- glmnet(X, price, alpha = 1, penalty.factor = penalty_factors, lambda = best_lambda)
  
  coefficients_list[[i]] <- as.numeric(coef(fit))[-1]  # インターセプトを除く
}

# 結果1: 各変数が選ばれた回数の集計
selected_counts <- sapply(1:length(variables), function(var_index) {
  sum(sapply(coefficients_list, function(coefs) coefs[var_index] != 0))
})

# 結果2: 各変数の推定値の平均と標準偏差の計算
coef_stats <- t(sapply(1:length(variables), function(var_index) {
  coefs <- sapply(coefficients_list, function(coefs) coefs[var_index])
  coefs <- coefs[coefs != 0]  # ゼロを除外
  c(mean = mean(coefs), sd = sd(coefs))
}))

変数が選ばれた回数

selected_counts_df <- data.frame(Variable = variables, SelectedCount = selected_counts)
selected_counts_df

       Variable SelectedCount
1           mpg            50
2         rep78            50
3      headroom            50
4         trunk            14
5        weight            50
6        length            23
7          turn            50
8  displacement            50
9    gear_ratio            23
10      foreign            50

mpgは50回すべてで選ばれているので、うまく言っていることが確認できた。

各変数の推定値の平均と標準偏差

coef_stats_df <- data.frame(Variable = variables, Mean = coef_stats[, "mean"], SD = coef_stats[, "sd"])
coef_stats_df

       Variable        Mean          SD
1           mpg  -32.666970   6.4224423
2         rep78  170.231132   9.3563991
3      headroom -489.302250  56.8276089
4         trunk   34.969388  17.5331591
5        weight    2.400621   0.9040117
6        length  -36.981308  20.4223442
7          turn  -64.729485  30.9242399
8  displacement   13.427778   0.6517092
9    gear_ratio -156.180244  84.7626351
10      foreign 2929.969871 187.8901826