Penalized linear regression
引言
在统计学中, 输入变量(inputs)通常称作预测变量(predictors), 也可以被叫做自变量(independent variables)。输出变量(outputs)被称作响应变量(responses)，也可以被叫做因变量(dependent variables)。输入变量可以是测量得到或者预设的. 这些变量对一个或多个输出变量有影响. 即利用输入变量去预测输出的值。这样的过程称之为监督学习(supervised learning)。
输出变量的类型, 根据度量可以分为定量的(quantitative)与定性的(qualitative)。定性变量也被称为类别型 (categories)或者离散(discrete)型变量, 也被称作因子(factors)。对于两种类型的输出变量, 考虑使用输入变量去预测输出变量是有意义的。当我们预测定量的输出时被称为回归(regression), 当我们预测定性的输出时被称为分类(classification)。这两个任务有共同点, 特别地, 两者都可以看成是函数逼近。
输入变量也有各种类型, 除了定性和定量以外, 还有第三类有序分类 (ordered categorical), 如 小(small)、中 (medium) 和 大 (large), 在这些值之间存在顺序, 但是没有合适的度量概念.。
把输入变量用符号\(X\)来表示。如果\(X\)是一个向量，则它的组成部分可以用下标\(X_j\)来取出。定量的输出变量用\(Y\)来表示，对于定性的输出变量采用\(G\)来表示。我们使用大写字母\(X\)、\(Y\)、\(G\)来表示变量，对变量的观测值我们用小写字母来表示；因此，\(X\)的第\(i\)个观测值记作\(x_i\) (可以是标量可以是向量)。举个例子，\(N\)个\(P\)维输入向量\(x_i\), \(i=1,2,3,…，N\), 可以表示为\(N*P\)的矩阵\(X\)： \[X=\left[\begin{array}{cccc}x_{11} & x_{12} & \cdots & x_{1 p} \\x_{21} & x_{22} & \cdots & x_{2 p} \\\vdots & \vdots & \ddots & \vdots \\x_{N 1} & x_{N 2} & \cdots & x_{N p}\end{array}\right]\] 通常我们把维度用列表示, 观测值用行表示。因此, 根据上述矩阵, 每 一行都是一个输入向量, 每一列则是一个维度。
最后给出一个简单定义, 如下: 给定输入向量\(X\) , 对输出做出一个估计, 记为\(\hat{Y}\). 并且如果\(Y\)取值为实数，则\(\hat{Y}\) ̂取值也是实数。需要数据去构建预测规则，因此我们假设有一系列可用的测量值(\(x_i\),\(y_i\)),\(i=1,…,N\),这也称为训练数据（training data)，利用这些训练数据去构建我们的预测规则。
\(X\)代表投入变量，如果\(X\)是个向量，组成部分为\(X_j\)。默认向量\(x_i\)是列向量，行向量为\(x_i^T\)。给定输入向量\(X^T=(X_1,X_2,…,X_P )\)，通过下面的等式预测\(Y\)值：
\[\hat{y}=\hat{\beta}_{0}+\sum_{j=1}^{p} X_{j} \hat{\beta}_{j}\]
\(\hatβ_0\) 是截距项，也叫机器学习里的偏差（bias）。
Linear regression
给定输入向量\(X^T=(X_1, X_2,…,X_P )\),想要预测实值输出变量\(Y\)，线性回归模型有如下形式：
\[f(X)=\hat{\beta}_{0}+\sum_{j=1}^{p} X_{j} \hat{\beta}_{j}\] 线性模型要么假设回归函数\(E(Y│X)\)是线性的，要么假设它本身是一个合理的近似。这里\(β_j\)是未知的参数或系数，变量\(X_j\)可以有下列不同的来源:
定量输入变量的变换，比如对数变换，平方根变换或者平方变换 ;
基函数展开，比如\(X_2=X_1^2\)，\(X_3=X_1^3\),一个多项式的表示;
虚拟变量形式，\(X_j\)其中一个系数为1，其他为0;
变量的交叉项，比如：\(X_3=X_2*X_1\) 。
无论\(X_j\)是哪个来源，模型关于参数都是线性的。（因为\(β_j\)是我们要估计的变量，与\(X\)无关） 一般地，假设有一个用来估计参数\(β\) 的训练数据集合(x)的训练数据集合\((x_1,y_i ),…(x_N,y_N )\)。每个\(x_i=(x_i1 ,x_i2,…,x_ip )^T\)是第i种情形的特征测量值的向量。
最小二乘估计
在给定经典线性回归的假设下，最小二乘估计量是具有最小方差的线性无偏估计量（Gauss-Markov定理）。线性回归模型可以写成： \[Y=Xβ+U, E(U)=0,Cov(U)=σ^2I\] \[Y=\left[\begin{array}{c}Y_{1} \\Y_{2} \\\vdots \\Y_{n}\end{array}\right], X=\left[\begin{array}{ccccc}1 & X_{11} & X_{12} & \cdots & X_{1 P} \\1 & X_{21} & X_{22} & \cdots & X_{2 p} \\\vdots & \vdots & \vdots & \ddots & \vdots \\1 & X_{n 1} & X_{n 2} & \cdots & X_{n p}\end{array}\right], \quad \beta=\left[\begin{array}{c}\beta_{0} \\\beta_{1} \\\vdots \\\beta_{p}\end{array}\right], U=\left[\begin{array}{c}u_{1} \\u_{2} \\\vdots \\u_{n}\end{array}\right]\] 其中，列为变量、行为样本点，经典线性回归的假定为：
（1） 线性性：模型对参数而言是线性的，\(Y=Xβ+U\)
（2） 外生性：\(E(U|X)=0\)
（3） 同方差：\(Var(U|X)=σ^2⋅I\)
（4） 独立性：\(rank(X)=p+1\)
取参数\(β=(β_0,β_1,…,β_p )^T\)使得下式的残差平方和最小。 \[\begin{array}{l} \operatorname{RSS}(\beta)=\sum_{i=1}^{N}\left(y_{i}-f\left(x_{i}\right)\right)^{2} \\=\sum_{i=1}^{N}\left(y_{i}-\beta_{0}-\sum_{j=1}^{p} x_{i j} \beta_{j}\right)^{2} \\=(Y-X \beta)^{T}(Y-X \beta) \\\frac{\partial \mathrm{RSS}}{\partial \beta}=-2 \mathrm{X}^{\mathrm{T}}(Y-\mathrm{X} \beta) \\\frac{\partial^{2} \mathrm{RSS}}{\partial \beta \partial \beta^{\mathrm{T}}}=2 \mathrm{X}^{\mathrm{T}} X \end{array} \] \[X^T (Y-Xβ)=0\]\[\hat{β} =(X^T X)^{-1} X^T Y\]
其中，\(y ̂_i=\hat{f}(x_i )\)，矩阵\(H=X(X^T X)^{-1}X^T\),时称之为”帽子”矩阵。估计值 可以看成y在子空间中的正交投影，帽子矩阵H计算了正交投影，因此也被成为投影矩阵。
OLS估计量的性质（最优线性无偏估计量，BLUE）
（1）\(\hat{\beta}\)是关于Y的线性估计量
（2）\(\hat{\beta}\)是无偏估计量
\(E(\hat{\beta}_{OLS})=E(E(\hat{\beta}_OLS│X))=E(E((X^T X)^{-1} X^T Y│X))=E((X^T X)^{-1} X^T E(Y|X))=E((X^T X)^{-1} X^T Xβ)=β\)
（3）有效性：最小方差 观测值\(y_i\)不相关，且有固定的方差\(σ^2\)，并且 是固定的（非随机的），最小二乘法的参数估计的最小二乘法的协方差矩阵 \[Var(\hatβ)=(X^T X)^{-1}σ^2\]
Guass-Markov 定理说明参数\(β\)的最小二乘估计在所有的线性无偏估计中有最小的方差。
均方误差准则
设未知参数是\(\theta\)，\(\hat\theta\)是它的某种估计\(\hat\theta\) 的均方误差定义如下： \[ MSE(\hat\theta)=E(\hat\theta-θ)^2=Var(\hat\theta)+[E(\hat\theta)-θ]^2 \]
第一项为方差，第二项为偏差平方和，均方误差要求两项都比较小。 假设\(X\)已中心标准化，且\(rank(X)=p\)，Gauss-Markov 定理表明最小二乘估计在所有无偏线性估计中有最小的方差，由于\(\hat\beta\)是\(\beta\)的无偏估计量，因此第二项方差项为0。 \[\hatβ=(X^T X)^{-1}X^TY \\=(X^T X)^{-1}X^T (Xβ+U) \\=β+(X^T X)^{-1} X^T U \] \[ MSE(\hat\beta)=Var(\hat\beta) \\=E(\hat\beta-E(\hat\beta)^2) \\=E(U^TX)((X^TX)^{-1})^{T}(X^TX)^{-1}X^TU \\=\sigma^2tr(X(X^TX)^{-1})^{T}((X^TX)^{-1}X^T) \\=\sigma^2tr(X^TX)^{-1} \]
\((X^TX)\)可逆，计\(λ_1，λ_2，…,λ_p>0\)是\((X^TX)\)的特征根。所以： \[MSE(\hat\beta)=\sigma^2\sum_{i=1}^pλ_i^{-1}\]
由上式可知，如果\((X^TX)\)存在一个非常小的特征根，那么\(MSE(\hat\beta)\)就会很大。尽管高斯-马尔科夫定理了保证了\(\sigma^2\sum_{i=1}^pλ_i^{-1}\)在线性无偏估计类中很小，但是它本身却很大，这种情况下，\(\hat\beta\)不再是\(\beta\)的一个良好估计。 估计量的期望预测误差为： \[ E(Y_0-\hat{f}(x_0)^2=\sigma^2+E(x_0^T \hat{\beta}-f(x_0 ))^2=\sigma^2+MSE(\hat{f}(x_0)) \]
在任意输入\(x_0\)处，预测值为\(\hat{f}(x_0)=x_0^T \beta\)，预测误差的估计值和均方误差只有常数值\(\sigma^2\)的差别， \(\sigma^2\)代表了新观测\(y_0\)的方差。
OLS存在的问题
\(\hatβ=(X^T X)^{-1}X^TY\)，该公式成立的条件是\((X^T X)\ne0\) ，也就是\((X^T X)\)能求逆，而当变量之间存在多重共线性或者特征数大于样本数时，\(X\)不是满秩矩阵，使得\((X^T X)\) 的结果趋近于0，回归系数估计的方差就很大，估计值不稳定。这也就是普通最小二乘法的局限性。
\((X^T X)\)存在一个非常小的特征根意味着什么？此时LS估计在均方误差准则小表现不好，并且\(X\)的列向量之间存在近似线性关系。
证明：记\(X_i，i=1，2，…，\) \(p\)是\(X\)的\(p\)列向量，若\(c\)对应\((X^TX)\)特征根\(\lambda\)的标准正交特征向量，且\(\lambda\approx0\)，那么： \[(X^T X)c=c\lambda\approx0\rightarrow c^T (X^T X)c=\lambda\approx0\rightarrow Xc\approx0 \\c_1 X_1+c_2 X_2+⋯c_p X_p\approx0\]
可以采用子集选择、降维、特征收缩-正则化来解决。本节课重点关注特征收缩-正则化的理解。（正则化就是通过增加一些限制项来缩减系数，达到降低模型复杂度的目的。）
Generalized Linear Models
线性模型存在以下局限：
第一，响应变量（response variable，\(Y_i\) ）不一定是服从正态分布的。数据可以是离散的，例如0/1变量，基数变量（count variable），序数变量（ordered variable），离散非序数变量（unordered categorical or nominal variable）等。
第二，\(Y\)的方差很有可能不是常数，会依赖于均值。一个例子就是0/1变量。很多时候我们会考虑使用伯努利分布（Bernoulli Distribution）来建模0/1变量， 假设\(Y_i~Bern(π_i )\)，其中\(π_i=P(Y_i=1│x_i )\)是真实参数，则有\(E(Y_i )=π_i\)，\(Var(Y_i )=π_i (1-π_i)\)，容易看出，方差会根据每一个数据的\(π_i\)不同而不同。
第三，存在值域问题。比如，对于Logistic回归模型，响应变量取0或1，对于OLS拟合的模型，响应变量预测值可能会位于0，1区间之外。
第四，数据是离散的情况下，线性回归是无法解释系数的。比如说对于一个分类变量（即数字仅仅代表所属的类别，而没有明确的数值含义），如果用正常的线性回归进行建模，事实上是不知道系数是什么意思的，因为在其他变量不变的情况下，从类别1变到类别2并不代表”增长”。 当响应变量为正态分布且有常数方差这一假设不适用时，对响应变量进行数据变换通常是解决响应变量非正态与方差不相等的非常有效的方法，加权最小二乘（需要知道权重）可能也是处理非常数方差问题的有效方法。更广泛使用的方法是广义线性模型（generalized linear model，GLM）。
Nelder and Wedderburn(1972)引入了广义线性模型（Generalized Linear Models，GLM），可以处理回归模型与实验设计模型、单位化常规正态理论的线性回归模型，单位化非线性模型。在广义线性模型的框架下，因变量不再要求连续、正态，能够对正态分布、二项分布、泊松分布、Gamma分布等随机因变量进行建模，OLS是一种特殊的GLM。
广义线性模型的参数估计一般不能用最小二乘估计，常用加权最小二乘法或最大似然法估计（MLE），假定响应变量的分布已知，各回归系数β需用迭代方法求解。参数估计也是基于大样本性质或者渐近理论。
Logistic回归模型
有二值响应变量的模型
首先考虑一种特殊情况：响应变量只有两个可能的取值，0和1，比如，要么成功要么失败。假设响应变量\(y_i\)为伯努利随机变量，概率分布为：
| \(y_i\) | 概率 |
|---|---|
| 1 | \(P(y_i=1)=\pi_i\) |
| 0 | \(P(y_i=0)=1-\pi_i\) |
由于\(E(\epsilon_i )\)，\(E(y_i)=π_i\)。一般情况下，当响应变量为二值变量时，会存在相当多的经验性证据表明响应变量函数的形状应当是非线性的，通常要使用S形的单调递增（或单调递减）函数。这一函数称为Logistic响应变量函数，形式为： \[E(y)=exp(x^{'}\beta)/(1+exp(x^{'}\beta)))\] 将这个函数线性化，将响应变量函数的均值函数的结构化部分定义出来，令： \[\eta=x^{'}{\beta}\] 为线性预测项，\(\eta\)通过以下变换定义： \[\eta=ln((π )/(1-π))\]
这一变换通常称为Logit变换，而变换中的比率\((π )/(1-π)\)称为优势比.有时将Logit变换称为对数优势比.
参数估计
Logistic回归模型的一般形式为： \[ y_i=E(y_i)+\epsilon_i \]\[ E(y)=\pi_i=exp(x^{'}\beta)/(1+exp(x^{'}\beta))) \]
用最大似然方法估计线性预测\(x_i^{'}\beta\)中参数。 每个样本观测值都服从伯努利分布，所以每个样本观测值的概率分布为： \[ f_i(y_i)=π_i^{y_i}(1-π_i)^{1-y_i} (i=1,2,⋯,n) \] 每个观测值yi的取值为0或1，由于观测值是独立的，所以似然函数就是: \[ L\left(y_1, y_2, \cdots, y_n, \boldsymbol{\beta}\right)=\prod_{i=1}^n f_i\left(y_i\right)=\prod_{i=1}^n \pi_i^{y_i}\left(1-\pi_i\right)^{1-y_i} \] 对数似然函数： \[ \ln L\left(y_1, y_2, \cdots, y_n, \beta\right)=\ln \prod_{i=1}^n f_i\left(y_i\right)=\sum_{i=1}^n\left[y_i \ln \left(\frac{\pi}{1-\pi_i}\right)\right]+\sum_{i=1}^n \ln \left(1-\pi_i\right) \]现在由于\(1-π_i=[1+exp(x_i^{'}\beta)]\)且\(η_i=ln[π_i/(1-π_i )]=x_i^{'}\beta\)，所以对似然函数可以写为 \[ lnL(y,\beta)=\sum_{i=1}^n y_i x^{'}_i\beta+\sum_{i=1}^n ln[1+exp( x^{'}_i\beta)] \]
通常在Logistic回归模型中，在\(x\) 变量的每个水平上都会有重复观测值或重复试验。这一点在试验设计时会经常出现，令 \(y_i\)表示观测到第i个观测值取 1 的个数，且\(η_i\) 为在每个观测值上的试验个数。然后对数似然函数变为: \[lnL(y,\beta)=\sum_{i=1}^n y_i ln(\pi_i)+\sum_{i=1}^n (n_i-y_i) ln(1-\pi_i)\] 令\(\hat\beta\)为模型参数的最终估计值，可以先使用迭代，然后再用加权最小二乘（IRLS）来求出MLE。如果模型假设正确，那么可以证明下式渐进成立：
\[E(\hat\beta)=\beta,且Var(\hat\beta)=(X^{'} VX)^{-1}\]式中：矩阵\(V\)为\(n\times n\)对角矩阵，其包含了主对角线上的每个观测值的方差的估计值，\(V\)的第i个对角元素为 \[V_ii=n_i π ̂_i (1-π ̂_i)\] 线性预测项的估计值为\(\hat\eta_i=x_i^{'} β\)，所以写出Logistic回归模型的拟合值为: \[\hat y_i=\hat\pi_i=exp(\hat\eta_i)/(1+exp(\hat\eta_i)))=exp(x^{'}\beta)/(1+exp(x^{'}\beta))) \]
GLM模型设定和基本假定
主要思想
GLM是一般线性模型的推广，它使因变量的总体均值通过一个非线性连接函数而依赖于线性预测值，允许响应概率分布为指数分布族中的任何一员。许多广泛应用的统计模型都属于广义线性模型，如常用于研究二元分类响应变量的Logistic回归、Poisson回归和负二项回归模型等。
GLM建模思想为：
假设因变量服从某个随机分布，如正态分布、二项分布
根据上述的假设分布构建因变量的转换形式
对转换后的随机变量进行线性拟合
基本假定
数据\(Y_1,Y_2,…，Y_n\)是独立分布的。
因变量\(Y_i\)不一定是正态分布，但是它通常假设是一个指数族的分布。比如二项分布、泊松分布、多项式分布等。
GLM不假设响应变量和解释变量之间存在线性关系，但它确实假设转换后的预期响应\(T(Y_i)\)与解释变量\(X\)之间存在线性关系。比如，对于二元逻辑回归\(logit(π)=β_0+β_1 x\)。
解释变量可以是一些原始变量的非线性变换。
方差的同质性不需要被满足。事实上，在许多情况下，考虑到模型结构，这甚至是不可能的。
误差是独立的，但不是正态分布的。
模型设定
在GLM中，回归形式可以如下： \[y_i=\beta_0+\beta_1 x_{1i}+⋯+\beta_p x_{pi}+\epsilon_i\] 响应变量\(y_i\)是由解释变量的线性函数加上一个误差项拟合。 具体来说，一个广义线性模型包含以下三个部分：
（1）随机成分（Random Component）：响应变量\(y_i\)的分布来自指数族分布的成员，这是模型中唯一的随机成分。指数族：大多数常用的统计分布是指数分布家族的成员，如正态分布、二项分布、泊松分布、指数分布、Gamma分布等，指数族分布的成员都有一般性的形式： \[f(y_i;\theta_i,\phi))=exp{(y_i \theta_i-b(\theta_i ))/(a(\phi))+h(y_i,\phi))}\] \(\phi\)是尺度化参数，\(\theta\)是自然位置参数。
（2）系统成分（Systematic Component）：模型中的解释变量\(（x_1,x_2,…，x_k)\)的线性预测算子（ linear predictor）： $\(\eta_i=\beta_0+\beta_1 x_{1i}+⋯+\beta_p x_{pi}\)$
（3）连接函数（Link Function）：\(\eta\)或者\(g(\mu)\)是单调可微的，指定了随机成分和系统成分之间的联系。它代表响应变量y预期值与解释变量的线性预测值之间的关系。 \[g[E(Y_i)]=g(\mu_i)=\mu_i\]
如果以概率论的方式解释回归（regression）这一过程，我们可以把通过给定的自变量\(x\)和相关的线性参数\(\theta\)估计因变量\(y\)的过程，理解为求解条件概率\(Pr(y|x,θ)\)的过程。也就是在给定\(η=θ^T x\)的条件下，求解因变量\(y\)的概率分布曲线。然后，计算这个概率分布的期望值\(E(y|x,θ)\)作为\(y\)的估计值，同时，这个概率分布的方差\(var(y|x,θ\))作为估计值的方差。
参数估计
GLM的均值和方差为： \[ \begin{aligned} & \mathrm{E}\left(y_i\right)=\mu_i=b^{\prime}\left(\theta_i\right) \\ & \operatorname{Var}\left(y_i\right)=a(\phi) b^{\prime \prime}\left(\theta_i\right) \end{aligned} \]
参数估计：用最大似然估计指数族GLM的参数。 对于样本\(y_1,…,y_n\)来说，似然估计为： \[L(\beta)=\prod_{i=1}^nf(y_i;\theta_i,\phi) \]
\[log L(\beta)=\sum_{i-1}^n\frac{y_i \theta_i-b(\theta_i)}{a(\phi)}+h(y_i, \phi)\]
Score function：对于\(\beta=\left(\beta_0, \beta_1, \cdots, \beta_K\right)\)
\[ S(\beta ; \phi)=\frac{\partial logL(\beta)}{\partial \beta}=\sum_{i=1}^n\left(\frac{\partial \mu_i(\beta)}{\partial \beta}\right)\left(\operatorname{Var}\left(y_i\right)\right)^{-1}\left(y_i-\mu_i(\boldsymbol{\beta})\right) \] 满足 \(g\left[\mu_i(\beta)\right]=\mathbf{x}_i^{\prime} \beta\)， 其中\(g[\mu_i(\beta)]=x_i^{'}\beta\) 令\(S(\beta;\phi)=0\)，得到\(\beta\)的信息矩阵（Information matrix）： \[ I(\beta;\phi)==\sum_{i=1}^n\frac{\partial \mu_i(\beta)}{\partial \beta}(\operatorname{Var}(y_i)^{-1})(\frac{\partial \mu_i(\beta)}{\partial \beta}) \]
Link function
满足单调可微的\(g(\mu)\)都可以是连接函数。
- Simple Linear Regression
SLR模拟连续响应变量的平均值如何取决于一组解释变量，其中\(i\)指的是每个观测值： \[\mu_i=\beta_0+\beta_1 x\]
①随机成分：\(Y\)从正态分布，均值为\(\mu\)，方差为\(\sigma^2\) ；②系统成分：x是解释变量（可以是连续的，可以是离散的），在参数\(\beta_0+\beta_1 x\)中是线性的，这可以扩展到多元线性回归。此外，解释变量本身也可以被转换，它们必须与参数系数以线性方式结合起来； ③连接函数：恒等连接（the identity link），\(η=g(E(Y_i ))=E(Y_i )\)是最简单的连接函数。
- Binary Logistic Regression
二元逻辑回归的模型是：二元响应变量的 “成功”几率如何取决于一组解释变量。 \[logit(π_i )=log (\frac{π_i}{1-π_i}) =\beta_0+\beta_1 x \]
①随机成分：响应变量的分布被假定为二项式，有一次试验和成功的概率，\(E(Y)=π\)； ②系统成分：\(x\)是解释变量（可以是连续的，可以是离散的），在参数中是线性的； ③连接函数：对数-赔率或对数连接(the log-odds or logit link)，\(\eta=g(π)=log (\frac{π_i}{1-π_i} )\)是最简单的连接函数。
- Poisson Regression
泊松回归模拟了离散(计数)响应变量的均值如何依赖于一组解释变量： \[log(\lambda_i)=\beta_0+\beta_1 x\]
①随机成分：\(Y\)服从泊松分布，均值为 \(\lambda\) ； ②系统成分：\(x\)是解释变量（可以是连续的，可以是离散的），在参数中是线性的； ③连接函数：对数连接(the log link)，\(\eta=g(\pi)。=log (\pi_i)\)
Canonical link
有一些连接函数叫规范连接（Canonical link）：对于指数族来说，如果\(\theta_i=g(\mu_i)\)，连接函数\(g(⋅)\)就叫规范连接。 在GLM中，响应变量服从指数族分布，所以有： \[g(\mu_i)=g(b^{'}(\theta_i ))=\beta_0+\beta_1 x_{1i}+⋯+\beta_p x_{pi}\] 定义规范连接： \[g=(b^{'})^{-1}\] \[g(\mu_i)=\theta_i=\beta_0+\beta_1 x_{1i}+⋯+\beta_p x_{pi}\] Canonical link可以简化参数估计： \[S(\beta;\phi)=\sum_{i=1}^nx_i(y_i-\mu_i(\beta))\]
常见的连接函数及对应分布
| 分布 | \(\theta\) | \(b(\theta)\) | \(\phi\) | \(E(y)=b^{'}(\theta)\) | \(Var(y)=b^{''}(\theta)\phi\) |
|---|---|---|---|---|---|
| 正态分布 | \(\mu\) | \(\theta^2/2\) | \(\sigma^2\) | \(\mu=\theta\) | \(\sigma^2\) |
| 逆高斯分布 | \(\frac{1}{\mu^2}\) | \(-(-2\theta)^{1/2}\) | \(\sigma^2\) | \(\mu=1/\sqrt(\theta)\) | \(\mu^3\sigma^2\) |
| 伽马分布 | \(1/\mu\) | \(-ln(\theta)\) | \(1/\gamma\) | \(\mu=1/\theta\) | \(\mu^3\gamma\) |
| 二项分布 | \(ln(p/(1-p)\) | \(ln(1+e^\theta)\) | 1 | \(p=\frac{e^\theta}{1+e^\theta}\) | \(p(1-p)\) |
| 泊松分布 | \(ln\lambda\) | \(e^\theta\) | 1 | \(\lambda=e^\theta\) | \(\lambda\) |
| 负二项分布 | \(ln\lambda\) | \(e^\theta\) | \(k\) | \(\lambda=e^\theta\) | \(\lambda+k\lambda^2\) |
R包glm简介
R提供了拟合广义线性模型的函数glm(),其调用格式为glm(formula, family = gaussian, data, weights, subset,na.action, start = NULL, etastart, mustart, offset,control = list(…), model = TRUE, method = “glm.fit”,x = FALSE, y = TRUE, contrasts = NULL, …)。 其中，formula为拟合公式，与函数lm()中的参数formula用法相同;最重要的参数是family,用于指定分布族，包括正态分布(gaussian)、二项分布(binomial)、泊松分布(poisson)和伪伽马分布(Gamma)，分布族还可以通过选项link来指定连接函数，默认值为family=gaussian(link=identity),而二项分布默认的连接函数是logit，即family=binomial(link=logit);data指定数据集;offset指定线性函数的常数部分，通常反映已知信息;control用于对待估参数的范围进行设置。
GLM的优点
与传统OLS回归相比，GLM的优点：
响应变量不需要转换为正态分布
链接的选择与随机成分的选择是分开的，使我们在建模时有更大的灵活性。
GLM模型通过最大似然估计进行拟合，因此似然函数和参数估计受益于渐进正态分布和齐次分布。
用于逻辑和泊松回归模型的所有参数估计和假设检验也适用于GLMs；例如，Wald和Likelihood比率检验、残差、置信区间等。
Shrinkage Methods
在线性回归模型中，其参数估计公式为\(\hat\beta=(X^T X)^{-1} X^TY\)，当(X^T X) 不可逆时，无法求出\(\beta\)，另外如果\(|X^TX |\)趋于0，会使得回归系数趋向于无穷大，此时得到的回归系数是无意义的。解决这类问题可以使用岭回归 (Ridge Regression)和套索回归(Lasso Regression))，主要针对自变量之间存在多重共线性或者自变量个数多于样本量的情况。
此外，多元线性回归模型中，如果所有特征一起考虑，容易造成过拟合使测试数据误差方差过大，最小二乘估计经常有小偏差大方差。预测精确性有时可以通过收缩或者令某些系数为0来提高，通过这些方法我们牺牲一点偏差来降低预测值的方差，因此可能提高整个预测的精确性。惩罚项越大意味着模型越简单，越来越多的特征系数被压缩到0，当惩罚项无限大的时候，只剩一个常数项，此时bias大variance小；惩罚项越小意味着模型越复杂，当惩罚项为0时，Lasso和OLS是一样的。
如何理解正则化？正则化可以降低模型复杂度，减少模型过拟合的风险，使模型权重之间的方差较小，在回归问题中的解释是：当特征之间存在高度相关关系的时候，假设有两个特征高度负相关，那么不带正则化的回归问题可能会赋予二者近似相等的很大权重，这样加权起来的结果仍然较小，但是由于权重很大，就导致了过拟合问题。
Ridge Regression
基本假定和模型参数推导
岭回归 (Ridge regression) 主要是针对过拟合问题，根据回归系数的大小加上惩罚因子对它们进行收缩．岭回归的系数使得带惩罚的残差平方和最小。岭回归的极值求解函数为（Penalty form）： \[ \hat{\beta}^{\text {Ridge }}=\underset{\beta}{\operatorname{argmin}}\left\{\sum_{i=1}^N\left(y_i-\beta_0-\sum_{j=1}^p x_{i j} \beta_j\right)^2+\lambda \sum_{j=1}^p \beta_j^2\right\} \]
\(\lambda \sum_{j=1}^p\beta_j^2\)为惩罚函数，保证了 值不会变得很大。这里\(\lambda\)≥0，是控制收缩程度的参数：λ值越大，收缩的程度越大。每个系数都向零收缩。 岭回归问题可以等价地写成（Constrained form）：（第二个条件：所有系数\(\beta\)平方和<某个阈值\(t\)，让\(\beta\)不会太大，\(\lambda\)和\(t\)是一一对应关系） \[ \hat\beta^{\text {Ridge }}=\underset{\beta}{\operatorname{argmin}}\sum_{i=1}^N(y_i-\beta_0-\sum_{j=1}^px_{ij}\beta_j)^2 \] subject to \[\sum_{j=1}^p\beta_j^2≤t\] 转为矩阵向量形式： \[J(B)=\sum(Y-X\beta)^2+\sum\lambda\beta^2\] 参数推导过程如下（闭式解）： \[ J(\beta)=(Y-X\beta)^T (Y-X\beta)+λ\beta^T β \] \[\frac{\partial J(\beta)}{\partial\beta}=0\] \[\hat\beta ^{Ridge}=(X^T X+\lambda I)^{-1} X^T Y\]
假设满足一个条件\(X^T X=I\)，即任意两个特征之间正交，那么有\(\hat\beta^{OLS}=X^T Y\)，那么\(\hat\beta ^{Ridge}=\hat\beta^{OLS}/(1+\lambda)\)。 L2范数惩罚项的加入使得\((X^T X+\lambda I)\) 满秩，保证了可逆，但是也由于惩罚项的加入，使得回归系数\(\beta\)的估计不再是无偏估计。所以岭回归是以放弃无偏性、降低精度为代价的回归方法。 如果\(\lambda>0\)，有：
岭回归使用了单位矩阵乘以常数\(\lambda\)，单位矩阵\(I\)的对角线上全是1，其余元素全是0，像一条山岭一样，这也是岭回归名称的由来。
岭回归的性质：
（1） 当岭参数为0，得到最小二乘解。
（2） 岭估计是最小二乘估计的线性变换。
（3） 当岭参数\(\lambda\)趋向更大时，岭回归系数\(\beta\)估计趋向于0。
（4） 岭回归中的\(\hat\beta(\lambda)\) 是回归参数β的有偏估计。它的结果是使得残差平和变大，但是会使系数检验变好。
（5） 对于回归系数向量来说，对任意\(\lambda>0\),\(||\hat\beta(\lambda)||\ne0\) 有偏估计回归系数向量长度<无偏估计回归系数向量长度， 即： \[||\hat\beta(\lambda)||<||\hat\beta||\]
（6）存在某一个\(\lambda\)，使得它所对应的MSE（估计向量的均方误差）<最小二乘法对应估计向量的MSE。即存在\(\lambda>0\)，使得,\(MSE(\hat\beta(\lambda))<MSE(\hat\beta)\)
岭回归的几何解释
以二维为例，残差平方和RSS可以表示为β1和β2的一个二次函数，数学上可以用一个抛物面表示，即图4.1所示。岭回归目标函数最小化问题可以等价为： 加入约束项，约束项对应着投影为β1和β2平面上的一个圆，即图4.2所示.
该圆柱与抛物面的交点对应的\(\beta1\)、\(\beta2\)值，即为满足约束项条件下的能取得的最小的\(\beta1\)和\(\beta2\)，从\(\beta1\beta2\)平面理解，即为抛物面等高线在水平面的投影和圆的交点，如下图所示：
岭回归的贝叶斯解释
正则化可以通过贝叶斯角度来进行解释。在贝叶斯学派中，任何事物都是有概率的，任何随机变量包括参数都要假设服从一定的分布。从贝叶斯的角度来看，正则化等价于对模型参数引入先验分布。Ridge回归可以解释为线性回归，其系数先验地服从正态分布（Normal prior distributions），而Lasso也可以被其系数服从拉普拉斯先验分布（ Laplace prior distributions）。
（1）回顾贝叶斯公式(Bayes’ Theorem)： \[P(A|B)=\frac{P(B|A)P(A)}{P(B)}\]
（2）理解先验分布和后验分布
贝叶斯理论假设回归系数向量β具有先验分布（Prior distribution） : \(P(\beta)=P(\beta|X)\)
给定参数\(\beta\) 以及自变量\(x_i\)后，数据\(Y\)的似然函数可以写成：\(f(Y│X,\beta)\)
将先验分布乘以似然，得到 “后验分布”（Posterior distribution）: \(P(\beta|X,Y)∝f(Y|X,\beta)P(\beta|X)\)
（3）区分最大似然估计MLE和最大后验概率估计MAP
对于\(P(x|\theta)\)，假设\(x\)表示具体的数据，\(\theta\)表示模型的参数。 如果\(\theta\)是已知确定的，\(x\)是变量，这个函数叫做概率函数(Probability function)，它描述对于不同的样本点\(x\)，其出现概率是多少。 如果\(x\)是已知确定的，\(\theta\)是变量，这个函数叫做似然函数(Likelihood function), 它描述对于不同的模型参数，出现\(x\) 这个样本点的概率是多少。 MLE是求参数θ,使似然函数P(x_0│θ) 最大，MAP则是想求\(\theta\)使\(P(x_0│\theta)P(\theta)\)最大，求得的\(\theta\)不单单让似然函数大，\(\theta\)自己出现的先验概率也得大。MAP是假设\(\theta\)服从先验分布（这是一个主观概率），然后利用已有的数据\(x\)更新\(\theta\)。
对于线性回归问题，假设给定参数\(\beta\)，并且服从分布\(p(\beta)\), 给定自变量\(x_i\)，观察到的值\(y_i\)服从高斯分布，即\(p(y_i│x_i,\beta)=N(\beta^T x_i,\delta)\)，那么回归问题就变成最大化后验概率（MAP）。
\[\beta^*=\underset{\beta}{\operatorname{argmin}}\ln(\prod_{i=1}^np(y_i\mid x_i,\beta)p(\beta)) \\=\underset{\beta}{\operatorname{argmin}}\sum_{i=1}^n(\ln p(y_i\mid x_i,\beta)+\ln p(\beta)) \]
证明过程： 引入参数\(\theta\)的先验分布，假设\(\theta^{(j)}\sim N（0,\delta^2）\) \[\begin{aligned}l(\theta) & =\ln L(\theta)=\ln \left(\prod_{i=1}^{N} p\left(y^{(i)} \mid x^{(i)} ; \theta\right) \prod_{j=1}^{M} p\left(\theta^{(j)}\right)\right) \\& =\ln \left(\prod_{j=1}^{N} e^{-\frac{\left.y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}} \prod_{i=1}^{M} \frac{1}{\sqrt{2 \pi} \delta} e^{-\frac{\left(\theta^{(j)}\right)^{2}}{2 \delta^{2}}}\right) \\& =\sum_{i=1} \ln \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}}+\sum_{j=1}^{M} \ln \frac{1}{\sqrt{2 \pi} \delta} e^{-\frac{\left(\theta^{(j)}\right)^{2}}{2 \delta^{2}}} \\& =N \ln \frac{1}{\sqrt{2 \pi} \sigma}+M \ln \frac{1}{\sqrt{2 \pi} \delta}-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{N}\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}-\frac{1}{2 \delta^{2}} \sum_{j=1}^{M}\left(\theta^{(j)}\right)^{2}\end{aligned} \]
等价于： \(\min _\theta J(\theta)=\sum_{i=1}^N\left(y^{(i)}-\theta^T x^{(i)}\right)^2+\eta\|\theta\|_2^2\)
上述的目标函数即Ridge损失函数。
参数的选择
随着模型复杂度的提升，在训练集上的效果就越好，即模型的偏差（预测值和真实值的差异）就越小；但是同时模型的方差（回归系数的方差）就越大。对于岭回归的岭参数\(\lambda\)而言，随着\(\lambda\)的增大， 就越大， 越小，模型的方差就越小；而\(\lambda\)越大使得\(\beta\)的估计值更加偏离真实值，模型的偏差就越大。所以岭回归的关键是找到一个合理的\(\lambda\)值来平衡模型的方差和偏差。
（1）岭迹法
当岭参数\(\lambda\)在\((\lambda,+\infty)\)内变化时，在平面坐标系上把\(\lambda-\beta(\lambda)\)- 画出来，画出来的曲线为岭迹。根据岭迹图我们可以选择合适\(\lambda\)值，称为岭迹法。
在\(\lambda\)很小时，通常各\(\beta\)系数取值较大；而如果\(\lambda\)=0，则跟普通意义的多元线性回归的最小二乘解完全一样；当\(\lambda\)略有增大，则各\(\beta\)系数取值迅速减小，即从不稳定趋于稳定。一般通过观察，选取喇叭口附近的λ值，此时各\(\beta\)值已趋于稳定，但总的RSS又不是很大。此外，通过岭迹图也可以观察变量是否有多重共线性，上图类似喇叭形状的岭迹图，一般都存在多重共线性。
岭迹图法与传统的基于残差的方法完全不同，它提供了一种全新的分析问题的思路，这是一种直观的易于理解的方法，这对于本文研究自变量之间的相互作用是有帮助的，可以说采用岭迹图法确定岭估计 \(k\) 值是定量分析与定性分析的一个很好的结合。但同时，岭迹图分析方法也存在着明显的缺点，比如采用岭迹图分析方法确定的岭参数 \(k\) 在一定程度上存在主观人为性，并且缺少可靠的理论依据。
(2)方差膨胀因子法
通过对于平均方差因子\(\widehat{VIF}\) 的计算，可以掌握变量之间存在多重共线性问题的程度。一般\(\widehat{VIF}>0\)，当 的时候，认为存在多重共线性问题。下面计算岭估计\(\hat\beta(\lambda)\)的协方差矩阵。 \[ \begin{aligned} & \mathrm{D}(\hat{\beta}(k))=\operatorname{cov}(\hat{\beta}(k), \hat{\beta}(k)) \\ & =\operatorname{cov}\left((X X+k I)^{-1} X^{\prime} y,(X X+k I)^{-1} X^{\prime} y\right) \\ & =\left(X^{\prime} X+k I\right)^{-1} X^{\prime} \operatorname{cov}(y, y) X\left(X^{\prime} X+k I\right)^{-1} \\ & =\sigma^2(X X+k I)^{-1} X X(X X+k I)^{-1} \\ & =\sigma^2 c(k) \end{aligned} \]\[ c(k)=(X X+k I)^{-1} X X(X X+k I)^{-1} \]
其对角线元素\(c(k)\)为岭估计的方差膨胀因子。从上面的分析中可以看出，当岭参数\(\lambda\)增大时，\(c(k)\)减小；当\(\lambda\)减小时\(c(k)\)增大。从大量的试验中得到经验，如果应用方差膨胀因子选择岭参数\(\lambda\)值，那么选择的\(\lambda\)是使得所有方差膨胀因子 \(c(k)≤ 10\)的 \(\lambda\) 值。
(3)最小化均方差预测误差（交叉验证法）
K折交叉验证，是说将样本数据随机分为K个等分。将第1个子样本作为”验证集”（validation set）而保留不用，而使用其余 K-1 个子样本作为”训练集”（training set）来估计此模型，再以此预测第 1 个子样本，并计算第1个子样本的”均方预测误差”（Mean Squared Prediction Error）。其次，将第2个子样本作为验证集，而使用其余K-1个子样本作为训练集来预测第2个子样本，并计算第2个子样本的MSPE。以此类推，将所有子样本的MSPE加总，即可得整个样本的 MSPE。最后，选择调整参数，使得整个样本的MSPE最小，故具有最佳的预测能力。这个方法最常用。
岭回归模型的优缺点
岭回归: 专用于共线性数据分析；当预测变量的数量超过观测变量的数量；不能用来筛选变量
岭回归模型的优点：可以提高预测精度，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的拟合要强于最小二乘法。对于有些矩阵，矩阵中某个元素的一个很小的变动，会引起最后计算结果误差很大，这种矩阵称为”病态矩阵”。） 岭回归模型的局限：但是因为λ不可能为无穷大，二次项求偏导时总会保留变量本身，不能使任何变量的系数等于零，很难满足减少变量个数的要求。
Lasso Regression
基本假定和模型参数推导
LASSO是个缩写，全称为The Least Absolute Shrinkage and Selection Operator ，1996年由Tibshrani发表于统计期刊Journal of the Royal Statistical Society。(Series B)，像岭回归一样是个收缩方法。
Lasso模型的基本假设：稀疏性假设。简单来说尽管世界如此复杂，但有用的信息却非常有限。套用到常见的统计学模型中，假如，我们考虑一个线性回归模型，有一个因变量\(Y\)，但有成百上千的自变量\(X\)。我们假设，只有有限个\(X\)的回归系数不为0，但其余的都是0。找到其中重要的\(X\)，对我们理解数据有重要的意义。
lasso 估计定义如下(拉格朗日形式)： \[\hat{\beta}^{\text {Ridge }}=\underset{\beta}{\operatorname{argmin}}\left\{\sum_{i=1}^N\left(y_i-\beta_0-\sum_{j=1}^p x_{i j} \beta_j\right)^2+\lambda \sum_{j=1}^p \mid \beta_j\mid\right\} . \]
Lasso回归问题可以等价地写成（Constrained form）： \[\hat\beta^{\text {Lasso }}=\underset{\beta}{\operatorname{argmin}}\sum_{i=1}^N(y_i-\beta_0-\sum_{j=1}^px_{ij}\beta_j)^2\] subject to \[\sum_{j=1}^p\mid \beta_j\mid≤t\]
两个式子是等价的，t越小，\(\lambda\)越大，对估计参数的压缩作用就越强。当我们对这个目标函数求最小时，一些不那么重要的自变量的系数将被压缩为0，从而达到筛选变量的作用。 参数推导(闭式解)： \[ J(\beta)=(Y-X\beta)^T (Y-X\beta)+\lambda\mid\beta\mid \\=ESS(\beta)+\lambda l_i(\beta) \] 首先求导得： \[2(X^T X\beta-X^T Y)+\lambda\frac{\partial\mid\beta\mid}{\partial\beta}=0\]
上面对一范数进行求导，要引入次梯度(sub gradient)，并且还要假设一个条件\(X^T X=I\),即任意两个特征之间正交，那么有\(β\beta^{OLS}=X^T Y\),\(\beta^{Lasso}=\beta^{OLS}-\frac{\lambda}{2}\frac{\partial\mid\beta\mid}{\partial\beta}\)，此时分情况考虑：
1）对于\(\beta_j^{OLS}<\frac{\lambda}{2}\),\(\beta_j=\beta^{OLS}-\frac{\lambda}{2}\frac{\partial\mid\beta\mid}{\partial\beta}\)，由于\(\frac{\partial\mid\beta\mid}{\partial\beta_j}<1\),所以\(\beta_j\>0\),此时，\(\beta_j=\beta_j^{OLS}-\frac{\lambda}{2}\)
2）对于\(\beta_j^{OLS}>-\frac{\lambda}{2}\),\(\beta_j=\beta^{OLS}-\frac{\lambda}{2}\frac{\partial\mid\beta\mid}{\partial\beta}\)，由于\(\frac{\partial\mid\beta\mid}{\partial\beta_j}\ge-1\),所以\(\beta_j<0\),此时，\(\beta_j=\beta_j^{OLS}+\frac{\lambda}{2}\)
3）对于\(-\frac{\lambda}{2}<\beta_j^{OLS}<\frac{\lambda}{2}\),用反证法，假设\(\beta_j>0\),\(\frac{\partial\mid\beta\mid}{\partial\beta_j}=1\),所以\(\beta_j=\beta_j^{OLS}-\frac{n\lambda}{2}\le0\)
综合，\(\beta_j^{Lasso}=sign(\beta_j^{OLS})（\mid \beta_j^{OLS} \mid-\frac{\lambda}{2})_+\),这是Lasso的闭式解，Lasso将OLS得到的系数的绝对值进行衰减\(\lambda/2\)，如果小于0了，就将其变为0，即稀疏化了系数。
Lasso的几何理解
以二维的情况为例，如下图所示，图中椭圆表示上式第一项在不同的取值时的图像椭圆上的点取值相同），越靠近椭圆中心越优。而图中以原点为中心的正方形则表示满足一式第二行的限制条件的点集，所以我们也只能选取落入该正方形的点。最终Lasso的估计值为椭圆和下面矩形的交点，除非椭圆与矩形正好相切在矩形的某条边上，否则交点将落在矩形的顶点上，这时某参数的估计值将被压缩到0，即该变量已被剔除出模型。
Lasso的贝叶斯理解
对于先验部分\(ln p(\beta)\)，假设参数\(\beta\)服从均值为0的拉普拉斯分布，那么可以得Lasso Regression： \[p_{l1} (\beta_j )=\frac{1}{2\alpha} exp(\frac{-\mid\beta_j \mid}{\alpha})\]
证明过程： 拉普拉斯(Laplace distribution)概率密度函数为： \[f(x\mid\mu,b)=\frac{1}{2b} exp(\frac{-(\mid x-\mu\mid)}{b})\]
引入参数\(\theta\)的先验分布，假设\(\theta^{(j)} \sim La（\lambda,0）\) \[ \begin{aligned}l(\theta) & =\ln L(\theta)=\ln \left(\prod_{i=1}^{N} p\left(y^{(i)} \mid x^{(i)} ; \theta\right) \prod_{j=1}^{M} p\left(\theta^{(j)}\right)\right) \\& =\ln \left(\prod_{i=1}^{N} \frac{1}{\sqrt{2 \pi} \sigma}^{-\frac{\left.y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}} \prod_{j=1}^{M} \frac{1}{2\lambda} e^{-\frac{\left|\theta^{(j)}\right|}{2 \delta^{2}}}\right) \\& =\sum_{i=1} \ln \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}}+\sum_{j=1}^{M} \ln \frac{1}{2\lambda} e^{-\frac{\left|\theta^{(j)}\right|}{\lambda}} \\& =N \ln \frac{1}{\sqrt{2 \pi} \sigma}+M \ln \frac{1}{2\lambda}-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{N}\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}-\frac{1}{\lambda} \sum_{j=1}^{M}\left|\theta^{(j)}\right|\end{aligned} \]
等价于： \(\min _\theta J(\theta)=\sum_{i=1}^N\left(y^{(i)}-\theta^T x^{(i)}\right)^2+\zeta\mid\mid \theta\mid\mid_1\)
上述的目标函数即Lasso损失函数。
参数的选择
通过下图了解一下Lasso筛选变量的动态过程，更直观地了解”压缩估计参数”的意思： 横轴为\(s=t/\sum_1^p |\hat\beta_j|\)，t越小表示限制越强，即\(s\)越小，图中垂直于横轴的虚线也更靠左，而图中其余曲线则代表在指定\(s\)的条件下，该变量系数的估计值。当我们确定\(s\)值后，各估计参数的值即为\(x=s\)的直线与图中其他曲线交点的横坐标。我们可以看到，当\(s\)越小，各估计参数相应的也被压缩得更小，而当s达到一定值以后，一部分不重要的变量将被压缩为0，代表该变量已被剔除出模型，图中从右至左不断下降的曲线如同被不断减小的\(s\)一步一步压缩，直到压缩为0。
在实际中，往往不会使用\(s\)，而是选择用\(\lambda\)作为压缩参数， \(\lambda\)越大，一般对应越严格的限制。有许多学者研究如何选取λ，一般采用最小化均方差预测误差（\(K\)折交叉验证）的方法，其他常见使用的准则有CV（Cross Validation）, AIC, BIC等。当存在真模型的假设下，常常选取BIC准则。
Lasso模型的优缺点
Lasso回归: 解决共线性问题；进行特征筛选； 当预测变量的数量超过观测变量的数量的时候。 Lasso模型的优缺点：Lasso相比于普通最小二乘估计，可以在变量众多的时候快速有效地提取出重要变量，产生稀疏估计，简化模型。但是， Lasso产生的是有偏估计。
Tibshirani 在他的原作中提出了两方面好处：
（1）虽然最小二乘估计是无偏估计，但是在变量过多的情况下，往往带有较大的方差。MSE（Mean Square Error）由偏差和方差两部分组成，Lasso虽然是有偏估计，但是在引入一定的偏差的同时，可能可以大量降低估计的方差，从而降低整体的MSE，如果我们拥有的样本信息是有限的，那么我们想要用有限的信息去估计过多的系数，此时信息很可能会出现不够用的情况，所以筛选变量提高估计效果十分必要。
（2）对于最后得到的回归方程，在估计出每一个放入模型的自变量的系数后，能够更好地解释它。当得到稀疏估计之后，能够明确得知哪些因素对Y有显著的影响。
岭回归和 lasso回归的对比
（1）在\(X\)为正交情况下的\(\beta_j\)的估计值 表格为在\(X\)为正交情况下\(（X^T X=I）\)的 的估计值，\(\lambda\)是通过对应的选择的常数。符号标记变量的符号（±1），而且 \(x+\)记\(x\)的正数部分。
| Ridge regression | \(\hat\beta ^{Ridge}=\hat\beta^{OLS}/(1+\lambda)\) |
| Lasso regression | \(\beta_j^{Lasso}=sign(\beta_j^{OLS})（\mid\beta_j^{OLS} \mid-\frac{\lambda}{2})_+\) |
在图形中，估计值由红色虚线来显示．灰色的 45∘直线作为参照显示了无约束的估计。
左图，岭回归做等比例的收缩。右图，lasso 通过常数因子\(\lambda\)变换每个系数，在 0 处截去。这也称作”软阈限”。
（2）在\(X\)为非正交情况下的\(\beta_j\)的估计值
当\(X\)为非正交的情形，当只有两个参数时，下图描绘了 lasso（左）和岭回归（右）的残差平方和的等高线。岭回归的约束区域为圆盘\(\beta_1^2+\beta_2^2≤t\)，lasso 的约束区域为菱形\(\mid\beta_1\mid+\mid\beta_2\mid≤t\)。两种方式都寻找当椭圆等高线达到约束区域的第一个点。与圆盘不同，菱形 (diamond) 有角；如果解出现在角上，则有一个参数\(\beta_j=0\)。当\(p>2\)，菱形变成了偏菱形 (rhomboid)，而且有许多角，平坦的边和面；对于参数估计有更多的可能为0。
Lasso限制各系数绝对值之和，而ridge regression则是限制各系数的平方和, 由于可行域从矩形变成的圆，其交点将为椭圆与圆的切点，且难以刚好落在坐标轴上，这也使得ridge regression很多时候并不能将多余变量的系数压缩为0。
Lasso损失函数等值线与正则化约束交点——使其中一个系数为0，此时实现了特征选择。Ridge既要满足损失函数尽可能小，又要满足约束条件，只能在切点处。因为若想让Ridge交在系数为0处，需要使损失函数的等高线向外扩张，此时不满足损失函数最小条件。
除了Lasso对应的L1惩罚以及ridge regression对应的L2惩罚，还有许多著名的惩罚函数。比如，在牺牲了惩罚函数凸性的情况下，一些non-concave的惩罚函数（如SCAD、MCP惩罚等）能够获得渐近无偏的估计，与之而来的是较高的计算成本。
（3）岭回归与Lasso回归，并没有绝对的好坏。
从预测的角度看，如果真实模型（或数据生成过程）确实是稀疏的，则 Lasso 一般表现更优。但如果真实模型并非稀疏，则岭回归可能预测效果优于 Lasso。 另一方面，从模型易于解释（interpretability）的角度，则 Lasso 显然是赢家，因为岭回归一般不具有变量筛选的功能。由于经济学通常强调模型易于解释，故 Lasso 在经济学中的应用更加广泛。
当只有一小部分预测变量是真实有效的而其他预测变量系数非常小或者等于0时， LASSO更好。当响应变量是很多(甚至全部的)预测变量的函数，Ridge 更好
Elastic Net
无论对于岭回归还是lasso回归，本质都是通过调节\(\lambda\)来实现模型误差vs方差的平衡调整。可以把岭回归和 lasso 一般化，并且可以看成是贝叶斯估计，\(q=1\)对应 lasso，\(q=2\)对应岭回归。 \[ \hat{\beta}^{\text {Ridge }}=\underset{\beta}{\operatorname{argmin}}\left\{\sum_{i=1}^N\left(y_i-\beta_0-\sum_{j=1}^p x_{i j} \beta_j\right)^2+\lambda \sum_{j=1}^p \mid \beta_j\mid^q\right\} . \]
\(q∈(1,2)\)表明在lasso和岭回归之间进行权衡，为了便于计算，Zou and Hastie (2005)引入弹性网惩罚（Elastic Net），它同时将 2-范数与 1-范数引入目标函数的惩罚项，在岭回归与套索估计的性质之间作了折衷，弹性网惩罚像 lasso 一样选择变量，同时像岭回归一样收缩相关变量的系数，原式变为： \[\hat{\beta}^{\text {Elastic net }}=\underset{\beta}{\operatorname{argmin}}\left\{\sum_{i=1}^N\left(y_i-\beta_0-\sum_{j=1}^p x_{i j} \beta_j\right)^2+\lambda \sum_{j=1}^p (\alpha\beta_j^2+(1-\alpha)\mid\beta_j\mid)\right\} ,q>0\]
在正交的假设下，参数估计结果为： \[\beta_j^{Elastic net}=sign(\beta_j^{OLS})（\frac{\mid \beta_j^{OLS} \mid-\frac{\lambda}{2}}{1+\lambda_2})_+\] \[\hat\beta ^{Ridge}=\hat\beta^{OLS}/(1+\lambda_1)\] \[\beta_j^{Lasso}=sign(\beta_j^{OLS})（\mid \beta_j^{OLS} \mid-\frac{\lambda_2}{2})_+\]
在图4.10中，圆圈黑实线为Elastic Net对应的约束，圆圈虚线为Ridge对应的约束，菱形虚线为Lasso对应的约束。
在图4.11中，黑实线为Elastic Net对应的系数估计值，黑虚线为Ridge对应的系数估计值，灰虚线为Lasso对应的系数估计值，可以比对图4.7
论文中的运用
具体运用方向
时间序列分析：LASSO和RIDGE也可以用于时间序列分析，例如对经济数据进行建模和预测。这些方法可以在数据集中存在大量变量和噪声的情况下提高模型的鲁棒性和准确性。
经济增长与收敛：LASSO回归和RIDGE回归被用来研究国家、地区或行业之间的经济增长和收敛现象，识别影响经济增长和收敛的因素。
贸易和投资流动：这两种方法被用来分析贸易和投资流动的决定因素，例如跨国公司的投资决策、国际经济政策制定等。
经济预测：这些方法可以用来预测宏观经济指标，例如通货膨胀率、失业率等，并且比传统的回归方法更准确。
金融风险管理：LASSO回归和RIDGE回归被用来识别金融市场中的风险因素，保险精算、股票市场波动预测等。
数据挖掘和机器学习： 这些方法在大量数据下进行特征选择和变量筛选，用于判断变量的重要性，提高模型准确性。
论文应用
Lasso 回归在经济学中的应用至少需要注意以下两方面的问题。
第一、作为收缩估计量，Lasso 是有偏的，而经济学家向来不喜欢有偏估计。解决方法之一为所谓 “Post Lasso” 估计量，即仅使用 Lasso 进行变量筛选，然后扔掉 Lasso 的回归系数，再对筛选出来的变量进行 OLS 回归。
第二、作为变量筛选算子（selection operator），Lasso 并不能保证避免 “遗漏变量偏差”（omitted variable bias）。比如，假设解释变量包含一个我们感兴趣的处理变量（treatment variable）以及诸多控制变量（control variables）。如果直接使用 Lasso 估计此方程，并进行控制变量的选择，则可能忽略对处理变量有影响的变量（由于这些变量可能与处理变量高度相关，故在回归方程包含处理变量的情况下，它们的作用可能被 Lasso 忽略），导致遗漏变量偏差。
为此，Belloni, Chernozhukov and Hansen (2014, REStudy) 提出了更为稳健的 “Post Double Lasso” 估计量，即将被解释变量与处理变量分别对所有控制变量进行 Lasso 回归，然后对这两个 Lasso 回归（即所谓 “Double Lasso”）所得的非零控制变量取并集（union）之后，再代入原方程进行 OLS 回归，以避免遗漏变量偏差。
另外，在存在很多工具变量的情况下，Belloni, Chernozhukov, Chen and Hansen (2012, Econometrica) 将 Lasso 应用于2SLS 的第一阶段回归，以获得最优的工具变量组合，这是 Lasso 在工具变量法方面的重要应用。
示例：《Gender Stereotyping in Academia: Evidence from Economics Job Market Rumors Forum》文章主要是通过文本挖掘的方法，利用Logistic LASSO、DID、Topic Analysis等模型，发现在EJMR网络论坛上，讨论一旦涉及到女性，其内容就会从专业探讨转变为私人话题。这是一篇颇具挑衅意味的论文，它探讨了一个活跃用户上万的匿名经济学家网站上的性别歧视现象。该篇文章在美国经济学界引起轰动，不仅被提名为2017年最有份量的经济学研究报告之一，还让美国经济学会还因此采取了政策干预。
作者爬取了EJMR论坛中自2014到2016年超过131913条帖子，共计1143416条的用户回复。作者设计了一个倾向得分模型来预测一个职位与一万个最常用单词的出现次数相关的性别，该模型有两个目的：第一，解决包含Level 1词语的重复帖子的情况；第二，找出对性别预测能力最强的单词。该模型的自变量为是否是女性，因变量是关于词频的词语矩阵。变量达10000左右，故采用LASSO的方法进行变量选择。 设\(Wi\)表示单词出现的频率，假设后验概率为: \[ P(Female_i=1\mid W_i)=\frac{exp(\theta_0+W_i^{'}\theta)}{1+exp(\theta_0+W_i^{'}\theta)} \] \[ P(Female_i=0\mid W_i)=\frac{1}{1+exp(\theta_0+W_i^{'}\theta)} \] 每个观察值的似然函数：
\[ P(Female_i=0\mid W_i)=P(Female_i=1\mid W_i)^{Female_i}\times P(Female_i=0\mid W_i)^{(1-Female_i)} \] 假设观测值是独立的，则N个观测值的对数似然函数为: \[ \begin{array}{c}l_{N}(\theta)=\log \left(\Pi_{i=1}^{N} P\left(\text { Female }_{i} \mid W_{i}\right)\right) \\=\sum_{i=1}^{N} \text { Female }_{i}\left(\theta_{0}+W_{i}^{\prime} \theta\right)-\log \left(1+\exp \left(\theta_{0}+W_{i}^{\prime} \theta\right)\right) \\\hat{\theta}_{\lambda}=\operatorname{argmin}_{\theta}\left(-l_{N}(\theta)\right)+\lambda\|\theta\|_{1} \\\text { where }\|\theta\|_{1}=\sum_{j \geq 1}\left|\theta^{j}\right|\end{array} \] Logistic-Lasso模型将对性别具有最强预测能力的词分出来。估计值\(\hat\theta_\lambda\) 是有偏差的，但模型的方差会减少，并倾向于产生更准确的预测结果。
作者取了 75% 的数据作为训练集，并使用5折交叉验证选择出了最优的调节参数。 剩下的数据用做训练集。将模型运用到 26,002 个重复项, 其中9,044个归类为Female = 1 其他的为Female = 0.6088个词在变量选择中系数被压缩为0，即这些词并未起到甄别性别的作用。
R的具体运用（Glmnet包）
glmnet功能强大，算法速度快。可以拟合多元线性回归模型、定制族广义线性回归模型和Ridge、lasso、弹性网回归模型，在该软件包中还拥有预测、绘图、交叉验证的方法。 glmnet算法采用循环坐标下降法，在其他参数不变的情况下，依次优化目标函数，并重复循环，直到收敛。该R包还利用强规则对活动集进行有效约束。 cv.glmnet使用k折交叉验证拟合模型，predict对数据进行预测（分类/回归），coef用于提取指定lambda时特征的系数。
加载包和数据
## Loading required package: glmnet## Loading required package: Matrix## Loaded glmnet 4.1-7回归方法的选择： alpha=0 岭回归 不进行变量的选择 alpha=1 lasso回归 进行变量的选择 其他[0,1]间取值则是Elastic net回归。
拟合模型
fit = glmnet(x,y,alpha = 1,family = "gaussian",standardize=TRUE)
# 这里alpha = 1，所以为lasso回归
##glmnet主要参数说明：
# alpha： elasticnet mixing parameter,0≤α≤ 1，
# 回归方法的选择： alpha=0 岭回归,alpha=1 lasso回归, alpha在[0,1]间取值则是Elastic net回归。
# nlambda: 默认100
# lambda: 用户提供的lambda序列。默认是让程序根据nlambda和lambda计算的lambda序列
# standardize：数据是否要标准化处理，默认为TRUE
# intercept，线性模型的截距，默认为TRUE
# family:模型中使用的误差分布和链接函数,取值范围:"gaussian", "binomial", "poisson", "multinomial","cox", "mgaussian",默认为gaussian
#fit = glmnet(x,y,alpha = 0,family = "gaussian",standardize=TRUE)
# Ridge回归
#fit = glmnet(x,y,alpha = 0.3,family = "gaussian",standardize=TRUE)
# elasticnet
## 自定义lambdas
# lambdas <- 10^seq(2, -2, by = -.1) # 100和0.01之间
# fit = glmnet(x,y,alpha = 1,family = "gaussian",
# lambda=lambdas,standardize=TRUE
class(fit)## [1] "elnet" "glmnet"##
## Call: glmnet(x = x, y = y, family = "gaussian", alpha = 1, standardize = TRUE)
##
## Df %Dev Lambda
## 1 0 0.00 1.63100
## 2 2 5.53 1.48600
## 3 2 14.59 1.35400
## 4 2 22.11 1.23400
## 5 2 28.36 1.12400
## 6 2 33.54 1.02400
## 7 4 39.04 0.93320
## 8 5 45.60 0.85030
## 9 5 51.54 0.77470
## 10 6 57.35 0.70590
## 11 6 62.55 0.64320
## 12 6 66.87 0.58610
## 13 6 70.46 0.53400
## 14 6 73.44 0.48660
## 15 7 76.21 0.44330
## 16 7 78.57 0.40400
## 17 7 80.53 0.36810
## 18 7 82.15 0.33540
## 19 7 83.50 0.30560
## 20 7 84.62 0.27840
## 21 7 85.55 0.25370
## 22 7 86.33 0.23120
## 23 8 87.06 0.21060
## 24 8 87.69 0.19190
## 25 8 88.21 0.17490
## 26 8 88.65 0.15930
## 27 8 89.01 0.14520
## 28 8 89.31 0.13230
## 29 8 89.56 0.12050
## 30 8 89.76 0.10980
## 31 9 89.94 0.10010
## 32 9 90.10 0.09117
## 33 9 90.23 0.08307
## 34 9 90.34 0.07569
## 35 10 90.43 0.06897
## 36 11 90.53 0.06284
## 37 11 90.62 0.05726
## 38 12 90.70 0.05217
## 39 15 90.78 0.04754
## 40 16 90.86 0.04331
## 41 16 90.93 0.03947
## 42 16 90.98 0.03596
## 43 17 91.03 0.03277
## 44 17 91.07 0.02985
## 45 18 91.11 0.02720
## 46 18 91.14 0.02479
## 47 19 91.17 0.02258
## 48 19 91.20 0.02058
## 49 19 91.22 0.01875
## 50 19 91.24 0.01708
## 51 19 91.25 0.01557
## 52 19 91.26 0.01418
## 53 19 91.27 0.01292
## 54 19 91.28 0.01178
## 55 19 91.29 0.01073
## 56 19 91.29 0.00978
## 57 19 91.30 0.00891
## 58 19 91.30 0.00812
## 59 19 91.31 0.00739
## 60 19 91.31 0.00674
## 61 19 91.31 0.00614
## 62 20 91.31 0.00559
## 63 20 91.31 0.00510
## 64 20 91.31 0.00464
## 65 20 91.32 0.00423
## 66 20 91.32 0.00386
## 67 20 91.32 0.00351输出结果解读:
Df: 自由度
%dev: 表示模型的可解释偏差，值越大表明该模型包括了越多样本的信息
Lambda: Lambda对应的值，可能小于指定的lambda数
查看模型参数
## [1] FALSE## 21 x 67 sparse Matrix of class "dgCMatrix"## [[ suppressing 67 column names 's0', 's1', 's2' ... ]]##
## (Intercept) 0.6607581 0.631235043 0.5874616 0.5475769 0.5112354 0.4781224
## V1 . 0.139264992 0.2698292 0.3887945 0.4971912 0.5959583
## V2 . . . . . .
## V3 . . . . . .
## V4 . . . . . .
## V5 . . . . . .
## V6 . . . . . .
## V7 . . . . . .
## V8 . . . . . .
## V9 . . . . . .
## V10 . . . . . .
## V11 . . . . . .
## V12 . . . . . .
## V13 . . . . . .
## V14 . -0.005878595 -0.1299063 -0.2429157 -0.3458857 -0.4397081
## V15 . . . . . .
## V16 . . . . . .
## V17 . . . . . .
## V18 . . . . . .
## V19 . . . . . .
## V20 . . . . . .
##
## (Intercept) 0.44579828 0.41340406 0.38515161 0.35554873 0.3268488
## V1 0.67960556 0.74364368 0.79672086 0.84797353 0.8959613
## V2 . . . . .
## V3 . . . 0.04902347 0.1156263
## V4 . . . . .
## V5 -0.03264358 -0.11207321 -0.18116746 -0.24689473 -0.3080361
## V6 . 0.01347456 0.06180401 0.10589979 0.1461096
## V7 . . . . .
## V8 . . . . .
## V9 . . . . .
## V10 . . . . .
## V11 . . . . .
## V12 . . . . .
## V13 . . . . .
## V14 -0.52122516 -0.58808817 -0.64543588 -0.69202504 -0.7319390
## V15 . . . . .
## V16 . . . . .
## V17 . . . . .
## V18 . . . . .
## V19 . . . . .
## V20 -0.03836926 -0.13486770 -0.22196207 -0.30436207 -0.3808057
##
## (Intercept) 0.3006986 0.2768715 0.2551611 0.23897850 0.22492518
## V1 0.9396872 0.9795286 1.0158306 1.05340431 1.08852453
## V2 . . . . .
## V3 0.1763105 0.2316038 0.2819849 0.32601558 0.36575238
## V4 . . . . .
## V5 -0.3637424 -0.4144999 -0.4607483 -0.50174373 -0.53883929
## V6 0.1827461 0.2161280 0.2465443 0.27480316 0.30064502
## V7 . . . . .
## V8 . . . 0.03113765 0.06549654
## V9 . . . . .
## V10 . . . . .
## V11 . . . . .
## V12 . . . . .
## V13 . . . . .
## V14 -0.7683077 -0.8014456 -0.8316396 -0.85716075 -0.88003892
## V15 . . . . .
## V16 . . . . .
## V17 . . . . .
## V18 . . . . .
## V19 . . . . .
## V20 -0.4504576 -0.5139218 -0.5717480 -0.61944973 -0.66194509
##
## (Intercept) 0.21212042 0.2004532 0.1898225 0.1801361 0.1713103 0.1632684
## V1 1.12051680 1.1496670 1.1762275 1.2004285 1.2224795 1.2425717
## V2 . . . . . .
## V3 0.40196705 0.4349645 0.4650306 0.4924256 0.5173870 0.5401292
## V4 . . . . . .
## V5 -0.57265298 -0.6034628 -0.6315355 -0.6571143 -0.6804208 -0.7016532
## V6 0.32419490 0.3456527 0.3652042 0.3830188 0.3992508 0.4140384
## V7 . . . . . .
## V8 0.09680196 0.1253263 0.1513166 0.1749980 0.1965756 0.2162341
## V9 . . . . . .
## V10 . . . . . .
## V11 . . . . . .
## V12 . . . . . .
## V13 . . . . . .
## V14 -0.90088249 -0.9198744 -0.9371791 -0.9529465 -0.9673131 -0.9804053
## V15 . . . . . .
## V16 . . . . . .
## V17 . . . . . .
## V18 . . . . . .
## V19 . . . . . .
## V20 -0.70066911 -0.7359530 -0.7681023 -0.7973956 -0.8240866 -0.8484063
##
## (Intercept) 0.16018043 0.1586314 0.15721186 0.1559184 0.15473991
## V1 1.25617136 1.2670788 1.27705027 1.2861359 1.29441444
## V2 . . . . .
## V3 0.56123386 0.5805823 0.59821036 0.6142724 0.62890756
## V4 . . . . .
## V5 -0.72008611 -0.7365962 -0.75164619 -0.7653592 -0.77785401
## V6 0.43008350 0.4454661 0.45947857 0.4722462 0.48387954
## V7 . . . . .
## V8 0.23282182 0.2474957 0.26088201 0.2730791 0.28419264
## V9 . . . . .
## V10 . . . . .
## V11 0.01749238 0.0386113 0.05783249 0.0753461 0.09130386
## V12 . . . . .
## V13 . . . . .
## V14 -0.99271796 -1.0040712 -1.01440902 -1.0238284 -1.03241106
## V15 . . . . .
## V16 . . . . .
## V17 . . . . .
## V18 . . . . .
## V19 . . . . .
## V20 -0.87452081 -0.8995116 -0.92227009 -0.9430068 -0.96190123
##
## (Intercept) 0.1536661 0.1526877 0.1517961 0.150933802 0.1501130
## V1 1.3019575 1.3088305 1.3150929 1.320563622 1.3253733
## V2 . . . . .
## V3 0.6422426 0.6543929 0.6654639 0.675052179 0.6833693
## V4 . . . . .
## V5 -0.7892388 -0.7996122 -0.8090641 -0.817360974 -0.8246020
## V6 0.4944794 0.5041376 0.5129379 0.521380737 0.5293940
## V7 . . . 0.004773491 0.0127738
## V8 0.2943189 0.3035455 0.3119525 0.319370009 0.3259468
## V9 . . . . .
## V10 . . . . .
## V11 0.1058440 0.1190924 0.1311638 0.142425610 0.1528706
## V12 . . . . .
## V13 . . . . .
## V14 -1.0402312 -1.0473567 -1.0538491 -1.059938983 -1.0656291
## V15 . . . . .
## V16 . . . . .
## V17 . . . . .
## V18 . . . . .
## V19 . . . . .
## V20 -0.9791172 -0.9948037 -1.0090967 -1.021794791 -1.0330999
##
## (Intercept) 0.14936467 0.14867414 0.147927056 0.145832036 0.14291659
## V1 1.32975267 1.33377821 1.337393911 1.340981414 1.34431995
## V2 . . . . .
## V3 0.69096092 0.69787701 0.704086481 0.708347140 0.71047211
## V4 . . . . .
## V5 -0.83122558 -0.83726751 -0.842853150 -0.848087765 -0.85182173
## V6 0.53669611 0.54334327 0.549403480 0.554823782 0.55960366
## V7 0.02005438 0.02668633 0.032703914 0.038519738 0.04475230
## V8 0.33193760 0.33741131 0.342430521 0.347221319 0.35117414
## V9 . . . . .
## V10 . . 0.001206608 0.010034050 0.01985732
## V11 0.16239419 0.17105029 0.178995989 0.186365265 0.19335515
## V12 . . . . .
## V13 . . . . .
## V14 -1.07081121 -1.07552680 -1.079993473 -1.086006902 -1.09350480
## V15 . . . -0.004444319 -0.01564233
## V16 . . . . .
## V17 . . . . .
## V18 . . . . .
## V19 . . . . .
## V20 -1.04340741 -1.05278699 -1.061382444 -1.069942845 -1.07874251
##
## (Intercept) 0.140076081 0.1374558396 0.1357322446 0.133073342 0.130676693
## V1 1.347383065 1.3500302965 1.3518343496 1.354221283 1.356361532
## V2 . 0.0003305368 0.0023648525 0.004540620 0.006504739
## V3 0.712555760 0.7149749797 0.7183628415 0.722517758 0.726148026
## V4 . . 0.0006067578 0.005717923 0.010273295
## V5 -0.855449422 -0.8591237651 -0.8627891179 -0.865619702 -0.868213547
## V6 0.564158207 0.5688257495 0.5734103636 0.576730059 0.579794269
## V7 0.050401187 0.0557935292 0.0616912891 0.067072276 0.072003466
## V8 0.354841505 0.3581403172 0.3606567338 0.363243052 0.365606825
## V9 . . . . .
## V10 0.028831014 0.0373622716 0.0457807585 0.053800086 0.061116762
## V11 0.199511972 0.2045194703 0.2080907599 0.210806709 0.213281313
## V12 . -0.0012547268 -0.0074697523 -0.013596900 -0.019177870
## V13 -0.001983046 -0.0072984582 -0.0111422626 -0.013991250 -0.016593423
## V14 -1.100192172 -1.1058992949 -1.1105322648 -1.115442340 -1.119919073
## V15 -0.025611555 -0.0348312100 -0.0451759986 -0.054376448 -0.062794808
## V16 . . . . .
## V17 . . . . .
## V18 . 0.0016051118 0.0074282515 0.011935155 0.016047237
## V19 . . . . .
## V20 -1.086745008 -1.0938819631 -1.0998730251 -1.105325880 -1.110285879
##
## (Intercept) 0.128660890 0.127100757 0.125658148 0.124459484 0.123083478
## V1 1.358328625 1.360077239 1.361826578 1.363516610 1.365012473
## V2 0.008183930 0.009356286 0.010663829 0.012117261 0.013394881
## V3 0.729395934 0.732747892 0.735687801 0.738489919 0.740982347
## V4 0.014458167 0.018153437 0.021807371 0.025450530 0.029051896
## V5 -0.870720166 -0.873442109 -0.875954408 -0.878483306 -0.880831285
## V6 0.582517165 0.585058350 0.587415384 0.589775898 0.592245889
## V7 0.076231973 0.079617935 0.083063167 0.086552751 0.089833929
## V8 0.367970084 0.370419245 0.372681284 0.374828053 0.377057526
## V9 . . . . -0.002621182
## V10 0.068030038 0.074233016 0.079723615 0.084204349 0.088799149
## V11 0.215812809 0.218374201 0.220998092 0.223955022 0.226508751
## V12 -0.024238268 -0.028646041 -0.032605455 -0.035873158 -0.038821242
## V13 -0.019278242 -0.022187336 -0.024764099 -0.027176713 -0.029215714
## V14 -1.124155296 -1.127787677 -1.131209037 -1.134226168 -1.136938041
## V15 -0.070162005 -0.076174415 -0.082038466 -0.087625136 -0.092981469
## V16 -0.002157957 -0.006135461 -0.009839126 -0.013415303 -0.016753758
## V17 . . -0.001922159 -0.006547775 -0.010856550
## V18 0.019755628 0.023150539 0.026194762 0.028945231 0.031379958
## V19 . . . . .
## V20 -1.114758360 -1.118794293 -1.122369457 -1.125537624 -1.127752032
##
## (Intercept) 0.121759609 0.120547496 0.11944328 0.11846462 0.11754879
## V1 1.366450256 1.367748542 1.36893056 1.36998470 1.37096264
## V2 0.014397675 0.015372031 0.01626067 0.01707702 0.01781339
## V3 0.743465148 0.745636509 0.74761351 0.74946962 0.75110911
## V4 0.032374342 0.035423370 0.03820166 0.04068249 0.04298924
## V5 -0.883160151 -0.885211320 -0.88707883 -0.88881016 -0.89035981
## V6 0.594632426 0.596777968 0.59873244 0.60050176 0.60212802
## V7 0.092794881 0.095510754 0.09798555 0.10021054 0.10226623
## V8 0.379289498 0.381282051 0.38309650 0.38471632 0.38622371
## V9 -0.005743726 -0.008552895 -0.01111123 -0.01336043 -0.01548785
## V10 0.092981099 0.096836620 0.10035051 0.10350507 0.10642372
## V11 0.228692663 0.230716422 0.23256140 0.23424727 0.23578173
## V12 -0.041612372 -0.044116969 -0.04639775 -0.04841234 -0.05030762
## V13 -0.031065122 -0.032739020 -0.03426420 -0.03568535 -0.03695091
## V14 -1.139223438 -1.141373902 -1.14333452 -1.14507381 -1.14670284
## V15 -0.097764379 -0.102178942 -0.10620215 -0.10981091 -0.11315326
## V16 -0.019873867 -0.022697069 -0.02526906 -0.02760956 -0.02974533
## V17 -0.014836101 -0.018454102 -0.02175038 -0.02474603 -0.02748296
## V18 0.033577799 0.035573518 0.03739207 0.03906579 0.04057628
## V19 . . . . .
## V20 -1.129497976 -1.131121615 -1.13260210 -1.13398890 -1.13521781
##
## (Intercept) 0.11671126 0.11594788 0.11525230 0.11461851 0.11404103
## V1 1.37185860 1.37267542 1.37341971 1.37409788 1.37471580
## V2 0.01848505 0.01909714 0.01965486 0.02016303 0.02062606
## V3 0.75260011 0.75395848 0.75519616 0.75632389 0.75735144
## V4 0.04509504 0.04701410 0.04876269 0.05035594 0.05180766
## V5 -0.89177002 -0.89305486 -0.89422555 -0.89529224 -0.89626417
## V6 0.60360748 0.60495526 0.60618329 0.60730223 0.60832176
## V7 0.10414080 0.10584892 0.10740530 0.10882341 0.11011555
## V8 0.38759898 0.38885221 0.38999412 0.39103458 0.39198261
## V9 -0.01742982 -0.01919950 -0.02081199 -0.02228122 -0.02361993
## V10 0.10908576 0.11151147 0.11372169 0.11573557 0.11757053
## V11 0.23717704 0.23844813 0.23960628 0.24066154 0.24162306
## V12 -0.05203815 -0.05361518 -0.05505213 -0.05636143 -0.05755441
## V13 -0.03810207 -0.03915083 -0.04010641 -0.04097710 -0.04177044
## V14 -1.14818987 -1.14954498 -1.15077972 -1.15190477 -1.15292987
## V15 -0.11620138 -0.11897887 -0.12150962 -0.12381555 -0.12591662
## V16 -0.03169119 -0.03346417 -0.03507964 -0.03655159 -0.03789278
## V17 -0.02997710 -0.03224969 -0.03432040 -0.03620715 -0.03792629
## V18 0.04195060 0.04320266 0.04434347 0.04538293 0.04633006
## V19 . . . . .
## V20 -1.13633467 -1.13735209 -1.13827911 -1.13912377 -1.13989340
##
## (Intercept) 0.11351485 0.11303541 0.11259857 0.11220053 0.1118690740
## V1 1.37527883 1.37579184 1.37625928 1.37668519 1.3770732640
## V2 0.02104795 0.02143237 0.02178263 0.02210178 0.0223925760
## V3 0.75828770 0.75914078 0.75991808 0.76062633 0.7612716571
## V4 0.05313040 0.05433564 0.05543381 0.05643442 0.0573461383
## V5 -0.89714976 -0.89795667 -0.89869190 -0.89936181 -0.8999722126
## V6 0.60925072 0.61009715 0.61086839 0.61157111 0.6122114061
## V7 0.11129289 0.11236564 0.11334309 0.11423371 0.1150452063
## V8 0.39284642 0.39363350 0.39435065 0.39500409 0.3955994813
## V9 -0.02483972 -0.02595114 -0.02696382 -0.02788655 -0.0287272973
## V10 0.11924248 0.12076591 0.12215399 0.12341876 0.1245711739
## V11 0.24249916 0.24329742 0.24402477 0.24468751 0.2452913700
## V12 -0.05864141 -0.05963185 -0.06053430 -0.06135657 -0.0621058001
## V13 -0.04249331 -0.04315195 -0.04375209 -0.04429891 -0.0447971478
## V14 -1.15386390 -1.15471496 -1.15549041 -1.15619698 -1.1568407715
## V15 -0.12783104 -0.12957539 -0.13116478 -0.13261297 -0.1339325016
## V16 -0.03911483 -0.04022831 -0.04124287 -0.04216730 -0.0430096118
## V17 -0.03949271 -0.04091996 -0.04222043 -0.04340536 -0.0444850328
## V18 0.04719304 0.04797936 0.04869582 0.04934863 0.0499434520
## V19 . . . . -0.0002826171
## V20 -1.14059465 -1.14123361 -1.14181581 -1.14234628 -1.1428571574
##
## (Intercept) 0.1116391211 0.111416020 0.111208836 0.111018972 0.110845721
## V1 1.3774624706 1.377778464 1.378068980 1.378335220 1.378578220
## V2 0.0226628898 0.022876591 0.023067319 0.023240539 0.023398270
## V3 0.7618229689 0.762333115 0.762792114 0.763209604 0.763589908
## V4 0.0581650485 0.058924124 0.059619334 0.060253956 0.060832496
## V5 -0.9005302148 -0.901019251 -0.901460720 -0.901862151 -0.902227796
## V6 0.6126935719 0.613199654 0.613661389 0.614081490 0.614464085
## V7 0.1158639519 0.116625929 0.117323876 0.117960550 0.118540773
## V8 0.3960464454 0.396486054 0.396890604 0.397260052 0.397596878
## V9 -0.0293237454 -0.029954538 -0.030538991 -0.031073136 -0.031560145
## V10 0.1255623472 0.126525647 0.127412702 0.128222375 0.128960349
## V11 0.2458135745 0.246332061 0.246801359 0.247227761 0.247615990
## V12 -0.0627355658 -0.063361864 -0.063941712 -0.064471794 -0.064955124
## V13 -0.0452205653 -0.045596759 -0.045935249 -0.046242852 -0.046522983
## V14 -1.1574445473 -1.158021349 -1.158552963 -1.159038292 -1.159480668
## V15 -0.1350226279 -0.136107335 -0.137103471 -0.138012175 -0.138840304
## V16 -0.0437592265 -0.044453331 -0.045085698 -0.045661882 -0.046186890
## V17 -0.0455062126 -0.046430870 -0.047272446 -0.048039238 -0.048737920
## V18 0.0505994498 0.051177394 0.051702567 0.052180547 0.052615915
## V19 -0.0008442286 -0.001340022 -0.001791685 -0.002203174 -0.002578088
## V20 -1.1433791483 -1.143843763 -1.144262012 -1.144641845 -1.144987654## 21 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 0.11476735
## V1 1.37393862
## V2 0.02004369
## V3 0.75605906
## V4 0.04998179
## V5 -0.89504175
## V6 0.60703946
## V7 0.10849039
## V8 0.39079024
## V9 -0.02193619
## V10 0.11526263
## V11 0.24041373
## V12 -0.05605396
## V13 -0.04077263
## V14 -1.15164057
## V15 -0.12327403
## V16 -0.03620592
## V17 -0.03576407
## V18 0.04513883
## V19 .
## V20 -1.13892541## 21 x 3 sparse Matrix of class "dgCMatrix"
## s1 s2 s3
## (Intercept) 0.6607581 0.46953687 0.150928072
## V1 . 0.61817569 1.320597195
## V2 . . .
## V3 . . 0.675110234
## V4 . . .
## V5 . -0.00867040 -0.817411518
## V6 . . 0.521436671
## V7 . . 0.004829335
## V8 . . 0.319415917
## V9 . . .
## V10 . . .
## V11 . . 0.142498519
## V12 . . .
## V13 . . .
## V14 . -0.46135970 -1.059978702
## V15 . . .
## V16 . . .
## V17 . . .
## V18 . . .
## V19 . . .
## V20 . -0.01019119 -1.021873704预测
模型拟合过程作图
交叉验证确定模型的最佳lambda（超参数）
##
## Call: cv.glmnet(x = x, y = y)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.07569 34 1.075 0.1483 9
## 1se 0.15933 26 1.210 0.1770 8
返回cv.glmnet对象，包含交叉验证的所有fits要素统计。
lambda.min:：最小的平均交叉验证误差
lambda.1se:：最正则化模型的λ值
对于图形来说，x轴代表经过log(lambda)，y轴代表模型的误差，cv.glmnet会自动选择使误差最小的lambda（左侧的虚线），最小的lambda值保存在cvfit$lambda.min中
Cox 回归
COX比例风险模型，主要用于带有时间的生存结局的影响因素研究，或评价某个临床治疗措施对患者生存的影响。
其他回归实例
## 当y仅包含两种分类时可以使用二分类模型，
data(BinomialExample)
x <- BinomialExample$x
y <- BinomialExample$y
fit = glmnet(x,y,alpha = 1,family = "binomial")
plot(fit,label = TRUE) 
#指定lambda在0.05和0.01时预测新样本的类别，type = "class"指定输出值为类别
predict(fit, newx = x[1:5,], type = "class", s = c(0.05, 0.01))## s1 s2
## [1,] "0" "0"
## [2,] "1" "1"
## [3,] "1" "1"
## [4,] "0" "0"
## [5,] "1" "1"#交叉验证
#cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class")
cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "auc")
plot(cvfit)
## 当y为多分类时
data(MultinomialExample)
x <- MultinomialExample$x
y <- MultinomialExample$y
#拟合模型
fit = glmnet(x, y, family = "multinomial", type.multinomial = "grouped")
predict(fit, newx = x[1:10,], type = "class")## s0 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17
## [1,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [2,] "3" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [3,] "3" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [4,] "3" "3" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "3" "3"
## [5,] "3" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [6,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [7,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [8,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "1" "1" "1" "1" "1" "1"
## [9,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "1" "1" "1" "1" "1"
## [10,] "3" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## s18 s19 s20 s21 s22 s23 s24 s25 s26 s27 s28 s29 s30 s31 s32 s33 s34 s35
## [1,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [2,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [3,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [4,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "1" "1" "1"
## [5,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "1" "1" "1" "1" "1" "1" "1" "1"
## [6,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [7,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [8,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [9,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [10,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## s36 s37 s38 s39 s40 s41 s42 s43 s44 s45 s46 s47 s48 s49 s50 s51 s52 s53
## [1,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [2,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [3,] "2" "2" "2" "2" "2" "2" "2" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [4,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [5,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [6,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [7,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [8,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [9,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [10,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## s54 s55 s56 s57 s58 s59 s60 s61 s62 s63 s64 s65 s66 s67 s68 s69 s70 s71
## [1,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [2,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
## [3,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [4,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [5,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [6,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [7,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
## [8,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [9,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [10,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"## Warning: executing %dopar% sequentially: no parallel backend registered#cvfit = cv.glmnet(x, y, family = "multinomial", type.measure = "class")
# Only deviance, class, mse, mae available as type.measure for Multinomial models;
#绘图展示
plot(fit, xvar = "lambda", label = TRUE, type.coef = "2norm")
## 1
## [1,] "3"
## [2,] "2"
## [3,] "2"
## [4,] "3"
## [5,] "1"
## [6,] "3"
## [7,] "3"
## [8,] "1"
## [9,] "1"
## [10,] "2"总结：
当构建回归模型时，y为连续变量，指定family = “gaussian/mgaussian”构建广义线性模型，可以进行如Lasso回归、Ridge回归、Elastic net等。在使用交叉验证cv.glmnet时损失函数选用type.measure = “mse”使用均方误差，使用predict预测新样本时指定type = “response”输出预测值。
当我们构建分类模型时，y为离散变量代表分类数据，指定family = “binomial/multinomial”以及type.multinomial = “grouped”构建逻辑回归模型，在使用交叉验证cv.glmnet时损失函数选用type.measure = “class/auc”使用错配误差或auc（二分类），使用predict预测新样本时指定type = “class”输出预测的类。
当我们构建Cox回归模型时，y为离散变量代表分类数据，指定family = “cox”构建Cox回归模型。
有关glmnet的更多介绍，可以参见https://cran.r-project.org/web/packages/glmnet/glmnet.pdf。
参考文献
(Hastie et al. 2009; Hoerl and Kennard 1970; McCullagh and Nelder 1989; Montgomery, Peck, and Vining 2021; Tibshirani 1996; Wu 2017)