The variance of causal effect estimators for binary v-structures

1 Introduction

As graphical representations of multivariate probability distributions, Bayesian networks are widely used statistical models with an underlying directed acyclic graph (DAG) structure. When taking DAGs to represent causal diagrams [1,2, 3,4], we may use a machinery based on the “do” calculus [5] to estimate potential intervention effects of any variable on any other. Different graphical criteria exist to identify valid adjustment sets, among which the back-door criterion [6] is probably the most well known, and with more generalised strategies developed more recently [7,8].

A valid adjustment set Z for the effect of X on Y is such that for any probability distribution p compatible with the underlying graphical structure, the probability distribution of Y after intervening on X (setting it to some value) satisfies [9]

For linear Gaussian models, the marginalisation can be simply estimated by regressing Y on X and Z and extracting the coefficient of X , hence the naming of “adjustment” sets. This also holds for linear non-Gaussian causal models [10].

The set of parents of X always satisfies the back-door criterion and is therefore a valid adjustment set, but there may be many more depending on the graphical structure of the DAG [8]. Although all valid adjustment sets provide consistent estimators of the causal effects, for finite-sized data, different adjustment sets can lead to different numerical estimates, and with different precisions.

In evaluating the variance of different estimators, the remarkable result that the asymptotically optimal adjustment set can be determined solely based on graphical criteria regardless of the edge coefficients has recently been obtained [10]. Even more recently, this has been extended to non-parametric estimators [11], and the asymptotically optimal set has been further characterised [12].

To explore the precision of causal estimators for non-linear models, we consider the simplest such case: a DAG with three nodes of binary variables organised in a v-structure with the outcome Y of interest as a collider with parents Z and X (Figure 1) and with the latter being the exposure whose effect we wish to estimate. For binary data and relatively small networks, one can explicitly marginalise over the remaining nodes in the DAG and its parameters [13] to derive interventional distributions as follows:

and estimate causal effects from them.

Figure 1

A v-structure on three nodes.

In the simple case of a v-structure (as in Figure 1), there is no confounding of the effect of X on Y (there are no common parents) so that the empty set constitutes a valid (and minimal) adjustment set, and the interventional distribution is simply

A valid expression for computing the total causal effect of X on Y , in accordance with equation (1), is then

where we used the subscript R for raw, to highlight the fact that the formula only involves raw (or observed) conditional probabilities of Y on X , which in this simple scenario are sufficient to identify the desired causal effects.

However, by definition, the conditional distribution of Y on X is also

Therefore, another valid expression for the total causal effect of X on Y is

where we used the subscript M to highlight the fact that the formula derives from explicitly marginalising Z out from the joint distribution p ( Y , Z ∣ X ) . In contrast, one could interpret the formula based on raw conditionals as performing the marginalisation implicitly (with the observations already providing a marginalised sample).

A more general way of understanding equation (6) is by observing that in the case of the v-structure Z also constitutes a valid adjustment set (albeit not a minimal one). Then starting from the joint interventional distribution p ( Y , Z ∣ do ( X ) ) , the interventional distribution of Y when intervening on X is also

with the latter equality justified by structural and invariance properties and also in agreement with the standardisation formula in equation (1).

Since we see that adjustment by Z is valid, but not necessary, it is natural to ask whether the two estimators differ in terms of precision. To answer the question we compute the variance, for finite sample sizes, of the two different estimators corresponding to the implicit and explicit marginalisation as outlined earlier.

It is instructive to also consider the DAG with the edge from Z → Y deleted. Since Y would then be independent of Z , the marginalisation would reduce to the raw conditionals. The estimator using raw conditionals is therefore the same whether the edge from Z → Y is present, while the approach using marginalisation would give different estimates for the two cases. Intuitively we would expect that the extent by which estimates differ will depend on the strength of the relationship between Y and Z . The underlying rationale is similar to that for the standard practice of adjusting for baseline covariates in models of the outcome in randomised controlled trials [14], where the prognostic factors Z and the (randomised) treatment X can be seen as forming a v-structure with the outcome Y as the collider.

The questions of whether adjusting for baseline covariates is justified or even desirable is indeed a recurrent one with a long history in the context of randomised controlled trials [15,16, 17,18,19, 20,21]. An extensive body of literature exists with several proposals to exploit covariate adjustment in order to build more efficient estimators [22,23,24, 25,26,27]. Health authority guidelines also typically recommend adjustment for the sake of improving precision [28,29,30]; however, in the case of non-linear models, special care must be taken to account for the potential non-collapsibility of effect measures [31,32,33] and that the estimand may change depending on the method used for adjustment [34]. Recent work [26] examines the performance of covariate-adjusted estimators in clinical trials with binary, ordinal, or time-to-event analysis. For binary outcomes, they present simulation results where the baseline covariate represents categories of age and the estimand is the risk difference, for a magnitude of the prognostic value of the covariate mimicking the association found in observational data. Our study considers the question from a theoretical and causal diagram perspective, in a slightly simplified scenario where the treatment, outcome, and prognostic factors are all binary variables, but where we evaluate how things change with the strength of the prognostic value of the covariate.

Therefore, we consider the v-structure since it provides the simplest example of a causal diagram where there is a choice between different adjustment sets. If we add an edge in the graph of Figure 1 connecting X and Z we end up with no choice about adjustment sets: in particular, if we add an edge from X → Z , then Z is not a valid adjustment set and the empty set is the only choice; conversely, if we add an edge from Z → X , then Z is a confounder (a common parent), and it must be adjusted for, making it the only valid adjustment set with the empty set no longer valid.

2 Causal estimates for a binary v-structure

For both causal estimators, we will use the maximum likelihood estimates of probabilities from the observed data. We consider the DAG in Figure 1 with the following probability tables:

When we generate data, as a collection of N binary vectors, from the DAG in Figure 1, instead of forward sampling along the topological order for this small example, we can sample directly from a multinomial with probabilities

(9) If we represent with N i the number of sampled binary vectors indexed by i = 4 X + 2 Z + Y , then the estimator of F from the raw conditionals is simply

By using the marginalisation, we would have the following estimator:

with the terms separated for later ease

These estimators, as they rely on observed data frequencies, are non-parametric and fit in the recent general framework for arbitrary graphs [11]. The key advance of our derivation with respect to their result is that we consider terms beyond the leading-order asymptotics and compute the variance of the estimators for arbitrary sample sizes, which further enables us to perform more detailed asymptotic analyses.

2.4 Relative difference in variances

To explore the difference in variances, we focus mainly on the case where the effect of Z on Y is the same for each X to remove one parameter degree of freedom from the seven in total. In this case, we write the probabilities as follows:

(29) p Y , 0 = q 0 − C , p Y , 1 = q 0 + C , p Y , 2 = q 1 − C , p Y , 3 = q 1 + C ,

where C is a measure of the effect of Z on Y (the same for each X ) and the causal effect of X on Y is q 1 − q 0 .

We plot the relative difference in variances of the two estimators, Δ = V [ M ] − V [ R ] V [ R ] . In Figure 2, we leave p X free, set p Z = 2 3 and set q 0 = 1 3 , q 1 = 2 3 and plot Δ for N = 100 and N = 400 . In the plot for N = 400 , we also scaled C by dividing by 2. The behaviour and rescaled plots are very similar, suggesting a N − 1 2 scaling.

Figure 2

The relative difference in variance of the two estimators for two sample sizes. (a) N = 100 and (b) N = 400 .

With our general result, we can further allow for an interaction between Z and X and have a different effect of Z on Y for each X . To account for the interaction strength, we introduce parameter D , which we multiply by p Z or ( 1 − p Z ) in the following parameterisation:

(30) p Y , 0 = q 0 − C − D p Z , p Y , 1 = q 0 + C + D ( 1 − p Z ) , p Y , 2 = q 1 − C + D p Z , p Y , 3 = q 1 + C − D ( 1 − p Z ) ,

so that the causal effect of X on Y remains unchanged as q 1 − q 0 . The effect of Z on Y in the probability tables is then 2 C + D for X = 0 and 2 C − D for X = 1 , making D a measure of the interaction effect between Z and X . In Figure 3, we again leave p X and C free, set p Z = 2 3 , q 0 = 1 3 , q 1 = 2 3 , N = 100 and plot Δ for D = 1 8 and D = 1 4 . We observe, as in Figure 2, a central region where R is the more efficient estimator and tails for larger C , where M is the better estimator, but now with a rotation dependent on D and p X .

Figure 3

The relative difference in variance of the two estimators for two interaction strengths. (a) D = 1 8 and (b) D = 1 4 .

3 Asymptotic behaviour

To examine the asymptotic behaviour of the causal effect estimators in more detail, we return here to the setting of equation (29) with the same effect of Z on Y for each X (or with D = 0 ), and we expand the hypergeometric function as in Appendix C. We treat the general case with interactions ( D ≠ 0 ) in Appendix D, while without interactions ( D = 0 ) we obtain the following formula for the variance of R :

and for the variance of M ,

To extract the asymptotic behaviour of the difference in variances of the two estimators, we consider C ∼ N − 1 2 to obtain

with root

so that

Note that although we used the scaling C ∼ N − 1 2 to extract this result, it holds more generally. For example, for fixed C ≠ 0 , it is trivial to see that C > C ∗ for some N and so that V [ M ] will become lower than V [ R ] in the limit N → ∞ . The asymptotically optimal adjustment set therefore uses marginalisation rather than raw conditioning, in line with previous results [10,11] from the leading-order asymptotics. For fixed C = 0 , however, raw conditioning would be better. It is exactly by treating subleading terms, as we do here, that we can examine where the transition occurs and how it depends on the coefficients. For weaker effects of the edge from Z → Y , with C ≲ N − 1 2 , the raw conditional can give a more precise estimate of the causal effect of X on Y .

The result in equation (35) therefore shows that the optimal adjustment set, for the effect of X on Y , does not depend solely on graphical criteria, but crucially on the relative effect size of the influence of the other node Z on the outcome Y . The result details, for large but finite sample sizes, the crossover in efficiency from R to M as the relative effect of C increases. Therefore, starting from plausible parameter values, we can use equations (34) and (35) to guide our choice of adjustment.

Since for many practical purposes the sample size may be determined by other considerations and we can never take the limit N → ∞ , the asymptotic regime developed here accounting for the relative scale of effects compared to the sample size is most relevant. Although derived just for the binary v-structure, this is a counterexample showing that recent leading-order aysmptotic results [11] cannot directly extend outside their particular asymptotic limit.

4 Implications for causal discovery

With a larger sample size, we may be able to detect and quantify smaller causal effects. Therefore, we wish to get a feeling for the strength of the edge Z → Y we would detect from the data, or equivalently for which values of C we would infer the presence of the edge. To do so, we calculate the expected difference in maximised log-likelihoods when including the edge compared to a DAG with the edge deleted:

where the 1 2 comes from Wilk’s theorem [36] for the additional parameter when maximising all the probabilities relative to evaluating with the restriction C = 0 .

The change is AIC is then

There is therefore an asymptotic regime where the edge is strong enough to detect on average using the AIC, but the estimator from raw conditionals that does not use the edge has lower variance:

which follows from the Cauchy–Schwarz inequality. The regime only vanishes when p Z = 1 2 and q 1 ( 1 − q 1 ) ( 1 − p X ) 2 = q 0 ( 1 − q 0 ) p X 2 , and the two bounds become equal. Utilising the BIC instead ( E [ Δ BIC ] = E [ Δ AIC ] + log ( N ) ) leads to a large regime where we would not detect the edge on average, but where the estimator using marginalisation that does rely on the edge has lower variance.

5 Discussion

To evaluate the precision of different estimators targeting the same causal effect in causal diagrams, we considered the simple case of a v-structure for binary data and explicitly computed the variance of the two different estimators for the effect of X on the collider Y , with Z as the other parent.

The results involve combinations of hypergeometric functions, suggesting that exact results for larger DAGs may be rather complex. Which estimator has the lower variance depends, among other parameters, on the relative strength of the edge from Z to Y . In general, estimating the causal effect through marginalisation offers better performance in the presence of a stronger direct effect of Z on Y . When the direct effect is weaker instead, ignoring the edge and estimating the causal effect through the raw conditionals provides higher precision.

In light of our results, it is instructive to go back to the parallel with clinical trials. Even if the simple case of three nodes with only binary variables does not cover the more general and realistic scenarios encountered in randomised clinical trials, the example is still enlightening to show that taking a gain in precision for granted under any circumstances may not always be justified. Asymptotic results suggest that the optimal adjustment set in the sense of efficiency should include the prognostic factor Z . In real-world situations, with finite sample sizes, whether adjusting actually improves efficiency may depend on the relative strength of the prognostic value of Z on the outcome Y with respect to the sample size. This would seem to match the intuition that including stronger prognostic factors will indeed benefit precision, while it may be counterproductive to adjust for weaker ones.

By examining the asymptotic regime of large sample sizes, we could confirm the intuition that for edge strengths statistically detectable by the AIC, accounting for the edge in the estimation should generally lead to lower variance. Conversely, that the presence of statistically non-detectable edges should be ignored to achieve a lower variance.

Most importantly, we could also discover an asymptotic regime where raw conditional estimates, ignoring the edge, were more precise in the presence of statistically detectable edges. One way to appreciate the practical relevance of these findings is by observing that we can expect ranges of causal strengths, which become statistically detectable from data before we can gain precision by accounting for them in the estimation. Our detailed asymptotic analysis for the v-structure goes beyond the leading-order asymptotic result where the optimal estimator does not depend on the edge coefficients [10,11].

Outside the asymptotic regime, for finite sample sizes, the gain in precision when using marginalisation and, thus explicitly accounting for the edge presence, appears to be linked to its strength. Although the example considered here is the simplest non-trivial DAG, this finding further supports the idea that learning the full structure of the graph, beyond simply identifying a valid adjustment set, may benefit the precision of causal inference. The practical limitation with observational data is that we can only learn structures up to an equivalence class, so that we need to consider the possible range of causal effects across the whole class [37] or implement Bayesian model averaging across DAGs [13].

If we use a more stringent criterion to decide about the presence of edges, such as the BIC, for example, which implements a stronger penalisation with respect to the AIC, we may end up missing edges too weak to detect on average, but whose presence would improve the precision of the causal estimation through marginalisation. In other words, for moderately weak direct effects, the selection of suitable adjustment sets may be relatively sensitive to the choice of the score. Analogously, we may expect that optimal causal estimation may also be sensitive to the choice of learning algorithm, whether constraint-based search [38,39], score-based search [40] or Bayesian sampling [41,42,43]. Quantifying the extent by which the structure learning affects causal estimation constitutes an interesting line of further investigation.

Finally, a limitation of the current study is its relatively reduced scope, treating a very simple situation with only three binary variables. It would be interesting to see if and how results generalise in the presence of multiple and possibly non-binary covariates as well as to non-binary outcomes. Unfortunately, extensions of our approach to higher dimensions are technically challenging, especially because the number of parameters grows exponentially for larger graphs with more edges. Similarly, extending to non-binary variables requires re-framing the problem in a relatively different setting. Nevertheless, we feel that even this simple example is of practical value to highlight that paying attention to the relative strength of the prognostic value with respect to the sample size is critical to determining the extent of gain in precision. This is echoed in simulations [26], with the discussion of ref. [44] further elaborating on the relative importance of the strength of the covariate associations with the outcomes and the sample size. The conclusions we draw from our exact calculations for a non-parametric estimator of the risk difference appear to match finite sample effects in simulation analyses [44] very closely. Furthermore, it seems natural to expect similar considerations applying to different estimands, such as the risk ratio and the odds ratio, though technically more challenging to treat. Notably, simulation studies provide for a valuable alternative to gain insights about more realistic scenarios, which are technically intractable, as nicely illustrated in ref. [26].