【數學建模】偏最小二乘回歸分析（PLSR）

PLSR的基本原理與推導，我在這篇博客中有講過。

 

0.
偏最小二乘回歸集成了多元線性回歸、主成分分析和典型相關分析的優點，在建模中是一個更好的選擇，並且MATLAB提供了完整的實現，應用時主要的問題是：

• 注意檢驗，各種檢驗參數：有關回歸的檢驗以及有關多元分析的檢驗
• 系數眾多，容易混淆
• 要清楚原理才能寫好論文
• 注意matlab函數plsregress的眾多返回值
• 例如累計貢獻度，建模時最好列出表格

1.

問題：

自變量組 X = [x1,x2…xn]  （n組自變量）

因變量組 Y = [y1,y2,…yp]   （p組因變量）

考慮到X、Y內部之間的多重相關性，可以使用PLSR建立Y對X的多元回歸模型。這是一種多對多回歸的模型。

 

偏最小二乘回歸的實現步驟：

1. X、Y標准化。若考慮標准化的不對等特性，考慮實現對應分析。
2. 求相關系數矩陣。可以把X、Y統一放到一個增廣矩陣中，實現求列向量之間的相關系數矩陣（corrcoef實現無需標准化，直接使用原始數據）
3. 求主成分對。（求出自變量與因變量的成分，類似於典型相關分析）這里對數其實是min(n-1,p)。求出<u1,v1>、<u2,v2>… 實際上，u、v是原始變量標准化后的線性組合、即投影。
4. 計算貢獻率表格。計算前k個主成分u對原始變量X的貢獻率、v對Y的貢獻率（函數直接返回結果）。
5. 根據貢獻率表格，選取k個主成分對。一般累計貢獻率達到90%合適。
6. 求出原始變量X對這k個主成分u的回歸方程以及Y對u的（不是v！）回歸方程。
7. 根據6的結果，可以求出因變量組Y與自變量組X的回歸方程，但這其實是標准化了的（常數項一定是0），進一步可以還原為真實原始變量的回歸方程，這也是我們所要求得的。
8. 模型的解釋與檢驗。

• • 首先得進行一個回歸檢驗：判定系數R方的檢驗（接近於1)。計算每一個回歸方程的R方，可以列出表格。
  • 之后進行交叉有效性檢驗：交叉系數Qh方 = 1 – （PRESS(h) / SS(h-1)）。這是從主成分分析的角度的檢驗，即檢驗提取的k個主成分。（這個檢驗比較復雜，詳細看推導）

2.

MATLAB實現命令：

[XL,YL,XS,YS,BETA,PCTVAR,MSE,stats] = plsregress(X,Y,ncomp)

param:

　　X: 標准化后的原始X數據，每行一個數據組，每列是一項指標，即一個自變量

　　Y：標准化后的原始Y數據，每行一個數據組，每列是一項指標，即一個因變量

　　ncomp：選取的主成分對數

return：

　　XL：自變量的負荷量矩陣。維度是（自變量數*ncomp）。每行是原始數據X對主成分u的回歸表達式的系數

　　YL：因變量的負荷量矩陣。維度是（自變量數*ncomp）。每行是原始數據Y對主成分u的回歸表達式的系數

 

　　XS：對應於主成分u的得分矩陣(得分說的是主成分的值)。每列是一個主成分得分向量。

　　　　如：每一列是一個主成分ui的值！列數是主成分數。

　　　　說明:主成分u1是個列向量.

　　YS：對應於主成分v的得分矩陣。每列是一個v對原始數據Y的線性組合的系數

 

　　BETA：最終的回歸表達式系數矩陣。每一列對應的，是一個yi對X的回歸表達式系數。

 

　　PCTVAR：兩行的矩陣。

　　　　第一行的每個元素代表着自變量提出主成分，相應主成分u的貢獻率。（特征值之比，詳細見主成分推導）

　　　　第二行的每個元素代表着因變量提出主成分，相應主成分v的貢獻率。這個貢獻率其實是主成分對原始變量的解釋能力大下。

　　MSE：兩行的矩陣。剩余標准差矩陣。第一行的第j個元素對應着自變量與他的前j-1個提出成分之間的剩余標准差。第二行對應因變量。

　　stats：返回4個值。結構體：stats。

　　    W — A p-by-ncomp matrix of PLS weights so that XS = X0*W.

　　　　    W = a\XS。 W每行是一個主成分得分向量的系數，如：

　　    T2 — The T2 statistic for each point in XS.

　　    Xresiduals — The predictor residuals, that is, X0-XS*XL'.

　　    Yresiduals — The response residuals, that is, Y0-XS*YL'.

 

 

3.

案例實現：

求Y對X的偏最小二乘回歸方程：

原始數據：

(前三列為X變量，后兩列為Y變量，共20組樣本。以下數據保存為pz.txt與matlab源文件同一文件夾下)

 

191 36 50 5 162 60  
 189 37 52 2 110 60
 193 38 58 12 101 101 
 162 35 62 12 105 37 
 189 35 46 13 155 58 
 182 36 56 4 101 42 
 211 38 56 8 101 38 
 167 34 60 6 125 40 
 176 31 74 15 200 40 
 154 33 56 17 251 250 
 169 34 50 17 120 38 
 166 33 52 13 210 115 
 154 34 64 14 215 105 
 247 46 50 1 50 50 
 193 36 46 6 70 31 
 202 37 62 12 210 120 
 176 37 54 4 60 25 
 157 32 52 11 230 80 
 156 33 54 15 225 73 
 138 33 68 2 110 43 

% PLSR 偏最小二乘

clc,clear
 ab0 = load('pz.txt');
 mu = mean(ab0);%均值
sig = std(ab0);% 標准差
rr = corrcoef(ab0) %相關系數矩陣
ab = zscore(ab0); %數據標准化
a = ab(:,[1:3]); %標准化的X
 b = ab(:,[4:end]); %標准化的Y
 % pls命令需要標准化變量
[XL,YL,XS,YS,BETA,PCTVAR,MSE,stats] = plsregress(a,b)
 contr = cumsum(PCTVAR,2) %每行累計求和，即計算累計貢獻率
XL
 YL
 XS
 YS
 xw = a\XS %自變量提出主成分系數，每列對應一個成分,這個就是stats.W
 yw = b\YS %因變量提出的主成分系數
ncomp = input('輸入主成分個數')
[XL2,YL2,XS2,YS2,BETA2,PCTVAR2,MSE2,stats2] = plsregress(a,b,ncomp)
 n = size(a,2);% n是自變量個數
m = size(b,2);
 %求原始數據回歸方程常數項
beta3(1,:) = mu(n+1:end) - mu(1:n)./sig(1:n)*BETA2([2:end],:).*sig(n+1:end);
 %求原始數據x1,x2...xn的系數，每一列是一個回歸方程
beta3([2:n+1],:) = (1./sig(1:n))'*sig(n+1:end).*BETA2([2:end],:)
 bar(BETA2','k') %畫直方圖

求解結果（部分）

假設采用2個主成分

ncomp =

     2

系數：
XL2 =

   -4.1306    0.0558
    -4.1933    1.0239
     2.2264    3.4441

YL2 =

    2.1191   -0.9714
     2.5809   -0.8398
     0.8869   -0.1877

主成分得分（每列一個主成分）：
XS2 =

   -0.1036   -0.2050
    -0.1241   -0.0577
    -0.1463    0.1807
     0.1110    0.2358
    -0.0785   -0.3927
    -0.0369    0.0249
    -0.2263    0.0263
     0.1199    0.0730
     0.2765    0.2263
     0.1874   -0.0577
     0.0588   -0.2428
     0.1198   -0.2420
     0.1913    0.2625
    -0.7077    0.2635
    -0.1327   -0.3375
    -0.1208    0.1803
    -0.0633    0.0707
     0.1933   -0.2712
     0.1690   -0.1291
     0.3131    0.3917

YS2 =

   -1.2834    0.1794
    -4.6311    1.3388
    -0.2845   -0.6256
    -1.2265    0.6851
     1.6002   -1.0788
    -4.5120    1.5408
    -2.9777   -0.0114
    -2.7548    1.5473
     3.9469   -0.4253
    10.4846   -2.6373
     1.4139   -0.6681
     4.8549   -1.1547
     5.2890   -1.0550
    -7.6800   -0.1989
    -5.1793    1.2090
     4.5405   -2.0460
    -6.4973    2.0374
     4.2728   -0.6046
     5.5489   -1.3537
    -4.9251    3.3215

 

標准化數據回歸方程系數（可以看到常數項系數是0）
BETA2 =

    0.0000    0.0000    0.0000
    -0.0773   -0.1380   -0.0603
    -0.4995   -0.5250   -0.1559
    -0.1323   -0.0855   -0.0072

 

貢獻率：
PCTVAR2 =

    0.6948    0.2265
     0.2094    0.0295

 

剩余標准差：
MSE2 =

    2.8500    0.8699    0.2242
     2.8500    2.2531    2.1689

stats2 =

             W: [3x2 double]
             T2: [20x1 double]
     Xresiduals: [20x3 double]
     Yresiduals: [20x3 double]

 

最終的回歸方程系數矩陣，每列一個方程:
beta3 =

   47.0375  612.7674  183.9130
    -0.0165   -0.3497   -0.1253
    -0.8246  -10.2576   -2.4964
    -0.0970   -0.7422   -0.0510

 

畫出回歸系數直方圖：

還可以用預測的方法做精度分析，在此略過。