andrey.resit

andrey.resit(data, *, alpha=0.05, knowledge=None)

Learn a DAG with RESIT, regression with subsequent independence tests (additive noise).

RESIT [1] models each variable as a function of its parents plus noise independent of its non-descendants, x_j = f_j(parents) + e_j, which identifies the DAG when the functions are nonlinear or the noise is not Gaussian [1] [3]. It finds a causal order from the sink up, taking as the next sink the variable whose additive-model residual is most independent of the others by an HSIC test [2], then removes each parent whose removal leaves the residual independent of the variable’s predecessors. It returns a DAG and a causal order.

Parameters:

• data (ArrayLike) – Observed data, one row per sample. Coerced to float64. A DataFrame’s column names become the structure’s labels.

• alpha (float) – The significance level of the edge-removal tests: a parent goes when the residual without it is independent of the variable’s predecessors at a p-value above alpha, in (0, 1). The causal order does not read it. 0.05 is the paper’s value.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. A variable is the next sink only when no remaining variable is in a later tier or is the effect of a required edge out of it. Forbidden edges, edges into an earlier tier, and pairs outside the allowed set are in no regression and never in the result; each test still reads every remaining variable or predecessor. The edge removal never removes a required edge.

Returns:

Wraps a dag GraphStructure; metadata["algorithm"] is "resit" and metadata["tie_break"] the rule ties followed. ordering is the causal order of the sink search, as a tuple of column indices, causes first.

Return type:

StructureOutput

Raises:

• TypeError – If knowledge is neither None nor a Knowledge.

• ValueError – If an argument has the wrong type, lies outside its range, or is not one of its choices; the message names the method, the argument, what it accepts, and the value. Also if data is not a usable 2-D matrix: not numeric, holding NaN or inf, with fewer than 6 rows or fewer than twice as many rows as columns, or with a constant column, a column of fewer than 3 distinct values, or a column that repeats another up to scale and shift, each named in the message. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

• FloatingPointError – If an HSIC test gives a non-finite value, rather than let it steer a choice.

See also

andrey.cam

Additive effects with Gaussian noise, scored by the likelihood.

andrey.anm

Decides the direction between two variables under an additive noise model.

andrey.pc, andrey.ges, which

Notes

When to use. RESIT suits i.i.d. continuous observational data whose effects are nonlinear and additive in the parents, with additive noise of any distribution and no unmeasured common causes [1]; it reads the noise through independence tests, not a likelihood, so the noise need not be Gaussian, but noise whose scale depends on the parents breaks the independence of the residuals. Each test of the causal order reads all the remaining variables jointly, so it is designed for a handful of variables [1]. It returns a DAG and a causal order. With linear effects and Gaussian noise only the Markov equivalence class is identified: use pc() or ges(); for more variables with Gaussian noise, use cam(), and with hidden variables, camuv().

The steps. The paper’s settings, in its two phases [1]:

1. Causal order. Each remaining variable is regressed on the other remaining ones and its residual tested against them jointly; the variable whose residual has the largest p-value is the next sink, and the variables still remaining are its possible parents.

2. Edge removal. Each variable’s possible parents are tried one at a time, in the rank order of Ties: a parent goes when the residual of the fit on the other current parents has a p-value above alpha against all the variable’s predecessors in the order. With one parent left, the residual without it is the variable itself.

Each regression is an additive model of penalized thin plate regression splines with ten basis functions per term and smoothing chosen by GCV [4], with ceil(n_samples / (3 n_parents)) functions per term, at least 3, when there are fewer than 30 rows per term, and at most one per distinct value. Each test is the HSIC test [2] with Gaussian kernels on the residual and on the other side’s joint squared distances, each bandwidth the median of the positive squared distances over every pair of rows, and the null approximated by a gamma law matched to its mean and variance.

The kernels read the columns centered and divided by their standard deviations, and each smooth term reads a variable through an affine map to a fixed interval, so rescaling or shifting a column changes no step; the centering also keeps a large offset from cancelling the distances’ digits. The fit draws no random numbers and runs on one BLAS thread.

Ties. p-values are compared as their logarithms rounded to 9 decimal places, as hc() compares its score changes per row, so two p-values that differ only in their rounding tie. The logarithms are computed directly, so p-values below the smallest double, about 1e-308, which strongly dependent residuals reach with a few thousand rows, keep their order. A tie goes to the variable with the lower rank computed from the data, as for hc(), and the edge removal tries the parents in that rank order; the fit runs on the columns in that rank order, so the result does not depend on the column order unless two columns are identical. The column-order guide says what the rank keeps under a new column order.

Missing values. A NaN in data raises ValueError: every fit reads every row of its variables, so a row with a missing value has no place in it. Drop the incomplete rows first, which suits values missing completely at random, or use pc() or fci(), whose tests delete rows per test.

Limitations. On additive-noise data with Gaussian or Gumbel noise, 5 to 20 variables, and 1,000 rows, standardized and not (the 0.13.0 release card), it made no more errors than the paper’s reference implementation in 12 of 12 settings, one of them after a rerun on fresh datasets, returned the same graph on 87 of 220 datasets, and fit 3.2 to 4.6 times faster on one core; against an implementation that regresses with random forests it made fewer errors in 11 of 12 settings and no more in the other, and fit 2.1 to 2.9 times faster. The accuracy claim holds up to 10 variables: at 20 variables, on Gaussian additive-noise data, no implementation clearly beats the empty graph (mean errors of 19.4 for this fit and 20.2 for the reference implementation against 20.0), so 20 variables is a range for speed. The edge removal tries a variable’s parents one at a time and the paper’s algorithm fixes no order, so the order is arbitrary and can move a graph by about one edge: over 40 datasets of 10 variables, 48 orders moved the mean errors (SHD) from 0.8 below to 0.8 above the column order’s, and returned its graph on 9 to 17 of the 40. Binary columns, and any column with fewer than 3 distinct values, are not supported: every regression fits a smooth term of each parent, and such a column raises ValueError before the fit. Measured on simulated additive-noise data with 1,000 rows from Erdos-Renyi graphs of mean degree 2 with standardized columns, on one core. It oriented 88% of the true edges correctly at 5 variables and 80% at 10, with median errors (SHD) of 1 and 5 against the empty graph’s 5 and 10. From 20 variables its median errors on these data come near the empty graph’s or exceed them, 18 against 20 at 20 variables, 31 against 30 at 30, and 91 against 50 at 50: each test of the causal order reads all the remaining variables at once and loses power, and an early wrong sink carries through the rest of the order. At 20 variables, over ten graphs per generator, it oriented 65.5% to 78% of the true edges correctly on four additive-noise generators, the most with Gumbel noise, with median errors of 18 to 19.5 against the empty graph’s 20; it oriented 46% on linear Gaussian data and 28% on location-scale data, whose noise scale changes with the parents. A fit took a median of 1.5 s at 5 variables, 45 s at 20, 3.3 minutes at 30, and 21 minutes at 50. The causal order fits about d^2 / 2 additive models of up to d - 1 terms; their cost grows with the terms up to about 35, past which each term’s basis shrinks to keep three rows per basis function, so a fit’s cost grows about as d^3.3 up to 50 variables and as d^2 beyond. A fit passes about a minute near 20 variables, 10 minutes near 40, and an hour near 75; at 100 variables, a fit stopped after an hour had run 3,565 of the 5,049 tests of its causal order, and the measured cost per test puts the whole fit near two hours. Its accuracy ends sooner: over five graphs per size, its median errors are at most half the empty graph’s up to 10 variables and 73% or more of it from 12 variables. Each HSIC test builds matrices with one row and one column per sample, so time and memory grow with the square of n_samples: at 5 variables a fit took 1.7 s and peaked at 0.17 GB with 1,000 rows, 25 s and 0.75 GB with 4,000, and 107 s and 2.6 GB with 8,000; on Sachs’s 7,466 rows and 11 variables it took 6.7 minutes and 2.3 GB, and kept all 55 possible edges, as every removal test rejected. The edge removal tests each parent against all the variable’s predecessors at once, and many of those p-values lie near alpha, so the order the parents are tried in changes the result: trying them in the data’s column order instead of by rank gave the same causal orders at 20 variables but graphs a median of 10 to 19.5 edges apart, at about the same accuracy. A cause with a long tail can reverse its edge: with a lognormal cause of log-sd 3 and an effect of its logarithm plus noise, on 500 rows, the true order came out in 0 of 5 fits, against 5 of 5 at log-sd 1 and 2; replacing the cause by its ranks before the fit gave 5 of 5 at every log-sd, while its logarithm, which makes that pair linear and Gaussian, did not. The smoothing of each fit is a local minimum of its GCV score, which can have several, so two variables of nearly equal p-value can be taken as the sink in either order. The numbers are a reference, not a hard limit: they change with the data generator, the density, the method’s settings, and the machine.

References

[1] (1,2,3,4,5)

Peters, J., Mooij, J. M., Janzing, D., and Scholkopf, B. (2014). Causal discovery with continuous additive noise models. Journal of Machine Learning Research 15, 2009-2053.

[2] (1,2)

Gretton, A., Fukumizu, K., Teo, C. H., Song, L., Scholkopf, B., and Smola, A. J. (2007). A kernel statistical test of independence. NeurIPS 20, 585-592.

[3]

Hoyer, P. O., Janzing, D., Mooij, J. M., Peters, J., and Scholkopf, B. (2009). Nonlinear causal discovery with additive noise models. NeurIPS 21, 689-696.

[4]

Wood, S. N. (2003). Thin plate regression splines. Journal of the Royal Statistical Society B 65(1), 95-114.

Examples

>>> import warnings
>>> import numpy as np
>>> from andrey import resit
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-2, 2, 300)
>>> y = np.sin(2 * x) + 0.3 * rng.standard_normal(300)
>>> with warnings.catch_warnings():
...     warnings.simplefilter("ignore")
...     out = resit(np.column_stack([x, y]))
>>> out.structure.kind, out.ordering
('dag', (0, 1))