{ "cells": [ { "cell_type": "markdown", "id": "rcd-0", "metadata": {}, "source": "# Hidden common causes (RCD)\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/rcd.ipynb) {download}`Download this notebook <../../../examples/rcd.ipynb>`\n\nWhen two variables share a cause that is not in the data, most methods either join them with a\nwrong edge or miss the dependence. RCD, repetitive causal discovery, assumes linear mechanisms with\nnon-Gaussian noise. It regresses each variable on candidate ancestors and tests the residuals for\nindependence, which finds the direct causes and the pairs that share a hidden cause, marked `<->`." }, { "cell_type": "markdown", "id": "rcd-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": "rcd-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": "rcd-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": "rcd-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": "rcd-5", "metadata": {}, "source": [ "## Fit RCD\n", "\n", "`andrey.rcd` returns a PAG over the observed variables: `->` for a direct cause and `<->` for a pair\n", "with a hidden common cause. `metadata[\"confounded_pairs\"]` lists the `<->` pairs by column index." ] }, { "cell_type": "code", "execution_count": 3, "id": "rcd-6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "confounded pairs: [[0, 1]] = [['x0', 'x1']]\n", "x0 <-> x1\n", "x0 -> x3\n", "x1 -> x2\n" ] }, { "data": { "text/html": [ "x0x1x3x2" ], "text/plain": [ "'x0x1x3x2'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "out = andrey.rcd(X)\n", "pairs = [[X.columns[i] for i in pair] for pair in out.metadata[\"confounded_pairs\"]]\n", "print(\"confounded pairs:\", out.metadata[\"confounded_pairs\"], \"=\", pairs)\n", "for source, target, kind in out.structure.oriented_edges():\n", " print(f\"{source} {'<->' if kind == 'bidirected' else '->'} {target}\")\n", "display(andrey.viz.draw(out.structure))" ] }, { "cell_type": "markdown", "id": "rcd-7", "metadata": {}, "source": [ "## Edge weights\n", "\n", "`weighted_adjacency[i, j]` is the coefficient of variable `i` in the equation for variable `j`. A\n", "pair marked `<->` has no weight: its dependence comes from the hidden cause." ] }, { "cell_type": "code", "execution_count": 4, "id": "rcd-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": "rcd-9", "metadata": {}, "source": [ "## Check over ten datasets\n", "\n", "The same model, drawn with ten seeds. For each, the table says whether RCD found the hidden pair\n", "`x0 <-> x1` and both direct edges, and how many edges it reported." ] }, { "cell_type": "code", "execution_count": 5, "id": "rcd-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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
hidden pairdirect edgesedgessame as the truth
seed
0TrueTrue3True
1TrueTrue3True
2TrueTrue3True
3TrueTrue3True
4TrueTrue3True
5TrueTrue3True
6TrueTrue3True
7TrueTrue3True
8TrueTrue3True
9TrueTrue3True
\n", "
" ], "text/plain": [ " hidden pair direct edges edges same as the truth\n", "seed \n", "0 True True 3 True\n", "1 True True 3 True\n", "2 True True 3 True\n", "3 True True 3 True\n", "4 True True 3 True\n", "5 True True 3 True\n", "6 True True 3 True\n", "7 True True 3 True\n", "8 True True 3 True\n", "9 True True 3 True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def named_edges(structure) -> set:\n", " \"\"\"The edges by name; a bidirected pair in alphabetical order, whatever the column order.\"\"\"\n", " return {\n", " (*sorted((a, b)), kind) if kind == \"bidirected\" else (a, b, kind)\n", " for a, b, kind in structure.oriented_edges()\n", " }\n", "\n", "\n", "expected = named_edges(truth)\n", "rows = {}\n", "for seed in range(10):\n", " found = named_edges(andrey.rcd(sample(seed)).structure)\n", " rows[seed] = {\n", " \"hidden pair\": (\"x0\", \"x1\", \"bidirected\") in found,\n", " \"direct edges\": {(\"x0\", \"x3\", \"directed\"), (\"x1\", \"x2\", \"directed\")} <= found,\n", " \"edges\": len(found),\n", " \"same as the truth\": found == expected,\n", " }\n", "pd.DataFrame(rows).T.rename_axis(\"seed\")" ] }, { "cell_type": "markdown", "id": "rcd-11", "metadata": {}, "source": [ "## When to use it\n", "\n", "RCD relies on linear mechanisms and non-Gaussian noise, and it needs many rows: a weak hidden\n", "cause leaves only a small trace in the residuals it tests. Without hidden causes, DirectLiNGAM\n", "returns a full DAG; with Gaussian noise, FCI returns the PAG the conditional independences support.\n", "The [0.2.0 release card](https://andrey.abel.ai/docs/benchmarks/0.2.0.html) measures RCD on larger\n", "graphs. Change the seed, the hidden effect (`1.5`), or `n` and rerun to see how RCD 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 }