LASSO线性回归模型

LASSO线性回归模型

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

1. 研究背景

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

$$mi{n}_{\beta \in R^p}L(\beta )=mi{n}_{\beta \in R^p}\{\frac{1}{2n}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_1\}(1)$$

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

2. 通过坐标下降法求解

坐标下降法

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

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

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

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\ne j}^p{x}_{ik}{\beta }_k$。

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

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

$$\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 <|\frac{1}{n}\sum _{i=1}^n(y_i-{\tilde{y}}_i^{(-j)})x_{ij}|$时，有：

$$\hat{\beta }=\frac{1}{n}\sum _{i=1}^n(y_i-{\tilde{y}}_i^{(-j)})x_{ij}+\lambda$$

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

$$\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);   #模型筛选出来的变量所在位置