andrey.jpcmciplus

andrey.jpcmciplus(data, *, time_contexts=None, space_contexts=None, time_dummy=None, tau_max=1, alpha=0.01, fdr_method='fdr_bh', collider_rule='majority', indep_test=None, seed=None, max_depth=None, max_p_global=None, knowledge=None)

Learn the links of time series from several datasets with shared contexts with J-PCMCI+.

J-PCMCI+ [1] extends PCMCI+ [2] to datasets of the same variables, such as sites, subjects, or regimes, whose differences come from context variables. A time context varies over time and is the same in every dataset at each step; a space context is fixed within a dataset. A context can be measured or not: a one-hot dummy for the dataset and one for the time step stand in for the contexts nobody measured, so that a hidden context does not link the variables it drives. The method finds the links from the contexts into the system variables, then the links between the system variables, and returns the graph over the system variables.

Parameters:

• data (Sequence[ArrayLike]) – Two or more datasets over the same columns, one row per time step; their lengths may differ. Each is coerced to float64. DataFrames must name their columns alike; the names of the system variables become the structure’s labels. The columns that time_contexts and space_contexts do not list are the system variables.

• time_contexts (Sequence[str | int] | None) – The columns that hold time contexts: column labels (strings) or 0-based positions. Each must hold the same values in every dataset at each step, counted from each dataset’s first row. None lists none.

• space_contexts (Sequence[str | int] | None) – The columns that hold space contexts, listed as in time_contexts. Each must be constant within a dataset. None lists none.

• time_dummy (bool | None) – Whether the fit adds the dummy of the time step, which stands in for hidden time contexts. It needs datasets of one length, aligned in time, with more than one step after the first 2 * tau_max. None, the default, adds it when the datasets allow it; True requires it and raises ValueError when they do not; False leaves it out, for datasets that are not aligned in time or to save its one column per step. The dummy of the dataset is always used. metadata["time_dummy"] says whether the fit used it.

• tau_max (int) – The largest lag tested, at least 1: links at lags 0 to tau_max are tested.

• alpha (float) – Significance level of every test, in (0, 1): a link is removed when a test gives it a p-value above alpha. The default is pcmciplus()’s; the published experiments of the method use 0.05 [1].

• fdr_method (str) – How the p-values of the lagged links between system variables are corrected for the number of links tested. "fdr_bh" replaces them with their Benjamini-Hochberg adjusted values [3] over all those links together, as in pcmciplus(), and "none" keeps them as tested. The links do not change: the skeleton passes decide them at alpha, so fdr_method only changes metadata["p_values"].

• collider_rule (str) – How an unshielded triple a - b - c with a same-step collider candidate b - c between system variables is judged, as in pcmciplus(); a can be a context or a dummy. The default, "majority", gives orientations that do not depend on the order of the variables.

• indep_test (str | None) – The conditional-independence test behind every p-value, a test of continuous data, as in pcmciplus(): "fisherz", "student_t", and "normal_scores" test partial correlation, and "kci", "rcot", "cmiknn", and "gpdc" detect nonlinear dependence at a far higher cost per test. The dummies enter a test as numeric columns, so a categorical or mixed-data test raises ValueError. None picks "fisherz". metadata["indep_test"] names the test the fit used.

• 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 of the lagged phase (PC1 as in pcmci()), at least 0; None sets no limit. It does not cap the later passes, which add the parents of both ends to every conditioning set; max_p_global caps those.

• max_p_global (int | None) – The most conditions a test of the skeleton passes and the collider phase takes beyond the default conditions, at least 0: every such test conditions on the parents of both ends, the default conditions, and on a subset of at most max_p_global same-step neighbors; the cap also bounds the subsets the collider vote reads. None sets no limit.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge, over the system variables, read as links between lagged variables as in pcmciplus(): required and forbidden links, the pairs allowed to be linked, and temporal tiers, all honored. A position counts the system variables. The contexts bring their own knowledge: no system variable causes a context, and no two contexts are linked.

Returns:

Wraps a temporal TemporalStructure over the system variables, as pcmciplus() returns it: structure.lag(0) is a cpdag and structure.lag(k), k >= 1, a digraph. metadata["p_values"], metadata["statistics"], metadata["parents"], metadata["conflicts"], and metadata["ambiguous_triples"] read as there, over the system variables by position. metadata["context_links"] lists the links from contexts into system variables as [context, lag, variable], a space context at lag 0, with a context by name or position in data and a variable by name or position among the system variables. metadata["index_edges"] lists the links from the dummy of the dataset as [variable, "directed"] pairs: the variables that change with the dataset in a way the measured contexts do not explain. metadata["time_index_edges"] does the same for the dummy of the time step, and metadata["time_dummy"] says whether the fit used it. A dummy link is not a causal link. There is no metadata["changing_variables"]: the dummy links stand in for the hidden contexts, and the links of the measured ones are in metadata["context_links"]. ordering is None, and metadata["algorithm"] is "jpcmciplus". With a randomized indep_test, metadata["seed"] is the seed used. metadata["tie_break"] is the tie-break rule the fit ran under, andrey.config.tie_break.

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 time_contexts or space_contexts 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. Also if data holds fewer than two datasets, a dataset is not a 2-D table, the datasets differ in their columns, or a dataset has at most 2 * tau_max rows; if time_dummy=True and the datasets differ in length or have one step after the first 2 * tau_max; if a test conditions on so many columns of lagged values, contexts, and dummy levels that the rows leave it no degree of freedom; if a column is both a time and a space context or a context names a column data lacks, or fewer than two system variables remain; if a time context differs between datasets at some step, a space context varies within a dataset, or a space context repeats another up to scale and shift, as any two do over two datasets; or if data holds NaN. With "fisherz", "student_t", or "normal_scores", also if a variable repeats another up to scale and shift, at the same step or up to 2 * tau_max steps apart, or a test finds a system variable in the span of its conditions, as with one that sums others. Also if knowledge names a column that is not a system variable, such as a context, contradicts itself, requires a link without a lag or beyond tau_max, or directs same-step links into a directed cycle over the links the fit kept.

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.

Notes

When to use. J-PCMCI+ suits several datasets of the same variables, each a stationary time series, that differ through context variables: no system variable causes a context, no hidden variable causes two system variables except a context, and no hidden context stands between a measured context and a system variable or causes both [1]. Every context must be a function of the dataset or of the time step that does not tell them all apart. The method suits few system variables: its context and dummy passes keep every same-step link between system variables, so their tests grow exponentially with those variables, and the dummies add a column per dataset and per step. It returns the temporal structure of pcmciplus() over the system variables and lists the context and dummy links in metadata. For one series, use pcmciplus(), or lpcmci() when hidden common causes are not contexts; for i.i.d. rows from several environments, cdnod() or fci_jci().

The method runs PCMCI+’s lagged phase on the system variables and time contexts, then the skeleton in three passes: the links from the contexts into the system, given the lagged parents; the links from the dummies into the system, given those and the context parents; and the links between system variables, given those and the dummy parents. Links from contexts and dummies point into the system, and PCMCI+’s collider and rule phases use them to orient the same-step links between system variables.

Every test reads the datasets’ lagged rows stacked, so a test reads M * (T - 2 * tau_max) rows for M datasets of T steps, and no lag window crosses from one dataset into the next. A dummy has one level per dataset or per step. A test of a dummy tests each level against the rest and multiplies the smallest p-value by the number of levels, Bonferroni’s correction, so a dummy link needs one level that stands out; two levels are one test. Conditioning on a dummy takes all its levels but one, which centers each variable within each dataset or step, and a context is fixed given its dummy, so it leaves the conditions.

Missing values. A NaN in data raises ValueError: the stacked design reads every step. Fill in the missing steps first, as dropping one would join the steps around it, or use pcmci(), pcmciplus(), or lpcmci(), whose tests delete rows per test.

Limitations. Measured on the linear model of its authors’ experiment script [1], with both dummies at tau_max=2 and alpha=0.05, on one core: 5 system variables with lagged links only, 2 time and 2 space contexts, every context observed or one of each kind, and 10 datasets of 10 to 200 steps or 5 to 200 datasets of 10 steps. Its mean errors (endpoint SHD over the links between system variables at lags 0 to 2) were below the empty graph’s 20 in every setting: 3.1 to 7.7 with every context observed; with one of each kind observed, where no hidden context confounds an observed one and a system variable, 3.6 to 7.2 with 10 to 200 datasets and 13.0 with 5 datasets of 10 steps. The median fit took 0.017 to 0.17 s. With 10 and 15 system variables and 10 datasets of 100 steps, its mean errors were 8.1 and 18.8 against the empty graph’s 40 and 60. With "gpdc" on the authors’ nonlinear model of 2 system variables and 20 steps, the median fit took 4.8 s with 5 datasets and 385 s with 50. It assumes faithfulness, so a dummy link it misses can leave the links a hidden context adds. A dataset effect far larger than the noise can hide its dummy, as the lagged values stand in for it: with 10 datasets of 100 steps, the dummy link was found for effects of 1 and 2 noise standard deviations and missed for 5, which left false lagged links. The lagged phase tests each candidate without conditions first, so an autocorrelated target can hide a context link from it for good: on data from sample_context_series() with 10 datasets, at alpha=0.01 and tau_max=2, it found 82% of the context links at 50 steps and 92% at 100. With few steps against many datasets, false positives stay high as the datasets grow, a finite-sample bias of the kind fixed-effects panels have [1]. A time context that the lags of the time contexts determine, such as a sine and cosine of one period, or a sine or a quadratic trend alone from tau_max=2, cannot be told apart from them: each of its links a test cannot read stays, so context_links can name other lags of the context that has the effect, and other time contexts that determine it at some lag, including ones with no effect. Only the variables the context links reach can be trusted, not the context or the lag that reaches them. Following the text of [1], it tests no same-step link between two contexts and reports no link between contexts: the lagged parents of each time context, which the lagged phase selects among the time contexts, serve only as conditions. A dummy link needs one level that stands out, so the time dummy misses a smooth effect spread over the steps, and the links it adds stay: list such a context if it is known. Under the "majority" and "conservative" rules, a context cannot orient a link into a variable that has the context’s dummy as a parent, so hidden contexts leave same-step links undirected: on 10 datasets of 100 steps, with one measured context and a hidden one on every variable, of 300 links 5 were oriented correctly, 5 wrongly, and 290 not at all. Cost grows fastest with the system variables: on one core and 10 datasets of 100 steps, the median fit with both dummies took 0.3 s at 10 variables and 11.5 s at 15 in the measurement above; without the time dummy, a fit took 0.1 s at 10, 1.5 s at 15, and 47 s at 20 and ran past 300 s at 30, and with it, a fit ran out of 12 GB at 20. The time dummy adds a column per step to the conditions of every test that takes it: at 5 variables and 1,600 steps a fit took 7 s with it and 0.05 s without. The nonlinear tests read the dummies’ columns as any other. 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,6)

Gunther, W., Ninad, U., and Runge, J. (2023). Causal discovery for time series from multiple datasets with latent contexts. Proceedings of the 39th Conference on Uncertainty in Artificial Intelligence (UAI), PMLR 216, 766-776.

[2]

Runge, J. (2020). Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), PMLR 124, 1388-1397.

[3]

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57(1), 289-300.

Examples

>>> import numpy as np
>>> from andrey import jpcmciplus
>>> rng = np.random.default_rng(0)
>>> def dataset(level):
...     x = level + rng.standard_normal(200)
...     y = np.r_[0.0, 0.8 * x[:-1]] + rng.standard_normal(200)
...     return np.column_stack([x, y])
>>> out = jpcmciplus([dataset(level) for level in (0, 3, 6, 9)], time_dummy=False)
>>> out.structure.lag(1).oriented_edges()
[(0, 1, 'directed')]
>>> [variable for variable, _ in out.metadata["index_edges"]]
[0]