{ "cells": [ { "cell_type": "markdown", "id": "bottom_up_parce_lingam-0", "metadata": {}, "source": [ "# Hidden common causes (BottomUpParceLiNGAM)\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/bottom_up_parce_lingam.ipynb)\n", "\n", "BottomUpParceLiNGAM assumes linear mechanisms with non-Gaussian noise and finds a causal order from\n", "the bottom up: at each step it looks for a sink, a variable whose regression residual on the others\n", "is independent of them. A hidden common cause leaves the variables it confounds without a sink, so\n", "the method stops there and leaves them unordered, joined by `o-o`, rather than guess an order." ] }, { "cell_type": "markdown", "id": "bottom_up_parce_lingam-1", "metadata": {}, "source": [ "## Start\n", "\n", "Choose a CPU runtime and run this cell first. In Colab it installs `andrey-core==0.9.0` from PyPI;\n", "local kernels use the installed environment." ] }, { "cell_type": "code", "execution_count": 1, "id": "bottom_up_parce_lingam-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.9.0\" pandas" ] }, { "cell_type": "markdown", "id": "bottom_up_parce_lingam-3", "metadata": {}, "source": [ "## Data\n", "\n", "A linear model with uniform noise over five variables, one of them hidden: `h` drives `x0` and\n", "`x1`, `x0` drives `x3`, and `x1` drives `x2`. The data keep only `x0` to `x3`, so `x0` and `x1`\n", "depend on each other with no edge between them. The true graph over the observed variables marks\n", "that pair `<->`. The columns come in an order drawn from the seed, so their order says nothing\n", "about the causal order." ] }, { "cell_type": "code", "execution_count": 2, "id": "bottom_up_parce_lingam-4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2000 observations, 4 variables, columns ['x0', 'x1', 'x3', 'x2']\n" ] }, { "data": { "text/html": [ "x0x1x2x3" ], "text/plain": [ "'x0x1x2x3'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import pandas as pd\n", "from IPython.display import display\n", "\n", "import andrey\n", "import andrey.viz\n", "from andrey.core import ARROW, TAIL\n", "\n", "COLUMNS = [\"x0\", \"x1\", \"x2\", \"x3\"]\n", "EDGES = [(\"h\", \"x0\"), (\"h\", \"x1\"), (\"x0\", \"x3\"), (\"x1\", \"x2\")] # h has no column\n", "\n", "\n", "def sample(seed: int, n: int = 2000, shuffle: bool = True) -> pd.DataFrame:\n", " \"\"\"Draw n rows of x0..x3, the columns in an order drawn from the seed; h is left out.\"\"\"\n", " rng = np.random.default_rng(seed)\n", " h, e0, e1, e2, e3 = rng.uniform(-1, 1, (5, n))\n", " x0, x1 = 1.5 * h + e0, 1.5 * h + e1\n", " frame = pd.DataFrame({\"x0\": x0, \"x1\": x1, \"x2\": 0.6 * x1 + e2, \"x3\": 0.6 * x0 + e3})\n", " return frame[rng.permutation(frame.columns)] if shuffle else frame\n", "\n", "\n", "# marks[i, j] is the mark at x_i of the edge between x_i and x_j.\n", "marks = np.zeros((4, 4), dtype=np.int8)\n", "marks[0, 1] = marks[1, 0] = ARROW # x0 <-> x1\n", "marks[0, 3], marks[3, 0] = TAIL, ARROW # x0 -> x3\n", "marks[1, 2], marks[2, 1] = TAIL, ARROW # x1 -> x2\n", "truth = andrey.GraphStructure.from_numpy(marks, kind=\"pag\", labels=tuple(COLUMNS))\n", "\n", "X = sample(seed=0)\n", "print(f\"{X.shape[0]} observations, {X.shape[1]} variables, columns {list(X.columns)}\")\n", "display(andrey.viz.draw(truth))" ] }, { "cell_type": "markdown", "id": "bottom_up_parce_lingam-5", "metadata": {}, "source": [ "## Fit BottomUpParceLiNGAM\n", "\n", "`andrey.bottom_up_parce_lingam` returns a PAG: `->` for a direct cause and `o-o` between variables\n", "it leaves unordered. `metadata[\"causal_order\"]` gives the order, the unordered variables nested in\n", "one list. `alpha` is the significance level of the sink test, Bonferroni-corrected across the\n", "remaining variables; at the default 0.1 a true sink is sometimes rejected by chance on this small\n", "model, so this notebook uses 0.01." ] }, { "cell_type": "code", "execution_count": 3, "id": "bottom_up_parce_lingam-6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "causal order: [[0, 1], 2, 3] over columns ['x0', 'x1', 'x3', 'x2']\n", "x0 o-o x1\n", "x0 -> x3\n", "x1 -> x2\n" ] }, { "data": { "text/html": [ "x0x1x3x2" ], "text/plain": [ "'x0x1x3x2'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "out = andrey.bottom_up_parce_lingam(X, alpha=0.01)\n", "print(\"causal order:\", out.metadata[\"causal_order\"], \"over columns\", list(X.columns))\n", "for source, target, kind in out.structure.oriented_edges():\n", " print(f\"{source} {'->' if kind == 'directed' else 'o-o'} {target}\")\n", "display(andrey.viz.draw(out.structure))" ] }, { "cell_type": "markdown", "id": "bottom_up_parce_lingam-7", "metadata": {}, "source": [ "## Edge weights\n", "\n", "`weighted_adjacency[i, j]` is the coefficient of variable `i` in the equation for variable `j`,\n", "from an adaptive-Lasso regression of each variable on those before it in the order." ] }, { "cell_type": "code", "execution_count": 4, "id": "bottom_up_parce_lingam-8", "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", "
x0x1x2x3
x00.00.00.000.6
x10.00.00.610.0
x20.00.00.000.0
x30.00.00.000.0
\n", "
" ], "text/plain": [ " x0 x1 x2 x3\n", "x0 0.0 0.0 0.00 0.6\n", "x1 0.0 0.0 0.61 0.0\n", "x2 0.0 0.0 0.00 0.0\n", "x3 0.0 0.0 0.00 0.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weights = pd.DataFrame(out.weighted_adjacency, index=X.columns, columns=X.columns)\n", "weights.loc[COLUMNS, COLUMNS].round(2)" ] }, { "cell_type": "markdown", "id": "bottom_up_parce_lingam-9", "metadata": {}, "source": [ "## Check over ten datasets\n", "\n", "The same model, drawn with ten seeds. For each, the table says whether the method left exactly\n", "`x0` and `x1` unordered, whether it found both direct edges, and its causal order." ] }, { "cell_type": "code", "execution_count": 5, "id": "bottom_up_parce_lingam-10", "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0, x1 unordereddirect edgescausal order
seed
0TrueTrue[['x0', 'x1'], 'x3', 'x2']
1TrueTrue[['x0', 'x1'], 'x2', 'x3']
2TrueTrue[['x1', 'x0'], 'x2', 'x3']
3TrueTrue[['x1', 'x0'], 'x3', 'x2']
4TrueTrue[['x0', 'x1'], 'x2', 'x3']
5TrueTrue[['x1', 'x0'], 'x2', 'x3']
6TrueTrue[['x1', 'x0'], 'x3', 'x2']
7TrueTrue[['x0', 'x1'], 'x2', 'x3']
8TrueTrue[['x0', 'x1'], 'x2', 'x3']
9TrueTrue[['x1', 'x0'], 'x2', 'x3']
\n", "
" ], "text/plain": [ " x0, x1 unordered direct edges causal order\n", "seed \n", "0 True True [['x0', 'x1'], 'x3', 'x2']\n", "1 True True [['x0', 'x1'], 'x2', 'x3']\n", "2 True True [['x1', 'x0'], 'x2', 'x3']\n", "3 True True [['x1', 'x0'], 'x3', 'x2']\n", "4 True True [['x0', 'x1'], 'x2', 'x3']\n", "5 True True [['x1', 'x0'], 'x2', 'x3']\n", "6 True True [['x1', 'x0'], 'x3', 'x2']\n", "7 True True [['x0', 'x1'], 'x2', 'x3']\n", "8 True True [['x0', 'x1'], 'x2', 'x3']\n", "9 True True [['x1', 'x0'], 'x2', 'x3']" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def by_name(order: list, columns) -> list:\n", " \"\"\"The causal order with each column index replaced by its name.\"\"\"\n", " return [\n", " [columns[i] for i in item] if isinstance(item, list) else columns[item] for item in order\n", " ]\n", "\n", "\n", "rows = {}\n", "for seed in range(10):\n", " X_seed = sample(seed)\n", " fit = andrey.bottom_up_parce_lingam(X_seed, alpha=0.01)\n", " order = by_name(fit.metadata[\"causal_order\"], X_seed.columns)\n", " blocks = [set(item) for item in order if isinstance(item, list)]\n", " found = set(fit.structure.oriented_edges())\n", " rows[seed] = {\n", " \"x0, x1 unordered\": blocks == [{\"x0\", \"x1\"}],\n", " \"direct edges\": {(\"x0\", \"x3\", \"directed\"), (\"x1\", \"x2\", \"directed\")} <= found,\n", " \"causal order\": str(order),\n", " }\n", "pd.DataFrame(rows).T.rename_axis(\"seed\")" ] }, { "cell_type": "markdown", "id": "bottom_up_parce_lingam-11", "metadata": {}, "source": [ "## When to use it\n", "\n", "BottomUpParceLiNGAM relies on linear mechanisms and non-Gaussian noise. Its unordered block holds\n", "the confounded variables and every variable above them in the order, so on a larger graph it can\n", "join pairs that share no cause, and edges inside the block stay unoriented. RCD marks only the\n", "pairs it attributes to a hidden cause; without hidden causes, DirectLiNGAM returns a full DAG. The\n", "[0.2.0 release card](https://andrey.abel.ai/docs/benchmarks/0.2.0.html) measures it on larger\n", "graphs. Change the seed, `alpha`, or `n` and rerun to see how it responds." ] } ], "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.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }