andrey.pc

andrey.pc(data, *, alpha=0.05, indep_test='fisherz', collider_order=None, categorical=None, max_categories=10, max_depth=None, knowledge=None)

Learn a CPDAG from observational data with the PC algorithm.

PC [1] [2] uses conditional-independence tests to remove edges from a complete graph. These tests check whether variables are independent given other variables. It then orients colliders, where arrows meet at a variable, and applies Meek’s rules [4]. The result is a CPDAG, representing an equivalence class of graphs whose unresolved edges stay undirected.

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 conditional-independence test. Smaller values give sparser graphs.

• indep_test (str) – Conditional-independence test. "fisherz" and "student_t" test partial correlation, for continuous data: "fisherz" applies Fisher’s z-transform and the normal distribution; "student_t" applies the Student-t distribution, which is exact for Gaussian data at any sample size. "chisq" (Pearson’s chi-square) and "gsq" (G-square, the likelihood ratio) test categorical data, coded by whole numbers. "degen_gauss" (degenerate Gaussian) and "cond_gauss" (conditional Gaussian) test mixed data, continuous and categorical columns, by a likelihood ratio. Four more tests detect nonlinear dependence in continuous data, at a far higher cost per test: "kci" is the kernel test; "rcot" its random-feature approximation; "cmiknn" a nearest-neighbor conditional mutual information test; and "gpdc" Gaussian-process regression with a distance correlation test. See Notes.

• collider_order (str | None) – The order in which colliders are oriented, which decides between two colliders that would orient one edge in opposite directions. "p_value" orients first the collider whose endpoints were separated with the largest p-value, so the colliders do not depend on the column order. "column" orients them in column order and keeps the first of two that conflict. None follows andrey.config.tie_break: "p_value" under "rank", "column" under "position". See Notes.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data test: column labels (strings) or 0-based positions. None infers them: a DataFrame’s category, bool, and text columns, and every column of at most max_categories distinct whole numbers, named in a DataTypeWarning. [] marks none. The other tests take no column types and reject it. See Data types in the getting-started guide.

• max_categories (int) – The most distinct whole numbers a column may hold to be inferred as categorical, at least 2. With a continuous test, columns within it get a one-time DataTypeWarning that they may be category codes.

• max_depth (int | None) – The maximum conditioning-set size of every conditional-independence test, at least 0: 0 runs only unconditional tests, and None sets no limit. A smaller cap runs fewer tests and can keep edges that a larger conditioning set would remove.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. The adjacency search starts from the allowed pairs and never removes a required pair; the tests can remove any other allowed pair. A pair the knowledge makes non-adjacent is tested for a separating set, so colliders over it are oriented as usual; when the search finds none within max_depth, no rule reads its unshielded triples as colliders or non-colliders. The result directs every edge the knowledge decides, then Meek’s rules direct the edges that follow, giving the maximally oriented PDAG (MPDAG) Knowledge describes; knowledge that orients the edges the search kept into a directed cycle raises ValueError.

Returns:

Wraps a cpdag GraphStructure; metadata["algorithm"] is "PC". Directed marks (->) are oriented edges, undirected (--) the ones the equivalence class leaves open. Under a mixed-data setting, metadata["categorical"] lists the categorical columns it read, by name or 0-based position.

Return type:

StructureOutput

Raises:

• NotImplementedError – If indep_test names no available test.

• TypeError – If categorical is not a list of column names and positions, or max_categories is not an int. Also if knowledge is neither None nor a Knowledge, or a given max_depth is not an int.

• ValueError – If data has too few samples for a test or a singular correlation matrix, if alpha is not in (0, 1), if max_categories is less than 2, with "chisq" or "gsq", if data holds a value that is not a whole number, or, with "degen_gauss" or "cond_gauss", if categorical names a column data lacks, or if a categorical column holds a number that is not whole or a missing value. Also if categorical is given with another test, or if collider_order is not "p_value" or "column". With "fisherz" or "student_t", also if a column repeats another up to scale and shift, before any test runs; one message names every group of copies and says to keep one column per group. Also if a given max_depth is negative, or if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

• PerformanceWarning – When a conditioning pass of size two or greater requires at least ten million CI tests; each message warns once per process. The warning does not limit the search.

• DataTypeWarning – When a mixed-data test infers the categorical columns, or finds none and so runs as the linear-Gaussian form, or when a continuous setting meets columns that look like category codes; each message warns once per process.

See also

andrey.fci

Allows unmeasured common causes and selection bias.

andrey.direct_lingam

Orients every edge when the noise is non-Gaussian.

Notes

When to use it. PC suits data with no unmeasured common causes. Its default Fisher-Z test assumes continuous data, where each variable depends linearly on its causes with roughly Gaussian noise; for categorical data, use "chisq" or "gsq", for mixed data "degen_gauss" or "cond_gauss", and for nonlinear effects in continuous data "kci", "rcot", "cmiknn", or "gpdc", at a far higher cost per test. It is fastest on sparse graphs, where every variable has few neighbors [5].

When not to. When unmeasured common causes are plausible, use fci(). When the noise is non-Gaussian and every edge should be oriented, use direct_lingam().

Choosing alpha. A test removes an edge when it cannot reject independence at alpha. A larger alpha therefore keeps more edges, true and false. On the Sachs data in the getting-started guide, raising alpha from 0.05 to 0.2 found more known links, but a smaller share of the links found were known. No value is right in general: it trades missed edges against false ones, so a result is best reported with how the graph changes with alpha.

How it works. Edge removal is order-independent (PC-stable) [3]. At each conditioning-set size, every test uses the adjacencies from the start of that size. Edges found independent are removed together. A removed pair’s separating set combines all conditioning sets that separated it. Two colliders conflict when they would orient one edge in opposite directions, so the order in which they are oriented can change the result [3]. A collider is skipped when it would give an edge an arrowhead at both ends, so of two conflicting colliders the one oriented first is kept. With collider_order="p_value" (the default), each collider x -> y <- z ranks by the largest p-value among the conditioning sets that separated x and z, and colliders are oriented from the highest rank to the lowest. When two conflicting colliders rank equal to nine decimals, the edge they dispute stays undirected, and Meek’s rules do not orient it. The colliders therefore do not depend on the column order. This adapts the conflict rule of Ramsey [10], which ranks a collider by separating sets that additional tests choose; Andrey ranks it by the tests that removed the edge, so orienting runs no new test. With collider_order="column", colliders are oriented by middle node, then by endpoint pair, in column order, so the result can depend on the column order. After the colliders, Meek’s rules run in rounds: each round collects every orientation they imply on the graph as it stands and applies them together, except that an edge two of them would direct both ways stays undirected, and so do orientations that would together close a directed cycle. The rounds repeat until one orients nothing, so Meek’s rules do not depend on the column order. Such conflicts arise only when no DAG on the skeleton has exactly the colliders found, as test errors can cause; otherwise the rounds give the CPDAG that firing the rules one at a time gives. Under andrey.config(tie_break="position"), the rules fire one at a time in column order instead, and of two firings that would orient one edge in opposite directions the first decides. The column-order guide lists what reordering the columns does to each method.

Limitations. The benchmark numbers below were measured with collider_order="column"; the default is not benchmarked yet. On simulated linear-Gaussian data from sparse random graphs (mean degree 2), from 20 to 400 variables with ten rows per variable, its median errors (SHD) stayed below the empty graph’s at every size, 357 against 400 at 400 variables. The median fit took 0.0076 s at 20 variables and 0.82 s at 400, on 16 cores. On the 223-variable ANDES network with 1,000 simulated rows, its median errors passed the empty graph’s (372 against 338). It is practical to roughly 1,000 variables on sparse data with ten rows per variable; runtime grows combinatorially with the largest neighborhood, so sparser graphs can stay feasible beyond that. On simulated mixed data (conditional linear Gaussian networks on sparse random graphs, the first half of the variables categorical with 2 to 5 categories), from 10 to 40 variables with 2,000 rows, its median errors stayed below the empty graph’s at every size: at 40 variables, 18.5 with "degen_gauss" and 9.5 with "cond_gauss", against 40. On simulated categorical networks with 2 or 3 categories per variable, of the same sizes and rows, they were 14 with "chisq" and 16 with "gsq" at 40 variables. The median fit at 40 variables took 0.38 s with "degen_gauss", 0.42 s with "cond_gauss", and 0.18 s with either categorical test, on 4 cores.

Categorical data. "chisq" and "gsq" take each distinct value of a column as one category, so the codes need not start at 0 or be consecutive; a DataFrame’s text, bool, and category columns are read the same way. A test compares the counts of x and y with independence within each configuration of the conditioning set that occurs in the data [2]. Each configuration adds (r_x - 1) * (r_y - 1) degrees of freedom, counting only the values of x and y present in it [6]; with none, the p-value is 1. A test needs several rows per cell of its table: with many categories or a large conditioning set, most cells are empty, and the p-value is unreliable.

Mixed data. "degen_gauss" and "cond_gauss" test data with continuous and categorical columns. Each distinct value of a categorical column is one category, and a DataFrame column of text or categories needs no coding. "degen_gauss", the degenerate Gaussian test [7], replaces a categorical column of k categories by k - 1 indicator columns and treats every column as Gaussian, so it suits effects that are linear in the continuous variables and add up across categories. "cond_gauss", the conditional Gaussian test [8], fits the continuous variables in each joint category of the categorical ones with a Gaussian of that category’s own mean and covariance, so a category can change a variance as well as a mean. Each category’s covariance is shrunk toward that of all rows by one pseudo-row [9], so a joint category with few rows takes mostly the covariance of all rows; to tell its own apart, a test needs several rows in each joint category of its categorical variables. Each compares the likelihoods of the models with and without the dependence, a chi-square statistic on the difference in their numbers of parameters. With no categorical column, both are the likelihood-ratio form of the partial-correlation test, which agrees with "fisherz" as the rows grow.

References

[1]

Spirtes, P., and Glymour, C. (1991). An algorithm for fast recovery of sparse causal graphs. Social Science Computer Review 9(1), 62-72.

[2] (1,2)

Spirtes, P., Glymour, C., and Scheines, R. (2000). Causation, Prediction, and Search, 2nd ed. MIT Press.

[3] (1,2)

Colombo, D., and Maathuis, M. H. (2014). Order-independent constraint-based causal structure learning. JMLR 15, 3921-3962.

[4]

Meek, C. (1995). Causal inference and causal explanation with background knowledge. UAI, 403-410.

[5]

Kalisch, M., and Buhlmann, P. (2007). Estimating high-dimensional directed acyclic graphs with the PC-algorithm. JMLR 8, 613-636.

[6]

Tsamardinos, I., Brown, L. E., and Aliferis, C. F. (2006). The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning 65(1), 31-78.

[7]

Andrews, B., Ramsey, J., and Cooper, G. F. (2019). Learning high-dimensional directed acyclic graphs with mixed data-types. Proceedings of Machine Learning Research 104, 4-21.

[8]

Andrews, B., Ramsey, J., and Cooper, G. F. (2018). Scoring Bayesian networks of mixed variables. International Journal of Data Science and Analytics 6(1), 3-18.

[9]

Friedman, J. H. (1989). Regularized discriminant analysis. Journal of the American Statistical Association 84(405), 165-175.

[10]

Ramsey, J. (2016). Improving accuracy and scalability of the PC algorithm by maximizing P-value. arXiv:1610.00378, Section 5.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from andrey import pc
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = rng.standard_normal(500)
>>> z = x + y + 0.5 * rng.standard_normal(500)
>>> pc(pd.DataFrame({"x": x, "y": y, "z": z})).structure.oriented_edges()
[('x', 'z', 'directed'), ('y', 'z', 'directed')]

On categorical data, with a chi-square test:

>>> a = rng.integers(0, 2, 2000)
>>> b = rng.integers(0, 2, 2000)
>>> c = np.where(rng.random(2000) < 0.8, a + b, rng.integers(0, 3, 2000))
>>> data = pd.DataFrame({"a": a, "b": b, "c": c})
>>> pc(data, indep_test="chisq").structure.oriented_edges()
[('a', 'c', 'directed'), ('b', 'c', 'directed')]

On mixed data, with a category column and the degenerate Gaussian test:

>>> g = rng.integers(0, 3, 2000)
>>> u = rng.standard_normal(2000)
>>> v = np.array([-1.0, 0.0, 1.0])[g] + u + 0.5 * rng.standard_normal(2000)
>>> data = pd.DataFrame({"g": pd.Categorical(g), "u": u, "v": v})
>>> pc(data, indep_test="degen_gauss").structure.oriented_edges()
[('g', 'v', 'directed'), ('u', 'v', 'directed')]