Discovery on a known network: grading the graph¶
Every other case study handed you the graph. Someone had drawn the arrows, and the only job left was to read an effect off them. That is a fair description of the easy half of causal work. The hard half is the question it skips: where does the graph come from when nobody can draw it?
Discovery learns the structure from the data. You feed it a table and it returns the arrows it can defend, with the undecided ones left marked rather than guessed. The trouble with trying it on real data is that you can never grade the answer. Nobody knows the true graph of a real business, so a discovered graph is a claim you have to take partly on faith.
The fix is the same move from why causality: simulate a world whose truth you wrote yourself, then grade the method against the answer key. Here the world is a textbook Bayesian network. The Asia chest-clinic network (Lauritzen and Spiegelhalter, 1988) has eight variables and eight known edges, so when discovery returns a graph we can score it edge by edge. That is the whole point of running it here. The number you get is not a guess about a method that might work; it is a measurement of one that did.
from ergodic.datasets import load_asia
ds = load_asia() # n=3000 by default; the draw is seeded and deterministic
print(ds.target["source"])
print("nodes:", ds.target["nodes"])
print("true edges:", ds.target["n_edges"])
ds.data
Lauritzen and Spiegelhalter (1988)
nodes: ('asia', 'tub', 'smoke', 'lung', 'bronc', 'either', 'xray', 'dysp')
true edges: 8
TabularData(schema=Schema(8 columns: 8 categorical), n_rows=3000)
Eight binary variables, one row per simulated patient. A recent visit to Asia, tuberculosis,
whether they smoke, lung cancer, bronchitis, a logical "either tuberculosis or lung cancer",
a chest x-ray, and dyspnoea (shortness of breath). The true graph wires them into a small
clinic: smoking drives lung cancer and bronchitis, a visit to Asia raises the tuberculosis
risk, tuberculosis and lung cancer feed the either node, and either plus bronchitis
explain the x-ray and the breathing trouble.
Before discovery, look at the marginals. How common is each condition?
df = ds.data.to_pandas()
rates = (df == "yes").mean().sort_values()
print(rates.round(4))
tub 0.0093 asia 0.0103 lung 0.0530 either 0.0623 xray 0.1033 dysp 0.4350 bronc 0.4613 smoke 0.5003 dtype: float64
Smoking is a coin flip and bronchitis is close behind, but the diseases are rare. A visit to
Asia sits near 1 percent, and tuberculosis is rarer still. Hold onto that. A cause that fires
in one row out of a hundred leaves a faint trace in the data, and a faint trace is hard for
any test to separate from noise. The arrow asia -> tub connects the two rarest variables in
the table, so if discovery is going to miss an edge, that is the one to bet on.
Run discovery and grade it¶
The truth is a DAG built straight from the known edges. The estimate comes from one call to
discover with the PC algorithm, which is the constraint-based method the
discovery guide walks through. It returns a CPDAG:
arrows where the data decides direction, dashes where several directions explain the data
equally well.
from ergodic import dag, discover
truth = dag(list(ds.target["edges"]))
result = discover(ds.data, method="pc")
print(result.graph.kind.value)
sorted(e.glyph() for e in result.graph.edges)
cpdag
['bronc --- dysp', 'bronc --> smoke', 'either <-- lung', 'either <-- tub', 'lung --> smoke']
To grade it, compare skeletons. The skeleton is the graph with directions thrown away, just the set of pairs that share an edge. Scoring the skeleton asks the first question that matters: did discovery find the right connections at all? Direction is a separate question, and on a CPDAG some directions are genuinely undecidable, so skeleton recall is the honest headline.
def skeleton(graph):
return {frozenset((x, y)) for x in graph.nodes for y in graph.neighbors(x)}
truth_skeleton = skeleton(truth)
found_skeleton = skeleton(result.graph)
recall = len(truth_skeleton & found_skeleton) / len(truth_skeleton)
precision = len(truth_skeleton & found_skeleton) / max(1, len(found_skeleton))
print(f"recall {recall:.3f} (found edges that are real)")
print(f"precision {precision:.3f} (real edges among those found)")
recall 0.625 (found edges that are real) precision 1.000 (real edges among those found)
Recall is 0.625 and precision is 1.0. Read those two numbers together. Precision 1.0 says every edge PC drew is a true edge: it invented nothing, added no phantom link between variables that only correlate. That is the trait you want from a discovery method on data you cannot check, because a false edge is a confounder you will adjust for that was never there. Recall 0.625 says it found five of the eight real edges and left three on the table.
Which three? Name them, because the misses are not random.
missed = truth_skeleton - found_skeleton
for pair in sorted(tuple(sorted(p)) for p in missed):
print(" - ".join(pair))
asia - tub dysp - either either - xray
asia - tub, either - xray, and either - dysp. Two of the three touch either, and the
third is the rare-cause edge we flagged from the marginals. Both misses have an honest reason.
The asia -> tub edge links the two rarest variables in the table. Tuberculosis appears in
under 1 percent of patients and a visit to Asia in about 1 percent, so the slice where Asia
moves the tuberculosis rate is a handful of rows in three thousand. The signal is
real but tiny, and PC's independence test cannot yet tell it from sampling noise. This is not
a flaw in the algorithm. It is the data being honest that it has not seen enough.
The either misses come from how that node is built. In the network either is a logical or:
it is "yes" exactly when tuberculosis or lung cancer is present, with no randomness of its own.
A deterministic node is a strange object for a conditional-independence test. Once you know
tuberculosis and lung cancer, you know either perfectly, which scrambles the independence
pattern its children depend on, and the edges from either to the x-ray and to dyspnoea get
lost. Determinism breaks an assumption these tests quietly rely on, so the honest result is a
gap, not a confident wrong arrow.
Here is the recovered graph. The five edges PC kept, drawn over the clinic.
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch
POS = {
"asia": (0.0, 2.0),
"tub": (0.0, 1.0),
"smoke": (3.0, 2.0),
"lung": (2.4, 1.0),
"bronc": (3.6, 1.0),
"either": (1.2, 1.0),
"xray": (0.6, 0.0),
"dysp": (2.4, 0.0),
}
def draw_graph(graph, pos, title, missed_pairs=()):
fig, ax = plt.subplots(figsize=(7, 4.2))
for a, b in missed_pairs:
ax.plot(
[pos[a][0], pos[b][0]],
[pos[a][1], pos[b][1]],
color="#d1495b",
ls=":",
lw=1.3,
zorder=1,
)
for edge in graph.edges:
u, mid, v = edge.glyph().split()
if mid == "<--":
u, v, mid = v, u, "-->"
style = "-" if mid == "---" else "-|>"
ax.add_patch(
FancyArrowPatch(
pos[u],
pos[v],
arrowstyle=style,
mutation_scale=16,
shrinkA=30,
shrinkB=30,
lw=1.6,
color="#33415c",
zorder=2,
)
)
for name, (x, y) in pos.items():
ax.text(
x,
y,
name,
ha="center",
va="center",
fontsize=9,
bbox={"boxstyle": "round,pad=0.4", "fc": "#eef1f6", "ec": "#33415c"},
zorder=3,
)
ax.set_xlim(-0.6, 4.2)
ax.set_ylim(-0.6, 2.6)
ax.set_title(title, fontsize=10)
ax.axis("off")
fig.tight_layout()
draw_graph(
result.graph,
POS,
"Asia at n=3000: five solid edges, three missed (dotted red)",
missed_pairs=[tuple(sorted(p)) for p in missed],
)
The funnel that survives is the spine of the clinic. Smoking, lung cancer, bronchitis, and the
two either parents all sit where they should. The dotted red lines are the three edges
discovery could not yet justify, and every one of them runs into a rare variable or the
deterministic node. The graph is partial, and it is partial in exactly the places the data is
thin.
More data sharpens the picture¶
The rare-cause miss was a sample-size problem, so give the algorithm more samples. The loader takes a row count, and nothing else about the world changes.
ds_big = load_asia(n=8000)
result_big = discover(ds_big.data, method="pc")
found_big = skeleton(result_big.graph)
recall_big = len(truth_skeleton & found_big) / len(truth_skeleton)
precision_big = len(truth_skeleton & found_big) / max(1, len(found_big))
print(f"n=8000 recall {recall_big:.3f} precision {precision_big:.3f}")
missed_big = truth_skeleton - found_big
print("still missed:", sorted(" - ".join(sorted(p)) for p in missed_big))
n=8000 recall 0.750 precision 1.000 still missed: ['dysp - either', 'either - xray']
Recall climbs from 0.625 to 0.75, and precision stays pinned at 1.0. The edge that came back
is asia -> tub: with 8000 rows instead of 3000, the faint trace of a rare cause finally
clears the test's bar, and PC draws it. No knob was turned and no hint was added. The same
method, fed more of the same data, recovered an edge it had honestly declined before.
What is still missing are the two either edges, and they stay missing because their cause is
structural, not statistical. More rows do not undo determinism. That distinction is the useful
lesson here. A miss from thin data heals with sample size; a miss from a broken assumption does
not, and it tells you something real about the world rather than the dataset.
draw_graph(result_big.graph, POS, "Asia at n=8000: asia -> tub returns")
A clean network for contrast¶
Asia recovers partly because Asia is hard: rare diseases and a logical gate. A simpler network shows what discovery does when the world cooperates. The burglar-alarm network (Pearl, 1988) has five variables and four edges. A burglary or an earthquake can trip an alarm, and the alarm is what makes a neighbor call.
from ergodic.datasets import load_earthquake
eq = load_earthquake()
print("nodes:", eq.target["nodes"])
print("true edges:", eq.target["edges"])
res = discover(eq.data, method="pc")
eq_truth = skeleton(dag(list(eq.target["edges"])))
eq_found = skeleton(res.graph)
print(f"recall {len(eq_truth & eq_found) / len(eq_truth):.3f}")
sorted(e.glyph() for e in res.graph.edges)
nodes: ('Burglary', 'Earthquake', 'Alarm', 'JohnCalls', 'MaryCalls')
true edges: ('Burglary -> Alarm', 'Earthquake -> Alarm', 'Alarm -> JohnCalls', 'Alarm -> MaryCalls')
recall 1.000
['Alarm --> JohnCalls', 'Alarm --> MaryCalls', 'Alarm <-- Burglary', 'Alarm <-- Earthquake']
Full skeleton, recall 1.0, and the arrows are oriented too. The interesting one is the alarm. Both a burglary and an earthquake point into it, and that meeting of two arrowheads at one node is a collider, the single structure observational data can orient on its own. Burglaries and earthquakes are independent until the alarm goes off; learning the alarm fired makes them dependent, because one cause explains away the other. That signature is what fixes the arrows, and you can read it straight off the recovered graph.
from ergodic import Mark
at_alarm_b = res.graph.mark_at("Burglary", "Alarm", "Alarm")
at_alarm_e = res.graph.mark_at("Earthquake", "Alarm", "Alarm")
print("Burglary -> Alarm, mark at Alarm:", at_alarm_b)
print("Earthquake -> Alarm, mark at Alarm:", at_alarm_e)
print("both arrowheads into Alarm:", at_alarm_b is Mark.ARROW and at_alarm_e is Mark.ARROW)
Burglary -> Alarm, mark at Alarm: Mark.ARROW Earthquake -> Alarm, mark at Alarm: Mark.ARROW both arrowheads into Alarm: True
Both endpoints at the alarm are arrowheads, so PC recovered the v-structure exactly. The contrast with Asia is the message. A small clean network with strong signals recovers in full, while a harder one recovers in part and is plain about which part. Same algorithm, two honesties, and neither of them pretends to know more than the data supports.
What this means for real data¶
On a real table you do not get an answer key, so you will not get a recall number. What you keep is the behavior these known networks just demonstrated. PC adds no edges it cannot support, it recovers the structure where the signal is strong, and it leaves gaps where a cause is rare or an assumption is strained. A partial graph that is honest about its gaps beats a complete graph that guessed.
A discovered graph is not the end of the line. It feeds straight into the rest of the library:
identify_effect reads a CPDAG and estimate_effect runs off it, so the arrows you just
learned become the effects you wanted, with the bookkeeping shown in the
discovery guide. For the next step toward real data,
the protein-signaling network of Sachs et al. (2005) is on the roadmap: a biology benchmark
where the true graph is known from intervention experiments, so discovery can be graded
against bench science instead of a simulator.