Hidden common causes (RCD)#

Download this notebook

When two variables share a cause that is not in the data, most methods either join them with a wrong edge or miss the dependence. RCD, repetitive causal discovery, assumes linear mechanisms with non-Gaussian noise. It regresses each variable on candidate ancestors and tests the residuals for independence, which finds the direct causes and the pairs that share a hidden cause, marked <->.

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 RCD#

andrey.rcd returns a PAG over the observed variables: -> for a direct cause and <-> for a pair with a hidden common cause. metadata["confounded_pairs"] lists the <-> pairs by column index.

out = andrey.rcd(X)
pairs = [[X.columns[i] for i in pair] for pair in out.metadata["confounded_pairs"]]
print("confounded pairs:", out.metadata["confounded_pairs"], "=", pairs)
for source, target, kind in out.structure.oriented_edges():
    print(f"{source} {'<->' if kind == 'bidirected' else '->'} {target}")
display(andrey.viz.draw(out.structure))

confounded pairs: [[0, 1]] = [['x0', 'x1']]
x0 <-> x1
x0 -> x3
x1 -> x2

Edge weights#

weighted_adjacency[i, j] is the coefficient of variable i in the equation for variable j. A pair marked <-> has no weight: its dependence comes from the hidden cause.

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 RCD found the hidden pair x0 <-> x1 and both direct edges, and how many edges it reported.

def named_edges(structure) -> set:
    """The edges by name; a bidirected pair in alphabetical order, whatever the column order."""
    return {
        (*sorted((a, b)), kind) if kind == "bidirected" else (a, b, kind)
        for a, b, kind in structure.oriented_edges()
    }


expected = named_edges(truth)
rows = {}
for seed in range(10):
    found = named_edges(andrey.rcd(sample(seed)).structure)
    rows[seed] = {
        "hidden pair": ("x0", "x1", "bidirected") in found,
        "direct edges": {("x0", "x3", "directed"), ("x1", "x2", "directed")} <= found,
        "edges": len(found),
        "same as the truth": found == expected,
    }
pd.DataFrame(rows).T.rename_axis("seed")

	hidden pair	direct edges	edges	same as the truth
seed
0	True	True	3	True
1	True	True	3	True
2	True	True	3	True
3	True	True	3	True
4	True	True	3	True
5	True	True	3	True
6	True	True	3	True
7	True	True	3	True
8	True	True	3	True
9	True	True	3	True

When to use it#

RCD relies on linear mechanisms and non-Gaussian noise, and it needs many rows: a weak hidden cause leaves only a small trace in the residuals it tests. Without hidden causes, DirectLiNGAM returns a full DAG; with Gaussian noise, FCI returns the PAG the conditional independences support. The 0.2.0 release card measures RCD on larger graphs. Change the seed, the hidden effect (1.5), or n and rerun to see how RCD responds.