andrey.gnies

andrey.gnies(data, env=None, *, targets=None, penalty_discount=1.0, ebic_gamma=1.0, max_parents=None, knowledge=None)

Learn an interventional essential graph and its unknown targets with greedy search (GnIES).

GnIES learns from data pooled from several environments that share one linear model with Gaussian noise and differ in the noise variances of some variables, the targets, which it need not be told [1]. It learns the targets together with the graph’s interventional equivalence class, so every edge at a target gets its direction from the data.

Parameters:

• data (ArrayLike) – Observed data, one row per sample. 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 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, when given, maps those labels to their known variables.

• env (ArrayLike | str | None) – The environment of each row: one label per row, strings or numbers, not both, or the name of a column of the DataFrame data that holds them. GnIES needs two or more environments, each with at least two rows, and at least as many rows as variables once each environment’s mean is taken out: n_samples - n_environments >= n_variables. One of them may be observational, but none has to be. None takes the environments from targets, as a column or a 0/1 matrix.

• targets (Mapping | ArrayLike | str | None) – The variables known to be intervened on, by column label or 0-based position. None, the setting GnIES is built for, leaves every environment’s targets unknown, and the search learns them. A mapping {label: [variables]} lists the known targets of each label of env; None, or a label the mapping leaves out, marks an environment whose targets are unknown, and [] adds none. When the labels of env are numbers, a text key names the label it spells, "1" the label 1. GnIES fits one set of targets for all environments: the variables the mapping lists stay targets, and while some environment’s targets are unknown the search may add more. When every environment’s targets are known (a mapping that gives each a list, a string naming a column of the DataFrame data whose cells give each row’s targets, or a 0/1 array that is True where the row’s value of the variable was intervened on), the targets are fixed and only the graph is searched.

• penalty_discount (float) – Multiplier on the BIC penalty of log(n_samples) per parameter: per parent, and per extra noise variance of a target, n_environments - 1 each. n_samples counts every row. Larger values give sparser graphs and fewer targets; 1.0 with ebic_gamma=0.0 is the standard BIC. Must be finite and non-negative.

• ebic_gamma (float) – Weight on the extended-BIC term: each variable’s score adds 2 * ebic_gamma * log C(n_variables - 1, n_parents), which grows with the number of parent sets a variable can choose from. Larger values give sparser graphs; 0.0 is the standard BIC. The term charges parents only, not targets. Must be finite and non-negative. See Notes.

• max_parents (int | None) – A cap on every node’s parents, at least 0: no DAG the returned graph represents gives a node more than max_parents parents. Each search over graphs skips a step that would give a node more for the next best step, so a cap can stop it short of the lowest-scoring graph. None sets no cap.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored in each search over graphs. No move adds a forbidden edge or removes a required one, and each 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, the interventional essential graph of the learned targets: directed marks (->) are edges every graph in the class orients alike, among them every edge at a target, and undirected (--) the ones it leaves open. metadata["algorithm"] is "gnies"; metadata["score"] the score of the result, lower is better (see Notes); metadata["changing_variables"] lists the targets, known and learned, by label or 0-based position, in column order: the variables whose noise variance differs between some two environments. A shift in a variable’s mean is not among them, since each environment’s mean is taken out, and neither is a change made alike in every environment. With env, metadata["environments"] holds its labels in sorted order, so targets=dict.fromkeys(out.metadata["environments"], out.metadata["changing_variables"]) gives every environment the learned targets as known: a refit on it searches the graph alone and returns the same graph and score. Without env, targets gives every environment’s targets, which are then known.

Return type:

StructureOutput

Raises:

• TypeError – If 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 data is not a 2-D numeric matrix with at least two rows, no NaN or inf, and no constant column; if a column is declared categorical by a DataFrame’s category, bool, or text type; if the rows come from fewer than two environments, an environment has fewer than two rows, n_samples - n_environments is below the variable count, or a column is constant within one environment; if a column repeats another up to scale and shift or is a linear function of other columns, within the environments (one message names every group of copies and says to keep one column per group). Also if data holds each environment’s rows apart with env or with targets other than a mapping, or its tables differ in kind or columns, or one has no rows; if env does not hold one label per row, misses a label, or mixes strings and numbers; if a targets mapping comes without env, names a label env lacks, or names one label by two keys; if a column or an array of targets gives two rows of one environment different targets, or an array does not have data’s shape or holds a value other than 0 or 1; if targets names a variable data lacks; or if env or targets names a column data lacks. Also if a given max_parents is below the parents some graph of the required edges of knowledge gives a node; if the required edges give a known target more parents than the smallest environment’s rows minus 2; or 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 the fit reads as numbers. Each message warns once per process.

See also

andrey.gies

Learns from environments whose targets are known and replaced by interventions.

andrey.cdnod

Learns from environments without a model of how they differ.

Notes

When to use. GnIES suits independent rows from two or more environments, with or without an observational one, of one linear structural equation model with Gaussian noise, in which every environment keeps the same weights and differs only in the noise variances of some variables, its targets, which it learns without being told them [1]. It assumes no unmeasured common causes and a faithful distribution. It returns the learned targets and the interventional essential graph, which directs every edge at a target, so the environments can orient edges that a CPDAG leaves undirected. Its authors simulate graphs of 10 variables and present it as most useful with many environments of few rows each; measured with 100 or more rows per environment, it learned graphs of 10 to 100 variables with errors far below the empty graph’s, a fit taking about a second at 10 variables and three minutes at 100 (see Limitations). When an intervention replaces a variable’s mechanism and the targets are known, use gies(); when how the environments differ is not linear and Gaussian, cdnod(); and on observational data alone, ges().

How it works. Each environment’s columns are centered by their own means. For a set of targets, a variable outside it is scored by the Gaussian likelihood of its regression on its parents with one noise variance over every environment, and a target with one set of coefficients and a noise variance per environment, fitted by alternating weighted least squares from the regression over every environment. The search keeps one set of targets for all environments [1]. The equivalence class of a set of targets holds the graphs with the same adjacencies and v-structures that orient every edge at a target alike [4], and all of them get the same score [1]. For each set it visits, a greedy search over these classes [2], which extends greedy equivalence search [3], adds, removes, and turns edges from the empty graph while the score improves. The outer search starts from the known targets, adds the variable whose set scores best while that improves the score, then removes added variables the same way. Steps 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. With every environment’s targets known, one search over graphs runs, with the union of the targets.

The score. metadata["score"] is the sum of each variable’s term. Outside the targets it is n_samples * log(residual variance) + penalty_discount * log(n_samples) * n_parents; at a target, sum_e n_e * log(residual variance in e) + penalty_discount * log(n_samples) * (n_parents + n_environments - 1), with n_e the rows of environment e. Each term also adds 2 * ebic_gamma * log C(n_variables - 1, n_parents). Variances are divided by the row count.

Choosing the penalty. The paper’s score is the standard BIC, penalty_discount=1.0 with ebic_gamma=0.0 [1]. The default adds the extended-BIC term [5], as ges() and gies() do, which describe the choice; with 10 variables it adds 2 * log 9 = 4.4 to the first parent’s penalty, against log(n_samples), 8.3 with 4,000 rows. On the simulated data of the Limitations, paired per setting, the default made fewer errors than the standard BIC on 20 datasets in 7 of the 21 settings of 10 variables on the SHD and in 1 on the target error, and more in none; on 20 more datasets of those 7 settings, the SHD difference held in 3 and the target-error difference in none. With 50 and 100 variables its median SHD was 1 and 2.5, against 16.5 and 54.5 at the standard BIC.

Data scaling. Rescaling or shifting a column alike in every row does not change the result, nor does shifting one environment’s rows. Do not standardize each environment on its own: that rescales the noise variances the method reads, and changes the result.

Limitations. Measured with the targets unknown on simulated data from the model above: Erdos-Renyi graphs of 10 variables at mean degree 2.7 or 4, two to ten environments with 10 to 1,000 rows each, with or without an observational one, one or three targets per environment, 20 datasets per setting, one core per fit. Against the paper’s reference implementation, the 95% interval of the mean difference in errors lay above 0 in none of the 21 settings, on the SHD against the interventional essential graph or on the targets learned wrongly or missed, and below 0 on the SHD in 2 settings with 10 rows per environment and on the targets in 1 of them, on 20 more datasets too; at the standard BIC the two returned the same graph and targets on 412 of the 420 datasets. With 100 or more rows per environment its errors were below the empty graph’s on both counts in every setting. With 10 rows per environment its target errors were 1.85 to 3 with 3 true targets and 4.35 and 7.7 with 9, and its SHD was not reliably below the empty graph’s at mean degree 4 (20 against 20) or with the targets’ noise variances set to 3 to 4 instead of 5 to 10 (11.9 against 14; the 95% interval of the difference reached above 0 on both sets of 20 datasets). Outside its model, with 1,000 rows per environment: after hard interventions, which replace a target’s mechanism, its mean SHD was 3.75 against the empty graph’s 14; a shifted mean alone is invisible to it, so it learned almost no targets (3.1 target errors with 3 true targets) and returned close to the observational class (SHD 2.55 against it); after a target’s incoming weights were halved it beat the empty graph on neither count (SHD 11.6 against 14, 3.05 target errors with 3 true targets); with Laplace noise its SHD was 1.65, and with additive nonlinear relations 7.55. With 1,000 rows in each of four environments, three of them with one target, its median SHD at 20, 30, 50, and 100 variables was 0, 0, 1, and 2.5, against the empty graph’s 27, 40, 68, and 135, and its median fit took 3.5, 9.0, 30, and 176 s. At 10 variables its median fit took 0.3 to 1.3 s per setting, a median 0.03 times the paper’s reference implementation’s time on the same core; at 20 variables the reference implementation took a median 388 s, and at 30 it finished 2 of 3 fits, in 2,735 and 3,327 s, within a cap of 3,600 s. Each set of targets it visits costs one search over graphs from the empty graph: each forward step scores every variable not yet a target, and each backward step every target, about (n_targets + 1) * n_variables - n_targets * (n_targets - 1) / 2 searches per fit, with n_targets the targets it returns. Given every environment’s targets, one search over graphs runs, and it can stop at a graph that scores worse than one the target search reaches: with the true targets given, its mean SHD was higher than when it learned them in 17 of the 21 settings. A target’s coefficients are shared across environments and its variances are not, so a target with k parents needs at least k + 2 rows in every environment: no parent set larger than the smallest environment’s row count minus 2 is selected for a target, and a variable that the required edges of knowledge give more parents than that is not learned as one. The numbers are a reference, not a hard limit: they change with the data generator, the density, the method’s settings, and the machine.

References

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

Gamella, J. L., Taeb, A., Heinze-Deml, C., and Buhlmann, P. (2022). Characterization and greedy learning of Gaussian structural causal models under unknown interventions. arXiv:2211.14897.

[2]

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

[3]

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

[4]

Tian, J., and Pearl, J. (2001). Causal discovery from changes. Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, 512-521.

[5]

Chen, J., and Chen, Z. (2008). Extended Bayesian information criteria for model selection with large model spaces. Biometrika 95(3), 759-771.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from andrey import gnies
>>> rng = np.random.default_rng(0)
>>> env = np.repeat(["observed", "noisy y"], 500)
>>> x = rng.standard_normal(1000)
>>> y = 0.8 * x + rng.standard_normal(1000) * np.where(env == "noisy y", 3.0, 1.0)
>>> out = gnies(pd.DataFrame({"x": x, "y": y}), env)
>>> out.structure.oriented_edges()
[('x', 'y', 'directed')]
>>> out.metadata["changing_variables"]
['y']