关于数据挖掘:高维数据惩罚回归方法主成分回归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

向下滑动查看后果▼

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 值的减少而变动。

点击题目查阅往期内容

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 值的序列

向下滑动查看后果▼

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

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

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

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

向下滑动查看后果▼

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

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

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

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

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

ridge_preds <- predict
## 计算 MSE

向下滑动查看后果▼

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

向下滑动查看后果▼

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