Examples

Markov check: do the data contradict a graph?#

Open in Colab Download this notebook

A graph implies conditional independences: each variable is independent of every earlier variable in a causal order that is not its parent, given its parents. andrey.markov_check tests each one with a conditional-independence test and asks whether the p-values look uniform, as they do when every independence holds. A small uniformity p-value rejects the graph. Here the check runs on data simulated on the ALARM network, whose true graph we know, to show what it catches, what it misses, and how it can choose among graphs.

Start#

Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.13.0 from PyPI; local kernels use the installed environment.

import os
import sys

IN_COLAB = "google.colab" in sys.modules

if IN_COLAB:
    os.environ["ANDREY_DEVICE"] = "cpu"
    !pip install --quiet "andrey-core[viz]==0.13.0" pandas

import pandas as pd

import andrey
import andrey.viz
from andrey import GraphStructure
from andrey.data import load_dataset
from andrey.metrics import shd, to_cpdag

ALARM is a published network of 37 variables and 46 arcs. load_dataset simulates rows on it from a linear-Gaussian model. A graph fitted to some rows can pass a check on those rows because the fit made its independences true there, so the check needs rows the graph was not learned from. We learn on the first 1,000 rows and check on the other 1,000.

alarm = load_dataset("alarm", n=2000, seed=1)
data, truth = alarm.data, alarm.graph
learn, held = data.iloc[:1000], data.iloc[1000:]
true_cpdag = to_cpdag(truth)
print(data.shape[1], "variables,", len(truth.oriented_edges()), "true edges")

37 variables, 46 true edges

Check a graph#

The true graph is the control: its independences hold, so the check should not reject it. The result prints the verdict, the number of independences tested, and the uniformity statistic.

check = andrey.markov_check(truth, held)
print(check)

Not rejected: the p-values of the 620 independences the graph implies are consistent with uniform (Anderson-Darling p = 0.656 > 0.05); 5.6% have p < 0.05, about 5% expected. A pass does not show the edges are right: a graph with extra edges also passes.
  620 tests | Anderson-Darling A^2 = 0.592 | each test on 500 rows | indep_test fisherz

A damaged graph is rejected#

Remove the edge CO -> BP from the true graph. The graph now claims that BP is independent of the causes of CO, such as LVFAILURE, given TPR, which the data contradict. table() lists every tested independence with its p-value, so the smallest ones show where the graph is wrong:

def with_pair(graph, source, target, edge):
    """Return a copy of the DAG with `source -> target` added, or removed when not `edge`."""
    marks = graph.to_numpy().copy()
    i, j = graph.labels.index(source), graph.labels.index(target)
    marks[i, j], marks[j, i] = (1, 2) if edge else (0, 0)  # tail at source, arrowhead at target
    return GraphStructure.from_numpy(marks, kind="dag", labels=graph.labels)


damaged = with_pair(truth, "CO", "BP", edge=False)
check = andrey.markov_check(damaged, held)
print(check)
pd.DataFrame(check.table()).nsmallest(3, "p_value")

Rejected: the p-values of the 621 independences the graph implies are not uniform (Anderson-Darling p = 9.66e-07 <= 0.05); 8.5% have p < 0.05, about 5% expected.
  621 tests | Anderson-Darling A^2 = 48.4 | each test on 500 rows | indep_test fisherz
x y given p_value
185 BP INSUFFANESTH [TPR] 0.0
188 BP LVFAILURE [TPR] 0.0
197 BP STROKEVOLUME [TPR] 0.0

What a pass does not show#

The check rejects a missing edge only when the data can see the dependence. PULMEMBOLUS -> PAP is weak enough here that a graph without it still passes. A graph with an extra edge passes too, because an extra edge removes an independence from the graph’s claims and leaves the rest true. So a pass says the data do not contradict the graph, not that its edges are right:

variants = {
    "without PULMEMBOLUS -> PAP": with_pair(truth, "PULMEMBOLUS", "PAP", edge=False),
    "plus HYPOVOLEMIA -> CVP": with_pair(truth, "HYPOVOLEMIA", "CVP", edge=True),
}
rows = []
for name, graph in variants.items():
    check = andrey.markov_check(graph, held)
    rows.append(
        {
            "graph": name,
            "edges": len(graph.oriented_edges()),
            "tests": check.n_tests,
            "uniformity p": round(check.uniformity_p_value, 3),
            "verdict": "not rejected" if check.passed else "rejected",
        }
    )
pd.DataFrame(rows)
graph edges tests uniformity p verdict
0 without PULMEMBOLUS -> PAP 45 621 0.216 not rejected
1 plus HYPOVOLEMIA -> CVP 47 619 0.903 not rejected

Choosing a graph: the sparsest one that passes#

When the truth is unknown, Ramsey, Andrews, and Spirtes (2024) use the check to choose among the graphs a method returns over a range of one setting: drop the graphs the data reject and keep the sparsest of the rest, the simplest explanation the data allow. GES’s penalty_discount multiplies the BIC penalty, and a larger value gives fewer edges. Fit GES on the learn rows over five values and check each fit on the held rows:

rows, fits = [], {}
for penalty in (0.5, 1, 2, 4, 8):
    fit = andrey.ges(learn, penalty_discount=penalty)
    check = andrey.markov_check(fit, held)
    fits[penalty] = fit
    rows.append(
        {
            "penalty_discount": penalty,
            "edges": len(fit.structure.oriented_edges()),
            "errors": shd(fit.structure, true_cpdag),
            "uniformity p": round(check.uniformity_p_value, 3),
            "verdict": "not rejected" if check.passed else "rejected",
        }
    )
sweep = pd.DataFrame(rows)
sweep
penalty_discount edges errors uniformity p verdict
0 0.5 59 18 0.267 not rejected
1 1.0 52 9 0.675 not rejected
2 2.0 51 8 0.149 not rejected
3 4.0 50 15 0.247 not rejected
4 8.0 50 15 0.247 not rejected

The data reject none of the five graphs, so the rule picks the sparsest, 50 edges at penalty_discount 4 (8 returns the same graph). It is not the most accurate: its 15 errors against the true CPDAG are more than the 8 of the 51-edge graph at 2 and the 9 of the 52-edge graph at 1. The check rules out graphs the data contradict, and among the rest fewer edges is a preference, not a measure of accuracy.

passing = sweep[sweep.verdict == "not rejected"]
sparsest = passing.loc[passing.edges.idxmin(), "penalty_discount"]
print("sparsest graph that passes: penalty_discount", sparsest)
andrey.viz.draw(fits[sparsest].structure)

sparsest graph that passes: penalty_discount 4.0
PAPHISTORYCOTPRPULMEMBOLUSBPHRBPHREKGEXPCO2VENTMACHLVEDVOLUMEANAPHYLAXISLVFAILUREPVSATMINVOLINSUFFANESTHVENTALVMINVOLSETVENTLUNGFIO2HRVENTTUBEDISCONNECTERRCAUTERKINKEDTUBESTROKEVOLUMEHYPOVOLEMIAERRLOWOUTPUTPCWPCATECHOLSAO2INTUBATIONSHUNTCVPHRSATARTCO2PRESS

PC: no graph passes#

The same sweep over PC’s alpha, the level of its independence tests. PC returns a CPDAG, which the check reads through a DAG of its class; a graph with no consistent DAG cannot be checked, and markov_check raises ValueError:

rows = []
for alpha in (0.001, 0.01, 0.05, 0.1, 0.3):
    fit = andrey.pc(learn, alpha=alpha)
    row = {
        "alpha": alpha,
        "edges": len(fit.structure.oriented_edges()),
        "errors": shd(fit.structure, true_cpdag),
    }
    try:
        check = andrey.markov_check(fit, held)
        row["uniformity p"] = round(check.uniformity_p_value, 3)
        row["verdict"] = "not rejected" if check.passed else "rejected"
    except ValueError:
        row["uniformity p"], row["verdict"] = None, "no DAG to check"
    rows.append(row)
pd.DataFrame(rows)
alpha edges errors uniformity p verdict
0 0.001 35 25 0.0 rejected
1 0.010 36 24 0.0 rejected
2 0.050 42 24 0.0 rejected
3 0.100 43 25 0.0 rejected
4 0.300 52 31 NaN no DAG to check

At every alpha the data reject the graph PC returns, or the graph has no DAG to check, so there is nothing to choose from. This is a finding about PC on these data and this sample size: PC’s graphs differ from the true CPDAG by many more errors than GES’s, and the check agrees that they do not fit the data. It does not say which edges to change; the smallest p-values in table() show where to look. andrey markov-check runs the same check from the shell.