Paper deep dive

NetworkNet: A Deep Neural Network Approach for Random Networks with Sparse Nodal Attributes and Complex Nodal Heterogeneity

Zhaoyu Xing, Xiufan Yu

Year: 2026Venue: arXiv preprintArea: stat.METype: PreprintEmbeddings: 71

Intelligence

Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%

Last extracted: 4/14/2026, 2:37:37 AM

Summary

NetworkNet is a unified deep neural network framework designed to model nodal heterogeneity in count-valued directed random networks. It integrates statistical network modeling with deep learning to estimate latent expansiveness and popularity functions while simultaneously performing high-dimensional attribute selection via hierarchical L1 regularization.

Entities (4)

NetworkNet · methodology · 100%Nodal Heterogeneity · concept · 98%Poisson Distribution · statistical-model · 95%Author-Citation Network · dataset · 90%

Relation Signals (3)

NetworkNet → models → Nodal Heterogeneity

confidence 95% · NetworkNet is a powerful unified approach that allows us to consistently model the nodal heterogeneity

NetworkNet → uses → Poisson Distribution

confidence 95% · NetworkNet models the count-valued random edges with Poisson distributions

NetworkNet → appliedto → Author-Citation Network

confidence 90% · We further apply NetworkNet to a large-scale author-citation network among statisticians

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Heterogeneous network data with rich nodal information become increasingly prevalent across multidisciplinary research, yet accurately modeling complex nodal heterogeneity and simultaneously selecting influential nodal attributes remains an open challenge. This problem is central to many applications in economics and sociology, when both nodal heterogeneity and high-dimensional individual characteristics highly affect network formation. We propose a statistically grounded, unified deep neural network approach for modeling nodal heterogeneity in random networks with high-dimensional nodal attributes, namely ``NetworkNet''. A key innovation of NetworkNet lies in a tailored neural architecture that explicitly parameterizes attribute-driven heterogeneity, and at the same time, embeds a scalable attribute selection mechanism. NetworkNet consistently estimates two types of latent heterogeneity functions, i.e., nodal expansiveness and popularity, while simultaneously performing data-driven attribute selection to extract influential nodal attributes. By unifying classical statistical network modeling with deep learning, NetworkNet delivers the expressive power of DNNs with methodological interpretability, algorithmic scalability, and statistical rigor with a non-asymptotic approximation error bound. Empirically, simulations demonstrate strong performance in both heterogeneity estimation and high-dimensional attribute selection. We further apply NetworkNet to a large-scale author-citation network among statisticians, revealing new insights into the dynamic evolution of research fields and scholarly impact.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

70,948 characters extracted from source content.

Expand or collapse full text

NetworkNet: A Deep Neural Network Approach for Random Networks with Sparse Nodal Attributes and Complex Nodal Heterogeneity Zhaoyu Xing Department of Applied and Computational Mathematics and Statistics University of Notre Dame, and Xiufan Yu Department of Applied and Computational Mathematics and Statistics University of Notre Dame April 14, 2026 Abstract Heterogeneous network data with rich nodal information become increasingly prevalent across multidisciplinary research, yet accurately modeling complex nodal heterogeneity and simultaneously selecting influential nodal attributes remains an open challenge. This problem is central to many applications in economics and sociology, when both nodal heterogeneity and high-dimensional individual character- istics highly affect network formation. We propose a statistically grounded, unified deep neural network approach for modeling nodal heterogeneity in random networks with high-dimensional nodal attributes, namely “NetworkNet”. A key innovation of NetworkNet lies in a tailored neural architecture that explicitly parameterizes attribute-driven heterogeneity, and at the same time, embeds a scalable attribute selec- tion mechanism. NetworkNet consistently estimates two types of latent heterogeneity functions, i.e., nodal expansiveness and popularity, while simultaneously performing data-driven attribute selection to extract influential nodal attributes. By unifying classical statistical network modeling with deep learning, NetworkNet delivers the expressive power of DNNs with methodological interpretability, algorithmic scalability, and statistical rigor with a non-asymptotic approximation error bound. Empirically, simulations demonstrate strong performance in both heterogeneity estimation and high-dimensional attribute selection. We further apply NetworkNet to a large-scale author-citation network among statisticians, revealing new insights into the dynamic evolution of research fields and scholarly impact. Keywords: Nodal Heterogeneity, High-Dimensional Nodal Attributes, Count-Valued Directed Networks, Deep Neural Network, Non-Asymptotic Bound, Consistency 1 arXiv:2604.11673v1 [stat.ME] 13 Apr 2026 1 Introduction Sociological and economic entities frequently interact through complex networks. Within these structures, firms trade goods, banks extend credit, and researchers disseminate knowledge through citations. A key characteristic of these networks is the underlying heterogeneity of the individuals (e.g., firm size, risk profile, research expertise) that shape both how actively they initiate interactions and how strongly they attract connections from others. Quantifying this individual-level heterogeneity is therefore a crucial step in modeling network formation, evaluating systemic evolutions, and designing effective policies. In practice, network data with rich nodal information are ubiquitous in modern science. Alongside rapidly expanding network scale, modern network datasets exhibit greater struc- tural complexity and growing granularity in edge information, evolving from simple binary undirected networks to weighted directed networks. Meanwhile, nodes themselves are often accompanied by rich attributes describing the individual entities, as seen in online social networks with user-specific features (Leskovec & Mcauley 2012, Xing, Wang, Wood & Zou 2025), academic collaboration networks annotated with scholars’ research interests (Gao et al. 2023, 2025), economic networks with different individual characteristics (Antinyan et al. 2020, Xing & Zhong 2024, Xing, Tan, Zhong & Shi 2025), and among others. Classical random network models primarily focus on the topological structure and treat all nodes as statistically equivalent (Erdős & R ̋enyi 1959, Watts & Strogatz 1998, Schweinberger 2011, Chatterjee & Diaconis 2013). In practice, however, real-world networks are formed among intrinsically heterogeneous entities, where influential node-level attributes can substantially impact the generative mechanisms of network connections. Modeling nodal heterogeneity, driven by the practical needs in network analysis, has thus become a critical yet challenging problem in network science. For example, in academic social networks, established researchers typically exert greater collaborative and academic 2 influence than early-career researchers who often possess fewer contributions concentrated in limited research areas. Overlooking such inherent individual nodal heterogeneity risks oversimplifying the underlying connection mechanisms and may lead to biased assessments of scholarly impact. As a classic application in social networks, academic collaboration and scholarly impact have long been central topics. Nevertheless, previous findings on scholarly impact analysis are solely based on network structures, and do not fully account for the effect of nodal heterogeneity (Lehmann et al. 2003, Moody 2004, Che et al. 2025). Although some approaches attempt to capture nodal heterogeneity using local topological structures (Ma et al. 2025, Yan et al. 2016), they rely exclusively on the network structures and are unable to fully utilize the high-dimensional nodal attributes now widely available in modern empirical data. The availability of rich nodal information in modern academic network data opens new opportunities for characterizing heterogeneity beyond network topology alone. A large academic network dataset (Gao et al. 2023) was recently released, containing count-valued edges and high-dimensional nodal attributes for 168,171 scholars based on 97,436 papers in representative statistical journals over a four-decade span. This dataset enables empirical investigations of attribute-driven nodal heterogeneity within an edge-generative random network, and provides a unique basis for studying how individual-level attributes affect the scholarly impact through heterogeneous network generative processes. This paper is empirically motivated by both the practical need for principled quantification of scholarly impact using this newly released public dataset of academic networks with informative count-valued edges and rich nodal attributes, and the lack of methodologies for modeling nodal heterogeneity with high-dimensional nodal attributes. We seek to fill this methodological gap and address the following fundamental question: 3 How to better quantify the nodal heterogeneity with high-dimensional nodal attributes? Prior works in network modeling reveal a clear methodological divide between statistics and deep learning perspectives: statistical network models principally focus on characterizing randomness in network topology without modeling the attribute-driven nodal heterogeneity, while deep learning approaches typically analyze observed networks as fixed input struc- tures without considering the nodal heterogeneity in the generating process of random networks. To our knowledge, no method yet adequately unifies the strengths of statistical and deep learning techniques for modeling nodal heterogeneity using high-dimensional nodal attributes. Specifically, related statistical network models treat the networks as random variables, and largely focus on the generating processes of network structures as well as the corresponding network distributions. Among them, the exponential random graph models (ERGMs) cap- ture structural randomness through their topology factors (Hunter et al. 2008, Schweinberger et al. 2020). However, ERGMs are computationally intensive and face inferential challenges, which become even more demanding for extensions to count-valued networks (Krivitsky 2012, Huang & Butts 2024), and their classical formulations assume node exchangeability that limits their ability to model nodal heterogeneity or incorporate influential node-specific variation. Stochastic block models capture the random interactions in a network with estimated community information of nodes (Lei & Rinaldo 2015, Lei & Lin 2023, Zhang & Amini 2023). Although they allow the different probabilities for connection among different communities, they usually fail to capture the individual nodal heterogeneity as the nodes within one community are often assumed to be identical. Degree-based models, such as the beta-model and its variants (Yan et al. 2016, Graham 2017, Stein et al. 2025), utilize the degree heterogeneity to capture the randomness of network structures where latent parame- ters are estimated to describe the nodal heterogeneity. However, previously proposed classic 4 beta-models focus on degrees of binary networks and model the network structures only. While few recent extensions can incorporate covariates to capture the nodal heterogeneity (Wang, Cai, Niu & Li 2024, Yan et al. 2019), they are specially designed for edge-wise features and assume restrictive linear relationships. Commonly collected nodal attributes are not directly used to capture the nodal heterogeneity in random network models. Despite the clear practical need, a general statistical method that can model complex count-valued networks with high-dimensional nodal attributes and non-linear effects remains an open challenge. Parallel to statistical network models, deep learning fields have recently made significant strides in network analysis, particularly in tasks such as link prediction and node classification (Yuan & Qu 2023, Abbe 2018). However, almost all of these related methods treat networks as fixed given structures that serve as inputs of the algorithms without considering the randomness in network generating process. A growing body of work has shown that the neural network approaches for network data can be highly sensitive to random noise and attacks in data (Zügner et al. 2018, Wang, Liu, Liu, Wang, Medya & Yu 2024, Goodfellow et al. 2015, Szegedy et al. 2013). Overlooking the inherent stochasticity of networks is one of the main reasons for overfitting problems while modeling network data. Besides, most deep learning-based graph models function as black boxes that offer limited interpretability within the theoretical context of network science, despite their strong predictive performance in practice. While some “explainer” models are proposed for graph neural networks (Vu & Thai 2020, Ying et al. 2019), they focus mainly on estimating the feature importance for specific learning tasks but do not evaluate their influence on network distribution. Thus, without accounting for the probabilistic generating nature of random networks, these deep learning-only approaches are unable to effectively capture and interpret the nodal heterogeneity. 5 Despite rapid methodological advances in statistics and deep learning techniques for network data, there remains a lack of methods capable of explicitly modeling nodal heterogeneity using high-dimensional nodal attributes in random networks, leaving a gap in the literature. To close this critical gap, we propose a novel framework for modeling nodal heterogeneity in count-valued directed random networks with high-dimensional nodal attributes, namely “NetworkNet”. NetworkNet is a powerful unified approach that allows us to consistently model the nodal heterogeneity utilizing both large-scale count-valued networks and high- dimensional nodal information, and effectively identify the influential subset of nodal attributes. As a synthesis of statistical and deep learning methods, NetworkNet models the count-valued random edges with Poisson distributions, where the heterogeneous rates for edge generation are determined by both the sender’s and receiver’s nodal heterogeneity. Influential nodal attributes are selected from the set of high-dimensional nodal attributes and used to capture complex nodal heterogeneity through specially designed deep neural networks (DNN). For the architectures of NetworkNet, we introduce two skip layers with corresponding constraints based on hierarchicalL 1 regularization to accurately estimate the nodal heterogeneity and select nodal attributes simultaneously. Theoretically, we show that the NetworkNet estimator can approximate the latent nodal heterogeneity functions consistently and achieve the non-asymptotic upper bound of the error. This work contributes to the network literature by addressing several critical challenges in a unified framework. First, we introduce a nodal attribute selection mechanism designed for the network context, enabling us to sift through high-dimensional nodal information and identify meaningful attributes. Second, we propose a principled statistical model for count-valued networks that directly targets the distribution of edge counts, effectively accommodating nodal heterogeneity in the interaction intensity. Third and most notably, we synthesize statistical methods with modern deep learning techniques by embedding flexible 6 neural components within a coherent probabilistic framework, bridging statistical rigor and deep learning expressiveness. This allows us to leverage the representational power of DNNs for complex nodal information while retaining a statistical-model-based treatment of edge counts and a clear link between nodal attributes and edge intensities. Together, NetworkNet offers a powerful tool that is simultaneously distributionally appropriate for count-valued edges, scalable to high-dimensional nodal attributes, and capable of integrating predictive flexibility with interpretable statistical learning. Furthermore, NetworkNet is not limited to count-valued network data, but offers a flexible and unified modeling framework that can be readily extended to capture nodal heterogeneity in a broad class of network types, such as undirected networks, binary networks, and continuously weighted networks. Empirically, the proposed general approach NetworkNet can be applied broadly to analyze not only academic collaboration networks, and also other practical network structures with nodal attributes, such as the social networks with personal characteristics, business networks with companies’ features and international trading networks with macroeconomic indexes. The rest of this paper is organized as follows. Section 2 formally introduces the NetworkNet model to estimate the nodal heterogeneity and select relevant nodal attributes with an alternating optimization algorithm. We demonstrate the dominance of NetworkNet in both nodal heterogeneity estimation and nodal attribute selection with comprehensive simulation studies in Section 3. In Section 4, we utilize the proposed NetworkNet to conduct an in-depth analysis of the academic author-citation networks over a four-decade span, investigating the evolution of scholarly impact and uncovering data-driven insights. Section 5 concludes with a discussion. Proofs of all theoretical results and additional details of empirical analysis are included in the Supplementary Materials. 7 2 Methods To introduce our approach to model nodal heterogeneity in random networks with high- dimensional nodal attributes clearly, we first introduce the necessary notation and basic model settings. Then we introduce the framework of NetworkNet, followed by the iterative optimization algorithms. Subsequently, theoretical results, including a non-asymptotic bound, are provided. 2.1 Notation We consider a count-valued directed network with n vertices and corresponding adjacency matrixA= (A ij ) n×n ,i,j ∈[n], whereA ij ∈Nrepresents the count-valued weight of the edge between the nodesiandj, and [n] =1,· ,n. For any constanta ∈R,⌊a⌋ represents the largest integer that is smaller than or equal toa. Besides the network structures, ap-dimensional vector of nodal covariates associated with each node is denoted byx i ∈R p ,i ∈[n]. These covariates may include both continuous and categorical attributes. The complete set of nodal attributes for the network is represented by the matrix X= [x 1 ,x 2 ,...,x n ] ′ ∈R n×p , which characterizes the nodal heterogeneity. For a function η(x i )∈F,i∈[n], theL n,2 norm defined on classFis∥η(·)∥ n,2 = [ P n i=1 η(x i ) 2 /n ] 1/2 . For any compact spaceX ⊆[a,b] p with reference measureμ, define∥η∥ L 2 =[ R X η(x) 2 dμ(x)] 1/2 for a measurable functionη:R p →R. For the pairη= (α,β)∈G(F) from a function class G(F), define∥η∥ L 2 = ∥α∥ 2 L 2 +∥β∥ 2 L 2 1/2 . ⌈r⌉ represents a ceiling function that indicates the smallest integer that is greater than or equal tor. For a vectorθ= (θ 1 ,· ,θ n ), ∥θ∥ 1 = P n i=1 |θ i |represents theL 1 -norm of the vectorθ. We denote scalars and functions using lowercase letters, and denote matrices and vectors by bold uppercase and bold lowercase letters, respectively. 8 2.2 NetworkNet Motivated by the author-citation networks, where citation interactions are initiated at irregular time points and citation counts represent the aggregation of these discrete random events, we model the count-valued network data via a set of Poisson distributions with latent parameters as A ij |λ ij ∼ Poisson(λ ij ),(1) where the non-negative rate parameterλ ij represents the expected number of interactions from nodeito nodejwithin a time range. For academic author-citation networks, the number of citations coming from different publications can only be count-valued data, which is modeled by(1)with Poisson distribution, and the general strength of the connection between scholars is captured by the Poisson parameterλ, which is assumed to be determined by both the expansiveness of origin node i and the popularity of the destination node j. The high-dimensional nodal attributes give a straightforward description about individual characteristics that are usually strongly related to the network generating process. To better model the nodal heterogeneity in the generating process of random graphs, we introduce two types of latent node-specific components based on high-dimensional nodal attributes, the expansiveness functionα(·) capturing the nodal inherent heterogeneity in generating connections to the other nodes, and the popularity functionβ(·) capturing the nodal inherent heterogeneity in attracting connections from the other nodes. Both of them are functions of nodal attributes from a general function class and capture two types of the nodal heterogeneity. Then we model the rateλ ij of the count-valued edgesA ij in(1)with both the expansiveness of the edge-sender and popularity of the edge-receiver as log(λ ij ) = [α(x i ) + β(x j )]/Z n ,(2) whereZ n is a scaling constant that ensures that expected edge counts remain well-behaved as the network grows, analogous to the sparsity-inducing scaling commonly adopted in 9 large-network asymptotics. As a constant shift can lead to model identifiability issues, we set P n i=1 α(x i ) = P n j=1 β(x j ), while other constraints can also be used. The flexible expansiveness functionα(·) and popularity functionβ(·) in our model settings can capture complex non-linear relationships between nodal attributes and nodal hetero- geneity that is directly related to the probability distribution of random networks. Modeling heterogeneity with high-dimensional and sparse nodal attributes in random network models is an open challenge. To address this problem, we embed flexible neural techniques within the coherent probabilistic framework to model the nodal heterogeneity in large-scale random networks, and select the high-dimensional nodal attributes simultaneously. Specifically, we approximate latent nodal heterogeneity functionsα(·) andβ(·) with deep neural networks from the following class H =f ≡ f(x,W,θ) : x7→ θ ′ x + h W (x),(3) whereθrepresents the coefficient vector for the skip layer andh W (·) represents a ReLU-DNN with weightsW. BothWandθare parameters to be estimated. To select the influential nodal attributes while modeling nodal heterogeneity within a probabilistic framework, we formulate a regularized optimization problem minimizing a composite objective function with the negative log-likelihood of the network under the Poisson assumption, additional regularization terms and constraints: min θ α ,θ β ,W α ,W β X i̸=j L ij (θ α ,θ β ,W α ,W β ) + λ 1 ∥θ α ∥ 1 + λ 2 ∥θ β ∥ 1 ,(4) s.t. ∥W (1) α,j ∥ ∞ ≤ M|θ α,j |,∥W (1) β,j ∥ ∞ ≤ M|θ β,j | and n X i=1 α(x i ) = n X j=1 β(x j ), where two neural networksf(X i ;W α ,θ α ) andg(X j ;W β ,θ β ) with skip-layers are used to approximate latent nodal heterogeneity functionsα(·) andβ(·), and the empirical loss 10 function is L ij (θ α ,θ β ,W α ,W β ) = exp [f(x i ;W α ,θ α )/Z n + g(x j ;W β ,θ β )/Z n ] −A ij [f(x i ;W α ,θ α )/Z n + g(x j ;W β ,θ β )/Z n ](5) with observed networkA=A ij and nodal attributesX. TheL 1 penalty is applied to shrink the weights of the skip layers with tuning-parametersλ 1 ,λ 2 >0, andθ α andθ β serve as the constraints for the first-layer weightsW (1) α andW (1) β in DNNs. When ( ˆ θ α ) j = 0, then the corresponding featurejis effectively removed from the subset of influential nodal attributes for nodal expansiveness, thus achieving data-driven nodal attributes selection. Figure 1: NetworkNet Methodology The architecture of NetworkNet is shown in Figure 1. The algorithm starts with the inputs of the count-valued networks with high-dimensional nodal attributes on the left-hand side. Each dimension of the nodal attributes is treated as the input of two separate parts of the DNNs in dash boxes, where the weights of blue and green linkages are trained within deep neural networks to approximate the latent expansiveness functionα(·) and popularity functionβ(·), respectively. Meanwhile, two skip layers are introduced to select nodal 11 attributes in each iterating update in Algorithm 1, where a hierarchical proximal operator described in Algorithm 3 is used for specially designed constraints in(4). As NetworkNet models nodal heterogeneity with random network models with high-dimensional nodal attributes, we have the additional Network Layer marked in yellow demonstrates that the observed count-valued network structures are included in the empirical loss function defined in(5). As illustrated in Figure 1, the architecture of NetworkNet is specifically designed to unify advanced principles from statistical random network models and deep learning. This synergistic approach offers a robust methodology capable of explicitly modeling complex nodal heterogeneity within random networks characterized by high-dimensional attributes. 2.3 Algorithms We propose an efficient and stable alternative optimization algorithm to effectively train the specially designed deep neural networks in Figure 1 with the objective function in(4). The core idea is to utilize the hierarchical proximal operator to shrink the first-layer weights of neural networks, and decouple the approximations of expansiveness functionα(·) and popularity function β(·) into two iterative steps. Algorithm 1 NetworkNet Algorithm Require:Adjacency matrixA ∈N n×n , nodal attributes matrixX ∈R n×p ; Hyper- parameters λ 1 ,λ 2 ,γ,M, learning rate ρ, tolerance ε 1: Initialization of ( ˆ α (0) , ˆ β (0) ) 2: for t = 0 to T max − 1 do 3: ( ˆ β (t+1) , (W β ,θ β ))← Update(g,X,A, ˆ α (t) ,λ 2 ,γ,ρ,M) for β 4: ( ˆ α (t+1) , (W α ,θ α ))← Update(f,X,A, ˆ β (t+1) ,λ 1 ,γ,ρ,M) for α 5: if ∥ ˆ α (t+1) − ˆ α (t) ∥ 2 < ε and ∥ ˆ β (t+1) − ˆ β (t) ∥ 2 < ε then 6:break 7: end if 8: end for 9: ˆ α← ˆ α (t+1) , ˆ β ← ˆ β (t+1) 10: S α ←k :|(θ α ) k | > 0, S β ←k :|(θ β ) k | > 0 11: return Estimated parameters ˆ α, ˆ β; Selected feature sets S α ,S β Specifically, we model two types of nodal heterogeneity via two separate neural networks 12 f(·;W α ,θ α ) andg(·;W β ,θ β ) from the classHwith skip-layers. Letα= [α(x 1 ),...,α(x n )] ⊤ andβ= [β(x 1 ),...,β(x n )] ⊤ represent the true nodal heterogeneity parameter vectors, and denote the estimated nodal heterogeneity parameter vectors as ˆ α= [f(x 1 ),· ,f(x n )] ⊤ and ˆ β= [g(x 1 ),· ,g(x n )] ⊤ . The iterative optimization procedures are straightforward as demonstrated in Algorithm 1: we update the network parameters (θ β ,W β ) forg(·) holding the ˆ α fixed; and then update the network parameters (θ α ,W α ) forf(·) holding the ˆ βfixed. The details of optimization in these iterative updates are given in Algorithm 2. As the hierarchy constraints are separable over the features but the identifiability constraint is related to both two updates, we optimize the objective in(4)by constrained proximal gradient descent, which is given in Algorithm 3. Algorithm 2 The Update Procedure (Steps 3 & 4 in Algorithm 1) Require: f,g,X,A,λ fixed ,λ,γ,ρ,M, let W (1) ∈R p×d 1 be the first layer weights 1: Let (W,θ) be the parameters of network f. 2: if to update α then 3: Define loss function for α: L(W,θ) = X i̸=j h e L α − A ij L α i + λ∥θ∥ 1 + γ " X k f W (x k )− X k λ fixed,k # 2 where L α = f W (x i ) + λ fixed,j . 4: else 5: Define loss function for β: L(W,θ) = X i̸=j h e L β − A ij L β i + λ∥θ∥ 1 + γ " X k λ fixed,k − X k g W (x k ) # 2 where L β = λ fixed,i + g W (x j ). 6: end if 7: for k = 1 to T max do 8: Gradient descent: (W,θ)← (W,θ)− ρ∇L(W,θ) 9: Proximal step: (W (1) ,θ)← HierarchicalProx(W (1) ,θ,ρλ,M) 10: end for 11: ˆ λ new ← [f(x 1 ),...,f(x n )] ⊤ 12: Return: ˆ λ new , (W,θ) There are three hyper-parameters in the NetworkNet optimization: the regularization 13 Algorithm 3 Hierarchical Prox Procedure (Step 9 in Algorithm 2) Require: Parameters W , θ,τ, M 1: for k = 1 to p do 2: θ new,k ← sign(θ k )· max(|θ k |− τ, 0) 3: Let w k,· be the k-th row of W (1) , then update w k,·,new ← w k,· min 1, M ·|θ new,k | ∥w k,· ∥ 2 ! 4: end for 5: return (W (1) new ,θ new ) parametersλ 1 andλ 2 , and the balance parameterMthat governs the trade-off between the neural-net components and skip layers. We employ the High-dimensional Bayesian Information Criterion (HBIC, Wang et al. 2013) for hyperparameter tuning in practical applications with large-scale networks, where parallel computing is applicable. 2.4 Theoretical Results The nodal heterogeneity functionsα(·) andβ(·) of nodal attributesXcan be complex and highly non-linear. Below, we consider them belonging to a general class of functions that are Hölder smooth and derive the non-asymptotic upper bound of the NetworkNet approximation. We prove that, with mild assumptions, the proposed NetworkNet can consistently approximate the latent nodal heterogeneity functions and give consistent estimations of individual expansiveness and popularity. Definition 1 (Hölder smooth function). A functionf:R s →Ris said to beγ-Hölder smooth if all ⌊γ⌋-order partial derivatives of f exist and satisfies the Hölder condition f (⌊γ⌋) (x)− f (⌊γ⌋) (y) ≤ a∥x− y∥ γ−⌊γ⌋ 2 for some constant a. Definition 2 (Hölder smooth compositional function). Consider the following collection of parameters withJ ∈N + : dimension vectors k = (k 1 ,...,k J+1 ) ⊤ ∈N J+1 + , intrinsic 14 dimension vectors s = (s 1 ,...,s J ) ⊤ ∈N J + , smoothness parametersγ= (γ 1 ,...,γ J ) ⊤ ∈R J + , and domain bounds a = (a 1 ,...,a J ) ⊤ ∈R J and b = (b 1 ,...,b J ) ⊤ ∈R J . A function f: [a 1 ,b 1 ] k 1 →R k J+1 belongs to the classFof Hölder smooth compositional functions if it admits the representationf(x) = (e J ◦·◦ e 2 ◦ e 1 )(x), ∀x∈[a 1 ,b 1 ] k 1 ,where for each u= 1,...,J, the mape u : [a u ,b u ] k u → [a u+1 ,b u+1 ] k u+1 is constructed component-wise: each componente uv : [a u ,b u ] s u →[a u+1 ,b u+1 ] is aγ u -Hölder smooth function that depends on at most s u of the k u input coordinates. Assumption 1 (Smoothness of nodal heterogeneity functions). DenoteFas the class of γ u -Hölder smooth compositional functions in terms of nodal attributesx u ,∀u ∈[p]. The true nodal heterogeneity function α 0 ,β 0 ∈F. Denoteη 0 = (α 0 ,β 0 )∈ G(F) as the true nodal heterogeneity function with both nodal expansiveness and nodal popularity, whereG(F) :=(α,β) :α∈F,β ∈F. Assumption 1 requires that the complex latent nodal heterogeneity functions are relatively smooth in terms of the nodal attributes. The mild assumption gives a wide class of functions to model the relationship between nodal attributes and nodal heterogeneity (Liu et al. 2020, Kohler & Langer 2021, Schmidt-Hieber 2020). Assumption 2 (Sparsity of nodal attributes). Assume the subset of true relevant nodal attributes satisfies |A α | +|A β | = r α + r β ≪ p where r α + r β = o( √ n log 2 (n)) and p = O n 3 log 3 (n) . For network data with rich nodal information, we usually consider that the high-dimensional true nodal attributes are sparse. As both the true relevant features and dimensions of nodal attributes can increase with sample size, and the fraction of true nodal attributes is (r α +r β )/p=o n −5/2 log −1 (n) , Assumption 2 gives a mild requirement about the sparsity of nodal attributes. 15 Assumption 3 (DNN technical conditions). There exists anL-layer ReLU-DNN classH with layer widths d v L v=2 and minimum width W = min v=2,·,L d v such that 1. L > J + 216( P J u=1 ⌈γ 2 u ⌉ + 1); 2. W ≥ 81d L max 1≤u≤J r u (⌈γ u ⌉ + s u + 2) s u +1 3 s u +1 ; 3. L max 2≤v≤L d v ≍ L min 2≤v≤L d v ≍ n ̃s 2(2 ̃γ+ ̃s) , where ̃γ= ̃γ ̃u and ̃s=s ̃u , in which the index ̃u is defined by ̃u = arg min 1≤u≤J ̃γ u s u with ̃γ u = γ u Q J v=u+1 (γ v ∧ 1). Assumption 3 is a commonly used technical assumption for DNNs (Schmidt-Hieber 2020, Liu et al. 2020). With the width and length of neural networks increasing with sample size, NetworkNet can approximate the Hölder nodal heterogeneity functions based on nodal attribute matrix X ∈R n×p with the non-asymptotic bound given in Theorem 1. Theorem 1. Denote the estimated nodal heterogeneity functions asˆη= (ˆα, ˆ β) and the true nodal heterogeneity functions asη 0 = (α 0 ,β 0 ). Under Assumption 1-3 and??in the Supplements, one has sup η 0 ∈G(F) ∥ˆη− η 0 ∥ n,2 = O p n − ̃γ 2 ̃γ+ ̃s log 2 (n) Theorem 1 shows the consistency in estimation of nodal heterogeneity functions, where the convergence rate achieves the classical nonparametric regression rate (Schmidt-Hieber 2020). Based on consistent approximation of nodal heterogeneity functions, NetworkNet can give a consistent estimation for each node about its expansiveness and popularity. Furthermore, we explore the general approximation error bound of expansiveness and popularity function over the compact attribute space. DenoteA=A α ∪A β represent the active subspace, and letr=r α +r β , then|A| ≤ r ≪ pbased on Assumption 2. With the following mild assumptions, we provide a general non-asymptotic bound for nodal heterogeneity approximation in Theorem 2. 16 Assumption 4 (Design regularity on active subspace). The fixed design pointsx 1 ,...,x n satisfy the following with respect to Borel reference measure μ and active subspace A: (i) LetX A =x A :x∈Xbe the projection onto the active subspace. The fill distance h n,A = sup z∈X A min 1≤i≤n ∥z− (x i ) A ∥ 2 satisfies h n,A = O(n −1/r ). (i) There exists C μ > 0 (independent of n) such that the Voronoi cells V i = n x∈X :∥(x) A − (x i ) A ∥ 2 ≤∥(x) A − (x j ) A ∥ 2 , ∀j ̸= i o satisfy max 1≤i≤n μ(V i )≤ C μ /n. Assumption 5 (Lipschitz regularity of the neural network class). There exists a sequence b n >0 such that everyf ∈ Hsatisfying the constraints in optimization problem (4) is b n -Lipschitz on X A : |f(x)− f(x ′ )|≤ b n ∥(x) A − (x ′ ) A ∥ 2 , ∀x,x ′ ∈X, with b n = O n ρ / log 2 (n) , where ρ = ̃γ/(2 ̃γ + ̃s). It is worth noting that Assumption 4 is related to unknown active setAonly, where the raten −1/r is mild and commonly-used (Wendland 2004, Narcowich et al. 2005, Pronzato & Zhigljavsky 2023) fors-dimensional space. In practice, it requires that the nodal attributes are reasonably spread over the low-dimensional active subspace, while the nodal attributes neither cluster nor leave large gaps. Theorem 2. Denote the estimated nodal heterogeneity functions asˆη= (ˆα, ˆ β) and the true functions asη 0 = (α 0 ,β 0 ). Under Assumptions in Theorem 1 and the additional Assumptions 4 and 5, one has sup η 0 ∈G(F) ∥ˆη− η 0 ∥ L 2 = O p n − ̃γ 2 ̃γ+ ̃s log 2 (n) + n τ log 2 (n) ! ,(6) where τ = ̃γ/(2 ̃γ + ̃s)− 1/r. We explain the conclusion in Theorem 2 in several regimes in the following Corollary 1. 17 Corollary 1. Letρ= ̃γ/(2 ̃γ+ ̃s). Under the assumptions of Theorem 2, the following regimes hold. (a)Ifτ <−ρ, or equivalently 2 ̃γ r <2 ̃γ+ ̃s,thenn τ / log 2 (n) =o(n −ρ log 2 (n)). Theorem 2 and Theorem 1 share the same rate. (b) If −ρ < τ < 0, or equivalently 2 ̃γ + ̃s < 2 ̃γ r < 2(2 ̃γ + ̃s), then ∥ˆη− η 0 ∥ L 2 = O p n τ log 2 (n) ! = O p n ρ−1/r log 2 (n) ! . Remark 1. The condition 2 ̃γ r < 2 ̃γ + ̃s, equivalently r < 1 + ̃s/(2 ̃γ), is a joint condition on the smoothness parameters ( ̃γ, ̃s) and the total number of active featuresr. For a fixedr, the condition 2 ̃γ r < 2 ̃γ + ̃s requires ̃s to be large relative to ̃γ. 3 Simulations To evaluate the numerical performance of NetworkNet on both estimation of nodal het- erogeneity and nodal attribute selection, we generate the random count-valued directed network from model(1)with sizen= 100. The nodal attribute matrixX ∈R n×p is then generated independently from a uniform distribution,x ij ∼ Unif(−1,1),i ∈[n],j ∈[p] with dimensionp= 1000. As shown in(2), the complex nodal heterogeneity, including both expansiveness and popularity, is determined by high-dimensional nodal attributes. To fully show the effectiveness of NetworkNet in modeling complex nodal heterogeneity, we consider the following two different settings to simulate the expansiveness functionα(·) and popularity function β(·): a linear setting α(x i ) = 5 X k=1 x i,k , β(x j ) = 10 X k=6 x j,k (7) 18 and a nonlinear setting α(x i ) = 5 [|x i,1 | +|x i,2 | + log(|x i,3 |) + log(|x i,4 | +|x i,5 |)] β(x j ) = 5 [|x j,6 | +|x j,7 | + log(|x j,8 |) + log(|x j,9 | +|x j,10 |)](8) As shown in equations(7)and(8), we set the true nodal attributes forαasx 1 tox 5 , and true nodal attributes forβare fromx 6 tox 10 . Correspondingly, the indicator of the true subset ofα-related nodal attributes isA α =1,2,3,4,5, and indicator of the true subset of β-related nodal attributes is A β =6, 7, 8, 9, 10 in all cases. We generater= 1,2,· ,100 replications of each of the simulation settings above. Denote the estimated nodal expansiveness and popularity for nodei∈[n] inr-th round of replications asˆα (r) i and ˆ β (r) i , respectively. Denote the estimated subsets ofα-related features andβ- related features as ˆ A (r) α and ˆ A (r) β respectively. With known true nodal heterogeneity (α 0 ,β 0 ) and true subset of nodal attributesA α ,A β , we evaluate and compare different approaches with the following classic evaluations: To assess the performance in estimating the nodal heterogeneity for each method, we compute the average root mean squared error for both αandβfor comparison based on the simulation results; to evaluate the performance in selecting relevant nodal attributes, we consider the precision, true-positive-rate and F1 value for both expansiveness-related selection and popularity-related selection. The definitions of these evaluations are shown in Table 1. As there are no specialized approaches in previous literature to model the nodal hetero- geneity in a count-valued random networks with high-dimensional nodal attributes that can also select the nodal attributes simultaneously, we consider the following 2 two-stage combinations of classic statistical and deep learning algorithms as competitors: (1) The classic Lasso estimator (Tibshirani 1996) based on maximum likelihood estimates (denoted as “MLE+Lasso”) and (2) The DNN estimator Lassonet (Lemhadri et al. 2021) based on maximum likelihood estimates (denoted as “MLE+LassoNet”). For both of the two-stage 19 Table 1: The definitions of the simulation performance measures to assess the performance of estimation of α(·) and β(·), and selection of nodal attributes. MeasuresNotation Formula Estimation Average root mean squared error for α(·)α-RMSE 100 −1 P r [ 1 n P n i=1 (ˆα (r) i − α i,0 ) 2 ] 1/2 Average root mean squared error for β(·)β-RMSE 100 −1 P r [ 1 n P n i=1 ( ˆ β (r) i − β i,0 ) 2 ] 1/2 Selection Average precision for α-related attributesα-Precision 100 −1 P r | ˆ A (r) α ∩A α |/| ˆ A (r) α | Average precision for β-related attributesβ-Precision 100 −1 P r | ˆ A (r) β ∩A β |/| ˆ A (r) β | Average true positive rate of α-related attributes α-TPR100 −1 P r | ˆ A (r) α ∩A α |/|A α | Average true positive rate of β-related attributes β-TPR100 −1 P r | ˆ A (r) β ∩A β |/|A β | Average F1 value for α-related attributesα-F1100 −1 P 100 r=1 2·| ˆ A (r) α ∩A α |/(| ˆ A (r) α | +|A α |) Average F1 value for β-related attributesβ-F1100 −1 P 100 r=1 2·| ˆ A (r) β ∩A β |/(| ˆ A (r) β | +|A β |) Note: The parameter r refers to the r-th replication of the simulations. methods, the maximum likelihood estimates for the two parameters of each individual node are numerically approximated based on model(1)with observed network structures. Then the estimated vectors of nodal expansiveness and popularity serve as the responses in both Lasso and Lassonet algorithms for nodal attribute selection. For all simulation cases, the network architectures are set to be (32,16) and training epochs are set to be 1000 with random initialization. The tuning parameters are chosen by high-dimensional Bayesian information criteria, and the optimal models are used to approximate the nodal heterogeneity function and estimate the nodal heterogeneity ˆ α and ˆ β. Tables 2 and 3 summarize the simulation results under linear and non-linear settings, respectively, and demonstrate the superior performance of the proposed method, NetworkNet, in both nodal heterogeneity estimation and nodal attribute selection. For the estimation of nodal heterogeneity, we use RMSE to quantify the distance between estimated nodal heterogeneity parameters ( ˆ α, ˆ β ) and true expansivenessαand popularity β. In all simulation settings, the maximum likelihood estimator has the largest average RMSE as it uses only information about network structure. Meanwhile, all the other three approaches benefit from the additional nodal attribute information. Specifically, the 20 LassoNet method, which includes Lasso as a special case, gives a smaller RMSE compared to Lasso-based maximum likelihood estimations when the latent nodal heterogeneity is non-linear. The proposed NetworkNet gives the estimations of heterogeneity parameters with much lower RMSE in both linear and nonlinear nodal heterogeneity settings. For the selection of nodal attributes, we evaluate the performance by comparing to the known true subset. As the maximum likelihood estimation does not perform selection, all the evaluations of attribute selection for MLE are blank (/). As the simulation results show, the proposed NetworkNet can select the nodal attributes accurately with precision, TPR, and F1 score close to 1. While both the two-stage methods, MLE+Lasso and MLE+Lassonet, cannot select the true nodal attributes, they usually have good performance in classic independent data. One of the reasons is that both two-stage estimates are highly dependent on the first stage MLE, and cannot capture the complex nodal heterogeneity based on nodal attributes in practice. However, the proposed NetworkNet, combining the statistical model of random networks and deep neural network techniques, can effectively model the nodal heterogeneity while accurately selecting nodal attributes. Table 2: Performance of estimation and nodal attribute selection for linear nodal hetero- geneity NetworkNet MLE+Lasso MLE+LassoNetMLE α RMSE 0.5870 (0.2623) 1.3196 (0.0835) 1.4785 (0.0941) 1.3758 (0.0833) Precision 0.9983 (0.0167) 0.0050 (0.0286) 0.0298 (0.0602)/ TPR0.9680 (0.1497) 0.0060 (0.0343) 0.0480 (0.0990)/ F10.9737 (0.1218) 0.0055 (0.0312) 0.0365 (0.0739)/ β RMSE 0.6163 (0.1142) 1.3081 (0.0876) 1.2672 (0.0896) 1.3654 (0.0885) Precision 0.7179 (0.1783) 0.0033 (0.0235) 0.0121 (0.0563)/ TPR1.0000 (0.0000) 0.0040 (0.0281) 0.0100 (0.0438)/ F10.8230 (0.1252) 0.0036 (0.0256) 0.0106 (0.0471)/ Note: The evaluations of attribute selection for MLE are blank (/) as the MLE method does not perform selection. For nonlinear cases, the complex nodal heterogeneity enlarges the advantages of the proposed NetworkNet. Specifically, the MLE estimate, which is based only on network structure 21 Table 3: Performance of estimation and nodal attribute selection for non-linear nodal heterogeneity NetworkNet MLE+Lasso MLE+LassoNetMLE α RMSE 1.8608 (0.3217) 4.0095 (0.3761) 3.6234 (0.4250) 5.2026 (0.4876) Precision 0.9105 (0.1791) 0.0131 (0.0468) 0.0503 (0.0646)/ TPR0.9880 (0.0844) 0.0180 (0.0642) 0.1620 (0.1943)/ F10.9327 (0.1414) 0.0151 (0.0539) 0.0737 (0.0870)/ β RMSE 2.0015 (0.4917) 4.2462 (0.4712) 3.7348 (0.3889) 5.5397 (0.7502) Precision 0.8291 (0.2872) 0.0150 (0.0480) 0.0375 (0.0760)/ TPR0.9300 (0.2564) 0.0220 (0.0690) 0.0540 (0.1058)/ F10.8662 (0.2651) 0.0177 (0.0562) 0.0426 (0.0832)/ Note: The evaluations of attribute selection for MLE are blank (/) as the MLE method does not perform selection. information, is less accurate compared to linear cases, which makes the MLE-based two-stage methods perform poorly on both estimation and nodal attribute selection. Meanwhile, NetworkNet approximates the nodal heterogeneity with the random network model and deep neural networks, which enables the capture of true nonlinear relationships and the accurate selection of the true subset of nodal attributes. As shown in the results in Table 3, NetworkNet gives a dominant performance in both nodal heterogeneity estimation and nodal attribute selection for non-linear settings. 4 Empirical Analysis Academic networks have been a longstanding focus in network science. In the statistics community, the increasing availability of collaboration and citation networks among statis- ticians (Ji & Jin 2016, Ji et al. 2022) has spurred increasing interest in network centrality, community structures, research patterns, and disciplinary trends in the field (Gao et al. 2025, Hayes & Rohe 2025). In this section, we study a related but distinct empirical question, emphasizing the dynamic evolution of influential research fields and their shifting prominence over time. 22 With the powerful specialized framework of NetworkNet, we are equipped to more precisely quantify the scholars’ nodal heterogeneity, using count-valued author-citation networks and incorporating individuals’ research interests as high-dimensional nodal attributes. Our empirical analysis reveals the influential research fields, selected by NetworkNet, that have largely affected the nodal heterogeneity through the expansiveness function and popularity function, respectively. Such nodal heterogeneity subsequently informs the latent interaction intensity that governs the generative mechanism of count-valued citation edges in author–citation networks. 4.1 Dataset and Pre-processing The recently released dataset of large-scale academic network (Gao et al. 2023) consists of information about 97,436 publications authored by 168,171 scholars over a four-decade span in representative statistical journals. In addition to the citation relationships, it records rich nodal information for each scholar. The availability of these nodal attributes enables deep exploration of scholarly interactions through the random networks. Furthermore, by selecting the influential nodal attributes within each decade, the dataset facilitates the analysis of the evolution and scholarly impact of various research topics over time. We first pre-process the raw data extracted from publications and construct the author- citation networks. As shown in Figure 2, we start from the journal publications data where three parts of information are collected: authorship, reference lists, and also keywords. After the data cleaning, we construct the academic author-citation networks by combining the reference lists and the authorship, where the vertices represent different scholars and the weighted edges between them indicate the aggregated citation counts from one scholar to another. As the citation count increases once a new paper is published and the corresponding reference list is available, the edges in author-citation networks are naturally count-valued 23 Figure 2: The construction of the academic author-citation networks with directions from the author to the cited scholar. Meanwhile, for any scholars in the author-citation networks, we aggregate the keywords from all their publications within a time period as high-dimensional nodal attributes, which directly represent the distributions of both their research interests and contributions. Figure 3: Example of count-valued author-citation networks with nodal attributes With the pre-processing and construction steps shown in Figure 2 and by splitting the 24 publication-based raw data between 1981 and 2020 by decade, we construct four decade- specific author-citation networks with count-valued adjacency matrices and high-dimensional nodal attribute matrices. In the subsequent descriptions, we denote them as ‘80-90’, ‘90-00’, ‘00-10’, and ‘10-20’, respectively. Figure 3 illustrates the specific data type associated with each decade, and the basic information about the four author-citation networks is summarized in Table 4. For each of the author-citation networks withnscholars, we use adjacency matrixA ∈N n×n with non-negative integers to summarize the observed network structure, where elementA ij represents the aggregated citation counts of scholarj contributed by the publications written by authori. Nodal research profiles are encoded through a nodal attribute matrixX ∈R n×p , wherepis the number of unique research keywords observed across all scholars, andx ik represents the appearance count ofk-th keyword in authori’s publications. Overall, we observe a rapid growth in the number of scholars and research keywords across the four decades. This expansion signifies the vibrant growth of the statistical discipline, accompanied by a progressive broadening and refinement of its core research areas. Table 4: Network information of the four author-citation networks 80-90 90-00 00-10 10-20 Publication counts11279 18212 25495 38432 Author counts296749988283 14347 Keyword counts393 23761 36400 57927 Edges counts78640 211668 447434 632485 Network density0.009 0.008 0.007 0.003 Average path length 3.290 3.113 3.106 3.336 Diameter12131413 Clustering coefficient 0.185 0.174 0.152 0.116 Average degree26.505 42.351 54.018 44.085 25 4.2 Empirical Findings The directed edge in author-citation networks indicates the citation relationship. Within this framework, the expansiveness functionα(·) describes the inherent characteristics of a scholar to cite the publications by diverse scholars and thereby generating outgoing citations broadly across the community. Conversely, the popularity functionβ(·) captures a scholar’s inherent attractiveness to receive more citations, which is also widely accepted as a measure of academic contributions. To analyze the evolution of influential research fields in statistics, we use the proposed NetworkNet to model the nodal heterogeneity in four decade-specific author-citation net- works, respectively. We select 20 influential nodal attributes for expansiveness function and popularity function, respectively, by tuning the regularization parameters to control model sparsity. After fitting the NetworkNet models, we evaluate the contributions of selected nodal attributes by their SHAP values (Lundberg & Lee 2017), a model-agnostic nodal attribute importance metric widely adopted in the machine learning community (Arenas et al. 2023, Wang, Liang, Hancock & Khoshgoftaar 2024). Since the magnitudes of SHAP values are not comparable across the four author-citation networks, within-network SHAP rankings are used instead. Specifically, attributes are ranked by their SHAP values in each decade, and temporal trends in these ranks are analyzed for major statistical research fields across the four decades. The resulting ranks of all selected nodal attributes are reported in Table 6 for nodal expansiveness and Table 7 for nodal popularity, respectively. The decade-specific SHAP values of influential research fields are shown in Table??-??of the supplementary material. As the selected top 20 keywords are all influential (top 0.1% for the latter three decades), it is worth noting that the repeated selection across the four decades indicates the sustained importance of specific research fields. Table 5 summarizes the consistently selected influential 26 research fields across all four decades, along with their rankings. As shown, bootstrap and consistency emerge as long-standing themes of substantial theoretical and practical scope, characterized by considerable depth and breadth. Bootstrap methods offer great potential for general problems and often provide a tractable way for sufficiently general applications (Diciccio & Romano 1988, Horowitz 2019). Likewise, consistency, as a foundational concept in statistics, always serves as a fundamental part of theoretical guarantees (Lin et al. 2014). Contributions to these fields have demonstrably influenced both the expansiveness and popularity of scholarly work over time. Furthermore, maximum likelihood remains a method of enduring academic relevance due to its foundational role in statistical modeling and its widespread use in machine learning (Kirch et al. 2025) This is also evidenced by its consistent presence among the top 20 popularity-related nodal attributes across all four periods under study. This finding indicates that research contributions specifically focusing on maximum likelihood estimation significantly impact a scholar’s citation attractiveness within the academic author-citation network. Table 5: The ranks of the consistently selected influential research fields in all four decades Attributes80-90 90-00 00-10 10-20 Expansiveness-related Attributes bootstrap2126 consistency77710 Popularity-related Attributes bootstrap2124 consistency66815 maximum likelihood estimation16151910 maximum likelihood1731019 Note: ‘80-90’, ‘90-00’, ‘00-10’ and ‘10-20’ represent the four decades “1981-1990”, “1991-2000”, “2001-2010” and “2011-2020”, respectively. Table 6 shows the rank changes of selected expansiveness-related nodal attributes over the four decades. Since expansiveness characterizes a node’s inherent tendency to connect with other nodes in a network, in the author-citation networks, it represents a scholar’s proclivity to generate outgoing citations broadly across the community. Scholars with higher 27 expansiveness tend to cite publications by a wide and diverse set of scholars, producing more outward edges and larger out-degrees. Meanwhile, research fields with upward trends in the rank indicate a growing influence on scholarly expansiveness that reflects broader citation behavior. This suggests that these fields are progressively evolving into broadly applicable methodological toolkits with strong extensibility and applicability. For instance, the field of Bayesian inference developed rapidly after the success of tackling the computing challenges with MCMC around 1990 (Gelfand & Smith 1990, Brooks et al. 2011), and has become one of the primary areas in statistics (Ji et al. 2022). This trajectory aligns with the results in Table 6 that Bayesian inference was not selected among influential attributes during the 1981–1990 period, but its relative position rose largely after 1990 from rank 18th to 11th. The increasing trend in its rank reflects the deepened research into Bayesian methods and their increased interdisciplinary integration with various fields. Conversely, research fields exhibiting a declining rank over time indicate a reduced role in shaping nodal expansiveness, implying that scholars in these areas increasingly cite work from a narrower or more specialized author set. This trend suggests that such areas are reaching methodological maturity and have formed a relatively stable scholarly community. For example, likelihood ratio test, a long-standing statistical method, was a leading contributor to expansive citation behavior in earlier decades (Moreira 2003, Arnold et al. 2024), supported by its frequent use across diverse models and its strong theoretical properties (Kent 1982, Hogg et al. 1977). Consistent with this, empirical results in Table 6 show a strong expansiveness effect in the 1981–2000 interval. In later years, as its core theory became standardized and incorporated into textbooks, its relative contribution to generating broad citation links declined, becoming non-selected after 2000. Similar downward expansiveness-related rankings are observed for other mature methodologies, such as logistic regression. 28 Table 6: The ranks of the selected influential research fields for nodal expansiveness in the four decades Nodal expansiveness(α)-related fields 80-90 90-00 00-10 10-20 simulation120 bootstrap2126 reliability3 logistic regression4815 mean squared error5 empirical bayes6 consistency77710 poisson regression8 numerical integration9 successive difference analysis10 maximum likelihood estimation11208 repeated measures12 analysis of variance13 likelihood ratio test1415 regression15 kurtosis16 nonparametric regression1723 significance level18 2-way layout19 hierarchical models20 maximum likelihood31219 asymptotic normality443 robustness5615 EM algorithm652 survival analysis91313 missing data101120 Gibbs sampling11 density estimation12 order statistics1310 Markov Chain Monte Carlo1414 smoothing16 gibbs sampler17 Bayesian inference181411 efficiency19 model selection85 longitudinal data914 measurement error16 random effects17 variable selection181 empirical likelihood197 quantile regression9 lasso12 sparsity16 functional data analysis17 causal inference18 Note: ‘80-90’, ‘90-00’, ‘00-10’ and ‘10-20’ represent the four decades “1981-1990”, “1991-2000”, “2001-2010” and “2011-2020”, respectively. 29 Table 7: The ranks of the selected influential research fields for nodal popularity in the four decades Nodal popularity(β)-related fields 80-90 90-00 00-10 10-20 simulation1 bootstrap2124 reliability3 poisson regression4 empirical bayes5 consistency66815 likelihood ratio test718 longitudinal data8116 mean squared error9 successive difference analysis10 analysis of variance11 local optimality12 parameter estimation13 aligned ranks14 normal distribution15 maximum likelihood estimation16151910 maximum likelihood1731019 growth curve18 logistic regression19816 order statistics20147 nonparametric regression23 asymptotic normality443 EM algorithm552 robustness7613 density estimation9 missing data101420 Gibbs sampling11 survival analysis121314 Markov Chain Monte Carlo1315 smoothing16 gibbs sampler17 efficiency19 Bayesian inference20129 model selection97 measurement error15 variable selection171 time series18 random effects20 empirical likelihood8 quantile regression11 lasso12 sparsity16 functional data analysis17 causal inference18 Note: ‘80-90’, ‘90-00’, ‘00-10’ and ‘10-20’ represent the four decades “1981-1990”, “1991-2000”, “2001-2010” and “2011-2020”, respectively. 30 Table 7 presents the selected influential research fields associated with nodal popularity and their ranks over time. As nodal popularity characterizes a node’s inherent characteristics to attract connections from other nodes, in the author-citation network, it directly reflects the perceived influence and citation pull of a scholar’s publications. Scholars with higher popularity typically generate more inward edges and exhibit larger in-degrees. Therefore, research fields with rising rankings indicate increasing contribution to enhancing scholars’ citation pull and collaborative influence. These areas are typically influential topics under- going fast development, marked by broad methodological development, strong theoretical support, and growing interdisciplinary reach. For example, model selection and variable selection gained tremendous momentum in the 2000s with the arrival of the Big Data era (Fan & Lv 2008, Li et al. 2012), and have become increasingly popular with a series of influential works proposed during the early 21st century (Fan & Li 2001, Fan et al. 2020). This trajectory aligns with the results in Table 7 that model selection and variable selection were not selected before 2000, yet they have become prominent as one of the most popular research fields since the beginning of the 21st century, culminating in the first-place ranking for variable selection in 2011–2020. Similar influential popularity-related topics with rising rankings include Markov Chain Monte Carlo, EM algorithm, and Bayesian inference. Likewise, declining rankings among popularity-related research fields indicate that publi- cations in these fields become less appealing for future citations. This pattern potentially implies a state of field maturity, where foundational frameworks are already widely estab- lished and additional publications will have diminishing benefits in terms of the connection in author-citation networks. For instance, analysis of variance (ANOVA) is one of the most classic statistical topics that has been extensively formalized and disseminated, including comprehensive coverage in many standard textbooks (Scheffe 1999, Casella & Berger 2024). As shown in the table, while it attracted active citation engagement in the 1981–1990 period 31 with rank 10th, it was no longer identified as a top contributor to nodal popularity after 1990. Another similar case is the growth curve model, a type of generalized multivariate analysis-of-variance (GMANOVA), which ranked 18th in 1980-1990. While it remains an important topic today, scholarly attention has been shifted toward derivative research directions of GMANOVA, such as survival analysis and longitudinal data analysis (Das 2008). 5 Conclusion and Discussion In this paper, we propose NetworkNet, a deep neural network approach for identifying the complex nodal heterogeneity in count-valued network data with high-dimensional nodal attributes. By combining the advantages of both statistical random network models and deep learning techniques, NetworkNet can effectively model the nodal heterogeneity and simultaneously select the influential subset of nodal attributes using a specially designed DNN with two skip layers. Theoretically, we prove the consistency of estimation for nodal heterogeneity and provide the specific non-asymptotic bound for approximation error. Em- pirically, the dominant performances of NetworkNet in both nodal heterogeneity estimation and nodal attribute selection are shown with comprehensive simulations. By investigating author-citation networks among statisticians with the newly proposed NetworkNet, we select the influential research fields in four decades and provide several insightful findings about the evolution of scientific fields and scholarly impact. While our primary focus has been on count-valued data motivated by author-citation networks, the methodology can be readily extended to other network structures, including binary and continuously weighted networks, by tailoring the likelihood-based loss function. The compelling empirical performance of NetworkNet highlights the value of integrating statistical network modeling with deep learning techniques, and opens promising avenues 32 for future research. For instance, extending the model to signed networks with negative edge weights introduces new challenges, particularly regarding the interpretation of selected attributes; Applying this framework to complex economic systems (such as firm covariates in trade networks or balance sheets in interbank lending) can offer a valuable tool for investigating network formation and systemic risk. Furthermore, exploring the interaction effects among selected influential nodal attributes remains an unresolved challenge. Overall, as a general framework to model complex random network data with high- dimensional nodal attributes via neural networks, NetworkNet provides a powerful tool to capture nodal heterogeneity and select influential nodal attributes. Its flexibility and strong empirical performance make it well-suited for a wide range of interdisciplinary applications in economic, social, and biological network studies. References Abbe, E. (2018), ‘Community detection and stochastic block models: recent developments’, Journal of Machine Learning Research 18(177), 1–86. Antinyan, A., Horváth, G. & Jia, M. (2020), ‘Positional concerns and social network structure: An experiment’, European Economic Review 129, 103547. Arenas, M., Barceló, P., Bertossi, L. & Monet, M. (2023), ‘On the complexity of shap-score- based explanations: Tractability via knowledge compilation and non-approximability results’, Journal of Machine Learning Research 24(63), 1–58. Arnold, B. C., Sengupta, A. & Boca, G. (2024), Likelihood Ratio Tests and Related Topics in Multivariate Analysis, Springer. Brooks, S., Gelman, A., Jones, G. & Meng, X.-L. (2011), Handbook of markov chain monte carlo, CRC press. Casella, G. & Berger, R. (2024), Statistical inference, Chapman and Hall/CRC. Chatterjee, S. & Diaconis, P. (2013), ‘Estimating and understanding exponential random graph models’, The Annals of Statistics p. 2428–2461. 33 Che, C., Chen, Y., Xing, Z. & Zhong, W. (2025), ‘Flat: Fused lasso regression with adaptive minimum spanning tree with applications on thermohaline circulation’, Preprint arXiv: 2507.09800 . Das, S. (2008), Generalized Linear Models and Beyond: An Innovative Approach from Bayesian Perspective, University of Connecticut. Diciccio, T. J. & Romano, J. P. (1988), ‘A review of bootstrap confidence intervals’, Journal of the Royal Statistical Society: Series B (Methodological) 50(3), 338–354. Erdős, P. & R ̋enyi, A. (1959), ‘On random graphs’, Publicationes Mathematicae 6, 290–297. Fan, J. & Li, R. (2001), ‘Variable selection via nonconcave penalized likelihood and its oracle properties’, Journal of the American Statistical Association 96(456), 1348–1360. Fan, J., Li, R., Zhang, C. H. & Zou, H. (2020), Statistical Foundations of Data Science, Chapman and Hall/CRC. Fan, J. & Lv, J. (2008), ‘Sure independence screening for ultrahigh dimensional feature space’, Journal of the Royal Statistical Society Series B: Statistical Methodology 70(5), 849–911. Gao, T., Liu, J., Pan, R. & Sun, A. (2025), ‘BASIC: Bipartite assisted spectral-clustering for identifying communities in large-scale networks’, Preprint arXiv:2503.06889 . Gao, T., Zhang, Y., Pan, R. & Wang, H. (2023), ‘Large-scale multi-layer academic networks derived from statistical publications’, Preprint arXiv:2308.11287 . Gelfand, A. E. & Smith, A. F. (1990), ‘Sampling-based approaches to calculating marginal densities’, Journal of the American Statistical Association 85(410), 398–409. Goodfellow, I. J., Shlens, J. & Szegedy, C. (2015), ‘Explaining and harnessing adversarial examples’, Stat 1050, 20. Graham, B. S. (2017), ‘An econometric model of network formation with degree heterogene- ity’, Econometrica 85(4), 1033–1063. Hayes, A. & Rohe, K. (2025), ‘Co-factor analysis of citation networks’, Journal of Compu- tational and Graphical Statistics 34(2), 448–461. Hogg, R. V., Tanis, E. A. & Zimmerman, D. L. (1977), Probability and Statistical Inference, Vol. 993, Macmillan New York. Horowitz, J. L. (2019), ‘Bootstrap methods in econometrics’, Annual Review of Economics 11, 193–224. 34 Huang, P. & Butts, C. T. (2024), ‘Parameter estimation procedures for exponential-family random graph models on count-valued networks: A comparative simulation study’, Social Networks 76, 51–67. Hunter, D. R., Handcock, M. S., Butts, C. T., Goodreau, S. M. & Morris, M. (2008), ‘ERGM: A package to fit, simulate and diagnose exponential-family models for networks’, Journal of Statistical Software 24(3). Ji, P. & Jin, J. (2016), ‘Coauthorship and citation networks for statisticians’, Annals of Applied Statistics 10(4), 1779 – 1812. Ji, P., Jin, J., Ke, Z. T. & Li, W. (2022), ‘Co-citation and co-authorship networks of statisticians’, Journal of Business & Economic Statistics 40(2), 469–485. Kent, J. T. (1982), ‘Robust properties of likelihood ratio tests’, Biometrika 69(1), 19–27. Kirch, C., Lahiri, S., Binder, H., Brannath, W., Cribben, I., Dette, H., Doebler, P., Feng, O., Gandy, A., Greven, S., Hammer, B., Harmeling, S., Hotz, T., Kauermann, G., Krause, J., Krempl, G., Nieto-Reyes, A., Okhrin, O., Ombao, H. & Lederer, J. (2025), ‘Challenges and opportunities for statistics in the era of data science’, Harvard Data Science Review 7(2). Kohler, M. & Langer, S. (2021), ‘On the rate of convergence of fully connected deep neural network regression estimates’, Annals of Statistics 49(4), 2231–2249. Krivitsky, P. N. (2012), ‘Exponential-family random graph models for valued networks’, Electronic Journal of Statistics 6, 1100. Lehmann, S., Lautrup, B. & Jackson, A. D. (2003), ‘Citation networks in high energy physics’, Physical Review E 68(2), 026113. Lei, J. & Lin, K. Z. (2023), ‘Bias-adjusted spectral clustering in multi-layer stochastic block models’, Journal of the American Statistical Association 118(544), 2433–2445. Lei, J. & Rinaldo, A. (2015), ‘Consistency of spectral clustering in stochastic block models’, Annals of Statistics p. 215–237. Lemhadri, I., Ruan, F., Abraham, L. & Tibshirani, R. (2021), ‘Lassonet: A neural network with feature sparsity’, Journal of Machine Learning Research 22(127), 1–29. Leskovec, J. & Mcauley, J. (2012), Learning to discover social circles in ego networks, in F. Pereira, C. Burges, L. Bottou & K. Weinberger, eds, ‘Advances in Neural Information 35 Processing Systems’, Vol. 25. Li, R., Zhong, W. & Zhu, L. (2012), ‘Feature screening via distance correlation learning’, Journal of the American Statistical Association 107(499), 1129–1139. Lin, X., Genest, C., Banks, D. L., Molenberghs, G., Scott, D. W. & Wang, J. L. (2014), Past, present, and future of statistical science, CRC Press. Liu, R., Shang, Z. & Cheng, G. (2020), ‘On deep instrumental variables estimate’, Preprint arXiv: 2004.14954 p. 1–57. Lundberg, S. M. & Lee, S. I. (2017), ‘A unified approach to interpreting model predictions’, Advances in Neural Information Processing Systems 30. Ma, Y., Lan, W., Leng, C., Li, T. & Wang, H. (2025), ‘Supervised centrality via sparse network influence regression: An application to the 2021 henan floods’ social network’, Annals of Applied Statistics 19(2), 1734–1752. Moody, J. (2004), ‘The structure of a social science collaboration network: Disciplinary cohesion from 1963 to 1999’, American Sociological Review 69(2), 213–238. Moreira, M. J. (2003), ‘A conditional likelihood ratio test for structural models’, Economet- rica 71(4), 1027–1048. Narcowich, F., Ward, J. & Wendland, H. (2005), ‘Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting’, Mathematics of Computation 74(250), 743–763. Pronzato, L. & Zhigljavsky, A. (2023), ‘Quasi-uniform designs with optimal and near-optimal uniformity constant’, Journal of Approximation Theory 294, 105931. Scheffe, H. (1999), The analysis of variance, John Wiley & Sons. Schmidt-Hieber, J. (2020), ‘Nonparametric regression using deep neural networks with relu activation function’, Annals of Statistics 48(4), 1875. Schweinberger, M. (2011), ‘Instability, sensitivity, and degeneracy of discrete exponential families’, Journal of the American Statistical Association 106(496), 1361–1370. Schweinberger, M., Krivitsky, P. N., Butts, C. T. & Stewart, J. R. (2020), ‘Exponential- family models of random graphs’, Statistical Science 35(4), 627–662. Stein, S., Feng, R. & Leng, C. (2025), ‘A sparse beta regression model for network analysis’, Journal of the American Statistical Association 120(550), 1281–1293. 36 Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I. & Fergus, R. (2013), ‘Intriguing properties of neural networks’, Preprint arXiv:1312.6199 . Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society Series B: Statistical Methodology 58(1), 267–288. Vu, M. & Thai, M. T. (2020), ‘Pgm-explainer: Probabilistic graphical model explanations for graph neural networks’, Advances in Neural Information Processing Systems 33, 12225– 12235. Wang, F., Liu, Y., Liu, K., Wang, Y., Medya, S. & Yu, P. S. (2024), ‘Uncertainty in graph neural networks: A survey’, Preprint arXiv:2403.07185 . Wang, H., Liang, Q., Hancock, J. T. & Khoshgoftaar, T. M. (2024), ‘Feature selection strategies: a comparative analysis of shap-value and importance-based methods’, Journal of Big Data 11(1), 44. Wang, J., Cai, X., Niu, X. & Li, R. (2024), ‘Variable selection for high-dimensional nodal attributes in social networks with degree heterogeneity’, Journal of the American Statistical Association 119(546), 1322–1335. Wang, L., Kim, Y. & Li, R. (2013), ‘Calibrating non-convex penalized regression in ultra-high dimension’, Annals of Statistics 41(5), 2505. Watts, D. J. & Strogatz, S. H. (1998), ‘Collective dynamics of ‘small-world’networks’, Nature 393(6684), 440–442. Wendland, H. (2004), Scattered Data Approximation, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press. Xing, Z., Tan, H., Zhong, W. & Shi, L. (2025), ‘Calms: Constrained adaptive lasso with multi-directional signals for latent networks reconstruction’, Neurocomputing p. 129545. Xing, Z., Wang, Y. X. R., Wood, A. T. A. & Zou, T. (2025), ‘Regularization and selection in a directed network model with nodal homophily and nodal effects’. URL: https://arxiv.org/abs/2504.04622 Xing, Z. & Zhong, W. (2024), ‘Palms: Parallel adaptive lasso with multi-directional signals for latent networks reconstruction’, Preprint arXiv:2411.11464 . Yan, T., Jiang, B., Fienberg, S. E. & Leng, C. (2019), ‘Statistical inference in a di- rected network model with covariates’, Journal of the American Statistical Association 37 114(526), 857–868. Yan, T., Leng, C. & Zhu, J. (2016), ‘Asymptotics in directed exponential random graph models with an increasing bi-degree sequence’, Annals of Statistics 44(1), 31–57. Ying, Z., Bourgeois, D., You, J., Zitnik, M. & Leskovec, J. (2019), ‘Gnnexplainer: Generating explanations for graph neural networks’, Advances in Neural Information Processing Systems 32. Yuan, Y. & Qu, A. (2023), ‘High-order joint embedding for multi-level link prediction’, Journal of the American Statistical Association 118(543), 1692–1706. Zhang, L. & Amini, A. A. (2023), ‘Adjusted chi-square test for degree-corrected block models’, Annals of Statistics 51(6), 2366–2385. Zügner, D., Akbarnejad, A. & Günnemann, S. (2018), Adversarial attacks on neural networks for graph data, in ‘Proceedings of International Conference on Knowledge Discovery & Data Mining’, p. 2847–2856. 38

Instant research discovery