偏最小二乘回归（Partial least square regression，PLSR）

发表于 更新于 分类于

如题。

PLSR，其概念可以追溯到PCA分析。主成分分析（PCA），是一种高维数据降维的方法，通过对数据进行正交主成分的查找，以寻找能够解释数据最大差异的特征组合。对于拥有成百上千维度的数据，使用PCA，一般可以降低到几个主成分（PCs），从而最终实现降维分析。

同样的，对于拥有成百上千维度的数据，如果要进行回归建模，那么高特征维度将会带来许多问题。一种解决方法是事先进行特征选择；另一种则是PCR，即主成分分析回归，先PCA找到各个主成分，然后对这些降维之后的主成分做回归分析，从而获得最终的结果。

或许可以这样简单地理解： PLSR=多元线性回归+PCA+相关性分析

虽然PLSR的数学原理与PCR不同，但实际解决的都是类似的问题，即如何对高维数据做回归。当数据量小，甚至比变量维数还小，而相关性又比较大时使用，这个方法甚至优于主成分回归。

模型介绍

设有 $m$ 个自变量 ${x_1,…,x_m}$ ， $p$ 个因变量 ${y_1,…,y_p}$ 。为了研究因变量和自变量的统计关系,我们观测了 $n$ 个样本点,由此构成了自变量与因变量的数据表 $X$ 和 $Y$ ：

$$
X=\left[\begin{array}{c}
x_{11} & x_{12} & \cdots & x_{1m} \\
x_{21} & x_{22} & \cdots & x_{2m} \\
\vdots & \vdots & \ddots & \vdots \\
x_{n1} & x_{n2} & \cdots & x_{nm} \\
\end{array}\right]
$$

$$
Y=\left[\begin{array}{c}
y_{11} & y_{12} & \cdots & y_{1p} \\
y_{21} & y_{22} & \cdots & y_{2p} \\
\vdots & \vdots & \ddots & \vdots \\
y_{n1} & y_{n2} & \cdots & y_{np} \\
\end{array}\right]
$$

其中，数据表的每一行代表一个样本点，每一列代表一个维度的变量。

偏最小二乘的一般多元底层模型是

$$
X = T P^T + E
$$

$$
Y = U Q^T + F
$$

其中

  • $X$ 是一个 $n\times m$ 的预测矩阵， $Y$ 是一个 $n\times p$ 的响应矩阵；
  • $T$ 和 $U$ 是 $n\times l$ 的矩阵，分别为 $X$ 的投影（X得分矩阵， X score）和 $Y$ 的投影（Y得分矩阵, Y score）；
  • $P$ 和 $Q$ 分别是 $m\times l$ 和 $p\times l$ 的正交载荷矩阵(“加载矩阵”，loading matrices)，
  • 矩阵 $E$ 和 $F$ 是错误项( error terms)，假设是独立同分布的随机正态变量。

偏最小二乘回归的目的是对 $X$ 和 $Y$ 分解来最大化 $T$ 和 $U$ 之间的协方差。具体来说：

  • 偏最小二乘回归首先分别在 $X$ 与 $Y$ 中提取出成分 $t_1$ 和 $u_1$ (也就是说, $t_1$ 是 $x_1,x_2,…,x_m$ 的线形组合, $u_1$ 是 $y_1,y_2,…,y_p$ 的线形组合，且满足 (1) $t_1$ 和 $u_1$ 应尽可能大地携带他们各自数据表中的变异信息; (2) $t_1$ 与 $u_1$ 的相关程度能够达到最大。这两个要求表明， $t_1$ 和 $u_1$ 应尽可能好的代表数据表 $X$ 和 $Y$ ,同时自变量的成分 $t_1$ 对因变量的成分 $u_1$ 又有最强的解释能力）。
  • 在第一个成分 $t_1$ 和 $u_1$ 被提取后，偏最小二乘回归分别实施 $X$ 对 $t_1$ 的回归以及 $Y$ 对 $u_1$ 的回归。
  • 如果回归方程已经达到满意的精度，则算法终止；否则,将利用 $X$ 被 $t_1$ 解释后的残余信息以及 $Y$ 被 $t_2$ 解释后的残余信息进行第二轮的成分提取。
  • 如此往复，直到能达到一个较满意的精度为止。
  • 若最终对 $X$ 共提取了 $l$ 个成分 $t_1,t_2,…,t_l$ （注意， $T={t_1,t_2,…,t_l}$ 就是前面我们提到的X得分矩阵；同理还有 $U={u_1,u_2,…,u_l}$ 是Y的得分矩阵），偏最小二乘回归将通过实施 $y_k$ 对 $t_1,t_2,…,t_l$ 的回归,然后再表达成 $y_k$ 关于原变量 $X_1,X_2,…,X_m$ 的回归方程, $k=1,2,…,l$ 。

循环迭代的三个向量

前文中的描述听起来有一些复杂。我们可以简化一下计算的逻辑，PLSR的数学步骤如下：

1
2
3
4
5
6
7
初始化：设置初始权重向量。
循环迭代：
    计算权重向量和得分向量。
    根据得分向量计算加载向量。
    更新权重向量。
停止条件：当达到预定的主成分数或满足某种收敛准则时停止。
模型评估：使用得到的模型进行预测和验证。

其中，“循环迭代”步骤是非常关键的部分，涉及到权重向量、得分向量和加载向量的计算。

下面我们详细解释一下权重向量、得分向量和加载向量这三个概念。

权重向量 (Weight Vector, w)

定义:

  • 权重向量是用于从原始自变量矩阵 $X$ 中提取得分向量 $t$ 的一组系数。
  • 在每次迭代中，权重向量指向 $X$ 的方向，该方向能够最大化 $X$ 和因变量 $Y$ 之间的协方差。

计算:

  • 初始权重向量通常是随机生成的。
  • 对于每次迭代 $i$，权重向量 $w_i$ 可以通过以下公式计算：

$$
w_i = \frac{X^T u_i}{|X^T u_i|}
$$

其中 $u_i$ 是残差向量 $e_i$（即 $Y - T B_Y$）的标准化版本。

得分向量 (Score Vector, t)

定义:

  • 得分向量是从原始自变量矩阵 $X$ 中提取出的低维表示。
  • 它们通常被看作是原始自变量的线性组合。

计算:

  • 得分向量 $t$ 是由 $X$ 与权重向量 $w$ 的点积获得的。

$$
t_i = X w_i
$$

加载向量 (Loading Vector, p)

定义:

  • 加载向量是得分向量 $t$ 和原始自变量 $X$ 之间的关系。
  • 加载向量指示了原始自变量与得分向量之间的相关性。

计算:

  • 加载向量 $p$ 可以通过以下公式计算：

$$
p_i = \frac{X^T t_i}{t_i^T t_i}
$$

三个向量的作用

  1. 权重向量 $w$ : 用于确定从原始数据中抽取得分向量的方向。指示了哪些变量对构建得分向量贡献最大。
  2. 得分向量 $t$ : 代表了从原始数据中提取的新变量。有助于减少数据维度，并且这些新变量与因变量 $Y$ 有较高的相关性。
  3. 加载向量 $p$ : 描述了得分向量与原始变量之间的关系。用于解释得分向量中每个变量的贡献程度。

基于scikit-learn库的PLSR计算

参考 sklearn - PLSRegression

在sklearn库中提供了PLSR的接口。使用方法很简单：

1
2
3
pls = PLSRegression(n_components=2)  # 创建 PLS 回归模型. n_components参数用于设定主成分数量
pls.fit(X_train, y_train) # 训练模型
y_pred = pls.predict(X_test) # 预测

一个代码实例

下面的代码块是一个实例。注意，sklearn的计算结果是一个对象，其包含了多个与PLS有关的属性，这些属性的理解如下：

  • x_scores_ or y_scores_ 是投影矩阵或者称之为得分矩阵（ $T$ or $U$ ），维度是 (n_samples,n_components) ，可以看作是降维以后的数据。
  • x_loadings_ or y_loadings_ 是加载矩阵（ $P$ or $Q$ ），维度是 (n_features,n_components) ， 在经过 $X = T P^T + E$ 和 $Y = U Q^T + F$ 的计算以后，得到的维度是 (n_samples,n_features) ，正好可以与原始矩阵一一对应。
  • x_rotations_ or y_rotations_ 是投影矩阵，维度为 (n_features,n_components) （与 loadings 矩阵一致），可以用于对X和Y的转换。（但不是简单的 X @ rotation -> scores ，似乎需要一些更复杂的操作）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
from sklearn.metrics import mean_squared_error
from sklearn.datasets import load_diabetes
from sklearn.cross_decomposition import PLSRegression
data=load_diabetes()
X=data.data  [0:100,0:10] # 训练数据的x，100个sample，10个feature，维度为(100,10)
X1 = np.zeros((X.shape[0],X.shape[1]*2))
X1[:,0:X.shape[1]] = X
X1[:,X.shape[1]:2*X.shape[1]]=X
X=X1 ## 此处做了一些改动，将数据集的各个特征复制了一份，以探索PLSR对于共线性特征的敏感性。
y=data.target[0:100]      # 训练数据的y，维度为(100,1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1) # 划分训练集和测试集
pls = PLSRegression(n_components=3)  # 创建 PLS 回归模型. n_components参数用于设定主成分数量
pls.fit(X_train, y_train) # 训练模型
y_pred = pls.predict(X_test) # 预测
mse = mean_squared_error(y_test, y_pred) # 计算均方误差
print("Mean Squared Error:", mse)
# 输出得分矩阵，即降维之后的X与Y
print("x_scores=",pls.x_scores_)
print("y_scores=",pls.y_scores_)
# 输出模型参数
print("x_loadings=",pls.x_loadings_)
print("y_loadings=",pls.y_loadings_)
print("x_rotation=",pls.x_rotations_)
print("y_rotation=",pls.y_rotations_)
print("x_weights_=",pls.x_weights_)
print("y_weights_=",pls.y_weights_)
print("coef_=",pls.coef_)
print("intercept_=",pls.intercept_)

上述代码块的输出如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
Mean Squared Error: 4441.943241667042
x_scores= [[ 1.98007095e-01  1.98610957e+00  9.12949852e-01]
 [ 2.21235201e-01  8.31216685e-01 -3.21104673e+00]
 ...
 [ 7.34663215e-01 -2.47641178e+00  7.53128615e-02]]
y_scores= [[-1.97439282e+00  2.51439028e+00  1.44606920e+00]
 [-2.79218866e+00  3.48781254e+00  7.27193215e+00]
 ...
 [ 8.48171117e+00 -8.96662807e+00 -1.77657480e+01]]
x_loadings= [[ 0.18945635  0.1960172   0.20756259]
 [ 0.15410153  0.35718437  0.25366944]
 ...
 [ 0.23590504  0.05698879  0.39514026]]
y_loadings= [[ 0.27374659 -0.23651343 -0.08640353]]
x_rotation= [[ 0.0475368   0.21231553  0.23576214]
 [-0.07582693  0.4342168   0.2386384 ]
 ...
 [ 0.14320436  0.07246786  0.43404028]]
y_rotation= [[ 1.97877605 -1.70963632 -0.62456754]]
x_weights_= [[ 0.0475368   0.24038099  0.18479015]
 [-0.07582693  0.38944901  0.13439312]
 ...
 [ 0.14320436  0.15701488  0.41664245]]
y_weights_= [[ 0.27374659 -0.23651343 -0.08640353]]
coef_= [[-3.6004284  -9.00994304  8.93538224  2.69884403  1.6118585  -4.9012268
  -1.02837883  0.9324192  16.95604055 -0.96559852 -3.6004284  -9.00994304
   8.93538224  2.69884403  1.6118585  -4.9012268  -1.02837883  0.9324192
  16.95604055 -0.96559852]]
intercept_= [130.8]

特征权重的探讨：如何确定在PLSR中每一个主成分都由哪些特征组成

在 PLSR中，x_loadings_（加载矩阵）是一个维度为 (n_features, n_components) 的矩阵，其中每一行对应一个特征，每一列对应一个主成分。然而加载矩阵中的元素含义并非特征权重，而是一个反映了原始特征与主成分之间的关系的数值。但我们可以使用下面的方法，通过对加载矩阵做一些操作，获取每一个主成分中各个特征的权重。

对于加载矩阵来说，矩阵中的每个元素 $p_{ij}$ 表示第 i 个特征与第 j 个主成分之间的相关性系数。如果 $p_{ij}$ 的值较大（无论是正值还是负值），则表示该特征对相应的主成分有较大的影响（加载矩阵中的元素可以是正数也可以是负数，这取决于特征与主成分之间的相关性方向）。我们可以使用下面的代码，提取每个主成分的特征组成以及这些特征的权重占比：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
# 假设我们已经有了一个训练好的 PLSRegression 模型 `pls_model` ，我们需要对 `pls_model`进行操作
# pls_model = PLSRegression(n_components=5).fit(X_train, y_train)
# 仍然以上面的那个数据集为例，因此这里令pls_model = pls
pls_model = pls
# 获取加载向量
x_loadings = pls_model.x_loadings_
# 计算每个主成分的特征权重占比
n_components = x_loadings.shape[1]
feature_weights = [] # feature_weights矩阵即为每个主成分的特征权重占比矩阵，维度为 (n_components,n_features)
for i in range(n_components):
    component_loadings = x_loadings[:, i]
    total_weight = np.abs(component_loadings).sum()
    weights = np.abs(component_loadings) / total_weight
    feature_weights.append(weights)

# 打印每个主成分的特征权重占比。请注意，feature1和feature11、feature2和feature12等的权重都是一致的，说明模型可以很好的捕获共线性特征。
for i, weights in enumerate(feature_weights):
    print(f"Component {i+1}:")
    count = 0 # 一个用于统计权重总和的变量。
    for j, weight in enumerate(weights):
        print(f"  Feature {j+1}: {weight:.4f}")
        count += weight
    print("total=",count)

仍以上面的那个数据集为例，我们可以看一看输出结果。事实上，前面在训练模型时，我用了一点点小心机，在代码块的第5-9行，将数据集的各个特征又复制了一份，以探索PLSR对于共线性特征的敏感性（也就是说，原先数据集有10个特征，经过处理现在有20个特征，多出来的10个特征和原先的10个特征是完全共线性的）。从下面的输出结果可以看出， feature1 和 feature11 、 feature2 和 feature12 等的权重都是一致的，说明模型可以很好的捕获共线性特征。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
Component 1:
  Feature 1: 0.0384
  Feature 2: 0.0312
  Feature 3: 0.0729
  Feature 4: 0.0523
  Feature 5: 0.0508
  Feature 6: 0.0366
  Feature 7: 0.0315
  Feature 8: 0.0544
  Feature 9: 0.0839
  Feature 10: 0.0478
  Feature 11: 0.0384
  Feature 12: 0.0312
  Feature 13: 0.0729
  Feature 14: 0.0523
  Feature 15: 0.0508
  Feature 16: 0.0366
  Feature 17: 0.0315
  Feature 18: 0.0544
  Feature 19: 0.0839
  Feature 20: 0.0478
total= 1.0
Component 2:
  Feature 1: 0.0493
  Feature 2: 0.0898
  Feature 3: 0.0210
  Feature 4: 0.0309
  Feature 5: 0.0683
  Feature 6: 0.1015
  Feature 7: 0.0282
  Feature 8: 0.0695
  Feature 9: 0.0272
  Feature 10: 0.0143
  Feature 11: 0.0493
  Feature 12: 0.0898
  Feature 13: 0.0210
  Feature 14: 0.0309
  Feature 15: 0.0683
  Feature 16: 0.1015
  Feature 17: 0.0282
  Feature 18: 0.0695
  Feature 19: 0.0272
  Feature 20: 0.0143
total= 1.0
Component 3:
  Feature 1: 0.0498
  Feature 2: 0.0609
  Feature 3: 0.0345
  Feature 4: 0.0196
  Feature 5: 0.0677
  Feature 6: 0.0770
  Feature 7: 0.0110
  Feature 8: 0.0617
  Feature 9: 0.0230
  Feature 10: 0.0949
  Feature 11: 0.0498
  Feature 12: 0.0609
  Feature 13: 0.0345
  Feature 14: 0.0196
  Feature 15: 0.0677
  Feature 16: 0.0770
  Feature 17: 0.0110
  Feature 18: 0.0617
  Feature 19: 0.0230
  Feature 20: 0.0949
total= 1.0

参考