Causal Reasoning — From Correlation to “Why”#
Security teams routinely confuse correlation with causation. “Organizations that deployed MFA had fewer breaches” doesn’t mean MFA caused the reduction — maybe mature programs deploy MFA AND do twenty other things that actually drive the outcome. Without causal reasoning, we credit the wrong controls, waste budget, and draw conclusions that fall apart under scrutiny.
Setup#
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from scipy.stats import beta as beta_dist
from decision_security.synth import make_rng
from decision_security.causal import (
backdoor_adjustment_set, parents, children, descendants, topological_sort
)
rng = make_rng(42)
plt.rcParams.update({
"font.family": "serif",
"font.size": 10,
"axes.labelsize": 11,
"axes.titlesize": 12,
"xtick.labelsize": 9,
"ytick.labelsize": 9,
"legend.fontsize": 9,
"figure.dpi": 150,
"axes.spines.top": False,
"axes.spines.right": False,
})
PRIMARY = "#1A1A1A"
ACCENT = "#E74C3C"
DARK_BG = "#34495E"
LIGHT_GRAY = "#95A5A6"
MED_GRAY = "#7F8C8D"
VERY_LIGHT = "#BDC3C7"
Simpson’s Paradox in Patch Rates#
Two business units report aggregate incident counts. Unit A looks worse overall. But when we stratify by risk tier, Unit A has better patch rates — and lower incident rates — in both tiers. The aggregate reversal happens because Unit A operates mostly high-risk systems.
data = pd.DataFrame([
{"unit": "A", "tier": "high-risk", "systems": 1000, "patch_rate": 0.80, "incident_rate": 0.15},
{"unit": "A", "tier": "low-risk", "systems": 200, "patch_rate": 0.90, "incident_rate": 0.02},
{"unit": "B", "tier": "high-risk", "systems": 200, "patch_rate": 0.70, "incident_rate": 0.18},
{"unit": "B", "tier": "low-risk", "systems": 1000, "patch_rate": 0.85, "incident_rate": 0.03},
])
data["incidents"] = (data["systems"] * data["incident_rate"]).astype(int)
data["patched"] = (data["systems"] * data["patch_rate"]).astype(int)
print("=== By risk tier (the truth) ===")
print(data[["unit", "tier", "systems", "patch_rate", "incident_rate"]].to_string(index=False))
agg = data.groupby("unit").agg(
total_systems=("systems", "sum"),
total_incidents=("incidents", "sum"),
total_patched=("patched", "sum")
).reset_index()
agg["agg_patch_rate"] = agg["total_patched"] / agg["total_systems"]
agg["agg_incident_rate"] = agg["total_incidents"] / agg["total_systems"]
print("\n=== Aggregate (the lie) ===")
print(agg[["unit", "total_systems", "agg_patch_rate", "agg_incident_rate"]].to_string(index=False))
print("\nUnit A looks worse in aggregate despite better rates in BOTH tiers.")
# Case mix is the confounder behind the reversal.
for u in ["A", "B"]:
sub = data[data["unit"] == u]
hi = int(sub[sub["tier"] == "high-risk"]["systems"].iloc[0])
tot = int(sub["systems"].sum())
print(f"Unit {u}: {hi / tot:.0%} of systems are high-risk")
=== By risk tier (the truth) ===
unit tier systems patch_rate incident_rate
A high-risk 1000 0.80 0.15
A low-risk 200 0.90 0.02
B high-risk 200 0.70 0.18
B low-risk 1000 0.85 0.03
=== Aggregate (the lie) ===
unit total_systems agg_patch_rate agg_incident_rate
A 1200 0.816667 0.128333
B 1200 0.825000 0.055000
Unit A looks worse in aggregate despite better rates in BOTH tiers.
Unit A: 83% of systems are high-risk
Unit B: 17% of systems are high-risk
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for ax, metric, label in zip(axes, ["patch_rate", "incident_rate"], ["Patch Rate", "Incident Rate"]):
x = np.arange(2)
width = 0.3
a_vals = data[data["unit"] == "A"][metric].values
b_vals = data[data["unit"] == "B"][metric].values
ax.bar(x - width/2, a_vals, width, label="Unit A", color=DARK_BG)
ax.bar(x + width/2, b_vals, width, label="Unit B", color=ACCENT)
ax.set_xticks(x)
ax.set_xticklabels(["High-Risk", "Low-Risk"])
ax.set_ylabel(label)
ax.set_title(f"{label} by Tier")
ax.legend()
fig.suptitle("Simpson's Paradox: Unit A is better in BOTH tiers but worse in aggregate",
fontsize=11, fontweight="bold", y=1.02)
fig.text(0.5, -0.04, "Case mix drives the reversal: Unit A is 83% high-risk systems, Unit B is 83% low-risk.",
ha="center", va="top", fontsize=8, style="italic", color=MED_GRAY)
plt.tight_layout()
plt.show()
DAGs: Drawing the Causal Structure#
A directed acyclic graph (DAG) encodes assumptions about what causes what. Drawing the arrows is the modeling decision — the math just follows. Here’s a simple security DAG: budget drives training, patching, and MFA; patching and MFA both affect incidents.
nodes = ["budget", "training", "patching", "mfa", "incidents"]
edges = [
("budget", "training"),
("budget", "patching"),
("budget", "mfa"),
("training", "patching"),
("patching", "incidents"),
("mfa", "incidents"),
]
order = topological_sort(nodes, edges)
print(f"Topological order: {order}")
print(f"Parents of patching: {parents(edges, 'patching')}")
print(f"Children of budget: {children(edges, 'budget')}")
print(f"Descendants of budget: {descendants(edges, 'budget')}")
Topological order: ['budget', 'training', 'mfa', 'patching', 'incidents']
Parents of patching: {'training', 'budget'}
Children of budget: {'patching', 'training', 'mfa'}
Descendants of budget: {'patching', 'training', 'mfa', 'incidents'}
pos = {
"budget": (0.5, 1.0),
"training": (0.15, 0.5),
"patching": (0.5, 0.5),
"mfa": (0.85, 0.5),
"incidents": (0.67, 0.0),
}
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.2, 1.2)
ax.set_aspect("equal")
ax.axis("off")
for node, (x, y) in pos.items():
circle = plt.Circle((x, y), 0.08, fc=DARK_BG, ec="white", linewidth=2, zorder=3)
ax.add_patch(circle)
ax.text(x, y, node, ha="center", va="center", fontsize=8,
fontweight="bold", color="white", zorder=4)
for u, v in edges:
x0, y0 = pos[u]
x1, y1 = pos[v]
dx, dy = x1 - x0, y1 - y0
dist = np.sqrt(dx**2 + dy**2)
shrink = 0.09 / dist
ax.annotate("",
xy=(x1 - dx * shrink, y1 - dy * shrink),
xytext=(x0 + dx * shrink, y0 + dy * shrink),
arrowprops=dict(arrowstyle="->", color=ACCENT, lw=1.8))
ax.set_title("Security DAG: Budget → Controls → Incidents", fontsize=12)
plt.tight_layout()
plt.show()
Confounders, Mediators, Colliders#
Three path types determine whether a statistical association reflects a causal relationship.
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
configs = [
("Confounder", {"budget": (0.5, 0.9), "patching": (0.1, 0.1), "incidents": (0.9, 0.1)},
[("budget", "patching"), ("budget", "incidents"), ("patching", "incidents")],
"Budget causes both better patching\nAND fewer incidents independently"),
("Mediator", {"budget": (0.1, 0.5), "training": (0.5, 0.5), "patching": (0.9, 0.5)},
[("budget", "training"), ("training", "patching")],
"Training partly mediates the\nbudget → patching path"),
("Collider", {"patching": (0.1, 0.9), "incidents": (0.9, 0.9), "audit": (0.5, 0.1)},
[("patching", "audit"), ("incidents", "audit")],
"Conditioning on audit findings\nopens a spurious path"),
]
for ax, (title, positions, arrows, note) in zip(axes, configs):
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.2, 1.15)
ax.set_aspect("equal")
ax.axis("off")
ax.set_title(title, fontsize=11, fontweight="bold")
for node, (x, y) in positions.items():
label = "audit\nfinding" if node == "audit" else node
circle = plt.Circle((x, y), 0.12, fc=DARK_BG, ec="white", lw=2, zorder=3)
ax.add_patch(circle)
ax.text(x, y, label, ha="center", va="center", fontsize=7,
fontweight="bold", color="white", zorder=4)
for u, v in arrows:
x0, y0 = positions[u]
x1, y1 = positions[v]
dx, dy = x1 - x0, y1 - y0
dist = np.sqrt(dx**2 + dy**2)
shrink = 0.13 / dist
ax.annotate("",
xy=(x1 - dx * shrink, y1 - dy * shrink),
xytext=(x0 + dx * shrink, y0 + dy * shrink),
arrowprops=dict(arrowstyle="->", color=ACCENT, lw=1.8))
ax.text(0.5, -0.15, note, ha="center", va="top", fontsize=7,
style="italic", color=MED_GRAY)
plt.tight_layout()
plt.show()
# Numerical collider demo: patching and incidents are INDEPENDENT by construction,
# but conditioning on the collider (audit finding) manufactures a correlation.
n_c = 5000
patch = rng.normal(0, 1, n_c)
inc = rng.normal(0, 1, n_c)
audit_finding = (inc - patch + rng.normal(0, 0.5, n_c)) > 1 # incidents and weak patching both raise the odds of a finding
r_all = np.corrcoef(patch, inc)[0, 1]
r_sub = np.corrcoef(patch[audit_finding], inc[audit_finding])[0, 1]
print(f"corr(patching, incidents) overall: {r_all:+.3f} (~0, independent)")
print(f"corr within the audit-finding subset (collider): {r_sub:+.3f} (spurious)")
corr(patching, incidents) overall: -0.034 (~0, independent)
corr within the audit-finding subset (collider): +0.469 (spurious)
The Backdoor Criterion#
To estimate the causal effect of patching on incidents, we need to block all
“backdoor paths” — non-causal paths that flow through common causes. The
library’s backdoor_adjustment_set proposes which variables to condition on (the
treatment’s parents, a heuristic); the next cell checks it against every backdoor path.
import itertools
candidates = {"budget", "training", "mfa"}
adj_set = backdoor_adjustment_set(edges, "patching", "incidents", candidates)
print(f"Adjustment set for patching -> incidents: {adj_set}")
# Check the result: list every backdoor path (a path that starts with an arrow
# INTO patching) and find the smallest sets that block all of them.
neighbors = {}
for u, v in edges:
neighbors.setdefault(u, set()).add(v)
neighbors.setdefault(v, set()).add(u)
def all_paths(path, target):
if path[-1] == target:
yield path
return
for nxt in sorted(neighbors[path[-1]] - set(path)):
yield from all_paths(path + [nxt], target)
def blocked(path, Z):
"""A path is blocked by Z if it has a non-collider in Z, or a collider
that is not in Z and has no descendant in Z."""
for a, b, c in zip(path, path[1:], path[2:]):
if (a, b) in edges and (c, b) in edges: # b is a collider
if b not in Z and not (descendants(edges, b) & Z):
return True
elif b in Z:
return True
return False
def arrows(path):
out = path[0]
for a, b in zip(path, path[1:]):
out += (" -> " if (a, b) in edges else " <- ") + b
return out
backdoor = [p for p in all_paths(["patching"], "incidents") if (p[1], "patching") in edges]
print("\nBackdoor paths:")
for p in backdoor:
print(" " + arrows(p))
allowed = sorted(candidates - descendants(edges, "patching"))
valid = [set(Z) for r in range(len(allowed) + 1) for Z in itertools.combinations(allowed, r)
if all(blocked(p, set(Z)) for p in backdoor)]
minimal = [Z for Z in valid if not any(V < Z for V in valid)]
print(f"\nPrinted set blocks every backdoor path: {all(blocked(p, adj_set) for p in backdoor)}")
print(f"Minimal valid sets: {[sorted(Z) for Z in minimal]}")
if adj_set in minimal:
print("The printed set is minimal.")
else:
print("The printed set is valid but not minimal; any one minimal set is enough.")
Adjustment set for patching -> incidents: {'training', 'budget'}
Backdoor paths:
patching <- budget -> mfa -> incidents
patching <- training <- budget -> mfa -> incidents
Printed set blocks every backdoor path: True
Minimal valid sets: [['budget'], ['mfa']]
The printed set is valid but not minimal; any one minimal set is enough.
Naive vs Adjusted: Why the Backdoor Matters#
We generate synthetic data where budget confounds patching and incidents. The naive correlation overstates the effect of patching because it picks up the budget confound. Stratifying by budget estimates the causal effect, which we can check against the true effect built into the simulation. Here do(patch > 0.7) means pushing every system’s patch rate above 0.7, spread as it is within each budget stratum, then averaging over strata.
n = 4000
budget = rng.choice(["low", "high"], n, p=[0.4, 0.6])
is_high = budget == "high"
BETA_PARAMS = {"high": (3, 2), "low": (2, 3)}
SLOPE = 0.12
# Budget leans patching up or down, but both Betas are spread out, so every
# budget stratum has systems at both patch levels (positivity).
patch_rate = np.where(is_high,
rng.beta(*BETA_PARAMS["high"], n),
rng.beta(*BETA_PARAMS["low"], n)
)
base_incident = np.where(is_high, 0.08, 0.25)
p_incident = np.clip(base_incident + SLOPE * (1 - patch_rate), 0, 1)
incident = rng.binomial(1, p_incident)
df = pd.DataFrame({"budget": budget, "patch_rate": patch_rate, "incident": incident})
bins = pd.cut(df["patch_rate"], bins=8)
naive = df.groupby(bins, observed=True)["incident"].mean()
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(range(len(naive)), naive.values, "o-", color=DARK_BG, markersize=5)
axes[0].set_xticks(range(len(naive)))
axes[0].set_xticklabels([f"{iv.mid:.2f}" for iv in naive.index], rotation=45, fontsize=7)
axes[0].set_xlabel("Patch Rate (binned)")
axes[0].set_ylabel("Incident Rate")
axes[0].set_title("Naive: patch rate vs incidents")
for label, color in [("high", DARK_BG), ("low", ACCENT)]:
sub = df[df["budget"] == label]
bins_s = pd.cut(sub["patch_rate"], bins=6)
strat = sub.groupby(bins_s, observed=True)["incident"].mean()
axes[1].plot(range(len(strat)), strat.values, "o-", color=color,
markersize=5, label=f"{label} budget")
axes[1].set_xlabel("Patch Rate (binned)")
axes[1].set_ylabel("Incident Rate")
axes[1].set_title("Stratified by budget (confounder held fixed)")
axes[1].legend()
fig.suptitle("Confounding makes patching look more effective than it is",
fontsize=11, fontweight="bold", y=1.02)
plt.tight_layout()
plt.show()
def adjusted_incidence(df, patch_hi=0.7, patch_lo=0.3):
"""P(incident | do(patch level)) via backdoor adjustment over budget.
Returns NaN if any stratum has no systems at a patch level: the effect is
then not identified, and dropping the stratum (renormalized or not) would
compare different stratum mixes, i.e. reintroduce confounding."""
strata = sorted(df["budget"].unique())
support = {}
for s in strata:
sub = df[df["budget"] == s]
support[s] = (int((sub["patch_rate"] > patch_hi).sum()),
int((sub["patch_rate"] <= patch_lo).sum()))
positivity_ok = all(hi > 0 and lo > 0 for hi, lo in support.values())
def do_level(mask_fn):
num, wsum = 0.0, 0.0
for s in strata:
sub = df[df["budget"] == s]
cell = sub[mask_fn(sub)]
if len(cell) == 0:
return np.nan
w = len(sub)
num += cell["incident"].mean() * w
wsum += w
return num / wsum
p_hi = do_level(lambda s: s["patch_rate"] > patch_hi)
p_lo = do_level(lambda s: s["patch_rate"] <= patch_lo)
return p_hi, p_lo, positivity_ok, support
p_do_high, p_do_low, positivity_ok, support = adjusted_incidence(df)
def true_effect(weights, patch_hi=0.7, patch_lo=0.3):
"""Exact do(low) - do(high) contrast from the generating process.
Within a stratum, E[p_incident | patch in range] = base + SLOPE * (1 - E[patch | range])
(the clip never binds), so the base cancels and only truncated Beta means matter."""
def mean_in(a, b, lo, hi):
mass = beta_dist(a, b).cdf(hi) - beta_dist(a, b).cdf(lo)
return a / (a + b) * (beta_dist(a + 1, b).cdf(hi) - beta_dist(a + 1, b).cdf(lo)) / mass
return sum(w * SLOPE * (mean_in(*BETA_PARAMS[s], patch_hi, 1) - mean_in(*BETA_PARAMS[s], 0, patch_lo))
for s, w in weights.items())
# Same stratum weights as the estimator: the sample share of each budget level.
true_ate = true_effect(df["budget"].value_counts(normalize=True).to_dict())
print("Positivity check (both patch levels needed in every stratum):")
for s, (n_hi, n_lo) in support.items():
print(f" budget={s:>4}: high-patch (>0.7) = {n_hi:>4}, low-patch (<=0.3) = {n_lo:>4}")
print(f" positivity satisfied: {positivity_ok}")
print()
print(f"Backdoor-adjusted P(incident | do(high patch)) = {p_do_high:.3f}")
print(f"Backdoor-adjusted P(incident | do(low patch)) = {p_do_low:.3f}")
print(f"Adjusted ATE (low - high) = {p_do_low - p_do_high:+.3f} (estimate; positive = patching reduces incidents)")
print(f"True ATE (low - high) = {true_ate:+.3f} (from the generating process)")
Positivity check (both patch levels needed in every stratum):
budget=high: high-patch (>0.7) = 864, low-patch (<=0.3) = 222
budget= low: high-patch (>0.7) = 126, low-patch (<=0.3) = 577
positivity satisfied: True
Backdoor-adjusted P(incident | do(high patch)) = 0.193
Backdoor-adjusted P(incident | do(low patch)) = 0.223
Adjusted ATE (low - high) = +0.030 (estimate; positive = patching reduces incidents)
True ATE (low - high) = +0.071 (from the generating process)
Observing vs Intervening#
The core insight of causal inference: P(incidents | observe patch=high) differs from P(incidents | do(patch=high)). Observing high patch rates includes the confound — high-budget orgs self-select into high patching. Intervening (forcing patch rates high for everyone) breaks the incoming arrows to patching and isolates its total effect (here also the direct effect, since nothing mediates patching → incidents).
high_patchers = df[df["patch_rate"] > 0.7]
low_patchers = df[df["patch_rate"] <= 0.3]
p_obs_high = high_patchers["incident"].mean()
p_obs_low = low_patchers["incident"].mean()
# do(): the backdoor-adjusted estimates from the previous section.
print("=== Observe vs Intervene ===")
print(f"P(incident | observe patch > 0.7) = {p_obs_high:.3f}")
print(f"P(incident | do(patch > 0.7)) = {p_do_high:.3f} (adjusted for budget)")
print()
print(f"P(incident | observe patch <= 0.3) = {p_obs_low:.3f}")
print(f"P(incident | do(patch <= 0.3)) = {p_do_low:.3f} (adjusted for budget)")
print()
print(f"Naive effect: {p_obs_low - p_obs_high:+.3f} (confounded)")
print(f"Causal effect: {p_do_low - p_do_high:+.3f} (adjusted estimate, with sampling noise)")
print(f"True effect: {true_ate:+.3f} (from the generating process)")
# Bootstrap the adjusted estimate to show how far sampling noise alone moves it.
boot = []
for i in range(200):
p_hi_b, p_lo_b, _, _ = adjusted_incidence(df.sample(len(df), replace=True, random_state=i))
boot.append(p_lo_b - p_hi_b)
boot_sd = np.nanstd(boot, ddof=1)
gap_sd = (p_do_low - p_do_high - true_ate) / boot_sd
print(f"Bootstrap SD: {boot_sd:.3f} -> estimate is {gap_sd:+.1f} SD from the true effect "
f"({'within' if abs(gap_sd) < 2 else 'outside'} typical sampling noise)")
if (p_do_low - p_do_high > 0) != (true_ate > 0):
print("WARNING: adjusted effect has the wrong sign; check positivity.")
=== Observe vs Intervene ===
P(incident | observe patch > 0.7) = 0.121
P(incident | do(patch > 0.7)) = 0.193 (adjusted for budget)
P(incident | observe patch <= 0.3) = 0.270
P(incident | do(patch <= 0.3)) = 0.223 (adjusted for budget)
Naive effect: +0.149 (confounded)
Causal effect: +0.030 (adjusted estimate, with sampling noise)
True effect: +0.071 (from the generating process)
Bootstrap SD: 0.026 -> estimate is -1.6 SD from the true effect (within typical sampling noise)
The gap between estimate and truth is part of the lesson. Adjustment removes the confounding, not the noise. The adjusted estimate is one draw from a sampling distribution centred on the true effect: across repeated samples it averages the true value, but any single sample can land well above or below it. A correct method on a finite sample still gives an uncertain answer, so report an adjusted effect with its uncertainty (a bootstrap interval or standard error), never as a single number.
When adjustment silently breaks. The estimate above is trustworthy only because the positivity check passed — every budget stratum contained systems at both patch levels. If a stratum is empty at one level (a positivity violation), the data say nothing about that level in that stratum. Dropping the stratum, with or without renormalizing the remaining weights, compares different stratum mixes at the two levels, so the estimate is confounded again and can flip sign.
Pitfalls#
Correlation between control deployment and fewer incidents almost always has confounders. Maturity, budget, and culture drive both control adoption and outcomes.
Conditioning on a collider creates phantom associations. “Among organizations that passed their audit” is a collider condition — it makes unrelated variables appear correlated because both contribute to passing.
You need the counterfactual. The question isn’t “did incidents drop after we deployed MFA?” It’s “what would have happened without MFA?” Without that comparison, post-hoc reasoning dominates.
DAGs encode assumptions, not facts. Two analysts can draw different DAGs for the same situation and reach opposite conclusions. The diagram should be debated, not the regression output.