andrey.camuv

andrey.camuv(data, *, alpha=0.01, max_depth=2, knowledge=None)

Learn a PAG of nonlinear additive effects with hidden variables by CAM-UV.

CAM-UV [1] models each variable, observed or hidden, as a sum of nonlinear functions of its direct causes plus independent noise. Where a hidden variable lies between two observed ones, as a common cause or on a causal path, their relation cannot be identified by regression and independence tests, and CAM-UV marks the pair instead of guessing. It finds each variable’s direct causes and those pairs.

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 every HSIC independence test, in (0, 1): of the sink tests, the tests that remove a candidate parent, and the tests that mark a pair. 0.01 is the paper’s value.

• max_depth (int) – The maximum conditioning-set size of a sink test, at least 1: the most variables tested together as a variable’s new candidate parents, so the search runs over sets of up to max_depth + 1 variables. Each sink test also regresses on the variable’s candidates found so far, however many. Larger values find parents that only show together and take longer. The default 2 is the paper’s setting, sets of at most 3 variables.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. A forbidden edge, an edge into an earlier tier, and a pair outside allowed never make a variable a candidate parent, so a variable is not tried as the sink of a set that holds a forbidden parent of it, and none is in the result. A required edge’s cause is a candidate parent from the start, and no step removes it; its pair is never marked. A pair forbidden in both directions is never marked; on a marked pair, a forbidden i -> j whose reverse is not forbidden, such as an edge into an earlier tier, puts an arrowhead at i, as in fci().

Returns:

A GraphStructure of kind pag: i -> j for a direct cause and i o-o j for a pair joined through a hidden variable, whose relation the method cannot determine, with an arrowhead where knowledge sets one. metadata["algorithm"] is "camuv", metadata["latent_pairs"] lists the marked pairs as sorted [i, j] lists, metadata["cycle_edges"] lists as [i, j] each edge i -> j dropped from a directed cycle, and metadata["tie_break"] is the rule ties followed.

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, a column that repeats another up to scale and shift, or a column that is a linear function of others, such as a total, 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.

Warns:

• ExperimentalWarning – On the first call in a process: the method has no published benchmark, and its API may change without deprecation.

• 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.cam, andrey.resit

andrey.rcd

Linear non-Gaussian effects with hidden common causes.

andrey.fci

Hidden variables with any effects, from conditional independence alone.

Notes

When to use. CAM-UV suits i.i.d. continuous observational data whose effects are nonlinear and additive in the causes, with independent noise of any distribution, where some causes may be unmeasured [1]. Its kernel independence tests compare every pair of rows, and its search tests sets of up to max_depth + 1 variables, so it suits a moderate number of variables and rows. It returns a PAG: i -> j for a direct cause and i o-o j for a pair that shares a hidden common cause or of which one causes the other through a hidden variable, which the method does not tell apart. With linear effects, a path through a hidden variable can be regressed away and the pairs are identified differently: use rcd() for non-Gaussian noise and fci() for Gaussian noise. Without hidden variables, cam() and resit() also give a causal order.

The steps. Every pair of variables is first tested for dependence; then come the paper’s three steps, with its settings [1]:

1. Candidate parents. Over every set of 2, then up to max_depth + 1 variables, each member dependent on every other member is regressed on the others and its candidate parents, and its residual is tested against the others’ residuals on their own candidates, jointly. The member of largest p-value is a sink of the set when that p-value is above alpha and its residual on its candidates alone is dependent on each other member’s, every p-value below alpha; the other members then become its candidates. After a sweep that adds a candidate, the search starts over from sets of 2. A set that holds a pair already related, one a candidate of the other, is skipped: the sink test identifies a sink only when each member’s candidates lie outside the set.

2. Parents. Each variable’s candidates are tried one at a time, in the rank order of Ties: a candidate goes when the variable’s residual on its other candidates is independent of the candidate’s residual on its own, at a p-value above alpha.

3. Hidden paths. A dependent pair with no edge whose residuals, each on its parents, are dependent at a p-value below alpha is marked o-o.

The two dependence screens come from the authors’ implementation, not from the paper’s algorithm, which tests every member of a set and every pair. Under the model they add no conclusion: two independent variables share no noise, so their residuals are independent too, and neither a sink test nor step 3 would accept them. With tests at alpha they skip most sink tests: on the simulated data of Limitations with 20 observed variables, a fit ran 12 times faster with them, and its directed errors were lower on 34 of 40 datasets and higher on 1. A parent whose dependence the test misses, such as a weak parent beside a strong one, is not found.

The paper does not say what to do when the parents form a directed cycle, which its model rules out and failed tests can leave: every edge on a directed cycle is dropped, except a required one, and step 3 tests the pair; rcd() raises instead when required edges close a cycle with the parents it found. metadata["cycle_edges"] lists the dropped edges, and an INFO log names them. On the simulated data of Limitations, the drop fired in 15 of 70 fits, and 74 of the 78 edges it dropped came back as o-o; 24 of the 78 were true edges. The child’s residual on its reduced parents also enters its other pairs, which can mark them too.

Each regression is an additive model of penalized thin plate regression splines with ten basis functions per term and smoothing chosen by GCV [3], with ceil(n_samples / (3 n_terms)) 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, each bandwidth the median of the positive squared distances over every pair of rows, a side of several residuals read through their joint distances, and the null approximated by a gamma law matched to its mean and variance.

Every step reads the columns divided by their standard deviations, so rescaling a column changes no step; the paper’s joint kernel reads raw residuals and so weighs each by its units, which is the same on standardized data. The kernels read the residuals centered, which keeps an 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 sets run in that rank order and the parents are tried in it; 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. Binary columns, and any column with fewer than 3 distinct values, are not supported: every regression fits a smooth term of each regressor, and such a column raises ValueError before the fit. Integer category codes are read as numbers, with a DataTypeWarning. Measured on simulated additive-noise data from Erdos-Renyi graphs of mean degree 2 with standardized columns, ten graphs per size, on one core. The comparison with the reference implementations of CAM-UV used two settings with hidden variables: one hidden variable per ten observed ones, at 10 to 50 observed variables and 1,000 rows, and one node in five hidden, chosen among those with two or more children, at 5 to 15 nodes and 2,000 rows. Against them it was on par with both at every measured size in both settings, each 95% interval of the paired difference in directed errors (SHD) including zero, with 30 graphs at 4 observed variables. It ran 2 to 16 times faster. Where hidden variables take most edges, as in the second setting, neither it nor the reference implementations beat the empty graph: at 12 observed variables each made more errors than the empty graph, 8.3 per graph for this method against 6.3. In the first setting its errors fell below the empty graph’s as the graphs grew, to 62% of them at 50 variables, then rose back to 96% at 300, over eight graphs, and passed them at 400, 317 against 315 on one graph; at 10 they were about the same, 7.0 against 7.1. At 50 observed variables it marked half again as many pairs o-o as the reference implementations, 164 against 112, at lower precision: 18% of its marks were pairs joined through a hidden variable, against 26%. With 20 observed variables, 1,000 rows, and 2 or 4 hidden ones, 23% and 42% of its marks were such pairs. With none hidden, on ten graphs of 20 variables, it marked 40 pairs with nonlinear effects and 19 with linear Gaussian ones, each a false mark; without the drop of the edges on directed cycles it would mark 22 and 16, but would find fewer hidden pairs, 18 against 24 with 4 hidden variables. On the linear Gaussian data, which identify only the Markov equivalence class, it oriented 53% of the true edges, with median errors of 14.5 against the empty graph’s 20. Each sink test reads the candidates found before it, so the order the sets are visited in changes the result: visiting them in the data’s column order instead of by rank gave graphs a median of 6 to 16.5 pairs apart, at about the same accuracy. A cause with a long tail loses its edge: with a lognormal cause of log-sd 2 or 3 and an effect of its logarithm plus noise, on 500 rows, it marked the pair o-o in 10 of 10 fits, and replacing the cause by its ranks found the edge in 9 of them. One value of 1e6 in the middle variable of a nonlinear chain of 800 rows turned both of its edges into o-o marks in 5 of 5 fits. 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 first setting a fit took 21 s at 20 variables, and its cost grows about as d^2: a fit passes about a minute near 40 variables, 10 minutes near 130, and an hour near 300. 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.5 s and 0.21 GB with 1,000 rows and 78 s and 3.6 GB with 8,000, and at 20 variables 6.1 minutes and 1.7 GB with 4,000. The kernels kept between tests stay within 1 GiB, so above about 11,600 rows none is kept.

References

[1] (1,2,3)

Maeda, T. N., and Shimizu, S. (2021). Causal additive models with unobserved variables. UAI, PMLR 161, 97-106.

[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]

Wood, S. N. (2004). Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association 99, 673-686.

Examples

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