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

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 200qr(X)$rank#> \[1\] 119XtX <- crossprod(X) # 更无效地计算t(X) %*% Xqr(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=4Vk <- pca$rotation\[, 1:k\] # 载荷矩阵Zk <- pca$x\[, 1:k\] # 分数矩阵# 在经典的线性回归中应用这些分数

 因为X和Y是中心化的，截距近似为0。

输入结果显示，PC1和PC4的估计值与0相差很大（在p<0.05），然而后果不能轻易解释，因为咱们没有对PC的间接解释。

2.2 应用软件包

PCR也能够间接在数据上进行（所以不用先手动进行PCA）。在应用这个函数时，你必须牢记几件事。

  1. 要应用的成分（PC）的数量是通过参数ncomp来确定
  2. 该函数容许你首先对预测因子进行标准化（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)再次成为满秩，并且是可逆的。能够应用两种不同的惩办项或正则化办法。

  1. L1正则化：这种正则化在预计方程中退出一个1‖‖1。该项将减少一个基于系数大小绝对值的惩办。这被Lasso回归所应用。

  1. 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 200qr(XtX_gammaI)$rank == 200 # #> \[1\] TRUE

2. 查看的逆值是否能够计算出来。

# 是的，能够被计算。XtX\_gammaI\_inv <- solve(XtX_gammaI)

3. 最初，计算。

## 计算岭估计值## 应用\`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值的减少而变动。

请留神，对于更高的值，系数估计值变得更靠近于0，显示了岭惩办的膨胀效应。

与PC回归的例子相似，咱们相当随便地抉择了=2和网格。咱们随后会看到，如何抉择，使预测误差最小。

6 练习: Lasso 回归

Lasso 回归也是惩罚性回归的一种模式，但咱们没有像最小二乘法和岭回归那样的^的剖析解。为了拟合一个Lasso 模型，咱们再次应用glmnet()函数。然而，这一次咱们应用的参数是=1

工作

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

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

你不用在这里提供一个自定义的（lambda）值序列，而是能够依附glmnet的默认行为，即依据数据抉择值的网格。

# 请留神，glmnet()函数能够主动提供伽马值# 默认状况下，它应用100个lambda值的序列

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

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回归对象，并像以前一样制作系数曲线图。 

咱们能够寻找能产生最佳成果的伽玛值。这里有两种可能性。

  1. lambda.min: 给出穿插验证最佳后果的值。
  2. lambda.1se:的最大值，使MSE在穿插验证的最佳后果的1个标准误差之内。

咱们将在这里应用lambda.min来拟合最终模型，并在测试数据上生成预测。请留神，咱们实际上不须要从新进行拟合，咱们只须要应用咱们现有的lasso_cv对象，它曾经蕴含了lambda值范畴的拟合模型。咱们能够应用predict函数并指定s参数（在这种状况下设置lambda）来对测试数据进行预测。 

2. 对岭回归做同样的解决。

请留神，咱们能够从CV后果中提取拟合的岭回归对象，并制作系数曲线图。

咱们能够寻找能产生最佳成果的伽玛值。这里有两种可能性。

  1. lambda.min: 给出穿插验证最佳后果的值。
  2. lambda.1se: 的最大值，使MSE在穿插验证的最佳后果的1个标准误差之内。

咱们在这里应用lambda.min来拟合最终的模型并在测试数据上生成预测。请留神，咱们实际上不须要从新进行拟合，咱们只须要应用咱们现有的ridge_cv对象，它曾经蕴含了lambda值范畴的拟合模型。咱们能够应用predict函数并指定s参数（在这种状况下凌乱地设置lambda）来对测试数据进行预测。

ridge_preds <- predict##计算MSE

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

模型

MSE

PCR

0.3655052

Lasso

0.3754368

Ridge

0.3066121

  1. 留神：R中的log()默认是自然对数（以e为底），咱们也会在文本中应用这个符号（比方下面图中的x轴题目）。这可能与你所习惯的符号（ln()）不同。要在R中取不同基数的对数，你能够指定log的基数=参数，或者应用函数log10(x)和log2(x)别离代表基数10和2︎

最受欢迎的见解 

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指标