API reference

andrey.pcmciplus#

andrey.pcmciplus(data, *, tau_max=1, pc_alpha=0.01, collider_rule='majority', indep_test='student_t')[source]#

Learn the lagged and same-step causal links of one multivariate time series with PCMCI+.

PCMCI+ [1] extends PCMCI to links within one time step. It first selects, for each variable, a few candidate parents among the lagged values of all variables. It then tests every lagged candidate and every pair of variables at the same step, given subsets of the same-step neighbors plus the selected lagged parents of both ends, and orients the same-step links with the collider and orientation rules of the PC algorithm. Lagged links always point forward in time. Use it on one series whose causal effects can act within a 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 0 to tau_max are tested.

  • pc_alpha (float) – Significance level of every test, in (0, 1): a link is removed when a test gives it a p-value above pc_alpha.

  • collider_rule (str) – How an unshielded triple a - b - c with a same-step link b - c is judged a collider a -> b <- c. "majority" and "conservative" test a and c again given every subset of their same-step neighbors and look at the subsets that separate them: "majority" calls the triple a collider when fewer than half hold b, and "conservative" when none does; a triple neither rule can decide is ambiguous and orients nothing. "none" calls it a collider when b is not in the one separating set the skeleton search found.

  • 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 0 to tau_max. structure.lag(0) is a cpdag of the same-step links: i -> j is oriented, and i -- j is a link whose direction the data cannot tell or whose orientations conflict; metadata["conflicts"] lists the latter as [i, j] pairs, i < j. structure.lag(k), k >= 1, 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][i][j] is the largest p-value among the skeleton-phase tests of the link between variable i at t - k and variable j at t (the first step’s tests for a lagged link that step removed), and metadata["statistics"][k][i][j] the partial correlation of that test; lag 0 is symmetric with 1 and 0 on the diagonal. metadata["parents"][j] lists the [i, k] pairs (variable i at lag k) the first step selected for variable j, strongest first, and metadata["ambiguous_triples"] the [i, k, m, j] triples (variable i at lag k, then m and j at lag 0) the collider rule left undecided. ordering is None, and metadata["algorithm"] is "pcmciplus".

Return type:

StructureOutput

Raises:

  • TypeError – If tau_max is not an int.

  • ValueError – If tau_max is below 1, if pc_alpha is not in (0, 1), if collider_rule or indep_test is not one of its choices, or if data is not a usable 2-D matrix with at least 2 * tau_max * (n_variables + 1) + n_variables + 1 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

The method runs in four phases. The lagged phase selects the lagged parents of each variable as PCMCI’s PC1 step does, with equal strengths ordered by each variable’s rank. The skeleton phase starts from those lagged links and from a link between every two variables at the same step. Round p = 0, 1, ... tests each linked pair given each subset of p same-step neighbors of the later end, with the lagged parents of both ends added, and removes the link at the first p-value above pc_alpha. A same-step pair is tested from both ends, first given the neighbors of its higher-numbered variable, and the neighbor sets stay fixed during a round, so the links found do not depend on the order of the variables. The collider phase orients the unshielded triples by collider_rule, and the rule phase applies Meek’s rules R1 to R3 to the same-step links until none applies. Each collider or rule pass collects its orientations before applying them, and marks a link that gets both orientations as a conflict. With the "majority" or "conservative" rule the orientations do not depend on the order of the variables either; "none" reads the one separating set the skeleton phase stored, which can change with the order.

All tests use the same time steps, the last n_samples - 2 * tau_max. The default test, "student_t", takes the p-value from the Student-t distribution with n - |S| - 2 degrees of freedom, for |S| conditions. 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.

Limitations. Measured on simulated linear-Gaussian series with strong autocorrelation and 1.5 cross-links per variable, 30% of them within a time step and the rest at lags 1 to 5, from 10 to 80 variables with 500 to 2,000 time steps, with tau_max=5. Its median errors (endpoint SHD) stayed well below the empty graph’s at every size. It assumes linear effects, a stationary series, and no hidden common causes. The skeleton and collider phases test every subset of a variable’s same-step neighbors, so the number of tests grows exponentially with the number of same-step links at one variable. A test conditions on up to 2 * n_variables * tau_max + n_variables - 2 values, so the series needs many more time steps than that.

References

Examples

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