Measuring Control Effectiveness — Observation vs Experiment#

We want to know whether a security control causes fewer incidents. But we can’t run randomized controlled trials on production systems — we can’t randomly deny half our business units EDR and see who gets breached. Most measurements we have are observational, and observational data carries structural biases that distort our conclusions.

This notebook walks through the core problems — selection bias, confounding, collider bias, time trends — and the techniques that extract causal signal from non-experimental data. Section 6 applies them to a public incident database.

Setup#

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from decision_security.synth import make_rng, sample
from decision_security.causal import backdoor_adjustment_set
from decision_security.montecarlo import simulate_aggregate_losses, make_lognormal_severity
from decision_security.bayes import beta_update

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"

1. The Selection Bias Problem#

Organizations that deploy EDR are not a random sample. They tend to be more security-mature overall — bigger budgets, better staffing, more process discipline. A naive comparison of incident rates between EDR-adopters and non-adopters confounds the effect of EDR with everything else mature organizations do.

EDR looks more effective than it actually is, because we’re crediting it for outcomes that maturity would have produced anyway.

n_orgs = 200
maturity = rng.uniform(0, 1, n_orgs)
has_edr = rng.random(n_orgs) < maturity                  # mature orgs adopt EDR more often
base_rate = rng.poisson(10 * (1 - 0.6 * maturity))       # maturity cuts incidents on its own

# EDR removes each incident with probability true_effect, so E[incidents] falls by exactly that.
true_effect = 0.30
incident_rate = np.where(has_edr, rng.binomial(base_rate, 1 - true_effect), base_rate)

df = pd.DataFrame({
    "maturity": maturity,
    "has_edr": has_edr,
    "incidents": incident_rate,
})

mean_edr = df[df["has_edr"]]["incidents"].mean()
mean_no_edr = df[~df["has_edr"]]["incidents"].mean()
naive_reduction = 1 - mean_edr / mean_no_edr

print("=== Naive comparison ===")
print("Estimand: relative reduction in E[incidents]")
print(f"Mean incidents WITH EDR:    {mean_edr:.1f}  (n={df['has_edr'].sum()})")
print(f"Mean incidents WITHOUT EDR: {mean_no_edr:.1f}  (n={(~df['has_edr']).sum()})")
print(f"Naive reduction estimate:   {naive_reduction:.0%}")
print(f"True causal effect of EDR:  {true_effect:.0%}")
=== Naive comparison ===
Estimand: relative reduction in E[incidents]
Mean incidents WITH EDR:    4.0  (n=96)
Mean incidents WITHOUT EDR: 8.2  (n=104)
Naive reduction estimate:   51%
True causal effect of EDR:  30%
fig, ax = plt.subplots(figsize=(8, 5))

edr_mask = df["has_edr"]
ax.scatter(df[edr_mask]["maturity"], df[edr_mask]["incidents"],
           c=ACCENT, alpha=0.6, s=30, label="Has EDR", edgecolors="white", linewidth=0.3)
ax.scatter(df[~edr_mask]["maturity"], df[~edr_mask]["incidents"],
           c=DARK_BG, alpha=0.6, s=30, label="No EDR", edgecolors="white", linewidth=0.3)

ax.set_xlabel("Security Maturity")
ax.set_ylabel("Annual Incidents")
ax.set_title("Selection bias: EDR adopters are already more mature")
ax.legend()
plt.tight_layout()
plt.show()
../_images/d26aeeff040ef0adcb327caa2c31b8923e7b0efe9ef708afc2075650ed8e5b32.png

Check it in your own environment.

  1. Ask who adopted the control first, and why. If it went to the strongest teams or the most exposed systems first, the raw comparison is biased one way or the other.

  2. Compare within groups that share that factor (maturity, exposure, asset criticality), not across the whole estate.

  3. Look at the trend before deployment: if adopters were already improving, a before-after drop is partly that trend (Section 4, difference-in-differences).

  4. Check how many teams sit in each group before trusting the number. A handful per group is anecdote.

2. Stratified Analysis — Reducing (Not Eliminating) the Bias#

If maturity confounds the EDR-incidents relationship, we can control for it by comparing EDR vs no-EDR within groups of similar maturity. Each quartile holds maturity roughly fixed, so the selection bias inside it is much smaller than in the naive pooled comparison.

But “roughly fixed” is not “fixed.” The naive estimate of 51% falls to about 38% after stratifying, closer to the true 30% but not on top of it. Some of the remaining gap is residual confounding: maturity still varies within each quartile. Most of it, in this sample, is noise; some quartile cells hold fewer than ten organizations. The within-quartile estimates (30%, 30%, 45%, 46%) bounce around for the same reason. With more organizations, stratification here lands within a few points of 30%. Finer adjustment (regression on maturity, matching) would shrink the gap further.

df["quartile"] = pd.qcut(df["maturity"], 4, labels=["Q1 (low)", "Q2", "Q3", "Q4 (high)"])

strat = df.groupby(["quartile", "has_edr"], observed=True)["incidents"].agg(["mean", "size"]).unstack()
strat = pd.DataFrame({
    "No EDR": strat[("mean", False)], "n (no EDR)": strat[("size", False)],
    "Has EDR": strat[("mean", True)], "n (EDR)": strat[("size", True)],
})
strat["Within-stratum reduction"] = 1 - strat["Has EDR"] / strat["No EDR"]
strat_mean = strat["Within-stratum reduction"].mean()

print("=== Stratified by maturity quartile ===")
print(strat.to_string(float_format="{:.1f}".format,
                      formatters={"Within-stratum reduction": "{:.0%}".format}))
print(f"\nNaive aggregate reduction:      {naive_reduction:.0%}")
print(f"Mean within-stratum reduction:  {strat_mean:.0%}")
print(f"True causal effect:             {true_effect:.0%}")
=== Stratified by maturity quartile ===
           No EDR  n (no EDR)  Has EDR  n (EDR) Within-stratum reduction
quartile                                                                
Q1 (low)      9.9          41      6.9        9                      30%
Q2            7.3          40      5.1       10                      30%
Q3            7.2          16      4.0       34                      45%
Q4 (high)     5.9           7      3.2       43                      46%

Naive aggregate reduction:      51%
Mean within-stratum reduction:  38%
True causal effect:             30%
fig, ax = plt.subplots(figsize=(8, 5))

quartiles = strat.index.tolist()
x = np.arange(len(quartiles))
width = 0.3

bars_no = ax.bar(x - width/2, strat["No EDR"].values, width,
                 label="No EDR", color=DARK_BG)
bars_yes = ax.bar(x + width/2, strat["Has EDR"].values, width,
                  label="Has EDR", color=ACCENT)

ax.set_xticks(x)
ax.set_xticklabels(quartiles)
ax.set_xlabel("Maturity Quartile")
ax.set_ylabel("Mean Annual Incidents")
ax.set_title(f"Stratifying shrinks the naive gap ({naive_reduction:.0%} → {strat_mean:.0%}; true {true_effect:.0%})")
ax.legend()

for i, q in enumerate(quartiles):
    reduction = strat.loc[q, "Within-stratum reduction"]
    y_pos = max(strat.loc[q, "No EDR"], strat.loc[q, "Has EDR"]) + 0.3
    ax.text(i, y_pos, f"{reduction:.0%}", ha="center", fontsize=8, color=LIGHT_GRAY)

plt.tight_layout()
plt.show()
../_images/b8d5f65a03e4b378704cd22a5766f27e747e2e23b7254125f3bc6e16b293a17f.png

3. The DAG Tells You What to Control For#

The causal DAG for this problem has three nodes: maturity causes both EDR adoption and lower incidents; EDR also directly reduces incidents. The backdoor path EDR <- maturity -> incidents creates a non-causal association. Conditioning on maturity blocks it.

But not all conditioning is safe. If we add a collider — say, detected incidents, which depend on both having EDR and having incidents — then conditioning on detected opens a spurious path and introduces bias that wasn’t there before.

edges_dag = [
    ("maturity", "edr"),
    ("maturity", "incidents"),
    ("edr", "incidents"),
]

adj = backdoor_adjustment_set(edges_dag, treatment="edr", outcome="incidents",
                              candidates={"maturity"})
print(f"Backdoor adjustment set: {adj}")
print(f"Conditioning on maturity blocks the backdoor path edr <- maturity -> incidents.")
Backdoor adjustment set: {'maturity'}
Conditioning on maturity blocks the backdoor path edr <- maturity -> incidents.
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# --- DAG 1: correct adjustment ---
ax = axes[0]
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.15, 1.15)
ax.set_aspect("equal")
ax.axis("off")
ax.set_title("Correct: condition on maturity", fontsize=11, fontweight="bold")

pos1 = {"maturity": (0.5, 1.0), "edr": (0.1, 0.2), "incidents": (0.9, 0.2)}
node_colors_1 = {"maturity": ACCENT, "edr": DARK_BG, "incidents": DARK_BG}
edges_draw_1 = [("maturity", "edr"), ("maturity", "incidents"), ("edr", "incidents")]

for node, (x, y) in pos1.items():
    fc = node_colors_1[node]
    circle = plt.Circle((x, y), 0.1, fc=fc, ec="white", lw=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_draw_1:
    x0, y0 = pos1[u]
    x1, y1 = pos1[v]
    dx, dy = x1 - x0, y1 - y0
    dist = np.sqrt(dx**2 + dy**2)
    shrink = 0.11 / dist
    style = "--" if u == "maturity" else "-"
    color = LIGHT_GRAY if u == "maturity" else ACCENT
    ax.annotate("",
        xy=(x1 - dx * shrink, y1 - dy * shrink),
        xytext=(x0 + dx * shrink, y0 + dy * shrink),
        arrowprops=dict(arrowstyle="->", color=color, lw=1.8,
                        linestyle=style))

ax.text(0.5, -0.1, "Maturity (red) = conditioned on.\nBackdoor path blocked.",
        ha="center", fontsize=8, style="italic", color=LIGHT_GRAY)

# --- DAG 2: collider bias ---
ax = axes[1]
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.15, 1.15)
ax.set_aspect("equal")
ax.axis("off")
ax.set_title("Wrong: condition on detected (collider)", fontsize=11, fontweight="bold")

pos2 = {"maturity": (0.5, 1.0), "edr": (0.1, 0.55), "incidents": (0.9, 0.55),
        "detected": (0.5, 0.1)}
node_colors_2 = {"maturity": DARK_BG, "edr": DARK_BG, "incidents": DARK_BG,
                 "detected": ACCENT}
edges_draw_2 = [("maturity", "edr"), ("maturity", "incidents"),
                ("edr", "incidents"), ("edr", "detected"), ("incidents", "detected")]

for node, (x, y) in pos2.items():
    fc = node_colors_2[node]
    circle = plt.Circle((x, y), 0.1, fc=fc, ec="white", lw=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_draw_2:
    x0, y0 = pos2[u]
    x1, y1 = pos2[v]
    dx, dy = x1 - x0, y1 - y0
    dist = np.sqrt(dx**2 + dy**2)
    shrink = 0.11 / dist
    color = ACCENT if v == "detected" else DARK_BG
    style = "-"
    ax.annotate("",
        xy=(x1 - dx * shrink, y1 - dy * shrink),
        xytext=(x0 + dx * shrink, y0 + dy * shrink),
        arrowprops=dict(arrowstyle="->", color=color, lw=1.8,
                        linestyle=style))

ax.text(0.5, -0.1, "Detected (red) = collider. Conditioning on it\nopens a spurious edr-incidents path.",
        ha="center", fontsize=8, style="italic", color=LIGHT_GRAY)

plt.tight_layout()
plt.show()
../_images/2e56e8e85f603ab8984f6daebc559111e67c1987bdec2c5a6fbf32a2dfe6fc87.png
# EDR and incidents are INDEPENDENT by construction, yet conditioning on `detected`
# manufactures an association.
rng_col = make_rng(3)
m = 8000
edr = rng_col.random(m) < 0.5
incidents_col = rng_col.poisson(9, m)
# A unit is detected if it runs EDR OR has many incidents. Detection is selective
# (~30% of units), so conditioning on it bites.
z = -2.5 + 0.5 * (incidents_col - 9) + 2.5 * edr
detected = rng_col.random(m) < 1 / (1 + np.exp(-z))

corr_all = np.corrcoef(edr.astype(float), incidents_col)[0, 1]
corr_detected = np.corrcoef(edr[detected].astype(float), incidents_col[detected])[0, 1]

print("=== Collider demo: conditioning on `detected` ===")
print(f"Detected fraction: {detected.mean():.0%}")
print(f"corr(edr, incidents), all units:     {corr_all:+.3f}   (independent by construction)")
print(f"corr(edr, incidents), detected only: {corr_detected:+.3f}   (spurious)")
print("\nAmong detected units, a non-EDR unit had to accumulate incidents to be seen,")
print("so EDR looks associated with fewer incidents even though it changed nothing.")
=== Collider demo: conditioning on `detected` ===
Detected fraction: 31%
corr(edr, incidents), all units:     -0.000   (independent by construction)
corr(edr, incidents), detected only: -0.249   (spurious)

Among detected units, a non-EDR unit had to accumulate incidents to be seen,
so EDR looks associated with fewer incidents even though it changed nothing.

4. Before-After vs Difference-in-Differences#

A common approach: measure the incident rate before and after deploying MFA, then attribute the drop to MFA. The problem is that incident rates trend downward over time anyway — other improvements, threat landscape shifts, regression to the mean.

Difference-in-differences (DiD) compares the change in the treatment group to the change in a control group that didn’t deploy MFA. The difference between those differences isolates the treatment effect.

DiD rests on one identifying assumption: parallel trends — absent the intervention, the treatment and control groups would have moved in parallel. If the MFA-adopting units were already improving faster than the control for unrelated reasons, DiD misattributes that divergence to MFA. The assumption is about the unobserved counterfactual, so it can be argued for (e.g. by showing parallel pre-period trends) but never fully proven.

months = np.arange(1, 25)
mfa_deploy_month = 12
true_mfa_effect = -2.0  # incidents per month

trend = -0.15 * months  # shared improving baseline
noise_treat = rng.normal(0, 0.8, 24)
noise_ctrl = rng.normal(0, 0.8, 24)

control = 12 + trend + noise_ctrl
mfa_effect = np.where(months > mfa_deploy_month, true_mfa_effect, 0)
treatment = 12 + trend + mfa_effect + noise_treat

before = treatment[:mfa_deploy_month].mean()
after = treatment[mfa_deploy_month:].mean()
naive_ba = after - before

treat_diff = naive_ba
ctrl_diff = control[mfa_deploy_month:].mean() - control[:mfa_deploy_month].mean()
did_estimate = treat_diff - ctrl_diff

print("=== Before-After vs DiD ===")
print(f"Treatment before:  {before:.1f}")
print(f"Treatment after:   {after:.1f}")
print(f"Naive before-after estimate: {naive_ba:.1f} incidents/month (includes the trend)")
print(f"\nTreatment change: {treat_diff:.1f}")
print(f"Control change:   {ctrl_diff:.1f}")
print(f"DiD estimate:     {did_estimate:.1f} incidents/month")
print(f"True MFA effect:  {true_mfa_effect:.1f} incidents/month")
=== Before-After vs DiD ===
Treatment before:  11.0
Treatment after:   7.2
Naive before-after estimate: -3.8 incidents/month (includes the trend)

Treatment change: -3.8
Control change:   -1.7
DiD estimate:     -2.2 incidents/month
True MFA effect:  -2.0 incidents/month
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Panel (a): time series
ax = axes[0]
ax.plot(months, treatment, "o-", color=ACCENT, markersize=4, label="Treatment (MFA deployed)")
ax.plot(months, control, "o-", color=DARK_BG, markersize=4, label="Control (no MFA)")
ax.axvline(mfa_deploy_month + 0.5, color=LIGHT_GRAY, ls="--", lw=1.5, label="MFA deployed")
ax.set_xlabel("Month")
ax.set_ylabel("Incidents / month")
ax.set_title("(a) Raw time series")
ax.legend(fontsize=8)

# Panel (b): estimates comparison
ax = axes[1]
labels = ["Naive\nbefore-after", "Difference-in-\ndifferences", "True\neffect"]
values = [naive_ba, did_estimate, true_mfa_effect]
colors = [LIGHT_GRAY, ACCENT, DARK_BG]
bars = ax.bar(range(3), values, color=colors, width=0.5)
ax.set_xticks(range(3))
ax.set_xticklabels(labels)
ax.set_ylabel("Estimated effect (incidents/month)")
ax.set_title("(b) Naive BA overstates the MFA effect")
ax.axhline(0, color=PRIMARY, lw=0.5)

for bar, val in zip(bars, values):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 0.15,
            f"{val:.1f}", ha="center", va="top", fontsize=9, fontweight="bold",
            color="white")

plt.tight_layout()
plt.show()
../_images/2ca7c0e2f33bf9e053d11a9bef5cd550564530eb6b615907d4717bc85133f76e.png

5. Bayesian Evidence Accumulation#

When you cannot run a controlled experiment, you can still accumulate observational evidence over time and update your beliefs about control effectiveness. Here we estimate the relative reduction in the incident rate from quarterly incident counts across matched business units, starting from a weak prior whose median is 0% effect (spanning harmful to helpful).

The value of the Bayesian view is not that it reaches an answer sooner. It is that at every sample size it reports a posterior that is wide when data is thin and narrow once evidence accumulates. You read your current uncertainty off the posterior instead of collapsing it into a significant/not-significant verdict.

from scipy.stats import beta as beta_dist

# We model incident COUNTS rather than "which unit won a matchup", so the posterior
# is about the rate reduction, not P(EDR unit has fewer incidents).
#
# Trick: with matched EDR / non-EDR pairs and equal exposure, the share of incidents
# landing in the EDR arm is theta = lambda_edr / (lambda_edr + lambda_ctrl). Each incident
# is a Bernoulli "which arm" draw, so a Beta prior on theta is conjugate, and
# reduction = (1 - 2*theta) / (1 - theta).

true_lambda_ctrl = 10.0                                  # incidents / BU-quarter, no EDR
true_reduction = 0.30
true_lambda_edr = true_lambda_ctrl * (1 - true_reduction)
n_bu_pairs = 8                                           # matched pairs observed per quarter

def reduction_from_theta(theta):
    return (1 - 2 * theta) / (1 - theta)

# Beta(5, 5): weak (10 pseudo-incidents), median theta 0.5 = median reduction 0%.
# The transform is nonlinear, so the prior MEAN reduction is negative.
a_prior, b_prior = 5.0, 5.0

a, b = a_prior, b_prior
history = [(a, b)]
for _ in range(12):
    edr_inc = rng.poisson(true_lambda_edr, n_bu_pairs).sum()
    ctrl_inc = rng.poisson(true_lambda_ctrl, n_bu_pairs).sum()
    a, b = beta_update(a, b, edr_inc, ctrl_inc)          # incidents in EDR arm = "successes"
    history.append((a, b))

# The transform is monotone, so pushing theta draws through it gives valid quantiles.
post_rng = np.random.default_rng(2024)
red_draws = reduction_from_theta(post_rng.beta(a, b, 50000))

print("=== Bayesian evidence accumulation (rate-ratio model) ===")
print("Estimand: relative reduction in E[incidents] = 1 - lambda_edr/lambda_ctrl")
print(f"Prior on theta: Beta({a_prior:.0f}, {b_prior:.0f})  ->  median reduction {reduction_from_theta(0.5):.0%}")
print(f"After 12 quarters: Beta({a:.0f}, {b:.0f})")
print(f"Posterior mean reduction:      {red_draws.mean():.1%}")
print(f"True reduction:                {true_reduction:.0%}")
print(f"Posterior 90% credible interval: "
      f"[{np.percentile(red_draws, 5):.1%}, {np.percentile(red_draws, 95):.1%}]")
=== Bayesian evidence accumulation (rate-ratio model) ===
Estimand: relative reduction in E[incidents] = 1 - lambda_edr/lambda_ctrl
Prior on theta: Beta(5, 5)  ->  median reduction 0%
After 12 quarters: Beta(662, 954)
Posterior mean reduction:      30.6%
True reduction:                30%
Posterior 90% credible interval: [24.6%, 36.2%]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

plot_rng = np.random.default_rng(7)

def reduction_draws(a_h, b_h, n=40000):
    return reduction_from_theta(plot_rng.beta(a_h, b_h, n))

# Panel (a): posterior of the reduction at several time steps
ax = axes[0]
steps_to_show = [0, 3, 6, 12]
alphas = [0.3, 0.5, 0.7, 1.0]
bins = np.linspace(-1.0, 0.7, 160)
centers = 0.5 * (bins[:-1] + bins[1:])
for step, alpha in zip(steps_to_show, alphas):
    a_s, b_s = history[step]
    draws = reduction_draws(a_s, b_s)
    # Normalize by ALL draws: the prior has mass below -100%, outside the plotted range.
    counts, _ = np.histogram(draws, bins=bins)
    dens = counts / (len(draws) * np.diff(bins))
    label = "Prior" if step == 0 else f"After Q{step}"
    color = LIGHT_GRAY if step == 0 else ACCENT
    ax.plot(centers, dens, color=color, alpha=alpha, lw=2, label=label)
    ax.fill_between(centers, dens, alpha=alpha * 0.15, color=color)

ax.axvline(true_reduction, color=DARK_BG, ls="--", lw=1.5, label=f"True reduction ({true_reduction:.0%})")
ax.axvline(0, color=PRIMARY, lw=0.6)
ax.set_xlabel("Relative reduction in E[incidents]")
ax.set_ylabel("Density")
ax.set_title("(a) Posterior concentrates on the true reduction")
ax.legend(fontsize=8)

# Panel (b): posterior mean reduction and 90% CI over time
ax = axes[1]
quarters = range(len(history))
means, ci_lo, ci_hi = [], [], []
for a_h, b_h in history:
    d = reduction_draws(a_h, b_h)
    means.append(d.mean())
    ci_lo.append(np.percentile(d, 5))
    ci_hi.append(np.percentile(d, 95))

ax.fill_between(quarters, ci_lo, ci_hi, alpha=0.2, color=ACCENT)
ax.plot(quarters, means, "o-", color=ACCENT, markersize=5, label="Posterior mean")
ax.axhline(true_reduction, color=DARK_BG, ls="--", lw=1.5, label=f"True reduction ({true_reduction:.0%})")
ax.axhline(0, color=PRIMARY, lw=0.6)
ax.set_xlabel("Quarter (0 = prior)")
ax.set_ylabel("Estimated reduction")
ax.set_title("(b) 90% credible interval narrows with evidence")
ax.legend(fontsize=8)
ax.set_xticks(range(0, 13, 2))

plt.tight_layout()
plt.show()
../_images/892c07bca7ae821d41243d98d30ae0e652d5c744618d1d96ebf5a1df4fb9f704.png
# Negative confounding: a WAF goes on the most-attacked apps, which have more incidents
# to begin with, so a naive comparison makes an effective control look HARMFUL.
# Stratifying on exposure recovers the sign and most of the benefit. Exposure still varies
# within each quartile, and the WAF apps sit at its high end, so some negative
# confounding survives.
rng_nc = make_rng(11)
n = 600
exposure = rng_nc.uniform(0, 1, n)
has_waf = rng_nc.random(n) < exposure
base = rng_nc.poisson(2 + 20 * exposure)                 # incidents climb steeply with exposure
true_effect_waf = 0.30
incidents_nc = np.where(has_waf, rng_nc.binomial(base, 1 - true_effect_waf), base)

naive_nc = 1 - incidents_nc[has_waf].mean() / incidents_nc[~has_waf].mean()

qx = pd.qcut(exposure, 4, labels=False)
strata_red = [
    1 - incidents_nc[has_waf & (qx == k)].mean() / incidents_nc[~has_waf & (qx == k)].mean()
    for k in range(4)
]

print("=== Negative confounding: control deployed where risk is highest ===")
print(f"Naive reduction estimate:       {naive_nc:+.0%}   "
      f"({'looks HARMFUL' if naive_nc < 0 else 'looks weaker than it is'})")
print(f"Mean within-stratum reduction:  {np.nanmean(strata_red):.0%}")
print(f"True causal effect of the WAF:  {true_effect_waf:.0%}")
=== Negative confounding: control deployed where risk is highest ===
Naive reduction estimate:       -33%   (looks HARMFUL)
Mean within-stratum reduction:  26%
True causal effect of the WAF:  30%

6. Real Data: An Incident Database Selected on the Outcome#

The VERIS Community Database (VCDB) is a public collection of security incidents that were reported in news stories, breach notifications, regulator filings. Most records note whether data was disclosed; only 27% of records note whether the affected data was stored encrypted.

The question: does encryption prevent disclosure? Compare the disclosure rate of incidents with encrypted data against incidents with unencrypted data. The snapshot in data/ keeps the result reproducible; load_vcdb(refresh=True) downloads the current database.

from _data import load_vcdb

vcdb = load_vcdb()
encrypted = vcdb["stored_encrypted"] & ~vcdb["stored_unencrypted"]
unencrypted = vcdb["stored_unencrypted"] & ~vcdb["stored_encrypted"]
known = vcdb[encrypted | unencrypted].assign(
    storage=np.where(encrypted[encrypted | unencrypted], "Stored encrypted", "Stored unencrypted"))

table = pd.crosstab(known["storage"], known["data_disclosure"])
table = table[[c for c in ["Yes", "Potentially", "No", "Unknown"] if c in table.columns]]
print(f"VCDB: {len(vcdb):,} incidents; {len(known):,} ({len(known) / len(vcdb):.0%}) "
      f"record whether stored data was encrypted\n")
print(table.to_string())

p_yes = (known["data_disclosure"] == "Yes").groupby(known["storage"]).mean()
print(f"\nConfirmed disclosure, encrypted:   {p_yes['Stored encrypted']:.0%}")
print(f"Confirmed disclosure, unencrypted: {p_yes['Stored unencrypted']:.0%}")
print(f"Naive reduction from encryption:   {1 - p_yes['Stored encrypted'] / p_yes['Stored unencrypted']:.0%}")

# data_disclosure exists only when VERIS codes a confidentiality section.
coded = vcdb["data_disclosure"].notna()
print(f"\nNo confidentiality section: {(~coded).sum()}")
print(f"data_disclosure = 'No': {(vcdb['data_disclosure'] == 'No').sum()} of {coded.sum():,} coded")

enc_actions = vcdb.loc[encrypted, "actions"].fillna("")
print(f"Encrypted records involving hacking or misuse: "
      f"{enc_actions.str.contains('hacking|misuse').sum()} of {len(enc_actions)}")
VCDB: 10,047 incidents; 2,731 (27%) record whether stored data was encrypted

data_disclosure      Yes  Potentially  No  Unknown
storage                                           
Stored encrypted      91           36   1        2
Stored unencrypted  2029          565   3        4

Confirmed disclosure, encrypted:   70%
Confirmed disclosure, unencrypted: 78%
Naive reduction from encryption:   10%

No confidentiality section: 675
data_disclosure = 'No': 83 of 9,372 coded
Encrypted records involving hacking or misuse: 88 of 130

Taken at face value, encryption barely helps.

Now look at the “No” column. Almost no incident whose confidentiality was coded has “no disclosure” (83 of 9,372). VERIS records the confidentiality fields, including encryption state, mainly when data was at stake. That is not a finding about encryption; it is how the database is built. An incident enters VCDB when it becomes public, and incidents mostly become public because data was, or may have been, disclosed. Many breach-notification laws exempt data that was encrypted, so a lost encrypted laptop rarely produces a filing or a news story and rarely becomes a record. Examples: most US state laws, HIPAA (only “unsecured” PHI) and GDPR Art. 34(3)(a).

In DAG terms, reported is a collider: encryption affects it (no notification needed) and disclosure affects it (notification required). The database conditions on reported = yes – the selection from Section 3, on real data. The encrypted incidents that do get in mostly involve hacking or insider misuse (88 of 130), where encryption at rest typically does not apply: an attacker or insider with valid credentials, or keys stored with the data.

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

shares = table.div(table.sum(axis=1), axis=0)
colors = {"Yes": ACCENT, "Potentially": LIGHT_GRAY, "No": DARK_BG, "Unknown": "#D6DBDF"}
left = np.zeros(len(shares))
ax = axes[0]
for col in shares.columns:
    ax.barh(shares.index, shares[col], left=left, color=colors[col], label=col, edgecolor="white")
    left += shares[col].values
ax.set_xlabel("Share of incidents")
ax.set_title("Disclosure outcome by storage (VCDB)")
ax.legend(title="data_disclosure", loc="lower right", fontsize=8)

per_year = vcdb.dropna(subset=["year"]).groupby("year").size()
per_year = per_year[per_year.index >= 2005]
ax = axes[1]
ax.bar(per_year.index.astype(int), per_year.values, color=DARK_BG, alpha=0.85)
ax.set_xlabel("Incident year")
ax.set_ylabel("Records in VCDB")
ax.set_title("Records per year track curation, not incident rates")

plt.tight_layout()
plt.show()

print("Left: almost no record has 'No' disclosure, whatever the storage.")
print("Right: record counts follow how much volunteers coded each year, not incidents.")
../_images/8ee7605032dbdfd2a8b5f2d7354cb4b4c2bf7677ecbb2f81a8686864b897df3e.png
Left: almost no record has 'No' disclosure, whatever the storage.
Right: record counts follow how much volunteers coded each year, not incidents.

What VCDB can and cannot tell you. It is a good source for how incidents that became public happened: actions, assets, attacker types. It cannot estimate the effect of a control that decides whether an incident becomes public, and its record counts are not incident rates. To measure encryption’s effect you need the denominator – all loss events, including the quiet ones – for example an organization’s own log of lost and stolen devices.

Data: VERIS Community Database (CC BY-SA 4.0), vz-risk/VCDB.

7. Pitfalls#

  • Selection bias is the default, not the exception. Organizations that deploy a control are systematically different from those that don’t.

  • Controls are deployed where risk is highest. This can make effective controls look ineffective (or even harmful) — the opposite of the EDR example. A firewall deployed at the most-attacked perimeter can still show more incidents than a quiet segment with no firewall. This is negative confounding, and it’s at least as dangerous as the positive kind.

  • Conditioning on post-treatment variables can introduce collider bias. Analyzing only “detected” incidents, or only incidents that triggered a response, conditions on a collider and distorts the estimated effect of the control.

  • Before-after designs confound treatment effects with time trends. Without a control group, you cannot separate the effect of the intervention from the secular trend.

  • Small observational samples can be worse than no data if they create false confidence. A single business unit’s before-after comparison, treated as conclusive, can anchor decisions in noise. The Bayesian approach represents uncertainty explicitly — but only helps if you respect the width of the posterior.

  • Public incident databases are selected on the outcome. A record exists because the incident became public, and controls like encryption decide whether it does. Comparing records with and without the control conditions on a collider: the database cannot show the incidents the control kept quiet.