1 为什么要做变量选择？
变量选择又叫模型选择，因为选出不同的变量自然就建立出不同的模型。近30年来，变量选择理论繁荣发展并逐渐渗透到各个领域中，成为大数据时代背后的中坚理论之一。
随着技术手段的发展和丰富，我们获取数据变得越来越容易，对于一个问题可能搜集到成百上千的变量。有人说，变量越多对一个问题的刻画就越细致，这是好事啊。那么，我们为什么在建立模型时要煞费苦心思考变量选择问题呢？
一方面，搜集到的变量中可能存在与目标完全无关的变量（冗余变量）。如果将冗余变量纳入模型之中肯定会得到一个错误的模型，用其进行预测自然就会做出误判，这明显与我们的目标相违背。因此，模型能够剔除不相关变量显得很重要。
另一方面，尽管有些变量已知和目标相关，但实际的影响微乎其微（类似线性模型中变量系数的绝对值近乎为零），把这些变量包含在模型中无疑增加了理解模型的难度。比如，建立一个含有50+变量与一个仅有5个变量的公司员工能力评估模型，如果两者在评估员工能力时效果相当，显然后者更易于理解员工能力的影响因素且更易于指导实际工作。
所以，研究变量选择是很有现实意义的。
2 早期的思考
原有的方法若可用，就省的重新造轮子了。然而，经典的最小二乘 (OLS) 没有变量筛选的能力。回忆OLS估计表达式可知，不论变量是否重要，OLS表达式
$$ \hat{\beta}^{ols} = (X^TX)^{-1}X^Ty $$
没有体现出变量是否重要的差异性。
因此，大家开始采用的是子集回归法 (subset regression) 而不是OLS方法。
子集回归的思想很朴素。给定数据集，我们从中选出最佳变量组合建立回归模型即可 (全子集回归)。不过，考虑变量的所有组合，这和变量数 $p$ 呈 $2^p$ 关系。这种指数级的关系使得全子集回归成本高昂、不切实际。因此，大家想到一种折衷的办法——考虑一部分变量组合，找到其中的最优模型。这种思路类似优化理论的局部最优解，不过也可以得到较为满意的结果。相关的方法有逐步回归 (stepwise regression)等。
子集回归有个严重缺陷：当原始数据中加入或者剔除某些变量后，子集回归得到的新模型可能发生巨大的变动。 比如，研究基因和某个疾病之间的关系。先搜集100个基因用子集回归得到3个重要的基因。后面研究后发现另外几个基因可能也和该疾病有关，但用子集回归重新建模后，原先的3个基因竟然没有出现在新模型中！这让人对该方法的可靠性产生怀疑。
于是，就轮到LASSO上场了。
3 LASSO模型
LASSO是Least Absolute Shrinkage and Selection Operator 的缩写。据LASSO的提出者Tibshirani本人说，LASSO是受到Breiman (1995) 提出的Nonnegative Garotte方法启发而得到的。NG方法首先建立原始数据的OLS回归，然后对得到的系数施加一个带有非负权重的约束惩罚，使得某些影响较小的系数被直接压缩到零，从而实现变量选择的目的。
$$ \begin{split} NG:\quad&\min \sum _{i=1}^n (y_i - \sum _{j=1}^p t_j \hat{\beta} _j^{ols} x _{ij})^2 \newline &t_j\geq 0,\quad\sum _{j=1}^p t_j \leq s \end{split} $$
NG方法基于OLS，因而深受OLS影响。换句话说，OLS做不好时，NG也好不到哪里去。因此，Tibshirani抛弃了“两步估计”，直接给系数添加了一个绝对值约束，这就是LASSO模型：
$$ \begin{split} LASSO:\quad &\min \sum _{i=1}^n (y_i - \sum _{j=1}^p \hat{\beta}_j x _{ij})^2 \newline &\sum _{j=1}^p \vert \hat{\beta}_j \vert \leq s \end{split} $$
两个模型形式确实相似。不过，LASSO和Hoerl与Kernnard在1970揭示的岭回归 (Ridge Regression) 更为相似，仅仅是罚项由绝对值变成平方项：
$$ \begin{split} Ridge:\quad &\min \sum _{i=1}^n (y_i - \sum _{j=1}^p \hat{\beta}_j x _{ij})^2 \newline &\sum _{j=1}^p \hat{\beta}_j^2 \leq s \end{split} $$
因此，这么一个看似极其微小的改动，究竟神奇在哪里呢？
- 参数s
模型的参数s是需要事先给定的一个非负常数。以LASSO模型为例，如果 $s=0$，那么显然每个模型的系数都是零；如果 $s$ 趋于正无穷，那么约束相当于没加，结果就是OLS估计。当 $s$ 取一个不太大的数时，比如OLS系数估计的绝对值之和的一半。显然此时LASSO无法得到OLS估计，它需要降低某些系数的绝对值以满足约束条件。所以，参数 $s$ 实际控制了系数的压缩程度。
OLS是无偏估计，但在变量数较多时估计量的MSE较大 (估计的系数不可靠)。MSE由估计量的方差和偏度组成，因此牺牲无偏性可以换来估计量方差显著减小，从而得到一个更加可靠的估计。
- 罚函数
LASSO约束系数从而提升估计量可靠性，你用这个思路分析岭回归也有同样的结论。那么，我们再来看两者罚函数形式的区别。考虑二变量情形，画出LASSO和Ridge模型的几何示意图有
两者目标函数都是二次损失，对应一个椭圆。LASSO罚函数采用的是绝对值优化，可行域的边界有很多“尖点” (不可导点)。因此，目标函数优先和尖点触碰，也就形成了稀疏解 (某些系数为零的解)。 然而Ridge罚函数边界光滑，故而难以形成稀疏解。
于是，我们可以看到LASSO的美妙之处：罚函数的小小改动就可以同时提升模型的性能、估计系数以及实现变量选择。
LASSO最早发表在统计“四大天王”期刊之一的JRSSB上。原文没有一个理论证明，但仍然影响深远。可见，创新不见得一定要平地起高楼，站在巨人的肩膀上完全可以有卓越的突破。
4 LASSO模型的求解
4.1 设计阵正交
一般来说，求解LASSO这种带约束的凸优化问题，常常先考虑相应的拉格朗日函数。对于LASSO，拉格朗日函数为
$$ LASSO:\quad \min \sum_{i=1}^n (y_i - \sum_{j=1}^p \hat{\beta}_j x_{ij})^2 + \lambda \sum_{j=1}^p \vert \hat{\beta}_j \vert $$
其中 $\lambda$ 是非负的拉格朗日乘数。如果我们假设数据矩阵 $X$ 是一个列正交矩阵，我们就能得到LASSO的显示解为
$$ \begin{split} \hat{\beta} _j^{LASSO} &= sign(\hat{\beta} _j^{ols})\left(\hat{\beta} _j^{ols} - \frac{\lambda}{2} \right) _+ \newline (x) _+ &= \begin{cases}x,&x\geq 0 \newline 0,&x < 0 \end{cases} \end{split} $$
其中sign()表示取符号。在X列正交时，我们顺便求出Ridge的系数估计表达式为
$$ \hat{\beta}_j^{Ridge} = \frac{1}{1+\lambda}\hat{\beta}_j^{ols} $$
可见，LASSO的解其实就是对最小二乘的解做了一个压缩：对于较小的那部分系数，直接将其压缩到零；对于较大的系数，作了一个平移压缩， 如左图中橙色的线所示。Ridge估计就是对OLS估计做了一个整体的放缩，并没有直接压缩为零的过程， 如右图中蓝色线所示。
因此，我们在这里也能看到罚函数的微小差异带来的巨大不同。绝对值惩罚在零点的不光滑性是LASSO产生稀疏性的本质。
4.1 LARS算法
实际中，数据矩阵X一般不列正交，此时LASSO没有显示解。
LASSO模型罚函数的奇异性像一柄双刃剑：它可以同时实现系数估计和变量选择，但也让基于梯度的算法失效而令其难以有效计算。事实上，LASSO问世初期没有得到人们关注的最大原因就是不好算！
据说，Tibshrani当时为了寻找一个有效的算法，曾四处拜访名家无果。最后，Tibshirani找到了自己的恩师Efron，大师出手一举摆平了此事。这就是著名的最小角回归 (Least Angle Regression, LARS)。
有了LARS算法，LASSO像开了挂一般，迅速得到了蓬勃的发展。确实，统计模型往往来源于实际问题，如果模型仅仅停留在理论的空中楼阁，往往难以拥有生命力。当然，LARS不仅是一个算法，重要的是其从理论上揭示稀疏性 (Sparsity) 的本质。有兴趣的朋友可以拜读一下大师这篇90多页的大作！
4.2 坐标下降法
LARS算法将LASSO的计算复杂度降到相当于OLS的复杂度，这基本就将其计算效率提升到极致。不过，这个记录后来被Friedman打破了。Friedman采用的是坐标下降法 (Coordinate Decent Algorithm, CD)。
我们思考这么一个问题：如果一个函数在某点处各个维度取到了最小值，那么该点是不是函数的最小值点？
假如该结论成立，那么全局优化的过程就可以从每个维度进行一维优化即可。这显然极大地提升了模型的求解速度。
可惜，上述问题在一般情形下的答案是不一定。幸运的是，在LASSO模型中这个结论是成立的！因此，我们优化LASSO可以从每个维度 (几何上看就是每个坐标) 进行优化，采用迭代计算收敛到最优解即可。这就大大简化了模型的求解过程。
考虑第 $j$ 个系数，我们根据KKT条件得到该出的最优解满足
$$ \begin{split} \hat{\beta}_j &= \text{sign}(x_j^T r_j) \frac{(x _j^T r _j - \frac{\lambda}{2}) _+}{\Vert x _j \Vert _2^2}\newline r_j &= y - \sum _{k\neq j}^p x_k \hat{\beta}_k \end{split} $$
于是，我们便可以利用上式来迭代求解LASSO模型。
4.3 正则化参数的选取
LASSO模型中的参数 $s$ 控制系数的压缩程度，拉格朗日函数中则由参数 $\lambda$ 扮演该角色。通常，我们称这种参数为正则化参数，它们需要事先给定。
正则化参数给大了，系数压缩的过偏，影响模型的表现；正则化参数给小了，系数压缩不明显，提升模型表现的作用不大。
那么问题来了：怎么设定一个最优参数？
通常，我们首先选定一系列的正则化候选参数；然后针对每个参数计算某一准则 (如AIC、BIC等) 的结果；最后根据准则确定最优参数的大小。
5 R与LASSO
在R中，我们可以使用lars进行LASSO求解。使用前需要事先安装lars程序包，它由Efron和Hastie编写。此外，我们还可以使用glmnet求解LASSO模型。使用之前需要事先安装glmnet程序包，它由Tibshirani、Friedman等编写的。它们的基本用法如下
# lars求解lasso的基本用法
fit <- lars(x, y, type = "lasso", intercept = T)
# glmnet求解lasso的基本用法
fit <- glmnet(x, y, alpha = 1, lambda = NULL, intercept = T)
predict.glmnet(fit, newx)
其中x为数据矩阵，y是响应变量，lambda是正则化参数。当alpha=1时，glmnet计算LASSO模型；当alpha=0时，glmnet计算Ridge模型；当0<alpha<1时，glmnet计算Elastic Net模型，一个LASSO的改进模型。此外，模型的截距项参数intercept默认是T，改成F即为无截距。
glmnet本身优化了模型的计算。对于给定的数据集，glmnet自动计算一系列lambda (默认100个) 对应的系数估计，且每次计算都会利用到上次估计的结果。因此，计算一系列lambda的速度比单独给lambda赋值计算的总耗时要快的多。如果你想用CV (cross validation) 的方式寻找最优参数，也可以自己手动设置参数或者使用cv.glmnet函数。
6 模拟分析
我们在R当中产生一个含有10个变量和1000个观测的数据矩阵，系数为 (1,0,1,0,…,1,0) 且模型无截距，具体代码如下
set.seed(123)
samplen <- 1000
samplep <- 10
x <- matrix(rnorm(samplep*samplen), ncol = samplep)
b <- rep(c(1,0),5)
y <- x%*%b + rnorm(samplen)
6.1 使用lars求解
使用lars函数求解模型，并用plot函数画出系数估计的路径图
library(lars)
fit <- lars(x, y, type = "lasso", intercept = F)
plot(fit)
图中的竖线表示lars算法的迭代次数，非零值对应的即是被选入的变量，图的上刻度数字表示本次迭代选入的变量数，图的右刻度数字表示最终选入的变量名。
我们可以使用summary函数查看每次迭代的具体情况，可以按照最小 $C_p$ 值得到最优的模型。
看来 $C_p$ 准则得到的模型和真实模型有点出入。实际估计的系数可以看到，非零系数基本接近1，零系数基本接近零，但 $C_p$ 准则在本案例中稀疏的能力不佳。
> fit$beta[which.min(fit$Cp),]
[1] 0.98459718 0.04174324 1.02200983 0.05431813 1.02213893 0.00000000
[7] 0.95807588 0.01790684 1.02019194 0.04779414
6.2使用glmnet求解
使用glmnet求解，再用print函数输出结果：每个lambda对应的模型选入的变量数Df和模型解释的偏差%Dev (R方)。
> library(glmnet)
> fit <- glmnet(x, y, intercept = F)
> print(fit)
Call: glmnet(x = x, y = y, intercept = F)
Df %Dev Lambda
1 0 0.00000 1.23800
2 1 0.04079 1.12800
3 3 0.08998 1.02800
4 5 0.18710 0.93620
5 5 0.29830 0.85310
6 5 0.39060 0.77730
7 5 0.46730 0.70820
8 5 0.53090 0.64530
9 5 0.58370 0.58800
10 5 0.62760 0.53570
···
选取Df=5和最大%Dev对应得lambda值，利用coef函数得到相应得系数估计结果。
> coef(fit, s = 0.06305)
11 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) .
V1 0.9243286
V2 .
V3 0.9650517
V4 .
V5 0.9656673
V6 .
V7 0.8985474
V8 .
V9 0.9700144
V10 .
看来结果不错。不过，我们利用cv.glmnet筛选最优参数发现，最优模型也是对应9个变量。
> cv_fit <- cv.glmnet(x, y)
> op_lam <- cv_fit$lambda[which.min(cv_fit$cvm)]
> coef(fit, s = op_lam)
11 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) .
V1 0.97767537
V2 0.03679953
V3 1.01549521
V4 0.04761603
V5 1.01535413
V6 .
V7 0.95157838
V8 0.01017659
V9 1.01429797
V10 0.04150092
看来我随便生成的这个数据把这两个准则给难住了。