andrey.igsp

andrey.igsp(data, env=None, *, targets=None, alpha=None, invariance_alpha=None, indep_test=None, invariance_test=None, n_starts=10, seed=None, categorical=None, knowledge=None)

Learn an interventional essential graph with IGSP from settings with known targets.

IGSP (interventional greedy sparsest permutation) learns from data pooled from an observational environment and interventional ones, the settings, whose targets, the variables each one intervened on, are known [1] [2]. It searches variable orders for the DAG with the fewest edges and, among those, the fewest edges the settings contradict: each order gives the DAG of the edges its conditional-independence tests keep on the observational rows, and invariance tests between each setting and the observational rows say which directions the settings support. It returns the DAG’s interventional essential graph under the targets.

Parameters:

• data (ArrayLike) – Observed data, one row per sample, from every environment. Coerced to float64. A DataFrame’s column names become the structure’s labels. A column that env or targets names is not a variable: it leaves the data before the fit, and categorical, knowledge, targets, and the result’s positions count the other columns. A list of DataFrames or 2-D arrays gives each environment’s rows apart, labeled by position 0, 1, ..., and a mapping labels each by its key: the rows are stacked in that order, env stays None, and targets maps those labels to their variables.

• env (ArrayLike | str | None) – The environment each row comes from: one label per row, strings or numbers, not both, or the name of a column of the DataFrame data that holds them. None takes the environments from targets: its distinct values.

• targets (Mapping | ArrayLike | str | None) – The variables intervened on in each environment, named by column label or 0-based position. A mapping {label: [variables]} takes each label of env to a list of its targets; [] marks an observational environment, and every label of env needs a list: None, or a label the mapping leaves out, raises ValueError, as a key that is not a label of env does. A string names a column of the DataFrame data whose cells give each row’s targets: a label or a position, a list of them, or an empty cell (None, NaN, "", []) for an observational row. An array is True (or 1) where the row’s value of the variable was set by an intervention. With a column or an array, every row of one environment of env must give the same targets; without env, rows with the same targets are one environment. At least one environment must be observational and one not, and every environment needs at least 3 rows. See Notes.

• alpha (float | None) – Significance level of the conditional-independence tests, in (0, 1): two variables are independent when a test’s p-value exceeds it. Smaller values keep more edges out. None uses one over the number of observational rows, 0.001 at 1,000 rows, so the level falls as the rows grow. Also the level of the invariance tests when invariance_alpha is None. metadata["alpha"] gives the level used. See Notes for the choice of the default.

• invariance_alpha (float | None) – Significance level of the invariance tests, in (0, 1): a variable’s distribution given a set of others differs between a setting and the observational rows when the test gives a p-value at most this level. None uses the level of the conditional-independence tests, metadata["alpha"].

• indep_test (str | None) – The conditional-independence test, on the observational rows. For continuous data, on the partial correlation: "fisherz", "student_t", and "normal_scores"; for categorical data coded by whole numbers: "chisq" and "gsq"; for mixed data, with the categorical columns named by categorical: "degen_gauss" and "cond_gauss"; for nonlinear effects in continuous data, each far slower: "kci", "rcot", "cmiknn", and "gpdc", as for pc(). None picks from the declared column types: "gsq" when every column is categorical, "cond_gauss" when some are, else "fisherz". metadata["indep_test"] names the test the fit used.

• invariance_test (str | None) – The invariance test, which asks whether a variable’s distribution given a set of others is the same in a setting as in the observational rows. "gaussian" compares the linear regressions of the two groups of rows: their coefficients and intercept by the Chow test and their residual variances by an F-test, each at half the level. "indicator" runs indep_test between the variable and a 0/1 indicator of the setting’s rows given the others. None picks "gaussian" under "fisherz", "student_t", and "normal_scores", else "indicator". "gaussian" reads every column as numbers, so with a categorical or mixed-data indep_test it raises ValueError. metadata["invariance_test"] names the test the fit used. See Notes.

• n_starts (int) – Number of searches, at least 1, each from a minimum-degree order of the variables: the first breaks ties by each column’s rank, and each later one draws them, and the order in which it visits moves, from seed and its own index. Returns the best-scoring search, the one that found the fewest edges, then the fewest contradicted ones; a later one replaces an earlier one only when it is strictly better. See Notes for the choice of the default.

• 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 searches after the first, and for the random draws of "rcot", "cmiknn", and "gpdc". With n_starts=1 and another test, the result does not depend on it.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data test: column labels (strings) or 0-based positions. None takes a DataFrame’s category, bool, and text columns, never a column by its values. [] marks none. A categorical test reads every column as categorical, and a continuous test raises ValueError on a declared column, named here or by its type.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. Every DAG the search visits has the required edges and none of the forbidden ones, and every order it visits puts the cause of each required edge, and each earlier tier, first. 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, the interventional essential graph of the DAG found under the observational environment and each setting’s targets: directed marks (->) are the edges that every DAG of the class orients alike, undirected (--) the ones it leaves open. metadata["algorithm"] is "igsp"; metadata["score"] is the DAG’s number of edges, lower is better, and metadata["n_contradicting"] its number of I-contradicting edges, the second part of the objective; metadata["indep_test"] and metadata["invariance_test"] name the tests used, and metadata["alpha"] is the level of the conditional-independence tests; metadata["seed"] is the seed used when the fit draws random numbers. Under a mixed-data test, metadata["categorical"] lists the categorical columns it read.

Return type:

StructureOutput

Raises:

• TypeError – If an argument has the wrong type, such as a string for a number; the message names the method, the argument, what it accepts, and the value. Also if categorical is not a list of column names and positions, or knowledge is neither None nor a Knowledge.

• ValueError – If an argument lies outside its range or is not one of its choices; the message names the method, the argument, what it accepts, and the value. If env or targets is malformed, or leaves an environment’s targets unknown, as gies() lists; if no environment is observational, or every one is, as with env and targets both None; if an environment has fewer than 3 rows; if invariance_test="gaussian" comes with a categorical or mixed-data indep_test; if data holds a missing or infinite value, a column constant on every row or on the observational rows, or does not suit the test, as for pc(). With a partial-correlation test, also if a column repeats another up to scale and shift on the observational rows; one message names every group of copies. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When a column holds at most 10 distinct whole numbers and so may hold category codes, which a continuous test reads as numbers. Each message warns once per process.

See also

andrey.gies

Learns the same kind of graph from the same inputs by a score search.

andrey.ut_igsp

Learns the targets too, when some are unknown.

andrey.grasp

Searches variable orders on observational data.

Notes

When to use. IGSP suits continuous data from an observational environment and settings whose targets are all known, with linear effects, Gaussian noise, and no unmeasured common causes [1] [2]. Its authors claim it for perfect interventions, which cut a variable from its causes, and for imperfect ones, which change its mechanism, such as a smaller weight on each cause or a perfect intervention that succeeds on some rows only. Their simulations use 10 or 20 variables with one or two settings of 4 to 10 targets and 1,000 to 100,000 rows per setting [1], or up to 20 single-variable or 45 two-variable settings [2]. It returns the interventional essential graph, the same kind of graph gies() returns from the same inputs; the authors show a 7-variable DAG where the score search of GIES is not consistent and IGSP is. Measured in the authors’ settings against their reference implementation, both at alpha=1e-3, the default’s level at 1,000 observational rows, it was worse in none of 114 comparisons and better in 10, its errors lower on two sets of datasets; all 10 are in settings of two or more targets, where that implementation’s test departs from the paper, and against a copy corrected to the paper it was worse in none and better in none. At its default, against the same implementation at 1e-3, it was worse in none of the 114 and had fewer errors in 17 on the second set (see Choosing alpha and n_starts). On the same data, gies() and ut_igsp() told every target often had a lower SHD (see Limitations). When some targets are unknown, use ut_igsp(); with no observational environment, use gies().

Choosing alpha and n_starts. The defaults come from simulated data in the authors’ settings [1] [2]: linear effects and Gaussian noise on 7 to 20 variables, with one or two settings of up to 10 targets or up to 45 settings of one or two targets, perfect, inhibiting, or imperfect, and 1,000 to 100,000 rows per setting: 87 groups of 20 datasets, one fit each, one set of datasets with no rerun. A candidate’s SHD counts as lower or higher than another’s in a group when the 95% bootstrap interval of the mean difference, paired by dataset, lies wholly below or above 0. Among the levels the authors swept, alpha=1e-3 was worse in none of the groups and lower than 1e-4 in 4 groups, 5e-3 in 16, 1e-2 in 20, and 0.05 in 83, by 0.2 to 10.7 edges; 5e-4 tied with it in every group. On 20 variables at mean degree 3 with 8 perfect single-variable settings and 1,000 or 10,000 rows per setting, outside these groups, no level was lower than 1e-3. At alpha=1e-3, 10 starts and 5 were each worse in none of the groups and tied with each other in every group; against one start, 10 were lower in 21 groups and 5 in 19, by 0.4 to 1.7 edges, at a median of 0.10 s per fit against 0.07 s for 5. On fresh datasets, 5 fits each, against the reference implementation at 1e-3, alpha=1e-3 and one over the observational rows, both with 10 starts, were each worse in none of 114 comparisons and had fewer errors in 17 each when interval ends that are floating-point residue of 0 are read as 0, or in 18 and 17 as computed; a rule written before the fits gives such a tie to one over the observational rows, the default of ut_igsp(). It equals 1e-3 at 1,000 observational rows, and in the 50 groups at 10,000 and 100,000 rows per setting it had fewer SHD errors than 1e-3 in 1, by 0.28 edges, and more in none. Of its 17, 15 are against the reference implementation, 12 of them in settings of two or more targets, by 0.3 to 3.58 edges, and 2 against the corrected copy.

How it works. A variable order gives a DAG, its minimal I-MAP: i -> j for i before j unless the conditional-independence test finds them independent given every other variable before j, on the observational rows (all environments whose targets are [], pooled). An edge i -> j is I-contradicting when the settings say it points the wrong way [2]: with single-variable settings, when j’s distribution is the same under every setting that targets i alone, or i’s changes under one that targets j alone; with larger target sets, when, given some set of j’s other neighbors, j’s distribution is the same under every setting that targets i and not j, or, given every set of i’s other neighbors, i’s changes under some setting that targets j and not i. The objective is the number of edges, then the number of I-contradicting edges. The search starts from the order of a minimum-degree elimination of the graph of pairs dependent given all other variables, and moves by reversing I-covered edges: i -> j where j’s parents are i and i’s parents, and j’s distribution does not change under a setting that targets i alone. A depth-first search visits the reversals of I-contradicting edges first and backtracks at a dead end, up to 5 reversals away; it restarts from a DAG with fewer edges whenever it finds one, and keeps the DAG with the lowest objective it saw [1] [2]. The result is the best of n_starts such searches, returned as its interventional essential graph under the targets [3].

Invariance tests. Under "gaussian", a variable given a set of others is invariant when the linear regressions with an intercept agree in the two groups of rows: the Chow test of equal coefficients and intercept and the F-test of equal residual variances [4], exact under Gaussian noise, each at half the level. With no other variable it compares the means and variances, so it sees a perfect intervention that keeps the mean. A group with no more rows than the regression has parameters leaves nothing to compare, and the variable reads as invariant. Under "normal_scores" the regressions read every column’s normal scores over all rows. Under "indicator", a variable is invariant when indep_test finds it independent of a 0/1 indicator of the setting’s rows given the others, on the rows of the setting and the observational ones [5]; under "fisherz" and "student_t" this sees a change of the mean alone, which is why "gaussian" is their default.

Environments and targets. Every environment whose targets are [] is observational, and their rows are pooled. Every other environment is a setting, compared with the observational rows alone, never with another setting. With env, two environments with the same targets stay two settings; without env, rows with the same targets are one setting. A setting’s targets are used as given: a target whose distribution the intervention leaves unchanged orients nothing.

Column order. Each column gets a tie rank computed from the data on every row, as for hc(), and the search runs on the columns sorted by it: the first search breaks ties between orders and between moves by the rank, and the later ones shuffle them. Reordering the columns relabels the result, and the labels of the environments and the order of the rows change nothing for the tests that draw no random numbers. The column-order guide says what the rank keeps.

Missing values. A NaN in data raises ValueError: every invariance test regresses or tests a variable on all of its conditioning variables in each environment. Drop the incomplete rows first.

Limitations. Measured on simulated data in the authors’ settings, the 87 groups above with 5 fits per dataset, one core per fit, against the reference implementation with its 5 runs at alpha 1e-4, 1e-3, and 1e-2, and in the 27 groups with settings of two or more targets also against a copy of it corrected to the paper’s test. In every group its SHD was below the empty graph’s and every fit finished. At the three levels, with 10 starts and with 5, it was worse in none. At alpha=1e-3 with 10 starts, the 10 comparisons it was better in were all against the uncorrected implementation; the corrected copy tied it in 25 of the 27 groups, and in the other 2 its SHD was lower on one set of datasets only. At alpha=1e-4 with 5 starts, on 10 variables with 2 settings of 4 targets and 1,000 rows per setting, its SHD was 2.51 against the corrected copy’s 2.12, an interval whose lower end is 0, and 3.07 against 2.76, higher, on fresh datasets. There the objective often ties the true DAG with a wrong one, and the first start, the minimum-degree order with its ties broken by the column ranks, is the same in every fit, so the five fits repeat one result where the reference implementation’s random starts vary; at 10 starts that group ties. Unless a level is named, the numbers below are at alpha=1e-3, the default’s level at 1,000 observational rows. On the same data, gies() had a lower SHD in 35 of the 87 groups, by 0.35 to 2.63, and a higher one in 6, and ut_igsp() told every target a lower one in 27, by 0.38 to 1.95, and a higher one in none. A perfect intervention that leaves its target’s distribution unchanged, such as cutting the causes of a variable that has none while keeping its noise, orients nothing, and the search then reverses some edges out of that target: 823 of 3,960 such edges over the 87 groups, against 838 for the reference implementation and 16 for gies(). Outside the authors’ setting, one fit per dataset: with 100 rows per setting on 10 variables its SHD was 4.65 against the reference implementation’s 4.05 and the empty graph’s 8; at mean degree 4 on 20 variables it was 11.8 against 12.65 and 40; with one unknown target per setting, 0.85 and 0.95 against 0.05 and 0.15 for ut_igsp() and 15 for the empty graph. With additive nonlinear effects on 10 variables its SHD was 2.05 and 2.3 with the Gaussian tests and 5.0 and 1.15 with indep_test="kci", at 200 and 500 rows per setting, against 8, a kernel fit taking 1.7 s and 13.8 s; on categorical data with indep_test="gsq", 6.4 at 1,000 rows per setting and 0.9 at 10,000, against 10. With 10 single-variable settings of 1,000 rows, its SHD was 3.9 at 50 variables, 15.3 at 100, and 42.9 at 200, against 38, 75, and 150, a median fit taking 1.6 s, 33 s, and 726 s and 0.5 GB of memory at 100; at 400 no fit finished within an hour, while gies() reached 0.2, 0.1, and 0.3 at 50 to 200 variables in under 2 s. Its fits take longer than the reference implementation’s on one core of the same processor, a median of 1.8 times at 1,000 rows per setting and 4.3 times at 10,000, its time growing with the number of rows; where the search dominates it is faster, 1.6 s against 4.3 s at 50 variables. With settings of more than one target, each edge’s test runs over the subsets of a neighborhood, so dense graphs are slow. The search’s cost grows with the number of I-covered edges to the power of its depth, and each start tests every pair given all other variables, so time grows fast with the number of variables. Unmeasured common causes are outside its model and not measured. One level serves both tests unless invariance_alpha is given. The numbers are a reference, not a hard limit: they change with the data generator, the density, the method’s settings, and the machine. In bootstrap(), the averaged graph is made a CPDAG without the targets, so an edge that only the settings orient comes back undirected there; each fit and the edge frequencies keep its direction.

References

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

Wang, Y., Solus, L., Yang, K. D., and Uhler, C. (2017). Permutation-based causal inference algorithms with interventions. NeurIPS 30. https://proceedings.neurips.cc/paper/2017/hash/275d7fb2fd45098ad5c3ece2ed4a2824-Abstract.html

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

Yang, K. D., Katcoff, A., and Uhler, C. (2018). Characterizing and learning equivalence classes of causal DAGs under interventions. ICML, PMLR 80, 5537-5546. https://proceedings.mlr.press/v80/yang18a.html

[3]

Hauser, A., and Buhlmann, P. (2012). Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. JMLR 13, 2409-2464.

[4]

Chow, G. C. (1960). Tests of equality between sets of coefficients in two linear regressions. Econometrica 28(3), 591-605.

[5]

Heinze-Deml, C., Peters, J., and Meinshausen, N. (2018). Invariant causal prediction for nonlinear models. Journal of Causal Inference 6(2).

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from andrey import igsp
>>> rng = np.random.default_rng(0)
>>> env = np.repeat(["observed", "set y"], 2000)
>>> x = rng.standard_normal(4000)
>>> y = np.where(env == "set y", rng.standard_normal(4000), x + rng.standard_normal(4000))
>>> out = igsp(pd.DataFrame({"x": x, "y": y}), env, targets={"observed": [], "set y": ["y"]})
>>> out.structure.oriented_edges()
[('x', 'y', 'directed')]