Network classification with applications to brain connectomics

TL;DR

This paper introduces a novel graph classification method that combines edge weights and network structure to analyze brain connectivity data, aiding in distinguishing brain disorders.

Contribution

The paper develops a new computationally efficient classification approach that promotes sparsity in nodes and edges, improving interpretability of brain network differences.

Findings

Effective in classifying schizophrenia from fMRI data

Produces interpretable models highlighting key brain regions

Outperforms existing methods in accuracy and sparsity

Abstract

While statistical analysis of a single network has received a lot of attention in recent years, with a focus on social networks, analysis of a sample of networks presents its own challenges which require a different set of analytic tools. Here we study the problem of classification of networks with labeled nodes, motivated by applications in neuroimaging. Brain networks are constructed from imaging data to represent functional connectivity between regions of the brain, and previous work has shown the potential of such networks to distinguish between various brain disorders, giving rise to a network classification problem. Existing approaches tend to either treat all edge weights as a long vector, ignoring the network structure, or focus on graph topology as represented by summary measures while ignoring the edge weights. Our goal is to design a classification method that uses both the individual edge information and the network structure of the data in a computationally efficient way, and that can produce a parsimonious and interpretable representation of differences in brain connectivity patterns between classes. We propose a graph classification method that uses edge weights as predictors but incorporates the network nature of the data via penalties that promote sparsity in the number of nodes, in addition to the usual sparsity penalties that encourage selection of edges. We implement the method via efficient convex optimization and provide a detailed analysis of data from two fMRI studies of schizophrenia.

Tables

Table 1. Table 1 : Cross-validated accuracy (average and standard errors over 10 folds) for different methods.

	Classification accuracy % (s.e.)
Method	COBRE	UMich
With variable selection
Our method (GC)	92.7 (2.6)	85.9 (3.6)
Elastic net	89.5 (1.8)	82.6 (4.7)
SVM-L1	87.9 (2.2)	86.2 (4.3)
Signal-subgraph	86.1 (3.3)	82.4 (3.3)
DLDA	84.6 (3.3)	73.4 (3.9)
Lasso	80.1 (5.6)	60.9 (5.6)
No variable selection
SVM	93.5 (2.1)	89.8 (2.5)
Ridge penalty	91 (2.6)	80.9 (3.5)
Random forest	74.2 (2.6)	82.1 (3.9)
Network summaries	61.4 (3.1)	65 (7.2)

Table 2. Table 2 : Edges with the top 15 largest selection probabilities from stability selection. The first column shows the pair of nodes making the edge, the second column the brain systems the nodes belong to in the Power parcellation, and the third column the fitted coefficient of the edge

	COBRE			UMich
	Edge	Systems	Coefficient	Edge	Systems	Coefficient
1	(208, 85)	(9, -1)	-0.187	(110, 207)	(5, 9)	-0.013
2	(260, 11)	(12, -1)	0.183	(255, 113)	(1, 5)	0.014
3	(194, 140)	(8, -1)	0.136	(33, 218)	(1, 9)	0.016
4	(52, 186)	(3, 8)	- 0.1	(46, 225)	(2, 10)	0.013
5	(160, 239)	(7, 11)	-0.082	(43, 90)	(2, 5)	-0.013
6	(120, 116)	(5, 5)	0.099	(23, 225)	(1, 10)	0.012
7	(57, 129)	(3, 5)	-0.128	(66, 118)	(4, 5)	-0.013
8	(24, 114)	(1, 5)	-0.148	(26, 145)	(1, 7)	0.013
9	(81, 179)	(5, 8)	-0.129	(186, 254)	(8, -1)	0.012
10	(193, 140)	(8, -1)	0.153	(15, 134)	(1, 6)	0.011
11	(178, 234)	(8, 10)	0.146	(76, 207)	(5, 9)	-0.012
12	(18, 194)	(1, 8)	0.116	(65, 84)	(4, -1)	-0.012
13	(215, 207)	(9, 9)	-0.076	(26, 122)	(1, 5)	0.012
14	(90, 224)	(5, 10)	0.123	(33, 145)	(1, 7)	0.012
15	(112, 253)	(5, -1)	0.136	(36, 224)	(1, 10)	0.011

Table 3. Table 3 : Classification accuracy (cross-validation average and standard error) of the classifier fitted on one dataset and evaluated on the other. The intercept (the mean) is fitted on the test data and the accuracy is estimated using 10-fold cross-validation on the test data.

	Test data
Training data	COBRE	UMich
COBRE	92.7 (2.6)	73.5 (3.4)
UMich	78.3 (3.0)	85.9 (3.6)

Table 4. Table 4 : Summary of the two datasets.

Dataset	# nodes	Status	#	Sex M/F	Age mean (s.d.)
COBRE	263	Schizophrenic	54	48/6	35.4 (13.1)
		Control	70	48/22	35.1 (11.5)
UMich	264	Schizophrenic	39	29/10	40.7 (11.5)
		Control	40	28/12	36.8 (12.3)

Equations

$$\hat{B}=\underset{B\in \mathcal{B}}{\mathrm{argmin}}\left\{\ell (B)+\Omega (B)\right\},$$

$$\ell (B)=\frac{1}{n}\sum _{k=1}^n\tilde{\ell }(Y_k,A^{(k)};B)$$

$$\tilde{\ell }\left(Y_k,A^{(k)};B,b\right)=\log \left\{1+\exp \left(-Y_k⟨A^{(k)},B⟩+b\right)\right\}.$$

$${\Omega }_{\lambda ,\rho }(B)=\lambda \left(\sum _{i=1}^N\Vert B_{(i)}{\Vert }_2+\rho \Vert B{\Vert }_1\right),$$

$$\tilde{\mathcal{B}}=\left\{B{\in }^{N\times N},\text{diag}(B)=0\right\}.$$

$${\Psi }_{\lambda ,\rho }(B)=\lambda \left(\sum _{i=1}^N\sqrt{\sum _j\left(B_{ij}^2+B_{ji}^2\right)}+\rho \Vert B{\Vert }_1\right),$$

$${\tilde{\Psi }}_{\lambda ,\rho }(B^{(1)},B^{(2)})=\lambda \left[\sum _{i=1}^N\sqrt{\sum _j\{(B_{ij}^{(1)}{)}^2+(B_{ji}^{(2)}{)}^2\}}+\rho (\Vert B^{(1)}{\Vert }_1+\Vert B^{(2)}{\Vert }_1)\right],$$

$$(\hat{B},\hat{b})=\underset{B\in \mathcal{B},b\in }{\mathrm{argmin}}$$

$$\lambda \left(\sum _{i=1}^N\Vert B_i{\Vert }_2+\rho \Vert B{\Vert }_1\right)\}$$

$$W^{(k)}$$

$$B^{(k)}$$

$$=\underset{B\in \mathcal{B}}{\mathrm{argmin}}\left\{\frac{1}{2}{\Vert B-\left(W^{(k)}-t^{(k)}\nabla \ell (W^{(k)})\right)\Vert }_2^2+t^{(k)}\Omega (B)\right\},$$

$$\ell (B^{(k)})\le \ell (W^{(k)})+⟨\nabla \ell (W^{(k)}),B^{(k)}-W^{(k)}⟩+\frac{1}{2t^{(k)}}\Vert B^{(k)}-W^{(k)}{\Vert }_2^2$$

$$\underset{\tilde{B},Q,R}{\min }$$

$$\tilde{B}=Q,\tilde{B}=R,\tilde{B}={\tilde{B}}^T,\text{diag}(\tilde{B})=0.$$

$${\mathcal{L}}_{\mu }=$$

$$⟨V,\tilde{B}-R⟩\frac{\mu }{2}\left(\Vert \tilde{B}-Q{\Vert }_2^2+\Vert \tilde{B}-R{\Vert }_2^2+\Vert \tilde{B}-{\tilde{B}}^T{\Vert }_2^2\right).$$

$$\underset{\tilde{B},Q,\tilde{Q},R}{\min }$$

$$\tilde{B}=Q,\tilde{B}=\tilde{Q},\tilde{B}=R,\tilde{B}={\tilde{B}}^T,\text{diag}(\tilde{B})=0.$$

$$Y^{(k)}$$

$$Z_{ij}^{(k)}$$

$$B_{(i)}^{(k)}$$

$$\mathcal{M}=\left\{B{\in }^{N\times N}|B_{ij}=0\text{ for all }i\in {\mathcal{G}}^C\text{ or }j\in {\mathcal{G}}^C,B=B^T\right\}$$

$${\nabla }_{(i,j),(k,l)}^2\ell (B)=\frac{{\partial }^2\ell (B)}{\partial B_{ij}\partial B_{kl}},$$

$${\left(H_{(i,j),\mathcal{G}}\right)}_{k,l}=\text{Tr}\left\{\left({\nabla }_{(i,j),(\mathcal{G},\mathcal{G})}^2\ell (B^{\star })\right){\Lambda }_{(k,l),(\cdot ,\cdot )}\right\},1\le k,l\le G,,$$

$$\sum _{i,j}{\Delta }_{i,j}\text{Tr}\left\{\left({\nabla }_{(i,j),(\cdot ,\cdot )}^2\ell (B)\right)\Delta \right\}$$

$$\Vert {\nabla }^2\ell (B)-{\nabla }^2\ell (B^{\star }){\Vert }_2$$

$$\underset{i\in {\mathcal{G}}^C}{\max }{\Vert {\left(\sum _{k=1}^G\Vert (H_{(i,j),\mathcal{G}}{)}_{k\cdot }{\Vert }_2\right)}_{j=1}^N\Vert }_2=1-\tau <0.$$

$$\underset{i,j}{\max }\mathbb{P}\left(\Vert {\nabla }_{ij}\tilde{\ell }\left(Y,A;B^{\star }\right){\Vert }_{\infty }>t\right)\le 2\exp \left(-t^2/{\sigma }^2\right).$$

$$\Vert \hat{B}-B^{\star }{\Vert }_2=O_P\left(\sigma N\sqrt{\frac{\log N}{n}}\right).$$

$$\rho >\frac{1}{\tau }-\frac{1}{\sqrt{G}},$$

$$\Vert \hat{B}-B^{\star }{\Vert }_2$$

$$\mathbb{P}\left(\hat{\mathcal{G}}\subseteq \mathcal{G}\right)$$

$$\mu_{kl}=\left\{\begin{array}[]{cl}0.3,&\text{if $k=l$},\\ 0.1&\text{if $k\neq l$},\end{array}\right.$$

$$\text{EFPR}(\mathcal{M},\eta )=\frac{\left|\hat{\mathcal{E}}(\mathcal{M},\eta )\cap {\mathcal{E}}^C\right|}{\left|{\mathcal{E}}^C\right|},\text{ETPR}(\mathcal{M},\eta )=\frac{\left|\hat{\mathcal{E}}(\mathcal{M},\eta )\cap \mathcal{E}\right|}{\left|\mathcal{E}\right|}.$$

Full text

graph classification, high-dimensional data, variable selection, fMRI data,

arXiv:1701.08140

\startlocaldefs

\endlocaldefs

,

and

1 Introduction

Network data analysis has received a lot of attention in recent literature, especially unsupervised analysis of a single network which is thought of as generated from an exchangeable random graph model, for example Bickel and Chen, (2009); Le et al., (2017); Zhang and Zhou, (2016); Gao et al., (2017) and many others. This setting is applicable to a number of real life scenarios, such as social networks, but there are situations where the network nodes are labeled and therefore not exchangeable, and/or more than one network is available for analysis, which have received relatively less attention. Here we focus on the setting motivated by brain connectomics studies, where a sample of networks is available from multiple populations of interest (for example, mentally ill patients and healthy controls). In this setting, each unit in the population (e.g., a patient) is represented by their own network, and the nodes (brain regions of interest) are labeled and shared across all networks through a registration process that maps all individual brains onto a common atlas. There are many classical statistical inference questions one can ask in this setting, for example, how to compare different populations (Tang et al., 2017b, ; Tang et al., 2017a, ). The question we focus on in this paper is a classification problem: given a training sample of networks with labeled nodes drawn from multiple classes, the goal is to learn the rules for predicting the class of a given network, and just as importantly, interpret these rules.

Network methods are a popular tool in the neuroscience literature (Bullmore and Sporns,, 2009; Bullmore and Bassett,, 2011). A brain network represents connectivity between different locations of an individual’s brain. How connectivity is defined varies with the type of imaging technology used and the conditions under which data were collected. In this paper, we focus on functional connectivity, which is a measure of statistical association between each pair of locations in the brain, constructed from functional magnetic resonance imaging (fMRI) data, although the methods we develop are applicable to any sample of weighted networks with labeled nodes. In fMRI studies, BOLD (blood oxygen-level dependent) signal, a known correlate of underlying neural activity, is measured at a sequence of time points at many spatial locations in the brain, known as voxels, resulting in a 4-dimensional data array, with three spatial dimensions and a time index. Brain networks constructed from fMRI data have been successfully used for various tasks, such as differentiating between certain illnesses, or between types of external stimuli (Bullmore and Sporns,, 2009), and contain enough information to identify individual subjects (Finn et al.,, 2015). Extensive statistical literature has focused on the analysis of raw fMRI data (Lindquist,, 2008; Zhou et al.,, 2013; Zhang et al.,, 2016), usually aiming to characterize brain activation patterns obtained from task-based fMRI experiments. In this paper, we focus on resting-state fMRI data, where no particular task is performed and subjects are free to think about anything they want. Thus registering the time dimension across different subjects is not possible. The connectivity network approach, which averages over the time dimension in computing a measure of dependence between different voxels, is thus a natural choice, and has been widely used with multiple types of neuroimaging data.

Two different datasets are analyzed in this paper, both of resting state fMRI studies on schizophrenic patients and healthy controls. One dataset, COBRE (54 schizophrenics and 70 controls), is publicly available (Aine et al.,, 2017); another, which we will refer to as UMich data (39 schizophrenics and 40 controls), was collected internally in the last author’s lab. Having two datasets on the same disease allows us to cross-check models trained on one of them for classification on the other to check the robustness of our approach. The raw data arrays undergo pre-processing and registration steps, discussed in detail in the Appendix C, along with additional details on data collection. To construct a brain network from fMRI measurements, a set of nodes is chosen, typically corresponding to regions of interests (ROIs) from some predefined parcellation. In our analysis we use the parcellation of Power et al., (2011) (see Figure 1), which consists of 264 ROIs divided into 14 functional brain systems (in the COBRE data, node 75 is missing). A connectivity measure is then computed for every pair of nodes, resulting in an adjacency matrix of size $264\times 264$. Many choices of connectivity measures are available (Smith et al.,, 2013); perhaps the most commonly used one is the Pearson correlation coefficient between locations, computed by averaging over the time dimension. It has been argued that partial correlations are a better measure of connectivity (Varoquaux and Craddock,, 2013; Narayan et al.,, 2015), but the choice depends on the final goal of analysis. In this paper we follow the vast majority of the connectomics literature and measure connectivity on each individual by using marginal correlations between the corresponding time series (see Figure 2). The correlations are then further rank-transformed and standardized; see Appendix C for details. These transformations are intended to deal with subject-to-subject variability and the global signal regression issue (Gotts et al.,, 2013), and although they lose some information, we observed that on our datasets classification based on standardized ranks of marginal correlations outperformed classification based on other connectivity measures, such as marginal correlations. The methods we develop here are applicable to networks that encode any type of connectivity measure.

The problem of graph classification has been studied previously in other contexts, with a substantial literature motivated by the problem of classification of chemical compounds (Srinivasan et al.,, 1996; Helma et al.,, 2001), where graphs represent the compound’s molecular structure. This setting is very different, with small networks of about 20 nodes on average, binary or categorical edges recorded with no noise, and different nodes corresponding to different networks, (Ketkar et al.,, 2009). Classification methods for chemical compounds are usually based on finding certain discriminative patterns in the graphs, like subgraphs or paths (Inokuchi et al.,, 2000; Gonzalez et al.,, 2000), and using them as features for training a standard classification method (Deshpande et al.,, 2005; Kudo et al.,, 2004; Fei and Huan,, 2010). Computationally, finding these patterns is only possible on small binary networks.

Another type of methods are based on graph kernels (Gärtner et al.,, 2003; Vishwanathan et al.,, 2010), which define a similarity measure between two networks. These kernels combined with support vector machines (SVMs) have been successfully used on small networks (Kashima et al.,, 2003; Borgwardt et al.,, 2005), but the curse of dimensionality makes local kernel methods unsuitable for large scale networks (Bengio and Monperrus,, 2005). On our datasets, graph kernel methods did not perform better than random guessing.

In the context of classifying large-scale brain networks, two main approaches have been followed. One approach is to reduce the network to its global summary measures such as the average degree, clustering coefficient, or average path length (Bullmore and Sporns,, 2009), and use those measures as features for training a classification method. Previous studies have reported significant differences on some of these network measures for groups of patients with certain brain diseases compared with healthy controls (Supekar et al.,, 2008; Liu et al.,, 2008), suggesting their usefulness as diagnostic biomarkers. However, global summary statistics collapse all local network information, which can harm the accuracy of classification and do not allow to identify local differences. In our data analysis, a method based on the network measures suggested in Prasad et al., (2015) performed poorly for classification (see Section 6).

An alternative approach to classification of large networks is to treat edge weights as a “bag of features”, vectorizing the unique elements of the adjacency matrix and ignoring the network nature of the data. This approach can leverage many existing classification methods for vectors, and provides an interpretation at the edge level if variable selection is applied (Richiardi et al.,, 2011; Craddock et al.,, 2009; Zhang et al.,, 2012). Spatial correlation between edges connecting neighboring nodes can be incorporated (Watanabe et al.,, 2014; Scott et al.,, 2015), although the effectiveness of this regularization will depend on the parcellation used to define nodes. Alternatively, an individual test can be used for each edge to find significant differences between two populations, with a multiple testing correction and without constructing a classifier at all (Narayan et al.,, 2015). While these methods can deliver good predictions, their interpretability is limited to individual edge selection, which is less scientifically interesting than identifying differentiating nodes or regions, and they cannot account for network structure.

Taking the network structure into account can have benefits for both testing and classification settings. Some methods perform inference over groups of edges based on the community assignments of the nodes to which they are incident. For example, Sripada et al., 2014a ; Sripada et al., 2014b introduced Network Contingency Analysis which begins with massive univariate testing at each edge, and then counts the number of superthreshold connections in each cell, a group of edges that connect nodes in two functional systems. Nonparametric methods are then used to conduct inference on the count statistic for each cell, with multiple comparison correction for inference at the cell level. Power can be improved by applying a network-based multiple testing dependence correction (Zalesky et al.,, 2010). For classification, better interpretability and potentially accuracy can be obtained if we focus on understanding which brain regions or interactions between them are responsible for the differences. In somewhat related work, Vogelstein et al., (2013) proposed to look for a minimal set of nodes which best explains the difference, though that requires solving a combinatorial problem. Hypothesis testing on a type of graph average has also been proposed (Ginestet et al.,, 2017). Bayesian nonparametrics approaches for modeling populations of networks allow to test for local edge differences between the groups (Durante and Dunson,, 2018), but are computationally feasible only for small networks.

Our goal in this paper is to develop a high-dimensional network classifier that uses all the individual edge weights but also respects the network structure of the data and produces more interpretable results. To achieve this goal, we use structured sparsity penalties to incorporate the network information by penalizing both the number of edges and the number of nodes selected. Although our main application here is classification of brain connectivity networks, our methods are applicable to any weighted graphs with labeled nodes, and to general prediction problems, not just classification.

The rest of this paper is organized as follows. In Section 2, we introduce our classifier and the structured penalties. In Section 3 we show how to efficiently solve the resulting convex optimization problem by a proximal algorithm, each step of which is a further optimization problem which we solve by the alternating direction method of multipliers (ADMM). The performance of our method is evaluated and compared with other methods using simulations in Section 5. In Section 6, we analyze two brain connectivity datasets, each containing schizophrenic patients and healthy controls, and show that our regularization framework leads to state-of-the-art accuracy while providing interpretable results, some of which are consistent with previous findings and some are new. We conclude with a brief discussion in Section 7.

2 A framework for node selection in graph classification

2.1 A penalized graph classification approach

We start from setting up notation. All graphs we consider are defined on the same set of $N$ labeled nodes. A graph can be represented with its adjacency matrix $A\in^{N\times N}$. We focus on graphs that are undirected ($A_{ij}=A_{ji}$) and contain no self-loops ($A_{ii}=0$). These assumptions are not required for the derivations below, but they match the neuroimaging setting and simplify notation. Our goal is predicting a class label $Y$ from the graph adjacency matrix $A$; in this paper we focus on the binary classification problem where $Y$ takes values $\{-1,1\}$, although extensions from binary to multi-class classification or real-valued responses are straightforward. Throughout this paper, we use $\|\cdot\|_{p}$ to denote the entry-wise $\ell_{p}$ norm, i.e., for a matrix $A\in^{m\times n}$, $\|A\|_{p}=\left(\sum_{i=1}^{m}\sum_{j=1}^{n}|A_{ij}|^{p}\right)^{1/p}$.

A standard general approach is to construct a linear classifier, which predicts the response $Y$ from a linear combination of the elements of $A$, $\langle A,B\rangle=\text{Tr}\left(B^{T}A\right)$, where we arrange the coefficients in a matrix $B\in^{N\times N}$ to emphasize the network nature of the predictors. We focus on linear classifiers here because variable selection is at least as important as prediction itself in the neuroimaging application, and setting coefficients to 0 is a natural way to achieve this. The coefficients are typically estimated from training data by minimizing an objective consisting of a loss function plus a penalty. The penalty can be used to regularize the problem to make the estimator well-defined in high-dimensional problems, to select important predictors, and to impose structure, and many such penalties have been proposed, starting from the lasso (Tibshirani,, 1996). Our focus is on designing a classifier in this framework that respects and utilizes the network nature of the predictors. In brain networks in particular, neuroscientists believe that edges are organized in subnetworks, also called brain systems (Power et al.,, 2011), that carry out specific functions, and certain subnetworks are important for prediction (Bullmore and Sporns,, 2009), although different studies tend to implicate different regions (Fornito et al.,, 2012). Thus we aim to find nodes or subnetworks with good discriminative power, and hence select both the most informative nodes and edges.

Although the methods we develop here can be used on small networks, our main focus is on the more challenging case of medium to large brain networks. In brain connectivity studies dealing with multiple subjects, while raw images may contain hundreds of thousands of voxels, they are commonly down-sampled according to a parcellation scheme with a coarser resolution, usually resulting in networks with hundreds or thousands of nodes representing ROIs (see for example Kiar et al., (2018)). This coarser resolution is essential for registration, as aligning different brains at a high resolution is much harder, but it still results in hundreds of thousands or millions of edges which serve as predictor variables. Given the typical data sizes in this area of application, we focus on methods based on convex formulations, which allow for efficient and scalable implementations with convergence guarantees.

Let $\{(A^{(1)},Y_{1}),\ldots,(A^{(n)},Y_{n})\}$ be the training sample of undirected adjacency matrices with their class labels, and let $Y=(Y_{1},\ldots,Y_{n})$. A generic linear classifier described above is computed by finding the coefficients $B$ defined by

[TABLE]

where $\mathcal{B}=\left\{B\in^{N\times N}:B=B^{T},\text{diag}(B)=0\right\}$, $\Omega$ is a penalty, and

[TABLE]

is a loss function evaluated on the training data. Our methodology can accommodate different choices of loss functions that can extend beyond classification problems (e.g., least squares or generalized linear models). The optimization algorithm presented in Section 3 can work with any convex and continuously differentiable loss function, and further assumptions are required for consistency (see Section 4). In this paper, for the purpose of classification we use the logistic loss function in the simulations and data analysis, which is defined as

[TABLE]

The threshold $b$ is an additional parameter to be estimated.

To capture structural assumptions on important predictive edges, we focus on convex structured sparsity penalties (Bach et al.,, 2012) that encourage a small number of active nodes, by which we mean nodes attached to at least one edge with a non-zero coefficient. One approach to finding a set of such nodes was proposed by Vogelstein et al., (2013), who called it a signal-subgraph, and proposed finding the minimal set of nodes (called signal vertices) which together are incident to all selected edges (but not every node connected to a selected edge is a signal vertex). Finding this set is a combinatorial optimization problem, and the set is not always uniquely defined. Instead, we focus on convex formulations that allow for efficient computation and encourage small active node sets indirectly.

Other convex penalties have been used for fMRI data as a way to enforce spatial smoothness in the solution (Grosenick et al.,, 2013; Xin et al.,, 2014; Hu and Allen,, 2015). These methods assume that voxels are equally spaced in the brain, and neighboring voxels are highly correlated. In particular, Watanabe et al., (2014) proposed penalties for brain network classification using these spatial assumptions. Here, instead of enforcing a spatial regularization directly, we aim for a regularization that can be applied to any type of network data, and in particular to brain networks with coarse and/or uneven parcellations where enforcing spatial smoothness may not work as well. In any case, the flexibility of convex optimization algorithms allows one to easily incorporate additional spatially-informed penalties if needed.

2.2 Selecting nodes and edges through group lasso

To reflect the network structure of the predictors, we use a penalty that promotes a sparse classifier not only in the number of edges used, but also in the number of nodes. The group lasso penalty (Yuan and Lin,, 2006) is designed to eliminate a group of variables simultaneously. Here we penalize the number of active nodes by treating all edges connected to one node as a group. Then eliminating this group (a row of coefficients in the matrix $B$) is equivalent to de-activating a node. The group penalty is defined as

[TABLE]

where $B_{(i)}$ denotes the vector of edge weights incident to the $i$-th node (or equivalently, the $i$-th row or column of $B$), and $\lambda,\rho\geq 0$ are tuning parameters. Note that the constraint $B=B^{T}$ makes the groups overlap, since a coefficient $B_{ij}$ belongs to groups associated with the nodes $i$ and $j$, and therefore, the edge between nodes $i$ and $j$ would be selected only if both nodes are activated. The second term in the penalty $\rho\|B\|_{1}$ acts as the usual lasso penalty to promote sparsity inside the group (Friedman et al., 2010a, ), allowing to select a subset of edges for an active node. Due to the overlap in the groups, this lasso penalty is usually necessary in order to produce sparse solution (see Proposition 1). The constraint $\text{diag}(B)=0$ in (2.1) is automatically enforced with this formulation.

Remark 1*.*

An alternative to the constraint in the problem (2.1) is to optimize over the set

[TABLE]

Without the symmetry constraint and assuming undirected graphs, the penalty (2.2) is equivalent to the overlapping group lasso formulation of Jacob et al., (2009). This formulation has some advantages. Since it gives group lasso without overlaps, the lasso penalty $\rho\|B\|_{1}$ is not required to obtain sparse solutions, and more efficient optimization algorithms exist for this case. This approach would loosely correspond to the idea of selecting signal nodes as in Vogelstein et al., (2013), in the sense that an edge can be selected if at least one of its nodes is selected, and the second node could be inactive. The downside is that each edge now corresponds to two different coefficients $B_{ij}$ and $B_{ji}$, the problem encountered by all variable selection methods that ignore symmetry, such as Meinshausen and Bühlmann, (2006). The standard solution for this problem, as suggested by Jacob et al., (2009), is to take the average of the coefficients. Intuitively, one would expect that the formulation using $\mathcal{B}$ would be better when the significant edges are incident to a small set of nodes, since both nodes have to be active for an edge to be selected, while using $\tilde{\mathcal{B}}$ may be better when for some nodes most of their edges are significant, creating “significant hubs”. Since in our application we are primarily looking for discriminative brain subnetworks, we focus on the symmetrically constrained formulation for the rest of the paper. We also found that in practice this second formulation results in less accurate classifiers for the neuroimaging data discussed in Section 6.

Remark 2*.*

The analogue to (2.2) for directed graphs would assign coefficients $B_{ij}$ and $B_{ij}$ to the same group, resulting in the penalty

[TABLE]

where $B\in\tilde{\mathcal{B}}$. Alternatively, we can also use the formulation of Remark 1, by replicating the variables and estimating two matrices of coefficients, say $B^{(1)}$ and $B^{(2)}$, with the penalty

[TABLE]

with $B^{(1)},B^{(2)}\in\tilde{\mathcal{B}}$, and set the coefficients matrix to $B=\left(B^{(1)}+B^{(2)}\right)/2$. This formulation will again not directly select subnetworks as discussed in Remark 1.

Finally, for numerical stability we add an extra ridge penalty term $\frac{\gamma}{2}\|B\|_{F}^{2}=\frac{\gamma}{2}\text{Tr}\left(B^{T}B\right)$, with $\gamma$ a small and fixed constant. There are several benefits of combining ridge and lasso penalties (see for example Zou and Hastie, (2005)). The parameter $\gamma$ can be potentially considered as an additional tuning parameter, but here we only use a small fixed constant $\gamma$ in order to avoid numerically degenerate solutions. In practice, the results are not sensitive to the exact value of $\gamma$.

Putting everything together, to fit our graph classifier, we solve the problem

[TABLE]

for given values of $\lambda$ and $\rho$, which will be chosen by cross-validation.

3 The optimization algorithm

Our optimization algorithm to solve the problem (2.4) combines two common approaches to convex optimization: proximal algorithms and alternating direction method of multipliers (ADMM). We use an accelerated version of the proximal algorithm (Beck and Teboulle,, 2009) to solve the main problem (2.4). In each step, we need to calculate a proximal operator, which is a further convex optimization problem solved with the ADMM algorithm.

The main optimization difficulty comes from the overlapping groups. Some algorithms have been proposed for this case, including a subgradient descent method (Duchi and Singer,, 2009), which has a slow rate of convergence, or proximal algorithms based on smoothing the original problem (Yuan et al.,, 2011; Chen et al.,, 2012). Although smoothing yields fast algorithms, it is not clear that the sparsity pattern is preserved with those approximations. We follow an approach similar to Yuan et al., (2011) and Chen et al., (2012), but solve the proximal operator for the penalty (2.2) directly using the ADMM method. This can potentially give a more accurate sparsity pattern, and the flexibility of the algorithm allows for additional penalties if desired, such as spatial smoothing similar to Watanabe et al., (2014) (see Remark 3).

The main problem (2.1) is solved with a proximal algorithm (see Parikh and Boyd, (2013)). Recall that the proximal operator for a function $f$ is defined as $\operatorname*{prox}_{f}(v)=\operatorname*{arg\,min}_{x}\{f(x)+\frac{1}{2}\|x-v\|_{2}^{2}\}$. Starting with an initial value $B^{(0)}\in^{N\times N}$, a proximal algorithm solves the optimization problem (2.1) by iteratively calculating the proximal operator of $\Omega=\Omega_{\lambda,\rho}$ for a descent direction of the differentiable loss function $\ell$. We use an accelerated version of the algorithm (Beck and Teboulle,, 2009), which for each $k=2,\ldots,$ until convergence, performs the updates

[TABLE]

where $\nabla\ell(W)\in^{N\times N}$ is the gradient of the loss function $\ell$ at $W$ and $\{t^{(k)}\}$ is a sequence of positive values. If $\nabla\ell$ is Lipschitz continuous, with $L$ its Lipschitz constant, the sequence of values $\ell(B^{(k)})+\Omega(B^{(k)})$ converges to the optimal value at rate $O(1/k^{2})$ if $t^{(k)}\in[0,1/L)$. The value of $t^{(k)}$ can be chosen using a backtracking search (Beck and Teboulle,, 2009), which decreases this value until the condition

[TABLE]

is satisfied. This procedure ensures that step sizes $\{t^{(k)}\}$ become smaller as the algorithm progresses, until $t^{(k)}<1/L$. In practice, $L$ might be large, which can make the algorithm slow to converge. It has been observed in other sparse high-dimensional problems that search strategies for $t^{(k)}$ which allow for $t^{(k)}>1/L$ when appropriate can actually speed up convergence (Scheinberg et al.,, 2014; Hastie et al.,, 2015). We use a strategy of this type, allowing $t^{(k)}$ to increase by a factor of $\alpha\geq 1$ if the relative improvement in the loss function on iteration $k$ becomes small. We observed that this strategy significantly reduces the number of iterations until convergence. The entire procedure is summarized in Algorithm 1 on the Appendix A.

The logistic loss function of (2.4) has an extra parameter $b$. Rather than including it as an unpenalized coefficient for a constant covariate, we use block coordinate descent and solve for $b$ separately. This is convenient because the threshold $b$ and the matrix of coefficients $B$ may not be on the same scale. Thus, $b$ can be updated by solving $b^{(k+1)}=\operatorname*{arg\,min}_{b\in}\ell(B^{(k)},b)$, which is easy to compute via Newton’s method.

The proximal algorithm requires solving the proximal operator (3.2), which has no closed form solution for the penalty (2.2) under the symmetry constraint. Strategies based on smoothing this penalty have been proposed (Yuan et al.,, 2011; Chen et al.,, 2012). However, to allow for variable selection which results from non-differentiability of the penalty, we aim to solve the proximal operator directly using ADMM (see Boyd et al., (2011) for a review). Note that if the symmetric constraint is relaxed as in Remark 1, the proximal operator has a closed form solution (see Remark 4).

The ADMM works by introducing additional constraints and performing coordinate descent in the corresponding augmented Lagrangian function. Setting $Z=W^{(k)}-t^{(k)}\nabla\ell(W^{(k)})$ and $t=t^{(k)}$, and introducing the variables $Q,R\in^{N\times N}$, we can formulate (3.2) as a convex optimization problem

[TABLE]

The ADMM algorithm introduces the multipliers $U,V\in^{N\times N}$ and a penalty parameter $\mu>0$ to perform gradient descent on the Lagrangian of (3.4), given by $\mathcal{L}_{\mu}=\mathcal{L}_{\mu}(\tilde{B},Q,R,U,V)$ as

[TABLE]

The value $\mu$ controls the gap between dual and primal feasibility. In practice, we observed that setting $\mu=0.1$ gives a good balance between primal and dual feasibility, although other self-tuning methods are available (Parikh and Boyd,, 2013). This function is optimized by coordinate descent, with each variable updated to minimize the value of $\mathcal{L}_{\mu}$ while all the other variables are fixed. This update has a closed form; the detailed steps of the ADMM are shown in Algorithm 2 on the Appendix A. These steps are performed until the algorithm converges within tolerance $\epsilon^{\text{ADMM}}>0$. Note that ADMM will be performed in each iteration of the algorithm to solve (2.4) and thus tolerance $\epsilon^{\text{ADMM}}$ can be decreased as the algorithm progresses. On the other hand, performing only one iteration of algorithm (2) gives a similar algorithm to the one of Chen et al., (2012).

Remark 3*.*

The ADMM makes it very easy to incorporate additional penalties. If $\Psi$ is a new penalty, we can rewrite (3.4) by introducing an additional parameter $\tilde{Q}$ so it becomes

[TABLE]

We can obtain the Lagrangian formulation (3.5) in a similar manner, and include new parameters in the ADMM updates, which can be performed efficiently as long as the proximal operator of $\Psi$ has a closed form solution. This is in fact the case for some other penalties of interest, such as the GraphNet penalty (Grosenick et al.,, 2013; Watanabe et al.,, 2014), which can incorporate spatial location information.

Remark 4*.*

The alternative formulation for the graph penalty given in Remark 1 corresponds to standard sparse group lasso (Friedman et al., 2010a, ). In particular, we can still employ the proximal algorithms (3.1) and (3.2), but instead optimize over the set $\tilde{\mathcal{B}}$. Without the symmetric constraint on $B$, the overlap in the group lasso penalty disappears, and this vastly simplifies the problem. Using Theorem 1 of Yuan et al., (2011), the update for $B^{(k)}$ has a closed form solution given by

[TABLE]

4 Theory

In this section, we show that the solution of the penalized problem (2.2) can recover the correct subgraph corresponding to the set of non-zero coefficients, and give its rates of convergence. The theory is a consequence of the results of Lee et al., (2015) for establishing model selection consistency of regularized M-estimators under geometric decomposability (see Appendix for details). We present explicit conditions for our penalty to work well, which depend on the data as well as the tuning parameters.

Let $B^{\star}\subset^{N\times N}$ be the unknown parameter we seek to estimate, and we assume there is a set of active nodes $\mathcal{G}\subset[N]$ with $|\mathcal{G}|=G$, so that $B^{\star}_{ij}=0$ if $i\in\mathcal{G}^{C}$ or $j\in\mathcal{G}^{C}$. We allow some edge weights inside the subgraph defined by $\mathcal{G}$ to be zero, but we focus on whether the set $\mathcal{G}$ is correctly estimated by the set $\hat{\mathcal{G}}$ of active nodes in $\hat{B}$. Denote by $\mathcal{M}\subseteq^{N\times N}$ the set of matrices where the only non-zero coefficients appear in the active subgraph, that is,

[TABLE]

There are two main assumptions on the loss function $\ell$ required for consistent selection in high-dimensional models (Lee et al.,, 2015). The first assumption is on the convexity of the loss function around $B^{\star}$, while the second assumption bounds the size of the entries in the loss Hessian between the variables in the active subgraph and the rest. Let the loss Hessian $\nabla^{2}\ell(B^{\star})\in^{N\times N}\otimes^{N\times N}$ be defined by

[TABLE]

and define the matrix $H_{(i,j),\mathcal{G}}\in^{G\times G}$ with $(i,j)\in\left(\mathcal{G}\times\mathcal{G}\right)^{C}$ such that

[TABLE]

where $\Lambda\in^{G\times G}\otimes^{G\times G}$ is a tensor such that $\text{Mat}(\Lambda)$ is a pseudoinverse of $\text{Mat}\left(\nabla^{2}_{(\mathcal{G},\mathcal{G}),(\mathcal{G},\mathcal{G})}\ell(B^{\star})\right)$, and Mat is the operation that unfolds the entries of a tensor $\Lambda$ into a $G^{2}\times G^{2}$ matrix. The matrix $H_{(i,j),\mathcal{G}}$ measures how well the variable corresponding to the edge $(i,j)$ can be represented by the variables in the active subgraph.

Assumption 1* (Restricted Strong Convexity).*

There exists a set $C\subset^{N\times N}$ with $B^{\star}\in C$, and constants $m>0$, $\tilde{L}<\infty$ such that

[TABLE]

Assumption 2* (Irrepresentability).*

There exists a constant $0<\tau<1$ such that

[TABLE]

This version of the irrepresentability condition corresponds to the one usually employed in group lasso penalties (Bach,, 2008), but as we will see later, due to overlaps in the groups it further requires a lower bound on $\rho$ to work for model selection.

The first two assumptions are stated directly as a function of the loss for a fixed design case, but they can be substituted with bounds in probability for the case of random designs. In order to obtain rates of convergence, we do require a distributional assumption on the first derivative of the loss. This assumption can be substituted with a bound on $\max_{i}\|\nabla\ell(B^{\star})_{(i)}\|_{2}$ (see Lemma B.6 in the Appendix).

Assumption 3* (Sub-Gaussian score function).*

Each pair in the sample $(A,Y)$ is independent and comes from a distribution such that the entries of the matrix $\nabla\tilde{\ell}(Y,A;B^{\star})$ are subgaussian. That is, for all $t>0$ there is a constant $\sigma^{2}>0$ such that

[TABLE]

With these assumptions, we establish consistency and correct model selection. The proof is given in the Appendix..

Proposition 1*.*

Suppose Assumptions 1 and 3 hold.

1. (a)

Setting the penalty parameters as $\lambda=c_{1}\sigma\sqrt{\frac{\log N}{n}}\min\left\{\frac{\sqrt{N}}{1+\rho},\frac{1}{\rho}\right\}$ and $\rho\geq 0$ for some constant $c_{1}>0$, with probability at least $1-2/N$ the optimal solution of (2.4) is unique and satisfies

[TABLE] 2. (b)

Suppose Assumption 2 also holds. If $n>c_{2}G^{2}\sigma^{2}\log N$ for a constant $c_{2}>0$, setting the penalty parameters as $\lambda=c_{3}\sigma\sqrt{\frac{\log N}{n}}\min\left\{\frac{\sqrt{N}}{1+\rho},\frac{1}{\rho}\right\}$ for some constant $c_{3}>0$, and

[TABLE]

then

[TABLE]

The part of the penalty associated with $\rho$ causes the solution to be sparse. Due to the overlap in the groups, a small value of $\rho$ will usually not result in zeros in the solution of the problem (2.4). The lower bound on $\rho$ in (4.4) ensures that the irrepresentability condition of Lee et al., (2015) holds (see Lemma 2 in the Appendix).

Proposition 1 ensures that, with high probability, all edges estimated to have non-zero weights are contained in the active subgraph, and quantifies the error in estimating the entries of $B^{\ast}$. To ensure that all active nodes are recovered, at least one edge corresponding to each active node needs to have a non-zero weight. A similar result can be obtained to guarantee recovery of all active nodes under a stronger assumption that the magnitude of the non-zero entries of $B^{\star}$ is bounded below by $|B_{ij}^{\star}|>c_{4}G^{2}\lambda$ for a constant $c_{4}$. The condition in part (b) requires a sample size that grows faster than the size of the active subgraph, which in practice can be much smaller than the size of the network, making the method suitable for our applications in which the sample size is limited and the number of nodes is large.

5 Numerical results on simulated networks

In this section, we evaluate the performance of our method using synthetic networks. We are interested in assessing the ability of the method to correctly identify predictive edges, its classification accuracy, and comparisons to benchmarks. We compare the different methods’ edge selection performance in simulations using area under the curve (AUC). Brain connectomic networks are characterized by organization of nodes into communities (Bullmore and Sporns,, 2009), in which nodes within the same community tend to have stronger connections than nodes belonging to different communities. In order to generate synthetic networks that mimic this property, we introduce community structure using the stochastic block model (SBM) (Holland et al.,, 1983). Before generating edges, we assign each node a community label, $C_{i}$, where $C_{i}\in[K]$ for each $i\in[N]$. The node assignments are the same for all networks in the population. Given the community labels, network edges are generated independently from a distribution that only depends on the community labels of the nodes associated with each edge. Since fMRI networks are real-valued networks, we generate edge weights from a Gaussian distribution, rather than the standard Bernoulli distribution normally used with the SBM. Specifically, we draw each $A_{ij}$ independently from $N(\mu_{C_{i}C_{j}},\sigma^{2})$, with $\mu\in\mathbb{R}^{K\times K}$ defined by

[TABLE]

and $\sigma^{2}=0.18$. These values were chosen to approximately match the distribution of edge weights in our datasets (see Section 6). We set the number of nodes $N=300$, with $K=12$ communities of size $25$ each. We work with undirected networks, so the adjacency matrices are symmetric, with 44,850 distinct edges. Although our method is able to scale to larger networks, this moderately sized setting is already highly computationally demanding for many of the comparison benchmarks.

To set up two different class distributions, we select a set of active nodes $\mathcal{G}$ first, defined by the nodes corresponding to some communities selected at random. Then, we alter a set of differentiating edges $\mathcal{E}$ selected at random from $\mathcal{G}\times\mathcal{G}$ with probability $p$. For each edge $(i,j)\in\mathcal{E}$, the distribution in class $Y=-1$ is $N(\mu_{C_{i}C_{j}},\sigma^{2})$, while the distribution in class $Y=1$ is changed to $N(0.2,\sigma^{2})$. Figure 3 shows example expected adjacency matrices for each class. We then generate 50 networks from each class, resulting in a sample size of $n=100$. We vary $G=|\mathcal{G}|$ by changing the number of communities selected, and the value of $p$, to study the effect of the number of active nodes and the density of differentiating edges inside a subgraph. The number of communities selected is set to 1 ($|\mathcal{G}|=25$), 2 ($|\mathcal{G}|=50$) and all communities ($|\mathcal{G}|=300$); note that in the last scenario all nodes are active and hence the node structure is not informative at all.

Since we are interested is identifying predictive edges and nodes, we use the AUC of the receiver operating characteristic (ROC) curve, for both edge and node selection. For each method, we calculate the ROC curve by changing its corresponding sparsity parameter to vary the number of edges selected . For a selection method $\mathcal{M}$ and a sparsity parameter $\eta$ let $\hat{\mathcal{E}}(\mathcal{M},\eta)$ be the set of edges selected by $\mathcal{M}$, and $\hat{\mathcal{G}}(\mathcal{M},\eta)$ the set of active nodes corresponding to $\hat{\mathcal{E}}(\mathcal{M},\eta)$. We calculate the edge false positive rate (EFPR) and edge true positive rate (ETPR) as

[TABLE]

The node FPR and TPR are calculated similarly.

We also evaluate the prediction accuracy of the methods. For each method, we use 5-fold cross-validation to select the best tuning parameter using the training data, and then compute the test error on a different dataset simulated under the same settings. The AUC and test errors reported are averaged over 30 replications.

Methods for benchmark comparisons on simulated networks were selected based on their good performance on real data (see Section 6). For our method (GC), we vary the parameter $\rho$ and compare results for two different values of $\lambda$, $.05$ (GC1) and $10^{-4}$ (GC2). For unstructured regularized logistic regression, we use the elastic net (Friedman et al., 2010b, ), with a fixed $\alpha=0.02$ (ENet). The performance of elastic net is not very sensitive to different values of $\alpha$, but the number of variables that the method is able to select with large values is limited (including the case of $\alpha=1$ that corresponds to the lasso). A support vector machine with $\ell_{1}$ penalty (Zhu et al.,, 2004; Becker et al.,, 2009) is also included (SVML1) for comparison, and additionally we evaluate the classification error of the original support vector machines (SVM) (Cortes and Vapnik,, 1995). For both SVMs, we use linear kernels, which performed better than nonlinear ones. We also consider an independent screening method for variable selection based on the two sample $t$-statistic (T-stat). Finally, we also compare with the signal-subgraph method (SS) (Vogelstein et al.,, 2013) which is the only other method that takes into account the network structure of the predictor variables. Note that the signal subgraph is designed for binary networks, so in order to apply it we thresholded each edge at the population mean. For each method, we fit 10 different tuning parameters to change the sparsity of the solution.

Figure 4 shows the values of the average AUC for selecting edges (top) and nodes (bottom). For $G=25$ and $50$, as the proportion of differentiating edges in the active subgraph increases, methods that take into account network structure (GC1, GC2 and SS) slightly improve their edge AUC, since enforcing node selection also results in better edge selection, while the edge AUC remains constant for unstructured methods (ENet, T-stat and SVML1). On node selection, all methods improve the node AUC as the fraction of significant edges increases, but GC and SS have the largest gains. A similar trend is observed in classification error shown in Figure 5. All methods improve as the proportion of differentiating edges increases, but our method has the best performance overall. Our method performed the best with the larger value of $\lambda$ (GC1) on variable selection, particularly when the set of active nodes is smaller, but both values of $\lambda$ give very good classification performance. In the last scenario ($G=300$), all nodes are active so the node AUC is undefined, and the node structure is not informative at all. Although the performance of our method is no longer the best, it performs comparably to state of the art methods that do not use network structure.

In terms of computing time, since there are many contributing factors including the software choice for implementation and the tuning parameters, a fair comparison is difficult. We can roughly say that elastic net is the fastest, taking about a minute to run a cross-validation instance, while our method takes about 10 minutes on average, and the signal-subgraph takes more than an hour.

6 Application to schizophrenia data

We analyze the performance of the classifier on the two different brain fMRI datasets previously described in Section 1. The code of our classifier and the processed connectomics datasets are available at https://github.com/jesusdaniel/graphclass.

6.1 Classification results

First, we evaluate our method’s classification accuracy. We use a nested 10-fold cross-validation to choose tuning parameters and estimate the test accuracy. The classifier is trained for a range of values of $\lambda$ and $\rho$, with $\lambda\in\left\{10^{-7},10^{-6.5},\ldots,10^{-2}\right\}$ and $\rho\in\left\{10^{-3},10^{-2.5},\ldots,10^{2}\right\}$. The value of $\gamma$ in (2.4) is set to $10^{-5}$; we observed that setting $\gamma$ to a small value speeds up convergence without affecting the accuracy or sparsity of the solution. Figure 6 shows the average cross-validated accuracy, sparsity (fraction of zero coefficients) and node sparsity (fraction of inactive nodes), as a heat map over the grid of tuning parameter values. We observe that $\lambda$ has little influence on sparsity, which is primarily controlled by $\rho$. Moreover, as Proposition 1 suggests, values of $\rho<1$ do not result in node selection. As expected, accuracy generally decreases as the solution becomes sparser, which is not uncommon in high-dimensional settings (Hastie et al.,, 2015). However, we can still achieve excellent accuracy with a substantially reduced set of features. In the COBRE dataset, the best accuracy is obtained with only 1886 edges (5.4%) but almost all nodes are active (260). On the UMich data, 29733 edges (85.6%) achieve the best performance, and all nodes are active. Choosing parameters by cross-validation often tends to include too many noise variables (Meinshausen,, 2007), as we also observed in simulations. A commonly used technique to report solutions that still achieve good accuracy with a substantially reduced set of features is the so-called “one-standard-error rule” (Hastie et al.,, 2015), in which one selects the most parsimonious classifier with cross-validation accuracy at most one standard error away from the best cross-validation accuracy. Figure 7 shows the solutions for each dataset obtained by this rule. Nodes are ordered by brain systems (see Figure 1). The fitted solution for COBRE has 549 non-zero coefficients (1.56%) and 217 active nodes (82.5%), while the UMich solution has 11748 non-zero entries (33.8%), and all nodes are active. Note that when many variables are selected, the magnitude of the coefficients becomes small due to the grouping effect of the penalty (Zou and Hastie,, 2005).

We also compared our method to benchmarks (Table 1), using the same methods as in the previous section and training and evaluating all methods using with the same nested 10-fold cross-validation. For SVM, we tested different kernels, including graph aware kernels (Gärtner et al.,, 2003), but in most cases local kernel methods were no better than random guessing. We additionally included random forests and a method based on global and local network summaries previously proposed as features for classifying brain data (Prasad et al.,, 2015). For the latter, because our dataset is much larger, we only considered global and node features proposed in Prasad et al., (2015), which resulted in about 30,000 features per individual, and omitted edge features. Watanabe et al., (2014) evaluated their classifiers on a different parcellation of the COBRE data, and we do not include their methods since they are based on the assumption of equally spaced nodes and cannot be directly applied to our data. Their reported accuracy of 71.9% and 73.5% for the COBRE data is substantially lower than our method, although the results are not directly comparable.

Results in Table 1 show that most methods performed better on the COBRE dataset than on the UMich dataset, which can be partially explained by the different sample sizes and possibly noise levels. Besides differences in sample size and demographic characteristics (Table 4), the COBRE dataset is more homogeneous as it was collected using identical acquisition parameters, whereas the UMich dataset was pooled across five different experiments spanning seven years.

Our method performs very well on both datasets, particularly among methods that do variable selection. SVMs, which use the hinge loss, perform well too, and generally outperform methods using the logistic loss. Our penalty can be combined with any loss, so we could also include our penalty combined with hinge loss which might potentially improve classification accuracy, but we do not pursue this direction, for two reasons: one, our method is close to SVM + L1 as it is (better on COBRE, slightly worse on UMich but the difference is within noise levels), and because solutions based on logistic loss are generally considered more stable and preferable for variable selection (Hastie et al.,, 2015). In Figure 8, we plot cross-validated classification accuracy of these methods as a function of the number of variables selected. For the COBRE data, as we have observed before, good accuracy can be achieved with a fairly small number of edges, and the noisier UMich data requires more edges. In all cases, our method uses fewer nodes than the others, as it is designed to do so.

Ultimately, assessing significance of the selected variables is necessary, which is in general a difficult task in high-dimensional settings and an active area of research (see for example Meinshausen and Bühlmann, (2010); Van de Geer et al., (2014); Lockhart et al., (2014); Lee et al., (2016)). In brain connectomics, it is particularly challenging to identify significant variables because of small sample sizes (Button et al.,, 2013). Here we employ stability selection (Meinshausen and Bühlmann,, 2010) which can be shown to control a type of false discovery rate by employing many rounds of random data splitting and calculating the probability of each variable being selected. Some versions of this method have been theoretically studied, and upper bounds on the expected number of variables with a low selection probability that are included in the final solution (i.e., errors) have been derived under mild conditions (Meinshausen and Bühlmann,, 2010; Shah and Samworth,, 2013). We implemented the version of stability selection proposed by Shah and Samworth, (2013), with values of $\lambda$ and $\rho$ obtained by cross-validation on the COBRE data, and by the “one standard error rule” on the UMich dataset, since stability selection is most relevant to sparse solutions. However, one of the advantages of stability selection is that it is not sensitive to the initial choice of tuning parameters, and changing tuning parameters only slightly alters the ordering of variables with the largest selection probabilities.

The edges with the 15 largest selection probabilities are reported in Table 2. Using the results of Shah and Samworth, (2013) (equation 8), we estimated that the expected number of falsely selected variables (variables with a probability of selection smaller than the estimated) is bounded by 6.1 for the COBRE dataset and 9.7 for the UMich data, which also suggests that results on the UMich data might be less reliable. While the two datasets yield somewhat different patterns of edge selection, it is notable that the default mode network (5) was often selected in both. This network has been consistently implicated in schizophrenia (Whitfield-Gabrieli et al.,, 2009; Öngür et al.,, 2010; Peeters et al.,, 2015), as well as other psychiatric disorders, possibly as a general marker of psychopathology (Broyd et al.,, 2009; Menon,, 2011). In the COBRE dataset, edges were also selected from the fronto-parietal task control region (8), previously linked to schizophrenia (Bunney and Bunney,, 2000; Fornito et al.,, 2012). These results coincide with the findings of Watanabe et al., (2014) on a different parcellation of the same data, which is an encouraging indication of robustness to the exact choice of node locations. Some of the variables with the highest estimated selection probabilities appear in the uncertain system (-1), in particular in the cell connecting it with salience system (9), which suggests that alternative parcellations that better characterize these regions may offer a better account of the schizophrenia-related changes. Additionally, sensory/somatomotor hand region (1) and salience system (9) also stand out in the UMich data, and these are networks that have also been implicated in schizophrenia (Dong et al.,, 2017).

While results in Table 2 do not fully coincide on the two datasets, there are clear commonalities. Table 3 compares classification accuracy when the classifier is trained on one dataset and tested on the other (with the exception of the intercept, since the datasets are not centered in the same way, which is fitted on a part of the test data, and the test error is then computed via 10-fold cross-validation). While the accuracy is lower than when the same dataset is used for training and testing, as one would expect, it is still reasonably good and in fact better than some of the benchmark methods even when they train and test on the same data. We again observe that the COBRE dataset is easier to classify.

Figure 9 shows the active nodes in the COBRE dataset (marked in green), corresponding to the endpoints of the edges listed in Table 2. We also identified a set of 25 nodes that are not selected in any of the sparse solutions with cross-validation accuracy within one standard error from the best solution (marked in purple). These consistently inactive nodes are mostly clustered in two anatomically coherent regions.

7 Discussion

We have presented a method for classifying graphs with labeled nodes, motivated by brain connectomics but generally applicable to any setting with such graphs. The distinct feature of our method is that it is graph-aware, aiming to select a sparse set of both edges and nodes, but it is general in the sense that it does not rely on the spatial structure of the brain. The method is computationally efficient since the regularization we use is convex, and the solution is implemented with efficient optimization algorithms. These properties guarantee fast convergence to the solution, making the methods scalable to networks with thousands of nodes, which is enough to deal with many of the brain atlases usually employed in neuroimaging (see for example Kiar et al., (2018)). Statistically, the rate of convergence depends on the number of active nodes only, not the total number of nodes, which allows for accurate results with even moderate sample sizes if the active node set is small.

The results we obtained on the schizophrenia data are generally in agreement with previous studies. In particular, the default mode network has been consistently implicated in schizophrenia and many other psychiatric disorders (Öngür et al.,, 2010; Broyd et al.,, 2009). While different subnetworks were implicated by the two different datasets, we are still able to predict the disease status fairly accurately by training on one dataset and testing on the other. The differences between the two datasets may reflect real differences in samples collected at different sites and in different experiments, as significant pathophysiological heterogeneity occurs for all psychiatric diagnoses, or they may simply reflect type 2 errors.

Our methods work very generally with a sample of networks with labeled nodes and associated responses. The many pre-processing steps inevitable when dealing with fMRI data always add some uncertainty, and pre-processing decisions can potentially affect downstream conclusions. We aimed to somewhat mitigate this by using ranks, which are more robust and showed a slightly better performance on our datasets. Another option, when practical, is to compare multiple pre-processing pipelines, and/or multiple measures of connectivity, to further validate results. Our method’s independence of these particular choices and its computational efficiency make it an attractive option for such comparisons.

Acknowledgments

This research was supported in part by NSF grant DMS-1521551, ONR grant N000141612910, and a Dana Foundation grant to E. Levina, NSF training grant DMS-1646108 support for D. Kessler, as well as by computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. S. F. Taylor’s research is supported by the National Institute of Mental Health (R01MH064148, R21MH086701, R21MH101676), the Boledovich Schizophrenia Research Fund and University of Michigan Clinical Translational Science Award (UL1RR024986).

Appendix A Optimization algorithm details

The optimization procedure for solving the penalized prediction problem introduced above consists in a proximal algorithm, and the steps are detailed in Algorithm 1. Each step requires to solve a further convex optimization problem via ADMM. The exact solution of the steps of this method is shown in Algorithm 2.

Appendix B Proofs

Here we prove the bounds on Frobenius norm error and probability of support selection in Proposition 1, following the framework of Lee et al., (2015) based on geometrical decomposability. A penalty $\Omega$ is geometrically decomposable if it can be written as

[TABLE]

for all $B$, with $A,I$ closed convex sets, $E$ a subspace, and $h_{C}$ the support function on $C$ defined as $h_{C}\left(B\right)=\sup\left\{\left\langle Y,B\right\rangle\left|Y\in C\right.\right\}$ .

The proof proceeds in the following steps. Lemma 1 shows that an equivalent form of our penalty (2.2) is geometrically decomposable, allowing us to use the framework of Lee et al., (2015). Lemma B.6 shows the Assumption 2 together with a lower bound on $\rho$ imply that the irrepresentability assumption of Lee et al., (2015) holds. Assumption 2 is directly on the entries of the loss Hessian, which simplifies the very general form of the assumption in Lee et al., (2015). Lemma B.6 gives a bound on the entries of the loss gradient under the sub-gaussianity assumption 3. Lemma 4 gives explicit bounds for the compatibility constants that appear on Theorem 1 of Lee et al., (2015). Finally, we combine these results to prove Proposition 1.

Without loss of generality, to simplify notation we assume that $\mathcal{G}=\{1,\ldots,G\}$, that is, the active subgraph is in the first $G$ rows of the matrix.

Lemma 1**.**

The penalty (2.2) can be written as geometrically decomposable.

Proof of Lemma 1.

We use an equivalent formulation of the penalty in which every coefficient is penalized only once. Let $B^{\prime},B^{\prime\prime}\in^{N\times N}$ be matrices such that the upper triangular part of $B^{\prime\prime}$ and the diagonals of $B^{\prime}$ and $B^{\prime\prime}$ are zero. Define

[TABLE]

and $E=\{(B^{\prime},B^{\prime\prime})\in^{N\times 2N}\ :\ B^{\prime}=B^{\prime T},B^{\prime\prime}_{ij}=B^{\prime}_{ij},\text{ for }i<j\text{ and }B^{\prime\prime}_{ij}=0\text{ for }i\geq j\}$. Denote by $R$ the transformation from N×N to N×2N that replicates entries appropriately,

[TABLE]

Therefore, for any $B\in^{N\times N}$, we can uniquely define $RB=(B^{\prime},B^{\prime\prime})$ such that $\Omega(B)=\tilde{\Omega}(B^{\prime},B^{\prime\prime})$. We then show that $\tilde{\Omega}$ is geometrically decomposable. Moreover, for any $(B^{\prime},B^{\prime\prime})\in E$ we can define $R^{-1}$, so the penalties $\Omega$ and $\tilde{\Omega}$ on $E$ are equivalent. Define the sets $A,I\subset^{N\times 2N}$ such that

[TABLE]

Letting $\left\langle Y,(B^{\prime},B^{\prime\prime})\right\rangle=\text{Tr}(Y^{\prime}B^{\prime T})+\text{Tr}(Y^{\prime\prime}B^{\prime\prime T})$, combining the arguments of Lee et al., (2015) for lasso and group lasso penalties,

[TABLE]

Hence, $\Omega$ can be written as a geometrically decomposable penalty

[TABLE]

∎

We introduce some notation in order to state the irrepresentability condition of Lee et al., (2015). For a set $F\subset^{N\times 2N}$ and $Y\in^{N\times 2N}$, denote by $\gamma_{F}\left(Y\right)=\inf\left\{\lambda>0:Y\in F\right\}$ the gauge function on $C$. Thus,

[TABLE]

where $\textbf{1}_{I}(B)=0$ if $B\in I$ and $\infty$ otherwise. Define

[TABLE]

for $Z\in^{N\times 2N}$. Let $\tilde{\mathcal{M}}=E\cap\text{span}(I)^{\perp}$ be the set of matrices with correct support in the extended space N×2N, similarly to $\mathcal{M}$ in (4.1). Denote by $P_{M}$ and $P_{M^{\perp}}$ the projections onto $\tilde{\mathcal{M}}$ and $\tilde{\mathcal{M}}^{\perp}$. Define the function $\mathcal{H}(Z):^{N\times N}\rightarrow^{N\times N}$ as

[TABLE]

where $H_{(i,j),\mathcal{G}}$ is the matrix defined in (4.2). The Irrepresentability Assumption 3.2 of Lee et al., (2015) requires the existence of $0<\tilde{\tau}<1$ such that

[TABLE]

For a support function $h$, denote by $\partial h(M)=\bigcup_{Y\in M}\partial h(Y)$ the set of subdifferentials of $h$ in $M$. Note that $\partial h_{A}(M)=A$, since $0\in M$ and $\partial h_{A}(0)=A$.

Lemma 2**.**

If Assumption 2 holds and $\rho>\frac{1}{\tau}-\frac{1}{\sqrt{G}}$, then there exists $0<\tilde{\tau}<1$ such that (B.2) holds.

Proof of Lemma 2..

Since $V$ is sublinear (Lemma 3.3 of Lee et al., (2015)),

[TABLE]

To bound the first term, note that $E^{\perp}=\{(Z^{\prime},Z^{\prime\prime})|Z_{ij}^{\prime}+Z_{ji}^{\prime}+Z_{ij}^{\prime\prime}=0,j<i\}$. Hence,

[TABLE]

Therefore,

[TABLE]

which implies that

[TABLE]

Let $Z=(Z^{\prime},Z^{\prime\prime})\in A$. Without loss of generality, assume that $Z^{\prime}_{\mathcal{G},\mathcal{G}}=Z^{\prime\prime}_{\mathcal{G},\mathcal{G}}=0$ (note that these entries do not change $V(P_{M^{\perp}}Z)$). Therefore, $P_{M^{\perp}}Z=Z$. Hence,

[TABLE]

The last bound from $|Z^{\prime}_{ji}|\leq 1$ and no longer depends on $Z$. It is easy to see that the minimum is attained when, for each $i>G$,

[TABLE]

and therefore

[TABLE]

Moreover, if $Z^{\ast}\in A$ is defined such that $(Z^{\ast})^{(1)}_{G+1,i}=1$ for $i=1,\ldots,G$ and [math] elsewhere, then $V(Z^{\ast})$ achieves this bound, which shows that $\rho>1-\frac{1}{\sqrt{G}}$ is a necessary condition for the irrepresentability to hold, even in the case where the entries of the Hessian that denote the information between active and inactive edges is zero. Therefore, plugging the bounds (B.4) and (B.5) into equation (B.3), we obtain (B.2) holds as long as $1-\tau+\frac{\sqrt{G}}{1+\rho\sqrt{G}}<1$, which implies that $\rho>\frac{1}{\tau}-\frac{1}{\sqrt{G}}$. ∎

The next lemma establishes a bound on the dual norm of $\Omega$ of the loss gradient function under a sub-Gaussian assumption. Let $\Omega^{\ast}$ denote the dual norm of $\Omega$, so $\Omega^{\ast}(B)=\sup\left\{\left\langle Y,B\right\rangle\left|\ Y\in C,\Omega(Y)\leq 1\right.\right\}$.

Lemma 3**.**

Under Assumption 3,

[TABLE]

Proof of Lemma B.6.

By Hoeffding’s inequality for sub-Gaussian variables, for all $j,k$ and $t>0$,

[TABLE]

Note that $(1+\rho)\sum_{i=1}^{N}\|B_{(i)}\|_{2}\leq\Omega(B)$. Let $\Phi(B)=\frac{1}{1+\rho}\max_{i=1,\ldots N}\|B_{(i)}\|_{2}$. Thus,

[TABLE]

In a similar manner, $\rho\|B\|_{1}\leq\Omega(B)$. Setting $\Xi(B)=\frac{1}{\rho}\|B\|_{\infty}$, we have

[TABLE]

Using (B.7) and setting $\Lambda=\nabla\ell(B^{\star})$,

[TABLE]

the last inequality obtained by arguments similar to Lemma 4.3 of Lee et al., (2015). In the same way, we can also bound the previous quantity using (B.8) by

[TABLE]

Combining (B) and (B), we can obtain equation (B.6). ∎

For a semi-norm $\Psi:^{N\times N}\rightarrow$, let $\kappa_{\Omega}$ be the compatibility constant between $\Psi$ and the Frobenius norm, defined as

[TABLE]

and let $\kappa_{\text{IC}}$ be the compatibility constant between the irrepresentable term and the dual norm $\Omega^{\ast}$ given by

[TABLE]

Lemma 4**.**

The following bounds on the compatibility constants hold:

[TABLE]

Proof of Lemma 4..

Note that $\Omega(Y)$ is maximized on $\{Y:\|Y\|_{2}\leq 1\}$ when all entries of $Y$ have magnitude equal to $\frac{1}{\sqrt{G(G-1)}}$. Therefore

[TABLE]

Similarly, (B.7) implies

[TABLE]

Finally,

[TABLE]

∎

Proof of Proposition 1.

Part (a). Since $\hat{B}$ minimizes the problem (2.4),

[TABLE]

Rearranging the terms, using Assumption 1, by the triangle inequality and the generalized Cauchy-Schwarz inequality,

[TABLE]

By the argument for computing $\kappa_{\Omega}$ in (4), $\Omega(Y)\leq\{\sqrt{N}+\rho\sqrt{N(N-1)}\}\|Y\|_{2}$. Rearranging the terms in (B.11),

[TABLE]

For any $\rho$, setting $\lambda=2\sqrt{\frac{\sigma^{2}\log N}{n}}\min\left\{\frac{\sqrt{N}}{1+\rho},\frac{1}{\rho}\right\}$, by Lemma B.6, with probability at least $1-2/N$,

[TABLE]

Part (b). Lemma 1 gives a geometric decomposition of the penalty. Therefore, we can directly use Theorem 3.1 of Lee et al., (2015), since Lemma 2 also ensures that their irrepresentability condition holds. Thus,

[TABLE]

and $\hat{\mathcal{G}}\subseteq\mathcal{G}$ as long as

[TABLE]

Setting

[TABLE]

using a similar argument than (B.12), the left hand size of (B.13) holds with probability at least $1-2/N$. The right hand side of (B.13) holds as long as the sample size satisfies

[TABLE]

with $C(L,m,\tau,\kappa_{\Omega},\kappa_{\Omega^{\ast}},\kappa_{IC})>0$ a positive constant. Therefore, claims (4.5) and (4.6) follow.

∎

Appendix C Data aquisition and pre-processing

C.1 Subjects and imaging

The COBRE data

Raw anatomic and functional scans from 146 subjects (72 psychosis patients and 74 healthy control subjects) were downloaded from a public database (http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html). Four subjects coded as ambidextrous (2 patients, 2 controls) were excluded to yield 70 psychosis patients and 72 controls for analysis. To enter the COBRE dataset, subjects had a diagnosis of either schizophrenia or schizoaffective disorder and were without histories of neurological disorder, mental retardation, severe head trauma with more than 5 minutes loss of consciousness and substance abuse/dependence within the last 12 months.

In the primary sample, two schizophrenic (SCZ) subjects and one healthy control (HC) subject had insufficient voxels in the cerebrospinal fluid (CSF) segmentation on the CSF, and they were dropped from additional analyses. Two additional SCZ subjects were excluded for scrub ratios (see discussion of scrubbing routine in fMRI Data Analysis) greater than 0.6, leaving 38 SCZ subjects and 42 HC subjects for the analysis. In the replication sample, 15 psychosis patients and two control subjects were excluded for scrub ratios greater than 0.6; one patient was excluded with incomplete data, leaving 54 SCZ and 69 HC subjects for analysis (see Table 4).

A full description of the imaging parameters for the COBRE dataset is available online at the link provided above and in Aine et al., (2017).

The UMich data

Subjects were selected from experiments conducted by Professor Stephan F. Taylor at the University of Michigan between 2004 and 2011 for task-based fMRI studies, which included resting state scans. Thirty-nine stable outpatients were selected with DSM-IV schizophrenia or schizoaffective disorder (SCZ) (American Psychiatric Association,, 1994). Forty healthy comparison (HC) subjects, without a lifetime history of Axis I psychiatric disorders (First et al.,, 1995), were selected to approximate the age range, gender distribution and family education level of the patients. Prior to initial data collection, all subjects gave written, informed consent to participate in the protocol approved by the University of Michigan institutional review board (IRBMED).

MRI scanning occurred on a GE 3T Signa scanner (LX [8.3] release, General Electric Healthcare, Buckinghamshire, United Kingdom). Functional images were acquired with a T2*-weighted, reverse spiral acquisition sequence (gradient recalled echo, TE=30 msec, FA=90 degrees, field of view=22 cm, 40 slices, 3.0mm thick/0mm skip, equivalent to 64 x 64 voxel grid – yielding isotropic voxels 3 mm on edge). Because the data were acquired across different experiments, acquisition parameters differed slightly in the aggregate sample: 240 volumes at TR=1500 msec (11 SCZ, 10 HC), 180 volumes at TR=2000 msec (17 SCZ, 16 HC) and 240 volumes at 2000 msec (14 SCZ, 17 HC). Acquisitions were acquired in the resting state with eyes open and fixated on a large ‘plus’ sign projected on a monitor.

C.2 Pre-processing

We first performed standard pre-processing steps. All scans were slice-time corrected and realigned to the 10th image acquired during a scanning session (Jenkinson et al.,, 2002). Subsequent processing was performed with the Statistical Parametric Mapping SPM8 package (Wellcome Institute of Cognitive Neurology, London). Anatomic normalization was done with the VBM8 toolbox in SPM8, using the high resolution structural scans obtained for both datasets. Normalizing warps were applied to the co-registered, functional volumes, which were re-sliced and smoothed with an 8 mm isotropic Gaussian smoothing kernel. To assess and manage movement, we calculated the frame-wise displacement (FD) (Power et al.,, 2012), for all 6 parameters of rotation and translation. We used a scrubbing routine to censor any frame with FD $>$ 0.5 mm from the regression analysis described below, yielding a scrub ratio for each subject. Three-compartment segmentation of the high-resolution structural image from the VBM8 normalization was applied to the functional time series to extract cerebral spinal volume (CSF) and white matter (WM) compartments, which were then subjected to a principal component analysis to identify the top 5 components in each (Behzadi et al.,, 2007), which should correspond to heart rate and respiratory effects on global signal (Chai et al.,, 2012). Multiple regressions were applied to the time series to remove the following nuisance effects: linear trend, 6 motion parameters, their temporal derivatives, the quadratics of these 12 parameters, 5 components from the PCA of CSF, 5 components of PCA of WM, followed by band pass filtering from 0.01 – 0.1 Hz, and then motion scrubbing. For each 4D data set, time courses were then extracted from 10 mm diameter spheres based on the 264 sets of coordinates identified by Power et al., (2011). From these time series, a cross-correlation matrix of Pearson r-values was obtained and Fisher’s R-to-Z transformation was applied for each of the 264 nodes with every other node (for COBRE dataset, node 75 is missing). Finally, for each individual, edge weights were assigned to be ranks of these score, with edge scores ranked separately for each subject, and then these values were centered and standardized across the individuals. Ranks have been used previously in brain connectomic studies to reduce the effect of potential outliers (Yan et al.,, 2013); although some information is lost with the rank transformation, we observed that while ranks do not increase the classification accuracy significantly, they tend to produce sparser solutions with a similar accuracy to Pearson correlations.