Example of causal discovery
Suppose we have this DAG:
We define a layout which we will use for plotting the graphs in this vignette:
We can create data from a linear Gaussian model corresponding to the
above DAG using generate_dag_data(). We generate 10,000
samples with a fixed random seed from the DAG above, where we let the
edge coefficients be sampled with absolute values between 0.1 and 0.9
and assigned random signs, and where the standard deviation of the
additive Gaussian noise at each node is sampled from a log-uniform
distribution between 0.3 and 2.
data_linear <- generate_dag_data(
cg,
n = 10000,
seed = 1405,
coef_range = c(0.1, 0.9),
error_sd = c(0.3, 2)
)
head(data_linear)
#> # A tibble: 6 × 5
#> Z X3 X1 X2 Y
#> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.632 0.297 0.773 -0.256 -0.730
#> 2 -1.67 -1.10 -1.20 1.10 -0.722
#> 3 -0.214 0.713 1.09 0.0264 0.261
#> 4 -1.61 -0.171 -1.21 -0.0346 0.0746
#> 5 1.57 0.0517 0.402 -0.874 1.40
#> 6 -1.08 0.178 -1.09 0.189 0.0591The R code used to generate the data is stored as an attribute of the data frame:
attr(data_linear, "generating_model")
#> $dgp
#> $dgp$X3
#> rnorm(n, sd = 0.95)
#>
#> $dgp$X2
#> X3 * 0.159 + rnorm(n, sd = 0.507)
#>
#> $dgp$Z
#> rnorm(n, sd = 1.794)
#>
#> $dgp$X1
#> Z * 0.3 + rnorm(n, sd = 1.634)
#>
#> $dgp$Y
#> X1 * 0.61 + X2 * 0.31 + rnorm(n, sd = 1.216)The goal is now to recover the causal structure from the data alone.
We will use the PC algorithm to learn the causal structure. One of its key assumptions is causal sufficiency, meaning that there are no unobserved confounders. We can use the PC algorithm from any of the “tetrad”, “pcalg”, or “bnlearn” engines to learn the structure from the data. The motivation for having multiple engines is that they may implement different algorithms, tests, or scoring methods, or the same algorithm may be implemented differently across engines.
Below, we set up the PC algorithm using Fisher’s Z test, a
significance level of 0.05, and pcalg as the engine. To do
so, we first define the PC method using the pc() function,
and then pass it to the disco() function along with the
data.
pc_pcalg <- pc(engine = "pcalg", test = "fisher_z", alpha = 0.05)
pc_result_pcalg <- disco(data = data_linear, method = pc_pcalg)We can visualize the results using the plot() function,
where we explicitly provide the layout defined above so graphs are
easier to compare.
The PC algorithm recovers the correct causal structure up to Markov equivalence, represented as a CPDAG. A CPDAG represents the equivalence class of DAGs that encode the same conditional independencies.
The first notable feature of this plot is that some edges are
directed, while others are undirected. For example, the edge from
X1 to Y is directed, indicating a causal
effect of X1 on Y, but not in the reverse
direction. In contrast, the edge between X2 and
X3 is undirected, indicating that the data alone do not
provide sufficient information to determine the causal direction. Both
orientations X2 %-->% X3 and X3 %-->% X2
are compatible with the observed conditional independencies.
We demonstrate the non-identifiability of the causal direction
between X2 and X3 by reversing the direction
of this edge in the data-generating process above and applying the PC
algorithm to the resulting data set.
cg_reverse <- caugi::caugi(
Z %-->% X1,
X2 %-->% X3,
X1 %-->% Y,
X2 %-->% Y
)
data_linear_reverse <- generate_dag_data(
cg_reverse,
n = 10000,
seed = 1405,
coef_range = c(0.1, 0.9),
error_sd = c(0.3, 2)
)
pc_result_reverse <- disco(data = data_linear_reverse, method = pc_pcalg)
plot(pc_result_reverse, layout = layout, main = "PC (pcalg) reversed")
We learn the same causal structure as before, demonstrating that the
direction of influence between X2 and X3
cannot be determined from the data alone.
