Node-Level Mechanism-Shift Detection in Linear Non-Gaussian Causal Networks: A Bootstrap-Based Framework with Application to the Human Gut Microbiome

1. Introduction

The microbiome consists of a characteristic community of microorganisms, including their genes and environmental context, that occupies a given habitat  [1]. The microbiome plays a crucial role in human health. In particular, the gut microbiome has been linked to metabolic, immune, and neurological disease states [2]. A common research question in this field is to define how the abundance of specific taxa varies between distinct cohorts, such as healthy versus diseased or pre- and post-intervention [3]. The primary objective is to determine whether any detected difference reflects a genuine change in a taxon’s ecological role.

Most methods developed for this task are associational, relying on correlation-based or co-occurrence networks that are constructed separately for each group and then compared, either edge-by-edge or via global topological summaries. NetCoMi [4] formalized this workflow in a general-purpose R package, providing compositionally aware association measures, such as correlation estimators (e.g., SparCC [5]), partial correlation estimators [6], and proportionality [7], together with sparsification and formal statistical tests for whether individual edges, or the network as a whole, differ between two groups. Liu et al. [8] distilled this into a broader methodological guide for comparing co-occurrence networks constructed from different sample groups or with different construction pipelines, and Wen et al. [9] extended multi-network comparison to multi-factor, multi-treatment, and transkingdom settings. However, despite the maturity of these methods, the principle correlation is not causation remains [10]. Ultimately, the relationships identified through association-network comparisons indicate association, rather than true causal dependencies.

This reflects a growing need within the field to transition from purely associational to causal approaches [11,12]. Recent studies have begun to address this challenge directly by adapting established causal discovery algorithms to the specific characteristics of microbiome data. Sazal et al. [13] combined Bayesian networks with co-occurrence networks into a signed Bayesian network to infer directional relationships and colonization order within microbial communities, and Corander et al. [14] extended standard causal graphs to labeled DAGs that encode context-specific independence, applying the framework to identify the causal role of Clostridioides difficile in driving diarrhea risk. Transitioning from structural representation to quantitative inference, Sazal et al. [15] later used interventional calculus to estimate causal effects between taxa directly from observational abundance data. As an alternative, functional causal models leverage the statistical properties of observational data to infer directionality of relationships. Among these, DirectLiNGAM [16] is particularly well suited to taxon-level abundance data, since it guarantees identification of the causal ordering in a finite number of steps under linearity, non-Gaussian noise, and acyclicity-assumptions with reasonable justification once compositional abundances are log-ratio transformed (Section 3.6.1). MultiGroupDirectLiNGAM [17] performs a joint estimation across several datasets under the assumption that they have a common causal order, leveraging shared information across datasets to improve estimation accuracy compared to independent fitting. However, in complex biological systems, maintaining a strictly conserved causal order is rarely a realistic assumption. Consequently, imposing a shared structure can hide the exact ecological and structural shifts that a comparative analysis aims to detect. While standard DirectLiNGAM has been employed to estimate causal networks within isolated gut-microbiome cohorts [18,19,20,21], these applications do not address whether the underlying causal architecture actually changes between different groups. To bridge this gap and evaluate structural differences without forcing a shared causal ordering, one must compare independently estimated networks. Toward this end, Ma and Tourani [22] defined a resampling-based estimator for quantifying, edge by edge, how confidently an individual causal relationship differs between two learned DAGs. However, this estimator operates strictly at the edge level. It tests for isolated changes rather than assessing whether a taxon’s full regulatory mechanism (its complete set of causal parents) has shifted. In microbiome data, evaluating the entire regulatory neighborhood is often more biologically meaningful, as a taxon’s ecological role is shaped by its broader network context rather than any single interaction [23,24].

In this work, we address this gap by extending the edge-level network-comparison framework of Ma and Tourani [22] to enable node-level mechanism-shift detection under DirectLiNGAM. DirectLiNGAM is well suited to this setting because it is a computationally fast algorithm capable of operating in high-dimensional regimes and relatively small sample sizes.

The prevalence of these small datasets, coupled with frequent sample size imbalances between cohorts, motivates our choice to evaluate three distinct estimators. First, we formalize a mechanism-perturbation ground truth where an entire node’s parent set, rather than a single edge, is rewired between two systems (Section 2.2.2). We then evaluate the Naive, Bootstrap, and RSBS (Relative sample size Bootstrap Stability) estimators. The RSBS estimator is specifically designed to remain reliable when the compared datasets differ substantially in size. This addresses a critical limitation of previous associational methods, as severe sample size imbalances are the rule rather than the exception when pooling independent microbiome cohorts. We validate this framework across an extensive set of synthetic experiments spanning both equal and heterogeneous sample sizes. To ensure these findings translate to data with authentic compositional structure and non-idealized noise, we also perform a semi-synthetic validation using a batch-corrected human gut-microbiome cohort (Section 3.6). Finally, we provide the full benchmark pipeline as an open-source Python package, including the three estimators and the underlying generators for synthetic graphs and mechanism perturbations, allowing other researchers to construct their own benchmarks with known ground truths.

This study reveals that resampling-based estimation does not uniformly outperform a single-fit Naive baseline. Specifically, standard Bootstrap estimation is the superior choice when the compared datasets are of equal size. Conversely, RSBS is the only estimator that reliably outperforms the Naive approach under sample size imbalance, an advantage that becomes more pronounced as network size grows. This pattern holds in the batch-corrected gut microbiome cohort, proving that the advantages of this approach generalize beyond idealized synthetic graphs, handling real compositional and technical noise.

The remainder of this paper is organized as follows. Section 2 introduces the LiNGAM framework, DirectLiNGAM, and our proposed node-level mechanism-shift extension. Section 3 details the experimental results on fully synthetic graphs under equal and heterogeneous sample sizes, as well as the semi-synthetic microbiome validation. Finally, Section 4 and Section 5 discuss the broader implications of this work and summarize the key contributions.

2. Materials and Methods

2.1. Background: The LiNGAM Framework & DirectLiNGAM

Our approach adopts the linear, non-Gaussian acyclic model (LiNGAM) framework, originally introduced by Shimizu et al. [25], for causal discovery from purely observational data. Rather than a single algorithm, LiNGAM is a general framework that includes a variety of specific causal discovery approaches.

Among these, we focus on DirectLiNGAM [16], which relies on four core assumptions about the data-generating process: the relationships are linear, the causal structure forms a directed acyclic graph (DAG), the noise terms are non-Gaussian, and there are no unobserved confounders.

Formally, the process can be represented as a directed acyclic graph (DAG) over p variables $\{x_1,\dots ,x_p\}$, encoded by a $p\times p$ adjacency matrix $\mathbf{B}=\left\{b_{ij}\right\}$, where $b_{ij}$ quantifies the direct causal effect of $x_j$ on $x_i$. Let $k\left(i\right)$ denote the position of $x_i$ in an unknown causal ordering, such that no variable appearing later in the ordering has a directed path to any earlier variable; under this ordering, an edge $x_j\to x_i$ is permitted only if $k\left(j\right)<k\left(i\right)$. Every variable is assumed to have zero mean, the data-generating process is linear, and it can be written as follows:

$$x_i=\sum _{k\left(j\right)<k\left(i\right)}{b}_{ij}{x}_j+e_i,$$

(1)

where $e_i$ are continuous random variables with non-Gaussian distributions, zero mean, and nonzero variance, and are mutually independent, so that there are no latent confounding variables. Equation (1) can be written in compact form as in Equation (2):

$$\mathbf{x}=\mathbf{B}\mathbf{x}+\mathbf{e},$$

(2)

where $\mathbf{x}$ is a p-dimensional random vector and $\mathbf{e}$ is the p-dimensional error vector; $\mathbf{B}$ can always be brought to strictly lower-triangular form by simultaneous, equal row and column permutation, as a consequence of the acyclicity assumption [26].

As an illustration, consider a three-variable causal chain $x_1\to x_2\to x_3$, with structural equations

$$x_1=e_1,x_2=b_{21}{x}_1+e_2,x_3=b_{32}{x}_2+e_3,$$

(3)

so that $x_1$ is an exogenous variable (a root node, determined only by its noise term), while $x_2$ and $x_3$ are endogenous. The corresponding adjacency matrix,

$$\mathbf{B}=\left(\begin{array}{ccc}0& 0& 0\\ b_{21}& 0& 0\\ 0& b_{32}& 0\end{array}\right),$$

is strictly lower triangular under the causal ordering $k=(1,2,3)$. Given observations only, the estimation problem is to recover $\mathbf{B}$ (equivalently, the ordering $k(·)$ together with the regression weights). Early approaches, such as ICA-LiNGAM [25], estimate the mixing matrix$\mathbf{A}={(I-\mathbf{B})}^{-1}$ via independent component analysis, but suffer from the absence of guaranteed convergence to the correct solution from a poor initialization and from the need for non-trivial parameter selection. DirectLiNGAM instead estimates the causal ordering directly, by sequentially identifying exogenous variables using an independence measure. Once the ordering is estimated, $\mathbf{B}$ is constructed by regressing each variable on its predecessors, through an adaptive LASSO.

2.1.1. DirectLiNGAM as a Causal Discovery Algorithm Reference

Throughout this paper we write M for a generic causal discovery algorithm: a procedure that takes a dataset D generated from the linear structural equation model (SEM) of Equation (2) and returns an estimated (weighted) adjacency matrix $\widehat{\mathbf{B}}$. M is treated as a black box by the Naive, Bootstrap, and RSBS estimators introduced in Section 2.2.3, which act only on ${\widehat{\mathbf{B}}}^{\left(s\right)}$ (Equations (6)–(10)) and never on the internal mechanics of how M produced it. We adopt DirectLiNGAM as the concrete instantiation of M used throughout this study.

2.2. Formal Framework

2.2.1. Problem Formulation and Notation

We extend the network-comparison framework of Ma and Tourani [22], applied to the LiNGAM structural equation model of Section 2.1, from edge-level difference estimation to node-level mechanism-shift detection between a pair of systems.

Let $V=\{x_1,\dots ,x_p\}$. A system $s\in \{1,2\}$ is defined by its own weighted, strictly lower-triangular matrix ${\mathbf{B}}^{\left(s\right)}\in {\mathbb{R}}^{p\times p}$, generating data via the LiNGAM structural equation model of Equation (2), ${\mathbf{x}}^{\left(s\right)}={\mathbf{B}}^{\left(s\right)}{\mathbf{x}}^{\left(s\right)}+{\mathbf{e}}^{\left(s\right)}$, i.e., ${\mathbf{x}}^{\left(s\right)}={\left(I-{\mathbf{B}}^{\left(s\right)}\right)}^{-1}{\mathbf{e}}^{\left(s\right)}$. For any node k, its set of causal parents in system s is denoted by ${\mathrm{Pa}}^{\left(s\right)}\left(k\right)=\{j:B_{kj}^{\left(s\right)}\ne 0\}$. Both systems share the same true causal order $\pi$ over $\{1,\dots ,p\}$: ${\mathbf{B}}^{\left(2\right)}$ is always constructed from ${\mathbf{B}}^{\left(1\right)}$ by rewiring entries within existing rows (Section 2.2.2), never by altering which index precedes which. Graphs, edge weights, and noise terms are generated following the standard DirectLiNGAM simulation protocol of Shimizu et al. [16], with the sole difference that we generate a pair of related systems, ${\mathbf{B}}^{\left(1\right)}$ and ${\mathbf{B}}^{\left(2\right)}$, rather than a single one.

Baseline graphs are generated with each candidate edge $(k,j)$, $j<k$, included independently with probability $p_e$ (edge density). Prior to estimation, node labels are blinded independently for each system: a random permutation ${\rho }^{\left(s\right)}$ of $\{1,\dots ,p\}$ is applied to the columns of $D^{\left(s\right)}$, so the causal discovery algorithm M observes variables in an order unrelated to $\pi$, exactly as M must recover in practice. Let ${\widehat{\mathbf{B}}}^{\left(s\right)}$ denote the weighted adjacency matrix returned by M in this permuted index space. Before any comparison, each ${\widehat{\mathbf{B}}}^{\left(s\right)}$ is re-aligned to the shared ground-truth order by inverting its blinding permutation, mapping row and column k back to index ${\left({\rho }^{\left(s\right)}\right)}^{-1}\left(k\right)$; all downstream quantities are defined on these realigned matrices.

2.2.2. Mechanism-Perturbation Ground Truth

Rather than perturbing individual edges independently anywhere in the graph, as in the general formulation of $E_i-E_j$ in [22] (Eq. 1), our extension constrains perturbations to act on an entire node’s mechanism (its full parent set), which is the natural unit of comparison once edge-level differences are aggregated row-wise (Section 2.2.3).

Given baseline ${\mathbf{B}}^{\left(1\right)}$, a target set $T\subset \{1,\dots ,p\}$ with $\left|T\right|=n_{\mathrm{pert}}$ is drawn from the subset of nodes k satisfying $|{\mathrm{Pa}}^{\left(1\right)}\left(k\right)|\ge n_{\mathrm{pos}}$ (enough existing parents to remove) and $|\{1,\dots ,k-1\}\setminus {\mathrm{Pa}}^{\left(1\right)}\left(k\right)|\ge n_{\mathrm{pos}}$ (enough eligible non-parent predecessors to add). For each $k\in T$, ${\mathbf{B}}^{\left(2\right)}$ is obtained from ${\mathbf{B}}^{\left(1\right)}$ by removing $n_{\mathrm{pos}}$ uniformly sampled entries from row k and adding $n_{\mathrm{pos}}$ new nonzero entries at uniformly sampled positions in row k not already in ${\mathrm{Pa}}^{\left(1\right)}\left(k\right)$, with weights redrawn under the same protocol as Section 2.2.1. For $k\notin T$, row k of ${\mathbf{B}}^{\left(2\right)}$ is identical to row k of ${\mathbf{B}}^{\left(1\right)}$, so the structural equation for every non-target node is exactly preserved across systems. The ground truth is the binary vector

$$\mathbf{y}\in {\{0,1\}}^p,y_k=\mathbf{1}[k\in T],$$

(4)

which is what the node-level estimators in Section 2.2.3 are evaluated against.

2.2.3. Estimators

For a chosen discovery algorithm M (DirectLiNGAM throughout, Section 2.1.1), let $D^{\left(1\right)},D^{\left(2\right)}$ be datasets of size $N_1,N_2$ simulated from systems 1 and 2 via Equation (2). We extend the three estimators of [22] from the constraint-based, unweighted setting to the weighted linear-SEM setting as follows.

Naive.

A single fit of M is obtained on each dataset, and the estimated adjacency is binarized:

$$\begin{array}{cc}\hfill {\widehat{B}}_{kj}^{\left(s\right),\mathrm{naive}}& =\mathbf{1}\left[\left|{\widehat{B}}_{kj}^{\left(s\right)}\right|>0\right],\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\Delta }^{\mathrm{naive}}& ={\widehat{\mathbf{B}}}^{\left(1\right),\mathrm{naive}}-{\widehat{\mathbf{B}}}^{\left(2\right),\mathrm{naive}}\in {\{-1,0,1\}}^{p\times p},\hfill \end{array}$$

(6)

the direct linear-SEM analogue of Eq. 3 in [22].

2.2.3.2. Bootstrap.

M is refit on $n_{\mathrm{sampling}}$ bootstrap resamples (with replacement) of each dataset, yielding the edge-selection frequency

$$\begin{array}{cc}\hfill {\widehat{B}}_{kj}^{\left(s\right),\mathrm{BS}}& ={\displaystyle \frac{1}{n_{\mathrm{sampling}}}}{\displaystyle \sum _{b=1}^{n_{\mathrm{sampling}}}}\mathbf{1}\left[\mathrm{edge}(j,k)\mathrm{recovered}\mathrm{in}\mathrm{resample}b\right]\in [0,1],\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\Delta }^{\mathrm{BS}}& ={\widehat{\mathbf{B}}}^{\left(1\right),\mathrm{BS}}-{\widehat{\mathbf{B}}}^{\left(2\right),\mathrm{BS}},\hfill \end{array}$$

(8)

matching Eq. 4 in [22].

RSBS (Relative sample size Bootstrap Stability).

This is our implementation’s realization of the equal-sample size resampling estimator of [22] (Eq. 5). Let $n_{min}=min(N_1,N_2)$. The larger of $D^{\left(1\right)},D^{\left(2\right)}$ is first reduced to a single subsample of size $n_{min}$, drawn once without replacement; the smaller dataset is left unchanged (it is already of size $n_{min}$). The standard bootstrap procedure of Equation (7) (with the same $n_{\mathrm{sampling}}$) is then applied to both size-$n_{min}$ datasets:

$$\begin{array}{cc}\hfill {\widehat{\mathbf{B}}}^{\left(s\right),\mathrm{RSBS}}& ={\widehat{\mathbf{B}}}^{\left(s^{\prime }\right),\mathrm{BS}}\mathrm{evaluated}\mathrm{on}\mathrm{the}{n}_{min}-\mathrm{sized}\mathrm{version}\mathrm{of}{D}^{\left(s\right)},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\Delta }^{\mathrm{RSBS}}& ={\widehat{\mathbf{B}}}^{\left(1\right),\mathrm{RSBS}}-{\widehat{\mathbf{B}}}^{\left(2\right),\mathrm{RSBS}}.\hfill \end{array}$$

(10)

This is a practical, single-draw simplification of the general equal-sample size resampling estimator described above, which allows for the down-sampling itself to be repeated across resamples; our implementation fixes one down-sampled subset (controlled by a random seed) and bootstraps only on top of it, trading a small amount of additional variance in the down-sampling step for a substantially lower computational cost, since only one M-fit per bootstrap replicate is required instead of two nested resampling loops.

2.2.4. From Edge-Level Differences to Node-Level Mechanism Scores

The central extension of this work is the reduction ${\Delta }^{•}\in {\mathbb{R}}^{p\times p}\mapsto {\mathbf{s}}^{•}\in {\mathbb{R}}_{\ge 0}^p$ of the $p\times p$ edge-difference matrix (for $•\in \{\mathrm{naive},\mathrm{BS},\mathrm{RSBS}\}$) to a p-dimensional node score vector, aligned with the mechanism-perturbation ground truth $\mathbf{y}$ of Equation (4). Row k of ${\Delta }^{•}$ collects the estimated differences over all candidate parents of node k; summing its absolute value therefore yields a single scalar quantifying how much node k’s mechanism is estimated to have changed between the two systems:

$$s_k^{•}={\displaystyle \sum _{j=1}^p}\left|{\Delta }_{kj}^{•}\right|,k=1,\dots ,p.$$

(11)

Two properties of Equation (11) are worth making explicit. First, taking the absolute value symmetrizes the estimator: unlike the directional definition $E_i-E_j=1$ in [22] (Eq. 1), which only credits edges lost from system 1 to system 2, $s_k^{•}$ accumulates contributions from both edge deletions (${\Delta }_{kj}^{•}>0$) and edge additions (${\Delta }_{kj}^{•}<0$) at node k, matching the symmetric construction of the mechanism perturbation itself (Section 2.2.2), which always removes and adds $n_{\mathrm{pos}}$ parents. Second, because row k of ${\mathbf{B}}^{\left(s\right)}$ contains exactly ${\mathrm{Pa}}^{\left(s\right)}\left(k\right)$, restricting the sum to row k (rather than column k, or both) isolates changes to node k’s incoming edges, i.e., its structural mechanism, as opposed to changes in how node k influences its descendants; this is the node-level quantity that corresponds by construction to $y_k$ in Equation (4). The vector ${\mathbf{s}}^{•}=(s_1^{•},\dots ,s_p^{•})$ is evaluated against $\mathbf{y}$ as a binary classification problem. Since ${\Delta }_{kj}^{•}$ is confined to the strictly lower or upper-triangular structure of ${\mathbf{B}}^{\left(s\right)}$, each $s_k^{•}$ is a nonnegative integer, or real value score with no fixed upper bound other than the row’s number of eligible entries; $s^{\mathrm{naive}}$ takes integer values (a sum of $\{0,1\}$ terms) and is treated as a ranking score exactly as $s^{\mathrm{BS}},s^{\mathrm{RSBS}}$ are (a sum of $[0,1]$-valued terms), from which AUCROC and AUPR are computed directly, and $F_1$, precision, and recall are computed at the score-maximizing threshold, $\widehat{\tau }=arg{max}_{\tau \in \left\{s_k^{•}\right\}}{F}_1(\mathbf{y},\mathbf{1}[{\mathbf{s}}^{•}\ge \tau ])$, searched over the set of unique observed scores (Section 3.1) as in [22].

2.3. Software Implementation

The full benchmark pipeline used to produce the synthetic results in Section 3, from graph generation through node-level scoring, is provided as an open-source Python package named causal-comparator. The complete source code is publicly available in an open repository (see the Data Availability Statement for access details).

Concretely, causal-comparator exposes the discovery algorithm M as a model_class parameter (defaulting to lingam.DirectLiNGAM), ensuring that the architecture remains modular. Beyond the three node-level estimators themselves, the software developed for this study also implements the synthetic-graph and mechanism-perturbation generators described in Section 2.2.1 and Section 2.2.2. By providing these underlying generators alongside the estimators, this framework allows researchers to construct custom benchmark datasets with known ground truths across arbitrary network sizes, edge densities, sample size regimes, and perturbation strengths.

3. Results

3.1. Experimental Design

We evaluated the Naive, Bootstrap, and RSBS estimators of Section 2 under the mechanism-perturbation ground truth of Section 2.2.2. Concretely, for a pair of systems $(G_i,G_j)$ we generated $G_i$ as a random strictly lower-triangular weight matrix $B_i$ (edge density/probability $p_e$, coefficients drawn from $[-1.5,-0.5]\cup [0.5,1.5]$), and obtained $G_j$ by perturbing the mechanism (parent set) of a subset of $n_{\mathrm{pert}}$ target nodes: for each target, $n_{\mathrm{pos}}$ existing parents were removed and $n_{\mathrm{pos}}$ new parents were added, leaving all other nodes unchanged. Data were simulated from each system, the node labels were randomly shuffled prior to discovery to blind the causal order (as DirectLiNGAM must recover it), and the estimated adjacency/probability matrices were re-aligned to the ground-truth order before scoring. Each of the three estimators produces an edge-level difference matrix $\widehat{\Delta }\in {[-1,1]}^{p\times p}$. We converted this into a node-level score, $s_k={\sum }_l\left|{\widehat{\Delta }}_{kl}\right|$, and evaluated it against the ground-truth binary indicator of whether node k’s mechanism was perturbed, using AUCROC, AUPR, and threshold-optimized $F_1$, precision, and recall. Thresholding was performed by grid search over the unique score values, as in [22]. We ran two complementary experiments. In the equal-sample size experiment, $N_1=N_2\in \{100,200,300\}$, and only the Naive and Bootstrap estimators were evaluated, since RSBS is defined only when the two sample sizes differ (Section 2.2.3). In the heterogeneous-sample size experiment, the primary and secondary dataset sizes were $N_1\in \{400,800\}$ and $N_2\in \{100,200,400\}$ respectively (yielding sample size ratios $r=N_2/N_1\in \{0.125,0.25,0.5,1.0\}$), and all three estimators were compared. Both experiments spanned network dimensions $n_{\mathrm{nodes}}$ (corresponding to p in Section 2.2.1) $\in \{15,20,25\}$, edge density $p_e\in \{0.15,0.20,0.25\}$, $n_{\mathrm{pert}}\in \{2,3,4\}$ (the number of target nodes with altered mechanisms), and $n_{\mathrm{pos}}\in \{1,2\}$ (the number of parents swapped per target node), with 50 independently sampled seeds per configuration. For a given seed, the same pair of graphs and datasets was reused across all three estimators, allowing paired comparisons. Bootstrap and RSBS estimators utilized $n_{\mathrm{sampling}}=100$ resamples per dataset. In total, this yielded $9,900$ discovery runs for the equal sample size experiment and $29,700$ for the heterogeneous one ($39,600$ overall).

All summary tables in this section (Table 1, Table 2, Table 3, Table 4 and Table 5) aggregate over all three network sizes ($n_{\mathrm{nodes}}\in \{15,20,25\}$) unless stated otherwise. Figure 1, Figure 2 and Figure 3 and Figure 4, Figure 5 and Figure 6 show results for a representative network size of $n_{\mathrm{nodes}}=15$. We chose to plot a single size rather than overlaying all three, as doing so would obscure the trends in sample size and perturbation load that these figures aim to isolate. The corresponding 20 and 25 node panels are visually consistent with the trends shown and are provided as Figures S1–S6 (equal-sample size regime) and Figures S7–S12 (heterogeneous-sample size regime) in the Supplementary Materials. The size-dependent trends themselves are analyzed separately and quantitatively in Section 3.3 (Table 2), which does use all three network sizes.

3.2. Naive versus Bootstrap Estimation under Equal Sample Sizes

Table 1 summarizes performance across all $4,950$ runs per method in the equal-sample size experiment. Bootstrap estimation outperformed the Naive method on every ranking and classification metric except recall: mean AUPR rose from $0.687$ to $0.825$ and mean AUCROC from $0.894$ to $0.941$. A paired comparison (same seed, same graphs, same data, differing only in estimator) confirmed this is not noise: Bootstrap exceeded Naive by $+0.139$ AUPR ($t=39.7$, $p<{10}^{-16}$), $+0.047$ AUCROC ($t=27.6$, $p<{10}^{-16}$), and $+0.094$$F_1$ ($t=35.6$, $p<{10}^{-16}$). Bootstrap’s estimates are also noticeably more stable than Naive’s. The standard deviation drops from $0.296$ to $0.244$ for AUPR, $0.142$ to $0.103$ for AUCROC, $0.233$ to $0.193$ for $F_1$, and $0.291$ to $0.237$ for precision (Table 1). The one exception is recall, where Bootstrap’s SD ($0.164$) is marginally higher than Naive’s ($0.153$). Naive only exceeded Bootstrap on recall ($0.945$ vs. $0.926$; $t=-6.5$, $p<{10}^{-10}$), at the cost of a 14-point drop in precision ($0.699$ vs. $0.837$), see Section 3.5.

Figure 1 shows, for the 15-node condition, mean AUPR as a function of sample size for each edge density; the corresponding aggregate-over-all-sizes numbers are shown in Table 1– Table 2. Both methods improve monotonically with sample size, but the ordering never reverses: averaging across the three density panels shown, Bootstrap leads Naive by $0.18$ AUPR at $N_1=N_2=100$ ($0.842$ vs. $0.665$), narrowing to a $0.05$ AUPR gap at $N_1=N_2=300$ ($0.993$ vs. $0.943$). This same non-reversing ordering holds once the 20 and 25-node conditions are pooled in, at lower absolute AUPR (Table 1), because larger networks are harder to discover from a fixed, small sample (Table 2). Edge density has only a minor effect on the gap; both curves shift slightly as the graph becomes sparser, consistent with fewer high-effect-size parents to anchor discovery.

Figure 2 stress-tests the estimators (15-node condition) as the fraction of perturbed nodes increases (holding $N_1=N_2=100$). Bootstrap’s advantage over Naive is stable across densities and perturbation fractions, indicating the benefit of resampling is not an artifact of a specific graph topology or perturbation load; the same qualitative pattern, at correspondingly lower absolute AUPR, holds in the 20 and 25-node conditions (Table 2).

Figure 3 examines the effect of perturbation strength $n_{\mathrm{pos}}$ (parents swapped per target node) on precision and recall, again for the 15-node condition at the highest edge density tested ($p_e=0.25$). Increasing $n_{\mathrm{pos}}$ from 1 to 2 raises precision for both methods (Naive: $0.793\to 0.893$; Bootstrap: $0.909\to 0.965$), consistent with a larger effect size making the perturbation easier to detect, while recall remains high ($>0.96$) for both methods.

3.3. Effect of Network Size

Table 2 reports AUCROC and AUPR as the number of nodes grows from 15 to 25, with the (small) sample sizes used in our grid held fixed. Performance degrades for every method as dimensionality increases, as expected given the growing number of parameters DirectLiNGAM must estimate; but the size of the drop is not uniform across estimators. For Naive, AUPR falls from $0.835$ (15 nodes) to $0.538$ (25 nodes), a 30-point drop; for Bootstrap, it falls from $0.938$ to $0.706$, a smaller relative decline that leaves Bootstrap’s advantage over Naive larger at 25 nodes ($+0.168$ AUPR) than at 15 ($+0.104$ AUPR). The same pattern holds in the heterogeneous experiment (Section 3.4), where RSBS’s margin over Naive widens from $+0.054$ AUPR at 15 nodes to $+0.070$ at 25 nodes. Resampling-based confidence therefore becomes more valuable as the discovery problem gets harder.

3.4. Naive versus Bootstrap versus RSBS under Heterogeneous Sample Sizes

When the two datasets differ substantially in size, the ranking among estimators changes in an important way: plain Bootstrap resampling is no longer reliably better than the Naive method. Over all $9,900$ runs of the heterogeneous experiment (Table 3), Bootstrap’s mean AUCROC ($0.907$) and $F_1$ ($0.802$) were lower than Naive’s ($0.927$ and $0.824$), and the paired difference was significant (AUCROC: $-0.020$, $t\left(9899\right)=-18.2$, $p<{10}^{-16}$; $F_1$: $-0.022$, $t=-11.9$, $p<{10}^{-16}$), though the AUPR difference was not distinguishable from zero ($-0.003$, $t=-1.2$, $p=0.24$). RSBS, designed specifically for this regime, was the only estimator that dominated Naive on every metric except recall (AUPR $+0.066$, $t=28.6$; AUCROC $+0.013$, $t=12.8$; $F_1$$+0.034$, $t=19.6$; all $p<{10}^{-16}$) and also outperformed plain Bootstrap on every metric except recall (AUPR $+0.069$, AUCROC $+0.034$, $F_1$$+0.056$; all $p<{10}^{-16}$). Under the same evaluation framework, RSBS achieved or shared the top performance in $97.5\%$ of the 198 heterogeneous parameter combinations, versus $82.3\%$ for Bootstrap and $68.2\%$ for Naive.

Figure 4 shows, for the 15-node condition, mean AUPR against the larger sample size $N_1$, marginalized over $N_2$. RSBS (top curve) leads throughout; Naive (bottom curve at low density, middle at higher density) and Bootstrap cross over depending on density, reflecting the non-uniform Bootstrap-vs-Naive result above.

Figure 5 repeats the topological stress test of Figure 2 in the heterogeneous regime. RSBS sits above both Naive and Bootstrap at every density and every perturbation load tested, consistent with the overall ranking in Table 3.

Figure 6 shows precision and recall vs. perturbation strength for the 15-node condition at the highest edge density tested ($p_e=0.25$). As in the equal-sample size experiment, increasing $n_{\mathrm{pos}}$ from 1 to 2 raises precision for all three methods (Naive: $0.852\to 0.929$; Bootstrap: $0.881\to 0.965$; RSBS: $0.913\to 0.964$), and RSBS’s precision advantage over Bootstrap is concentrated at $n_{\mathrm{pos}}=1$ (the harder, single-parent-swap condition), where RSBS also has the edge over Naive.

Table 4 breaks the heterogeneous experiment down by sample size ratio $r=N_2/N_1$. RSBS’s advantage is largest exactly where it is designed to help: at the most imbalanced ratio ($r=0.125$, i.e., $N_1=800$, $N_2=100$), RSBS beats Naive by $+0.076$ AUPR and Bootstrap by $+0.079$ AUPR (paired $t=10.9$ and $11.1$, both $p<{10}^{-16}$; see Table 5), while plain Bootstrap is statistically indistinguishable from Naive on AUPR at this ratio ($-0.003$, $t=-0.5$, $p=0.61$) despite being worse on AUCROC ($-0.030$, $t=-8.5$, $p<{10}^{-16}$). As $r\to 1$ (equal sample sizes), the three methods converge, and RSBS becomes numerically indistinguishable from Bootstrap (AUCROC $0.995$ vs. $0.995$; AUPR $0.972$ vs. $0.972$ at $r=1$)—exactly as expected, since the equal-sample size resampling step in Equation 10 performs no downsampling when $N_1=N_2$ and RSBS reduces to standard bootstrap by construction. This internal consistency check gives us additional confidence in the correctness of our LiNGAM adaptation of the CausalComparator estimators.

3.5. The Precision–Recall Trade-off across Both Regimes

A consistent pattern emerges across both experiments: the Naive method achieves the highest recall of the three estimators but the lowest precision, in essentially every configuration we tested (Table 1 and Table 3). Because Naive assigns a hard $0/1$ edge label from a single DirectLiNGAM fit per dataset, it has no mechanism to discount edges that are only weakly or inconsistently supported by the data; borderline edges are as likely to flip the node-level score as strongly supported ones. This over-flags a wider set of nodes as perturbed, inflating recall while depressing precision. In the equal-sample size experiment, Naive’s recall ($0.945$) exceeds Bootstrap’s ($0.926$) by only 2 points, while its precision ($0.699$) trails Bootstrap’s ($0.837$) by 14 points. In the heterogeneous experiment, Naive’s recall ($0.966$) exceeds RSBS’s ($0.927$) by 4 points, while its precision ($0.772$) trails RSBS’s ($0.844$) by 7 points; the gap widens further under the most imbalanced sample sizes tested ($r=0.125$: Naive precision $0.582$ vs. RSBS $0.676$, a 9-point gap against a 8-point recall deficit). In every configuration we examined, the precision gain for the resampling-based estimators exceeded their recall cost, which is why they dominate on $F_1$, AUPR, and AUCROC despite never leading on recall alone. This trade-off is most pronounced exactly where discovery is hardest—at the smallest sample sizes and under the greatest sample size imbalance (Table 4, $r=0.125$)—suggesting that the confidence information carried by resampling is most useful precisely when a single point-estimate discovery run is least trustworthy.

3.6. Semi-Synthetic Benchmark Derived from Human Gut Microbiome Data

3.6.1. Cohort Construction

To assess whether the findings of Section 3.2, Section 3.3, Section 3.4 and Section 3.5 transfer beyond fully synthetic graphs, we applied the same three estimators to a real gut-microbiome cohort drawn from curatedMetagenomicData (version 3.16.1) [27], a curated, uniformly reprocessed collection of publicly available human microbiome shotgun-sequencing studies. We restricted to disease-free adult stool samples ($\mathrm{disease}=\mathrm{healthy}$, $\mathrm{body}\text{\_}\mathrm{site}=\mathrm{stool}$, $\mathrm{age}\text{\_}\mathrm{category}=\mathrm{adult}$), excluded samples with recorded current antibiotic use, pregnancy, or lactation, applied a minimum sequencing-depth floor of $1,000,000$ reads [28,29], and retained one sample per subject (earliest visit, in studies with repeated sampling), to avoid pseudo-replication from longitudinally resampled individuals. This yielded $6,871$ samples spanning multiple independent studies. Species-level relative abundances were aggregated to the family level; families present in fewer than $50\%$ of samples or averaging less than $0.01\%$ relative abundance were dropped [30,31], leaving 24 families.

Because samples originate from different studies, we corrected for study-level technical batch effects with MMUPHin [32] (adjust_batch, batch variable = study, applied on the relative-abundance scale prior to log-ratio transformation), then applied a centered log-ratio (CLR) transform [33,34] to remove the compositional (sum-to-one) constraint inherent to relative-abundance data, which otherwise induces spurious negative correlations among taxa unrelated to any real biological association. A quick diagnostic check (a one-way, Euclidean-distance analogue of PERMANOVA on the CLR-transformed family profile) confirmed that study identity explained $20.6\%$ of total compositional variance prior to correction, supporting a non-negligible batch effect worth correcting for. Because the full 24-family CLR composition sums to exactly zero for every sample by construction, its coordinates are linearly dependent, so the resulting covariance structure has rank at most $p-1$ for p parts rather than full rank [33], a well-known obstacle in compositional data analysis that also affects CLR-based microbial network inference specifically [6]. Any full-rank parameterization of a p-part composition must be $(p-1)$-dimensional; standard constructions include the isometric log-ratio basis and Aitchison’s additive log-ratio, which fixes one part as reference and is equivalent, up to a linear transformation, to dropping one CLR coordinate [33,35]. We adopted the latter in its simplest form, dropping the single highest-mean-CLR family (Lachnospiraceae) prior to modeling, leaving 23 usable, well-posed variables.

3.6.2. Semi-Synthetic Ground Truth Construction

A fundamental obstacle in validating causal discovery or causal inference methods on gut microbiome data, and on microbiome data more generally, is the near-total absence of biological ground truth. Communities with fully known interaction sets are rare, and existing validated interaction lists are incomplete rather than exhaustive. No such resource exists at the scale this cohort operates at (thousands of samples, dozens of taxa). As a result, evaluation of microbial network inference is still carried out largely in silico, using synthetic or semi-synthetic data rather than communities with fully known interactions [36]. Fully synthetic graphs (Section 2.2.1) let us control the ground truth exactly, but say nothing about whether the estimators remain informative once applied to data with real compositional structure and real non-idealized noise. We therefore constructed a semi-synthetic ground truth. It inherits this cohort’s actual compositional structure and residual noise, but includes a known, injected structural change against which detection performance can still be measured directly. For each node count $p\in \{15,20,23\}$, treated independently, we selected the p most abundant families from the pool of 23 usable families described in Section 3.6.1 and fit a separate real DirectLiNGAM model (Section 2.1.1) to the full $6,871$-sample cohort restricted to just those p families. Fitting the baseline on the full cohort, and only subsequently drawing smaller simulated samples from it for evaluation, follows the same real-data design used by Ma and Tourani [22] in their own real-world evaluation, where causal graphs were likewise estimated on the full sample size of each real dataset before subsampling for the comparison experiment. This produces three independent baseline fits, one per node count, each with its own causal order $\pi$, adjacency matrix ${\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}$, and residuals $\mathbf{e}=\mathbf{x}-\mathbf{x}{{\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}}^{\top }$ implied by that fit.

We ranked families by mean CLR value, rather than sampling them at random, for two reasons. First, low-abundance taxa are disproportionately affected by zero-inflation and detection-limit noise, which degrades the reliability of any network or causal structure built on them [36]; restricting to the families with the highest mean CLR value at each p favors the more reliably measured taxa, consistent with the prevalence/abundance filter already applied during cohort construction (Section 3.6.1). Second, ranking by mean CLR value produces nested panels across node counts (15-family ⊂ 20-family ⊂ 23-family), so that increasing p consistently means adding progressively rarer, harder-to-estimate families.

${\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}$ was then perturbed into ${\mathbf{B}}^{\left(2\right)}$ using exactly the mechanism-perturbation procedure of Section 2.2.2 (Equation (4)): $n_{\mathrm{pert}}$ target nodes have $n_{\mathrm{pos}}$ parent edges removed and $n_{\mathrm{pos}}$ new parent edges added, with new-edge weights built by drawing a magnitude from the empirical distribution of $|{\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}|$’s nonzero weights and assigning an independent random sign, so injected edges have plausible magnitudes and can be positive or negative, just like the real network’s own edges. What differs from the fully synthetic setup of Section 2.2.1 is where the baseline graph, causal order, and noise come from, and how data are subsequently simulated: instead of a randomly generated graph and fresh i.i.d. noise from a specified non-Gaussian distribution, the baseline ${\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}$ and $\pi$ are estimated directly from the real cohort, and each of the two simulated datasets is generated as ${\mathbf{x}}^{\left(s\right)}={(I-{\mathbf{B}}^{\left(s\right)})}^{-1}{\mathbf{e}}^{\left(s\right)}$, where ${\mathbf{e}}^{\left(s\right)}$ is built by bootstrap-resampling rows of the real residual matrix $\mathbf{e}$. Resampling whole residual rows, rather than resampling each family’s residual independently, preserves whatever real cross-family residual dependence the linear fit did not absorb; using the empirical residual distribution as-is, rather than fitting a parametric noise model to it, avoids introducing a further, potentially misspecified, distributional assumption on top of the real data.

Therefore, it is worthwhile to mention that ${\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}$ is itself a DirectLiNGAM fit, not independently verified ground truth about the true causal structure of the human gut microbiome. Moreover, this design does not test whether DirectLiNGAM recovers the true biological network. It tests whether the Naive, Bootstrap, and RSBS estimators can reliably detect a controlled, known structural change injected on top of a real, data-derived baseline, under real noise, real edge scale, and a real causal order, rather than a synthetic one.

3.6.3. Experimental Design

For every $(p,n_{\mathrm{pert}},n_{\mathrm{pos}})$ combination, twenty independent replicates were run, each drawing a new random target set T, simulating $N_1=300$ samples from ${\mathbf{B}}_{\mathrm{real}}^{\left(1\right)}$ and $N_2=150$ from ${\mathbf{B}}^{\left(2\right)}$ (sample size ratio $r=N_2/N_1=0.5$, held fixed rather than swept, unlike the synthetic heterogeneous experiment of Section 3.4), and scoring all three estimators with 100 bootstrap resamples each. We tested $p\in \{15,20,23\}$ (the last using the full 23-family panel), $n_{\mathrm{pert}}\in \{2,3,4\}$ (equivalently $8.7\%$–$26.7\%$ of nodes perturbed, depending on p), and $n_{\mathrm{pos}}\in \{2,4\}$ parent edges swapped per perturbed node. This gives 18 parameter combinations $\times 20$ repeats $=360$ runs per estimator, $1,080$ discovery-comparison runs in total. Metrics (AUCROC, AUPR, threshold-optimized $F_1$, precision, recall) are computed exactly as in Section 2.2.4.

3.6.4. Results

Table 6 reports overall performance, split by perturbation intensity $n_{\mathrm{pos}}$. Pooled across all 18 conditions (360 paired runs), both resampling-based estimators outperformed Naive by a statistically significant margin: Bootstrap led Naive by $+0.032$ AUCROC ($t=3.50$, $p=4.7\times {10}^{-4}$), $+0.054$ AUPR ($t=5.71$, $p=1.1\times {10}^{-8}$), and $+0.050$$F_1$ ($t=6.59$, $p=4.3\times {10}^{-11}$); RSBS led Naive by $+0.030$ AUCROC ($t=2.79$, $p=5.3\times {10}^{-3}$), $+0.062$ AUPR ($t=5.60$, $p=2.1\times {10}^{-8}$), and $+0.053$$F_1$ ($t=5.92$, $p=3.2\times {10}^{-9}$). Bootstrap and RSBS, in contrast, were not distinguishable from each other pooled over all conditions ($\Delta$AUCROC $=-0.001$, $p=0.90$; $\Delta$AUPR $=+0.008$, $p=0.44$), at the single, moderate sample size ratio ($r=0.5$) used throughout this cohort. This is a narrower separation between the two resampling-based estimators than in the synthetic heterogeneous experiment at the same ratio (Section 3.4, Table 4, where RSBS still led Bootstrap by $+0.062$ AUPR at $r=0.5$), consistent with the additional biological and technical noise present in real, compositional microbiome data relative to a clean synthetic linear SEM.

Figure 7 shows mean AUPR against node count p, faceted by $n_{\mathrm{pos}}$. Absolute performance is markedly lower than on the fully synthetic graphs of Section 3.2, Section 3.3 and Section 3.4 (AUPR in the $0.22$–$0.50$ range here, versus $0.54$–$0.99$ there), which is expected: real, batch-corrected, compositional microbiome data carries genuine biological and technical noise that a clean synthetic linear SEM does not, and should not be expected to match. More informative than the absolute level is the trend across p: Bootstrap’s advantage over the other two estimators, clearest at $p=15$ (AUCROC $0.665$ vs. RSBS $0.626$ vs. Naive $0.592$), erodes as p grows, and by $p=23$ the ranking has reversed on both ranking metrics — RSBS leads (AUCROC $0.609$, AUPR $0.310$), with Bootstrap ($0.569$/$0.245$) statistically indistinguishable from Naive ($0.575$/$0.237$) on AUCROC. This mirrors, within the microbiome-derived semi-synthetic benchmark, the qualitative pattern already established synthetically in Section 3.3: resampling-based confidence becomes more valuable, not less, as the discovery problem gets harder, and here RSBS specifically is the estimator that keeps that promise as network size grows.

Figure 8 repeats the topological stress test of Figure 2 and Figure 5 on the real cohort, faceted by p instead of edge density (the real cohort’s structural axis). The same crossover is visible panel-by-panel: at $p=15$, Bootstrap leads at every perturbation load tested; by $p=23$, RSBS leads at every perturbation load tested, with the gap most pronounced at the lowest perturbation fraction shown ($8.7\%$).

Table 7 reports the paired comparison at the hardest condition tested (largest node count, strongest perturbation: $p=23$, $n_{\mathrm{pos}}=4$). RSBS’s advantage over Naive is significant on AUPR ($+0.074$, $t=2.75$, $p=5.9\times {10}^{-3}$) and $F_1$ ($+0.061$, $t=2.78$, $p=5.5\times {10}^{-3}$) even at this reduced, condition-specific sample size ($n=60$ paired repeats); its advantage over Bootstrap on AUPR is directionally consistent but only borderline significant at this single condition ($+0.050$, $t=1.94$, $p=0.053$), while Bootstrap is statistically indistinguishable from Naive on every metric at this hardest condition, consistent with the crossover already visible in Figure 7.

Finally, Figure 9 examines the effect of perturbation strength $n_{\mathrm{pos}}$ itself on precision and recall, pooled over all three node counts. Increasing $n_{\mathrm{pos}}$ from 2 to 4 raised precision for every method (Naive: $0.341\to 0.368$; Bootstrap: $0.399\to 0.460$; RSBS: $0.413\to 0.457$), while recall stayed essentially flat for all three ($0.75$–$0.80$ throughout). This is the same qualitative pattern as Figure 3 and Figure 6: a larger injected structural change is easier to distinguish from background noise without any corresponding cost in the estimators’ sensitivity, and its replication here, on a real, batch-corrected, compositional cohort rather than a synthetic graph, is itself a validity check on the estimators’ behavior: despite the added difficulty of real, compositional, batch-corrected data, effect size still behaves the way it should.

3.7. Summary of Findings

• Under equal sample sizes, Bootstrap dominated the Naive estimator on AUCROC, AUPR, $F_1$, and precision across every network size, edge density, and perturbation setting tested (best-or-tied in $100\%$ of 99 parameter combinations, vs. $50.5\%$ for Naive).
• Under heterogeneous sample sizes, plain Bootstrap resampling lost its advantage over the Naive method (and was significantly worse on AUCROC and $F_1$), confirming that resampling alone does not compensate for unequal discovery power between datasets of very different sizes.
• RSBS (equal-sample size resampling) was the most robust estimator overall in the heterogeneous experiment, best-or-tied in $97.5\%$ of 198 parameter combinations, with its largest advantage precisely at the most imbalanced sample size ratio tested ($N_2/N_1=0.125$).
• RSBS converges to plain Bootstrap as the sample size ratio approaches 1, as it must by construction, which we take as an internal validity check on our LiNGAM implementation of the three estimators.
• Performance for all three estimators degraded as the number of nodes grew from 15 to 25 at fixed (small) sample sizes, but the relative benefit of resampling-based estimation over the Naive method increased with network size, rather than diminishing.
• Across both experiments, the Naive method consistently achieved the highest recall but the lowest precision; Bootstrap and RSBS traded a small amount of recall for a substantially larger gain in precision, which drove their advantage on $F_1$, AUPR, and AUCROC.
• Because gut microbiome data lack an established causal ground truth and are compositional, with CLR-transformed abundances rank-deficient by construction, evaluating the estimators on real data required a purpose-built, semi-synthetic, real-data-grounded validation (Section 3.6.2) rather than a direct comparison against known structure.
• On this real, batch-corrected microbiome cohort, both resampling-based estimators reproduced the synthetic experiments’ qualitative ranking, outperforming Naive by a modest but statistically significant margin (pooled across node counts and perturbation intensities), though Bootstrap and RSBS were not distinguishable from each other at the single, moderate sample size ratio tested, a narrower separation than observed in the synthetic heterogeneous experiment at the same ratio.

4. Discussion

As the need for mechanistic approaches to discover ecological shifts in microbial communities grows, transitioning from purely associational to causal methods has become essential. Traditional comparative methods rely on correlational summaries. However, these approaches cannot reliably distinguish true mechanistic changes from purely associational differences. To address this limitation, we extended the network comparison framework of Ma and Tourani [22] to detect mechanism shifts at the node level under DirectLiNGAM. We evaluated three estimators, Naive, Bootstrap, and RSBS, that are designed to handle the small and often imbalanced sample sizes typical of microbiome data.

Our empirical evaluation across fully synthetic networks and a semi-synthetic gut microbiome cohort reveals that the optimal estimation strategy depends fundamentally on the sample size. When comparing datasets of equal size, standard Bootstrap estimation consistently dominated the Naive baseline across all tested network sizes, edge densities, and perturbation settings. However, this advantage deteriorated under heterogeneous sample sizes because plain Bootstrap failed to compensate for unequal discovery power. In these imbalanced scenarios, the RSBS estimator was the only method that reliably outperformed the Naive approach. We specifically designed RSBS to equalize sample sizes prior to resampling. As the network size increased and the discovery problem became more difficult, the relative benefit of estimation based on resampling grew. The underlying mechanism for this performance gain is rooted in a trade-off between precision and recall. While the Naive method achieves high recall, it suffers from poor precision by overestimating weakly supported edges. Both Bootstrap and RSBS trade a marginal reduction in recall for a substantial increase in precision. This property proved robust even when subjected to the authentic compositional and technical noise of a real microbiome cohort corrected for batch effects.

Despite these promising results, we must acknowledge several limitations. First, an absolute and independently verified causal network for the human gut microbiome does not currently exist. Therefore, our real-world validation relied on a semi-synthetic ground truth. While our findings confirm that the estimators can detect structural shifts against a baseline derived from data, they do not guarantee that DirectLiNGAM perfectly recovers the true underlying biological network. Second, the framework operates under the assumption of no unobserved confounders. In microbial ecosystems, unmeasured variables like host metabolites or environmental factors can act as hidden common causes that drive the apparent relationship between two species [15]. Third, DirectLiNGAM relies on linearity and acyclicity. In contrast, true microbial interactions frequently feature complex nonlinear dynamics and feedback loops that violate the directed acyclic graph (DAG) assumption [37]. Finally, our current synthetic validation was computationally constrained to a moderate number of nodes. Capturing microbiome dynamics at the species level naturally involves significantly higher dimensionality.

Future work will focus on addressing these dimensional limitations by extending the synthetic benchmarking to accommodate more than 50 variables. This expansion is a necessary prerequisite before applying the pipeline at the species level. It will require expanded computational resources and optimized hardware configurations. Ultimately, we aim to extend the framework along with the Python package to evaluate different classes of structural perturbations in other fields. For example, we could model gene knockout experiments in broader systems biology applications.

5. Conclusions

In this study, we extended the network comparison framework of Ma and Tourani [22] to detect node-level mechanism shifts under DirectLiNGAM. We rigorously evaluated three estimators (Naive, Bootstrap, and RSBS) across nearly 40,000 synthetic discovery runs and a semi-synthetic, batch-corrected human gut microbiome cohort. Our findings demonstrate that estimation based on resampling consistently outperforms a single-fit Naive baseline. Specifically, standard Bootstrap is optimal for equal sample sizes, while RSBS is the only estimator that reliably handles imbalanced cohorts.

The complete benchmark pipeline, named causal-comparator, is openly available as a Python package hosted in a public GitHub repository. We designed this framework specifically to accommodate custom structural perturbations and dynamic synthetic data generation. As a result, researchers can easily adapt and extend the pipeline to generate their own benchmarks with known ground truths. Because of its flexible architecture, causal-comparator is not restricted to the microbiome field, but serves as a general-purpose tool for causal mechanism-shift detection across broader systems biology applications.