andrey.nogam#

andrey.nogam(data, *, alpha=0.001, knowledge=None)[source]#

Learn a DAG with NoGAM, a causal order from the score and residuals (additive noise).

NoGAM [1] models each variable as a nonlinear function of its parents plus noise of any distribution, x_j = f_j(parents) + e_j. Under that model the entry j of the score of the joint density, d log p(x) / d x_j, is a function of the variable’s noise alone exactly when j has no children, and that noise is what a regression of x_j on all the other variables leaves [1]. NoGAM estimates the score at the rows by a kernel estimator built on Stein’s identity [2], predicts each entry from its variable’s regression residual, takes the variable predicted with the smallest error as the next leaf, repeats on the others, and prunes the complete DAG of the order by the smooth-term tests of CAM [3]. 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) – Significance level of the pruning test, in (0, 1): a parent stays when its smooth term’s p-value is at most alpha. The causal order does not read it. 0.001 is the value of the paper’s experiments.
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 leaf 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 pruning fit and never in the result; the score estimate and the residual regressions still read every remaining variable. The pruning never removes a required edge.

Returns:

Wraps a dag GraphStructure; metadata["algorithm"] is "nogam". ordering is the causal order of the leaf search, as a tuple of column indices, causes first. metadata["tie_break"] is the tie-break rule the fit ran under, andrey.config.tie_break.

Return type:

StructureOutput

Raises:

TypeError – If an argument has the wrong type, such as a string for a number; the message names the method, the argument, what it accepts, and the value. Also if knowledge is neither None nor a Knowledge.
ValueError – If an argument 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 5 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; or if more than half of the pairs of rows are equal on the remaining variables at a step, which leaves the score estimate’s kernel no bandwidth. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When a column holds at most 10 distinct whole numbers and so may hold category codes, which the fit reads as numbers. Each message warns once per process.

Notes

When to use. NoGAM suits i.i.d. continuous observational data whose effects are nonlinear, with additive noise of any distribution and no unmeasured common causes [1]; under that model the DAG is identified [4], and it returns the DAG with a causal order. It does not assume Gaussian noise, as score_matching() and cam() do, and its order search needs no independence test, as resit()’s does. The paper shows it on 10 and 20 variables and 1,000 rows, with Gaussian, Laplace, Gumbel, Beta, exponential, and gamma noise. With linear effects and Gaussian noise only the Markov equivalence class is identified: use pc() or ges().

The steps. The paper’s two phases [1]:

Causal order. On the remaining variables, three estimates are made at every row. The score grad log p(x), by the Stein gradient estimator [2] with the Gaussian kernel exp(-||x - y||^2 / (2 s^2)), the bandwidth s the median distance between the rows. Each variable’s residual, from its kernel ridge regression on all the other remaining variables. And each score entry’s prediction from its own variable’s residual, by another kernel ridge regression. The variable whose score entry has the smallest mean squared prediction error is the next leaf; it is removed, and the estimates are made again on the others. The causal order is the leaves in reverse.
Pruning. Each variable is fit on all the variables before it in the order, an additive model of penalized thin plate regression splines with ten basis functions per parent and smoothing chosen by GCV [5], 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. A parent stays when its smooth term’s p-value is at most alpha, the pruning of cam() [3].

The kernel ridge regressions take the settings of the paper’s experiments: the kernel exp(-gamma ||x - y||^2) with gamma = 0.1, the ridge term 0.01, no intercept, and predictions by K-fold cross-validation, each row predicted by the model fit on the other folds. gamma is the value the paper states; the authors’ implementation sets it to the ridge term’s value, 0.01.

The choices the paper leaves open. It names no number of folds: there are 5, each a block of consecutive rows, so the order of the rows can change the result. It states the cross-validation for the score entries; the residuals are predicted across the same folds. It states no ridge term for the score estimate: it is 0.001, the value in the authors’ implementation, added to the kernel matrix as the estimator’s definition writes it, a matrix with ones on its diagonal; the median distance is taken over the pairs of distinct rows. The errors are compared as estimated, not divided by the entries’ variances.

The score estimate and the kernel ridge regressions read the columns’ units, which the model does not. The columns are therefore centered and divided by their standard deviations before the order search, so rescaling or shifting a column changes no step; the pruning reads each variable through an affine map to a fixed interval. The paper’s experiments run on the columns as simulated. The fit draws no random numbers and runs on one BLAS thread.

Ties. The errors are compared rounded to 9 decimal places, as hc() compares its score changes per row, so two errors that differ only in their rounding tie. A tie goes to the variable with the lower rank computed from the data, as for hc(); 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. The numbers below were measured on the paper’s simulation with standardized columns and 1,000 rows, on one core, 20 graphs per setting: Erdos-Renyi graphs with as many edges as variables or four times as many, each link one draw of a Gaussian process, with Gaussian, Laplace, or Gumbel noise, at 10 and 20 variables. With as many edges as variables the order put 90% to 92% of the true edges in the true direction at 10 variables and 81% to 86% at 20, and the DAG was 1.2 to 1.7 and 4.2 to 6.0 edges from the true one (SHD), where the empty graph is 10 and 20 away; with four times as many edges the order put 98% to 99% and 92% to 95% in the true direction, and the SHD was 19.6 to 21.1 and 48.8 to 49.4 against the empty graph’s 40 and 80. On the graphs with as many edges as variables cam() was at least as accurate on the same data, under every noise: its order put 90% to 98% of the true edges in the true direction, and its SHD was 0.7 to 0.9 at 10 variables and 1.4 to 3.3 at 20. The paper’s tables are on data as simulated, and the order search reads standardized columns, where it orders fewer edges of a sparse graph correctly: in a separate check on ten graphs of 10 variables with as many edges, the same search put 99% of the true edges in the true direction under Gaussian noise and 98% under Laplace noise when it read the columns as simulated, and 87% and 91% on standardized columns. The search compares the variables’ prediction errors in absolute terms on standardized columns, and a sink’s error grows with its variance-to-noise ratio (its variance over the variance of its noise term), since the same error in its residual costs more where the noise is a small part of the variable. Where links have a strong linear part that ratio is largest at the sinks: on additive-noise data with links z + sin z (10 and 20 variables, as many edges as variables; Gaussian, Laplace, and Gumbel noise) the order put 47% to 62% of the true edges in the true direction at 1,000 rows, and 4,000 rows raised it only from 56% to 62% under Gaussian noise at 10 variables. Use cam() on such data (76% to 94% on the same datasets). A fit took 7.5 to 9.7 s at 10 variables and 28 to 40 s at 20, within 0.16 GB; it was not measured in graphs of 50 variables, where a single fit took 194 s. Each step fits ten kernel ridge regressions per remaining variable, each on four fifths of the rows, so the order search’s time grows with the square of n_variables and the cube of n_samples, and its memory with the square of n_samples. The pruning starts from the complete DAG of the order, so its last fits have nearly n_variables smooth terms each, and it needs at least twice as many rows as columns. Binary columns, and any column with fewer than 3 distinct values, are not supported: the pruning fits a smooth term of each parent, and such a column raises ValueError before the fit. 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

Examples

>>> import numpy as np
>>> from andrey import nogam
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-2, 2, 300)
>>> y = np.sin(2 * x) + 0.3 * rng.uniform(-1, 1, 300)
>>> out = nogam(np.column_stack([x, y]))
>>> out.structure.kind, out.ordering
('dag', (0, 1))