這個練習的目的是根據微陣列實驗中測量的200個基因的表達水平預測TRIM32的表達水平。為此，需要從構建中心化數據開始。我們將其存儲在兩個矩陣X和Y中。

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

Y <- scale(tri, center = TRUE)

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

1.2 奇異性詛咒

我們首先假設預測因子和結果已經中心化，因此截距為0。我們會看到通常的回歸模型。

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

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

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

#> [1] 120 200

qr(X)$rank

#> [1] 119

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

qr(XtX)$rank

#> [1] 119

# 嘗試用solve進行求解。

solve(XtX)

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

2 主成分回歸

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

2.1對主成分PC的經典線性回歸

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

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

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

Zk <- pca$x[, 1:k] # 分數矩陣

# 在經典的線性回歸中使用這些分數

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

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

2.2 使用軟體包

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

• 要使用的成分（PC）的數量是通過參數ncomp來確定
• 該函數允許你首先對預測因子進行標準化（set scale = TRUE）和中心化（set center = TRUE）（在這裡的例子中，XX已經被中心化和標準化了）。

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

#X已經被標準化和中心化了

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

在後面的階段，我們將研究如何選擇預測誤差最小的成分數。

3 嶺回歸、Lasso 和彈性網Elastic Nets

嶺回歸、Lasso 回歸和彈性網Elastic Nets都是密切相關的技術，基於同樣的想法：在估計函數中加入一個懲罰項，使(XTX)再次成為滿秩，並且是可逆的。可以使用兩種不同的懲罰項或正則化方法。

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

• L2正則化：這種正則化在估計方程中增加了一個項γ2‖β‖22。這個懲罰項是基於係數大小的平方。這被嶺回歸所使用。

彈性網結合了兩種類型的正則化。它是通過引入一個α混合參數來實現的，該參數本質上是將L1和L2規範結合在一個加權平均中。

4 練習:嶺回歸的驗證

在最小平方回歸中，估計函數的最小化 可以得到解。

對於嶺回歸所使用的懲罰性最小二乘法準則，你要最小化，可以得到解。

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

脊參數γ將係數縮減為0，γ=0相當於OLS（無縮減），γ=+∞相當於將所有β^設置為0。最佳參數位於兩者之間，需要由用戶進行調整。

習題

使用R解決以下練習。

gamma <- 2 #

# 計算懲罰矩陣

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

dim(XtX_gammaI)

#> [1] 200 200

qr(XtX_gammaI)$rank == 200 #

#> [1] TRUE

向下滑動查看結果

# 是的，可以被計算。

XtX_gammaI_inv <- solve(XtX_gammaI)

向下滑動查看結果

## 計算嶺β估計值

## 使用`drop`來刪除維度並創建向量

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

#> [1] 200

我們現在已經手動計算了嶺回歸的估計值。

向下滑動查看結果

5 用glmnet進行嶺回歸和套索lasso回歸

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

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

示範：嶺回歸

讓我們進行嶺回歸，以便用200個基因探針數據預測TRIM32基因的表達水平。我們可以從使用γ值為2開始。

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

#看一下前10個係數

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

我們創建一個γ值的網格，也就是作為glmnet函數的輸入值的範圍。請注意，這個函數的lambda參數可以採用一個值的向量作為輸入，允許用相同的輸入數據但不同的超參數來擬合多個模型。

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

# 繪製係數與對數 lambda序列的對比圖!

plot(ridge_mod_grid)

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

這張圖被稱為係數曲線圖，每條彩線代表回歸模型中的一個係數β^，並顯示它們如何隨著γ（對數）1值的增加而變化。

點擊標題查閱往期內容

r語言中對LASSO回歸，Ridge嶺回歸和彈性網絡Elastic Net模型實現

左右滑動查看更多

01

02

03

04

請注意，對於更高的γ值，係數估計值變得更接近於0，顯示了嶺懲罰的收縮效應。

與PC回歸的例子類似，我們相當隨意地選擇了γ=2和網格。我們隨後會看到，如何選擇γ，使預測誤差最小。

6 練習: Lasso 回歸

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

任務

你不必在這裡提供一個自定義的γ（lambda）值序列，而是可以依靠glmnet的默認行為，即根據數據選擇γ值的網格。

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

# 默認情況下，它使用100個lambda值的序列

向下滑動查看結果

plot(lasso_model

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

向下滑動查看結果

7 預測模型的評估和超參數的調整

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

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

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

set.seed(1)

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

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

# 訓練數據

trainX <- X[trainID, ]

trainY <- Y[trainID]

# 測試數據

testX <- X[-trainID, ]

testY <- Y[-trainID]

為了使以後的模型擬合更容易一些，我們還將創建2個數據框，將訓練和測試數據的因變量和預測因素結合起來。

7.1 模型評估

我們對我們的模型的樣本外誤差感興趣，即我們的模型在未見過的數據上的表現如何。 這將使我們能夠比較不同類別的模型。對於連續結果，我們將使用平均平方誤差（MSE）（或其平方根版本，RMSE）。

該評估使我們能夠在數據上比較不同類型模型的性能，例如PC主成分回歸、嶺回歸和套索lasso回歸。然而，我們仍然需要通過選擇最佳的超參數（PC回歸的PC數和lasso和山脊的γ數）來找到這些類別中的最佳模型。為此，我們將在訓練集上使用k-fold交叉驗證。

7.2 調整超參數

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

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

8 例子: PC回歸的評估

我們從PC回歸開始，使用k-fold交叉驗證尋找使MSE最小的最佳PC數。然後，我們使用這個最優的PC數來訓練最終模型，並在測試數據上對其進行評估。

8.1 用k-fold交叉驗證來調整主成分的數量

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

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

## 為可重複性設置種子，kCV是一個隨機的過程!

set.seed(123)

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

summary(pcr_cv)

我們可以繪製每個成分數量的預測均方根誤差（RMSEP），如下所示。

plot(pcr_cv, plottype = "validation")

選擇最佳的成分數。這裡我們使用 "one-sigma "方法，它返回RMSE在絕對最小值的一個標準誤差內的最低成分數。

plot(pcr, method = "onesigma")

這個結果告訴我們，我們模型的最佳成分數是13。

8.2 對測試數據進行驗證

我們現在使用最佳成分數來訓練最終的PCR模型。然後通過對測試數據進行預測並計算MSE來驗證這個模型。

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

#平均平方誤差

## obs: 觀測值; pred: 預測值

MSE <- function(obs, pred)

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

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

最後，我們將我們的因變量（TRIM32基因表達）的預測值與我們測試集的實際觀察值進行對比。

plot(pcr_model, line = TRUE)

9 練習：評估和比較預測模型

lasso_cv

#>

請注意，我們可以從CV結果中提取擬合的 lasso回歸對象，並像以前一樣製作係數曲線圖。

我們可以尋找能產生最佳效果的伽瑪值。這裡有兩種可能性。

• lambda.min: 給出交叉驗證最佳結果的γ值。
• lambda.1se:γ的最大值，使MSE在交叉驗證的最佳結果的1個標準誤差之內。

我們將在這裡使用lambda.min來擬合最終模型，並在測試數據上生成預測。請注意，我們實際上不需要重新進行擬合，我們只需要使用我們現有的lasso_cv對象，它已經包含了lambda值範圍的擬合模型。我們可以使用predict函數並指定s參數（在這種情況下設置lambda）來對測試數據進行預測。

向下滑動查看結果

請注意，我們可以從CV結果中提取擬合的嶺回歸對象，並製作係數曲線圖。

我們可以尋找能產生最佳效果的伽瑪值。這裡有兩種可能性。

• lambda.min: 給出交叉驗證最佳結果的γ值。
• lambda.1se: γ的最大值，使MSE在交叉驗證的最佳結果的1個標準誤差之內。

我們在這裡使用lambda.min來擬合最終的模型並在測試數據上生成預測。請注意，我們實際上不需要重新進行擬合，我們只需要使用我們現有的ridge_cv對象，它已經包含了lambda值範圍的擬合模型。我們可以使用predict函數並指定s參數（在這種情況下混亂地設置lambda）來對測試數據進行預測。

ridge_preds <- predict

##計算MSE

向下滑動查看結果

模型MSEPCR0.3655052Lasso0.3754368Ridge0.3066121向下滑動查看結果

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

本文摘選 《 R語言高維數據懲罰回歸方法：主成分回歸PCR、嶺回歸、lasso、彈性網絡elastic net分析基因數據（含練習題） 》 ，點擊「閱讀原文」獲取全文完整資料。

點擊標題查閱往期內容

Python高維變量選擇:SCAD平滑剪切絕對偏差懲罰、Lasso懲罰函數比較

R語言懲罰logistic邏輯回歸（LASSO,嶺回歸）高維變量選擇的分類模型案例

R使用LASSO回歸預測股票收益

廣義線性模型glm泊松回歸的lasso、彈性網絡分類預測學生考試成績數據和交叉驗證

貝葉斯分位數回歸、lasso和自適應lasso貝葉斯分位數回歸分析免疫球蛋白、前列腺癌數據

R語言RSTAN MCMC：NUTS採樣算法用LASSO 構建貝葉斯線性回歸模型分析職業聲望數據

r語言中對LASSO回歸，Ridge嶺回歸和彈性網絡Elastic Net模型實現

R語言高維數據懲罰回歸方法：主成分回歸PCR、嶺回歸、lasso、彈性網絡elastic net分析基因數據（含練習題）

Python中LARS和Lasso回歸之最小角算法Lars分析波士頓住房數據實例

R語言Bootstrap的嶺回歸和自適應LASSO回歸可視化

R語言Lasso回歸模型變量選擇和糖尿病發展預測模型R語言實現貝葉斯分位數回歸、lasso和自適應lasso貝葉斯分位數回歸分析

基於R語言實現LASSO回歸分析

R語言用LASSO，adaptive LASSO預測通貨膨脹時間序列

R語言自適應LASSO 多項式回歸、二元邏輯回歸和嶺回歸應用分析

R語言懲罰logistic邏輯回歸（LASSO,嶺回歸）高維變量選擇的分類模型案例

Python中的Lasso回歸之最小角算法LARS

r語言中對LASSO回歸，Ridge嶺回歸和彈性網絡Elastic Net模型實現

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

R語言實現LASSO回歸——自己編寫LASSO回歸算法

R使用LASSO回歸預測股票收益

python使用LASSO回歸預測股票收益Python中LARS和Lasso回歸之最小角算法Lars分析波士頓住房數據實例

R語言Bootstrap的嶺回歸和自適應LASSO回歸可視化

R語言Lasso回歸模型變量選擇和糖尿病發展預測模型R語言實現貝葉斯分位數回歸、lasso和自適應lasso貝葉斯分位數回歸分析

基於R語言實現LASSO回歸分析

R語言用LASSO，adaptive LASSO預測通貨膨脹時間序列

R語言自適應LASSO 多項式回歸、二元邏輯回歸和嶺回歸應用分析

R語言懲罰logistic邏輯回歸（LASSO,嶺回歸）高維變量選擇的分類模型案例

Python中的Lasso回歸之最小角算法LARS

r語言中對LASSO回歸，Ridge嶺回歸和彈性網絡Elastic Net模型實現

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

R語言實現LASSO回歸——自己編寫LASSO回歸算法

R使用LASSO回歸預測股票收益

python使用LASSO回歸預測股票收益