{ "cells": [ { "cell_type": "markdown", "id": "bootstrap-0", "metadata": {}, "source": [ "# Bootstrap: how stable is each edge?\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/bootstrap.ipynb)\n", "\n", "A method fitted to one sample returns one graph, and another sample from the same source would\n", "give a somewhat different one. `andrey.bootstrap` refits the method on many resamples of the rows\n", "and counts how often each edge appears. Here PC runs on data simulated on the ALARM network, whose\n", "true graph we know, to show how to read those counts, what the averaged graph is, and which edges\n", "to trust." ] }, { "cell_type": "markdown", "id": "bootstrap-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": "bootstrap-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": "bootstrap-3", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import andrey\n", "import andrey.viz\n", "from andrey.data import load_dataset\n", "from andrey.metrics import shd, to_cpdag" ] }, { "cell_type": "markdown", "id": "bootstrap-4", "metadata": {}, "source": [ "ALARM is a published network of 37 variables and 46 arcs. `load_dataset` simulates 1,000 rows\n", "on it from a linear-Gaussian model and lists the columns in an order drawn from the seed. From\n", "these data PC can recover at best the CPDAG of the true graph, so that is what we score against." ] }, { "cell_type": "code", "execution_count": 3, "id": "bootstrap-5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 rows, 37 variables, 46 true edges\n" ] } ], "source": [ "alarm = load_dataset(\"alarm\", n=1000, seed=1)\n", "data, truth = alarm.data, alarm.graph\n", "true_cpdag = to_cpdag(truth)\n", "upper = np.triu_indices(truth.n_nodes, 1) # each pair of variables once\n", "true_pair = (truth.to_numpy() != 0)[upper]\n", "print(data.shape[0], \"rows,\", data.shape[1], \"variables,\", int(true_pair.sum()), \"true edges\")" ] }, { "cell_type": "markdown", "id": "bootstrap-6", "metadata": {}, "source": [ "## Edge frequencies rank the edges\n", "\n", "`bootstrap` fits PC on 50 resamples, each half of the rows drawn without replacement, the default.\n", "`frequency(\"absent\")` holds the share of the fits in which each pair has no edge, so one minus it\n", "is the share in which the pair is adjacent. Here are the pairs grouped by that share:" ] }, { "cell_type": "code", "execution_count": 4, "id": "bootstrap-7", "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", "
adjacent inpairstrue pairs
0more than 10%, up to 30%62
1more than 30%, up to 50%43
2more than 50%, up to 70%00
3more than 70%, up to 90%22
4more than 90%, up to 100%3231
\n", "
" ], "text/plain": [ " adjacent in pairs true pairs\n", "0 more than 10%, up to 30% 6 2\n", "1 more than 30%, up to 50% 4 3\n", "2 more than 50%, up to 70% 0 0\n", "3 more than 70%, up to 90% 2 2\n", "4 more than 90%, up to 100% 32 31" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "out = andrey.bootstrap(\"pc\", data, n_resamples=50, seed=1)\n", "adjacent = (1 - out.frequency(\"absent\"))[upper]\n", "\n", "bins = [0.1, 0.3, 0.5, 0.7, 0.9, 1.0]\n", "rows = []\n", "for low, high in zip(bins[:-1], bins[1:]):\n", " inside = (adjacent > low) & (adjacent <= high)\n", " rows.append(\n", " {\n", " \"adjacent in\": f\"more than {low:.0%}, up to {high:.0%}\",\n", " \"pairs\": int(inside.sum()),\n", " \"true pairs\": int((inside & true_pair).sum()),\n", " }\n", " )\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "bootstrap-8", "metadata": {}, "source": [ "The pairs adjacent in more than 90% of the fits are nearly all true edges. The shares in between\n", "are not probabilities: most pairs adjacent in 30% to 50% of the fits are true edges here, not 30%\n", "to 50% of them. A frequency says how stable an edge is when the rows change, which ranks the\n", "edges; how often the stable ones are right depends on the method, the data, and the sample\n", "size." ] }, { "cell_type": "markdown", "id": "bootstrap-9", "metadata": {}, "source": [ "## The averaged graph\n", "\n", "`out.structure` is the averaged graph. Each pair takes its most frequent state in the fits, with\n", "no edge counted as one of the states. PC returns CPDAGs, and a vote of CPDAGs is often not a\n", "CPDAG, so the vote is then replaced by the CPDAG of a DAG with the same adjacencies, v-structures,\n", "and directed edges. The result is a valid CPDAG, here with fewer errors than PC on all the rows:" ] }, { "cell_type": "code", "execution_count": 5, "id": "bootstrap-10", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "single fit: 24 errors\n", "averaged graph: 21 errors\n", "kind: cpdag | built from: vote\n", "the CPDAG of one of its DAGs: True\n" ] } ], "source": [ "single = andrey.pc(data)\n", "print(\"single fit: \", shd(single.structure, true_cpdag), \"errors\")\n", "print(\"averaged graph:\", shd(out.structure, true_cpdag), \"errors\")\n", "print(\"kind:\", out.structure.kind, \"| built from:\", out.metadata[\"averaged_graph\"])\n", "dag = out.structure.to_dag(tie_break=\"label\")\n", "print(\n", " \"the CPDAG of one of its DAGs:\",\n", " np.array_equal(to_cpdag(dag).to_numpy(), out.structure.to_numpy()),\n", ")" ] }, { "cell_type": "code", "execution_count": 6, "id": "bootstrap-11", "metadata": {}, "outputs": [ { "data": { "text/html": [ "PAPHISTORYCOTPRPULMEMBOLUSBPHRBPHREKGEXPCO2VENTMACHLVEDVOLUMEANAPHYLAXISLVFAILUREPVSATMINVOLINSUFFANESTHVENTALVMINVOLSETVENTLUNGFIO2HRVENTTUBEDISCONNECTERRCAUTERKINKEDTUBESTROKEVOLUMEHYPOVOLEMIAERRLOWOUTPUTPCWPCATECHOLSAO2INTUBATIONSHUNTCVPHRSATARTCO2PRESS" ], "text/plain": [ "'PAPHISTORYCOTPRPULMEMBOLUSBPHRBPHREKGEXPCO2VENTMACHLVEDVOLUMEANAPHYLAXISLVFAILUREPVSATMINVOLINSUFFANESTHVENTALVMINVOLSETVENTLUNGFIO2HRVENTTUBEDISCONNECTERRCAUTERKINKEDTUBESTROKEVOLUMEHYPOVOLEMIAERRLOWOUTPUTPCWPCATECHOLSAO2INTUBATIONSHUNTCVPHRSATARTCO2PRESS'" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "andrey.viz.draw(out.structure)" ] }, { "cell_type": "markdown", "id": "bootstrap-12", "metadata": {}, "source": [ "## Subsampling or the bootstrap\n", "\n", "`resampling=\"bootstrap\"` draws as many rows as the data have, with replacement, instead of half of\n", "them without. A repeated row makes a test see more evidence than the resample holds, so weak and\n", "false edges survive in more of the fits. `sample_fraction` sets the share of the rows in a\n", "subsample: a larger share gives fits closer to the fit on all the rows, so the fits agree more and\n", "their frequencies say less about how much the result varies." ] }, { "cell_type": "code", "execution_count": 7, "id": "bootstrap-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", "
resamplingerrorsabove 90%: true of all50% to 90%: true of allaveraged graph
0subsample, 50% of the rows2131 of 322 of 2vote
1subsample, 80% of the rows1533 of 344 of 5vote
2bootstrap1932 of 334 of 5nearest_resample
\n", "
" ], "text/plain": [ " resampling errors above 90%: true of all \\\n", "0 subsample, 50% of the rows 21 31 of 32 \n", "1 subsample, 80% of the rows 15 33 of 34 \n", "2 bootstrap 19 32 of 33 \n", "\n", " 50% to 90%: true of all averaged graph \n", "0 2 of 2 vote \n", "1 4 of 5 vote \n", "2 4 of 5 nearest_resample " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "runs = {\n", " \"subsample, 50% of the rows\": {},\n", " \"subsample, 80% of the rows\": {\"sample_fraction\": 0.8},\n", " \"bootstrap\": {\"resampling\": \"bootstrap\"},\n", "}\n", "results, raised = {}, {}\n", "rows = []\n", "for name, settings in runs.items():\n", " with warnings.catch_warnings(record=True) as caught:\n", " warnings.simplefilter(\"always\")\n", " result = andrey.bootstrap(\"pc\", data, n_resamples=50, seed=1, **settings)\n", " results[name] = result\n", " raised[name] = [str(w.message) for w in caught if issubclass(w.category, andrey.AndreyWarning)]\n", " share = (1 - result.frequency(\"absent\"))[upper]\n", " above, middle = share > 0.9, (share > 0.5) & (share <= 0.9)\n", " rows.append(\n", " {\n", " \"resampling\": name,\n", " \"errors\": shd(result.structure, true_cpdag),\n", " \"above 90%: true of all\": f\"{(above & true_pair).sum()} of {above.sum()}\",\n", " \"50% to 90%: true of all\": f\"{(middle & true_pair).sum()} of {middle.sum()}\",\n", " \"averaged graph\": result.metadata[\"averaged_graph\"],\n", " }\n", " )\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "bootstrap-14", "metadata": {}, "source": [ "## When the vote has no valid graph\n", "\n", "With the bootstrap, the vote of these fits has no consistent extension: no DAG keeps all its\n", "adjacencies, v-structures, and directed edges, so no CPDAG holds the vote. `bootstrap` then\n", "returns the resample graph nearest the others among those a DAG extends, the fit whose state of\n", "each pair the most fits share. It is a graph PC itself returned, and a warning and\n", "`metadata[\"averaged_graph\"]` say so:" ] }, { "cell_type": "code", "execution_count": 8, "id": "bootstrap-15", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "the vote of the fits has no valid cpdag: PDAG admits no consistent DAG extension: no DAG has its skeleton, its v-structures, and its directed edges, because the undirected edges among nodes [2, 3, 4, 5, 11, 13, 15, 29, 30, 31, 32] cannot be oriented without adding a v-structure or closing a directed cycle. The averaged graph is resample 25, the valid fit nearest the others; strict=True raises instead\n", "built from: nearest_resample | kind: cpdag\n" ] } ], "source": [ "boot = results[\"bootstrap\"]\n", "print(\"\\n\".join(raised[\"bootstrap\"]))\n", "print(\"built from:\", boot.metadata[\"averaged_graph\"], \"| kind:\", boot.structure.kind)" ] }, { "cell_type": "markdown", "id": "bootstrap-16", "metadata": {}, "source": [ "`strict=True` raises instead, so you can decide from the frequencies, which keep the whole vote\n", "either way:" ] }, { "cell_type": "code", "execution_count": 9, "id": "bootstrap-17", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ValueError: the vote of the fits has no valid cpdag: PDAG admits no consistent DAG extension: no DAG has its skeleton, its v-structures, and its directed edges, because the undirected edges among nodes [2, 3, 4, 5, 11, 13, 15, 29, 30, 31, 32] cannot be oriented without adding a v-structure or closing a directed cycle. strict=False returns the fit nearest the others instead\n" ] } ], "source": [ "try:\n", " andrey.bootstrap(\"pc\", data, n_resamples=50, seed=1, resampling=\"bootstrap\", strict=True)\n", "except ValueError as error:\n", " print(\"ValueError:\", error)" ] }, { "cell_type": "markdown", "id": "cbd690bd", "metadata": {}, "source": [ "## Choosing the settings\n", "\n", "The defaults, 100 resamples of half of the rows each, suit most uses. Each setting trades one thing\n", "for another; here is how to choose it." ] }, { "cell_type": "markdown", "id": "7c7be4da", "metadata": {}, "source": [ "### How many resamples\n", "\n", "A frequency is a share of `n_resamples` fits, so it carries its own sampling noise: for a pair\n", "adjacent in a share `p` of all possible fits, about `sqrt(p (1 - p) / n_resamples)`, 0.10 at 25\n", "resamples and 0.05 at 100 for a pair near one half. Five seeds at each count show it on the pairs\n", "whose frequency lies between 20% and 80%:" ] }, { "cell_type": "code", "execution_count": 10, "id": "0ea86bc4", "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", "
resamplespairs between 20% and 80%spread over seeds (SD)expected SDpairs above 90%, per seed
02570.0900.090[32, 32, 32, 33, 32]
15070.0700.065[32, 32, 32, 33, 32]
210070.0450.046[32, 32, 32, 32, 32]
320070.0290.032[32, 32, 32, 32, 32]
\n", "
" ], "text/plain": [ " resamples pairs between 20% and 80% spread over seeds (SD) expected SD \\\n", "0 25 7 0.090 0.090 \n", "1 50 7 0.070 0.065 \n", "2 100 7 0.045 0.046 \n", "3 200 7 0.029 0.032 \n", "\n", " pairs above 90%, per seed \n", "0 [32, 32, 32, 33, 32] \n", "1 [32, 32, 32, 33, 32] \n", "2 [32, 32, 32, 32, 32] \n", "3 [32, 32, 32, 32, 32] " ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows = []\n", "for n_resamples in (25, 50, 100, 200):\n", " fits = [andrey.bootstrap(\"pc\", data, n_resamples=n_resamples, seed=s) for s in range(1, 6)]\n", " shares = np.array([(1 - fit.frequency(\"absent\"))[upper] for fit in fits])\n", " mean = shares.mean(axis=0)\n", " middle = (mean > 0.2) & (mean < 0.8)\n", " rows.append(\n", " {\n", " \"resamples\": n_resamples,\n", " \"pairs between 20% and 80%\": int(middle.sum()),\n", " \"spread over seeds (SD)\": shares[:, middle].std(axis=0, ddof=1).mean(),\n", " \"expected SD\": np.sqrt(mean[middle] * (1 - mean[middle]) / n_resamples).mean(),\n", " \"pairs above 90%, per seed\": [int((s > 0.9).sum()) for s in shares],\n", " }\n", " )\n", "pd.DataFrame(rows).round(3)" ] }, { "cell_type": "markdown", "id": "b6db3bb1", "metadata": {}, "source": [ "The spread follows the formula. The pairs above 90% barely move at any count, so 100 resamples\n", "rank the edges well. When a decision rests on one frequency near a threshold, such as whether a\n", "pair passes 0.5, use more: at 200 the spread is near 0.03. The time grows in proportion." ] }, { "cell_type": "markdown", "id": "a738a5e5", "metadata": {}, "source": [ "### How large a subsample\n", "\n", "`sample_fraction` sets how many rows each fit sees. A larger share gives fits closer to the fit on\n", "all the rows, so more pairs reach each level, true and false; but the resamples then share most of\n", "their rows, so the frequencies overstate how stable an edge is. Half of the rows keeps the fits\n", "apart and performs as well as the best correction of the bootstrap (Banerjee, Andrews, and\n", "Kummerfeld 2025). Here are the shares on ALARM at 200, 1,000, and 5,000 rows:" ] }, { "cell_type": "code", "execution_count": 11, "id": "e2de4d64", "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rowsshare of rowsabove 90%true above 90%above 50%true above 50%
020050%23232424
120080%25252928
2100050%32313433
3100080%34333937
4500050%35353838
5500080%38384139
\n", "
" ], "text/plain": [ " rows share of rows above 90% true above 90% above 50% true above 50%\n", "0 200 50% 23 23 24 24\n", "1 200 80% 25 25 29 28\n", "2 1000 50% 32 31 34 33\n", "3 1000 80% 34 33 39 37\n", "4 5000 50% 35 35 38 38\n", "5 5000 80% 38 38 41 39" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows = []\n", "for n_rows in (200, 1000, 5000):\n", " sample = load_dataset(\"alarm\", n=n_rows, seed=1).data\n", " for name, settings in {\"50%\": {}, \"80%\": {\"sample_fraction\": 0.8}}.items():\n", " result = andrey.bootstrap(\"pc\", sample, n_resamples=50, seed=1, **settings)\n", " adjacent = (1 - result.frequency(\"absent\"))[upper]\n", " rows.append(\n", " {\n", " \"rows\": n_rows,\n", " \"share of rows\": name,\n", " \"above 90%\": int((adjacent > 0.9).sum()),\n", " \"true above 90%\": int(((adjacent > 0.9) & true_pair).sum()),\n", " \"above 50%\": int((adjacent > 0.5).sum()),\n", " \"true above 50%\": int(((adjacent > 0.5) & true_pair).sum()),\n", " }\n", " )\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "c0651184", "metadata": {}, "source": [ "With few rows, a half-size fit loses the weak edges, so fewer pairs pass each level. Keep the\n", "default to rank edges by stability; a larger share finds more of the weak edges on small data, at\n", "the cost of frequencies that read higher. Whatever the share, compare frequencies only between runs\n", "with the same settings." ] }, { "cell_type": "markdown", "id": "6c41755c", "metadata": {}, "source": [ "### Method settings and workers\n", "\n", "Every keyword that `bootstrap` does not take goes to each fit: `alpha`, `indep_test`, `knowledge`,\n", "or any other argument of the method. A lower `alpha` asks each test for more evidence, so fewer\n", "pairs are adjacent in each fit:" ] }, { "cell_type": "code", "execution_count": 12, "id": "4019c0df", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "alpha 0.05: 32 pairs above 90%, 31 of them true\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "alpha 0.01: 30 pairs above 90%, 29 of them true\n" ] } ], "source": [ "for alpha in (0.05, 0.01):\n", " result = andrey.bootstrap(\"pc\", data, n_resamples=50, seed=1, alpha=alpha)\n", " adjacent = (1 - result.frequency(\"absent\"))[upper]\n", " print(f\"alpha {alpha}: {int((adjacent > 0.9).sum())} pairs above 90%,\", end=\" \")\n", " print(f\"{int(((adjacent > 0.9) & true_pair).sum())} of them true\")" ] }, { "cell_type": "markdown", "id": "6e7a5e3c", "metadata": {}, "source": [ "`andrey.config(n_jobs=4)` runs the fits in four processes. Starting the processes takes time, so it\n", "pays off when each fit takes a second or more; a PC fit on ALARM takes a few milliseconds, and one\n", "process is faster. In a script, run it under `if __name__ == \"__main__\":`, since each process\n", "imports the script again." ] }, { "cell_type": "code", "execution_count": 13, "id": "d3c5bfaf", "metadata": { "tags": [ "output-varies-by-cpu" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 process(es): 0.6 s\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "4 process(es): 2.0 s\n" ] } ], "source": [ "import time\n", "\n", "for n_jobs in (1, 4):\n", " with andrey.config(n_jobs=n_jobs):\n", " start = time.perf_counter()\n", " andrey.bootstrap(\"pc\", data, n_resamples=100, seed=1)\n", " print(f\"{n_jobs} process(es): {time.perf_counter() - start:.1f} s\")" ] }, { "cell_type": "markdown", "id": "1c8ce5ed", "metadata": {}, "source": [ "### Time series\n", "\n", "A time-series method resamples time, not independent rows. A subsample is one stretch of\n", "consecutive steps, half of the series by default, so every lag inside it is intact.\n", "`resampling=\"bootstrap\"` draws a series of the same length from blocks of consecutive steps, the\n", "moving block bootstrap, with `block_length` steps each, by default `max(ceil(n ** (1/3)), 10 *\n", "tau_max)`. Here PCMCI looks for the links at lag 1 in 1,000 steps of four variables:" ] }, { "cell_type": "code", "execution_count": 14, "id": "6b028712", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10 true links at lag 1\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "subsample, block length None:\n", " 0 -> 0 at lag 1: 1.00 (true)\n", " 0 -> 1 at lag 1: 1.00 (true)\n", " 1 -> 0 at lag 1: 1.00 (true)\n", " 1 -> 1 at lag 1: 1.00 (true)\n", " 2 -> 1 at lag 1: 1.00 (true)\n", " 2 -> 3 at lag 1: 1.00 (true)\n", " 2 -> 2 at lag 1: 0.98 (true)\n", " 1 -> 3 at lag 1: 0.84 (true)\n", " 3 -> 1 at lag 1: 0.20 (false)\n", " 0 -> 3 at lag 1: 0.12 (true)\n", "\n", "block bootstrap, block length 10:\n", " 0 -> 0 at lag 1: 1.00 (true)\n", " 1 -> 0 at lag 1: 1.00 (true)\n", " 1 -> 1 at lag 1: 1.00 (true)\n", " 2 -> 1 at lag 1: 1.00 (true)\n", " 2 -> 2 at lag 1: 1.00 (true)\n", " 2 -> 3 at lag 1: 1.00 (true)\n", " 1 -> 3 at lag 1: 0.96 (true)\n", " 0 -> 1 at lag 1: 0.92 (true)\n", " 0 -> 3 at lag 1: 0.64 (true)\n", " 3 -> 3 at lag 1: 0.44 (true)\n", " 3 -> 2 at lag 1: 0.16 (false)\n", " 3 -> 1 at lag 1: 0.14 (false)\n" ] } ], "source": [ "from andrey.data import sample_scm\n", "\n", "series = sample_scm(\n", " graph=\"erdos_renyi\",\n", " functional=\"var\",\n", " noise=\"uniform\",\n", " d=4,\n", " n=1000,\n", " seed=1,\n", " density=1.0,\n", " scale=\"raw\",\n", ")\n", "true_links = {(s, t) for s, t, _ in series.graph.lag(1).oriented_edges()}\n", "print(len(true_links), \"true links at lag 1\")\n", "for name, settings in {\"subsample\": {}, \"block bootstrap\": {\"resampling\": \"bootstrap\"}}.items():\n", " result = andrey.bootstrap(\"pcmci\", series.data, n_resamples=50, seed=1, tau_max=1, **settings)\n", " print(f\"\\n{name}, block length {result.metadata['block_length']}:\")\n", " for edge in result.table(min_frequency=0.1):\n", " link = (edge[\"source\"], edge[\"target\"])\n", " verdict = \"true\" if link in true_links else \"false\"\n", " print(f\" {link[0]} -> {link[1]} at lag 1: {edge['frequency']:.2f} ({verdict})\")" ] }, { "cell_type": "markdown", "id": "fb7998e6", "metadata": {}, "source": [ "Half of the series, 500 steps, misses the two weakest true links: 0 -> 3 comes back in 12% of the\n", "fits and 3 -> 3 in none. The block bootstrap's full-length resamples find them in 64% and 44% of\n", "its fits. In both, a false link appears in a fifth of the fits or fewer. On a short series, or when\n", "the weak links matter, a larger `sample_fraction` or the block bootstrap keeps them; a block must\n", "be long enough to hold the lags, which the default ensures." ] }, { "cell_type": "markdown", "id": "bootstrap-18", "metadata": {}, "source": [ "## Which edges to trust, and how to report them\n", "\n", "Trust the edges most fits agree on, and report each one with its frequency. Here are the averaged\n", "graph's most and least stable edges, each with the share of the fits in which the pair is adjacent\n", "and the share with this edge:" ] }, { "cell_type": "code", "execution_count": 15, "id": "bootstrap-19", "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
edgeadjacentthis edgetrue pair
0PAP -- PULMEMBOLUS1.000.90True
1HISTORY -- LVFAILURE1.001.00True
2CO -> BP1.000.60True
4TPR -> BP1.000.60True
26VENTTUBE -> VENTLUNG0.980.82True
32INTUBATION -> SHUNT0.980.92True
3CO -- STROKEVOLUME0.800.68True
14LVFAILURE -> LVEDVOLUME0.800.78True
\n", "
" ], "text/plain": [ " edge adjacent this edge true pair\n", "0 PAP -- PULMEMBOLUS 1.00 0.90 True\n", "1 HISTORY -- LVFAILURE 1.00 1.00 True\n", "2 CO -> BP 1.00 0.60 True\n", "4 TPR -> BP 1.00 0.60 True\n", "26 VENTTUBE -> VENTLUNG 0.98 0.82 True\n", "32 INTUBATION -> SHUNT 0.98 0.92 True\n", "3 CO -- STROKEVOLUME 0.80 0.68 True\n", "14 LVFAILURE -> LVEDVOLUME 0.80 0.78 True" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "position = {name: k for k, name in enumerate(data.columns)}\n", "glyph = {\"directed\": \"->\", \"undirected\": \"--\"}\n", "pairs = 1 - out.frequency(\"absent\")\n", "report = pd.DataFrame(\n", " {\n", " \"edge\": f\"{source} {glyph[kind]} {target}\",\n", " \"adjacent\": pairs[position[source], position[target]],\n", " \"this edge\": out.frequency(glyph[kind])[position[source], position[target]],\n", " \"true pair\": bool(truth.to_numpy()[position[source], position[target]]),\n", " }\n", " for source, target, kind in out.structure.oriented_edges()\n", ").sort_values(\"adjacent\", ascending=False, kind=\"stable\")\n", "pd.concat([report.head(4), report.tail(4)])" ] }, { "cell_type": "markdown", "id": "bootstrap-20", "metadata": {}, "source": [ "`threshold` asks more of the averaged graph: `out.average(0.9)` keeps only the pairs adjacent in\n", "more than 90% of the fits, without refitting, which here drops true edges as well as false ones.\n", "With each edge, report what produced the frequencies:\n", "the method and its settings, the number of resamples, the resampling and its share of rows, the\n", "seed, the threshold, and whether the averaged graph is the vote or the nearest resample." ] }, { "cell_type": "code", "execution_count": 16, "id": "bootstrap-21", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "32 edges, 25 errors\n", "{'method': 'pc', 'n_resamples': 50, 'threshold': 0.5, 'resampling': 'subsample', 'sample_fraction': 0.5, 'seed': 1, 'averaged_graph': 'vote'}\n" ] } ], "source": [ "strong = out.average(0.9)\n", "print(len(strong.oriented_edges()), \"edges,\", shd(strong, true_cpdag), \"errors\")\n", "print(\n", " {\n", " \"method\": out.metadata[\"method\"],\n", " \"n_resamples\": out.n_resamples,\n", " \"threshold\": out.threshold,\n", " **{\n", " key: out.metadata[key]\n", " for key in (\"resampling\", \"sample_fraction\", \"seed\", \"averaged_graph\")\n", " },\n", " }\n", ")" ] } ], "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 }