A connectome is a map of the structural and/or functional connections in the brain. This information-rich representation has the potential to transform our understanding of the relationship between patterns in brain connectivity and neurological processes, disorders, and diseases. However, existing computational techniques used to analyze connectomes are often insufficient for interrogating multi-subject connectomics datasets. Several methods are either solely designed to analyze single connectomes, or leverage heuristic graph invariants that ignore the complete topology of connections between brain regions. To enable more rigorous comparative connectomics analysis, we introduce robust and interpretable statistical methods motivated by recent theoretical advances in random graph models. These methods enable simultaneous analysis of multiple connectomes across different scales of network topology, facilitating the discovery of hierarchical brain structures that vary in relation with phenotypic profiles. We validated these methods through extensive simulation studies, as well as synthetic and real-data experiments. Using a set of high-resolution connectomes obtained from genetically distinct mouse strains (including the BTBR mouse -- a standard model of autism -- and three behavioral wild-types), we show that these methods uncover valuable latent information in multi-subject connectomics data and yield novel insights into the connective correlates of neurological phenotypes.