A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Authors: Yufei Huang
Changhu Wang
Junjie Tang
Weichi Wu
Ruibin Xi
Venue: ICML 2025 poster (2025)
Keywords: High-dimensional graphical model
Maximum marginal likelihood estimation
Marginal recoverability
Mixture model
Gene regulatory network
Single-cell RNA sequencing
Primary area: probabilistic_methods->graphical_models

Abstract: Traditional network inference methods, such as Gaussian Graphical Models, which are built on continuity and homogeneity, face challenges when modeling discrete data and heterogeneous frameworks. Furthermore, under high-dimensionality, the parameter estimation of such models can be hindered by the notorious intractability of high-dimensional integrals. In this paper, we introduce a new and flexible device for graphical models, which accommodates diverse data types, including Gaussian, Poisson log-normal, and latent Gaussian copula models. The new device is driven by a new marginally recoverable parametric family, which can be effectively estimated without evaluating the high-dimensional integration in high-dimensional settings thanks to the marginal recoverability. We further introduce a mixture of marginally recoverable models to capture ubiquitous heterogeneous structures. We show the validity of the desirable properties of the models and the effective estimation methods, and demonstrate their advantages over the state-of-the-art network inference methods via extensive simulation studies and a gene regulatory network analysis of real single-cell RNA sequencing data.

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Yufei Huang" 1
Changhu Wang 23 Junjie Tang? Weichi Wu Ruibin X123
Abstract
ingly focused on high-dimensional settings, commonly en- countered in applications such as microarray experiments (Amaratunga et al., 2014) and single-cell RNA sequencing (scRNA-seq) studies. However, parameter estimation in such settings poses significant challenges due to the compu- tational intractability of nigh-dimensional integrals, limiting the efficiency of existing methods. This highlights the crit- ical need for developing efficient estimation frameworks tailored to high-dimensional graphical models.
Traditional network inference methods, such as Gaussian Graphical Models, which are built continuity and homogeneity, face challenges when modeling discrete data and heteroge- neous frameworks. Furthermore, under high- dimensionality, the parameter estimation of such models can be hindered by the notorious in- tractability of high-dimensional integrals. In this paper, we introduce new and flexible device for graphical models, which accommodates di- verse data types, including Gaussian, Poisson log-normal, and latent Gaussian copula models. The new device is driven by a new marginally recoverable parametric family, which can be ef- fectively estimated without evaluating the high- dimensional integration in high-dimensional set- tings thanks to the marginal recoverability. We further introduce mixture of marginally recover- able models to capture ubiquitous heterogencous structures. We show the validity of the desirable properties of the models and the effective esti- mation methods, and demonstrate their advan- tages over the tfateoof-the-ar network inference methods via extensive simulation studies and a gene regulatory network analysis of real single- cell RNA sequencing data.
Due to well-established theoretical properties and high inter- pretability, Gaussian graphical model (GGM) (Meinshausen & Bühlmann, 2006) has been extensively studied recent years. Notable methods include -i-pelalized log likelihood maximization (Yuan & Lin, 2007; Banerjee et al., 2008; Friedman et al., 2008), penalized regression (Meinshausen & Bühlmann, 2006; Peng et al., 2009), adaptive thresholding (Cai & Liu, 201 1), and D-trace loss (a smooth convex loss function) (Zhang & Zou, 2014), among others. To address the limitations the Gaussian assumption, semiparametric Gaussian copula models (Liu al., 2009; Xue & Zou, 2012: Liu al., 2012) have been developed to handle continuous data, using monotonic univariate transformations. These methods, based on Gaussian and Gaussian copula models, are usually efficient but are only applicable to continuous data and cannot handle discrete data.
In many applications, particularly in genomics studies, dis- crete data are ubiquitous. The observed data are considered to be generated from the discretization of underlying la- tent variables (Skrondal & Rabe-Hesketh, 2007). which naturally leads to the use of hierarchical models for data modeling. For example, in scRNA-seq data, gene expres- sion levels are count data, often featuring numerous zero values (Islam et al., 2014; Zheng et al., 2017). To effectively model the underlying dependencies in count data, the mul- tivariate Poisson log-normal (PLN) distribution (Aitchison & Ho, 1989) is plausible and popular, since it can capture the conditional dependencies while accommodating over- dispersion in the marginal distributions. However, compared to GGM, maximizing the likelihood of the PLNN model is more challenging due to the high-dimensionall integrals in- volved, without known closed-forn ffoem solution. To address this issue, Choi et al. (2017) used Laplace's 's method for
1. Introduction
Graphical models (Lauriizen,, 1996) are widely used to ex- plore network structures and identify complex interactions between random variables. Recent research has increas-
"Equal contribution 'Center for Data Seience, Peking Univer- sity, Beijing, China "School of Mathematical Sciences, Peking University, Beijing, China 'Center for Statistical Science, Peking University. Beijing, China 'Department of Statistics and Data Sci- ence, Tsinghua University, Beijing, China. Correspondence to: Ruibin Xi <ruibinxie @ math, math.pku.edu.comm pku.edu.en>, Weichi Wu <wuwe- ichi @ www.:tsiinngnuu....om
Proceedings of the 4nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
likelihood approximation and applied Alternating Direction Method of Multipliers (ADMM) to compute the penalized maximum likelihood estimator. Chiquet et (2.19) in- troduced a variational approximation method to infer the network of the PLN model. However, these methods rely on approximations of the likelihood function and lack theoreti- cal guarantees.
networks, showeasing its practical utility. All code is avail- able at http::::t//... http//gtttuuupt wtt.c/oo///iii. htmmm/iiDaaa//// XiD.stababb//b............... EMMMLE.
Our work balances diversity, efficiency, and heterogeneity, whereas existing methods typically address at most two of these aspects. For instance, Fan et al. (2017) proposed a latent Gaussian copula model for mixed data, combining continuous and binary variables, and introduced a gener- alized rank-based method. While this approach ensures diversity and efficiency, it does not account for heterogene- ity as it cannot handle mixture models.
In practice, samples often arise from mixed populations exhibiting heterogeneous patterns. Traditional graphical models usually assume that samples come from single population with a shared network, which limits their ability to capture data heterogeneity. To address this, mixture mod- els are widely utilized. For example, in scRNA- data, samples comprise single cells from different cell types, each with its own gene regulatory network. One of the arguably most important tools for studying the mixture and count data with different reglatory networks is the mixture Poisson log- normal (MPLN). Silva et al.(2019) employed the MPLN model for clustering count data and used Expectation- Maximization (EM) algorithms with MCMC steps to maxi- mize the computationally intractable oogllikelihood. Tang et al. a2.24) adopted a variational inference approach for inferring cell-type-specifii gene regulatory networks. How- ever, both EM algorithms with MCMC and variational in- ference are computationally expensive.
2. Flexible Marginally Recoverable Family
In this section, we propose the marginally recoverable para- metric family to address diversity. First, we introduce some notations. For p-dimensional vector with the i-th entry ag [[j = (a,, az) represent the two-dimensional sub- vector for S. j<k S p. For a p X p symmetric matrix A with the (i. entry oij, the submatrix Ajj, is defined as:
Ojj ajk
azk
(1)
Okk
where 1 j S p. This submatrix is constructed by selecting the entries in the j-th and k-th rows and columns of A. Let M, denote the set of p x p symmetric positive semi-definite matrices.
To address the aforementioned challenges, this paper makes the following contributions with theoretical guarantees.
Our work is inspired by the properties of the Gaussian dis- tribution (Lauritzen, 1996). For a p-dimensional random vector X (1,, E), any of its two-dimensional marginal distributions X(j.A] (2 Eg,a) for 1 S < k S p. The parameters of the p-dimensional distribution can be fully characterized by the parameters of all two-dimensional marginal distributions. Consequently, the model parameters can be estimated independently through marginal likeli- hoods. Inspired by this property, we define the marginally recoverable parametric family as follows.
1) To accommodate diverse data types, we propose a novel class of distributions, the marginally recoverable para- metric family, which unifies existing models and is not restricted to single data type. This family includes the Gaussian, PLN, and latent Gaussian copula models introduced by Fan et al. (2017) for binary data.
2) To overcome the computational intractability of high- dimensional integration, we develop an efficient param- eter estimation framework for this family based on the Maximum Marginal Likelihood Estimator (MMLE). Thanks to marginal recoverabilityy this method sim- plifies nigh-dimensional integrals into multiple low- dimensional integral computations.
Definition 2.1 (Marginally recoverable parametric family). Let {Ha] be a sequence of distribution functions, where for each d 2 1, any p-dimensional marginal distribution (1 S p S d) of Ha belongs to the family
Hpp = (Hp Σ) μ E E Mp).
3) To capture heterogeneous structures, we extend the framework to a mixture of marginal recoverable mod- els by integrating the EM algorithm to update the MMLE to EM-MMLE.
For any p-dimensional random vector X = (XI...., Xp) with p 2 2 such that X Hp(m, Σ) where P (ui,..., :pp) TT and Σ = [ojn
we say that X or H, is marginally recoverable if the following conditions hold:
We establish the consistency of MMLE for covariance matri- ces and networks under mild conditions. Simulations show our method outperforms existing ones for mixture count and binary data. Furthermore, we apply the EM-MMLE to real scRNA-seg seq data to infer seleeyfp--peecffi gene regulatory
. For 1 $ X, ~ Hı Tzz).
For 1sj k S p. X [J,1 ~ H2( E(g.)))
2

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
The marginally recoverable parametric family includes many of the most common distributions. For instance, we have the following remark.
prominent example is DNA nucleotide data (Abbasy et al., 2012). More concretely, genes that exhibit higher levels of expression are represented as 1, whereas genes with lower levels of expression are represented as O. The latent Gaus- sian copula model for binary data is proposed by Fan et al. (2017).
Remark 2.2. Elliptical distributions, including the Gaus- sian distribution and the multivariate t-distribution, are marginally recoverable.
The following proposition elucidates that hierarchical mod- els with marginally recoverable inner layers also satisfy the marginally recoverable condition in Definition 2.1.
Example 2.6 (Latent Gaussian copula model for binary data). Let X = (XI,..., be a random p-vector and Y= (Y,,..,, ..))))))) represents p-dimensional binary vari- ables. Y is sampled from the latent Gaussian copula model for binary data, if
Let Q()) be a distribution function char- acterized by a single parameter A. If is marginally re- coverable and Y X ([==Q((X)) then Y is also marginally recoverable.
Proposition 2.3.
Y, 1 > Cj),
(3)
Gaussian copula a widely used semiparametric model, overcoming the drawback of the Gaussian model's reliance on exact normality (Liu et al., 2009). It is easy to verify that the Gaussian copula model is marginally recoverable.
X ~ NPN(O, f),
where Cj is a constant, I(-) is the indicator functionn and Jj for any 1 1 p.
Definition 2.4 (Gaussian copula model). Let X be ran- dom p-vector. X is sampled from the Gaussian copula model, if there exists monotonic transformation f such that S(X)= = (Xi) (X.)...., (Xp)) N2(,, Then we denote X ~ NPN(H,E,f)
In the latent Gaussian copula model for binary data, each binary component Yj, which takes values of 0 or 1, is gener- ated from a latent continuous random variable X, truncated at an unknown threshold value Cj.
Despite its flexibility, the Gaussian copula model cannot be directly applied to discrete data. For practical applications in network inference, we introduce two hierarchical models designed for count data, including binary data, both which are marginally recoverable according to Proposition 2.3.
Examples 2.5 and 2.6 can be regarded as two specific in- stances of marginally recoverable distributions. In both cases, Σ represents the covariance matrix of latent vari- ables. In Example 2.5, μ4 is the mean of the Gaussian copula distributions in the inner layer, while in Example 2.6, μ rep- resents the threshold values. The inverse of the covariance matrix Σ, denoted as θ reveals the network. Specifically, in Example 2.5 and 2.6, X; and X3 are independent given the remaining variables if and only if O, 0. Conse- quently, inferring the graph structure can be achieved by estimating O However, the likelihood function in hierarchi- cal models, such as Example 2.5 and Example 2.6, involves high-dimensional integrations that pose significant compu- tational challenges, since these integrations seldom have closed-form solutions, making the computation of max- imum likelihood estimation infeasible. Fortunately, the computational issue can be circumvented by leveraging the properties of marginally recoverable distributions and sim- plifying the problem of estimating the parameters of a high- dimensional distribution into multiple lower-dimensional parameter estimation problems.
Example 2.5 (Latent Gaussian copula model for count data). Let X = (XI,...,X) Xp) and Y = (Y1...., .p)) be two random p-vectors. Y is sampled from the latent Gaussian copula model for count data, if
Y X Poisson (S exp (X,)),
1=1 X ~ NPN(H, E.
This distribution is commonly used to model genomic data (Sarkar & Stephens, 2021; Sinclair & Hooker, 2019). When X ~ (p,, E), we say that Y follows the PLN distribution, denoted as PTNN (S: μ. Σ). The PLNN model is widely used for single-cell RNA sequencing data, influenza-like illness dataset and purchase frequency counts (Silva et al., 2019; 9; Wu et al., 2018: 8: Trinh al., 2014). In scRNA-seq data, X denotes the underlying expression levels of genes and S represents the sequencing depth of the cell, which can be estimated by the sum of UMI counts across all genes (Sarkar& & Stephens, 2021; Hafemeister & Satija, 2019).
3. Efficient Estimation
In this section, we introduce an efficient estimation frame- work designed to estimate parameters, particularly Σ, which is associated with the network structure in H,, E). Then, we extend the estimation framework to accommodate mix- ture distributions for inferring networks.
The binary data type is an important special class of count data and is often observed in genetic and genomic studies. A
3

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
3.1. The Maximum Marginal Likelihood Estimator
that the distribution Hp is marginally recoverable. A G- component mixture of marginally recoverable distributions can be expressed as H (T,, 2t) = ->g To Tg Eg) where ! = [ MG: Ei.... Eg}.
Let hyp denote the density function of the distribution Hp- Suppose that Y;, for i ,,....;7 n, are random p- dimensional vectors sampled from Hp (He, Σ). When p is large, computing the maximum likelihood estimator is in- feasible due to the intractable high-dimensional integration. By leveraging the properties of the marginally recoverable component as defined in Definition 2.1, we propose an effi- cient estimation framework based on MMLE for parameter estimation. This framework is inspired by the pairwise like- lihood method (Cox& & Reid, 2004; Varin, 2008), which com- putes only the likelihoods of pairs of observations. There- fore, this method significantly reduces the computational cost compared to conventional likelihood. This reduction in computation relies on evaluating a limited number of sets of two -ro-dimensional integrals instead of computing the full high-dimensional likelihood integral, i.e.,
However, due to the unknown sample categories, maximiz- ing the marginal ligllikelhhood of mixed distributions is computationally intractable practice. Therefore, we pro- pose the EM algorithm to update the MMLE for data from mixed populations and name the estimator EM-MMLE.
Assume that Y;, for i = 1. Rr are p-dimensional random vectors generated from HP 12). To indicate cluster membership, we introduce the indicator variable Z, for the i-th sample, which follows a multinomial distri- bution with proportion parameter π.. For S < k S p. let Øy { TT. Mgss.A. E g1b 9== = 1, ...,, G. containing the unknown parameters of Y (.... In the 1-th iteration, the ex-
pected complete log-likelihood Q(0,a, (((-1) is taken as the optimization objective, which is defined as
p-1 1 d
Spair(p. E:y) = II II -=1-=
h2(us, be( M.N) Di (4)
Q jk, ( = E
p(Z vo "y Dogr(Y (Y, 2z0xx)
Unfortunately, the pairwise maximum likelihood estimator (PMLE) is often inconsistent (Varin et al., 2011). However, for marginally recoverable distributions, the PMLE is con- sistent because we can estimate the parameters from the two-dimensional marginals.
We update (1) = arg maxe,, Q( (1 The itera- jk tion terminates when the change in the optimization objec- tive between consecutive steps is negligible. The framework of EM-MMLE is summarized in Algorithm 1.
Note that when μa and تأزوه .P) p) are known, maximizing Lpai E; is equivalent to maximizing h22 h255, M&s M[j.A) ()) for 1 S j<k S p. Motivated by this and to avoid redundant computation of μ and دزj- we first estimate Hjy and وزj by maximizing the one- dimensional marginal log-likelioood. Then, Jys is estimated by maximizing the two-dimensional marginal log- likelihoo.. The MMLE Σ- (ō is derived follows:
A key application of the mixture distributions HT( is to infer networks from its parameters { (gg 9 1...., G}. Therefore, with EM-MMLE ES (g ... G), we apply the D-trace method to estimate the sparse precision matrix + Σ g i.c.
θ = arg min tr (0 tr + +A11 e 2
θ 11l off
(6)
where θ - indicates that 0, is positive semi-definite, and 00 = Σijz @gijl This optimization problem is efficiently solved using an alternating direction method as described in Zhang & Zou (2014). To ensure the convexity of the loss function, the D-trace approach requires the input covariance matrix estimator to be positive semi-definite. To satisfy this requirement, we propose the Eg:
رزj = arg max 0z
log hi (Ysji My 03))
(5)
oi4 = arg max
log h2
Y Yifjiaii (1j E13.0)..
where in the optimization for dzk: Jj and Okk in E [iM are fixed to وزو and δb.. respectively.
ΣS Σn + 1| - E,
3.2. Mixture for Heterogeneity
where I is the identity matrix and ΣS = arg min ||
The mixture model surpasses the limitations of the single- component model and offers additional flexibility mod- eling complex data. In the following, we focus on the mix- ture model within the marginally recoverable family. Let (M...... ) denote the mixing proportions, where 1= Mg 1 and Fg > for g = 1,...,C G. Suppose
Eg boo is the projection of Σn onto the space of positive semi-definite matrices, which can be solved by a splitting conic solver (Fu et al., 2020).
The selection of the tuning parameter in Equation (6) is achieved by minimizing the approximate Bayesian informa-
4

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Algorithm 1 Framework of EM-MMLE
the necessary conditions.
Input: Observed data YI,.., ,, the number of pop- ulations G, the maximum iteration number T and the convergence threshold e1.
Let 41 (y; h (y: represent the one-dimensional and two-dimensional marginal density functions of Hp (Z), respectively, with parameters W1 E O1 and W2 E O2, where O1, O2 are parameter spaces in finite- dimensional Euclidean space, The Hellinger distance be- tween two densities pl and p2 is defined as: d(pi.P2)
Output: E,,... ....... Σg.
for j = 1 to p- 1 do
for k = j + to p do
while t < T and I. > ez do
S (p) - P2 1/22 dv 1/2 =||01 1/2 - P2
for i = to n do
1/2
for g = 1 to G do
P(0) = h22
Y Mi
(1-1) 11--D) g[ink] [
Condition 4.1 (Lower boundedness condition). For k = 1,2, and for any Wx E Ob there ex- ists a positive constant e such that c||Wx= - w/ll2 S ( (y; hid (y;wk)).
x(0)) gifk
1) pic)
Σ 11
end for end for
Condition 4.2 (Upper boundedness condition). For k 1,2, and for any Wk, w's E Ork there exist measurable function m(y) and a constant e. succ that f m² my))dv h/ 00. and I (y;wk) - hhk
for g = 1 to G do
A n Di=1
(0)
losk
²²
(y: (2r
S
Update Pai and a and by maximizing
m(y) ||W> - w'1l12
(t) Ei=1 'gi8 log h2 (Viz: pggs ggs))-
Update Ag and ((1)
Based on these boundedness conditions, we establish a the- orem that elucidates the convergence rate of MMLE Σ.
glik by maximizing
200) EL= 'oik
log h1 (Yuxi Myke s.k)).
Theorem 4.3 (Rate of convergence from MMLE Σg), For the mixture marginally recoverable model, assume that Conditions 4.1 and 4.2 hold. Then, there exists a con- stant such that for any 1 S S G and > 0,
(1) gsk
Update
by
maximizing (t) "biA) Egg..))-
(1))) Ei=1 giik
log h2(2)/2
end for
8 (Q(4:/0 (0) jk: (1))) Q(0%/ ",00""000000 o0-"))
pr ( [ E. l 1 26) S 5Gpp exp (-ane?).
where 8(z,8) = 1x - yl/y.
Theorem 4.3 shows that if p < exp(én) for some constant c, or in other words, if p tends infinity not faster than exponential of n, then Eg is consistent estimator of Eg-
end while end for end for
The precision matrix θ g ΣØ¹ encodes the network struc- ture. Specifically, an edge exists between vertices i and j in the g-th group if O * 0; otherwise, Ogi = 0 indicates edge. With the rate of convergence for Eg. apply the Theorem Xiao et al. (2022) to each Σ, and then estabilsh the sign consistency of the estimator θg.
tion criterion
-I F
+ +(1gglo logn) /n.
(8)
2
where |/0ylo denotes the number of nonzero entries in θ
Theorem 4.4 (Sign consistency for θ₂). Assume that for each g I....,G G, the true precision matrix O, satisfies the conditions described in the Appendix A.2. Then, for some > 2, pr (vec(Og)a; =0) > 1 p2- if is sufficiently large (depending on 1J, see Appendix A.2), where vec( +gy) denotes the vector formed by stacking the columns of go and veel @ @g)cg represents the subvector indexed by g ((i.j) 0gr = 0).
4. Theoretical Properties
In this section, we establish theoretical results concerning the convergence rate of MMLE and the application to net- work estimation. We focus on the mixture of marginally recoverable distributions H.M (IT, 2),, which reduces to the single-component model discussed in Section 3.1 when G = 1. In theory, under the assumptions that the mixing proportions Ty and means M4 (g- = I,....C G) are known, we can estimate (Eg.9 = 1,.... G} gg 9 I.....G! G} using the MMLE. Before presenting the theoretical results, we first provide
Theorem 4.4 demonstrates that the estimator θ recovers all zeros and nonzeros in O, with probability 1 2-p2-9 The proof of Theorem 4.4 is similar to that of Theorem 3 in Xiao et al. (2022) and is omitted here.
Notably, many mixture marginally recoverable distributions
5

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
satisfy both Conditions 4.1 and 4.2. A crucial example the MPLN distribution, well-suited for gene regulatory network inference in scRNA-seq studies. Unlike single-model ap- proaches requiring prior knowledge of cell type labels, the MPIN model handles network inference without such prior knowledge. Additionally, it accounts for overdispersion in scRNA-seq data and supports both positive and negative correlations (Silva et al., 2019; Tunaru, 2002). Tang al. (2024) proposed VMPLN, a variational inference-basee al- gorithm for estimating the precision matrices of MPLN, but it is time-consuming and lacks theoretical guarantees. In contrast, our estimation method is supported by theoretical guarantees. Under mild conditions, we show the following:
Frequents
ENENL-MLLE
-10
VYPPIN
1) The MPLN distribution satisfies Conditions 4.1 and 4.2. This nontrivial proof is provided in Appendix A.3.4.
Figure I. Mean networks predicted by EM-MMLE, PLNet, VM- PLN, and Glasso for the banded graph with p 100, 3000, high mixing, and a high zero-pras pero-proportion rate. False edges are colored in red and true edges are in dark blue.
2) The binary data model in Example 2.6 also satisfies both conditions, with the proof detailed in Appendix A.4.
are provided in Tables 3-5 in the Appendix B.3. As ex- pected, the AUPR decreases with an increase in the number of genes or higher zero-proportion levels. Among the four methods, EM-MMLE proves to be the most robust, consis- tently outperforming PINet, VMPLN, and Glasso in AUPR across all simulated scenarios. This advantage is particu- larly evident in high-dimensional settings, situations with higher population mixing levels, when zero-proportion rates are high. As the sample size increases, the AUPR of EM-MMLE improves significantly, aligning well with theoretical expectations.
Thus, Theorem 4.3 and Theorem 4.4 demonstrate broad applicability.
5. Simulation
To evaluate the performance of our method, we conduct simulations on mixed count data and binary data.
5.1. Simulation for Mixed Count Data
We generate simulation data following the MPLN distribu- tion and compare EM-MMLE with the available network inference methods including PLNet (Xiao et al., 2022), VM- PLN (Tang al., 2024), and Glasso (Friedman et al., 2008). EM-MMLE and VMPLN are developed to directly esti- mate the precision matrices from the MPLN model, using K-means clustering results as initial values. For PLNet and Glasso, samples are clustered using the K-means algorithm, and then the network is inferred separately for each cluster.
Moreover, the superiority of EM-MMLE over PLNet be- comes more evident with increased population mixing, For example, in simulations with blocked random graph (n 3000, p 300, low pero-proportion), EM-MMLE achieves mean AUPRs of 0.84 and 0.94 in high and low mixing scenarios, respectively, approximatelyy 15% and 8% higher than PLNet's AUPRs (0.73 and 0.87).
Additionally, compared to VMPLN, another network infer- ence method based on the MPLN model, EM-MMLE shows superior performance, especially in higher-dimensional set- tings. For instance, in random graph simulations (n 3000, low mixing, low zero- pero-proportion), EM-MMLE achieves mean AUPRs of 0.96 (p = 300) and 0.99 (p = 100), outper- forming VMPLN by 28% and 4%, respectively.
Simulation Settings. The number of populations G is set 3, and the proportion parameter и is set as (1/3,1/3, 1/3). The synthetic datasets vary across different network struc- tures (random, blocked random, banded, and scale-free), dimensions (p= 100,300), sample sizes (n = 1800, 3000), population mixing levels (low, middle, and high), and zero- proportion levels that represent the proportion of zeros in the count matrix (low and high). In each scenario, we inde- pendently repeat simulations 50 times. Details of the data generation process are provided in Appendix B.1.
To further assess the performance of EM-MMLE, we vi- sualize the average network predicted by the four methods across 50 simulations. We computed the relative frequency matrix F to capture the accuracy of network recovery. For Og f 0, Flj represents the proportion of simulations in which the edge was correctly recovered. Conversely, when Oi = 0, Fly denotes the negative proportion of simulations in which an edge between nodes and j was incorrectly
Performance Comparison. Table 1 shows the area un- der the precision-recal curve (AUPR) of each estimator for random graphs, calculated by varying the penalty pa- rameters, while results for the other three graph structures
6

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Table 1. Comparisons of FM-MMI.E with Pl,Net, VMPLN and Glasso in terms of AUPR on simulation results for random graphs generated by the MPLN model. The results are averages over 50 replicates with standard deviations in brackets.
Zero-proportion
Low
Dimension
p = 100 Middle
= 300 Middle
Mixing degree
Low
High
Low
High
1800
EM-MMLE
0.95 (0.013)
0.94 (0.012)
0.91 (0.018)
0.86 (0.024)
0.81 (0.022)
0.72 (0.04)
PLNet
0.89 (0.03)
0.86 (0.046)
0.81 (0.047)
0.74 (0.061)
0.67 (0.066)
0.57 (0.085)
VMPLN
0.9 (0.023)
0.9 (0.022)
0.89 (0.037)
0.67 (0.026)
0.66 (0.017)
0.64 (0.022)
Glasso
0.83 (0.036)
0.8 (0.032)
0.76 (0.047)
0.72 (0.028)
0.69 (0.027)
0.61 (0.034)
3000
EM-MMLE
0.99 (0.006)
0.98 (0.006)
0.98 (0.008)
0.96 (0.008)
0.94 (0.011)
0.87 (0.019)
PLNet
0.95 (0.047)
0.95 (0.033)
0.93 (0.03)
0.91 (0.052)
0.87 (0.063)
0.75 (0.052)
VMPLN
0.95 (0.028)
0.95 (0.021)
0.94 (0.016)
0.75 (0.021)
0.74 (0.013)
0.72 (0.019)
Glasso
0.89 (0.04)
0.88 (0.033)
0.85 (0.035)
0.82 (0.021)
0.79 (0.019)
0.72 (0.024)
propootoootion
High
Dimension Mixing degree
p = 100 Middle
= 300 Middle
Low
High
Low
High
1800
EM-MMLE
0.83 (0.026)
0.8 (0.037)
0.75 (0.03)
0.62 (0.041)
0.55 (0.028)
0.49 (0.03)
PLNet
0.72 (0.073)
0.69 (0.059)
0.61 (0.069)
0.55 (0.064)
0.47 (0.04)
0.41 (0.031)
VMPLN
0.77 (0.037)
0,75 (0.034)
0.72 (0.033)
0.47 (0.034)
0.45 (0.025)
0.44 (0.021)
Glasso
0.55 (0.045)
0.49 (0.046)
0.45 (0.045)
0,47 (0.032)
0,44 (0.028)
0,41 (0.026)
3000
EM-MMLE
0.95 (0.01)
0.94 (0.013)
0.92 (0.017)
0.83 (0.028)
0.79 (0.025)
0.71 (0.037)
PLNet
0.89 (0.037)
0.87 (0.048)
0.81 (0.065)
0.77 (0.048)
0.71 (0.039)
0.62 (0.048)
VMPLN
0,86 (0.038)
0.84 (0.042)
0.81 (0.047)
0.56 (0.031)
0.55 (0.026)
0.54 (0.018)
Glasso
0.63 (0.04)
0.6 (0.06)
0.55 (0.063)
0.58 (0.029)
0.58 (0.033)
0.53 (0.043)
predicted. Figure 2 shows the results for the banded graph with p = 100, n = 3000, high mixing, and a high zero- proportion rate. EM-MMLE closely aligns with the true net- work structure, detecting more true edges while maintaining the lowest false positive rate compared other methods.
its input, where ΣA = and X = A Σ. X,
1 Σ (X, -X) (x, xX
Simulation Settings. We set p = 50, 200 and evaluate the performance for three sample sizes: n = 200, 500, 3000. Each simulation scenario is repeated 100 times. The data generative process detailed in Appendix B.2
5.2. Simulation for Binary Data
Performance Comparison. Table 2 presents the average estimation errors of Σ - Σ and θ O as measured by the Frobenius norm. It is seen that the estimation errors of MMLE and Fan's 's method remain nearly identical across different dimensions and both are significantly lower than those of AMLE. When n is small, MMLE demonstrates higher accuracy in estimating Σ compared to Fan's 's method. Notably, Fan's method can only handle model (3), whereas MMLE offers greater generalizability. Compared to the benchmark Oracle estimator, the results in Table 2 indicates that MMLE has almost no information loss at n = 200 and n = 500, and only moderate information loss in the high-dimensional environment at n 3000.
To evaluate the performance of MMLE in estimating the inverse correlation matrix + for binary data, we follow similar data-generating procedure as described in Fan al. (2017). Specifically, we generate simulation data Y = (Yi,.... YD)", where Y3- T (X; C2) for j ... with X (0, E).
We then compare the performance of MMLE with three estimation methods: the rank-based estimator by Fan et al.(2017). AMLE (d' Aspremont et al., 2008), and an Ora- cle estimator, which utilizes the Pearson correlation of the latent variable X within the D-trace loss to benchmark the information loss of estimators derived from observed data Y. The rank-based estimator by Fan et al.(2017) is designed to estimate the correlation matrix in latent Gaussian copula model. To estimate the precision matrix, we apply D-trace to the correlation matrix from this estimator to assess its accu- racy in estimating O AMLE is graphical lasso estimator that takes the modified sample covariance matrix ΣA as
6. Application to a scRNA-seq Dataset
In this section, we evaluate the performance of EM-MMLE for gene regulatory network inference using real scRNA- seq dataset (Zheng et al., 2017), comprising 6,952 cells
7

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Table Average estimation MML.,, Fan et al.'s method, Oracle, and AML. measured by the Frohenius binary data. The results are averages over 100 replicates with standard deviations in brackets.
p
50
200
T0
200
500
3000
200
500
3000
Σ Σ
Oracle
3.599 (0.02)
2.173 (0.01)
0.909 (0.01)
14.153 (0.03)
8.911 (0.02)
3.629 (0.01)
MMLE
6.343 (0.15)
3.978 (0.07)
1.597 (0.03)
25.84 (0.14)
15.869 (0.11)
6.529 (0.03)
Fan al.
6.375 (0.15)
3.986 (0.07)
1.598 (0.03)
25.978 (0. 14)
15.901 (0.11)
6,531 (0.03)
AMLE
31.254 (0.11)
31.247 (0.06)
31.036 (0.03)
125.278 (0.31)
123.72 (0.36)
124.367 (0.08)
O
Oracle
2.327 (0.03)
2.335 (0.03)
1.12 (0.08)
2.365 (0.02)
2.342 (0.03)
1.045 (0.03)
MMLE
2.328 (0.03)
2.334 (0.03)
1.335 (0.06)
2.365 (0.02)
2.363 (0.03)
2.332 (0.04)
Fan al.
2.327 (0.03)
2.333 (0.03)
1.336 (0.09)
2.365 (0.02)
2.363 (0.03)
2.33 (0.04)
AMLE
2.684 (0.34)
2.834 (0.03)
2.198 (0.03)
2.365 (0.02)
3.564 (0.02)
2.871 (0.04)
heatmap of AUPRC ratios for six methods (four evaluated in the simulation and two additional single-cell gene regula- tory network inference methods: PPCOR (Kim, 2015) and GENIE3 (Huynh-Thu al., 2010)). EM-MMLE achieves the highest AUPRC ratios in most cases and consistently outperforms the random predictor across all cell types.
EM MMLE
14.6
9.6
3.9
3.4
VMPLN
2.1
5.1
7.6
AUPRC ratio High
PLNet
8.9
1.3
3.5
1.9
From a biological perspective, some interesting association patternss are identified by EM-MMLE (Figure 3). It reveals an association between the genes MYBL2 and TK1, which is predicted across the regulatory networks of three cell types. This finding is supported by the literature (Qiu et al., 2022), which links these genes N glycan biosynthesis and p53 signaling pathways. Additionally, EM-MMLE predicts an association between the genes MYBL2 and BIRC5, sup- ported by et al. (2021). Interestingly. we also identified cell-type-specific hub gene ID3 in the gene regulatory network of CD4+ T cells. Gene ontology analysis of ID3 target genes (Figure 4) reveals enrichment terms related to antigen processing and presentation, including exogenous and endogenous peptide antigens via MHC class II. The rela- tionship between CD4+ T cells and MHC class II is central to adaptive immune responses in the immune system.
PPCOR
1.3
16.2
Glasso
4.7
Low
GENIE3
2.4
B cells
Monocytes
CDB+T cells CD4+T cells
Figure 2. AUPRC ratios of network inference algorithms on the scRNA-seg dataset. Algorithms ordered decreasing median AUPRC ratios. The color in each cell is proportional to the corre- sponding value (scaled between 0 and with values below those of a random predictor shown as grey squares). The actual values are displayed the boxes.
across four cell types. The dataset includes batches, sequenced by 3' and 5' scRNA-seq technologies. One batch is used to construct silver standards based on public regu- latory network databases (Appendix C.1), while the other batch is reserved for algorithm testing. We infer regulatory interactions for the top 300 highly variable genes identi- fied by Seurat (Stuart et al., 2019) and assess the results by comparison these silver standards.
7. Discussion
In this paper, introduce and generic family of graphical models, the marginally recoverable family, which is not limited to specific models or data types. We propose MMLE for parameter estimation with theoretical guaran- tees, avoiding high-dimensional integration, and introduce EM-MMLE to handle mixture distributions and capture het- erogeneous structures. The effectiveness of our method is demonstrated through several specific distributions.
To compare network predictions fairly, we use the AUPRC ratio, a metric from a previous benchmark (Pratapa et al., 2020) (see Appendix C.2). Algorithms are compared at a fixed network density (5%), with AUPRC ratios calcu- lated relative the silver standards. Figure presents the
To facilitate intuitive network representation, we impose certain restrictions on the form of the parameters within
8

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
the marginally recoverable family and require that the one- dimensional and two-dimensional marginal distributions satisfy specific conditions. However, motivated by the un- derlying principles of the family, we can provide more general definition that extends to distributions determined by higher moments.
for sequence count data. Nature Precedings, pp. 1-1, 2010.
Banerjee, O., El Ghaoui, L., and d' A'Apremoot, A. Model selection through sparse maximum likelihood estimation for multivariate gaussian binary data. The Journal of Machine Learning Research, 9:485-516, 2008
Let X be p-dimensional random variable, and let d p be fixed integer. We say that X is d-marginally recoverable if any d-dimensional marginal distribution of X belongs to the same distribution family, and the parameters of all d- dimensional marginal distributions collectively characterize the parameters the full distribution.
Cai, T. and Liu, W. Adaptive thresholding for sparse variance matrix estimation. Journal of the American Statistical Association, 10 106(494):672-6688 684, 2011.
Chiquet, J., Robin, S., and Mariadassou, M. Variational inference for sparse network reconstruction from count data. In International Conference on Machine Learning, pp. 1 1162-1171. PMLR, 2019.
This generalization improves the flexibility of the marginally recoverable family and broadens its applicability to more general parameter estimation problems. Future work will explore additional distributions within this framework and apply the proposed estimation method to diverse real-world applications.
Choi, Y., Coram, M., Peng, J., and Tang, H. A poisson log- normal model for constructing gene covariation network using RNA- data. Journal of Computational Biology, 24(7):721-731, 2017.
Acknowledgements
Cox, D. R. and Reid, N. A note on pscudolikelihood con- structed from marginal densities. Biometrika, 91(3):729- 737, 2004.
We thank the anonymous reviewers for their valuable com- ments and Zihao Chen for helpful discussions. This work was supported by the National Key R&D Program of China (2024YFF0507404 to R.X.), the National Natural Science Foundation of China (12425110, 10, 12371286 to R.X., 12271287 to W.W.), and the Sino-Russian Mathematics Center. Part of the analysis was performed on the high- performance computing platform of the Center for Life Sciences (Peking University).
d' Aspremont, A., Banerjee, O., and El Ghaoui, L. First- order methods for sparse covariance selection. SIAM Journal on Matrix Analysis and Applications, 30(1):56- 66, 2008.
Fan, J., Liu, H., Ning, Y., and Zou, H. High dimensional semiparametric latent graphical model for mixed data. Journal of the Royal Statistical Society Series B: Statisti- cal Methodology, 79(2):00-442) 2017.
Impact Statement
Friedman, J., Hastie, T., and Tibshirani, R. Sparse inverse
This paper presents work whose goal is to advance the field of machine learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
covariance estimation with the graphical lasso. Biostatis- tics, 9(3):432 441, 2008
Fu, A., Narasimhan, B., and Boyd. S. CVXR: An package for disciplined convex optimization. Journal of Statistical Software, 94(14):1 2020.
References
Hafemeister, C. and Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regular- ized negative binomial regression. Genome biology, 20 (1):296, 2019.
Abbasy, M. R., Nikfard, P., Ordi, A., and Torkaman, M. R. DNA base data hiding algorithm. International Journal of New Computer Architectures and their Applications (LJNCAA), 2(1):183-192 2012.
Han, H., Cho, J.-W., Lee, S., Yun, A., Kim, H., Bae, D., Yang, S., Kim, C. Y., Lee, M., Kim, E., et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic acids re- search, (OGD1))D330-DSSS0 2018.
Aitchison, J. and Ho, C. The multivariate Poisson-log nor- mal distribution. Biometrika, 7644)5643-653, 1989.
Amaratunga, D., Cabrera, J., and Shkedy, Z. Explo- ration and analysis of DNA microurray and other high- dimensional data. John Wiley & Sons, 2014.
Hu, H., Miao, Y.R... Jia, L.-H., Yu, Q.-Y, Zhang, Q., and Guo, A.Y AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription fac- tors. Nucleic acids research, 47(D1):D33-DDS) 2019.
Anders, S. and Huber, W. Differential expression analysis

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Huynh-Thu, V. A., Intthum, A., Wehenkel, L., and Geurts, Inferring regulatory networks from expression data using tree-based methods. PloS one, 5(9)e 2010.
Qiu, C.-G., Shen, B., and Sun, X.-Q. Significant biomarkers identification associated with cutaneous squamous cell carcinoma progression. International Journal of General Medicine, pp. 2347-2366 2022.
Islam, S., Zeisel, A., Joost, S., La Manno, G., Zajac, P., Kasper, M., Lönnerberg, P., and Linnarsson, S. Quantita- tive single-cell rna-seq with unique molecular identifiers. Nature methods, 11(2):663-166 2014.
Sarkar, and Stephens, M. Separating measurement and expression models clarifies confusion in single-cell RNA sequencing analysis. Nature genetics, 53(6))770-777, 2021.
Kim, S. ppcor: an package for a fast calculation to semi-partial correlation coefficients. Communications for statistical applications and methods, 22(()):655, 2015.
Shao, J. Mathematical statistics. Springer Science & Busi- ness Media, 2003.
Lachmann, A., Xu, H., Krishnan, J., Berger, S. I., Mazloom, A. R., and Ma' ayan, A. ChEA: transcription factor reg- ulation inferred from integrating genome-wide ChIP-X experiments. Bioinformatics, 26(19):2438-244 2010.
Silva, A., Rothstein, S. J., McNicholas, P. D., and Subedi, S. A multivariate Poisson-lov normal mixture model for clustering transcriptome sequencing data. BMC bioinfor- matics, 20:1-11, 2019.
Lauritzen, S. Graphical models. Clarendon Press, 1996.
Sinclair, D. and Hooker, G. Sparse inverse covariance esti- mation for high- thgohhhruughpu microRNA sequencing data in the Poisson log-normal graphical model. Journal of Sta- tistical Computation and Simulation, 89(16):3105-3117, 2019.
Li, X., Zhang, X., Wu, C.-C., P.P. Fu, Y--.,, Xie, L.-H., Sun, S.,S., Zhou, Y.-Y., and Zhu, B.-L. The role of MYB proto-oncogene like 2 in tamoxifen resistance in breast cancer. Journal of Molecular Histology, 52:21-30, 2021.
Skrondal, A. and Rabe-Hesketh, S. Latent variable mod- elling: A survey. Scandinavian Journal of Statistics, 34 (4):712-745, 2007.
Liu, H., Lafferty, J., and Wasserman, L. The nonparanor- mal: semiparametric estimation high dimensional undi- rected graphs. Journal of Machine Learning Research, 10(10), 2009.
Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M., Hao, Y., Stoeckius, M., Smibert, P., and Satija, R. Comprehensive integration of single-cell data. cell, 177())1888-1902, 2019.
Liu, H., Han, F. Yuan, M., Lafferty, J., and Wasserman, L. figh-dimensiooal semiparametric Gaussian copula graphical models. 2012.
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., Simonovic, M., Doncheva, N. T. Mor- ris, J. H., Bork, P., et al. STRING vl1: protein-protein association networks with increased coverage, support- ing functional discovery in genome-wide experimental datasets. Nucleic acids research, 47DDI):D077-D613, 2019.
Liu, Z.-P., Wu, C., Miao, H., and Wu, H. RegNet- work: an integrated database of transcriptional and post- transcriptional regulatory networks human and mouse. Database, 2015 b015:bav095, 2015.
Meinshausen, N. and Bühlmann, P. High-dimensional graphs and variable selection with the lasso. 2006.
Tang, J., Wang, C., Xiao, F., and Xi, R. Single-cell gene regulatory network analysis for mixed cell populations. Quantitative Biology, 2024.
Oki, S., Ohta, T., Shioi, G., Hatanaka, H., Ogasawara, O., Okuda, Y., Kawaji, H., Nakaki, R., Sese, J., and Meno, C. ChIP-Atlas: a data-mining suite powered by full integration of public ChIP-seq data. EMBO reports, 19 (12): 2018.
Trinh, G., Rungie, C., Wright, M., Driesener, C., and Dawes, J. Predicting future purchases with the Poisson log- normal model. Marketing Letters, 25:219-2344 2014.
Peng, J.. Wang. P., Zhou, N., and Zhu, J. Partial correlation estimation by joint sparse regression models. Journal of the American Statistical Association, 104(486):735-46, 2009.
Tunaru, R. Hicrarchical Bayesian models for multiple count data. Austrian Journal of statistics, 31(2&:3) 2002.
Varin, C. On composite marginal likelihoods. Asta advances in stutistical analysis, 92(1):1-22 2008
Pratapa, A., Jalihal, A. P., Law, J. N., Bharadwaj. A., and Murali, T. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature methods, ((22)1147-155 2020.
Varin, C., Reid, N., and Firth, D. An overview of composite likelihood methods. Statistica Sinica, pp. 5-42, 2011.
10

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Wong, W. H. and Shen, X. Probability inequalities for likelihood ratios and convergence rates of sieve MI.Es. The Annals of Statistics, pp. 339-362, 1995.
Wu, H., Deng, X., and Ramakrishnan, N. Sparse estimation of multivariate Poisson log-normal models from count data. Statistical Analysis and Data Mining: The ASA
Data Science Journal, 11(2):66-77, 2018.
Xiao, F., Tang, J., Fang, H., and Xi, Estimating graphical models for count data with applications single-cell gene network. Advances in Neural Information Processing Systems, 35:29038-29053 2022.
Xu, H., Baroukh, C., Dannenfelser, R., Chen, E. Y., Tan, C. M., Kou, Y., Kim, Y. E., Lemischka, R., and Ma'ay Ma'ayan, A. ESCAPE: database for integrating high- content published data collected from human and mouse embryonic stem cells. Database, 201 2013:bat045, 2013.
Xue. L. and Zou, H. Regularized rank-based estimation of high -igheninnni nonparanormal graphical models. 2012.
Yakowitz, and Spragins, D. On identifiability of finite mixtures. The Annals of Mathematical Statistics, -901)220--214 1968.
Yuan, M. and Lin, Y. Model selection and estimation in the Gaussian graphical model. Biometrika, 94(1):19-35, 2007.
Zhang, Q., Liu, W., Zhang, H.-M., Xie, G.-Y., Miao, Y.-R., Xia, M., and Guo, A.-Y. hTFiarget: comprehensive database for regulations human transcription factors and their targets. Genomics, Proteomics and Bioinformat- ics, 18(2):120-125 2020.
Zhang, T. and Zou, H. Sparse precision matrix estimation via lasso penalized D-trace loss. Biometrika, 101(1): 103-120, 2014.
Zheng, G. X., Terry, J. M., Belgrader, P., Ryvkin, P., Bent, Z. W., Wilson, R., Ziraldo, S. B., Wheeler, T. D., Mc- Dermott, G. P., Zhu, J., et al. Massively parallel digital transcriptional profiling of single cells. Nature communi- cations, 3(1):14049, 2017.
Zhou, K.-R., Liu. S., Sun, W.-J., Zheng, L.-L., Zhou, H., Yang, J..H., and Qu, L.-H. ChIPBase v2. O: decoding transcriptional regulatory networks of non-coding RNAs and protein-coding genes from ChIP-seq data. Nucleic acids research, pp. gkw965, 2016.

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
A. Theoretical Results and Proofs
A.1. Proofs for Theorem 4.3
In this section, we first provide some lemmas concerning the Maximum Likelihood Estimator (MLE) and subsequently prove Theorem 4.3.
Recall that h (1; h. (y; denote the one-dimensional and two-dimensional marginal density functions of H (T, 2), respectively, with parameters w1 E On and w2 E O2. where O1, O2 are parameter spaces in finite-iimensional Euclidean space. For simplicity, we omit the subscripts and use hM (y: w) to represent either the me-dimensionall or two-dimensional marginal density function of HD (TT, S2,, where w O is k-dimensional parameter vector. Define
F= {TrM w) : EO.
For any u > 0, if there exists finite set {V5 1f);; -, } such that |1, -11)) S u for = l,..., N and for any hM w) F. there exists j such that f/ < (y; s 1 we say that {(f - -1, N) is Hellinger u-bracketing of F, and N is the size of the Hellinger u-bracketing. Let Nr, be the set of sizes of all Hellinger u-bracketings. We define the bracketing Hellinger metric entropy of F as
H), min log N. NEN.
A.1.1. LEMMAS
Lemma A.1. There exist positive constants C1. C22 C3, such that, for any > if
/2x
12
u//22 {or F h (y:w)) S v2s}) du 12
(9)
32
for all s 2 e, then for the MLE is of the true parameter w using n independent samples, we have
pr (11001/2 (y: - 1 (y: (x: La 2c) S 5exp(-cincc
Lemma A. is a local version of Theorems and 2 from Wong Shen (1995).
Lemma A.2. Let is be the MLE of the true parameter w restricted on O. Under Condition 4. I and Condition 4.2, there exists a positive constant e independent with parameters, such that, for any > O, we have
pr (1 wll2 2) 5exp(-cnee).
Proof. To prove Lemma A.2, we adopt a proof framework similar to that of Lemma S11 in Xiao et al. (2022). First we will show that there exist postive constants c4, cs. such that,
H uu. pu E F d(pr, h (y: S <2s)) s ca log(ess/u).
(10)
According to the lower boundedness condition 4.1, the parameters whose index in {pI E F d(pi, hM (y:w)) S /25} can be covered by Fure {w"l || 2اl2 v2s/C1,, }. where Cj is constant. Using the upper boundedness condition 4.2, for any w, w' O, there exist measurable function m(y) and constant C2 such that m² (y)dv m C² co, and 112- It is straightforward to verify that can cover the set Fu.. by h 1/2 (y:w) - h 1/2 (y: w/) I m(y) ||ا w'll
most (2v2C2s/Cru) balls, while each ball has radius of u/2C2. For any ball B with radius u/ u/2C2, we define the centre of B Wa, then |w wollz S u/2C, for any w/ B. Then, we have 1/2 (y:w²) h 1/2 ( (y:wo)) S mly)u/2C2z 1/2
1/2
Therefore, we can select the minimum and maximum density within each ball as follows:
1 max { 11 1/2 (y:
(123 mly)u/2C2. 0). s1/2 1,1 (y:wo) + mly)u//222.
12

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity Consequently, have d( fu) s (J m² (y)u
u, thus can derive (10).
Next, to apply Lemma A.I, we aim to prove the existence of positive constants C2 and c' such that (9) holds for any
$ 2 c'n-1/2 The case that < d/2 will be discussed later. According the Cauchy inequality (Ix² f1/2 dx) " -
(b- a) la¹ f dr.
(1 11 (w/o2, (VV2w//o
d(pi, N S v20}) du (y:w)) (x:ww) v28) du
( H
S(V2s - 22121
fni EF: d(pi, (u/osffo
22/25
(11)
/
22
Sv2s
(u/o2, Pi EF: :(pi, (y; (y:w)) s v22)) du
2/120
SVDs
log((2css/ujdd
22/29
Thus, it suffices to show there exist positive constants C2 and d' such that
1 /20
121/2 128
log log(c2055/1) du S ns²,
(12)
for all 2 -1/2
Let C2 2 %5
After calculating the left hand of (12), we have
[ J33/28
log( log(c2055/u)dde log(e2(5 (0022ss)(((2)) 52/2²) + (u u log u| 1/2
<- 2--2ssogg 20logl 22) 5² /22 (log(() log(s²/28)) 122)) (V2s - 22/22) log(c22s) VV2 (13)
= Vlog(c205 22)s s² 1/2a log(2>e/s)
Note that v2s 2 82/28 implies log(2%e/s) > 0. Therefore, if nss 2 v2 1og(e203e// i.e. 2 1/22 where d' = (V2log(c20ze// 1/2 inequality (12) will hold.
Then using Lemma A. I, there exists a positive constant ei such that
pr (11:012 (y: (0) - h A 1/2 (y: (y:w)/ La
20) s 5exp(-cenee)
for any 2 1/2 Notice that for << c'n we can have constant cp satisfy
pr (11001// (y: (o) - 1 (y:w) 2 20) s (exp(=cune?)
Choosing constant min(co. )) we have for any <> 0,
pr (1
(y: - 1, (y; [x;]) I11 La
20) < 5exp(-ene?
Noting that w, w E O, we apply the lower boundedness condition 4.1 to obtain
{ {d(h M (y; ()) h (y; < } C {1اا wll2 < </O
where C1 is constant. Thus.
pr (1|@ wll2 < <</Ci) 2 pr (d(h (y; / (y: <e) 2 1 - 5exp(-cne's
which implics that
pr (1|ا 8ااس こ - 5exp(-Cne()
for a constant C' and all é' > 0. Then we finish our proof.
13

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
A.1.2. PROOF OF THEOREM 4.3
Proof. Recall the definition of MMLE in Equation (5), Apply Lemma A.2 to Y, and Y(. for any S jck S p. we have pr (lāg5k Jaskl 2 6) S 5xp (-Cne) for any 1 S j S k S p. S S G and 6 > O. Thus, we have pr - 2 20) & 5Gp² exp (-cne?) for any 1 S g S G and é > O with a constant e and finish the proof.
A.2. Detailed Restatement of Theorem 4.4
For a vector a with the i-th entry (s, let ||a||| = Ei lail. For a matrix A, let Xmax and Amin (A) be the largest and smallest eigenvalues of A, (ec(A) denote the vector formed by stacking the columns of A, || = Ei, [(A # 0) denote the number of nonzero entries, || I = Eiiy [Ayl and ||All,, = max, (E, Assl).
For each g 1, .... G, we introduce the following notation. Define Gg = {(i, g) * 0} as the set of positions corresponding to nonzero elements in Og and G {(t,j))@gi = 0}. Let d denote the maximum node degree in θ g. define - = r(2gg)
| g lo. and min (,,,,, 1r be the minimal absolute value of nonzero elements of Og. We
( @ I I+I& & 2g)/2, where A @ B is the Kronecker product. For two subsets Ti and T2 of {(i,)) I S i.j S p). we define I(A) .T be the submatrix of ((A) whose rows and columns indexed by Ti and T2. respectively. Other notations are as follows,
Fiin,a (T (Ta,,,,(T,...) (Ta,,c,)"'(
7 1 max (K))EG;
kr o,, 11
1,000 .kk = 1g
Recall that the precision matrix estimator G is defined as:
G arg min e,
((t=)))t(((,) + +1g1(
(14)
We now provide detailed restatement of Theorem 4.4, aligning the conditions with those in Theorem 3 of Xiao et al. (2022).
For each g= = 1, ..... G, there exist constants 0g and byr such that for some 1) > 2, if the true precision matrix A g satisfies the following conditions:
n >>>>> b1 (nlog + log ag) max [12dkr; 127 (ksk?? + kr), {12 (kyaki + kf) + 5ck?) 0-1
min { 81/2 d + 11 {127 (ksskit + kt) + 5dk}) X (0))
and
1/2 12 (nloog + 1/2 Xy 12y - (kg + kr) 1/2
then
pr (vec(
3 g).5 = =0) - 2-
where vcc(@y)og represents the subvector indexed by G (i. Agis + 0}.
The proof of Theorem 4.4 is similar that of Theorem 3 in Xiao al. (2022) and is omitted here.
A.3. Theoretical Results and Proofs for the MPLN Model
Recall the definition of the PLN distribution from the manuscript:
y x Poisson {S exp (Iz)),
(15)
j=1
X ~ N(H, E),
14

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
where S denotes the total sequencing reads, which be estimated as the sum of counts per cell or by other methods (Anders & Huber, 2010; Hafemeister & Satija, 2019). Without loss of generality, we assume that S is known and constant. For simplicity, we set 5 = 1 in the proof.
We hereafter represent the PLN distribution as PLN(O) (+), where e = Σ E-' For the MPLN distribution, defined as ES=1 T PLN g Hg), we assume that the true means My and proportions my 1, ...... G) are known.
We impose the following conditions:
Condition A.3. There exist positive constants m, M, and I, such that maxigjAso ga S and m Amin (o)) S
Amax (Og) I M for g = ==....,,G.
Condition A.4. The parameters vech(O@))))) (g 1, ,G) are different from each other.
Based on these conditions, we establish the following theoretical results for the MPLN model:
Theorem A.5. Under Conditions A.3 and A.4, the MPLN model is identifiable, and its Fisher information matrix is positive definite.
Theorem A.6. Under Conditions A.3 and A.4, the MPLN model satisfies Conditions 4.1 and 4.2.
Theorem A.5 establishes the identifiability of the MPLN model and the positive definiteness of its Fisher information matrix, ensuring that the model behaves well under relatively mild conditions. Building on Theorem A.5, we demonstrate Theorem A.6, which represents key theoretical contribution of this study.
In this section, we first introduce some notations, followed by the presentation of several lemmas. Finally, we provide proofs for Theorem A.5 and Theorem A.6.
A. NOTATION
We define two vectorization operators, vech and vech2. For symmetric matrix A (asklisiAs E RPXP, vech(A) defined as
vech(A) (A11, 012, ......., 01ps 022, 023.....
Ogp: agp:.... aip i\p 1), G(p 1)px app)
and (A) is
vecho A) (A11, 2012, 2013, 201ps 022, 2023.
202p. ............. a(p-1)(p-1)))) 20 10(p-11pp l)p. App) )pp)
Note that vech and vech2 only differ at i # elements.
We represent the MPLN distribution as MPLN (v, ( where v = (IT. 28) and v, = vech(
Define &p = [ bijective. the following discussion, defining either &p or v implicitly determines the other.
pT.... .T)
where 9 = vech Eg). Since Eg is positive definite, the mapping between &o and is
Additionally, we define the bounded sets:
D. = (v -VI _
......,)
max 19yjal S 1', m S Amin (0g) S Amax (0g) S M,g 1,...,. G}. 1SSASP
and
O (1p (yTT...,,pp) T
max la S I, S Amin (Eg) Amar (2g) S M 1,.... 1Sj.ksp
We assume the true parameter n² is interior point of D.M and the true parameter up* is restricted on O.M.
p
We denote the density of the PLN distribution PLN(O μ) by ply: O. 4)) or p(y: E, μ)), and the density of the MPLN distribution MPLN(U, {1r (11)) by p (y; (1n))= or p (y: 4, (2))) In the following sections, we always write h(y, x) = 11 exp(x,U)) exp( exx(xx))) Given a single sample i, we write the likeiiild function of the PLN at Y: as
((O) y:) log 1,
log det((O)) exp (- 1 (- (x- 4) TO(x wi) h(y,,))dx C(y),
15

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
where Cly) is term that depends only y. Also, define (,)) Σi=1 ((O,);) the log-likelihood in the PLN. For the MPLN, its log-likelihood function Y; is
(v) = log
Σ 9=1
T Typ y. H4)
/ det (0₂)² exp (-1 -- (x- Mgg) ,²)T0,( uz)) M(xx,))xx O(y). "g
log
The log-likelihood of the MPLN model is (x(v)= ((v,y.). If we define
L) = (de((0)) (-5(x- (-5(x-1g)
(gg T0)(( - -11)) hiy,, x)dx,
and (v,)) = Eg= (v v). then ((v, = log ( + C(y.). Note that the function Lg (v y) is proportional to the density ply; θ Hz).
For the PLN model, denote the derivative (the score function) and the second order derivative (the Hessian matrix) of its log-likelihood as
@ y) vvech ())
8 @
S(O.y)
H) = Ôvech + (aachh@)@vechiee @e)
(16)
For the MPLN model, we can similarly define its score function StM y), its Hessian matrix H (v, y), and its Fisher information matrix -D("").
D(v) Ev- Eu-(H (v. y)).
(17)
Finally, we denote No as the set all p-dimensional non-negative integer vector. For a vector a (a),.... ... ap) we denote ||a||2 = Ei i az as its L2-norm and ||a|| = max, layl as its norm. For matrix A, we denote || ||2 as its largest singular value of A. Given ( and μl, we define operator T that maps functions in functions in y,
T(f)= / exp 1-2(x
1 (x n1) T0(x u)) f(x)h(y, x)dx.
In particular,
T(1) = f (-1 n1) T) - 10) h(y, x)dx,
where 1 denotes the constant function 1(x) = 1.
Definition A.7 (Good vector). We call vector (E1, di, Ea, aG) E R2G good vector if there exists g such that (Eg', )gg) A (g)) for all of + g. We call the index y good index with respect to &.
A.3.2. LEMMAS
Lemma A.8. Let ~ PLN(O, μ1). For any n. E N, we define
[1
n = 0, n >
((n,v) y11 -1)... +1
Then, for n = -m1, ",p)" ,Rp)T
we have
E
onnj 9/)
exp (n T # + θ 'n/2). mnn2))
16

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Proof. By the property of conditional expectation, have
E
11 66n,, us) = Ex
II
o(n, Uj)
1=1
From the moments of the Poisson distribution, we have
I w))l = II exp(ngay))
By
Further, since x ~ N(u,, ! we have
Ex
II exp(n)2z)
Ex(exp(n
x)) (x)) exp (n mo+nee o 'n/2).
and the conclusion follows.
Lemma A.9. Let & = Ea, aa) be a good vector with good index 8, satisfying 0g > 0fror g = 1,, G and
a = (a).-.. ag) T. for any 2
c
Σ 9=1
0g exp(E,= + + = =0,
then 0.
Proof. Without loss of generality, we assume that (Eg. og) (g 1,... G) are first sorted by og in increasing order, and for the same gg. they are further sorted by &g in increasing order. We say that (Eg. og) and (E.,a,) are equivalent (g) (Es, (.)) ). By this equivalence relationship, {(Eg, 0 can be partitioned into Q groups (Q Z 1). Let S, be the index set of the q-th group. We have
a exp(E)= + =
4=1 JES,
for all 2 E N. Dividing exp(5a= + oG22) on both sides of the above equation, we get
ΣΣu aj exp(5)> o,2² - (c2=-062)) + @=1 JES,
Σ jESq
=0
(18)
for all z N. have Ejesa aj O. By mathematical induction, we have Ejes, a; = for q I,..., Q. Since & is good vector with good index (E) itself forms a group, and hence a, = 0.
By the choice of &G, the first summation Equation (18) converges to zero when a goes infinity. So,
Lemma A.10. For any O, let M; C. RP, 1, .... nn be linear proper subspaces of R.. Then, there exists non-negative integer vector E Nº such that y & U M,
Proof. We prove by mathematical induction. The conclusion clearly holds for n = 1. Now we assume that Lemma A.10 holds for n and we aim to prove that it also holds for n + 1. Note that by linear algebra, U M; is proper subspace of IRP for all n.
By induction hypothesis, we can take a No AU" M, If Ma+1, we have a &f Upti M, and the proof finished. Thus, we only need to consider a E MM+1-- Similarly, we can take B E NP NU=-1 M; and ( E MM. For any + 1, we can prove that there at most one k such that a + Ki M₁. In fact, there are k k2 such that b & h2 and a + 3 E M;, a + 3 E M;, then B M,, which contradictory to the fact that 3 E No Urtt M, Furthermore. there is no i
such that a + kB M1. Otherwise, there exists k N such that a + KB M1. Then, we have E Mi, which is also contradiction. So could find most positive integers for k such that a + h E Upt! M,.
Since there are infinitely many non negaaivee numbers, we prove that there exists N such that a + kB o U Mi, and Lemma 10 is proved.
17

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Lemma A.11. Assume that A; E Rexp (i I,...,n) are non-zero matrix. For any > 0, let N₁ {x:x] A 0, E Ro) (i = -1ff=1.... , be proper subspaces of RR, Then, there exists a non-negative integer vector ? such that 7 $ N
Proof. We prove by mathematical induction. The conclusion clearly holds for 1. Now assume that Lemma A.11 holds for n and we aim to prove that it also holds for n + 1.
By induction hypothesis, we can take ,3 E NP U Ng. 3 & Mi, we have 3 E U Ni and the proof is finished. Thus, we only need to consider 3 Ni. Next, we can take a such that a & Ni. if + 1, then there are at most two integers E N satisfying the quadratic equation (a k3)) A,,a + kB3 = 0 + 1) because of 3 A,, o O.
1, then there
most integers k E N satisfying the quadratic equation (a + kB) A (a kB3) 0 because of ( A, 0 and x T A a A 0. follows that there are at most 2n + positive integers E N such that a + kB E Up+1 Nr. Since there are infinitely many non-negative numbers, we prove that there exists k E N such that a + kB # U N1, and Lemma A.11 proved.
Lemma A.12. Assume that A; E RPXP (i = n. are n. non-zero matrices. For any n > 0, let N {x: x" A, = 0, E } (i=1...., be n proper subspaces of RR. For any m > 0, let M; C CR, , =1, m be m linear proper subspaces. Then, there exists a non-negative integer vector y such that 7 & (U M.) U (G)
Proof. By the proof of Lemma A.10, there exists a & U. M,, We can assume a E Uf= M₁ Otherwise, we complete the proof. By the proof of Lemma A. 1,, there exists 7 & NG. On the one hand, since y E N" \ U NY, there are most integers & N satisfying the quadratic equation (ax kry) A,(a kry 0. On the other hand, for any M, = I,..., there is at most one integer such that KB E M, Otherwise, if there exist k o l2 and i satisfying
a KIY E M k2 E M,,
then we have (A:2 - k2 E M₁. It follows that a M₁, which is contradictory to the fact that a & U Mr Hence, there are at most 2m + n positive integers & E N such that a + ky E (U) M.) (U=_, N,). Since there are infinitely many non-negative integers, there exists an integer is such that ky & (U' M.) U N and Lemma A.12 is proved.
Lemma A.13. Let y and be random variables defined in for the PIN model, and let ((n,y) y) be the same in Lemma A.8. Let (M1, np) and T Rexp be P x p matrix. We have
E
p(n My tr (T (x =)((xx - 41) 11)
n +
TO nxp exp (n
u θ n/2).
Proof. By Lemma A.8
I (n,, n,) tr (TO-1) tr (TO-T)ess exp (n + o
E
Similar to the proof of Lemma A.8, by the moment generating function of the normal distribution, we have
E II p(ny, (T(x 1)(x 10)T)
Ex (exp(n T x) tr (T(x 4i(( - 1n
ttr (TO-1) + n θ -o ') exp (n T
r + n T o 'n/))).
Lemma A.13 follows from the above two equations.
18

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Lemma A.14. For any (Rg vech((,g))) (g p (y: 11),
G, that are bounded and different from each other, G, PG) are linearly independent.
Proof. We prove this by mathematical induction. The independence for G = 1 is trivial. Now we assume that Lemma A. 14 holds for 1. For any (#g v (g = 1,... G) that are bounded and different from each other, if we can prove that there exist = (01, ... aG) and an index E (1...., G. such that
Σ agp (y: gg Hg) 0 and as 0,
(19)
then by induction, we have a = 0 and hence ply; , (g) (g = .,..., G) are linearly independent. So our goal is to prove that if 2 (yp(y; ; g (2)) 0, then we can always find an index such that 0s O.
Let n = (n1,.... mp) T be any non-negative integer vector. Then, for any positive integer 2, by Lemma A.8, there exists a polynomial function Pa (y) II d( Mj) such that, if (19) holds, then
Σ 0g Ey (y)) = 0 and Ey (p=(y)) = exp (an (anJ,o, + +22 o
where Eg represents taking expectation with respect PLN Hg)) Let
= (n p1. T n/2......... Mc n o n/2).
By Lemma A.9, if there exists an n such that & is a good vector with good index s, then a, = 0 and we complete the proof. If, on the other hand, c (n" f.. n O 1/22..... nn Hg. n T " n/2) is not good vector for any non-negative integer vector n. Therefore, for any n, there exists # such that n T 11 = n μ4, and n T = n Toe1_n 'n. Thus, n the solution to the linear equation x T (11 - p.) = 0 and the quadratic equation T (0, -
) -6--) x 0. We define My as the linear space consisting of solutions to the linear equation x (11 14g) = 0 (g * 1). We define N₂ as the
space consisting of solutions to the quadratic equation x (0r θ-1) x 0 (g + 1). Define Q = Uy (Mg
Thus, for any non-negative integer vector n, we have E Q. Since (H4 vech (g) T) i.... G, are different from each other, Mg nNg is a proper subspace of TRP. More precisely, if u4 M4 then we have θ, of O., which means Mg NNg = Ng. Let I = {g- : 11 = =My}- Then, we have
n ne = (M, c UU g=
SEI
9E2
This contradicts Lemma A.12. Therefore, there exists vector such that is good vector. Consequently, for any (nT voch(O,) (g 1, (.....G) G) that are different from each other, ply: Q₁, fi)........ ply: P mc) are linearly inde- pendent.
Lemma A.15. For any 10 OM, there exist two functions (. and (.) ofy such that fM (y) S ply: sp (4g1g S (y), and] f JM (y)ddy < oo, PM (y) K (y)ddy < co for any polynomial K(y) of
Proof. By Remark of Xiao et al. (2022).
A.3.3. PROOF OF THEOREM A.5
Proof. By Yakowitz. & Spragins (1968), under Conditions A.3 and A.4, the identifiability of the MPLN model is equivalent to the linear independence of the PLN components. By Lemma A. 14, the first conclusion holds.
To establish the second conclusion of Theorem A.5, it is necessary to derive the explicit formulas for the score functions and Fisher information matrices of the PLN and MPLN models. clear that the densities of the PLN and MPLN satisfy the
19

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
regularity conditions in the literature (Shao, 2003). Then, can calculate the score function and the Fisher information as follows.
The Hessian matrix 7t (@.y y) of the PLN model is a p(p + 11/2 x p(p + 1)/2 matrix. For notational convenience, let (i. j)
-p be the index of the element at position (2p - + 1 j of the score function S and let H(LD) = ,00,, /2 p + j' of the Hessian matrix.
(i S j, S j') be the element at row (2p + 1 - p and column (2p- - 1'
Recall the Equation (16) above. The score function of the PLN model can be written
exp (->(x- 11) @ (4) vecha ((x - 11)(x 11) M(y, x)dx
y(O,y) 2
(0- vechg vechg (o
2
J exp (- T 0 - 11) hiy, x)dx
Especially, at the true parameter @* have S(O'. Ex (vech2 (0- (x (x 11 ly). By applying the operator T, the score function element at position (2p - i + 1)i/2 - p can be rewritten
T
vech2 ((x - (0)=1))(-- 14) -1)) ((5))
((,) (+,y)
vecha (0- (i,))
2
(())
Using the operator T, the Hessian matrix can be written follows. Let == -
When = j, j'.
IT (xx==1))
1 T (13) T (xx-1077
1 Enw Evi Ew
Ho (0.
4 T (1)
4
T2 (1)
When i j'
T (xx - 1x1(x- - (x),(xx-1)) - 11)))
T ((x - 11)²) T - -x) (x),)
H( (0,y) = Eig Dgn + 2 (())
2
T (1)
When "hen j
T (k-M),(e - -M)), - w)²
Hanees (0.y) = -EN + 2
T T(1))
1 T ((x - -1))3) T (x=)))x(x- (11)):)))))))))))))))))))))))))
2
(2(1)
When i 1 #
T T((x 1xx(x - ());(x--)((x-1))
1)((x
H(.))) (0,y) = -(Ew Eng + Eeg +
T
T((x 1xx(x -m1)) T ((x (C--))(( - ()))()-1)))))))))))))))))))))))
(2(1)
Now we consider the score function and the Fisher information matrix of the MPIN model. The score function of the MPLN model can be written as
S =( ( TT
(vI.)) dvl
oLa(va.3)) =
... .. (v.y))
TG
ÔvG
Similarly, we use the triplet (g. i,)), with < j. to index S Recall that V, = voch(@,)) Let (g. (.j),, (g', 1,, i'), j. 1' < be the index of the element at the (g- - 1)p(p + 1)/2 (2p - i + 1)i/2 p + row and (g' - 1)p(p + 1)/2+
(2p 1)²//2 /2- p + j' column of HM (v,y). respectively.
20

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
When g =
a
H (ox
SM (e
d0pr''
Mg (g.g((993) g.))
(H
(u,)
(gg (9.
( (0.))(
(20)
I LMM(, (u,y)
When g + &g' we have
Hig (...
(v.))
80g'y (.i
0
SM
T Lig O g. y) Lo yy
(21)
y
Recall the definition (17) of D(v). Then, -D(v") = E (S (v" Sy) (v", ()) is the Fisher information matrix of the MPLN at v*.
To establish positive definiteness, we show that there exists no non non-zero vector = (t1, .. tc) such
that B(tTSM((> ,y)* Tt) 0. Assuming that there exists vector (t,, ; tt) satisfying
E(t (CTSSM 8M((", (u""r)) (v' yy",))" t) O, we aim prove that t O. Since
T (V" (u', y) x)
Ep(y; (1) (t (M' y))² 0,
and under the assumption that (y: (H4 10 1 0 for all y. it follows that t SM SM(v", yy 0 for any y.
Then, when v= = v², we have
@Lg(Vg.y) 0.
L yv
T
dvg
Since L (v,y) + O, we have Eg t ng
Bba(,ay) 02g
0. Let (n ;pp)
be any non-negative integer vector and
(yy) - p(zn,, ,/), E N. Then
log Lg Lg(Ug. y) Ly (V = 0.
kg
Ov,
(22)
Since Lg y) proportional to the density the PLN with parameters O, and Mg (22) can be rewritten as
Olog (v
Σ
ply: Gy O(y: My)
g kg
Ov,
9=1
Summing over we get
G
i log (Vg, ply: O 0g (y)) = 0.
t T
dv,
By Fubini 's Theorem, get
6
i log (y
9 n9
ply; θ 1p. = 0.
(23)
2v n
1=1
y
21

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Then, let T be the symmetric matrix such that vech(T, g) :... For fixed g. have
T I log (vn p(y: ② Mzg (y)) = Fy, (
Olog Lu No)
dvg
dvg
2 N (M(yg) tr (T TT 1)(xg 11) 1))
where y, follows the PLN distribution with parameters O and Mg and Xg N (Mg: 2ggg is the corresponding latent variable. By Lemma A.13, we get
Ey, (t d log (o).
2 (n T exp exp(an) Mg + +22n "='n"o''n 902n/2)
Ov,
Then, Equation (23) can be rewritten as, for all = E E.
Σ T nn θ T exp (an My + 2 " O.
In order to show T1 = O, similar to the proof of the first conclusion, We define Mg as the linear space consist- ing of solutions to the linear equation x² (m4 (g) = 0 (g * 1). We define Ng as the space consisting of solu- tions to the quadratic equation x (0, -,,)) x = 0 (g * 1). Define Q = - (Mg Ng).) For any n & Q (n 1 n OO, n/,...., n PG Og n/2) good vector with a good index 1. Since 0, we must have
o T
(24)
By Lemma A.12. if 0 O is not a zero matrix, then there exists an such that n & Q and n - O T n * 0. which is contradictory to (24). Hence, we must have T 0 = and thus T₁ = O. Similarly, we get T, = 0 for all g = I,...,G. G. follows that t = 0, and we compete the proof.
A.3.4. PROOF OF THEOREM A.6
Proof. For the MPLN model, the marginal distribution is also an MPLN distribution. Therefore, is unnecessary to separately verify the one-dimensional and two-dimensional marginal density funetions. Instead, it suffices to verify that the p-dimensional density function of the MPLN distribution satisfies the inequalities stated in Conditions 4.1 and 4.2. Specifically, there exist measurable function myy) and constants C₁, C2, such that f my² du Ci co, for any
1 P2 O.M, let P1 p (y: 1=))))) and p2 (y; 22, (1g ] 1==1), we have
(i) C2 |||p| -2212 S d()),PP).
1/2
1/2 P2 S m(y) || 22ll2
(ii
To establish this result, we follow a proof framework similar to that of Lemma S9 in Xiao et al. (2022). Here we denote the score function and Hessian matrix of the log- likeliooloo in the PLN model with respect to as
y) (
S(E,y)
H(E,y) = (vech(E))veer Σ)
For the MPLN model, we can similarly denote the score function and Hessian matrix with respect to yp as S*M (4, and H. y).
1/2 First, we verify that (i) is satisfied. We get the Taylor expansion of p2 on 411-
(pp (y: (
(14
1/2
1/2
-1/2
-2)22 (422 41) R((por)) 12) (25)
P2
= P/ +P1
0p
kp
= =p + Pi
1/2 M
(T (422- - /2) + (522 - -T) R(pp,.)) (422 -11)
22

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
where R((Po: y) (y: Po: (1g)g=1) [SM
(4o: () /4 HM (4ou v)/2] and Po is between 91 and 422
We define p. = p Py: ( 1c ) and ò = l22 - 411 then we have
1/2
/ (Dit2 -P2 /2/2) du
d(p1.P2)
71/2 81)) du
=ppp22 (M)))) 8 / (p! ()
/2
SooT 8M (491, (41, (1),)) 8dv+
(26)
/ 221 1/2 1/2 SM
8M (1)1)) 580 R((oo,y)ddve v)ddv+
1/2
f |pOO" R(( R(yoo:yddd
/2
(I+
where I, II, and III are defined in an obvious way.
Define the minimum eigenvalue of Fup, (49 y) (()))))))))))) is Amin (49)): Amin(()))
then
1 8 E,, (S y) v) 1) 8 2 Amin (nin(())11011
(27)
Since 101 E Old defined on a compact set, and Eye, (SM (411, y) is positive definite and continuous for L there exists positive constant Clow > O, such that Crow S Amin (4p1) for any 40 E O,,, then we have
I 2 Clow 1012.
(28)
For part II, we have
m S f pı (T),)))) T8 210 du du + / pa (P 3)
(29)
-E,, (()))) _8 R( + po (M 88 ((.
For notational convenience, let (1,j) be the index of the element at position (2p- -pp-+ i 1)1/2 - p p in S, while (g. i, is similar index for S. with i I j. Similarly, for S j, S J' (1,j) refers to the element at row 1)i/2- - p p+3 j and column (2p i' + 1)d' p j' H and (g. (g,, (",'' indexes the element row (g - 1)p(p 1)/2 (2p- i 1)i/2 p j and column (g' 1)p(p 1)/2 (2p - -' 1):' 12j p
(2p
1pp-i4 i
d'+1))''
p j' in H.
We now claim that for y ~ MPLN(U, (1g )-=1) there exist two polynomial functions Ki(y) and K2(y) sat- isfying [yK < and E[[2(y)] < xo, such that for any g, i.ij i V. J. y)l Kí(y) and H (o..s) (f) (p < K2(y).
To prove this, first consider the case g = g'. Since Ty TgLg(U,y)/LL (v.y) $ 1, then we have,
(M v) < (( (Zg. v)
and
H.M. (g (9 (p. < (4 (E +
Then, by applying Remark 3 in Xiao al. (2022), know that, there exist two polynomial functions K2(y), K2(y) such that (S ( < Kı(y) and HH (o (p. < K2(y) with EyK) < oo and E(K2(y)) co. The same proof can be applied to the g # g' case. This completes the proof of the claim.
23

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity Therefore, there exist two polynomial functions K1(y) and Kg(y), such that
SSM (4P1.y) - || (4.. |8]|2 S Ki(y) ||0||
(30)
and
8" R(40,, < || 1)1 |001 - [R( y) |101 S (2 1181
(31)
Using the fact that any MPLN distribution have any finite moments, we can derive that the first part of Equation (29) can be bounded by C' || with a constant C. Then using Lemma A.15, we know the second part of Equation (29) satisfies
f P.
f S (P), 88 1ipoo)00 du { ((
du S C' [1812
with constant C Then we can derive
II| S (C+ 1811
(32)
For III = f p. R((Po ((P0o))0) du, using (31) and similar technique in II, we have
m S C" ||81|
(33)
with a constant C".
Since the dominating functions fup K1, K2 are all independent from parameters, then the constants C, C', C" are all independent from parameters. Thus, there exist positive constants K and C22 both independent of the parameters, such that for any ||9; P2112 I K, we have C₁2 ||/p] 42112 S (P),
Next, we prove that (i) holds for any P1: P2 E O.d. Given that maxi |0 S 1, there exists positive constant K" such that ||ب; 2221l l1 K' for any ( P2 E OM Define the set W { 422 |||; P2lla K, 1 P2 E O.d By the identifiability of the MPLN model, we have dl (P1.P2) o 0 when P1 * P2. Since d (PI, P2) is continuous function
p J. and W is a closed and bounded domain, d (P1, P2) has a minimum value 4 > 0 on W. For any (41:492) E W, we have k ||+p 921l2 / K" S d ((1, Let C2 = min (C2, ki/K") K'). Then, we have (i) holds for any P1 92 E O.M.
Finally, we verify (ii) Letpi P1 ES 1= Pg and p2 = Eg= TpP2g. where Pi, ply: My and p2: = ply; Eg ng) According to the Lemma S10 in Xiao et al. (2022), there exist measurable functions my and constants Cg. such that mi(y)dv )dv Cg < and Pig have
1/2
1/2
< myly) ||4g 4 l for 1 S g I G. Using Cauchy- Seuchy-Schwarz inequality, we
P2g
G TaPigy
G
1/2 P1
1/2
P2
kg (pig + p2n) - 2
P2
g=1
G
Fg (Pig + p22)) 2
>
TgPig
1/2 1/2
P2g
(34)
G
2
Σ
2/2 - p2
Pig
9-1
< Σ my 2 (y) ||,,, -
g
9=1
Let my = mg (y), then f (y?(y)dv= = C2 < oo for a constant C1 and
1/2
1/2
< m my
Σ 1 4.112 = 353(-)3( m²(yy) |/1 40112-
p2
So P1 P2
1/2 1/2
< m(y) 1/:p -p'1l2-
24

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
A.4. Theoretical Results and Proofs for the Binary Data Model in Example 2.6
Let (y1 ........30 Bp)) 3p) represent p-dimensional binary variables and (Iiy..., Ep)) be p-dimensional random vector. Assume that y follows the latent Gaussian copula model for binary data defined in Example 2.6, hereafter referred to the binary model, then have
Vj = [2, > Cÿ) x ~ NP,(O, f)
(35)
where / (+) the indicator function and Σjj for any S
Let C _ The joint density function of y E {0, 1) given by:
1 dx
ply; C) (2T) det ( 1/ xEUly)
exp
(-5x
where Uyy) U1(11) x ... Up(Up), with:
(g(4,) [ff(C)), (-o,
if = 1
if リ = 0
for 1 S p.
Given E and C, we define an operator IC that maps functions in to functions in y.
K(g) [ exp
(-2x Σ xx xx g(x)dx.
(XEU()
In particular,
K((1) / | exp
1 1-xxI=11 dx.
where 1 denotes the constant function 1(x) =
Let + --- denote the precision matrix. Assuming that is known, we establish the following theoretical results for the binary model:
Condition A.16. There exist positive constants m and M, such that m S Amin ()) S Amax (e) S M.
Theorem A.17. Under Condition A.16, the binary model satisfies Conditions 4.1 and 4.2.
Proof. For the binary model, marginal distributions of any dimension remain consistent with the binary model framework. One-dimensional marginal density functions do not depend any parameters. Thus, focus the two-dimensional case, where y = (y1, 3/2
y I (xy > Cg) 1,2 NPN(0, Σ,f).
(36)
where x = (E1, F2) and
= o 6=
(37)
For notational simplicity, we denote the density function as hu 0r) According to the Condition A.16, we have -1 -1<
< C < for some constants and C. Define the bounded set D = -1 < <<s <
We only need to prove that there exist a measurable function myyy) and constants C1, Cg. such that my) du = C9 < <0, and for any σ1, 02 E D, let P1 les (y: 01) and P2 = h2 (y; 02), we have:
ii) C2lo1 - 02l s d(pr,p2).
25

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
1/2 1/2 (ii) Ip P2
S m - -2/-
First, we demonstrate that (i) satisfied, following the proof of Theorem A.6.
We denote the score function and Hessian matrix of the log-likelihood. Straightforward calculation shows that
S(o,y) = log oy)
do
(1+ (((()))) ((11))))))))) - x2))
02-02)3((1)
1-
H(aay) 8 log O)) do2 (1+02) + OKK(((] +2))) -20 + + +2²))
(38)
(1 (1-02)(((1)
(11+0") - ak{ + 2²))²
1+
+ (1-02)2
(1 -1-02 (1-02)+(3(1)
(6o + 200 - (1 31 K(r +
(1 (1-02)(((1)
1/2 We get the Taylor expansion of P2
01
1/2
=p
1/2 + pi / S(1 (02 01) /2 (02
(39)
P2
where R(on.y) = (y: 0o)
5) H(oo.y)/2) y)/2] and do is between 01 and 02. We define po h2 (y; oo) and
ò 02 01, then we have
1/2
/ DiS(on,y)88t /
1/2 1/2 v)ooo) de + I /ooloo 2pp Po pp
d()).P2)
/2
(40)
11+ 111 1/2/2
where I, II, and III are defined in an obvious way. Since 01 E D defined on a compact set, and Foz, ()) is positive definite and continuous for σı, there exists positive constant Clow 0, such that for any σi D,
12 Clow 2
(41)
For part II, we have
I < f pi [8oo,.))R(oo, du / (((o,,)) y)83| du.
(42)
-E., () / pol y)
Since takes finite number of values, according to Equation (38), there exist constants K₃ and K2, independent of the parameters, such that (S(a, < K1 and |H(a, y) < K2 for any E D. Similarly, we have fyg (y; o) S K3, where K3 is constant.
Then we can derive
|I| S 8
(43)
with constant C².
For part III pollooo.y)&e 8 dv, using similar technique in 11, we have
|T S 1"1
(44)
26

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
with constant C".
Since the constants K1, K2, C², ,, are all independent from parameters, there exist positive constants K and C22 both independent of the parameters, such that C2 02| S d (pi. P2) for any |0i -02l S K.
Next, we prove that (i) holds for any 01, 02 E 'D. Given that 1 < c - o - C < 1, there exists a positive constant K" such that | - 02l S K' for any 01, 01,02 E D. Define W {(0) 02) 01 02l 2 K, 11, 02 E D). By the identifiability of the binary model, pp) + 0 when pl + p2- Since p2) is continuous and W is compact, attains a minimum value L on W. For any (2) W, k 01 - 02l/K' - d(pi, p2). Let C22 min / k/KK). Thus, (i) holds for all 01, 02 D.
Finally, we verify (ii). For any 01, D, let pl ha (y: and p2 h2 (2) Similarly, we have Taylor expansion of 1/2 P2 on 01.
P2 = Pi + p.
1/2
1122
1pp
(02 -0)
000
where o² is between 01 and d2, and p. h22 (y;
Let 1) = 02 01, we have,
1/2
1/2
1 1/2 y)nn)))nl= S 2
1 11/21 Ki
P2 Pi
2"
So we finish the proof.
B. Simulation
B.1. Details of the Data Generation Process for the MPLN Model
The network structures are generated based the following procedures.
Random Graph: Pairs of nodes are connected with probability 0.016 and the nonzero edges are randomly set as 0.3 or -0.3.
Blocked Random Graph: The blocked random graph comprises 5 independent groups of randomly connected nodes. Pairs of nodes within the same group are connected with probability 0.016, while nodes belonging to separate groups are not related. The nonzero edges are randomly set as 0.3 or -0.3.
Banded Graph: Pairs (i,j) of nodes are connected if | - j| < 2, = j. All nonzero edges are set as 0.3.
Scale-free Graph: A feele-free network follows a power law which suggests that the central node have more connections. The scale-free graphs are generated with power and the nonzero edges are randomly set as 0.3 or -0.3.
For each simulation dataset, we independently generate the precision matrix for each of the three latent normal distributions based on one of four graph structures. To ensure positive definiteness, the diagonal elements of the precision matrix are set to 1 plus small positive value. The mean vectors Hg (Agl...... Mppp) for each latent normal distribution (g 1, 2, 3) are then generated as follows: the first pa elements of f4, are independently sampled from {v1, (vi + 02) w2), and the remaining p- Pa elements are shared across μ1, H2: f3g, independently sampled from (Vy, ). The values of (V1, U2, v3 va) are set to (1.9, -0 0.6, 0.4, -0.6) for the low zero-proportion case (about 20% zeros) and (1.1, - 1.4, - -0.4, - 1.44 for the high zero-proportion case (about 40% zeros). We vary Pa to control the mixing degree of the three populations.
The scaling factors S- (S1. ... . ST) are independently sampled from distribution defined as S. = exp (2,), where Zi N(log(log 10), 0.05²) for i 1, 1,...,"! n. Using these model parameters, the observed expression values YI,..., Y,, generated from the MPLN model. We then compute the Adjusted Rand Index (ARI) between the true population labels and the labels obtained from K-means clustering of the normalized data
Y: + 1 )
Y, = log(
S,
27

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Table Comparisons of EM-MMLE with Pl.Net, VMPLN and Glasso in terms of AUPR on simulation results for blocked random graphs generated by the MPLN model. The results are averages over 50 replicates with standard deviations in brackets.
Zero-proportion
Low
Dimension
p 100 Middle
300 Middle
Mixing degree
Low
High
Low
High
1800
EM-MMLE
0.95 (0.011)
0.93 (0.017)
0.9 (0.016)
0.82 (0.028)
0.78 (0.023)
0.67 (0.043)
PLNet
0.86 (0.05)
0.79 (0.084)
0.79 (0.074)
0.7 (0.033)
0.65 (0.056)
0.52 (0.054)
VMPLN
0.91 (0.028)
0.88 (0.037)
0.89 (0.024)
0.62 (0.024)
0.62 (0.019)
0.6 (0. 025)
Glasso
0.82 (0.038)
0.78 (0.038)
0.75 (0.05)
0.66 (0.032)
0.64 (0.027)
0.56 (0.037)
3000
EM-MMLE
0.99 (0.005)
0.98 (0.005)
0.98 (0.006)
0.94 (0.014)
0.91 (0.013)
0.84 (0.047)
PI.Net
0.96 (0.02)
0.95 (0.041)
0.91 (0.057)
0.87 (0.062)
0.83 (0.048)
0.73 (0.068)
VMPLN
0.95 (0.03)
0.93 (0.036)
0.94 (0.021)
0.7 (0.021)
0.7 (0.02)
0.7 (0.024)
Glasso
0.89 (0.037)
0.87 (0.043)
0.85 (0.041)
0.75 (0.03)
0.74 (0.021)
0.68 (0.039)
Zero-proportion
High
Dimension
p p= 100 Middle
p 300 Middle
Mixing degree
Low
High
Low
High
1800
EM-MMLE
0.8 (0.029)
0.78 (0.028)
0.74 (0.035)
0.58 (0.048)
0.53 (0.048)
0.45 (0.045)
PLNet
0.7 (0.064)
0.67 (0.049)
0.63 (0.045)
0.49 (0.054)
0.44 (0.049)
0.35 (0.055)
VMPLN
0.76 (0.054)
0.74 (0.034)
0.71 (0.036)
0,43 (0.034)
0.42 (0.031)
0.41 (0.036)
Glasso
0.52 (0.045)
0.47 (0.036)
0.45 (0.048)
0.44 (0.037)
0.42 (0.031)
0.36 (0.038)
3000
EM-MMLE
0.95 (0.012)
0.93 (0.014)
0.91 (0.018)
0.79 (0.037)
0.75 (0.037)
0.68 (0.039)
PLNet
0.87 (0.061)
0.86 (0.034)
0.8 (0. (0,062)
0.73 (0.057)
0.67 (0.057)
0.6 (0.047)
VMPLN
0.86 (0.048)
0.84 (0.043)
0.83 (0.046)
0.53 (0.028)
0.53 (0.023)
0.51 (0.019)
Glasso
0.62 (0.058)
0.59 (0.041)
0.54 (0.067)
0.55 (0.034)
0.54 (0.028)
0.49 (0.029)
where 3, Σラ Y,,/10 f 1.2,...,mo n. The ARI values are adjusted as follows: low-level mixing data correspond an ARI in (0.9, 1], middle-level mixing data to an ARI in (0.75, 0.85], and high-level mixing data to ARI in (0.6, 0.7].
B.2. Details of the Data Generation Process for Binary Data
The threshold parameter Cj is sampled from uniform distribution on [-1, 1]. We construct the inverse correlation matrix Q such that ,, =1 1 and θ f 0zzzk for j * k. Here, we set 01 0.15 ensure the positive definiteness of O and zys: Bernoulli random variable with success probability
400
|t, tz 112
Pik:
p(p 1) exp
20²
where ty and t& are independently drawn from a bivariate uniform distribution on j0, 1]. The constant 0z is adjusted to generate approximately 200 edges in the graph, and the covariance matrix Σ is rescaled so that all diagonal elements equal
B.3. Additional Results
Tables 3-5 present the AUPR performance of EM-MMLE, PLNet, VMPLN, and Glasso for blocked random graphs, banded graphs, and scale-free graphs, respectively.
C. Real Data Analysis
C.1. Silver Standard Construction for Benchmarking on scRNA-seq Data
The databases utilized in this study include STRING (Szklarezyk et al., 2019), HumanTHDB (Hu et al., 2019), hTFtaaget (Zhang et al., 2020), ChEA (Lachmann et al., 2010), ChIP-Atlas (Oki et al., 2018), ChIPBase (Zhou et al., 2016), ESCAPE (Xu et al., 2013), TRRUST (Han al., 2018), and RegNetwork (Liu al., 2015).
28

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Table Comparisons of EM-MMLE with PLNet, VMPLN and Glasso in terms of AUPR on simulation results for banded graphs generated by the MPLN model. The results are averages over 50 replicates with standard deviations in brackets.
Zero-proportion Dimension
Low
p 100 Middle
p 300 Middle
Mixing degree
Low
High
Low
High
1800
EM-MMLE
0.87 (0.017)
0.88 (0.017)
0.85 (0.022)
0.87 (0.01)
0.83 (0.017)
0.76 (0.021)
PLNet
0.8 (0.019)
0.77 (0.03)
0.72 (0.017)
0.73 (0.027)
0.67 (0.052)
(0.031)
VMPLN
0.71 (0.02)
0.71 (0.021)
0.71 (0.019)
0.71 (0.012)
0.71 (0.014)
0.69 (0.016)
Glasso
0.65 (0.037)
0.63 (0.034)
0.6 (0.03)
0.72 (0.018)
0.7 (0.018)
0.65 (0.021)
3000
EM-MMLE
0.93 (0.015)
0.93 (0.016)
0.93 (0.017)
0.95 (0.008)
0.93 (0.01)
0.88 (0.019)
Pl.Net
0.88 (0.014)
0.85 (0.023)
0.81 (0.03)
0.85 (0.037)
0.82 (0.025)
0.75 (0.052)
VMPLN
0.76 (0.024)
0.75 (0.024)
0.76 (0.025)
0.77 (0.014)
0.77 (0.015)
0.77 (0.012)
Glasso
0.69 (0.028)
0.67 (0.034)
0.66 (0.044)
0.79 (0.018)
0.78 (0.019)
0.74 (0.02)
Zero-proportion
High
Dimension
p 100 Middle
p 300 Middle
Mixing degree
Low
High
Low
High
1800
EM-MMLE
0.79 (0.019)
0.79 (0.023)
0.77 (0.025)
0.74 (0.015)
0.7 (0.018)
0.64 (0.022)
PLNet
0.72 (0.028)
0.69 (0.031)
0.65 (0.035)
0.62 (0,033)
0.59 (0.031)
0.53 (0.03)
VMPLN
0.62 (0.021)
0.61 (0.034)
0.61 (0.028)
0.59 (0.019)
0.57 (0.017)
0.57 (0.021)
Glasso
0.5 (0.023)
0.48 (0.035)
0.46 (0.037)
0.55 (0.015)
0.54 (0.019)
0.5 (0.025)
3000
EM-MMLE
0.88 (0.012)
0.89 (0.022)
0.87 (0.019)
0.87 (0.01)
0.84 (0.015)
0.79 (0.015)
PLNet
0.81 (0.023)
0.79 (0.033)
0.75 (0.03)
0.79 (0.034)
0.74 (0.039)
0.69 (0.036)
VMPLN
0.67 (0.025)
0.67 (0.036)
0.66 (0.029)
0.66 (0.015)
0.65 (0.019)
0.65 (0.017)
Glasso
0.53 (0.028)
0.52 (0.042)
0.49 (0.037)
0.64 (0.018)
0.62 (0.024)
0.59 (0.024)
Table Comparisons of EM-MMLE with PLNet, VMPLN and Glasso terms of AUPR on simulation results for scale-free graphs generated by the MPLN model. The results are averages over 50 replicates with standard deviations in brackets.
Zero-proportion
Low
Dimension
p 100 Middle
300 Middle
Mixing degree
Low
High
Low
High
1800
EM-MMLE
0.74 (0.036)
0.71 (0.035)
0.64 (0.05)
0.64 (0.034)
0.58 (0.04)
0.5 (0.04)
PLNet
0.68 (0.044)
0.6 (0.042)
0.54 (0.043)
0.57 (0.049)
0.49 (0.046)
0.43 (0.045)
VMPLN
0.54 (0.025)
0.53 (0.03)
0.52 (0.029)
0.48 (0.028)
0.48 (0.029)
0.43 (0.042)
Glasso
0.52 (0.038)
0.48 (0.037)
0.44 (0.04)
0.53 (0.035)
0.49 (0.037)
0.42 (0.033)
3000
EM-MMLE
0.84 (0.029)
0.85 (0.035)
0.78 (0.045)
0.83 (0.025)
0.78 (0.031)
0.68 (0.033)
PLNet
0.78 (0.045)
0.75 (0.053)
0.67 (0.042)
0.78 (0.042)
0.71 (0.05)
0.63 (0.04)
VMPLN
0.59 (0.027)
0.59 (0.025)
0.57 (0.034)
0.59 (0.021)
0.58 (0.018)
0.56 (0.023)
Glasso
0.56 (0.044)
0.56 (0.037)
0.5 (0.038)
0.66 (0.031)
0.62 (0.03)
0.54 (0.03)
Zero-proportion
High
Dimension
p 100
300 Middle
Mixing degree
Low
Middle
High
Low
High
1800
EM-MMLE
0.56 (0.037)
0.55 (0.033)
0.49 (0.047)
0.37 (0.039)
0.32 (0.04)
0.28 (0.043)
PLNet
0.5 (0.04)
0.47 (0.035)
0.41 (0.032)
0.33 (0.04)
0.29 (0.039)
0.25 (0.043)
VMPLN
0.42 (0.026)
0.42 (0.025)
0.4 (0.031)
0.3 (0.026)
0.28 (0.029)
0.27 (0.036)
Glasso
0.36 (0.032)
0.35 (0.033)
0.31 (0.031)
0.31 (0.032)
0.29 (0.034)
0.25 (0.035)
3000
EM-MMLE
0.7 (0.032)
0.7 (0.031)
0.65 (0.052)
0.58 (0.046)
0.53 (0.05)
0.47 (0.073)
PLNet
0,63 (0.035)
0.61 (0.034)
0,53 (0.043)
0.54 (0.042)
0.49 (0.047)
0.42 (0.067)
VMPLN
0,46 (0.02)
0.46 (0.023)
0.44 (0.02)
0.39 (0.025)
0.38 (0.028)
0.37 (0.043)
Glasso
0.4 (0.038)
0.39 (0.038)
0.35 (0.032)
0.42 (0.034)
0.41 (0.037)
0.35 (0.056)
29

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
Silver standards are derived from the 3' batch data, with gene pairs from public gene regulatory network databases considered potential regulatory relationships. Each identified regulatory relationship involves least transcription factor. For each cell type in the 3' batch, Spearman's p correlation is calculated between gene pairs with potential regulatory connections. Gene pairs showing significant Spearman's 's p correlations are designated as true regulatory relationships and included in the silver standard edge set for the respective cell type.
C.2. Description of the AUPRC Ratio
First, given network estimation ② from an algorithm, we define an edge score for each edge. For EM-MMLE, VMPLN, PLNet, and Glasso, the edge score for the edge (i,j) is defined as:
| (G,
ky
For PPCOR and GENIE3, the edge score is defined as the estimated connected weight for each edge.
Since the network inferred by the method contains connected edges with varying scores and unconnected edges, we calculate the area under the partial precision-recall curve (AUPRC) by applying different thresholds to the edge scores. To mitigate the impact of varying network densities across methods, define the AUPRC ratio as the ratio between the AUPRC of given method and the expected AUPRC of random network prediction. The precision of the random predictor is the edge density of the ground-truth network.
C.3. Supporting Figures
Figure 3 shows the gene regulatory networks inferred by EM-MMLE for the top 300 highly variable genes. ID3 is identified as cell-type-speciiee hub gene in the CD4+ T cell network. Figure 4 presents the top 10 biological processes from gene ontology analysis of ID3 target genes, highlighting the most significant biological processes associated with these genes.
30

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
092522
TRAE-2
Peum
WYLL
CISA
TEATS
JSST
S10099
UZ
S5410-11
VITE
12/6444
C377
TEVV
T8412-2
NIA-2MA
TREV12-3
STAN
SECTIC
ITACC
TUBA
HGA D0ka
HLA-D0A
12822-33
200
119/3-1
Tayn 1-11--2
NLA-008"
(a) cells
(b) Monocytes
KUC
ACTE
GERCIO
300W
DECL
12J7-44
K20/D-33
TURR
FLAI DNN)
STAIK
MISTINMN
KO00-13 VAGE
ONL
TUIKOR
KRVID-3
PTTOT
000
1197-2
TUSIFF
KCAM
VKST
ORR5
ICR00-12
CEAL
RAU
HIA CR
080-001
TRADOS
HISTHLC
KHHN
DESLE
1212-30
Cupars
COCa
IORO--O GUO-TE
CCS
CCUL2
CORVI
WAT
Uni2g1372
TEPG-1
TBW220K
TRUEQI-1
100g-11 OUCASS
(c) CD8+ T cells
(d) CD4+ cells
Figure 3. Inferred gene regulatory networks for four cell types. The size of each node represents its weighted node degree, while the edge width indicates correlation weight. Genes highlighted represent the transcription factors interest.
31

---

A Generic Family of Graphical Models: Diversity, Efficiency, and Heterogeneity
antigen processing presentation of exogenous
antigen
antigen processing and presentation of exogenous
peptide antigen via MHC
class I!
antigen processing and presentation of endogenous
antigen
antigen processing and presentation peptide antigen via MHC class antigen processing and presentation of peptide or polysaccharide antigen
Adjusted p-value
MHC c ass
2e-05
lymphocyte mediated immunity
antigen processing and presentation of exogenous peptide antigen
antigen processing and presentation
eukocyte mediated immunity
MHC class protein comp ex assembly
Count
Figure 4. Gene ontology enrichment analysis of ID3 target genes within gene regulatory network CD4+ Tells.
32


Reviews:
[{"review_id": "0mXqIpxilY", "reviewer": "Reviewer_7ySj", "summary": "The paper proposes a new, unified class of graphical models—called the marginally recoverable parametric family—designed to handle diverse data types (e.g., Gaussian, Poisson log‑normal, and latent Gaussian copula models) and heterogeneous structures via mixture modeling. The authors introduce an efficient maximum marginal likelihood estimator (MMLE) that avoids full high‑dimensional integration by leveraging low‑dimensional marginal computations, and extend it to EM‑MMLE for mixture contexts. The work is supported by extensive theory, simulations, and an application to single‑cell RNA‑seq gene regulatory network inference.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your time in reviewing our paper and your insightful comments. In the following response, we would like to address your major concern and provide additional clarification.\n\nFor your major concern in the claims and evidence section, as well as the similar concern raised in the theoretical claims and weaknesses sections:\n\n> The theoretical guarantees hinge on two stringent conditions (the lower-boundedness and upper-boundedness conditions in Conditions 4.1 and 4.2). Their restrictive nature may limit applicability in practical scenarios where the data do not closely conform to these idealized requirements.\n\nThank you for your comment. We appreciate your concern about the potential restrictiveness of Conditions 4.1 and 4.2 and their applicability to real-world data. We would like to clarify that Conditions 4.1 and 4.2 are, in fact, relatively mild. As demonstrated in Appendix A.3.4 and Appendix A.4, both the MPLN distribution and the latent Gaussian copula model for binary data satisfy these conditions.\n\nMore specifically, our proof establishes that a distribution satisfies Conditions 4.1 and 4.2 as long as the following assumptions hold:\n\n(A1) The Fisher information matrix is positive definite and continuous with respect to the parameters.\n\n(A2) Both the score function and the Hessian matrix are dominated by integrable functions.\n\n(A3) The parameter space is compact, and the model is identifiable.\n\nAmong these, (A1)-(A3) are standard regularity conditions, which are widely used to ensure the asymptotic efficiency of the maximum likelihood estimator, as shown in Shao's book [1].\n\nFor the other weakness:\n\n> More experiments showing robustness under deviations from the assumed conditions would be valuable.\n\nThank you for your comment. Due to rebuttal time constraints and the fact that we have not identified any marginally recoverable distributions that do not satisfy the conditions, we have not conducted additional experiments to demonstrate robustness under deviations from the assumed conditions at this time.\n\nReference\n\n[1] Shao, J. Mathematical statistics. Springer Science \\& Business Media, 2008.", "claims_and_evidence": "### Main Claims:\n\n- The proposed framework unifies and generalizes a range of graphical models while reducing computational complexity.\n\n- MMLE (and its EM extension) achieves consistent parameter estimation and sign recovery with nearly optimal convergence rates.\n\n### Major Concern:\n\nThe theoretical guarantees hinge on two stringent conditions (the lower‑boundedness and upper‑boundedness conditions in Conditions 4.1 and 4.2). Their restrictive nature may limit applicability in practical scenarios where the data do not closely conform to these idealized requirements.", "methods_and_evaluation_criteria": "The estimation strategy first maximizes one‑dimensional marginal likelihoods to estimate location and scale parameters, then refines covariance estimation by using two‑dimensional marginal likelihoods. This decomposition reduces the burden of high‑dimensional integration. The extension to mixture models is handled via an EM algorithm updating the MMLE.\n\nPerformance is measured using standard network inference metrics (e.g., AUPR, Frobenius norm error) and clustering accuracy (Adjusted Rand Index). The study covers various simulation settings—including different network structures, dimensions, mixing levels, and zero proportions—to assess robustness.", "theoretical_claims": "The paper provides detailed proofs—offered in the supplementary material—for key results such as convergence rates (Theorem 4.3) and sign consistency (Theorem 4.4).\n\nConcern: The reliance on Conditions 4.1 and 4.2 is a potential weakness. These assumptions, pivotal to the derived guarantees, appear rather strong and may not hold in many real‑world applications.", "experimental_designs_or_analyses": "### Simulation Studies:\n\nSimulations are performed on mixed count and binary data generated from realistic distributions (e.g., MPLN for count data and latent Gaussian copula for binary data).\n\n### Real‑Data Application:\n\nThe method is applied to scRNA‑seq data to infer gene regulatory networks, compared against silver standards from public databases.", "relation_to_broader_scientific_literature": "The work builds on several key strands of graphical model research—extending classical Gaussian approaches to handle count and binary data—and connects with recent work on latent Gaussian copula models.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "### Strengths:\n\n- Presents a versatile, unified framework that handles multiple data types and heterogeneous populations.\n\n- Develops an innovative estimation strategy (MMLE/EM‑MMLE) that significantly reduces computational complexity.\n\n- Offers comprehensive theory and extensive simulations alongside a compelling real‑data application in genomics.\n\n### Weaknesses:\n\n- The two core conditions (Conditions 4.1 and 4.2) are quite strong, potentially restricting the method’s applicability.\n\n- More experiments showing robustness under deviations from the assumed conditions would be valuable.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "See my previous comments.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FYmrJsFBzV", "reviewer": "Reviewer_fTvr", "summary": "This article introduces a new class of graphical models, referred to as the marginally recoverable parametric family, which aims to tackle challenges related to efficiency and heterogeneity in high-dimensional network inference tasks. The proposed family is rather flexible, with the joint distribution characterized by pairwise marginal distributions. To facilitate efficient parameter estimation, the paper presents a Maximum Marginal Likelihood Estimator, which is further extended to handle heterogeneous scenarios through mixture modeling. The method’s effectiveness is supported by theoretical consistency, simulations, and applications to real-world sequencing data.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your time in reviewing our paper and your insightful comments. In the following response, we would like to address your major concern and provide additional clarification.\n\nFor theoretical claims:\n\n> I believe that neither Theorem 4.2 nor Theorem 4.3 should be labeled as \"Theorems\". Instead, I would suggest calling them \"Lemmas,\" as they directly follow from other established theorems.\n\nThank you for your valuable suggestion. Since these two theorems are derived from established results, we acknowledge that referring to them as \"Lemmas\" is more appropriate. Nevertheless, we would like to emphasize that these results establish that once a connection between the parameter distance and the density distance is made (as specified in Conditions 4.1 and 4.2), we can obtain tail probability control for high-dimensional parameters. We believe that this provides new insights and a novel proof strategy. Furthermore, verifying that a specific model satisfies the necessary conditions (Conditions 4.1 and 4.2) is nontrivial, especially for complex models such as MPLN. In fact, a substantial portion of the appendix (see Appendix A.3) is dedicated to rigorously demonstrating that the MPLN model satisfies these conditions. \n\nFor weaknesses:\n\n> While the family is presented in a general form, it appears to only include elliptical copula distributions (and simple compositions thereof). This raises some concerns about the novelty of the approach, as elliptical copula families are already well-established in the literature. Additionally, elliptical copulas are quite limited in their expressive power, particularly when it comes to modeling tail dependencies, which either vanish asymptotically or are entirely absent in the special case of the Gaussian copula. In many applications of copulas, this limitation has led to the development of more complex copula families. Therefore, the claim on page 2 that \"The marginally recoverable parametric family includes many of the most common distributions\" might be somewhat misleading.\n\nThanks for your comment. In fact, the family includes distributions beyond elliptical copula distributions (and simple compositions thereof). For instance, the multinomial distribution belongs to the family but does not fall within the class of elliptical copula distributions. Specifically, let $\\mathbf{X}=\\left(X_{1}, \\ldots, X_{p}\\right)^\\top \\sim Multinomial(n,\\boldsymbol{\\theta})$, where $\\boldsymbol{\\theta}=(\\theta_1,\\dots,\\theta_p)^\\top$. Then, for $1 \\leq j < k \\leq p $, we have $X_j \\sim Multinomial(n, \\theta_{j}) $ and $\\left(X_j, X_k\\right)^\\top \\sim Multinomial(n, (\\theta_j,\\theta_k)^\\top ) $. This demonstrates that the multinomial distribution satisfies the definition of the family. \n\nFurthermore, the concept of the marginally recoverable family can be extended to a more general framework. We introduce the definition of a $d$-marginally recoverable family as follows:\n\nLet $\\mathbf{X}$ be a $p$-dimensional random variable, and let $d$ be a constant such that $d<p$. We say that $\\mathbf{X}$ is $d$-marginally recoverable if any $d$-dimensional marginal distribution of $\\mathbf{X}$ belongs to the same distribution family, and the parameters of all $d$-dimensional marginal distributions collectively characterize the parameters of the full distribution. \n\nIn this way, the marginally recoverable family becomes more flexible. Under this definition, the multinomial distribution mentioned earlier is actually 1-marginally recoverable, as only the parameters of the 1-dimensional marginal distributions are required to recover all the parameters of the full distribution. We will introduce this general definition in the Discussion section of a future revision and explore additional examples that align with this broader definition in our future work.\n\nFor modeling tail dependencies, we acknowledge that elliptical copulas do have limitations in this regard. That said, this concern is primarily relevant in financial applications. Our work, on the other hand, focuses on network inference and its applications to biological genomics data, where the elliptical copula family is commonly used, and tail dependencies are generally not a central issue. Exploring tail dependencies within our framework is an interesting direction for future research. The claim on page 2 could be revised to: \"The marginally recoverable parametric family includes distributions commonly used in the field of bioinformatics,\" which would be a more precise statement. We sincerely appreciate your constructive suggestion.", "claims_and_evidence": "The theoretical results are well-supported, though not particularly surprising, as they follow directly from well-established theorems. The evidence provided is both clear and convincing.\nWhile the marginally recoverable parametric family introduced is intriguing, I remain somewhat skeptical that it encompasses distributions beyond elliptical copulas and simple hierarchical compositions that do not alter the dependency structure.", "methods_and_evaluation_criteria": "On the experimental side, the model's performance is assessed through various simulation benchmarks and large-scale real-world applications. All evaluation procedures appear to be appropriate and well-executed.", "theoretical_claims": "The consistency proofs are mathematically sound and based on widely recognized statistical methods. However, I believe that neither Theorem 4.2 nor Theorem 4.3 should be labeled as \"Theorems\". Instead, I would suggest calling them \"Lemmas,\" as they directly follow from other established theorems.", "experimental_designs_or_analyses": "Experimental results show that the proposed method performs well on benchmark datasets. The simulation studies are carefully designed, addressing diverse experimental settings and data heterogeneity. The real-world applications are both relevant and well-motivated, with analyses supported by biological literature.", "relation_to_broader_scientific_literature": "The family introduced in this paper is versatile, making the proposed approach applicable to several problems related to graphical models. However, it is also relatively limited in its ability to model tail dependencies, which might be problematic in some applications.", "essential_references_not_discussed": "No.", "other_strengths_and_weaknesses": "The primary contribution of the paper lies in the introduction of the marginally recoverable family. While the family is presented in a general form, it appears to only include elliptical copula distributions (and simple compositions thereof). This raises some concerns about the novelty of the approach, as elliptical copula families are already well-established in the literature. Additionally, elliptical copulas are quite limited in their expressive power, particularly when it comes to modeling tail dependencies, which either vanish asymptotically or are entirely absent in the special case of the Gaussian copula. In many applications of copulas, this limitation has led to the development of more complex copula families. Therefore, the claim on page 2 that \"The marginally recoverable parametric family includes many of the most common distributions\" might be somewhat misleading.", "other_comments_or_suggestions": "No.", "questions_for_authors": "Please comment on my concerns regarding the \"real\" flexibility of the model beyond elliptical copulas.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0HhA4aCfCJ", "reviewer": "Reviewer_HnXb", "summary": "The paper introduces a novel family of graphical models, termed the marginally recoverable parametric family, to address diversity, efficiency, and heterogeneity in high-dimensional network inference. The paper proposes a Maximum Marginal Likelihood Estimator (MMLE) for efficient parameter estimation and extend it to EM-MMLE to handle heterogeneous membership via mixture modelling. The approach is demonstrated through theoretical consistency, simulations, and application to real single-cell RNA sequencing data.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your time in reviewing our paper and your positive feedback. Here we address your comments as follows.\n\nFor weakness:\n>Marginal recoverability may not extend to distributions with higher moments.\n\nThanks for your comment. To facilitate an intuitive network representation, we impose certain restrictions on the parameters of the family, typically related to the first and second moments. However, motivated by the core ideas behind this family, we can provide a more general definition that extends to distributions determined by higher moments. Specifically, we introduce the definition of a $d$-marginally recoverable family. Let $\\mathbf{X}$ be a $p$-dimensional random variable, and let $d$ be a constant such that $d<p$. We say that $\\mathbf{X}$ is $d$-marginally recoverable if any $d$-dimensional marginal distribution of $\\mathbf{X}$ belongs to the same distribution family, and the parameters of all $d$-dimensional marginal distributions collectively characterize the parameters of the full distribution. We will introduce this general definition in the Discussion section of our revision and explore additional distribution examples that align with it in our future work.\n\nFor suggestions:\n>1. Notation $d$ is used for both dimension and Hellinger distance in Condition 4.1 and 4.2.\n\nThank you for your helpful suggestions. We have corrected the duplicate symbols in the manuscript by replacing $d$ (used for dimension in Conditions 4.1 and 4.2) with $k$.\n\n>2. It is very hard to interpret the sample size requirement in Theorem 4.4. Some discussion of sample complexity depending on the dimension / sparsity will be appreciated.\n\nFor some $\\eta >2$, the result in Theorem 4.4 holds under the sample size condition $n>O(\\eta s^{1/2}\\log p)$, where $n$ is the sample size, $s$ denotes the sparsity level of the network (i.e., the number of nonzero entries), and $p$ denotes the dimension. Since $\\widehat{\\Theta}_g$ recovers all the zeros and nonzeros in $\\Theta_g$ with probability $1-p^{2-\\eta}$, increasing $\\eta$ raises this probability towards 1, but requires a larger sample size.\n\n>3. It would be helpful to provide an example with both continuous and discrete variables with explicit modeling formula.\n\nHere, we provide a specific example involving both continuous and discrete variables, with an explicit modeling formula. Assume that $\\mathbf{X}=(\\mathbf{X}_1,\\mathbf{X}_2)$, where $\\mathbf{X}_1$ is a random $a$-vector and $\\mathbf{X}_2$ is a random $b$-vector. Suppose that there is a random vector $\\mathbf{Z}_1=(Z_1,...,Z_a)$ such that $(\\mathbf{Z}_1,\\mathbf{X}_2)\\sim N(\\mu,\\Sigma)$ and $X_j=I(Z_j>C_j)$ for all $j=1, ..., a$, where $C_j$ is a constant. We can easily verify that $\\mathbf{X}$ is marginally recoverable. While this paper does not primarily focus on mixed-type data, our approach can be applied to parameter estimation as long as the variables meet the marginal recoverability definition.\n\nFor questions:\n>1. How do we justify Condition 4.1 and 4.2 intuitively?\n\nConditions 4.1 and 4.2 establish a relationship between the parameter distance and the distance between the marginal density functions. Intuitively, for any two parameters within the parameter space, if their Euclidean distance is small, the corresponding Hellinger distance between their marginal density functions will also be small. When verifying that specific distributions satisfy Conditions 1 and 2, we found that this is equivalent to verifying the following conditions:\n\n(A1) The Fisher information matrix is positive definite and continuous with respect to the parameters. \n\n(A2) Both the score function and the Hessian matrix are dominated by integrable functions. \n\n(A3) The parameter space is compact, and the model is identifiable. \n\nAmong these, (A1)-(A3) are standard regularity conditions. However, expressing the original conditions as (A1)-(A3) would appear somewhat redundant. Conditions 4.1 and 4.2, on the other hand, provide a more intuitive link between the distances of densities and parameters, making the formulation clearer and more accessible. Therefore, we use these boundedness conditions to present the core conditions of our theorem instead of (A1)-(A3).\n\n>2. Is there any theoretical guarantee for the consistency of membership assignment in the mixture modeling?\n\nFor the Gaussian mixture model, Chen and Zhang [1] derive minimax lower bounds for clustering with respect to the misclustering error rate. To the best of our knowledge, establishing theoretical guarantees for membership assignment remains a challenging problem. Due to time constraints, we have not addressed the consistency of membership assignment in general mixture models in this paper, but this could be explored in future work.\n\nReference\n\n[1] Chen, X. and Zhang, A. Y. Achieving optimal clustering in gaussian mixture models with anisotropic covariance structures. In The Thirty-eighth Annual\nConference on Neural Information Processing Systems, 2024.", "claims_and_evidence": "Most of the claims are theoretical statements and are verified via proofs.", "methods_and_evaluation_criteria": "The proposed marginally recoverable parametric family is well-motivated and mathematically rigorous. It is evaluated via theoretical consistency, simulation benchmarks (e.g. AUPRC) and application to real RNA-seq datasets.", "theoretical_claims": "The consistency proofs for MMLE (Theorem 4.3) and EM-MMLE (Theorem 4.4) appear mathematically sound. The proofs rely on Conditions 4.1 and 4.2, which ensure proper convergence.", "experimental_designs_or_analyses": "The simulation studies are well-structured, covering different aspects of experiment settings and data heterogeneity. The real-world application (scRNA-seq) is relevant and well-motivated.", "relation_to_broader_scientific_literature": "Graphical modelling for network analysis is widely used in scientific research, especially in genetics and epidemiology. Developing flexible and consistent methods is significant for the exploratory analysis.", "essential_references_not_discussed": "Not sure.", "other_strengths_and_weaknesses": "**Strength**:\n- This paper introduces a flexble class of graphical models with rigorous mathematical justification.\n- The advantage of this class is to avoid high-dimensional integration via marginal likelihood estimation.\n- The experiments show the proposal outperforms existing methods on benchmark datasets.\n\n**Weakness**:\n- Marginal recoverability may not extend to distributions with higher moments.", "other_comments_or_suggestions": "- Notation $d$ is used for both dimension and Hellinger distance in Condition 4.1 and 4.2.\n- It is very hard to interpret the sample size requirement in Theorem 4.4. Some discussion of sample complexity depending on the dimension / sparsity wil be appreciated.\n- One motivation of this paper is for mixed type data analysis. It would be helpful to provide an example with both continuous and discrete variables with explicit modelling formula.", "questions_for_authors": "- How do we justify Condition 4.1 and 4.2 intuitively? I can find parametric examples in Line 275-278.\n- Is there any theoretical guarantee for the consistency of membership assignment in the mixture modelling?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZD0PC9OA78", "reviewer": "Reviewer_oABz", "summary": "To infer high dimensional graphical models with a variety of data types, this paper introduces a marginally recoverable parametric family. The family is very flexible, and it includes Gaussian, PLN, and latent Gaussian copula for count and binary data. Within this family, the joint distribution can be characterized by all pairwise marginal distributions, and hence the dimensionality is reduced significantly. To further capture heterogeneous structures, they extend the method to mixture models. The proposed methods are evaluated via simulations and single-cell RNA sequencing data.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate your time in reviewing our paper and your insightful comments. In the following response, we would like to address your major concern and provide additional clarification.\n\nFor weakness:\n\n> The author(s) write the family in a very general way, but the examples are more or less all Gaussian related. Therefore, there are some concerns about novolties, as the defined family may be some extensions of Gaussian and we can directly use results.\n\nThanks for your comment. First, since the focus of this work is primarily on graphical network inference, with the application scenario being gene regulatory network inference using biological genomics data, the examples provided here are mainly based on Gaussian distribution models in hierarchical settings. However, the new concept introduced—the marginally recoverable family—includes distributions beyond just Gaussian-related distributions. For instance, the multinomial distribution belongs to this family but does not Gaussian-related. Let $\\mathbf{X}=\\left(X_{1}, \\ldots, X_{p}\\right)^\\top \\sim Multinomial(n,\\boldsymbol{\\theta})$, where $\\boldsymbol{\\theta}=(\\theta_1,\\dots,\\theta_p)^\\top$. Then, for $1 \\leq j < k \\leq p $, we have $X_j \\sim Multinomial(n, \\theta_{j}) $ and $\\left(X_j , X_k \\right)^\\top \\sim Multinomial(n, (\\theta_j, \\theta_k)^\\top) $. This demonstrates that the multinomial distribution satisfies the definition of the family. \n\nSecondly, our defined family is not an extension of the Gaussian distribution that can be addressed by existing methods. Even when the latent layer is modeled using a Gaussian distribution, distributions with hierarchical structures still pose significant challenges for parameter estimation. For example, as mentioned in the paper, there are few methods developed for estimating the precision matrix of the MPLN model in high-dimensional settings. Compared to Gaussian mixture models, the MPLN model presents a greater challenge during the parameter estimation process using the EM algorithm, as its expected complete log-likelihood involves high-dimensional integrals that are difficult to compute. Therefore, for such models, we cannot directly apply existing methods or results.\n\nFor other comments or suggestions:\n\n> Figure 1 is blank.\n\nThank you for your comment. We appreciate your attention to Figure 1. We would like to clarify that Figure 1 is not blank. It is possible that the white areas in the figure were confused with the background, which might have led to the impression that the figure is empty. To avoid this misunderstanding, we will add borders around each subplot in Figure 1 in the final version. We would like to clarify that in subplot (a), the white areas represent positions where the true network does not have edges. In subplots (b) to (e), the white areas at position $(i,j)$ indicate that the frequency of false positives (i.e., incorrectly predicting an edge where there is none) for that position is zero. We hope this explanation clarifies the meaning of the white areas.", "claims_and_evidence": "Yes, the claims are supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "Yes, the author study their methods by both simulations and single-cell RNA data. They evaluate the model performance using diverse criteria for different data types, e.g., AUPR for count data and 3 estimation methods for binary data. All these methods and criteria are standard and valid.", "theoretical_claims": "No, I didn't, since these are well-established statistical methods and I believe the authors did it correctly.", "experimental_designs_or_analyses": "Yes. All the simulations and application to single-cell RNA data is solid. For application in singel-cell RNA data, they further provide biological perspective with supporting biological literature, which make the analyses more meaningful.", "relation_to_broader_scientific_literature": "The family defined in this paper is general, which contains many widely-used distributions. Therefore, we can use the technique proposed here to models (not limited to graphical models) within this family, to deal with high dimensionality issue.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "Strength: the writing is clear and easy to follow. The application part provides biological insights.\n\nWeakness: The major contribution for this paper is proposing a new marginally recoverable family. The author(s) write the family in a very general way, but the examples are more or less all Gaussian related. Therefore, there are some concerns about novolties, as the defined family may be some extensions of Gaussian and we can directly use results. For mixture part, overlaiding the MMLE with EM is standard.", "other_comments_or_suggestions": "Figure 1 is blank.", "questions_for_authors": "NA", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: All reviewers assigned (mostly weakly) positive scores. They all agreed that the underlying research question of designing flexible models for handling heterogeneous data is interesting and relevant, that  the proposed framework of the "marginally recoverable family" of copula distributions has contains some novel inference ideas, and that the experiments conducted are rather convincing. On the more negative side, reviewers mentioned potential limitations of the marginally recoverable family concerning the restriction to Gaussian-like distributions and problems with handling tail dependencies in data, and some missing experiments that would  indicate how robust the proposed framework is when some of the theoretical assumptions are violated. Some of these concerns, however,  could be addressed (at least to some degree) in the rebuttal. In summary, I think that for this paper, the positive aspects outweigh the negative ones.
Decision: Accept (poster)