andrey.cam

andrey.cam(data, *, alpha=0.001, max_parents=None, knowledge=None)

Learn a DAG with CAM, the causal additive model (nonlinear additive effects).

CAM [1] models each variable as a sum of smooth functions of its parents plus Gaussian noise, x_j = sum_k f_jk(x_k) + e_j, which identifies the DAG when the functions are nonlinear [6] [7]. It preselects each variable’s possible parents by boosting, adds edges greedily by the drop in the log residual variance until the DAG is complete within them, and prunes the edges whose smooth terms are not significant. 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 pruning test: a parent stays when its smooth term’s p-value is at most alpha, in (0, 1). 0.001 is the paper’s value at every sample size. A larger value keeps more of the true edges a small sample tests weakly, and more false edges with them (see Notes).

• max_parents (int | None) – A cap on every node’s parents, at least 0: the order search adds no edge that gives a node more, so no node of the returned DAG has more than max_parents parents. A cap can stop the search short of the best-scoring DAG. None sets no cap, as in the paper.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. Forbidden edges, edges into an earlier tier, and pairs outside the allowed set are never possible parents; the boosting that preselects the possible parents reads the data alone. The search starts from the required edges and the pruning never removes one. The causal order puts the cause of each required edge, and every earlier tier, first.

Returns:

Wraps a dag GraphStructure; metadata["algorithm"] is "cam" and metadata["tie_break"] the rule ties followed. ordering is a causal order of the search’s DAG before pruning, 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, 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; if its required edges give a node more than max_parents parents; or if a column’s required parents need more coefficients than data has rows, naming the column.

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.

See also

andrey.anm

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

andrey.direct_lingam

Linear effects with non-Gaussian noise.

andrey.pc, andrey.ges, which

Notes

When to use. CAM suits i.i.d. continuous observational data whose effects are nonlinear and additive in the parents, with additive Gaussian noise and no unmeasured common causes [1]; under that model the DAG is identified [6] [7], and it returns the DAG with a causal order. Its preselection of each variable’s possible parents is designed for many variables with few parents each [1]. With linear effects and Gaussian noise only the Markov equivalence class is identified, and the orientation of the edges it leaves undirected is arbitrary: use pc() or ges(). Noise whose scale depends on the parents (location-scale data) can reverse edges, because the model’s noise variance is fixed.

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

1. Neighbor preselection. Each variable is regressed on all the others by componentwise L2 boosting [2], 100 iterations of step size 0.1, with one P-spline base learner per other variable [3]: a cubic B-spline basis on 20 equidistant interior knots over the variable’s range, a second-order difference penalty, and the penalty weight that gives the learner 4 degrees of freedom. The possible parents are the ten variables picked most often, among those picked at least three times. When more than ten qualify, the ones picked more often than the eleventh stay, so a tie at the tenth place keeps fewer than ten rather than choosing by column order. When more than ten tie at the top count, none is picked more often than the eleventh, and the variable has no possible parents.

2. Order search. From the empty graph, the search adds the edge i -> j whose fit most lowers the log residual variance of j, among the possible parents that close no cycle, whatever the sign of the change, and refits only j’s remaining candidates, until none is left. Each fit is an additive model of penalized thin plate regression splines with ten basis functions per parent and smoothing chosen by GCV [4]; with fewer than 30 rows per parent each term gets ceil(n_samples / (3 n_parents)) functions, at least 3, and a variable with fewer distinct values than that gets one per distinct value. The paper states no cap on a variable’s parents, so max_parents sets none by default; a candidate whose model would have as many coefficients as rows is never added.

3. Pruning. Each variable is fit on its parents in that DAG, and a parent stays when its smooth term’s p-value [5] is at most alpha, the paper’s rule.

The rescaling of a column changes no step: the base learners and the smooth terms read a variable through an affine map to a fixed interval, and each gain is a ratio of residual variances. The fit draws no random numbers and runs on one BLAS thread.

Ties. Gains are compared rounded to 9 decimal places, as hc() compares its score changes per row, so two gains that differ only in their rounding tie. A tie goes to the edge whose cause, then effect, has the lower rank computed from the data, as for hc(), and so do the variables the causal order could place next; 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, 10 to 100 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 14 of 16 settings and fewer in 2, returned the same graph on 278 of 320 datasets, and fit 6.5 to 14.7 times faster on one core. The reference implementation caps a variable’s parents at the smaller of the variables minus one and the rows over 20, 50 at 1,000 rows, which binds at 100 variables, where max_parents sets none by default. The numbers below were measured on simulated data with 1,000 rows from Erdos-Renyi graphs of mean degree 2 with standardized columns, ten graphs at 20 and at 100 variables, on one core. On additive-noise data with Gaussian noise it oriented 88.5% to 95% of the true edges correctly at 20 variables, depending on the generator, and 97.5% at 100, with median errors (SHD) of 2 to 2.5 and 8 against the empty graph’s 20 and 100; with Gumbel noise at 20 variables it reversed 22.5% of them. Location-scale data, whose noise scale changes with the parents, break the model’s fixed noise variance: it reversed 42% of the true edges at 20 variables and 32% at 100, with median errors of 11.5 and 55.5. Linear Gaussian data identify only the Markov equivalence class, and there the two directions of an edge fit equally, their gains differing by less than 1e-6, so the direction of an edge the class leaves undirected is arbitrary: it reversed 36% of the true edges at 20 variables and 32% at 100, with median errors of 11 and 68. The model has no hidden common causes: with 2 hidden variables among 22, over ten graphs of the Gaussian additive-noise generator, its median errors were 12 against the empty graph’s 16, mostly extra edges, and it joined 42 of the 49 pairs of observed variables that share only a hidden cause. Effects that interact are outside the model too: with x3 the product of x1 and x2 plus noise, on 1,000 rows and standard normal causes, whose product has no additive part, it found at most one of the two edges into x3 and returned an edge out of x3 in 6 of 10 fits; with causes uniform on (0, 2) it made x3 the cause of both, and added an edge between them, in 8 of 10. On the additive-noise data a fit took a median of 3 s at 20 variables and 21 s at 100, about 40% of it in the neighbor preselection and 55% in the order search. The smoothing of each fit is a local minimum of its GCV score, which can have several, so two edges of nearly equal gain can be added in either order. Binary columns, and any column with fewer than 3 distinct values, are not supported: a smooth term needs at least 3, so such a column raises ValueError. Categorical columns are not modeled: text labels raise ValueError, and integer codes are read as numbers, a smooth function of the code, with a DataTypeWarning. The numbers are a reference, not a hard limit: they change with the data generator, the density, the method’s settings, and the machine.

Large graphs and many rows. On the rescaled additive-noise data with 1,000 rows, over five graphs at 200 variables, three at 500, and two at 1,000 and at 2,000, its median errors stayed below the empty graph’s at every size but grew faster than the variables, from 8% of the empty graph’s at 100 variables to 37% at 2,000. Nearly all of them were extra edges, and on one graph each at 100 and 500 variables all but one lay between variables unrelated or more than three edges apart in the true graph: the preselection picks, among more variables, those that fit a variable’s residuals best by chance, and the pruning tests them on the same rows, so of the searched edges between such pairs it kept 1.1% and 5.4%. The preselection fits d - 1 base learners for each of the d variables, so its cost grows as d^2 and takes over the fit, 85% of it at 1,000 variables and 88% at 2,000; a fit’s cost grows about as d^1.6 up to 1,000 variables and d^1.9 from 1,000 to 2,000. A fit passes about a minute near 200 variables, 10 minutes near 800, and an hour near 2,000, where it took a median of 57 minutes and peaked at 0.64 GB of memory. With 20 variables, over three graphs each, a fit took a median of 5 s at 5,000 rows and 17.5 s at 20,000 and peaked at 140 MB and 214 MB, about 100 MB of it the imports; the median errors stayed at 2. Above 2,000 distinct values a column’s thin plate basis is built on a fixed subsample of 2,000 of them: at 5,000 and 20,000 rows that returned the same graph as knots on every value in 6 of 6 fits, in up to a fifth less time.

Small samples and skewed causes. With few rows per variable it keeps false edges: on independent columns it kept a median of 26 at 40 rows and 60 variables and 7.5 at 200 rows and 60 variables, against 0 at 500 rows and 20. The search keeps each variable’s best-fitting candidates and the pruning tests them on the same rows, and the smooth-term tests run slightly liberal at small samples: an unrelated parent among five reached a p-value of at most 0.01 in 3.3% of fits at 100 rows. A larger alpha trades false edges for true ones: at 100 rows, alpha=0.01 found 2 more of 100 true adjacencies on additive-noise graphs of 10 variables and 5 more false ones. A strongly skewed cause can be reversed: with a lognormal cause of log-scale 2 and an effect of log(x) plus noise, it reversed the edge in 6 of 10 fits. A variance-stabilizing transform such as the log of a positive column suits data whose effect stays nonlinear on the new scale: with an effect of sin(2 log(x) / 3) at log-scale 3, it found the edge in 6 of 10 fits on x and 10 of 10 on log(x). An effect that becomes linear, as log(x) itself does, leaves a linear Gaussian pair, whose direction no fit identifies.

References

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

Buhlmann, P., Peters, J., and Ernest, J. (2014). CAM: Causal additive models, high-dimensional order search and penalized regression. Annals of Statistics 42(6), 2526-2556.

[2]

Buhlmann, P., and Hothorn, T. (2007). Boosting algorithms: Regularization, prediction and model fitting. Statistical Science 22(4), 477-505.

[3]

Eilers, P. H. C., and Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science 11(2), 89-121.

[4]

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

[5]

Wood, S. N. (2013). On p-values for smooth components of an extended generalized additive model. Biometrika 100(1), 221-228.

[6] (1,2)

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.

[7] (1,2)

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.

Examples

>>> import warnings
>>> import numpy as np
>>> from andrey import cam
>>> 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 = cam(np.column_stack([x, y]))
>>> out.structure.kind, out.ordering
('dag', (0, 1))