Ready for real data#

Cover: a pantry shelf of mixed jars, cartons, and sacks, one jar amber, drawn in navy lines

Shu Wan · October 29, 2026 · Andrey 0.11.0

TL;DR#

Gaps: PC, FCI, and the other constraint-based methods read, for each test, the samples it can use, so they keep the third of the adults that dropping every incomplete one would lose.
Units: the same graph whether blood pressure is in mmHg or kPa.
Mixed types: declare a column as a category, and the default test becomes one for mixed data.
What you already know: background knowledge in every family, 27 of 34 methods. On the survey, it removes all 10 arrows that pointed into age or sex.
How sure each edge is: andrey.bootstrap refits on bootstrap samples and gives each edge’s frequency and an averaged graph.
New defaults: the extended BIC for ten score searches, which cuts GES’s median errors on our simulations from 87 to 25, and standardized simulated data.

In 0.11.0, we make Andrey ready for real-world chaos.

Real tables are messy. Values go missing. Columns come in different units. Some columns are categories, not numbers. You need to know which edges to trust, and you already know things the data cannot tell you.

Andrey now handles each of these. We take one real table, 5,569 adults from a national health survey, from raw columns to a causal graph, in one short script.

A real table#

The National Health and Nutrition Examination Survey (NHANES) measures the health of people across the United States. Its 2017-2018 cycle is public domain, published by the US Centers for Disease Control and Prevention. We join six of its files on each participant’s number and keep the adults, aged 20 and over, with eleven columns. The first five look like this:

age	sex	income	bmi	waist	systolic	diastolic	pulse	hba1c	cholesterol	hdl
66	female	0.82	31.7	101.8	`NaN`	`NaN`	52	6.2	157	60
66	female	1.63	23.7	88.2	`NaN`	`NaN`	68	6.2	209	88
75	female	0.41	38.9	113.0	120	66	74	6.3	176	65
56	male	5.00	21.3	86.6	108	68	62	5.7	238	72
67	male	2.65	23.5	99.7	104	70	90	5.6	184	48

Over all 5,569 adults:

Column	What it measures	Mean (SD)	Range	Missing
`age`	age in years	51.5 (17.8)	20 to 80	0%
`sex`	male or female, a category	51% female		0%
`income`	family income as a ratio to the poverty line	2.6 (1.6)	0 to 5	14%
`bmi`	body mass index, weight in kg over height in m squared	29.9 (7.4)	14.2 to 86.2	7%
`waist`	waist circumference, cm	100.8 (17.0)	57.9 to 169.5	11%
`systolic`	systolic blood pressure, mmHg	126.8 (19.8)	72 to 228	16%
`diastolic`	diastolic blood pressure, mmHg	72.1 (13.6)	0 to 136	16%
`pulse`	resting pulse, beats per minute	71.7 (11.7)	34 to 136	10%
`hba1c`	glycated hemoglobin, %	5.9 (1.1)	3.8 to 16.2	10%
`cholesterol`	total cholesterol, mg/dL	188.1 (41.3)	76 to 446	11%
`hdl`	HDL cholesterol, mg/dL	53.3 (15.7)	10 to 189	11%

Not everyone came to the exam, had their blood pressure read, or gave blood, and some did not report their income. Only 3,694 of the 5,569 adults have every value. We do not use the survey weights, so the graph describes this sample, not the US population.

The whole script is below; each section that follows runs one part of it.

example.py#

import numpy as np
import pandas as pd

import andrey
from andrey.metrics import shd

# NHANES 2017-2018 (CDC, public domain): six files joined on the participant's number, adults 20
# and over. Survey weights are not used.
URL = "https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/{}.xpt"
COLUMNS = {
    "RIDAGEYR": "age", "RIAGENDR": "sex", "INDFMPIR": "income", "BMXBMI": "bmi",
    "BMXWAIST": "waist", "BPXSY1": "systolic", "BPXDI1": "diastolic", "BPXPLS": "pulse",
    "LBXGH": "hba1c", "LBXTC": "cholesterol", "LBDHDD": "hdl",
}  # fmt: skip


def edges(fit):
    a = np.asarray(fit.structure.adjacency)
    return int(((a != 0) | (a.T != 0)).sum() // 2)


def into_age_or_sex(fit):
    oriented = fit.structure.oriented_edges()
    return sum(kind == "directed" and target in ("age", "sex") for _, target, kind in oriented)


if __name__ == "__main__":  # worker processes re-import this script
    files = ["DEMO_J", "BMX_J", "BPX_J", "GHB_J", "TCHOL_J", "HDL_J"]
    table = pd.read_sas(URL.format(files[0]), format="xport")
    for name in files[1:]:
        table = table.merge(pd.read_sas(URL.format(name), format="xport"), on="SEQN", how="left")
    data = table[table.RIDAGEYR >= 20][list(COLUMNS)].rename(columns=COLUMNS)
    # Sex is a category, so each test is one for mixed data.
    data["sex"] = data["sex"].map({1: "male", 2: "female"}).astype("category")
    print(f"{len(data):,} adults; {len(data.dropna()):,} with no value missing")

    # Gaps: each test reads the samples where its variables are observed.
    every_adult = andrey.pc(data)
    complete = andrey.pc(data.dropna())
    differ = shd(every_adult.structure, complete.structure)
    print(f"gaps: {edges(every_adult)} edges from every adult, {edges(complete)} from the")
    print(f"  adults with no value missing; the two differ in {differ} places")
    print(f"test: {every_adult.metadata['indep_test']}")

    # Units: blood pressure in kPa and cholesterol in mmol/L give the same graph.
    metric = data.assign(
        systolic=data["systolic"] * 0.133322,
        diastolic=data["diastolic"] * 0.133322,
        cholesterol=data["cholesterol"] * 0.02586,
        hdl=data["hdl"] * 0.02586,
    )
    print(f"units: {shd(andrey.pc(metric).structure, every_adult.structure)} differences")

    # What you already know: nothing causes age or sex.
    known = andrey.Knowledge(tiers=[["age", "sex"], [c for c in data if c not in ("age", "sex")]])
    guided = andrey.pc(data, knowledge=known)
    before, after = into_age_or_sex(every_adult), into_age_or_sex(guided)
    print(f"arrows into age or sex: {before} without it, {after} with it")

    # How sure each edge is: 100 fits on bootstrap samples, four at a time.
    with andrey.config(num_workers=4):
        boot = andrey.bootstrap("pc", data, n_boot=100, seed=0, knowledge=known)
    print("edges in at least half the fits:")
    MARK = {"directed": "->", "undirected": "--"}
    for edge in boot.table(min_frequency=0.5):
        mark = MARK.get(edge["type"], edge["type"])
        print(f"  {edge['source']} {mark} {edge['target']}: {edge['frequency']:.2f}")

    # Hand it on: the averaged graph as Graphviz DOT, which DoWhy reads.
    print(boot.structure.to_dot().splitlines()[0])

example.py#

5,569 adults; 3,694 with no value missing
gaps: 25 edges from every adult, 25 from the
  adults with no value missing; the two differ in 6 places
test: cond_gauss
units: 0 differences
arrows into age or sex: 10 without it, 0 with it
edges in at least half the fits:
  age -> systolic: 1.00
  age -> pulse: 1.00
  age -> hba1c: 1.00
  sex -> bmi: 1.00
  sex -> systolic: 1.00
  sex -> diastolic: 1.00
  sex -> pulse: 1.00
  sex -> hdl: 1.00
  age -> hdl: 0.95
  waist -> hdl: 0.86
  sex -> cholesterol: 0.71
  cholesterol -> hdl: 0.71
  hba1c -> hdl: 0.70
  age -> income: 0.66
  sex -> waist: 0.63
  hdl -> income: 0.58
  sex -> hba1c: 0.57
  bmi -> systolic: 0.56
  waist -> hba1c: 0.54
  systolic -> diastolic: 0.54
  age -> waist: 0.53
  bmi -> waist: 0.52
  hba1c -> pulse: 0.51
  pulse -> diastolic: 0.50
digraph {

Gaps#

The usual way to handle a gap is to drop every sample that has one, listwise deletion. Here that drops 1,875 adults, a third of the table, each missing at least one value.

PC, FCI, CD-NOD, IAMB, Inter-IAMB, MMPC, PCMCI, PCMCI+, and LPCMCI take the table as it is, with NaN where a value is missing. Each conditional-independence test reads the samples where all of its variables are observed, test-wise deletion. A test of age and HbA1c given waist keeps an adult whose blood pressure is missing, because it does not need it.

On the survey, the two ways agree on most of the graph. Both find 25 edges, and they differ in 6 places, only 2 of which are about whether an edge is there at all: sex and waist are adjacent only on the complete adults, and sex and cholesterol only on every adult. The other 4 are directions. Where the two disagree, the bootstrap below shows which edges hold up.

On wider tables, listwise deletion runs out of data. In our simulations of 20 variables and 1,000 samples, with values missing completely at random, dropping the incomplete samples leaves about 125 when a tenth of the values are missing, about 10 when a fifth are, and almost none when 30% are. Test-wise deletion keeps every sample and was never significantly worse than listwise deletion. PC’s errors rise only from 4.0 on the complete data to 6.5 with a tenth missing and 7.5 with 30%. The release card gives every size, rate, and test.

Units#

Blood pressure in mmHg or kPa, cholesterol in mg/dL or mmol/L: a unit is a choice of scale, and it should not change the causes. The script converts four columns to SI units and fits again. The graph is the same, with 0 differences.

Mixed types#

Sex is a category, not a number. In pandas, give it the category type, and Andrey reads it as one:

data["sex"] = data["sex"].map({1: "male", 2: "female"}).astype("category")

With no test named, PC now picks it from the column types you declare: gsq when every column is a category, the conditional Gaussian test cond_gauss when some are, and Fisher’s z when none are. Here it picks cond_gauss, and the fit records the test it used. Missing values work with it too.

A column of small whole numbers, such as 1 and 2 for sex, may also be a category. Andrey does not guess: it warns and reads the column as numbers until you declare it.

What you already know#

Nothing in a body causes a person’s age or sex. The data do not know that, and PC, left to itself, draws 10 arrows into age or sex, among them pulse into age and HDL cholesterol into sex.

Tell it with andrey.Knowledge: age and sex in the first tier, every other column after them, so no edge may point back into the first tier.

known = andrey.Knowledge(tiers=[["age", "sex"], [c for c in data if c not in ("age", "sex")]])

With it, 0 arrows point into age or sex. Knowledge also takes required, forbidden, and allowed edges, and 27 of the 34 methods accept it, in every family.

How sure is each edge#

One fit gives one graph. It does not say which edges would come back with a different sample of adults. andrey.bootstrap refits on 100 bootstrap samples, each drawn with replacement, and counts how often each edge appears. With andrey.config(num_workers=4), the fits run four at a time.

Of the edges in at least half the fits, eight appear in all of them: age into systolic pressure, pulse, and HbA1c, and sex into BMI, blood pressure, pulse, and HDL cholesterol. Then come age into HDL cholesterol (0.95) and waist into HDL cholesterol (0.86). Further down, near one half, are the edges a single fit may or may not draw, such as pulse and diastolic pressure (0.50).

Read the frequencies as a ranking. In our simulations, of the edges seen in more than nine fits in ten, 98% to 100% were true, and of those seen in fewer than one in ten, about 1% or fewer. The frequencies in between order the edges from likely to unlikely, but are not probabilities.

Hand it on#

The bootstrap also returns an averaged graph: the pairs adjacent in at least half the fits, each with its most frequent edge. Like every Andrey graph, it converts to NetworkX, a pcalg matrix, NumPy, or Graphviz DOT text, which DoWhy reads to estimate an effect. A CPDAG turns into one of its DAGs with to_dag().

New defaults#

The extended BIC. GES, GIES, XGES, hill climbing, Tabu search, BOSS, GRaSP, exact search, MMHC, and RSMAX2 score continuous data with the extended BIC (Chen and Chen, 2008). Its penalty grows with the number of variables, so a large graph keeps fewer false edges. Over 299 simulated datasets, GES’s median errors fall from 87 to 25, and BOSS’s from 69 to 15. GFCI keeps the standard BIC. Pass ebic_gamma=0.0 for the graphs of before.

Standardized simulated data. andrey.data.sample_scm and andrey.data.generate give every column a mean of 0 and a standard deviation of 1, so no method can read the causal order off the variances a simulation happens to give. Pass scale="raw" for the values of before.

One name for each argument. The level of a test is alpha in every method, the BIC’s weight is penalty_discount, the largest lag is tau_max, and the random seed is seed. The former names still work in this release, with a warning that gives the new one. The upgrade guide shows the common calls before and after, and the changelog lists every renamed argument, changed default, and changed output.

Limits#

Why values are missing matters. Test-wise deletion, like dropping incomplete samples, works when whether a value is missing has nothing to do with the values themselves. In the survey, that means people skipped the blood test for reasons unrelated to their health. If people with high blood sugar skipped it more often, the samples a test reads would no longer be a fair sample of the adults, and the test could find a dependence that is not there or miss one that is. FCI tolerates more of this than PC; the getting-started guide explains when.

The mixed-data tests are experimental. When a column is a category, Andrey picks the conditional Gaussian test. It runs on the survey, but we have not yet benchmarked it the way we have Fisher’s z. The release card’s results for missing values are on continuous data.

A frequency is not a probability. The bootstrap refits the method on resamples of the same 5,569 adults, so a frequency tells you how stable an edge is, not whether it is true. If the method makes the same mistake every time, the wrong edge appears in every fit and gets a frequency of 1. In our simulations, the most frequent edges were almost always true, so the frequencies are good for ranking edges. The averaged graph was not consistently more accurate than a single fit, so we do not yet recommend it over one.

The graph describes this sample and these columns. NHANES samples some groups more than others on purpose and gives each participant a weight to correct for it. We fit without the weights, so the graph describes these 5,569 adults, not the US population. PC also assumes that no variable left out of the table causes two that are in it. Diet and exercise are not in the table, yet they can affect both BMI and blood pressure; PC could then draw an edge between those two that is really the work of the missing variables. FCI allows for such hidden causes.

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.

The boundary Calibration