Hidden common causes (BottomUpParceLiNGAM)#

Download this notebook

BottomUpParceLiNGAM assumes linear mechanisms with non-Gaussian noise and finds a causal order from the bottom up: at each step it looks for a sink, a variable whose regression residual on the others is independent of them. A hidden common cause leaves the variables it confounds without a sink, so the method stops there and leaves them unordered, joined by o-o, rather than guess an order.

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

Data#

A linear model with uniform noise over five variables, one of them hidden: h drives x0 and x1, x0 drives x3, and x1 drives x2. The data keep only x0 to x3, so x0 and x1 depend on each other with no edge between them. The true graph over the observed variables marks that pair <->. The columns come in an order drawn from the seed, so their order says nothing about the causal order.

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

import andrey
import andrey.viz
from andrey.core import ARROW, TAIL

COLUMNS = ["x0", "x1", "x2", "x3"]
EDGES = [("h", "x0"), ("h", "x1"), ("x0", "x3"), ("x1", "x2")]  # h has no column


def sample(seed: int, n: int = 2000, shuffle: bool = True) -> pd.DataFrame:
    """Draw n rows of x0..x3, the columns in an order drawn from the seed; h is left out."""
    rng = np.random.default_rng(seed)
    h, e0, e1, e2, e3 = rng.uniform(-1, 1, (5, n))
    x0, x1 = 1.5 * h + e0, 1.5 * h + e1
    frame = pd.DataFrame({"x0": x0, "x1": x1, "x2": 0.6 * x1 + e2, "x3": 0.6 * x0 + e3})
    return frame[rng.permutation(frame.columns)] if shuffle else frame


# marks[i, j] is the mark at x_i of the edge between x_i and x_j.
marks = np.zeros((4, 4), dtype=np.int8)
marks[0, 1] = marks[1, 0] = ARROW  # x0 <-> x1
marks[0, 3], marks[3, 0] = TAIL, ARROW  # x0 -> x3
marks[1, 2], marks[2, 1] = TAIL, ARROW  # x1 -> x2
truth = andrey.GraphStructure.from_numpy(marks, kind="pag", labels=tuple(COLUMNS))

X = sample(seed=0)
print(f"{X.shape[0]} observations, {X.shape[1]} variables, columns {list(X.columns)}")
display(andrey.viz.draw(truth))

2000 observations, 4 variables, columns ['x0', 'x1', 'x3', 'x2']

Fit BottomUpParceLiNGAM#

andrey.bottom_up_parce_lingam returns a PAG: -> for a direct cause and o-o between variables it leaves unordered. metadata["causal_order"] gives the order, the unordered variables nested in one list. alpha is the significance level of the sink test, Bonferroni-corrected across the remaining variables; at the default 0.1 a true sink is sometimes rejected by chance on this small model, so this notebook uses 0.01.

out = andrey.bottom_up_parce_lingam(X, alpha=0.01)
print("causal order:", out.metadata["causal_order"], "over columns", list(X.columns))
for source, target, kind in out.structure.oriented_edges():
    print(f"{source} {'->' if kind == 'directed' else 'o-o'} {target}")
display(andrey.viz.draw(out.structure))

causal order: [[0, 1], 2, 3] over columns ['x0', 'x1', 'x3', 'x2']
x0 o-o x1
x0 -> x3
x1 -> x2

Edge weights#

weighted_adjacency[i, j] is the coefficient of variable i in the equation for variable j, from an adaptive-Lasso regression of each variable on those before it in the order.

weights = pd.DataFrame(out.weighted_adjacency, index=X.columns, columns=X.columns)
weights.loc[COLUMNS, COLUMNS].round(2)

	x2	x3
x0	0.00	0.6
x1	0.61	0.0
x2	0.00	0.0
x3	0.00	0.0

Check over ten datasets#

The same model, drawn with ten seeds. For each, the table says whether the method left exactly x0 and x1 unordered, whether it found both direct edges, and its causal order.

def by_name(order: list, columns) -> list:
    """The causal order with each column index replaced by its name."""
    return [
        [columns[i] for i in item] if isinstance(item, list) else columns[item] for item in order
    ]


rows = {}
for seed in range(10):
    X_seed = sample(seed)
    fit = andrey.bottom_up_parce_lingam(X_seed, alpha=0.01)
    order = by_name(fit.metadata["causal_order"], X_seed.columns)
    blocks = [set(item) for item in order if isinstance(item, list)]
    found = set(fit.structure.oriented_edges())
    rows[seed] = {
        "x0, x1 unordered": blocks == [{"x0", "x1"}],
        "direct edges": {("x0", "x3", "directed"), ("x1", "x2", "directed")} <= found,
        "causal order": str(order),
    }
pd.DataFrame(rows).T.rename_axis("seed")

	x0, x1 unordered	direct edges	causal order
seed
0	True	True	[['x0', 'x1'], 'x3', 'x2']
1	True	True	[['x0', 'x1'], 'x2', 'x3']
2	True	True	[['x1', 'x0'], 'x2', 'x3']
3	True	True	[['x1', 'x0'], 'x3', 'x2']
4	True	True	[['x0', 'x1'], 'x2', 'x3']
5	True	True	[['x1', 'x0'], 'x2', 'x3']
6	True	True	[['x1', 'x0'], 'x3', 'x2']
7	True	True	[['x0', 'x1'], 'x2', 'x3']
8	True	True	[['x0', 'x1'], 'x2', 'x3']
9	True	True	[['x1', 'x0'], 'x2', 'x3']

When to use it#

BottomUpParceLiNGAM relies on linear mechanisms and non-Gaussian noise. Its unordered block holds the confounded variables and every variable above them in the order, so on a larger graph it can join pairs that share no cause, and edges inside the block stay unoriented. RCD marks only the pairs it attributes to a hidden cause; without hidden causes, DirectLiNGAM returns a full DAG. The 0.2.0 release card measures it on larger graphs. Change the seed, alpha, or n and rerun to see how it responds.