andrey.loci

andrey.loci(data, *, alpha=None, criterion='independence', missing=None)

Test the direction of a variable pair with LOCI, the location-scale noise model.

LOCI [1] fits y = f(x) + g(x) n in each direction, a noise scale g that changes with the cause, by a Gaussian likelihood over B-spline features of the cause. Either an HSIC test [2] of the scaled residual against the proposed cause or the two likelihoods determine the edge.

Parameters:

• data (ArrayLike) – The two variables, one per column: column 0 is x and column 1 is y. Coerced to float64. A DataFrame’s column names become the structure’s labels.

• alpha (float | None) – Under criterion="independence": None chooses the direction with the larger p-value, or with the smaller test statistic when the p-values are equal; a significance level such as 0.05 draws an edge only when exactly one p-value exceeds it. Must be None under criterion="likelihood". See Notes.

• criterion (str) – "independence" tests whether each direction’s scaled residual is independent of its proposed cause; "likelihood" chooses the direction whose model has the larger likelihood. See Notes.

• missing (str | None) – How missing values (NaN) are handled. "testwise" deletes rows per test: each test reads only the rows where all of its variables are observed. None picks the method’s default, "testwise". On data with missing values, metadata["missing"] records the handling used. See Missing values. in Notes.

Returns:

Contains a two-node GraphStructure of kind dag with the chosen edge, or no edge. metadata["algorithm"] is "loci". metadata["log_likelihoods"] holds the mean log-likelihood of y given x under "forward" and of x given y under "backward". Under criterion="independence", metadata["p_values"] holds the independence p-value for x -> y under "forward" and for y -> x under "backward", and metadata["statistics"] the HSIC statistics of the same two tests, each divided by its mean under independence, under the same keys.

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.

• 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 does not have 2 columns, or if alpha is not None under criterion="likelihood".

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. LOCI decides the direction between two continuous variables observed i.i.d., under the location-scale noise model y = f(x) + g(x) n with g > 0, the noise n independent of the cause x, and no hidden common cause of the two [1]. The direction is identifiable, outside pathological cases, for noise of any distribution with a twice-differentiable log density and twice-differentiable f and g; the cause must have unbounded support, which excludes a uniform cause [1]. It is designed for one pair whose effect’s spread changes with the cause, which the additive noise model of anm() does not allow, and it returns a two-node DAG whose edge depends on criterion. For more than two variables, use cam(), whose Limitations give its errors under such noise. Keep the default criterion="independence" unless the noise is known to be Gaussian: it named the true direction on most pairs under every noise distribution measured, while criterion="likelihood" was ahead of it under Gaussian noise and below one pair in two under Gumbel noise (Limitations).

With criterion="independence", LOCI-H in [1], each direction’s scaled residual (y - mu(x)) / sigma(x) is tested against the proposed cause by HSIC with Gaussian kernels, each bandwidth the median distance between the rows, and a gamma approximation of the null distribution [2]. The p-values then decide as for anm(): with alpha=None, the larger p-value, compared through their logarithms, and on equal p-values the smaller statistic, with no edge when those are equal too; with a significance level, an edge only when exactly one direction’s p-value exceeds it. Mooij et al. (2016) [3] compare p-values to avoid choosing a level. The paper proposes this criterion for noise that is not Gaussian.

With criterion="likelihood", LOCI-M in [1], the edge points in the direction whose model has the larger mean log-likelihood, and equal ones give no edge. The likelihood identifies the direction consistently when the noise is Gaussian, as the fitted model assumes [1]. In the paper’s benchmarks it is ahead of the independence criterion on location-scale data and behind it on data outside that model.

Each column is standardized first, as in [1], so the result does not depend on the units of either variable; the Gaussian marginals of the two standardized causes are then equal, so the likelihoods compare the two conditional models. In each direction the effect given the cause is Gaussian with natural parameters eta_1 = Phi w_1 and eta_2 = -Phi w_2, where Phi is a B-spline basis of the cause with 25 knots, of degree 5, and w_2 >= 0, so the log-likelihood is jointly concave [1]. A Gaussian prior of precision 1e-6 on the weights makes it strictly concave [1]. Projected Newton steps [4] maximize it until the mean objective changes by less than 1e-6 in a full step, or for at most 100 steps; the paper takes a closed-form step in w_1 and L-BFGS steps in w_2 in turn, which approach the same maximum in many more steps. The fit draws no random numbers, so it is deterministic.

Missing values. A NaN in data marks a missing value. Both fits read only the rows where both columns are observed, which for a pair are the rows test-wise deletion keeps; metadata["rows_tested"] gives their count as both min and median. When values are missing completely at random, those rows are a random sample of all rows. Fewer than 2 such rows, or a column constant on them, raises ValueError.

Limitations. Measured on simulated cause-effect pairs of 500, 1,000, and 2,000 rows, 100 pairs per setting, on one core: location-scale pairs with Gaussian, Gumbel, or Laplace noise, additive-noise pairs with Gaussian noise, and, at 1,000 rows, the five synthetic sets of Tagasovska et al. (2020) [5]. With criterion="independence" it named the true direction on 79% to 86% of the pairs of each setting at 500 rows, 90% to 100% at 1,000, and 96% to 99% at 2,000, and the median fit took 0.07, 0.2, and 0.8 s. With criterion="likelihood" it named it on every location-scale pair with Gaussian noise at each size and on 85%, 93%, and 99% of the additive-noise pairs, in a median of at most 0.1 s. The model’s noise is Gaussian, so under other noise the likelihood criterion is misspecified: it named the true direction on 35% to 38% of the pairs under Gumbel noise at every size, fewer than a coin flip, and under Laplace noise on 37% at 500 rows, 44% to 47% at 1,000, and 57% to 59% at 2,000. On 20 location-scale pairs with Gaussian noise of more rows, with a cap of 1 h and 28 GB per fit, both criteria named the true direction on every pair they fitted; under criterion="independence" the median fit took 3.9 s and at most 0.8 GB at 4,000 rows, 15 s and 2.7 GB at 8,000, and 60 s and 10 GB at 16,000, and every fit at 32,000 rows ran out of the 28 GB, while under criterion="likelihood" it took 1.4 s and 0.2 GB at 32,000. Each HSIC test builds matrices with one row and one column per sample, so under criterion="independence" time and memory grow with the square of n_samples; the likelihood fit grows linearly with it. It decides one pair: the paper defines no procedure for more variables. A spline at either end of the cause’s range covers few rows, where the fitted variance can shrink toward zero as the likelihood grows, which only the prior bounds. Missing values are not benchmarked yet; the tests run them on pairs with values missing completely at random. 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,7,8,9)

Immer, A., Schultheiss, C., Vogt, J. E., Scholkopf, B., Buhlmann, P., and Marx, A. (2023). On the identifiability and estimation of causal location-scale noise models. ICML, PMLR 202.

[2] (1,2)

Gretton, A., Fukumizu, K., Teo, C. H., Song, L., Scholkopf, B., and Smola, A. J. (2008). A kernel statistical test of independence. NeurIPS 20, 585-592.

[3]

Mooij, J. M., Peters, J., Janzing, D., Zscheischler, J., and Scholkopf, B. (2016). Distinguishing cause from effect using observational data: methods and benchmarks. JMLR 17(32), 1-102.

[4]

Bertsekas, D. P. (1982). Projected Newton methods for optimization problems with simple constraints. SIAM Journal on Control and Optimization 20(2), 221-246.

[5]

Tagasovska, N., Chavez-Demoulin, V., and Vatter, T. (2020). Distinguishing cause from effect using quantiles: bivariate quantile causal discovery. ICML, PMLR 119, 9311-9323.

Examples

>>> import numpy as np
>>> from andrey import loci
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = np.sin(2 * x) + x + (0.6 + 0.5 * np.tanh(x)) * rng.standard_normal(500)
>>> out = loci(np.column_stack([x, y]))
>>> out.structure.oriented_edges()
[(0, 1, 'directed')]
>>> out = loci(np.column_stack([x, y]), criterion="likelihood")
>>> out.metadata["log_likelihoods"]["forward"] > out.metadata["log_likelihoods"]["backward"]
True