一元回归分析

一元回归分析的模型

1、因变量$Y_i$与自变量$x_i$存在确定的线性关系：$Y_i=\beta_0+\beta_1x_i+\epsilon_i$ ，其中$x_i$是人为挑选的数值，不是随机变量。$\epsilon_i$是引入的误差项。

2、给定$x_i$时$Y_i$的预测值$\hat{Y_i}$满足如下关系：$\hat{Y_i}=\beta_0+\beta_1x_i$

3、$b_0,b_1$是从现有样本得来的$\beta_0,\beta_1$的估计值。$\hat{y_i}=b_0+b_1x_i$

4、$(x_i,y_i)$是手上的第$i$个样本

$\beta_0,\beta_1$的点估计

求偏导进行计算

为了估计$\beta_0,\beta_1$的大小，我们希望找到一组合适的$b_0,b_1$，使得样本中的每个点到预测直线的距离的平方和最小，即：

令$L=\sum_i^n{(y_i-b_0-b_1x_i)^2}$，再分别对$b_0,b_1$求偏导：

联立后可以解得：

最小二乘法计算

$b_0,b_1$的值还可以用最小二乘法进行计算。通过引入矩阵计算，可以让计算机快速完成参数的估计。此外，最小二乘法在多因素回归分析、甚至是曲线拟合中仍然适用。

假设$n>k$，对于所拥有的n个样本：

可以写成矩阵乘法的形式：

接下来构建损失函数，并且希望找到合适的$\pmb{b}$，使得残差的平方和最小，即$\pmb{y}-\pmb{\hat{y}}$的模长最小：

令$L(\pmb{b})$对$\pmb{b}$的偏导为0（这里是函数矩阵的导数，可以参看矩阵分析相关书籍）：

由此可以解得：

得到的结果应当和直接求偏导解出的结果一致。

点估计参数的性质

1、$b_0,b_1$是随机变量$y_i$的线性函数：

由点估计结果得：

上述两式说明$b_0,b_1$是$y_1,y_2,...,y_n$的线性组合。

2、$b_0,b_1$是参数$\beta_0,\beta_1$的无偏估计：

上述两式说明$b_0,b_1$是$\beta_0,\beta_1$的无偏估计。此外，易得$\hat{y_i}$也是$\hat{Y_i}$的无偏估计。

3、$b_0,b_1$的方差：

4、$b_0,b_1$的分布：

由上述内容可知，$b_0,b_1$均为$y_i$的线性组合，而每个$y_i$服从方差为$\sigma^2$的正态分布，故$b_0,b_1$也服从正态分布，而且我们已经知道了$b_0,b_1$的期望与方差，故：

5、$b_0,b_1$的协方差：

由于$y_i$之间相互独立，即$Cov(y_i,y_j)=0,i≠j$，故上式等于：

6、$\hat{y_i}$的分布：

因为$\hat{y_i}=b_0+b_1x_i$，而$b_0,b_1$都是$y_i$的线性组合，故$\hat{y_i}$也是$y_i$的线性组合，故其服从正态分布。

又因为：

故$\hat{y_i}\sim N(\beta_0+\beta_1x_i,[\frac{1}{n}+\frac{(\bar{x}-x_i)^2}{SS_{xx}}]\sigma^2)$

线性回归的拟合优度检验

由决定系数$R^2$判断拟合的优劣

至此，我们已经找到了使样本点到预测距离的平方和最小的拟合参数$b_0$和$b_1$，接下来需要判断这种回归模式是好是坏，即拟合的这条直线是否能够，或在多大程度上能够代表真实的$x$与$y$之间的关系。我们计算$y_i$的总变异（方差），并进行一系列变形：

其中第三项：

又由于$L$对$b_0$和$b_1$的偏导都为0：

将( 3 )和( 4 )带入( 2 )，发现( 2 )式正好是0：

故将( 2 )带回至( 1 )：

我们可以记：

• SST 总变异 = $\sum_{i=1}^n{(y_i-\bar{y})^2}$，代表$y_i$的总体变异。
• SSE 残差平方和 = $\sum_{i=1}^n{(y_i-\hat{y_i})^2}$, 代表不能被回归模型所解释的变异。
• SSR 回归平方和 = $\sum_{i=1}^n{(\hat{y_i}-\bar{y})^2}$，代表能够用回归模型解释的变异。

三者满足关系式：

当SST一定的时候，SSR接近SST，则说明回归模型的解释性越好，故不妨构造判定系数$r^2$：

Pearson相关系数

又称简单相关系数，定义为：

根据表达式发现，可以直接使用原始数据一步算出相关系数，而不用经过回归过程。此外，我们发现决定系数$R^2$满足：

即：在一元线性回归模型中，皮尔森相关系数的平方恰好与决定系数$R^2$相等。

由假设检验的方式判断是否存在显著线性关系

1、类似ANOVA的F检验，检验$\beta_1$是否为零：

H0：$\beta_1=0$，即在预测$y$的值的时候，$x$并未提供任何有用的信息。

H1：$\beta_1≠0$

在虚无假设的条件下，$y_i=\beta_0+\epsilon_i$，其中$\epsilon_i$彼此相互独立，且都服从$N(0,\sigma^2)$，故$y_i$也彼此相互独立，且都服从$N(\beta_0,\sigma^2)$。故$y_i$符合独立性、正态性、方差齐性。

这里不加证明地给出：$SSR/\sigma^2\sim \chi^2(1)$，$SSE/\sigma^2\sim \chi^2(n-2)$，则可以构建统计量F:

此时根据样本算出F，若$F＞F_{(1,n-2)}(0.05)$，即p<0.05，故认为$\beta_1≠0$，至少有一个$Y_i$的值与$x_i$有关。

2、t检验，检验$\beta_1$是否等于一个特定的值$\beta_1^*$

H0：$\beta_1=\beta_1^*$

H1：$\beta_1≠\beta_1^*$

因为$b_1\sim N(\beta_1^\star,\frac{\sigma^2}{SS_{xx}})$ ，故$\frac{b_1-b_1^\star}{\sqrt{\frac{\sigma^2}{SS_{xx}}}}\sim N(0,1)$，用s来代替$\sigma$得到：$\frac{b_1-b_1^*}{\sqrt{\frac{s^2}{SS_{xx}}}}\sim t(n-2)$

其中s为$\sigma$的估计值：$s^2=\frac{1}{n-2}\sum_{i=1}^n{e_i^2}=\frac{SSE}{d.f.}=MSE$。这里在算$e_i$的方差的时候分母是SSE的自由度（d.f = n-2），因为只有这样$s$才是$\sigma$的无偏估计量。

则构建统计量t:

若算出$|t|>t_{(n-2)}(0.05/2)$，则认为$\beta_1≠\beta_1^*$

3、Pearson相关系数$r$的假设检验：

H0：$r=0$

H1：$r≠0$

这里不加证明地给出：

若算出$|t|>t_{(n-2)}(0.05/2)$，则认为$r≠0$。

残差分析

定义残差$e_i=y_i-\hat{y_i}=y_i-b_0-b_1x_i$，则其具有如下性质：

又因为$e_i$也是$y_i$的线性组合，故$e_i\sim N(0,[1-\frac{1}{n}+\frac{(\bar{x}-x_i)^2}{SS_{xx}}]\sigma^2)$，服从正态分布。其中$|e_i|>2\sigma$的即标记为离群点，可以删除。

在实际应用中（不知道$\sigma$的真实值）可以使用学生化的残差：$SRE_i=\frac{e_i}{s\sqrt{1-\frac{1}{n}+\frac{(\bar{x}-x_i)^2}{SS_{xx}}}}$，若$|SRE_i>2(or\ 3)|$即视为离群点并剔除。

$\beta_0,\beta_1$的区间估计

在对$\beta_0,\beta_1$进行点估计后，还可以对其进行区间估计。即对每个$b_i$找到其置信下届$b_{i1}$和置信上届$b_{i2}$，使得：

因为$b_1\sim N(\beta_1,\frac{\sigma^2}{SS_{xx}})$，故：

故：

该式等价于：

故$\beta_1$的95%置信区间为$[b_1-t_{n-2}(0.05/2)\frac{s}{\sqrt{SS_{xx}}},b_1+t_{n-2}(0.05/2)\frac{s}{\sqrt{SS_{xx}}}]$

同理可得，$\beta_0$的95%置信区间为：$[b_0-t_{n-2}(0.05/2)\frac{s\sqrt{SS_{xx}+m\bar{x}^2}}{\sqrt{nSS_{xx}}},b_0+t_{n-2}(0.05/2)\frac{s\sqrt{SS_{xx}+m\bar{x}^2}}{\sqrt{nSS_{xx}}}]$

对于新的自变量$x$，对因变量$y$进行区间预测

因为$\hat{y}=b_0+b_1x$，故$\hat{y}\sim N(\beta_0+\beta_1x,[\frac{1}{n}+\frac{(\bar{x}-x)^2}{SS_{xx}}]\sigma^2)$

同时：

由于新的$y$与任何$y_i$都是独立的，故：

于是可以构建统计量t：

故$y$的95%置信区间为：$[y-t_{n-2}(0.05/2)\sqrt{[1+\frac{1}{n}+\frac{(\bar{x}-x)^2}{SS_{xx}}]}s,y+t_{n-2}(0.05/2)\sqrt{[1+\frac{1}{n}+\frac{(\bar{x}-x)^2}{SS_{xx}}]}s]$

即：$[b_0+b_1x-t_{n-2}(0.05/2)\sqrt{[1+\frac{1}{n}+\frac{(\bar{x}-x)^2}{SS_{xx}}]}MSE,b_0+b_1x+t_{n-2}(0.05/2)\sqrt{[1+\frac{1}{n}+\frac{(\bar{x}-x)^2}{SS_{xx}}]}MSE]$，该区间即为预测带。

回归分析实例

我们使用如下Python程序生成数据集，生成一百个随机的x值，并且令y=1.8x+2+e，其中e服从N(0,2)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

import numpy as np
import pandas as pd

# 设置随机数种子，以便结果可重复
np.random.seed(0)

# 生成X值，服从[3,9]的均匀分布
X = np.random.uniform(3, 9, 100)

# 生成Y值，Y=1.8X+2+e，其中e服从N(0,2)
e = np.random.normal(0, 2, 100)
Y = 1.8 * X + 2 + e

# 创建一个数据框
df = pd.DataFrame({'X': X, 'Y': Y})

# 将数据框保存为tsv文件
df.to_csv('output.tsv', sep='\t', index=False)

现在需要对这一百个点进行线性拟合，并判断拟合的优劣。

将数据导入origin软件

打开Origin2021软件，把数据粘贴进去：

进行线性回归

1、点击分析-拟合-线性拟合-打开对话框

2、输入数据选择B(Y)

3、输出量勾选上“置信区间的下限”与“置信区间的上限”

4、残差分析选择“学生化删除后”

5、拟合曲线图勾选上置信带与预测带，然后点击“确定”

解读程序报表

生成报表包含了拟合的一些相关信息：

1、拟合直线的截距$b_0$为2.50737，其95%CI为[1.10175,3.91299]，斜率为1.77898，其95%CI为[1.54808,2.00987]。

2、皮尔森相关系数$r$为0.83942，决定系数$R^2$为0.70462。一般来说，r在[0.8,1.0]为极强相关，故本例认为X和Y是极强相关。

解读拟合曲线图

图中较为重要的概念即置信带与预测带，其中置信带是指对于每一个x，观测到的y值的期望（均值）有95%的概率落在置信带中；预测带是指对于每一个x，观测到的y值有95%的概率落在预测带中。预测带一定比置信带要更宽些。

解读残差分析图

对于每个$y_i$，计算出残差$e_i=y_i-\hat{y_i}=y_i-b_0-b_1x_i$

根据回归分析的基本假设，$e_i$应当服从正态分布。这里我们可以检验$e_i$是否真的相互独立，且符合正态分布。如果成立，则说明结果的可解释性较好。左上图和左下图分别展示了残差与X、Y取值之间的关系，这里发现残差的值均匀的分布在0的两边，且不受X、Y的值影响；右上图为残差的频数分布直方图，应当类似正态分布的钟形曲线；右下图式残差的百分位图，应当与图中的蓝色曲线贴近。在本例中残差分析的结果良好，可以使用上述直线拟合。

多元回归分析

多元回归模型与参数的点估计

多元回归模型与一元回归模型类似，只不过是因变量$y$受到多个自变量$x$的调控：$Y_i=\beta_0+\beta_1x_{1i}+\beta_2x_{2i}+...+\beta_px_{pi}+\epsilon_i$

同样，$\beta_0,\beta_1,...,\beta_p$的估计值$b_0,b_1,...,b_p$可以用最小二乘法计算：

其中：

则有：

记$\pmb{H}=\pmb{x}(\pmb{x}^T\pmb{x})^{-1}\pmb{x}^T$，则有：

故$\pmb{H}$既对称矩阵，也是幂等矩阵，且$\pmb{H}$的迹为p+1。

残差分析

有了$\hat{\pmb{y}}$，还可以写出残差的表达式：

则 $ Cov(\pmb{e},\pmb{e})=(Cov(e_i,e_j))_{n \times n} $ ：

故 $ Cov(e_i,e_i)=Var(e_i)=\sigma^2(1-h_{ii}) $ ，其中 $ h_{ii} $ 是 $ \pmb{H} $ 主对角线上的元素。

此外，这里不加证明地给出，$s^2=\frac{1}{n-p-1}\sum_{i=1}^n{e_i}^2=\frac{1}{n-p-1}\pmb{e}^T\pmb{e}=\frac{1}{n-p-1}\pmb{y}^T(\pmb{I}-\pmb{H})\pmb{y}$是$\sigma^2$的一个无偏估计值。

点估计参数的性质

1、$b_i$是$y_i$的一个线性组合：

因为有最小二乘估计得：$\pmb{b}=(\pmb{x}^T\pmb{x})^{-1}\pmb{x}^T\pmb{y}$，即$\pmb{b}$是随机向量$\pmb{y}$的一个线性变换，故$b_i$是$y_i$的一个线性组合。这说明每个$b_i$都服从正态分布。

2、$\pmb{b}$是$\pmb{\beta}$的无偏估计：

3、$\pmb{b}$的协方差阵：

故$Var(b_i)$即为$\sigma^2(\pmb{x}^T\pmb{x})^{-1}$的主对角线上的元素。

拟合优度检验

由决定系数$R^2$判断拟合的优劣

类似的，决定系数$R^2$为：

将其开根号可以得到样本复回归系数$R$：

由假设检验的方式判断是否存在显著线性关系

1、方差分析

与一元回归类似。这里不加说明地给出方差分析表（ANOVA表格）：

H0：$\beta_0=\beta_1=...=\beta_p=0$

H1：至少有一个$\beta_i$不为0

2、对每个$\beta_j$进行检验：

$H_{0j}$：$\beta_j=0,j=1,2,...,p$

$H_{1j}$：$\beta_j \ne 0,j=1,2,...,p$

构建统计量：$t_j=\frac{b_j-\beta_j}{\sqrt{((\pmb{x}^T\pmb{x})^{-1}_{jj}s)}}\sim t(n-p-1)$

接下来对每个$\beta_j$进行假设检验，然后将不显著的$β_i$剔除后重新进行假设检验。

$\beta$的区间估计

与一元回归类似，$\beta_j$的95%CI为：

数据的标准化

数据标准化的第一步是将参数整体挪动到零点附近，这样拟合的结果就不会包含截距项，同时还能消除不同参数的量纲差异。第二步是让每个数据除以对应的标准差，使得标准化后的结果的方差为1。标准化的主要目的是使不同特征具有相同的尺度，从而提高模型的性能和收敛速度。这里的方差不用再除以n，因为每个标准化数据都要除以n，故除不除都无所谓。

则标准化后的回归方程为：

里面的系数同样可以使用最小二乘法轻易解出。

相关阵与偏相关系数

在一元回归模型中，简单相关系数（皮尔森相关系数，$r$）的平方与判定系数$R^2$相等，然而在多元回归模型中二者不再相等。其中判定系数是量化$y$与一组$x_i$之间的总的相关性，而简单相关系数量化了$y$与某一个特定的$x_i$之间的相关性。此外，简单相关系数还能反应不同自变量之间是否存在关联。

相关矩阵

构建相关矩阵$\pmb{r}=(r_{ij})_{p \times p}$，其中$r_{ij}$是第$i$个自变量和第$j$个自变量之间的简单相关系数：

其中$r_{ij}=\frac{\sum_{k=1}^n{(x_{ik}-\bar{x_i})(x_{jk}-\bar{x_j})}}{\sqrt{SS_{x_ix_i}SS_{x_j,x_j}}}$

如果$x_i$与$x_j(i \ne j)$是互相独立的，则$r_{ij}$应当是一个靠近0的数字，否则二者之间可能有一定的相关关系，对结果的判断可能会造成一定的影响。此外，我们发现相关矩阵还有一个简便的计算方式：

记标准化后的样本矩阵为$\pmb{x^\star}=(x_{ij}^\star)_{p \times p}=(\frac{x_{ij}-\bar{x_i}}{\sqrt{SS_{x_ix_i}}})_{p\times p}$，则发现：

即$(\pmb{x^\star})^T\pmb{x^\star}$恰好等于相关矩阵$\pmb{r}$。

偏决定系数与偏相关系数

偏决定系数与偏相关系数测量回归方程中已经包含若干个自变量时，再次引入一个新的自变量时，y的残差的相对减少量。它衡量的时某个自变量对y的残差的减少的边际贡献。

以$y$与$x_1$的偏决定系数为例：

其中$SSE(x_2,x_3,...,x_p)$即把样本中$x_{1i}$的部分去掉（假装没有这部分信息），然后算SSE。

对偏决定系数开根号就是偏相关系数。

小结

多元回归分析实例

数据准备

选取例3.1的数据，导入SPSS：

进行多元回归分析

1、点击分析-回归-线性，选取对应的自变量与因变量。方法选择后退法：

2、在【统计】中勾选：

估算值：输出回归系数的点估计结果

置信区间：输出回归系数的95%置信区间

协方差矩阵：输出系数的协方差矩阵$D(\pmb{\hat{\beta}})$

模型拟合：输出复决定系数$R^2$

描述：输出各变量的描述性统计结果，包括均值、标准差、个案数

部分相关性和偏相关性：输出因变量与各自变量之间的部分相关系数（part）和偏相关系数（partial）

3、在【保存】中勾选：

预测值-未标准化：添加未标准化的预测值

残差-未标准化：添加未标准化的残差

残差-学生化：添加学生化的残差

预测区间-平均值、单值

4、点击确定，分析输出结果

相关性表格可以看到变量间两两的简单相关系数以及显著性，其主对角线上的元素一定是1，此外，y和其它每个变量都应当有一定的相关性，而各自变量之间最好不要有较高的相关性（否则会导致多重共线性从而造成参数估计不准确）。

下表展示了使用后退法每步去除的变量，以及每次去除之后新模型的摘要。最终保留四个因变量（x2,x1,x3,x5），将其纳入模型。

ANOVA表反映了该模型的总体拟合优劣，对于最终的模型，P值小于0.001，即该模型有大于99.999%的概率反映了自变量与因变量之间的一定程度的线性关系。

系数表格记录了系数的估计值、标准误、标准化系数、t值、P值、系数的95%置信区间，y与该自变量的简单相关系数（Zero-order）,y与该自变量之间的偏相关系数（partial），y与该自变量之间的部分相关系数（part）。

系数相关性表格，记录了系数之间的协方差矩阵以及相关性矩阵。

残差分析

1、通过绘制散点图来判断残差是否均匀分布在0两侧，并且不与其它变量产生明显的相关性。

点击图形-旧对话框-散点图/点图-简单散点图，选取学生化残差作为y轴坐标，并随便选取一列作为x轴坐标：

这里可以发现学生化残差基本上均匀分布在0的两侧，符合基本假定。

2、检验原始的残差是否服从正态分布：

选取未标准化的残差绘制Q-Q图，检验其是否符合正态分布。选择分析-描述统计-Q-Q图：

发现点基本上落在直线上，故认为残差符合正态分布：

多重共线性的情形与处理

多重共线性对参数估计的影响

在多元回归模型的参数估计中，系数向量$\pmb{\hat{\beta}}$的协方差矩阵为：

当设计矩阵$X$的列向量之间存在完全多重共线性，则$|\pmb{X}^T\pmb{X}|=0$，无法通过最小二乘法求解出系数。而在实际问题中，当$|\pmb{X}^T\pmb{X}|≈0$时，系数估计值的方差就会很大，导致估计值与真实值相去甚远。

值得注意的是，多重共线性并不影响模型对样本的拟合效果，即使存在严重的多重共线性情况，模型也可能拥有极高的复决定系数。如果构建模型的目的就是为了根据新的数据进行预测，那么只要保证待预测的样本中共线性的关系保持不变，模型依然可以给出准确的预测值。然而，对于每一个系数来说其显著性差，方差很大，无法解释其实际意义。

多重共线性的判定方法

方差扩大因子 VIF

记$\pmb{C}=(\pmb{X^\star}^T\pmb{X^\star})^{-1}$，其中$\pmb{X^\star}$为中心标准化后的设计矩阵，则$\pmb{C}$的对角线元素$c_{jj}$即为$x_j$的方差扩大因子（$VIF_j$）。

可以证明：

$Var(\hat{\beta_j})=c_{jj}\sigma^2/SS_{x_jx_j}$

$c_{jj}=1/(1-R_j^2)$，其中$R_j$是自变量$x_j$与其余p-1个自变量的复决定系数。

则$VIF_j$的取值范围为$[1,+\infin)$。$VIF$值越大，则说明该自变量与其它自变量之间的线性相关程度就越强。一般经验认为当$VIF>10$时，该自变量与其它自变量之间存在严重的多重共线性问题。

在SPSS中，$Tol_j=1-R_j^2$为自变量$x_j$的容忍度，默认阈值为0.0001。当某个自变量与其余自变量之间的复决定系数＞0.9999，即容忍度＜0.0001时，SPSS会拒绝将该变量纳入模型。

特征根与条件数

由于矩阵的行列式等于其特征根的连乘，故当矩阵行列式接近0时，其特征根中至少有一个也近似为0。反之可以证明当矩阵至少有一个特征根接近0时，矩阵的列向量之间存在多重共线性。

记$\pmb{X}^T\pmb{X}$的最大的特征根为$\lambda_m$，则：

为特征根$\lambda_i$的条件数。

$k_i$体现了矩阵所有特征根之间的散布情况。条件数越大，说明该特征根相比于最大的特征根就越小，其值越接近0。

通常认为：

$0

$10≤k＜100$时，设计矩阵有较强的多重共线性；

$100≤k$时，设计矩阵存在严重的多重共线性。

经验判定法

1、增加或剔除一个变量，其它自变量的系数发生显著改变，则认为存在严重的多重共线性

2、主观上认为的某些重要的自变量在回归方程中没有通过显著性检验，则可能存在多重共线性

3、与因变量的简单相关系数很大，但是在模型中没有通过显著性检验，则可能存在多重共线性

4、回归系数与预期的差异很大，甚至正负号与预期不符，则可能存在多重共线性

5、自变量之间的简单相关系数很大，则可能存在多重共线性

6、某些重要自变量的回归系数的标准误很大，则可能存在多重共线性

消除多重共线性的方法

剔除变量法

可以将方差扩大因子最大者剔除，然后重新建立回归模型。如果仍存在者多重共线性的问题，则再剔除方差扩大因子最大者，直至回归方程中不存在严重的多重共线性问题为止。但是，倘若剔除了重要变量，则可能引起模型的设定误差。故在剔除变量时应当综合考虑系数的显著性检验、方差扩大因子以及生物学意义来判断。

增大样本量

当样本量增大时，$SS_{x_jx_j}$也同时增大，估计的参数的方差就会相应减小。在实验设计时尽可能使得n远大于p。

使用有偏估计

以有偏估计为代价，换取更小的方差。常用的有偏估计方法包括：

• 正则化：包括Lasso回归（L1正则化）、岭回归（L2正则化）

• 主成分回归

• 偏最小二乘回归

多重共线性诊断实例

在分析-回归-线性中勾选共线性诊断：

在输出的结果中可以看到拒绝阈、方差扩大因子和条件数。这里最大的条件数超过了100，说明该模型存在严重的多重共线性问题。

有偏估计

岭回归

岭回归的想法与推导

当设计矩阵呈现病态时，系数估计值的方差就会很大，与真实值可能相去甚远从而导致无法通过系数解读其生物学含义。岭回归（L2正则化）是解决多重共线性的一种方法，其以有偏估计为代价换来估计值更小的方差，从而导致估计值反而可能更加接近真实值。

与普通最小二乘法相比，岭回归即在目标函数的最后加上平方正则项，来约束每一个$\hat{\beta}$不至于太大：

故令：

可得：

该结果与普通最小二乘的结果多出了一项$\lambda\pmb{I}$，这可以看作：当$\pmb{X}^T\pmb{X}$呈病态、$|\pmb{X}^T\pmb{X}|≈0$时，给该矩阵人为加上一个正常数矩阵$\lambda\pmb{I}$，使得$|\pmb{X}^T\pmb{X}+\lambda\pmb{I}|$不再约等于0。则$\pmb{\hat{\beta}}(\lambda)$即为$\pmb{\beta}$的岭回归估计，$\lambda$为岭参数。实际上，$\pmb{\hat{\beta}}(\lambda)$是回归参数的一个估计族，估计值随着$\lambda$的变化而变化。当$\lambda$为0时，结果与普通最小二乘方法一致；当$\lambda$趋近于正无穷时，$\pmb{\beta}$趋近于$\pmb{0}$。

岭参数的选择

• 目测法

目测岭迹图，选取合适的岭参数，一般的原则为：

1、各回归系数的岭估计基本稳定

2、没有符号不合理的回归系数

3、没有不合理的绝对值

4、残差平方和不要增加太多

• 方差扩大因子法

与普通最小二乘类似，令$\pmb{c}(\lambda)=(\pmb{X}^T\pmb{X}+\lambda\pmb{I})^{-1}\pmb{X}^T\pmb{X}(\pmb{X}^T\pmb{X}+\lambda\pmb{I})^{-1}$，则其主对角线元素即为各因变量的方差扩大因子。此时选取合适的$\lambda$使得所有的方差扩大因子均≤10即可。

• 由残差平方和确定

由于使用岭估计时，在减少系数均方误差的同时增大了模型的残差平方和，故我们希望将残差平方和的增量限制在一定的范围之内。选取一个大于1的数c，要求$SSE(k)＜c\times SSE$

根据岭回归结果选择变量

首先将设计矩阵中心化并标准化，这样可以直接比较各个系数的大小。此时可以剔除掉标准化岭回归系数较为稳定且绝对值很小的变量。然后剔除岭回归系数不稳定，震动趋于0的变量。然而具体应该如何剔除并无一般规律可循，需要结合生物学意义来综合分析。

L1正则化与弹性网模型

L1正则化（Lasso回归）和弹性网模型的基本想法与岭回归一致，唯一的区别在于目标函数的正则项的次数。Lasso回归的正则项为$\lambda ||\pmb{\beta}||$，岭回归为$\lambda ||\pmb{\beta}||^2$，而弹性网模型为二者的线性组合：$\lambda ((1-\alpha)||\pmb{\beta}||^2+\alpha||\pmb{\beta}||)$

三者关系总结如下：

主成分回归

主成分回归的基本思想

在自变量较多的时候，通过正交旋转变换将原始的自变量转化为自变量的线性组合，并以其作为新的自变量。我们希望这些新的变量的其中几个能够最大程度的保留原始样本的信息，从而达到降维的效果。

主成分回归的算法简介

设有正交矩阵$\pmb{\mu}$，对$\pmb{X}$进行线性变换：

则：

其中$\pmb{\Sigma}$是随机向量$\pmb{X}$的协方差矩阵。我们将线性变换约束在以下原则之下：

• $\pmb{\mu}$是正交矩阵
• $Y_i$与$Y_j$不相关
• $Y_1$是$X_i$的所有满足第一个约束条件的线性组合中方差最大者；$Y_2$是与$Y_1$无关的、$X_i$的所有满足第一个约束条件的线性组合中方差最大者；以此类推。

则有如下结论：

设随机向量$\pmb{X}$的协方差矩阵为$\pmb{\Sigma}$，其特征根为$\lambda_1,\lambda_2,...,\lambda_p$，且$\lambda_1 \ge \lambda_2 \ge ...\ge \lambda_p$，$\pmb{\gamma_1},\pmb{\gamma_2},...,\pmb{\gamma_p}$是各个特征根对应的标准正交向量，则$Y_i=\pmb{\gamma_i}^T\pmb{X_i}$即为第$i$个主成分。

主成分的性质

1、$\pmb{Y}$的协方差阵$\pmb{\Lambda}$为对角阵，其主对角线上的元素为$\lambda_1,\lambda_2,...,\lambda_p$

2、记$\Sigma=(\sigma_{ij})_{p*p}$，则$\sum_i{\lambda_i}=\sum_i{\sigma_{ii}}$，即$Y_i$的方差之和等于$X_i$的方差之和。此外还有以下两个概念：

• 第k个主成分的方差贡献率$\alpha_k=\frac{\lambda_k}{\sum_i{\lambda_i}}$
• 前k个主成分的累计贡献率为$\frac{\sum_i^k{\lambda_i}}{\sum_i^p{\lambda_i}}$

3、因子负荷量$\rho(Y_k,X_i)=\mu_{ki}\sqrt{\lambda_k}/\sqrt{\sigma_{ii}}$即$Y_k$和$X_i$的简单相关系数，体现了第k个主成分中包含了多少的$X_i$成分。因子负荷量符合以下性质：

• $\sum_i{\rho^2(Y_k,X_i)\sigma_{ii}=\lambda_k}$
• $\sum_k{\rho^2(Y_k,X_i)=1}$

记$X_i$与前m个主成分的相关系数平方和为$Y_1,Y_2,...,Y_m$对$X_i$的方差贡献率，即$\nu_i=\frac{1}{\sigma_{ii}}\sum_k^m\lambda_k\mu_{ki}^2$，该值体现了前m个主成分中提取了原始变量中多少的信息。

偏最小二乘

偏最小二乘的基本思想

在主成分分析中，通过对设计矩阵X正交旋转变换，使得$X_i$的某些线性组合最大程度地保留了原始矩阵的信息。然而主成分分析完全从设计矩阵地角度选择合适的线性组合，然而，包含最多信息的线性组合并不都是与$Y$最相关的。偏最小二乘法解决了这一个矛盾，该算法在选取变量时综合考虑了设计矩阵以及观测值，主要选取那些与$Y$有相关性的变量进行最小二乘估计，故称为偏最小二乘。

偏最小二乘的算法

偏最小二乘算法可以概括为，选取$X_i$的某个线性组合$t$，将$Y$对$t$进行最小二乘估计。然后用$Y-\hat{Y}$算出残差，即在$Y$中将与$t$有关的那一部分减去。接着对每一个$X_i$，将$X_i$中与$t$有关的那一部分减去得到残存的$X_i$，然后将残差与残存$X_i$再进行同样的回归操作，直至$\pmb{X}^T\pmb{X}$的主对角线元素近似为0.

1、将$y$对每一个$x_i$进行回归。假设这里都是中心标准化后的数据：

2、将上述量进行加权求和得到$x_i$的一个线性组合，权重可以有很多种选择，为了方便起见一般选择$\omega_i=\pmb{x_i}^T\pmb{x_i}$：

则$t_1$的资料$\pmb{t_1}=\sum_i^p(\pmb{x_i}^T\pmb{y})\pmb{x_i}$

3、将$y$对$t_1$做回归分析：

得到残差$\pmb{y^{(1)}}=\pmb{y}-\pmb{\hat{y}}(t_1)$，则此时残差种便不再含有$t_1$的信息。

4、将每个$x_i$对$t_i$做最小二乘回归，并减去$t_i$的成分：

5、将$\pmb{y^{(1)}},\pmb{x_1^{(1)}},\pmb{x_2^{(1)}},...,\pmb{x_p^{(1)}}$作为新的原始资料并重复上述步骤，直到$\pmb{X}^T\pmb{X}$的主对角线元素近似为0。

6、最后将变量的线性组合还原回去即可得到偏最小二乘结果。

有偏估计实例

使用SPSS进行岭回归

进入SPSS的语法窗口，依次点击File-New-Syntax，在语法窗口录入以下命令：

1
2

INCLUDE'D:\SPSS\Samples\English\Ridge regression.sps'.
RIDGEREG DEP = y /ENTER    x1 x2 x3 x4 x5 x6 x7 x8 x9.

点击运行-全部，程序将生成不同的k值对应的判定系数、系数的岭估计值，以及岭迹图：

可以发现，参数的岭估计值在0.2左右变得相对稳定，并且判定系数下降的也不太多，此时可以调整k从0.1到0.3，步长为0.02进行岭参数估计：

1
2

INCLUDE'D:\SPSS\Samples\English\Ridge regression.sps'.
RIDGEREG DEP = y /ENTER x1 x2 x3 x4 x5 x6 x7 x8 x9 /START=0.1/STOP=0.3/INC=0.02.

得到如下结果：

最后综合考虑选取k=0.2：

1
2

INCLUDE'D:\SPSS\Samples\English\Ridge regression.sps'.
RIDGEREG DEP = y /ENTER x1 x2 x3 x4 x5 x6 x7 x8 x9 /k=0.2.

其中B是未标准化的参数，Beta是标准化后的参数。

使用SPSS进行主成分回归

点击分析-降维-因子，将所有的自变量选入变量框中，然后在【提取】中选择提取因子数为9（即与因变量的个数一致）；在【得分】中勾选‘保存至变量’：

点击【继续】得到输出结果：

一般而言，当累计贡献率＞80%的时候即可认为前k个主成分可以代表原始数据。这里选择前3或前4个主成分即可。将y与前四个主成分再做线性回归，得到如下结果：

此时可以看到方差扩大因子都变为了1，已经完全消除了多重共线性问题。该模型可以写成：

之后，再把每一个主成分还原为因变量的线性组合即可。

使用R进行偏最小二乘估计

使用SPSS进行偏最小二乘估计较为麻烦，需要安装相应的模块。这里使用R语言进行分析：

1
2
3
4
5
6
7

install.packages('pls') # 安装偏最小二乘包
data <- read.csv('./li3.1.csv') # 读取数据
data <- data.frame(scale(data)) # 中心标准化。偏最小二乘估计之前必须要中心标准化，之后还要手动将系数还原回去。

library('pls')
pls_result <- plsr(y~.,data=data,validation='LOO') # 建立偏最小二乘回归模型，LOO指使用留一交叉验证计RMSEP
summary(pls_result,what='all') # 输出回归结果

用调整的交叉验证法（adjCV）计算的RMSEP数值再分量取6的时候最小，可以选择$\alpha=6$，即进行六轮迭代。不过一般来说迭代3-4次基本已经足够：

1
2

pls_result_4 <- plsr(y~.,data=data,ncomp=4,validation='LOO',jackknife=TRUE)
coef(pls_result_4) # 输出回归系数

这个是标准化之后的系数估计值，还需要手动还原成原始的系数。

非线性回归

可化为线性回归的曲线回归

有些非线性模型可通过变量的代换，变为线性回归模型。例如假设$y=\beta_0+\beta_1e^{bx}+\epsilon$，只需要对样本中的每一个$x$，令$x'=e^{bx}$，再将$y$对$x'$做线性回归即可。

SPSS给出了一些常见的可线性化的曲线回归方程：

此外，多项式函数也是很常用的回归模型，并且可以轻易转为线性回归来处理。多项式模型可写作：

一般来说次数不超过三次，因为高次多项式模型一般很难解释其生物学含义。如果变量过多可以使用逐步引入变量的方式选择应当将哪些变量纳入模型中。

非线性最小二乘

分线性回归模型一般可记为：

则目标函数为：

即令：

即可求出$\pmb{\theta}$的值。一般用Newton迭代法求解。使用SPSS还可以给出近似的参数区间估计、显著性检验以及回归方程的显著性检验。在非线性回归中$SST=SSR+SSE$的关系式不再成立。定义非线性回归的相关指数$R^2=1-\frac{SSE}{SST}$