Pairwise direction#

Cover: an hourglass whose sand falls one way only, drawn in navy lines

Does x cause y, or does y cause x? With only two variables, independence tests cannot tell: both directions fit the data equally well. ANM and PNL instead use the shape of the relationship to infer direction. This release adds ANM. As of 0.4.0, we support both methods.

We explain how two variables can reveal their direction, run both tests on two pairs, and report benchmark results.

Direction from two variables#

PC and GES can orient an edge only through a third variable, such as a common effect. A pair has no third variable, so x -> y and y -> x explain the data equally well. These methods leave the edge undirected.

A model of the mechanism can break the tie. Suppose y = f(x) + e, where the noise e is independent of x. Regress y on x, and the residual looks independent of x. Regress x on y, and in general no residual is independent of y, unless f is linear and the noise is Gaussian.

Each test fits both directions and checks whether the residual is independent of the proposed cause. It uses HSIC, a kernel test that detects nonlinear dependence as well as correlation. The direction whose residual looks independent is the causal one.

The methods fit different models:

  • ANM (additive noise model) fits y = f(x) + e by Gaussian-process regression. The direction is generally identifiable when f is nonlinear or the noise is non-Gaussian.

  • PNL (post-nonlinear model) fits y = f2(f1(x) + e), where f2 is invertible. This is an additive noise model observed through a distortion, such as a sensor that saturates. ANM is the case where f2 changes nothing. PNL fits each transform as a sum of smooth basis functions.

Supported methods#

Before each release, we benchmark each supported method against another package and check its results on known data. ANM is new in this release. PNL was experimental until now.

Method

Model

Returns

ANM

y = f(x) + e

A two-node DAG, both p-values, and both test statistics

PNL

y = f2(f1(x) + e)

A two-node DAG, both p-values, and both test statistics

Choose a method for your data:

  • ANM when the effect is a smooth function of the cause plus independent noise. Its regression time grows with the cube of the number of rows, so it suits up to a few thousand rows.

  • PNL when the effect may pass through a distortion or you have more rows. In our benchmarks, its median fit took under 2 s at 5,000 rows.

Both assume that no hidden common cause links the pair. When one is plausible, use the methods in the 0.2.0 post.

Reading the result#

Each test returns a two-node DAG and four values in metadata:

  • pval_forward and pval_backward: the p-values of the independence tests for x -> y and y -> x. A large p-value means the residual looks independent of the proposed cause, so that direction fits.

  • stat_forward and stat_backward: the two HSIC statistics, each divided by its mean under independence. Smaller values mean weaker dependence.

By default, the edge points in the direction with the larger p-value. On a large sample, both p-values can round to zero. The smaller statistic then decides. This rule always selects a direction, even for unrelated variables.

With a significance level such as alpha=0.05, a direction fits when its p-value exceeds the level. The test draws an edge only when exactly one direction fits.

Example pairs#

In both pairs, x causes y, and e is Laplace noise:

  • Additive: y = x + 0.3 * x**3 + 0.5 * e. The noise adds to a function of x, so both models hold.

  • Saturating: a sensor reads y = tanh(x + 0.5 * e). The noise passes through the saturation, so only the post-nonlinear model holds.

example.py#
import numpy as np
import pandas as pd

import andrey

# x causes y in both pairs; e is noise.
rng = np.random.default_rng(0)
x = rng.laplace(size=500)
e = rng.laplace(size=500)

# The additive pair: y = x + 0.3 * x**3 + 0.5 * e.
additive = pd.DataFrame({"x": x, "y": x + 0.3 * x**3 + 0.5 * e})
# The saturating pair, which a sensor reads: y = tanh(x + 0.5 * e).
saturating = pd.DataFrame({"x": x, "y": np.tanh(x + 0.5 * e)})

print("Additive pair\n")
print(andrey.anm(additive), andrey.pnl(additive), sep="\n\n")
print("\nSaturating pair\n")
print(andrey.anm(saturating), andrey.pnl(saturating), sep="\n\n")

# For comparison: ANM on the saturating pair with a significance level.
print()
print(andrey.anm(saturating, alpha=0.05))

Two pairs, two models

The additive pair500 rows of y = x + 0.3 * x**3 + 0.5 * e
-200-1000100-8-6-4-20246xy

ANMnames x -> y, the true direction

x -> y0.029y -> x00.05

PNLnames x -> y, the true direction

x -> y0.072y -> x2.4e-090.05
The saturating pair500 rows of y = tanh(x + 0.5 * e)
-101-8-6-4-20246xy

ANMnames y -> x, reversed

x -> y0y -> x1.7e-070.05

PNLnames x -> y, the true direction

x -> y0.586y -> x2.1e-060.05
Each bar is one direction's p-value, on a scale from 0 to 1; the cobalt bar is the direction the test names, and the dashed line marks 0.05. On the additive pair, both tests name x -> y and reject y -> x. On the saturating pair, ANM rejects both directions, since the saturation is not additive noise, and names the one with the larger p-value; PNL keeps only x -> y.

On the additive pair, ANM and PNL both select the true direction, x -> y, and reject y -> x.

On the saturating pair, the distortion breaks ANM, and PNL models it:

  • ANM rejects both directions. Both p-values are near zero because the saturation is not additive noise. Its default rule still selects a direction, here the reversed one, y -> x. With alpha=0.05, ANM draws no edge (the last result below).

  • PNL models the saturation as its outer transform. Its p-value for x -> y is large. Its p-value for y -> x is near zero, so it selects the true direction.

The script prints each result in full:

example.py#
Additive pair

ANM  dag  |  2 nodes  |  1 edges (1 directed)

  x -> y

  pval_forward: 0.02928
  pval_backward: 0
  stat_forward: 1.725
  stat_backward: 11.33

PNL  dag  |  2 nodes  |  1 edges (1 directed)

  x -> y

  pval_forward: 0.07243
  pval_backward: 2.379e-09
  stat_forward: 1.525
  stat_backward: 4.287

Saturating pair

ANM  dag  |  2 nodes  |  1 edges (1 directed)

  y -> x

  pval_forward: 0
  pval_backward: 1.695e-07
  stat_forward: 21.66
  stat_backward: 4.322

PNL  dag  |  2 nodes  |  1 edges (1 directed)

  x -> y

  pval_forward: 0.5863
  pval_backward: 2.091e-06
  stat_forward: 0.8921
  stat_backward: 3.344

ANM  dag  |  2 nodes  |  0 edges (no edges)

  pval_forward: 0
  pval_backward: 1.695e-07
  stat_forward: 21.66
  stat_backward: 4.322

Benchmark results#

We compared both methods with causal-learn, a Python package that implements them, on simulated cause-effect pairs. We tested ANM on additive-noise pairs and PNL on post-nonlinear pairs.

Method

Compared with

Andrey’s speedup

Direction

ANM

causal-learn

about 2x

the same as causal-learn’s on every pair

PNL

causal-learn

up to 1,470x

true on more pairs than causal-learn’s

The 0.4.0 release card reports every measurement, including the sizes, machine, and settings.

Speedup over causal-learnRelease card
Method
200 rows2.08x
500 rows1.88x
1,000 rows2.05x
2,000 rows2.38x
4,000 rows2.47x
8,000 rows2.39x
200 rows1,470x
500 rows1,310x
1,000 rows993x
2,000 rows696x
5,000 rows275x

Each bar is causal-learn's time divided by Andrey's on the same pair, the median over the pairs of that size, on 4 cores per fit. The bars are on a log scale, so equal steps are equal factors; each method's bars share its own scale.

Our ANM implementation fits the same Gaussian process and uses the same independence test as causal-learn, so both return the same graph on every pair.

The PNL implementations differ more in how they fit and test each direction. Ours selects the true direction on 94% of pairs; causal-learn’s selects it on 45%.

Pairs with the true direction

ANMThe same direction as causal-learn on every pair

Andrey and causal-learn

00.250.50.7512005001,0002,0004,0008,000rowsAndrey and causal-learn, 200 rows: 0.75Andrey and causal-learn, 500 rows: 0.8Andrey and causal-learn, 1,000 rows: 0.95Andrey and causal-learn, 2,000 rows: 0.95Andrey and causal-learn, 4,000 rows: 1Andrey and causal-learn, 8,000 rows: 1
PNLAndrey: every pair from 1,000 rows

Andreycausal-learn

00.250.50.7512005001,0002,0005,000rowscausal-learn, 200 rows: 0.5causal-learn, 500 rows: 0.5causal-learn, 1,000 rows: 0.5causal-learn, 2,000 rows: 0.5causal-learn, 5,000 rows: 0.25Andrey, 200 rows: 0.75Andrey, 500 rows: 0.95Andrey, 1,000 rows: 1Andrey, 2,000 rows: 1Andrey, 5,000 rows: 1
The share of pairs, per size, on which each package names the true direction; higher is better. Each size has 20 pairs, except 10 at 4,000 rows and 5 at 8,000 rows for ANM. The empty graph names no direction, so its share is 0 at every size.

At 200 rows, both methods select the reversed direction on a quarter of the pairs. A direction test selects a direction on every pair, right or wrong. Check both p-values before you trust it.

Other highlights#

  • Statistics break ties. On a large sample, both p-values can round to zero. PNL now selects the direction with the smaller HSIC statistic, relative to its mean under independence, instead of drawing no edge. ANM uses the same rule. Both statistics are in metadata.

  • ANM in any units. ANM standardizes each column first, so its answer does not depend on either variable’s units.

For details, read the changelog.

Try it#

pip install -U andrey-core

The pairwise notebook runs both tests on simulated pairs. You can open it in Colab, then pass your own two-column table to andrey.anm or andrey.pnl.

Help us build it#

Andrey is in alpha. The API may change between releases. When you try it on your own data, tell us what breaks or which method you need. We also welcome pull requests. The contributing guide explains how you can set up.