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

原文出處：拓端數(shù)據(jù)部落公眾號

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

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

例子

1. 


2. # 對數(shù)據(jù)

3. set.seed(123)

4. 


5. da1 <- marix(c(x1, x2, x3, x4, y), ncol = 5, row = F)

6. 


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

1. # 對于data1：（正確的順序?qū)⒏鶕?jù)模擬情況而改變）。

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. 


?

?

配對關(guān)系圖

pais(atix, ncol = 5, byrow = F

?

如果重復(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)

?

?

?數(shù)據(jù)集 2 的繪圖：

?

使用四個(gè) x 的均值作為單個(gè)變量來分析兩個(gè)數(shù)據(jù)集：?

1. xn1 <- (dt1[,1] + a1[,2] + at1[,3] + dt1[,4])/4

2. lm(data1[,5] ~ xn1)

3. lm(data2[,5] ~ xn2)

4. 


?

檢查一下X數(shù)據(jù)的PCA的載荷loading是什么。

1. # 幾乎所有的方差都在第一主成分解釋。

2. prnmp(dt1[,1:4])

3. 


1. # 第一個(gè)成分的載荷

2. picp(dta1[,1:4])$lads[,1]

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

1. c1 <- smry(lm(dta1[,5] ~ dta1[,1] + dta1[,2] + ata1[,3] +

2. dt1[,4]))$coficns[,1]

3. f <- summry(rm2)$cefets[,1]

我們現(xiàn)在模擬三種方法（完整模型、均值（=PCR）和單個(gè)變量）在 7000 次預(yù)測中的表現(xiàn)：

1. 


2. # 對預(yù)測進(jìn)行模擬。

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 = "平均預(yù)測誤差平方")

9. 


PCR 分析誤差最小。

示例：光譜類型數(shù)據(jù)

構(gòu)建一些人工光譜數(shù)據(jù)：（7 個(gè)觀測值，100 個(gè)波長）
?

1. 


2. 


3. # 光譜數(shù)據(jù)實(shí)例

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="")

標(biāo)準(zhǔn)化光譜：?

1. sptracs<-scale(spetra,scale=T,center=T)

2. matott(specrams),tye="n",

3. matlies(t(sectramcs))

?

1. 


2. # 用特征函數(shù)對相關(guān)矩陣做PCA。

3. pcaes <- eien(cor(spra))

4. ladigs <- pces$vectors[,1].

5. score <- peramcs%*%t(t(lodis1))

6. pred <- soes1 %*% loadings1

7. ## 1-PCA預(yù)測值轉(zhuǎn)換為原始尺度和平均值。

8. mtrx(repeasp, 7), nro=7, brw=T)

9. 


?在單個(gè)概覽圖中收集的所有圖：

1. par(mfrow = c(3, 3)

2. matlot(t(sectr)

?

PCR是什么？

? 數(shù)據(jù)情況：

? 用A 主成分t1、t2...?做MLR而不是所有（或部分）x。
? 多少個(gè)成分：通過交叉驗(yàn)證確定。

怎么做？

1. 探索數(shù)據(jù)
2. 進(jìn)行建模（選擇主成分?jǐn)?shù)量，考慮變量選擇）
3. 驗(yàn)證（殘差、異常值、影響等）
4. 迭代 2. 和 3。
5. 解釋、總結(jié)、報(bào)告。
6. 如果相關(guān)：預(yù)測未來值。

交叉驗(yàn)證

? 忽略一部分觀察值
? 在剩余（減少的）數(shù)據(jù)上擬合模型
? 預(yù)測模型遺漏的觀察值：y?i,val
? 對所有觀察值依次執(zhí)行此操作并計(jì)算總體模型性能：

（預(yù)測的均方根誤差）

最后：對所有選擇的分量（0、1、2、...、... ）進(jìn)行交叉驗(yàn)證并繪制模型性能

barplot(names.arg)

?

選擇最佳成分?jǐn)?shù)：
? 總體誤差最小的主成分。

重采樣

? 交叉驗(yàn)證 (CV)

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

? Bootstrapping
? 一個(gè)很好的通用方法：
– 將數(shù)據(jù)分成訓(xùn)練集和測試集。
– 在訓(xùn)練數(shù)據(jù)上使用交叉驗(yàn)證
– 檢查測試集上的模型性能
– 可能：重復(fù)所有這些多次（重復(fù)雙交叉驗(yàn)證）

交叉驗(yàn)證 - 原則

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

驗(yàn)證 - 存在于不同的級別

1. 分為 3 個(gè)：訓(xùn)練（50%）、驗(yàn)證（25%）和測試（25%）
2. 拆分為 2：校準(zhǔn)/訓(xùn)練 (67%) 和測試 (33%)?
訓(xùn)練中，CV/bootstrap ?更常用
3. 沒有 "固定分割"，而是通過CV/bootstrap反復(fù)分割，然后在每個(gè)訓(xùn)練組內(nèi)進(jìn)行CV。
4. 沒有分割，但使用（一級）CV/bootstrap。
5. 只對所有數(shù)據(jù)進(jìn)行擬合--并檢查誤差。

示例：汽車數(shù)據(jù)

1. 


2. # 例子：使用汽車數(shù)據(jù)。

3. # 將X矩陣定義為數(shù)據(jù)框中的一個(gè)矩陣。

4. mtas$X <- as.ix(mcas[, 2:11])

5. # 首先，我們考慮隨機(jī)選擇4個(gè)屬性作為測試集

6. mtcrs_EST<- mtcrs[tcars$rai == FASE,] 。

7. tcaTRAIN <- mtars[tcarstrai == TUE,] 。

8. 


現(xiàn)在所有的工作都在 訓(xùn)練數(shù)據(jù)集上進(jìn)行。

探索數(shù)據(jù)

我們之前已經(jīng)這樣做了，所以這里不再贅述

數(shù)據(jù)建模

使用pls軟件包以最大/大量的主成分運(yùn)行PCR。
?

1. 


2. # 使用pls軟件包，以最大/較大的成分?jǐn)?shù)運(yùn)行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. # 分?jǐn)?shù)。

3. scre(mod)

?

1. #負(fù)荷

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

?

我們選擇 3 個(gè)主成分：?

1. 


2. # 我們選擇4個(gè)成分

3. m <- ?ncmp = 3, data = mrs_TAI vdon = "LOO", akknie = RUE

4. 


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

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. # 讓我們也繪制一下殘差與每個(gè)輸入X的關(guān)系。

3. 


4. for ( i in 2:11){

5. plot(res~masAN[,i],type="n")

6. }

7. 


?

解釋/結(jié)論

現(xiàn)在讓我們看一下結(jié)果——“解釋/結(jié)論”：

1. 


2. # 現(xiàn)在我們來看看結(jié)果 - 4) "解釋/結(jié)論"

3. par(mfrw = c(2, 2))

4. # 繪制具有Jacknife不確定性的系數(shù)。

5. obfi <- red(mod3,, wich = "vltn)

6. abe(lm(ft[,2] ~ fit[,1])

7. plt(mo3, ses = TUE,)

1. # 最后是一些輸出

2. test(mo3, nm = 3)

預(yù)測

1. 


2. # 現(xiàn)在讓我們試著預(yù)測TEST集的4個(gè)數(shù)據(jù)點(diǎn)。

3. prdit(md3, nwaa =TEST)

4. plt(TEST$lgg, pes)

5. 


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

2. rmsep

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應(yīng)用案例

2.面板平滑轉(zhuǎn)移回歸(PSTR)分析案例實(shí)現(xiàn)

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

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

5.R語言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)

6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實(shí)現(xiàn)

7.在R語言中實(shí)現(xiàn)Logistic邏輯回歸

8.python用線性回歸預(yù)測股票價(jià)格

9.R語言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)