Feature Selection for Data Integration with Mixed Multi-view Data

TL;DR

This paper introduces B-RAIL, a novel feature selection method designed for high-dimensional multi-view data with mixed types, improving data integration and network inference in complex biological datasets.

Contribution

We develop B-RAIL, a practical, theoretically-guided feature selection method that effectively handles heterogeneity in multi-view data for sparse regression and graph analysis.

Findings

B-RAIL outperforms existing methods in simulations.

B-RAIL identifies known ovarian cancer biomarkers.

B-RAIL suggests new potential biomarkers.

Abstract

Data integration methods that analyze multiple sources of data simultaneously can often provide more holistic insights than can separate inquiries of each data source. Motivated by the advantages of data integration in the era of "big data", we investigate feature selection for high-dimensional multi-view data with mixed data types (e.g. continuous, binary, count-valued). This heterogeneity of multi-view data poses numerous challenges for existing feature selection methods. However, after critically examining these issues through empirical and theoretically-guided lenses, we develop a practical solution, the Block Randomized Adaptive Iterative Lasso (B-RAIL), which combines the strengths of the randomized Lasso, adaptive weighting schemes, and stability selection. B-RAIL serves as a versatile data integration method for sparse regression and graph selection, and we demonstrate the effectiveness of B-RAIL through extensive simulations and a case study to infer the ovarian cancer gene regulatory network. In this case study, B-RAIL successfully identifies well-known biomarkers associated with ovarian cancer and hints at novel candidates for future ovarian cancer research.

Equations

$$\underset{\alpha ,{\mathbf{\beta }}_1,\dots {\mathbf{\beta }}_K}{\mathrm{minimize}}-\frac{1}{n}\ell \left(y;\alpha 1_n+\sum _{k=1}^K{X}_k{\mathbf{\beta }}_k\right)\text{subject to }\sum _{k=1}^K||{\mathbf{\beta }}_k|{|}_0\le \nu .$$

$$\hat{\alpha },\hat{\mathbf{\beta }}=\underset{\alpha \in R,\mathbf{\beta }\in R^p}{\mathrm{argmin}}-\frac{1}{n}\ell \left(y;\alpha 1_n+X\mathbf{\beta }\right)+\lambda ||\mathbf{\beta }|{|}_1,$$

$${\lambda }_{k,j}^{(t)}=\{\begin{array}{ll}{\eta }_k^{(t)}& \text{if }{\hat{\beta }}_{k,j}^{(t-1)}\ne 0,\\ 2{\eta }_k^{(t)}& \text{otherwise}\end{array}$$

$${\eta }_k^{(t)}=\frac{{\Lambda }_{\max }({\hat{\mathbf{\Theta }}}^{(t-1)})}{{\Lambda }_{\max }({\mathbf{X}}^T\mathbf{X})}\frac{1}{\sqrt{n}}\Vert {\hat{\mathbf{\beta }}}_k^{(t-1)}{\Vert }_2\sqrt{\frac{\log (p_k)}{n}\Vert {\hat{\mathbf{\beta }}}_k^{(t-1)}{\Vert }_0}$$

$${\hat{\mathbf{\beta }}}_k^{(t)}(b)=\underset{\alpha ,\mathbf{\beta }}{\mathrm{argmin}}-\frac{1}{n}\ell \left({\mathbf{y}}^{\ast b};\alpha +{\mathbf{X}}_k^{\ast b}\mathbf{\beta }+{\mathbf{\Phi }}_k^{(t)}(b)\right)+\sum _{j=1}^{p_k}{\gamma }_j{\lambda }_{k,j}^{(t)}|{\beta }_j|$$

$${\hat{S}}_k^{(t)}=\left\{j:\frac{1}{B}\sum _{b=1}^B\text{mathbbm}1\left({\hat{\beta }}_{k,j}^{(t)}(b)\ne 0\right)\ge \tau \right\}.$$

$$[{\hat{\mathbf{\beta }}}_k^{(t)}{]}_{{\hat{S}}_k^{(t)}}=\underset{\alpha ,\mathbf{\beta }}{\mathrm{argmin}}-\frac{1}{n}\ell \left(\mathbf{y};\alpha +{\left[{\mathbf{X}}_k\right]}_{{\hat{S}}_k^{(t)}}\mathbf{\beta }+{\mathbf{\Phi }}_k^{(t)}\right)+ϵ\Vert \mathbf{\beta }{\Vert }_2^2,$$

$${\lambda }_{k,j}^{(t)}=\{\begin{array}{ll}\eta & \text{if }{\hat{\beta }}_{k,j}^{(t-1)}\ne 0,\\ 2\eta & \text{otherwise}\end{array}$$

$$\eta =\underset{\begin{array}{c}(a)\text{ domain}\\ \text{correction}\end{array}}{\underbrace{\frac{{\Lambda }_{\max }({\hat{\mathbf{\Theta }}}^{(t-1)})}{{\Lambda }_{\max }(X^{T}X)}}}\underset{\begin{array}{c}(b)\text{ signal}\\ \text{correction}\end{array}}{\underbrace{\frac{1}{\sqrt{n}}||{\hat{\mathbf{\beta }}}_k^{(t-1)}|{|}_2}}\underset{\begin{array}{c}(c)\text{ Lasso}\\ \text{penalty}\end{array}}{\underbrace{\sqrt{\frac{\log (p_k)}{n}\Vert {\hat{\mathbf{\beta }}}_k^{(t-1)}{\Vert }_0}}}.$$

$$\hat{\alpha },\hat{\mathbf{\beta }}=\underset{\begin{array}{c}{\alpha }_1,\dots ,{\alpha }_K\in R,\\ {\mathbf{\beta }}_1,\dots ,{\mathbf{\beta }}_K\in R^p\end{array}}{\mathrm{argmin}}-\frac{1}{n}\sum _{k=1}^K\ell \left(y;{\alpha }_k{1}_n+X_k{\mathbf{\beta }}_k\right)+\sum _{k=1}^K\sum _{j=1}^{p_k}{\lambda }_{k,j}|{\mathbf{\beta }}_{k,j}|,$$

$${\hat{\mathbf{\beta }}}_k^{(t)}=\underset{\alpha ,\mathbf{\beta }}{\mathrm{argmin}}-\frac{1}{n}\ell \left(\mathbf{y};\alpha +{\mathbf{X}}_k\mathbf{\beta }+{\mathbf{\Phi }}_k^{(t)}\right)+\sum _{j=1}^{p_k}{\lambda }_{k,j}|{\beta }_j|$$

Full text

data fusion,

multi-modal data,

integrative genomics,

variable selection,

Lasso/GLM Lasso,

stability selection,

mixed graphical models,

journalname

, ,

1 Introduction

As the amount of data grows in volume and variety, data integration, or the analysis of multiple sources of data simultaneously, is becoming increasingly necessary in numerous disciplines. For example, in genomics, scientists can gather data from many related, yet distinct sources including gene expression, miRNA expression, point mutations, and DNA methylation. Since all of these genomic sources interact within the same biological system, it can be advantageous to analyze them together via data integration. Ultimately, the abundance and diversity of information captured by integrated data offers an invaluable opportunity to gain a better and more holistic understanding of the phenomena at hand.

In this work, we aim to perform feature selection for a common family of integrated data sets called high-dimensional multi-view data. Multi-view data refers to data collected on the same set of samples, but with features from multiple sources of potentially mixed types (e.g. categorical, binary, count, proportion, continuous, and skewed continuous values). More formally, suppose we observe multi-view data with $K$ high-dimensional views (or sources), $\operatorname{\mathbf{X}}_{1}\in\operatorname{\mathbb{R}}^{n\times p_{1}},\dots,\operatorname{\mathbf{X}}_{K}\in\operatorname{\mathbb{R}}^{n\times p_{K}}$, which are measured on the same $n$ samples but with features of mixed types. We seek to recover a sparse set of features from each $\operatorname{\mathbf{X}}_{k}$ associated with the response $\operatorname{\mathbf{y}}\in\operatorname{\mathbb{R}}^{n}$ by considering:

[TABLE]

Here, $\bm{\beta}_{k}\in\operatorname{\mathbb{R}}^{p_{k}}$ are the coefficients associated with view $k$, $\nu>0$ is a tuning parameter which regulates the sparsity level, and $\ell()$ is the generalized linear model (GLM) log-likelihood associated with $\operatorname{\mathbf{y}}$. Note we not only consider continuous (Gaussian) responses, but also the broader class of GLMs including the Poisson (log-linear) and Bernoulli (logistic) families.

While there are many applications for multi-view feature selection in genomics, imaging, national security, economics, and other fields, major difficulties, stemming from the heterogeneity of features and how to appropriately integrate such differences, have prevented the successful use of multi-view feature selection in practice. To our knowledge, no one has proposed an effective practical solution to perform feature selection with multi-view data. A plethora of works have studied feature selection in the high-dimensional setting via the Lasso or GLM Lasso (Tibshirani, 1996; Yuan and Lin, 2007; Tibshirani et al., 2013; Zhao and Yu, 2006), and others have studied various data integration problems (Hall and Llinas, 1997; Shen, Olshen and Ladanyi, 2009; Acar, Kolda and Dunlavy, 2011). However, there is limited research at the intersection of the two fields.

The one area that touches on multi-view feature selection is in the context of mixed graphical models, which estimate sparse graphs between features in multi-view data (Yang et al., 2014a, b; Lee and Hastie, 2013; Cheng et al., 2013; Haslbeck and Waldorp, 2015). Using the node-wise neighborhood estimation approach of Meinshausen and Bühlmann (2006), mixed graphical models estimate the neighborhood of each node (i.e. feature) separately via a penalized regression model (typically based on the Lasso or GLM Lasso) and combine neighborhoods using an “AND” or “OR” rule. Though mixed graphical models perform well in idealized settings for which theoretical guarantees have been proven, we will demonstrate in Section 2 that there are severe limitations with these approaches in realistic settings with correlated, heterogeneous features, commonly found in multi-view data.

To facilitate more effective integrative analyses in practice, we investigate the under-studied problem of high-dimensional multi-view feature selection, and we propose a practical solution. Our work is the first to identify and critically examine the fundamental challenges of multi-view feature selection, and we leverage this deep understanding of the challenges to develop a new high-dimensional multi-view selection method, the Block Randomized Adaptive Iterative Lasso (B-RAIL). B-RAIL is a practical tool for multi-view feature selection with its roots grounded in theory, and it builds upon adaptive $\ell_{1}$ penalties, the randomized Lasso, and stability selection (Meinshausen and Bühlmann, 2010) to overcome the issues incurred by existing methods. Our method can be used for both regression and mixed graphical selection, thus lending itself to a host of important applications.

In Section 2, we investigate the major challenges of multi-view feature selection and highlight the literature gaps relating to these issues. We also show that the culmination of these challenges lead to poor feature recovery in existing Lasso-type methods and mixed graphical models. In Section 3, we introduce our proposed method, B-RAIL, which takes steps to address the challenges from Section 2. In Section 4, we showcase the strong empirical performance of B-RAIL through simulations and contrast it to existing methods. In Section 5, we further demonstrate the effectiveness of B-RAIL in a novel integrative genomics case study for ovarian cancer, and we provide concluding remarks in Section 6.

2 Challenges

Before introducing our proposed method, it is instructive to understand the challenges posed by feature selection in the multi-view setting. These challenges have been over-looked in previous methods and thus contribute to many of their shortcomings. In this section, we focus on the challenges faced by linear models with Lasso-type penalties due to their overwhelming popularity and desirable statistical properties (Tibshirani, 1996; Yuan and Lin, 2007; Meinshausen and Yu, 2009; Tibshirani et al., 2013; Zhao and Yu, 2006; Zhang and Huang, 2008). Given data $\operatorname{\mathbf{X}}\in\operatorname{\mathbb{R}}^{n\times p}$ and response $\operatorname{\mathbf{y}}\in\operatorname{\mathbb{R}}^{n}$, recall that the (GLM) Lasso solves

[TABLE]

where $\lambda>0$ is a regularization parameter and $\ell()$ is the GLM log-likelihood associated with the response. For clarity, we use the term “Lasso” to refer to the $\ell_{1}$-penalized model with continuous (Gaussian) responses, “GLM Lasso” to mean the $\ell_{1}$-penalized model with non-Gaussian GLM responses (e.g. binary, Poisson), and “Lasso-type” methods to mean either the Lasso or GLM Lasso with some form of $\ell_{1}$ penalty (e.g. a global penalty, separate penalties, adaptive penalties).

Our focus here is not on deriving new theoretical guarantees for the Lasso in multi-view settings. Rather, we highlight deep practical concerns, which are rooted in theory and commonly arise in feature selection for data integration. By identifying these practical challenges, we open up numerous avenues for future theoretical research and set the stage for the construction of a new method, which overcomes the identified issues.

2.1 Motivating Example

To first illustrate the current challenges and motivate the need for a solution, we present in Figure 3 the estimated graphs from common Lasso-type methods and our proposed method when applied to real ovarian cancer genomics data. Here, there are $n=293$ samples and $p=836$ features from three views: count-valued RNASeq data ($p_{\text{RNASeq}}=408$), continuous miRNA data ($p_{\text{miRNA}}=307$), and proportion-valued methylation data ($p_{\text{Methyl}}=301$) (refer to Section 5 for data collection and pre-processing details). As in several previous graphical models and mixed graphical models (Meinshausen and Bühlmann, 2006; Ravikumar et al., 2010; Jalali et al., 2011), we estimated the graphs using node-wise neighborhood selection. We then combined neighborhoods using the “AND” rule and applied stability selection (Meinshausen and Bühlmann, 2010; Liu, Roeder and Wasserman, 2010) with the threshold $0.98$ to select stable edges.

Figure 3 specifically compares three types of estimation schemas at each node: (a) GLM Lasso with one global penalty, (b) GLM Lasso with separate penalties for each view, and (c) our proposed B-RAIL algorithm (introduced in Section 3). The first two methods have been proposed in several mixed graphical models (Yang et al., 2014a; Chen, Witten and Shojaie, 2014; Haslbeck and Waldorp, 2015) and satisfy strong theoretical guarantees in idealized settings. In the real data example however, the Lasso-type methods are unstable (illustrated by the fewer edges), favor feature selection within one view, and select only a few edges between views. This overall instability indicates that the Lasso-type methods are not robust to small perturbations of the data and raises serious concerns about the reproducibility and reliability of the results (Yu, 2013). Our proposed B-RAIL algorithm, in contrast, avoids these issues and exhibits greater stability as well as balance, selecting a larger number of within and between block edges under the same thresholding value. We will later see through extensive simulations in Section 4 that the issues with existing Lasso-type methods observed here are recurring problems in very general multi-view scenarios.

To begin understanding why existing Lasso-type methods struggle in practice, we identify and study four major challenges of feature selection for high-dimensional multi-view data: 1) scaling, 2) ultra-high-dimensionality, 3) signal interference, and 4) domain-specific beta-min. These issues stem from a combination of domain differences, signal differences, and the high-dimensionality of each view. Together, these challenges can have a significant adverse effect on feature recovery for data integration. We next examine each of these challenges in greater detail.

2.2 Scaling

The first and most obvious challenge with integrative analyses revolves around scaling. That is, each view in a multi-view data set is often measured on a different scale, and it is unclear how to most effectively integrate such differences. Many believe that normalizing all features to mean [math] and variance $1$ remedies the scaling differences, but this is not always the case. Even after centering and scaling, data views remain distinct if they differ in ways beyond the first and second moments. This issue is especially problematic with binary and count-valued data blocks, two common types in multi-view data, since they are defined by much higher moments. We thus highly discourage using the ordinary (GLM) Lasso with a single penalty (2.1) on normalized multi-view data.

Now while one can use different regularization parameters for each view to help alleviate the scaling differences, this generates another set of issues that are complicated by the following challenges. We will revisit the scaling issue in light of these complications later in this section.

2.3 Ultra-High-Dimensionality

In addition to the scaling issue, performing exact feature selection with the Lasso is already difficult in the ordinary high-dimensional setting. For exact feature selection, the number of samples $n$ must be above a theoretical minimum known as the sample complexity. In the highly idealized scenario of an iid standard Gaussian design and a Gaussian response, Wainwright (2009) showed that the sample complexity scales at approximately $2s\log(p-s)$, where $p$ is the number of features, and $s$ is the number of non-zero features. This idealized lower bound can be difficult to attain in many applications including genomics, where typical values of $p=1000$ and $s=30$ demand $n\approx 400$ patients - a large and highly expensive study. We informally refer to the regime where $p\gg n\geq 2s\log(p-s)$ as “high-dimensional” and $n<2s\log(p-s)$ as “ultra-high-dimensional.” Roughly, the Lasso can never perform exact feature selection in the ultra-high-dimensional regime.

For non-Gaussian responses and more realistic designs such as correlated, heterogeneous views in multi-view data, the sample complexity is significantly higher than the idealized Gaussian bound (Chen, Witten and Shojaie, 2014; Ravikumar et al., 2010). As an example, the Poisson GLM’s sample complexity scales at approximately $s^{2}\log(p(\log p)^{2})$ (Yang et al., 2015), so if $p=1000$ and $s=30$, we require $n\approx 10,000$ samples. This problem is further exacerbated in multi-view settings since combining multiple high-dimensional views for data integration almost always results in an ultra-high-dimensional problem.

2.4 Signal Interference

The third challenge we identify stems from a problem with the Lasso known as shrinkage noise. Su, Bogdan and Candes (2015) showed that with high probability, no matter how strong the effect sizes, false discoveries appear early on the Lasso path due to pseudo-noise introduced by shrinkage in the high- and ultra-high-dimensional regimes. When the Lasso selects its first few features using large regularization parameters, the residuals still contain much of the signal associated with the selected features, and it is this extra noise which Su, Bogdan and Candes (2015) calls shrinkage noise.

In the multi-view context, shrinkage noise becomes a very complex and serious issue due to the different signals across blocks. Since the Lasso naturally selects features from the block with the highest signal first, the resulting shrinkage noise will mask the weaker signals from other blocks and compromise our ability to select from these weaker blocks. We refer to this adverse consequence of shrinkage noise as signal interference.

The problem of shrinkage noise has not been widely studied beyond the iid Gaussian design in Su, Bogdan and Candes (2015), but we provide strong empirical evidence in Figure 5 that confirms the existence of shrinkage noise and signal interference in non-Gaussian multi-view settings. In the case of an iid Gaussian and an iid binary block, shown in Figure 5, the Lasso achieves perfect recovery in the Gaussian block when the signal-to-noise ratio (SNR) in the binary block is [math], but as the SNR of the binary block increases, it interferes with our ability to recover the Gaussian features in the small sample scenario of $n=200$. This signal interference is especially disastrous in the GLM Lasso with binary responses, where support recovery in the Gaussian block tends to [math]. When we increase the sample size to $n=300$ however, there is no decline in the recovery of the Gaussian block in Figure 5(a). This agrees with the known result from Su, Bogdan and Candes (2015) that shrinkage noise occurs when the Lasso’s theoretical conditions are violated and in particular, when $n$ is not sufficiently large.

2.5 Domain-Specific Beta-min Condition

Finally, analogous to how signal differences can exacerbate the Lasso’s shrinkage noise issue, domain differences in multi-view problems can complicate the Lasso’s beta-min condition, which establishes a lower bound for the minimum amount of signal (i.e. SNR) required for feature recovery (Zhao and Yu, 2006; Meinshausen and Bühlmann, 2006; Bühlmann et al., 2013).

In Figure 6, we report our ability to recover the true features for a simple simulation with iid features from four data types (Gaussian, binary, uniform, and Poisson). In each subplot, the solid vertical lines indicate the minimum SNR needed to recover $99\%$ of all true features when $p=500$ (red) and $p=5000$ (blue). We observe that the minimum SNR requirement varies based upon the domain of the features and that the Gaussian predictors can tolerate the lowest SNR. These empirical results reveal that if two blocks have the same amount of signal but are from different domains, we can only recover the features that pass the minimum signal threshold dictated by the domains. Put concretely, if we were to perform feature selection on our simulated multi-view data with $p=500$ and SNR $=0.66$, we would be able to recover $99\%$ of the true features in the Gaussian block but only about $3/4$ of the true features in the binary block. This observed phenomenon agrees with previous work, which has shown that an increase in the sparsity in $\operatorname{\mathbf{X}}$ effectively reduces the SNR in the high-dimensional setting (Wang, Wainwright and Ramchandran, 2010). Beyond this however, the beta-min condition has been relatively unexplored for the GLM Lasso and domain differences and remains a ripe area for future theoretical work.

2.6 Additional Challenges

It is important to note that the four challenges above do not act independently from one another. In fact, the main source of difficulty with multi-view feature selection is arguably the interactions between challenges. For instance, consider the problem of selecting features from high-dimensional discrete blocks with weak signals. The ultra-high-dimensionality issue can exacerbate the already existing problem of signal interference, which can then worsen scaling issues, increase minimum SNR requirements, and amplify the overall difficulty of the problem.

In conjunction with these complex interactions, the need to select an appropriate amount of regularization $\lambda$ through model selection methods can also increase the difficulty of multi-view feature selection. We will compare three common selection methods, namely, stability selection (Meinshausen and Bühlmann, 2010; Liu, Roeder and Wasserman, 2010), cross validation (Stone, 1974; Allen, 1974; Shao, 1993), and extended BIC (Chen and Chen, 2012), and discuss their additional challenges in Section 4.

We lastly note that the majority of our discussion has been focused on the Lasso. Feature selection is even more challenging for the GLM Lasso. Chen, Witten and Shojaie (2014) investigated this for mixed graphical models and concluded that the predictors associated with Gaussian responses are easier to recover than those with responses from other exponential families. Specifically, Gaussian responses require fewer samples, allow for a wider range of tuning parameters, and generally have a higher probability of success.

2.7 Challenges with Existing Methods

Having identified a host of challenges, we return to address why common Lasso-type methods are not well-suited for multi-view feature selection.

To begin with its most simple form, the Lasso with a single global penalty (2.1) uses the same penalty for all views and does not alleviate the scaling issues or signal interference issues in multi-view data. The consequences of these problems are evident, especially in the case of non-Gaussian blocks with weak signals, in Figure 3(a), where fewer edges are selected within the proportion-valued methylation block.

By employing the Lasso with separate penalties for each data view, we can mitigate the issue of scaling. Nevertheless, model selection becomes more challenging with multiple penalties, and signal interference remains a driver of poor recovery. In fact, having a separate penalty for each view exacerbates signal interference and encourages selection from the block with the strongest signal and no selection from the blocks with weak signals. This signal interference is exemplified in Figure 3 by the extreme selection imbalance among views, with almost no selection in the miRNA block and heavy selection in the RNASeq block.

In the Adaptive Lasso (Zou, 2006), the amount of $\ell_{1}$-regularization associated with $\beta_{j}$ is typically $\lambda/|\hat{\beta}_{j}^{OLS}|^{\gamma}$ for some constants $\gamma,\lambda>0$. This adaptive penalty mitigates the scaling issue by adjusting for signal differences through $\hat{\beta}_{j}^{OLS}$, but it also encourages selection of features with higher signals and penalizes features with weaker signals. The Adaptive Lasso hence complicates signal interference by treating weaker signals as noise and results in little to no selection in the blocks with weak signals.

While the previous methods all struggle with signal interference, one simple way to reduce the signal interference between blocks is to perform separate Lassos for each data view. Since independently-estimated blocks cannot possibly interfere with one another, this method addresses both scaling and signal interference issues. It also avoids the problem of ultra-high-dimensionality. However, each view by itself usually does not contain sufficient information to explain much of the variability in the response, and we lose the advantages of data integration.

Beyond the Lasso-type methods, there are selection methods with non-convex penalties such as SCAD and MCP (Fan and Li, 2001; Zhang et al., 2010). These non-convex penalties tend to scale better than the Lasso-type penalties but are still not variable selection consistent in the ultra-high-dimensional regime, especially for non-Gaussian responses and highly correlated data. We investigate MCP/SCAD feature selection in Table 6 in the Appendix, but our primary focus in this paper is on the more commonly used Lasso-type penalties.

3 Block Randomized Adaptive Iterative Lasso

Driven by the many challenges and the lack of effective tools, we propose a new method for multi-view feature selection, the Block Randomized Adaptive Iterative Lasso (B-RAIL). For the sake of notation, suppose we observe the response vector $\operatorname{\mathbf{y}}$ and multi-view data $\operatorname{\mathbf{X}}=[\operatorname{\mathbf{X}}_{1},\dots,\operatorname{\mathbf{X}}_{K}]$ with $K$ views of potentially mixed types, $n$ samples, and $p$ total features. We will assume $p\gg n$ and typically $p_{k}\gg n$ for each view. Let $S$ denote the indices of the support, and let $[\operatorname{\mathbf{X}}]_{S}$ denote the columns of $\operatorname{\mathbf{X}}$ indexed by $S$. We will introduce B-RAIL in the context of regression and later discuss its extension to graph selection.

Under the regression framework, the goal of B-RAIL can be viewed as two-fold: 1) to select features from each view $\operatorname{\mathbf{X}}_{k}$ that are associated with the response $\operatorname{\mathbf{y}}$, and 2) to do so while avoiding the challenges discussed in Section 2. With this goal in mind, we briefly outline the B-RAIL algorithm in Algorithm 1 and summarize the key steps taken to overcome the current challenges.

At a high-level, B-RAIL iterates across the data blocks $k=1,\dots,K$ and estimates $\hat{\bm{\beta}}_{k}$ for each data block $\operatorname{\mathbf{X}}_{k}$ separately while holding all other blocks fixed. This iterative procedure is motivated by the advantages of performing separate Lassos - namely, that it mitigates the ultra-high-dimensionality and signal interference issues. Then, within each of the individual block estimations, B-RAIL first estimates the block’s support and subsequently, the coefficient values given the support. Here, B-RAIL leverages ideas from adaptive weighting schemes, stability selection, and the randomized Lasso in an attempt to reduce the scaling discrepancies and domain-specific beta-min issues.

We next provide the full B-RAIL algorithm in Algorithm 2 and proceed to discuss each step of the B-RAIL algorithm in greater detail.

3.1 Initialization

In our proposed B-RAIL algorithm, the coefficients are first initialized to a pre-specified sparsity level (e.g. $0.2p_{k}$ non-zero features per block) by fitting separate Lasso path regressions for each block. As long as the algorithm is initialized to an over-selection of the support of $\bm{\beta}$, we have found in all of our empirical simulations, the B-RAIL algorithm tends to perform well and is very robust to the exact choice of initialization.

After initializing $\bm{\beta}$, we must specify the order of the blocks to iterate over. This ordering can slightly alter the estimation results of B-RAIL as accurate estimation of previous blocks makes subsequent estimations of other blocks easier, but in most cases, we have found that the block ordering is not as important to B-RAIL’s performance as initializing the coefficients to an over-selection of the support. Nevertheless, for best practices, since previous Lasso results guarantee a high probability of support recovery when $n$ is sufficiently large compared to $p$, we advise estimating the blocks with the smallest $p$ first, especially if $p\leq n$. If dimensions of all the blocks are of similar sizes or much larger than $n$, we recommend starting with Gaussian blocks, which tend to have better support recovery than non-Gaussian blocks.

3.2 Estimating Support ($\hat{S}^{(t)}_{k}$)

After initialization, we repeatedly iterate across the $K$ data blocks and estimate the support of each block separately, holding the estimates of all other blocks fixed. This block-wise estimation avoids the ultra-high-dimensionality issue, and because shrinkage noise is mainly a problem in the ultra-high-dimensional regime, the signal interference issue is also mitigated as a direct bi-product.

Furthermore, to effectively handle correlated features in practice, we incorporate stability selection with the randomized Lasso (Meinshausen and Bühlmann, 2010) to estimate each block. As given by step 1(b) in Algorithm 2, we solve the Lasso $B$ times using the bootstrap and randomized penalty terms, and we threshold the stability score at $\tau$ (3.4) to select the most stable features. Though $\tau\in(0,1)$ is a user-specified hyperparameter, the B-RAIL algorithm is insensitive to choices of $\tau$ within reasonable ranges. This insensitivity to $\tau$ has also been observed in previous work on stability selection (Meinshausen and Bühlmann, 2010). Ultimately, by utilizing randomized penalties and stability selection when estimating the support of each block, B-RAIL leverages the key property that the randomized Lasso is feature selection consistent even when the Lasso’s irrepresentable condition is violated (Meinshausen and Bühlmann, 2010) and hence can effectively handle correlated features.

3.3 Adaptive Regularization ($\lambda$)

Now, looking more closely at the penalty term in (3.3) of the randomized Lasso, the penalty term includes a random weight $\gamma$ like the original randomized Lasso. However, in order to account for the scaling discrepancies, signal variability, and domain differences between blocks, we introduce a block-specific adaptive penalty $\lambda$ in (3.3) as well. For feature $j$ in block $k$, we define the adaptive weight

[TABLE]

where

[TABLE]

Here, $\widehat{\bm{\Theta}}^{(t-1)}$ is the Fisher information matrix corresponding to the GLM of the response $\operatorname{\mathbf{y}}$, and $\Lambda_{\max}$ denotes the maximum eigenvalue.

In this definition of $\lambda$, there are two moving parts. First, the multiplicative scheme in (3.1) encourages previously selected features to remain selected while still allowing all features to freely enter or exit the model. Secondly, $\eta$ accounts for the heterogeneity of multi-view data and helps to mitigate the challenges of Section 2.

Though the exact form of $\eta$ was derived experimentally, $\eta$ can be interpreted as the product of three factors, each of which is rooted in solid theoretical foundations. Namely, part (c) of (3.2) is closely related to the theoretical bound on the regularization parameter needed for selection consistency of the Lasso (Zhao and Yu, 2006; Meinshausen and Bühlmann, 2006). The ratio of eigenvalues in part (a) (i.e. the domain correction term) is motivated by the theoretical conditions imposed on the Fisher information matrix for exponential family distributions (Yang et al., 2015), and part (b) of (3.2) (i.e. the signal correction term) can be viewed as the average signal in block $k$ since $\frac{1}{\sqrt{n}}$ is derived from the theoretical sparsity level within each block (Bunea et al., 2007).

By constructing the adaptive penalty $\eta$ in this way, B-RAIL accounts for different block sizes through the $\frac{\log(p_{k})}{n}$ term and automatically penalizes non-Gaussian blocks less heavily than Gaussian blocks since $\Lambda_{\text{max}}(\hat{\Theta})$ is larger for Gaussian blocks. This helps to balance the inherently different beta-min conditions. In addition, because $||\hat{\bm{\beta}}_{k}||_{2}$ captures information about both the signal and scale of the $k^{th}$ block, $\eta$ addresses the scaling differences and penalizes the stronger signal blocks more heavily to allow for the possibility of selection from weaker blocks.

While in theory this specific combination of weights should correct for scaling and domain-specific beta-min differences across views, we reinforce our choice of $\eta$ through strong empirical results in Section 4. We also note that even if the form of $\eta$ is slightly misspecified, stability selection is known to be fairly robust to the exact amount of regularization, as long as the amount of regularization is within reason (Meinshausen and Bühlmann, 2010). Incorporating stability selection with the randomized Lasso thus serves as a built-in check within B-RAIL, which is advantageous in practice.

3.4 Coefficient Estimations

After estimating the block-wise support using the randomized Lasso with adaptive weights, we seek to estimate the coefficients of the support as accurately as possible since these values are used in future block estimations and iterations of B-RAIL. We hence refit a penalized regression model with a small ridge penalty (e.g., $\epsilon\approx 10^{-4}$) in (3.5) to avoid the known bias issues with the Lasso. The only reason to include the tiny ridge penalty is to ensure that we can still estimate coefficients when the selected support is greater than $n$. The exact choice of $\epsilon$ has a negligible impact in practice because it is chosen to be so small.

3.5 Convergence

We finally declare convergence of B-RAIL’s iterative block estimation procedure when the estimated support remains unchanged. Our empirical analysis indicates that B-RAIL has quick support convergence, and we provide one example of this fast convergence in Figure 14 in the Appendix. Using the ovarian cancer simulation (see Section 5) for three different responses (Gaussian, binary, and Poisson), we report that the average number of iterations until convergence is between 4 and 5 with the maximum number of iterations reaching 15 (over 100 runs). These ranges are similar for all designs, empirically demonstrating B-RAIL’s fast convergence.

Though convergence of the full-blown B-RAIL algorithm is currently limited to empirical analysis, B-RAIL can be viewed as a block coordinate descent algorithm, which can be studied theoretically under some simplifications to gain additional insights into the full B-RAIL algorithm. Namely, if we omit the adaptive regularization parameter and apply the ordinary (GLM) Lasso in each block of the algorithm (as detailed in Algorithm 3 in the Appendix), we call the resulting algorithm the blockwise (GLM) Lasso and discuss its convergence below.

Proposition 1**.**

Consider the optimization problem:

[TABLE]

where $\hat{\alpha}=[\hat{\alpha}_{1},\dots,\hat{\alpha}_{K}]$, $\hat{\bm{\beta}}=[\hat{\bm{\beta}}_{1},\dots,\hat{\bm{\beta}}_{K}]$, and $\lambda_{k,j}\geq 0$, and suppose that the objective function in (3.6) is bounded below. Then the blockwise (GLM) Lasso converges to a global solution of(3.6).

To prove Proposition 1, we leverage the block coordinate descent view of the blockwise Lasso and apply Theorem 4.1 from Tseng (2001) since the objective function is convex and separable with respect to each block $\bm{\beta}_{k}$. The detailed proof is provided in the Appendix.

The two main differences between the blockwise Lasso and B-RAIL are the adaptive regularization parameter $\eta$ and the use of stability selection with the randomized Lasso in each block’s update. If the estimated support from stability selection converges to a common support across the many block iterations in B-RAIL, then B-RAIL reduces to the blockwise Lasso algorithm, for which we have shown convergence. While there is empirical evidence to believe that stability selection, applied iteratively with the adaptive regularization parameter as in B-RAIL, converges to a common support, proving this theoretically is challenging due to the purely algorithmic and random nature of stability selection. In addition, existing optimization-theoretic frameworks cannot handle adaptive parameters ($\eta$) that are dependent on previous iterates ($\hat{\bm{\beta}}^{(t-1)}$). Developing such a framework to handle both of these issues is beyond the scope of this work, but we plan to further investigate it in the future. For now, due to these serious difficulties, we rely on our empirical analysis to demonstrate B-RAIL’s quick convergence.

3.6 B-RAIL Summary

While we have introduced B-RAIL under the regression framework, B-RAIL can be naturally extended to estimate mixed graphical models via a penalized node-wise regression approach (Meinshausen and Bühlmann, 2006). As in the motivating example in Section 2, we can use B-RAIL to estimate the neighborhood of each node separately via penalized regressions and then combine the neighborhoods using an“AND” or “OR” rule to obtain the graph.

In either the regression or graph selection setting, our B-RAIL algorithm deliberately takes steps to exploit the practical advantages of existing Lasso-type methods while avoiding the drawbacks described in Section 2. For instance, by performing iterative block-by-block estimations, B-RAIL inherits the advantages of performing separate Lassos and avoids the issue of ultra-high-dimensionality. This in turn reduces signal interference between blocks since shrinkage noise is only a concern when $n$ is not sufficiently large relative to $p$. Furthermore, we mitigate the scaling and beta-min problems by engineering adaptive $\ell_{1}$ penalties in B-RAIL to correct for domain and signal differences between blocks. In this construction, slightly weaker non-Gaussian blocks are penalized less heavily and thus not completely overshadowed by Gaussian blocks. Still, selecting an appropriate amount of $\ell_{1}$ regularization is challenging in practice, especially due to highly correlated data. B-RAIL thus incorporates randomized stability selection, which is known to be feature selection consistent under stronger and more complex dependencies than can be handled by the Lasso. This boosts the support estimation of correlated features, and together, with the previous components, B-RAIL effectively overcomes the many practical challenges of multi-view feature selection and lends itself to a plethora of data integration applications.

4 Numerical Studies

We next reinforce the theoretically-guided choices in our B-RAIL construction and demonstrate its effectiveness through extensive simulations. In these simulations, we evaluate B-RAIL against four common Lasso-based parametric methods: (i) Lasso with a global penalty for all blocks, (ii) Lasso with separate penalties for each block, (iii) separate Lassos for each block, and (iv) Adaptive Lasso. For the Adaptive Lasso, we use ridge weights as they are better adapted to handle correlated features. Moreover, to avoid biases from penalty selection methods, we use oracle information to select features in the Lasso-based models. That is, if $k$ is the number of true features in the simulation, we fit the full path of the Lasso and select the first $k$ features. In the case of the Lasso with separate penalties, we select the $k$ features with the largest number of true positives. We do not, however, use oracle information for B-RAIL. Instead, B-RAIL internally selects the number of features using stability selection with the threshold $\tau=0.8$ as outlined in Algorithm 2, and we set $\epsilon=0.001/p$.

To systematically compare these methods, we simulate from various designs of $\operatorname{\mathbf{X}}$ with three blocks - namely, a Gaussian $\operatorname{\mathbf{X}}_{1}$, Bernoulli $\operatorname{\mathbf{X}}_{2}$, and Poisson $\operatorname{\mathbf{X}}_{3}$ block - and various types of GLM responses $\operatorname{\mathbf{y}}$. Due to the popular use of the Gaussian, Bernoulli, and Poisson GLMs, we run simulations with responses $\operatorname{\mathbf{y}}$ from each of these families. For the Gaussian response, we fit the linear model $\operatorname{\mathbf{y}}=\operatorname{\mathbf{X}}\beta+\epsilon$, where $\epsilon\sim N(0,1)$. For the binary and Poisson responses, we use copula transformations (Nelsen, 1999) to simulate $\operatorname{\mathbf{y}}$.

In addition to these response models for $\operatorname{\mathbf{y}}$, we consider four simulation designs for $\operatorname{\mathbf{X}}$ to understand model behavior under different assumptions. The four simulations designs are: (i) iid features, (ii) independent features with non-constant variance, (iii) correlated features with covariance structure from a Block Directed Markov Random Field, and (iv) a real data inspired simulation with features from The Cancer Genome Atlas (TCGA) ovarian cancer study. We elaborate on each of these designs below.

Note in all of the simulations, we set the number of true features in each covariate block to 10, and the magnitudes of the true features are drawn from $Unif(4,10)$ with random sign assignment. However, for the Gaussian block in the non-iid simulations below, we artificially lowered the SNR since we know that recovering continuous features is easier than recovering non-continuous features (see Figure 6). Unless stated otherwise, we simulate $n=200$ samples and $p_{1}=p_{2}=p_{3}=300$ features. We also center and scale the design matrix $\operatorname{\mathbf{X}}$ before estimation.

iid Design.   For each of the three covariate blocks, we simulate $n=200$ samples from iid features. Here, $p_{1}=p_{2}=p_{3}=300$ for the high-dimensional design and $p_{1}=p_{2}=p_{3}=100$ for the low-dimensional design.

Heteroscedasticity Design.   In this design, we assume that the features are independent but have non-constant variance. For the Gaussian block, the entries in each column are simulated from the normal distribution $N(0,\sigma^{2})$, where $\sigma\sim\Gamma(3,0.6)$. In the Bernoulli block, each column is simulated independently with entries drawn from $Bern(p)$ with $p\sim Unif(0.2,0.8)$. Similarly, in the Poisson block, the mean $\lambda$ of each column is drawn from the Gamma distribution $\Gamma(4,0.6)$ (using the shape/scale parameterization).

Block Directed Graph Design.   We next drop the independence assumption and use a Block Directed Markov Random Field (BDMRF) (Yang et al., 2014a) graph to simulate correlated features. In this case, $\operatorname{\mathbf{X}}$ is simulated via Gibbs sampling with the partial ordering of the underlying mixed graph given by $P[X_{1},X_{2},X_{3}]=P[X_{1}|X_{2},X_{3}]P[X_{2}|X_{3}]P[X_{3}]$, where $P[X_{1}|X_{2},X_{3}]$ is a pairwise Gaussian conditional random field (CRF), $P[X_{2}|X_{3}]$ is a pairwise Ising CRF (Ravikumar et al., 2010), and $P[X_{3}]$ is a pairwise Poisson Markov Random Field (MRF) (Yang et al., 2013, 2012). We set high correlations for the Gaussian and Poisson blocks and low correlations for the binary block and between block structure.

Ovarian Cancer Inspired Simulation Design.   In an attempt to simulate data closest to real world scenarios, we take the continuous-valued miRNA data, proportion-valued methylation data, and the count-valued RNASeq data from The Cancer Genome Atlas (TCGA) ovarian cancer database (The Cancer Genome Atlas Research Network, 2011) to be our covariates. After merging and preprocessing the TCGA ovarian cancer data (refer to Section 5 for details), we arrive at $n=293$ samples and $p_{\text{RNASeq}}=408$, $p_{\text{miRNA}}=307$, and $p_{\text{Methyl}}=301$ features.

Under each of these simulation scenarios, we evaluate the performance of B-RAIL and the oracle Lasso-type methods by reporting the true positive rate (TPR) and false discovery proportion (FDP) for overall feature recovery and individual block recoveries. Due to the large number of features, we use FDP, defined as the number of false positives divided by total the number of recovered non-zero features, instead of the false discovery rate.

We summarize the results of our simulations with Gaussian responses in Table 1 and those with binary and Poisson responses in Table 2. Note that for the binary and Poisson responses, we show the block directed graph results here and provide the other simulation results in the Appendix. We also highlight in bold the TPR/FDP combination with the highest TPR*(1-FDP) value for overall recovery. In almost all scenarios, the results in Table 1 and Table 2 indicate that B-RAIL (with no oracle information) is able to achieve a higher TPR and lower FDP than its competitive Lasso-type methods with oracle information.

When oracle information is unavailable, model selection techniques can introduce additional errors and further complicate feature selection. Table 3 shows one such case and compares the block directed graph simulation performance of B-RAIL against the Lasso-type methods using 5-fold cross-validation, extended BIC, and stability selection to select the penalty parameters. We also include the oracle estimators for the same set of simulations to emphasize the large decrease in performance when the Lasso-type methods do not have oracle information. These simulations indicate that cross-validation tends to over-select the number of features in the model while extended BIC under-selects, and stability selection performs the best but pales in comparison to oracle selection. In contrast, B-RAIL when initialized to an over-selection using the pre-specified sparsity level of $0.2p_{k}$ outperforms the Lasso-type methods even when oracle selection is used for these competitive methods. Additional simulations, confirming the strong empirical performance of B-RAIL, are provided in the Appendix.

5 Case Study: Integrative Genomics of Ovarian Cancer

One promising practical application for our research on multi-view feature selection lies in integrative cancer genomics. Here, scientists seek to integrate data from multiple sources of high-throughput genomic data to more holistically model the genomic systems in cancer cells, leading to a better understanding of disease mechanisms and possible therapies.

In this case study, we seek to integrate three different types of genomic data to study how epigenetics and short RNAs influence the gene regulatory system in ovarian cancer. Specifically, we are interested in discovering miRNAs and CpG sites, which affect the gene expression of well-known oncogenes in ovarian cancer and hence can serve as potential drug targets for blocking or decreasing the expression of these oncogenes. Driven by this goal of discovering potential drug targets, we use our proposed B-RAIL method to estimate the integrative ovarian cancer gene regulatory network with the specific intention of identifying miRNAs and CpG sites that are directly linked to known oncogenes of ovarian cancer.

In this investigation, we integrate the following three data sets: (1) count-valued gene expression measured via RNASeq, (2) continuous (Gaussian) miRNA expression, and (3) proportion-valued DNA methylation data from The Cancer Genome Atlas (TCGA) ovarian cancer study, which is publicly available (The Cancer Genome Atlas Research Network, 2011). The TCGA data originally contained $19,990$ genes, $27,578$ CpG sites, and $799$ miRNAs but only $n=293$ common patients across all three data sets of interest. We hence reduced the number of features to manageable sizes by first filtering features according to their association with several important clinical outcomes - survival via a univariate cox model, chemo-resistance via a univariate logistic model, and recurrence via a univariate logistic model. In addition, we transformed the RNASeq data using the Kolmogorov-Smirnov Test ($\alpha=0.262$) to alleviate the problem of very large counts (up to 20,000). This preprocessing yielded $p_{1}=408$ genes, $p_{2}=301$ CpG sites, and $p_{3}=307$ miRNAs in the RNASeq, methylation, and miRNA data sets, respectively. Lastly, per the recommendation of scientists, we included 20 additional highly mutated genes that were experimentally identified as important in ovarian cancer, resulting in $p_{1}=428$ genes in the RNASeq data set.

To estimate the integrated ovarian cancer network, we fit a Block Directed Markov Random Field (BDMRF) model (Yang et al., 2014a) using B-RAIL to estimate the neighborhood of each node in the graph. Note that since miRNAs and methylation are both gene regulatory mechanisms, miRNAs and methylation can affect expression levels (measured via RNASeq), but the converse is not possible. To agree with this known physical mechanism, we set the partial ordering of the mixed graph underlying BDMRF as $P[X_{1},X_{2},X_{3}]=P[X_{1}|X_{2},X_{3}]P[X_{2}]P[X_{3}]$, where $P[X_{2}]$ is a pairwise Ising MRF for the proportion-valued methylation data, $P[X_{3}]$ is a pairwise Gaussian MRF for the continuous miRNA data, and $P[X_{1}|X_{2},X_{3}]$ is a pairwise Poisson CRF for the count-valued RNASeq data. However, we recall that only negative conditional dependencies are permitted in the Poisson MRF and CRF models. Since this constraint is unrealistic for genomics data, we fit a Sub-Linear Poisson CRF, in lieu of the usual Poisson CRF, to allow for both positive and negative conditional dependencies Yang et al. (2013). Under this specified BDMRF model, we employ node-wise neighborhood selection (Meinshausen and Bühlmann, 2006; Yang et al., 2015) using B-RAIL to learn the edge structure of the integrated network.

Our overall BRAIL-estimated network is presented in Figure 7, and in Figure 10, we more closely examine the relationships between the oncogenes, miRNAs, and CpG sites by zooming in on the sub-networks for the well-known oncogene, BRCA1, and its direct neighbor, miRNA23b. Both BRCA1 and miRNA23b are well-known biomarkers and have been implicated in several ovarian cancer studies (Antoniou et al., 2003; King et al., 2003; BRCA, 1994; Li et al., 2014; Yan et al., 2016; Geng et al., 2012). Moreover, miRNA23b is known to play a key role in p53 signaling (via TP53) (Boren et al., 2009), agreeing with the estimated edge between the TP53 oncogene and miRNA23b in Figure 10(b).

Aside from this link however, the estimated edges between genes, CpG sites, and miRNAs in Figure 7 are largely unexplored and unknown by researchers since B-RAIL is one of the first practical approaches for multi-view feature selection. Nevertheless, we can partially validate our B-RAIL-estimated network by highlighting the many genes with verified connections in the ovarian cancer and cancer proliferation/suppression literatures. In Figure 10, we circle this collection of implicated genes, which includes LDOC1, SGCB, and miRNA210 (Buchholtz et al., 2014; Obermayr et al., 2010; Giannakakis et al., 2008).

As we have noted, there is substantial evidence in the scientific literature, suggesting that our proposed B-RAIL algorithm successfully identified promising candidates as well as known biomarkers involved in ovarian cancer. By taking into account the relationships between genes, miRNAs, and CpG sites, our integrative analysis via B-RAIL leads to valuable insights beyond a single biomarker type and to novel discoveries of direct connections between miRNAs, CpG sites, and known oncogenes, which may aid the development of targeted drug therapies for ovarian cancer. This is the first integrative analysis of its kind, and future experiments studying the connections between known ovarian cancer oncogenes and candidate miRNAs and CpG sites would be of great value to validate our findings.

6 Discussion

Though we have primarily focused on applications to integrative genomics in this work, B-RAIL is not limited to this context. B-RAIL can be applied to any field that yields high-dimensional multi-view data, and with the rapid advances in technologies, we expect B-RAIL to have a growing and far-reaching impact in fields such as imaging genetics, national security, climate studies, spatial statistics, Internet data, marketing, and economics. B-RAIL is also a versatile tool that can be used to for any sparse regression or graph selection problem in this multi-view context.

In addition to developing an effective data integration tool for multi-view feature selection, our work addresses the many difficulties of performing multi-view feature selection in practice. These practical challenges were severely under-studied prior to this work, but we partially resolve this gap, identifying four root challenges which interact with one another to impede recovery. Throughout our investigation of these practical challenges, we provide strong empirical evidence of the existence as well as the adverse consequences of such challenges. However, the theoretical underpinnings of these issues are still unknown. Understanding exactly how challenges such as shrinkage noise and the beta-min condition are influenced by varying domains and signals would be of great benefit to the field of data integration as a whole. We also highlight that while the Lasso has been well-studied under Gaussianity and idealized assumptions, the increasing abundance of correlated non-Gaussian data in multi-view settings requires a greater push for theoretical studies on feature selection with heterogeneous data and the GLM Lasso.

Overall, we have demonstrated many challenges posed by multi-view feature selection, and in our investigation of these challenges, we opened up new avenues for future theoretical work to rigorously understand how the heterogeneity of multi-view data complicates feature selection. Driven by these challenges and the ineffectiveness of existing methods, we developed a practical solution to overcome the current challenges. Our method, B-RAIL, is one of the first practical tools for multi-view feature selection and is grounded in deep theoretical foundations. With its versatility and strong empirical performance, B-RAIL facilitates impactful integrative analyses across a broad spectrum of fields.

Acknowledgements

Y.B. acknowledges support from the NIH/NCI T32 CA096520 training program in Biostatistics for Cancer Research, grant number 5T32CA09652011. T.T. acknowledges support from the NSF Graduate Research Fellowship Program DGE-1752814 and ARO, grant number W911NF1710005. GA acknowledges support from NSF DMS-1554821, NSF NeuroNex-1707400, and NSF DMS-1264058. We also thank Zhandong Liu, Ying-Wooi Wan, and Matthew L. Anderson at Baylor College of Medicine for thoughtful discussions related to this work.

Appendix A B-RAIL Convergence

We investigate the convergence of the B-RAIL algorithm. Our empirical analysis indicates that B-RAIL has quick support convergence for all simulation scenarios. We demonstrate this convergence for one type of simulation in Figure 14. Here, we simulate data using predictors from the TCGA ovarian cancer data (see Section 5) for three types of responses: (a) continuous, (b) binary, and (c) counts. We report the number of iterations over 100 runs of the B-RAIL algorithm and denote the average number of iterations until convergence by the dashed vertical black line. We see that the average number of iterations is between 4 and 5 with the maximum number of iterations reaching 15. These ranges were similar for all simulation designs, confirming relatively fast convergence of B-RAIL. Furthermore, we also show the true positive rates and the total number of selected features in Figure 14 to highlight B-RAIL’s convergence to a relatively accurate solution.

We provide the blockwise Lasso algorithm in Algorithm 3 and prove its convergence below.

Proposition 2**.**

Suppose that the objective function in (3.6) is bounded below. Then the blockwise (GLM) Lasso converges to a global minimum of(3.6).

Proof.

Let $f(z_{1},\dots,z_{k})$ denote the objective function in (3.6), where $z_{k}=(\alpha_{k},\bm{\beta}_{k})$ for each $k=1,\dots,K$. Since the domain of the GLM log-likelihood $\ell$ is open and Gateaux-differentiable on its domain, then $f$ is regular at each point in the domain of $f$ by Lemma 3.1 in Tseng (2001). Note also that since the GLM negative log-likelihood is convex and the $\ell_{1}$ penalty term in (3.6) is convex and separable, then $f(z_{1},\dots,z_{K})$ is convex with respect to each block $z_{k}$ for $k=1,\dots,K$. Because $f$ is regular and convex with respect to each block $z_{k}$ ($k=1,\dots,K$), it follows that the blockwise (GLM) Lasso always converges to a stationary point of $f$ by Theorem 4.1 in Tseng (2001). By convexity of $f$, this implies that the blockwise (GLM) Lasso converges to a global minimum. ∎

Appendix B Additional B-RAIL Simulations

We provide the following figures and tables to supplement our simulations and to further support the strong empirical performance of B-RAIL.

To augment the motivating example in Figure 3, we provide comparisons of B-RAIL to two other mixed graphical selection methods from the mgm R package (Haslbeck and Waldorp, 2015) in Figure 17. mgm takes a node-wise neighborhood estimation approach, and for each node, mgm selects the Lasso regularization parameter using either the extended BIC or CV, fits a penalized GLM model, and applies additional thresholding to the estimated coefficients to remove noise. Here, we used the “AND” rule to combine estimated neighborhoods for all three graphs. (Note that we had to convert the proportion-valued methylation values into 0-1 binary values in order to comply with mgm package restrictions.)

From Figure 17, we see that mgm with the extended BIC selection criteria tends to under-select features while mgm with CV often over-selects. This agrees with our simulations and discussion of the Lasso-type model selection biases in Section 4. We also observe that like the Lasso-type methods in Figure 3, mgm with CV and extended BIC can result in imbalanced selection between the blocks.

Figure 23 duplicates the information in Table 1 but using boxplots for better visualization. We recall that these simulations compared B-RAIL and various Lasso-type methods (using oracle information) under four simulation designs with Gaussian responses (see Section 4 for further details). In almost all of these simulations, B-RAIL is able to achieve a higher TPR while maintaining a low FDP.

We next verify that the above simulation results are not heavily dependent on our choice of $p_{k}$ and $||\beta_{k}||_{0}$. In Figure 26, we ran the iid simulation design with Gaussian responses for different non-uniform values of $p_{k}$ and $||\beta_{k}||_{0}$. These results show that B-RAIL can successfully recover features from unequally sized blocks with different amounts of sparsity while other methods may struggle to account for biases introduced by the different $p_{k}$’s and $||\beta_{k}||_{0}$’s.

For Poisson and binary responses, Tables 4 and 5 provide the results from additional simulation designs to supplement Table 2. Here, the response $\operatorname{\mathbf{y}}$ and predictors $\operatorname{\mathbf{X}}$ were simulated according to the description in Section 4 with $n=200$, $p_{1}=p_{2}=p_{3}=300$, and $||\bm{\beta}_{k}||_{0}=10$ for each block.

For easier visualizations, we duplicate the results of Tables 4 and 5 using boxplots in Figure 33. As shown by the plots, B-RAIL is able to achieve higher TPR and maintain low FDP across various simulations for both binary and Poisson responses.

In Table 6, we compare B-RAIL to the non-convex penalties, MCP and SCAD, under the block directed graph simulation design with three different types of responses. We see that MCP performs well with Gaussian responses, and in fact, the MCP penalty can be used instead of a Lasso penalty in the B-RAIL algorithm for Gaussian responses. However, the B-RAIL algorithm outperforms MCP and SCAD for non-Gaussian responses. We thus chose to use a Lasso penalty when introducing the B-RAIL algorithm for consistency.

Figure 38 provides the same information as Table 3 but using boxplots for easier visualization. While the model selection techniques (i.e. CV, extended BIC, and stability selection) give lower values of $TPR*(1-FDP)$ than their oracle selection counterparts, B-RAIL outperforms even the oracle selection methods and yields the highest $TPR*(1-FDP)$.