{ "cells": [ { "cell_type": "markdown", "id": "markov_check-0", "metadata": {}, "source": [ "# Markov check: do the data contradict a graph?\n", "\n", "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/Abel-ai-lab/andrey/blob/main/examples/markov_check.ipynb)\n", "\n", "A graph implies conditional independences: each variable is independent of every earlier variable\n", "in a causal order that is not its parent, given its parents. `andrey.markov_check` tests each one\n", "with a conditional-independence test and asks whether the p-values look uniform, as they do when\n", "every independence holds. A small uniformity p-value rejects the graph. Here the check runs on data\n", "simulated on the ALARM network, whose true graph we know, to show what it catches, what it misses,\n", "and how it can choose among graphs." ] }, { "cell_type": "markdown", "id": "markov_check-1", "metadata": {}, "source": [ "## Start\n", "\n", "Choose a CPU runtime and run this cell first. In Colab it installs `andrey-core==0.13.0` from PyPI;\n", "local kernels use the installed environment." ] }, { "cell_type": "code", "execution_count": 1, "id": "markov_check-2", "metadata": { "tags": [ "colab-setup" ] }, "outputs": [], "source": [ "import os\n", "import sys\n", "\n", "IN_COLAB = \"google.colab\" in sys.modules\n", "\n", "if IN_COLAB:\n", " os.environ[\"ANDREY_DEVICE\"] = \"cpu\"\n", " !pip install --quiet \"andrey-core[viz]==0.13.0\" pandas" ] }, { "cell_type": "code", "execution_count": 2, "id": "markov_check-3", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "\n", "import andrey\n", "import andrey.viz\n", "from andrey import GraphStructure\n", "from andrey.data import load_dataset\n", "from andrey.metrics import shd, to_cpdag" ] }, { "cell_type": "markdown", "id": "markov_check-4", "metadata": {}, "source": [ "ALARM is a published network of 37 variables and 46 arcs. `load_dataset` simulates rows on it from\n", "a linear-Gaussian model. A graph fitted to some rows can pass a check on those rows because the fit\n", "made its independences true there, so the check needs rows the graph was not learned from. We learn\n", "on the first 1,000 rows and check on the other 1,000." ] }, { "cell_type": "code", "execution_count": 3, "id": "markov_check-5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "37 variables, 46 true edges\n" ] } ], "source": [ "alarm = load_dataset(\"alarm\", n=2000, seed=1)\n", "data, truth = alarm.data, alarm.graph\n", "learn, held = data.iloc[:1000], data.iloc[1000:]\n", "true_cpdag = to_cpdag(truth)\n", "print(data.shape[1], \"variables,\", len(truth.oriented_edges()), \"true edges\")" ] }, { "cell_type": "markdown", "id": "markov_check-6", "metadata": {}, "source": [ "## Check a graph\n", "\n", "The true graph is the control: its independences hold, so the check should not reject it. The\n", "result prints the verdict, the number of independences tested, and the uniformity statistic." ] }, { "cell_type": "code", "execution_count": 4, "id": "markov_check-7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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.\n", " 620 tests | Anderson-Darling A^2 = 0.592 | each test on 500 rows | indep_test fisherz\n" ] } ], "source": [ "check = andrey.markov_check(truth, held)\n", "print(check)" ] }, { "cell_type": "markdown", "id": "markov_check-8", "metadata": {}, "source": [ "## A damaged graph is rejected\n", "\n", "Remove the edge `CO -> BP` from the true graph. The graph now claims that `BP` is independent of\n", "the causes of `CO`, such as `LVFAILURE`, given `TPR`, which the data contradict. `table()` lists\n", "every tested independence with its p-value, so the smallest ones show where the graph is wrong:" ] }, { "cell_type": "code", "execution_count": 5, "id": "markov_check-9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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.\n", " 621 tests | Anderson-Darling A^2 = 48.4 | each test on 500 rows | indep_test fisherz\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
xygivenp_value
185BPINSUFFANESTH[TPR]0.0
188BPLVFAILURE[TPR]0.0
197BPSTROKEVOLUME[TPR]0.0
\n", "
" ], "text/plain": [ " x y given p_value\n", "185 BP INSUFFANESTH [TPR] 0.0\n", "188 BP LVFAILURE [TPR] 0.0\n", "197 BP STROKEVOLUME [TPR] 0.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def with_pair(graph, source, target, edge):\n", " \"\"\"Return a copy of the DAG with `source -> target` added, or removed when not `edge`.\"\"\"\n", " marks = graph.to_numpy().copy()\n", " i, j = graph.labels.index(source), graph.labels.index(target)\n", " marks[i, j], marks[j, i] = (1, 2) if edge else (0, 0) # tail at source, arrowhead at target\n", " return GraphStructure.from_numpy(marks, kind=\"dag\", labels=graph.labels)\n", "\n", "\n", "damaged = with_pair(truth, \"CO\", \"BP\", edge=False)\n", "check = andrey.markov_check(damaged, held)\n", "print(check)\n", "pd.DataFrame(check.table()).nsmallest(3, \"p_value\")" ] }, { "cell_type": "markdown", "id": "markov_check-10", "metadata": {}, "source": [ "## What a pass does not show\n", "\n", "The check rejects a missing edge only when the data can see the dependence. `PULMEMBOLUS -> PAP`\n", "is weak enough here that a graph without it still passes. A graph with an extra edge passes too,\n", "because an extra edge removes an independence from the graph's claims and leaves the rest true.\n", "So a pass says the data do not contradict the graph, not that its edges are right:" ] }, { "cell_type": "code", "execution_count": 6, "id": "markov_check-11", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
graphedgestestsuniformity pverdict
0without PULMEMBOLUS -> PAP456210.216not rejected
1plus HYPOVOLEMIA -> CVP476190.903not rejected
\n", "
" ], "text/plain": [ " graph edges tests uniformity p verdict\n", "0 without PULMEMBOLUS -> PAP 45 621 0.216 not rejected\n", "1 plus HYPOVOLEMIA -> CVP 47 619 0.903 not rejected" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "variants = {\n", " \"without PULMEMBOLUS -> PAP\": with_pair(truth, \"PULMEMBOLUS\", \"PAP\", edge=False),\n", " \"plus HYPOVOLEMIA -> CVP\": with_pair(truth, \"HYPOVOLEMIA\", \"CVP\", edge=True),\n", "}\n", "rows = []\n", "for name, graph in variants.items():\n", " check = andrey.markov_check(graph, held)\n", " rows.append(\n", " {\n", " \"graph\": name,\n", " \"edges\": len(graph.oriented_edges()),\n", " \"tests\": check.n_tests,\n", " \"uniformity p\": round(check.uniformity_p_value, 3),\n", " \"verdict\": \"not rejected\" if check.passed else \"rejected\",\n", " }\n", " )\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "markov_check-12", "metadata": {}, "source": [ "## Choosing a graph: the sparsest one that passes\n", "\n", "When the truth is unknown, Ramsey, Andrews, and Spirtes (2024) use the check to choose among the\n", "graphs a method returns over a range of one setting: drop the graphs the data reject and keep the\n", "sparsest of the rest, the simplest explanation the data allow. GES's `penalty_discount` multiplies\n", "the BIC penalty, and a larger value gives fewer edges. Fit GES on the learn rows over five values\n", "and check each fit on the held rows:" ] }, { "cell_type": "code", "execution_count": 7, "id": "markov_check-13", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
penalty_discountedgeserrorsuniformity pverdict
00.559180.267not rejected
11.05290.675not rejected
22.05180.149not rejected
34.050150.247not rejected
48.050150.247not rejected
\n", "
" ], "text/plain": [ " penalty_discount edges errors uniformity p verdict\n", "0 0.5 59 18 0.267 not rejected\n", "1 1.0 52 9 0.675 not rejected\n", "2 2.0 51 8 0.149 not rejected\n", "3 4.0 50 15 0.247 not rejected\n", "4 8.0 50 15 0.247 not rejected" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows, fits = [], {}\n", "for penalty in (0.5, 1, 2, 4, 8):\n", " fit = andrey.ges(learn, penalty_discount=penalty)\n", " check = andrey.markov_check(fit, held)\n", " fits[penalty] = fit\n", " rows.append(\n", " {\n", " \"penalty_discount\": penalty,\n", " \"edges\": len(fit.structure.oriented_edges()),\n", " \"errors\": shd(fit.structure, true_cpdag),\n", " \"uniformity p\": round(check.uniformity_p_value, 3),\n", " \"verdict\": \"not rejected\" if check.passed else \"rejected\",\n", " }\n", " )\n", "sweep = pd.DataFrame(rows)\n", "sweep" ] }, { "cell_type": "markdown", "id": "markov_check-14", "metadata": {}, "source": [ "The data reject none of the five graphs, so the rule picks the sparsest, 50 edges at\n", "`penalty_discount` 4 (8 returns the same graph). It is not the most accurate: its 15 errors against\n", "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.\n", "The check rules out graphs the data contradict, and among the rest fewer edges is a preference,\n", "not a measure of accuracy." ] }, { "cell_type": "code", "execution_count": 8, "id": "markov_check-15", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sparsest graph that passes: penalty_discount 4.0\n" ] }, { "data": { "text/html": [ "PAPHISTORYCOTPRPULMEMBOLUSBPHRBPHREKGEXPCO2VENTMACHLVEDVOLUMEANAPHYLAXISLVFAILUREPVSATMINVOLINSUFFANESTHVENTALVMINVOLSETVENTLUNGFIO2HRVENTTUBEDISCONNECTERRCAUTERKINKEDTUBESTROKEVOLUMEHYPOVOLEMIAERRLOWOUTPUTPCWPCATECHOLSAO2INTUBATIONSHUNTCVPHRSATARTCO2PRESS" ], "text/plain": [ "'PAPHISTORYCOTPRPULMEMBOLUSBPHRBPHREKGEXPCO2VENTMACHLVEDVOLUMEANAPHYLAXISLVFAILUREPVSATMINVOLINSUFFANESTHVENTALVMINVOLSETVENTLUNGFIO2HRVENTTUBEDISCONNECTERRCAUTERKINKEDTUBESTROKEVOLUMEHYPOVOLEMIAERRLOWOUTPUTPCWPCATECHOLSAO2INTUBATIONSHUNTCVPHRSATARTCO2PRESS'" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "passing = sweep[sweep.verdict == \"not rejected\"]\n", "sparsest = passing.loc[passing.edges.idxmin(), \"penalty_discount\"]\n", "print(\"sparsest graph that passes: penalty_discount\", sparsest)\n", "andrey.viz.draw(fits[sparsest].structure)" ] }, { "cell_type": "markdown", "id": "markov_check-16", "metadata": {}, "source": [ "## PC: no graph passes\n", "\n", "The same sweep over PC's `alpha`, the level of its independence tests. PC returns a CPDAG, which\n", "the check reads through a DAG of its class; a graph with no consistent DAG cannot be checked, and\n", "`markov_check` raises `ValueError`:" ] }, { "cell_type": "code", "execution_count": 9, "id": "markov_check-17", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
alphaedgeserrorsuniformity pverdict
00.00135250.0rejected
10.01036240.0rejected
20.05042240.0rejected
30.10043250.0rejected
40.3005231NaNno DAG to check
\n", "
" ], "text/plain": [ " alpha edges errors uniformity p verdict\n", "0 0.001 35 25 0.0 rejected\n", "1 0.010 36 24 0.0 rejected\n", "2 0.050 42 24 0.0 rejected\n", "3 0.100 43 25 0.0 rejected\n", "4 0.300 52 31 NaN no DAG to check" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows = []\n", "for alpha in (0.001, 0.01, 0.05, 0.1, 0.3):\n", " fit = andrey.pc(learn, alpha=alpha)\n", " row = {\n", " \"alpha\": alpha,\n", " \"edges\": len(fit.structure.oriented_edges()),\n", " \"errors\": shd(fit.structure, true_cpdag),\n", " }\n", " try:\n", " check = andrey.markov_check(fit, held)\n", " row[\"uniformity p\"] = round(check.uniformity_p_value, 3)\n", " row[\"verdict\"] = \"not rejected\" if check.passed else \"rejected\"\n", " except ValueError:\n", " row[\"uniformity p\"], row[\"verdict\"] = None, \"no DAG to check\"\n", " rows.append(row)\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "markov_check-18", "metadata": {}, "source": [ "At every `alpha` the data reject the graph PC returns, or the graph has no DAG to check, so there\n", "is nothing to choose from. This is a finding about PC on these data and this sample size: PC's\n", "graphs differ from the true CPDAG by many more errors than GES's, and the check agrees that they do\n", "not fit the data. It does not say which edges to change; the smallest p-values in `table()` show\n", "where to look. `andrey markov-check` runs the same check from the shell." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.4" } }, "nbformat": 4, "nbformat_minor": 5 }