Intro

On my second post I’ve stressed how important it is to use the correct adjustment set when trying to estimate a causal relationship between some treatment and exposure variables. The key to using the correct adjustment set is correctly identifying the underling DAG generating a dataset.

A few approaches I can think of to obtain a DAG for causal inference would be:

1. Use domain knowledge and theory
   
2. Pick one of a few candidate DAGs by comparing their fit to the data
   
3. Use algorithms to automatically “learn” the underlying DAG

In this post I’ll explore the last approach. Being able to learn the underlying model DAG from a dataset “automatically” would enable non subject matter experts to do Causal Inference in a very similar fashion to classic ML.

Existing methods for learning DAGs

R has tons of packages for DAG estimation. Fortunately, the causalDisco is a great shiny app that enables searching DAG estimation functions using different criteria.

In this post I’ll restrict the discussion to the scenario where there are no hidden confounders and where the dataset contains both discrete (i.e. categorical) and numeric variables. We’d also like our estimation function to output a DAG with all edges oriented. Below we can see the results of applying these conditions in the causalDisco app:

We get the following list: bnlearn::hc, bnlearn::tabu, bnlearn::mmhc, bnlearn::rsmax2, deal::autosearch and deal::heuristic.

After playing around a bit with the deal functions I found they were too manual for a quick algorithm comparison. I’ve thus decided to use only the bnlearn functions.

Estimated DAG dissimilarity measurements

How can we measure how different our estimated DAG is from the true DAG?

DAG estimation algorithms benchmarking

In this section we’ll use the simMixedDAG package to simulate datasets from a DAG and benchmark different algorithms on recovering it using the above measures.

The DAG and data generating process (DGP) in this section will be based on the General Social Survey data.

The data structure is shown below:

We’ll assume throughout the example that the following is the true DAG:

The DGP will be constructed using the non_parametric_dag_model function from simMixedDAG.

For every algorithm and sample size we simulate 200 datasets and measure the estimated DAG accuracy.

Below we can see the average simulation results for the edge distance measure for different samples sizes:

We can see that the algorithm performance is comparable up to around 27360 observations, after which mmhc and rsmax2 flatten out while hc and tabu approach perfect estimation of the DAG skeleton.

Below we can see the average simulation results for the Hamming distance measure:

We can see that pretty fast all algorithms hit a performance plateau at around 6.

When it comes to Structural Intervention Distance things look even worse:

Let’s revisit that last one with error bars just for mmhc:

We can see that the performance is essentially unchanged across different sample sizes (all confidence intervals overlap).

To get a sense of what different SID scores mean, we can calculate the probability of finding the correct adjustment set given a random selection of a treatment and control pair using the following equation:

\[P(\text{correct adj set)} = 1 – \frac{SID}{\#\{\text{possible treatment – exposure pairs}\}}\] In our case given there are 8 nodes we have \(\#\{\text{possible treatment – exposure pairs}\} = 8 \cdot 7 = 56\).

Below is the resulting graph:

Pretty abysmal.

Summary

While all algorithms converge on the correct DAG skeleton as the sample size increases, none of them is able to orient the edges correctly. So at least for the algorithms surveyed, benchmarked on our specific simulated dataset it would seem automatic DAG learning isn’t a very reliable way to find the correct adjustment set.

All is lost?

Despair not! For in my next blog post I’ll introduce a small package I’ve written to address the issue of DAG edge orientation.