LASSO回归变量选择¶

假设你在装修房子，有100个装饰品可选：

Ridge（L2）：每个装饰品都摆上，但都缩小了
  → 所有变量保留，系数都变小

LASSO（L1）：只选20个摆上，其他80个扔掉
  → 自动选变量，不重要的系数=0

ElasticNet：选30个摆上，但有些缩小
  → 介于两者之间，适合相关变量多的情况

# 安装glmnet（2026最新版）
install.packages("glmnet")  # 安装
library(glmnet)              # 加载

# 准备数据（以基因表达预测疾病为例）
# x: 预测变量矩阵（行=样本，列=基因）
# y: 响应变量（0/1的疾病状态）

# 读取数据
expr_data <- read.csv("expression_matrix.csv", row.names = 1)  # 表达矩阵

# 转为矩阵格式（glmnet要求）
x <- as.matrix(expr_data[, -1])  # 预测变量（去掉第一列如果是标签）
y <- expr_data[, 1]               # 响应变量

# 或者用model.matrix（处理分类变量更方便）
# x <- model.matrix(y ~ ., data = expr_data)[, -1]

# ===== 步骤1：交叉验证找最优lambda =====
set.seed(42)  # 设置随机种子（保证可重复）

cv_fit <- cv.glmnet(
  x = x,                    # 预测变量矩阵
  y = y,                    # 响应变量
  family = "binomial",      # 二分类（逻辑回归）；连续变量用"gaussian"
  alpha = 1,                # alpha=1 是LASSO
  nfolds = 10,              # 10折交叉验证
  type.measure = "auc"      # 用AUC作为评价指标；也可用"deviance"
)

# 画交叉验证图
plot(cv_fit)  # X轴=log(lambda)，Y轴=AUC，上方数字=非零系数个数
# 左虚线 = lambda.min（最优性能）
# 右虚线 = lambda.1se（简约模型）

# 查看两个关键lambda值
cat("lambda.min:", cv_fit$lambda.min, "\n")  # 最优lambda
cat("lambda.1se:", cv_fit$lambda.1se, "\n")  # 简约lambda

# ===== 步骤2：提取系数（选出的变量） =====
# 用lambda.1se（更简约，推荐发文章）
coef_1se <- coef(cv_fit, s = "lambda.1se")  # 提取系数
# 用lambda.min（更多变量，性能更好）
coef_min <- coef(cv_fit, s = "lambda.min")

# 找出非零系数的变量（即被选中的基因）
selected_genes <- rownames(coef_1se)[which(coef_1se != 0)]  # 非零系数
selected_genes <- selected_genes[selected_genes != "(Intercept)"]  # 去掉截距项
cat("选出的基因数:", length(selected_genes), "\n")
print(selected_genes)  # 打印选出的基因

# ===== 步骤3：用选出的变量重新建模 =====
# 方法1：直接用glmnet的预测
pred_prob <- predict(cv_fit, newx = x, s = "lambda.1se", type = "response")
# type="response"返回概率值（0-1之间）

# 方法2：用选出的变量建logistic回归
selected_data <- data.frame(y = y, x[, selected_genes])
final_model <- glm(y ~ ., data = selected_data, family = binomial)
summary(final_model)  # 查看系数和p值

# 在生存分析中用LASSO选变量（非常常见）
library(survival)  # 加载survival包

# 准备Surv对象
y_surv <- Surv(time = clinical$time, event = clinical$status)  # 生存时间和事件

# 交叉验证找最优lambda
cv_cox <- cv.glmnet(
  x = x,                    # 基因表达矩阵
  y = y_surv,               # Surv对象
  family = "cox",           # Cox回归
  alpha = 1,                # LASSO
  nfolds = 10,              # 10折CV
  type.measure = "C"        # C-index作为评价指标
)

# 画图和提取系数
plot(cv_cox)
coef_cox <- coef(cv_cox, s = "lambda.min")  # Cox常用lambda.min
selected_cox <- rownames(coef_cox)[which(coef_cox != 0)]
cat("选出的基因数:", length(selected_cox), "\n")

# 计算风险评分
risk_score <- predict(cv_cox, newx = x, s = "lambda.min", type = "link")
# risk_score越高，风险越大

# 当基因间高度共线性时，用ElasticNet
# 先找最优alpha
alpha_values <- seq(0.1, 0.9, by = 0.1)  # 测试多个alpha值
cv_results <- list()

for (a in alpha_values) {
  set.seed(42)
  cv_results[[as.character(a)]] <- cv.glmnet(
    x = x, y = y, family = "binomial",
    alpha = a, nfolds = 10, type.measure = "auc"
  )
}

# 找AUC最高的alpha
best_auc <- sapply(cv_results, function(fit) max(fit$cvm))
best_alpha <- alpha_values[which.max(best_auc)]
cat("最优alpha:", best_alpha, "\n")

# 用最优alpha重新拟合
cv_enet <- cv.glmnet(x, y, family = "binomial", alpha = best_alpha, nfolds = 10)

# 画系数随lambda变化的路径
fit <- glmnet(x, y, family = "binomial", alpha = 1)  # 拟合模型
plot(fit, xvar = "lambda", label = TRUE)  # 系数路径图
# X轴=log(lambda)，Y轴=系数值
# 每条线是一个基因，lambda增大时系数逐渐归零
abline(v = log(cv_fit$lambda.1se), lty = 2)  # 标记选择的lambda

# 画选出变量的系数条形图
library(ggplot2)

coef_df <- data.frame(
  Gene = selected_genes,
  Coefficient = as.numeric(coef_1se[selected_genes, ])
)
coef_df <- coef_df[order(abs(coef_df$Coefficient), decreasing = TRUE), ]

ggplot(coef_df, aes(x = reorder(Gene, abs(Coefficient)), y = Coefficient)) +
  geom_bar(stat = "identity", fill = ifelse(coef_df$Coefficient > 0, "red", "blue")) +
  coord_flip() +                    # 横向条形图
  labs(x = "Gene", y = "LASSO Coefficient",
       title = "Selected Genes by LASSO") +
  theme_minimal()

lambda.min: 交叉验证误差最小的λ
  优点：模型性能最好
  缺点：可能选太多变量（过拟合风险）
  适用：探索性分析、预测为主

lambda.1se: 误差在最小值±1SE范围内的最大λ
  优点：模型更简洁，泛化更好
  缺点：可能丢失一些有用变量
  适用：发文章、需要简洁模型

生信论文推荐：
  - 发现生物标志物 → lambda.1se（选变量更严格）
  - 构建预测模型 → lambda.min（保留更多信息）
  - 实际用时画两个λ的结果对比，讨论差异

# LASSO快速流程
as.matrix(x) → cv.glmnet(x, y, alpha=1) → plot() → coef(s="lambda.1se")

# 三种回归
alpha = 1:   LASSO（选变量）
alpha = 0:   Ridge（不选变量）
alpha = 0.5: ElasticNet（折中）

# family参数
"gaussian":  连续变量（线性回归）
"binomial":  二分类（逻辑回归）
"cox":       生存分析（Cox回归）
"multinomial": 多分类
"poisson":   计数数据

# 两个关键lambda
lambda.min: 性能最好
lambda.1se: 模型最简洁（推荐）

# 常用评价指标（type.measure）
"deviance": 偏差（默认）
"auc":      AUC（二分类推荐）
"C":        C-index（Cox推荐）
"mse":      均方误差（连续变量）