分类目录归档:回归分析

回归分析、线性模型、非线性模型、广义线性模型

共轭梯度法计算回归

共轭梯度示意图(图片来源:维基百科)
轮回眼 共轭梯度示意图(图片来源:维基百科

引子

之所以写这篇文章,是因为前几天统计之都的微信群里有同学提了一个问题,想要对一个很大的数据集做回归。然后大家纷纷给出了自己的建议,而我觉得共轭梯度算回归的方法跟这个背景比较契合,所以就正好写成一篇小文,与大家分享一下。

说到算回归,或许大家都会觉得这个问题太过简单了,如果用 $X$ 表示自变量矩阵,$y$ 表示因变量向量,那么回归系数的最小二乘解就是 $\hat{\beta}=(X’X)^{-1}X’y$。(本文完)



哎等等,别真走啊,我们的主角共轭梯度还没出场呢。前面的这个算系数的公式确实非常简洁、优雅、纯天然、不做作,但要往里面深究的话,还是有很多问题值得挖掘的。

最简单暴力的方法,就是从左向右,依次计算矩阵乘法,矩阵求逆,又一个矩阵乘法,最后是矩阵和向量的乘法。如果你就是这么算的,那么可以先默默地去面壁两分钟了。

更合理的方法,要么是对 $X’X$ 进行 Cholesky 分解,要么是对 $X$ 进行 QR 分解,它们基本上是现在算回归的软件中最常见的方法。关于暴力方法和矩阵分解方法的介绍和对比,可以参见这个B站上的视频。(什么?你问我这么严肃的话题为什么要放B站上?因为大部分时间都是在吐槽啊)

好,刚才去面壁的同学现在应该已经回来了,我们继续。前面这些通过矩阵运算求回归系数的方法,我们可以统称为直接法。叫这个名字,是因为它们都可以在确定数目的步骤内得到最终的结果。而与之相对的,则叫做迭代法,意思是通过不断更新已经得到的结果,来逐渐逼近真实的取值。打个比方,你想要知道一瓶82年的拉菲值多少钱,直接法就是去做调研,原料值多少,品牌值多少,加工费多少,运输费多少……然后加总起来得到最终的定价;而迭代法就是去问酒庄老板,你先随便蒙一个数,然后老板告诉你高了还是低了,反复循环,总能猜个八九不离十。

说到这里,你自然要问了,既然算回归的软件大都是用直接法,为什么还要考虑迭代法?莫非直接法有什么不好的地方?这就说到问题的点子上了。

继续阅读共轭梯度法计算回归

为什么我不是R方的粉丝

本文翻译自 John Myles White 的博客 Why I’m Not a Fan of R-Squared。翻译工作已经获得作者授权同意。

本文大意

人们通常喜欢用 $R^2$ 作为评判模型拟合好坏的标准。与 MSE MAD 不同,$R^2$ 不只是模型误差的函数,它的定义中还隐含了两个模型的比较:一个是当前被分析的模型,一个是所谓的常数模型,即只利用因变量均值进行预测的模型。基于此,$R^2$ 回答的是这样一个问题:“我的模型是否比一个常数模型更好?”,然而我们通常想要回答的是另一个完全不同的问题:“我的模型是否比真实的模型更差?

通过一些人为构造的例子我们可以很容易发现,对这两个问题的回答是不可互换的。我们可以构造一个这样的例子,其中我们的模型并不比常数模型好多少,但同时它也并不比真实的模型差多少。同样,我们也可以构造出另一个例子,使得我们的模型远比常数模型要好,但也远比真实模型要差。

与所有的模型比较方法一样,$R^2$ 不单是被比较模型的函数,它也是观测数据的函数。几乎对于所有的模型,都存在一个数据集,使得常数模型与真实模型之间是无法区分开的。具体来说,当使用一个模型区分效能很低的数据集时,$R^2$ 可以任意地向零趋近——即使我们对真实模型计算 $R^2$ 也是如此。因此,我们必须始终记住,$R^2$ 并不能告诉我们模型是否是对真实模型的一个良好近似:$R^2$ 只告诉我们,我们的模型在当前的数据下是否远比一个常数模型要好。 继续阅读为什么我不是R方的粉丝

R 中大型数据集的回归

众所周知,R 是一个依赖于内存的软件,就是说一般情况下,数据集都会被整个地复制到内存之中再被处理。对于小型或者中型的数据集,这样处理当然没有什么问题。但是对于大型的数据集,例如网上抓取的金融类型时间序列数据或者一些日志数据,这样做就有很多因为内存不足导致的问题了。
继续阅读R 中大型数据集的回归

分组最小角回归算法(group LARS)

继续前两篇博文中对于最小角回归(LARS)和lasso的介绍。在这篇文章中,我打算介绍一下分组最小角回归算法(Group LARS)。本文的主要观点均来自Ming Yuan和Yi Lin二人2006合作发表在JRSSB上的论文Model selection and estimation in regression with grouped variables.

首先,我想说明一下,为何要引入分组变量(grouped variable)的概念。举一个简单的例子,在可加的多项式模型中,每一项都是多项式。这个多项式有可能可以通过最初的变量的线性组合来表达。在进行这种类型的回归中,挑选重要的变量其实质是挑选重要的因子(factor),而因子则是最初的那些变量的线性组合。分组变量的回归问题,实际上就是我们一般所说的回归问题的推广。如果我们把每一个单独的变量都看成一个因子,那么这种情况下的回归就是一般意义下的回归。下面用公式更加直白的说明这个问题:

$Y=\sum_{j=1}^JX_j\beta_j+e$

其中$Y$是个$n$维向量,$e\~N_n(0,\sigma^2I)$.$X_j$是$n\times p_j$矩阵,代表的是第j个因子(factor,是变量variables的线性组合)。$\beta_j$是$p_j$维的系数向量。依然假定$Y$是中心化的,$X_j$是中心化并且正交化的($X_j’X_j=I$)。这个就是分组变量的回归模型。

在小弟之前的文章中介绍过最小角回归(LARS)算法,对于变量选择来说,这种算法是比较易于计算机实现的,而且比传统的向前逐步回归(forward selection)更加谨慎,但是没有逐段回归(stagewise)那么谨小慎微。对于分组变量的回归模型选择问题来说,直接套用LARS算法可能会挑选出一些不必要的因子进入回归模型中。因此在文章中,作者提出了分组最小角回归(group LARS).这种算法对最小角回归进行了推广,下面简单介绍一下这种算法。

推广是循序渐进的,因此,我们先来看看$p_j$均相等的情况(注意这里$p_j$未必等于1)。在传统的LARS算法中,我们是要保证入选变量和当前残差的角度都相同。如今,我们就需要保证入选的因子们(再次强调,这里是因子,是原始变量的线性组合)与当前残差的角度相同。直观地来定义一个角度$\theta(r,X_j)$,这个角度就是向量$r$与$X_j$中的变量所构成的空间的夹角,也就是$r$与它在$X_j$中的变量所构成空间的投影的夹角。(再废话两句,这个夹角越小,就意味着$r$与$X_j$中的变量构成的空间中的相关系数越大。)当$r$是当前残差的时候,我们就可以确定求解路径(solution path)为当前残差在$X_j$中的变量所构成的空间上的投影的方向。

上面说的可能还不够清楚,下面我分步骤具体介绍每一步的做法。

  1. 最开始时,回归变量集依然是空集,此时的残差就是$Y$,我们寻找和$Y$夹角最小的$X_j$.记它为$X_{j1}$.(这里要插一句,就是如何度量这个角度。用这个角的余弦值的平方即可:$cos^2(\theta(r,X_j))=\frac{\|X_j’r\|^2}{\|r\|^2}$,角度越小,这个值越大。)
  2. 确定求解路径为$Y$在$X_{j1}$中的变量所构成的空间上的投影的方向。然后在这个路径上前进,直到出现第二个因子$X_{j2}$,使得当前残差$r$与这两个因子的夹角同样小。也就是$\|X_{j1}’r\|^2=\|X_{j2}’r\|^2$。
  3. 计算新的求解路径,这个路径其实就是当前残差在$X_{j1}$和$X_{j2}$的变量所构成的空间中投影的方向。然后再沿着这条路径前进,直到第三个因子出现,使得当前残差和这三个因子的夹角都是同样小。

然后继续按照上面的步骤,计算出求解路径,然后继续前进。直到找到第四个满足条件的因子。

这些步骤和原始的LARS的进程是同样的。只不过在确定求解方向时有了变化。

上面介绍的都是$p_j$相同的情况,当$p_j$不同的时候,只要对上面的过程做一个小的改动即可。用$\frac{\|X_j’r\|^2}{p_j}$代替$\|X_j’r\|^2$,其他过程都不变。这种对于LARS算法的推广是比较直观的。而且M.Yuan和Y.Lin在文章中还指出,这种算法相对于正交化$X_j$的过程是稳定的,通过Gram-Schmidt方法进行正交化本来是依赖于选择的参数的,但是分组LARS确定的求解方向并不依赖该参数。

修正的LARS算法和lasso

在小弟的上一篇文章中,简单的介绍了LARS算法是怎么回事。主要参考的是Efron等人的经典文章least angle regression。在这篇文章中,还提到了一些有趣的看法,比如如何用LARS算法来求解lasso estimate和forward stagewise estimate。这种看法将我对于模型选择的认识提升了一个层次。在这个更高的层次下看回归的变量选择过程,似乎能有一些更加创新的想法。

lasso estimate的提出是Tibshirani在1996年JRSSB上的一篇文章Regression shrinkage and selection via lasso。所谓lasso,其全称是least absolute shrinkage and selection operator。其想法可以用如下的最优化问题来表述:

在限制了$\sum_{j=1}^J|\hat{\beta_j}|\leq t$的情况下,求使得残差平和$\|y-X\hat{\beta}\|^2$达到最小的回归系数的估值。

我们熟悉如何求解限制条件为等号时,回归方程的求解。也就是用lagrange乘子法求解。但是对于这种,限制条件是不等号的情况,该如何求解,则有两种想法。第一种,也是我比较倾向于的方法,是利用计算机程序,对$t$从$0$开始,不断慢慢增加它的值,然后对每个$t$,求限制条件为等号时候的回归系数的估计,从而可以以$t$的值为横轴,作出一系列的回归系数向量的估计值,这一系列的回归系数的估计值就是lasso estimation。另外一种想法,是借助与最优化问题中的KKT条件,用某个黑箱式的算法,求解。(本人对于最优化方面的东西实在是不很熟悉,故不在此弄斧,只求抛砖引玉,能有高手给出这种想法的具体介绍。)

lasso estimate具有shrinkage和selection两种功能,shrinkage这个不用多讲,本科期间学过回归分析的同学应该都知道岭估计会有shrinkage的功效,lasso也同样。关于selection功能,Tibshirani提出,当$t$值小到一定程度的时候,lasso estimate会使得某些回归系数的估值是$0$,这确实是起到了变量选择的作用。当$t$不断增大时,选入回归模型的变量会逐渐增多,当$t$增大到某个值时,所有变量都入选了回归模型,这个时候得到的回归模型的系数是通常意义下的最小二乘估计。从这个角度上来看,lasso也可以看做是一种逐步回归的过程。

在我的上一篇文章中,提到了Efron对于逐步回归的一种看法,就是在某个标准之下(比如LARS的标准就是要保证当前残差和已入选变量之间的相关系数相等,也就是当前残差在已入选变量的构成空间中的投影,是那些变量的角平分线)选择一条solution path,在这个solution path上proceed,不断吸收新的变量进入,然后调整solution path 继续proceed。那么对于求解lasso的算法,也有一个相应的对应。Efron提出了一种修正的LARS算法,可以用修正的LARS算法来求解所有的lasso estimates。下面我介绍一下这种修正的LARS算法。

首先假设我们已经完成了几步LARS steps。这时候,我们已经有了一个回归变量集,我们记这个回归变量集为$X_{\mathscr{A}}$。这个集合就对应着一个对于$Y$的估计,我们记为$\hat{\mu}_{\mathscr{A}}$。这个估值对应着一个lasso方法对于响应的估值(这里我认为LARS估值和lasso估值应该是一样的),lasso的估值,对应着回归系数的lasso估值,回归系数向量的lasso估值我们记为$\hat{\beta}$。

为了继续进行下一步,我们先给出一个向量的表达式,然后再解释一下它

$w_{A}=(1_{A}'(X_{A}’X_{A})^{-1}1_{A})^{-\frac{1}{2}}(X_{A}’X_{A})^{-1}1_{A}$.

$X_{A}w_{A}$就是LARS算法的在当前回归变量集下的solution path。那么我们可以把$w_{A}$作为$\beta$的proceed的path。Efron定义了一个向量$\hat{d}$,这个向量的元素是$s_jw_j$,其中$s_j$是入选变量$x_j$与当前残差的相关系数的符号,也是$\hat{\beta_j}$的符号。对于没有入选的变量,他们对应在$\hat{d}$中的元素为0。也就是对应着$\mu(r)=X\beta(r)$,我们有

$\beta_j(r)=\hat{\beta_j}+r\hat{d_j}$

将LARS的solution path对应到lasso estimate的path上,这种对应的想法非常值得借鉴。

很显然,$\beta_j(r)$会在$r_j=-\hat{\beta_j}/\hat{d_j}$处变号。那么对于我们已经有的lasso estimate$\beta(r)$,它中的元素会在最小的的那个大于$0$的$r_j$处变号。我们记之为$\bar{r}$。如果没有$r_j$大于$0$,那么$\bar{r}$就记为无穷大。

对于LARS本身而言,在已经有了如今的回归变量集和当前残差的基础上,我们就会有条solution path,在这个solution path上proceed的最大步记为$\hat{r}$.通过比较$\hat{r}$和$\bar{r}$就会有进一步的想法。Efron的文章证明了如果$\bar{r}$小于$\hat{r}$,则对应于LARS估计的那个$\beta_j(r)$不会成为一个lasso estimation。(这个是因为当前残差和对应变量的相关系数的符号一定是和该变量的系数符号一致才行)。在这种情况下,我们就不能继续在LARS的solution path上继续前进了,为了利用LARS算法求得lasso estimate,Efron提出把$\bar{r}$所对应的那个$r_j$所对应的$x_j$从回归变量中去掉。去掉之后再计算当前残差和当前这些变量集之间的相关系数,从而确定一条新的solution path,继续进行LARS step。这样进行下去,可以通过LARS算法得到所有的lasso estimate。

这个对于LARS的lasso修正算法,被Efron称作“one at a time”条件,也就是每一步都要增加或删掉一个变量。下图显示了用修正了的LARS算法求lasso estimate的过程。

这个图是Efron等人的文章中,对于一个实际数据进行回归得到的。该数据一共有10个变量。图的横轴,是所有回归系数估值的绝对值之和,这个值从$0$增加。左侧的纵轴,是回归系数的估值,右侧纵轴是这些回归系数对应的变量的下标。这个图中,我们可以看到每一个回归系数的path。可以看到第七个变量对应的回归系数在横轴快到3000的时候变为了0,说明到这一步时,该变量被删除掉,之后又被重新添加到了回归变量集中。

下面通过一个简单的模拟,对lars和lasso以及forward stagewise做一个简单的实现。其实在R中已经有了一个名为lars的包,可以实现上述三种回归。

首先,我要模拟的方程为

$y={x^3}_1+{x^2}_1+x_1+\frac{1}{3}{x^3}_2-{x^2}_2+\frac{2}{3}x_2+e$

其中$x_1$和$x_2$是服从二维联合正态分布,均值为零向量,$cov(x_1,x_2)=0.5$,$var(x_1)=var(x_2)=1$,$e$服从$N(0,9)$。我取了50次观测,然后分别通过lasso,lars,以及forward stagewise三种算法进行了回归,其变量的回归路径如下图。

简单的代码我直接贴在本文的最后。从这三个算法的图中,我们并看不出有特别的区别,只能看出一些细小的差别。至于要判断哪种算法更好,则应该因问题而异。也不是本文能够论述的问题了。

对于LARS算法的修正,还可以应用到计算forward stagewise的estimate中,在Efron的文章中也有介绍。他的这种看法,好似凌驾在整个回归变量选择过程之上,从一个更高的角度观察之,给出一种更为一般性的视角。这也就是大牛和一般人之间的差别。读Efron的文章,总有一种让人想要膜拜的冲动。对于模型选择方面的东西,值得挖掘的还很多。Tibshirani在最新的一篇综述性的文章中,给出了lasso的诞生到现今发展的一系列流程。感兴趣的读者,可以去看看这篇文章,在cos论坛上有。链接如下:

http://cos.name/cn/topic/104104

用lars算法做模拟的代码:

利用lars模拟