Time series (PCMCI, PCMCI+, LPCMCI, Granger, VARMA-LiNGAM, longitudinal LiNGAM)

Download this notebook

Six methods find causes in data measured over time. PCMCI and Granger find the links from earlier steps; PCMCI+ also finds links within one step, and LPCMCI allows series that are not measured. VARMA-LiNGAM estimates the weights of the links within one step and from the step before. Longitudinal LiNGAM takes a panel instead: the same units measured at a few time points. Each returns a TemporalStructure, one graph per lag.

Start

Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.9.0 from PyPI; local kernels use the installed environment.

import os
import sys

IN_COLAB = "google.colab" in sys.modules

if IN_COLAB:
    os.environ["ANDREY_DEVICE"] = "cpu"
    !pip install --quiet "andrey-core[viz]==0.9.0" pandas

import numpy as np
import pandas as pd
from IPython.display import display

import andrey
import andrey.viz
from andrey.data import sample_scm

Data

sample_scm draws one series of 1,000 time steps over 4 variables: a vector autoregression with one lag, links within each step from a random DAG, and uniform noise. The drawing merges the lags into one graph; a loop is a variable’s dependence on its own past.

series = sample_scm(
    graph="erdos_renyi",
    functional="var",
    noise="uniform",
    d=4,
    n=1000,
    seed=1,
    density=1.0,
)
X, truth = series.data, series.graph
print(f"{X.shape[0]} time steps of {X.shape[1]} variables")
print(truth)
display(andrey.viz.draw(truth.summary_graph()))

1000 time steps of 4 variables
temporal  |  4 nodes  |  2 lags

  lag 0  |  2 edges (2 directed)
    0 -> 1
    0 -> 3

  lag 1  |  10 edges (10 directed)
    0 -> 0
    0 -> 1
    1 -> 0
    0 -> 3
    1 -> 1
    2 -> 1
    1 -> 3
    2 -> 2
    2 -> 3
    3 -> 3

Links from the past: PCMCI and Granger

PCMCI tests each lagged link with partial correlations, given the strongest earlier causes it finds first; tau_max is the largest lag tested. Granger regresses each variable on the past maxlag steps of every variable with a cross-validated lasso. Neither tests links within one step.

print(andrey.pcmci(X, tau_max=1), end="\n\n")
print(andrey.granger(X, maxlag=1))

pcmci  temporal  |  4 nodes  |  1 lags

  lag 1  |  10 edges (10 directed)
    0 -> 0
    0 -> 1
    1 -> 0
    0 -> 3
    1 -> 1
    2 -> 1
    1 -> 3
    2 -> 2
    2 -> 3
    3 -> 3

  p_values: a list of 1, too long to print; read metadata['p_values']
  statistics: a list of 1, too long to print; read metadata['statistics']
  parents: [[[1, 1], [0, 1]], [[1, 1], [0, 1], [2, 1]], [[2, 1], [3, 1]], [[1, 1], [2, 1], [0, 1]]]

granger  temporal  |  4 nodes  |  1 lags

  lag 1  |  13 edges (13 directed)
    0 -> 0
    0 -> 1
    1 -> 0
    2 -> 0
    0 -> 3
    1 -> 1
    2 -> 1
    1 -> 3
    3 -> 1
    2 -> 2
    2 -> 3
    3 -> 2
    3 -> 3

Links within a step: PCMCI+ and LPCMCI

PCMCI+ adds the same-step links and orients them where the data allow; lag 0 is a CPDAG. LPCMCI allows hidden series, so it returns a partial ancestral graph (PAG): a circle marks an end the data leave open.

print(andrey.pcmciplus(X, tau_max=1), end="\n\n")
print(andrey.lpcmci(X, tau_max=1))

pcmciplus  temporal  |  4 nodes  |  2 lags

  lag 0  |  3 edges (3 directed)
    0 -> 1
    3 -> 0
    1 -> 3

  lag 1  |  7 edges (7 directed)
    0 -> 0
    0 -> 1
    1 -> 0
    1 -> 1
    2 -> 1
    1 -> 3
    2 -> 2

  p_values: a list of 2, too long to print; read metadata['p_values']
  statistics: a list of 2, too long to print; read metadata['statistics']
  parents: [[[1, 1], [0, 1]], [[1, 1], [0, 1], [2, 1]], [[2, 1]], [[1, 1], [2, 1]]]
  conflicts: []
  ambiguous_triples: [[0, 1, 1, 3]]

lpcmci  temporal  |  4 nodes  |  2 lags

  lag 0  |  2 edges (1 partially_directed, 1 directed)
    0 o-> 1
    0 -> 3

  lag 1  |  0 edges (no edges)

  links: a list of 2, too long to print; read metadata['links']
  p_values: a list of 2, too long to print; read metadata['p_values']
  statistics: a list of 2, too long to print; read metadata['statistics']
  parents: [[], [], [], [[0, 0]]]
  conflicts: []

Weights: VARMA-LiNGAM

VARMA-LiNGAM uses the non-Gaussian noise to order the variables within a step, and estimates the weight of every link. prune=True sets small weights to zero.

out = andrey.varma_lingam(X, order=(1, 1), prune=True)
print(out, end="\n\n")
print("lag 0 weights, row = cause:")
print(np.round(out.structure.lag_weights[0], 2))

varmalingam  temporal  |  4 nodes  |  2 lags

  lag 0  |  2 edges (2 directed)
    0 -> 1
    0 -> 3

  lag 1  |  7 edges (7 directed)
    0 -> 0
    1 -> 0
    1 -> 1
    2 -> 1
    1 -> 3
    2 -> 3
    3 -> 3

  order: 0 -> 3 -> 2 -> 1

lag 0 weights, row = cause:
[[ 0.    1.19  0.   -1.38]
 [ 0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.  ]]

Compare the links

Each method’s links at lags 0 and 1 beside the true ones, a variable’s dependence on its own past included. A same-step link counts once, whatever its direction; a lagged link counts with its direction. LPCMCI is left out: its lagged graphs keep only the links it orients as causes, so its PAG above is the answer to read.

def links(structure):
    found = set()
    for lag in (0, 1):
        if lag not in structure.lags:
            continue
        for i, j, _ in structure.lag(lag).oriented_edges():
            found.add((min(i, j), max(i, j), 0) if lag == 0 else (i, j, 1))
    return found


true_links = links(truth)
fits = {
    "PCMCI": andrey.pcmci(X, tau_max=1),
    "Granger": andrey.granger(X, maxlag=1),
    "PCMCI+": andrey.pcmciplus(X, tau_max=1),
    "VARMA-LiNGAM": andrey.varma_lingam(X, order=(1, 1), prune=True),
}
rows = {}
for name, fit in fits.items():
    found = links(fit.structure)
    rows[name] = {
        "true links found": f"{len(found & true_links)} of {len(true_links)}",
        "extra links": len(found - true_links),
    }
pd.DataFrame(rows).T

	true links found	extra links
PCMCI	10 of 12	0
Granger	10 of 12	3
PCMCI+	9 of 12	1
VARMA-LiNGAM	9 of 12	0

Panels: longitudinal LiNGAM

Longitudinal LiNGAM takes one table per time point, the same units in the same rows. Here 500 units are measured at 3 time points; within each, x0 -> x1 -> x2, and each variable depends on its own value at the time point before. The columns come in one order drawn from the seed.

CHAIN = [("x0", "x1"), ("x1", "x2")]  # the true edges within a time point, cause first


def sample_panel(seed: int, n: int = 500, shuffle: bool = True) -> list[pd.DataFrame]:
    """Draw n units at 3 time points, the columns in one order drawn from the seed."""
    rng = np.random.default_rng(seed)
    b0 = np.array([[0.0, 0.0, 0.0], [0.8, 0.0, 0.0], [0.0, -0.7, 0.0]])  # row = effect
    mix = np.linalg.inv(np.eye(3) - b0)
    panel = [rng.uniform(-1, 1, (n, 3)) @ mix.T]
    for _ in range(2):
        panel.append((0.5 * panel[-1] + rng.uniform(-1, 1, (n, 3))) @ mix.T)
    names = ["x0", "x1", "x2"]
    columns = list(rng.permutation(names)) if shuffle else names
    return [pd.DataFrame(x, columns=names)[columns] for x in panel]


out = andrey.longitudinal_lingam(sample_panel(seed=0), n_lags=1)
print(out, end="\n\n")
print("lag 0 at the last time point:", out.structure.lag(0).oriented_edges())

longitudinallingam  temporal  |  3 nodes  |  2 lags

  lag 0  |  2 edges (2 directed)
    x0 -> x1
    x1 -> x2

  lag 1  |  9 edges (9 directed)
    x0 -> x0
    x0 -> x1
    x1 -> x0
    x0 -> x2
    x2 -> x0
    x1 -> x1
    x1 -> x2
    x2 -> x1
    x2 -> x2

  causal_orders: [[], [0, 1, 2], [0, 1, 2]]

lag 0 at the last time point: [('x0', 'x1', 'directed'), ('x1', 'x2', 'directed')]

The release cards compare each method with tigramite, lingam, or causal-learn on the simulations of the methods’ papers.

Use your own data

Pass a table with one row per time step, such as a DataFrame whose column names become the graph’s labels, or change seed and n above and rerun.