跳至主要内容

GFLASSO：R 中的图引导融合 LASSO

通过本教程，使用 GFLASSO R 包探索图结构的多任务回归！
· 13分钟 阅读

使用 AI 探索

尽管机器学习领域伴随着越来越复杂的方法迅速发展，但对需要同时预测多个响应的高维多任务问题关注较少。本教程将展示图引导融合 LASSO（GFLASSO）在单一正则化的线性回归框架下预测多个响应的威力。

简介

在监督学习中，我们通常希望根据一组解释变量或预测变量，对一组样本或观测的因变量或响应变量进行预测。正则化方法通过引入惩罚项来防止高维数据的过拟合，尤其是在预测变量个数超过观测数的情况下。这些惩罚项被添加到目标函数中，使得对误差最小化贡献很小的非信息性预测变量的系数估计也被最小化。最小绝对收缩与选择算子（LASSO）[1] 就是一种这样的方法。

什么是 LASSO？

与普通最小二乘法（OLS）相比，LASSO 能够将系数估计（β）精确收缩到零，从而排除非信息性预测变量并执行特征选择，其形式为

$$argmin_\beta \sum_n(y_n-\hat{y_n})^2+\lambda\sum_{j}|\beta_{j}|$$

其中，n 和 j 分别表示任意给定的观测和预测变量。残差平方和（RSS），即 OLS 中唯一使用的项，也可等价写为代数形式 $RSS = \sum_n(y_n-\hat{y_n})^2 = (y-X\beta)^T.(y-X\beta)$。LASSO 的惩罚项为 λ∑j|βj|，即由 λ 加权的系数估计的 L1 范数。

为何使用图引导融合 LASSO（GFLASSO）？

如果您希望从同一组预测变量同时预测多个相关的响应，该怎么办？虽然为每个响应分别拟合多个独立的 LASSO 模型也能取得不错效果，但如果能按照响应之间的关联强度来协调这些预测，效果会更好。因为这种协调可以抵消包含噪声的响应特异性变化——这正是 GFLASSO 的关键优势。

原始论文给出了一个很好的示例，作者在其中解析了 543 名患者中 34 个遗传标记与 53 项哮喘性状之间的关联 [2]。

什么是 GFLASSO？

设X 为大小为 n × p 的矩阵，包含 n 个观测和 p 个预测变量；Y 为大小为 n × k 的矩阵，包含相同的 n 个观测和 k 个响应。比如，73 个国家中 1390 条不同电子产品的购买记录，用于预测所有 73 个国家对 50 部 Netflix 制作的评分。

适合建模两组高维数据的模型包括正交双向偏最小二乘（O2PLS）、典型相关分析（CCA）和共同惯量分析（CIA），这些方法都涉及矩阵分解 [3]。另外，由于这些模型基于潜变量（即基于原始预测变量的投影），计算效率的提升是以可解释性为代价的。

然而，这种权衡并非总是划算的。通过在统一的回归框架中，直接用 X 中选定的特征去预测 k 个单独的响应，并考虑响应之间的关系，可以扭转这种局面。

从数学上讲，GFLASSO 在上述 LASSO [1] 正则化的基础上，进一步利用由 k × k 相关矩阵量化的 Y 的图依赖结构（也就是您之前读到的“关联强度”）来构建模型。结果是，相似（或不相似）的响应会由相似（或不相似）的被选预测变量子集来解释。

更正式地，沿用原始论文 [2] 中的记号，GFLASSO 的目标函数为

function

其中，跨越所有 k 个响应，∑k(yk − X**βk)T.(yk − X**βk) 给出 RSS，而 λ∑k∑j|βj**k| 为借自 LASSO 的正则化惩罚项，由参数 λ 加权，并作用于每个预测变量 j 的系数 β。GFLASSO 的新颖之处在于

function

由 γ 加权的融合惩罚，它保证任一预测变量 j 在任意一对响应 m 与 l 上的系数 βj**m 与 βj**l 的绝对差，会随着它们的成对相关（可变换或不变换的 f(rm**l））越正相关而越小（反之越负相关而越大）。这种融合惩罚倾向于保留响应中的整体有意义变化，而非各自的噪声。当成对相关接近于零时，它不起作用，此时模型退化为纯 LASSO。所有 k 个响应的底层相关结构可表示为加权网络结构，默认采用绝对相关，f(rm**l)=|rm**l|，也可通过任意用户自定义函数变换，形成 GFLASSO 的变体，例如：

  1. 平方相关，f(rm**l)=rm**l2（加权）
  2. 阈值相关，$f(r_{ml}) = \begin{cases} 1, & \mbox{if } r_{ml} > \tau \\ 0, & \mbox{otherwise} \end{cases}$（不加权）

还有很大的创新空间。尽管与 1. 和默认的绝对相关 [2] 相比，方案 2. 的计算量小得多，但它需要预先设定一个阈值，例如 τ = 0.8。

总之，要拟合一个 GFLASSO 模型，您需要预测变量矩阵 X、响应矩阵 Y，以及描绘 Y 中所有响应两两间关联强度的相关矩阵。请注意，GFLASSO 得到的是一个 p × k 的 β 矩阵，不同于 LASSO 的（p × 1），该系数矩阵承载任一给定响应 k 与预测变量 j 的关联。

开始使用

Kris Sankaran 和我开发了一个实验性的R 包，实现了 GFLASSO，并提供交叉验证与绘图方法。我们最近使用 doParallel 包实现了多线程，大幅加快了交叉验证（CV）流程。

要在 R 中运行 GFLASSO，您需要安装 devtools，加载它，并从我的 GitHub 仓库安装 gflasso 包。演示将基于 bgsmtr 包中包含的数据集进行。我也建议安装 corrplot 和 pheatmap 来可视化结果。

# Install the packages if necessary
#install.packages("devtools")
#install.packages("bgsmtr")
#install.packages("corrplot")
#install.packages("pheatmap")
library(devtools)
library(bgsmtr)
library(corrplot)
library(pheatmap)
#install_github("monogenea/gflasso")
library(gflasso)

模拟

您可以轻松运行 cv_gflasso() 的帮助页面中概述的模拟。默认情况下，CV 在 5 折交叉验证的一次重复中，跨越 λ ∈ {0, 0.1, 0.2, ..., 0.9, 1} 与 γ ∈ {0, 0.1, 0.2, ..., 0.9, 1} 的所有可能组合（调参网格），计算均方根误差（RMSE）。

注意：用户自定义的误差函数同样适用！

除固有的统计假设与速度表现外，调参网格范围的选择在很大程度上取决于对 X 与 Y 所有列进行均值中心化与单位方差缩放，请务必事先完成此步骤。

在下面的示例中您不需要这样做，因为您将从标准正态分布中随机抽样。您可以尝试从一个不加权的相关网络推导融合惩罚，阈值为 r > 0.8：

?cv_gflasso
set.seed(100)
X <- matrix(rnorm(100 * 10), 100, 10)
u <- matrix(rnorm(10), 10, 1)
B <- u %*% t(u) + matrix(rnorm(10 * 10, 0, 0.1), 10, 10)
Y <- X %*% B + matrix(rnorm(100 * 10), 100, 10)
R <- ifelse(cor(Y) > .8, 1, 0)
system.time(testCV <- cv_gflasso(X, Y, R, nCores = 1))

## [1] 1.826146 1.819430 1.384595 1.420058 1.408619

##    user  system elapsed 
##  45.808   3.070  55.271

system.time(testCV <- cv_gflasso(X, Y, R, nCores = 2))

## [1] 1.595413 1.492953 1.469917 1.366832 1.441642

##    user  system elapsed 
##  22.231   1.698  26.441

cv_plot_gflasso(testCV)

cv mean rmse

在本次模拟中，使 RMSE 最小的最佳 λ（行）与 γ（列）分别为 0.3 与 0.8，确实捕捉到了所设定的关系。

提示：尝试用另一种度量重新运行该示例——决定系数（R2）。使用 R2 的一个关键优势是其取值范围为 0 到 1。

请记住，当您提供自定义的拟合优度函数 err_fun() 时，需要通过参数 err_opt 指定是最大化还是最小化所得到的度量。

下面的示例旨在最大化 R2，使用带平方相关系数的加权关联网络（即 f(rm**l)=rm**l2）。如果您有超过 2 个内核，也可以修改 nCores 参数加速！

# Write R2 function
R2 <- function(pred, y){
      cor(as.vector(pred), as.vector(y))**2
}

# X, u, B and Y are still in memory
R <- cor(Y)**2

# Change nCores, if you have more than 2, re-run CV
testCV <- cv_gflasso(X, Y, R, nCores = 2, err_fun = R2, err_opt = "max")

## [1] 0.6209191 0.7207394 0.7262781 0.7193907 0.6303187

cv_plot_gflasso(testCV)

cv mean r2

此时，最佳参数 λ 与 γ 分别为 0.6 和 0.3。

此外，请注意，cv_gflasso 对象是包含四个元素的单一列表：网格中所有单元的度量均值（$mean）与标准误（$SE），最佳的 λ 与 γ 参数（$optimal），以及拟合优度函数的名称（$err_fun）。当前示例的交叉验证模型明显同时倾向于稀疏性（λ）与融合（γ）。

最后，请记住，您可以微调其他参数，例如 Nesterov 梯度的收敛阈值 δ 和最大迭代次数，分别通过 delta_conv 与 iter_max 传递给 additionalOpts。这些将在下一个示例中使用。

使用 GFLASSO 确定 SNP—神经影像学的关联

为展示 GFLASSO 在相对高维问题中的简洁性与稳健性，接下来您将使用阿尔茨海默病神经影像学倡议（ADNI-1）数据库中的 bgsmtr_example_data 数据集进行建模。

这是一个包含 3 个元素的列表对象，属于 bgsmtr 包，其中包括 15 项结构性神经影像测度和 486 个单核苷酸多态性（SNP，遗传标记），样本量为 632 名受试者。重要的是，这 486 个 SNP 涵盖了被认为与阿尔茨海默病相关的 33 个基因。

您的任务是从 SNP 数据预测形态学神经影像学测度，并利用前者的相关结构。

让我们先整理数据并探索所有神经影像特征之间的相互依赖关系：

data(bgsmtr_example_data)
str(bgsmtr_example_data)

## List of 3
##  $ SNP_data     : int [1:486, 1:632] 2 0 2 0 0 0 0 1 0 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:486] "rs4305" "rs4309" "rs4311" "rs4329" ...
##   .. ..$ : chr [1:632] "V1" "V2" "V3" "V4" ...
##  $ SNP_groups   : chr [1:486] "ACE" "ACE" "ACE" "ACE" ...
##  $ BrainMeasures: num [1:15, 1:632] 116.5 4477.9 28631.9 34.1 -473.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:15] "Left_AmygVol.adj" "Left_CerebCtx.adj" "Left_CerebWM.adj" "Left_HippVol.adj" ...
##   .. ..$ : chr [1:632] "V1" "V2" "V3" "V4" ...

# Transpose, so that samples are distributed as rows, predictors / responses as columns
SNP <- t(bgsmtr_example_data$SNP_data)
BM <- t(bgsmtr_example_data$BrainMeasures)

# Define dependency structure
DS <- cor(BM)

# Plot correlation matrix of the 15 neuroimaging measures
corrplot(DS)

example

上图显示了神经影像特征之间的相互依赖关系。现在尝试对 GFLASSO 进行交叉验证（在笔记本电脑上可能需要几个小时！），并确定 SNP—神经影像学的关联。

注意：在下面的示例中，指定了收敛容差和最大迭代次数。欢迎您自行尝试不同的取值！

system.time(CV <- cv_gflasso(X = scale(SNP), Y = scale(BM), R = DS, nCores = 2,
                 additionalOpts = list(delta_conv = 1e-5, iter_max = 1e5)))

## [1] 1.550294 1.471637 1.470133 1.514425 1.504215

##      user    system   elapsed 
##  2611.492   323.441 51129.936

cv_plot_gflasso(CV)

example

通过将 GFLASSO 与纯 LASSO（γ = 0，第一列）、纯融合最小二乘（λ = 0，第一行）和 OLS（γ = 0 且 λ = 0，左上角单元）进行比较，我们可以得出结论：当前示例最好使用非零惩罚权重来建模，也即采用完整的 GFLASSO。取交叉验证的最优参数（λ = 0.7，γ = 1）来构建 GFLASSO 模型，并解读得到的系数矩阵：

gfMod <- gflasso(X = scale(SNP), Y = scale(BM), R = DS, opts = list(lambda = CV$optimal$lambda,
                                                                    gamma = CV$optimal$gamma, 
                                                                    delta_conv = 1e-5,
                                                                    iter_max = 1e5))
colnames(gfMod$B) <- colnames(BM)
pheatmap(gfMod$B, annotation_row = data.frame("Gene" = bgsmtr_example_data$SNP_groups,
                                              row.names = rownames(gfMod$B)),
         show_rownames = F)

example

上图显示了非常大比例的系数为零或近零。尽管在基因层面上（行注释，注意图例不完整）看不到明显的 SNP 聚类，但某些 SNP 与性状之间存在清晰的关联。

为确认是否存在非随机的预测机制，您可以在置换 X 或 Y 的取值后重复该过程。实验研究有助于从选定的 SNP 中阐明因果关系与作用机制。例如，影响蛋白序列与结构、进而干扰阿尔茨海默病中 β-淀粉样斑块清除的 SNP。

总结

GFLASSO 在建模多个响应时同时采用正则化与融合，从而有助于识别预测变量（X）与响应（Y）之间的关联。它最适用于观测很少而数据维度很高的情形，因为其速度比竞品方法更慢。例如，稀疏条件高斯图模型 [4] 和贝叶斯组稀疏多任务回归模型 [5] 可能主要因性能优势而更受青睐。尽管如此，GFLASSO 的可解释性很强。我最近在一个组学整合的研究中使用 GFLASSO 揭示了玉米中的新脂质基因 [6]。

参阅 DataCamp 的正则化教程：Ridge、Lasso 与 Elastic Net。

Kris 和我非常期待您的来信。该项目目前由 Kris 在仓库 krisrs1128/gflasso 维护，也可通过我的仓库获取，尽管变更较频繁：monogenea/gflasso。欢迎随时写信给我（[email protected]），非常感谢您的任何反馈。

祝您编码愉快！

参考文献

  1. Robert Tibshirani (1994). Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society, 58, 267-288.
  2. Seyoung Kim, Kyung-Ah Sohn, Eric P. Xing (2009). A multivariate regression approach to association analysis of a quantitative trait network. Bioinformatics, 25, 12:i204–i212.
  3. Chen Meng, Oana A. Zeleznik, Gerhard G. Thallinger, Bernhard Kuster, Amin M. Gholami, Aedín C. Culhane (2016). Dimension reduction techniques for the integrative analysis of multi-omics data. Briefings in Bioinformatics, 17, 4:628–641.
  4. Lingxue Zhang, Seyoung Kim (2014). Learning Gene Networks under SNP Perturbations Using eQTL Datasets. PLoS Comput Biol, 10, 2:e1003420.
  5. Keelin Greenlaw, Elena Szefer, Jinko Graham, Mary Lesperance, Farouk S. Nathoo (2017). A Bayesian group sparse multi-task regression model for imaging genetics. Bioinformatics, 33, 16:2513–2522.
  6. Francisco de Abreu e Lima, Kun Li, Weiwei Wen, Jianbing Yan, Zoran Nikoloski, Lothar Willmitzer, Yariv Brotman (2018). Unraveling the lipid metabolism in maize with time-resolved multi-omics data. The Plant Journal, 93, 6:1102-1115.

主题

R数据科学