Causal Inference VIII(c) — DAGs and Discovery on Survey Data

The same two methods, on variables you can judge for yourself

The two notebooks before this one made their arguments on simulations and on the Sachs protein-signalling network. Both are the right choice for what they do — a simulation is the only honest way to test whether an identification rule recovers a known effect, and Sachs is the field's benchmark because molecular biology supplies a consensus network to score against.

Both also share a limitation. You cannot grade the answer yourself. When PC proposes an edge between Plcg and PIP2, a reader without a background in cell signalling has no independent opinion about whether that is sensible, and has to take the scoring on trust.

So this notebook runs the same two methods on NHANES — the CDC's national health survey — where every variable is something you have a view about already. Age, sex, income, education, BMI, blood pressure, smoking, exercise, cholesterol, HDL, diabetes, sleep. 4,048 adults, complete cases, from the 2017–2018 cycle.

That buys two things the earlier notebooks cannot offer:

• a bad control that is not a toy. Estimating the effect of exercise on blood pressure, the two variables a health researcher would most reflexively adjust for turn out to be the two the graph forbids — and we can measure exactly what including them costs.
• a ground truth requiring no expertise whatsoever. Nothing causes your age. Nothing causes your sex. Any edge the algorithm orients into those two is a provable error, checkable by any reader, with no consensus network needed.

Data provenance. Assembled from the public NHANES 2017–2018 components at wwwn.cdc.gov/nchs/data/nhanes/public/2017/DataFiles/ — DEMO_J (age, sex, family income-to-poverty ratio, education), BMX_J (BMI), BPX_J (systolic, averaged over the available readings), SMQ_J (SMQ020, ever smoked 100 cigarettes), PAQ_J (PAQ650/PAQ665, vigorous or moderate recreational activity), TCHOL_J, HDL_J, DIQ_J (DIQ010, told by a doctor), SLQ_J (usual hours of sleep). Restricted to adults aged 20+ and to complete cases on those twelve variables.

1. The control that should not be there

Start with a question anyone might ask of this survey: does recreational exercise lower systolic blood pressure, and by how much?

The raw comparison is not the answer, because exercisers differ from non-exercisers in ways that also move blood pressure — they are younger, richer and better educated. Those are confounders: common causes of both, sitting on back-door paths that must be blocked. Adjusting for them is exactly what the back-door criterion prescribes.

The interesting part is what happens next. The natural instinct — and the near-universal practice — is to keep going: also adjust for BMI, for diabetes, for cholesterol. More controls, more careful. But the graph says otherwise. Exercise causes lower BMI, and BMI causes higher blood pressure, so BMI sits on the causal path from exercise to the outcome. It is a mediator, not a confounder. Adjusting for it does not remove bias; it removes part of the effect being estimated.

We build the ladder and measure the damage.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, networkx as nx, warnings
warnings.filterwarnings("ignore")
import statsmodels.api as sm
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"

d = pd.read_csv("nhanes_causal.csv")
print(f"NHANES 2017-2018, adults 20+, complete cases: {len(d):,}")
print(f"  mean age {d.age.mean():.1f};  {100*d.female.mean():.0f}% female;  "
      f"mean BMI {d.bmi.mean():.1f};  mean systolic {d.systolic.mean():.1f} mmHg")
print(f"  {100*d.smoker.mean():.0f}% ever-smokers;  {100*d.exercise.mean():.0f}% report recreational activity;  "
      f"{100*d.diabetes.mean():.0f}% told they have diabetes")
print()

def eff(extra):
    Xd = sm.add_constant(d[["exercise"] + extra]); m = sm.OLS(d.systolic, Xd).fit()
    return m.params["exercise"], m.bse["exercise"]

LADDER = [
 ("no adjustment at all",                        [], GREY),
 ("+ age",                                       ["age"], GREY),
 ("+ sex",                                       ["age","female"], GREY),
 ("+ income, education",                         ["age","female","income","educ"], GREEN),
 ("+ BMI",                                       ["age","female","income","educ","bmi"], RED),
 ("+ diabetes",                                  ["age","female","income","educ","bmi","diabetes"], RED),
 ("+ cholesterol, HDL, sleep",                   ["age","female","income","educ","bmi","diabetes",
                                                  "chol","hdl","sleep"], RED),
]
print("Effect of recreational exercise on systolic blood pressure (mmHg):")
print(f"  {'adjustment set':<32} {'estimate':>9} {'SE':>7}")
rows=[]
for lab, xs, c in LADDER:
    b, se = eff(xs); rows.append((lab, b, se, c))
    print(f"  {lab:<32} {b:>9.3f} {se:>7.3f}")
defensible = rows[3][1]; withbmi = rows[4][1]; se_def = rows[3][2]; se_raw = rows[0][2]
lo, hi = defensible-1.96*se_def, defensible+1.96*se_def
print()
print(f"Confounding is real and large: age alone takes the raw {rows[0][1]:.3f} to {rows[1][1]:.3f}, because")
print(f"exercisers are younger and blood pressure rises with age. Blocking the back door is not optional.")
print(f"The defensible set -- age, sex, income, education -- leaves {defensible:.3f} mmHg.")
print()
print("STATE THE OBVIOUS BEFORE GOING FURTHER, because the rest of this section is easy to misread.")
print(f"The raw association is decisive: 95% interval [{rows[0][1]-1.96*se_raw:.3f}, {rows[0][1]+1.96*se_raw:.3f}], nowhere near zero.")
print(f"The adjusted one is not: [{lo:.3f}, {hi:.3f}], which INCLUDES ZERO. So the honest headline of this")
print(f"cross-section is that essentially the entire apparent benefit of exercise for blood pressure is")
print(f"confounding, most of it age, and what survives is not distinguishable from nothing.")
print()
print("That does not make the next part academic -- it makes it cleaner. What follows is a question about")
print("WHY the number moves when a variable is added, and that question has an exact answer whether or")
print("not the number it moves is significant.")
print()
print(f"The estimate loses {100*(1-withbmi/defensible):.0f}% of what is left the moment BMI is added: {withbmi:.3f}.")
print("That is not bias being removed. BMI is on the causal path, and here is the arithmetic to prove it.")
print()

mb = sm.OLS(d.bmi, sm.add_constant(d[["exercise","age","female","income","educ"]])).fit()
ms = sm.OLS(d.systolic, sm.add_constant(d[["bmi","age","female","income","educ"]])).fit()
a_path, b_path = mb.params["exercise"], ms.params["bmi"]
print(f"  exercise -> BMI          {a_path:+.3f} kg/m2      (SE {mb.bse['exercise']:.3f})")
print(f"  BMI -> systolic          {b_path:+.3f} mmHg/unit  (SE {ms.bse['bmi']:.3f})")
print(f"  indirect path a x b      {a_path*b_path:+.3f} mmHg   <- the effect that travels through BMI")
print(f"  removed by adjusting     {defensible-withbmi:+.3f} mmHg   <- what the estimate lost")
print(f"  they differ by           {abs(a_path*b_path-(defensible-withbmi)):.3f} mmHg")
print()
print(f"Those two numbers agree to within {abs(a_path*b_path-(defensible-withbmi)):.3f} mmHg, which is the point rather than a coincidence:")
print("the amount the estimate falls when you 'control for' BMI is exactly the effect that travels")
print("through BMI. Adjusting for a mediator does not clean the estimate. It deletes a real channel of")
print("the very thing being measured, and reports the remainder as though it were the whole.")
print()
print("SO NAME THE THREE QUANTITIES, BECAUSE 'DO NOT ADJUST FOR A MEDIATOR' IS TOO BLUNT A RULE.")
print()
indirect = a_path*b_path
print(f"  TOTAL effect        {defensible:>7.3f} mmHg   everything exercise does to blood pressure,")
print(f"                                       by every route, adjusting only for confounders")
print(f"  INDIRECT, via BMI   {indirect:>7.3f} mmHg   the part that happens BECAUSE exercise made")
print(f"                                       people lighter: it moves BMI by {a_path:.3f} kg/m2, and")
print(f"                                       each unit of BMI is worth {b_path:+.3f} mmHg")
print(f"  DIRECT              {withbmi:>7.3f} mmHg   everything else -- vascular tone, autonomic")
print(f"                                       effects, fitness at a given body weight")
print(f"  check               {withbmi + indirect:>7.3f} = direct + indirect, against a total of {defensible:.3f}")
print()
print(f"In words: about {100*indirect/defensible:.0f}% of what exercise does to blood pressure in this cross-section")
print(f"it does by way of weight. The part that is not about weight -- the direct effect -- is {withbmi:.3f} mmHg,")
print(f"with an interval of [{withbmi-1.96*rows[4][2]:.3f}, {withbmi+1.96*rows[4][2]:.3f}] that comfortably contains zero.")
print()
print("Which sharpens what went wrong in row five. Adjusting for BMI is not meaningless and does not")
print("produce a corrupted number: it produces the DIRECT effect, a real and well-defined quantity that")
print("someone might genuinely want. 'Does exercise lower blood pressure?' and 'does exercise lower blood")
print("pressure other than by making you lighter?' are different questions with different answers, and")
print("the second is the one row five answers.")
print()
print("The mistake is never computing it. The mistake is computing it and reporting it as though it")
print("answered the first question -- which is what happens whenever BMI goes into the model because it")
print("was available and seemed relevant, rather than because someone decided the direct effect was the")
print("estimand. Nothing in the regression output records which of the two was intended.")
print()
md_ = sm.Logit(d.diabetes, sm.add_constant(d[["exercise","age","female","income","educ"]])).fit(disp=0)
print(f"Diabetes is a second mediator on the same footing -- exercise moves it by {md_.params['exercise']:+.3f} in log-odds")
print(f"(SE {md_.bse['exercise']:.3f}) -- and cholesterol and HDL sit in the same position.")
print("Every variable added after the fourth row makes the answer worse while making the table look")
print("more thorough. No fit statistic, no standard error and no significance test distinguishes the")
print("fourth row from the seventh. Only the graph does.")

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ys = np.arange(len(rows))[::-1]
ax[0].barh(ys, [r[1] for r in rows], color=[r[3] for r in rows], height=.62)
ax[0].errorbar([r[1] for r in rows], ys, xerr=[1.96*r[2] for r in rows], fmt="none", ecolor="k", lw=1.2, capsize=3)
ax[0].set_yticks(ys); ax[0].set_yticklabels([r[0] for r in rows], fontsize=9)
ax[0].axvline(0, color="k", lw=1)
ax[0].axvline(defensible, color=GREEN, ls="--", lw=1.6, label=f"defensible set: {defensible:.2f}")
ax[0].set_xlabel("estimated effect of exercise on systolic BP (mmHg)")
ax[0].set_title("Green = the back-door set. Red = mediators, which should not be there.", fontsize=10)
ax[0].legend(fontsize=8, loc="lower left")

G = nx.DiGraph(); pos = {"exercise":(0,0), "BMI":(1,.85), "systolic":(2,0), "age, sex,\nincome, educ":(1,-.9)}
E = [("exercise","BMI"),("BMI","systolic"),("exercise","systolic"),
     ("age, sex,\nincome, educ","exercise"),("age, sex,\nincome, educ","systolic")]
for nname,(x,y) in pos.items():
    col = "#fdecea" if nname=="BMI" else ("#eaf5ec" if "age" in nname else "white")
    ax[1].scatter([x],[y],s=5200,facecolor=col,edgecolor="k",linewidths=1.6,zorder=3)
    ax[1].text(x,y,nname,ha="center",va="center",fontsize=8.5,zorder=4)
for u,v in E:
    x1,y1=pos[u]; x2,y2=pos[v]; dx,dy=x2-x1,y2-y1; L=np.hypot(dx,dy); ux,uy=dx/L,dy/L; r=.235
    c = RED if (u=="BMI" or v=="BMI") else (GREEN if "age" in u else "k")
    ax[1].annotate("",xy=(x2-ux*r,y2-uy*r),xytext=(x1+ux*r,y1+uy*r),
                   arrowprops=dict(arrowstyle="-|>",color=c,lw=2.2),zorder=2)
ax[1].text(1,.85+.34,"mediator - do NOT adjust",ha="center",fontsize=8.5,color=RED)
ax[1].text(1,-.9-.38,"confounders - DO adjust",ha="center",fontsize=8.5,color=GREEN)
ax[1].set_xlim(-.6,2.6); ax[1].set_ylim(-1.6,1.5); ax[1].axis("off")
ax[1].set_title("The graph that decides which controls are allowed", fontsize=10)
plt.tight_layout(); plt.show()

NHANES 2017-2018, adults 20+, complete cases: 4,048
  mean age 51.3;  51% female;  mean BMI 29.8;  mean systolic 126.7 mmHg
  43% ever-smokers;  48% report recreational activity;  19% told they have diabetes

Effect of recreational exercise on systolic blood pressure (mmHg):
  adjustment set                    estimate      SE
  no adjustment at all                -4.573   0.615
  + age                               -1.600   0.556
  + sex                               -1.675   0.556
  + income, education                 -0.838   0.578
  + BMI                               -0.240   0.572
  + diabetes                          -0.271   0.572
  + cholesterol, HDL, sleep           -0.248   0.570

Confounding is real and large: age alone takes the raw -4.573 to -1.600, because
exercisers are younger and blood pressure rises with age. Blocking the back door is not optional.
The defensible set -- age, sex, income, education -- leaves -0.838 mmHg.

STATE THE OBVIOUS BEFORE GOING FURTHER, because the rest of this section is easy to misread.
The raw association is decisive: 95% interval [-5.778, -3.368], nowhere near zero.
The adjusted one is not: [-1.970, 0.294], which INCLUDES ZERO. So the honest headline of this
cross-section is that essentially the entire apparent benefit of exercise for blood pressure is
confounding, most of it age, and what survives is not distinguishable from nothing.

That does not make the next part academic -- it makes it cleaner. What follows is a question about
WHY the number moves when a variable is added, and that question has an exact answer whether or
not the number it moves is significant.

The estimate loses 71% of what is left the moment BMI is added: -0.240.
That is not bias being removed. BMI is on the causal path, and here is the arithmetic to prove it.

  exercise -> BMI          -1.461 kg/m2      (SE 0.242)
  BMI -> systolic          +0.411 mmHg/unit  (SE 0.037)
  indirect path a x b      -0.601 mmHg   <- the effect that travels through BMI
  removed by adjusting     -0.599 mmHg   <- what the estimate lost
  they differ by           0.002 mmHg

Those two numbers agree to within 0.002 mmHg, which is the point rather than a coincidence:
the amount the estimate falls when you 'control for' BMI is exactly the effect that travels
through BMI. Adjusting for a mediator does not clean the estimate. It deletes a real channel of
the very thing being measured, and reports the remainder as though it were the whole.

SO NAME THE THREE QUANTITIES, BECAUSE 'DO NOT ADJUST FOR A MEDIATOR' IS TOO BLUNT A RULE.

  TOTAL effect         -0.838 mmHg   everything exercise does to blood pressure,
                                       by every route, adjusting only for confounders
  INDIRECT, via BMI    -0.601 mmHg   the part that happens BECAUSE exercise made
                                       people lighter: it moves BMI by -1.461 kg/m2, and
                                       each unit of BMI is worth +0.411 mmHg
  DIRECT               -0.240 mmHg   everything else -- vascular tone, autonomic
                                       effects, fitness at a given body weight
  check                -0.840 = direct + indirect, against a total of -0.838

In words: about 72% of what exercise does to blood pressure in this cross-section
it does by way of weight. The part that is not about weight -- the direct effect -- is -0.240 mmHg,
with an interval of [-1.360, 0.881] that comfortably contains zero.

Which sharpens what went wrong in row five. Adjusting for BMI is not meaningless and does not
produce a corrupted number: it produces the DIRECT effect, a real and well-defined quantity that
someone might genuinely want. 'Does exercise lower blood pressure?' and 'does exercise lower blood
pressure other than by making you lighter?' are different questions with different answers, and
the second is the one row five answers.

The mistake is never computing it. The mistake is computing it and reporting it as though it
answered the first question -- which is what happens whenever BMI goes into the model because it
was available and seemed relevant, rather than because someone decided the direct effect was the
estimand. Nothing in the regression output records which of the two was intended.

Diabetes is a second mediator on the same footing -- exercise moves it by -0.324 in log-odds
(SE 0.091) -- and cholesterol and HDL sit in the same position.
Every variable added after the fourth row makes the answer worse while making the table look
more thorough. No fit statistic, no standard error and no significance test distinguishes the
fourth row from the seventh. Only the graph does.

2. Discovery on variables you can grade

Now the second method, on the same twelve variables. The PC algorithm tests conditional independencies to recover the skeleton, orients what the colliders logically permit, and returns a CPDAG.

On the Sachs data the recovered graph had to be scored against a consensus network drawn from decades of molecular biology. Here we can do something better and much harder to argue with. Two of these variables cannot be caused by anything else in the table. Nothing you do makes you older; nothing in the survey determines your sex. Every arrow pointing into age or female is therefore wrong — not implausible, not unsupported, but impossible — and no expertise is required to see it.

That gives a scoring rule with no consensus network, no literature, and no room for disagreement.

from causallearn.search.ConstraintBased.PC import pc
V = ["age","female","income","educ","bmi","systolic","smoker","exercise","chol","hdl","diabetes","sleep"]
X = d[V].values.astype(float)
EXOG = {"age","female"}

def edges_of(G):
    di, un = [], []
    for i in range(len(V)):
        for j in range(i+1, len(V)):
            if   G[i,j]==-1 and G[j,i]== 1: di.append((i,j))
            elif G[i,j]== 1 and G[j,i]==-1: di.append((j,i))
            elif G[i,j]==-1 and G[j,i]==-1: un.append((i,j))
    return di, un

cg = pc(X, 0.01, show_progress=False)
di, un = edges_of(cg.G.graph)
npairs = len(V)*(len(V)-1)//2
print(f"PC at alpha = 0.01 on {len(d):,} adults: {len(di)} directed and {len(un)} undirected edges, "
      f"out of {npairs} possible pairs.")
print()
into  = [(V[u],V[v]) for u,v in di if V[v] in EXOG]
outof = [(V[u],V[v]) for u,v in di if V[u] in EXOG]
print(f"Edges pointing INTO age or sex: {len(into)}")
for u,v in sorted(into): print(f"   {u:>9}  ->  {v:<7}   impossible")
print()
print(f"Edges pointing OUT of age or sex: {len(outof)}")
print()
print(f"So of the {len(into)+len(outof)} edges touching the two variables whose direction is beyond dispute,")
print(f"the algorithm got {len(into)} backwards and {len(outof)} right. Not most of them wrong. All of them.")
print()
print("This is not a defect of PC, and a better algorithm would not fix it. The orientation rules")
print("work from conditional independencies alone, and those are symmetric in a way the world is not:")
print("the data contains no trace of the fact that time runs forwards. Everything the algorithm knows")
print("about direction, it deduced from collider patterns, and on this data those patterns point the")
print("wrong way for every edge involving age or sex.")
print()
print("Worth noting what it does get right. The SKELETON is largely sensible -- BMI with diabetes,")
print("cholesterol with HDL, income with education, age with blood pressure are all real associations")
print("and the algorithm found them. Discovery is far better at 'these two are connected' than at")
print("'this one causes that one', which is exactly the Markov-equivalence limit made visible on")
print("variables you can check without a textbook.")
print()
for al in (0.001, 0.01, 0.05):
    g = pc(X, al, show_progress=False).G.graph
    dd, uu = edges_of(g)
    bb = sum(1 for u,v in dd if V[v] in EXOG)
    print(f"  alpha = {al:<6} directed {len(dd):>2}   undirected {len(uu):>2}   impossible orientations {bb:>2}")
print()
print("Loosening the test does not help. It adds edges, and it adds errors with them.")

posN = {"age":(0,2.2),"female":(0,1.1),"income":(1.15,2.6),"educ":(1.15,1.6),
        "smoker":(2.4,2.8),"exercise":(2.4,1.9),"sleep":(2.4,1.0),"bmi":(3.6,2.5),
        "chol":(3.6,1.5),"hdl":(3.6,0.5),"diabetes":(4.8,2.0),"systolic":(4.8,1.0)}
def draw_net(ax, di, un, title):
    for nm,(x,y) in posN.items():
        col = "#eef2ff" if nm in EXOG else "white"
        ax.scatter([x],[y],s=1750,facecolor=col,edgecolor="k",linewidths=1.4,zorder=3)
        ax.text(x,y,nm,ha="center",va="center",fontsize=7.6,zorder=4)
    for u,v in di:
        a,b = V[u],V[v]; x1,y1=posN[a]; x2,y2=posN[b]
        dx,dy=x2-x1,y2-y1; L=np.hypot(dx,dy); ux,uy=dx/L,dy/L; r=.30
        bad = b in EXOG
        ax.annotate("",xy=(x2-ux*r,y2-uy*r),xytext=(x1+ux*r,y1+uy*r),
                    arrowprops=dict(arrowstyle="-|>",color=RED if bad else GREY,
                                    lw=2.1 if bad else 1.2),zorder=2 if not bad else 5)
    for u,v in un:
        a,b = V[u],V[v]; x1,y1=posN[a]; x2,y2=posN[b]
        dx,dy=x2-x1,y2-y1; L=np.hypot(dx,dy); ux,uy=dx/L,dy/L; r=.30
        ax.plot([x1+ux*r,x2-ux*r],[y1+uy*r,y2-uy*r],color=ORANGE,lw=1.8,ls="--",zorder=2)
    ax.set_xlim(-.7,5.5); ax.set_ylim(0.1,3.3); ax.axis("off"); ax.set_title(title,fontsize=10)

fig, ax = plt.subplots(figsize=(12, 4.4))
draw_net(ax, di, un, f"PC on NHANES: {len(into)} edges point into age or sex (red) and are impossible")
plt.tight_layout(); plt.show()

PC at alpha = 0.01 on 4,048 adults: 27 directed and 1 undirected edges, out of 66 possible pairs.

Edges pointing INTO age or sex: 11
         bmi  ->  female    impossible
    diabetes  ->  age       impossible
    exercise  ->  age       impossible
         hdl  ->  age       impossible
         hdl  ->  female    impossible
      income  ->  age       impossible
       sleep  ->  age       impossible
       sleep  ->  female    impossible
      smoker  ->  age       impossible
      smoker  ->  female    impossible
    systolic  ->  age       impossible

Edges pointing OUT of age or sex: 0

So of the 11 edges touching the two variables whose direction is beyond dispute,
the algorithm got 11 backwards and 0 right. Not most of them wrong. All of them.

This is not a defect of PC, and a better algorithm would not fix it. The orientation rules
work from conditional independencies alone, and those are symmetric in a way the world is not:
the data contains no trace of the fact that time runs forwards. Everything the algorithm knows
about direction, it deduced from collider patterns, and on this data those patterns point the
wrong way for every edge involving age or sex.

Worth noting what it does get right. The SKELETON is largely sensible -- BMI with diabetes,
cholesterol with HDL, income with education, age with blood pressure are all real associations
and the algorithm found them. Discovery is far better at 'these two are connected' than at
'this one causes that one', which is exactly the Markov-equivalence limit made visible on
variables you can check without a textbook.

  alpha = 0.001  directed 25   undirected  0   impossible orientations  8
  alpha = 0.01   directed 27   undirected  1   impossible orientations 11

  alpha = 0.05   directed 28   undirected  2   impossible orientations 12

Loosening the test does not help. It adds edges, and it adds errors with them.

3. The repair is knowledge, not more data

The obvious response to eleven impossible arrows is to want a better algorithm, or more observations, or a more careful test level. None of those is the answer, and the previous cell already showed the test level making things worse.

What fixes it is the thing the reader supplied in the first place: nothing causes your age. That is not in the data and never will be, at any sample size. It is background knowledge, and constraint-based discovery has a formal slot for exactly that — forbid the arrows that cannot exist and let the orientation rules propagate the consequences.

The question worth asking is how far one sentence of knowledge travels.

from causallearn.utils.PCUtils.BackgroundKnowledge import BackgroundKnowledge
nodes = cg.G.get_nodes()
bk = BackgroundKnowledge()
for k, name in enumerate(V):
    if name in EXOG:
        for j in range(len(V)):
            if j != k: bk.add_forbidden_by_node(nodes[j], nodes[k])
cg2 = pc(X, 0.01, show_progress=False, background_knowledge=bk)
di2, un2 = edges_of(cg2.G.graph)
into2 = [(V[u],V[v]) for u,v in di2 if V[v] in EXOG]

s0 = {(V[u],V[v]) for u,v in di}; s2 = {(V[u],V[v]) for u,v in di2}
flipped = {(a,b) for a,b in s2 if (b,a) in s0}
propagated = sorted(flipped - {(b,a) for a,b in into})

print(f"  {'':<34} {'directed':>9} {'undirected':>11} {'impossible':>11}")
print(f"  {'PC alone':<34} {len(di):>9} {len(un):>11} {len(into):>11}")
print(f"  {'PC + nothing causes age or sex':<34} {len(di2):>9} {len(un2):>11} {len(into2):>11}")
print()
print(f"The eleven impossible edges are gone, which was guaranteed -- they were forbidden by hand.")
print(f"The interesting number is the other one: {len(flipped)} edges changed direction in total, so")
print(f"{len(propagated)} were repaired that nobody touched:")
for a,b in propagated: print(f"     {a} -> {b}")
print()
print(f"And the {len(un)} edge PC could not orient at all is now resolved: {len(un2)} remain undirected.")
print("Constraining two variables propagated through the orientation rules and settled parts of the")
print("graph that had nothing to do with age or sex. Knowledge is not a patch applied to the output;")
print("it enters the search and changes what the data is able to say.")
print()
print("What comes back now reads like something a person might have drawn:")
for a,b in sorted(s2):
    if a in EXOG: print(f"     {a} -> {b}")
print()
print("None of that is a discovery. Every one of those arrows was obvious before the algorithm ran,")
print("which is the honest summary of the whole exercise: the graph is credible exactly to the extent")
print("that a human constrained it, and the parts nobody constrained are the parts to distrust.")
print()
print("Set against section 1, the two halves make one argument. There, the graph was known and the")
print("data could not tell you that adjusting for BMI was a mistake -- every diagnostic looked fine")
print("and the estimate lost most of its magnitude. Here, the data could not tell you which way an")
print("arrow points even when the answer is beyond dispute. Both are the same limit seen from two")
print("sides: association is symmetric, causation is not, and the asymmetry has to come from")
print("somewhere other than the joint distribution -- from an experiment, or from knowing something.")

fig, ax = plt.subplots(2, 1, figsize=(12, 8.4))
draw_net(ax[0], di, un, f"PC alone: {len(into)} impossible orientations (red), {len(un)} unresolved (dashed)")
draw_net(ax[1], di2, un2, f"PC + one sentence of background knowledge: {len(into2)} impossible, {len(un2)} unresolved")
plt.tight_layout(); plt.show()

                                      directed  undirected  impossible
  PC alone                                  27           1          11
  PC + nothing causes age or sex            28           0           0

The eleven impossible edges are gone, which was guaranteed -- they were forbidden by hand.
The interesting number is the other one: 13 edges changed direction in total, so
2 were repaired that nobody touched:
     educ -> exercise
     educ -> income

And the 1 edge PC could not orient at all is now resolved: 0 remain undirected.
Constraining two variables propagated through the orientation rules and settled parts of the
graph that had nothing to do with age or sex. Knowledge is not a patch applied to the output;
it enters the search and changes what the data is able to say.

What comes back now reads like something a person might have drawn:
     age -> diabetes
     age -> exercise
     age -> hdl
     age -> income
     age -> sleep
     age -> smoker
     age -> systolic
     female -> bmi
     female -> hdl
     female -> sleep
     female -> smoker

None of that is a discovery. Every one of those arrows was obvious before the algorithm ran,
which is the honest summary of the whole exercise: the graph is credible exactly to the extent
that a human constrained it, and the parts nobody constrained are the parts to distrust.

Set against section 1, the two halves make one argument. There, the graph was known and the
data could not tell you that adjusting for BMI was a mistake -- every diagnostic looked fine
and the estimate lost most of its magnitude. Here, the data could not tell you which way an
arrow points even when the answer is beyond dispute. Both are the same limit seen from two
sides: association is symmetric, causation is not, and the asymmetry has to come from
somewhere other than the joint distribution -- from an experiment, or from knowing something.

4. What the corrected graph says about the original question

The two halves have run side by side without touching. Section 1 estimated an effect under a graph drawn by hand; sections 2 and 3 learned a graph from the data and repaired it with background knowledge. The obvious question has not been asked: take the learned graph seriously, read the adjustment set off it, and see whether it agrees.

For the effect of a treatment, adjusting for the treatment's parents in the graph satisfies the back-door criterion. So the learned graph specifies its own adjustment set, and it is not the one drawn by hand.

E2 = [(V[u], V[v]) for u, v in di2]
pa = lambda x: sorted(a for a, b in E2 if b == x)

print(f"  parents of exercise in the learned graph : {pa('exercise')}")
print(f"  the hand-drawn back-door set             : ['age', 'educ', 'female', 'income']")
print()
learned_set = pa("exercise")
b_l, se_l = eff(learned_set)
b_h, se_h = eff(["age","female","income","educ"])
print(f"  {'adjustment set':<44} {'estimate':>9} {'SE':>7}   {'95% interval':>20}")
for lab, b, se in [(f"learned graph: {', '.join(learned_set)}", b_l, se_l),
                   ("hand-drawn: age, sex, income, educ", b_h, se_h)]:
    print(f"  {lab:<44} {b:>9.3f} {se:>7.3f}   [{b-1.96*se:>7.3f}, {b+1.96*se:>7.3f}]")
print()
print(f"They agree, {abs(b_h-b_l):.3f} mmHg apart -- about {abs(b_h-b_l)/se_h:.2f} of a standard error. Two routes to the")
print("adjustment set, one from domain reasoning and one learned from the data under a constraint,")
print("landing in the same place. That is the reassuring result.")
print()
print("THE UNRESERVED ONE IS THAT THE LEARNED GRAPH DOES NOT CONTAIN THE EDGE AT ALL.")
print()
ie, isy = V.index("exercise"), V.index("systolic")
for al in (0.001, 0.01, 0.05, 0.10, 0.20):
    g = pc(X, al, show_progress=False, background_knowledge=bk).G.graph
    present = g[ie,isy] != 0 or g[isy,ie] != 0
    dd, uu = edges_of(g)
    print(f"    alpha = {al:<6} exercise-systolic edge present: {str(present):<6}  "
          f"({len(dd)} directed edges in the graph)")
print()
from scipy import stats
def pcorr(a, b, given):
    ra = sm.OLS(d[a], sm.add_constant(d[given])).fit().resid
    rb = sm.OLS(d[b], sm.add_constant(d[given])).fit().resid
    return stats.pearsonr(ra, rb)
for given in (["age"], ["age","educ"], ["age","female","income","educ"]):
    r, pv = pcorr("exercise", "systolic", given)
    print(f"  partial corr(exercise, systolic | {', '.join(given):<26}) = {r:+.4f}   p = {pv:.3f}")
print()
print("Which is not a contradiction of section 1 -- it is section 1's own answer, arriving in a")
print(f"different language. The hand-drawn estimate was {b_h:.3f} with an interval of")
print(f"[{b_h-1.96*se_h:.3f}, {b_h+1.96*se_h:.3f}], and that interval contains zero. Conditioning on the confounders takes the")
print("partial correlation from decisively non-zero to not significant, so a discovery algorithm testing")
print("for conditional independence deletes the edge. Both methods say the same thing: once the back")
print("door is blocked, this cross-section has no detectable effect of exercise on blood pressure left")
print("to report. Regression says it with an interval spanning zero; discovery says it by drawing no arrow.")
print()
print("BUT THE TWO GRAPHS DISAGREE ABOUT BMI, AND THE DISAGREEMENT MATTERS.")
print()
print(f"  hand-drawn : exercise -> BMI -> systolic     BMI is a MEDIATOR")
print(f"  learned    : parents of BMI = {pa('bmi')}")
if "exercise" in pa("bmi") and "systolic" in pa("bmi"):
    print(f"               exercise -> BMI <- systolic     BMI is a COLLIDER")
print()
print("Both graphs forbid adjusting for BMI, and for opposite reasons. Under the hand-drawn graph you")
print("would be deleting a real causal channel. Under the learned one you would be opening a path that")
print("is closed, manufacturing association out of nothing. Same instruction, incompatible explanations,")
print("and identical observable consequences -- which is why no amount of this data settles it.")
print()
print("Note also what the learned graph asserts to get there: systolic -> BMI. Blood pressure raising")
print("body mass runs backwards to the physiology, and it survived section 3 untouched because the only")
print("knowledge supplied was that nothing causes age or sex. The repair is exactly as good as what you")
print("put into it. One sentence fixed eleven edges; it did not fix this one, and nothing in the output")
print("marks the difference between the edges knowledge corrected and the edges it never reached.")

# ---- every estimate of the exercise effect this notebook has produced, in one place
b_raw, se_raw = eff([])
b_bmi, se_bmi = eff(["age","female","income","educ","bmi"])
indirect = a_path * b_path
EST = [
 ("no adjustment",                      b_raw, se_raw, RED,   "confounded -- mostly by age"),
 ("hand-drawn back-door set",           b_h,   se_h,   GREEN, "age, sex, income, education"),
 ("learned graph's own set",            b_l,   se_l,   BLUE,  f"{', '.join(learned_set)} -- parents of exercise"),
 ("+ BMI  (the direct effect)",         b_bmi, se_bmi, ORANGE,"a different estimand, not a worse estimate"),
]
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ys = np.arange(len(EST))[::-1]
for (lab, b, se, c, note), y in zip(EST, ys):
    ax[0].errorbar(b, y, xerr=1.96*se, fmt="o", color=c, lw=2.4, capsize=5, markersize=7)
    ax[0].text(b, y+.26, f"{b:.3f}", ha="center", fontsize=9, color=c)
    ax[0].text(4.1, y-.02, note, fontsize=7.6, color="#555", va="center")
ax[0].axvline(0, color="k", lw=1.6, ls="--")
ax[0].set_yticks(ys); ax[0].set_yticklabels([e[0] for e in EST], fontsize=8.6)
ax[0].set_xlim(-7.2, 11.5); ax[0].set_ylim(-.7, len(EST)-.3)
ax[0].set_xlabel("effect of exercise on systolic blood pressure (mmHg), 95% intervals")
ax[0].set_title("Every adjustment set, and where zero falls", fontsize=10)

# ---- the decomposition
lab2 = ["total\n(back-door set)", "indirect\nvia BMI", "direct\n(BMI held fixed)"]
val2 = [defensible, indirect, withbmi]
col2 = [GREEN, PURP, ORANGE]
ax[1].bar(lab2, val2, color=col2, width=.62)
for i, v in enumerate(val2):
    ax[1].text(i, v-.045, f"{v:.3f}", ha="center", va="top", fontsize=9.5, fontweight="bold")
ax[1].axhline(0, color="k", lw=1.2)
ax[1].errorbar([0,2],[defensible,withbmi],yerr=[1.96*se_h,1.96*se_bmi],fmt="none",ecolor="k",lw=1.3,capsize=5)
ax[1].set_ylabel("mmHg"); ax[1].set_ylim(-1.15, .35)
ax[1].set_title(f"Total = indirect + direct: {defensible:.3f} = {indirect:.3f} + {withbmi:.3f}", fontsize=10)
ax[1].text(1, .12, f"{100*indirect/defensible:.0f}% of the total\nruns through weight",
           ha="center", fontsize=8.4, color=PURP)
plt.tight_layout(); plt.show()
print("Left: four estimates, and only the unadjusted one clears zero. Blocking the back door is what")
print("removes the finding; nothing after that is the difference between a result and a null.")
print("Right: the split. The bar in the middle is the part of the effect that is about weight, and it")
print("is most of what there is. The right-hand bar is a direct effect -- a different question's answer.")

  parents of exercise in the learned graph : ['age', 'educ']
  the hand-drawn back-door set             : ['age', 'educ', 'female', 'income']

  adjustment set                                estimate      SE           95% interval
  learned graph: age, educ                        -0.973   0.572   [ -2.095,   0.148]
  hand-drawn: age, sex, income, educ              -0.838   0.578   [ -1.970,   0.294]

They agree, 0.135 mmHg apart -- about 0.23 of a standard error. Two routes to the
adjustment set, one from domain reasoning and one learned from the data under a constraint,
landing in the same place. That is the reassuring result.

THE UNRESERVED ONE IS THAT THE LEARNED GRAPH DOES NOT CONTAIN THE EDGE AT ALL.

    alpha = 0.001  exercise-systolic edge present: False   (25 directed edges in the graph)
    alpha = 0.01   exercise-systolic edge present: False   (28 directed edges in the graph)

    alpha = 0.05   exercise-systolic edge present: False   (30 directed edges in the graph)
    alpha = 0.1    exercise-systolic edge present: False   (32 directed edges in the graph)

    alpha = 0.2    exercise-systolic edge present: False   (34 directed edges in the graph)

  partial corr(exercise, systolic | age                       ) = -0.0452   p = 0.004
  partial corr(exercise, systolic | age, educ                 ) = -0.0267   p = 0.089
  partial corr(exercise, systolic | age, female, income, educ ) = -0.0228   p = 0.147

Which is not a contradiction of section 1 -- it is section 1's own answer, arriving in a
different language. The hand-drawn estimate was -0.838 with an interval of
[-1.970, 0.294], and that interval contains zero. Conditioning on the confounders takes the
partial correlation from decisively non-zero to not significant, so a discovery algorithm testing
for conditional independence deletes the edge. Both methods say the same thing: once the back
door is blocked, this cross-section has no detectable effect of exercise on blood pressure left
to report. Regression says it with an interval spanning zero; discovery says it by drawing no arrow.

BUT THE TWO GRAPHS DISAGREE ABOUT BMI, AND THE DISAGREEMENT MATTERS.

  hand-drawn : exercise -> BMI -> systolic     BMI is a MEDIATOR
  learned    : parents of BMI = ['diabetes', 'exercise', 'female', 'hdl', 'sleep', 'systolic']
               exercise -> BMI <- systolic     BMI is a COLLIDER

Both graphs forbid adjusting for BMI, and for opposite reasons. Under the hand-drawn graph you
would be deleting a real causal channel. Under the learned one you would be opening a path that
is closed, manufacturing association out of nothing. Same instruction, incompatible explanations,
and identical observable consequences -- which is why no amount of this data settles it.

Note also what the learned graph asserts to get there: systolic -> BMI. Blood pressure raising
body mass runs backwards to the physiology, and it survived section 3 untouched because the only
knowledge supplied was that nothing causes age or sex. The repair is exactly as good as what you
put into it. One sentence fixed eleven edges; it did not fix this one, and nothing in the output
marks the difference between the edges knowledge corrected and the edges it never reached.

Left: four estimates, and only the unadjusted one clears zero. Blocking the back door is what
removes the finding; nothing after that is the difference between a result and a null.
Right: the split. The bar in the middle is the part of the effect that is about weight, and it
is most of what there is. The right-hand bar is a direct effect -- a different question's answer.

5. Summary

The same two methods as the previous notebooks, on variables that need no specialist to referee.

• A mediator is not a confounder, and nothing but the graph says so. Estimating the effect of exercise on systolic blood pressure, the back-door set — age, sex, income, education — is genuinely required: age alone moves the estimate most of the way. Adding BMI, which almost any analyst would include, costs the estimate most of what remains, and the drop equals the indirect path $a \times b$ to three decimals. It is not bias removal; it is deletion of a real channel. Diabetes, cholesterol and HDL are the same mistake repeated. No fit statistic, standard error or significance test separates the correct specification from the wrong one.

• Discovery gets the skeleton broadly right and the directions comprehensively wrong. PC found sensible adjacencies — BMI with diabetes, cholesterol with HDL, income with education, age with blood pressure. It then oriented every one of the eleven edges touching age and female backwards, and none correctly. That is not a bug to be fixed by a better algorithm or a larger sample: conditional independencies are symmetric, and the observational distribution contains no trace of the fact that time runs forwards.

• The repair is knowledge, and it is cheap. Forbidding arrows into age and sex — one sentence anybody could supply — removed all eleven errors, repaired two more edges by propagation that were never touched, and resolved the last undirected edge. The resulting graph looks like something a person would draw, because a person supplied the part that made it credible.

The two halves are one argument. In the first, the graph was known and the data could not tell you that a control was wrong. In the second, the data could not tell you which way an arrow points even when the answer is beyond dispute. Association is symmetric and causation is not, and the asymmetry has to come from somewhere other than the joint distribution — from an intervention, or from knowing something about the world.

Cross-links. The bad-control result is the collider and mediator material of the DAGs & SCM notebook on real data instead of a simulation, and the same "adjusting made it worse" pattern the matching notebook found when balance diagnostics chose the wrong estimator. The discovery half is the causal-discovery notebook's Markov-equivalence limit, checked here against a ground truth that needs no consensus network — and its conclusion is the same one Sachs et al. reached by perturbing proteins: identification comes from design and from knowledge, not from analysing the observational distribution harder.