如何用R的lars()函数复现多数据集LASSO回归流程

原elasticnet实现代码 #

我原本使用elasticnet包的enet()函数对k个CSV格式数据集拟合LASSO回归，代码如下：

set.seed(150)
system.time(LASSO <- lapply(datasets, function(J) 
               elasticnet::enet(x = as.matrix(dplyr::select(J, 
                                         starts_with("X"))), 
               y = J$Y, lambda = 0, normalize = FALSE)))

提取系数的代码：

## 存储并输出LASSO选中的所有回归方程系数估计值
LASSO_Coeffs <- lapply(LASSO, 
                       function(i) predict(i, 
                                           x = as.matrix(dplyr::select(i, starts_with("X"))), 
                                           s = 0.1, mode = "fraction", 
                                           type = "coefficients")[["coefficients"]]) 

筛选正系数变量名称的代码：

IVs_Selected <- lapply(LASSO_Coeffs, function(i) names(i[i > 0])) 

使用lars包的等价实现 #

以下是用lars包的lars()函数复现上述流程的代码：

1. 拟合LASSO回归 #

set.seed(150)
# 加载lars包
library(lars)
system.time(LASSO_lars <- lapply(datasets, function(J) {
  # 提取自变量矩阵和因变量向量
  x <- as.matrix(dplyr::select(J, starts_with("X")))
  y <- J$Y
  # 拟合LASSO（type="lasso"），normalize=FALSE对应原代码设置
  lars(x = x, y = y, type = "lasso", normalize = FALSE)
}))

2. 提取指定lambda比例下的系数 #

lars的predict()函数中，s参数同样支持mode="fraction"（表示L1范数占最大L1范数的比例），对应原代码的s=0.1：

LASSO_Coeffs_lars <- lapply(LASSO_lars, function(model) {
  predict(model, s = 0.1, mode = "fraction", type = "coefficients")$coefficients
})

3. 筛选正系数的变量名称 #

这一步逻辑和原代码完全一致：

IVs_Selected_lars <- lapply(LASSO_Coeffs_lars, function(coeffs) {
  names(coeffs[coeffs > 0])
})

数据集加载参考代码（中文注释） #

我加载n个数据集到R并存储为datasets列表的完整代码如下：

# 生成数据集文件夹内所有文件的完整路径列表
folderpath <- "C:/Users/Spencer/Documents/EER Project/Data/0.5-5-1-1 to 0.5-6-10-500"
paths_list <- list.files(path = folderpath, full.names = T, recursive = T)

# 生成数据集名称列表（去掉文件扩展名）
DS_names_list <- basename(paths_list)
DS_names_list <- tools::file_path_sans_ext(DS_names_list)

# 批量读取数据集到R环境
system.time( datasets <- lapply(paths_list, fread) )

因为要一次性加载5000、10000或15000个数据集，且初始加载为字符类型，后续处理代码如下：

# 为所有数据集统一设置列名
datasets <- lapply(datasets, function(dataset_i) { 
  colnames(dataset_i) <- c("Y","X1","X2","X3","X4","X5","X6","X7","X8",
                           "X9","X10","X11","X12","X13","X14","X15",
                           "X16","X17","X18","X19","X20","X21","X22", 
                           "X23","X24","X25","X26","X27","X28","X29","X30")
  dataset_i })

# 提取结构工具变量（首行除Y外的列）
Structural_IVs <- lapply(datasets, function(j) {j[1, -1]})
Structural_Variables <- lapply(Structural_IVs, function(i) {names(i)[i == 1]})

# 移除前3行，转换所有列为数值型，再转为data.table格式
datasets <- lapply(datasets, function(i) {i[-1:-3, ]})
datasets <- lapply(datasets, \(X) { lapply(X, as.numeric) })
datasets <- lapply(datasets, function(i) { as.data.table(i) })

内容的提问来源于stack exchange，提问作者Marlen