偏最小二乘方法

第六章偏最小二乘方法偏最小二乘方法(PLS-PartialLeastSquares))是近年来发展起来的一种新的多元统计分析法,现已成功地应用于分析化学,如紫外光谱、气相色谱和电分析化学等等。该种方法，在化合物结构-活性/性质相关性研究中是一种非常有用的手段。如美国Tripos公司用于化合物三维构效关系研究的CoMFA(ComparativeMolecularFieldAnalysis)方法,其中，数据统计处理部分主要是PLS。在PLS方法中用的是替潜变量，其数学基础是主成分分析。替潜变量的个数一般少于原自变量的个数，所以PLS特别适用于自变量的个数多于试样个数的情况。在此种情况下，亦可运用主成分回归方法，但不能够运用一般的多元回归分析，因为一般多元回归分析要求试样的个数必须多于自变量的个数。§6.1多元线性回归（MLR）若自变量为m个，xj(j=1,2,…,m)，因变量为y，在y与xj间，我们可以建立一线性模型，即exbxbxbymm...2211(6.1a)exbyjmjj1(6.1b)ebxy(6.1c)在式中，bj为回归系数。在式（6.1）中仅有一个试样，若有n个试样，即为yi(i=1,2,…,n)，它的列向量形式为y，b与原来相同，矢量xj’为矩阵X的行，则：y=Xb+e若用图形表示，则为：y=XB+e1m11nnnm在此情况下，n为试样数，m为自变量数。有如下三种情况：(1)mn，即变量数多于试样数，对于b来说，则有无穷多个解。(2)m=n，变量数与试样数相等，若矩阵X满秩时，则矢量b有唯一解。但是，在实际工作中，这种情况是极少能碰到的。此时我们有：e=y–Xb=0(3)mn，变量数小于试样数，尽管我们得不到准确解b，但是可以使残差矢量e尽可能小而得到解，e=y–Xb这就是我们所熟知的最小二乘法。其解为：yXXXb1)((6.2)bybbx1)(在上边的叙述中，因变量为1个，而事实上可以有多个因变量。如有两个因变量y1和y2，我们可以简单地写成两个线性方程：y1=Xb1+e;y2=Xb2+e若用矩阵标表示，则：nnnnnnxxxxxxxxxX....................212222111211nnyyyyyyyyY212212121121......)(mmbbbbbbbbB212212211121......)(nneeeeeeeeE212212211121......)(由此得到Y=XB+E对于2-P个因变量的图形表示为：Y=XB+E2-p2-p2-pnmnmn最小二乘的解为：YXXXB1)((6.3)多元线性回归应用很广泛，因为在许多情况下该种方法具有良好的性能。但是，此种方法也有固有的缺点。假若体系的响应（即因变量）呈现线性，无干扰，无溶液间的相互作用，低噪声无共线性，则多元线性回归是一种非常好的方法。事实上，完全满足上述条件比较困难。当噪声较强，或干扰较严重时，有可能导致所得数学模型失真，如下例：1241576917621896821326310215275X2863129334572Y运用式(6.3)则可得B矩阵：050280080240410420480550710.........B所用数学模型有效性的量度可用Err：211)ˆ(ikKkIiikrryyEKkIiik112式中，yik为矩阵Y中第i行第k列的矩阵元，为由矩阵B所得的计算值，ik为前面所介绍的矩阵E的矩阵元。此例中，Err=0.49。若由于噪音使得X增广一列（注意：对于试样浓度的测定，它并不包含有用信息），即：5112415769741762189636821326391102152752X2863129334572Y由此得到的B矩阵为：0100301200302002402001904204201807102............B对于此模型，Err=0.07。它比前者为小，这就意味着对于矩阵Y，第二个数学模型比第个要更有效，这是一种假象。由于X中引入最后一列，使得B2中上部3*3部分与前边所提B不相等（B为真实模型）。由B2计算所得Y尽管误差要小，但其数学模型所描述的自变量与因变量间的关系并不真实。其原因主要为多元线性回归方法是采用整个X矩阵来建立数学模型，而并不顾及在X中的信息与真实模型相关与否。很显然，若所得结果偏离了其实际数学模型，则对于未知试样的预测也是错误的。为了克服多元线性回归的不足，在数学方法上引进了主成分回归方法（PCR）。§6.2主成分回归主成分回归可分为两步：测定主成分数，并由主成分分析将X矩阵降维；对于降维的X矩阵再进行线性回归分析。主成分分析的概念在前一章已经作了介绍。所谓主成分，它为一新的变量，而该新变量是原变量xij的线性组合。第一个主成分所能解释原变量的方差量最大，第二个次之，第三个再次之，等等。也就是说，主成分是一种线性组合，用它来表征原来变量时所产生的平方误差最小。运用主成分分析，原变量矩阵X可以表达为得分（即主成分）矩阵T，而T由X在本征矢量P上的投影所得。主成分与矩阵X的本征矢量一一对应，即T=XP。设矩阵X的阶为I*J，若T的阶与J相等，则主成分回归与多元线性回归所得结果相同，并不能显示出主成分回归的优越之处。选取的主成分数一般应该比J小，而删去那些不重要的主成分，因为这些主成分所包含的信息主要是噪声，由此所得的回归方程稳定性较好。另外，由X所定义的空间可以进一步来说明主成分回归与多元线性回归的区别。多元线性回归应用了由X的列所定义的全部空间，而主成分回归所占用的是一子空间。当X的J列中，有一列可为其它J—1列的线性组合时，则X可用J-1列的矩阵T来描述，而并不丢失信息。新的矩阵T定义了X的一个子空间。综合上述，X可由它的得分矩阵T来描述（由于删去与小的本征值相应的维，所以T的维小于X的维）：T=XP若用图形表示，则为：T=XPamannm由此可得多线性方程：Y=TB+E其解为：YTTTB1)(其图形表示为：Y=TB+Epppannna主成分分析可以解决共线问题，同时由于去掉了不太重要的主成分，因而可以削弱噪声（随机误差）所产生的影响。但是，由于主成分回归为二步法，若在第一步中消去的是有用的主成分，而保留的是噪声，则在第二步多元线性回归所得结果就将偏离真实的数学模型。§6.3偏最小二乘（PLS）§6.3.1基本原理为了叙述上的方便，我们首先引进“因子”的概念。一个因子为原来变量的线性组合，所以矩阵的某一主成分即为一因子，而某矩阵的诸主成分是彼此相互正交的，但因子不一定，因为一因子可由某一成分经坐标旋转而得。在主成分回归中，第一步，在矩阵X的本征矢量或因子数测试中，所处理的仅为X矩阵，而对于矩阵Y中信息并未考虑。事实上，Y中亦可能包含非有用的信息。所以很自然的一种想法是，在矩阵X因子的测试中应同时考虑矩阵Y的作用。偏最小二乘正是基于这种思想的一种回归方法。偏最小二乘和主成分分析很相似，其差别在于用于描述变量Y中因子的同时也用于描述变量X。为了实现这一点，在数学上是以矩阵Y的列去计算矩阵X的因子，与此同时，矩阵Y的因子则由矩阵X的列去预测。其数学模型为：EPTX（6.4）及FQUY(6.5)此处，T和U的矩阵元分别为X和Y的得分，而P和Q的矩阵元分别为X和Y的装载，E和F分别为运用偏最小二乘模型法去拟合X和Y所引进的误差。T=XP(主成分分析)TP’=XPP’PP’=IX=TP’(因子分析)在理想的情况下，X中误差的来源和Y中的误差的来源完全相同，即影响X与Y的因素相同。但实际上，X中误差与Y中误差并不相关，因而t≠u，但当两个矩阵同时用于确定因子时，则X和Y的因子具有如下关系：u=bt+e（6.6）式中b所表征的即为u和t间的内在关系。为了使因子T既可描述X矩阵，同时又可描述Y矩阵，则需采取折衷方案，即将T进行坐标旋转。显然，坐标旋转后的T因子对于X矩阵的表达已不再是最优的状况。如假设X矩阵和Y矩阵均为6*3，即行为6，列为3。在列空间，X和Y矩阵的行分别示于图6.1（上部）。PLS第一个因子（t和u）方向在各自的空间均可解释试样的最大偏差。若PLS模型是正确的，将t对u作图则可得一线性关系。事实上，PLS要将各自空间中的因子进行折衷以增加t对u的相关性（图6.1下部）。由于这种折衷才可使所得数学模型较好地同时描述X和Y。在行空间，情况与列空间类同。图6.1PLS处理的图形表示如有矩阵（见§6.2）：1241576917621896821326310215275X2863129334572Y数据的预处理为：每列减去相应列的平均值（mean-centered），PLS所得结果为：8.157.45.20t100.10u将t对u作图（图6.2）可显示出二者的线性关系，其斜率b=0.53。图6.2矩阵X的因子t对矩阵Y的因子u作图对于未知试样的预测，要应用X和Y的得分模型及相关性bi。若有L个因子，则bl为表达第l个因子相关性的系数，其步骤为：由未知试样的测定值x末通过校正模型（式（6.4）计算出t末，进而由（式6.6）及bl可计算未知试样的得分矢量u末，最后由校正模型(式6.5)得未知试样含量。EPTXFQUYu=bt+e(6.4)(6.5)(6.6)§6.3.2偏最小二乘算法1.校正模型的建立首先我们从一最简单的模型开始，然后给出偏最小二乘的完整算法。若仅有二矩阵块(block),即X块和Y块。对于X：(1)将某xj赋值给tstart，即tstart=xj；);/(/)2(uuXuttXtp;/)3(oldoldnewppp;/)4(ppXpt(5)比较步(2)和步(4)中的t，若二者相等，则停，否则转到(2)。对于Y：(1)将某yj赋值给ustart即ustart=yi;)/(/)2(ttYtuuYuq;/)3(oldoldnewqqq;/)4(qqYqu(5)比较步(2)和步(4)中的u，若二者相等，则停，否则到步(2)。在上述的算法中，X和Y是分别独立进行的，为了建立二者内在的相关性，则将得分t和u在步(2)中的位置相交换（上述算法中的括号内部分）；istartyu(1)令);/(/)2(uuXuuuXupoldoldnewoldoldnewppp/(/)3();/(/)4(XtppXpt;/)5(ttYtq;/)6(oldoldnewqqq;/)7(qqYqu(8)将步(4)中t与前一次迭代所得t相比较，若二者相等（有一定的舍入误差），则停，否则转入步(2)[若Y为一维，即仅一个变量，则跳过步(5)—(8)，并置q=1]。此算法一般收敛很快。所得到的为X和Y的经过旋转的主成分，即t不互相正交，其原因是在主成分计算中，运算的顺序发生了变化。因此，将权重w’（见上述运算中括号内等式）替代p’，并在收敛之后，再加入：ttXtp/以得到正交的t值。由，则可计算新的t，t=Xp/p，其实此即为oldoldnewppp/p.oldoldnewpttT的相互正交并非绝对必要，但当与主成分回归比较时，t正交的条件还是需要满足的。当预测时，需将ω’作同样的标准化处理：，否则，将引入误差。然后，t可用于内部的相关：（下脚意为对于h因子，大小为n*1），此处。oldoldnewhhhtbuˆhhhhtttub/其残差的计算分别为：0*