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

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

1 介紹

在本文中，我們將研究以下主題

• 證明為什么低維預(yù)測模型在高維中會失敗。

• 進(jìn)行主成分回歸（PCR）。

• 使用glmnet()進(jìn)行嶺回歸、lasso 和彈性網(wǎng)elastic net

• 對這些預(yù)測模型進(jìn)行評估

1.1?數(shù)據(jù)集

在本文中，我們將使用基因表達(dá)數(shù)據(jù)。這個數(shù)據(jù)集包含120個樣本的200個基因的基因表達(dá)數(shù)據(jù)。這些數(shù)據(jù)來源于哺乳動物眼組織樣本的微陣列實驗。

該數(shù)據(jù)集由兩個對象組成：

• genes: 一個120×200的矩陣，包含120個樣本（行）的200個基因的表達(dá)水平（列）。

• trim32: 一個含有120個TRIM32基因表達(dá)水平的向量。

1. ##查看剛剛加載的對象

2. str(genes)

?

這個練習(xí)的目的是根據(jù)微陣列實驗中測量的200個基因的表達(dá)水平預(yù)測TRIM32的表達(dá)水平。為此，需要從構(gòu)建中心化數(shù)據(jù)開始。我們將其存儲在兩個矩陣X和Y中。?

1. X <- scale(gen, center = TRUE, scale = TRUE)

2. Y <- scale(tri, center = TRUE)

請記住，標(biāo)準(zhǔn)化可以避免量綱上的差異，使一個變量（基因）在結(jié)果中具有更大的影響力。對于Y向量，這不是一個問題，因為我們討論的是一個單一的變量。不進(jìn)行標(biāo)準(zhǔn)化會使預(yù)測結(jié)果可解釋為 "偏離平均值"。

1.2?奇異性詛咒

我們首先假設(shè)預(yù)測因子和結(jié)果已經(jīng)中心化，因此截距為0。我們會看到通常的回歸模型。

我們的目標(biāo)是得到β的最小二乘估計值，由以下公式給出

其中p×p矩陣(XTX)-1是關(guān)鍵! 為了能夠計算出XTX的逆，它必須是滿秩p。我們檢查一下。?

1. dim(X) # 120 x 200, ?p > n!

2. #> [1] 120 200

3. qr(X)$rank

4. #> [1] 119

5. 


6. XtX <- crossprod(X) # 更有效地計算t(X) %*% X

7. qr(XtX)$rank

8. #> [1] 119

9. 


10. # ?嘗試用solve進(jìn)行求解。

11. solve(XtX)

??

我們意識到無法計算(XTX)-1，因為(XTX)的秩小于p，因此我們無法通過最小二乘法得到β^! 這通常被稱為奇異性問題。

2?主成分回歸

處理這種奇異性的第一個方法是使用主成分繞過它。由于min(n,p)=n=120，PCA將得到120個成分，每個成分是p=200個變量的線性組合。這120個PC包含了原始數(shù)據(jù)中的所有信息。我們也可以使用X的近似值，即只使用幾個（k<120）PC。因此，我們使用PCA作為減少維度的方法，同時盡可能多地保留觀測值之間的變化。一旦我們有了這些PC，我們就可以把它們作為線性回歸模型的變量。

2.1對主成分PC的經(jīng)典線性回歸

我們首先用prcomp計算數(shù)據(jù)的PCA。我們將使用一個任意的k=4個PC的截止點來說明對PC進(jìn)行回歸的過程。

1. k <- 4 #任意選擇k=4

2. Vk <- pca$rotation[, 1:k] # 載荷矩陣

3. Zk <- pca$x[, 1:k] # 分?jǐn)?shù)矩陣

4. 


5. # 在經(jīng)典的線性回歸中使用這些分?jǐn)?shù)

?由于X和Y是中心化的，截距近似為0。

輸出結(jié)果顯示，PC1和PC4的β估計值與0相差很大（在p<0.05），但是結(jié)果不能輕易解釋，因為我們沒有對PC的直接解釋。

2.2 使用軟件包

PCR也可以直接在數(shù)據(jù)上進(jìn)行（所以不必先手動進(jìn)行PCA）。在使用這個函數(shù)時，你必須牢記幾件事。

1. 要使用的成分（PC）的數(shù)量是通過參數(shù)ncomp來確定

2. 該函數(shù)允許你首先對預(yù)測因子進(jìn)行標(biāo)準(zhǔn)化（set scale = TRUE）和中心化（set center = TRUE）（在這里的例子中，XX已經(jīng)被中心化和標(biāo)準(zhǔn)化了）。

你可以用與使用lm()相同的方式使用pcr()函數(shù)。使用函數(shù)summary()可以很容易地檢查得出的擬合結(jié)果，但輸出結(jié)果看起來與你從lm得到的結(jié)果完全不同。

#X已經(jīng)被標(biāo)準(zhǔn)化和中心化了

??

首先，輸出顯示了數(shù)據(jù)維度和使用的擬合方法。在本例中，是基于SVD的主成分PC計算。summary()函數(shù)還提供了使用不同數(shù)量的成分在預(yù)測因子和響應(yīng)中解釋方差的百分比。例如，第一個PC只解釋了所有方差的61.22%，或預(yù)測因子中的信息，它解釋了結(jié)果中方差的62.9%。請注意，對于這兩種方法，主成分?jǐn)?shù)量的選擇都是任意選擇的，即4個。

在后面的階段，我們將研究如何選擇預(yù)測誤差最小的成分?jǐn)?shù)。

3 嶺回歸、Lasso?和彈性網(wǎng)Elastic Nets

嶺回歸、Lasso?回歸和彈性網(wǎng)Elastic Nets都是密切相關(guān)的技術(shù)，基于同樣的想法：在估計函數(shù)中加入一個懲罰項，使(XTX)再次成為滿秩，并且是可逆的?？梢允褂脙煞N不同的懲罰項或正則化方法。

1. L1正則化：這種正則化在估計方程中加入一個γ1‖β‖1。該項將增加一個基于系數(shù)大小絕對值的懲罰。這被Lasso回歸所使用。

1. L2正則化：這種正則化在估計方程中增加了一個項γ2‖β‖22。這個懲罰項是基于系數(shù)大小的平方。這被嶺回歸所使用。

彈性網(wǎng)結(jié)合了兩種類型的正則化。它是通過引入一個α混合參數(shù)來實現(xiàn)的，該參數(shù)本質(zhì)上是將L1和L2規(guī)范結(jié)合在一個加權(quán)平均中。?

4 練習(xí):嶺回歸的驗證

在最小平方回歸中，估計函數(shù)的最小化

?可以得到解

。

對于嶺回歸所使用的懲罰性最小二乘法準(zhǔn)則，你要最小化

，可以得到解

。?

其中II是p×p的識別矩陣。

脊參數(shù)γ將系數(shù)縮減為0，γ=0相當(dāng)于OLS（無縮減），γ=+∞相當(dāng)于將所有β^設(shè)置為0。最佳參數(shù)位于兩者之間，需要由用戶進(jìn)行調(diào)整。

習(xí)題

使用R解決以下練習(xí)。

1. 驗證

?秩為200,對于任何一個?

?.

1. 


2. gamma <- 2 #

3. 


4. # 計算懲罰矩陣

5. XtX_gammaI <- XtX + (gamma * diag(p))

6. dim(XtX_gammaI)

7. #> [1] 200 200

8. qr(XtX_gammaI)$rank == 200 #

9. #> [1] TRUE

2. 檢查

的逆值是否可以計算出來。

1. # 是的，可以被計算。

2. XtX_gammaI_inv <- solve(XtX_gammaI)

3. 最后，計算

。

1. ## 計算嶺β估計值

2. ## 使用`drop`來刪除維度并創(chuàng)建向量

3. length(ridge_betas) # 每個基因都有一個

4. #> [1] 200

?我們現(xiàn)在已經(jīng)手動計算了嶺回歸的估計值。?

5?用glmnet進(jìn)行嶺回歸和套索lasso回歸

glmnet允許你擬合所有三種類型的回歸。使用哪種類型，可以通過指定alpha參數(shù)來決定。對于嶺回歸，你將alpha設(shè)置為0，而對于套索lasso回歸，你將alpha設(shè)置為1。其他介于0和1之間的α值將適合一種彈性網(wǎng)的形式。這個函數(shù)的語法與其他的模型擬合函數(shù)略有不同。你必須傳遞一個x矩陣以及一個y向量。

控制懲罰 "強度 "的gamma值可以通過參數(shù)lambda傳遞。函數(shù)glmnet()還可以進(jìn)行搜索，來找到最佳的擬合伽馬值。這可以通過向參數(shù)lambda傳遞多個值來實現(xiàn)。如果不提供，glmnet將根據(jù)數(shù)據(jù)自己生成一個數(shù)值范圍，而數(shù)值的數(shù)量可以用nlambda參數(shù)控制。這通常是使用glmnet的推薦方式，詳見glmnet。?

示范：嶺回歸?

讓我們進(jìn)行嶺回歸，以便用200個基因探針數(shù)據(jù)預(yù)測TRIM32基因的表達(dá)水平。我們可以從使用γ值為2開始。

1. glmnet(X, Y, alpha = 0, lambda = gamma)

2. 


3. #看一下前10個系數(shù)

第一個系數(shù)是截距，基本上也是0。但γ的值為2可能不是最好的選擇，所以讓我們看看系數(shù)在γ的不同值下如何變化。

我們創(chuàng)建一個γ值的網(wǎng)格，也就是作為glmnet函數(shù)的輸入值的范圍。請注意，這個函數(shù)的lambda參數(shù)可以采用一個值的向量作為輸入，允許用相同的輸入數(shù)據(jù)但不同的超參數(shù)來擬合多個模型。?

1. grid <- seq(1, 1000, by = 10) ?# 1到1000，步驟為10

2. # 繪制系數(shù)與對數(shù) lambda序列的對比圖!

3. plot(ridge_mod_grid)

4. # 在gamma = 2處添加一條垂直線

這張圖被稱為系數(shù)曲線圖，每條彩線代表回歸模型中的一個系數(shù)β^，并顯示它們?nèi)绾坞S著γ（對數(shù)）1值的增加而變化。

請注意，對于更高的γ值，系數(shù)估計值變得更接近于0，顯示了嶺懲罰的收縮效應(yīng)。

與PC回歸的例子類似，我們相當(dāng)隨意地選擇了γ=2和網(wǎng)格。我們隨后會看到，如何選擇γ，使預(yù)測誤差最小。

6?練習(xí): Lasso 回歸

Lasso 回歸也是懲罰性回歸的一種形式，但我們沒有像最小二乘法和嶺回歸那樣的β^的分析解。為了擬合一個Lasso 模型，我們再次使用glmnet()函數(shù)。然而，這一次我們使用的參數(shù)是α=1

任務(wù)

1. 驗證設(shè)置α=1確實對應(yīng)于使用第3節(jié)的方程進(jìn)行套索回歸。

2. 用glmnet函數(shù)進(jìn)行Lasso 套索回歸，Y為因變量，X為預(yù)測因子。

你不必在這里提供一個自定義的γ（lambda）值序列，而是可以依靠glmnet的默認(rèn)行為，即根據(jù)數(shù)據(jù)選擇γ值的網(wǎng)格。

1. # 請注意，glmnet()函數(shù)可以自動提供伽馬值

2. # 默認(rèn)情況下，它使用100個lambda值的序列

3. 繪制系數(shù)曲線圖并進(jìn)行解釋。

plot(lasso_model

請注意，非零系數(shù)的數(shù)量顯示在圖的頂部。在lasso回歸的情況下，與嶺回歸相比，正則化要不那么平滑，一些系數(shù)在較高的γ值下會增加，然后急劇下降到0。與嶺回歸相反，lasso最終將所有系數(shù)縮減為0。?

7?預(yù)測模型的評估和超參數(shù)的調(diào)整

首先，我們將把我們的原始數(shù)據(jù)分成訓(xùn)練集和測試集來驗證我們的模型。訓(xùn)練集將被用來訓(xùn)練模型和調(diào)整超參數(shù)，而測試集將被用來評估我們最終模型的樣本外性能。如果我們使用相同的數(shù)據(jù)來擬合和測試模型，我們會得到有偏見的結(jié)果。

在開始之前，我們使用set.seed()函數(shù)來為R的隨機數(shù)生成器設(shè)置一個種子，這樣我們就能得到與下面所示完全相同的結(jié)果。一般來說，在進(jìn)行交叉驗證等包含隨機性元素的分析時，設(shè)置一個隨機種子是很好的做法，這樣所得到的結(jié)果就可以在以后的時間里重現(xiàn)。

我們首先使用sample()函數(shù)將樣本集分成兩個子集，從原來的120個觀測值中隨機選擇80個觀測值的子集。我們把這些觀測值稱為訓(xùn)練集。其余的觀察值將被用作測試集。

1. set.seed(1)

2. # 從X的行中隨機抽取80個ID（共120個）。

3. trainID <- sample(nrow(X), 80)

5. # 訓(xùn)練數(shù)據(jù)

6. trainX <- X[trainID, ]

7. trainY <- Y[trainID]

9. # 測試數(shù)據(jù)

10. testX <- X[-trainID, ]

11. testY <- Y[-trainID]

為了使以后的模型擬合更容易一些，我們還將創(chuàng)建2個數(shù)據(jù)框，將訓(xùn)練和測試數(shù)據(jù)的因變量和預(yù)測因素結(jié)合起來。?

7.1?模型評估

我們對我們的模型的樣本外誤差感興趣，即我們的模型在未見過的數(shù)據(jù)上的表現(xiàn)如何。?這將使我們能夠比較不同類別的模型。對于連續(xù)結(jié)果，我們將使用平均平方誤差（MSE）（或其平方根版本，RMSE）。

該評估使我們能夠在數(shù)據(jù)上比較不同類型模型的性能，例如PC主成分回歸、嶺回歸和套索lasso回歸。然而，我們?nèi)匀恍枰ㄟ^選擇最佳的超參數(shù)（PC回歸的PC數(shù)和lasso和山脊的γ數(shù)）來找到這些類別中的最佳模型。為此，我們將在訓(xùn)練集上使用k-fold交叉驗證。

7.2?調(diào)整超參數(shù)

測試集只用于評估最終模型。為了實現(xiàn)這個最終模型，我們需要找到最佳的超參數(shù)，即對未見過的數(shù)據(jù)最能概括模型的超參數(shù)。我們可以通過在訓(xùn)練數(shù)據(jù)上使用k倍交叉驗證（CVk）來估計這一點。

對于任何廣義線性模型，CVk估計值都可以用cv.glm()函數(shù)自動計算出來。

8?例子: PC回歸的評估

我們從PC回歸開始，使用k-fold交叉驗證尋找使MSE最小的最佳PC數(shù)。然后，我們使用這個最優(yōu)的PC數(shù)來訓(xùn)練最終模型，并在測試數(shù)據(jù)上對其進(jìn)行評估。

8.1?用k-fold交叉驗證來調(diào)整主成分的數(shù)量

方便的是，pcr函數(shù)有一個k-fold交叉驗證的實現(xiàn)。我們只需要設(shè)置validation = CV和segments = 20就可以用PC回歸進(jìn)行20折交叉驗證。如果我們不指定ncomp，pcr將選擇可用于CV的最大數(shù)量的PC。

請注意，我們的訓(xùn)練數(shù)據(jù)trainX由80個觀測值（行）組成。如果我們執(zhí)行20折的CV，這意味著我們將把數(shù)據(jù)分成20組，所以每組由4個觀測值組成。在每個CV周期中，有一個組將被排除，模型將在剩余的組上進(jìn)行訓(xùn)練。這使得我們在每個CV周期有76個訓(xùn)練觀測值，所以可以用于線性回歸的最大成分?jǐn)?shù)是75。

1. ## 為可重復(fù)性設(shè)置種子，kCV是一個隨機的過程!

2. set.seed(123)

4. ##Y ~ . "符號的意思是：用數(shù)據(jù)中的每個其他變量來擬合Y。

5. summary(pcr_cv)

我們可以繪制每個成分?jǐn)?shù)量的預(yù)測均方根誤差（RMSEP），如下所示。

plot(pcr_cv, plottype = "validation")

?

?

?

選擇最佳的成分?jǐn)?shù)。這里我們使用 "one-sigma "方法，它返回RMSE在絕對最小值的一個標(biāo)準(zhǔn)誤差內(nèi)的最低成分?jǐn)?shù)。

plot(pcr, method = "onesigma")

這個結(jié)果告訴我們，我們模型的最佳成分?jǐn)?shù)是13。

8.2?對測試數(shù)據(jù)進(jìn)行驗證

我們現(xiàn)在使用最佳成分?jǐn)?shù)來訓(xùn)練最終的PCR模型。然后通過對測試數(shù)據(jù)進(jìn)行預(yù)測并計算MSE來驗證這個模型。

我們定義了一個自定義函數(shù)來計算MSE。請注意，可以一次性完成預(yù)測和MSE計算。但是我們自己的函數(shù)在后面的lasso和ridge嶺回歸中會派上用場。

1. #平均平方誤差

2. ## obs: 觀測值; pred: 預(yù)測值

3. MSE <- function(obs, pred)

1. pcr_preds <- predict(model, newdata = test_data, ncomp = optimal_ncomp)

2. 
   

這個值本身并不能告訴我們什么，但是我們可以用它來比較我們的PCR模型和其他類型的模型。

最后，我們將我們的因變量（TRIM32基因表達(dá)）的預(yù)測值與我們測試集的實際觀察值進(jìn)行對比。

plot(pcr_model, line = TRUE)

9 練習(xí)：評估和比較預(yù)測模型

1. 對訓(xùn)練數(shù)據(jù)（trainX, trainY）進(jìn)行20次交叉驗證的lasso回歸。繪制結(jié)果并選擇最佳的λ（γ）參數(shù)。用選定的lambda擬合一個最終模型，并在測試數(shù)據(jù)上驗證它。

1. 


2. lasso_cv

3. #>

請注意，我們可以從CV結(jié)果中提取擬合的 lasso回歸對象，并像以前一樣制作系數(shù)曲線圖。?

我們可以尋找能產(chǎn)生最佳效果的伽瑪值。這里有兩種可能性。

1. lambda.min: 給出交叉驗證最佳結(jié)果的γ值。

2. lambda.1se:γ的最大值，使MSE在交叉驗證的最佳結(jié)果的1個標(biāo)準(zhǔn)誤差之內(nèi)。

我們將在這里使用lambda.min來擬合最終模型，并在測試數(shù)據(jù)上生成預(yù)測。請注意，我們實際上不需要重新進(jìn)行擬合，我們只需要使用我們現(xiàn)有的lasso_cv對象，它已經(jīng)包含了lambda值范圍的擬合模型。我們可以使用predict函數(shù)并指定s參數(shù)（在這種情況下設(shè)置lambda）來對測試數(shù)據(jù)進(jìn)行預(yù)測。?

2. 對嶺回歸做同樣的處理。

請注意，我們可以從CV結(jié)果中提取擬合的嶺回歸對象，并制作系數(shù)曲線圖。

我們可以尋找能產(chǎn)生最佳效果的伽瑪值。這里有兩種可能性。

1. lambda.min: 給出交叉驗證最佳結(jié)果的γ值。

2. lambda.1se: γ的最大值，使MSE在交叉驗證的最佳結(jié)果的1個標(biāo)準(zhǔn)誤差之內(nèi)。

我們在這里使用lambda.min來擬合最終的模型并在測試數(shù)據(jù)上生成預(yù)測。請注意，我們實際上不需要重新進(jìn)行擬合，我們只需要使用我們現(xiàn)有的ridge_cv對象，它已經(jīng)包含了lambda值范圍的擬合模型。我們可以使用predict函數(shù)并指定s參數(shù)（在這種情況下混亂地設(shè)置lambda）來對測試數(shù)據(jù)進(jìn)行預(yù)測。

1. ridge_preds <- predict

2. ##計算MSE

3. 在所考慮的模型（PCR、lasso、嶺回歸）中，哪一個表現(xiàn)最好？

模型MSEPCR0.3655052Lasso0.3754368Ridge0.3066121

1. 注意：R中的log()默認(rèn)是自然對數(shù)（以e為底），我們也會在文本中使用這個符號（比如上面圖中的x軸標(biāo)題）。這可能與你所習(xí)慣的符號（ln()）不同。要在R中取不同基數(shù)的對數(shù)，你可以指定log的基數(shù)=參數(shù)，或者使用函數(shù)log10(x)和log2(x)分別代表基數(shù)10和2?

最受歡迎的見解?

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

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

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

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

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

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

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

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

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