Examples

Nonlinear additive-noise models (CAM and RESIT)#

Open in Colab Download this notebook

CAM, the causal additive model, learns a DAG from data in which every variable is a nonlinear function of its parents plus independent noise, x_j = f_j(parents) + e_j. When the functions are nonlinear, the direction of every edge is identifiable, so CAM returns a DAG instead of the equivalence class that PC and GES return. It fits a smooth function of each parent, adds edges greedily by the likelihood, and drops the parents whose smooth terms are not significant. It also returns a causal order.

RESIT, regression with subsequent independence tests, fits the same kind of model with noise of any distribution, for a handful of variables.

Start#

Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.13.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.13.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
from andrey.metrics import shd

Data#

sample_scm draws six variables on an Erdos-Renyi graph, each the sum of a nonlinear function of its parents and Gaussian noise, with 1,000 rows and standardized columns. The columns come in an order drawn from the seed, so their order says nothing about the causal order.

def draw(seed: int, d: int = 6, noise: str = "gaussian"):
    """Draw 1,000 rows of d variables with additive noise; return the data and the true DAG."""
    sample = sample_scm(
        graph="erdos_renyi",
        functional="additive_noise",
        noise=noise,
        d=d,
        n=1000,
        seed=seed,
        scale="standardize",
        density=2.0,
    )
    return sample.data, sample.graph


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

1000 observations of 6 variables
012345

Fit CAM#

CAM returns a DAG and ordering, a causal order as column indices, causes first: the order of its search before pruning.

cam = andrey.cam(X)
print(f"causal order {cam.ordering}")
for source, target, _ in cam.structure.oriented_edges():
    print(f"  {source} -> {target}")
display(andrey.viz.draw(cam.structure))

causal order (5, 2, 4, 3, 0, 1)
  0 -> 1
  3 -> 0
  3 -> 1
  4 -> 1
  5 -> 2
  4 -> 3
012345

Compare with the truth#

The structural Hamming distance (SHD) counts the edges to add, remove, or reverse to reach the true graph; the empty graph’s SHD is the number of true edges.

rows = {
    "CAM": {
        "edges": len(cam.structure.oriented_edges()),
        "SHD": shd(cam.structure, truth),
    },
    "empty graph": {"edges": 0, "SHD": len(truth.oriented_edges())},
}
pd.DataFrame(rows).T
edges SHD
CAM 6 3
empty graph 0 6

Check over ten datasets#

The same model, drawn with ten seeds. The table gives CAM’s SHD, with the empty graph’s beside it.

rows = {}
for seed in range(10):
    data, graph = draw(seed)
    rows[seed] = {
        "CAM": shd(andrey.cam(data).structure, graph),
        "empty graph": len(graph.oriented_edges()),
    }
table = pd.DataFrame(rows).T
table.index.name = "seed"
display(table)
print("median SHD:", table.median().to_dict())
CAM empty graph
seed
0 3 6
1 0 6
2 0 6
3 0 6
4 1 6
5 0 6
6 2 6
7 3 6
8 0 6
9 6 6 
median SHD: {'CAM': 0.5, 'empty graph': 6.0}

Rescaling changes nothing#

CAM does not depend on a column’s scale or location, so multiplying or shifting a column leaves its graph unchanged. Here every column is rescaled and shifted differently.

scale = np.array([0.01, 1.0, 5.0, 50.0, 300.0, 1e4])
unit_free = X * scale + 7 * scale
same = np.array_equal(andrey.cam(unit_free).structure.to_numpy(), cam.structure.to_numpy())
print(f"same graph on the rescaled data: {same}")

same graph on the rescaled data: True

RESIT#

RESIT fits the same model with noise of any distribution. It finds a causal order from the sink up: the next sink is the variable whose additive-model residual an independence test finds most independent of the remaining variables. Then it removes each parent whose removal leaves the residual independent. Each independence test reads all the remaining variables at once, so RESIT suits a handful of variables: its accuracy claim holds up to 10, and it takes seconds at 5 and minutes at 30. Here it runs on five variables with Gumbel noise.

X5, truth5 = draw(seed=0, d=5, noise="gumbel")
resit = andrey.resit(X5)
print(f"causal order {resit.ordering}")
for source, target, _ in resit.structure.oriented_edges():
    print(f"  {source} -> {target}")
print(f"SHD {shd(resit.structure, truth5)}, empty graph {len(truth5.oriented_edges())}")

causal order (1, 2, 0, 3, 4)
  1 -> 0
  1 -> 2
  2 -> 3
  2 -> 4
  3 -> 4
SHD 0, empty graph 5

Check RESIT over ten datasets#

The same model with ten seeds: RESIT’s SHD, with the empty graph’s beside it.

rows = {}
for seed in range(10):
    data, graph = draw(seed, d=5, noise="gumbel")
    rows[seed] = {
        "RESIT": shd(andrey.resit(data).structure, graph),
        "empty graph": len(graph.oriented_edges()),
    }
table = pd.DataFrame(rows).T
table.index.name = "seed"
display(table)
print("median SHD:", table.median().to_dict())
RESIT empty graph
seed
0 0 5
1 0 5
2 0 5
3 0 5
4 0 5
5 0 5
6 2 5
7 0 5
8 0 5
9 0 5 
median SHD: {'RESIT': 0.0, 'empty graph': 5.0}

When to use it#

CAM needs effects that are nonlinear and additive in the parents, with additive noise and no hidden common causes. With linear effects and Gaussian noise the data identify only the equivalence class and the direction of the other edges is arbitrary, so use PC or GES; noise whose scale depends on the parents can reverse edges. The 0.13.0 release card compares CAM with CAM 1.0 from 10 to 100 variables. Change seed, d, or the noise ("gumbel") and rerun to see how it responds.

RESIT needs the same kind of data with noise of any distribution, and no more than about ten variables: from 20 variables its errors come near the empty graph’s on Gaussian noise and a fit takes about a minute, so use CAM for more variables. The release card compares RESIT with the authors’ code and with lingam from 5 to 20 variables.