Network modeling in neuroimaging holds promise in probing the interrelationships among brain regions and potential clinical applications. Two types of matrix-based analysis (MBA) are usually seen in neuroimaging connectomics: one is the functional attribute matrix (FAM) of, for example, correlations, that measures the similarity of BOLD response patterns among a list of predefined regions of interest (ROIs). Another type of MBA involves the structural attribute matrix (SAM), e.g., describing the properties of white matter between any pair of gray-matter regions such as fractional anisotropy, mean diffusivity, axial and radial diffusivity. There are different methods that have been developed or adopted to summarize such matrices across subjects, including general linear models (GLMs) and various versions of graph theoretic analysis. We argue that these types of modeling strategies tend to be "inefficient" in statistical inferences and have many pitfalls, such as having strong dependence on arbitrary thresholding under conventional statistical frameworks.Here we offer an alternative approach that integrates the analyses of all the regions, region pairs (RPs) and subjects into one framework, called Bayesian multilevel (BML) modeling. In this approach, the intricate relationships across regions as well as across RPs are quantitatively characterized. This integrative approach avoids the multiple testing issue that typically plagues the conventional statistical analysis in neuroimaging, and it provides a principled way to quantify both the effect and its uncertainty at each region as well as for each RP. As a result, a unique feature of BML is that the effect at each region and the corresponding uncertainty can be estimated, revealing the relative strength or importance of each region; in addition, the effect at each RP is obtained along with its uncertainty as statistical evidence. Most importantly, the BML approach can be scrutinized for consistency through validation and comparisons with alternative assumptions or models. We demonstrate the BML methodology with a real dataset with 16 ROIs from 41 subjects, and compare it to the conventional GLM approach in terms of model efficiency, performance and inferences. Furthermore, we emphasize the notion of full results reporting through "highlighting," instead of through the common practice of "hiding." The associated program will be available as part of the AFNI suite for general use.