拓端tecdat：R語言主成分回歸（PCR）、 多元線性回歸特征降維分析光譜數據和汽車油耗、性能數據

原文鏈接：http://tecdat.cn/?p=24152

原文出處：拓端數據部落公眾號

什么是PCR？（PCR = PCA + MLR）

• PCR是處理許多 x 變量的回歸技術
• 給定 Y 和 X 數據：
• 在 X 矩陣上進行 PCA
– 定義新變量：主成分（分數）
• 在 多元線性回歸(MLR)  中使用這些新變量中的一些來建模/預測 Y
• Y 可能是單變量或多變量。 

例子

1.  
    
2.  
   # 對數據
3.  
   set.seed(123)
4.  
    
5.  
   da1 <- marix(c(x1, x2, x3, x4, y), ncol = 5, row = F)
6.  
    

多元線性回歸和逐步剔除變量，手動：

1.  
   # 對於data1：（正確的順序將根據模擬情況而改變）。
2.  
   lm(y ~ x1 + x2 + x3 + x4)
3.  
    
4.  
   lm(y ~ x2 + x3 + x4)
5.  
    
6.  
    
7.  
   lm(y ~ x2 + x3)
8.  
    
9.  
    
10.  
    lm(y ~ x3)
11.  
     

 

 

配對關系圖

pais(atix, ncol = 5, byrow = F

 

如果重復： 

1.  
   # 對於data2:
2.  
    
3.  
   lm(y ~ x1 + x2 + x3 + x4)
4.  
    
5.  
   lm(y ~ x1 + x2 + x4)
6.  
    
7.  
    
8.  
    
9.  
   lm(y ~ x2 + x4)
10.  
     
11.  
    lm(y ~ x2)

 

 

 數據集 2 的繪圖：

 

使用四個 x 的均值作為單個變量來分析兩個數據集： 

1.  
   xn1 <- (dt1[,1] + a1[,2] + at1[,3] + dt1[,4])/4
2.  
   lm(data1[,5] ~ xn1)
3.  
   lm(data2[,5] ~ xn2)
4.  
    

 

檢查一下X數據的PCA的載荷loading是什么。

1.  
   # 幾乎所有的方差都在第一主成分解釋。
2.  
   prnmp(dt1[,1:4])
3.  
    

1.  
   # 第一個成分的載荷
2.  
   picp(dta1[,1:4])$lads[,1]

它們幾乎相同，以至於第一個主成分本質上是四個變量的平均值。讓我們保存一些預測的 beta 系數 - 一組來自數據 1 的完整集和一組來自均值分析的：

1.  
   c1 <- smry(lm(dta1[,5] ~ dta1[,1] + dta1[,2] + ata1[,3] +
2.  
   dt1[,4]))$coficns[,1]
3.  
   f <- summry(rm2)$cefets[,1]

我們現在模擬三種方法（完整模型、均值（=PCR）和單個變量）在 7000 次預測中的表現：

1.  
    
2.  
   # 對預測進行模擬。
3.  
   誤差<- 0.2
4.  
    
5.  
   xn <- (x1 + x2 + x3 + x4)/4
6.  
   yt2 <- cf[1] + cf[2] * xn
7.  
   yht3 <- cf[1] + cf[2] * x3
8.  
   bro(c(um((y-hat)^2)/7000 min = "平均預測誤差平方")
9.  
    

PCR 分析誤差最小。

示例：光譜類型數據

構建一些人工光譜數據：（7 個觀測值，100 個波長）
 

1.  
    
2.  
    
3.  
   # 光譜數據實例
4.  
    
5.  
   mapot(t(spcra) )
6.  
   mtlnes(t(spcra))
7.  
    
8.  
    

 

平均光譜表明： 

1.  
   mtpot(t(secra))
2.  
   malies(t(spcta))
3.  
   mnp <- apply(spcra, 2, mean)
4.  
   lines(1:100, mnp, lwd = 2)

 

平均中心光譜： 

1.  
   spcamc<-scae(spcta,scale=F)
2.  
   plot(t(spermc),tpe="")

標准化光譜： 

1.  
   sptracs<-scale(spetra,scale=T,center=T)
2.  
   matott(specrams),tye="n",
3.  
   matlies(t(sectramcs))

 

1.  
    
2.  
   # 用特征函數對相關矩陣做PCA。
3.  
   pcaes <- eien(cor(spra))
4.  
   ladigs <- pces$vectors[,1].
5.  
   score <- peramcs%*%t(t(lodis1))
6.  
   pred <- soes1 %*% loadings1
7.  
   ## 1-PCA預測值轉換為原始尺度和平均值。
8.  
   mtrx(repeasp, 7), nro=7, brw=T)
9.  
    

 在單個概覽圖中收集的所有圖：

1.  
   par(mfrow = c(3, 3)
2.  
   matlot(t(sectr)

 

PCR是什么？

• 數據情況：

• 用A 主成分t1、t2... 做MLR而不是所有（或部分）x。
• 多少個成分：通過交叉驗證確定。

怎么做？

1. 探索數據
2. 進行建模（選擇主成分數量，考慮變量選擇）
3. 驗證（殘差、異常值、影響等）
4. 迭代 2. 和 3。
5. 解釋、總結、報告。
6. 如果相關：預測未來值。

交叉驗證

• 忽略一部分觀察值
• 在剩余（減少的）數據上擬合模型
• 預測模型遺漏的觀察值：yˆi,val
• 對所有觀察值依次執行此操作並計算總體模型性能：

（預測的均方根誤差）

最后：對所有選擇的分量（0、1、2、...、... ）進行交叉驗證並繪制模型性能

barplot(names.arg)

 

選擇最佳成分數：
• 總體誤差最小的主成分。

重采樣

• 交叉驗證 (CV)

•留一法(Leave-One-Out,簡稱LOO)

• Bootstrapping
• 一個很好的通用方法：
– 將數據分成訓練集和測試集。
– 在訓練數據上使用交叉驗證
– 檢查測試集上的模型性能
– 可能：重復所有這些多次（重復雙交叉驗證）

交叉驗證 - 原則

• 最小化預期預測誤差：
平方預測誤差 = Bias2 +方差
• 包括“許多”PC主成分：低偏差，但高方差
• 包括“很少”PC 主成分：高偏差，但低方差
• 選擇最佳折衷！

驗證 - 存在於不同的級別

1. 分為 3 個：訓練（50%）、驗證（25%）和測試（25%）
2. 拆分為 2：校准/訓練 (67%) 和測試 (33%) 
訓練中，CV/bootstrap •更常用
3. 沒有 "固定分割"，而是通過CV/bootstrap反復分割，然后在每個訓練組內進行CV。
4. 沒有分割，但使用（一級）CV/bootstrap。
5. 只對所有數據進行擬合--並檢查誤差。

示例：汽車數據

1.  
    
2.  
   # 例子：使用汽車數據。
3.  
   # 將X矩陣定義為數據框中的一個矩陣。
4.  
   mtas$X <- as.ix(mcas[, 2:11])
5.  
   # 首先，我們考慮隨機選擇4個屬性作為測試集
6.  
   mtcrs_EST<- mtcrs[tcars$rai == FASE,] 。
7.  
   tcaTRAIN <- mtars[tcarstrai == TUE,] 。
8.  
    

現在所有的工作都在 訓練數據集上進行。

探索數據

我們之前已經這樣做了，所以這里不再贅述

數據建模

使用pls軟件包以最大/大量的主成分運行PCR。
 

1.  
    
2.  
   # 使用pls軟件包，以最大/較大的成分數運行PCR。
3.  
   pls(lomg ~ X , ncop = 10, dta = marsTRAN,
4.  
   aliaon="LOO")
5.  
    

 初始圖集：

1.  
    
2.  
   # 初始化的繪圖集。
3.  
   par(mfrow = c(2, 2)
4.  
   plot(mod)

 

主成分的選擇： 

1.  
    
2.  
   # 主成分的選擇。
3.  
   # 分段的CV會得到什么。
4.  
   modseCV <- pcr(lomg ~ X , ncp = 10, dta = marTIN
5.  
   vai ="CV"
6.  
   )
7.  
   # 初始圖集。
8.  
   par(mfrow = c(1, 2))
9.  
   plot(odsC, "vadaion")

 

讓我們看看更多的主成分： 

1.  
   # 讓我們看看更多的主成分。
2.  
   # 分數。
3.  
   scre(mod)

 

1.  
   #負荷
2.  
   loading(md,cms = 1:4)

 

我們選擇 3 個主成分： 

1.  
    
2.  
   # 我們選擇4個成分
3.  
   m <- ncmp = 3, data = mrs_TAI vdon = "LOO", akknie = RUE
4.  
    

然后： 驗證：
讓我們驗證更多：使用 3 個主成分。我們從中獲取預測的殘差，因此這些是（CV）驗證版本！

1.  
    
2.  
    
3.  
   oit <- ppo(mod3, whih = "litin")
4.  
   plot(obft[,2], Rsds)
5.  
   # 為了繪制殘差與X-杠桿的對比，我們需要找到X-杠桿。
6.  
   # 然后找到杠桿值作為Hat矩陣的對角線。
7.  
   # 基於擬合的X值。
8.  
   Xf <- sors(md3)
9.  
   plot(lvge, abs(Rsidals))
10.  
    text(leage, abs(Reuls))
11.  
     

 

1.  
    
2.  
   # 讓我們也繪制一下殘差與每個輸入X的關系。
3.  
    
4.  
   for ( i in 2:11){
5.  
   plot(res~masAN[,i],type="n")
6.  
   }
7.  
    

 

解釋/結論

現在讓我們看一下結果——“解釋/結論”：

1.  
    
2.  
   # 現在我們來看看結果 - 4) "解釋/結論"
3.  
   par(mfrw = c(2, 2))
4.  
   # 繪制具有Jacknife不確定性的系數。
5.  
   obfi <- red(mod3,, wich = "vltn)
6.  
   abe(lm(ft[,2] ~ fit[,1])
7.  
   plt(mo3, ses = TUE,)
8.  
    
9.  
    

1.  
   # 最后是一些輸出
2.  
   test(mo3, nm = 3)

預測

1.  
    
2.  
   # 現在讓我們試着預測TEST集的4個數據點。
3.  
   prdit(md3, nwaa =TEST)
4.  
   plt(TEST$lgg, pes)
5.  
    

1.  
   rmsep <- sqrt(men(log - prd)^2))
2.  
   rmsep

---

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應用案例

2.面板平滑轉移回歸(PSTR)分析案例實現

3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）

4.R語言泊松Poisson回歸模型分析案例

5.R語言回歸中的Hosmer-Lemeshow擬合優度檢驗

6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實現

7.在R語言中實現Logistic邏輯回歸

8.python用線性回歸預測股票價格

9.R語言如何在生存分析與Cox回歸中計算IDI，NRI指標