andrey.ges

andrey.ges(data, *, score_func=None, lambda_value=1.0, max_parents=None, categorical=None, max_categories=10, knowledge=None)

Learn a CPDAG from observational data with greedy equivalence search (GES).

GES searches CPDAGs, which represent equivalence classes of graphs. It adds edges while they improve a score, the BIC (Bayesian information criterion) by default, then removes edges in a backward sweep [1].

Parameters:

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

• score_func (str | None) – Local score objective: "local_score_BIC", the linear-Gaussian BIC, for continuous data. "local_score_BDeu" is the BDeu score [3] and "local_score_BIC_categorical" the multinomial BIC [7], for categorical data coded by whole numbers. "local_score_degen_gauss" (degenerate Gaussian) [4] and "local_score_cond_gauss" (conditional Gaussian) [5] are BIC scores for mixed data, continuous and categorical columns. None reads the column types from categorical, or infers them, and picks "local_score_BDeu" when every column is categorical or "local_score_cond_gauss" when some are, each with a DataTypeWarning; otherwise it picks "local_score_BIC". Pass "local_score_BIC" to read whole numbers, such as counts, as numbers. See Notes.

• lambda_value (float) – Multiplier on the BIC penalty of log(n_samples) per parameter: per parent for the linear-Gaussian BIC. Larger values give sparser graphs; 1.0 is the standard BIC. Must be finite and non-negative. BDeu ignores it; for the multinomial BIC it weighs log(n) per free parameter of a variable’s conditional probability table.

• max_parents (int | None) – A cap on every node’s parents, at least 0: no DAG the returned CPDAG represents gives a node more than max_parents parents, and some DAG reaches each node’s count. A move whose CPDAG would break the cap, in either phase, is skipped for the next best move. A cap can stop the search short of the lowest-scoring graph, which GES otherwise finds. None sets no cap; the forward phase then keeps its gate of n_variables / 2. See Notes.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data score: 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 scores 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 score, columns within it get a one-time DataTypeWarning that they may be category codes. With score_func=None, the columns it infers pick the default score.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. No move adds a forbidden edge or removes a required one, and the search starts from the required edges. The result directs every edge the knowledge decides, then Meek’s rules direct the edges that follow, as Knowledge describes.

Returns:

Wraps a cpdag GraphStructure; metadata["algorithm"] is "GES" and metadata["score"] the score of the result; lower is better. For BIC, it is the sum over variables of n * log(residual variance) + lambda_value * log(n) * n_parents; for BDeu, -2 times the log of the structure prior times the marginal likelihood; for the multinomial BIC and the mixed-data scores, see Notes. 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 score_func names no available score.

• 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_parents is not an int.

• ValueError – If lambda_value is negative or non-finite, if a given max_parents is negative or below the parents some graph of the required edges of knowledge gives a node, if max_categories is less than 2, with a categorical score, if data holds a number that is not whole, or, with a mixed-data score, 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 score. With the BIC score, also if a column repeats another up to scale and shift given more than two rows or, with more rows than columns, is a linear function of other columns. One message names every group of copies and says to keep one column per group. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When score_func is None and the columns are all categorical or mixed, so the search scores with BDeu or the conditional Gaussian BIC; when a mixed-data score 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. GES suits data with no unmeasured common causes. Its default BIC score assumes continuous data, where each variable depends linearly on its causes with Gaussian noise; in that setting it is a fast, accurate default. When every column is categorical, the default is "local_score_BDeu", and on mixed data "local_score_cond_gauss"; "local_score_BIC_categorical" and "local_score_degen_gauss" also score those data.

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 lambda_value. A larger value raises the penalty per parent, giving a sparser graph: fewer false edges and more missed ones. 1.0 is the standard BIC. With many variables and about 10 * n_variables samples, lambda_value=2.0 gave more accurate graphs than 1.0 on the published benchmarks. FGES, a fast GES, multiplied the BIC penalty by 4 when searching continuous data on graphs of up to a million variables [2].

How it works. Without max_parents, a variable that already has more than n_variables / 2 parents gains no new parents in the forward phase; the backward phase has no such gate. With max_parents, that gate gives way to the cap. One move can add several parents to a variable at once, and rebuilding the CPDAG after a move can direct or undirect other edges, so the search tests the CPDAG after each move and takes the best move whose CPDAG keeps every variable within the cap. A variable’s count is its directed parents plus the largest set of its undirected neighbors that are adjacent to one another: a graph of the class directs such a set into it and no larger one, since two neighbors that are not adjacent would add a collider. A cap of 1 therefore allows a chain or a tree of undirected edges but not a collider. Moves whose score changes are equal after rounding are chosen by each column’s rank, as in hc(); the column-order guide says what the rank keeps under a new column order and on another machine.

Categorical data. "local_score_BDeu" takes 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. It scores a variable given its parents by the marginal likelihood of its counts in each parent configuration, under a Dirichlet prior with an equivalent sample size of 1 spread evenly over the cells [3]. A structure prior with one expected parent per variable multiplies it, from three variables up. Equivalent graphs score the same, as GES requires. Moves are scored in this process, without the worker processes the BIC search can use. "local_score_BIC_categorical" reads the columns the same way and scores a variable by -2 times the log-likelihood of its counts in each parent configuration, plus lambda_value * log(n) for each of the q * (r - 1) free parameters, with r its categories and q its parents’ joint categories [7]. Both scores are on the scale of the BIC: -2 times a log score.

Mixed data. "local_score_degen_gauss" and "local_score_cond_gauss" score data with continuous and categorical columns; the categorical columns are read as pc() describes. Each is a BIC: -2 times the fitted log-likelihood of a variable given its parents, plus lambda_value * log(n) per parameter; the degenerate Gaussian’s leaves out a constant per variable, which does not change which graph scores best. Equivalent graphs score the same. The degenerate Gaussian score [4] replaces a categorical column of k categories by k - 1 indicator columns and scores every column as Gaussian, by the residual covariance of a linear regression on the parents’ columns; on continuous columns only, it equals "local_score_BIC". The conditional Gaussian score [5] fits the continuous variables of a variable and its parents in each joint category of their categorical ones, with a Gaussian of that category’s own mean and a covariance shrunk toward that of all rows by one pseudo-row [6], and divides by the same fit of the parents alone, so a categorical variable can have continuous parents. It counts parameters for every joint category, so a parent with many categories costs many parameters; on continuous columns only, it differs from "local_score_BIC" by a constant per variable, so it finds the same graph. Moves are scored in this process.

Limitations. Measured on simulated linear-Gaussian data from sparse random graphs (mean degree 2), from 20 to 200 variables with ten rows per variable and a forward-phase gate of 4 parents rather than n_variables / 2; the cap max_parents sets is not benchmarked yet. At lambda_value=1.0 its median errors (SHD) stayed below the empty graph’s up to 150 variables (141 against 150) and passed it at 200 (206 against 200), where lambda_value=2.0 brought them to 43. The median fit took 0.15 s at 20 variables and 40 s at 150, on 128 cores. On the 223-variable ANDES network with 1,000 simulated rows and the default limit, its median errors passed the empty graph’s (612 against 338). With fewer rows than columns the BIC score is unreliable. A parent set of n_samples - 1 or more variables fits its child exactly, so the search never selects one, nor any parent set that leaves its child a residual variance of at most 1e-12 times its variance. 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 at 40 variables were 8 with "local_score_degen_gauss" and 6.5 with "local_score_cond_gauss", against 40, and the median fit took 0.75 s and 1.2 s on 4 cores. On simulated categorical networks with 2 or 3 categories per variable, of the same sizes and rows, "local_score_BDeu" gave 6 against 40 at 40 variables, in 0.60 s. The multinomial BIC is not benchmarked with GES yet: its tests check it against the formula and run GES on small fixed categorical networks.

References

[1]

Chickering, D. M. (2002). Optimal structure identification with greedy search. JMLR 3, 507-554.

[2]

Ramsey, J., Glymour, M., Sanchez-Romero, R., and Glymour, C. (2017). A million variables and more: the Fast Greedy Equivalence Search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and Analytics 3(2), 121-129.

[3] (1,2)

Heckerman, D., Geiger, D., and Chickering, D. M. (1995). Learning Bayesian networks: the combination of knowledge and statistical data. Machine Learning 20(3), 197-243.

[4] (1,2)

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.

[5] (1,2)

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.

[6]

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

[7] (1,2)

Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics 6(2), 461-464.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from andrey import ges
>>> 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)
>>> ges(pd.DataFrame({"x": x, "y": y, "z": z})).structure.oriented_edges()
[('x', 'z', 'directed'), ('y', 'z', 'directed')]

On categorical data, with BDeu:

>>> 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})
>>> ges(data, score_func="local_score_BDeu").structure.oriented_edges()
[('a', 'c', 'directed'), ('b', 'c', 'directed')]

On mixed data, with a category column and the conditional Gaussian score:

>>> 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})
>>> ges(data, score_func="local_score_cond_gauss").structure.oriented_edges()
[('g', 'v', 'directed'), ('u', 'v', 'directed')]