NotebooksDiscovery on a known network: grading the graph

Static notebook preview · outputs are cleared. Download the notebook to run the cells.

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.

python

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

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?

python

df = ds.data.to_pandas()
rates = (df == "yes").mean().sort_values()
print(rates.round(4))

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.

python

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)

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.

python

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 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.

python

missed = truth_skeleton - found_skeleton
for pair in sorted(tuple(sorted(p)) for p in missed):
    print(" - ".join(pair))

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.

python

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.

python

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))

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.

python

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.

python

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)

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.

python

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)

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.

Carrying this back to a real table#

On real data you never get an answer key, so you never get a recall number either. What carries over is the behavior these known networks just showed you. PC adds no edges it can't 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's 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.

Suggest an edit