GFLASSO: Graph-Guided Fused LASSO en R

Aunque el campo del aprendizaje automático avanza a gran velocidad con técnicas cada vez más sofisticadas, se ha prestado poca atención a los problemas multitarea de alta dimensión que exigen predecir varias respuestas de forma simultánea. Este tutorial te mostrará el potencial de Graph-Guided Fused LASSO (GFLASSO) para predecir múltiples respuestas dentro de un único marco de regresión lineal regularizada.

Introducción

En el aprendizaje supervisado, el objetivo habitual es predecir una variable dependiente o de respuesta a partir de un conjunto de variables explicativas o predictores sobre una serie de muestras u observaciones. Los métodos de regularización introducen penalizaciones que evitan el sobreajuste en datos de alta dimensión, especialmente cuando el número de predictores supera al de observaciones. Estas penalizaciones se añaden a la función objetivo para que las estimaciones de los coeficientes de los predictores no informativos (los que apenas contribuyen a minimizar el error) se reduzcan por sí mismas. El least absolute shrinkage and selection operator (LASSO) [1] es uno de estos métodos.

¿Qué es LASSO?

Frente a los mínimos cuadrados ordinarios (OLS), LASSO puede encoger estimaciones de coeficientes (β) exactamente a cero, descartando así predictores no informativos y realizando selección de variables, mediante

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

donde n y j denotan, respectivamente, una observación y un predictor. La suma de cuadrados de residuos (RSS), único término usado en OLS, puede escribirse de forma equivalente como $RSS = \\sum_n(y_n-\\hat{y_n})^2 = (y-X\\beta)^T.(y-X\\beta)$. La penalización de LASSO es λ∑_j|β_j|, la norma L1 de los coeficientes ponderada por λ.

¿Por qué Graph-Guided Fused LASSO (GFLASSO)?

¿Y si necesitas predecir varias respuestas relacionadas a la vez a partir de un conjunto común de predictores? Podrías ajustar múltiples modelos LASSO independientes, uno por respuesta, pero te irá mejor si coordinas esas predicciones según la fuerza de asociación entre respuestas. Esta coordinación elimina la variación específica de cada respuesta que incluye ruido: esa es la clave de GFLASSO.

Un buen ejemplo aparece en el artículo original, donde los autores estudian las asociaciones entre 34 marcadores genéticos y 53 rasgos de asma en 543 pacientes [2].

¿Qué es GFLASSO?

Sea X una matriz de tamaño n × p, con n observaciones y p predictores, y Y una matriz de tamaño n × k, con las mismas n observaciones y k respuestas; por ejemplo, 1390 compras de electrónica distintas en 73 países para predecir las valoraciones de 50 producciones de Netflix en esos 73 países.

Entre los modelos adecuados para relacionar pares de conjuntos de datos de alta dimensión están Partial Least Squares bidireccional ortogonal (O2PLS), Canonical Correlation Analysis (CCA) y Co-Inertia Analysis (CIA), todos basados en descomposición matricial [3]. Además, como se basan en variables latentes (proyecciones de los predictores originales), la eficiencia computacional se logra a costa de la interpretabilidad.

Sin embargo, este intercambio no siempre compensa y puede evitarse con la predicción directa de k respuestas individuales a partir de características seleccionadas de X, dentro de un marco de regresión unificado que tiene en cuenta las relaciones entre respuestas.

En términos matemáticos, GFLASSO toma la regularización de LASSO [1] comentada arriba y construye el modelo sobre la estructura de dependencia en grafo subyacente en Y, cuantificada por la matriz de correlación k × k (la \"fuerza de asociación\" que mencionábamos). Como resultado, respuestas similares (o disimilares) se explicarán mediante subconjuntos similares (o disimilares) de predictores seleccionados.

Más formalmente, y siguiendo la notación del manuscrito original [2], la función objetivo de GFLASSO es

donde, sobre las k respuestas, ∑_k(y_k − X**β_k)^T.(y_k − X**β_k) aporta el RSS y λ∑_k∑_j|β_j**k| la penalización de regularización heredada de LASSO, ponderada por λ y aplicada a los coeficientes β de cada predictor j. La novedad de GFLASSO reside en

la penalización de fusión, ponderada por γ, que hace que la diferencia absoluta entre los coeficientes β_j**m y β_j**l, de cualquier predictor j y par de respuestas m y l, sea menor (o mayor) cuanto más positiva (o más negativa) sea su correlación por pares, transformada o no, f(r_m**l). Esta penalización de fusión favorece la variación globalmente significativa en las respuestas frente al ruido propio de cada una. Cuando la correlación por pares es cercana a cero, no actúa, y te quedas con un LASSO puro. Esta estructura de correlación para las k respuestas, que puede representarse como una red ponderada, por defecto es la correlación absoluta, f(r_m**l)=|r_m**l|, pero puede transformarse para crear variantes de GFLASSO con cualquier función definida por el usuario, por ejemplo:

\n
1. Correlación al cuadrado, f(r_m**l)=r_m**l² (ponderada)
\n
2. Correlación umbralizada, $f(r_{ml}) = \\begin{cases} 1, & \\mbox{if } r_{ml} > \\tau \\\\ 0, & \\mbox{otherwise} \\end{cases}$ (no ponderada)
\n

con mucho margen para innovar. Aunque la opción 2 es bastante menos costosa computacionalmente que la 1 y que la correlación absoluta por defecto [2], requiere fijar un umbral previo, por ejemplo, τ = 0.8.

En resumen, para ajustar un modelo GFLASSO necesitarás una matriz de predictores X, una matriz de respuestas Y y una matriz de correlación que represente la fuerza de asociación entre todos los pares de respuestas en Y. Ten en cuenta que GFLASSO produce una matriz p × k de β, a diferencia de LASSO (p × 1), y que esta matriz de coeficientes contiene las asociaciones entre cada respuesta k y predictor j.

Primeros pasos

Kris Sankaran y yo hemos estado trabajando en un paquete de R experimental que implementa GFLASSO junto con métodos de validación cruzada y de visualización. Recientemente hemos incorporado multihilo con el paquete doParallel, acelerando mucho las rutinas de validación cruzada (CV).

Para ejecutar GFLASSO en R necesitarás instalar devtools, cargarlo e instalar el paquete gflasso desde mi repositorio de GitHub. La demostración se hará con un conjunto de datos incluido en el paquete bgsmtr. También te recomiendo instalar corrplot y pheatmap para visualizar los resultados.

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

Simulación

Puedes ejecutar fácilmente la simulación descrita en la ayuda de la función de CV cv_gflasso(). Por defecto, la CV calcula el error cuadrático medio de raíz (RMSE) a lo largo de una única repetición de una CV de 5 particiones, sobre todos los pares posibles entre λ ∈ {0, 0.1, 0.2, ..., 0.9, 1} y γ ∈ {0, 0.1, 0.2, ..., 0.9, 1}, la malla de ajuste.

Nota: ¡también funcionarán funciones de error definidas por el usuario!

Más allá de los supuestos estadísticos y del rendimiento, la elección de los rangos de la malla de ajuste depende en gran medida de centrar a media y escalar a varianza unitaria todas las columnas de X y Y, así que asegúrate de hacerlo antes.

En el siguiente ejemplo no hará falta, ya que generarás muestras aleatorias de una distribución normal estándar. Puedes probar a derivar la penalización de fusión a partir de una red de correlación no ponderada, con un umbral de r > 0.8:

?cv_gflasso\nset.seed(100)\nX <- matrix(rnorm(100 * 10), 100, 10)\nu <- matrix(rnorm(10), 10, 1)\nB <- u %*% t(u) + matrix(rnorm(10 * 10, 0, 0.1), 10, 10)\nY <- X %*% B + matrix(rnorm(100 * 10), 100, 10)\nR <- ifelse(cor(Y) > .8, 1, 0)\nsystem.time(testCV <- cv_gflasso(X, Y, R, nCores = 1))\n\n## [1] 1.826146 1.819430 1.384595 1.420058 1.408619\n\n##    user  system elapsed \n##  45.808   3.070  55.271\n\nsystem.time(testCV <- cv_gflasso(X, Y, R, nCores = 2))\n\n## [1] 1.595413 1.492953 1.469917 1.366832 1.441642\n\n##    user  system elapsed \n##  22.231   1.698  26.441\n\ncv_plot_gflasso(testCV)

Los valores óptimos de λ (filas) y γ (columnas) que minimizan el RMSE en esta simulación, 0.3 y 0.8 respectivamente, capturan las relaciones impuestas.

Truco: prueba a repetir este ejemplo con otra métrica, el coeficiente de determinación (R²). Una ventaja clave de usar R² es que varía entre 0 y 1.

Recuerda que, si proporcionas una función de bondad de ajuste personalizada err_fun(), debes indicar si quieres maximizar o minimizar la métrica resultante mediante el argumento err_opt.

El siguiente ejemplo busca maximizar R², usando una red de asociación ponderada con coeficientes de correlación al cuadrado (es decir, f(r_m**l)=r_m**l²). Si tienes más de 2 núcleos, cambia el argumento nCores y dale más impulso.

# Write R2 function\nR2 <- function(pred, y){\n      cor(as.vector(pred), as.vector(y))**2\n}\n\n# X, u, B and Y are still in memory\nR <- cor(Y)**2\n\n# Change nCores, if you have more than 2, re-run CV\ntestCV <- cv_gflasso(X, Y, R, nCores = 2, err_fun = R2, err_opt = \"max\")\n\n## [1] 0.6209191 0.7207394 0.7262781 0.7193907 0.6303187\n\ncv_plot_gflasso(testCV)\n

Los parámetros óptimos λ y γ ahora son 0.6 y 0.3, respectivamente.

Además, ten en cuenta que los objetos cv_gflasso son listas con cuatro elementos: la media ($mean) y el error estándar ($SE) de la métrica sobre todas las celdas de la malla, los parámetros óptimos λ y γ ($optimal) y el nombre de la función de bondad de ajuste ($err_fun). El modelo validado del ejemplo presente favorece claramente tanto la esparsidad (λ) como la fusión (γ).

Por último, recuerda que puedes afinar parámetros adicionales, como el umbral de convergencia del gradiente de Nesterov δ y el número máximo de iteraciones, pasando delta_conv e iter_max a additionalOpts, respectivamente. Los usaremos en el siguiente ejemplo.

Determinación de asociaciones SNP–neuroimagen con GFLASSO

Para demostrar la sencillez y robustez de GFLASSO en un problema relativamente de alta dimensión, a continuación modelarás los conjuntos de datos bgsmtr_example_data obtenidos de la base de datos Alzheimer’s Disease Neuroimaging Initiative (ADNI-1).

Es un objeto lista de 3 elementos, parte del paquete bgsmtr, que contiene 15 medidas estructurales de neuroimagen y 486 polimorfismos de nucleótido único (SNPs, marcadores genéticos) determinados en una muestra de 632 sujetos. Importante: los 486 SNPs cubren 33 genes considerados asociados a la enfermedad de Alzheimer.

Tu tarea es predecir las medidas morfológicas de neuroimagen a partir de los datos de SNP, aprovechando la estructura de correlación de las primeras.

Empecemos organizando los datos y explorando las interdependencias entre todas las características de neuroimagen:

data(bgsmtr_example_data)\nstr(bgsmtr_example_data)\n\n## List of 3\n##  $ SNP_data     : int [1:486, 1:632] 2 0 2 0 0 0 0 1 0 1 ...\n##   ..- attr(*, \"dimnames\")=List of 2\n##   .. ..$ : chr [1:486] \"rs4305\" \"rs4309\" \"rs4311\" \"rs4329\" ...\n##   .. ..$ : chr [1:632] \"V1\" \"V2\" \"V3\" \"V4\" ...\n##  $ SNP_groups   : chr [1:486] \"ACE\" \"ACE\" \"ACE\" \"ACE\" ...\n##  $ BrainMeasures: num [1:15, 1:632] 116.5 4477.9 28631.9 34.1 -473.4 ...\n##   ..- attr(*, \"dimnames\")=List of 2\n##   .. ..$ : chr [1:15] \"Left_AmygVol.adj\" \"Left_CerebCtx.adj\" \"Left_CerebWM.adj\" \"Left_HippVol.adj\" ...\n##   .. ..$ : chr [1:632] \"V1\" \"V2\" \"V3\" \"V4\" ...\n\n# Transpose, so that samples are distributed as rows, predictors / responses as columns\nSNP <- t(bgsmtr_example_data$SNP_data)\nBM <- t(bgsmtr_example_data$BrainMeasures)\n\n# Define dependency structure\nDS <- cor(BM)\n\n# Plot correlation matrix of the 15 neuroimaging measures\ncorrplot(DS)\n

La figura anterior muestra las interdependencias entre las características de neuroimagen. Ahora valida mediante CV el GFLASSO (¡en un portátil puede tardar un par de horas!) y determina las asociaciones SNP–neuroimagen.

Nota: en el ejemplo de abajo se especifican la tolerancia de convergencia y el número máximo de iteraciones. ¡No dudes en probar valores distintos!

system.time(CV <- cv_gflasso(X = scale(SNP), Y = scale(BM), R = DS, nCores = 2,\n                 additionalOpts = list(delta_conv = 1e-5, iter_max = 1e5)))\n\n## [1] 1.550294 1.471637 1.470133 1.514425 1.504215\n\n##      user    system   elapsed \n##  2611.492   323.441 51129.936\n\ncv_plot_gflasso(CV)\n

Comparando GFLASSO con modelos LASSO puros (γ = 0, primera columna), una fusión de mínimos cuadrados pura (λ = 0, primera fila) y OLS (γ = 0 y λ = 0, celda superior izquierda), podemos concluir que este ejemplo se modela mejor con penalizaciones no nulas y, por tanto, con el GFLASSO completo. Toma los parámetros óptimos de CV (λ = 0.7 y γ = 1) para construir un modelo GFLASSO e interpretar la matriz de coeficientes resultante:

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

La figura muestra una proporción muy alta de coeficientes con valores cero o cercanos a cero. Aunque no hay un agrupamiento evidente de SNPs por genes (ver la anotación por filas; la leyenda está incompleta), sí se aprecian asociaciones claras entre ciertos SNPs y rasgos.

Para comprobar si existe un mecanismo predictivo no aleatorio, podrías repetir el procedimiento tras permutar los valores de X o de Y. Experimentos posteriores podrían ayudar a dilucidar causalidad y mecanismos a partir de SNPs seleccionados. Por ejemplo, SNPs que afectan a la secuencia y estructura de proteínas, alterando la eliminación de las placas de β-Amiloide implicadas en la enfermedad de Alzheimer.

Resumen

GFLASSO combina regularización y fusión al modelar múltiples respuestas, lo que facilita identificar asociaciones entre predictores (X) y respuestas (Y). Es especialmente útil con datos de alta dimensión y pocas observaciones, ya que es más lento que otros métodos. Por ejemplo, los modelos gráficos gaussianos condicionales dispersos [4] y el modelo bayesiano de regresión multitarea con grupos dispersos [5] pueden preferirse sobre todo por rendimiento. Aun así, GFLASSO es muy interpretable. Recientemente utilicé GFLASSO en un enfoque integrativo ómico para descubrir nuevos genes lipídicos en maíz [6].

Echa un vistazo al tutorial de regularización: Ridge, Lasso y Elastic Net de DataCamp.

A Kris y a mí nos encantará saber de ti. Este proyecto lo mantiene actualmente Kris en el repositorio krisrs1128/gflasso y también está disponible en el mío, aunque cambia con frecuencia, en monogenea/gflasso. Escríbeme cuando quieras (francisco.lima278@gmail.com), todo feedback es bienvenido.

¡Feliz programación!

Referencias

\n
1. Robert Tibshirani (1994). Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society, 58, 267-288.
\n
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.
\n
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.
\n
4. Lingxue Zhang, Seyoung Kim (2014). Learning Gene Networks under SNP Perturbations Using eQTL Datasets. PLoS Comput Biol, 10, 2:e1003420.
\n
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.
\n
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.
\n