{'fraction_non_alphanumeric': 0.03561246436252689, 'fraction_numerical': 0.016155654479067673, 'mean_word_length': 4.581658291457287, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High Dimensional Bayesian Network Classification with Network Global-Local Shrinkage Priors

This article proposes a novel Bayesian classification framework for networks with labeled nodes. While literature on statistical modeling of network data typically involves analysis of a single network, the recent emergence of complex data in several biological applications, including brain imaging studies, presents a need to devise a network classifier for subjects. This article considers an application from a brain connectome study, where the overarching goal is to classify subjects into two separate groups based on their brain network data, along with identifying influential regions of interest (ROIs) (referred to as nodes). Existing approaches either treat all edge weights as a long vector or summarize the network information with a few summary measures. Both these approaches ignore the full network structure, may lead to less desirable inference in small samples and are not designed to identify significant network nodes. We propose a novel binary logistic regression framework with the network as the predictor and a binary response, the network predictor coefficient being modeled using a novel class global-local shrinkage priors. The framework is able to accurately detect nodes and edges in the network influencing the classification. Our framework is implemented using an efficient Markov Chain Monte Carlo algorithm. Theoretically, we show asymptotically optimal classification for the proposed framework when the number of network edges grows faster than the sample size. The framework is empirically validated by extensive simulation studies and analysis of a brain connectome data.


Introduction
Of late, the statistical literature has paid substantial attention to the unsupervised analysis of a single network, thought to be generated from a variety of classic models, including random graph models Erdos and Rényi (1960), exponential random graph models Frank and Strauss (1986), social space models Hoff et al. (2002); Hoff (2005Hoff ( , 2009) and stochastic block models Nowicki and Snijders (2001). These models have found prominence in social networking applications where the nodes in the network are exchangeable. However, there are pertinent biological and physiological applications in which network nodes are labeled and a network is available corresponding to every individual. Section 6 presents one such example from a brain connectome study, where brain networks are available for multiple individuals who are classified as subjects with high or low IQ (Intelligence Quotient). In this study, the human brain has been divided according to the Desikan Atlas Desikan et al. (2006) that identifies 34 cortical regions of interest (ROIs) both in the left and the right hemispheres of the human brain, implying 68 cortical ROIs in all. A brain network for each subject is represented by a symmetric adjacency matrix whose rows and columns are labeled corresponding to different ROIs (shared among networks corresponding to all individuals) and entries correspond to estimates of the number of fibers connecting pairs of brain regions. The scientific goal in this setting pertains to developing a predictive rule for classifying a new subject as having low or high IQ based on his/her observed brain network with labeled nodes. Additionally, it is of specific interest for neuroscientists to identify influential brain regions (nodes in the brain network) and significant connections between different brain regions predictive of IQ. Guha and Rodriguez (2020) discuss the network regression problem with a continuous response and an undirected network predictor. However, there are pertinent biological and physiological studies where a network along with a binary response is obtained for each subject. The goal of these studies is usually to classify the networks according to the binary response and predict the associated binary response from a network. We refer to this problem as the network or graph classification problem. Additionally, Guha and Rodriguez (2020) focus on a specific network shrinkage prior, whereas this article generalizes the inference to a class of network global-local shrinkage priors, which includes the prior specification in Guha and Rodriguez (2020) as a special case.
Earlier literature on network or graph classification has been substantially motivated by the problem of classification of chemical compounds Srinivasan et al. (1996), Helma et al. (2001), where a graph represents a compound's molecular structure. In such analyses, certain discriminative patterns in a graph are identified and used as features for training a standard classification method Deshpande et al. (2005), Fei and Huan (2010). Another type of method is based on graph kernels Vishwanathan et al. (2010), which defines a similarity measure between two networks. Both of these approaches are computationally feasible only for small networks, do not account for uncertainty, and do not facilitate influential network node identification. When the number of network nodes is moderately large, a common approach to network classification is to use a few summary measures (average degree, clustering coefficient, or average path length) from the network and then apply statistical procedures in the context of standard classification methods (see, for e.g., Bullmore and Sporns (2009) and references therein). These procedures have been recently employed in exploring the relationship between the brain network and neuropsychiatric diseases, such as Parkinson's Olde Dubbelink et al. (2013) and Alzheimer's Daianu et al. (2013), but the analyses are sensitive to the chosen network topological measures, with substantially different results obtained for different types of summary statistics. Indeed, global summary statistics collapse all local network information, which can affect the accuracy of classification. Furthermore, identification of the impact of specific nodes on the response, which is of clear interest in our setting, is not feasible. As with network regression problems, an alternate approach proceeds to vectorize the network predictor and treat edge weights together as a long vector followed by developing a high dimensional regression model with this long vector of edge weights as predictors Richiardi et al. (2011); Craddock et al. (2009);Zhang et al. (2012).
This approach can take advantage of the recent developments in high dimensional binary regression, consisting of both penalized optimization Tibshirani (1996) and Bayesian shrinkage Park and Casella (2008); Carvalho et al. (2010); Armagan et al. (2013a) perspectives.
However, as mentioned in Guha and Rodriguez (2020), this treats the links of the network as exchangeable, ignoring the fact that coefficients involving common nodes can be expected to be correlated a priori. In a related work, Vogelstein et al. (2013) propose to look for a minimal set of nodes which best explains the difference between two groups of networks.
This requires solving a combinatorial problem. Again, Durante and Dunson (2017) propose a high dimensional Bayesian tensor factorization model for a population of networks that allows to test for local edge differences between two groups of subjects. Both of these approaches tend to focus mainly on classification and are not designed to detect important nodes and edges impacting the response.
Our goal in this article is to develop a high-dimensional Bayesian network classifier that additionally infers on influential nodes and edges impacting classification. To achieve this goal, we formulate a high dimensional logistic network regression model with the binary response regressed on the network predictor corresponding to each subject. The network predictor coefficient is assigned a prior from the class of Bayesian network global-local shrinkage priors discussed in this article. The proposed prior imparts low-rank and near sparse structures a priori on the network predictor coefficient. The low-rank structure of the coefficient is designed to address the transitivity effect on the network predictor coefficient and captures the effect of network edge coefficients on classification due to the interaction between nodes.
On the other hand, the near sparse structure accounts for the residual effect due to edges.
One important contribution of this article is a careful study of the asymptotic properties of the proposed binary network classification (BNC) framework. In particular, we focus on consistency properties for the posterior distribution of the BNC framework using a specific network global-local shrinkage prior, namely the Bayesian Network Lasso prior. Theory of posterior contraction for high dimensional regression models has gained traction lately, though the literature is less developed in shrinkage priors compared to point-mass priors.
For example, Castillo et al. (2012) and Belitser and Nurushev (2015) have established posterior concentration and variable selection properties for certain point-mass priors in the normal-means models. The latter article also establishes asymptotically nominal coverage of Bayesian credible sets. Results on posterior concentration and variable selection in high dimensional linear models are also established by Castillo et al. (2015) and Martin et al. (2017) for certain point-mass priors. In contrast, literature on posterior contraction properties for high dimensional Bayesian shrinkage priors is relatively limited. To this end, Armagan et al. (2013b) were the first to show posterior consistency in the ordinary linear regression model with shrinkage priors for low-dimensional settings under the assumption that the number of covariates does not exceed the number of observations. Using direct calculations, Van Der Pas et al. (2014) show that the posterior based on the ordinary horseshoe prior concentrates at the optimal rate for normal-mean problems. Recently, Song and Liang (2017) considers a general class of continuous shrinkage priors and obtains posterior contraction rates in ordinary high dimensional linear regression models. In the same vein, Wei and Ghosal (2017) offers analysis of posterior concentration for logistic regression models with shrinkage priors on coefficients. While Wei and Ghosal (2017) are the first to delineate a theoretical approach for ordinary high dimensional binary classification models with shrinkage priors, the study of posterior contraction properties for more structured binary network classification problems in the Bayesian paradigm has not appeared in the literature. In fact, developing the theory for Bayesian network classification with the Bayesian Network Lasso prior proposed in this article is faced with two major challenges. First, the novel Bayesian Network Lasso prior imparts a more complex prior structure (incorporating a low-rank structure in the prior mean of edge coefficients, as described in Guha and Rodriguez (2020) than that in Wei and Ghosal (2017), introducing additional theoretical challenges. Second, we aim at proving a challenging but practically desirable result of asymptotically optimal classification when the number of edges in the network predictor grows at a super-linear rate as a function of the sample size. Both of these present obstacles which we overcome in this work. The theoretical results provide insights on how the number of nodes in the network predictor, or the sparsity in the true network predictor coefficients should vary with sample size n to achieve asymptotically optimal classification. We must mention that developing a similar theory for the Bayesian Network Horseshoe prior proposed in this article faces more challenges due to complex prior structure in parameters. We plan to tackle that problem as part of future work.
Section 2 develops the model and the prior distributions. Section 3 discusses theoretical developments justifying the asymptotically desirable prediction from the proposed model. Section 4 details posterior computation. Results from various simulation experiments and a brain connectome data analysis have been presented in Sections 5 and 6 respectively. Finally, Section 7 concludes the article with a brief discussion of the proposed methodology.

Model Formulation
In the context of network classification, we propose the high dimensional logistic regression model of the binary response y i ∈ {0, 1} on the undirected network predictor A i as where Γ is a V × V symmetric network coefficient matrix whose (k, l)th element is given by γ k,l /2, with γ k,k = 0, for all k = 1, ..., V .
Model (1) can be expressed in the form of a generalized linear model. To be more specific, A i , Γ F = 1≤k<l≤V a i,k,l γ k,l , so that ψ i = µ + 1≤k<l≤V a i,k,l γ k,l and the probability mass function of y i can be written as Note that, if x i = (a i,1,2 , ..., a i,(V −1),V ) ∈ R V (V −1)/2 is the collection of all upper triangular elements of A i , and γ = (γ 1,2 , ..., γ (V −1),V ) ∈ R V (V −1)/2 is the vector of corresponding upper triangular elements of 2Γ, then (1) can be written as .
( 3) Although the binary network regression model is proposed for the logit link, it assumes natural extension for any other link function. The next section describes a class of network global-local shrinkage priors on network coefficients.

Bayesian network global-local shrinkage prior on the network predictor coefficient
In this article, we propose the network global-local shrinkage prior given by, Note that this framework a priori centers γ k,l at a low-rank decomposition and controls the spread of the prior distribution of γ k,l using a global-local shrinkage prior. The formulation includes a wide variety of network shrinkage priors by choosing different functions H 1 (·) and H 2 (·). For example, Guha and Rodriguez (2020) have investigated a particular class of such prior distributions, obtained by choosing H 1 (σ) = δ 1 (σ), where δ 1 (σ) is the Dirac-delta function that is defined as δ 1 (σ) = 1 if σ = 1, and 0 otherwise; and H 2 (s 2 k,l ) as an exponential density, referred to as the Network Lasso prior. To show the generality of (4), we additionally investigate performance of (4) in binary regression with s k,l ∼ C + (0, 1) and σ ∼ C + (0, 1).
The resulting prior is referred to as the Network Horseshoe prior. The rest of the hierarchy on λ r 's, u k 's follows as in Guha and Rodriguez (2020).

Posterior Contraction of the Binary Network Classification Model
This section establishes convergence results for (1) with γ k,l 's following the Bayesian Network Lasso shrinkage prior. From the hierarchical specification given in (4), the Bayesian Network Lasso shrinkage prior is given by γ k,l |s k,l ∼ N (u k Λu l , s 2 k,l ), s 2 k,l ∼ Exp(θ n /2). For the theoretical study, a common practice is to fix θ n as a function of n Armagan et al. (2013a).
Our theoretical investigations will also fix θ n (the exact expression is given in Condition (F) in the next subsection) with the fixed values specified later.
Here we consider an asymptotic setting in which the number of nodes in the network predictor, V n , grows with the sample size n. This paradigm attempts to capture the fact that the number of elements in A i , given by V 2 n can be substantially larger than sample size. Since model (1)  Let y n = (y 1 , ..., y n ) . Using the superscript (0) to indicate true parameters, the true data generating model is given by where Γ (0) is the true network coefficient. Let γ (0) be the vectorized upper triangular part of 2 is the vector of all γ (0) 2,k,l , k < l, and we denote the number of nonzero elements of γ For any > 0, define A n = γ : 1 as a neighborhood around the true density. Further suppose π n (·) and Π n (·) are the prior and posterior densities of γ with n observations, so that where p γ (y n ) denotes the likelihood of the ndimensional response vector y n .

Main Results
To show the posterior contraction results, we follow Wei and Ghosal (2017) and Armagan et al. (2013a), with substantial modifications required due to the nature of our proposed network lasso prior distribution. In proving the results, we make a couple of simplifications.
It is assumed that the dimension R of u k is fixed and is the same as R 0 , the dimension of u k . Consequently, effective dimensionality is not required to be estimated, and hence Λ = I is a non-random matrix. Additionally, we assume M to be non-random and M = I.
We emphasize that both these assumptions are not essential for the posterior contraction rate result to be true, and are only introduced for simplifying calculations.
The following theorem shows contraction of the posterior asymptotically under mild sufficient conditions on V n , s 0 2,n . The proof of the theorem is provided in Appendix F.
Theorem 3.1 Under assumptions (A)-(F) for the Bayesian Network Lasso prior on γ,
While the full posterior distributions for the parameters are not in closed forms, they mostly belong to the standard families. Hence drawing posterior samples using MCMC can be readily implemented. Appendix D and Appendix E describe full conditional distributions of parameters for Bayesian Network Lasso and Network Horseshoe priors on γ, respectively.
Let Ω (1) , ..., Ω (L) , Γ (1) , ..., Γ (L) and µ (1) , ..., µ (L) be the L post burn-in MCMC samples for Ω, Γ and µ respectively after suitable thinning. To classify a newly observed network M * as a member of one of the two groups, we compute S (l) = exp(µ (1) + M * ,Γ (l) ) 1+exp(µ (1) + M * ,Γ (l) ) for l = 1, ..., L. M * is classified as a member of group 'low' or 'high' if 1 L L l=1 S (l) is less than or greater than 0.5, respectively. To judge sensitivity to the choice of the cut-off, the simulation section Node k is recognized to be influential in the classification process if 1 are the L post burn-in MCMC samples of ξ k . Again, one of the goals of the proposed framework is to identify influential network edges impacting the response. We employ the algorithm described in Appendix C to identify influential edges. The algorithm takes care of multiplicity correction by controlling the false discovery rate (FDR) at 5% level. To study all competitors under various data generation schemes, we simulate the response from (1) given by where Γ 0 is a symmetric matrix with zero diagonal entries. The intercept µ 0 is fixed at 2 in all simulation scenarios. We consider two different schemes of generating the network A i , referred to as Simulation 1 and Simulation 2, respectively.
Simulation 1. In Simulation 1, the network edges (i.e., the elements of the matrix A i ) are simulated from N(0, 1). Thus, Simulation 1 assumes that the network predictor follows an Erdos-Renyi graph.
Simulation 2. In Simulation 2, the network predictor A i corresponding to the ith sample is generated from a stochastic blockmodel. Here nodes in a simulated network are organized into communities so that nodes in the same community tend to have stronger connections than nodes belonging to different communities. This simulation scenario simulates networks which closely mimic brain connectome networks Bullmore and Sporns (2009). To simulate networks with such community structures, we assign each node a community label, f k ∈ {1, 2, ..., 3}, k = 1, ..., V . The node assignments are the same for all networks in the population. Given the community labels, the (k, k )th element of A is simulated from , where m k,l = 0.5 when k = l. When k = l, i.e., the concerned edges connect nodes belonging to different clusters, we sample a fixed number of edge locations randomly and simulate the values from N (0, 1), assigning the values at the remaining locations to be 0. We set σ 2 0 = 1 and the three clusters with 8, 9 and 8 nodes respectively, in the three communities. We note that the network predictors are simulated from a stochastic blockmodel in Simulation 2 which also ensures transitivity in the network predictor.
Simulating the network predictor coefficient Γ 0 . In both Simulations 1 and 2, the network predictor coefficient Γ 0 is constructed as the sum of two matrices Γ 0,1 and Γ 0,2 . We provide the details of constructing the two matrices as below.
In both Simulations 1 and 2, we draw V latent variables u k,0 , each of dimension R g , from a mixture distribution given by where δ 0 is the Dirac-delta function and π is the probability of any u k,0 being nonzero.
Define a symmetric matrix Γ 0,1 whose (k, l)th element is given by u k,0 u l,0 2 , k < l and = 0 if k = l. Note that if u k,0 is zero, then the kth node has no contribution to the mean function in (7), i.e., the kth node becomes non-influential in predicting the response. Since (1 − π) is the probability of a node being inactive, it is referred to as the node sparsity parameter in the context of the data generation mechanism under Simulations 1 and 2. All elements of u m,g are taken to be 0.5 and u s,g is taken to be 1.
We also construct another symmetric sparse matrix Γ 0,2 to add additional edge effects corresponding to edges connecting a few randomly selected nodes. Let π 2 be the proportion of nonzero elements of Γ 0,2 , set randomly at either 0.05 or 0.1. We randomly choose π 2 proportion of locations from the set of all (k, l). The nonzero entries are drawn using one of the three following strategies: Strategy 1: Nonzero entries are simulated from N(1, 0.1).
Strategy 3: All nonzero entries are fixed at 0.5.
The quantity (1 − π 2 ) is referred to as the residual edge sparsity.
For a comprehensive picture of Simulation 1 and Simulation 2, we consider 4 different cases each in both simulations as summarized in Table 1 and 2 respectively. In each of these cases, the network predictor coefficient and the response are generated by changing the node sparsity (1 − π), the residual edge sparsity (1 − π 2 ) and the true dimension R g of the latent variables u k,0 's. The table also presents the maximum fitted dimension R of the latent variables u k for the logistic regression model (2). Note that the various cases also allow model mis-specification with unequal choices of R and R g .
As competitors, we use generic variable selection and shrinkage methods that treat edges between nodes together as a long predictor vector to run high dimensional regression, thereby ignoring the relational nature of the predictor. More specifically, we use Lasso Tibshirani (1996), which is a popular penalized optimization scheme, and the Bayesian Lasso (BLasso for short) Park and Casella (2008) and Bayesian Horseshoe (BHS for short) priors Carvalho Cases R g R Node Residual Edge Strategy Sparsity (1 − π) Sparsity (1 − π 2 ) Case -1 2 2 0.5 0.95 Strategy 1 Case -2 3 5 0.6 0.95 Strategy 1 Case -3 2 5 0.5 0.90 Strategy 2 Case -4 2 5 0.4 0.90 Strategy 3 Table 1: Table presents different cases for Simulation 1 . The true dimension R g is the dimension of vector object u k,0 using which data has been generated. The maximum dimension R is the dimension of vector object u k using which the model has been fitted. Node sparsity and residual edge sparsity are described in the text.
Cases  The true dimension R g is the dimension of vector object u k,0 using which data has been generated. The maximum dimension R is the dimension of vector object u k using which the model has been fitted. Node sparsity and residual edge sparsity are described in the text. to implement the frequentist Lasso, while we write our own codes for BLasso and BHS.
A comparison with these methods will indicate any relative advantage of exploiting the structure of the network predictor. Additionally, we compare our methods to a frequentist approach that develops network classification in the presence of a network predictor and a binary response Relión et al. (2017). We refer to this approach as Relión.
All Bayesian competitors are allowed to draw 50, 000 MCMC samples, out of which the first 30, 000 are discarded as burn-ins. Convergence is assessed by comparing different in Park and Casella (2008). Sensitivity to the choice of hyper-parameters is discussed later, both for simulation studies and for the real data analysis.

Identification of Influential Nodes
Figures 1 and 2 show the posterior probability of the k-th node being detected as influential, i.e., P (ξ k = 1|Data), by BNLC and BNHC for each node and each case within Simulations 1 and 2, respectively. Some interesting observations emerge from the results. We find that both methods work well with lower node sparsity and higher residual edge sparsity.
Decreasing the residual edge sparsity and increasing the node sparsity have adverse effects on the performance. In general, BNLC shows relatively better performance than BNHC in cases with higher node sparsity and/or lower residual edge sparsity. We provide a brief discussion below to support these observations.
For BNHC, case 2 exhibits a few false positives, and the separation of posterior probabilities for truly active and truly inactive nodes is much more stark in case 1 than in case 2.
BNLC does a better job of node identification than BNHC in case 2. Residual edge effect does have an impact on the probabilities, which is evident by comparing cases 1 and 3. For BNHC, case 3 (Simulation 1) displays poor performance with a higher number of both false positives and false negatives. Performance of BNLC appears to be better than BNHC in case 3. Fixing the residual edge sparsity and increasing the node sparsity has a negative impact on node identification, as seen by comparing performances in cases 3 and 4 (Simulation 1). For Simulation 2, both competitors perform quite well in cases 1 a nd 2. Again, case 3 (Simulation 2) represents a higher node sparsity, so that both BNHC and BNLC do not perform well in this case. Similar to Simulation 1, BNHC shows inferior performance to BNLC in case 3. While BNHC offers a few false positives and false negatives in case 4 (Simulation 2), the performance appears to be much better than in case 3. Notice that case 3 has both higher node sparsity and residual edge sparsity than case 4. While they have opposing effects, it appears that higher node sparsity demonstrates more of an adverse effect here compared to a small perturbation in the residual edge sparsity. Recall that Relión et al.
(2017) is the only other competitor which is designed to detect influential nodes. It detects all nodes to be influential in all simulation cases.

Identification of Influential Edges
We apply the algorithm with a mixture of skewed t-distributions described in Appendix C to detect influential edges from the post burn-in MCMC samples of the edge coefficients using a threshold of t = 0.05. The proposed approach controls FDR below a threshold of 0.05 to account for multiplicity correction. Tables 3 and 4 provide the true positive rates (TPR) and false positive rates (FPR) in detecting important edges for Simulations 1 and 2 for the competitors, respectively. It is observed that when node sparsity is moderate and residual edge sparsity is high (cases 1 and 2), both BNLC and BNHC offer moderate performance in terms of identifying true positives, and include very few false positives. In these cases, BNHC generally exhibits a little higher FPR than BNLC. In the case of high node sparsity (e.g., case 3, Simulation 2) both these methods unfortunately show much lower true positive rates.
Again, lower edge sparsity (case 3, Simulation 1) has almost no effect on FPR of BNLC,   The results in Tables 3 and 4 indicate higher number of edges identified as influential

Estimation of Edge Coefficients and Classification Accuracy
The

Sensitivity to the choice of Hyperparameters
To assess how sensitive the inferences from BNLC and BNHC are, we analyze BNLC and BNHC with different combinations of hyperparameters. Specifically for BNLC, we use the       Thus, we classify the subjects in our study as belonging to the low IQ group if y = 0, and the high IQ group if y = 1.
Along with FSIQ measurements, brain connectome information for n = 114 subjects is gathered using weighted diffusion tensor imaging (DTI). DTI is a brain imaging technique that enables measurement of the restricted diffusion of water in tissue in order to produce neural tract images. The brain imaging data we use has been pre-processed using the NDMG pre-processing pipeline Kiar et al. (2016); Kiar et al. (2017a); Kiar et al. (2017b). In the context of DTI, the human brain is divided according to the Desikan atlas Desikan et al. (2006), which identifies 34 cortical regions of interest (ROIs) both in the left and right hemispheres of the human brain, implying 68 cortical ROIs in all. Similar to Guha and Rodriguez (2020), this results in a brain network of a 68 × 68 matrix for each individual.
Our scientific goals in this setting include identification of brain regions or network nodes significantly related to FSIQ and classification of a subject into the low IQ or high IQ group based on his/her brain connectome information.
Identical prior distributions for all the parameters as in the simulation studies have been used. BNLC and BNHC are both fitted with R = 4, which is found to be sufficient for this study. Further, Guha and Rodriguez (2020) show robust inference as long as the chosen R is bigger than the effective dimensionality of the latent variables. Similar to article Guha and Rodriguez (2020) All inference is based on the remaining 20, 000 post burn-in iterates appropriately thinned.

Findings from the Brain Connectome Application
As in simulation studies, we put our emphasis on identifying influential The ROIs are mainly detected in the temporal, frontal, parietal and cingulate lobes in both hemispheres. Figure 8 plots the estimated posterior probability of an ROI being detected as influential by the BNLC and BNHC models. Notably, there are 29 ROIs identified by both BNLC and BNHC, given in Table 10.  Stuss et al. (1985); Razumnikova (2007); Miller and Milner (1985); Kolb and Milner (1981). Our method also finds a significant association between FSIQ and the left inferior parietal lobule, the left precuneus and the supramarginal gyri in both the hemispheres, in the parietal lobe, regions also found to be significantly related to FSIQ by Yoon et al. (2017).
We additionally look into ROIs which are detected by only of the two methods (lets say, BNLC), and report the posterior probabilities of these ROIs being active under the other method (i.e., BNHC). Figure 7 shows the posterior probabilities of nodes being active under the 'other' method as discussed above. It is observed that the nodes selected by BNHC but not by BNLC have probabilities not very far from 0.5 under BNLC, which says that BNLC is not enough confident to exclude these nodes from the set of influential nodes. However, most of the nodes selected by BNLC but not by BNHC show smaller probabilities of being influential under BNHC. Perhaps, BNLC is more conservative in including nodes in the set of influential nodes, which is responsible for the discrepancy between the number of identified nodes by BNHC and BNLC.
As described earlier, we identify influential edges connecting pairs of influential nodes using the algorithm described in Appendix C. Figure 9 presents the influential edges (among all edges connecting pairs of influential nodes) identified by the BNLC and BNHC models.
Note that BNLC and BNHC identify 142 and 291 edges as being influential out of 38 2 and 48 2 possibilities, respectively. Since a different number of nodes are detected as influential by BNHC and BNLC, to make a fair comparison, we consider the 29 nodes detected as influential by both methods, and use our algorithm to find the number of influential edges among these 29 2 possibilities for both BNLC and BNHC. The numbers turn out to be 96 and 184, respectively. We note that there are a few nodes which are identified as influential by either BNHC or BNLC, but none of the edges connecting these nodes are found to be influential. As an example, although the frontal pole and the temporal pole in the left hemisphere are identified as influential nodes by BNLC, none of the edges connecting these two nodes turn out to be influential. This phenomenon may be due to the use of the FDR in the edge selection procedure, which finds edges that are most likely to be active while controlling for false discoveries. Hence, not identifying an edge does not necessarily mean that the edge is not active, it just means that there are others that satisfy the criteria better. Similar to simulation studies, we dig deeper to analyze the discrepancy in the number of influential edges identified by BNLC and BNHC. Specifically, we rank the 29 2 = 406 edges connecting the nodes found to be influential by both BNLC and BNHC, according to the absolute values of their posterior means. Table 11 shows between 23-74% intersections.

Nodes Selected by BNLC, but not by BNHC
To examine the predictive ability of the Bayesian network classification model, we report the area under curve (AUC) of the ROC curve for BNLC and BNHC, along with all competing methods. The AUCs are computed using a 10-fold cross validation approach. The AUC estimates presented in Table 12 indicate better performance of both BNLC and BNHC, with BNLC slightly outperforming. Frequentist Binary Lasso turns out to be the next best performer, while BLasso and BHS perform very similar to a random classifier. Finally, the effective dimensionality of the model is investigated for both BNLC and BNHC, and they turn out to be 2.17 and 2, respectively.     Table 12: Predictive performance of Bayesian Network Classification (BNC) vis-a-vis competitors in terms of Area Under Curve (AUC) of the ROC. AUC has been calculated in each case using 10-fold cross validation.

Sensitivity to the choice of hyperparameters
We have already discussed how the hyperparameters are chosen for the simulation studies and data analysis. To assess how sensitive the inferences from BNLC and BNHC are, we analyze BNLC and BNHC with different combinations of hyperparameters. Specifically for BNLC, we use the five different combinations (i)-(v) given in Section 5.5, and three different combinations (i)'-(iii)' for BNHC also mentioned in Section 5.5. We report performances on the number of influential nodes identified. We also find the number of influential edges connecting influential nodes.
rh-transversetemporal rh-frontalpole rh-supramarginal rh-superiortemporal rh-superiorparietal rh-rostralmiddlefrontal rh-rostralanteriorcingulate rh-precuneus rh-precentral rh-postcentral rh-pericalcarine rh-parstriangularis rh-parsorbitalis rh-parsopercularis rh-paracentral rh-parahippocampal rh-middletemporal rh-lingual rh-lateraloccipital rh-inferiortemporal rh-entorhinal rh-caudalmiddlefrontal rh-caudalanteriorcingulate lh-transversetemporal lh-frontalpole lh-supramarginal lh-superiortemporal lh-superiorfrontal lh-rostralmiddlefrontal lh-rostralanteriorcingulate lh-posteriorcingulate lh-postcentral lh-parstriangularis lh-parsorbitalis lh-parsopercularis lh-paracentral lh-middletemporal lh-lingual lh-lateralorbitofrontal lh-lateraloccipital lh-isthmuscingulate lh-inferiortemporal lh-inferiorparietal lh-fusiform lh-entorhinal lh-cuneus lh-corpuscallosum lh-bankssts (b) BNHC Figure 9: Plot showing whether an edge connecting two influential nodes is influential or not. Note that the map is a M × M symmetric matrix, where M denotes the number of influential nodes, and each cell denotes an edge connecting the corresponding pair of nodes. The axis labels are the abbreviated names of the influential ROIs in the left (starting with 'lh -') and the right (starting with 'rh -') hemispheres of the brain. Full names of the ROIs can be obtained from the widely available Desikan brain atlas. A white cell represents an influential edge, while red cell represents a non-influential edge.     Table 14 presents the number of edges detected as influential, as well as the number of intersecting edges with the original analysis. Again, due to the high dimensionality of the problem, the variation in the number of identified edges with different choices of hyperparameters is expected, though the variation turns out not to be very significant.
Finally, to check sensitivity to the choice of R on the performance of BNLC and BNHC, we run the data analysis for BNHC and BNLC with R = 8 and R = 10, and report the posterior mean of the effective dimensionality, along with AUC. Table 15 reports the posterior mean of effective dimensionality, which shows very moderate increase with increasing R. However, increasing R seems to have almost no effect on AUC.

Conclusion
We develop a binary Bayesian network regression model that enables classifying multiple networks with "labeled nodes" into two groups, identifies influential network nodes and predicts the class in which a newly observed network belongs. Our contribution lies in carefully constructing a class of network global-local shrinkage priors on the network predictor coefficient while recognizing the latent network structure in the predictor variable. In particular, we investigate two specific network shrinkage priors from this general class, leading to two network classifiers BNLC and BNHC. Our extensive simulation study shows competitive performance p(ω i ) ∼ P G(1, 0), λ r ∼ Ber(π r ), π r ∼ Beta(1, r η ), η > 1.
The full conditional distributions of the model parameters are given below.
where GIG denotes the generalized inverse Gaussian distribution.

Appendix C
Similar to the assumptions made by Wei and Ghosal (2017) in their proof of posterior consistency for binary logistic regression, we prove our results assuming that the centering parameter µ = 0 in both the true and the data generating models. We note that the main structure of the proof will remain unchanged with this assumption and the result proved in this article can be trivially extended to the setting with nonzero µ.
We begin by defining some notations. In the proof, Π(·) will be used to denote the generic probability notation. We define the notation of the log-likelihood function by We also introduce the function C y n ,n (·) to quantify the curvature of w γ,n (y n ) around γ (0) , where ∇w γ (0) ,n (y n ) is the derivative of w γ (0) ,n (y n ) w.r.t. γ, evaluated at γ (0) . Also the likelihood p γ (y n ) can be written using the above notations as p γ (y n ) = n i=1 exp(w γ,n (y i )). The notations E γ (·) and E γ (0) (·) have been reserved to denote expectation w.r.t the distribution of y n |γ and y n |γ (0) respectively. The proof of Theorem 3.1 relies in part on the existence of exponentially consistent sequence of tests.
Proof The construction of the test is provided in the proof of Theorem 2 and Lemma 4 in Ghosal et al. (2006).
We also state another result which will be subsequently used in the proof.

Proof of Theorem 3.1
Suppose E n = y : ||∇w γ (0) ,n (y)|| ∞ ≤ 2 √ nq n . Then the probability of the vector y n belonging to the set E n is given by, where the last step follows from the Hoeffding inequality. Note that as n → ∞, q n → ∞, hence P γ (0) (y n ∈ E n ) → 1. Hence, in the subsequent proof we can assume without loss of generality that y n ∈ E n . It can be observed that Π n (A c n ) = A c n p γ (y n )π n (γ) p γ (y n )π n (γ) = A c n pγ (y n ) p γ (0) (y n ) π n (γ) pγ (y n ) p γ (0) (y n ) π n (γ) where Φ n is the exponentially consistent sequence of tests given in Lemma 8.1. The above equation is true as N n /D n ≤ 1. This is in turn true as both are integrals of the same nonnegative functions, D n is the integral of that function over the entire set of possible γ's, while N n is the integral over a subset A c n . In proving Theorem 3.1, we will proceed in three steps as following. (a) Step 1 shows that Φ n → 0, as n → ∞, almost surely. (b) Step 2 shows that exp(hn/2)(1 − Φ n )N n → 0, as n → ∞, almost surely.