{ "cells": [ { "cell_type": "markdown", "id": "time_series-0", "metadata": {}, "source": [ "# Time series (PCMCI, PCMCI+, LPCMCI, Granger, VARMA-LiNGAM, longitudinal LiNGAM)\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/time_series.ipynb)\n", "\n", "Six methods find causes in data measured over time. PCMCI and Granger find the links from earlier\n", "steps; PCMCI+ also finds links within one step, and LPCMCI allows series that are not measured.\n", "VARMA-LiNGAM estimates the weights of the links within one step and from the step before.\n", "Longitudinal LiNGAM takes a panel instead: the same units measured at a few time points. Each\n", "returns a `TemporalStructure`, one graph per lag." ] }, { "cell_type": "markdown", "id": "time_series-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": "time_series-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": "code", "execution_count": 2, "id": "time_series-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" ] }, { "cell_type": "markdown", "id": "time_series-4", "metadata": {}, "source": [ "## Data\n", "\n", "`sample_scm` draws one series of 1,000 time steps over 4 variables: a vector autoregression with\n", "one lag, links within each step from a random DAG, and uniform noise. The drawing merges the lags\n", "into one graph; a loop is a variable's dependence on its own past." ] }, { "cell_type": "code", "execution_count": 3, "id": "time_series-5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 time steps of 4 variables\n", "temporal | 4 nodes | 2 lags\n", "\n", " lag 0 | 2 edges (2 directed)\n", " 0 -> 1\n", " 0 -> 3\n", "\n", " lag 1 | 10 edges (10 directed)\n", " 0 -> 0\n", " 0 -> 1\n", " 1 -> 0\n", " 0 -> 3\n", " 1 -> 1\n", " 2 -> 1\n", " 1 -> 3\n", " 2 -> 2\n", " 2 -> 3\n", " 3 -> 3\n" ] }, { "data": { "text/html": [ "0123" ], "text/plain": [ "'0123'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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", ")\n", "X, truth = series.data, series.graph\n", "print(f\"{X.shape[0]} time steps of {X.shape[1]} variables\")\n", "print(truth)\n", "display(andrey.viz.draw(truth.summary_graph()))" ] }, { "cell_type": "markdown", "id": "time_series-6", "metadata": {}, "source": [ "## Links from the past: PCMCI and Granger\n", "\n", "PCMCI tests each lagged link with partial correlations, given the strongest earlier causes it\n", "finds first; `tau_max` is the largest lag tested. Granger regresses each variable on the past\n", "`maxlag` steps of every variable with a cross-validated lasso. Neither tests links within one step." ] }, { "cell_type": "code", "execution_count": 4, "id": "time_series-7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "pcmci temporal | 4 nodes | 1 lags\n", "\n", " lag 1 | 10 edges (10 directed)\n", " 0 -> 0\n", " 0 -> 1\n", " 1 -> 0\n", " 0 -> 3\n", " 1 -> 1\n", " 2 -> 1\n", " 1 -> 3\n", " 2 -> 2\n", " 2 -> 3\n", " 3 -> 3\n", "\n", " p_values: a list of 1, too long to print; read metadata['p_values']\n", " statistics: a list of 1, too long to print; read metadata['statistics']\n", " parents: [[[1, 1], [0, 1]], [[1, 1], [0, 1], [2, 1]], [[2, 1], [3, 1]], [[1, 1], [2, 1], [0, 1]]]\n", "\n", "granger temporal | 4 nodes | 1 lags\n", "\n", " lag 1 | 13 edges (13 directed)\n", " 0 -> 0\n", " 0 -> 1\n", " 1 -> 0\n", " 2 -> 0\n", " 0 -> 3\n", " 1 -> 1\n", " 2 -> 1\n", " 1 -> 3\n", " 3 -> 1\n", " 2 -> 2\n", " 2 -> 3\n", " 3 -> 2\n", " 3 -> 3\n" ] } ], "source": [ "print(andrey.pcmci(X, tau_max=1), end=\"\\n\\n\")\n", "print(andrey.granger(X, maxlag=1))" ] }, { "cell_type": "markdown", "id": "time_series-8", "metadata": {}, "source": [ "## Links within a step: PCMCI+ and LPCMCI\n", "\n", "PCMCI+ adds the same-step links and orients them where the data allow; lag 0 is a CPDAG. LPCMCI\n", "allows hidden series, so it returns a partial ancestral graph (PAG): a circle marks an end the data\n", "leave open." ] }, { "cell_type": "code", "execution_count": 5, "id": "time_series-9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "pcmciplus temporal | 4 nodes | 2 lags\n", "\n", " lag 0 | 3 edges (3 directed)\n", " 0 -> 1\n", " 3 -> 0\n", " 1 -> 3\n", "\n", " lag 1 | 7 edges (7 directed)\n", " 0 -> 0\n", " 0 -> 1\n", " 1 -> 0\n", " 1 -> 1\n", " 2 -> 1\n", " 1 -> 3\n", " 2 -> 2\n", "\n", " p_values: a list of 2, too long to print; read metadata['p_values']\n", " statistics: a list of 2, too long to print; read metadata['statistics']\n", " parents: [[[1, 1], [0, 1]], [[1, 1], [0, 1], [2, 1]], [[2, 1]], [[1, 1], [2, 1]]]\n", " conflicts: []\n", " ambiguous_triples: [[0, 1, 1, 3]]\n", "\n", "lpcmci temporal | 4 nodes | 2 lags\n", "\n", " lag 0 | 2 edges (1 partially_directed, 1 directed)\n", " 0 o-> 1\n", " 0 -> 3\n", "\n", " lag 1 | 0 edges (no edges)\n", "\n", " links: a list of 2, too long to print; read metadata['links']\n", " p_values: a list of 2, too long to print; read metadata['p_values']\n", " statistics: a list of 2, too long to print; read metadata['statistics']\n", " parents: [[], [], [], [[0, 0]]]\n", " conflicts: []\n" ] } ], "source": [ "print(andrey.pcmciplus(X, tau_max=1), end=\"\\n\\n\")\n", "print(andrey.lpcmci(X, tau_max=1))" ] }, { "cell_type": "markdown", "id": "time_series-10", "metadata": {}, "source": [ "## Weights: VARMA-LiNGAM\n", "\n", "VARMA-LiNGAM uses the non-Gaussian noise to order the variables within a step, and estimates the\n", "weight of every link. `prune=True` sets small weights to zero." ] }, { "cell_type": "code", "execution_count": 6, "id": "time_series-11", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "varmalingam temporal | 4 nodes | 2 lags\n", "\n", " lag 0 | 2 edges (2 directed)\n", " 0 -> 1\n", " 0 -> 3\n", "\n", " lag 1 | 7 edges (7 directed)\n", " 0 -> 0\n", " 1 -> 0\n", " 1 -> 1\n", " 2 -> 1\n", " 1 -> 3\n", " 2 -> 3\n", " 3 -> 3\n", "\n", " order: 0 -> 3 -> 2 -> 1\n", "\n", "lag 0 weights, row = cause:\n", "[[ 0. 1.19 0. -1.38]\n", " [ 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. ]]\n" ] } ], "source": [ "out = andrey.varma_lingam(X, order=(1, 1), prune=True)\n", "print(out, end=\"\\n\\n\")\n", "print(\"lag 0 weights, row = cause:\")\n", "print(np.round(out.structure.lag_weights[0], 2))" ] }, { "cell_type": "markdown", "id": "time_series-12", "metadata": {}, "source": [ "## Compare the links\n", "\n", "Each method's links at lags 0 and 1 beside the true ones, a variable's dependence on its own past\n", "included. A same-step link counts once, whatever its direction; a lagged link counts with its\n", "direction. LPCMCI is left out: its lagged graphs keep only the links it orients as causes, so its\n", "PAG above is the answer to read." ] }, { "cell_type": "code", "execution_count": 7, "id": "time_series-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", "
true links foundextra links
PCMCI10 of 120
Granger10 of 123
PCMCI+9 of 121
VARMA-LiNGAM9 of 120
\n", "
" ], "text/plain": [ " true links found extra links\n", "PCMCI 10 of 12 0\n", "Granger 10 of 12 3\n", "PCMCI+ 9 of 12 1\n", "VARMA-LiNGAM 9 of 12 0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def links(structure):\n", " found = set()\n", " for lag in (0, 1):\n", " if lag not in structure.lags:\n", " continue\n", " for i, j, _ in structure.lag(lag).oriented_edges():\n", " found.add((min(i, j), max(i, j), 0) if lag == 0 else (i, j, 1))\n", " return found\n", "\n", "\n", "true_links = links(truth)\n", "fits = {\n", " \"PCMCI\": andrey.pcmci(X, tau_max=1),\n", " \"Granger\": andrey.granger(X, maxlag=1),\n", " \"PCMCI+\": andrey.pcmciplus(X, tau_max=1),\n", " \"VARMA-LiNGAM\": andrey.varma_lingam(X, order=(1, 1), prune=True),\n", "}\n", "rows = {}\n", "for name, fit in fits.items():\n", " found = links(fit.structure)\n", " rows[name] = {\n", " \"true links found\": f\"{len(found & true_links)} of {len(true_links)}\",\n", " \"extra links\": len(found - true_links),\n", " }\n", "pd.DataFrame(rows).T" ] }, { "cell_type": "markdown", "id": "time_series-14", "metadata": {}, "source": [ "## Panels: longitudinal LiNGAM\n", "\n", "Longitudinal LiNGAM takes one table per time point, the same units in the same rows. Here 500 units\n", "are measured at 3 time points; within each, `x0 -> x1 -> x2`, and each variable depends on its own\n", "value at the time point before. The columns come in one order drawn from the seed." ] }, { "cell_type": "code", "execution_count": 8, "id": "time_series-15", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "longitudinallingam temporal | 3 nodes | 2 lags\n", "\n", " lag 0 | 2 edges (2 directed)\n", " x0 -> x1\n", " x1 -> x2\n", "\n", " lag 1 | 9 edges (9 directed)\n", " x0 -> x0\n", " x0 -> x1\n", " x1 -> x0\n", " x0 -> x2\n", " x2 -> x0\n", " x1 -> x1\n", " x1 -> x2\n", " x2 -> x1\n", " x2 -> x2\n", "\n", " causal_orders: [[], [0, 1, 2], [0, 1, 2]]\n", "\n", "lag 0 at the last time point: [('x0', 'x1', 'directed'), ('x1', 'x2', 'directed')]\n" ] } ], "source": [ "CHAIN = [(\"x0\", \"x1\"), (\"x1\", \"x2\")] # the true edges within a time point, cause first\n", "\n", "\n", "def sample_panel(seed: int, n: int = 500, shuffle: bool = True) -> list[pd.DataFrame]:\n", " \"\"\"Draw n units at 3 time points, the columns in one order drawn from the seed.\"\"\"\n", " rng = np.random.default_rng(seed)\n", " b0 = np.array([[0.0, 0.0, 0.0], [0.8, 0.0, 0.0], [0.0, -0.7, 0.0]]) # row = effect\n", " mix = np.linalg.inv(np.eye(3) - b0)\n", " panel = [rng.uniform(-1, 1, (n, 3)) @ mix.T]\n", " for _ in range(2):\n", " panel.append((0.5 * panel[-1] + rng.uniform(-1, 1, (n, 3))) @ mix.T)\n", " names = [\"x0\", \"x1\", \"x2\"]\n", " columns = list(rng.permutation(names)) if shuffle else names\n", " return [pd.DataFrame(x, columns=names)[columns] for x in panel]\n", "\n", "\n", "out = andrey.longitudinal_lingam(sample_panel(seed=0), n_lags=1)\n", "print(out, end=\"\\n\\n\")\n", "print(\"lag 0 at the last time point:\", out.structure.lag(0).oriented_edges())" ] }, { "cell_type": "markdown", "id": "time_series-16", "metadata": {}, "source": [ "The [release cards](https://andrey.abel.ai/docs/benchmarks.html) compare each method with tigramite,\n", "lingam, or causal-learn on the simulations of the methods' papers.\n", "\n", "## Use your own data\n", "\n", "Pass a table with one row per time step, such as a DataFrame whose column names become the graph's\n", "labels, or change `seed` and `n` above and rerun." ] } ], "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 }