Unifying inference on brain network variations in neurological diseases: The Alzheimer's case

There is growing interest in understanding how the structural interconnections among brain regions change with the occurrence of neurological diseases. Diffusion weighted MRI imaging has allowed researchers to non-invasively estimate a network of structural cortical connections made by white matter tracts, but current statistical methods for relating such networks to the presence or absence of a disease cannot exploit this rich network information. Standard practice considers each edge independently or summarizes the network with a few simple features. We enable dramatic gains in biological insight via a novel unifying methodology for inference on brain network variations associated to the occurrence of neurological diseases. The key of this approach is to define a probabilistic generative mechanism directly on the space of network configurations via dependent mixtures of low-rank factorizations, which efficiently exploit network information and allow the probability mass function for the brain network-valued random variable to vary flexibly across the group of patients characterized by a specific neurological disease and the one comprising age-matched cognitively healthy individuals.


Introduction
There is fundamental interest in understanding the relationship between brain connectivity structure and neurodegenerative disorders, such as Alzheimer's [1], Parkinson's [2] and other dementias [3]. These diseases are mainly found in aged populations and affect the normal functions of the central and peripheral nervous system causing, among others, muscle weakness, loss of coordination and cognitive impairment.
Alarming prevalence projections of dementia cases by the World Health Organization in 2006 [4], and the rapid devel-opment of brain imaging technologies in recent years, have stimulated intensive research aimed at understanding how the brain structure is compromised with specific neurological diseases. This is key to improving diagnosis as well as providing increasingly targeted therapies.
Most of the related literature has focused on the modular paradigm, which considers brain regions as specialized actors in specific cognitive functions [5]. Under this paradigm, inference focuses on multivariate data ψ i = (ψ i1 , . . . , ψ iV ) T , where ψ iv is the activity level in region v (v = 1, . . . ,V ) for individual i (i = 1, . . . , n), often measured via functional magnetic resonance imaging (fMRI). Let y i denote an indicator of whether individual i is in the control group (y i = 0) or has a specific disorder of interest (y i = 1). Inference then proceeds by testing for global and local variations in the brain activity vector across these two groups.
In contrast to functional correlations, diffusion weighted magnetic resonance imaging (dw-MRI) can provide information on the structural connection network in the brain by approximating the diffusion of water molecules across tissues. Water diffuses equally in all directions in grey matter, but is constrained and does not similarly diffuse across white matter fibers. This allows for the reconstruction of white matter tracts using tractography methods for dw-MRI [6] and mapping the patterns of connections formed from the white matter fiber bundles. Using this information, we can obtain measurements Adjacency matrix A i representing the brain network of the i-th subject. Black corresponds to edges, white to nonedges. In our Alzheimer's study, brain anatomical regions are defined by the Desikan atlas [9] for a total of V = 68 nodes equally divided between left and right hemisphere.
of the brain's network of cortical connectivity by seeing if certain anatomical regions are connected by white matter.
In this article, we focus on brain network data corresponding to a V × V symmetric binary adjacency matrix A i for individual i having elements A i[vu] = A i [uv] = 1 if there is one or more white matter tracts connecting brain regions v (v = 2, . . . ,V ) and u (u = 1, . . . , v − 1) in individual i (i = 1, . . . , n) and A i[vu] = A i [uv] = 0 otherwise. Refer to Fig. 1 for an illustrative example. Such connectome data provide substantial insights into the sources of complex cognitive processes relative to brain activity data analyzed in the modular paradigm [7]. Cognitive network theory views the studying of structural connectivity networks as a key to learning the complexity of information processing mechanisms and understanding the effect of neurological disorders on these cognitive processes [8].
Much of today's literature is focused on analytic methods for understanding localized brain activity data ψ i , yet methodologies for analyzing brain network data A i is still in its infancy. Our main aim is to develop techniques to assess whether and how a network-valued random variable generating structural brain networks A i (i = 1, . . . , n) varies across diagnostic groups. In particular, it is of interest to test for global variation in the overall brain network structure across groups, while identifying specific local variations to understand if and which brain connections are compromised by a specific neurological disease of interest.
There has been an increasing attention in the literature toward methods for addressing these aims; see [10] and the references cited therein for an overview. Common prac-tice proceeds by reducing A i to summary measures θ i = (θ i1 , . . . , θ ip ) T (i = 1, . . . , n) and then applying standard multivariate analyses such as MANOVA (see e.g. [11]) to assess whether these network measures change with the occurrence of a disease. Summary statistics are commonly chosen to represent global network characteristics, such as the average path length, network density, transitivity and k-core [12]. These measures provide a useful simplification of a complex problem, but cannot characterize the entire network structure and hence may fail to detect important relationships between brain network and neurological disorders, leading to inconsistent results in the literature; see [13] for a review of inconsistencies when relating brain networks to creative reasoning.
An alternative class of methods avoid reducing brain networks to summary statistics by performing separate and independent tests to assess which brain connections are compromised by a neurological disease. As there are V (V − 1)/2 edges in the V brain regions under study, with V = 68 using the Desikan atlas [9], the number of tests is substantial and requires adjustments of the significance threshold to control for multiplicity [14]. Such mass-univariate approaches do not exploit network information, leading to low power [15], and substantially underestimating the number of connections compromised by a specific neurological disorder. Recent proposals try to gain power by replacing the false discovery rate (FDR) control [16], with thresholding procedures that account for the network structure [17]. However, such approaches require careful interpretation, while being highly computationally intensive and complex.
We propose a fundamentally new approach based on defining a generative probabilistic model for the brain network data. In particular, the probability mass function (pmf) for the network-valued random variable is assigned a mixture model, allocating individuals to subpopulations in terms of their brain network structure. Within a subpopulation, the edge probabilities are related to a latent similarity measure via a logistic mapping. The similarity matrix is then factorized as the sum of a common component and a subpopulationspecific one that arises from embedding the brain regions in a low dimensional latent space that accounts for the network structure. Using this flexible mixture model as scaffolding, we develop an efficient testing method by allowing the mixture weights to vary across case and control groups. This induces a highly efficient Bayesian testing procedure that adjusts automatically for multiple comparisons in drawing inferences on brain structural differences across groups.
In the rest of the paper, we describe the model formulation, providing insights on theoretical properties, estimation procedures and inference techniques. An application to assess variations in the brain architecture with Alzheimer's disease illustrates the benefits of the proposed procedure, while providing novel insights into the impact of the disease on brain structural connections.

Dependent Mixture of Low-rank Factorizations
A network-valued random variable A generating symmetric V × V binary adjacency matrices can be uniquely characterized by its lower triangular matrix elements L (A ) = (A 21 , A 31 , . . . , This number becomes intractable and massively larger than the sample size n even in using coarse brain regions. In the motivating Alzheimer's disease study, there are V = 68 brain regions. This implies that, in the absence of constraints, we need to estimtae 2 68(68−1)/2 − 1 = 2 2278 − 1 free parameters characterizing p L (A ) . Clearly no studies will ever have this many subjects, and hence it is necessary to reduce dimensionality and make the problem tractable. However, in reducing dimension, it is important to maintain flexibility in characterizing the structure underlying brain networks.
We propose to reduce dimensionality, while maintaining flexibility in characterizing p L (A ) , by using a hierarchical latent space representation. The idea is to assign each brain region a coordinate in a lower dimensional Euclidean space; such models have been used effectively in social network contexts [18]. Our contribution is the generalization to a latent space random effects model which defines the population distribution of network-valued data, while characterizing individual differences in the architecture of interconnections in the brain. The random effects are modeled as a mixture of lowrank factorizations. This induces clustering of individuals in terms of their brain structure, and facilitates inferences on differences in case and controls by allowing the mixing weights to vary across these groups. We first describe the low-rank factorization structure, which represents the key building block to reduce dimensionality and borrow network information.

Low-rank factorization
the vector of probabilities of edges between each pair of brain regions, our probabilistic low-rank factorization generates brain networks in two main steps displayed in Fig. 2 and Fig.  3. The first constructs the edge probability vector π exploiting network information as follows: 1. Embed each region v in an R-dimensional Euclidean space (in Fig. 2, R = 2), with X a V × R matrix of latent coordinates; X vr is the rth coordinate for brain region v.
2. The rth column of X corresponds to the coordinates for the different brain regions in dimension r. These latent Weighted Latent Coordinates  coordinates are sampled from a standard normal distribution, and then multiplied by √ λ r , with λ r measuring the overall importance of the rth dimension of the latent space on brain network structure. In Fig. 2 λ 1 = 1 while λ 2 = 0, meaning that the second coordinate doesn't play any role in generating the network.

Weighted Latent Coordinates
3. Focusing on the lth pair of brain regions, corresponding to nodes v and u, v > u, construct the similarity measure S l via dot product of their weighted latent coordinate vectors S l = ∑ R r=1 λ r X vr X ur , for each l = 1, . . . ,V (V − 1)/2. 4. Define π l by mapping each similarity measure S l ∈ ℜ into the probability space via the one-to-one continuous increasing logistic mapping, so that The second step generates networks L (A i ) by sampling their edges L (A i ) l (l = 1, . . . ,V (V − 1)/2) from conditionally independent Bernoulli random variables given connection probabilities π l = Pr{L (A i ) l = 1} ∈ (0, 1). Although the mechanism generates networks with conditionally independent edges given π, the shared dependence on a common set of node-specific latent coordinates induced by the dot product representation of S = (S 1 , . . . , S V (V −1)/2 ) T ∈ ℜ V (V −1)/2 can define arbitrarily rich dependence structures. This is illustrated in Fig. 3, where the networks preserve a two-block structure induced by factorization in Fig. 2. The low-rank factorization allows dimensionality reduction from V (V − 1)/2 edge probabilities to V × R latent coordinates and R weights. The formulation facilitates adaptive collapsing on lower dimensional models by shrinking the weights λ r towards 0 as r increases.

Mixture of low-rank factorizations
Under the proposed model, the probability assigned to config- where Bern(a l ; π l ) = π a l l (1 − π l ) 1−a l and π is generated according to Fig. 2. Although providing a key building block to reduce dimensionality and include network information, the low-rank factorization is insufficiently flexible in assuming independence across pairs of brain regions in the occurrence of connections conditionally on π. To improve flexibility, we generalize the model to mix together H subpopulations.
Let G i ∈ {1, . . . , H} index the subpopulation that the ith brain network A i is generated from, with p G (h) = Pr(G = h) = ν h the probability of being drawn from that subpopulation, and π (h) the connections probability vector specific to that subpopulation. Then, we have the following generative process: 1. Allocate the ith individual to a subpopulation by sampling G i according to p G .
2. Given G i = h and the corresponding π (h) , generate A i by sampling its edges L (A i ) l (l = 1, . . . ,V (V − 1)/2) from conditionally independent Bernoulli random variables given the connection probabilities specific to subpopulation h, π . By marginalizing out the subpopulation indicator G i , this hierarchical probabilistic generative process leads to a mixture representation for p L (A ) of the form Equation (2) is much more flexible than (1). Statistical properties of representation (2) have been studied in [19], showing the full flexibility of the mixture of low-rank factorizations in representing every possible pmf p L (A ) . They additionally slightly modify the steps generating the class-specific edge probability vectors π (h) = {1 + exp(−S (h) )} −1 (h = 1, . . . , H) to allow inference on shared versus class-specific components of variability in the brain connectivity structure. In particular, instead of defining the entries in the latent similarity vector for the hth subpopulation for each l = 1, . . . ,V (V − 1)/2, with Z l ∈ ℜ a latent similarity measure for the lth pair of regions common to all individuals and D (h) l ∈ ℜ a deviation specific to subpopulation h. The common similarity Z l in (3) is left unstructured as it can be estimated borrowing information across all individuals, while D (h) l is constructed according to steps 1 − 3 of the low-rank factorization mechanism for each h = 1, . . . , H exploiting the network structure to cope with less information in the data about class-specific deviations.

Inferences on differences across groups
When the focus is on inference on changes in brain networks with the occurrence of a disease, the model described above needs to be generalized. Although we treat Y as a random variable through a prospective likelihood, the method we propose is valid for inference on differences across groups in brain network structure also for case control studies that sample under a retrospective design.
The conditional pmf for the brain network given the presence or absence of the disease, p L (A )|y , y ∈ {0, 1}, is characterized via a slight modification of the mixture of lowrank network factorization proposed in the previous section. In particular, we simply allow the proportions of individuals in each subpopulation to vary across groups by letting p G|y (h) = Pr(G = h | Y = y) = ν hy , while holding the subpopulation specific edge probability vectors π (h) fixed. A graphical representation of the final model is provided in Fig.  4. Replacing Pr(G = h) = ν h with Pr(G = h | Y = y) = ν hy (h = 1, . . . , H) in equation (2), leads to the dependent mixture representation for control (y = 0) and case (y = 1) group. According to theoretical properties in [20], factorizations (4) − (5), with each π (h) constructed as in equation (3) are fully flexible in characterizing any joint pmf p Y ,L (A ) , while reducing the dimensionality in requiring estimation of H low-rank factorization mechanisms, rather than 2 V (V −1)/2 − 1 configuration probabilities for each group y.
Under this formulation, testing of the null hypothesis of no differences in the distribution of brain networks between the case and control groups versus the alternative of some difference can be mathematically expressed as

Local inferences on whether individual edge probabilities vary with the presence of a neurological disorder can instead be based on Cramer's V [21]
for l = 1, . . . ,V (V − 1)/2, where p Y ,L (A ) l (y, a l ) = Pr{Y = y, L (A ) l = a l } and p L (A ) l (a l ) = Pr{L (A ) l = a l }. The coefficient ρ l has values in the interval [0, 1], with ρ l = 0 meaning that p Y ,L (A ) l = p Y p L (A ) l , so that the probability of a connection between the lth pair of brain regions is the same for cases and controls.

Bayesian inference
Inference for the model in Fig. 4 is performed under a Bayesian paradigm. Bayesian inference characterizes uncertainty and facilitates adaptive choice of R and H through carefully specified priors. Hence, the focus is on the posterior distribution As the joint p Y ,L (A ) is defined via equations (3) − (5), the posterior distribution in equation (8) can be easily obtained as a function of the posteriors for parameters in (3) − (5). Inference proceeds by updating the independent priors for the quan- . . , H and ν y = (ν 1y , . . . , ν Hy ) ∼ Π ν , y ∈ {0, 1} given brain network data and case-control status, to obtain posterior samples via MCMC methods. We compute (8) as a function of these posteriors via (3) − (5).
Prior distributions are carefully defined to induce a prior Π on the joint pmf p Y ,L (A ) with simple posterior computation procedures, allowance for testing and flexibility. These aims are accomplished by letting p Y (0) = 1− p Y (1) ∼ Beta(a 0 , a 1 ), while choosing Normal priors for the entries in Z and standard Gaussians for the elements in X (h) . To learn the dimensions of the latent spaces and penalize high dimensional representations, a multiplicative inverse gamma is defined for λ (h) . This choice for Π λ favors shrinkage effects with elements in λ (h) stochastically decreasing towards 0 as r increases; see [19] for details.
Prior Π ν is instead defined to incorporate global hypothesis testing in equation (6), while allowing automatic deletion of redundant classes. This is accomplished by introducing an hypothesis indicator T ∈ {0, 1}, with T = 0 for H 0 and T = 1 for H 1 . If T = 1 the latent class assignment mechanism is different between case and controls and two independent Dirichlet priors priors are assigned to ν 0 and ν 1 , respectively. If T = 0, we set ν 0 = ν 1 = u ∼ Dir(1/H, . . . , 1/H) as p G is the same in the two groups, under H 0 . Hence in assessing evidence in favor of the alternative, one can rely on the posterior probability, Pr{H 1 | y 1 , L (A 1 ), . . . , y n , L (A n )} = Pr{T = 1 | y 1 , L (A 1 ), . . . , y n , L (A n )}. This quantity can be easily computed in the MCMC algorithm for posterior computation; see [20] for details. Beside this benefit for global hypothesis testing, deletion of redundant classes can also be easily accomplished by choosing small values for the hyperparameters in the Dirichlet priors [22].
Simulation studies in [20] highlight the good performance of the proposed methodology in accurately estimating the joint pmf p Y ,L (A ) and in accurately identifying differences between case and controls in the brain network structure, while identifying local variations in each edge probability. Across multiple scenarios, the proposed local testing procedure has a Type I error of 0.00044, Type II of 0.0587 and FDR of 0.0023. Independent screening via separate two-sided Fisher's exact tests (see e.g. [24]) with FDR control has Type I error of 0.0036, Type II of 0.5983 and FDR of 0.0387. These improvements are also exhibited in global testing with our procedure having both Type I and II errors of 0.01. In contrast MANOVA test on network features has Type I error of 0.10 and Type II error of 0.87.

Application to Alzheimer's Disease
According to the Centers for Disease Control and Prevention (CDC), Alzheimer's disease (AD) is the most common form of dementia and the sixth leading cause of death in the United States. Unlike cancer and heart disease death rates, which are expected to decline, the growth of elderly population in the age range most commonly affected by dementia is leading to an increase of the death rates due to AD [25]. This has strongly motivated intensive research aimed at finding the sources of AD in the human brain to develop increasingly refined diagnosis and prognosis procedures and improve therapy.
Current understanding of variations in brain behavior across AD is mostly available via early neuropathological studies (e.g. [26]), and contributions analyzing joint or local changes in the activity of each region under the modular paradigm (e.g. [27]). More recent proposals shift increasingly away from the above approach towards studying brain activity networks via changes of the covariance in activity across brain regions for AD and controls (e.g. [28]). However, functional connectivity matrices estimated from fMRI data do not reflect the underlying axonal pathways that can give rise to changes in function, and often require caution in interpreting the results [8]. This has motivated an increasing interest in structural connectivity matrices estimated from diffusion scans. Early studies on these data proceed by assessing variations of global brain network measures or region-specific connectivity statistics across AD and controls (e.g. [29]). As previously noted, these methods may fail in flexibly characterizing the richness of the brain network structure, leading to inconsistent results.
To address these issues, we apply our methodology to brain networks derived from the Alzheimer's Disease Neuroimaging Initiative (ADNI), providing consistent new insights which contribute to solving the ongoing mystery behind the mechanisms of AD in the human brain. Data (y i , A i ) are available for n = 92 individuals, with 50 in the control group (y i = 0) and 42 age-matched patients with AD (y i = 1). Each adjacency matrix A i represents the brain network of the i-th individual, measuring the structural connectivity among V = 68 lateralized brain regions defined by the Desikan atlas [9] as part of FreeSurfer [30]. In particular A i has elements A i[vu] = A i [uv] = 1 if at least one white matter tract connects brain regions v and u in individual i and A i[vu] = A i[uv] = 0, otherwise. These structural networks, also known as connectomes, represent estimates of the axonal-fiber pathways connecting the different regions. Connectomes considered in this application have been estimated as in [29] via recently developed pipelines, which efficiently exploit structural MRI data to obtain a parcellation of the brain in anatomical regions, and dw-MRI images to recover the fiber streamlines connecting each pair of brain regions; see [29] for details. Posterior computation and inference is performed considering the same settings as in the application to creativity in [20].

Changes in Brain Network with Alzheimer's
The global testing procedure in (6) strongly favors the hypothesis of association between brain structural connectivity and AD diagnosis withPr{H 1 | y 1 , L (A 1 ), . . . , y n , L (A n )} > 0.99. This confirms findings in [29] highlighting significant variations in brain network summary measures when comparing AD patients with cognitively healthy controls. As expected the estimated significant differences between the edge probabilities in AD group and control group in Fig.  5 show an overall less connected brain network for the AD group compared to controls, in line with [29] and literature on AD. The main differences appear in terms of intra-hemispheric connections in the left hemisphere, while fewer local differences are found also in terms of inter-hemispheric connections and right intra-hemispheric. This major role of the left hemisphere agrees with [29,27].
The agreement with previous studies highlights the consistency of our methodology, which has the additional benefit of providing inference not only on the scale of the network summary measures but in terms of variations of the entire pmf for the brain network-valued random variable representing brain interconnections. This rules out the issue of conflicting conclusions when different network statistics are considered, while also avoiding ad-hoc choices when defining certain summary measures. Recalling for example [29] one may obtain different results when considering an order for the k-core different from 18. An additional benefit of our approach, as outlined in the simulation study, is that local testing automatically controls for multiplicity, while out-performing frequentist competitors controlling for FDR, in terms of power. Recalling the application to AD, this leads to a procedure which can more easily identify connections significantly varying between control and AD subjects. This is evident when comparing Fig. 5 to results in Fig. 1 in [29] learning less significant local differences. This result may be related to the Figure 5: For connections significantly varying in AD according to (7), weighted network visualization with weights given by the estimated Pr{L . . ,V (V − 1)/2). Edges colors range from red to green as the corresponding difference goes from −1 to 1. Solid lines refer to inter-hemispheric connections and dashed lines to intrahemispheric connections. Frontal lobe regions (circles), Insular (ellipse), Limbic (square), Temporal (triangle), Parietal (diamond), Occipital (rectangle) according to classification of Desikan atlas in anatomical lobes [31].
use of a region-specific network statistic which displays low variations across case and controls as well as the choice of an overly conservative level for the FDR and the less power related to mass-univariate local testing procedures.
Our approach doesn't rely on the choice of network summary measures and automatically controls for multiplicity, overcoming previous issues while strongly gaining power. As a result we learn more connections significantly varying between control and AD groups. This provides interesting new insights according to Fig. 6, which displays for each region v (v = 1, . . . ,V ) the total number of connections among v and the remaining V − 1 regions significantly varying between controls and AD group under our local testing procedure (7) with ε = 0.1. To highlight the roles of higher level brain systems, regions are grouped in anatomical lobes according to [31] and in hemispheres.
Results in Fig. 6 highlight the connectivity breakdown for regions in the left hemisphere while providing new insights with respect to [29]. In particular we learn the major role of regions in the left limbic lobe consistently with initial neuropathological studies [26,32] and more recent empirical findings via MRI [33,34] highlighting the key role of the limbic system in memory, attention and executive functioning, while focusing on this lobe as one of the areas mainly affected by AD. Significant changes are also found in the connectivity of the other anatomical lobes such as temporal, parietal and occipital, consistent with [35,36,37,38]. According to Fig. 6 the regions mostly affected by AD in terms of connectivity behavior are the left isthmus of the cingulate (10L), left parahippocampal (16L), left posterior cingulate (23L), left fusiform (7L) and left precuneus (25L). These results provide a unifying answer to different insights arising from several studies, typically focusing on the activity of a subset of regions. Parahippocampal atrophy is found in [39] and [40]; [41] highlights abnormal connectivity in hippocampus and posterior cingulate, while [42] learn reduced functional activity in hippocampus and precuneus, with the latter showing atrophy also in [43]. Metabolic reduction in the posterior cingulate is studied in [44] and [45]. Reduced functional connectivity in the fusiform is found in [46] and [28] via fMRI. Fewer studies are available on the role of the isthmus of the cingulate with only a recent work of [47] trying a first attempt in this direction. We provide a unifying vision, consistent with previous literature, while highlighting the role of the isthmus. This region represents a bridge between the parahippocampal and the posterior cingulate, two critical regions extensively explored in the literature in terms of atrophy and metabolic reduction in AD subjects. Hence a reduced metabolic activity and increased atrophy of parahippocampal and the posterior cingulate, may be related to a disruption of the circuits from the left cingulate isthmus.
Besides providing unifying novel results on brain network variations in AD, our methodology also represents the unique ability to assess evidence of AD according to the subject's full brain network structure. In fact, under our framework, the probability Pr{Y = 1 | L (A i )} that a subject i has AD, conditionally on his brain's structural connectivity network A i is simply equal to where L (A i ) = a i is the network configuration of the i-th subject and p L (A )|y (a i ), (y = 0, 1) can be easily computed from (5). Current classification procedures exploit either region activity vectors ψ i [48] or network summary statistics vectors θ i [49,50], rather than the whole brain network L (A i ), to predict y i . We evaluate our procedure in (9) in terms of insample and out-of-sample classification performance. In the first case, we compute (9) for each subject after considering all data in model estimation. Out-of-sample classification is instead performed by training the model on 69 subjects and predicting the AD status via (9) on the remaining one fourth of the individuals, with the training and test samples randomly selected. Our methodology provides an overall good classification performance, with an area under the ROC curve of 0.91 for in-sample classification and 0.83 for out-of-sample. The accuracy is instead 87% in the former, and 75% in the latter. These results out perform [48], and [49] when summary statistics θ i are extracted from undirected brain networks, while providing similar performance to [50]. It is important to note that [50] utilizes substantially more information in considering both weighted and flow connectivity networks for a total of 298,600 network summary measures, rather than only binary connections encoding presence or absence of fibers.

Discussion
Brain connectivity plays a key role in brain function and dysfunction. Modern magnetic resonance imaging technologies, combined with state-of-the-art data processing algorithms, have made it possible to reliably measure brain structural connectivity networks non-invasively. Understanding brain networks is necessary for developing improved diagnosis and treatment strategies for neurological disorders, but a deep understanding of the relationship between the network of structural interconnections in the brain and such disorders remains still elusive. The statistical methodology presented in this paper defines the first-ever probabilistic generative mechanism to draw tractable and efficient inference directly on the probability mass function associated to a network-valued random variable, rather than on network summary measures or multivariate activity data. In allowing the brain network data to be appropriately analyzed as network-valued, these methods enable substantial improvements in accurately detecting group differences, isolating specific aspects of the network that vary across neurological disorders, and enhancing performance of predictive models as outlined in the application to Alzheimer's disorder.