using CausalGraphs
using CausalEstimate
import CausalGraphs: identify, make_graph, to_mermaid, markov_pillow, is_fix, is_p_fix, ID_algorithm, estimate_causal, estimate_id
import CausalEstimate: estimate, confint, pvalue
using DataFrames
using Random
using Statistics
using LoggingThe potential-outcomes chapter asked when an observed data set can tell us about \(E[Y(1)-Y(0)]\). Graphs give one way to state the assumptions behind the answer. A graph records which variables are allowed to cause other variables, and whether some observed variables may share unobserved common causes. Once the graph is stated, identification is a graph problem. Estimation comes after that: estimate the functional that the graph identified.
I use CausalGraphs.jl for the graph and identification steps, then CausalEstimate.jl for backdoor/a-fixable effects, and CausalGraphs.estimate_causal directly for p-fixable and nested-fixable effects. The order is the important part: graph first, identified functional second, estimator third.
A directed acyclic graph, or DAG, has directed arrows such as X -> A. In causal work, an arrow means a direct causal relation is allowed by the model. Absence of an arrow is also a substantive assumption.
An acyclic directed mixed graph, or ADMG, adds bidirected arrows such as A <-> Y. A bidirected edge represents an unobserved common cause. This is useful because many econometric applications have variables we cannot measure: ability, preferences, latent health, firm quality, and so on.
The graph is not an estimator. It is a compact statement of assumptions. The workflow is:
Start with a simple observational study. A is treatment, Y is the outcome, and X is an observed confounder:
g_dag = make_graph(
vertices = [:X, :A, :Y],
di_edges = [(:X, :A), (:X, :Y), (:A, :Y)],
)ADMG([:X, :A, :Y], [(:X, :A), (:X, :Y), (:A, :Y)], Tuple{Symbol, Symbol}[], Dict{Symbol, Vector{Symbol}}(), Dict{Symbol, Bool}(:A => 0, :X => 0, :Y => 0))println("```{mermaid}")
println(to_mermaid(g_dag; direction = "LR"))
println("```")flowchart LR n1_X["X"] n2_A["A"] n3_Y["Y"] n1_X --> n2_A n1_X --> n3_Y n2_A --> n3_Y classDef fixed stroke:#111827,stroke-width:2px linkStyle 0 stroke:#2563eb,stroke-width:1.8px linkStyle 1 stroke:#2563eb,stroke-width:1.8px linkStyle 2 stroke:#2563eb,stroke-width:1.8px
The graph says X affects both treatment and outcome. It also says there is no unobserved common cause between A and Y after we condition on X. In this case the treatment is a-fixable, which corresponds to backdoor-style identification.
id_dag = identify(g_dag, :A, :Y)
(strategy = id_dag.strategy,
markov_pillow = markov_pillow(g_dag, :A; treatment = :A))(strategy = :a_fixable, markov_pillow = [:X])The Markov pillow is the adjustment set used by the estimator. Under this graph,
\[ E[Y(a)] = E_X\{E(Y \mid A=a, X)\}. \tag{3.1}\]
So the graph has translated a causal query into an observed-data functional.
Here is a small simulated data set where X confounds the treatment-outcome association.
Random.seed!(2026)
n = 400
X = randn(n)
pA = logistic.(-0.2 .+ 0.8 .* X)
A = Float64.(rand(n) .< pA)
Y = 1 .+ 1.5 .* A .+ 0.7 .* X .+ 0.3 .* randn(n)
df_dag = DataFrame(X = X, A = A, Y = Y)
first(df_dag, 5)| Row | X | A | Y |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | -0.91813 | 1.0 | 1.76741 |
| 2 | 0.313593 | 1.0 | 2.77638 |
| 3 | -1.23885 | 0.0 | 0.130173 |
| 4 | 0.546413 | 0.0 | 1.43362 |
| 5 | 0.114729 | 0.0 | 0.906974 |
Because the graph is a-fixable, GraphID routes to AIPW backdoor adjustment. The call below estimates the average contrast \(E[Y(1)-Y(0)]\).
dag_res = estimate(
ATE(outcome = :Y, treatment = :A),
GraphID(graph = g_dag),
AIPW(crossfit = 5),
df_dag,
)
effect_table("Backdoor ACE", dag_res)| Row | estimand | estimate | lower_95 | upper_95 | se |
|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Float64 | |
| 1 | Backdoor ACE | 1.465 | 1.265 | 1.666 | 0.1023 |
The graph chose the adjustment set (Markov pillow of A) and the estimator used that set. Access it via dag_res.components[:adjustment_set]. The DGP sets the coefficient on A to 1.5, so the point estimate recovers the true ACE within sampling noise.
crossfit = 5 is the package default, and the number of folds matters more than it looks. With crossfit = 2 each nuisance is fit on only half the sample, and the extra noise in the out-of-fold predictions propagates into both the point estimate and the influence-function variance — fewer folds means a noisier estimate reported with a wider interval, not a cheap approximation to the same answer. Two folds is the smallest value the package allows; it is not a sensible default. The next chapter measures the size of this effect on real data.
Now suppose treatment and outcome share an unobserved common cause. A simple adjustment argument no longer works. The front-door ADMG has a mediator M that transmits the effect of A to Y, while A and Y are hidden-confounded:
g_fd = make_graph(
vertices = [:A, :M, :Y],
di_edges = [(:A, :M), (:M, :Y)],
bi_edges = [(:A, :Y)],
)ADMG([:A, :M, :Y], [(:A, :M), (:M, :Y)], [(:A, :Y)], Dict{Symbol, Vector{Symbol}}(), Dict{Symbol, Bool}(:M => 0, :A => 0, :Y => 0))println("```{mermaid}")
println(to_mermaid(g_fd; direction = "LR"))
println("```")flowchart LR n1_A["A"] n2_M["M"] n3_Y["Y"] n1_A --> n2_M n2_M --> n3_Y n1_A <--> n3_Y classDef fixed stroke:#111827,stroke-width:2px linkStyle 0 stroke:#2563eb,stroke-width:1.8px linkStyle 1 stroke:#2563eb,stroke-width:1.8px linkStyle 2 stroke:#dc2626,stroke-width:1.8px
The bidirected edge A <-> Y means that ordinary backdoor adjustment is not available. But the graph is p-fixable, the ADMG generalization that includes front-door identification.
id_fd = identify(g_fd, :A, :Y)
(strategy = id_fd.strategy,
a_fixable = is_fix(g_fd, :A),
p_fixable = is_p_fix(g_fd, :A))(strategy = :p_fixable, a_fixable = false, p_fixable = true)The intuition is that M carries the causal effect of A, and M is not itself hidden-confounded with A. The graph lets us use mediator information to recover the causal effect even though A and Y are confounded.
Simulate a front-door setting. The latent U is not included in the observed data; it is represented in the graph by A <-> Y.
Random.seed!(2027)
U = randn(n)
A_fd = Float64.(rand(n) .< logistic.(-0.1 .+ 0.8 .* U))
M = Float64.(rand(n) .< logistic.(-0.4 .+ 1.2 .* A_fd))
Y_fd = Float64.(rand(n) .< logistic.(-1 .+ 0.9 .* M .+ 0.7 .* U))
df_fd = DataFrame(A = A_fd, M = M, Y = Y_fd)
first(df_fd, 5)| Row | A | M | Y |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 1.0 | 0.0 | 0.0 |
| 2 | 0.0 | 0.0 | 0.0 |
| 3 | 0.0 | 0.0 | 0.0 |
| 4 | 0.0 | 0.0 | 1.0 |
| 5 | 0.0 | 1.0 | 1.0 |
The DGP has a known truth. Under \(do(A=a)\) the mediator is \(M \sim \text{Bern}(\Lambda(-0.4 + 1.2a))\) while \(U\) keeps its marginal \(N(0,1)\) distribution, so
\[ \tau = \big[\Lambda(0.8) - \Lambda(-0.4)\big]\; \big[g(-0.1) - g(-1.0)\big], \qquad g(c) = E_U\big[\Lambda(c + 0.7U)\big], \tag{3.2}\]
which evaluates to \(0.2887 \times 0.1894 = 0.0547\). The naive contrast \(E[Y\mid A=1]-E[Y\mid A=0]\) equals \(0.158\) in the population this DGP defines — almost three times the causal effect, because \(U\) raises both the chance of treatment and the outcome. That gap is what the front-door step removes. The table below reports the sample version of the naive contrast alongside the front-door estimate.
Since identify() returns :p_fixable, this is a front-door effect. CausalEstimate.jl currently handles backdoor (a-fixable) graphs; for p-fixable effects we call CausalGraphs.estimate_causal directly, which routes to NPS TMLE.
fd_res = with_logger(NullLogger()) do
CausalGraphs.estimate_causal(
a = [1.0, 0.0],
data = df_fd,
graph = g_fd,
treatment = :A,
outcome = :Y,
n_iter = 80,
cvg_criteria = 0.02,
truncate_lower = 0.01,
truncate_upper = 0.99,
)
end
naive_fd = mean(df_fd.Y[df_fd.A .== 1]) - mean(df_fd.Y[df_fd.A .== 0])
vcat(
DataFrame(estimand = ["Naive E[Y|A=1] - E[Y|A=0]"], estimate = [r4(naive_fd)],
lower_95 = [missing], upper_95 = [missing], se = [missing]),
effect_table("Front-door ACE", fd_res[:TMLE]),
DataFrame(estimand = ["Truth"], estimate = [0.0547],
lower_95 = [missing], upper_95 = [missing], se = [missing]);
cols = :union,
)| Row | estimand | estimate | lower_95 | upper_95 | se |
|---|---|---|---|---|---|
| String | Float64 | Float64? | Float64? | Missing | |
| 1 | Naive E[Y|A=1] - E[Y|A=0] | 0.1014 | missing | missing | missing |
| 2 | Front-door ACE | 0.05239 | 0.02208 | 0.0827 | missing |
| 3 | Truth | 0.0547 | missing | missing | missing |
The naive contrast is far above the truth; the front-door estimate brackets it.
Backdoor/a-fixable, p-fixable/front-door and nested-fixable are useful because they are easy to estimate. ADMGs can be more general. The Pearl-Shpitser ID algorithm asks the broader question:
Can P(Y | do(A)) be written using only the observed joint law P(V)?The algorithm works recursively:
For the front-door graph, the general ID algorithm also succeeds:
id_alg_fd = ID_algorithm(g_fd, :A, :Y)
(identified = id_alg_fd.identified,
expression = string(id_alg_fd.expression))(identified = true, expression = "sum_{M} P(M | A) * sum_{A} P(A) * P(Y | A, M)")The printed expression reuses A as the inner summation index; read it as \(\sum_M P(M\mid A)\sum_{A'} P(Y\mid A', M)P(A')\), where the inner \(A'\) ranges over the marginal distribution of treatment and is distinct from the \(A\) being intervened on. The expression is symbolic. It is not yet an estimate. For finite-support variables, CausalGraphs.estimate_causal evaluates the symbolic functional by empirical probabilities and attaches a finite-support EIF standard error via the :IDPlugin key.
id_res = with_logger(NullLogger()) do
estimate_id(
a = [1.0, 0.0],
data = df_fd,
graph = g_fd,
treatment = :A,
outcome = :Y,
)
end
effect_table("ID plug-in ACE", id_res[:IDPlugin];
se = id_res[:IDPlugin].standard_error)| Row | estimand | estimate | lower_95 | upper_95 | se |
|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Float64 | |
| 1 | ID plug-in ACE | 0.05346 | 0.02278 | 0.08413 | 0.01565 |
This plug-in estimator is useful for discrete examples and diagnostics. It is not a universal TMLE for every continuous ID functional.
The bow graph has both a direct causal edge and unobserved confounding between the same two variables:
g_bow = make_graph(
vertices = [:A, :Y],
di_edges = [(:A, :Y)],
bi_edges = [(:A, :Y)],
)ADMG([:A, :Y], [(:A, :Y)], [(:A, :Y)], Dict{Symbol, Vector{Symbol}}(), Dict{Symbol, Bool}(:A => 0, :Y => 0))println("```{mermaid}")
println(to_mermaid(g_bow; direction = "LR"))
println("```")flowchart LR n1_A["A"] n2_Y["Y"] n1_A --> n2_Y n1_A <--> n2_Y classDef fixed stroke:#111827,stroke-width:2px linkStyle 0 stroke:#2563eb,stroke-width:1.8px linkStyle 1 stroke:#dc2626,stroke-width:1.8px
The effect of A on Y is not identified in this ADMG. The observed association mixes the direct causal effect with the latent common cause.
id_bow = identify(g_bow, :A, :Y)
id_alg_bow = ID_algorithm(g_bow, :A, :Y)
(route = id_bow.strategy,
id_algorithm_identified = id_alg_bow.identified,
hedge = id_alg_bow.hedge)(route = :not_identified, id_algorithm_identified = false, hedge = (S = [:Y], F = [:A, :Y], treatment = [:A], outcome = [:Y]))For applied work, a useful workflow is:
identify(graph, treatment, outcome);estimate(ATE(...), GraphID(graph=...), AIPW(...), df); for p-fixable or nested-fixable effects, call CausalGraphs.estimate_causal(...);