GFLASSO: LASSO Hợp nhất có hướng dẫn bởi đồ thị trong R

Mặc dù lĩnh vực học máy đang tiến rất nhanh với sự ra đời của những kỹ thuật ngày càng tinh vi, các bài toán đa nhiệm chiều cao đòi hỏi dự đoán đồng thời nhiều biến phản hồi lại ít được chú ý. Hướng dẫn này sẽ cho bạn thấy sức mạnh của Graph-Guided Fused LASSO (GFLASSO) trong việc dự đoán nhiều phản hồi trong cùng một khung hồi quy tuyến tính được chuẩn hóa.

Giới thiệu

Trong học có giám sát, mục tiêu thường là dự đoán một biến phụ thuộc hay biến phản hồi từ một tập biến giải thích hay biến dự đoán trên một tập mẫu hay quan sát. Các phương pháp chuẩn hóa đưa vào các hình phạt nhằm ngăn quá khớp dữ liệu có chiều cao, đặc biệt khi số biến dự đoán vượt quá số quan sát. Những hình phạt này được cộng vào hàm mục tiêu để các ước lượng hệ số của những biến dự đoán kém thông tin (đóng góp ít vào việc giảm thiểu sai số) cũng bị giảm thiểu. Toán tử co rút và chọn lựa tối thiểu (LASSO) [1] là một phương pháp như vậy.

LASSO là gì?

So với bình phương tối thiểu thông thường (OLS), LASSO có khả năng co rút các ước lượng hệ số (β) về đúng 0, qua đó loại bỏ các biến dự đoán không hữu ích và thực hiện chọn đặc trưng, thông qua

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

trong đó n và j lần lượt ký hiệu một quan sát và một biến dự đoán bất kỳ. Tổng bình phương phần dư (RSS), hạng tử duy nhất dùng trong OLS, có thể được viết tương đương ở dạng đại số $RSS = \sum_n(y_n-\hat{y_n})^2 = (y-X\beta)^T.(y-X\beta)$. Hình phạt LASSO là λ∑_j|β_j|, là chuẩn L1 của các ước lượng hệ số, được trọng số bởi λ.

Tại sao là Graph-Guided Fused LASSO (GFLASSO)?

Nếu bạn cần dự đoán đồng thời nhiều phản hồi có liên quan, từ cùng một tập biến dự đoán thì sao? Dù bạn có thể đạt kết quả ổn với nhiều mô hình LASSO độc lập, mỗi mô hình cho một phản hồi, nhưng bạn sẽ làm tốt hơn nếu phối hợp các dự đoán theo sức mạnh liên hệ giữa các phản hồi. Bởi vì sự phối hợp đó khử đi biến thiên đặc thù theo từng phản hồi có chứa nhiễu - chính là điểm mạnh của GFLASSO.

Một ví dụ hay được nêu trong bài báo gốc, nơi các tác giả phân giải liên hệ giữa 34 chỉ thị di truyền và 53 đặc tính hen suyễn ở 543 bệnh nhân [2].

GFLASSO là gì?

Gọi X là ma trận kích thước n × p, với n quan sát và p biến dự đoán, và Y là ma trận kích thước n × k, với cùng n quan sát và k phản hồi, ví dụ, 1390 bản ghi mua sắm đồ điện tử khác nhau ở 73 quốc gia, nhằm dự đoán xếp hạng của 50 tác phẩm Netflix trên cả 73 quốc gia.

Các mô hình phù hợp để mô hình hóa cặp dữ liệu chiều cao bao gồm Partial Least Squares hai chiều trực giao (O2PLS), Phân tích Tương quan Canonical (CCA) và Co-Inertia Analysis (CIA), tất cả đều liên quan đến phân rã ma trận [3]. Thêm nữa, vì các mô hình này dựa trên biến ẩn (tức là các phép chiếu dựa trên biến dự đoán gốc), hiệu quả tính toán phải trả giá bằng tính dễ diễn giải.

Tuy nhiên, sự đánh đổi này không phải lúc nào cũng xứng đáng, và có thể đảo ngược bằng cách dự đoán trực tiếp k phản hồi riêng lẻ từ các đặc trưng được chọn trong X, trong một khung hồi quy hợp nhất có xét đến quan hệ giữa các phản hồi.

Về mặt toán học, GFLASSO kế thừa chuẩn hóa của LASSO [1] như đã bàn ở trên và xây dựng mô hình dựa trên cấu trúc phụ thuộc dạng đồ thị nằm dưới Y, được định lượng bởi ma trận tương quan k × k (chính là “sức mạnh liên hệ” mà bạn đã đọc). Hệ quả là các phản hồi giống nhau (hoặc khác nhau) sẽ được giải thích bởi một tập con biến dự đoán được chọn giống nhau (hoặc khác nhau).

Chính xác hơn, theo ký hiệu dùng trong bản thảo gốc [2], hàm mục tiêu của GFLASSO là

trong đó, trên toàn bộ k phản hồi, ∑_k(y_k − X**β_k)^T.(y_k − X**β_k) cho RSS và λ∑_k∑_j|β_j**k| là hình phạt chuẩn hóa vay mượn từ LASSO, được trọng số bởi tham số λ và tác động lên các hệ số β của từng biến dự đoán j. Điểm mới của GFLASSO nằm ở

hình phạt hợp nhất (fusion) được trọng số bởi γ, đảm bảo độ chênh tuyệt đối giữa các hệ số β_j**m và β_j**l, từ bất kỳ biến dự đoán j và cặp phản hồi m, l, sẽ nhỏ hơn (tương ứng lớn hơn) khi tương quan cặp của chúng dương hơn (tương ứng âm hơn), đã biến đổi hay chưa, f(r_m**l). Hình phạt hợp nhất này ưu tiên biến thiên có ý nghĩa toàn cục trong các phản hồi hơn là nhiễu riêng lẻ của từng phản hồi. Khi tương quan cặp xấp xỉ 0, nó không tác động, khi đó bạn còn lại LASSO thuần túy. Cấu trúc tương quan nền của toàn bộ k phản hồi, có thể biểu diễn như một mạng có trọng số, mặc định là tương quan tuyệt đối, f(r_m**l)=|r_m**l| nhưng có thể biến đổi để tạo biến thể GFLASSO theo bất kỳ hàm do người dùng chỉ định, chẳng hạn

1. Tương quan bình phương, f(r_m**l)=r_m**l² (có trọng số)
2. Tương quan ngưỡng, $f(r_{ml}) = \begin{cases} 1, & \mbox{if } r_{ml} > \tau \\ 0, & \mbox{otherwise} \end{cases}$ (không trọng số)

và còn nhiều dư địa để sáng tạo. Mặc dù cách 2. ít tốn kém tính toán hơn nhiều so với 1. và so với mặc định tương quan tuyệt đối [2], nó đòi hỏi một ngưỡng định sẵn, ví dụ, τ = 0.8.

Tóm lại, để khớp một mô hình GFLASSO, bạn sẽ cần ma trận biến dự đoán X, ma trận phản hồi Y và một ma trận tương quan mô tả sức mạnh liên hệ giữa mọi cặp phản hồi trong Y. Lưu ý rằng GFLASSO cho ra một ma trận p × k của β, khác với LASSO (p × 1), và ma trận hệ số này mang các liên hệ giữa bất kỳ phản hồi k nào và biến dự đoán j.

Bắt đầu

Kris Sankaran và tôi đã phát triển một gói R thử nghiệm triển khai GFLASSO cùng các phương pháp đánh giá chéo và trực quan hóa. Gần đây chúng tôi đã triển khai đa luồng bằng gói doParallel, giúp các quy trình đánh giá chéo (CV) nhanh hơn nhiều.

Để chạy GFLASSO trong R bạn cần cài devtools, nạp nó và cài gói gflasso từ kho GitHub của tôi. Phần minh họa sẽ thực hiện trên một bộ dữ liệu nằm trong gói bgsmtr. Tôi cũng khuyến nghị cài corrplot và pheatmap để trực quan hóa kết quả.

# 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)Được hỗ trợ bởi

Mô phỏng

Bạn có thể dễ dàng chạy mô phỏng được phác thảo trong trang trợ giúp cho hàm CV cv_gflasso(). Theo mặc định, CV tính căn sai số bình phương trung bình (RMSE) qua một lần lặp của CV 5-fold, trên tất cả các cặp có thể giữa λ ∈ {0, 0.1, 0.2, ..., 0.9, 1} và γ ∈ {0, 0.1, 0.2, ..., 0.9, 1}, tức lưới tinh chỉnh.

Lưu ý rằng các hàm sai số do người dùng cung cấp cũng hoạt động!

Ngoài các giả định thống kê vốn có và hiệu năng tốc độ, lựa chọn miền giá trị của lưới tinh chỉnh phụ thuộc nhiều vào việc chuẩn hóa trung bình và chia tỷ lệ phương sai đơn vị cho tất cả các cột trong X và Y, vì vậy hãy đảm bảo bạn làm điều đó trước.

Trong ví dụ sau bạn sẽ không cần, vì bạn sẽ lấy mẫu ngẫu nhiên từ phân phối chuẩn tắc. Bạn có thể thử suy ra hình phạt hợp nhất từ một mạng tương quan không trọng số, với ngưỡng 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)Được hỗ trợ bởi

Giá trị tối ưu của λ (hàng) và γ (cột) tối thiểu hóa RMSE trong mô phỏng này, lần lượt là 0.3 và 0.8, phản ánh đúng các quan hệ đã áp đặt.

Mẹo: Hãy thử chạy lại ví dụ này với một thước đo khác, hệ số xác định (R²). Một ưu điểm chính của R² là có miền từ 0 đến 1.

Hãy nhớ rằng khi bạn cung cấp một hàm độ phù hợp tự viết err_fun(), bạn phải xác định là tối đa hóa hay tối thiểu hóa thước đo kết quả bằng đối số err_opt.

Ví dụ sau nhằm tối đa hóa R², sử dụng mạng liên hệ có trọng số với hệ số tương quan bình phương (tức f(r_m**l)=r_m**l²). Nếu bạn có hơn 2 lõi, hãy đổi đối số nCores để tăng tốc!

# 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)
Được hỗ trợ bởi

Các tham số tối ưu λ và γ giờ là 0.6 và 0.3, tương ứng.

Cũng lưu ý rằng đối tượng cv_gflasso là một danh sách đơn với bốn phần tử: giá trị trung bình ($mean) và sai số chuẩn ($SE) của thước đo trên tất cả các ô của lưới, các tham số λ và γ tối ưu ($optimal) và tên của hàm độ phù hợp ($err_fun). Mô hình được đánh giá chéo trong ví dụ này rõ ràng ưa chuộng cả thưa (sparsity, λ) và hợp nhất (fusion, γ).

Cuối cùng, hãy nhớ bạn có thể tinh chỉnh thêm các tham số như ngưỡng hội tụ theo gradient Nesterov δ và số vòng lặp tối đa bằng cách truyền delta_conv và iter_max vào additionalOpts, tương ứng. Những tham số này sẽ được dùng trong ví dụ sau.

Xác định liên hệ SNP - hình ảnh thần kinh với GFLASSO

Để minh họa sự đơn giản và vững vàng của GFLASSO trong một bài toán tương đối nhiều chiều, tiếp theo bạn sẽ mô hình hóa các tập dữ liệu bgsmtr_example_data lấy từ cơ sở dữ liệu Alzheimer’s Disease Neuroimaging Initiative (ADNI-1).

Đây là một đối tượng danh sách gồm 3 phần, thuộc gói bgsmtr, bao gồm 15 thước đo cấu trúc hình ảnh thần kinh và 486 đa hình đơn nucleotide (SNP, chỉ thị di truyền) xác định từ mẫu 632 đối tượng. Quan trọng là 486 SNP này bao phủ 33 gen được xem là liên quan đến bệnh Alzheimer.

Nhiệm vụ của bạn là dự đoán các thước đo hình thái học, hình ảnh thần kinh từ dữ liệu SNP, khai thác cấu trúc tương quan của nhóm thước đo này.

Hãy bắt đầu bằng cách sắp xếp dữ liệu và khám phá phụ thuộc giữa tất cả các đặc trưng hình ảnh thần kinh:

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)
Được hỗ trợ bởi

Hình trên cho thấy các phụ thuộc giữa các đặc trưng hình ảnh thần kinh. Bây giờ hãy thử đánh giá chéo GFLASSO (có thể mất đến vài giờ trên laptop!) và xác định các liên hệ SNP - hình ảnh thần kinh.

Lưu ý rằng trong ví dụ dưới đây, đã chỉ định dung sai hội tụ và số vòng lặp tối đa. Bạn cứ thoải mái thử các giá trị khác!

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)
Được hỗ trợ bởi

Khi thử thách GFLASSO bằng các mô hình LASSO thuần (γ = 0, cột đầu), bình phương tối thiểu hợp nhất thuần (λ = 0, hàng đầu) và OLS (γ = 0 và λ = 0, ô trên cùng bên trái), ta có thể kết luận ví dụ hiện tại được mô hình hóa tốt nhất với trọng số hình phạt khác 0 và do đó, với GFLASSO đầy đủ. Lấy các tham số CV tối ưu (λ = 0.7 và γ = 1) để xây dựng mô hình GFLASSO và diễn giải ma trận hệ số thu được:

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)Được hỗ trợ bởi

Hình trên mô tả một tỷ lệ rất lớn các hệ số bằng 0 hoặc xấp xỉ 0. Dù không có sự gom cụm rõ ràng của các SNP theo gen (xem chú giải theo hàng, lưu ý chú giải chưa đầy đủ), vẫn có các liên hệ rõ ràng giữa một số SNP và đặc tính.

Để xác nhận sự tồn tại của một cơ chế dự đoán không ngẫu nhiên, bạn có thể lặp lại quy trình sau khi hoán vị các giá trị trong X hoặc Y. Thực nghiệm có thể giúp làm sáng tỏ quan hệ nhân quả và cơ chế từ các SNP được chọn. Ví dụ, các SNP ảnh hưởng đến trình tự và cấu trúc protein, làm gián đoạn việc dọn dẹp mảng bám β-Amyloid vốn là nền tảng của bệnh Alzheimer.

Tổng kết

GFLASSO sử dụng cả chuẩn hóa và hợp nhất khi mô hình hóa nhiều phản hồi, nhờ đó hỗ trợ nhận diện các liên hệ giữa biến dự đoán (X) và phản hồi (Y). Nó phù hợp nhất khi xử lý dữ liệu nhiều chiều với rất ít quan sát, vì chậm hơn đáng kể so với các phương pháp cạnh tranh. Ví dụ, các mô hình đồ thị Gaussian điều kiện thưa [4] và mô hình hồi quy đa nhiệm thưa theo nhóm Bayesian [5] có thể được ưa chuộng chủ yếu vì hiệu năng. Tuy nhiên, GFLASSO có tính diễn giải cao. Gần đây tôi đã dùng GFLASSO trong một cách tiếp cận tích hợp omics để khám phá các gen lipid mới ở ngô [6].

Hãy xem thêm Hướng dẫn Chuẩn hóa: Ridge, Lasso và Elastic Net của DataCamp.

Kris và tôi rất mong nhận phản hồi từ bạn. Dự án hiện được Kris duy trì tại kho krisrs1128/gflasso và thông qua kho của tôi, dù thay đổi thường xuyên, tại monogenea/gflasso. Hãy viết cho tôi bất cứ lúc nào ([email protected]), mọi góp ý đều được trân trọng.

Chúc bạn mã hóa vui vẻ!

Tài liệu tham khảo

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.