21 From Graph to Estimate

Code

library(dagitty)
library(igraph)
library(causaleffect)
library(ggplot2)

The previous chapters covered discovery: algorithms that learn a graph from data. But a graph is not an estimate. After we have a graph, we still need to ask two questions:

Identification: is the target effect a function of the observed-data distribution?
Estimation: once the identifying formula is known, estimate it from the data.

In R, dagitty handles identification via the backdoor and front-door criteria; causaleffect implements the full Pearl–Shpitser ID algorithm for acyclic directed mixed graphs (ADMGs) with hidden confounders. Estimation is done directly from the identifying formula — a step that is explicit in R, making the connection between the formula and the estimator transparent.

21.1 Backdoor: standard adjustment

Start with the simple case. $X$ is a common cause of treatment $A$ and outcome $Y$, and $X$ is observed. Backdoor adjustment applies.

Code

# X is a common cause of A and Y
g_bd <- dagitty("dag { X -> A; X -> Y; A -> Y }")
exposures(g_bd) <- "A"
outcomes(g_bd)  <- "Y"

sets <- adjustmentSets(g_bd, exposure = "A", outcome = "Y")
cat("Adjustment sets:\n")

Adjustment sets:

Code

print(sets)

{ X }

dagitty finds {X} as the minimal adjustment set. Adjusting on $X$ blocks the only open backdoor path $A \leftarrow X \rightarrow Y$.

We simulate data matching that graph. We draw $n = 2000$ observations with $X \sim N(0,1)$, treatment $A \sim \text{Bernoulli}(\Lambda(X))$, and outcome $Y = 2A + X + 0.5\varepsilon$ with $\varepsilon \sim N(0,1)$. The effect is constant, so the true ACE is 2. We then estimate it with AIPW using the adjustment set the graph chose, and take a 95% interval from 500 bootstrap resamples.

Code

set.seed(1)
n     <- 2000
X     <- rnorm(n)
A     <- as.numeric(runif(n) < 1 / (1 + exp(-X)))  # logistic: A depends on X
Y     <- 2 * A + X + 0.5 * rnorm(n)                 # true ACE = 2
data_bd <- data.frame(X = X, A = A, Y = Y)

# AIPW (doubly-robust) estimator for the backdoor-adjusted ACE
# Outcome model: E[Y | A, X]
mu1 <- predict(lm(Y ~ A + X, data = data_bd), newdata = transform(data_bd, A = 1))
mu0 <- predict(lm(Y ~ A + X, data = data_bd), newdata = transform(data_bd, A = 0))

# Propensity score: P(A = 1 | X)
ps  <- predict(glm(A ~ X, data = data_bd, family = binomial), type = "response")

# AIPW influence function
aipw <- mean(
  (data_bd$A / ps) * (Y - mu1) - ((1 - data_bd$A) / (1 - ps)) * (Y - mu0) +
  mu1 - mu0
)

cat(sprintf("Backdoor AIPW ACE (true = 2.0): %.3f\n", aipw))

Backdoor AIPW ACE (true = 2.0): 1.990

Code

# Bootstrap 95% CI
set.seed(99)
boot_bd <- replicate(500, {
  idx  <- sample.int(n, replace = TRUE)
  bd   <- data_bd[idx, ]
  m1   <- predict(lm(Y ~ A + X, data = bd), newdata = transform(bd, A = 1))
  m0   <- predict(lm(Y ~ A + X, data = bd), newdata = transform(bd, A = 0))
  p    <- predict(glm(A ~ X, data = bd, family = binomial), type = "response")
  mean((bd$A / p) * (bd$Y - m1) - ((1 - bd$A) / (1 - p)) * (bd$Y - m0) + m1 - m0)
})
ci_bd <- quantile(boot_bd, c(0.025, 0.975))
cat(sprintf("95%% CI: [%.3f, %.3f]\n", ci_bd[1], ci_bd[2]))

95% CI: [1.942, 2.043]

AIPW returns 1.990 with a 95% bootstrap interval of $[1.942, 2.043]$, covering the true 2. The AIPW estimator is consistent if either the outcome model or the propensity score is correctly specified; here both are, since the DGP is linear in the outcome and logistic in the treatment.

21.2 Front-door: p-fixable effects

Now suppose $A$ and $Y$ have an unmeasured common cause $U$, but there is a measured mediator $M$ on the path $A \rightarrow M \rightarrow Y$. Backdoor adjustment fails because $U$ is unobserved. The front-door formula can still identify the effect.

In dagitty, the bidirected edge $A \leftrightarrow Y$ (representing the latent $U$) appears as A <-> Y. For causaleffect, the same edge is encoded as two directed edges A -> Y and Y -> A both marked with description = "U" in the igraph object.

Code

# dagitty for display and adjustment-set check
g_fd <- dagitty("dag { A -> M; M -> Y; A <-> Y }")
exposures(g_fd) <- "A"
outcomes(g_fd)  <- "Y"

# An empty dagitty adjustment-set list prints nothing at all, so report the
# count rather than calling print() on it and appearing to produce no output.
sets_fd <- adjustmentSets(g_fd, exposure = "A", outcome = "Y")
cat("Valid backdoor adjustment sets:", length(sets_fd), "\n")

Valid backdoor adjustment sets: 0

Code

# causaleffect for ID algorithm: bidirected A <-> Y = two directed edges + description "U"
g_fd_ig <- graph_from_literal(A -+ M, M -+ Y, A -+ Y, Y -+ A)
# Mark the A<->Y bidirected edge by endpoint name rather than a hardcoded
# edge index: igraph happens to sort edges by source vertex today (giving
# order A->M(1), A->Y(2), M->Y(3), Y->A(4)), but that ordering isn't part of
# igraph's documented contract, so indices like c(2,4) can silently break if
# igraph's internals or this graph's construction ever changes.
edge_ends <- ends(g_fd_ig, E(g_fd_ig))
confounded <- (edge_ends[, 1] == "A" & edge_ends[, 2] == "Y") |
              (edge_ends[, 1] == "Y" & edge_ends[, 2] == "A")
g_fd_ig <- set_edge_attr(g_fd_ig, "description", index = which(confounded), value = "U")

expr_fd <- causal.effect(y = "Y", x = "A", G = g_fd_ig, simp = TRUE)
cat("\nIdentifying expression:\n")


Identifying expression:

Code

cat(expr_fd, "\n")

\sum_{M}P(M|A)\left(\sum_{A}P(Y|A,M)P(A)\right)

There are no valid backdoor adjustment sets, as expected: $U$ is unobserved, so no set of measured variables blocks the path it opens. And the expression causal.effect returns is the front-door formula, $\sum_M P(M \mid A) \sum_{A'} P(Y \mid A', M)\, P(A')$.

We now simulate data for it. We draw $n = 3000$ observations with an unmeasured confounder $U \sim N(0,1)$, treatment $A \sim \text{Bernoulli}(\Lambda(U))$, mediator $M \sim \text{Bernoulli}(\Lambda(0.5 + 2A))$, and outcome $Y = 1.5M + 2U + 0.5\varepsilon$ with $\varepsilon \sim N(0,1)$. The mediator depends on $A$ alone and the outcome depends on $M$ and $U$ but not directly on $A$, which is what the front-door criterion requires.

The truth follows from the DGP. Intervening on $A$ changes $Y$ only through $M$, so the ACE is $1.5 \times [P(M=1 \mid A=1) - P(M=1 \mid A=0)] = 1.5 \times [\Lambda(2.5) - \Lambda(0.5)] = 1.5 \times 0.302 = 0.453$.

Code

set.seed(2)
n     <- 3000
U     <- rnorm(n)                                     # unmeasured confounder
A_fd  <- as.numeric(runif(n) < 1 / (1 + exp(-U)))     # A depends on U
M_fd  <- as.numeric(runif(n) < 1 / (1 + exp(-(0.5 + 2 * A_fd))))  # M depends on A only
Y_fd  <- 1.5 * M_fd + 2 * U + 0.5 * rnorm(n)          # Y depends on M and U

# True ACE: E[Y | do(A=1)] - E[Y | do(A=0)]
# = 1.5 * (E[M | A=1] - E[M | A=0])
true_EMA1  <- 1 / (1 + exp(-(0.5 + 2)))   # scalar P(M=1 | A=1)
true_EMA0  <- 1 / (1 + exp(-0.5))         # scalar P(M=1 | A=0)
true_ace_fd <- 1.5 * (true_EMA1 - true_EMA0)
cat(sprintf("True front-door ACE = %.3f\n", true_ace_fd))

True front-door ACE = 0.453

Code

data_fd <- data.frame(A = A_fd, M = M_fd, Y = Y_fd)

Estimate the front-door functional using the plug-in formula:

Code

front_door_ace <- function(data) {
  pA <- mean(data$A)
  pM_given_A <- prop.table(table(data$A, data$M), margin = 1)
  pY_given_AM <- with(data, tapply(Y, list(A, M), mean))

  ey_do <- function(a) {
    sum(sapply(c(0, 1), function(m) {
      p_m_given_a <- pM_given_A[as.character(a), as.character(m)]
      inner <- sum(sapply(c(0, 1), function(ap)
        pY_given_AM[as.character(ap), as.character(m)] *
          ifelse(ap == 1, pA, 1 - pA)))
      p_m_given_a * inner
    }))
  }
  ey_do(1) - ey_do(0)
}

est_fd <- front_door_ace(data_fd)

set.seed(77)
boot_fd <- replicate(500, {
  idx <- sample.int(nrow(data_fd), replace = TRUE)
  front_door_ace(data_fd[idx, ])
})
ci_fd  <- quantile(boot_fd, c(0.025, 0.975))

cat(sprintf("Front-door ACE: %.3f  [95%% CI: %.3f, %.3f]\n",
            est_fd, ci_fd[1], ci_fd[2]))

Front-door ACE: 0.499  [95% CI: 0.432, 0.577]

Code

# Naïve mean difference for comparison
naive_ace <- mean(Y_fd[A_fd == 1]) - mean(Y_fd[A_fd == 0])
cat(sprintf("Naïve mean difference (biased!): %.3f\n", naive_ace))

Naïve mean difference (biased!): 2.083

Code

cat(sprintf("True ACE:                        %.3f\n", true_ace_fd))

True ACE:                        0.453

The front-door estimate is 0.499 with a 95% bootstrap interval of $[0.432, 0.577]$, which covers the true 0.453. The naive mean difference is 2.083 — more than four times the truth, and outside the front-door interval by a factor of four. The formula recovers the effect even though $U$ is never measured; the naive comparison is picking up the path through $U$, which enters the outcome with coefficient 2 and the treatment through the logit.

The bias here is far larger than in the front-door example of the graphs chapter, where the naive contrast was 0.14 against a truth of 0.055. The difference is the strength of the confounding: $U$ enters $Y$ with coefficient 2 here against 0.7 there. Nothing about the front-door formula changes; what changes is how much work it is doing.

21.3 When identification fails

Not every graph identifies the effect. If $A$ and $Y$ have a hidden common cause and there is no front-door variable, the effect is not identified. The graph below is the bow: a directed edge $A \to Y$ plus a bidirected $A \leftrightarrow Y$, and nothing else. Remove the mediator $M$ from the front-door graph and this is what is left.

Code

# A → Y with unmeasured confounding (A <-> Y) and no mediator: the "bow"
# graph needs THREE igraph edges — the directed A → Y plus the
# bidirected pair (encoded as two opposite edges marked "U").
# graph_from_literal() silently collapses a duplicated A -+ Y, so we
# build the multigraph with make_graph().
g_unident_ig <- make_graph(c("A","Y",  "A","Y",  "Y","A"), directed = TRUE)
# make_graph() preserves the literal edge order given above (unlike
# graph_from_literal(), which re-sorts), so index 2:3 is reliable as long as
# the c(...) vector above isn't reordered -- mark by endpoint pair anyway so
# a future edit to that vector can't silently mismark the real A->Y edge.
edge_ends <- ends(g_unident_ig, E(g_unident_ig))
# The first A->Y edge is the real causal edge; the second (a duplicate,
# per `duplicated()`) together with Y->A forms the bidirected U pair.
u_edges <- which(((edge_ends[,1]=="A" & edge_ends[,2]=="Y") & duplicated(edge_ends)) |
                  (edge_ends[,1]=="Y" & edge_ends[,2]=="A"))
g_unident_ig <- set_edge_attr(g_unident_ig, "description",
                               index = u_edges, value = "U")

result <- tryCatch(
  causal.effect(y = "Y", x = "A", G = g_unident_ig, simp = TRUE),
  error = function(e) conditionMessage(e)
)
cat("Result from causal.effect:\n", result, "\n")

Result from causal.effect:
 Not identifiable.

causaleffect raises an error, printed above as “Not identifiable.”, which the tryCatch captures. This is a useful failure. The algorithm did not return a formula that would quietly estimate the wrong thing; it reported that no function of the observed distribution equals the causal effect under this graph. No amount of data helps, because the obstacle is the graph and not the sample size.

The honest response in applied work is to report the non-identifiability, then either:

Find additional variables that close the unmeasured-confounding paths
Apply sensitivity analysis to bound the effect under maintained assumptions
Find an instrumental variable for a LATE-style identification (IV chapter)

Saying that the hidden confounding is probably small is sensitivity analysis, not identification.

21.4 End-to-end: discover, identify, estimate

Discovery algorithms return graphs that can be passed to the identification step:

Code

library(pcalg)

# Step 1: discover the graph (chapters on causal discovery)
suffStat <- list(C = cor(data_obs), n = nrow(data_obs))
pc_fit   <- pc(suffStat, indepTest = gaussCItest,
               alpha = 0.05, p = ncol(data_obs))

# Step 2: check for backdoor adjustment set using dagitty
# (convert pcalg CPDAG → dagitty for convenience).
# NOTE the t(): as(graphNEL, "matrix") is from→to, but pcalg2dagitty
# expects the amat.cpdag coding, which is its transpose — without it
# every directed edge is silently reversed.
g_dag <- pcalg2dagitty(t(as(pc_fit@graph, "matrix")),
                        colnames(data_obs), type = "cpdag")
sets  <- adjustmentSets(g_dag, exposure = "A", outcome = "Y")

# Step 3: if an adjustment set is valid for every DAG in the class, estimate
#         it with AIPW
# Step 4: if the CPDAG is still unresolved, do NOT hand it to
#         causaleffect::causal.effect() -- that function expects a specified
#         DAG/ADMG, and different orientations in the class can imply
#         different effects, different valid adjustment sets, or a different
#         identification status. Instead either
#           (a) orient the remaining edges with background knowledge and then
#               run ordinary identification, or
#           (b) report what the whole class permits: pcalg::ida() for
#               linear-Gaussian total effects, otherwise enumerate the
#               compatible DAGs and keep only conclusions that hold across
#               all of them.
# Step 5: estimate the returned functional by plug-in or regression

Step 4 is where this workflow most often goes wrong, so it is worth stating plainly: a CPDAG is an equivalence class, not a graph. An empty adjustment-set result for a CPDAG is not permission to pass one arbitrary orientation to the ID algorithm. The causal discovery chapter sets out what to do instead.

In applied work, I would not trust the discovered graph blindly. A more reasonable workflow is:

Hypothesise a graph based on subject-matter knowledge
Discover a data-driven graph using PC/GES/FCI
Compare the two — disagreement points to substantive questions
Identify under both graphs separately
Estimate under both, and report the range as part of the result

If the estimate changes a lot across plausible graphs, the data alone is not settling the question.

21.5 What changes after identification

Backdoor adjustment is not always available because many realistic DAGs have hidden confounders. Front-door, nested, and general ID algorithms may still identify the effect, although their formulas quickly become difficult to derive by hand.

The identification result also determines the estimator. AIPW for a backdoor adjustment uses a different formula from a front-door plug-in estimator. Keeping these steps together helps us avoid applying an estimator that does not match the graph. It also forces us to state which variables support an unconfoundedness claim instead of leaving the claim implicit.

21.6 Summary

A graph is not an estimate. Identification is the bridge.
dagitty handles backdoor and front-door checks; causaleffect handles general ADMG identification.
Backdoor effects can be estimated with AIPW.
Front-door effects require estimating the front-door functional.
If the effect is not identified, report that instead of forcing a biased regression estimate.

# From Graph to Estimate ```{r} #| include: false library(dagitty) library(igraph) library(causaleffect) library(ggplot2) library(dplyr) ``` ```{r} #| eval: false library(dagitty) library(igraph) library(causaleffect) library(ggplot2) ``` The previous chapters covered discovery: algorithms that learn a graph from data. But a graph is not an estimate. After we have a graph, we still need to ask two questions: 1. **Identification**: is the target effect a function of the observed-data distribution? 2. **Estimation**: once the identifying formula is known, estimate it from the data. In R, `dagitty` handles identification via the backdoor and front-door criteria; `causaleffect` implements the full Pearl–Shpitser ID algorithm for acyclic directed mixed graphs (ADMGs) with hidden confounders. Estimation is done directly from the identifying formula — a step that is explicit in R, making the connection between the formula and the estimator transparent. ## Backdoor: standard adjustment Start with the simple case. $X$ is a common cause of treatment $A$ and outcome $Y$, and $X$ is observed. Backdoor adjustment applies. ```{r} #| label: backdoor-graph #| cache: true # X is a common cause of A and Y g_bd <- dagitty("dag { X -> A; X -> Y; A -> Y }") exposures(g_bd) <- "A" outcomes(g_bd) <- "Y" sets <- adjustmentSets(g_bd, exposure = "A", outcome = "Y") cat("Adjustment sets:\n") print(sets) ``` `dagitty` finds `{X}` as the minimal adjustment set. Adjusting on $X$ blocks the only open backdoor path $A \leftarrow X \rightarrow Y$. We simulate data matching that graph. We draw $n = 2000$ observations with $X \sim N(0,1)$, treatment $A \sim \text{Bernoulli}(\Lambda(X))$, and outcome $Y = 2A + X + 0.5\varepsilon$ with $\varepsilon \sim N(0,1)$. The effect is constant, so the true ACE is 2. We then estimate it with AIPW using the adjustment set the graph chose, and take a 95% interval from 500 bootstrap resamples. ```{r} #| label: backdoor-sim #| cache: true set.seed(1) n <- 2000 X <- rnorm(n) A <- as.numeric(runif(n) < 1 / (1 + exp(-X))) # logistic: A depends on X Y <- 2 * A + X + 0.5 * rnorm(n) # true ACE = 2 data_bd <- data.frame(X = X, A = A, Y = Y) # AIPW (doubly-robust) estimator for the backdoor-adjusted ACE # Outcome model: E[Y | A, X] mu1 <- predict(lm(Y ~ A + X, data = data_bd), newdata = transform(data_bd, A = 1)) mu0 <- predict(lm(Y ~ A + X, data = data_bd), newdata = transform(data_bd, A = 0)) # Propensity score: P(A = 1 | X) ps <- predict(glm(A ~ X, data = data_bd, family = binomial), type = "response") # AIPW influence function aipw <- mean( (data_bd$A / ps) * (Y - mu1) - ((1 - data_bd$A) / (1 - ps)) * (Y - mu0) + mu1 - mu0 ) cat(sprintf("Backdoor AIPW ACE (true = 2.0): %.3f\n", aipw)) # Bootstrap 95% CI set.seed(99) boot_bd <- replicate(500, { idx <- sample.int(n, replace = TRUE) bd <- data_bd[idx, ] m1 <- predict(lm(Y ~ A + X, data = bd), newdata = transform(bd, A = 1)) m0 <- predict(lm(Y ~ A + X, data = bd), newdata = transform(bd, A = 0)) p <- predict(glm(A ~ X, data = bd, family = binomial), type = "response") mean((bd$A / p) * (bd$Y - m1) - ((1 - bd$A) / (1 - p)) * (bd$Y - m0) + m1 - m0) }) ci_bd <- quantile(boot_bd, c(0.025, 0.975)) cat(sprintf("95%% CI: [%.3f, %.3f]\n", ci_bd[1], ci_bd[2])) ``` AIPW returns 1.990 with a 95% bootstrap interval of $[1.942, 2.043]$, covering the true 2. The AIPW estimator is consistent if either the outcome model or the propensity score is correctly specified; here both are, since the DGP is linear in the outcome and logistic in the treatment. ## Front-door: p-fixable effects Now suppose $A$ and $Y$ have an unmeasured common cause $U$, but there is a measured mediator $M$ on the path $A \rightarrow M \rightarrow Y$. Backdoor adjustment fails because $U$ is unobserved. The front-door formula can still identify the effect. In `dagitty`, the bidirected edge $A \leftrightarrow Y$ (representing the latent $U$) appears as `A <-> Y`. For `causaleffect`, the same edge is encoded as two directed edges `A -> Y` and `Y -> A` both marked with `description = "U"` in the igraph object. ```{r} #| label: frontdoor-graph #| cache: true # dagitty for display and adjustment-set check g_fd <- dagitty("dag { A -> M; M -> Y; A <-> Y }") exposures(g_fd) <- "A" outcomes(g_fd) <- "Y" # An empty dagitty adjustment-set list prints nothing at all, so report the # count rather than calling print() on it and appearing to produce no output. sets_fd <- adjustmentSets(g_fd, exposure = "A", outcome = "Y") cat("Valid backdoor adjustment sets:", length(sets_fd), "\n") # causaleffect for ID algorithm: bidirected A <-> Y = two directed edges + description "U" g_fd_ig <- graph_from_literal(A -+ M, M -+ Y, A -+ Y, Y -+ A) # Mark the A<->Y bidirected edge by endpoint name rather than a hardcoded # edge index: igraph happens to sort edges by source vertex today (giving # order A->M(1), A->Y(2), M->Y(3), Y->A(4)), but that ordering isn't part of # igraph's documented contract, so indices like c(2,4) can silently break if # igraph's internals or this graph's construction ever changes. edge_ends <- ends(g_fd_ig, E(g_fd_ig)) confounded <- (edge_ends[, 1] == "A" & edge_ends[, 2] == "Y") | (edge_ends[, 1] == "Y" & edge_ends[, 2] == "A") g_fd_ig <- set_edge_attr(g_fd_ig, "description", index = which(confounded), value = "U") expr_fd <- causal.effect(y = "Y", x = "A", G = g_fd_ig, simp = TRUE) cat("\nIdentifying expression:\n") cat(expr_fd, "\n") ``` There are no valid backdoor adjustment sets, as expected: $U$ is unobserved, so no set of measured variables blocks the path it opens. And the expression `causal.effect` returns is the front-door formula, $\sum_M P(M \mid A) \sum_{A'} P(Y \mid A', M)\, P(A')$. We now simulate data for it. We draw $n = 3000$ observations with an unmeasured confounder $U \sim N(0,1)$, treatment $A \sim \text{Bernoulli}(\Lambda(U))$, mediator $M \sim \text{Bernoulli}(\Lambda(0.5 + 2A))$, and outcome $Y = 1.5M + 2U + 0.5\varepsilon$ with $\varepsilon \sim N(0,1)$. The mediator depends on $A$ alone and the outcome depends on $M$ and $U$ but not directly on $A$, which is what the front-door criterion requires. The truth follows from the DGP. Intervening on $A$ changes $Y$ only through $M$, so the ACE is $1.5 \times [P(M=1 \mid A=1) - P(M=1 \mid A=0)] = 1.5 \times [\Lambda(2.5) - \Lambda(0.5)] = 1.5 \times 0.302 = 0.453$. ```{r} #| label: frontdoor-sim #| cache: true set.seed(2) n <- 3000 U <- rnorm(n) # unmeasured confounder A_fd <- as.numeric(runif(n) < 1 / (1 + exp(-U))) # A depends on U M_fd <- as.numeric(runif(n) < 1 / (1 + exp(-(0.5 + 2 * A_fd)))) # M depends on A only Y_fd <- 1.5 * M_fd + 2 * U + 0.5 * rnorm(n) # Y depends on M and U # True ACE: E[Y | do(A=1)] - E[Y | do(A=0)] # = 1.5 * (E[M | A=1] - E[M | A=0]) true_EMA1 <- 1 / (1 + exp(-(0.5 + 2))) # scalar P(M=1 | A=1) true_EMA0 <- 1 / (1 + exp(-0.5)) # scalar P(M=1 | A=0) true_ace_fd <- 1.5 * (true_EMA1 - true_EMA0) cat(sprintf("True front-door ACE = %.3f\n", true_ace_fd)) data_fd <- data.frame(A = A_fd, M = M_fd, Y = Y_fd) ``` Estimate the front-door functional using the plug-in formula: ```{r} #| label: frontdoor-est #| cache: true front_door_ace <- function(data) { pA <- mean(data$A) pM_given_A <- prop.table(table(data$A, data$M), margin = 1) pY_given_AM <- with(data, tapply(Y, list(A, M), mean)) ey_do <- function(a) { sum(sapply(c(0, 1), function(m) { p_m_given_a <- pM_given_A[as.character(a), as.character(m)] inner <- sum(sapply(c(0, 1), function(ap) pY_given_AM[as.character(ap), as.character(m)] * ifelse(ap == 1, pA, 1 - pA))) p_m_given_a * inner })) } ey_do(1) - ey_do(0) } est_fd <- front_door_ace(data_fd) set.seed(77) boot_fd <- replicate(500, { idx <- sample.int(nrow(data_fd), replace = TRUE) front_door_ace(data_fd[idx, ]) }) ci_fd <- quantile(boot_fd, c(0.025, 0.975)) cat(sprintf("Front-door ACE: %.3f [95%% CI: %.3f, %.3f]\n", est_fd, ci_fd[1], ci_fd[2])) # Naïve mean difference for comparison naive_ace <- mean(Y_fd[A_fd == 1]) - mean(Y_fd[A_fd == 0]) cat(sprintf("Naïve mean difference (biased!): %.3f\n", naive_ace)) cat(sprintf("True ACE: %.3f\n", true_ace_fd)) ``` The front-door estimate is 0.499 with a 95% bootstrap interval of $[0.432, 0.577]$, which covers the true 0.453. The naive mean difference is 2.083 --- more than four times the truth, and outside the front-door interval by a factor of four. The formula recovers the effect even though $U$ is never measured; the naive comparison is picking up the path through $U$, which enters the outcome with coefficient 2 and the treatment through the logit. The bias here is far larger than in the front-door example of the [graphs chapter](graphs-identification-estimation.qmd), where the naive contrast was 0.14 against a truth of 0.055. The difference is the strength of the confounding: $U$ enters $Y$ with coefficient 2 here against 0.7 there. Nothing about the front-door formula changes; what changes is how much work it is doing. ## When identification fails Not every graph identifies the effect. If $A$ and $Y$ have a hidden common cause and there is no front-door variable, the effect is not identified. The graph below is the **bow**: a directed edge $A \to Y$ plus a bidirected $A \leftrightarrow Y$, and nothing else. Remove the mediator $M$ from the front-door graph and this is what is left. ```{r} #| label: not-identified #| cache: true # A → Y with unmeasured confounding (A <-> Y) and no mediator: the "bow" # graph needs THREE igraph edges — the directed A → Y plus the # bidirected pair (encoded as two opposite edges marked "U"). # graph_from_literal() silently collapses a duplicated A -+ Y, so we # build the multigraph with make_graph(). g_unident_ig <- make_graph(c("A","Y", "A","Y", "Y","A"), directed = TRUE) # make_graph() preserves the literal edge order given above (unlike # graph_from_literal(), which re-sorts), so index 2:3 is reliable as long as # the c(...) vector above isn't reordered -- mark by endpoint pair anyway so # a future edit to that vector can't silently mismark the real A->Y edge. edge_ends <- ends(g_unident_ig, E(g_unident_ig)) # The first A->Y edge is the real causal edge; the second (a duplicate, # per `duplicated()`) together with Y->A forms the bidirected U pair. u_edges <- which(((edge_ends[,1]=="A" & edge_ends[,2]=="Y") & duplicated(edge_ends)) | (edge_ends[,1]=="Y" & edge_ends[,2]=="A")) g_unident_ig <- set_edge_attr(g_unident_ig, "description", index = u_edges, value = "U") result <- tryCatch( causal.effect(y = "Y", x = "A", G = g_unident_ig, simp = TRUE), error = function(e) conditionMessage(e) ) cat("Result from causal.effect:\n", result, "\n") ``` `causaleffect` raises an error, printed above as "Not identifiable.", which the `tryCatch` captures. This is a useful failure. The algorithm did not return a formula that would quietly estimate the wrong thing; it reported that no function of the observed distribution equals the causal effect under this graph. No amount of data helps, because the obstacle is the graph and not the sample size. The honest response in applied work is to report the non-identifiability, then either: 1. Find additional variables that close the unmeasured-confounding paths 2. Apply sensitivity analysis to bound the effect under maintained assumptions 3. Find an instrumental variable for a LATE-style identification (IV chapter) Saying that the hidden confounding is probably small is sensitivity analysis, not identification. ## End-to-end: discover, identify, estimate Discovery algorithms return graphs that can be passed to the identification step: ```{r} #| eval: false #| echo: true library(pcalg) # Step 1: discover the graph (chapters on causal discovery) suffStat <- list(C = cor(data_obs), n = nrow(data_obs)) pc_fit <- pc(suffStat, indepTest = gaussCItest, alpha = 0.05, p = ncol(data_obs)) # Step 2: check for backdoor adjustment set using dagitty # (convert pcalg CPDAG → dagitty for convenience). # NOTE the t(): as(graphNEL, "matrix") is from→to, but pcalg2dagitty # expects the amat.cpdag coding, which is its transpose — without it # every directed edge is silently reversed. g_dag <- pcalg2dagitty(t(as(pc_fit@graph, "matrix")), colnames(data_obs), type = "cpdag") sets <- adjustmentSets(g_dag, exposure = "A", outcome = "Y") # Step 3: if an adjustment set is valid for every DAG in the class, estimate # it with AIPW # Step 4: if the CPDAG is still unresolved, do NOT hand it to # causaleffect::causal.effect() -- that function expects a specified # DAG/ADMG, and different orientations in the class can imply # different effects, different valid adjustment sets, or a different # identification status. Instead either # (a) orient the remaining edges with background knowledge and then # run ordinary identification, or # (b) report what the whole class permits: pcalg::ida() for # linear-Gaussian total effects, otherwise enumerate the # compatible DAGs and keep only conclusions that hold across # all of them. # Step 5: estimate the returned functional by plug-in or regression ``` Step 4 is where this workflow most often goes wrong, so it is worth stating plainly: a CPDAG is an equivalence class, not a graph. An empty adjustment-set result for a CPDAG is not permission to pass one arbitrary orientation to the ID algorithm. The [causal discovery chapter](causal-discovery.qmd) sets out what to do instead. In applied work, I would not trust the discovered graph blindly. A more reasonable workflow is: 1. **Hypothesise** a graph based on subject-matter knowledge 2. **Discover** a data-driven graph using PC/GES/FCI 3. **Compare** the two — disagreement points to substantive questions 4. **Identify** under both graphs separately 5. **Estimate** under both, and report the range as part of the result If the estimate changes a lot across plausible graphs, the data alone is not settling the question. ## What changes after identification Backdoor adjustment is not always available because many realistic DAGs have hidden confounders. Front-door, nested, and general ID algorithms may still identify the effect, although their formulas quickly become difficult to derive by hand. The identification result also determines the estimator. AIPW for a backdoor adjustment uses a different formula from a front-door plug-in estimator. Keeping these steps together helps us avoid applying an estimator that does not match the graph. It also forces us to state which variables support an unconfoundedness claim instead of leaving the claim implicit. ## Summary - A graph is not an estimate. Identification is the bridge. - `dagitty` handles backdoor and front-door checks; `causaleffect` handles general ADMG identification. - Backdoor effects can be estimated with AIPW. - Front-door effects require estimating the front-door functional. - If the effect is not identified, report that instead of forcing a biased regression estimate.