Sensitivity Analysis: How Robust Is Your Estimate to an Unmeasured Confounder?

Last Time on the 90-Day Policy

Maya — the new data scientist at Blue Harbor Hotels — just handed her first real result to Mr. Delgado, the VP of Revenue: for every extra day of lead_time (the time between booking and arrival), the probability of a guest cancelling rises by about 0.3% (a coefficient of 0.003). She built the DAG (Directed Acyclic Graph — a picture of what causes what), identified the estimand (the mathematical recipe for the effect), estimated her number, and even ran refutation tests — attempts to break her own result on purpose. The number survived.

Then Delgado asked the question that ended the meeting: “What if guests who book far in advance are just ‘flaky’ people who change their minds often — and we don’t have a ‘flakiness’ column in our database?”

That’s the Ghost in the Machine: the unmeasured confounder. It’s the variable you never collected, the one that might sit behind both your treatment (lead time) and your outcome (cancellation). If a ghost like that is real, Maya’s estimate isn’t just slightly off. It could be completely wrong.

So Maya learns to stress test her model — to check whether hidden variables can undo her conclusions. You’ll learn the same.

What you’ll learn here:

• What an unmeasured confounder is and why it’s so dangerous
• What sensitivity analysis actually does (and what it can’t do)
• How to run a hidden-confounder stress test with DoWhy — with a plot built from real refuter output, not a mocked-up line
• How to read the tipping point
• What the E-value is, computed correctly on a risk-ratio scale, and how to use it to talk to busy stakeholders

Prerequisites: Tutorials 1–3, or comfort with DoWhy’s CausalModel, estimands, and refutation. Every new term is defined when it first appears.

---

1. The Ghost in the Machine: The Unmeasured Confounder

In Part 3, Maya built a model assuming she saw everything important. Let’s rebuild that baseline together — but this time, honestly acknowledge what she might be missing.

Note on the data: This tutorial reuses the exact canonical dataset from Part 2/3 — 3,000 rows, np.random.seed(0), the same lead_time/is_repeat_guest/market_segment/special_requests/adr/is_cancelled columns, the same 0.003/day true lead_time effect. We are not rebuilding it from scratch with a different seed or a different effect size. We extend it with exactly one new column: loyalty_score, the hidden “Loyalty/Affinity” factor this whole tutorial is about.

A loyal guest books early (short lead time) and is much less likely to cancel. If we don’t measure “Loyalty,” our model might wrongly attribute the lower cancellations to lead time itself — when really the kind of person who books early is also the kind of person who stays. For that story to be worth stress-testing, loyalty_score has to actually pull on both lead_time and is_cancelled — the same shape as Part 2’s is_repeat_guest → lead_time arrow. So we wire it in as a real cause, not a bystander column that happens to sit in the same dataframe.

A quick vocabulary check, since these words carry the whole tutorial:

• Unmeasured confounder: a confounder (a variable that affects both treatment and outcome) that you didn’t collect in your data. You can’t see it. But it’s there.
• Bias: a systematic error in your estimate — a reason your number is wrong in a consistent direction, not just random sloppiness.

import pandas as pd
import numpy as np
from dowhy import CausalModel

np.random.seed(0)
n = 3000

is_repeat_guest = np.random.binomial(1, 0.3, n)
market_segment = np.random.choice(
    ["Direct", "Corporate", "Online TA", "Groups"], n,
    p=[0.2, 0.25, 0.45, 0.1]
)

# NEW in this tutorial: loyalty_score, the "ghost" confounder. Drawn right
# after the other guest-level traits, before lead_time, so it can feed
# into lead_time below -- loyal guests book earlier AND cancel less, which
# is what makes it a real confounder instead of an inert extra column.
loyalty_score = np.random.normal(0, 1, n)

# lead_time: same Part 2/3 shape (segment + is_repeat_guest), now also
# pulled shorter by loyalty_score.
segment_lead_scale = np.select(
    [market_segment == "Corporate", market_segment == "Groups"],
    [35, 90],
    default=60,
)
lead_time_scale = np.clip(
    segment_lead_scale - 45 * is_repeat_guest - 10 * loyalty_score, 10, None
)
lead_time = np.random.exponential(scale=lead_time_scale).astype(int)

# special_requests: unchanged mediator from Part 2 (lead_time -> special_requests -> is_cancelled).
segment_request_base = np.select(
    [market_segment == "Corporate", market_segment == "Groups"],
    [0.3, 1.1],
    default=0.7,
)
special_requests_rate = segment_request_base + 0.01 * lead_time
special_requests = np.random.poisson(special_requests_rate, n)

# adr: unchanged collider from Part 2. Still excluded from every graph below.
adr_base = np.where(market_segment == "Corporate", 120, 95)
adr = adr_base + lead_time * 0.1 + np.random.normal(0, 15, n)

# cancel_prob: same Part 2/3 structure, with loyalty_score's direct pull
# added on top. 0.003 is the TRUE per-day lead_time effect -- the same
# number Part 2 and Part 3 use. Nothing about the treatment effect changed;
# we only added one more thing that ALSO moves the outcome.
segment_cancel_shift = np.select(
    [market_segment == "Corporate", market_segment == "Groups"],
    [-0.05, 0.08],
    default=0.0,
)
cancel_prob = (
    0.2
    + 0.003 * lead_time
    - 0.3 * is_repeat_guest
    - 0.03 * special_requests
    - 0.05 * loyalty_score
    + segment_cancel_shift
    + np.random.normal(0, 0.08, n)
)
is_cancelled = np.random.binomial(1, np.clip(cancel_prob, 0, 1))

df_full = pd.DataFrame({
    "lead_time": lead_time,
    "is_repeat_guest": is_repeat_guest,
    "market_segment": market_segment,
    "special_requests": special_requests,
    "adr": adr.round(2),
    "loyalty_score": loyalty_score,   # kept ONLY so we can check our work later
    "is_cancelled": is_cancelled,
})

# The dataframe Maya actually has. loyalty_score does not appear in it --
# exactly like the real world, where nobody collected it.
df = df_full.drop(columns=["loyalty_score"])

print(f"Cancellation rate: {df['is_cancelled'].mean():.1%}")
print(f"loyalty_score vs lead_time correlation:  {df_full['loyalty_score'].corr(df_full['lead_time']):.4f}")
print(f"loyalty_score vs is_cancelled correlation: {df_full['loyalty_score'].corr(df_full['is_cancelled']):.4f}")

Running this prints:

Cancellation rate: 22.6%
loyalty_score vs lead_time correlation:  -0.1455
loyalty_score vs is_cancelled correlation: -0.1413

Both correlations are real and non-trivial — loyalty_score genuinely pulls on both lead_time (negatively: loyal guests book closer in) and is_cancelled (negatively: loyal guests cancel less). That’s what makes it worth worrying about. A “confounder” that happened to be statistically independent of the treatment wouldn’t confound anything — it’d just be noise sitting in the dataframe.

# The model Maya actually has: is_repeat_guest, market_segment, and
# special_requests are in her graph (same as Part 2/3). loyalty_score is
# NOT -- she doesn't know it exists.
model_biased = CausalModel(
    data=df,
    treatment="lead_time",
    outcome="is_cancelled",
    common_causes=["is_repeat_guest", "market_segment", "special_requests"],
)
estimand = model_biased.identify_effect(proceed_when_unidentifiable=True)
estimate = model_biased.estimate_effect(estimand, method_name="backdoor.linear_regression")
print(f"Biased estimate: {estimate.value:.6f}")

# The oracle model: same data, but WITH loyalty_score included -- only
# possible here because we're the ones who generated the ghost. This is
# not a model Maya could ever run for real; it's how we, the tutorial
# authors, check that the ghost actually did something.
model_oracle = CausalModel(
    data=df_full,
    treatment="lead_time",
    outcome="is_cancelled",
    common_causes=["is_repeat_guest", "market_segment", "special_requests", "loyalty_score"],
)
estimand_o = model_oracle.identify_effect(proceed_when_unidentifiable=True)
estimate_o = model_oracle.estimate_effect(estimand_o, method_name="backdoor.linear_regression")
print(f"Oracle estimate: {estimate_o.value:.6f}")

Biased estimate: 0.002638
Oracle estimate: 0.002509

What this actually means: Maya’s real-world model — blind to loyalty_score, exactly as she’d be in production — puts the effect at 0.002638. The oracle model, which cheats by including the ghost, puts it at 0.002509. The gap (0.000128) is small but real: omitting loyalty_score overstates the per-day effect by about 5.1%. That’s the mechanical proof this dataset’s ghost isn’t decorative — leaving it out of the graph measurably moves the number, the same way Part 2’s Scenario A (missing is_repeat_guest) did. It’s a smaller bias than Part 2’s ~25% swing, because loyalty_score’s pull, while real, is more modest than is_repeat_guest’s — which is exactly the kind of thing sensitivity analysis is built to quantify precisely, instead of eyeballing.

Both numbers sit close to the true 0.003 — the estimation procedure (linear regression on a finite, noisy sample) doesn’t recover the exact structural coefficient even in the oracle case, the same pattern Part 2 showed. What sensitivity analysis cares about isn’t hitting 0.003 exactly; it’s the gap between “what I can measure” and “what I’m missing” — and here, that gap is real, small, and about to get measured properly.

2. What is Sensitivity Analysis Actually Doing?

Think of sensitivity analysis as a stress test for a bridge. You don’t know if a category-5 hurricane will ever hit, but you calculate how much wind the structure can take before it fails. That number is what lets you sleep at night.

We aren’t proving the hidden variable exists. We’re playing a game of “what if.” If a hidden factor explained some share of the leftover variation in cancellations, would Maya’s result still hold?

Here’s the hard part: you can’t fix the bias when the data is missing. You can only measure your vulnerability to it. What we’re after is the tipping point — the moment where the hidden variable grows strong enough that the estimated effect vanishes or reverses.

That’s a humble position. Maya can’t prove the ghost isn’t there. She can, however, measure exactly how big it would have to be to eat her result — and if that size is unrealistic (or, as we’ll see, bigger than the ghost we actually know is sitting in this data), she can defend her number in front of Delgado.

3. Running the Stress Test with DoWhy

DoWhy handles this with the refute_estimate method from Part 3. The specific refuter here is add_unobserved_common_cause. It simulates a “ghost” variable with a given strength and re-estimates the effect with that ghost partialled in.

# Run the sensitivity analysis over the same effect-strength grid on both axes
refute = model_biased.refute_estimate(
    estimand,
    estimate,
    method_name="add_unobserved_common_cause",
    confounders_effect_on_treatment="linear",
    confounders_effect_on_outcome="linear",
    effect_strength_on_treatment=[0.01, 0.05, 0.1],
    effect_strength_on_outcome=[0.01, 0.05, 0.1],
)
print(refute)

Two parameters control this test:

• effect_strength_on_treatment — how hard the ghost variable pushes on lead_time, in raw regression-coefficient units.
• effect_strength_on_outcome — how hard the ghost variable pushes on is_cancelled, same units.

Running the code above prints, among other things, a New effect range for the 3×3 grid of strengths:

New effect:(0.0026078259413437194, 0.002667985423156452)

At these (small) strengths, the estimate barely moves — it stays inside 0.00261–0.00267, close to the biased 0.002638 baseline, and nowhere near zero. That’s a genuinely reassuring — and genuinely computed — result. But notice something about the range’s own width: 0.01 to 0.1 is a fairly narrow slice of “how strong could a ghost be.” It tells us the model survives small ghosts. It doesn’t yet tell us the actual breaking point. For that, we need a wider sweep and a method that doesn’t depend on redrawing a random ghost each time.

4. Reading the ‘Tipping Point’ Plot

A full sensitivity suite in the literature ends with a contour plot — a 2-D chart where the X-axis is how strongly the ghost pushes on the treatment, the Y-axis is how strongly it pushes on the outcome, and the shading traces “equal-estimate” contours across the two strengths. DoWhy ships a deterministic version of exactly this, based on Cinelli & Hazlett’s partial-R² sensitivity framework (simulated_method_name="linear-partial-R2") — instead of resimulating a random confounder each call (which is what Section 3’s plot would do if you swept it over many points, and why that sweep would come out noisy and non-monotonic), it computes the bias analytically as a function of how much residual variance the hypothetical confounder explains in the treatment and the outcome. That’s what lets us actually plot real numbers instead of a schematic.

We benchmark the hypothetical ghost’s strength against is_repeat_guest — our strongest known confounder — and ask: “what if the ghost were 1×, 2×, 3×… as strongly associated with lead time and cancellation (in partial-R² terms) as is_repeat_guest already is?”

import statsmodels.api as sm
import matplotlib.pyplot as plt

# The deterministic partial-R2 method needs a fully numeric design matrix,
# so we one-hot encode market_segment before running it.
df_num = pd.concat(
    [df.drop(columns=["market_segment"]),
     pd.get_dummies(df["market_segment"], prefix="seg", drop_first=True).astype(float)],
    axis=1,
)
seg_cols = [c for c in df_num.columns if c.startswith("seg_")]
common_causes = ["is_repeat_guest", "special_requests"] + seg_cols

model_n = CausalModel(data=df_num, treatment="lead_time", outcome="is_cancelled", common_causes=common_causes)
estimand_n = model_n.identify_effect(proceed_when_unidentifiable=True)
estimate_n = model_n.estimate_effect(estimand_n, method_name="backdoor.linear_regression")

tipping = model_n.refute_estimate(
    estimand_n, estimate_n,
    method_name="add_unobserved_common_cause",
    simulated_method_name="linear-partial-R2",
    benchmark_common_causes=["is_repeat_guest"],
    effect_fraction_on_treatment=[1, 2, 3, 4, 5, 6],
    effect_fraction_on_outcome=[1, 2, 3, 4, 5, 6],
)
print(tipping.benchmarking_results[["r2tu_w", "r2yu_tw", "bias_adjusted_estimate"]])

multiplier = [1, 2, 3, 4, 5, 6]
new_effects = tipping.benchmarking_results["bias_adjusted_estimate"].tolist()
plt.plot(multiplier, new_effects, marker='o')
plt.axhline(0, color='red', linestyle='--')
plt.xlabel("Ghost strength, as a multiple of is_repeat_guest's own strength")
plt.ylabel("Bias-adjusted causal effect")
plt.title("Real Tipping-Point Sweep (linear-partial-R2)")
plt.show()

This is the real returned data — not a straight line produced by subtraction:

   multiplier  r2tu_w   r2yu_tw   bias_adjusted_estimate
            1  0.0986   0.0721    0.001900
            2  0.1972   0.1473    0.001058
            3  0.2957   0.2268    0.000074
            4  0.3943   0.3126   -0.001110
            5  0.4929   0.4081   -0.002594
            6  0.5915   0.5196   -0.004567

How to read this plot, using these actual numbers:

• Each dot is a real re-estimated effect, computed by DoWhy’s sensitivity analyzer at that ghost strength — not estimate.value - s for some made-up s.
• The red dashed line is zero. The curve crosses it between the 3× row (0.000074 — a hair’s breadth above zero) and the 4× row (-0.001110, already negative). Interpolating, the tipping point lands at roughly 3.06× is_repeat_guest’s partial-R² strength.
• The plot also reports a single summary number, DoWhy’s robustness value: 27.11%. That’s the share of residual variance a confounder would need to explain in both lead_time and is_cancelled, symmetrically, to fully explain away the effect — a different (stricter, benchmark-free) framing of the same tipping point. At the 5% significance level, a confounder explaining more than 24.48% of residual variance on both sides would already make the estimate statistically indistinguishable from zero.

Does the sensitivity analysis actually catch our real ghost? We know the true answer here, because we built loyalty_score ourselves. Its actual partial R² — computed the same way, controlling for the same observed covariates — is:

X_controls = sm.add_constant(df_num[common_causes].astype(float))
resid_treat = sm.OLS(df_num["lead_time"].astype(float), X_controls).fit().resid
resid_loy_t = sm.OLS(df_full["loyalty_score"], X_controls).fit().resid
r2tu_actual = np.corrcoef(resid_treat, resid_loy_t)[0, 1] ** 2

X_controls_y = sm.add_constant(df_num[["lead_time"] + common_causes].astype(float))
resid_out = sm.OLS(df_num["is_cancelled"].astype(float), X_controls_y).fit().resid
resid_loy_y = sm.OLS(df_full["loyalty_score"], X_controls_y).fit().resid
r2yu_actual = np.corrcoef(resid_out, resid_loy_y)[0, 1] ** 2

print(f"loyalty_score partial R2 with treatment: {r2tu_actual:.4f}")
print(f"loyalty_score partial R2 with outcome:   {r2yu_actual:.4f}")

loyalty_score partial R2 with treatment: 0.0203
loyalty_score partial R2 with outcome:   0.0115

0.0203 and 0.0115 — only about 20% and 16% of is_repeat_guest’s own strength on each axis (row “1×” above is 0.0986/0.0721), and nowhere close to the 27.11% robustness value needed to erase the effect on either axis. This is the real, numeric confirmation the earlier finding demanded: the ghost we built genuinely confounds (Section 1’s correlations, and the 5.1% bias from Section 1’s biased-vs-oracle comparison, prove that), and the sensitivity analysis correctly places it well inside “survivable” territory rather than at the tipping point. A stronger ghost — one actually comparable to is_repeat_guest — would be a different story, and that’s exactly what the sweep above shows happening between 2× and 4×.

5. The ‘E-Value’: A Shortcut for the Busy Data Scientist

Sometimes your boss doesn’t want a tipping-point sweep. They want one number.

Here’s what it is: the E-value is the minimum risk-ratio strength a hidden confounder would need — with both the treatment and the outcome — to fully explain away an observed effect. It’s defined on a ratio scale, which matters: you cannot hand it a raw probability difference like 0.003 and get a sensible answer. It has to be built from a risk ratio.

So, to compute a real E-value for Maya’s result, we first have to turn her linear per-day effect into a risk ratio. We use the same threshold the whole series has been built around: Blue Harbor’s 90-day booking-fee policy. Compare the predicted cancellation probability at lead_time = 90 against lead_time = 0, using the biased estimate’s slope:

p0 = df["is_cancelled"].mean()
beta = estimate.value          # 0.002638, the biased per-day estimate
delta_days = 90                # the tutorial's own 90-day policy threshold

p1 = p0 + beta * delta_days
RR = p1 / p0

def e_value(rr):
    if rr < 1:
        rr = 1 / rr
    return rr + np.sqrt(rr * (rr - 1))

E = e_value(RR)
print(f"p0={p0:.4f}  p1={p1:.4f}  RR={RR:.4f}  E-value={E:.4f}")

p0=0.2260  p1=0.4634  RR=2.0504  E-value=3.5179

Booking 90+ days out roughly doubles the predicted cancellation probability relative to booking with no lead time at all (RR = 2.05), on top of everything is_repeat_guest, market_segment, and special_requests already explain. Plugging that risk ratio into the standard E-value formula (E = RR + √(RR·(RR−1))) gives E = 3.52.

What does 3.52 mean? An unmeasured confounder would need to be associated with both “booked 90+ days out” and “cancelled” by a risk ratio of roughly 3.5×ish each — over and above everything already in the graph — to fully explain away the observed policy-threshold effect. Running the same calculation on the oracle model (which already includes loyalty_score) gives RR = 1.9993, E = 3.4127 — nearly identical, because, as Section 4 just showed numerically, loyalty_score’s real strength is well below what it would take to move this number much.

• E-value near 1: Fragile. A weak hidden factor could undo your result.
• E-value of 3.52, like Maya’s: Moderately robust. Meaningfully sturdier than a toy “1.2,” well short of an extreme “5.0” — a real number in between, not a hand-picked bookend.
• E-value of 5+: Robust. A hidden factor would have to be dramatically stronger than your measured variables to flip the result.

The intuition comes down to this: how tough does the ghost have to be? A ghost that only needs to be 1.2× as strong as your best measured variable is easy to picture — and that’s bad. 3.5× is a genuinely demanding bar to clear. 5×+ is hard to picture — and that’s good.

6. Conclusion: When to Trust Your Data

Sensitivity analysis doesn’t find the truth. It finds the breaking point. Before Maya moves this hotel model into production, she runs her three-question checklist:

1. Is the estimate robust? The default-strength add_unobserved_common_cause sweep (Section 3) barely moved the number (0.00261–0.00267). The wider, deterministic tipping-point sweep (Section 4) shows the estimate doesn’t cross zero until a hypothetical ghost reaches roughly 3× is_repeat_guest’s strength — or, in the robustness-value framing, until it explains 27.11% of residual variance on both sides. Neither number is a hair-trigger.
2. Is the tipping point realistic? Here’s the part that used to be hand-waved: we actually know the tipping point isn’t close, because we know the real ghost. loyalty_score’s partial R² (0.0203 on treatment, 0.0115 on outcome) sits at only 16–20% of is_repeat_guest’s own strength — nowhere near the 27.11% robustness value. If Maya had a real, uncollected “loyalty” signal even this strong, her lead-time effect would still hold. A genuinely fragile result — a robustness value in the single digits — would be a to-do list, not a conclusion.
3. Can I explain this? Maya can now tell Delgado, precisely: “Even a hidden factor three times as strongly tied to booking-and-cancelling as repeat-guest status wouldn’t erase this effect. And the E-value — 3.5 — says a hidden confounder would need a risk ratio north of 3× with both lead time and cancellation to undo it. The one candidate ghost we can name, guest loyalty, isn’t nearly that strong.” That’s real trust — the kind that survives a board meeting, because every number in it came out of code that actually ran.

Maya’s result survives the stress test. Not with a vague asterisk — with a specific number attached to exactly how big an unknown ghost would have to be, and a specific comparison showing the one ghost she can name falls well short of that bar.

Then Delgado pushes once more. “Maya, what if the ghost isn’t a maybe? What if we know there’s an unmeasured factor — like a guest’s stress level — and we simply can’t measure it, no matter what? Your stress test tells me how big it would have to be. It doesn’t tell me the answer anyway.”

He’s right. That’s the problem the next tutorial solves — with a method called an Instrumental Variable. A side door when the front one is locked.

In this tutorial, you learned:

• What an unmeasured confounder is and why it threatens every causal claim
• Why sensitivity analysis measures your vulnerability, not the truth
• How to run add_unobserved_common_cause in DoWhy, and why its default random-simulation grid is too narrow to show a real tipping point
• How to build a tipping-point curve from DoWhy’s deterministic linear-partial-R2 sensitivity analyzer, benchmarked against a known confounder — real numbers, not a subtraction line
• How to compute a proper E-value on the risk-ratio scale it’s actually defined on, and use it to communicate robustness to a stakeholder

Next up: Tutorial 5 — Instrumental Variables: when you can’t measure the confounder directly.

Check Your Understanding

Questions below move from simple recall up to open-ended design, roughly following Bloom’s Taxonomy.

Remember What does the “E-value” measure, and why does it need to be computed from a risk ratio rather than applied directly to a raw probability-scale effect like 0.003?

Understand In your own words, explain why sensitivity analysis “doesn’t find the truth; it finds the breaking point”—what can it tell you, and what can it never tell you about whether a hidden confounder actually exists?

Apply Using the tutorial’s numbers: the biased estimate is 0.002638, and the tipping-point sweep shows it takes a hidden confounder with roughly 3× the partial-R² strength of is_repeat_guest to push the effect to zero. Is that closer to the “very fragile” or “very robust” end of the spectrum the tutorial describes? What would make you less confident in that answer?

Analyze The tutorial builds a synthetic loyalty_score that it deliberately excludes from the model to demonstrate the effect of a hidden confounder — and then, uniquely, goes back and checks its actual partial R² against the tipping-point threshold. Walk through why this “secretly knowing the true confounder” setup is a good teaching tool but can’t tell us what to do in a real analysis where we don’t know what we’re missing.

Evaluate Section 6’s checklist item 3 says Maya has “built real trust” if she can tell Delgado her effect survives a confounder “three times as strongly tied to booking-and-cancelling as repeat-guest status.” Critique this benchmark: what makes “3× as strong as a specific, named observed confounder” a more defensible bar than an arbitrary-sounding round number like “an E-value of 2”?

Create Design a sensitivity-analysis stress test for a different claim: “customers who read at least one blog post before signing up churn less.” Name a plausible unmeasured confounder for this claim, and describe what observed variable you’d benchmark it against and what effect_fraction_on_treatment / effect_fraction_on_outcome multipliers you’d want to test to see if the claim survives.

---

Each card below turns this article into something you build, not just read. Pick the card that matches your role (or do more than one). Every 🔬 and 🛠 card ships a stub file in projects/causal-p04/ — clone it, fill in the # TODO markers, and grade yourself against the rubric.