井出草平の研究ノート

トップ > Lasso, Ridge, EN > ロジスティックLasso回帰の安定性選択(Stability Selection)[R]

ロジスティックLasso回帰の安定性選択(Stability Selection)[R]

Lasso, Ridge, EN R 計量

安定性選択(Stability Selection)について。

ides.hatenablog.com

データは説明はこちら。

ides.hatenablog.com

データの作成と前処理

library(dplyr)
url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.cleveland.data"
column_names <- c("age", "sex", "cp", "trestbps", "chol", "fbs", "restecg", "thalach", 
                  "exang", "oldpeak", "slope", "ca", "thal", "num")
heart_data <- read.csv(url, header = FALSE, col.names = column_names, na.strings = "?")

heart_data <- heart_data %>%
  mutate(class = ifelse(num > 0, 1, 0)) %>%  # num > 0 を心臓病ありとする
  select(-num) %>%  # 使用しない変数を除去
  na.omit()  # 欠損値を除去

# デザインマトリックスと応答ベクトルの作成
x <- as.matrix(heart_data %>% select(-class))
y <- heart_data$class

stabsパッケージを利用

args.fitfun = list(family = "binomial")と書くことでglmnetパッケージにわたる引数を指定できる。

library(glmnet)
library(stabs)
library(dplyr)

# 再現性のためのシード設定
set.seed(1234)  

# 安定性選択の実行
stab_glmnet_model <- stabsel(
  x = x,
  y = y,
  fitfun = glmnet.lasso,
  cutoff = 0.75,
  PFER = 1,
  B = 100,
  sampling.type = "SS",
  assumption = "r-concave",
  args.fitfun = list(family = "binomial")
)

# 結果の表示
print(stab_glmnet_model)
plot(stab_glmnet_model, main = "Stability Selection with Lasso (glmnet)")

結果。

 Stability Selection with r-concavity assumption

Selected variables:
  ca thal 
  12   13 

Selection probabilities:
     age trestbps     chol      fbs  restecg    slope      sex       cp  oldpeak    exang  thalach       ca     thal 
   0.000    0.000    0.000    0.000    0.005    0.030    0.035    0.605    0.605    0.630    0.720    0.945    1.000 

---
Cutoff: 0.75; q: 5; PFER (*):  0.976 
   (*) or expected number of low selection probability variables
PFER (specified upper bound):  1 
PFER corresponds to signif. level 0.0751 (without multiplicity adjustment)

hdiパッケージを利用

こちらは、args.model.selector = list(family = "binomial")と書くことでglmnetパッケージにわたる引数を指定できる。

library(glmnet)
library(hdi)
library(dplyr)


# 安定性選択の実行
fit.stab <- stability(
  x, 
  y,
  EV = 1,  # PFER
  threshold = 0.75,  # cutoff
  B = 100,  # サブサンプル数
  args.model.selector = list(family = "binomial")
)

# 結果の表示
print(fit.stab)
print(fit.stab$freq)

結果。

Selected predictors:
--------------------
thal 
  13 
--------------------
Expected number of false positives controlled at level 1 
     age      sex       cp trestbps     chol      fbs  restecg  thalach    exang  oldpeak    slope       ca     thal 
    0.00     0.00     0.23     0.00     0.00     0.00     0.00     0.30     0.24     0.23     0.00     0.73     0.94

hdiパッケージではsampling.type = "SS"やassumption = "r-concave"といった推定方法の指定ができないので、値が少し異なる。

プロット

hdiパッケージにはプロット機能がないので、ggplotでプロットを描く。

library(ggplot2)

# 選択頻度のデータフレームを作成
freq_df <- data.frame(
  Predictor = names(fit.stab$freq),
  Frequency = fit.stab$freq
)

# 頻度で降順にソート
freq_df <- freq_df %>%
  arrange(desc(Frequency))

# プロットの作成
ggplot(freq_df, aes(x = Frequency, y = reorder(Predictor, Frequency))) +
  geom_point(color = "red", size = 3) +
  geom_segment(aes(x = 0, xend = Frequency, y = Predictor, yend = Predictor), color = "grey") +
  geom_vline(xintercept = 0.75, linetype = "dashed", color = "red") +  # カットオフの縦線を追加
  labs(title = "Stability Selection with Lasso (glmnet)",
       x = expression(hat(pi)),
       y = "") +
  theme_minimal()