{ "cells": [ { "cell_type": "markdown", "id": "functional-0", "metadata": {}, "source": "# Nonlinear additive-noise models (CAM and RESIT)\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/functional.ipynb) {download}`Download this notebook <../../../examples/functional.ipynb>`\n\nCAM, the causal additive model, learns a DAG from data in which every variable is a nonlinear\nfunction of its parents plus independent noise, `x_j = f_j(parents) + e_j`. When the functions are\nnonlinear, the direction of every edge is identifiable, so CAM returns a DAG instead of the\nequivalence class that PC and GES return. It fits a smooth function of each parent, adds edges\ngreedily by the likelihood, and drops the parents whose smooth terms are not significant. It also\nreturns a causal order.\n\nRESIT, regression with subsequent independence tests, fits the same kind of model with noise of\nany distribution, for a handful of variables." }, { "cell_type": "markdown", "id": "functional-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": "functional-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": "functional-3", "metadata": {}, "outputs": [], "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.data import sample_scm\n", "from andrey.metrics import shd" ] }, { "cell_type": "markdown", "id": "functional-4", "metadata": {}, "source": [ "## Data\n", "\n", "`sample_scm` draws six variables on an Erdos-Renyi graph, each the sum of a nonlinear function of\n", "its parents and Gaussian noise, with 1,000 rows and standardized columns. The columns come in an\n", "order drawn from the seed, so their order says nothing about the causal order." ] }, { "cell_type": "code", "execution_count": 3, "id": "functional-5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 observations of 6 variables\n" ] }, { "data": { "text/html": [ "012345" ], "text/plain": [ "'012345'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def draw(seed: int, d: int = 6, noise: str = \"gaussian\"):\n", " \"\"\"Draw 1,000 rows of d variables with additive noise; return the data and the true DAG.\"\"\"\n", " sample = sample_scm(\n", " graph=\"erdos_renyi\",\n", " functional=\"additive_noise\",\n", " noise=noise,\n", " d=d,\n", " n=1000,\n", " seed=seed,\n", " scale=\"standardize\",\n", " density=2.0,\n", " )\n", " return sample.data, sample.graph\n", "\n", "\n", "X, truth = draw(seed=0)\n", "print(f\"{X.shape[0]} observations of {X.shape[1]} variables\")\n", "display(andrey.viz.draw(truth))" ] }, { "cell_type": "markdown", "id": "functional-6", "metadata": {}, "source": [ "## Fit CAM\n", "\n", "CAM returns a DAG and `ordering`, a causal order as column indices, causes first: the order of its\n", "search before pruning." ] }, { "cell_type": "code", "execution_count": 4, "id": "functional-7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "causal order (5, 2, 4, 3, 0, 1)\n", " 0 -> 1\n", " 3 -> 0\n", " 3 -> 1\n", " 4 -> 1\n", " 5 -> 2\n", " 4 -> 3\n" ] }, { "data": { "text/html": [ "012345" ], "text/plain": [ "'012345'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cam = andrey.cam(X)\n", "print(f\"causal order {cam.ordering}\")\n", "for source, target, _ in cam.structure.oriented_edges():\n", " print(f\" {source} -> {target}\")\n", "display(andrey.viz.draw(cam.structure))" ] }, { "cell_type": "markdown", "id": "functional-8", "metadata": {}, "source": [ "## Compare with the truth\n", "\n", "The structural Hamming distance (SHD) counts the edges to add, remove, or reverse to reach the true\n", "graph; the empty graph's SHD is the number of true edges." ] }, { "cell_type": "code", "execution_count": 5, "id": "functional-9", "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", "
edgesSHD
CAM63
empty graph06
\n", "
" ], "text/plain": [ " edges SHD\n", "CAM 6 3\n", "empty graph 0 6" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows = {\n", " \"CAM\": {\n", " \"edges\": len(cam.structure.oriented_edges()),\n", " \"SHD\": shd(cam.structure, truth),\n", " },\n", " \"empty graph\": {\"edges\": 0, \"SHD\": len(truth.oriented_edges())},\n", "}\n", "pd.DataFrame(rows).T" ] }, { "cell_type": "markdown", "id": "functional-10", "metadata": {}, "source": [ "## Check over ten datasets\n", "\n", "The same model, drawn with ten seeds. The table gives CAM's SHD, with the empty graph's beside it." ] }, { "cell_type": "code", "execution_count": 6, "id": "functional-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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CAMempty graph
seed
036
106
206
306
416
506
626
736
806
966
\n", "
" ], "text/plain": [ " CAM empty graph\n", "seed \n", "0 3 6\n", "1 0 6\n", "2 0 6\n", "3 0 6\n", "4 1 6\n", "5 0 6\n", "6 2 6\n", "7 3 6\n", "8 0 6\n", "9 6 6" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "median SHD: {'CAM': 0.5, 'empty graph': 6.0}\n" ] } ], "source": [ "rows = {}\n", "for seed in range(10):\n", " data, graph = draw(seed)\n", " rows[seed] = {\n", " \"CAM\": shd(andrey.cam(data).structure, graph),\n", " \"empty graph\": len(graph.oriented_edges()),\n", " }\n", "table = pd.DataFrame(rows).T\n", "table.index.name = \"seed\"\n", "display(table)\n", "print(\"median SHD:\", table.median().to_dict())" ] }, { "cell_type": "markdown", "id": "functional-12", "metadata": {}, "source": [ "## Rescaling changes nothing\n", "\n", "CAM does not depend on a column's scale or location, so multiplying or shifting a column leaves its\n", "graph unchanged. Here every column is rescaled and shifted differently." ] }, { "cell_type": "code", "execution_count": 7, "id": "functional-13", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "same graph on the rescaled data: True\n" ] } ], "source": [ "scale = np.array([0.01, 1.0, 5.0, 50.0, 300.0, 1e4])\n", "unit_free = X * scale + 7 * scale\n", "same = np.array_equal(andrey.cam(unit_free).structure.to_numpy(), cam.structure.to_numpy())\n", "print(f\"same graph on the rescaled data: {same}\")" ] }, { "cell_type": "markdown", "id": "functional-115", "metadata": {}, "source": [ "## RESIT\n", "\n", "RESIT fits the same model with noise of any distribution. It finds a causal order from the sink up:\n", "the next sink is the variable whose additive-model residual an independence test finds most\n", "independent of the remaining variables. Then it removes each parent whose removal leaves the residual\n", "independent. Each independence test reads all the remaining variables at once, so RESIT suits a\n", "handful of variables: its accuracy claim holds up to 10, and it takes seconds at 5 and minutes at 30.\n", "Here it runs on five variables with Gumbel noise." ] }, { "cell_type": "code", "execution_count": 8, "id": "functional-116", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "causal order (1, 2, 0, 3, 4)\n", " 1 -> 0\n", " 1 -> 2\n", " 2 -> 3\n", " 2 -> 4\n", " 3 -> 4\n", "SHD 0, empty graph 5\n" ] } ], "source": [ "X5, truth5 = draw(seed=0, d=5, noise=\"gumbel\")\n", "resit = andrey.resit(X5)\n", "print(f\"causal order {resit.ordering}\")\n", "for source, target, _ in resit.structure.oriented_edges():\n", " print(f\" {source} -> {target}\")\n", "print(f\"SHD {shd(resit.structure, truth5)}, empty graph {len(truth5.oriented_edges())}\")" ] }, { "cell_type": "markdown", "id": "functional-117", "metadata": {}, "source": [ "## Check RESIT over ten datasets\n", "\n", "The same model with ten seeds: RESIT's SHD, with the empty graph's beside it." ] }, { "cell_type": "code", "execution_count": 9, "id": "functional-118", "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", "
RESITempty graph
seed
005
105
205
305
405
505
625
705
805
905
\n", "
" ], "text/plain": [ " RESIT empty graph\n", "seed \n", "0 0 5\n", "1 0 5\n", "2 0 5\n", "3 0 5\n", "4 0 5\n", "5 0 5\n", "6 2 5\n", "7 0 5\n", "8 0 5\n", "9 0 5" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "median SHD: {'RESIT': 0.0, 'empty graph': 5.0}\n" ] } ], "source": [ "rows = {}\n", "for seed in range(10):\n", " data, graph = draw(seed, d=5, noise=\"gumbel\")\n", " rows[seed] = {\n", " \"RESIT\": shd(andrey.resit(data).structure, graph),\n", " \"empty graph\": len(graph.oriented_edges()),\n", " }\n", "table = pd.DataFrame(rows).T\n", "table.index.name = \"seed\"\n", "display(table)\n", "print(\"median SHD:\", table.median().to_dict())" ] }, { "cell_type": "markdown", "id": "functional-14", "metadata": {}, "source": [ "## When to use it\n", "\n", "CAM needs effects that are nonlinear and additive in the parents, with additive noise and no hidden\n", "common causes. With linear effects and Gaussian noise the data identify only the equivalence class\n", "and the direction of the other edges is arbitrary, so use PC or GES; noise whose scale depends on\n", "the parents can reverse edges. The [0.13.0 release card](https://andrey.abel.ai/docs/benchmarks/0.13.0.html#cam)\n", "compares CAM with CAM 1.0 from 10 to 100 variables. Change `seed`, `d`, or the noise (`\"gumbel\"`)\n", "and rerun to see how it responds.\n", "\n", "RESIT needs the same kind of data with noise of any distribution, and no more than about ten\n", "variables: from 20 variables its errors come near the empty graph's on Gaussian noise and a fit\n", "takes about a minute, so use CAM for more variables. The release card compares RESIT with the\n", "authors' code and with lingam from 5 to 20 variables." ] } ], "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 }