偏最小二乘回归(PLSR)

# 1、作用

偏最小二乘回归（Partial Least Squares Regression, PLSR）是一种统计方法，用于建模两个数据集之间的线性关系，当这两个数据集都包含多个变量且变量间存在多重共线性时尤其有用。PLSR通过提取两组变量（自变量集和因变量集）的潜在成分（latent components）来建立模型，这些潜在成分能够最大程度地解释自变量和因变量之间的协方差结构。这种方法不仅解决了自变量之间的多重共线性问题，还允许在自变量和因变量之间存在复杂关系时进行有效的预测。
主要用途：
预测建模：使用一组自变量（预测变量）来预测一组因变量（响应变量）。
变量选择：通过成分的重要性分析，可以帮助识别哪些自变量对因变量的预测最为重要。
数据结构分析：分析两组变量之间的复杂关系，特别是当直接的关系难以通过简单的线性回归模型描述时。

# 2、输入输出描述

输入：一组预测变量和一组响应变量。
输出：偏最小二乘回归结果，包括自变量累计投影重要性、成分矩阵表和模型系数结果等。

# 3、案例示例

案例：某康复俱乐部对20名中年人测量了三个生理指标：体重 x1, 腰围 x2, 脉搏 x3；三个训练指标：单杠 y1, 弯曲 y2, 跳高 y3。用偏最小二乘回归建立由三个生理指标分别预测三个训练指标的回归模型。

# 4、案例数据

# 5、案例操作

step1：新建项目；
Step2：上传数据；
Step3：选择对应数据打开后进行预览，确认无误后点击开始分析； step4：选择【偏最小二乘回归（PLSR）】；
step5：查看对应的数据数据格式，按【偏最小二乘回归（PLSR）】要求输入多个定量变量X和多个定量变量Y；
step6：选择是否自动确定最大主成分数量；
step7：点击【开始分析】，完成全部操作。

# 6、输出结果分析

输出结果1：因子方差解释情况表

潜在因子	X方差	累计的X方差	Y方差	累计的Y方差（R²）	调整后的R²
1	0.641	0.641	0.209	0.209	0.166
2	0.139	0.78	0.029	0.239	0.149
3	0.075	0.855	0.038	0.277	0.141

图表说明：上表展示潜在因子的信息综合解释能力。其中，累计的X方差代表对自变量信息的提取，累计的Y²（R²）代表对因变量信息的提取，可以以此为依据确定参数最大主成分数量。
分析：因子对方差解释情况表的结果显示，前【3】个潜在因子就可解释自变量80%的信息，全部潜在因子也不能解释因变量80%的信息。

输出结果2：自变量VIP（累积投影重要性）汇总表

变量	因子1	因子2	因子3
体重	1.022	0.998	1.019
腰围	1.336	1.298	1.22
脉搏	0.414	0.565	0.689

图表说明：上表展示VIP（累积投影重要性）的情况，它表示成分个数不同时，X对于Y的解释重要性力度，也可以用于参考最大主成分数量。其中对于VIP很大（大于1）的自变量，它在解释潜在因子（从而在解释因变量）时作用相对更大一些。

输出结果3：自变量VIP（累积投影重要性）图

图表说明：上图将VIP（累积投影重要性）的情况可视化的展示出来。

输出结果4：成分矩阵表

变量	因子1	因子2	因子3
体重	0.59	-0.368	0.935
腰围	0.771	0.7	-0.802
脉搏	-0.239	0.636	0.223
弯曲	-7.602	10.434	1.913
单杠	6.778	-11.376	-1.071
跳高	-1.068	-3.185	-3.01

图表说明：由上表可得到主成分分析降维后的成分矩阵表。

输出结果5：因子载荷系数表

变量	因子1	因子2	因子3
体重	0.666	0.02	0.657
腰围	0.676	0.355	-0.287
脉搏	-0.359	1.194	0.697
弯曲	-0.416	-0.291	0.455
单杠	-0.342	-0.336	0.477
跳高	-0.143	-0.065	-0.213

图表说明：上表为因子载荷系数表，可以分析到每个因子中显变量的重要性。

输出结果6：模型系数结果表

	弯曲	单杠	跳高
常数	145.55	9.45	70.3
体重	17.966	1.947	-13.28
腰围	-55.673	-4.662	0.749
脉搏	1.005	-0.137	-2.802
R²	0.34	0.436	0.054

图表说明：展示了本次PLSR模型结果，主要包括模型的系数，用于分析自变量X对于因变量Y的影响关系情况。
分析：模型的标准化公式为：
单杠 = 9.45+1.947 * 体重(标准化)-4.662 * 腰围(标准化)-0.137 * 脉搏(标准化)
弯曲 = 145.55+17.966 * 体重(标准化)-55.673 * 腰围(标准化)+1.005 * 脉搏(标准化)
跳高 = 70.3-13.28 * 体重(标准化)+0.749 * 腰围(标准化)-2.802 * 脉搏(标准化)
注：变量必须进行标准化才允许使用该公式计算

输出结果7：拟合图

# 7、注意事项

• 预测变量可以是连续的定量变量或定类变量
• 响应变量应当是连续变量
• 预测变量无需固定

# 8、模型理论

# 1.简介

偏最小二乘回归 (PLSR) 是一种多因变量对多自变量的回归建模方法，是最小二乘方法的推广。用于解决两组多重相关变量间的相互依赖关系，并研究用一组变量（自变量或预测变量）去预测另一组变量（因变量或响应变量）。
当两组变量的个数很多，且线性相关，而观测数据的数量又较少时，适合用偏最小二乘回归建立的模型。
偏最小二乘回归有如下优点：

1. 提供了一种多因变量对多自变量的回归建模方法；
2. 有效地解决变量之间的多重共线性问题；
3. 适合当样本点数量少于自变量个数时进行回归建模；
4. 最终模型中含原有所有自变量，回归系数容易解释。
5. 计算简单、预测精度高，易于定性解释。

# 2.原理

典型相关分析（CCA）对特征的处理方式比较粗糙，用的是线性回归来，因此会存在线性回归的一些缺点。
多元线性回归的缺点： 当自变量的数量大于样本量的时候，解不出$\theta$ ,$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$ ，设$X_{n\ast k}$,当 k>n 时，$(X^{T}X{)}_{k\ast k}$的秩为n，不是满秩的，所以没有逆矩阵 。
PCA 的缺点：PCA只考虑了自变量的方差，然后选取方差最大的几个正交变量，可以用于解决共线性问题(计量)，没有考虑自变量对因变量的贡献。
偏最小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线性回归分析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回归模型外，还可以同时完成一些类似于主成分分析和典型相关分析的研究内容，提供更丰富、深入的一些信息。
对于 P 个因变量$y_1,...,y_p$与 m 个自变量$x_1,...,x_m$的回归问题。首先在自变量集中提出第一成分$u_1$($u_1$是$x_1,...,x_m$的线性组合，且尽可能多地提取原自变量集中的变异信息)；同时在因变量集中也提取第一成分$v_1$，并要求$u_1$与$v_1$相关程度达到最大。然后建立因变量$y_1,...,y_p$与$u_1$的回归，重复这个过程直到提取到足够的指定的成分。

# 3.适用情况

偏最小二乘回归是集【主成分分析】，【典型相关分析】和【多元线性回归分析】3种分析方法的优点于一身的方法。

• 多元线性回归的缺点： 当自变量的数量大于样本量的时候，解不出$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$，设$X_{n\ast k}$,当 k>n 时，$(X^{T}X{)}_{k\ast k}$ 的秩为n，不是满秩的，所以没有逆矩阵Rank(A,B)<=Rank(B)

• PCA的缺点：PCA只考虑了自变量的方差，然后选取方差最大的几个正交变量，可以用于解决共线性问题(计量)，没有考虑自变量对因变量的贡献

• PLS：偏最小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的个数很多，且都存在多重相关性，而样本又比较少的时候

# 4.计算步骤

先将 X 与 Y 标准化:

$A=\left[\begin{array}{ccc}{x}_{11}& ...& x_{1m}\\ ...& & ...\\ x_{n1}& ...& x_{nm}\end{array}\right],B=\left[\begin{array}{ccc}{y}_{11}& ...& y_{1m}\\ ...& & ...\\ y_{n1}& ...& y_{nm}\end{array}\right]$
第一步：
分别提取两组 X 与 Y 变量的第一对成分$u_1$和$v_1$，并使他们的相关性达到最大。假设从两组变量中第一成分为$u_1$和$v_1$，$u_1$是自变量集$X=[x_1,...,x_m{]}^T$的线性组合，$v_1$是自变量集$Y=[y_1,...,y_p{]}^T$的线性组合：
$\{\begin{array}{r}{u}_1={\rho }_1^{T}X\\ v_1={\gamma }_1^{T}Y\end{array}$ ，$(u_{m\ast m}={\rho }_{n\ast m}^T\cdot X_{n\ast m},u_1\mathrm{为}to{p}_1\mathrm{主}\mathrm{成}\mathrm{分},v\mathrm{相}\mathrm{同})$

为了回归分析的需要，要求：
(1)$u_1$ 和 $v_1$各自尽可能多地提取所在变量组的变异信息；
(2))$u_1$ 和 $v_1$的相关程度达到最大。

第二步：计算${\rho }_1$与${\gamma }_1$

1. 最大化协方差，使得)$u_1$ 和 $v_1$的相关程度达到最大，可以用得分向量)${\hat{u}}_1$ 和 ${\hat{v}}_1$的内积来计算：

   $max<{\hat{u}}_1,{\hat{v}}_1>={\rho }_1^{T}AB{\gamma }_1$
   $s.t.=\{\begin{array}{r}{\rho }_1^T{\rho }_1=1\\ {\gamma }_1^T{\gamma }_1=1\end{array}$

2. 采用拉格朗日乘数法，问题化为求单位向量${\rho }_1$与${\gamma }_1$，使${\gamma }_1$达到最大，问题求解只需计算$M=A^{T}B{B}^{T}A$的特征值与特征向量，且M的最大特征值为${\theta }_1^2$，相应的特征向量就是所要求解的${\rho }_1$，进而也能得到${\gamma }_1,{\gamma }_1=\frac{1}{{\theta }_1}{B}^{T}A{\rho }_1$

第三步
由两组变量集的标准化观察数据矩阵 X 和 Y，可以计算第一对成分的得分向量，记为${\hat{u}}_1$ 和 ${\hat{v}}_1$，${\hat{u}}_1=A{\rho }_1$，${\hat{v}}_1=B{\gamma }_1$

1. 建立<$y_1,...,y_p$对$u_1$的回归及$x_1,...,x_m$对$u_1$的回归，假定回归模型:
   $\{\begin{array}{r}A={\hat{u}}_1{\sigma }_1^T+A_1\\ B={\hat{u}}_1{\tau }_1^T+B_1\end{array}$

其中，${\sigma }_1^T=[{\sigma }_1,...,{\sigma }_m]$,${\tau }_1^T=[{\tau }_1,...,{\tau }_m]$分别是多对一回归模型中的参数向量，$A_1,B_1$是残差阵。

2. 回归系数${\sigma }_1,{\tau }_1$的最小二乘估计为

   $\{\begin{array}{r}{\sigma }_1=\frac{A^T{\hat{u}}_1}{||{\hat{u}}_1|{|}^2}\\ {\tau }_1=\frac{B^T{\hat{u}}_1}{||{\hat{u}}_1|{|}^2}\end{array}$

3. 用残差阵$A_1\mathrm{和}{B}_1$代替$A,B$，重复以上步骤，直到残差阵中元素的绝对值近似为0，每进行一次得到一个${\sigma }_1\mathrm{和}{\tau }_1$。

第四步
重复上面的步骤，得到 r 个成分

$\{\begin{array}{r}A={\hat{u}}_1{\sigma }_1^T+...+{\hat{u}}_r{\sigma }_r^T+A_r\\ B={\hat{u}}_1{\tau }_1^T+...+{\hat{u}}_r{\tau }_r^T+B_r\end{array}$
将$u_1={\rho }_1^{T}X$代入$Y={\hat{u}}_1{\tau }_1^T+...+{\hat{u}}_r{\tau }_r^T$，即得P个因变量的偏最小二乘回归方程式$y_j=c_{j1}{x}_1+...+c_{jm}{x}_m$,$j=1,2,...,p$。

第五步：交叉有效性检验
应该提取多个成分，可以使用交叉有效性检验，每次舍去第 i 个观察数据，对余下的 n-1个观测数据用最小二乘回归方法，并考虑抽取 h(h<=r) 个成分后拟合的回归式，然后把舍去的自变量组第 j 个观测数据代入所拟合的回归方程式，得到$y_j(j=1,2,...,p)$在第 i 观测点上的预测值为。对$i=1,2,...,n$重复以上的验证，即得抽取 h 个成分时第 j 个因变量 $y_j(j=1,2,...,p)$的预测误差平方和为 $PRES{S}_j(h)=\sum _{i=1}^n(b_{ij}-b_{\hat{i}j}(h){)}^2,j=1,2,...,p$ $Y=[y_1,...,y_p{]}^T$ 的预测误差平方和为$PRESS(h)=\sum _{j=1}^{p}PRES{S}_j(h)$。 另外，再采用所有的样本点，拟合含 h 个成分的回归方程。这时，记第 i 个样本点的预测值为${\hat{b}}_{ij}(h)$，则可定义$y_i$的误差平方和为$SS(h)=\sum _{j=1}^{p}S{S}_j(h)$，当 PRESS (h) 达到最小值时，对应的 h 即为所求的成分个数。通常，总有 PRESS(h)>SS(h)，而 SS(h)<SS(h-1)。因此在提取成分时，总是希望$\frac{PRESS(h)}{SS(h-1)}$越小越好，一般可以设定阈值为0.05，判定规则为：$\frac{PRESS(h)}{SS(h-1)}<=(1-0.05{)}^2$时，新成分对回归改善是有帮助的，因此，可以定义交叉有效性$Q_h^2=1-\frac{PRESS(h)}{SS(h-1)}$,在每一步计算结束前，计算交叉有效性，在第 h 步有 $Q_h^2<1-{0.095}^2$，则模型到达精度，可以停止提取成分。

# 9、手推步骤

Step 1. 数据标准化
原始数据（前3行示例）：

体重($X_1$)	腰围($X_2$)	脉搏($X_3$)	单杠($Y_1$)	弯曲($Y_2$)	跳高($Y_3$)
191	36	50	5	162	60
189	37	52	2	110	60
193	38	58	12	101	101

计算均值和标准差：

${\mu }_{X_1}=\frac{191+189+193}{3}=191,{\sigma }_{X_1}=\sqrt{\frac{(191-191{)}^2+(189-191{)}^2+(193-191{)}^2}{3}}=1.63$

标准化后：

$X{1}_{std}=\frac{191-191}{1.63}=0,\frac{189-191}{1.63}=-1.23,\frac{193-191}{1.63}=1.23$

标准化矩阵 $A$（X）和 $B$（Y）：

$A=\left[\begin{array}{ccc}0& -0.87& -1.34\\ -1.23& 0& -0.45\\ 1.23& 1.15& 1.79\end{array}\right],B=\left[\begin{array}{ccc}-0.57& 1.15& -0.67\\ -1.43& -0.92& -0.67\\ 1.43& -1.38& 1.79\end{array}\right]$

Step 2：第一主成分计算
(1) 计算 $M=A^{T}B{B}^{T}A$

$A^{T}B=\left[\begin{array}{ccc}0& -1.23& 1.23\\ -0.87& 0& 1.15\\ -1.34& -0.45& 1.79\end{array}\right]\left[\begin{array}{ccc}-0.57& 1.15& -0.67\\ -1.43& -0.92& -0.67\\ 1.43& -1.38& 1.79\end{array}\right]=\left[\begin{array}{ccc}1.75& -3.38& 3.03\\ 1.64& -0.60& 1.15\\ 3.18& -3.12& 4.21\end{array}\right]$

$M=(A^{T}B)(A^{T}B{)}^T=\left[\begin{array}{ccc}1.75& -3.38& 3.03\\ 1.64& -0.60& 1.15\\ 3.18& -3.12& 4.21\end{array}\right]\left[\begin{array}{ccc}1.75& 1.64& 3.18\\ -3.38& -0.60& -3.12\\ 3.03& 1.15& 4.21\end{array}\right]=\left[\begin{array}{ccc}23.2& 12.4& 30.1\\ 12.4& 4.0& 13.2\\ 30.1& 13.2& 41.3\end{array}\right]$

(2) 求最大特征值和特征向量
解 $|M-\lambda I|=0$：

$\text{特征值：}{\lambda }_1=64.1,{\lambda }_2=3.8,{\lambda }_3=0.6$

最大特征值 ${\theta }_1^2=64.1$，对应特征向量：

${\rho }_1=[0.59,0.77,-0.24{]}^T(\text{标准化后})$

(3) 计算 ${\gamma }_1$

${\gamma }_1=\frac{B^{T}A{\rho }_1}{{\theta }_1}=\frac{[-0.34,-0.42,-0.14{]}^T}{8.01}=[-0.04,-0.05,-0.02{]}^T$

(4) 计算得分向量 $u_1$ 和 $v_1$

$u_1=A{\rho }_1=\left[\begin{array}{c}0\times 0.59+(-0.87)\times 0.77+(-1.34)\times (-0.24)\\ (-1.23)\times 0.59+0\times 0.77+(-0.45)\times (-0.24)\\ 1.23\times 0.59+1.15\times 0.77+1.79\times (-0.24)\end{array}\right]=\left[\begin{array}{c}-0.42\\ -0.67\\ 1.09\end{array}\right]$

Step 3：回归残差计算 (1) 回归系数 ${\sigma }_1$ 和 ${\tau }_1$

${\sigma }_1=\frac{A^T{u}_1}{\parallel u_1{\parallel }^2}=\frac{[-1.23,0.87,0.24{]}^T}{1.82}=[-0.68,0.48,0.13{]}^T$

${\tau }_1=\frac{B^T{u}_1}{\parallel u_1{\parallel }^2}=\frac{[-0.34,-0.42,-0.14{]}^T}{1.82}=[-0.19,-0.23,-0.08{]}^T$

(2) 更新残差矩阵 $A_1$ 和 $B_1$

$A_1=A-u_1{\sigma }_1^T=\left[\begin{array}{ccc}0& -0.87& -1.34\\ -1.23& 0& -0.45\\ 1.23& 1.15& 1.79\end{array}\right]-\left[\begin{array}{c}-0.42\\ -0.67\\ 1.09\end{array}\right][-0.68,0.48,0.13]=\left[\begin{array}{ccc}-0.29& -0.67& -1.28\\ -0.77& 0.32& -0.36\\ 1.97& 0.63& 1.65\end{array}\right]$

Step 4：后续主成分提取
重复上述步骤计算 ${\rho }_2,{\gamma }_2,u_2,v_2$ 和 ${\rho }_3,{\gamma }_3,u_3,v_3$，最后求得总体数据的各项值,最终可以得到：

$W=\left[\begin{array}{ccc}0.59& -0.37& 0.94\\ 0.77& 0.70& -0.80\\ -0.24& 0.64& 0.22\end{array}\right],C=\left[\begin{array}{ccc}-0.34& -0.34& 0.48\\ -0.42& -0.29& 0.46\\ -0.14& -0.07& -0.21\end{array}\right]$

Step 5:回归系数矩阵 $B$

$B=W(P^{T}W{)}^{-1}{C}^T=\left[\begin{array}{ccc}1.95& 17.97& -13.28\\ -4.66& -55.67& 0.75\\ -0.14& 1.01& -2.80\end{array}\right]$

最终回归方程:

$\{\begin{array}{l}\text{单杠}=9.45+1.95\cdot {\text{体重}}_{std}-4.66\cdot {\text{腰围}}_{std}-0.14\cdot {\text{脉搏}}_{std}\\ \text{弯曲}=145.55+17.97\cdot {\text{体重}}_{std}-55.67\cdot {\text{腰围}}_{std}+1.01\cdot {\text{脉搏}}_{std}\\ \text{跳高}=70.30-13.28\cdot {\text{体重}}_{std}+0.75\cdot {\text{腰围}}_{std}-2.80\cdot {\text{脉搏}}_{std}\end{array}$

# 10、参考文献

[1] Scientific Platform Serving for Statistics Professional 2021. SPSSPRO. (Version 1.0.11)[Online Application Software]. Retrieved from https://www.spsspro.com.
[2] Jong S D , Phatak A . Partial least squares regression[C]// International Workshop on Recent Advances in Total Least Squares Techniques & Errors-in-variables Modeling. Society for Industrial and Applied Mathematics, 1997.