{ "cells": [ { "cell_type": "markdown", "id": "lingam-0", "metadata": {}, "source": [ "# Linear non-Gaussian models (DirectLiNGAM / ICA-LiNGAM / MultiGroupDirectLiNGAM / GIN)\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/lingam.ipynb)\n", "\n", "When the mechanisms are linear and the noise is non-Gaussian, the data identify a single DAG, not\n", "only its equivalence class. DirectLiNGAM finds a causal order one variable at a time with a\n", "pairwise independence measure. ICA-LiNGAM runs independent component analysis (FastICA) and reads\n", "the order from the unmixing matrix. Both return a fully oriented DAG, a causal order, and edge\n", "weights. MultiGroupDirectLiNGAM fits several datasets that share one causal order, and GIN finds\n", "the latent variables that observed indicators measure." ] }, { "cell_type": "markdown", "id": "7fb27b941602401d91542211134fc71a", "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": "acae54e37e7d407bbb7b55eff062a284", "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": "code", "execution_count": 2, "id": "lingam-1", "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.metrics as metrics\n", "import andrey.viz\n", "from andrey.data import benchmarks" ] }, { "cell_type": "markdown", "id": "lingam-2", "metadata": {}, "source": [ "## Data\n", "\n", "`lingam_sf` draws a scale-free graph with linear mechanisms and uniform noise: here 8 variables\n", "and 2,000 samples." ] }, { "cell_type": "code", "execution_count": 3, "id": "lingam-3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2000 observations, 8 variables, seed 7\n" ] }, { "data": { "text/html": [ "01234567" ], "text/plain": [ "'01234567'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "seed = 7\n", "dataset = benchmarks.scm(\"lingam_sf\", d=8).sample(n=2000, seed=seed, scale=\"standardize\")\n", "X, truth = dataset.data, dataset.observed_graph\n", "print(f\"{X.shape[0]} observations, {X.shape[1]} variables, seed {seed}\")\n", "display(andrey.viz.draw(truth))" ] }, { "cell_type": "markdown", "id": "lingam-4", "metadata": {}, "source": [ "## Fit both methods\n", "\n", "DirectLiNGAM is deterministic. ICA-LiNGAM starts FastICA from a random point; `random_state` fixes\n", "it, so a run repeats exactly." ] }, { "cell_type": "code", "execution_count": 4, "id": "lingam-5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DirectLiNGAM: dag, causal order (7, 5, 1, 3, 6, 0, 4, 2)\n" ] }, { "data": { "text/html": [ "01234567" ], "text/plain": [ "'01234567'" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "ICA-LiNGAM: dag, causal order (5, 7, 1, 3, 6, 0, 2, 4)\n" ] }, { "data": { "text/html": [ "01234567" ], "text/plain": [ "'01234567'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "results = {\n", " \"DirectLiNGAM\": andrey.direct_lingam(X),\n", " \"ICA-LiNGAM\": andrey.ica_lingam(X, random_state=0),\n", "}\n", "for name, result in results.items():\n", " print(f\"{name}: {result.structure.kind}, causal order {result.ordering}\")\n", " display(andrey.viz.draw(result.structure))" ] }, { "cell_type": "markdown", "id": "lingam-6", "metadata": {}, "source": [ "## Edge weights\n", "\n", "`weighted_adjacency[i, j]` is the coefficient of variable `i` in the equation for variable `j`, so\n", "each nonzero entry is an edge `i -> j`." ] }, { "cell_type": "code", "execution_count": 5, "id": "lingam-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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01234567
00.000.00-0.720.000.510.00.000.0
10.000.000.000.94-0.710.00.000.0
20.000.000.000.000.000.00.000.0
30.000.000.000.000.000.00.760.0
40.000.000.000.000.000.00.000.0
50.00-0.720.000.380.000.00.660.0
60.790.000.470.000.000.00.000.0
70.000.00-0.770.380.000.00.000.0
\n", "
" ], "text/plain": [ " 0 1 2 3 4 5 6 7\n", "0 0.00 0.00 -0.72 0.00 0.51 0.0 0.00 0.0\n", "1 0.00 0.00 0.00 0.94 -0.71 0.0 0.00 0.0\n", "2 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.0\n", "3 0.00 0.00 0.00 0.00 0.00 0.0 0.76 0.0\n", "4 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.0\n", "5 0.00 -0.72 0.00 0.38 0.00 0.0 0.66 0.0\n", "6 0.79 0.00 0.47 0.00 0.00 0.0 0.00 0.0\n", "7 0.00 0.00 -0.77 0.38 0.00 0.0 0.00 0.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.DataFrame(results[\"ICA-LiNGAM\"].weighted_adjacency).round(2)" ] }, { "cell_type": "markdown", "id": "lingam-8", "metadata": {}, "source": [ "## Check accuracy over ten datasets\n", "\n", "SHD counts errors against the true DAG; lower is better. The empty graph's SHD is the number of\n", "true edges. The table covers ten seeds and their median." ] }, { "cell_type": "code", "execution_count": 6, "id": "lingam-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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
seed0123456789median
Empty graph1212121212121212121212.0
DirectLiNGAM00000000000.0
ICA-LiNGAM08003000000.0
\n", "
" ], "text/plain": [ "seed 0 1 2 3 4 5 6 7 8 9 median\n", "Empty graph 12 12 12 12 12 12 12 12 12 12 12.0\n", "DirectLiNGAM 0 0 0 0 0 0 0 0 0 0 0.0\n", "ICA-LiNGAM 0 8 0 0 3 0 0 0 0 0 0.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "methods = {\n", " \"DirectLiNGAM\": lambda X: andrey.direct_lingam(X),\n", " \"ICA-LiNGAM\": lambda X: andrey.ica_lingam(X, random_state=0),\n", "}\n", "rows = {}\n", "for s in range(10):\n", " draw = benchmarks.scm(\"lingam_sf\", d=8).sample(n=2000, seed=s, scale=\"standardize\")\n", " d = draw.observed_graph.n_nodes\n", " empty = andrey.GraphStructure.from_numpy(np.zeros((d, d), dtype=np.uint8), kind=\"dag\")\n", " rows[s] = {\"Empty graph\": metrics.score(empty, draw.observed_graph)[\"shd\"]}\n", " for name, fit in methods.items():\n", " rows[s][name] = metrics.score(fit(draw.data), draw.observed_graph)[\"shd\"]\n", "shd = pd.DataFrame(rows)\n", "shd.columns.name = \"seed\"\n", "shd[\"median\"] = shd.median(axis=1)\n", "shd" ] }, { "cell_type": "markdown", "id": "lingam-10", "metadata": {}, "source": [ "## When to use them\n", "\n", "Both methods rely on non-Gaussian noise. With Gaussian noise the order is not identified; PC, GES,\n", "BOSS, or GRaSP return the CPDAG that the data do support. ICA-LiNGAM needs many samples: at\n", "`n = 10d`, the size the [release cards](https://andrey.abel.ai/docs/benchmarks.html) use, it scores\n", "worse than the empty graph, and this notebook uses `n = 250d`. Change `seed`, `d`, or `n` and rerun\n", "to see how each method responds." ] }, { "cell_type": "markdown", "id": "lingam-11", "metadata": {}, "source": [ "## Several datasets: MultiGroupDirectLiNGAM\n", "\n", "`multi_group_direct_lingam` takes datasets over the same variables that share one causal order\n", "and returns one result per dataset, each with its own weights. Here two groups follow the chain\n", "`x0 -> x1 -> x2`, with different weights on `x0 -> x1`, their columns in one order drawn from the\n", "seed." ] }, { "cell_type": "code", "execution_count": 7, "id": "lingam-12", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "group 0: [('x0', 'x1', 'directed'), ('x1', 'x2', 'directed')], weight 0.74\n", "group 1: [('x0', 'x1', 'directed'), ('x1', 'x2', 'directed')], weight -1.17\n" ] } ], "source": [ "CHAIN = [(\"x0\", \"x1\"), (\"x1\", \"x2\")] # the true edges, cause first\n", "\n", "\n", "def sample_groups(seed: int, n: int = 1000, shuffle: bool = True) -> list[pd.DataFrame]:\n", " \"\"\"Draw n rows of the chain per weight on x0 -> x1, all in one column order from the seed.\"\"\"\n", " rng = np.random.default_rng(seed)\n", " groups = []\n", " for weight in (0.8, -1.2):\n", " x0 = rng.uniform(-1, 1, n)\n", " x1 = weight * x0 + rng.uniform(-1, 1, n)\n", " x2 = 0.6 * x1 + rng.uniform(-1, 1, n)\n", " groups.append(pd.DataFrame({\"x0\": x0, \"x1\": x1, \"x2\": x2}))\n", " columns = rng.permutation(groups[0].columns) if shuffle else groups[0].columns\n", " return [frame[columns] for frame in groups]\n", "\n", "\n", "groups = sample_groups(seed=3)\n", "for g, out in enumerate(andrey.multi_group_direct_lingam(groups)):\n", " weights = pd.DataFrame(\n", " out.weighted_adjacency, index=groups[g].columns, columns=groups[g].columns\n", " )\n", " print(f\"group {g}: {out.structure.oriented_edges()}, weight {weights.loc['x0', 'x1']:.2f}\")" ] }, { "cell_type": "markdown", "id": "lingam-13", "metadata": {}, "source": [ "## Latent variables: GIN\n", "\n", "`gin` groups observed indicators under the latent variables that cause them and orders the\n", "latents. Here `a1` to `a3` measure one latent and `b1` to `b3` another, which the first causes. The\n", "indicators come in an order drawn from the seed." ] }, { "cell_type": "code", "execution_count": 8, "id": "lingam-14", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[('L1', 'a1', 'directed'), ('L2', 'b1', 'directed'), ('L2', 'b2', 'directed'), ('L1', 'a2', 'directed'), ('L1', 'a3', 'directed'), ('L2', 'b3', 'directed'), ('L1', 'L2', 'directed')]\n" ] }, { "data": { "text/html": [ "a1b1b2a2a3b3L1L2" ], "text/plain": [ "'a1b1b2a2a3b3L1L2'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "INDICATORS = [(\"l1\", \"l2\")] # the latents l1 -> l2 have no column\n", "INDICATORS += [(\"l1\", f\"a{i}\") for i in (1, 2, 3)] + [(\"l2\", f\"b{i}\") for i in (1, 2, 3)]\n", "\n", "\n", "def sample_indicators(seed: int, n: int = 1000, shuffle: bool = True) -> pd.DataFrame:\n", " \"\"\"Draw n rows of the six indicators of l1 -> l2, in an order drawn from the seed.\"\"\"\n", " rng = np.random.default_rng(seed)\n", " l1 = rng.uniform(-1, 1, n)\n", " l2 = 0.9 * l1 + rng.uniform(-1, 1, n)\n", " loadings = {\"a1\": (l1, 1.0), \"a2\": (l1, 0.8), \"a3\": (l1, -1.2)}\n", " loadings |= {\"b1\": (l2, 1.0), \"b2\": (l2, 1.1), \"b3\": (l2, -0.7)}\n", " frame = pd.DataFrame(\n", " {name: w * latent + 0.3 * rng.uniform(-1, 1, n) for name, (latent, w) in loadings.items()}\n", " )\n", " return frame[rng.permutation(frame.columns)] if shuffle else frame\n", "\n", "\n", "gin_out = andrey.gin(sample_indicators(seed=5))\n", "print(gin_out.structure.oriented_edges())\n", "display(andrey.viz.draw(gin_out.structure))" ] } ], "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 }