andrey.gfci

andrey.gfci(data, *, score_func='local_score_BIC', lambda_value=1.0, indep_test='fisherz', alpha=0.05, collider_rule='sepsets', categorical=None, max_categories=10, possible_dsep=None, max_depth=None, knowledge=None)

Learn a PAG from observational data with GFCI, a score-then-FCI hybrid.

GFCI [1] starts from the adjacencies of a score-based GES search, removes those that conditional-independence tests reject, and orients the endpoints with FCI’s rules to account for possible latent confounding. Its scored skeleton is usually more accurate than FCI’s test-only one.

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) – Local score driving the GES skeleton phase, as for ges(): "local_score_BIC", the linear-Gaussian BIC; "local_score_BDeu" or "local_score_BIC_categorical" for categorical data, coded by whole numbers; or "local_score_degen_gauss" or "local_score_cond_gauss", the degenerate or conditional Gaussian BIC for mixed data.

• lambda_value (float) – Multiplier on the BIC penalty of log(n_samples) per parameter in the GES phase. Larger values give sparser graphs; 1.0 is the standard BIC. Must be finite and non-negative. BDeu ignores it.

• indep_test (str) – Conditional-independence test for the FCI phase. "fisherz" (Fisher’s z) and "student_t" (Student-t) test partial correlation, for continuous data. "chisq" (Pearson’s chi-square) and "gsq" (G-square, the likelihood ratio) test categorical data, coded by whole numbers, and "degen_gauss" and "cond_gauss" mixed data. "kci", "rcot", "cmiknn", and "gpdc" also detect nonlinear dependence in continuous data, at a far higher cost per test. Each as for pc().

• alpha (float) – Significance level of the conditional-independence test. Smaller values give sparser graphs.

• collider_rule (str) – How GFCI decides whether z is a collider on a triple x - z - y where the tests removed the GES edge between x and y; on the other triples whose ends are not adjacent, the GES CPDAG decides (see Notes). "sepsets" makes it a collider when z is not in the recorded separating set of x and y. "majority" tests x and y given every subset of x’s neighbors and every subset of y’s neighbors, and makes it a collider when fewer than half of the separating subsets (p >= alpha) contain z [3]; at exactly half the triple stays unoriented. Majority runs 2**degree(x) + 2**degree(y) tests per triple.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for the mixed-data score and test: DataFrame column names 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 settings 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 setting, columns within it get a one-time DataTypeWarning that they may be category codes.

• possible_dsep (str | None) – The order of GFCI’s possible-d-separation step [2] [4], which removes edges the adjacency search kept, one at a time, so each removal changes what the step tests after it. "order_free" visits the pairs and nodes by a rank computed from each column’s own data, so the result does not depend on the column order unless two columns are identical. "fast" visits them in column order. None follows andrey.config.tie_break: "order_free" under "rank", "fast" under "position". See Notes and Limitations.

• 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. It caps the tests of the adjacency search, the separating sets for the pairs GES leaves non-adjacent, the possible-d-separation step, the majority rule, and the discriminating-path rule, but not the GES phase.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. The GES phase honors it as ges() does, and the test phase starts from the adjacencies GES found and never removes a required pair. In the PAG, a forbidden a -> b on an adjacent pair puts an arrowhead at a and can leave a circle at b, and a required one is a -> b.

Returns:

Wraps a pag GraphStructure; metadata["algorithm"] is "GFCI". Endpoints carry arrowheads, tails, or circle marks (o); a bidirected edge (<->) marks a probable latent common cause. 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, or indep_test 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 alpha is not in (0, 1), if lambda_value is negative or non-finite, if collider_rule or possible_dsep is not a supported choice, if score_func and indep_test read different data (continuous, categorical, or mixed), or if data or categorical does not suit them, as for pc(). 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 a given max_depth is negative, or if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When the mixed-data score and 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.

Notes

GFCI first runs ges() and keeps its CPDAG. FCI’s adjacency search then starts from the CPDAG’s adjacencies instead of a complete graph: it removes an edge x - y at the smallest size at which a subset of x’s or y’s neighbors makes the two independent, and records the union of those subsets as their separating set. On a triple x - z - y whose ends are not adjacent, z is a collider when the CPDAG has x -> z <- y; when the CPDAG had the edge x - y, the separating set decides, as collider_rule describes. FCI’s possible-d-separation step can remove more edges, after which the colliders are decided again. No step adds an edge, so the PAG’s adjacencies are a subset of the GES CPDAG’s. Zhang’s rules R1-R10 [2] then orient the endpoints, in rounds as in fci(); an endpoint no rule resolves stays a circle. A rule that needs the separating set of a pair GES left non-adjacent uses the union of the smallest subsets of the two endpoints’ CPDAG neighbors that make the pair independent, which does not depend on the column order. Under andrey.config(tie_break="position"), the rules fire one at a time in column order, such a pair’s separating set is the first subset found in column order, and possible_dsep defaults to "fast".

Rule R5 turns a cycle of circle-circle edges into undirected edges only when the whole cycle, closing edge included, is uncovered, that is, no two nodes two steps apart on it are adjacent. It considers every circle path, up to a fixed search bound, and rule R8 checks each pair of endpoints in both orders. GFCI shares FCI’s possible-d-separation step [2] [4], in the order possible_dsep selects, and its separating sets; the Notes of fci() describe both and the reasons for keeping them.

Limitations. Measured with possible_dsep="fast" on simulated linear-Gaussian data from random graphs with mean degree 2 or 4 and one variable in ten or one in five hidden, from 20 to 800 variables counting the hidden ones, with 10 or 50 rows per variable, on 4 cores. Its median endpoint errors, each mismatched tail, arrowhead, or circle, stayed below the empty graph’s wherever its fits finished (197 against 444 at 200 variables with mean degree 2, one variable in ten hidden, and ten rows per variable). With those settings the median fit took 0.53 s at 100 variables, 2.6 s at 200, and 132 s at 800. Every fit at 800 variables finished within 7,200 s in four of the eight settings; with mean degree 4, one variable in five hidden, and 50 rows per variable, none did, and three of five fits at 400 ran past 3,600 s. The default, "order_free", took 1.00 to 1.02 times as long as "fast" up to 200 variables and 0.87 to 1.03 times at 400 and 800, and its median endpoint errors were within 2% of "fast"’s. On some graphs from 400 variables its order needs more tests than the column order: on one with 400 variables, mean degree 2, one variable in five hidden, and 50 rows per variable, it took 407 s against 120 to 336 s for "fast" over five column orders. possible_dsep="fast" gives the column-order mode. It runs a full ges() search before any test, so it takes at least as long as GES on the same data. Its BIC score and Fisher-Z tests assume linear relations with Gaussian noise. 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 GES search never selects one, nor any parent set that leaves its child a residual variance of at most 1e-12 times its variance. Rule R5 has the search bound described under fci(). With the mixed-data score and test, it is experimental and not benchmarked yet; its tests run it on 5 variables with 3,000 rows.

References

[1]

Ogarrio, J. M., Spirtes, P., and Ramsey, J. (2016). A hybrid causal search algorithm for latent variable models. PGM, PMLR 52, 368-379.

[2] (1,2,3)

Zhang, J. (2008). On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence 172(16-17), 1873-1896.

[3]

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

[4] (1,2)

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

Examples

>>> import numpy as np
>>> from andrey import gfci
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = x + 0.1 * rng.standard_normal(500)
>>> gfci(np.column_stack([x, y])).structure.kind
'pag'