Calibration#
TL;DR#
Skewed columns: the new
normal_scorestest reads only ranks, so it keeps its 5% level where Fisher’s z rejects a true independence up to 56% of the time.Sparse tables: the G-square test’s false rejections fall from 21.7% to 3.8%.
Many lags at once: PCMCI and Granger control the false discovery rate by default. At 80 variables, PCMCI’s extra links fall from 616 to 13.
Check a graph:
andrey.markov_checktests the independences a graph implies against data it was not learned from.How stable: the bootstrap refits on half-size subsamples, and its averaged graph is a valid graph of the method’s kind, or a warning says why it cannot be.
Three more supported methods: Inter-IAMB, MMHC, and RSMAX2, none worse than bnlearn in any of 72 comparisons.
In 0.12.0, we make Andrey’s results something you can check.
A causal search is a long chain of statistical tests. Each test at level 0.05 promises to call a true independence a dependence only 5% of the time. When that promise slips, the search draws edges that are not there, and nothing in the graph tells you so.
We found three places where the promise slipped, and fixed them. Then we added a way to test the result itself against the data.
A skewed-data example#
We simulate ten variables on a known graph and take the exponential of each column, so every column is log-normal, as many lab values are. We learn on one half of the 4,000 samples and check on the other. The whole script is below; each section that follows runs one part of it.
import numpy as np
import pandas as pd
import andrey
from andrey.core.structure import GraphStructure
from andrey.data import sample_scm
from andrey.metrics import shd
# Ten variables on a known graph, each column the exponential of a Gaussian one: log-normal, as
# many lab values are.
truth = sample_scm(d=10, n=4000, seed=2, density=3.0)
names = [f"x{i}" for i in range(10)]
data = np.exp(pd.DataFrame(truth.data, columns=names))
graph = GraphStructure.from_edges(truth.graph.to_edges(), n_nodes=10, kind="dag", labels=names)
learn, held_out = data.iloc[:2000], data.iloc[2000:] # learn on one half, check on the other
if __name__ == "__main__": # worker processes re-import this script
# The test: Fisher's z assumes Gaussian columns; normal_scores reads only their ranks.
for test in ["fisherz", "normal_scores"]:
fit = andrey.pc(learn, indep_test=test)
print(f"{test}: {shd(fit.structure, graph)} errors against the true graph")
# The check: does the held-out half contradict a graph?
for test in ["fisherz", "normal_scores"]:
check = andrey.markov_check(graph, held_out, indep_test=test)
print(
f"true graph, checked with {test}: p = {check.uniformity_p_value:.3f}, "
f"{'passed' if check.passed else 'rejected'}"
)
edges = graph.to_edges()
missing_one = GraphStructure.from_edges(edges[1:], n_nodes=10, kind="dag", labels=names)
check = andrey.markov_check(missing_one, held_out, indep_test="normal_scores")
print(
f"one edge missing: p = {check.uniformity_p_value:.3f}, "
f"{'passed' if check.passed else 'rejected'}"
)
# How stable: 100 fits on half-size subsamples, four at a time, and their averaged graph.
with andrey.config(n_jobs=4):
boot = andrey.bootstrap("pc", learn, n_resamples=100, seed=0, indep_test="normal_scores")
pairs = {frozenset(e[:2]) for e in graph.oriented_edges()}
stable = boot.table(min_frequency=0.9)
true = sum(frozenset((e["source"], e["target"])) in pairs for e in stable)
print(f"edges in at least 9 fits in 10: {len(stable)}, {true} of them in the true graph")
averaged = boot.structure
print(f"averaged graph: {boot.metadata['averaged_graph']}, a {averaged.kind}")
# Ask the graph: what screens x0 off from the rest?
print(f"Markov blanket of x0: {averaged.markov_blanket('x0')}")
fisherz: 11 errors against the true graph
normal_scores: 5 errors against the true graph
true graph, checked with fisherz: p = 0.000, rejected
true graph, checked with normal_scores: p = 0.288, passed
one edge missing: p = 0.000, rejected
edges in at least 9 fits in 10: 7, 7 of them in the true graph
averaged graph: vote, a cpdag
Markov blanket of x0: ['x2', 'x6']
Rank-based test for skewed columns#
Fisher’s z test assumes Gaussian columns. On skewed ones, it finds dependences that are not there.
In our simulations, with two conditioning variables and columns passed through exp, a cube, or
a heavy-tailed map, it rejected a true independence in 24% to 31% of datasets at 200 samples and
52% to 56% at 2,000. More samples made it worse.
indep_test="normal_scores" replaces each column with the normal scores of its ranks, then runs
Fisher’s z. Any strictly increasing change to a column leaves its ranks, and so the p-value, the
same. In the same simulations it rejected 4% to 6% of the true independences, and it is as fast as
Fisher’s z. In the script, PC makes 11 errors with Fisher’s z and 5 with normal scores.
Every search with a conditional-independence test takes it, including PCMCI and its relatives for time series.
Calibrated tests for sparse tables#
The G-square and chi-square tests for categorical data split the table by every combination of the conditioning variables. With several of them, most cells are nearly empty, and the textbook degrees of freedom overstate what the table holds. Over 3 to 5 categories, 0 to 3 conditioning variables, and 500 to 5,000 samples, the G-square test rejected 21.7% of true independences at level 0.05, and the chi-square test 11.1%. With 5 categories and 3 conditioning variables, it was 83% and 48%.
Both tests now take the degrees of freedom from the statistic’s exact mean when one variable is shuffled within each configuration (Tsamardinos and Borboudakis, 2010). The rates are 3.8% and 4.9% overall, and 0.8% and 3.2% in the hardest cell. PC with the G-square test, the default for categorical data, recovered graphs nearer the truth than with the chi-square test in 3 of the 16 settings we measured, and no farther in the other 13.
False discovery rate control#
PCMCI tests every variable at every lag. At 80 variables and two lags, that is thousands of tests,
and 5% of thousands is hundreds of false links. PCMCI, PCMCI+, and Granger now take fdr_method,
and its default is "fdr_bh": PCMCI and Granger keep a link when its Benjamini-Hochberg adjusted
p-value is at most alpha, and PCMCI+ reports the adjusted p-values.
On simulated series of 80 variables with 1,000 steps, PCMCI’s errors fall from 779 to 252, against
481 for the empty graph: the extra links fall from 616 to 13, at the cost of more missed links.
Granger now decides each link by an F-test at alpha, where it kept the nonzero coefficients of a
lasso, and its median errors on our quality checks fall from 13.5 to 6.5.
Markov checker#
andrey.markov_check asks whether the data contradict a graph. A graph implies independences:
each variable is independent of the earlier variables that are not its parents, given its parents.
The check tests each one and asks whether the p-values look uniform, as they do when every
independence holds (Ramsey, Andrews, and Spirtes, 2024).
Check on samples the graph was not learned from: a graph fitted to the same samples can pass because the fit made its independences true there. In the script, the true graph passes with normal scores and fails with Fisher’s z, and a graph missing one edge fails. A check is only as good as its test.
Each test reads its own random half of the rows, so a true graph is rejected at close to the
stated level: 6% to 9% of the time at alpha=0.05 in our simulations, against 9% to 21% when
every test reads every row. A graph missing one edge was rejected 96% to 99% of the time at 1,000
samples. andrey markov-check GRAPH DATA runs the same check from the command line.
Bootstrap by subsampling#
andrey.bootstrap now refits the method on half the rows, drawn without replacement, instead of a
resample of all the rows with repeats. A repeated row looks like more evidence than the resample
holds; subsampling half of the rows performs as well as the best correction
(Banerjee, Andrews, and Kummerfeld, 2025), and each fit is
faster.
The averaged graph is now a valid graph of the method’s kind: a CPDAG for PC and GES, a PAG for FCI,
and a DAG without a cycle for a DAG method. On simulated data with 20 and 50 variables, it had fewer
errors than the averaged graph of 0.11.0 in 11 of 24 settings and more in none. With strict=True,
the bootstrap raises when no valid graph holds the fits’ vote; otherwise it returns the valid fit
nearest the others, and warns in the rare case that no fit is valid either. In the script, all seven
edges that appear in at least 9 fits in 10 are in the true graph.
The bootstrap example shows how to choose the settings: how many resamples, how large a subsample, and how to resample a time series.
Graph queries#
A graph now answers questions: parents, children, ancestors, descendants,
markov_blanket, d_separated for a DAG or CPDAG, and m_separated for a PAG or MAG. A CPDAG
or PAG answers as every graph in its class does. In the script, the Markov blanket of x0 in the
averaged graph is x2 and x6, as in the true graph.
Three newly supported methods#
Inter-IAMB, MMHC, and RSMAX2 move from experimental to supported. Against bnlearn 5.2.1, on 12 settings of 20 to 100 variables and 1,000 and 5,000 samples, none of them makes significantly more errors in any setting, for each graph they return. Each is faster than bnlearn: the median ratio of bnlearn’s time to Andrey’s is 3.8 for Inter-IAMB’s CPDAG, 6.6 for MMHC, and 7.8 for RSMAX2.
IAMB, Inter-IAMB, and MMPC take structure and return a CPDAG, the skeleton, or the moral graph.
IAMB and Inter-IAMB return a CPDAG by default, as bnlearn’s do.
Other changes#
BOSS, GRaSP, hill climbing, and Tabu search take n_starts and run several searches, keeping the
best-scoring one. DirectLiNGAM and its relatives take measure="skew" and measure="tanh", faster
approximations of the default pairwise measure. The former argument names that 0.11.0 kept for one
release no longer work, and a few more are renamed, with a warning that names the new one. The
upgrade guide shows the calls before and after, and the
changelog lists every change.
Limits#
Ranks see only monotone dependence. Normal scores find a relation that rises or falls. A
relation like y = x**2 for a symmetric x has no rank correlation, and the test misses it; the
kernel tests find it.
A pass is not proof. A graph with extra edges also passes the Markov check, since an extra edge removes an independence without making another one false. A pass says the data do not contradict the graph, not that its edges are right. The check also depends on its test: with Fisher’s z on skewed data, it rejects even the true graph.
Fewer false links means more missed ones. The false discovery rate control keeps far fewer
false links, and misses more true ones. Pass fdr_method="none" for the uncorrected tests.
The calibrated G-square test is more cautious, and slower. It calls fewer dependences, so at small sample sizes PC can miss weak edges it found before: in 2 of the settings we measured, its graphs were farther from the truth than in 0.11.0. It also takes about three times as long.
Try it#
pip install -U andrey-core
Help us build it#
Andrey is in alpha. The API may change between releases. When you try it on your own data, tell us what breaks or which method you need. We also welcome pull requests. The contributing guide explains how you can set up the development environment.