Empirical Bayesian Approaches for Robust Constraint-based Causal Discovery under Insufficient Data

1 Introduction

Learning causal relations has been a fundamental and widelyinvestigated topic. The causal relations are captured by a directed acyclic graph (DAG), and a directed link in DAG captures cause-effect relation between two variables connected by the link. Specifically, a directed link from node X to node Y indicates the cause-effect relation between cause variable X and effect variable Y . Causal discovery aims at learning a DAG capturing causal-effect relationships among a set of random variables from observational data. Existing causal discovery methods focus on learning a DAG with high confidence from sufficient data samples. Not much attention, however, has been paid to performance improvement of causal discovery under limited data. Such work is important, as even in the era of big data, there are still domains in which the availability of data is very limited. For example, in biological or clinical disciplines, data can be severely insufficient either because of high cost or lack of cases from which data is collected [Mukherjee and Speed, 2007]. Furthermore, even for applications with a vast amount of data, the data may not adequately cover all possible states of the nodes, leading to insufficient data for certain states. For example, the observed data under the absence of earthquake is adequate, while the observed data under the occurrence of earthquake is limited, due to the fact that earthquake rarely happens in nature.

∗Corresponding author

Constraint-based causal discovery methods apply independence tests to determine a DAG from observational data and can be performed globally or locally. Global approaches aim at learning cause-effect relationships among all random variables, such as PC-stable [Colombo and Maathuis, 2014], and Sepset consistent PC (SC-PC) [Li et al., 2019]. Global causal discovery methods discussed above learn DAGs that are in the same markov equivalent class of ground truth DAG. Further tests under certain assumptions about the graph or data distribution are needed to resolve the causal ambiguity [Glymour et al., 2019]. In this paper, we focus on learning markov equivalent DAGs. In contrast to global approaches, local approaches identify the direct causes and effects of a target variable, represented by a causal Markov Blanket [Gao and Ji, 2015; Yang et al., 2021]. A causal Markov Blanket captures local relationships of a target variable by identifying its parents, children, and spouses. For both global and local approaches, the main challenge of constraint-based causal discovery methods is that their performance highly depends on the accuracy of the independence test. Independence test error, even one mistake in independence decision, can propagate throughout the graph, causing a sequence of errors and resulting in an erroneous DAG with incorrect orientations [Spirtes, 2010]. Hence, to perform a robust constraintbased causal discovery, it is crucial to improve the robustness of the independence test.

To improve the causal discovery performance under insufficient data, we propose to introduce Bayesian approaches to independence tests for accurate and efficient constraint-based causal discovery. Specifically, two Bayesian-augmented frequentist independence tests are proposed, whereby we use Bayesian approach to reliably estimate, under low data regime, independence test statistics used by frequentist independence tests. For both MI estimation (Sec.3.1) and hypothesis likelihood estimation (Sec.3.2), we employ Bayesian inference to calculate statistics by considering the entire parameter space instead of using a point estimate one. Given the estimated Bayesian statistics, we follow the standard frequentist framework to perform independence test. The proposed Bayesian-augmented independence tests are then applied to improve the constraint-based causal structure learning. We evaluate both local and global causal discovery performance with proposed independence tests on benchmark datasets and compare them to state-of-the-art methods. We empirically demonstrate the effectiveness of the proposed Bayesian approaches in improving both the accuracy and efficiency of the local and global causal discovery under insufficient data.

2 RelatedWork

To handle causal discovery under insufficient data, some methods downsize the problem domain to sub-domains. Rohekar et al., [2020] approximated the structure by performing independence tests with a small fixed size of the condition set. The structure was then refined by iteratively increasing the condition set. A similar idea was explored in the Recursive Autonomy Identification (RAI) method [Yehezkel and Lerner, 2009]. Related works along this line always assume that there exist sufficient data for sub-domains. Besides, Claassen and Heskes [2012] estimated the posterior distribution of the independence hypothesis between two variables, based on which reliability was quantified. The causal discovery was then processed in decreasing order of reliability. Rohekar et al., [2018] estimated posterior distribution of DAG through bootstrap samples. The negative effect from independence tests error was minimized through model averaging.

Some causal discovery methods address the limited data issue by directly improving the independence test [Marx and Vreeken, 2018]. A Bayesian-augmented frequntist independence test based on Bayes Factor (BF) was proposed [Natori et al., 2017] whereby Bayesian parameter estimate is employed in computing BF while the value of BF is then applied to a frequentist independence test. The proposed independence test is incorporated into RAI, achieving competitive DAG learning performance. However, a threshold is required in [Natori et al., 2017] and the selection of threshold can be heuristic. Instead, we propose to formulate the Bayes Factor into a well-defined statistical independence test without requiring threshold tuning.

In addition, different approaches have been proposed for robust independence tests under insufficient data. These methods, however, are not aimed at improving the causal discovery performance. Seok and Seon Kang [2015] improved the estimation of mutual information (MI) by partitioning the whole sample space into sub-regions. For better MI estimation under limited data, Bayesian approaches have been widely considered [Hutter, 2002; Archer et al., 2013]. Another category of recent independence test techniques are focused on developing non-parametric methods to improve efficiency, such as CCIT [Sen et al., 2017] and RCIT [Strobl et al., 2019]. These works assume the availability of sufficient data and are mainly focused on continuous variables, while we are focused on discrete ones.

3 ProposedMethods

We consider two types of independence test: MI based and statistical testing based independence tests. We introduce Bayesian approaches to improve both types of independence tests through a Bayesian-augmented frequentist framework. Particularly, for MI based approach, we employ empirical Bayesian approach for better MI estimation under limited data. For statistical testing based approach, we consider the empirical Bayesian estimation of hypothesis likelihood and formulate it into χ2 statistical independence test, providing an accurate p-value under limited data.

3.1 Bayesian Approach for Mutual Information Based Independence Test

The mutual information (MI) of two discrete random variables X and Y is defined as MI(X; Y ) = ∑Kx i=1 ∑Ky j=1 P(xi, yj) log P(xi,yj) P(xi)P(yj) , where Kx and Ky denote the total number of possible states of X and Y respectively. P(xi, yj), P(xi), and P(yj) represent the joint probability of (X, Y ), and the marginal probabilities of X and Y respectively. By definition, MI(X; Y ) = 0 if and only if X and Y are independent. In practice, the true MI is unknown, and the estimated MI is always larger than zero. In the following, we denote the probability distribution parameters as θ, i.e., P(xi) = θi, P(yj) = θj and P(xi, yj) = θij. Conventionally, MLE is employed to estimate θ from data as ˆ θ = argmaxθ p(D|θ), where P(D|θ) is the likelihood of parameter θ given the data D. MI is then estimated as MI = MI(X; Y | ˆ θ). When data is insufficient, MLE is not reliable [Geweke and Singleton, 1980] and MI tends to be overestimated. Instead, the full Bayesian MI is estimated from data over the entire parameter and hyper-parameter space, i.e.,

where α is the hyper-parameter for symmetric Dirichlet prior of θ. The full Bayesian MI is the expected MI over the joint posterior distribution of the parameters and hyper-parameter, i.e., p(θ, α|D). The integration over hyper-parameter α can be computationally challenging [Archer et al., 2013]. Instead of marginalizing out α, we propose to maximize it out. Particularly, we approximate the integration over the hyperparameter space by its mode α∗ that maximizes a posterior (MAP) of α, i.e., α∗ = argmaxα p(α|D). By assuming uniform distribution of p(α), we have α∗ = argmaxα p(D|α). The likelihood p(D|α) can be computed as,

where K is the number of states for the random variable, ni is the number of samples for state i, and N = ∑Ki ni. P(D|α)

follows Polya distribution and Γ(x) is the gamma function. We solve for α∗ with a fixed-point update [Minka, 2000]. Given α∗ , the full Bayesian method is converted to the empirical Bayesian method, and we have the proposed empirical Bayesian MI ˆ MI eB defined as,

with a closed-form solution:

where ψ(x) is the digamma function. ni and nj are the number of samples for X = i and Y = j respectively, and nij is the number of samples for (X, Y ) = (i, j). Given the estimated MI, we compare it against a pre-defined threshold for independence test. If MI is smaller than the threshold, two random variables will be declared to be independent, and dependent otherwise.

3.2 Bayesian Approach for Statistical Testing Based Independence Test

We now introduce our proposed Bayesian approach to improve the statistical testing based independence test. We firstly consider a standard independence test, G test [McDonald, 2009], which is a likelihood ratio test with null hypothesis assuming two random variables are independent. G test is a widely used statistical test. As the same with other statistical tests, G test doesn't require threshold tuning and the significance level is set to be 5% by default. The formula for the statistic G reads as G = -2 ∑Kx i=1 ∑Ky j=1 nij ln ˆ θi ˆ θj ˆ θij where ˆ θ = argmaxθ P(D|θ). Samples D = {Dn}Nn=1 are i.i.d given parameter ˆ θ and the statistic G follows asymptotic χ2df=(Kx-1)(Ky-1) distribution, based on which a statistical test can be performed. As MLE parameter estimates are not reliable under insufficient data, leading to inaccurate estimation of the likelihood of hypothesis, we instead consider the empirical Bayesian estimation. Specifically, we employ the Bayes Factor (BF) [Kass and Raftery, 1995] which defines the ratio of expected likelihoods of null hypothesis(H0) and that of the alternative hypothesis(H1) over all possible parameter settings with the posterior distributions of parameters under null and alternate hypothesis respectively,

where α0 and α1 are the hyper-parameters for the symmetric Dirichlet prior under null and alternate hypothesis respectively. Both hypothesis likelihoods P(D|H0, α0 ) and P (D|H1, α1 ) can be analytically solved, and BF can be computed. However, BF can't be directly applied to a statistical test because samples D = {Dn}Nn=1 are not i.i.d given hyper-parameter α and BF no longer follows the χ2 distribution under the null hypothesis. Instead, we propose to approximate P(D|α) by a multinomial distribution and calculate modified parameters of multinomial distribution with α taken into account, as both capture the distributions for integer random variables, i.e.,

where K is the total number of states, and N = ∑Ki=1 ni is the total number samples with ni being the number of samples for state i. ˜ θi are the modified parameters of the multinomial distribution. ˜ θi = g(ni,α) g(N,Kα) with g(ni, α) = ani + bα

P(D|α) (defined in Eq. 2) into Eq. 6, it is clear that to satisfy Eq. 6, we must have ni ln g(ni, α) = ln Γ(ni +α) - ln Γ(α). Given {ni}Ki=1 and α, we can construct a system of K such equations through which we can solve for Λ∗ , i.e.,

exp( 1 ni (lnΓ(ni + α) - lnΓ(α))). Given Λ∗ , we have ˜ θi as ˜ θi = g(ni,α) g(N,Kα) = a∗ ni+b∗α a∗N+b∗Kα . P(D| ˜ θ) can well approximate P(D|α). Our proposed approximation is different from the method provided in [Minka, 2000] where the Polya distribution P(D|α) is interpreted as a multinomial distribution with modified counts ˜ni. In addition, our proposed estimation can better approximate the Polya distribution given the symmetric Dirichlet prior compared to [Minka, 2000]. We then approximate the hypothesis likelihood under null and alternative hypothesis respectively and obtain a modified Bayes Factor ˜ BF

We obtain the statistic BFchi2 for the statistical test as,

The statistic BFchi2 asymptomatically follows the χ2df=(Kx-1)(Ky-1) distribution. If p-value is smaller than the significance level, we reject the null hypothesis and accept the alternative hypothesis. It is worth noting that BF can be directly applied for a frequentist independence test where a pre-defined threshold η is required [Natori et al., 2017]. The value of the threshold is unconstrained [Kass and Raftery, 1995] making it hard to be properly selected. Instead, our approach only requires a significant level for independence test which is usually set to be 5% by default.

4 Experiments

We evaluate both the local and global constraint-based causal discovery performance on benchmark datasets. Through exhaustive experiments, we show that our approaches can significantly improve causal discovery performance in terms of both accuracy and efficiency over state-of-the-art methods. Besides, we compare proposed independence tests to stateof-the-art independence tests to further show the effectiveness of the proposed methods.

Experiment Settings. We employ six benchmark datasets 1 that are widely used for causal discovery evaluation: CHILD, INSURANCE, ALARM, HAILFINDER, CHILD3 and CHILD5. The causal discovery performance is evaluated in terms of both accuracy and efficiency. For accuracy, we employ the structural hamming distance (SHD) [Tsamardinos et al., 2006a]. SHD computes the number of extra and incorrect (missing and reverse) edges in the learned causal structure compared to the ground truth one. For efficiency, we consider the number of conducted independence test. We perform evaluation on a number of small sized datasets. These small sample sizes are chosen to mimic insufficient data scenario through significantly small number of samples per configuration. For each sample size, we repeat 10 runs and report the averaged performance over 10 runs. In addition, we report standard derivation of SHD. Experiments are performed on a laptop with a 8-Core Intel Core i9 processor with CPU only.

1 https://www.bnlearn.com/bnrepository/.

4.1 Local Constraint-based Causal Discovery

For the local causal discovery, we employ Causal Markov Blanket (CMB) [Gao and Ji, 2015], which is the state-ofthe-art method. CMB employs constraint-based approach and performs conditional independence test using MI to identify the CMB of a target node. We incorporate the proposed independence tests into CMB and compare to the original CMB. We denote cI eB as the CMB with empirical Bayesian MI estimation and cBFchi2 as CMB with BFchi2 independence test. SHD is 0 if learned CMB is identical to the ground truth CMB.

From Table 1, we can see that both cI eB and cBFchi2 outperform the CMB on all datasets in terms of both accuracy and efficiency under insufficient data. The number of performed independence test reduces dramatically. On ALARM dataset, cIeB only performs 2210 independence tests on average, while CMB requires 17907 tests on average. The proposed methods improve the accuracy significantly. On INSURANCE dataset, cBFchi2 improves the averaged SHD by 3.92 compared to CMB. From the results we can see that, by introducing Bayesian approaches, both the accuracy and the efficiency can be improved. Comparing the performance between the two proposed methods, cBFchi2 achieves overall better accuracy, and cI eB is more efficient with the fewest number of independence test on all datasets.

It is worth noting that the number of independence test increases with reduced samples in CMB, but decreases with the proposed methods. The reason is that under insufficient data, MLE will lead to an overestimated MI. Hence, conventional MI based independence test is likely to declare dependence when data size is small, resulting in a large number of independence test. As the sample size increases, the incorrect dependency declarations will be corrected and the number of independence tests will decrease. On the other hand, our methods are more accurate and show a preference of independence under insufficient data, resulting a small number of performed independence test.

4.2 Global Constraint-based Causal Discovery

Majority of global causal discovery algorithms are under causal sufficiency assumption, whereby all random variables are observed in data and there is no latent variable. However, causal sufficiency assumption can be violated since the real data may fail to capture the values for all the variables, leaving some variables to be latent. To address this issue, several recent causal discovery methods [Ramsey et al., 2012; Colombo et al., 2012] have been developed to identify latent common confounders of the observed variables. In our evaluations, we mainly focus on standard algorithms that are under causal sufficiency assumptions. We firstly employ RAI [Yehezkel and Lerner, 2009] as our baseline and compare to two state-of-the-art methods. Then, to demonstrate that our proposed methods can consistently improve causal discovery performance, we consider wellknown DAG learning algorithms: PC [Spirtes et al., 2000] and MMHC [Tsamardinos et al., 2006b] as two additional baselines. In the end, we consider the algorithms without causal sufficiency assumption to demonstrate that our proposed methods can be applied to different causal discovery methods, independent of the existence of latent confounders.

Global causal discovery with causal sufficiency assumption. We employ RAI as our baseline algorithm and incorporate the proposed independence tests. We denote rI eB as the RAI with empirical Bayesian MI estimation, and rBFchi2 as RAI with BFchi2 independence test. We compare our approaches to two state-of-the-art methods: RAIBF method [Natori et al., 2017] and PC-stable [Colombo and Maathuis, 2014]. SC-PC can't be performed under insufficient data smoothly, and thus we exclude this method for comparison. SHD is 0 if the learned DAG and the ground truth DAG belong to the same equivalence class.

From Table 2, we can see that rBFchi2 outperforms RAIBF and PC-stable on almost all datasets in terms of both accuracy and efficiency. rI eB also achieves overall better accuracy and significantly improves efficiency. For example, on CHILD3, rBFchi2 improves the SHD by 21.7 and 27.6 compared to RAI-BF and PC-stable. In terms of efficiency, on HAILFINDER, rI eB only performs 2215 independence tests in average, while RAI-BF requires 10212 tests in average. Comparing between the two proposed methods, rBFchi2 achieves better accuracy and rI eB achieves better efficiency. With the proposed methods, the number of independence tests decreases due to the reduced samples for all datasets, which is consistent with the conclusion we have from the local causal discovery. In addition, both RAI-BF and PC-stable show a preference of independence under insufficient data, leading to the decreased number of independence tests with reduced number of samples.

Since BFchi2 essentially is only an approximate of original BF, BF with the optimal threshold should outperform BFchi2

in principle. However, selecting the optimal threshold for BF can be challenging and incorrect thresholds can lead to inferior causal discovery performance. Instead of fixing the threshold of RAI-BF with its default value, we consider the optimal performance of RAI-BF with tuned thresholds for comparison. According to the experimental results, RAI-BF with the optimally tuned threshold at best achieves comparable performance compared to rBFchi2 in terms of both accuracy and efficiency, which is expected. While rBFchi2 only requires a fixed significance level (5% by default) without additional tuning process.

To further show that our proposed methods can consistently improve the causal discovery performance, we consider another two widely used DAG learning algorithms: PC [Spirtes et al., 2000] and MMHC [Tsamardinos et al., 2006b]. We incorporate the proposed methods into PC and MMHC for evaluation. According to the experimental results, our proposed methods can consistently improve the DAG learning performance, particularly with PC. For example, on ALARM, PC with BFchi2 achieves averaged SHD 40.5, while PC only achieves averaged SHD 58.2. Overall, BFchi2 achieves better accuracy and I eB achieves better efficiency with both PC and MMHC on different datasets.

Global causal discovery without causal sufficiency assumption. To demonstrate that our robust independent tests can also be applied to causal discovery without causal sufficiency assumption, we employ the conservative FCI (cFCI) method [Ramsey et al., 2012] as our baseline. cFCI is considered as the state-of-the-art causal discovery algorithm that identifies latent confounders. We denote cI eB as the cFCI

with empirical Bayesian MI estimation, and cBFchi2 as cFCI with BFchi2 independence test. We compare our approaches to cFCI with default g2 statistical based independence test 2 . As we can see from Table 3, cI eB achieves best efficiency by performing the smallest number of independence tests. In terms of accuracy, cBFchi2 achieves overall better performance. The consistent performance improvement further demonstrates that the proposed independence test can improve the causal discovery performance under insufficient data, independent of the existence of latent confounders.

4.3 Bayesian Approaches for Independence Tests

To compare the proposed independence tests to state-of-theart methods, we firstly perform a direct evaluation of proposed independence tests on synthetic data, and we then compare to state-of-the-art methods in terms of causal discovery performance on benchmark datasets. On synthetic data, we compare to three state-of-the-art independence tests: adaptive partition [Seok and Seon Kang, 2015], empirical Bayesian with fixed α [Hutter, 2002] and full Bayesian method [Archer et al., 2013]. We evaluate the performance in terms of both accuracy and efficiency. Experimental results show that the proposed methods achieve better accuracy with significantly improved efficiency. More importantly, we compare proposed independence tests to two state-of-the-art methods: adaptive partition and empirical Bayesian with fixed α methods in terms of causal discovery performance on benchmark datasets. Because the full Bayesian method is of high computational complexity, making it impractical to be applied to constraint-based causal discovery, we exclude the comparison to this method. We incorporate the adaptive partition method and the empirical Bayesian with fixed α method to RAI (denoted as rI AdP and rI eBFix respectively). As shown in Table 4, our methods achieve overall better accuracy than rI AdP and rI eBFix on different datasets. For example, on CHILD3, rBFchi2 achieves averaged SHD 53.2, significantly better than rI AdP which achieves averaged SHD 81.6. In terms of efficiency evaluation, rI eB also achieves competitive efficiency. Overall, our proposed methods outperform other SoTA independence tests in terms of causal discovery performance. On HAILFINDER, because rI AdP tends to declare independence, the learned DAG contains fewer false positive edges compared to other methods and thus its averaged SHD is the best.

2 https://github.com/striantafillou/causal-graphs.

5 Conclusion

In this paper, we introduce Bayesian methods for robust constraint-based causal discovery under insufficient data. Two Bayesian-augmented frequentist independence tests are proposed for reliable statistic estimation under a frequentist independence test framework. Specifically, we propose: 1) an effective empirical Bayesian method for accurate estimation of mutual information under limited data; 2) a Bayesian statistical testing method for independence test by formulating the Bayes Factor into the well-defined χ2 statistical test. We apply the proposed methods to both local and global causal discovery algorithms and evaluate their performance against state-of-the-art methods on different benchmark datasets. The experiments show that, by introducing Bayesian approaches, the proposed methods not only outperform the competing methods in terms of accuracy, but also improve efficiency significantly.

Acknowledgments

This work is supported in part by a DARPA grant FA875017-2-0132, and in part by the Rensselaer-IBM AI Research Collaboration (http://airc.rpi.edu), part of the IBM AI Horizons Network (http://ibm.biz/AIHorizons).