E-values and Robustness Values

62.1 Partial R² and the robustness value (40 min)

For regression estimates, Cinelli & Hazlett (2020) describe a hidden confounder Z by two partial R² values:

• R²(D~Z | X): the share of the treatment's variation, after the controls, that Z explains;
• R²(Y~Z | D, X): the share of the outcome's variation, after treatment and controls, that Z explains.

Everything follows from the regression output you already have. The robustness value RV_q is the smallest strength (equal on both sides) a confounder needs to reduce the estimate by a fraction q:

f = q × |t| / √dof          RV_q = ½ (√(f⁴ + 4f²) − f²)

A confounder weaker than RV on both sides can't do it; one stronger on both sides can. To ask "could a confounder push the estimate below a threshold c?", use t = (estimate − c) / SE with q = 1.

Create causalkit/sensitivity.py:

"""Sensitivity to an omitted confounder in OLS (Cinelli & Hazlett, 2020) and E-values (VanderWeele & Ding, 2017).

The language: a hidden confounder Z is described by two partial R2 values,
    r2dz: share of the treatment's residual variance (after the controls) that Z explains,
    r2yz: share of the outcome's residual variance (after treatment and controls) that Z explains.
"""
import numpy as np


def partial_r2(t: float, dof: float) -> float:
    """Partial R2 of a regressor with t-statistic t: t^2 / (t^2 + dof)."""
    return t**2 / (t**2 + dof)


def robustness_value(t: float, dof: float, q: float = 1.0) -> float:
    """Minimum strength (same partial R2 with treatment and outcome) a confounder needs to reduce the
    estimate by a fraction q (q = 1: to zero). Pass t = (estimate - target) / se for other targets with q = 1."""
    f = q * abs(t) / np.sqrt(dof)
    return 0.5 * (np.sqrt(f**4 + 4 * f**2) - f**2)


def bias(se: float, dof: float, r2yz: float, r2dz: float) -> float:
    """Largest possible bias of the OLS estimate from a confounder with strengths (r2yz, r2dz)."""
    return se * np.sqrt(dof * r2yz * r2dz / (1 - r2dz))


def benchmark(r2yxj: float, r2dxj: float, k: float = 1.0) -> tuple[float, float]:
    """Strengths (r2yz, r2dz) of a confounder k times as strong as observed covariate Xj.
    r2dxj: partial R2 of Xj with the treatment given the other controls;
    r2yxj: partial R2 of Xj with the outcome given the treatment and the other controls."""
    r2dz = k * r2dxj / (1 - r2dxj)
    r2zxj = k * r2dxj**2 / ((1 - k * r2dxj) * (1 - r2dxj))
    r2yz = ((np.sqrt(k) + np.sqrt(r2zxj)) / np.sqrt(1 - r2zxj)) ** 2 * (r2yxj / (1 - r2yxj))
    return min(r2yz, 1.0), min(r2dz, 1.0)


def e_value(rr: float) -> float:
    """E-value: the risk ratio a confounder needs with BOTH treatment and outcome to explain away rr."""
    rr = rr if rr >= 1 else 1 / rr
    return rr + np.sqrt(rr * (rr - 1))

62.2 Back to the loyalty programme (50 min)

Module 4's question: does joining the loyalty programme raise next-year spending by more than $150, the break-even (Day 40)? World A has no hidden confounder; world B has the unmeasured brand affinity. Create day62/robustness_value.py:

# day62/robustness_value.py   (python -m day62.robustness_value)
"""
Back to Module 4's loyalty programme. In world B, an unmeasured brand affinity drives both joining and
spending. The robustness value (Cinelli & Hazlett, 2020) asks: how strongly would a hidden confounder have
to be associated with BOTH joining and spending (as a share of their remaining variation) to change the
conclusion? The business line from Day 40: joining pays off if it raises spending by more than $150.
"""
import numpy as np  # noqa: F401  (used inside the formula: np.sqrt)
import statsmodels.formula.api as smf

from causalkit.loyalty import truth, world
from causalkit.sensitivity import partial_r2, robustness_value

FORMULA = ("spend ~ joined + pre_spend + I(pre_spend**2) + pre_orders + engaged + tenure_years"
           " + np.sqrt(tenure_years) + age + mobile")
BREAK_EVEN = 150.0

if __name__ == "__main__":
    for title, hidden in [("world A (no hidden confounder)", 0.0), ("world B (hidden brand affinity)", 1.0)]:
        full = world(hidden=hidden)
        fit = smf.ols(FORMULA, full.to_pandas()).fit()
        est, se, t, dof = fit.params["joined"], fit.bse["joined"], fit.tvalues["joined"], fit.df_resid
        print(f"{title}: OLS {est:.1f} (SE {se:.1f}); true ATE {truth(full)['ATE']:.1f}")
        print(f"   partial R2 of joining with spend: {partial_r2(t, dof):.1%}")
        print(f"   robustness value to explain the effect away entirely:   {robustness_value(t, dof):.1%}")
        print(f"   robustness value to push it below the $150 break-even:  {robustness_value((est - BREAK_EVEN) / se, dof):.1%}")

python -m day62.robustness_value

world A (no hidden confounder): OLS 151.7 (SE 1.9); true ATE 158.4
   partial R2 of joining with spend: 23.4%
   robustness value to explain the effect away entirely:   42.1%
   robustness value to push it below the $150 break-even:  0.6%
world B (hidden brand affinity): OLS 226.1 (SE 2.3); true ATE 158.4
   partial R2 of joining with spend: 32.4%
   robustness value to explain the effect away entirely:   49.3%
   robustness value to push it below the $150 break-even:  20.7%

Read world B like a skeptical reader would, without the truth column: the estimate is $226, and a confounder would need 20.7% partial R² with both joining and spending to push it below $150. That sounds robust. Is a confounder that strong plausible? The robustness value can't answer that. Benchmarks can.

62.3 Benchmarks a CFO can follow (45 min)

Engagement (opening marketing emails) is a confounder we did measure. How strong is it, and what would a hidden confounder 1, 2 or 3 times as strong do? Create day62/benchmarks.py:

# day62/benchmarks.py   (python -m day62.benchmarks)
"""
'A confounder with partial R2 of 21%' still means nothing to a CFO. Benchmark it: how strong is ENGAGEMENT
(opening marketing emails), a confounder we did measure? Then ask what a hidden confounder 1, 2 or 3 times
as strong as engagement would do. Finally, because it's a simulator, measure the real hidden affinity.
"""
import numpy as np  # noqa: F401  (used inside the formula: np.sqrt)
import polars as pl
import statsmodels.formula.api as smf

from causalkit.loyalty import world
from causalkit.sensitivity import benchmark, bias, partial_r2
from day62.robustness_value import BREAK_EVEN, FORMULA

df = world(hidden=1.0).to_pandas()
y_fit = smf.ols(FORMULA, df).fit()
d_fit = smf.ols(FORMULA.replace("spend ~ joined +", "joined ~"), df).fit()
est, se, dof = y_fit.params["joined"], y_fit.bse["joined"], y_fit.df_resid
r2yx = partial_r2(y_fit.tvalues["engaged"], dof)
r2dx = partial_r2(d_fit.tvalues["engaged"], d_fit.df_resid)
print(f"Engagement: partial R2 {r2yx:.1%} with spend, {r2dx:.1%} with joining")

rows = []
for k in [1, 2, 3]:
    r2yz, r2dz = benchmark(r2yx, r2dx, k)
    adj = est - bias(se, dof, r2yz, r2dz)
    rows.append({"hidden confounder": f"{k}x as strong as engagement", "r2 with spend": round(r2yz, 3),
                 "r2 with joining": round(r2dz, 3), "adjusted estimate": round(adj, 1),
                 "still above $150?": "yes" if adj > BREAK_EVEN else "no"})
long_y = smf.ols(FORMULA + " + _affinity", df).fit()
long_d = smf.ols(FORMULA.replace("spend ~ joined +", "joined ~") + " + _affinity", df).fit()
rows.append({"hidden confounder": "the real brand affinity", "r2 with spend": round(partial_r2(long_y.tvalues['_affinity'], long_y.df_resid), 3),
             "r2 with joining": round(partial_r2(long_d.tvalues['_affinity'], long_d.df_resid), 3),
             "adjusted estimate": round(long_y.params["joined"], 1),
             "still above $150?": "yes" if long_y.params["joined"] > BREAK_EVEN else "no"})
with pl.Config(tbl_width_chars=120, fmt_str_lengths=40):
    print(pl.DataFrame(rows))
print(f"Unadjusted OLS: {est:.1f}")

python -m day62.benchmarks

Engagement: partial R2 11.6% with spend, 2.8% with joining
│ hidden confounder          ┆ r2 with spend ┆ r2 with joining ┆ adjusted estimate ┆ still above $150? │
│ 1x as strong as engagement ┆ 0.139         ┆ 0.029           ┆ 205.1             ┆ yes               │
│ 2x as strong as engagement ┆ 0.278         ┆ 0.058           ┆ 183.4             ┆ yes               │
│ 3x as strong as engagement ┆ 0.418         ┆ 0.086           ┆ 161.1             ┆ yes               │
│ the real brand affinity    ┆ 0.319         ┆ 0.132           ┆ 153.9             ┆ yes               │
Unadjusted OLS: 226.1

The statement for the CFO: "Even a hidden factor three times as strong as email engagement would leave the effect above break-even ($161)." And that was true: the real brand affinity is a bit over twice engagement's strength with spend and about four and a half times with joining. Adjusting for it gives $153.9, above $150 by a whisker.

Notice what the benchmarks don't say: the estimate moved from $226 to $154. The conclusion (above $150) survived; the number did not. If the decision depended on the size of the effect (how much to invest, say), the sensitivity analysis would be a warning, not a reassurance.

62.4 E-values for yes/no outcomes (50 min)

Concept. The E-value (VanderWeele & Ding, 2017) is the simplest sensitivity statistic of all, for risk ratios: the minimum strength of association, on the risk-ratio scale, that a hidden confounder would need with both treatment and outcome to explain the estimate away.

E-value = RR + √(RR × (RR − 1))

No modelling of the confounder at all, which makes it popular in medicine. Here the yes/no outcome is "high spender next year" (over $1,000). Create day62/e_value.py (about 30 seconds; the bootstrap refits the model 400 times):

# day62/e_value.py   (python -m day62.e_value)
"""
E-values (VanderWeele & Ding, 2017) for a yes/no outcome: does joining make a customer more likely to be a
'high spender' (over $1,000 next year)? The E-value is the risk ratio a hidden confounder would need with
BOTH joining and high spending to explain the observed risk ratio away. No model of the confounder needed.
The risk ratio is standardised (g-computation, Day 31) from a logistic model with squares and interactions;
its 95% CI comes from 200 bootstrap samples.
"""
import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

from causalkit.loyalty import world
from causalkit.sensitivity import e_value

COLS = ["joined", "pre_spend", "pre_orders", "engaged", "tenure_years", "age", "mobile"]


def standardised_rr(df):
    """Logistic model with all pairwise interactions and squares; predict everyone joined vs not joined."""
    model = make_pipeline(StandardScaler(), PolynomialFeatures(2), LogisticRegression(C=1e4, max_iter=2000))
    model.fit(df[COLS], df.high)
    p1 = model.predict_proba(df[COLS].assign(joined=1))[:, 1].mean()
    p0 = model.predict_proba(df[COLS].assign(joined=0))[:, 1].mean()
    return p1 / p0


for title, hidden in [("world A", 0.0), ("world B", 1.0)]:
    df = world(hidden=hidden).to_pandas().assign(high=lambda d: (d.spend > 1000).astype(int))
    rr = standardised_rr(df)
    rng = np.random.default_rng(0)
    boots = [standardised_rr(df.iloc[rng.integers(0, len(df), len(df))]) for _ in range(200)]
    lo = np.percentile(boots, 2.5)
    true_rr = (df._y1 > 1000).mean() / (df._y0 > 1000).mean()
    print(f"{title}: risk ratio {rr:.2f} (95% CI from {lo:.2f}); true risk ratio {true_rr:.2f}")
    print(f"   E-value: {e_value(rr):.2f} for the estimate, {e_value(lo):.2f} for the CI's lower end")

python -m day62.e_value

world A: risk ratio 1.36 (95% CI from 1.33); true risk ratio 1.38
   E-value: 2.05 for the estimate, 2.00 for the CI's lower end
world B: risk ratio 1.59 (95% CI from 1.55); true risk ratio 1.36
   E-value: 2.56 for the estimate, 2.47 for the CI's lower end

In world B: "a hidden confounder would need to make customers about 2.5 times as likely to join, and 2.5 times as likely to be high spenders, to explain away the effect". Report the E-value for the CI's lower end too: it tells you what it takes to make the effect statistically indistinguishable from zero.

E-values are easy to compute and to misread. A large E-value doesn't mean confounding is absent; and because E-values only address "explain away entirely", they can't answer threshold questions like break-even. Use them for yes/no outcomes in a world that thinks in risk ratios; use robustness values and benchmarks for decisions.

Try it. In world B, what strength would a confounder need to halve the effect, or cut it by a quarter? Use robustness_value(t, dof, q). Create day62/halving.py.

# day62/halving.py   (python -m day62.halving)
import numpy as np  # noqa: F401  (used inside the formula: np.sqrt)
import statsmodels.formula.api as smf

from causalkit.loyalty import world
from causalkit.sensitivity import robustness_value
from day62.robustness_value import FORMULA

fit = smf.ols(FORMULA, world(hidden=1.0).to_pandas()).fit()
t, dof = fit.tvalues["joined"], fit.df_resid
for q in [1.0, 0.5, 0.25]:
    print(f"RV to reduce the estimate by {q:.0%}: {robustness_value(t, dof, q):.1%}")

RV to reduce the estimate by 100%: 49.3%
RV to reduce the estimate by 50%: 29.1%
RV to reduce the estimate by 25%: 15.9%

62.5 Day 62 mini-project: the sensitivity contour plot (65 min)

Cinelli & Hazlett's signature figure: adjusted estimates for every hypothetical confounder, with benchmarks marked. Create day62/contour.py:

# day62/contour.py   (python -m day62.contour)
"""
Mini-project: the sensitivity contour plot (Cinelli & Hazlett's signature figure). Every point is a
hypothetical confounder (its strength with joining, its strength with spending); the contour lines show the
adjusted estimate. Mark the break-even line, the engagement benchmarks, and the real affinity.
"""
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.formula.api as smf

from causalkit.loyalty import world
from causalkit.sensitivity import benchmark, bias, partial_r2
from day62.robustness_value import BREAK_EVEN, FORMULA

df = world(hidden=1.0).to_pandas()
y_fit = smf.ols(FORMULA, df).fit()
d_fit = smf.ols(FORMULA.replace("spend ~ joined +", "joined ~"), df).fit()
est, se, dof = y_fit.params["joined"], y_fit.bse["joined"], y_fit.df_resid

r2d = np.linspace(0.001, 0.3, 120)
r2y = np.linspace(0.001, 0.6, 120)
D, Y = np.meshgrid(r2d, r2y)
adjusted = est - bias(se, dof, Y, D)

fig, ax = plt.subplots(figsize=(6.5, 5.5))
cs = ax.contour(D, Y, adjusted, levels=[0, 50, 100, 200], colors="#94a3b8")
ax.clabel(cs, fmt="%.0f")
be = ax.contour(D, Y, adjusted, levels=[BREAK_EVEN], colors="#dc2626", linewidths=2)
ax.clabel(be, fmt="break-even %.0f")
r2yx = partial_r2(y_fit.tvalues["engaged"], dof)
r2dx = partial_r2(d_fit.tvalues["engaged"], d_fit.df_resid)
for k in [1, 2, 3]:
    ry, rd = benchmark(r2yx, r2dx, k)
    ax.plot(rd, ry, "s", color="#2563eb")
    ax.annotate(f"{k}x engagement", (rd, ry), textcoords="offset points", xytext=(6, -10), fontsize=8)
long_y = smf.ols(FORMULA + " + _affinity", df).fit()
long_d = smf.ols(FORMULA.replace("spend ~ joined +", "joined ~") + " + _affinity", df).fit()
ta = (partial_r2(long_d.tvalues["_affinity"], long_d.df_resid), partial_r2(long_y.tvalues["_affinity"], long_y.df_resid))
ax.plot(*ta, "*", ms=14, color="#d97706")
ax.annotate("real brand affinity", ta, textcoords="offset points", xytext=(6, 6), fontsize=8)
ax.set_xlabel("partial R2 of the confounder with joining")
ax.set_ylabel("partial R2 of the confounder with spending")
ax.set_title(f"Adjusted effect of joining (unadjusted {est:.0f})")
fig.tight_layout()
fig.savefig("figures/day62_contour.png", dpi=120)
print("Saved figures/day62_contour.png")
print(f"Adjusted estimate at the real affinity's strength (worst case): {est - bias(se, dof, ta[1], ta[0]):.1f}")

python -m day62.contour

Saved figures/day62_contour.png
Adjusted estimate at the real affinity's strength (worst case): 153.9

Open the figure. The red line is break-even; the blue squares are 1×, 2× and 3× engagement; the star is the real affinity, just on the profitable side of the red line. For an OLS estimate the bias formula is exact, so the worst-case adjustment at the affinity's strength reproduces the long regression ($153.9) to the decimal.

Put this figure in any report that rests on "we controlled for everything important". It shows the reader exactly how much "everything" has to cover.

62.6 Save your work with Git (20 min)

git add .
git commit -m "Day 62: partial R2, robustness values, benchmarks, E-values, contour plot"
git push

Day 62 checklist

• I can compute a robustness value from a t-statistic, for zero and for a decision threshold
• I can benchmark a hidden confounder as k times an observed covariate
• I can compute and interpret an E-value for a risk ratio, and its limits
• I distinguish "the conclusion survives" from "the number survives"
• I can draw a sensitivity contour plot with benchmarks