API reference

andrey.pcmci#

andrey.pcmci(data, *, tau_max=1, pc_alpha=0.2, alpha_level=0.05, indep_test='student_t')[source]#

Learn the lagged causal links of one multivariate time series with PCMCI.

PCMCI [1] works in two steps. The PC1 step selects, for each variable, a few candidate parents among the lagged values of all variables. The momentary conditional independence (MCI) step then tests each lagged link given the selected parents of both its ends. Conditioning on the cause’s own past removes the autocorrelation that makes a plain lagged-correlation test report false links. Use it on one series whose causal effects take at least one time step.

Parameters:
  • data (ArrayLike) – A multivariate time series, one row per time step. Coerced to float64. A DataFrame’s column names become the structure’s labels.

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

  • pc_alpha (float) – Significance level of the PC1 step, in (0, 1): a candidate parent is dropped when a test gives it a p-value above pc_alpha. A larger value keeps more parents as conditions for the MCI step.

  • alpha_level (float) – Significance level of the MCI step, in (0, 1): a link is in the graph when its p-value is at most alpha_level.

  • indep_test (str) – The partial-correlation test behind every p-value. Both report the same partial correlation; "student_t" takes its p-value from the Student-t distribution and "fisherz" from Fisher’s z-transform.

Returns:

Wraps a temporal TemporalStructure with lags 1 to tau_max and no lag 0. structure.lag(k) is a digraph whose edge i -> j is the link from variable i at time t - k to variable j at time t; a self-loop is a variable’s dependence on its own past. The structure has no lag_weights. metadata["p_values"][k - 1][i][j] and metadata["statistics"][k - 1][i][j] hold the MCI p-value and partial correlation of that link, for every pair, significant or not. metadata["parents"][j] lists the [i, k] pairs (variable i at lag k) that PC1 selected for variable j, strongest first. ordering is None, and metadata["algorithm"] is "pcmci".

Return type:

StructureOutput

Raises:
  • TypeError – If tau_max is not an int.

  • ValueError – If tau_max is below 1, if pc_alpha or alpha_level is not in (0, 1), if indep_test is neither "student_t" nor "fisherz", or if data is not a usable 2-D matrix with at least 2 * tau_max * (n_variables + 1) + 3 rows, the length at which every test keeps a degree of freedom. Also if a variable repeats another up to scale and shift, at the same step or up to 2 * tau_max steps apart, before any test runs; one message names every group of copies, a lagged value as 'x' at lag k, and says to keep one column per group.

Notes

Every test is a partial-correlation test. The default, "student_t", takes the p-value from the Student-t distribution with n - |S| - 2 degrees of freedom, for |S| conditions, as the authors’ reference implementation does. "fisherz" uses Fisher’s z-transform instead; its p-values differ slightly, so a link close to alpha_level can flip. The partial correlation is the correlation of the least-squares residuals of the two tested values given the conditions. It is computed from the correlation matrix, and from the residuals themselves when the conditions are collinear to working precision, as lagged copies of a strongly autocorrelated series can be.

All tests use the same time steps, the last n_samples - 2 * tau_max, where every lag a test can condition on exists: the MCI test of a link at lag k conditions on the cause’s parents, which reach back k + tau_max steps. PC1 runs in rounds p = 0, 1, ...: each round tests every remaining candidate given the p strongest other candidates, one conditioning set per candidate, then drops those that test independent. A candidate’s strength is the smallest absolute partial correlation it has had so far. Equal strengths are ordered by each variable’s rank, as in hc(); the column-order guide says what the rank keeps.

Limitations. Measured on simulated linear-Gaussian series with an autodependency in every variable and one cross-link per variable at lag 1 or 2, from 10 to 160 variables with 500 to 2,000 time steps, with tau_max=2. It found nearly every true link at every size, but its false links grew with the number of variables: its median errors (lag-resolved SHD) stayed below the empty graph’s at 10 variables, came close to it at 20, and exceeded it from 40 on. It assumes linear effects, a stationary series, and no hidden common causes. It does not test same-step (lag 0) links: an effect within one time step is missed and can show up as a false lagged link. An MCI test conditions on up to 2 * n_variables * tau_max - 1 lagged values, so the series needs many more time steps than that. The MCI step alone runs n_variables ** 2 * tau_max tests.

References

Examples

>>> import numpy as np
>>> from andrey import pcmci
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = np.r_[0.0, 0.8 * x[:-1]] + rng.standard_normal(500)
>>> out = pcmci(np.column_stack([x, y]), tau_max=2, alpha_level=0.01)
>>> out.structure.lag(1).oriented_edges()
[(0, 1, 'directed')]
>>> out.structure.lag(2).oriented_edges()
[]