高维数据惩罚回归方法：主成分回归PCR岭回归lasso弹性网络elastic net分析基因数据|附代码数据

Posted 大数据部落

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了高维数据惩罚回归方法：主成分回归PCR岭回归lasso弹性网络elastic net分析基因数据|附代码数据相关的知识，希望对你有一定的参考价值。

全文链接：http://tecdat.cn/?p=23378

最近我们被客户要求撰写关于高维数据惩罚回归方法的研究报告，包括一些图形和统计输出。

在本文中，我们将使用基因表达数据。这个数据集包含120个样本的200个基因的基因表达数据。这些数据来源于哺乳动物眼组织样本的微阵列实验

1 介绍

在本文中，我们将研究以下主题

  • 证明为什么低维预测模型在高维中会失败。
  • 进行主成分回归（PCR）。
  • 使用glmnet()进行岭回归、lasso 和弹性网elastic net
  • 对这些预测模型进行评估

1.1 数据集

在本文中，我们将使用基因表达数据。这个数据集包含120个样本的200个基因的基因表达数据。这些数据来源于哺乳动物眼组织样本的微阵列实验。

该数据集由两个对象组成：

  • genes: 一个120×200的矩阵，包含120个样本（行）的200个基因的表达水平（列）。
  • trim32: 一个含有120个TRIM32基因表达水平的向量。
##查看刚刚加载的对象
str(genes)

这个练习的目的是根据微阵列实验中测量的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解决以下练习。

  1. 验证 秩为200,对于任何一个  .

gamma <- 2 # 

# 计算惩罚矩阵
XtX_gammaI <- XtX + (gamma * diag(p))
dim(XtX_gammaI)
#> [1] 200 200
qr(XtX_gammaI)$rank == 200 # 
#> [1] TRUE

向下滑动查看结果▼

  1. 检查的逆值是否可以计算出来。

# 是的，可以被计算。
XtX_gammaI_inv <- solve(XtX_gammaI)

向下滑动查看结果▼

  1. 最后，计算。

## 计算岭β估计值
## 使用`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

任务

  1. 验证设置α=1确实对应于使用第3节的方程进行套索回归。

  2. 用glmnet函数进行Lasso 套索回归，Y为因变量，X为预测因子。

你不必在这里提供一个自定义的γ（lambda）值序列，而是可以依靠glmnet的默认行为，即根据数据选择γ值的网格。

# 请注意，glmnet()函数可以自动提供伽马值
# 默认情况下，它使用100个lambda值的序列

向下滑动查看结果▼

  1. 绘制系数曲线图并进行解释。

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 练习：评估和比较预测模型

  1. 对训练数据（trainX, trainY）进行20次交叉验证的lasso回归。绘制结果并选择最佳的λ（γ）参数。用选定的lambda拟合一个最终模型，并在测试数据上验证它。

lasso_cv
#>

请注意，我们可以从CV结果中提取拟合的 lasso回归对象，并像以前一样制作系数曲线图。 

我们可以寻找能产生最佳效果的伽玛值。这里有两种可能性。

  • lambda.min: 给出交叉验证最佳结果的γ值。
  • lambda.1se:γ的最大值，使MSE在交叉验证的最佳结果的1个标准误差之内。

我们将在这里使用lambda.min来拟合最终模型，并在测试数据上生成预测。请注意，我们实际上不需要重新进行拟合，我们只需要使用我们现有的lasso_cv对象，它已经包含了lambda值范围的拟合模型。我们可以使用predict函数并指定s参数（在这种情况下设置lambda）来对测试数据进行预测。 

向下滑动查看结果▼

  1. 对岭回归做同样的处理。

请注意，我们可以从CV结果中提取拟合的岭回归对象，并制作系数曲线图。

我们可以寻找能产生最佳效果的伽玛值。这里有两种可能性。

  • lambda.min: 给出交叉验证最佳结果的γ值。
  • lambda.1se: γ的最大值，使MSE在交叉验证的最佳结果的1个标准误差之内。

我们在这里使用lambda.min来拟合最终的模型并在测试数据上生成预测。请注意，我们实际上不需要重新进行拟合，我们只需要使用我们现有的ridge_cv对象，它已经包含了lambda值范围的拟合模型。我们可以使用predict函数并指定s参数（在这种情况下混乱地设置lambda）来对测试数据进行预测。

ridge_preds <- predict
##计算MSE

向下滑动查看结果▼

  1. 在所考虑的模型（PCR、lasso、岭回归）中，哪一个表现最好？

模型MSE
PCR 0.3655052
Lasso 0.3754368
Ridge 0.3066121

向下滑动查看结果▼

  • 注意：R中的log()默认是自然对数（以e为底），我们也会在文本中使用这个符号（比如上面图中的x轴标题）。这可能与你所习惯的符号（ln()）不同。要在R中取不同基数的对数，你可以指定log的基数=参数，或者使用函数log10(x)和log2(x)分别代表基数10和2︎

回归分析06：回归参数的估计

目录

Chapter 6：回归参数的估计(4)

3.8 岭估计

3.8.1 岭估计的定义和性质

当自变量之间具有多重共线性时，岭估计是一种为了克服最小二乘估计的方差较大的问题而提出的改进的最小二乘估计方法。

岭估计的主要思想为：多重共线性下的设计矩阵 \\(X\\) 是病态的，即 \\(\\left|X\'X\\right|\\approx0\\) ，从而使得 \\(\\left(X\'X\\right)^{-1}\\) 接近奇异。为避免这一现象，给 \\(X\'X\\) 加上一个正常数对角矩阵 \\(kI\\ (k>0)\\) ，使得矩阵 \\(\\left(X\'X+kI\\right)^{-1}\\) 接近奇异的可能性要比 \\(\\left(X\'X\\right)^{-1}\\) 接近奇异的可能性小得多。因此用

\\[\\hat\\beta(k)=\\left(X\'X+kI\\right)^{-1}X\'Y \\]

作为未知参数 \\(\\beta\\) 的估计，可以得到比最小二乘估计 \\(\\hat\\beta\\) 更加稳定的估计。

岭估计：对给定的 \\(k>0\\) ，称 \\(\\hat\\beta(k)=\\left(X\'X+kI\\right)^{-1}X\'Y\\) 为回归系数 \\(\\beta\\) 的岭估计。由岭估计所建立的回归方程称为岭回归方程，称 \\(k\\) 为岭参数。对于 \\(\\hat\\beta(k)\\) 的分量 \\(\\hat\\beta_j(k)\\) ，把在平面直角坐标系中 \\(\\hat\\beta_j(k)\\) 关于 \\(k\\) 的变化所表现出来的曲线称为岭迹。

  • 岭估计 \\(\\hat\\beta(k)\\) 是一个关于 \\(k\\) 的估计类。当 \\(k=0\\) 时，\\(\\hat\\beta(k)\\) 就是通常的最小二乘估计。
  • 在进行岭估计之前，需要消除量纲的影响，故假设自变量与因变量均已标准化，故这里的所讨论的设计矩阵 \\(X\\) 均是 \\(n\\times p\\) 的矩阵。

性质 1：岭估计 \\(\\hat\\beta(k)\\) 是 \\(\\beta\\) 的有偏估计，即对 \\(\\forall k>0\\) ，\\({\\rm E}\\left[\\hat\\beta(k)\\right]\\neq\\beta\\) 。

当自变量之间存在多重共线性时，最小二乘估计虽然保持偏差部分为 \\(0\\) ，但它的方差部分却很大，最终导致它的均方误差很大。岭估计的引入就是一种牺牲无偏性，换取方差的大幅度减少，从而降低均方误差的方法。

性质 2：岭估计 \\(\\hat\\beta(k)\\) 是最小二乘估计 \\(\\hat\\beta\\) 的一个线性变换。

只需注意到

\\[\\begin{aligned} \\hat\\beta(k)&=\\left(X\'X+kI\\right)^{-1}X\'Y \\\\ \\\\ &=\\left(X\'X+kI\\right)^{-1}X\'X\\left(X\'X\\right)^{-1}XY \\\\ \\\\ &=\\left(X\'X+kI\\right)^{-1}X\'X\\hat\\beta \\ . \\end{aligned} \\]

性质 3：对任意的 \\(k>0\\) ，若 \\(\\left\\|\\hat\\beta\\right\\|\\neq0\\) ，则总有 \\(\\left\\|\\hat\\beta(k)\\right\\|<\\left\\|\\hat\\beta\\right\\|\\) 。即岭估计是把最小二乘估计 \\(\\hat\\beta\\) 向原点作适度的压缩而得到的，岭估计是一个压缩的有偏估计。

考虑多元线性回归模型 \\(Y=X\\beta+e\\) ，令

\\[Z=XP \\ , \\quad \\alpha=P\'\\beta \\ , \\]

其中 \\(P\\) 为正交矩阵，满足

\\[P\'X\'XP=\\Lambda={\\rm diag}\\left(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p\\right) \\ , \\]

这里 \\(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p>0\\) 为 \\(X\'X\\) 的特征根。将多元线性回归模型写为

\\[Y=Z\\alpha+e \\ , \\quad {\\rm E}(e)=0 \\ , \\quad {\\rm Cov}(e)=\\sigma^2I_n \\ . \\]

我们将上述模型称为线性回归模型的典则形式，称 \\(\\alpha\\) 为典则回归系数。

注意到 \\(Z\'Z=P\'X\'XP=\\Lambda\\) ，所以

\\[\\hat\\alpha=\\left(Z\'Z\\right)^{-1}Z\'Y=\\Lambda^{-1}Z\'Y \\ . \\]

而又因为

\\[\\hat\\beta=\\left(X\'X\\right)^{-1}X\'Y=P\\Lambda^{-1}P\'X\'Y=P\\Lambda^{-1}Z\'Y=P\\hat\\alpha \\ . \\]

它们相应的岭估计为

\\[\\begin{aligned} \\hat\\alpha(k)&=\\left(Z\'Z+kI\\right)^{-1}Z\'Y=\\left(\\Lambda+kI\\right)^{-1}Z\'Y \\ . \\\\ \\\\ \\hat\\beta(k)&=\\left(X\'X+kI\\right)^{-1}X\'Y \\\\ \\\\ &=PP\'\\left(X\'X+kI\\right)^{-1}PP\'X\'Y \\\\ \\\\ &=P\\hat\\alpha(k) \\ . \\end{aligned} \\]

因此有

\\[\\left\\|\\hat\\beta(k)\\right\\|=\\left\\|\\hat\\alpha(k)\\right\\|=\\left\\|\\left(\\Lambda+kI\\right)^{-1}\\Lambda\\hat\\alpha\\right\\|<\\left\\|\\hat\\alpha\\right\\|=\\left\\|\\hat\\beta\\right\\| \\ . \\]

容易证明，典则回归系数的最小二乘估计（或岭估计）和原回归系数的最小二乘估计（或岭估计）具有相同的均方误差：

\\[{\\rm MSE}(\\hat\\alpha)={\\rm MSE}(\\hat\\beta) \\ , \\quad {\\rm MSE}(\\hat\\alpha(k))={\\rm MSE}(\\hat\\beta(k)) \\ . \\]

定理 3.8.1 (岭估计存在定理)：存在 \\(k>0\\) ，使得在均方误差意义下，岭估计优于最小二乘估计，即

\\[{\\rm MSE}(\\hat\\beta(k))<{\\rm MSE}(\\hat\\beta) \\ . \\]

由岭估计的性质 3 可知，只需证存在 \\(k>0\\) ，使得

\\[{\\rm MSE}(\\hat\\alpha(k))<{\\rm MSE}(\\hat\\alpha) \\ . \\]

记 \\(f(k)={\\rm MSE}(\\hat\\alpha(k)),\\,k\\geq0\\) 。注意 \\(f(0)={\\rm MSE}(\\hat\\alpha)\\) 。只需证明 \\(f(k)\\) 在 \\([0,\\infty)\\) 上是连续函数且 \\(f\'(0)<0\\) ，则必存在一个较小的 \\(k>0\\) 使得上述不等式成立。

注意到

\\[\\begin{aligned} {\\rm E}\\left[\\hat\\alpha(k)\\right]&=\\left(\\Lambda+kI\\right)^{-1}Z\'{\\rm E}(Y) \\\\ \\\\ &=\\left(\\Lambda+kI\\right)^{-1}Z\'Z\\alpha \\\\ \\\\ &=\\left(\\Lambda+kI\\right)^{-1}\\Lambda\\alpha \\ , \\\\ \\\\ {\\rm Cov}\\left[\\hat\\alpha(k)\\right]&=\\sigma^2\\left(\\Lambda+kI\\right)^{-1}Z\'Z\\left(\\Lambda+kI\\right)^{-1} \\\\ \\\\ &=\\sigma^2\\left(\\Lambda+kI\\right)^{-1}\\Lambda\\left(\\Lambda+kI\\right)^{-1} \\ . \\end{aligned} \\]

所以

\\[\\begin{aligned} f(k)&={\\rm MSE}(\\hat\\alpha(k))={\\rm tr}\\left[{\\rm Cov}\\left[\\hat\\alpha(k)\\right]\\right]+\\left\\|{\\rm E}\\left[\\hat\\alpha(k)\\right]-\\alpha\\right\\|^2 \\\\ \\\\ &=\\sigma^2\\sum_{j=1}^p\\frac{\\lambda_j}{\\left(\\lambda_j+k\\right)^2}+k^2\\sum_{j=1}^p\\frac{\\alpha_j^2}{\\left(\\lambda_j+k\\right)^2} \\\\ \\\\ &\\xlongequal{def}f_1(k)+f_2(k) \\ . \\end{aligned} \\]

显然 \\(f(k)\\) 是 \\([0,\\infty)\\) 上的连续函数，又因为

\\[f_1\'(k)=-2\\sigma^2\\sum_{j=1}^p\\frac{\\lambda_j}{\\left(\\lambda_j+k\\right)^3} \\ , \\quad f_1\'(0)=-2\\sigma^2\\sum_{j=1}^p\\frac{1}{\\lambda_j^2}<0 \\ , \\]

以及

\\[f_2\'(k)=2k\\sum_{j=1}^p\\frac{\\lambda_j\\alpha_j^2}{\\left(\\lambda_j+k\\right)^3} \\ , \\quad f_2\'(0)=0 \\ , \\]

所以

\\[f\'(0)=f_1\'(0)+f_2\'(0)=-2\\sigma^2\\sum_{j=1}^p\\frac{1}{\\lambda_j^2}<0 \\ . \\]

岭估计的存在性定理从理论上证明了存在某个岭估计优于最小二乘估计，但要找出这个岭参数 \\(k\\) 是不容易的。这个解依赖于未知参数 \\(\\alpha_i,\\,i=1,2,\\cdots,p\\) 和 \\(\\sigma^2\\) ，所以不可能从解方程的角度获得。因此，我们需要提出从其他途径选择岭参数 \\(k\\) 的方法。

3.8.2 岭参数的选择方法

(1) Hoerl-Kennard 公式

Hoerl 和 Kennard 提出的选择岭参数 \\(k\\) 的公式为

\\[\\hat k=\\frac{\\hat\\sigma^2}{\\max\\limits_i\\hat\\alpha_i^2} \\ . \\]

注意到，理论上的最优岭参数是下列方程的解：令 \\(f\'(k)=0\\) ，则有

\\[f\'(k)=-2\\sigma^2\\sum_{i=1}^p\\frac{\\lambda_i}{\\left(\\lambda_i+k\\right)^3}+2k\\sum_{i=1}^p\\frac{\\lambda_i\\alpha_i^2}{\\left(\\lambda_i+k\\right)^3}=2\\sum_{i=1}^p\\frac{\\lambda_i\\left(k\\alpha_i^2-\\sigma^2\\right)}{\\left(\\lambda_i+k\\right)^3}=0 \\ . \\]

若 \\(k\\alpha_i^2-\\sigma^2<0\\) 对 \\(i=1,2,\\cdots,p\\) 都成立，则 \\(f\'(k)<0\\) ，于是取

\\[k^*=\\frac{\\hat\\sigma^2}{\\max\\limits_i\\hat\\alpha_i^2} \\ , \\]

当 \\(0<k<k^*\\) 时，\\(f\'(k)<0\\) 恒成立，因而 \\(f(k)\\) 在 \\((0,k^*)\\) 上是单调递减函数。再由 \\(f(k)\\) 是 \\((0,k^*)\\) 上的连续函数得到 \\(f(k^*)<f(0)\\) 。用 \\(\\hat\\alpha_i\\) 和 \\(\\hat\\sigma^2\\) 代替 \\(\\alpha_i\\) 和 \\(\\sigma^2\\) 即可得到我们需要的岭参数 \\(\\hat k\\) 。

(2) 岭迹法

将 \\(\\hat\\beta_1(k),\\hat\\beta_2(k),\\cdots,\\hat\\beta_p(k)\\) 的岭迹画在一张图上，根据岭迹的变化趋势选择岭参数 \\(k\\) 。以下是几条选择岭参数 \\(k\\) 的准则：

  • 各回归系数的岭估计大致趋于稳定；
  • 用最小二乘估计时符号不合理的回归系数，其岭估计的符号变得合理；
  • 残差平方和不要上升太多；

一般情况下，我们选择使得各条岭迹均开始趋于稳定的最小的 \\(k\\) 值。

3.8.3 岭估计的几何意义

前面已经证明，岭估计 \\(\\hat\\beta(k)\\) 是最小二乘估计 \\(\\hat\\beta\\) 的一种压缩。如果我们已经有了 \\(\\hat\\beta\\) ，希望将它的长度压缩到原来的 \\(c\\) 倍 \\((0<c<1)\\) ，并使残差平方和上升尽可能小，可以证明，这样的估计就是岭估计。

设 \\(b\\) 为 \\(\\beta\\) 的任意估计，对应的残差平方和为

\\[\\begin{aligned} {\\rm RSS}(b)&=\\left\\|Y-Xb\\right\\|^2 \\\\ \\\\ &=\\left\\|Y-X\\hat\\beta+X(\\hat\\beta-b)\\right\\|^2 \\\\ \\\\ &=\\left\\|Y-X\\hat\\beta\\right\\|^2+(\\hat\\beta-b)\'X\'X(\\hat\\beta-b) \\ . \\end{aligned} \\]

所以，将 \\(\\hat\\beta\\) 的长度压缩到原来的 \\(c\\) 倍，且使残差平方和上升最小，等价于求解下列极值问题：

\\[\\begin{aligned} \\min_b \\quad &(b-\\hat\\beta)\'X\'X(b-\\hat\\beta) \\ , \\\\ {\\rm s.t.}\\quad & \\|b\\|^2=\\left\\|c\\hat\\beta\\right\\|^2 \\ . \\end{aligned} \\]

设 \\(P\\) 为正交矩阵，满足

\\[P\'X\'XP=\\Lambda={\\rm diag}\\left(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p\\right) \\ , \\]

其中 \\(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p>0\\) 为 \\(X\'X\\) 的特征根。记

\\[\\alpha=P\'\\beta \\ , \\quad d=P\'b \\ , \\quad \\hat\\alpha=P\'\\hat\\beta \\ . \\]

显然上述极值问题等价于

\\[\\begin{aligned} \\min_d \\quad &(d-\\hat\\alpha)\'\\Lambda(d-\\hat\\alpha) \\ , \\\\ {\\rm s.t.}\\quad & \\|d\\|^2=\\left\\|c\\hat\\alpha\\right\\|^2 \\ . \\end{aligned} \\]

用 Lagrange 乘子法，构造辅助函数

\\[F(d,k)=(d-\\hat\\alpha)\'\\Lambda(d-\\hat\\alpha)+k\\left(d\'d-c^2\\left\\|\\hat\\alpha\\right\\|^2\\right) \\ , \\]

其中 \\(k\\ (k\\neq0)\\) 为 Lagrange 乘子。对上式关于 \\(d\\) 求导，并令其等于 \\(0\\) ，则有

\\[\\frac{\\partial F(d,k)}{\\partial d}=2(\\Lambda+kI)d-2\\Lambda\\hat\\alpha=0 \\ . \\]

解得

\\[d=(\\Lambda+kI)^{-1}\\Lambda\\hat\\alpha \\ , \\quad \\Longrightarrow \\quad b=\\left(X\'X+kI\\right)^{-1}X\'Y \\ . \\]

下面只需证 \\(k>0\\) 。将 \\(d\\) 的解析表达式代入目标函数中，记为 \\(Q(k)\\) ，于是

\\[\\begin{aligned} Q(k)&=(d-\\hat\\alpha)\'\\Lambda(d-\\hat\\alpha) \\\\ \\\\ &=\\hat\\alpha\'\\left[\\left((\\Lambda+kI)^{-1}\\Lambda-I\\right)\'\\Lambda\\left((\\Lambda+kI)^{-1}\\Lambda-I\\right)\\right]\\hat\\alpha \\\\ \\\\ &=k^2\\hat\\alpha\'{\\rm diag}\\left(\\frac{\\lambda_1}{\\left(\\lambda_1+k\\right)^2},\\frac{\\lambda_2}{\\left(\\lambda_2+k\\right)^2},\\cdots,\\frac{\\lambda_p}{\\left(\\lambda_p+k\\right)^2}\\right)\\hat\\alpha \\\\ \\\\ &=k^2\\sum_{i=1}^p\\frac{\\lambda_i\\hat\\alpha_i^2}{\\left(\\lambda_i+k\\right)^2} \\ . \\end{aligned} \\]

由于 \\(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p>0\\) ，所以对 \\(k>0\\) ，都有 \\(\\left(\\lambda_i+k\\right)^2>(\\lambda_i-k)^2\\) ，所以 \\(Q(k)<Q(-k)\\) 。这说明 \\(Q(k)\\) 的极小值不会在 \\((-\\infty,0)\\) 上取到，所以该极值问题取到极值点时，一定有 \\(k>0\\) 。

从几何上说，约束条件 \\(\\|b\\|^2=\\left\\|c\\hat\\beta\\right\\|^2\\xlongequal{def}h^2\\) 是一个中心在原点，半径为 \\(h\\) 的球面。对目标函数作椭球

\\[(b-\\hat\\beta)\'X\'X(b-\\hat\\beta)=\\delta^2 \\ , \\]

因为 \\(0<c<1\\) ，所以 \\(\\hat\\beta\\) 在 \\(\\|b\\|^2=\\left\\|c\\hat\\beta\\right\\|^2=h^2\\) 的球面之外。故总可以找到 \\(\\delta>0\\) 使得球 \\(\\|b\\|^2=h^2\\) 和目标函数的椭球相切于某点 \\(\\tilde\\beta\\) 。这个 \\(\\tilde\\beta\\) 就是上述极值问题的解，也就是岭估计 \\(\\hat\\beta(k)\\) 。

以二维向量为例，该极值问题如下图所示

容易看出，该极值问题与如下优化问题等价：

\\[\\min_\\beta\\ \\|Y-X\\beta\\|^2+\\lambda\\|\\beta\\|^2 \\ , \\]

其中 \\(\\lambda\\geq0\\) 称为调节参数。在这个优化问题中，\\(\\|Y-X\\beta\\|^2\\) 称为损失函数，\\(\\|\\beta\\|^2\\) 称为惩罚函数，参数 \\(\\lambda\\) 用于在损失和惩罚之间控制权衡。

3.9 主成分估计

3.9.1 主成分估计的过程

主成分估计的基本思想：

  1. 首先借助于正交变换将回归自变量变为对应的主成分，要求主成分的观测向量是正交的，且某些观测向量近似为 \\(0\\) 向量；
  2. 从所有的主成分中删去观测向量近似为 \\(0\\) 的主成分，起到消除多重共线性和降维的双重作用；
  3. 将保留下来的主成分作为新的回归自变量建立回归模型，用最小二乘估计模型中的回归系数，得到主成分回归方程。
  4. 基于得到的主成分回归方程，将它们转换为原始变量，即可得到原来的回归方程。

为了消除量纲的影响，假设自变量与因变量均已标准化。考虑回归模型

\\[Y=X\\beta+e \\ , \\quad {\\rm E}(e)=0 \\ , \\quad {\\rm Cov}(e)=\\sigma^2I_n \\ , \\]

记 \\(\\lambda_1\\geq\\lambda_2\\geq\\cdots\\geq\\lambda_p>0\\) 为 \\(X\'X\\) 的特征根，\\(\\phi_1,\\phi_2,\\cdots,\\phi_p\\) 为对应的标准正交化特征向量，则有

\\[\\Phi=(\\phi_1,\\phi_2,\\cdots,\\phi_p) \\]

为 \\(p\\times p\\) 正交矩阵，且有

\\[\\Phi\'X\'X\\Phi={\\rm diag}\\left(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p\\right)\\xlongequal{def}\\Lambda \\ . \\]

再记 \\(Z=X\\Phi,\\,\\alpha=\\Phi\'\\beta\\) ，则线性回归模型可以改写为

\\[Y=Z\\alpha+e \\ , \\quad {\\rm E}(e)=0 \\ , \\quad {\\rm Cov}(e)=\\sigma^2I_n \\ . \\]

主成分的性质：任意两个的主成分的观测向量都是互不相关的，且第 \\(j\\) 个主成分的偏差平方和

\\[\\sum_{i=1}^n\\left(z_{ij}-\\bar{z}_j\\right)^2=\\lambda_j \\ . \\]

因为 \\(Z\'Z=\\Phi\'X\'X\\Phi=\\Lambda={\\rm diag}\\left(\\lambda_1,\\lambda_2,\\cdots,\\lambda_p\\right)\\) ，所以

\\[z_j\'z_k=0 \\ , \\quad \\forall j\\neq k \\ , \\]

且 \\(z_j\'z_j=\\lambda_j,\\,j=1,2,\\cdots,p\\) 。又因为 \\(X\\) 是标准化设计矩阵，所以

\\[\\bar{z}_j=\\frac1n\\sum_{i=1}^nz_{ij}=\\frac1n\\sum_{i=1}^n\\sum_{k=1}^p\\phi_{kj}x_{ik}=\\frac{1}{n}\\sum_{k=1}^p\\phi_{kj}\\sum_{i=1}^nx_{ik}=0 \\ . \\]

因此有

\\[\\sum_{i=1}^n\\left(z_{ij}-\\bar{z}_j\\right)^2=\\sum_{i=1}^nz_{ij}^2=z_j\'z_j=\\lambda_j \\ , \\quad j=1,2,\\cdots,p \\ . \\]

由此可知 \\(\\lambda_j\\) 度量了第 \\(j\\) 个主成分 \\(z_j\\) 取值变动的大小。因为 \\(\\lambda_1\\geq\\lambda_2\\geq\\cdots\\geq\\lambda_p>0\\) ，所以我们称 \\(z_1\\) 为第一主成分，称 \\(z_2\\) 为第二主成分，以此类推。这 \\(p\\) 个主成分的观测向量是正交的。

由主成分的性质可知，\\(z_1\\) 对因变量的解释能力最强，\\(z_2\\) 次之，\\(z_p\\) 最弱。若设计矩阵 \\(X\\) 是病态矩阵，则存在一些 \\(X\'X\\) 的特征根接近于 \\(0\\) 。

不妨设 \\(\\lambda_{r+1},\\lambda_{r+2},\\cdots,\\lambda_p\\approx0\\) ，此时后面的 \\(p-r\\) 个主成分的取值变动就很小，且均在 \\(0\\) 附近取值，所以这 \\(p-r\\) 个主成分对因变量的影响可以忽略，可将它们从回归模型中剔除。

剩下的主成分 \\(z_1,z_2,\\cdots,z_r\\) 就不存在多重共线性问题，用最小二乘法对剩下的 \\(r\\) 个主成分做回归即可。

我们用分块的方式建立回归方程：对 \\(\\Lambda,\\alpha,Z,\\Phi\\) 进行分块

\\[\\Lambda=\\begin{pmatrix} \\Lambda_1 & 0 \\\\ 0 & \\Lambda_2 \\end{pmatrix} \\ , \\quad \\alpha=\\begin{pmatrix} \\alpha_1 \\\\ \\alpha_2 \\end{pmatrix} \\ , \\quad Z=\\left(\\begin{array}{c:c}Z_1&Z_2\\end{array}\\right) \\ , \\quad \\Phi=\\left(\\begin{array}{c:c}\\Phi_1& \\Phi_2\\end{array}\\right) \\ , \\]

其中 \\(\\Lambda_1\\) 为 \\(r\\times r\\) 矩阵，\\(\\alpha_1\\) 为 \\(r\\times1\\) 向量，\\(Z_1\\) 为 \\(n\\times r\\) 矩阵，\\(\\Phi_1\\) 为 \\(p\\times r\\) 矩阵。因为 \\(Z_2\\) 近似是 \\(0\\) 矩阵，所以剔除 \\(Z_2\\alpha_2\\) ，模型可以写为

\\[Y\\approx Z_1\\alpha_1+e \\ , \\quad {\\rm E}(e)=0 \\ , \\quad {\\rm Cov}(e)=\\sigma^2I_n \\ . \\]

这里 \\(Z_1\\) 不是病态矩阵，因为 \\(Z_1\'Z_1\\) 的特征根为 \\(\\lambda_1,\\lambda_2,\\cdots,\\lambda_r\\) 均远离 \\(0\\) ，所以可直接利用最小二乘法求得 \\(\\alpha_1\\) 的最小二乘估计

\\[\\hat\\alpha_1=\\left(Z_1\'Z_1\\right)^{-1}Z_1\'Y=\\Lambda_1^{-1}Z_1\'Y \\ . \\]

前面我们从模型中剔除了后面的 \\(p-r\\) 个主成分，这相当于我们用 \\(\\hat\\alpha_2=0\\) 来估计 \\(\\alpha_2\\) ，利用 \\(\\beta=\\Phi\\alpha\\) 可以得到 \\(\\beta\\) 的主成分估计为

\\[\\begin{aligned} \\tilde\\beta=\\Phi\\hat\\alpha=\\left(\\begin{array}{c:c}\\Phi_1& \\Phi_2\\end{array}\\right)\\begin{pmatrix} \\hat\\alpha_1 \\\\ \\hat\\alpha_2 \\end{pmatrix}=\\Phi_1\\hat\\alpha_1=\\Phi_1\\Lambda_1^{-1}Z_1\'Y=\\Phi_1\\Lambda_1^{-1}\\Phi_1\'X\'Y \\ . \\end{aligned} \\]

相应的主成分回归方程为 \\(\\hat{Y}=X\\tilde{\\beta}\\) 。

3.9.2 主成分估计的性质

上述主成分估计的过程可以概括为以下三步：

  1. 做正交变换 \\(Z=X\\Phi\\) ，获得新的自变量，称为主成分；
  2. 做回归自变量选择，提出特征根比较小的主成分；
  3. 用标准化后的 \\(y\\) 对剩余的主成分做回归，得到最小二乘估计和主成分回归方程，再将这个回归方程转换为关于原始变量的回归方程。

经过上述过程得到的估计量 \\(\\tilde\\beta\\) 称为 \\(\\beta\\) 的主成分估计，下面我们来研究 \\(\\tilde\\beta\\) 的统计性质。

性质 1：主成分估计 \\(\\tilde\\beta=\\Phi_1\\Phi_1\'\\hat\\beta\\) 是最小二乘估计的一个线性变换。

根据下列关系

\\[\\Phi_1\'\\Phi_1=I_r \\ , \\quad \\Phi_1\'\\Phi_2=0 \\ , \\quad X\'X=\\Phi\\Lambda\\Phi\'=\\Phi_1\\Lambda_1\\Phi_1\'+\\Phi_2\\Lambda_2\\Phi_2\' \\ , \\]

可以得到

\\[\\begin{aligned} \\tilde\\beta&=\\Phi_1\\Lambda_1^{-1}\\Phi_1\'X\'Y \\\\ \\\\ &=\\Phi_1\\Lambda_1^{-1}\\Phi_1\'X\'X\\hat\\beta \\\\ \\\\ &=\\Phi_1\\Lambda_1^{-1}\\Phi_1\'\\Phi_1\\Lambda_1\\Phi_1\'\\hat\\beta+\\Phi_1\\Lambda_1^{-1}\\Phi_1\'\\Phi_2\\Lambda_2\\Phi_2\'\\hat\\beta \\\\ \\\\ &=\\Phi_1\\Phi_1\'\\hat\\beta \\ . \\end{aligned} \\]

性质 2：\\({\\rm E}(\\tilde{\\beta})=\\Phi_1\\Phi_1\'\\beta\\) ，即只要 \\(r<p\\) ，主成分估计就是有偏估计。

性质 3：\\(\\left\\|\\tilde\\beta\\right\\|\\leq\\left\\|\\hat\\beta\\right\\|\\) ，即主成分估计是压缩估计。

构造 \\(p\\times p\\) 矩阵 \\(\\tilde I={\\rm diag}\\left(I_r,0\\right)\\) ，则由 \\(\\Phi\\) 的定义可知

\\[\\Phi_1\\Phi_1\'=\\Phi\\tilde{I}\\Phi\' \\ . \\]

从而有

\\[\\left\\|\\tilde\\beta\\right\\|=\\left\\|\\Phi\\tilde{I}\\Phi\'\\hat\\beta\\right\\|=\\left\\|\\Phi\\right\\|\\left\\|\\tilde{I}\\Phi\'\\hat\\beta\\right\\|=1\\times \\left\\|\\tilde{I}\\Phi\'\\hat\\beta\\right\\|\\leq \\left\\|\\Phi\'\\hat\\beta\\right\\|=\\left\\|\\hat\\beta\\right\\| \\ . \\]

定理 3.9.1：当原始自变量存在足够严重的多重共线性时，适当选择保留的主成分个数可使主成分估计比最小二乘估计拥有较小的均方误差，即

\\[{\\rm MSE}(\\tilde\\beta)<{\\rm MSE}(\\hat\\beta) \\ . \\]

假设 \\(X\'X\\) 的后 \\(p-r\\) 个特征根 \\(\\lambda_{r+1},\\cdots,\\lambda_p\\) 很接近于 \\(0\\) ，不难看出

\\[\\begin{aligned} {\\rm MSE}(\\tilde\\beta)&={\\rm MSE}\\begin{pmatrix} \\hat\\alpha_1 \\\\ 0 \\end{pmatrix} \\\\ \\\\ &={\\rm tr}\\left[{\\rm Cov}\\begin{pmatrix} \\hat\\alpha_1 \\\\ 0 \\end{pmatrix} \\right]+\\left\\|{\\rm E}\\begin{pmatrix} \\hat\\alpha_1 \\\\ 0 \\end{pmatrix} -\\alpha\\right\\|^2 \\\\ \\\\ &=\\sigma^2{\\rm tr}\\left(\\Lambda_1^{-1}\\right)+\\left\\|\\alpha_2\\right\\|^2 \\ . \\end{aligned} \\]

因为

\\[{\\rm MSE}(\\hat\\beta)=\\sigma^2{\\rm tr}\\left(\\Lambda^{-1}\\right)=\\sigma^2{\\rm tr}\\left(\\Lambda_1^{-1}\\right)+\\sigma^2{\\rm tr}\\left(\\Lambda_2^{-1}\\right) \\ , \\]

所以

\\[{\\rm MSE}(\\tilde\\beta)={\\rm MSE}(\\hat\\beta)+\\left(\\left\\|\\alpha_2\\right\\|^2-\\sigma^2{\\rm tr}\\left(\\Lambda_2^{-1}\\right)\\right) \\ . \\]

于是

\\[{\\rm MSE}(\\tilde\\beta)<{\\rm MSE}(\\hat\\beta) \\quad \\iff \\quad \\left\\|\\alpha_2\\right\\|^2<\\sigma^2{\\rm tr}\\left(\\Lambda_2^{-1}\\right)=\\sigma^2\\sum_{i=r+1}^p\\frac{1}{\\lambda_i} \\ . \\]

当多重共线性足够严重的时候，\\(\\lambda_{r+1},\\cdots,\\lambda_{p}\\) 中的某一个可以充分接近于 \\(0\\) ，因此上式右端可以足够大使得该不等式成立。

因为 \\(\\alpha_2=\\Phi_2\'\\beta\\) ，我们可以将上述不等式写为

\\[\\left(\\frac{\\beta}{\\sigma}\\right)\'\\Phi_2\\Phi_2\'\\left(\\frac{\\beta}{\\sigma}\\right)<{\\rm tr}\\left(\\Lambda_2^{-1}\\right) \\ , \\]

这就是说，当 \\(\\beta\\) 和 \\(\\sigma\\) 满足该不等式时，主成分估计才比最小二乘估计拥有较小的均方误差。如果将 \\(\\beta/\\sigma\\) 视为参数空间中的一个参数，则上述不等式表示一个中心在原点的椭球，有以下两个结论：

  • 对给定的参数 \\(\\beta\\) 和 \\(\\sigma^2\\) ，当 \\(X\'X\\) 的后 \\(p-r\\) 个特征根比较小时，主成分估计比最小二乘估计拥有较小的均方误差；
  • 对给定的\\(X\'X\\) ，即固定的 \\(\\Lambda_2\\) ，对绝对值相对较小的 \\(\\beta/\\sigma\\) ，主成分估计比最小二乘估计拥有较小的均方误差。

最后我们给出两条主成分个数 \\(r\\) 的选取准则：

  1. 略去特征根接近于 \\(0\\) 的主成分；
  2. 选择 \\(r\\) 使得累计贡献率（前 \\(r\\) 个特征根之和在 \\(p\\) 个特征根之和中所占的比例）达到预想给定的值。

根据经验，在实际操作中，我们一般选择最小的 \\(r\\) 使得

\\[\\frac{\\sum_{i=1}^r\\lambda_i}{\\sum_{i=1}^p\\lambda_i}>0.85 \\ . \\]

以上是关于高维数据惩罚回归方法：主成分回归PCR岭回归lasso弹性网络elastic net分析基因数据|附代码数据的主要内容，如果未能解决你的问题，请参考以下文章

主成分回归的一般步骤是怎样的

回归分析06：回归参数的估计

Matlab主成分回归（pcr）分析中的常数项

sklearn help之岭回归 ridge regression

从最小二乘到岭回归(Ridge Regression)的深刻理解

岭回归和Lasso回归有啥区别？