head(x, 3)

head(y, 3)

# 必要なパッケージの読み込み
library(glmnet)

# 乱数の固定
set.seed(1234)

# Lasso回帰を使って「効きそうな説明変数」を自動選抜
cv_fit <- cv.glmnet(x, y, alpha = 1, standardize = FALSE)
lasso_fit <- glmnet(x, y, alpha = 1, standardize = FALSE, lambda = cv_fit$lambda.min)

# Lasso回帰の結果、係数が0でない列名を抽出
lasso_coefs <- as.matrix(coef(lasso_fit))
active_vars <- which(lasso_coefs[-1, 1] != 0)
print(paste("選ばれた説明変数:", names(active_vars)))

# 選択された説明変数で回帰分析
selected_x <- x[, active_vars, drop = FALSE]
model <- lm(y ~ selected_x)
summary(model)

set.seed(1234)
sim_p <- replicate(200, {
  xs <- matrix(rnorm(n * p), n, p)
  ys <- rnorm(n)
  cvs <- cv.glmnet(xs, ys, alpha = 1)
  fs <- glmnet(xs, ys, alpha = 1, lambda = cvs$lambda.min)
  act <- which(as.matrix(coef(fs))[-1, 1] != 0)
  if(length(act) > 0) {
    m_sim <- lm(ys ~ xs[, act, drop=FALSE])
    summary(m_sim)$coefficients[-1, 4]
  } else {
    NULL
  }
})
hist(unlist(sim_p), breaks = 20, col = "#e63946", 
     main = "変数選択後に回帰分析を行った場合のP値の分布",
     xlab = "P値", ylab = "頻度")

set.seed(1234)
sim_p <- replicate(200, {
  xs <- matrix(rnorm(n * p), n, p)
  ys <- rnorm(n)
  m_sim <- lm(ys ~ xs)
  summary(m_sim)$coefficients[-1, 4]
})
hist(unlist(sim_p), breaks = 20, col = "#e63946", 
     main = "変数選択をせずに回帰分析を行った場合のP値の分布",
     xlab = "P値", ylab = "頻度")

# 簡易シミュレーション比較
set.seed(1234)
n_sim <- 100
n_samples <- 60
p_vars <- 10
signal <- 0.4

power_res <- replicate(n_sim, {
  # データ生成
  X <- matrix(rnorm(n_samples * p_vars), n_samples, p_vars)
  # 第1変数だけ真の効果がある
  Y <- signal * X[, 1] + rnorm(n_samples)

  # 1. 選択的推論 (全データ)
  fit_all <- glmnet(X, Y, alpha = 1, lambda = 0.1)
  # 簡易的な評価（実際はもっと複雑）
  out_si <- tryCatch(fixedLassoInf(X, Y, beta = as.vector(coef(fit_all, s=0.1)[-1]), lambda = 0.1*n_samples), error=function(e) NULL)
  p_si <- if(!is.null(out_si) && out_si$vars[1] == 1) out_si$pv[1] else 1
  
  # 2. データ分割 (半分)
  idx <- 1:(n_samples/2)
  fit_split <- glmnet(X[idx,], Y[idx], alpha = 1, lambda = 0.1)
  selected <- which(as.vector(coef(fit_split, s=0.1)[-1]) != 0)
  p_split <- 1
  if(1 %in% selected) {
    m_test <- lm(Y[-idx] ~ X[-idx, 1])
    p_split <- summary(m_test)$coefficients[2, 4]
  }

  c(SI = p_si < 0.05, Split = p_split < 0.05)
})

power_df <- data.frame(
  Method = c("選択的推論", "データ分割"),
  Power = rowMeans(power_res, na.rm = TRUE)
)

ggplot(power_df, aes(x = Method, y = Power, fill = Method)) +
  geom_bar(stat = "identity") +
  theme_minimal(base_size = 16) +
  scale_fill_manual(values = c("#a8dadc", "#457b9d")) +
  labs(title = "真の効果を検出できる確率 (Power)", y = "検出力") +
  ylim(0, 1)