本帖最后由 ysq 于 2016-3-7 20:23 编辑
岭回归和lasso---回归的拓展，目的：解决回归中的重大疑难问题“排除多重共线性，进行变量的选择”。
一、回归问题之数学背景
多元线性回归的最小二乘解（无偏估计）
将原先很复杂回归问题，写成矩阵形式，以简化回归公式和对回归问题的描述。
n个样本，p个变量，X，y已知。对数据中心化、标准化处理后，可以去掉截距项。
矩阵形式的多元线性模型为
，求解β，使得误差项ε能达到最低。
残差ε为
，残差平方和RSS为file:///C:\Users\Administrator\AppData\Roaming\Tencent\Users\277418979\QQ\WinTemp\RichOle\I_)0))93PL7IQRA87Q]S{6U.png
∴多元线性回归问题变为求解β，从而使残差平方和极小值问题（关于系数向量β的二次函数极值问题）。
求解方法：配方法（最小二乘法）和求偏导法。
残差向量的几何意义：响应y向量到由p个x向量组成的超平面的距离向量。
残差平方和几何意义：残差向量长度的平方。
可以证明β的最小二乘估计为
X矩阵不存在广义逆（即奇异性）的情况：
1）X本身存在线性相关关系（即多重共线性），即非满秩矩阵。
当采样值误差造成本身线性相关的样本矩阵仍然可以求出逆阵时，此时的逆阵非常不稳定，所求的解也没有什么意义。
2）当变量比样本多，即p>n时.
这时，回归系数会变得很大，无法求解。在统计学上，可证明β的最小二乘解为无偏估计，即多次得到的采样值X而计算出来的多个系数估计值向量
的平均值将无限接近于真实值向量β。
岭回归和lasso回归属于有偏估计。
二、岭回归（Ridge Regression，RR）
1962年由Heer首先提出，1970年后他不肯纳德合作进一步发展了该方法。
先对数据做标准化，为了记号方便，标准化后癿样本数据学习集仍然用X表示。
思路：在原先的β的最小二乘估计中加一个小扰动λI，是原先无法求广义逆的情况变成可以求出其广义逆，使得问题稳定并得以求解。
（式1）
上式中，此时随着λ而变化，是λ的函数，称为关于岭参数λ的岭回归估计。
此时极值问题的求解函数也不再是原来的残差平方和表达式了，而变成了下式：
(式2)
对上式用偏导数求极值，结果就是 （式1）。后面这一项
称为惩罚函数，它保证了β值不会变的很大。岭参数λ不同，岭回归系数也会不同。
在数学上可以证明，（式2）代表的极值问题和下式（式3）的极值问题完全一致：
（式3）
但（式3）可以更好地描述岭回归的几何意义。
（式3）中第一行与多元线性回归的残差平方和Q(β)相同，第二行是对第一行的一个约束项，也就是说所有系数β平方和应该<某个阈值t。
（式2）中的λ和（式3）中的t是一一对应的关系，也就是说两个问题是等价的。
以两个变量为例，解释岭回归的几何意义1）没有约束项时
系数β1和β2已经经过标准化。残差平方和RSS可以表示为β1和β2的一个二次函数，数学上可以用一个抛物面表示。
图1
2）岭回归时
约束项为β12 β22≤t，对应着投影为β1和β2平面上的一个圆，即下图中的圆柱。
图2
该圆柱与抛物面的交点对应的β1、β2值，即为满足约束项条件下的能取得的最小的β1和β2.
从β1β2平面理解，即为抛物面等高线在水平面的投影和圆的交点，如下图所示
图3
可见岭回归解与原先的最小二乘解是有一定距离的。
岭回归性质
1）当岭参数为0，得到最小二乘解。
2）当岭参数λ趋向更大时，岭回归系数β估计趋向于0。---为了使（式2）约束项很小，才能实现（式2）整体岭函数达到极小。
是回归参数β的
有偏估计。它的结果是使得残差平和变大，但是会使系数检验变好，即
R语言summary结果中变量后的*变多。
4）在认为岭参数λ是与y无关的常数时，是最小二乘估计的一个线性变换，也是y的线性函数。
但在实际应用中，由于λ总是要通过数据确定，因此λ也依赖于y、因此从本质上说，并非的线性变换，也非y的线性函数。
5）对于回归系数向量来说，有偏估计回归系数向量长度<无偏估计回归系数向量长度，，即比理想值要短。
6）存在某一个λ，使得它所对应的的MSE（估计向量的均方误差）<最小二乘法对应估计向量的的MSE。
即 存在λ>0，使得
。
理解 ： 和β平均值有偏差，但不能排除局部可以找到一个比更接近β。
岭迹图
是λ的函数，岭迹图的横坐标为λ，纵坐标为β(λ)。而β(λ)是一个向量，由β1(λ)、β2(λ)、...等很多分量组成，每一个分量都是λ的函数，将每一个分量分别用一条线。
当不存在奇异性时，岭迹应是稳定地逐渐趋向于0
岭迹图作用：
1）观察λ最佳取值；
2）观察变量是否有多重共线性；
可见，在λ很小时，通常各β系数取值较大；而如果λ=0，则跟普通意义的多元线性回归的最小二乘解完全一样；当λ略有增大，则各β系数取值迅速减小，即从不稳定趋于稳定。
上图类似喇叭形状的岭迹图，一般都存在多重共线性。
λ的选择：一般通过观察，选取喇叭口附近的值，此时各β值已趋于稳定，但总的RSS又不是很大。
选择变量：删除那些β取值一直趋于0的变量。
注意：用岭迹图筛选变量并非十分靠谱。
岭参数的一般选择原则
选择λ值，使到
1）各回归系数癿岭估计基本稳定；
2）用最小二乘估计时符号不合理癿回归系数，其岭估计的符号变得合理；
3）回归系数没有不合乎实际意义癿绝对值；
4）残差平方和增大不太多。 一般λ越大，系数β会出现稳定的假象，但是残差平方和也会更大。
取λ的方法比较多，但是结果差异较大。这是岭回归的弱点之一。
岭回归选择变量的原则（不靠谱，仅供参考）
1）在岭回归中设计矩阵X已经中心化和标准化了，这样可以直接比较标准化岭回归系数癿大小。可以剔除掉标准化岭回归系数比较稳定且绝对值很小癿自变量。
2）随着λ的增加，回归系数不稳定，震动趋于零的自变量也可以剔除。
3）如果依照上述去掉变量的原则，有若干个回归系数不稳定，究竟去掉几个，去掉哪几个，这幵无一般原则可循，这需根据去掉某个变量后重新进行岭回归分析的效果来确定。
用R做岭回归
library(MASS)#岭回归在MASS包中。
longley #内置数据集，有关国民经济情况的数据，以多重共线性较强著称
summary(fm1<-lm(Employed~.,data=longley)) #最小二乘估计的多元线性回归
#结果可见，R^2很高，但是系数检验不是非常理想
names(longley)[1]<-'y'
lm.ridge(y~.,longley) #此时，仍为线性回归
plot(lm.ridge(y~.,longley,lambda=seq(0,0.1,0.001))) #加了参数lambda的描述后才画出响应的岭迹图
#由于lambda趋于0时，出现了不稳定的情况，所以可以断定变量中存在多重共线性
select(lm.ridge(y~.,longley,lambda=seq(0,0.1,0.001))) #用select函数可算lambda值，结果给出了3种方法算的的lambda的估计值
modified HKB estimator is 0.006836982
modified L-W estimator is 0.05267247
smallest value of GCV at 0.006
#以上结果通常取GCV估计，或者观察大多数方法趋近哪个值。
R的ridge包
###在R语言中已下架，没有安装成功
下载地址 https://cran.r-project.org/src/contrib/Archive/ridge/?C=D;O=A
岭回归缺陷
1.主要靠目测选择岭参数
2.计算岭参数时，各种方法结果差异较大
所以一般认为，岭迹图只能看多重共线性，却很难做变量筛选。
三、LASSO
Tibshirani(1996)提出了Lasso(The Least Absolute Shrinkage and Selectionatoroperator)
算法，这里 Absolute 指的绝对值。
Shrinkage收缩的含义：
类似于图3，即系数收缩在一定区域内（比如圆内）。
主要思想：
通过构造一个一阶惩罚函数获得一个精炼的模型；通过最终确定一些指标（变量）癿系数为零（岭回归估计系数等于0癿机会微乎其微，造成筛选变量困难），解释力很强。
擅长处理具有多重共线性癿数据，筛选变量，与岭回归一样是有偏估计。
LASSO vs 岭回归
1）对于岭回归的两个表达式：
一方面可以将其变成一个最小二乘问题
另一方面可以将它解释成一个带约束项的系数优化问题。
λ增大的过程就是t减小的过程，该图也说明了岭回归系数估计值为什么通常不为0，因为随着抛物面的扩展，它与约束圆的交点可能在圆周上的任意位置，除非交点恰好位于某个坐标轴或坐标平面上，否则大多数情况交点对应的系数值都不为零。再加上λ的选择应使椭球面和圆周的交点恰好在一个坐标平面上，更增加了求解λ的难度。
2）对于LASSO回归
（式3）---第一个表达式与岭回归的区别在于惩罚项上，其惩罚函数是一个β的一阶函数。
（式4）---第二个表达式将上式也转化为一个约束项。t和λ也是一一对应的，λ增大的过程就是t减小的过程。
由于方框的顶点更容易交于抛物面，也就是lasso更易求解，而该顶点对应的很多系数为0，也就是起到了筛选变量的目的。
回归过程对比：
左图为岭回归，右图为lasso回归
横轴越往左，自由度越小（即圆或方框在收缩的过程），λ越大，系数（即纵轴）会越趋于0。但是岭回归没有系数真正为0，但lasso的不断有系数变为0.
四、更一般化的模型
不同q对应的约束域形状
五、弹性网
Zouand Hastie (2005)提出elasticnet，介于岭回归和lasso回归之间，现在被认为是处理多重共线性和变量筛选最好的收缩方法，而且损失的精度不会太多。
六、LAR（最小角回归）Lasso回归中第一个表达式（式3）用偏导求极值时，存在部分点不可导的情况（如方框的尖点），如何解决？
Least Angel Regression ，参考书The Elements of Statistical Learning .pdf
Efron于2004年提出癿一种变量选择癿方法，类似于向前逐步回归(Forward Stepwise)癿形式，最初用于解决传统的线性回归问题，有清晰的几何意义。
是lasso regression癿一种高效解法。
向前逐步回归(Forward Stepwise)丌同点在于，Forward Stepwise每次都是根据选择癿变量子集，完全拟合出线性模型，计算出RSS，再设计统计量（如AIC）对较高癿模型复杂度作出惩罚，而LAR是每次先找出和因变量相关度最高癿那个变量, 再沿着LSE癿方向一点点调整这个predictor癿系数，在这个过程中，这个变量和残差癿相关系数会逐渐减小，等到这个相关性没那么显著癿时候，就要选进新癿相关性最高癿变量，然后重新沿着LSE癿方向进行变动。而到最后，所有变量都被选中，就和LSE相同了。
左图为LAR逐步加上变量的过程（从左往右看），右图为LASSO变量逐渐淘汰的收缩过程（从右往左看）。
对比两幅图，非常类似。所以可以用LAR方法来计算LASSO，该方法完全是线性解决方法，没有迭代的过程。
相关系数的几何意义
设变量y=[y1,y2,...yn]; 变量x=[x1,x2,...,xn].
其相关系数为
，其中cov--协方差、var---方差。
如果对x和y进行中心化、标准化，则var(y)=var(x)=1,相关系数变为x1y1 x2y2 .... xnyn，即为向量x和y的内积=||x||*||y||*cos θ，其中θ为x和y的夹角。而对于标准化和中心化后的x和y，则有||x||=||y||=1，所以此时x和y的内积就是它们夹角的余弦。
---如果x和y向量很像，几乎重合，则夹角θ=0，也就是相关系数=内积=1，此时称为高度相关.
---如果x和y相关程度很低，则表现出来的x和y向量相互垂直，相关系数=0.
---- 如果相关系数=-1，标明x和y呈180°，即负相关。
LAR算法及几何意义
参考书The Elements of Statistical Learning .pdf的74页。LAR和Lasso的区别以及LAR解Lasso的修正
参考书The Elements of Statistical Learning .pdf的76页。
LAR过程图解
有6个变量，最先加入与残差向量相关系数最大的浅蓝色v2，在v2变化过程中，相关系数越变越小，直到等于深蓝色的v6，于是加入v6，沿着v2与v6的最小角方向（即向量角分线方向）前进，此后v2和v6与残差向量的相关系数是共同变化的，即两者合并变化，使得相关系数越来越小，直到加入黑色v4为止，三个变量一起变化，...，一直打到最小二乘解为止，此时残差向量与所有变量的相关系数都为0，即与他们都垂直。
横坐标L1 Length表示：从原点开始走了多长距离，就是绝对值距离，L1范数。
R语言中对LAR的实现
install.packages('lars') #lars包
longley #用longley数据集，它是一个著名的多重共线性例子
w=as.matrix(longley) #将数据集转换为一个矩阵
laa=lars(w[,2:7],w[,1]) #w的2:7列为自变量，第1列为因变量
laa #显示LAR回归过程
Call:
lars(x = w[, 2:7], y = w[, 1])
R-squared: 0.993
Sequence of LASSO moves:
GNP Year Armed.Forces Unemployed Employed Population Year Employed Employed Year Employed Employed
Var 1 5 3 2 6 4 -5 -6 6 5 -6 6
Step 1 2 3 4 5 6 7 8 9 10 11 12
#第一步加入1号变量，第二步加入5号变量，...，第七步5号变量被去掉，第八步6号变量被去掉，...
plot(laa) #画lasso回归过程图
summary(laa)
LARS/LASSO
Call: lars(x = w[, 2:7], y = w[, 1])
Df Rss Cp
0 1 1746.86 1210.0561
1 2 1439.51 996.6871
2 3 32.31 12.6400
3 4 23.18 8.2425
4 5 22.91 10.0505
5 6 22.63 11.8595
6 7 18.04 10.6409
7 6 14.74 6.3262
8 5 13.54 3.4848
9 6 13.27 5.2974
10 7 13.01 7.1189
11 6 12.93 5.0624
12 7 12.84 7.0000
#以上结果显示了每一步的残差平方和RSS和多重共线性指标Cp（Mallows's Cp http://en.wikipedia.org/wiki/Mallows%27_Cp）
#Cp越小，多重共线性越小，因此结果以第八步为准，即只剩下第1、2、3、4个变量。
