Epidemiological Causal Graph Identification: Challenges, Identifiability and Algorithms
arXiv:2609.20676v1 [cs.LG] 17 Sep 2026
Sambit Mishra∗ , Yingying Wang† , Christine K. Johnson† , and Urbashi Mitra∗ ∗ University of Southern California, † University of California, Davis
Abstract—Causal discovery from observational data is fundamental to statistics and machine learning, yet determining causal direction without interventions necessitates structural assumptions. Existing identifiability research primarily focuses on continuous variables under additive noise models, often neglecting mixed datasets containing ordinal scales, counts, and continuous measurements. This paper investigates causal discovery in Directed Acyclic Graphs (DAGs) where nodes follow either an ordinal distribution (via an ordered logit model) or a regular one-parameter exponential family distribution. We prove that the edge direction between an ordinal and an exponential family node is distributionally identifiable for generic parameter values. Our findings generalize previous Ordinal-Poisson results to the broader exponential family. Computationally, we introduce a score-based exhaustive search and a masked continuous optimization framework using DAGMA for larger graphs. Numerical results validate the theory, recovering edge orientations within a Markov equivalence class that are unidentifiable under classical structural equation models. Index Terms—Causal Discovery, Directed Acyclic Graphs, Identifiability, Exponential Family, Ordinal Distribution, Structure Learning
I. I NTRODUCTION Causal reasoning distinguishes the variables that produce an outcome from those that merely correlate with it. Causal directed acyclic graphs (DAGs) [1] compactly encode these relationships, using nodes for random variables and directed edges for causal influence. We focus on discovery from observations, where the key challenge is that multiple DAGs can induce the same joint distribution, forming a Markov equivalence class (MEC) [2]. Constraint-based algorithms such as PC [3] recover the MEC but cannot orient edges outside Vstructures; going beyond the MEC requires model assumptions that break this symmetry. Identifiability has been shown for linear Gaussian structural equation models (SEMs) with equal noise variances [4], linear non-Gaussian SEMs [5], and nonlinear additive noise models [6]. This prior work focuses on continuous valued random variables and additive noise models, a mold that real observational datasets rarely fit. Our current work is motivated by causal discovery in epidemiology, where ordinal, binary, and counting data appear alongside continuous-valued measurements [7], [8], [9]. Forcing these into a Gaussian additive-noise model introduces mis-specifications that invalidate both identifiability guarantees and downstream inferences. A key contribution is the structured statistical model (SSM), which generalizes structural equation models and appears to induce stronger identifiability properties. We leverage prior
identifiability work for specific discrete families: Poisson DAGs [10], extended to the generalized hypergeometric family in [11], and ordinal DAGs [12]. We generalize our prior work in [13], which considered mixed ordinal and Poisson nodes. We also adapt key algorithms to the SSM. Since the DAG space grows super-exponentially with the number of nodes [14], methods such as NOTEARS [15], GOLEM [16], DAGMA [17], and SURE-Ridge [18] reformulate DAG learning as a continuous optimization problem with a smooth acyclicity characterization or regress-and-threshold style DAG learning, but focus on linear or additive noise models. The contributions of this work are as follows: 1) We introduce the structured statistical model (SSM), which generalizes structural equation models. 2) We prove that the edge between an ordinal node and a regular one-parameter exponential family node in a bivariate SSM is distributionally identifiable for generic parameter values, subsuming [13]. 3) We develop a masked DAGMA procedure with joint cutpoint and weight optimization for graphs beyond the reach of exhaustive search. 4) We validate the framework on three-node and 50-node experiments, where the normalized structural Hamming distance (SHD) converges to zero with sample size, including correct orientation of edges within the MEC. II. BACKGROUND We consider a causal directed acyclic graph (DAG) with d nodes, denoted by G = (V, E), where V = {1, . . . , d} is the vertex set and E = {(i, j) : i → j} is the edge set. The graph G is represented by a deterministic weighted adjacency matrix WG ∈ Rd×d whose (i, j) entry is nonzero if and only if (i, j) ∈ E. Each node i is associated with a random variable Xi and parent set P (i) = {k : (k, i) ∈ E}, and we write x = (X1 , . . ., Xd ). For each node, the conditional distribution p Xi | xP(i) belongs to a known parametric family whose parameters are deterministic functions of the parent variables. This generalizes the structural equation models typically seen in causal inference; we call it our structured statistical model (SSM). The observed data matrix X ∈ Rn×d contains n independent realizations of x. A. Ordinal Nodes Let i ∈ V correspond to an ordinal random variable Xi with finite ordered support Xi = {1, . . . , s}, s ≥ 2. If P (i) = ∅, Xi is categorical with p (Xi = x) = πi,x . Otherwise, we
adopt a cumulative-link (ordinal regression) model. Let wi := [WG ]P(i), i denote the causal coefficients from the parents of node i, with parent vector xP(i) , and introduce strictly ordered cutpoints −∞ = γi,0 < γi,1 < · · · < γi,s−1 < γi,s = +∞. The conditional probability mass function is p Xi = x | xP(i) = f γi,x − wiT xP(i) − f γi,x−1 − wiT xP(i) . (1)
sufficient statistic T (y) is one of y, log (y), or y k with fixed k > 0.
Throughout this work, we adopt the ordered logit model [19], f (u) = σ (u) = 1/ (1 + e−u ); location shifts are absorbed into the cutpoints.
is an affine function of T (y) on the support of Y .
B. Exponential-Family Nodes Let i ∈ V correspond to a random variable Xi whose conditional distribution belongs to a regular one-parameter exponential family with known sufficient statistic Ti : Xi → R and natural parameter space Ω ⊆ R. The support Xi ⊆ R may be discrete or continuous. If P (i) = ∅, the natural parameter is a fixed ηi ∈ Ω. Otherwise, we model it as a deterministic function of a linear predictor: wi := [WG ]P(i), i , ηi xP(i) = gi wiT xP(i) , where gi : R → Ω is a known, strictly monotone link function. The conditional distribution becomes pXi |xP(i) (x) = , Hi (x) exp ηi xP(i) Ti (x) − ai ηi xP(i)
(2)
where Hi (x) > 0 is the base measure and ai (·) is the logpartition function. III. I DENTIFIABILITY We begin with a two-node system with random variables X and Y . We study two competing causal models MX→Y and MY →X , with the edge weight of the single directed edge between X and Y being w ̸= 0. We assume that X is an ordinal random variable with support X = {1, . . . , s}, and that the conditional distribution of Y belongs to a regular one-parameter exponential family with known sufficient statistic T (y). Under MX→Y , X is a root node with probabilities πx and Y | X = x has natural parameter η (x) = g (wx), so that pX,Y (x, y; MX→Y ) = πx H (y) exp(η (x) T (y) − a (η (x))) . (3) Under MY →X , Y is a root node with natural parameter η and X | Y = y follows the ordered logit with cutpoints γ and linear predictor wy, so that pX,Y (x, y; MY →X ) = H (y) exp(ηT (y) − a (η)) (σ (γx − wy) − σ (γx−1 − wy)) . (4) We use the following assumptions throughout the work. Assumption 1. (i) X has at least three ordinal levels, s ≥ 3; (ii) the support of Y contains at least three distinct points for discrete Y and extends to +∞ for continuous Y ; (iii) the
Our key results, all under Assumption 1, are as follows. Lemma 1. For MX→Y and any distinct u, v ∈ X , pX|Y (u | y) log pX|Y (v | y)
Proof sketch. By Bayes’ rule and (3), the base measure H (y) cancels and the log-ratio equals (ηu − ηv ) T (y) + log (πu /πv ) − (a (ηu ) − a (ηv )), which is affine in T (y). Lemma 2. Under MY →X with the ordered logit cumulativelink model, for any distinct categories u, v ∈ X , the posterior ratio pX|Y (u | y) R (y; u, v) ≜ pX|Y (v | y) is a strictly positive rational function of ewy on the support of Y . Its leading-order behavior as ewy → 0+ and ewy → ∞ depends only on whether u or v lies at the boundary of X : if both are interior, R tends to a positive constant at both ends; if exactly one is a boundary category, R ∼ c e±wy at one end and tends to a positive constant at the other; if both are boundary, R ∼ c e±wy at both ends. The sign of the exponent is determined by which side carries the boundary category, and all constants are strictly positive and finite, depending only on the cutpoints γ. Proof sketch. With t = ewy , the ordered logit gives pX|Y (1 | y) = eγ1 / (eγ1 + t), pX|Y (s | y) = t/ (eγs−1 + t), and pX|Y (x | y) = (eγx − eγx−1 ) t/ [(eγx + t) (eγx−1 + t)] for interior x; boundary and interior categories have different leading orders in t at each end, and forming ratios yields the three cases. Lemma 3. There do not exist parameters (γ, w) for which pX|Y (u | y) log pX|Y (v | y) is affine in T (y), on a continuous interval of y, for any distinct u, v ∈ X . Furthermore, for discrete Y , this affine relationship cannot hold for generic parameters (γ, w). Proof sketch. Suppose R (y; u, v) = eβ eαT (y) on the support; if Y is continuous, real-analyticity extends this to a half-line. For T (y) = y, this reads R = eβ tr with r = α/w; matching leading-order exponents from Lemma 2 at t → 0+ and t → ∞ forces r ∈ {−1, 0, +1}. r = 0 forces coinciding cutpoints, contradicting u ̸= v; r = ±1 forces {u, v} = {1, s} and then eγ1 = eγs−1 , contradicting s ≥ 3. For T (y) = logy or y k (k ̸= 1), log R would grow as Θ (log y) or Θ y k , incompatible with the O (1) or Θ (y) growth of log R under Lemma 2, forcing α = 0 and reducing to the constant case. For discrete Y , three support points yield a collinearity condition Φ (γ, w) = 0 with Φ real-analytic and not identically zero, so it fails outside a Lebesgue-null set.
TABLE I R EGULAR O NE -PARAMETER E XPONENTIAL FAMILY D ISTRIBUTIONS Distribution
Natural Parameter η
Link Function η = g wiT xP(i)
Suff. Stat. T (x)
Support X
Exponential Poisson Gaussian (fixed variance σ 2 )
−λ, (λ > 0) log λ, (λ > 0) µ , (µ ∈ R) σ2
η = −ewi xP(i) η = wiT xP(i) η = wiT xP(i)
x x x
[0, ∞) {0, 1, . . . } R
Gamma (fixed shape α) Binomial (fixed trials m ≥ 3)
−β, (β > 0) p log 1−p , (p ∈ (0, 1))
η = −ewi xP(i) η = wiT xP(i)
x x
(0, ∞) {0, . . . , m}
Pascal (fixed successes r)
log(1 − p),
η = − log(1 + ewi xP(i) )
Gamma (fixed rate β)
(p ∈ (0, 1))
α − 1, (α > 0)
Theorem 1. The models MX→Y and MY →X are distributionally identifiable (generically identifiable when Y is discrete). Proof sketch. Suppose forward parameters (η, π) and reverse parameters (η, γ, w) yield the same joint, hence the same posteriors pX|Y (x | y). By Lemma 1, every log-posterior ratio is affine in T (y), which Lemma 3 rules out for w ̸= 0 (generically, for discrete Y ). Hence no such reverse parameters exist. IV. A LGORITHMS FOR M IXED DAG D ISCOVERY Algorithm 1 Masked DAGMA with m-NLL Score for Large DAG Discovery Require: Data matrix X, number of nodes d ≥ 6, bipartite partition (Vord , Vexp ), initial central path coefficient µ(0) (e.g., 1), decay factor α ∈ (0, 1) (e.g., 0.1), ℓ1 parameter λ > 0 (e.g., 0.01), log-det parameter τ > 0 (e.g., 1), number of iterations T , threshold ω > 0 (e.g., 0.3) ∗ ∗ induced by West Ensure: Estimated DAG Gest d×d 1: Construct bipartite mask M ∈ {0, 1} where [M]i,j = 0 if nodes i and j belong to the same type, and [M]i,j = 1 otherwise n o (0) (0) 2: Initialize West = 0d×d and γ i =0 i∈Vord 3: for t = 0, 1, 2, . . . T − 1 do n o (t) (t) , solve 4: Starting at West and γ i i∈Vord
n o (t+1) (t+1) West , γ i = arg min
W,{γ i }
where
µ(t) Q (W, {γ i } ; X) + h (W ⊙ M)
Q (W, {γ i } ; X) = e W ⊙ M, {γ i }i∈Vord ; X + λ∥W ⊙ M∥1
Set µ(t+1) = αµ(t) 6: end for (T ) (T ) 7: Apply bipartite mask to final iterate: West ← West ⊙M (T ) ∗ 8: Threshold small ientries: [West ]i,j ← [West ]i,j · h (T ) 1 |[West ]i,j | ≥ ω ∗ 9: return Gest 5:
Let V = Vord ∪ Vexp be partitioned into ordinal and exponential family nodes, with edges permitted only between
T
T
T
wiT xP(i)
η=e
−1
x
{0, 1, . . . }
log (x)
(0, ∞)
nodes of different types, so G is bipartite; we encode this with d×d the mask M ∈ {0, 1} , where [M]i,j = 0 if i and j are of the same type and 1 otherwise. For a candidate graph Gest with weighted adjacency West and parent sets PGest (i), we score the fit by the mean negative log-likelihood (m-NLL) e West , {γ i }i∈Vord ; X n d 1 XX ≜− ln p Xj,i | Xj,PGest (i) ; West , γ i , (5) n j=1 i=1 where γ i are the cutpoints of ordinal node i; (5) decomposes over nodes. We use two strategies depending on d. A. Exhaustive search for small graphs For d ≤ 5 we enumerate the feasible set U ′ (d) of d×d all acyclic binary matrices West,B ∈ {0, 1} satisfying West,B ⊙(1 − M) = 0, which always contains the true graph, and solve ∗ Gest = arg min′ e West , {γ i }i∈Vord ; X . (6) Gest ∈U (d)
For each candidate skeleton, the edge weights of exponential family nodes and the joint weights and cutpoints of ordinal nodes are obtained by per-node conditional maximum likelihood via L-BFGS-B, with the cutpoints reparametrized through softplus increments to enforce strict ordering. Exhaustive search returns the exact minimizer of (6), so any recovery error is attributable to finite-sample noise alone. B. Masked continuous optimization for large graphs Since |U ′ (d)| grows super-exponentially in d, for d ≥ 6 we embed structure learning into a continuous optimization over weighted adjacency matrices using DAGMA [17] as the backbone. The bipartite constraint is enforced throughout by optimizing over the masked parameters W ⊙ M, since the ℓ1 penalty alone does not guarantee a bipartite solution: ∗ West = arg
min
W∈Rd×d , {γ i }
subject to
e (W ⊙ M, {γ i } ; X) + λ||W ⊙ M||1
h (W ⊙ M) = 0,
(7)
where λ > 0 controls sparsity and h (W) = − log det (τ Id − W ◦ W) + d log τ , τ > 0, vanishes if and only if W is acyclic. The cutpoints are optimized jointly with W, since the skeleton is itself a product
of the optimization. Denoting the penalized objective by Q (W, {γ i } ; X) ≜ e (W ⊙ M, {γ i } ; X) + λ||W ⊙ M||1 , we solve (7) by the central path method of [17], with the mask applied before every loss, gradient, and acyclicity evaluation, as summarized in Algorithm 1. V. N UMERICAL R ESULTS AND D ISCUSSION We validate the identifiability result using the algorithms of Section IV: exhaustive search on the d = 3 case, the smallest setting containing distinct Markov equivalence classes, and masked DAGMA on d = 50 bipartite DAGs. Classically, graphs within a single MEC are not identifiable; Theorem 1 predicts that the SSM breaks this symmetry for ordinalexponential family edges. We measure recovery by the structural Hamming distance (SHD), which counts all errors in the presence/absence and direction of edges, normalized by the number of possible directed edges (nSHD). We consider four ground truth DAGs: G1 : X1 → X2 → X3 ,
G2 : X1 → X2 ← X3 ,
G3 : X1 ← X2 → X3 ,
G4 : X1 ← X2 ← X3 .
G1 , G3 , and G4 share a MEC and are classically unidentifiable from observations alone, whereas Theorem 1 predicts that all four are distinguishable under the SSM. Here X2 is the exponential family node, drawn from each of the seven distributions in Table I in turn, and X1 , X3 are ordinal with s = 4 categories. Figure 1 shows nSHD vs the number of samples for these four graphs, with one curve per distribution of X2 , all edge weights in {0, 1}, and N ranging from 1 to 1000. For all 4 DAGs and all seven distributions, the nSHD approaches 0 as the number of samples increases, empirically demonstrating the identifiability of edge directions between ordinal and exponential-family nodes. Figure 2 presents the results for the 50 node bipartite DAG case, which uses the masked DAGMA algorithm, where ER-k denotes Erdős–Rényi graphs with kd expected edges. We consider eight distribution modes: seven homogeneous modes in which all exponential family nodes share one distribution from Table I, and a mixed mode in which each exponential family node is assigned one of the seven uniformly at random. Across both ER-2 and ER4 settings, the nSHD decreases sharply with increasing N , approaching zero by moderate sample sizes for most modes, confirming that the masked DAGMA estimator recovers both the sparsity pattern and the edge orientations as the data grows. VI. C ONCLUSIONS We establish identifiability for causal discovery in DAGs whose nodes follow either ordinal or regular one-parameter exponential family distributions: the direction of the edge between an ordinal node and an exponential family node in a bivariate SSM is uniquely determined from the joint distribution for generic parameter values, subsuming the OrdinalPoisson result of [13]. We develop an exhaustive search for small DAGs and a masked DAGMA optimization for larger DAGs, and our numerical results confirm the decay of nSHD
to zero with sample size for three-node and 50-node bipartite DAGs, signaling recovery of the true DAG even within the MEC. R EFERENCES [1] J. Pearl, Causality, 2nd ed. Cambridge University Press, 2009. [2] J. Peters, D. Janzing, and B. Schölkopf, Elements of causal inference: foundations and learning algorithms. The MIT press, 2017. [3] P. Spirtes, C. Glymour, and R. Scheines, Causation, prediction, and search. The MIT press, 2001. [4] J. Peters and P. Bühlmann, “Identifiability of gaussian structural equation models with equal error variances,” Biometrika, vol. 101, no. 1, pp. 219– 228, 03 2014. [5] S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen, “A linear non-gaussian acyclic model for causal discovery,” Journal of Machine Learning Research, vol. 7, no. 72, pp. 2003–2030, 2006. [6] J. Peters, J. M. Mooij, D. Janzing, and B. Schölkopf, “Causal discovery with continuous additive noise models,” Journal of Machine Learning Research, vol. 15, no. 58, pp. 2009–2053, 2014. [7] K. J. Olival, P. R. Hosseini, C. Zambrana-Torrelio, N. Ross, T. L. Bogich, and P. Daszak, “Host and viral traits predict zoonotic spillover from mammals,” Nature, vol. 546, pp. 646 – 650, 2017. [8] C. K. Johnson, P. L. Hitchens, P. S. Pandit, J. Rushmore, T. S. Evans, C. C. W. Young, and M. M. Doyle, “Global shifts in mammalian population trends reveal key predictors of virus spillover risk,” Proceedings of the Royal Society B: Biological Sciences, vol. 287, no. 1924, p. 20192736, 04 2020. [9] C. K. Johnson, P. L. Hitchens, T. S. Evans, T. Goldstein, K. Thomas, A. Clements, D. O. Joly, N. D. Wolfe, P. Daszak, W. B. Karesh, and J. A. K. Mazet, “Spillover and pandemic properties of zoonotic viruses with high host plasticity,” Scientific Reports, vol. 5, 2015. [10] G. Park and G. Raskutti, “Learning large-scale poisson DAG models based on overdispersion scoring,” in Advances in Neural Information Processing Systems, vol. 28, 2015. [11] G. Park and H. Park, “Identifiability of generalized hypergeometric distribution (GHD) directed acyclic graphical models,” in Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, 16–18 Apr 2019, pp. 158–166. [12] Y. Ni and B. Mallick, “Ordinal causal discovery,” in Proceedings of the 38th Conference on Uncertainty in Artificial Intelligence, 01–05 Aug 2022, pp. 1530–1540. [13] J. Shaska, Y. Wang, C. K. Johnson, and U. Mitra, “Ordinal-Poisson causal discovery,” in Proceedings of the 61st Annual Allerton Conference on Communication, Control, and Computing, 2025. [14] D. M. Chickering, “Optimal structure identification with greedy search,” Journal of Machine Learning Research, vol. 3, no. Nov, pp. 507–554, 2002. [15] X. Zheng, B. Aragam, P. K. Ravikumar, and E. Xing, “DAGs with NO TEARS: Continuous optimization for structure learning,” in Advances in Neural Information Processing Systems, vol. 31, 2018. [16] I. Ng, A. Ghassami, and K. Zhang, “On the role of sparsity and dag constraints for learning linear dags,” in Advances in Neural Information Processing Systems, vol. 33, 2020, pp. 17 943–17 954. [17] K. Bello, B. Aragam, and P. Ravikumar, “DAGMA: Learning DAGs via M-matrices and a log-determinant acyclicity characterization,” in Advances in Neural Information Processing Systems, vol. 35, 2022, pp. 8226–8239. [18] S. Mishra and U. Mitra, “Causal discovery in equal variance linear gaussian DAGs via SURE-tuned ridge regression,” arXiv preprint arXiv:2608.17132, 2026. [19] P. McCullagh, “Regression models for ordinal data,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 42, no. 2, pp. 109–127, 01 1980.
7 U X W K G3 X1 → X2 → X3
X2 a 3 R L V V R Q X2 a ( [ S R Q H Q W L D O X2 a * D X V V L D Q
X2 a * D P P D α X2 a % L Q R P L D O
7 U X W K G4 X1 → X2 ← X3
X2 a 3 D V F D O X2 a * D P P D β
X2 a 3 R L V V R Q X2 a ( [ S R Q H Q W L D O X2 a * D X V V L D Q
X2 a * D P P D α X2 a % L Q R P L D O
X2 a 3 D V F D O X2 a * D P P D β
D Y H U D J H 6 + '
D Y H U D J H 6 + '
Q X P E H U R I V D P S O H V N
(a) G1 : X1 → X2 → X3 7 U X W K G5 X1 ← X2 → X3
X2 a 3 R L V V R Q X2 a ( [ S R Q H Q W L D O X2 a * D X V V L D Q
Q X P E H U R I V D P S O H V N
(b) G2 : X1 → X2 ← X3
X2 a * D P P D α X2 a % L Q R P L D O
7 U X W K G6 X1 ← X2 ← X3
X2 a 3 D V F D O X2 a * D P P D β
X2 a 3 R L V V R Q X2 a ( [ S R Q H Q W L D O X2 a * D X V V L D Q
X2 a * D P P D α X2 a % L Q R P L D O
X2 a 3 D V F D O X2 a * D P P D β
D Y H U D J H 6 + '
D Y H U D J H 6 + '
Q X P E H U R I V D P S O H V N
(c) G3 : X1 ← X2 → X3
$ O O 3 R L V V R Q $ O O ( [ S R Q H Q W L D O $ O O * D X V V L D Q $ O O * D P P D α
$ O O % L Q R P L D O $ O O 3 D V F D O $ O O * D P P D β ) X O O 0 L [ D O O
Q X P E H U R I V D P S O H V N
(d) G4 : X1 ← X2 ← X3
Fig. 1. Normalized SHD versus number of samples for the four three-node graphs with edge weights in {0, 1}.
( 5 % L S D U W L W H ' $ * $ O O 3 R L V V R Q $ O O ( [ S R Q H Q W L D O $ O O * D X V V L D Q $ O O * D P P D α
( 5 % L S D U W L W H ' $ *
$ O O % L Q R P L D O $ O O 3 D V F D O $ O O * D P P D β ) X O O 0 L [ D O O
Q R U P D O L ] H G 6 + '
Q R U P D O L ] H G 6 + '
Q X P E H U R I V D P S O H V N
(a) ER-2
Q X P E H U R I V D P S O H V N
(b) ER-4
Fig. 2. Normalized SHD versus number of samples for d = 50 node bipartite DAGs recovered by masked DAGMA.