正則化項付き線形回帰は真の偏回帰係数を推定しているのか？
最近、正則化項付き線形回帰についてちょっと調べてます。
それで以下の記事が気になりました。
どちらも人工データを用いて、真の偏回帰係数を正則化項付き線形回帰で推定できるか？というシミュレーションをされています。
これは非常に興味深いので自分でもやってみようと思います。
先の記事はどちらもPythonを使われてましたが、私はR言語でやってみます。
試すのは以下の5つの手法です。
- 線形回帰
- Ridge回帰
- LASSO回帰
- 適応的LASSO回帰
- Elastic net回帰
確認したいのは真の偏回帰係数に対する推定された偏回帰係数の分布です、
準備
まず下準備として、必要なパッケージの呼び出しと、必要な関数の定義をします。
パッケージの呼び出し。
if (!require("pacman")) {install.packages("pacman")}
pacman::p_load(conflicted, tidyverse, glmnetUtils)λをAICcで決める関数。RidgeとLASSO用。
AICc_glmnet <- function(glmnet_model) { tLL <- glmnet_model$nulldev - deviance(glmnet_model) k <- glmnet_model$df n <- glmnet_model$nobs df <- tibble( lambda = glmnet_model$lambda, AICc = -tLL + 2 * k + 2 * k * (k + 1) / (n - k - 1) ) best_AICc_lambda <- df|> slice_min(AICc) |> pull(lambda) return( list( best_AICc_lambda = best_AICc_lambda, AICc = df$AICc ) ) }
αとλをAICcで決める関数。Elastic net用。
get_ic_best_alpha_lambda <- function(cva_model) { n <- cva_model$alpha |> length() ics <- map( 1:n, \(list_no) cva_model |> pluck("modlist", list_no, "glmnet.fit") |> AICc_glmnet() ) df_ic <- map( 1:n, \(list_no) data.frame( best_alpah = cva_model |> pluck("alpha", list_no), best_lambda = cva_model |> pluck("modlist", list_no, "lambda"), AICc = ics |> pluck(list_no, "AICc") ) ) |> list_rbind() df_ic_best <- list( AICc_best_alpha = df_ic |> slice_min(AICc) |> pull(best_alpah), AICc_best_lambda = df_ic |> slice_min(AICc) |> pull(best_lambda) ) return(df_ic_best) } || そして、シミュレーションを回す関数。 指定した説明変数の数とサンプルサイズのダミーデータを生成し、各手法に当てはめてそれぞれの偏回帰係数を出力します。 >|r| # sim_f(N, K, seed) # N: サンプルサイズ # K: 変数の数 # seed: 再現のための乱数種 sim_f <- function(N, K, seed) { set.seed(seed) X <- MASS::mvrnorm(N, mu = rep(0, K), Sigma = diag(K)) Y <- X %*% beta + rnorm(N) df <- X |> as.data.frame() |> bind_cols(data.frame(Y = Y)) # lm coef_lm <- lm(Y ~ ., data = df) |> broom::tidy() |> slice(-1) |> select(term, estimate) # Ridge ridge_best_lambda <- glmnet(Y ~ ., data = df, alpha = 0) |> AICc_glmnet() coef_ridge <- glmnet(Y ~ ., data = df, alpha = 0, lambda = ridge_best_lambda$best_AICc_lambda) |> broom::tidy() |> slice(-1) |> select(term, estimate) # LASSO lasso_best_lambda <- glmnet(Y ~ ., data = df, alpha = 1) |> AICc_glmnet() coef_lasso <- glmnet(Y ~ ., data = df, alpha = 1, lambda = lasso_best_lambda$best_AICc_lambda) |> broom::tidy() |> slice(-1) |> select(term, estimate) # 適応的LASSO coef_alasso <- glmnet(Y ~ ., data = df, alpha = 1, lambda = 5 / nrow(df)^(3 / 4), penalty.factor = 1 / abs(coef_lm$estimate) ) |> broom::tidy() |> slice(-1) |> select(term, estimate) # EN EN_best_alpha_lambda <- cva.glmnet(Y ~ ., data = df) |> get_ic_best_alpha_lambda() coef_EN <- glmnet(Y ~ ., data = df, alpha = EN_best_alpha_lambda$AICc_best_alpha, lambda = EN_best_alpha_lambda$AICc_best_lambda) |> broom::tidy() |> slice(-1) |> select(term, estimate) # res res <- coef_lm |> left_join(coef_ridge, by = join_by(term)) |> left_join(coef_lasso, by = join_by(term)) |> left_join(coef_alasso, by = join_by(term)) |> left_join(coef_EN, by = join_by(term)) |> rename( lm = estimate.x, ridge = estimate.y, lasso = estimate.x.x, alasso = estimate.y.y, EN = estimate ) return(res) }
これで準備完了です。
素直なデータセット
まずは素直なデータセットで試してみます。
説明変数5個、サンプルサイズ100です。
それぞれの真の偏回帰係数は-1, -0.5, 0, 0.5, 1です。
線形回帰は真の偏回帰係数を中心に分布していて、バラツキも小さくいい感じです。
正則化項付き回帰はいずれも過小推定の傾向があり、真の偏回帰係数をうまく推定できているとは言えそうにありません。
K <- 5 N <- 100 I <- 300 # 試行回数 beta <- seq(-1, 1, length.out = K) res <- map( 1:I, \(seed) sim_f(N, K, seed), .progress = TRUE ) |> list_rbind() |> pivot_longer(!term) |> mutate(value = replace_na(value, 0)) res |> ggplot(aes(x = term, y = value, fill = name)) + geom_violin() + geom_boxplot(fill="white", width = 0.2) + geom_hline(yintercept = beta, linetype = 2, color = "darkgray") + facet_wrap(vars(name)) + theme(legend.position = 'none') + labs(x = "説明変数", y = "偏回帰係数", title = "素直なデータ", subtitle = "300回試行")
サンプルサイズが少ないケース
今度は説明変数の数が多く、サンプルサイズが少ないケースを考えます。
ただ、サンプルサイズが説明変数の数より小さくなると線形回帰の解が出ないので、サンプルサイズを説明変数の数＋１にしています。
説明変数101個、サンプルサイズ102個です。
真の偏回帰係数がそれぞれ-1, -0.5, 0, 0.5, 1である説明変数について推定値を確認します。
線形回帰は真の値を中心に分布していますがとにかくバラツキが大きいので、個々の分析においては外れている可能性が大きそうです。
Ridgeは推定値が安定していますが、ほとんどゼロに近く過小推定の傾向にあります。
LASSO, エラスティックネットも安定しているがやや過小評価の傾向が見られます。
適応的LASSOが比較的ましですがやはり若干過小評価の傾向にはあります。
K <- 101 N <- 102 I <- 300 beta <- seq(-1, 1, length.out = K) res <- map( 1:I, \(seed) sim_f(N, K, seed), .progress = TRUE ) |> list_rbind() |> pivot_longer(!term) |> mutate(value = replace_na(value, 0)) res |> dplyr::filter(term %in% c("V1", "V26", "V51", "V76", "V101")) |> mutate(term = factor(term, levels=c("V1", "V26", "V51", "V76", "V101"))) |> ggplot(aes(x = term, y = value, fill = name)) + geom_violin() + geom_boxplot(fill="white", width = 0.2) + geom_hline(yintercept = c(-1, -0.5, 0, 0.5, 1.0), linetype = 2, color = "darkgray") + facet_wrap(vars(name)) + theme(legend.position = 'none') + coord_cartesian(ylim = c(-2,2))+ labs(x = "説明変数", y = "偏回帰係数", title = "サンプルサイズが少ないデータ", subtitle = "300回試行")
といったところで、正則化項付き線形回帰は予測モデルとしては使えても、それぞれの偏回帰係数を影響力の大きさとして解釈するのはやはり難しそうですね。
参考文献はこちら。
