Examples

Linear non-Gaussian models (DirectLiNGAM / ICA-LiNGAM / MultiGroupDirectLiNGAM / GIN)#

Open in Colab Download this notebook

When the mechanisms are linear and the noise is non-Gaussian, the data identify a single DAG, not only its equivalence class. DirectLiNGAM finds a causal order one variable at a time with a pairwise independence measure. ICA-LiNGAM runs independent component analysis (FastICA) and reads the order from the unmixing matrix. Both return a fully oriented DAG, a causal order, and edge weights. MultiGroupDirectLiNGAM fits several datasets that share one causal order, and GIN finds the latent variables that observed indicators measure.

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.metrics as metrics
import andrey.viz
from andrey.data import benchmarks

Data#

lingam_sf draws a scale-free graph with linear mechanisms and uniform noise: here 8 variables and 2,000 samples.

seed = 7
dataset = benchmarks.scm("lingam_sf", d=8).sample(n=2000, seed=seed, scale="standardize")
X, truth = dataset.data, dataset.observed_graph
print(f"{X.shape[0]} observations, {X.shape[1]} variables, seed {seed}")
display(andrey.viz.draw(truth))

2000 observations, 8 variables, seed 7
01234567

Fit both methods#

DirectLiNGAM is deterministic. ICA-LiNGAM starts FastICA from a random point; random_state fixes it, so a run repeats exactly.

results = {
    "DirectLiNGAM": andrey.direct_lingam(X),
    "ICA-LiNGAM": andrey.ica_lingam(X, random_state=0),
}
for name, result in results.items():
    print(f"{name}: {result.structure.kind}, causal order {result.ordering}")
    display(andrey.viz.draw(result.structure))

DirectLiNGAM: dag, causal order (7, 5, 1, 3, 6, 0, 4, 2)
01234567
ICA-LiNGAM: dag, causal order (5, 7, 1, 3, 6, 0, 2, 4)
01234567

Edge weights#

weighted_adjacency[i, j] is the coefficient of variable i in the equation for variable j, so each nonzero entry is an edge i -> j.

pd.DataFrame(results["ICA-LiNGAM"].weighted_adjacency).round(2)
0 1 2 3 4 5 6 7
0 0.00 0.00 -0.72 0.00 0.51 0.0 0.00 0.0
1 0.00 0.00 0.00 0.94 -0.71 0.0 0.00 0.0
2 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.0
3 0.00 0.00 0.00 0.00 0.00 0.0 0.76 0.0
4 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.0
5 0.00 -0.72 0.00 0.38 0.00 0.0 0.66 0.0
6 0.79 0.00 0.47 0.00 0.00 0.0 0.00 0.0
7 0.00 0.00 -0.77 0.38 0.00 0.0 0.00 0.0

Check accuracy over ten datasets#

SHD counts errors against the true DAG; lower is better. The empty graph’s SHD is the number of true edges. The table covers ten seeds and their median.

methods = {
    "DirectLiNGAM": lambda X: andrey.direct_lingam(X),
    "ICA-LiNGAM": lambda X: andrey.ica_lingam(X, random_state=0),
}
rows = {}
for s in range(10):
    draw = benchmarks.scm("lingam_sf", d=8).sample(n=2000, seed=s, scale="standardize")
    d = draw.observed_graph.n_nodes
    empty = andrey.GraphStructure.from_numpy(np.zeros((d, d), dtype=np.uint8), kind="dag")
    rows[s] = {"Empty graph": metrics.score(empty, draw.observed_graph)["shd"]}
    for name, fit in methods.items():
        rows[s][name] = metrics.score(fit(draw.data), draw.observed_graph)["shd"]
shd = pd.DataFrame(rows)
shd.columns.name = "seed"
shd["median"] = shd.median(axis=1)
shd
seed 0 1 2 3 4 5 6 7 8 9 median
Empty graph 12 12 12 12 12 12 12 12 12 12 12.0
DirectLiNGAM 0 0 0 0 0 0 0 0 0 0 0.0
ICA-LiNGAM 0 8 0 0 3 0 0 0 0 0 0.0

When to use them#

Both methods rely on non-Gaussian noise. With Gaussian noise the order is not identified; PC, GES, BOSS, or GRaSP return the CPDAG that the data do support. ICA-LiNGAM needs many samples: at n = 10d, the size the release cards use, it scores worse than the empty graph, and this notebook uses n = 250d. Change seed, d, or n and rerun to see how each method responds.

Several datasets: MultiGroupDirectLiNGAM#

multi_group_direct_lingam takes datasets over the same variables that share one causal order and returns one result per dataset, each with its own weights. Here two groups follow the chain x0 -> x1 -> x2, with different weights on x0 -> x1, their columns in one order drawn from the seed.

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


def sample_groups(seed: int, n: int = 1000, shuffle: bool = True) -> list[pd.DataFrame]:
    """Draw n rows of the chain per weight on x0 -> x1, all in one column order from the seed."""
    rng = np.random.default_rng(seed)
    groups = []
    for weight in (0.8, -1.2):
        x0 = rng.uniform(-1, 1, n)
        x1 = weight * x0 + rng.uniform(-1, 1, n)
        x2 = 0.6 * x1 + rng.uniform(-1, 1, n)
        groups.append(pd.DataFrame({"x0": x0, "x1": x1, "x2": x2}))
    columns = rng.permutation(groups[0].columns) if shuffle else groups[0].columns
    return [frame[columns] for frame in groups]


groups = sample_groups(seed=3)
for g, out in enumerate(andrey.multi_group_direct_lingam(groups)):
    weights = pd.DataFrame(
        out.weighted_adjacency, index=groups[g].columns, columns=groups[g].columns
    )
    print(f"group {g}: {out.structure.oriented_edges()}, weight {weights.loc['x0', 'x1']:.2f}")

group 0: [('x0', 'x1', 'directed'), ('x1', 'x2', 'directed')], weight 0.74
group 1: [('x0', 'x1', 'directed'), ('x1', 'x2', 'directed')], weight -1.17

Latent variables: GIN#

gin groups observed indicators under the latent variables that cause them and orders the latents. Here a1 to a3 measure one latent and b1 to b3 another, which the first causes. The indicators come in an order drawn from the seed.

INDICATORS = [("l1", "l2")]  # the latents l1 -> l2 have no column
INDICATORS += [("l1", f"a{i}") for i in (1, 2, 3)] + [("l2", f"b{i}") for i in (1, 2, 3)]


def sample_indicators(seed: int, n: int = 1000, shuffle: bool = True) -> pd.DataFrame:
    """Draw n rows of the six indicators of l1 -> l2, in an order drawn from the seed."""
    rng = np.random.default_rng(seed)
    l1 = rng.uniform(-1, 1, n)
    l2 = 0.9 * l1 + rng.uniform(-1, 1, n)
    loadings = {"a1": (l1, 1.0), "a2": (l1, 0.8), "a3": (l1, -1.2)}
    loadings |= {"b1": (l2, 1.0), "b2": (l2, 1.1), "b3": (l2, -0.7)}
    frame = pd.DataFrame(
        {name: w * latent + 0.3 * rng.uniform(-1, 1, n) for name, (latent, w) in loadings.items()}
    )
    return frame[rng.permutation(frame.columns)] if shuffle else frame


gin_out = andrey.gin(sample_indicators(seed=5))
print(gin_out.structure.oriented_edges())
display(andrey.viz.draw(gin_out.structure))

[('L1', 'a1', 'directed'), ('L2', 'b1', 'directed'), ('L2', 'b2', 'directed'), ('L1', 'a2', 'directed'), ('L1', 'a3', 'directed'), ('L2', 'b3', 'directed'), ('L1', 'L2', 'directed')]
a1b1b2a2a3b3L1L2