Bootstrap: how stable is each edge?#
A method fitted to one sample returns one graph, and another sample from the same source would
give a somewhat different one. andrey.bootstrap refits the method on many resamples of the rows
and counts how often each edge appears. Here PC runs on data simulated on the ALARM network, whose
true graph we know, to show how to read those counts, what the averaged graph is, and which edges
to trust.
Start#
Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.13.0 from PyPI;
local kernels use the installed environment.
import os
import sys
IN_COLAB = "google.colab" in sys.modules
if IN_COLAB:
os.environ["ANDREY_DEVICE"] = "cpu"
!pip install --quiet "andrey-core[viz]==0.13.0" pandas
import warnings
import numpy as np
import pandas as pd
import andrey
import andrey.viz
from andrey.data import load_dataset
from andrey.metrics import shd, to_cpdag
ALARM is a published network of 37 variables and 46 arcs. load_dataset simulates 1,000 rows
on it from a linear-Gaussian model and lists the columns in an order drawn from the seed. From
these data PC can recover at best the CPDAG of the true graph, so that is what we score against.
alarm = load_dataset("alarm", n=1000, seed=1)
data, truth = alarm.data, alarm.graph
true_cpdag = to_cpdag(truth)
upper = np.triu_indices(truth.n_nodes, 1) # each pair of variables once
true_pair = (truth.to_numpy() != 0)[upper]
print(data.shape[0], "rows,", data.shape[1], "variables,", int(true_pair.sum()), "true edges")
1000 rows, 37 variables, 46 true edges
Edge frequencies rank the edges#
bootstrap fits PC on 50 resamples, each half of the rows drawn without replacement, the default.
frequency("absent") holds the share of the fits in which each pair has no edge, so one minus it
is the share in which the pair is adjacent. Here are the pairs grouped by that share:
out = andrey.bootstrap("pc", data, n_resamples=50, seed=1)
adjacent = (1 - out.frequency("absent"))[upper]
bins = [0.1, 0.3, 0.5, 0.7, 0.9, 1.0]
rows = []
for low, high in zip(bins[:-1], bins[1:]):
inside = (adjacent > low) & (adjacent <= high)
rows.append(
{
"adjacent in": f"more than {low:.0%}, up to {high:.0%}",
"pairs": int(inside.sum()),
"true pairs": int((inside & true_pair).sum()),
}
)
pd.DataFrame(rows)
| adjacent in | pairs | true pairs | |
|---|---|---|---|
| 0 | more than 10%, up to 30% | 6 | 2 |
| 1 | more than 30%, up to 50% | 4 | 3 |
| 2 | more than 50%, up to 70% | 0 | 0 |
| 3 | more than 70%, up to 90% | 2 | 2 |
| 4 | more than 90%, up to 100% | 32 | 31 |
The pairs adjacent in more than 90% of the fits are nearly all true edges. The shares in between are not probabilities: most pairs adjacent in 30% to 50% of the fits are true edges here, not 30% to 50% of them. A frequency says how stable an edge is when the rows change, which ranks the edges; how often the stable ones are right depends on the method, the data, and the sample size.
The averaged graph#
out.structure is the averaged graph. Each pair takes its most frequent state in the fits, with
no edge counted as one of the states. PC returns CPDAGs, and a vote of CPDAGs is often not a
CPDAG, so the vote is then replaced by the CPDAG of a DAG with the same adjacencies, v-structures,
and directed edges. The result is a valid CPDAG, here with fewer errors than PC on all the rows:
single = andrey.pc(data)
print("single fit: ", shd(single.structure, true_cpdag), "errors")
print("averaged graph:", shd(out.structure, true_cpdag), "errors")
print("kind:", out.structure.kind, "| built from:", out.metadata["averaged_graph"])
dag = out.structure.to_dag(tie_break="label")
print(
"the CPDAG of one of its DAGs:",
np.array_equal(to_cpdag(dag).to_numpy(), out.structure.to_numpy()),
)
single fit: 24 errors
averaged graph: 21 errors
kind: cpdag | built from: vote
the CPDAG of one of its DAGs: True
andrey.viz.draw(out.structure)
Subsampling or the bootstrap#
resampling="bootstrap" draws as many rows as the data have, with replacement, instead of half of
them without. A repeated row makes a test see more evidence than the resample holds, so weak and
false edges survive in more of the fits. sample_fraction sets the share of the rows in a
subsample: a larger share gives fits closer to the fit on all the rows, so the fits agree more and
their frequencies say less about how much the result varies.
runs = {
"subsample, 50% of the rows": {},
"subsample, 80% of the rows": {"sample_fraction": 0.8},
"bootstrap": {"resampling": "bootstrap"},
}
results, raised = {}, {}
rows = []
for name, settings in runs.items():
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
result = andrey.bootstrap("pc", data, n_resamples=50, seed=1, **settings)
results[name] = result
raised[name] = [str(w.message) for w in caught if issubclass(w.category, andrey.AndreyWarning)]
share = (1 - result.frequency("absent"))[upper]
above, middle = share > 0.9, (share > 0.5) & (share <= 0.9)
rows.append(
{
"resampling": name,
"errors": shd(result.structure, true_cpdag),
"above 90%: true of all": f"{(above & true_pair).sum()} of {above.sum()}",
"50% to 90%: true of all": f"{(middle & true_pair).sum()} of {middle.sum()}",
"averaged graph": result.metadata["averaged_graph"],
}
)
pd.DataFrame(rows)
| resampling | errors | above 90%: true of all | 50% to 90%: true of all | averaged graph | |
|---|---|---|---|---|---|
| 0 | subsample, 50% of the rows | 21 | 31 of 32 | 2 of 2 | vote |
| 1 | subsample, 80% of the rows | 15 | 33 of 34 | 4 of 5 | vote |
| 2 | bootstrap | 19 | 32 of 33 | 4 of 5 | nearest_resample |
When the vote has no valid graph#
With the bootstrap, the vote of these fits has no consistent extension: no DAG keeps all its
adjacencies, v-structures, and directed edges, so no CPDAG holds the vote. bootstrap then
returns the resample graph nearest the others among those a DAG extends, the fit whose state of
each pair the most fits share. It is a graph PC itself returned, and a warning and
metadata["averaged_graph"] say so:
boot = results["bootstrap"]
print("\n".join(raised["bootstrap"]))
print("built from:", boot.metadata["averaged_graph"], "| kind:", boot.structure.kind)
the vote of the fits has no valid cpdag: PDAG admits no consistent DAG extension: no DAG has its skeleton, its v-structures, and its directed edges, because the undirected edges among nodes [2, 3, 4, 5, 11, 13, 15, 29, 30, 31, 32] cannot be oriented without adding a v-structure or closing a directed cycle. The averaged graph is resample 25, the valid fit nearest the others; strict=True raises instead
built from: nearest_resample | kind: cpdag
strict=True raises instead, so you can decide from the frequencies, which keep the whole vote
either way:
try:
andrey.bootstrap("pc", data, n_resamples=50, seed=1, resampling="bootstrap", strict=True)
except ValueError as error:
print("ValueError:", error)
ValueError: the vote of the fits has no valid cpdag: PDAG admits no consistent DAG extension: no DAG has its skeleton, its v-structures, and its directed edges, because the undirected edges among nodes [2, 3, 4, 5, 11, 13, 15, 29, 30, 31, 32] cannot be oriented without adding a v-structure or closing a directed cycle. strict=False returns the fit nearest the others instead
Choosing the settings#
The defaults, 100 resamples of half of the rows each, suit most uses. Each setting trades one thing for another; here is how to choose it.
How many resamples#
A frequency is a share of n_resamples fits, so it carries its own sampling noise: for a pair
adjacent in a share p of all possible fits, about sqrt(p (1 - p) / n_resamples), 0.10 at 25
resamples and 0.05 at 100 for a pair near one half. Five seeds at each count show it on the pairs
whose frequency lies between 20% and 80%:
rows = []
for n_resamples in (25, 50, 100, 200):
fits = [andrey.bootstrap("pc", data, n_resamples=n_resamples, seed=s) for s in range(1, 6)]
shares = np.array([(1 - fit.frequency("absent"))[upper] for fit in fits])
mean = shares.mean(axis=0)
middle = (mean > 0.2) & (mean < 0.8)
rows.append(
{
"resamples": n_resamples,
"pairs between 20% and 80%": int(middle.sum()),
"spread over seeds (SD)": shares[:, middle].std(axis=0, ddof=1).mean(),
"expected SD": np.sqrt(mean[middle] * (1 - mean[middle]) / n_resamples).mean(),
"pairs above 90%, per seed": [int((s > 0.9).sum()) for s in shares],
}
)
pd.DataFrame(rows).round(3)
| resamples | pairs between 20% and 80% | spread over seeds (SD) | expected SD | pairs above 90%, per seed | |
|---|---|---|---|---|---|
| 0 | 25 | 7 | 0.090 | 0.090 | [32, 32, 32, 33, 32] |
| 1 | 50 | 7 | 0.070 | 0.065 | [32, 32, 32, 33, 32] |
| 2 | 100 | 7 | 0.045 | 0.046 | [32, 32, 32, 32, 32] |
| 3 | 200 | 7 | 0.029 | 0.032 | [32, 32, 32, 32, 32] |
The spread follows the formula. The pairs above 90% barely move at any count, so 100 resamples rank the edges well. When a decision rests on one frequency near a threshold, such as whether a pair passes 0.5, use more: at 200 the spread is near 0.03. The time grows in proportion.
How large a subsample#
sample_fraction sets how many rows each fit sees. A larger share gives fits closer to the fit on
all the rows, so more pairs reach each level, true and false; but the resamples then share most of
their rows, so the frequencies overstate how stable an edge is. Half of the rows keeps the fits
apart and performs as well as the best correction of the bootstrap (Banerjee, Andrews, and
Kummerfeld 2025). Here are the shares on ALARM at 200, 1,000, and 5,000 rows:
rows = []
for n_rows in (200, 1000, 5000):
sample = load_dataset("alarm", n=n_rows, seed=1).data
for name, settings in {"50%": {}, "80%": {"sample_fraction": 0.8}}.items():
result = andrey.bootstrap("pc", sample, n_resamples=50, seed=1, **settings)
adjacent = (1 - result.frequency("absent"))[upper]
rows.append(
{
"rows": n_rows,
"share of rows": name,
"above 90%": int((adjacent > 0.9).sum()),
"true above 90%": int(((adjacent > 0.9) & true_pair).sum()),
"above 50%": int((adjacent > 0.5).sum()),
"true above 50%": int(((adjacent > 0.5) & true_pair).sum()),
}
)
pd.DataFrame(rows)
| rows | share of rows | above 90% | true above 90% | above 50% | true above 50% | |
|---|---|---|---|---|---|---|
| 0 | 200 | 50% | 23 | 23 | 24 | 24 |
| 1 | 200 | 80% | 25 | 25 | 29 | 28 |
| 2 | 1000 | 50% | 32 | 31 | 34 | 33 |
| 3 | 1000 | 80% | 34 | 33 | 39 | 37 |
| 4 | 5000 | 50% | 35 | 35 | 38 | 38 |
| 5 | 5000 | 80% | 38 | 38 | 41 | 39 |
With few rows, a half-size fit loses the weak edges, so fewer pairs pass each level. Keep the default to rank edges by stability; a larger share finds more of the weak edges on small data, at the cost of frequencies that read higher. Whatever the share, compare frequencies only between runs with the same settings.
Method settings and workers#
Every keyword that bootstrap does not take goes to each fit: alpha, indep_test, knowledge,
or any other argument of the method. A lower alpha asks each test for more evidence, so fewer
pairs are adjacent in each fit:
for alpha in (0.05, 0.01):
result = andrey.bootstrap("pc", data, n_resamples=50, seed=1, alpha=alpha)
adjacent = (1 - result.frequency("absent"))[upper]
print(f"alpha {alpha}: {int((adjacent > 0.9).sum())} pairs above 90%,", end=" ")
print(f"{int(((adjacent > 0.9) & true_pair).sum())} of them true")
alpha 0.05: 32 pairs above 90%, 31 of them true
alpha 0.01: 30 pairs above 90%, 29 of them true
andrey.config(n_jobs=4) runs the fits in four processes. Starting the processes takes time, so it
pays off when each fit takes a second or more; a PC fit on ALARM takes a few milliseconds, and one
process is faster. In a script, run it under if __name__ == "__main__":, since each process
imports the script again.
import time
for n_jobs in (1, 4):
with andrey.config(n_jobs=n_jobs):
start = time.perf_counter()
andrey.bootstrap("pc", data, n_resamples=100, seed=1)
print(f"{n_jobs} process(es): {time.perf_counter() - start:.1f} s")
1 process(es): 0.6 s
4 process(es): 2.0 s
Time series#
A time-series method resamples time, not independent rows. A subsample is one stretch of
consecutive steps, half of the series by default, so every lag inside it is intact.
resampling="bootstrap" draws a series of the same length from blocks of consecutive steps, the
moving block bootstrap, with block_length steps each, by default max(ceil(n ** (1/3)), 10 * tau_max). Here PCMCI looks for the links at lag 1 in 1,000 steps of four variables:
from andrey.data import sample_scm
series = sample_scm(
graph="erdos_renyi",
functional="var",
noise="uniform",
d=4,
n=1000,
seed=1,
density=1.0,
scale="raw",
)
true_links = {(s, t) for s, t, _ in series.graph.lag(1).oriented_edges()}
print(len(true_links), "true links at lag 1")
for name, settings in {"subsample": {}, "block bootstrap": {"resampling": "bootstrap"}}.items():
result = andrey.bootstrap("pcmci", series.data, n_resamples=50, seed=1, tau_max=1, **settings)
print(f"\n{name}, block length {result.metadata['block_length']}:")
for edge in result.table(min_frequency=0.1):
link = (edge["source"], edge["target"])
verdict = "true" if link in true_links else "false"
print(f" {link[0]} -> {link[1]} at lag 1: {edge['frequency']:.2f} ({verdict})")
10 true links at lag 1
subsample, block length None:
0 -> 0 at lag 1: 1.00 (true)
0 -> 1 at lag 1: 1.00 (true)
1 -> 0 at lag 1: 1.00 (true)
1 -> 1 at lag 1: 1.00 (true)
2 -> 1 at lag 1: 1.00 (true)
2 -> 3 at lag 1: 1.00 (true)
2 -> 2 at lag 1: 0.98 (true)
1 -> 3 at lag 1: 0.84 (true)
3 -> 1 at lag 1: 0.20 (false)
0 -> 3 at lag 1: 0.12 (true)
block bootstrap, block length 10:
0 -> 0 at lag 1: 1.00 (true)
1 -> 0 at lag 1: 1.00 (true)
1 -> 1 at lag 1: 1.00 (true)
2 -> 1 at lag 1: 1.00 (true)
2 -> 2 at lag 1: 1.00 (true)
2 -> 3 at lag 1: 1.00 (true)
1 -> 3 at lag 1: 0.96 (true)
0 -> 1 at lag 1: 0.92 (true)
0 -> 3 at lag 1: 0.64 (true)
3 -> 3 at lag 1: 0.44 (true)
3 -> 2 at lag 1: 0.16 (false)
3 -> 1 at lag 1: 0.14 (false)
Half of the series, 500 steps, misses the two weakest true links: 0 -> 3 comes back in 12% of the
fits and 3 -> 3 in none. The block bootstrap’s full-length resamples find them in 64% and 44% of
its fits. In both, a false link appears in a fifth of the fits or fewer. On a short series, or when
the weak links matter, a larger sample_fraction or the block bootstrap keeps them; a block must
be long enough to hold the lags, which the default ensures.
Which edges to trust, and how to report them#
Trust the edges most fits agree on, and report each one with its frequency. Here are the averaged graph’s most and least stable edges, each with the share of the fits in which the pair is adjacent and the share with this edge:
position = {name: k for k, name in enumerate(data.columns)}
glyph = {"directed": "->", "undirected": "--"}
pairs = 1 - out.frequency("absent")
report = pd.DataFrame(
{
"edge": f"{source} {glyph[kind]} {target}",
"adjacent": pairs[position[source], position[target]],
"this edge": out.frequency(glyph[kind])[position[source], position[target]],
"true pair": bool(truth.to_numpy()[position[source], position[target]]),
}
for source, target, kind in out.structure.oriented_edges()
).sort_values("adjacent", ascending=False, kind="stable")
pd.concat([report.head(4), report.tail(4)])
| edge | adjacent | this edge | true pair | |
|---|---|---|---|---|
| 0 | PAP -- PULMEMBOLUS | 1.00 | 0.90 | True |
| 1 | HISTORY -- LVFAILURE | 1.00 | 1.00 | True |
| 2 | CO -> BP | 1.00 | 0.60 | True |
| 4 | TPR -> BP | 1.00 | 0.60 | True |
| 26 | VENTTUBE -> VENTLUNG | 0.98 | 0.82 | True |
| 32 | INTUBATION -> SHUNT | 0.98 | 0.92 | True |
| 3 | CO -- STROKEVOLUME | 0.80 | 0.68 | True |
| 14 | LVFAILURE -> LVEDVOLUME | 0.80 | 0.78 | True |
threshold asks more of the averaged graph: out.average(0.9) keeps only the pairs adjacent in
more than 90% of the fits, without refitting, which here drops true edges as well as false ones.
With each edge, report what produced the frequencies:
the method and its settings, the number of resamples, the resampling and its share of rows, the
seed, the threshold, and whether the averaged graph is the vote or the nearest resample.
strong = out.average(0.9)
print(len(strong.oriented_edges()), "edges,", shd(strong, true_cpdag), "errors")
print(
{
"method": out.metadata["method"],
"n_resamples": out.n_resamples,
"threshold": out.threshold,
**{
key: out.metadata[key]
for key in ("resampling", "sample_fraction", "seed", "averaged_graph")
},
}
)
32 edges, 25 errors
{'method': 'pc', 'n_resamples': 50, 'threshold': 0.5, 'resampling': 'subsample', 'sample_fraction': 0.5, 'seed': 1, 'averaged_graph': 'vote'}