LASSO线性回归模型

LASSO线性回归模型

LASSO是1996年由Tibshirani提出的一种惩罚方法，可以同时进行变量选择和参数估计，适用于高维数据。
特点：稀疏性，不具有无偏性和一致性，不具有Oracle属性

1. 研究背景

例如研究基因对某个生物表征的影响，假定共有p个基因的n次观测值（p>>n)，因变量是连续型变量。我们想要筛选与因变量强相关的基因，直接建立线性回归模型会造成过拟合，可以通过对系数施加惩罚，如L1范数，即增加LASSO惩罚项，将大量冗余变量的系数收缩到0，以此达到变量筛选的目的。
具有LASSO惩罚项的线性回归模型要求解的目标函数为：

\[min_{\beta\in R^{p}}L(\beta)=min_{\beta\in R^{p}} \left\{\frac{1}{2n}\left\Vert Y-X\beta\right\Vert_{2}^{2}+\lambda\Vert \beta\Vert_{1}\right\} \ \ \ (1) \]

其中，第一项的L2范数是用误差平方和表示的损失函数，第二项是LASSO惩罚项，使得所有系数绝对值之和小于某个较小的正常数。

2. 通过坐标下降法求解

坐标下降法

基本思想：在每步迭代中沿一个坐标的方向进行搜索，通过循环使用不同的坐标方向来达到目标函数的局部极小值，属于非梯度优化方法。适用于线性回归模型求解，可广泛用于不可导的优化问题。

具体过程：
假设目标函数为\(f(x_1,x_2,⋯,x_p )\)

1. 选取\(x_2,x_3,⋯,x_p\)的初始值；
2. 在每轮迭代中：
   固定\(x_2,x_3,⋯,x_p\)，将\(x_1\)当作自变量，最小化\(f(x_1)\)，得到\(x_1^{*}\)；
   将\(x_1^{∗}\)代入目标函数，同时固定\(x_3,⋯,x_p\)，最小化\(f(x_2)\)，得到\(x_2^{∗}\)；
   将\(x_1^{∗},x_2^{∗}\)代入目标函数，同时固定\(x_4,⋯,x_p\)，最小化\(f(x_3)\)，得到\(x_3^{∗}\)；
   ……
   得到本轮迭代的一组值\(x_1^{∗},x_2^{∗},⋯,x_p^{∗}\)；
   若满足迭代终止条件则终止迭代，否则进行下一轮迭代。

LASSO线性回归求解过程

将目标函数(1)对\(\beta_j\)求导，令导数等于零，即：

\[\frac{\partial L(\beta)}{\partial \beta}=-\frac{1}{n}\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}+\frac{1}{n}\sum_{i=1}^n x_{ij}^2\beta_j+\lambda sign{\beta_j}=0, \]

其中\(\tilde{y}_i^{(-j)}=\sum_{k=1,k\neq j}^p x_{ik}\beta_k\)。

求解上式时分情况讨论，可以得到\(\beta\)的估计值\(\hat{\beta}\)：

当\(\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}>0\)，且\(\lambda<\left|\frac{1}{n}\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}\right|\)时，有：

\[\hat{\beta}=\frac{1}{n}\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}-\lambda \]

当\(\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}<0\)，且\(\lambda<\left|\frac{1}{n}\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}\right|\)时，有：

\[\hat{\beta}=\frac{1}{n}\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}+\lambda \]

当\(\lambda\ge\left|\frac{1}{n}\sum_{i=1}^n (y_i-\tilde{y}_i^{(-j)})x_{ij}\right|\)时，有：

\[\hat{\beta}=0 \]

该解析解也可以用软阈值函数表示

3. Matlab代码

通过glmnet包求解

clear;
clc;
load data;   #加载数据，包括数据矩阵x，因变量y
x=zscore(x);
fold = 10;    #10折交叉验证
m = 10;       #参数lambda的个数
options = glmnetSet;
options.alpha = 1;
options.nlambda = m;
fit = glmnet(x,y,'gaussian',options);
CVerr = cvglmnet(x,y,fold,[],'response','gaussian',glmnetSet,0);
beta_path = fit.beta;    #得到对应于所有lambda的系数矩阵
intercept_path = fit.a0;   #截距项
lambda_path = fit.lambda;   #lambda的取值
opt_index = find(CVerr.glmnetOptions.lambda == CVerr.lambda_min); #最优lambda所在位置
beta = beta_path(:,opt_index);
selectedgene = find(beta~=0);   #模型筛选出来的变量所在位置