May 29, 2026
For many real-world systems, causal ground truth is difficult to obtain, making claims about causal effects hard to assess. We develop methods for evaluating collections of \(\binom{n}{2}\) bivariate causal statements over a set of \(n\) variables. In the setting of acyclic linear statements, any such collection can be extended to a unique multivariate causal model, but we argue that this induced model is implausible if it imposes substantial additional confounding to explain observed correlations. We introduce a compatibility score that quantifies this notion of plausibility, notably without relying on the faithfulness assumption. Additionally, we define an incompatibility score for purely graphical bivariate causal statements, based on global consistency constraints that are derived from acyclicity and faithfulness assumptions. We give theoretical and empirical evidence that both scores can successfully distinguish correct from incorrect causal statements in generic settings. Moreover, we demonstrate the practical applicability of our methods by analyzing causal claims made by large language models. Our work aims to provide a foundation for assessing the reliability of causal information derived from human experts or artificial intelligence in settings where alternative forms of validation are unavailable.
Causal inference in many domains such as medicine or economics still relies heavily on expert knowledge [1]–[3]. Such knowledge is often expressed in the form of pairwise cause-effect statements. Indeed, there is empirical evidence that humans find it easier to assess bivariate causal relations than to infer entire multivariate causal graphs [4]. A similar pattern can be observed when querying large language models (LLMs) for causal relationships [5], [6], an approach that has attracted substantial recent attention [7]–[9]. While LLMs have shown promising performance on some causal reasoning tasks, there are virtually no guarantees regarding the quality of their outputs. Unlike most causal discovery algorithms, which at least come with well-understood consistency guarantees under strong assumptions, most LLMs are not specifically adapted or trained for accurate causal inference [10]. Causal statements from humans can be subject to a variety of well-known biases and are not always based on objective data [11]. This raises a fundamental question: how can we ever trust human- or AI-generated causal statements?
In many real-world settings, obtaining causal ground truth through experimentation is either infeasible or prohibitively expensive, leaving no straightforward answer. The few existing approaches for evaluating causal statements in the absence of ground truth are mostly based on defining and measuring different notions of consistency between an estimated causal graph and the data [12]–[14]. In contrast, [15] consider the compatibility of multiple estimated causal graphs across different subsets of variables to evaluate causal graph learning algorithms. In this paper, we focus on the setting, where only bivariate causal statements are available, and develop notions of compatibility for evaluating and falsifying such statements.
We pursue this approach at two levels of causal information. Given the prevalence of linearity assumptions in causal inference [16], we first focus on falsifying bivariate causal statements that quantify total pairwise causal effects by a linear causal coefficient. In the absence of faithfulness, we show that for any acyclic list of such statements, there exists a unique multivariate causal model that perfectly fits the data and whose marginalizations exactly coincide with the given statements. As a consequence, there are no hard compatibility constraints in the sense of previous work for such lists, beyond consistency with an acyclic causal ordering. We therefore propose to measure compatibility based on the plausibility of the unique induced multivariate causal model. We postulate that a plausible multivariate causal model should exhibit less confounding than its bivariate marginal models, as correlations induced by observed causal backdoor paths in the multivariate model become unobserved confounding in the marginals. We formalize this intuition mathematically, leading to a measure of model plausibility. Theoretically and empirically, we demonstrate that this measure successfully distinguishes correct from incorrect bivariate causal statements. We further apply our approach to causal statements obtained from LLMs for socioeconomic indicators and show that the resulting scores vary substantially across different models and strongly correlate with general model performance.
Second, we consider the case where only qualitative bivariate causal statements are available, indicating the presence or absence of a causal effect and of confounding. Under the faithfulness assumption, such statements are subject to hard graphical compatibility constraints, following [15]. In this setting, we develop an incompatibility score that measures the number of violations of these constraints. In synthetic experiments we investigate how the score reflects the true number of errors in a list of bivariate causal statements and demonstrate our approach for statements generated by LLMs.
In this section, the causal models we consider are linear structural equation models (SEMs). A linear SEM for a vector of observable random variables \({\boldsymbol{X}}= (X_1, \dots X_n)\) is specified by a noise vector of random variables \({\boldsymbol{N}}= (N_1, \dots, N_n)\) and a matrix of causal coefficients \(\Gamma \in \mathbb{R}^{n \times n}\). We assume acyclicity and require \(\Gamma\) to be strictly lower-triangular. Then, the model is given by \[\begin{align} {\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}. \label{eq:SEM} \end{align}\tag{1}\] Since \(\Gamma\) is lower-triangular, \((I - \Gamma)\) is invertible, so \({\boldsymbol{X}}= (I - \Gamma)^{-1} {\boldsymbol{N}}\) is the unique solution to equation 1 . We do not assume the entries of the noise vector to be independent, therefore allowing for correlations representing unmeasured confounding.
The causal interpretation of SEMs is that they define the distribution of \({\boldsymbol{X}}\) under arbitrary interventions (or experiments). An intervention sets variables \(X_j\) for \(j \in S \subseteq \{1, \dots, n\}\) to fixed values \(x_j \in \mathbb{R}\). Then, the remaining variables can be found from solving the equation \[\begin{align} {\boldsymbol{X}}= I_{\overline{S}}\Gamma {\boldsymbol{X}}+ I_{\overline{S}} {\boldsymbol{N}}+ {\boldsymbol{x}}_S.\label{eq:intervention} \end{align}\tag{2}\] Here, \(I_{\overline{S}}\) is the diagonal matrix that has \(1\)’s on the diagonal for indices in the complement \(\overline{S}\) of \(S\), and zeros otherwise, and \({\boldsymbol{x}}_S\) has entries \(x_j\) for \(j \in S\) and zeros elsewhere. We denote the resulting interventional probability distribution of the variables solving equation 2 by \(P({\boldsymbol{X}}\mid \mathop{\mathrm{do}}{\boldsymbol{X}}_S = {\boldsymbol{x}}_S)\). Linear SEMs are closed under marginalization in the following sense:
Lemma 1 (marginalized SEM). Let \({\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}\) be a linear SEM and let \({\boldsymbol{Y}}, {\boldsymbol{Z}}\) be a partition of \({\boldsymbol{X}}\). Then, there exists a noise vector \(\tilde{{\boldsymbol{N}}}\) such that the SEM \[\begin{align} {\boldsymbol{Y}}= \left(\Gamma_{{\boldsymbol{Y}}{\boldsymbol{Y}}} + \Gamma_{{\boldsymbol{Y}}{\boldsymbol{Z}}} (I - \Gamma_{{\boldsymbol{Z}}{\boldsymbol{Z}}})^{-1} \Gamma_{{\boldsymbol{Z}}{\boldsymbol{Y}}}\right) {\boldsymbol{Y}}+ \tilde{{\boldsymbol{N}}}. \end{align}\] induces the same observational and interventional distributions as the original model, marginalized to \({\boldsymbol{Y}}\).
Here, \(\Gamma_{{\boldsymbol{Y}}{\boldsymbol{Z}}}\) denotes the submatrix of \(\Gamma\) whose rows match the indices of variables in \({\boldsymbol{Y}}\) and whose columns match the indices of variables in \({\boldsymbol{Z}}\). For a proof of Lemma 1, see for instance Lemma 7 of [17].
Definition 1 (Bivariate causal statements). A linear bivariate causal statement for a pair \((X_i, X_j)\) specifies
a causal direction, e.g. \(i \to j\);
a linear causal coefficient \(\alpha_{ij} \in \mathbb{R}\).
Our goal is to evaluate a complete list of \(\binom{n}{2}\) linear bivariate causal statements for a fixed set of variables \(\{X_1, \dots, X_n\}\) based on their mutual compatibility. After potentially renaming variables, we assume that all proposed causal directions are compatible with the ordering of \(X_1, \dots, X_n\) given by their indices, that is, each statement with a nonzero coefficient says \(i \to j\) for \(i < j\). Note that this ordering need not be a correct causal ordering. We further assume that the joint probability distribution of \((X_1, \dots, X_n)\) is known. Then, a linear bivariate causal statement assigning the causal coefficient \(\alpha_{ij}\) to \((X_i, X_j)\) equivalently proposes the linear bivariate SEM \[\label{eq:bivariate95sem} X_j = \alpha_{ij} X_i + \tilde{N}_{ij},\tag{3}\] where the distribution of \(\tilde{N}_{ij}\) is just defined as the distribution of \(X_j - \alpha_{ij} X_i\). Crucially, when all proposed causal directions respect a common ordering, there are no further hard compatibility constraints: any complete list of linear bivariate causal statements induces a unique multivariate linear SEM whose pairwise marginals agree with the proposed bivariate SEMs. We encode a list of \(\binom{n}{2}\) bivariate statements by the unit lower-triangular matrix \(A\in\mathbb{R}^{n\times n}\) with entries \[\begin{align} A_{ii} = 1, \qquad A_{ij} = \alpha_{ji} \;\;(i>j), \qquad A_{ij}=0 \;\;(i<j). \end{align}\]
Lemma 2 (Existence of a unique compatible SEM). Let \(A\) be a unit lower-triangular matrix of linear bivariate causal statements for a vector \({\boldsymbol{X}}\) of observed variables and define \(\Gamma = I - A^{-1}\). Then, the multivariate structural equation model \[\begin{align} {\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}, \end{align}\] where the distribution of \({\boldsymbol{N}}\) is defined as the distribution of \((I-\Gamma) {\boldsymbol{X}}\), is the unique linear SEM whose pairwise marginal submodels are \[\begin{align} X_j = A_{ji} X_i + \tilde{N}_{ij}. \end{align}\]
The proof of Lemma 2 only uses basic linear algebra and we provide it in Appendix 6.
Given that any acyclic collection of linear bivariate causal statements is compatible with a unique multivariate SEM, we propose to evaluate such collections based on the plausibility of the induced multivariate model. Clearly, there is no unique way to define plausibility. In this work, we suggest a necessary criterion for plausibility based on the following postulate:
Assumption 1 (Confounding postulate). A generic multivariate causal model should not exhibit more confounding than its pairwise marginal models.
Here, confounding corresponds to correlation between variables that cannot be explained by observed causal effects, and we will formalize this notion in the next section. Intuitively, marginalizing a multivariate causal model to a pairwise model drastically decreases the number of observed causal effects and therefore the amount of confounding in the marginal model should increase. This intuition only fails when the multivariate model is specifically designed so that the correlation induced by its additional observed causal pathways cancels with the multivariate confounding (see Figure 1). Ruling out such cancellations is reminiscent of the standard faithfulness condition (see Appendix 5), but our postulate is logically independent from it. While some near-cancellations of causal pathways can randomly appear [18], a lot more fine-tuning is required for the global confounding of a large multivariate model to be smaller than that of all its pairwise marginals. We provide a theoretical foundation for this claim in Section 2.6 by showing that random causal models satisfy Assumption 1 in expectation, and in Section 2.7, we demonstrate empirically that this is also true with high probability.
To formalize Assumption 1, we need to quantify the amount of confounding in a causal model. There is no universally accepted confounding measure (see [19] for some options). In principle, our approach of falsifying linear bivariate causal statements based on Assumption 1 can be implemented with different choices of confounding measures. We simply need one that enables a fair comparison between a multivariate model and its pairwise marginals. It turns out that quantifying confounding based on squared parts of pairwise covariances between variables that are not explained by observable causal pathways works well for our purpose. Let \(\Sigma\) denote the covariance matrix of a vector of observable variables \({\boldsymbol{X}}\). In a linear bivariate SEM \(X_{j} = \alpha_{ij} X_i + \tilde{N}_{ij}\), the covariance decomposes as \[\begin{align} \Sigma_{ij} = \alpha_{ij} \Sigma_{ii} + \mathop{\mathrm{Cov}}(X_i, \tilde{N}_{ij}). \end{align}\] The first term describes the covariance due to the causal effect from \(X_i\) to \(X_j\), and the second is due to confounding. This motivates the following definition:
Definition 2 (amount of confounding in bivariate models). Let \(\Sigma\) be a covariance matrix and let \(\alpha_{ij}\) be a proposed bivariate causal coefficient from \(X_i\) to \(X_j\). We define the bivariate amount of confounding for the pair \((i,j)\) as \[\begin{align} \mathcal{C}_{ij}^{\text{biv}}(\Sigma, \alpha_{ij}) =\left(\Sigma_{ij} - \alpha_{ij}\Sigma_{ii}\right)^2. \end{align}\]
Here, we use squaring to get a confounding score which is non-negative, and equals zero if and only if the noise term for \(X_j\) is uncorrelated with \(X_i\). Up to normalization, this score has been studied before [20], [21], but only in the bivariate setting. We now extend this idea to multivariate models. Consider a linear SEM \({\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}\) on \(n\) variables with covariance matrix \(\Sigma = \mathop{\mathrm{Cov}}({\boldsymbol{X}})\). For an ordered tuple \(P = (t_1, \dots, t_m)\) with \(1 \leq t_0 < \dots < t_m \leq n\), we define \[\begin{align} \Gamma^P = \prod_{s = 0}^{m-1} \Gamma_{t_{s+1}, t_s} \end{align}\] with the convention that the empty tuple \(P_0 = ()\) gives \(\Gamma^{P_0} = 1\). The tuple \(P\) can be thought of as specifying a directed causal path \(X_{t_0} \to \dots \to X_{t_m}\) and \(\Gamma^P\) is the product of causal coefficients along the path. By Wright’s path tracing rules [22], the covariance between \(X_i\) and \(X_j\) decomposes as
\[\begin{align} \Sigma_{ij} &= \sum_{k \leq i} \Sigma_{kk} \cdot \sum_{\substack{P_1: k \leadsto i, P_2: k \leadsto j,\\ P_1 \cap P_2 = \emptyset}} \Gamma^{P_1} \Gamma ^{P_2} \notag\\ &+ \sum_{\ell \neq k \leq j} \mathop{\mathrm{Cov}}(N_\ell, N_k) \cdot \sum_{\substack{P_1: \ell \leadsto i, P_2: k \leadsto j \\ P_1 \cap P_2 = \emptyset}} \Gamma^{P_1} \Gamma ^{P_2} \label{eq:paths} \end{align}\tag{4}\]
Here, each \(P: k \leadsto i\) runs over all directed paths from \(k\) to \(i\), and two paths are disjoint if they share no vertices, except possibly endpoints. The first term in equation 4 quantifies the contribution of all direct causal paths and fully observed back-door paths to the covariance of \(X_i\) and \(X_j\), whereas the second term is the part of the covariance that cannot be explained by observed causal effects but is due to correlation of the noise variables.
Definition 3 (amount of confounding in multivariate models). Let \(\Sigma\) be a covariance matrix and let \(\Gamma\) be the lower-triangular causal coefficient matrix of a multivariate linear SEM. For \(i < j\), we define the multivariate amount of confounding for the pair \((i,j)\) as \[\begin{align} &\mathcal{C}_{ij}^{\text{mult}}(\Sigma, \Gamma) =\left(\Sigma_{ij} - \sum_{k \leq i} \Sigma_{kk}\sum_{\substack{P_1: k \leadsto i, P_2: k \leadsto j,\\ P_1 \cap P_2 = \emptyset}} \Gamma^{P_1} \Gamma ^{P_2}\right)^2 \end{align}\]
For two variables, this definition coincides exactly with Definition 2.
We can now quantify the compatibility of a list of linear bivariate causal statements given by a matrix \(A\) as the degree to which Assumption 1 holds for the unique induced multivariate model with causal coefficient matrix \(\Gamma = I - A^{-1}\) (see Lemma 2).
Definition 4 (compatibility score). For a list of \(n\) observable variables, let \(\Sigma\) be the covariance matrix and \(A\) be a unit lower-triangular matrix of linear bivariate causal statements. We define the compatibility score: \[\begin{align} \mathop{\mathrm{comp}}(\Sigma, A) = \sum_{1 \leq i < j \leq n} \mathcal{C}_{ij}^\text{biv}(\Sigma, A_{ji}) - \mathcal{C}_{ij}^\text{mult}(\Sigma, I - A^{-1}). \end{align}\]
The exact size of the compatibility score depends on the scaling of the variables. Therefore, in practice, we compute the score after standardizing all variables to unit variance, equivalently replacing \(\Sigma\) by the corresponding correlation matrix. Unless the entries of \(A\) are already standardized causal effects, the matrix \(A\) should be transformed under the same rescaling. However, regardless of standardization, the compatibility score is negative if and only if Assumption 1 is violated (see Figure 1 for an example). In the following sections, we provide evidence that true causal statements typically do not induce models violating Assumption 1. Hence, we view a negative compatibility score as evidence for falsification of the corresponding causal statements, while a positive compatibility score is inconclusive, since the causal statements could simply induce a plausible but wrong causal model.
As a first justification for Assumption 1, we show that the expected compatibility score of true bivariate causal statements is positive for all distributions over causal models that satisfy only a few natural properties. We consider causal models given by linear SEMs with centered Gaussian noise. Such a model is fully specified by the strictly lower-triangular matrix \(\Gamma\) of causal coefficients and a symmetric positive definite noise covariance matrix \(\Sigma_{{\boldsymbol{N}}}\). Note that the covariance matrix \(\Sigma_{\boldsymbol{X}}\) of the observable variables \({\boldsymbol{X}}\), that is needed to compute compatibility scores, is given by \[\begin{align} \Sigma_{\boldsymbol{X}}= (I-\Gamma)^{-1}\Sigma_{\boldsymbol{N}}(I-\Gamma)^{-\top}. \end{align}\]
Assumption 2. We consider probability distributions over \((\Gamma, \Sigma_{\boldsymbol{N}})\) satisfying the following properties:
(Unbiasedness) For all \(i,j\), we have \(\mathop{\mathrm{\mathbb{E}}}[\Gamma_{ij}] = 0\);
(Independence of causal mechanisms 1) The noise covariance matrix and all entries of \(\Gamma\) are mutually independent.
(Non-degeneracy) \(\Sigma_{\boldsymbol{N}}\succ 0\) almost surely and \(\mathop{\mathrm{Var}}(\Gamma_{ij}) > 0\) for all \(i > j\).
Note that these assumptions still allow for arbitrary choices of Gaussian noise, and arbitrary symmetric probability density functions for each causal coefficient.
Theorem 1. For \(n \geq 3\), let \((\Gamma, \Sigma_{\boldsymbol{N}})\) specify a random \(n\)-dimensional linear Gaussian SEM drawn from a distribution that satisfies Assumption 2. Define \(A = (I - \Gamma)^{-1}\) to be the matrix of marginal bivariate causal effects. Then, \[\begin{align} \mathbb{E}_{(\Gamma, \Sigma_{\boldsymbol{N}})}\left[\mathop{\mathrm{comp}}(\Sigma_{\boldsymbol{X}}, A)\right] > 0. \end{align}\]
The proof of Theorem 1 decomposes \(\mathop{\mathrm{comp}}(\Sigma_{\boldsymbol{X}}, A)\) into two components: a squared term, which is strictly positive by the non-degeneracy assumption, and a mixed product of terms whose expectation vanishes due to unbiasedness and independence of the causal coefficients. A full proof is given in Appendix 6.
In practice, we usually do not have access to the precise covariance matrix \(\Sigma\) of the vector \({\boldsymbol{X}}\) of observed variables that is needed to compute the compatibility score of a list of causal statements \(A\). Instead, we use samples \(X^{(1)}, \dots, X^{(N)}\) to compute an empirical covariance matrix \[\begin{align} \widehat{\Sigma} = \frac{1}{N}\sum_{r=1}^N X^{(r)}X^{(r)\top}. \end{align}\] Let us define \[\begin{align} V &:= \max_i \Sigma_{ii}, \qquad a := \max_{i<j} |A_{ji}|, \\ b &:= \max_{i<j} \sum_{k<i} \left|\sum_{\substack{ P_1:k\leadsto i,\;P_2:k\leadsto j\\ P_1\cap P_2=\emptyset }} \Gamma^{P_1}\Gamma^{P_2}\right|, \end{align}\] where \(\Gamma = I - A^{-1}\). Here, \(a\) and \(b\) are parameters that only depend on the matrix \(A\) of bivariate causal statements. While \(a\) simply captures the largest proclaimed total causal effect between two variables, the parameter \(b\) essentially quantifies the largest proclaimed total contribution of causal back-door paths to the covariance between two variables. The following result states that polynomially many samples in \(a,b, V\) and the number of variables \(n\) suffice to approximate the true compatibility score well by using an empirical estimate of the covariance matrix. The proof is based on standard concentration results and can be found in Appendix 6.
Theorem 2. Fix \(0 < {\varepsilon}, \delta < 1\), let \({\boldsymbol{X}}\) be an \(n\)-dimensional centered Gaussian vector with covariance matrix \(\Sigma\) and suppose that \(\widehat{\Sigma}\) is estimated from \[\begin{align} N \geq C \frac{n^4 (1+a+b)^4V^4}{{\varepsilon}^2} \log \frac{n}{\delta} \end{align}\] many iid samples, where \(C\) is a universal constant. Then, with probability at least \(1- \delta\), we have \[\begin{align} \left| \mathop{\mathrm{comp}}(\widehat{\Sigma},A)-\mathop{\mathrm{comp}}(\Sigma,A)\right| \leq \epsilon. \end{align}\]



Figure 2: Percentage of positive compatibility scores for lists of bivariate causal statements on synthetic linear models with increasing amount of error. Each plot contains four different curves for four variations of a model parameter. The fixed model parameters in the first plot are \(n=10\), \(p=0.5\), in the second \(m=3\), \(p=0.5\), and in the third \(n=10\), \(m=3\)..
The source code for all our experiments can be found under https://github.com/ejahn17/compatibility-scores. First, we validate our approach using synthetic experiments. We sample random causal ground truth in the form of linear Gaussian SEMs with \(n\) observed variables, \(m\) hidden variables and sparsity parameter \(p\), which sets causal coefficients to zero with probability \(1-p\) (see Appendix 8.1). Starting from the true bivariate causal statements implied by the ground truth, we generate lists of bivariate statements of decreasing quality by adding independent Gaussian noise with increasing variance \(\sigma\) to the causal coefficients. We then compute compatibility scores for these lists (Figure 2). Each data point in Figure 2 averages over 50 draws of different causal models and 20 draws of the random noise per model. Across all tested combinations of model parameters, the fraction of statements with positive compatibility scores strictly decreases as the error increases, which implies that our compatibility score successfully distinguishes between correct and incorrect statements on average. Our results also provide empirical evidence for the validity of Assumption 1, as most of the true sets of bivariate causal statements have positive compatibility scores. We further observe that the score is robust to variations in the number of hidden variables and that it becomes more successful at distinguishing between correct and incorrect causal statements as the number of observed variables and the density of the true causal model increases. This is consistent with the intuition that both factors increase the number of observed back-door paths in the multivariate model, which are unobserved in the bivariate marginals and therefore increase the difference of the confounding scores. We also examine the differences between compatibility scores that are computed with an empirical covariance matrix based on different numbers of samples and scores computed with the true covariance matrix (see Figure 3). The results suggest that the empirical scores converge faster to the true scores in practice than the worst-case analysis in Theorem 2 suggests. Moreover, less than \(100\) samples are sufficient for the signs of the empirical and the true scores to match in more than \(90\%\) of the test cases. See Appendix 8.4 for additional plots for different parameter combinations.
In a second set of experiments, we compute compatibility scores for lists of linear bivariate causal statements produced by different large language models (LLMs) that we accessed via Amazon Bedrock. For the causal system, we consider seven variables from the gapminder dataset [25]: population density, literacy rate, average daily income, access to sanitation, percentage of adults that smoke, happiness score, life expectancy (see Appendix 8.2 for details). The dataset consists of country-level data points for 179 different countries over the years 1950-2025. We prompt LLMs to choose a causal ordering and estimate total linear causal effects for each pair of variables, based on variable descriptions and the empirical correlation matrix (see Appendix 8.3). We compute compatibility scores averaged over 15 independent runs for each model and compare them to a random baseline, where causal coefficients are drawn from a centered normal distribution with variance matching the variance of the LLM outputs (Figure 4). The results show that our compatibility score can demonstrate systematic differences across models, with higher-capacity models tending to achieve higher scores. Even though the random baseline scores positively in our experiment, many LLMs still receive negative scores, in which case our approach provides evidence for the incorrectness of the generated causal statements.
When quantitative causal statements are unavailable, it can still be desirable to have qualitative causal statements about the existence and direction of causal effects. In this section, we consider acyclic directed mixed graphs (ADMGs) as our causal models. An ADMG is a mixed graph on a set of variables \(V\), that is, a graph that can have both directed edges and bidirected edges, with the restriction that the set of directed edges is acyclic. ADMGs have been widely studied as graphical causal models for causal structures with confounding [26], [27]. Informally, a directed edge \(v \to w\) represents a direct causal effect from \(v\) to \(w\), and a bidirected edge \(v \leftrightarrow w\) indicates the existence of an unobserved causal backdoor path, i.e. confounding, between \(v\) and \(w\). The fact that ADMGs are capable of representing unobserved variables makes them closed under marginalization.
Definition 5. Let \(v,w\) be vertices in an ADMG \(G\). A path between \(v,w\) is a confounding path if both \(v\) and \(w\) are adjacent to an arrowhead of the path and no intermediate vertex is adjacent to two arrowheads (e.g. \(v \leftrightarrow \to \to w\) or \(v \leftarrow \leftrightarrow \to w\)).
Definition 6 (marginal ADMG, see [28]). Let \(G\) be an ADMG on the vertex set \(V \cup L\), where \(V\) and \(L\) are disjoint. The marginal ADMG \(H\) on \(V\) contains all edges of \(G\) within \(V\), and additionally
a directed edge \(v \to w\) if there exists a directed path from \(v\) to \(w\) in \(G\) with all intermediate vertices in \(L\);
a bidirected edge \(v \leftrightarrow w\) if there exists a confounding path from \(v\) to \(w\) in \(G\) with all intermediate vertices in \(L\).
This definition is motivated by the fact that if a normal distribution is Markov and faithful (see Appendix 5) with respect to an ADMG \(G\) on \(V \cup L\), then its marginal distribution on \(V\) is Markov and faithful with respect to the marginal ADMG \(H\) on \(V\) (see Corollaries 7.2 and 7.3 in [29]). We now focus on evaluating lists of bivariate graphical causal statements.
Definition 7 (bivariate statements and statement graph). A bivariate graphical causal statement for two variables \(X_i, X_j\) specifies an ADMG on \(X_i, X_j\), i.e. the existence and direction of a direct causal effect and the existence of confounding. Given a complete list of \(\binom{n}{2}\) such statements for \(n\) variables, we define the statement graph \(\mathop{\mathrm{\mathcal{G}}}\) to be the mixed graph obtained by taking the union of all bivariate ADMGs.
[15] already introduced the concept of graphical compatibility for ADMGs.
Definition 8 (graphical compatibility). A collection of ADMGs on subsets of a common vertex set \(V\) is graphically compatible if they coincide with the marginalizations of a single large ADMG on \(V\).
Specifically for complete lists of bivariate ADMGs, we have the following characterization of graphical compatibility:
Lemma 3. A list of \(\binom{n}{2}\) graphical bivariate causal statements with statement graph \(\mathop{\mathrm{\mathcal{G}}}\) is graphically compatible if and only if
the directed part of \(\mathop{\mathrm{\mathcal{G}}}\) is acyclic;
the directed part of \(\mathop{\mathrm{\mathcal{G}}}\) is transitively closed, i.e. if there is a directed path from \(X_i\) to \(X_j\), then there must be an edge from \(X_i\) to \(X_j\);
if there is a confounding path between \(X_i\) and \(X_j\), then there must be a bidirected edge between \(X_i\) and \(X_j\).
The proof is straightforward and can be found in Appendix 6.
Let \(d(\mathop{\mathrm{\mathcal{G}}}, \mathop{\mathrm{\mathcal{G}}}^*)\) denote the Hamming distance between two mixed graphs, that is the number of directed and bidirected edge deletions and additions needed to transform \(\mathop{\mathrm{\mathcal{G}}}\) into \(\mathop{\mathrm{\mathcal{G}}}^*\).
Definition 9. For a statement graph \(\mathop{\mathrm{\mathcal{G}}}\) summarizing a list of \(\binom{n}{2}\) bivariate graphical causal statements, we define the incompatibility score \[\begin{align} \mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}}) = \min_{\mathop{\mathrm{\mathcal{G}}}^*} d(\mathop{\mathrm{\mathcal{G}}}, \mathcal{G}^* ), \end{align}\] where the minimum ranges over all mixed graphs \(\mathop{\mathrm{\mathcal{G}}}^*\) that satisfy the three properties stated in Lemma 3.
There is an important distinction from the graphical incompatibility score defined by [15]. While their paper assesses compatibility with a specific large ADMG, that is inferred by the same causal discovery algorithm as the respective small ADMGs, we do not assume access to an inferred multivariate ADMG. Instead, we implicitly optimize over the multivariate ADMG, by comparing \(\mathop{\mathrm{\mathcal{G}}}\) with the closest possible compatible set of statements. However, this generalization makes the problem of computing \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) \(\mathsf{NP}\)-hard. This is because it contains the following \(\mathsf{NP}\)-hard problem.
Acyclic Transitivity Editing: Given an acyclic directed graph \(G\), find the minimum number of edge deletions and additions to make it transitively closed. See [30] for a hardness proof.
Indeed, computing \(\mathop{\mathrm{incomp}}(G)\) for an acyclic directed graph after adding all possible bidirected edges is equivalent to solving the transitivity editing problem. Still, this does not mean any attempt at computing \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) is fruitless. Intuitively, the \(\mathsf{NP}\)-hardness comes from regimes where \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) is large. Indeed, a brute-force algorithm can decide in polynomial time if it is smaller than a given constant. But when \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) is large, a precise computation is unnecessary for our purpose: any evidence for \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) exceeding a certain threshold implies that the graph \(\mathop{\mathrm{\mathcal{G}}}\) is highly incompatible, which is already useful for falsification. Motivated by this insight, we present a heuristic algorithm for approximating \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\).
The exact solution for the minimization problem in Definition 9 depends on the complex interaction between all three properties of Lemma 3. As a first heuristic, we decouple these constraints and address them sequentially. Given a mixed graph \(\mathop{\mathrm{\mathcal{G}}}\) with directed part \(D = D(\mathop{\mathrm{\mathcal{G}}})\), we proceed as follows:
First, we find a small set of directed edges whose deletion makes \(D\) acyclic. This is exactly the Minimum Feedback Arc problem, for which there is a heuristic algorithm called GreedyFAS, which runs in linear time and has good empirical performance [31], [32]. The algorithm iteratively constructs a vertex ordering by picking sources and sinks whenever possible, otherwise, it greedily picks the vertex with maximal difference between out- and in-degree instead of a source. Its output set \(E_{\mathrm{cycles}}\) is the set of edges that is not consistent with the obtained ordering (see Appendix 7 for more details).
Set \(D' = D \setminus E_{\mathrm{cycles}}\) to be the directed acyclic graph obtained in step 1. We design a greedy, heuristic algorithm to solve the Acyclic Transitivity Editing problem for \(D'\). Let the transitive closure \(tc(D')\) be the smallest directed graph that contains \(D'\) and is transitively closed. Our algorithm GreedyTE iteratively deletes the edge that reduces the Hamming distance from the current graph to its transitive closure as much as possible, resulting in the deletion of \(E_\mathrm{del} \subseteq E(D')\). Once no more reduction is possible, the algorithm adds the set of edges \(E_\mathrm{add}\) that is missing to get to \(tc(D' \setminus E_\mathrm{del})\) (see Appendix 7).
After steps 1 and 2, the directed part \(D\) of the statement graph \(\mathop{\mathrm{\mathcal{G}}}\) has been transformed into a DAG \(D' = (D \setminus (E_{\mathrm{cycles}} \cup E_{\mathrm{del}})) \cup E_{\mathrm{add}}\) that is transitively closed. Now, we keep \(D'\) fixed and only operate on the bidirected edges of \(\mathop{\mathrm{\mathcal{G}}}\). Define the confounding path closure \(cpc(G)\) of a mixed graph \(G\) to be the graph obtained from \(G\) by adding bidirected edges between each pair of vertices that is connected by a confounding path. We then apply a greedy procedure analogous to Step 2, where the transitive closure is replaced by \(cpc(G)\), to identify bidirected edge sets \(B_\mathrm{del}\) and \(B_\mathrm{add}\) whose deletion and addition makes the resulting graph satisfy property [prop:cpc] of Lemma 3. Again, the full algorithm is given in Appendix 7.
We report \[\begin{align} c(\mathop{\mathrm{\mathcal{G}}}) := |E_{\mathrm{cycles}}| + |E_{\mathrm{add}}| + |E_{\mathrm{del}}| + |B_{\mathrm{add}}| + |B_{\mathrm{del}}| \label{eq:heuristic-score} \end{align}\tag{5}\] as a heuristic approximation for \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\). This approximation satisfies the following properties:
Lemma 4. For any statement graph \(\mathop{\mathrm{\mathcal{G}}}\), we have
\(c(\mathop{\mathrm{\mathcal{G}}}) \geq \mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\);
\(c(\mathop{\mathrm{\mathcal{G}}}) = 0 \iff \mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}}) = 0\).
Hence, whenever \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) is either zero or it is large, this reflects correctly in our heuristic approximation. We give a proof of Lemma 4 in Appendix 6.



Figure 5: Average incompatibility scores for lists of bivariate causal statements on synthetic graphical models with increasing numbers of errors. The fixed parameters in the first plot are \(n=10\), \(p=0.3\), in the second \(m=3\), \(p=0.3\), and in the third \(n=10\), \(m=3\)..
Again, we start with synthetic experiments and sample causal ground truth with \(n\) observed variables, \(m\) hidden variables and sparsity parameter \(p\) (see Appendix 8.1). For each list of correct bivariate causal statements, we inject a fixed number of random errors, by randomly selecting an edge out of all possible directed and bidirected edges and deleting or inserting it, depending on whether it was present or absent. Then we compute incompatibility scores using the heuristic approach of Section 3.3, see Figure 5. As expected, incompatibility scores increase on average monotonically with the number of injected errors. This suggests that the incompatibility score can serve as a tool for comparing the quality of different lists of graphical bivariate causal statements for the same causal graph. For small or sparse graphs, the scores also approximate the true numbers of errors quite closely on average, and they are robust against variations in the number of hidden confounders. However, for larger and denser models, the scores tend to overestimate the true error count, which indicates that the approximation quality of our heuristic algorithm decreases as the number of injected errors is an upper bound for the exact incompatibility score.
We also apply our approach to graphical bivariate causal statements obtained from LLMs, using the same dataset and models as in Section 2.7 (see Appendix 8.3 for prompt details). The results are shown in Figure 6. Among all derived statement graphs, only one achieves \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}}) = 0\), and it is complete, meaning that it only needs to be acyclic to be compatible. Indeed, several low-scoring LLM outputs are extremely dense statement graphs, which achieve low incompatibility scores, but are not very informative2. This is why we also report incompatibility scores for statement graphs filtered by a capped edge density, see Figure 7. Among this subset of more informative statement graphs, the results again reveal a correlation between low incompatibility scores and higher general model capacity. We also provide a more detailed breakdown of incompatibility scores in Appendix 8.4.
We developed two complementary approaches to assess bivariate causal statements based on their mutual compatibility. For graphical statements, we approximate the minimum number of statements that need to change to give rise to a single consistent multivariate causal graph under the assumptions of acyclicity and faithfulness. In the absence of faithfulness and when quantitative linear causal statements are available, we measure compatibility by evaluating how well the unique induced linear multivariate model satisfies our plausibility criterion, namely that it should exhibit less confounding than its bivariate marginals. Importantly, a positive compatibility score does not imply correctness of the causal statements, but a negative score provides evidence for internal inconsistencies to some degree.
There is potential for future work that extends our framework to cases where the collection of bivariate causal statements may be incomplete, or accompanied by statements over small subsets of variables. Especially graphical compatibility seems amenable to such a generalization, while extending compatibility for linear statements may require new ideas, as the existence of a unique composite multivariate model can no longer be guaranteed.
Another interesting direction would be to test the application of our compatibility scores for causal inference in the setting where incomplete lists of bivariate causal statements are already available. The missing statements could then be obtained by maximizing the compatibility score (or minimizing the incompatibility score).
Finally, compatibility scores can provide a meaningful quantitative comparison between collections of causal claims in the absence of causal ground truth. Thus, our methods could guide the automated evaluation of causal statements, for example as part of post-training, filtering, or validation pipelines for LLMs. This could be a step towards obtaining more reliable causal information from AI models that can effectively support human decision-making.
This paper describes an approach for assessing the reliability of AI-generated causal statements. Our goal is to enable future improvements of their accuracy. To mitigate the risk of over-interpreting compatibility scores, we want to emphasize that a high score does not provide a proof of correctness of the respective causal statements, it only indicates that they pass a sanity check.
Below, we summarize how causal graphical models such as ADMGs can be linked to probability distributions (see [29]):
Definition 10 (Collider). Given a path in an ADMG \(G\), a collider is an interior vertex \(x\) on the path such that the edges preceding and succeeding \(x\) have arrowheads pointing to \(x\) (that is, \(\to x \leftarrow\) or \(\leftrightarrow x \leftarrow\) or \(\to x \leftrightarrow\) or \(\leftrightarrow x \leftrightarrow\)).
Definition 11 (\(m\)-separation). Let \(G\) be an ADMG with disjoint vertex sets \(X,Y, Z\). Then, \(X\) and \(Y\) are said to be \(m\)-connected given \(Z\) if there exists a path connecting a vertex in \(X\) with a vertex in \(Y\) such that
each non-collider on the path is not in \(Z\);
for each collider on the path, there exists a directed path to a vertex in \(Z\).
If \(X\) and \(Y\) are not \(m\)-connected given \(Z\), they are said to be \(m\)-separated given \(Z\).
Definition 12 (Markov condition). A probability distribution \(P\) on \(V\) satisfies the Markov condition with respect to an ADMG \(G\) with vertex set \(V\) if for any disjoint subsets \(X,Y, Z\), we have \[\begin{align} \text{X and Y are m-separated given Z in G} \implies \text{X and Y are independent conditioned on Z in P.} \end{align}\]
Definition 13 (Faithfulness condition). A probability distribution \(P\) on \(V\) satisfies the faithfulness condition with respect to an ADMG \(G\) with vertex set \(V\) if for any disjoint subsets \(X,Y, Z\), we have \[\begin{align} \text{X and Y are independent conditioned on Z in P} \implies \text{X and Y are m-separated given Z in G.} \end{align}\]
The Markov and faithfulness conditions are standard assumptions for the task of inferring a causal graph from an empirical probability distribution. For linear SEMs \({\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}\) with Gaussian noise, one can show that the distribution of \({\boldsymbol{X}}\) always satisfies the Markov condition with respect to the corresponding ADMG on \({\boldsymbol{X}}\) that is defined by setting \(X_i \to X_j\) if and only if \(\Gamma_{ji} \neq 0\) and \(X_i \leftrightarrow X_j\) if and only if \(\mathop{\mathrm{Cov}}(N_i, N_j) \neq 0\). In contrast, the faithfulness condition is not always satisfied, even for linear Gaussian models (which is the case when the contributions of multiple \(m\)-connecting paths to the covariance between two variables exactly cancel), and we do not need this condition for our approach developed in Section 2. However, the faithfulness condition is needed to ensure that a causal graph inferred from a marginal probability distribution coincides with the marginal causal graph that is inferred from the full probability distribution. The soundness of Definition 8 relies on this fact, and therefore our approach in Section 3 implicitly relies on the faithfulness assumption.
Lemma 5 (Existence of a unique compatible SEM). Let \(A\) be a unit lower-triangular matrix of linear bivariate causal statements for a vector \({\boldsymbol{X}}\) of observed variables and define \(\Gamma = I - A^{-1}\). Then, the multivariate structural equation model \[\begin{align} {\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}, \end{align}\] where the distribution of \({\boldsymbol{N}}\) is defined as the distribution of \((I-\Gamma) {\boldsymbol{X}}\), is the unique linear SEM whose pairwise marginal submodels are \[\begin{align} X_j = A_{ji} X_i + \tilde{N}_{ij}. \end{align}\]
Proof. Consider an arbitrary SEM \({\boldsymbol{X}}= \Gamma {\boldsymbol{X}}+ {\boldsymbol{N}}\). By Lemma 1, marginalizing to the pair \((X_i, X_j)\) with \(i < j\) results in the SEM \[\begin{align} \begin{pmatrix} X_i \\X_j \end{pmatrix} = (\Gamma_{YY} + \Gamma_{YZ} (I - \Gamma_{ZZ})^{-1} \Gamma_{ZY}) \begin{pmatrix} X_i \\X_j \end{pmatrix} + \tilde{{\boldsymbol{N}}}, \end{align}\] where \(Y = {i,j}\) and \(Z = \{1, \dots, n\} \setminus \{i,j\}\). Since \(\Gamma\) is strictly lower-triangular, we have that \[\begin{align} \Gamma_{YY} + \Gamma_{YZ} (I - \Gamma_{ZZ})^{-1} \Gamma_{ZY} = \begin{pmatrix} 0 & 0\\ r & 0 \end{pmatrix}, \end{align}\] for some entry \(r \in \mathbb{R}\). On the other hand, block matrix inversion of \((I-\Gamma)\) using the Schur complement tells us that \[\begin{align} (I - \Gamma)^{-1}_{YY} = (I- \Gamma_{YY} - \Gamma_{YZ}(I-\Gamma_{ZZ})^{-1}\Gamma_{ZY})^{-1} = \begin{pmatrix} 1 & 0\\ -r & 1 \end{pmatrix}^{-1} = \begin{pmatrix} 1 & 0\\ r & 1 \end{pmatrix}. \end{align}\] This implies \(r = (I - \Gamma)^{-1}_{ji}\), so the pairwise marginal submodels of an SEM with causal coefficients \(\Gamma\) are given by \[\begin{align} X_j = (I - \Gamma)^{-1}_{ji} X_i + \tilde{N}_{ij}. \end{align}\] Requiring that all these submodels align with \(X_j = \alpha_{ij} X_i + \tilde{N}_{ij}\) is therefore equivalent to \((I- \Gamma)^{-1} = A\) or \(\Gamma = I - A^{-1}\), which completes the proof. ◻
Theorem 3. For \(n \geq 3\), let \((\Gamma, \Sigma_{\boldsymbol{N}})\) specify a random \(n\)-dimensional linear Gaussian SEM drawn from a distribution that satisfies Assumption 2. Define \(A = (I - \Gamma)^{-1}\) to be the matrix of marginal bivariate causal effects. Then, \[\begin{align} \mathbb{E}_{(\Gamma, \Sigma_{\boldsymbol{N}})}\left[\mathop{\mathrm{comp}}(\Sigma_{\boldsymbol{X}}, A)\right] > 0. \end{align}\]
Proof. For \(1 \leq i < j \leq n\), we define \[\begin{align} B_{ij} := \sum_{k < i} \Sigma_{{\boldsymbol{X}}, kk} \cdot \sum_{\substack{P_1: k \leadsto i \\ P_2: k \leadsto j\\ P_1 \cap P_2 = \emptyset}} \Gamma^{P_1} \Gamma ^{P_2} \end{align}\] to be the part of the covariance between \(X_i, X_j\) coming from back-door paths through observed variables and \[\begin{align} {\varepsilon}_{ij} := \sum_{\ell \neq k \leq j} \Sigma_{{\boldsymbol{N}}, lk} \cdot \sum_{\substack{P_1: \ell \leadsto i\\ P_2: k \leadsto j \\ P_1 \cap P_2 = \emptyset}} \Gamma^{P_1} \Gamma ^{P_2} \end{align}\] to be the part of the covariance between \(X_i, X_j\) coming from back-door paths that are not fully observed. First, we show that \(B_{ij}\) and \({\varepsilon}_{ij}\) are uncorrelated. For \(i=1\), we have \(B_{ij} = 0\), so this is trivially true. For \(i >1\), we get \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[B_{ij} {\varepsilon}_{ij}] = \sum_{k<i} \sum_{\ell \neq m \leq j} \sum_{\substack{P_1: k \leadsto i \\ P_2: k \leadsto j\\ P_1 \cap P_2 = \emptyset}}\sum_{\substack{Q_1: \ell \leadsto i\\ Q_2: m \leadsto j \\ Q_1 \cap Q_2 = \emptyset}} \mathop{\mathrm{\mathbb{E}}}\left[\Sigma_{{\boldsymbol{X}}, kk} \cdot \Sigma_{{\boldsymbol{N}}, lm} \cdot \Gamma^{P_1}\Gamma^{P_2}\Gamma^{Q_1}\Gamma^{Q_2}\right] \label{eq:correlation} \end{align}\tag{6}\] Fix a summand in the expression above, determined by specific choices of \(k,\ell, m, P_1, P_2, Q_1, Q_2\). Since \(Q_1, Q_2\) are disjoint paths with all different endpoints, they cannot fully contain the paths \(P_1, P_2\) as those intersect in the common endpoint \(k\). Hence, there must exist an entry \(\Gamma_{rs}\) with \(r >k\) that only appears once in the product \(\Gamma^{P_1}\Gamma^{P_2}\Gamma^{Q_1}\Gamma^{Q_2}\). Note that \(\Sigma_{{\boldsymbol{X}}, kk}\) only depends on \(\Sigma_{\boldsymbol{N}}\) and entries \(\Gamma_{pq}\) with \(p \leq k\). By part [part:independence] of Assumption 2, the entry \(\Gamma_{rs}\) is therefore independent from all other factors in the product \(\Sigma_{{\boldsymbol{X}}, kk} \cdot \Sigma_{{\boldsymbol{N}}, lm} \cdot \Gamma^{P_1}\Gamma^{P_2}\Gamma^{Q_1}\Gamma^{Q_2}\). This implies \[\begin{align} \mathop{\mathrm{\mathbb{E}}}\left[\Sigma_{{\boldsymbol{X}}, kk} \cdot \Sigma_{{\boldsymbol{N}}, lm} \cdot \Gamma^{P_1}\Gamma^{P_2}\Gamma^{Q_1}\Gamma^{Q_2}\right] = \mathop{\mathrm{\mathbb{E}}}\left[\Gamma_{rs}\right] \cdot \mathop{\mathrm{\mathbb{E}}}\left[\Sigma_{{\boldsymbol{X}}, kk} \cdot \Sigma_{{\boldsymbol{N}}, lm} \cdot \Gamma^{P_1}\Gamma^{P_2}\Gamma^{Q_1}\Gamma^{Q_2} \cdot \frac{1}{\Gamma_{rs}}\right] = 0, \end{align}\] by part [part:bias] of Assumption 2. Since this argument holds for all summands of equation 6 , we get \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[B_{ij} {\varepsilon}_{ij}] = 0. \end{align}\] Now, recall equation 4 stating that \[\begin{align} \Sigma_{{\boldsymbol{X}}, ij} = \Sigma_{{\boldsymbol{X}}, ii} \cdot A_{ji} + B_{ij} + {\varepsilon}_{ij}. \end{align}\] Therefore, we get for the bivariate confounding between \(X_i\) and \(X_j\) \[\begin{align} \mathcal{C}_{ij}^\text{biv}(\Sigma_{\boldsymbol{X}}, A_{ji}) = (B_{ij} + {\varepsilon}_{ij})^2 \end{align}\] and the multivariate confounding is given by \[\begin{align} \mathcal{C}_{ij}^\text{mult}(\Sigma_{\boldsymbol{X}}, \Gamma) = {\varepsilon}_{ij}^2. \end{align}\] Hence, we can write \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[\mathop{\mathrm{comp}}(\Sigma_{\boldsymbol{X}}, A)] = \mathop{\mathrm{\mathbb{E}}}\left[\sum_{i < j} (B_{ij} + {\varepsilon}_{ij})^2 - {\varepsilon}_{ij}^2\right] = \sum_{i<j} \left(\mathop{\mathrm{\mathbb{E}}}\left[B_{ij}^2 \right] + 2 \mathop{\mathrm{\mathbb{E}}}\left[B_{ij} {\varepsilon}_{ij}\right] \right) = \sum_{i<j} \mathop{\mathrm{\mathbb{E}}}[B_{ij}^2] \geq 0. \end{align}\] For \(i = 2, j=3\), we get from part [part:independence] of Assumption 2 \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[B_{23}^2] = \mathop{\mathrm{\mathbb{E}}}[\Sigma_{{\boldsymbol{X}}, 11}^2 \cdot \Gamma_{21}^2 \cdot \Gamma_{31}^2] = \mathop{\mathrm{\mathbb{E}}}[\Sigma_{{\boldsymbol{N}}, 11}^2 \cdot \Gamma_{21}^2 \cdot \Gamma_{31}^2] = \mathop{\mathrm{\mathbb{E}}}[\Sigma_{{\boldsymbol{N}}, 11}^2] \cdot \mathop{\mathrm{\mathbb{E}}}[\Gamma_{21}^2] \cdot \mathop{\mathrm{\mathbb{E}}}[\Gamma_{31}^2]. \end{align}\] By part [part:nondegenerate] of Assumption 2, we have that \(\Sigma_{{\boldsymbol{N}}, 11}^2 > 0\) almost surely, and \(\mathop{\mathrm{\mathbb{E}}}[\Gamma_{21}^2] = \mathop{\mathrm{Var}}(\Gamma_{21}) > 0, \mathop{\mathrm{\mathbb{E}}}[\Gamma_{31}^2] = \mathop{\mathrm{Var}}(\Gamma_{31}) > 0\). This implies \(\mathop{\mathrm{\mathbb{E}}}[B_{23}^2] > 0\) and therefore also \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[\mathop{\mathrm{comp}}(\Sigma_{\boldsymbol{X}}, A)] > 0. \end{align}\] ◻
Theorem 4. Fix \(0 < {\varepsilon}, \delta < 1\), let \({\boldsymbol{X}}\) be an \(n\)-dimensional centered Gaussian vector with covariance matrix \(\Sigma\) and suppose that \(\widehat{\Sigma}\) is estimated from \[\begin{align} N \geq C \frac{n^4 (1+a+b)^4V^4}{{\varepsilon}^2} \log \frac{n}{\delta} \end{align}\] many iid samples, where \(C\) is a universal constant. Then, with probability at least \(1- \delta\), we have \[\begin{align} \left| \mathop{\mathrm{comp}}(\widehat{\Sigma},A)-\mathop{\mathrm{comp}}(\Sigma,A)\right| \leq \epsilon. \end{align}\]
Proof. Since \(X\) is centered Gaussian and \(\max_i \Sigma_{ii}\leq V\), a standard concentration result (see [34], Lemma 1) shows that there exists a universal constant \(c>0\) such that, for all \(0<\gamma\leq V\), \[\begin{align} \Pr\left[ \max_{p,q} |\widehat{\Sigma}_{pq}-\Sigma_{pq}|>\gamma \right] \leq 4n^2\exp\left( -cN\frac{\gamma^2}{V^2} \right). \end{align}\]
Let \(\mathcal{E}_\gamma\) denote the event \[\begin{align} \mathcal{E}_\gamma := \left\{\max_{p,q}|\widehat{\Sigma}_{pq}-\Sigma_{pq}| \leq \gamma \right\}. \end{align}\]
Conditional on \(\mathcal{E}_\gamma\), we prove a bound on \(| \mathop{\mathrm{comp}}(\widehat{\Sigma},A)-\mathop{\mathrm{comp}}(\Sigma,A)|\). For \(i<j\) and fixed \(A\), define the bivariate residual \[\begin{align} B_{ij}(\Sigma) := \Sigma_{ij}-A_{ji}\Sigma_{ii}, \end{align}\] and the multivariate residual \[\begin{align} M_{ij}(\Sigma) := \Sigma_{ij}-\sum_{k\leq i}\Sigma_{kk} \sum_{\substack{ P_1:k\leadsto i,\;P_2:k\leadsto j\\ P_1\cap P_2=\emptyset }}\Gamma^{P_1}\Gamma^{P_2}, \end{align}\] where \(\Gamma = I - A^{-1}\). The term corresponding to \(k=i\) equals \(A_{ji}\Sigma_{ii}\), because it sums over all directed paths from \(i\) to \(j\). Hence \[\begin{align} M_{ij}(\Sigma) = \Sigma_{ij}-A_{ji}\Sigma_{ii}-\sum_{k<i}\Sigma_{kk} \sum_{\substack{ P_1:k\leadsto i,\;P_2:k\leadsto j\\ P_1\cap P_2=\emptyset }}\Gamma^{P_1}\Gamma^{P_2}. \end{align}\] Let \(K := 1 + a + b\) with \(a,b\) as defined in Section 2.6. Then, \(B_{ij}\) and \(M_{ij}\) are linear functions of the entries of \(\Sigma\) with coefficient \(\ell_1\)-norm at most \(K\). On the event \(\mathcal{E}_\gamma\), this implies \[\begin{align} |B_{ij}(\widehat{\Sigma})-B_{ij}(\Sigma)|\leq K\gamma, \qquad |M_{ij}(\widehat{\Sigma})-M_{ij}(\Sigma)|\leq K\gamma. \end{align}\] Since \(|\Sigma_{pq}|\leq V\) for all \(p,q\) by Cauchy–Schwarz, we also have \[\begin{align} |B_{ij}(\Sigma)|,\;|M_{ij}(\Sigma)| \leq KV, \end{align}\] and on \(\mathcal{E}_\gamma\), assuming \(\gamma \leq V\), we have \(|\widehat{\Sigma}_{pq}|\leq 2V\), and therefore \[\begin{align} |B_{ij}(\widehat{\Sigma})|,\;|M_{ij}(\widehat{\Sigma})|\leq 2KV. \end{align}\] Using \(|x^2-y^2|=|x-y||x+y|\), we get \[\begin{align} \left|\mathcal{C}_{ij}^\text{biv}(\widehat{\Sigma}, A_{ji})-\mathcal{C}_{ij}^\text{biv}(\Sigma, A_{ji})\right| = \left|B_{ij}(\widehat{\Sigma})^2 - B_{ij}(\Sigma)^2\right| \leq K \gamma \cdot 3KV = 3K^2V \gamma, \end{align}\] and similarly, \[\begin{align} \left|\mathcal{C}_{ij}^\text{mult}(\widehat{\Sigma}, \Gamma) - \mathcal{C}_{ij}^\text{mult}(\Sigma, \Gamma)\right|= \left|M_{ij}(\widehat{\Sigma})^2 - M_{ij}(\Sigma)^2\right| \leq 3K^2V \gamma. \end{align}\] Thus, for each pair \(i<j\), we obtain \[\begin{align} \left|\left(\mathcal{C}_{ij}^\text{biv}(\widehat{\Sigma}, A_{ji}) - \mathcal{C}_{ij}^\text{mult}(\widehat{\Sigma}, \Gamma)\right) - \left(\mathcal{C}_{ij}^\text{biv}(\Sigma, A_{ji}) - \mathcal{C}_{ij}^\text{mult}(\Sigma, \Gamma)\right)\right| \leq 6K^2V\gamma. \end{align}\] Summing over at most \(n^2/2\) pairs gives \[\begin{align} \left| \mathop{\mathrm{comp}}(\widehat{\Sigma},A)-\mathop{\mathrm{comp}}(\Sigma,A)\right| \leq 3n^2K^2V\gamma = {\varepsilon}, \end{align}\] where the last step follows by choosing \(\gamma = \frac{\epsilon}{3n^2K^2V}\). Finally, the Gaussian covariance concentration bound shows that \(\mathcal{E}_\gamma\) holds with probability at least \(1-\delta\) whenever \[\begin{align} N \geq C'\frac{V^2}{\gamma^2} \log \frac{n}{\delta} = C \frac{n^4K^4V^4}{{\varepsilon}^2}\log\frac{n}{\delta} \end{align}\] for some sufficiently large constant \(C\). This completes the proof. ◻
Lemma 6. A list of \(\binom{n}{2}\) graphical bivariate causal statements with statement graph \(\mathop{\mathrm{\mathcal{G}}}\) is graphically compatible if and only if
the directed part of \(\mathop{\mathrm{\mathcal{G}}}\) is acyclic;
the directed part of \(\mathop{\mathrm{\mathcal{G}}}\) is transitively closed, i.e. if there is a directed path from \(X_i\) to \(X_j\), then there must be an edge from \(X_i\) to \(X_j\);
if there is a confounding path between \(X_i\) and \(X_j\), then there must be a bidirected edge between \(X_i\) and \(X_j\).
Proof. First, suppose the statement graph \(\mathop{\mathrm{\mathcal{G}}}\) satisfies the three conditions of Lemma 3. By the first condition, \(\mathop{\mathrm{\mathcal{G}}}\) is an ADMG itself. By the second condition, all its directed edges coincide with the directed edges of its pairwise
marginalizations and by the third condition, all its bidirected edges coincide with the bidirected edges of its pairwise marginalizations. Hence, \(\mathop{\mathrm{\mathcal{G}}}\) itself is a multivariate ADMG that is
compatible with all bivariate statements.
Secondly, suppose the statement graph \(\mathop{\mathrm{\mathcal{G}}}\) is graphically compatible, so there exists a multivariate ADMG \(G\) whose pairwise marginalizations coincide with the
edges of \(\mathop{\mathrm{\mathcal{G}}}\). Assume that all vertices are ordered consistently with respect to the directed edges of \(G\). Then, also the directed edges of the pairwise
marginalizations of \(G\) (which are obtained from directed paths), must be consistent with the vertex ordering. Hence, the directed part of \(\mathop{\mathrm{\mathcal{G}}}\) is acyclic,
showing that the first condition of Lemma 3 holds. Next, suppose there is a directed path from \(X_i\) to \(X_j\) in \(\mathop{\mathrm{\mathcal{G}}}\). Every directed edge \(X_k \to X_\ell\) in the path corresponds to either a directed edge \(X_k \to X_\ell\) in \(G\) or a directed path from \(X_k\) to \(X_\ell\). Concatenating all these directed paths gives a directed
path from \(X_i\) to \(X_j\) in \(G\). Hence, marginalizing \(G\) onto the vertex set \(\{X_i,
X_j\}\) results in an edge \(X_i \to X_j\), which proves that \(\mathop{\mathrm{\mathcal{G}}}\) satisfies the second property of Lemma 3. Finally, assume that there is a confounding path between \(X_i, X_j\) in \(\mathop{\mathrm{\mathcal{G}}}\). As before,
each directed edge on this path must correspond to a directed path in \(G\) and if there is a bidirected edge on the path, then it must correspond to a confounding path in \(G\).
Concatenating these directed paths and the confounding path, gives a confounding path between \(X_i\) and \(X_j\) in \(G\), which implies the existence of a
bidirected edge \(X_i \leftrightarrow X_j\) in the pairwise marginal ADMG. This shows that \(\mathop{\mathrm{\mathcal{G}}}\) also satisfies the third property of Lemma 3 and completes the proof. ◻
Lemma 7. For any statement graph \(\mathop{\mathrm{\mathcal{G}}}\), we have
\(c(\mathop{\mathrm{\mathcal{G}}}) \geq \mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}})\) ;
\(c(\mathop{\mathrm{\mathcal{G}}}) = 0 \iff \mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}}) = 0\).
Proof. For part [part-1] of Lemma 4, it suffices to show that \(\mathop{\mathrm{\mathcal{G}}}^* = \mathop{\mathrm{\mathcal{G}}}\setminus (E_{\text{cycles}}\cup E_{\text{del}} \cup B_{\text{del}}) \cup E_{\text{add}} \cup B_{\text{add}}\) is graphically compatible. Let \(\mathop{\mathrm{\mathcal{G}}}' := \mathop{\mathrm{\mathcal{G}}}\setminus (E_{\text{cycles}}\cup E_{\text{del}})\). By definition, \(B_{\text{add}}\) completes the confounding path closure of \(\mathop{\mathrm{\mathcal{G}}}' \setminus B_{\text{del}} \cup E_{\text{add}}\). Hence, \(\mathop{\mathrm{\mathcal{G}}}^*\) satisfies property [prop:cpc] of Lemma 3. Similarly, by definition of \(E_{\text{add}}\), the graph \(\mathop{\mathrm{\mathcal{G}}}' \cup E_{\text{add}}\) is the directed transitive closure of \(\mathop{\mathrm{\mathcal{G}}}'\). This property cannot be destroyed by adding or deleting bidirected edges, which implies that \(\mathop{\mathrm{\mathcal{G}}}^*\) is also transitively closed. Finally, we know that by definition of \(E_{\text{cycles}}\), the graph \(\mathop{\mathrm{\mathcal{G}}}\setminus E_{\text{cycles}}\) is acyclic. This implies that \(\mathop{\mathrm{\mathcal{G}}}'\) is also acyclic. Since the directed transitive closure of an acyclic graph is still acyclic, we know that \(\mathop{\mathrm{\mathcal{G}}}' \cup E_{\text{add}}\) is acyclic. Again, acyclicity cannot be destroyed by adding or removing bidirected edges, which implies that \(\mathop{\mathrm{\mathcal{G}}}^*\) is acyclic.
To show part [part-2] of Lemma 4, note that whenever \(\mathop{\mathrm{incomp}}(\mathop{\mathrm{\mathcal{G}}}) = 0\), i.e. \(\mathop{\mathrm{\mathcal{G}}}\) is graphically compatible, we must have \(E_{\text{cycles}} = E_{\text{add}} = E_{\text{del}} = B_{\text{add}} = B_{\text{del}} = \emptyset\), and hence \(c(\mathop{\mathrm{\mathcal{G}}}) = 0\). The converse follows from part [part-1]. ◻
In the description of Algorithm 8, \(d^{\mathrm{in}}_H(v)\) denotes the in-degree of a vertex \(v\) with respect to the graph \(H\) and similarly, \(d^{\mathrm{out}}_H(v)\) denotes the out-degree. In the description of Algorithm 9, \(tc(H)\) denotes the transitive closure of \(H\) and in the description of Algorithm 10, \(cpc(H)\) denotes the confounding path closure of \(H\) (see Section 3.3).
We begin by stating the details for sampling causal ground truth in our synthetic experiments. For sampling linear bivariate causal statements (for the results in Figures 2 and 3), we implement the following procedure:
Draw random causal coefficients for \(n+m\) variables from a standard normal distribution and error variances for each variable from an exponential distribution with variance \(1\);
Set each causal coefficient to \(0\) independently with probability \(1-p\);
Given the remaining causal coefficients, calculate the covariance matrix of all \(n+m\) variables, assuming independence between the error terms (resulting in a linear SEM without confounding);
Rescale error variances and causal coefficients so that each variable has unit variance;
Select \(m\) variables uniformly at random to be the hidden variables and marginalize the covariance matrix and the causal coefficients (using Lemma 1) to the remaining \(n\) variables.
Obtain the correct linear bivariate causal statements from the equation \(A = (I-\Gamma)^{-1}\) (see Lemma 2).
Perturb each bivariate causal statement by adding centered Gaussian noise with variance \(\sigma\).
Our procedure for graphical bivariate causal statements for the results in Figure 5 is the following:
Draw a uniformly random permutation on \(\{1, \dots, n+m\}\) to get a causal ordering;
Set each directed edge that is consistent with the ordering independently with probability \(p\) to get a causal DAG;
Choose the set of \(m\) hidden variables uniformly at random and marginalize the graph to the remaining \(n\) observed variables according to Definition 6;
Obtain the correct graphical bivariate causal statements by further marginalizing the causal graph to each possible pair of variables.
Add errors by choosing a possible directed or bidirected edge between all variables uniformly at random and adding it or deleting it from the statement graph depending on whether it was absent or present.
The data for the \(7\) country-level indicators that we consider is obtained from the gapminder dataset [25]. We use the following detailed variable descriptions for our experiments:
population density: Average number of people per square kilometer of land in the given country,
literacy rate: adult literacy rate is the percentage of people ages 15 and above who can, with understanding, read and write a short, simple statement on their everyday life,
daily income: mean daily household per capita income or consumption expenditure in constant international dollars,
sanitation access: percentage of people using at least basic sanitation services (improved sanitation facilities not shared with other households),
smoking: percentage of people over age 15 that smoke,
happiness score: national average response to a happiness survey, with scores ranging from 0 (worst) to 100 (best),
life expectancy: average life expectancy at birth in years.
| Population density | 1.000 | 0.109 | 0.708 | 0.104 | -0.018 | 0.078 | 0.128 |
| Literacy rate | 0.109 | 1.000 | 0.373 | 0.798 | 0.109 | 0.526 | 0.716 |
| Daily income | 0.708 | 0.373 | 1.000 | 0.381 | 0.019 | 0.745 | 0.424 |
| Sanitation access | 0.104 | 0.798 | 0.381 | 1.000 | 0.190 | 0.656 | 0.817 |
| Smoking | -0.018 | 0.109 | 0.019 | 0.190 | 1.000 | 0.103 | 0.096 |
| Happiness score | 0.078 | 0.526 | 0.745 | 0.656 | 0.103 | 1.000 | 0.737 |
| Life expectancy | 0.128 | 0.716 | 0.424 | 0.817 | 0.096 | 0.737 | 1.000 |
We show the empirical correlation matrix of these variables in Table 1. We selected these variables as they exhibit a variety of pairwise correlation strengths and causal structure: most notably, population density was chosen as a variable that likely only has weak direct causal relations with many of the other variables, but there is still strong confounding with daily income. On the other hand, there likely are strong causal pathways such as literacy rate \(\to\) daily income \(\to\) happiness score. While we accessed all data through https://www.gapminder.org/data/, the data for some of the variables is originally from different sources, shown in the table below. All data is made available through the CC-BY license.
Table 2 shows the large language models that we query. All models were accessed through Amazon Bedrock and we used the following parameters (except that the Claude Opus models do not allow for specifying top P):
Maximum number of output tokens: 2048
temperature: 0.6
top P: 0.7.
| Model name | Model ID | Provider |
|---|---|---|
| Claude Opus 4.5 | claude-opus-4-5-20251101-v1:0 | Anthropic |
| gpt-oss-120b | gpt-oss-120b-1:0 | OpenAI |
| gpt-oss-20b | gpt-oss-20b-1:0 | OpenAI |
| Gemma 3 27B IT | gemma-3-27b-it | |
| Gemma 3 4B IT | gemma-3-4b-it | |
| Kimi K2 Thinking | kimi-k2-thinking | Moonshot AI |
| Mistral Large 3 | mistral-large-3-675b-instruct | Mistral AI |
| Magistral Small 2509 | magistral-small-2509 | Mistral AI |
| Qwen3 235B A22B 2507 | qwen3-235b-a22b-2507-v1:0 | Qwen |
| Qwen3 Next 80B A3B | qwen3-next-80b-a3b | Qwen |
We use the following conversation to elicit linear bivariate statements from the LLMs that we test.
You are a causality expert, tasked to estimate standardized TOTAL causal effects between country development indicators. Return your answer in HTML format:
<answer>CAUSALCOEFFICIENT: <number ></answer>
For example: <answer>CAUSALCOEFFICIENT: 0.35</answer>or <answer>CAUSALCOEFFICIENT: -0.62</answer>
The causal coefficient quantifies the expected change of the effect variable in standard deviations, given an intervention that changes the cause variable by 1 standard deviation. It includes the effect of all direct causal pathways from the cause to the effect variable. Do not assume away confounding; use realistic domain knowledge. No other text.
I have observational data on 7 country-level development indicators.
Correlation matrix: {corr}
Variable descriptions: {varnames, vardesc}
Before we begin estimating causal coefficients, please determine a plausible causal ordering of these 7 variables (from root causes to downstream effects). Return the ordering as a comma-separated list of variable names inside <ordering>tags. Use the exact variable names shown above.
For example: <ordering>vara, varb, varc, ...</ordering>
<ordering>populationdensity, literacyrate, dailyincome, sanitationaccess, smoking, happinessscore, lifeexpectancy </ordering>
I will now ask you to estimate the total linear causal coefficient for several pairs of variables, following the ordering you provided. For each question, please provide your answer in the format: <answer>CAUSALCOEFFICIENT: <number></answer>
Estimate the total linear causal coefficient for the causal effect of populationdensity on literacyrate.
<answer>CAUSALCOEFFICIENT: 0.1 </answer>
Estimate the total linear causal coefficient for the causal effect of population density on daily income.
<answer>CAUSALCOEFFICIENT: 0.2 </answer>
…
Here, the assistant answers are only examples. Moreover, {corr} is a placeholder for the correlation matrix in Table 1 that we provide as a string, {varnames} is a placeholder for the variable names and {vardesc} is a placeholder for the variable descriptions listed in Section 8.2. If all assistant answers are in the correct format, then we simply keep asking for the causal effects for all pairs of the seven variables in our dataset, according to the provided causal ordering. If a model answer is not in the correct format, then we send one of the following messages:
Please provide the causal ordering as a comma-separated list of the exact variable names inside <ordering>tags. The variable names are: {varnames}
Please provide your answer in the correct format: <answer>CAUSALCOEFFICIENT: <number></answer>
These messages are sent up to 5 times in a row until an answer in the correct format is obtained. In our experiments, this was sufficient to obtain answers in the correct format in most runs.
To obtain lists of graphical bivariate causal statements from LLMs, we use the following template for conversation:
You are a causality expert analyzing relationships between country-level development indicators.
For each pair of variables, you will assess:
1. Whether there is a total causal effect between them, and in which direction
2. Whether there is confounding (correlation not explained by the causal effect between the pair)
Guidelines:
- Answer YES for causal effect if you expect that an intervention on the cause variable would significantly change the effect variable
- Be conservative about confounding - only answer YES for confounding if MOST of the correlation cannot be explained by the causal effect
- Use realistic domain knowledge about socioeconomic factors
Return your answer in the following HTML format:
<causaleffect>
<exists>YES or NO</exists>
<direction>ATOB or BTOA or NONE</direction>
</causaleffect>
<confounding>
<exists>YES or NO</exists>
</confounding>
Where ATOB means the first variable causes the second, and BTOA means the second causes the first.
I have observational data on 7 country-level development indicators.
Correlation matrix: {corr}
Variable descriptions: {varnames, vardesc}
I will now ask you about each pair of variables. For each pair:
1. Assess if there is a causal effect
2. Assess if there is strong confounding
Consider the pair: "Population density" (A) and "Literacy rate" (B).
1. Is there a causal effect between them? If yes, in which direction?
2. Is there confounding (significant correlation not explained by the causal effect between these two)?
Provide your assessment in the required format.
<causaleffect>
<exists>NO</exists>
<direction>NONE</direction>
</causaleffect>
<confounding>
<exists>NO</exists>
</confounding>
Consider the pair: "Daily income" (A) and "Population density" (B).
1. Is there a causal effect between them? If yes, in which direction?
2. Is there confounding (significant correlation not explained by the causal effect between these two)?
Provide your assessment in the required format.
<causaleffect>
<exists>YES </exists>
<direction>BTOA</direction>
</causaleffect>
<confounding>
<exists>YES </exists>
</confounding>
…
If all assistant answers are in the correct format, we keep asking for the causal statements on all pairs of the seven variables in our dataset. However, we present each pair in a random order, since our graphical incompatibility score also tests the extent to which the causal statements are acyclic. If a model answer is not in the correct format, then we send the following message:
Please provide your answer in the correct format: <causaleffect>
<exists>YES or NO</exists>
<direction>ATOB or BTOA or NONE</direction>
</causaleffect>
<confounding>
<exists>YES or NO</exists>
</confounding>
Again, sending this message up to \(5\) times was sufficient to obtain valid model answers.
First, we present additional plots for the difference between compatibility scores based on empirical and true covariance matrices over different combinations of parameters, see Figures 11 and 12. The plots are obtained as follows:
sample \(50\) different linear Gaussian SEMs (see Appendix 8.1),
for each model, sample \(20\) different random perturbations with variance \(\sigma\) for the bivariate causal statements,
for each model-statement combination, draw iid samples \(X \sim \mathcal{N}(0, \Sigma_{\boldsymbol{X}})\) 10 different times and compute compatibility scores.
As the random error in the statements grows, the absolute value of the compatibility scores becomes larger, and more samples might be required to accurately approximate the true compatibility score. However, since only the sign matters for falsification, an accurate approximation is not always needed, and in fact, the signs converge very quickly across all tested combinations of model parameters.
Next, accompanying our results in Section 3.4, we provide a breakdown of the average incompatibility scores into the different components calculated by our heuristic algorithm (see Equation 5 ) in Figure 13. Notably, Gemma 3 27B IT is the only model that does not have inconsistencies with respect to its bidirected edges. However, a closer look (Figure 14) reveals that it is also the model that reports complete or almost complete graphs in all runs. Indeed, the third property of Lemma 3 is trivially satisfied when there is a bidirected edge between every single pair of variables. We tried to avoid this situation by formulating the prompt for graphical causal statements in a way that encourages models to be conservative with placing bidirected edges. However, there is still a large variance in the density of the statement graphs for different models. When focusing on the subset of statement graphs with density at most \(2/3\), the lowest achieved incompatibility score is \(4\) for a statement graph given by Claude Opus 4.6. For reference, we show this graph in Figure 15, but we remark that a low incompatibility score does not imply correctness of the graph.
The independence of causal mechanisms assumption states that the conditional distributions of each variable given its parents should not be informative of each other [23]. In our setting, this requires the rows of \(\Gamma\) and the noise covariance matrix to be chosen independently. Our assumption is therefore a slight strengthening (see [24] for a similar proposal).↩︎
For this reason, the quality of a set of causal statements should never be judged based on compatibility alone, but also account for informativeness and falsifiability [33]. Clearly, stating an unconfounded causal link is a more specific statement - and thus more falsifiable - than allowing confounding in addition.↩︎