andrey.mmpc

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

Learn a variable’s parents and children with max-min parents and children (MMPC).

MMPC [1] finds the variables adjacent to a target in the causal graph, its parents and children, without telling them apart, by conditional-independence tests. Every variable’s set together gives the skeleton mmhc() searches, which structure="cpdag" orients into a CPDAG by the collider tests 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. "skeleton": the adjacencies with no directions, each variable joined to its parents and children, the graph the sets give. "cpdag": the CPDAG, the causal graphs the tests cannot tell apart, whose directed edges they all share, which tests every unshielded triple for a collider and costs more. A variable’s parents and children hold no child’s other parents, so there is no moral graph. With target, the result is the target’s parents and children, 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. A test conditions on a subset of the current set; the set itself may grow larger. 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 set from the start, and no test removes it; a pair the knowledge makes non-adjacent, forbidden both ways or outside allowed, is never a candidate. 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: an undirected GraphStructure for "skeleton", whose edges have no arrowhead and join each pair in each other’s set, or a cpdag one for "cpdag", whose undirected edges the tests cannot orient. With target, an undirected one whose edges join target to each of its parents and children; metadata["target"] names it and metadata["parents_and_children"] lists the members in the order they joined, by column label or 0-based position. metadata["algorithm"] is "mmpc". 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.mmhc

Orients the skeleton by hill climbing.

andrey.iamb

The Markov blanket, which adds the children’s other parents.

Notes

When to use. MMPC 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, each with few parents and children [1], since each candidate is tested given subsets of the members found so far. With target, it returns an undirected graph joining target to each of its parents and children, without telling the two apart; without, the graph structure names: the undirected skeleton by default, or the CPDAG. For the Markov blanket, which adds the children’s other parents, use iamb(); to orient the skeleton by a score search, use mmhc().

How it works. A candidate’s association with the target is its weakest over every conditioning set drawn from the current members: the largest p-value among those tests. The forward phase adds, one at a time, the candidate whose weakest association is the strongest, and drops for good every candidate that some conditioning set separates (a p-value above alpha). The backward phase visits the members in the order they joined and drops each one some subset of the other members separates. A member stays only when the target is in that member’s own set as well, the symmetry check of [1], which removes a descendant of the target that no subset of the target’s set separates from it. With exact tests, the result is the true set of parents and children [1].

Two candidates whose largest p-values agree to nine decimals are compared by their smallest strength, and candidates whose strengths also agree 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. Where the p-values of several candidates are zero to that precision, the candidate with the strongest statistic joins first, as the tie rule of [1] orders them.

From sets to a graph. Without a target, a pair is kept when each variable is in the other’s set, which gives the skeleton (structure="skeleton"). 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] [3], and a non-collider when one such set separates them. A variable’s blanket here is its set together with its members’ sets, before the symmetry check, which holds its children’s other parents. 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 Fisher-z tests at alpha 0.05, on one core: 20, 50, and 100 variables with 1,000 and 5,000 rows and mean degree 2 and 4, with the columns standardized and not, and 200 variables with 1,000 and 5,000 rows. Wherever fits finished, the mean skeleton errors stay below the empty graph’s; the skeleton is closest to it on unstandardized data at 100 variables, 1,000 rows, and mean degree 2, with 56.9 pairs wrong against 100, and the CPDAG with 78.1 against 100. The median skeleton fit took 0.01 to 1.6 s up to 100 variables and 0.3 to 3.3 s at 200 variables. The CPDAG fit took up to 30 s where it finished and reached the 600 s cap at 100 variables, mean degree 4: on 6 of 20 datasets at 1,000 rows and 19 of 20 at 5,000 on standardized data. At 200 variables and 5,000 rows, MMPC’s forward step often faces candidates whose dependence is too strong to rank by p-value, and it admits the strongest, which is a true neighbor 99.7% of the time. Borderline false edges between unrelated pairs, with a statistic just past the cutoff, then survive: about 2 more per graph than with an arbitrary order, which happens to admit a strongly correlated non-neighbor that removes them. At mean degree 2 the skeleton there makes 75.5 errors on average against 73.9 for bnlearn’s, and 86.9 against 84.2 on 10 fresh datasets. The continuous benchmark does not cover the categorical and mixed-data tests, which only the discrete-network lanes cover. Tested with exact d-separation on 15 variables. It tests each candidate against every subset of the current members, so its cost grows exponentially with the size of the set; the symmetry check runs the search once more for each member. The CPDAG tests every subset of the smaller blanket for each unshielded triple, and a collider takes all of them. These blankets reach two steps from a variable, so they are larger than a Markov blanket, and 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 collider tests took 1.6 to 2.2 s on four datasets and 277 s on a fifth, whose smaller blanket held 24 variables (38 million tests), and learning the sets 0.13 to 0.18 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,5,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.

[2] (1,2,3)

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 mmpc
>>> 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 = mmpc(np.column_stack([a, b, c]), target=0)
>>> out.metadata["parents_and_children"]
[2]