Causal Discovery on Dependent Mixed Data with Applications to Gene Regulatory Network Inference

Abstract: Causal discovery aims to infer causal relationships among variables from observational data, typically represented by a directed acyclic graph (DAG). Most existing methods assume independent and identically distributed observations, an assumption often violated in modern applications. In addition, many datasets contain a mixture of continuous and discrete variables, which further complicates causal modeling when dependence across samples is present. To address these challenges, we propose a de-correlation framework for causal discovery from dependent mixed data. Our approach formulates a structural equation model with latent variables that accommodates both continuous and discrete variables while allowing correlated Gaussian errors across units. We estimate the dependence structure among samples via a pairwise maximum likelihood estimator for the covariance matrix and develop an EM algorithm to impute latent variables underlying discrete observations. A de-correlation transformation of the recovered latent data enables the use of standard causal discovery algorithms to learn the underlying causal graph. Simulation studies demonstrate that the proposed method substantially improves causal graph recovery compared with applying standard methods directly to the original dependent data. We apply our approach to single-cell RNA sequencing data to infer gene regulatory networks governing embryonic stem cell differentiation. The inferred regulatory networks show significantly improved predictive likelihood on test data, and many high-confidence edges are supported by known regulatory interactions reported in the literature.

• The paper extends de-correlation-based causal discovery to mixed datasets with continuous and multi-level discrete variables by modeling discrete observations through latent Gaussian thresholds and correlated errors.
• The proposed consensus and averaging strategies improve CPDAG recovery across MMHC, PC, and Copula-PC, with reported F1 gains of up to 38% in high-dimensional, small-sample simulations.
• Applied to 859 single-cell samples and 51 marker genes, the method improved median normalized test log-likelihood from −1.5 to −0.55 and recovered biologically supported regulatory edges, including POU5F1–ZFP42 and LEFTY1–NODAL.

Causal discovery methods overwhelmingly assume that observational units are independent and identically distributed, yet this assumption is routinely violated in applications ranging from social network data to single-cell genomics. The paper under review extends a de-correlation framework for causal structure learning—previously developed for continuous data and binary data—to mixed datasets containing both continuous and multi-level discrete variables with correlated samples. The authors instantiate the framework on single-cell RNA sequencing (scRNA-seq) data to infer gene regulatory networks (GRNs), demonstrating both improved statistical fit and biologically supported edge recovery.

Model for dependent mixed data

The core modeling device is a latent-variable structural equation model (SEM). Each observed variable $X_j$ is either continuous or discrete with $C_j$ levels. A continuous latent variable $z_{ij}$ is generated by a linear SEM over a DAG $\mathcal{G}$:

$$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$

where discrete observations arise through a quantization mapping $x_{ij} = Q(z_{ij}; T_j)$ defined by ordered thresholds $\tau_{j,c}$, while continuous variables are observed directly as $x_{ij} = z_{ij}$. Fixing $\text{diag}(\Sigma)=1$ ensures identifiability of $\{\beta_j, T_j, \Sigma\}$. Between-unit dependence enters exclusively through the exogenous noise: correlated errors across rows of the data matrix induce dependence among units without altering the causal structure encoded across variables.

The authors argue this formulation is preferable to a multi-logit alternative for discrete variables on two grounds: parameter economy—a $C_j$0-category variable requires only $C_j$1 thresholds rather than $C_j$2 regression coefficients—and preservation of ordinal structure among levels, which matters for discretized expression states such as low/medium/high.

A key observation motivates the entire approach: the SEM can be rewritten entirely over the latent variables $C_j$3, so the same DAG governs the latent continuous process. If the sample-level dependence can be removed from $C_j$4, any standard i.i.d.-based DAG learning algorithm becomes applicable.

Pre-estimation of thresholds and covariance

Estimation proceeds in two phases. In the pre-estimation phase, an initial parent set is obtained via a mixed-data version of MMHC applied to the raw data. Thresholds $C_j$5 and coefficients $C_j$6 are then estimated jointly by blockwise coordinate ascent on the observed-data log-likelihood, initialized at $C_j$7 so that the likelihood factorizes across samples. The $C_j$8-step uses an EM algorithm with truncated-normal conditional expectations; the threshold step re-parameterizes successive differences to enforce ordering and is solved with L-BFGS-B. Convergence is reported to be fast.

Given $C_j$9 and $z_{ij}$0, the covariance matrix $z_{ij}$1 is estimated pairwise: for each unit pair $z_{ij}$2, a composite likelihood combines bivariate Gaussian CDF probabilities for discrete variable pairs with bivariate Gaussian densities for continuous pairs, and the scalar correlation $z_{ij}$3 is optimized univariately. A block-diagonal structure for $z_{ij}$4 is assumed, with blocks corresponding to clusters of dependent units; if the assembled matrix fails to be positive definite, negative eigenvalues are truncated and off-diagonals rescaled.

Two caveats bear directly on downstream performance. First, the block structure must be supplied or inferred—in simulations it is known, and in the scRNA-seq application it is derived from hierarchical clustering on background genes—so misspecification of blocks is an implicit assumption throughout. Second, covariance estimation accuracy depends on having many variables: the authors note that RMSE improves as $z_{ij}$5 grows because the pairwise likelihood effectively treats variables as replicates, which explains why gains are larger in high-dimensional settings.

De-correlation and structure learning

In the second phase, $z_{ij}$6 and $z_{ij}$7 are fixed. Continuous columns are de-correlated directly via $z_{ij}$8, where $z_{ij}$9 is the Cholesky factor of $\mathcal{G}$0. For discrete columns, an EM algorithm iteratively reconstructs latent data: the E-step approximates $\mathcal{G}$1 by Monte Carlo averaging over Gibbs samples from the $\mathcal{G}$2-dimensional truncated Gaussian (exploiting block structure to keep sampling tractable), and the M-step regresses the Cholesky-transformed latents on transformed parents with ridge regularization. Because Monte Carlo noise propagates into structure estimates, two aggregation strategies are proposed: Consensus, which runs a standard discovery algorithm on $\mathcal{G}$3 de-correlated datasets and retains edges appearing in at least half of the resulting CPDAGs, and Average, which applies the algorithm once to the element-wise average of the de-correlated datasets. Both are benchmarked against applying the same algorithm to the raw mixed data.

Simulation results

Simulations use random DAGs and eight bnlearn benchmark networks (Hailfinder through Link, spanning $\mathcal{G}$4 to $\mathcal{G}$5), with half the variables discretized at thresholds $\mathcal{G}$6 and block-diagonal covariance structures of Equal or Toeplitz type. Performance is measured by F1-score on CPDAG recovery, counting undirected edges in both directions so that misorientation is penalized.

Three findings stand out. First, de-correlation substantially removes dependence: average within-block correlation among continuous variables drops from 0.62 in the original data to 0.09 after transformation. Second, covariance estimation is robust to imperfect parent sets—with MMHC-estimated parents yielding F1-scores of 0.75 at $\mathcal{G}$7 but only 0.41 at $\mathcal{G}$8, the covariance RMSE is nonetheless similar whether true or estimated parents are used. Third, the consensus and average approaches outperform the baseline across MMHC, PC, and Copula-PC in nearly all settings, with the largest gains in high-dimensional, small-sample regimes ($\mathcal{G}$9): F1 improvements of 38% (Pigs), 24% (Link), 20% (Munin1), and 16% (Andes) are reported for specific method–network combinations. Copula-PC benefits most, consistent with its sensitivity to correlation-matrix sampling error at large $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$0. The consensus approach is best or comparable almost everywhere, while the average approach offers a cheaper alternative with comparable recovery.

An honest limitation visible here is that gains shrink when $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$1 is adequate relative to graph complexity—for example, MMHC on Hepar2 and Hailfinder at $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$2 shows near-parity between baseline and de-correlated approaches—so the framework's value concentrates precisely where dependence and dimensionality stress standard methods.

Application to GRN inference

The empirical application uses the Chu et al. ESC differentiation dataset (GSE75748): 859 cells across six developmental cell types (HFF controls excluded) and $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$3 lineage-informative marker genes. Genes are classified as continuous or discretized via BIC comparison of single-Gaussian versus Gaussian mixture fits combined with Shapiro-Wilk normality tests per component, yielding 31 continuous and 20 discrete variables. Cell-type block structure is inferred by hierarchical clustering on 2,000 background genes, cut into 100 clusters (90% of size ≤ 30).

Model evaluation uses 10-fold cross-validation with folds built from independent blocks, comparing three variants: Baseline (mixed-data MMHC on raw data), Consensus-Ident (the full pipeline with $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$4, isolating discretization handling alone), and Consensus (full de-correlation). Test log-likelihoods are computed with a covariance matrix estimated from 100 held-out background genes to avoid leakage. The result is stark: the consensus method achieves a median normalized test log-likelihood of −0.55 versus −1.5 for both alternatives, with non-overlapping distributions. Because Consensus-Ident shares everything except the de-correlation step, this gap attributes the improvement specifically to removing between-cell dependence rather than to the latent-variable treatment of discreteness.

For the final network, bootstrap resampling (50% of cells, 429 cells, with re-estimated block structure per replicate) yields confidence scores per edge; directed edges with confidence ≥ 0.5 are retained, and undirected edges are oriented when one direction dominates by a factor of 3 or more. The resulting graph contains 39 edges (23 directed, 16 undirected). Many high-confidence edges match known biology: POU5F1–ZFP42 (100%), NANOG–POU5F1 (99%), LEFTY1–NODAL (98.5%, a documented feedback loop), PRDM14 → NANOG (61.2%), and DNMT3B–POU5F1 (84.5%). Edges lacking literature support may reflect unknown interactions or indirect regulation routed through background genes not included in the graph—an interpretive caveat the authors acknowledge.

Limitations and open questions

The paper is explicit about several constraints. The framework targets purely observational data and cannot incorporate interventional data or prior causal constraints, both of which would aid identifiability. The SEM is linear in the latent variables; the authors sketch how a nonlinear extension could work—if $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$5, then Cholesky-decorrelated noise yields $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$6 with independent arguments—but leave the construction and estimation of such nonlinear de-correlated models open. Additionally, the block-diagonal assumption on $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$7 is imposed rather than tested, sparsity within blocks is not exploited, and the requirement that blocks be identifiable from auxiliary features (here, background genes) limits applicability when no such side information exists. Whether the consensus vote's majority rule remains reliable when individual CPDAG estimates are highly unstable at very large $$z_{ij} = \sum_{k \in PA_j} \beta_{kj} x_{ik} + \varepsilon_{ij}, \qquad \varepsilon_j \sim \mathcal{N}_n(0,\Sigma),$$8 is also not examined.

Conclusion

This work generalizes de-correlation-based causal discovery from continuous and binary data to general mixed data with dependent units, combining threshold-based latent Gaussian modeling, pairwise maximum likelihood covariance estimation, and EM-based latent imputation with Cholesky transformation. The central methodological claim—that dependence among units can be corrected upstream of any standard DAG learner, rather than requiring redesigned algorithms—is supported by consistent simulation gains across MMHC, PC, and Copula-PC, particularly in high-dimensional settings, and by a cross-validated likelihood improvement on scRNA-seq data attributable specifically to de-correlation. The bootstrap confidence scheme provides a practical mechanism for extracting literature-consistent regulatory interactions from equivalence-class output. Remaining questions concern nonlinear SEMs, integration of interventional knowledge, and robustness to block-structure misspecification.