andrey.inter_iamb

andrey.inter_iamb(data, *, target=None, structure='cpdag', alpha=0.05, indep_test=None, categorical=None, seed=None, max_depth=None, knowledge=None)

Learn a variable’s Markov blanket with interleaved IAMB (Inter-IAMB).

Inter-IAMB [1] is IAMB with its shrink step after every addition, so a member that becomes independent leaves before it enlarges the later conditioning sets. Without a target, every variable’s blanket is learned and joined into the graph structure names, by default the CPDAG of the Grow-Shrink structure learner [2] [3].

Parameters:

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

• target (str | int | None) – The variable whose set to learn: a column label (a string) or a 0-based position (an integer). None learns every variable’s set and keeps a pair when each variable is in the other’s set.

• structure (str) – The graph returned when target is None. "cpdag": the CPDAG, the causal graphs the tests cannot tell apart, whose directed edges they all share. "skeleton": the same adjacencies with no directions, each variable joined to its parents and children. "moral": the moral graph, each variable joined to every member of its Markov blanket, a child’s other parents included, for feature selection or as a superset of the skeleton; an edge there does not say that the two variables are adjacent. With target, the result is the target’s Markov blanket, whatever structure is.

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

• indep_test (str | None) – Conditional-independence test, any of pc()’s: "fisherz", "student_t", or "normal_scores" for continuous data, "chisq" or "gsq" for categorical data, "degen_gauss" or "cond_gauss" for mixed data, or a nonlinear test. None picks from the declared column types, as for pc(). metadata["indep_test"] names the test the fit used: the one picked from the declared column types when this is None.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data test, as for pc(). A continuous test raises ValueError on a declared column.

• seed (int | None) – None, the default, is the configured seed, ANDREY_SEED or 0 when it is unset, not a fresh random one, so a fit repeats. Seed for the random draws of indep_test, which only "rcot", "cmiknn", and "gpdc" make; the other tests ignore it. With those tests, metadata["seed"] records the seed used.

• 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 variables that a larger conditioning set would remove. Every test conditions on the whole current blanket, so the blanket stops growing once it holds more than max_depth members, counting those knowledge requires. A collider test conditions on the middle node as well, so with 0 no triple is decided and the CPDAG directs only the edges knowledge decides.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. Each end of a required edge is in the other’s blanket from the start, and no test removes it. A pair the knowledge makes non-adjacent, forbidden both ways or outside allowed, stays a candidate, since a child’s other parent belongs to the blanket without being adjacent, and is left out of the result. Under structure="cpdag", the result directs every edge the knowledge decides, by a tier or an edge forbidden one way, before the colliders, which never reverse it, and 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. The undirected outputs do not read tiers or edges forbidden one way.

Returns:

Without target, wraps the graph structure names: a cpdag GraphStructure for "cpdag", whose undirected edges the tests cannot orient, or an undirected one, whose edges have no arrowhead, for "skeleton" and "moral". With target, an undirected one whose edges join target to each member of its Markov blanket, which holds its parents, its children, and its children’s other parents; metadata["target"] names it and metadata["markov_blanket"] lists the members in the order they joined, by column label or 0-based position. metadata["algorithm"] is "inter_iamb". Under a mixed-data setting, metadata["categorical"] lists the categorical columns it read. With a randomized indep_test, metadata["seed"] is the seed used.

Return type:

StructureOutput

Raises:

• TypeError – If target is neither a column label nor a position, categorical is not a list of column names and positions, or 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 target names a column data lacks, or if data or categorical does not suit the test, as for pc(). With "fisherz" or "student_t", also if a column repeats another up to scale and shift. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists, or under structure="cpdag" orients the edges the search kept into a directed cycle.

Warns:

DataTypeWarning – As for pc(), when a mixed-data test finds no declared categorical column or a setting reads columns that look like category codes as numbers; once per process.

See also

andrey.iamb

Shrinks once, after the blanket stops growing.

Notes

When to use. Inter-IAMB suits i.i.d. observational data with no unmeasured common causes, in which every independence follows from the graph (faithfulness) [1]. Its default Fisher-Z test assumes continuous data with linear effects and Gaussian noise; the other tests read categorical or mixed data or detect nonlinear effects, as pc() describes. It is designed for many variables whose Markov blankets are small relative to the rows [1]: shrinking after each addition aims to keep the current blanket, and so the conditioning sets, small at every step. With target, it returns an undirected graph joining target to each member of its Markov blanket, which holds its parents, its children, and its children’s other parents; without, the graph structure names: the CPDAG by default, or the undirected skeleton or moral graph. For the parents and children alone, use mmpc(); when unmeasured common causes are plausible, use fci().

How it works. Each step adds the candidate most dependent on the target given the current blanket, among those whose p-value is at or below alpha, then visits the members in the order they joined and drops each one the test finds independent of the target given the other current members. The search stops when no candidate qualifies, or when the blanket after a step repeats one it held after an earlier step, which on sample data could otherwise make it cycle. With exact tests, the result is the true Markov blanket [1].

Candidates whose strengths agree to nine decimals go to the column whose rank, computed from the data, comes first, as for hc(), so the sets do not depend on the column order unless two columns are identical. The strength is the absolute partial correlation for "fisherz" and "student_t", which keeps its order where p-values underflow to zero, and the p-value for the other tests.

From blankets to a graph. Without a target, a pair is kept when each variable is in the other’s blanket, which gives the moral graph (structure="moral"). The skeleton drops each pair x - y that tests independent given some subset, the empty set included, of the smaller of x’s blanket without y and y’s blanket without x, or of either when both hold as many members [2] [3]: a child’s other parent is in the blanket without being adjacent, and such a subset separates the two. Every pair is tested against the moral graph, so no verdict depends on another. The CPDAG decides each triple x - z - y of the skeleton with x and y not adjacent: z is a collider, x -> z <- y, when x and y test dependent given z together with every subset of the smaller of their blankets without x, y, and z [2], and a non-collider when one such set separates them. The colliders are oriented as pc() orients them, the most clearly dependent first: by the smallest absolute partial correlation over the collider’s tests for "fisherz" and "student_t", and by the largest p-value, smallest first, for the other tests. When two colliders that rank equal to nine decimals dispute an edge, it stays undirected. Meek’s rules [4] then orient the edges that follow, in rounds, as for pc(), so the result does not depend on the column order; under andrey.config(tie_break="position") the colliders and the rules follow the column order instead. With exact tests, the result is the true CPDAG [2].

Missing values. A NaN in data marks a missing value. Each test reads only the rows where all of its variables are observed (test-wise deletion) [5] [6], so it keeps rows that deleting every incomplete row first would drop; metadata["rows_tested"] gives the fewest and the median rows a test read. When values are missing completely at random, the rows a test keeps are a random sample, so each test keeps its level [6]. When whether a value is missing depends on the variables, the rows a test keeps are a selected sample, and the tests can find dependence where the population has none and keep extra edges [6]. A test raises ValueError when it has fewer rows than it needs, |S| + 3 for "fisherz" and "student_t" and 2 for the other tests, or when a continuous column is constant on its rows.

Limitations. Measured on simulated linear-Gaussian data with 20, 50, and 100 variables, 1,000 and 5,000 rows, and mean degree 2 and 4, with Fisher-z tests at alpha 0.05, on one core. In every setting the mean errors of the CPDAG, the skeleton, and the moral graph stay below the empty graph’s; the CPDAG is closest at 20 variables, 1,000 rows, and mean degree 4, with 24 pairs wrong against 40 for the empty graph. A fit took 0.06 to 8.9 s, the longest at 100 variables, 5,000 rows, and mean degree 4. Tested with exact d-separation on 15 variables. It conditions on the whole current blanket, so each test needs more rows as the blanket grows. The skeleton and the CPDAG test every subset of the smaller blanket for each pair of the moral graph and each unshielded triple, and a collider takes all of them, so a triple whose smaller blanket holds m variables takes up to 2**m tests. On simulated linear-Gaussian data with 50 variables, 5,000 rows, and mean degree 4, the smaller blankets held at most 16 variables, the collider tests took 0.06 to 0.7 s on five datasets, and learning the blankets 0.9 to 1.3 s. Missing values are not benchmarked yet; the tests run them on up to 20 variables with up to a fifth of the values missing completely at random.

References

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

Tsamardinos, I., Aliferis, C. F., and Statnikov, A. (2003). Algorithms for large scale Markov blanket discovery. FLAIRS, 376-380.

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

Margaritis, D. (2003). Learning Bayesian network model structure from data. PhD thesis, Carnegie Mellon University, CMU-CS-03-153.

[3] (1,2)

Margaritis, D., and Thrun, S. (2000). Bayesian network induction via local neighborhoods. Advances in Neural Information Processing Systems 12, 505-511.

[4]

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

[5]

Strobl, E. V., Visweswaran, S., and Spirtes, P. L. (2018). Fast causal inference with non-random missingness by test-wise deletion. International Journal of Data Science and Analytics 6(1), 47-62.

[6] (1,2,3)

Tu, R., Zhang, C., Ackermann, P., Mohan, K., Kjellstrom, H., and Zhang, K. (2019). Causal discovery in the presence of missing data. Proceedings of Machine Learning Research 89, 1762-1770.

Examples

>>> import warnings
>>> import numpy as np
>>> from andrey import inter_iamb
>>> rng = np.random.default_rng(0)
>>> a, b = rng.standard_normal((2, 1000))
>>> c = a + b + 0.5 * rng.standard_normal(1000)
>>> with warnings.catch_warnings():
...     warnings.simplefilter("ignore")
...     out = inter_iamb(np.column_stack([a, b, c]), target=0)
>>> sorted(out.metadata["markov_blanket"])
[1, 2]