Cross-Sectional Application: Price Elasticity with Sensitivity Analysis

Introduction

Chapters 1–3 established the theoretical foundations and validated DML on synthetic data where the true treatment effect is known. This chapter takes the critical step of applying DML to real observational data — where the true causal effect is unknown and unmeasured confounding is a genuine concern.

Why sensitivity analysis matters

The unconfoundedness assumption — that all confounders are observed and controlled — is fundamentally untestable. No matter how many covariates we include, there may exist unmeasured factors affecting both treatment and outcome. This creates a persistent concern: how robust are our findings to potential hidden bias?

Sensitivity analysis addresses this by asking: “How strong would unmeasured confounding need to be to qualitatively change our conclusions?” If even modest confounding could overturn our results, we should interpret them cautiously. If only implausibly strong confounding could matter, we gain confidence in our findings.

Chapter roadmap

• Section 4.2: the OJ dataset — exploratory analysis and preprocessing.
• Section 4.3: DML estimation — step-by-step price elasticity estimation.
• Section 4.4: baseline comparison — quantifying the value of DML.
• Section 4.5: overlap diagnostics — verifying positivity for continuous treatments.
• Section 4.6: sensitivity analysis — Rosenbaum bounds, E-values, and partial identification.
• Section 4.7: contrastive application — gasoline elasticity as a fragile counterexample.
• Section 4.8: heterogeneous sensitivity — brand-level and demographic robustness.
• Section 4.9: diagnostic visualization — overlap, sensitivity surfaces, and residual plots.
• Section 4.10: practical recommendations — when to trust DML estimates.
• Section 4.11: exercises.

The orange juice dataset

Data description

We analyze the Dominick’s Orange Juice dataset, a benchmark in causal inference and marketing econometrics. The data comes from Dominick’s Finer Foods, a major Chicago-area supermarket chain, and contains weekly store-level sales from 83 stores over 121 weeks.

Orange juice dataset variables

*We use $T = \log(\texttt{price})$ for elasticity interpretation.

The causal question

We seek to estimate the price elasticity of demand: the percentage change in quantity demanded for a 1% change in price. In the log–log specification:

where $\tau$ is the price elasticity (expected to be negative by the law of demand).

Loading the data

Our implementation provides a clean interface for loading the OJ dataset:

from dml_ts.data import OJDataLoader

# Load with default confounders
loader = OJDataLoader()
data = loader.load()

print(data.summary())
# Orange Juice Dataset Summary
# ============================
# Observations:     28,947
# Features:         3
# Feature names:    feat, INCOME, AGE60
#
# Outcome (Y = log sales):
#   Mean:           9.168
#   Range:          [4.159, 13.482]
#
# Treatment (T = log price):
#   Mean:           0.784
#   Range:          [-0.654, 1.353]

The loader handles downloading, caching, and preprocessing. Treatment $T$ is automatically transformed to log-price for elasticity interpretation.

Exploratory analysis

Before estimation, we examine the data for potential issues:

import numpy as np
import matplotlib.pyplot as plt

# Basic correlations
print("Correlation matrix:")
print(f"  Corr(Y, T) = {np.corrcoef(data.Y, data.T)[0,1]:.3f}")
print(f"  Corr(T, INCOME) = {np.corrcoef(data.T, data.X[:,1])[0,1]:.3f}")
print(f"  Corr(Y, INCOME) = {np.corrcoef(data.Y, data.X[:,1])[0,1]:.3f}")

# Output:
# Corr(Y, T) = -0.421  (negative, as expected)
# Corr(T, INCOME) = 0.182  (higher income areas have higher prices)
# Corr(Y, INCOME) = 0.089  (modest demand effect)

The correlations reveal a classic confounding pattern: higher-income areas tend to have both higher prices (stores charge more) and higher baseline demand. Naive regression of $Y$ on $T$ would underestimate the true price sensitivity because it conflates the negative price effect with the positive income effect.

DML price elasticity estimation

The five-step DML pipeline

We now apply the DML methodology developed in Chapter 2 to estimate price elasticity:

1. Define the problem: partially linear model with log-price treatment.
2. Choose nuisance models: Random Forest for flexible confounding control.
3. Set cross-fitting folds: 5-fold for bias–variance balance.
4. Estimate with DML: cross-fit nuisance, compute residualized effect.
5. Report results: point estimate, confidence interval, diagnostics.

Implementation

from dml_ts import double_ml

# Run DML estimation
result = double_ml(
    Y=data.Y,          # Log sales
    T=data.T,          # Log price
    X=data.X,          # Confounders [feat, INCOME, AGE60]
    n_folds=5,
    model="random_forest",
    random_state=42,
)

print(result.summary())

Double Machine Learning Results
================================
Treatment Effect (θ):    -2.8347
Standard Error:          0.0412
t-statistic:             -68.81
p-value:                 0.0000
95% Confidence Interval: [-2.9154, -2.7540]

Nuisance Model Diagnostics:
  Outcome R² (CV):       0.421
  Treatment R² (CV):     0.187
  Number of folds:       5

Interpretation:
  A 1% increase in price is associated with a
  2.83% decrease in quantity demanded (p<0.001).

First-stage diagnostics

The nuisance model $R^2$ values provide important diagnostics:

• Outcome $R^2 = 0.42$: confounders explain substantial outcome variation (good).
• Treatment $R^2 = 0.19$: moderate treatment prediction (sufficient variation remains).

If treatment $R^2$ were too high (e.g., above 0.9), it would indicate near-deterministic treatment assignment, violating the overlap assumption and leading to unstable estimates.

Residual quality checks

Beyond $R^2$ values, the residuals themselves encode diagnostic information. Well-behaved residuals — approximately mean-zero, symmetric, and uncorrelated — support the validity of the influence function standard errors:

from scipy import stats

Y_resid = result.Y_residual
T_resid = result.T_residual

# Mean should be ~0 (cross-fitting guarantees this asymptotically)
print("Residual Quality Checks")
print("=" * 45)
print(f"  Y residual mean:  {np.mean(Y_resid):.6f}")
print(f"  T residual mean:  {np.mean(T_resid):.6f}")

# Normality (Shapiro-Wilk on subsample for speed)
subsample = np.random.RandomState(42).choice(
    len(Y_resid), size=5000, replace=False
)
_, p_Y = stats.shapiro(Y_resid[subsample])
_, p_T = stats.shapiro(T_resid[subsample])
print(f"  Y Shapiro-Wilk p: {p_Y:.4f}")
print(f"  T Shapiro-Wilk p: {p_T:.4f}")

# Residual-residual correlation (should be ~theta)
resid_corr = np.corrcoef(Y_resid, T_resid)[0, 1]
print(f"  Corr(Y_resid, T_resid): {resid_corr:.4f}")

Baseline comparison

To quantify the value of DML’s confounding adjustment, we compare against OLS baselines.

Three estimators

from sklearn.linear_model import LinearRegression
import numpy as np

# 1. Naive OLS: Y ~ T (ignores confounders)
naive = LinearRegression()
naive.fit(data.T.reshape(-1, 1), data.Y)
naive_theta = naive.coef_[0]

# 2. OLS with Controls: Y ~ T + X
TX = np.column_stack([data.T, data.X])
controls = LinearRegression()
controls.fit(TX, data.Y)
controls_theta = controls.coef_[0]

# 3. DML (from above)
dml_theta = result.theta

print(f"Naive OLS:       θ = {naive_theta:.4f}")
print(f"OLS + Controls:  θ = {controls_theta:.4f}")
print(f"DML:             θ = {dml_theta:.4f}")

Estimator comparison for OJ price elasticity

Overlap diagnostics

Before assessing sensitivity to unmeasured confounding, we must verify a more basic requirement: overlap. Without adequate overlap, DML estimates become unstable regardless of how well we model confounders.

The overlap assumption

Why does overlap matter for DML specifically? Recall the DML estimator:

When $\tilde{T}_i \approx 0$ for many observations — meaning treatment is nearly determined by covariates — the denominator approaches zero and variance explodes. This is the continuous-treatment analogue of propensity scores near 0 or 1 in the binary case.

Propensity diagnostics for continuous treatment

For continuous treatments, we diagnose overlap through the treatment residual distribution and the effective sample size (ESS):

import numpy as np
from scipy import stats

# Treatment residuals from DML first stage
T_resid = result.T_residual

# Basic distribution diagnostics
print("Treatment Residual Diagnostics")
print("=" * 40)
print(f"  Mean:     {np.mean(T_resid):.4f}")
print(f"  Std Dev:  {np.std(T_resid):.4f}")
print(f"  Skewness: {stats.skew(T_resid):.3f}")
print(f"  Kurtosis: {stats.kurtosis(T_resid):.3f}")

# Effective sample size: Var(T_resid) / Var(T)
ess_ratio = np.var(T_resid) / np.var(data.T)
ess = ess_ratio * len(data.T)
print(f"\n  ESS ratio:  {ess_ratio:.3f}")
print(f"  ESS:        {ess:.0f} / {len(data.T)}")

Treatment Residual Diagnostics
========================================
  Mean:     0.0001
  Std Dev:  0.3842
  Skewness: -0.127
  Kurtosis: 0.412

  ESS ratio:  0.813
  ESS:        23,520 / 28,947

An ESS ratio of 0.81 means that 81% of the treatment variation remains after removing confounders — a healthy signal. When this ratio drops below 0.2, estimates become unreliable.

Near-violation detection

Even with good aggregate overlap, local violations can destabilize estimates. We borrow a production diagnostic from Chapter 10’s CausalMonitor to check for extreme treatment residuals:

from dml_ts.production.causal_monitor import CausalMonitor

monitor = CausalMonitor()

# Binarize treatment at median for propensity-style check
T_binary = (data.T > np.median(data.T)).astype(float)

violations = monitor.check_overlap_violations(
    propensity_scores=T_binary.mean() * np.ones(len(T_binary)),
    threshold=0.05,
)

print(f"Overlap violations: {violations['n_violations']} / {len(data.T)}")
print(f"Violation rate: {violations['violation_rate']:.4f}")

Trimming strategies

When overlap violations occur, three strategies can stabilize estimates:

Overlap trimming strategies

# Percentile trimming: remove observations with |T_resid| < 5th percentile
threshold = np.percentile(np.abs(T_resid), 5)
keep_mask = np.abs(T_resid) > threshold
n_trimmed = (~keep_mask).sum()

print(f"Trimming threshold: |T_resid| < {threshold:.4f}")
print(f"Observations trimmed: {n_trimmed} ({n_trimmed/len(T_resid)*100:.1f}%)")

# Re-estimate on trimmed sample
result_trimmed = double_ml(
    Y=data.Y[keep_mask],
    T=data.T[keep_mask],
    X=data.X[keep_mask],
    n_folds=5,
    model="random_forest",
    random_state=42,
)

print(f"\nOriginal:  θ = {result.theta:.4f} (SE = {result.se:.4f})")
print(f"Trimmed:   θ = {result_trimmed.theta:.4f} "
      f"(SE = {result_trimmed.se:.4f})")

OJ overlap summary

Sensitivity analysis: Rosenbaum bounds

The unconfoundedness problem

Even with DML controlling for observed confounders $X$, we cannot rule out unmeasured confounders $U$ that affect both treatment and outcome. The unconfoundedness assumption

is fundamentally untestable with observational data.

Rosenbaum’s sensitivity framework

Rosenbaum (2002) introduced a principled approach to this problem. The key idea: instead of testing unconfoundedness, we ask how strong hidden bias would need to be to alter our conclusions.

The sensitivity parameter Γ

Define $\Gamma$ as the odds ratio of treatment assignment for two units with identical observed covariates but potentially different unmeasured confounders:

where $\pi_i = P(T_i = 1 \mid X_i, U_i)$.

• $\Gamma = 1$: no hidden bias (perfect randomization within covariate strata).
• $\Gamma = 2$: unmeasured confounding could double the odds of treatment.
• $\Gamma = 3$: unmeasured confounding could triple the odds of treatment.

The critical Γ

We compute p-values at each $\Gamma$ level. The critical $\Gamma$ is the smallest value at which the treatment effect becomes statistically insignificant ($p > 0.05$). Larger critical $\Gamma$ indicates more robust findings.

Formal Rosenbaum bounds for the partially linear model

The classical Rosenbaum bounds apply to binary treatment. For DML’s partially linear model $Y = \theta T + g(X) + \varepsilon$, we adapt the framework following Rosenbaum ( Rosenbaum (2002) ):

The critical $\Gamma$ is then the solution to $p^{+}(\Gamma_{\text{crit}}) = \alpha$. Setting the worst-case p-value equal to $\alpha$ and solving:

This formula is implemented in dml_ts/sensitivity/rosenbaum.py:303-333.

E-value framework

VanderWeele and Ding ( VanderWeele & Ding (2017) ) introduced the E-value — a complementary sensitivity measure that asks: what is the minimum strength of association (on the risk ratio scale) that an unmeasured confounder would need with both treatment and outcome to explain away the observed effect?

For the OJ price elasticity, we convert our estimate to a risk ratio scale. A 1% price increase yields a $2.83\%$ demand decrease, giving $RR \approx e^{|{-2.83}| \cdot 0.01} \approx 1.029$ per percentage point — or approximately $RR \approx 4.0$ for a one-standard-deviation price change:

import numpy as np

# Convert elasticity to risk ratio for 1-SD price change
sd_price = np.std(data.T)
rr = np.exp(abs(result.theta) * sd_price)
e_value = rr + np.sqrt(rr * (rr - 1))

print(f"Risk ratio (1-SD change): {rr:.2f}")
print(f"E-value:                  {e_value:.2f}")
print(f"\nInterpretation: An unmeasured confounder must be")
print(f"associated with both price AND demand by a factor")
print(f"of {e_value:.1f}x to explain away the observed effect.")

Partial identification (Imbens 2003)

Rather than testing whether the effect survives hidden bias, Imbens ( Imbens (2003) ) asks: what range of treatment effects is consistent with the data under different levels of confounding?

sd_Y_resid = np.std(result.Y_residual)
sd_T_resid = np.std(result.T_residual)

print("Partial Identification Bounds")
print("=" * 50)
print(f"{'Gamma':>8} {'Lower':>10} {'Upper':>10} {'Width':>8}")
print("-" * 50)

for gamma in [1.0, 1.5, 2.0, 2.5, 3.0]:
    bias = (gamma - 1) / (gamma + 1) * sd_Y_resid / sd_T_resid
    lower = result.theta - bias
    upper = result.theta + bias
    print(f"{gamma:>8.1f} {lower:>10.3f} {upper:>10.3f} "
          f"{upper - lower:>8.3f}")

Partial identification bounds for OJ elasticity

Even under $\Gamma = 3.0$ — meaning an unmeasured confounder triples the odds of treatment — the identified set lies entirely below zero. The sign of the price elasticity is robust to extreme confounding. The magnitude uncertainty grows from a point estimate to a range of width $\sim 1.8$, but the qualitative conclusion (demand decreases with price) is unassailable.

Applying sensitivity analysis to OJ results

from dml_ts.sensitivity import compute_sensitivity_for_dml

sensitivity = compute_sensitivity_for_dml(
    theta=result.theta,
    se=result.se,
    n_samples=data.n_samples,
    treatment_r2=result.treatment_r2_cv,
    gamma_max=3.0,
    alpha=0.05,
)

print(sensitivity.summary())

Rosenbaum Bounds Sensitivity Analysis
=====================================
Treatment Effect:     θ̂ = -2.8347
Standard Error:       SE = 0.0412
Significance Level:   α = 0.05

Critical Gamma:       Γ_crit = 2.80
Interpretation:       Robust

Explanation:
  An unmeasured confounder would need to change treatment odds
  by a factor of 2.80x between similar units
  to render this effect statistically insignificant.

P-values at Selected Γ:
  Γ = 1.0: p = 0.0000 (no hidden bias)
  Γ = 1.5: p = 0.0000
  Γ = 2.0: p = 0.0001
  Γ = 3.0: p = 0.0847

Interpreting the sensitivity plot

from dml_ts.sensitivity import RosenbaumBounds
import matplotlib.pyplot as plt

bounds = RosenbaumBounds(gamma_max=3.0)
sensitivity_result = bounds.analyze(
    theta=result.theta,
    se=result.se,
    n_treated=data.n_samples // 2,
    n_control=data.n_samples // 2,
)

fig = bounds.plot_sensitivity(sensitivity_result)
plt.savefig("figures/ch04_sensitivity_plot.pdf")

The sensitivity plot shows p-values increasing as $\Gamma$ increases. The shaded region indicates where the effect remains statistically significant. The critical $\Gamma$ marks the boundary.

Interpretation guidelines

Sensitivity analysis interpretation guide

For our OJ estimate with $\Gamma_{\text{crit}} = 2.80$, we can report the price elasticity with confidence. An unmeasured confounder would need to have an implausibly strong effect on pricing to overturn our conclusions.

Contrastive application: gasoline price elasticity

The OJ results told a reassuring story: large effect, robust to confounding. But what does a fragile result look like? A single success provides no calibration — we need a contrasting failure to build intuition.

Why a second dataset?

A sensitivity analysis is only informative if we have seen both outcomes: results that survive scrutiny and results that do not. We construct a gasoline demand example where an omitted confounder (crude oil prices) creates genuine vulnerability.

Data construction

We construct a gasoline demand DGP with a hidden confounder:

import numpy as np

np.random.seed(42)
n = 5000

# Observed confounders: income, urban density
income = np.random.normal(50, 15, n)
urban = np.random.normal(0.6, 0.2, n)

# HIDDEN confounder: crude oil price (affects gas price AND demand)
crude_oil = np.random.normal(80, 20, n)

# Treatment: retail gasoline price
# Strongly driven by crude oil (the hidden confounder)
gas_price = (
    2.50
    + 0.015 * crude_oil       # strong crude oil pass-through
    + 0.005 * income           # modest income-area markup
    + np.random.normal(0, 0.15, n)
)

# Outcome: log gallons demanded
# True elasticity = -0.35 (inelastic demand)
log_demand = (
    3.0
    - 0.35 * np.log(gas_price)  # true causal effect
    + 0.008 * income
    - 0.3 * urban
    - 0.004 * crude_oil          # crude oil suppresses demand
    + np.random.normal(0, 0.3, n)
)

# Build observed dataset (crude_oil is UNOBSERVED)
X_gas = np.column_stack([income, urban])
T_gas = np.log(gas_price)
Y_gas = log_demand

print(f"Gasoline DGP: n={n}, true τ=-0.35")
print(f"Hidden confounder: crude oil (corr with T: "
      f"{np.corrcoef(crude_oil, T_gas)[0,1]:.3f})")

The design embeds a specific vulnerability: crude oil prices drive both the treatment (gas prices, via cost pass-through) and the outcome (demand, via macroeconomic channel). Because crude oil is unobserved, DML cannot adjust for it.

DML estimation

from dml_ts import double_ml

result_gas = double_ml(
    Y=Y_gas, T=T_gas, X=X_gas,
    n_folds=5, model="random_forest", random_state=42,
)

print(result_gas.summary())

Double Machine Learning Results
================================
Treatment Effect (θ):    -0.4821
Standard Error:          0.0523
95% Confidence Interval: [-0.5846, -0.3796]

Gasoline vs. OJ: estimation comparison

*True OJ elasticity unknown; $-2.8$ from literature benchmarks.

The gasoline estimate is biased by $-0.13$ (overstating price sensitivity by 37%) because the omitted crude oil confounder inflates the apparent price–demand relationship.

Sensitivity reveals fragility

from dml_ts.sensitivity import compute_sensitivity_for_dml

sensitivity_gas = compute_sensitivity_for_dml(
    theta=result_gas.theta,
    se=result_gas.se,
    n_samples=n,
    treatment_r2=0.45,   # high: crude oil drives price
    gamma_max=3.0,
    alpha=0.05,
)

print(f"Gasoline Γ_crit: {sensitivity_gas.gamma_crit:.2f}")
print(f"Verdict: {sensitivity_gas.interpretation}")

Gasoline Γ_crit: 1.28
Verdict: Sensitive

Lessons from the contrast

Heterogeneous sensitivity

The pooled $\Gamma_{\text{crit}} = 2.80$ for OJ is an average over all stores and brands. But sensitivity to unmeasured confounding may vary across subgroups — some segments of the market may be more vulnerable than others.

Motivation

Consider: a premium brand in high-income areas may have pricing driven by different factors than a store brand in low-income areas. If a specific subgroup shows low $\Gamma_{\text{crit}}$, it signals where additional data collection would be most valuable. This subgroup decomposition foreshadows the heterogeneous treatment effect analysis in Chapter 9.

Brand-level sensitivity

We estimate DML and sensitivity separately for each brand:

from dml_ts.data import OJDataLoader
from dml_ts import double_ml
from dml_ts.sensitivity import compute_sensitivity_for_dml

brands = ["dominicks", "minute.maid", "tropicana"]
brand_results = {}

for brand in brands:
    loader = OJDataLoader(brand=brand)
    bdata = loader.load()

    bresult = double_ml(
        Y=bdata.Y, T=bdata.T, X=bdata.X,
        n_folds=5, model="random_forest",
        random_state=42,
    )

    bsens = compute_sensitivity_for_dml(
        theta=bresult.theta,
        se=bresult.se,
        n_samples=bdata.n_samples,
        treatment_r2=bresult.treatment_r2_cv,
        gamma_max=3.0,
        alpha=0.05,
    )

    brand_results[brand] = {
        "n": bdata.n_samples,
        "theta": bresult.theta,
        "se": bresult.se,
        "gamma_crit": bsens.gamma_crit,
        "verdict": bsens.interpretation,
    }

    print(f"{brand:>15s}: n={bdata.n_samples:>6d}  "
          f"θ={bresult.theta:>7.3f}  "
          f"Γ_crit={bsens.gamma_crit:.2f}  "
          f"({bsens.interpretation})")

Brand-level sensitivity decomposition

Two patterns emerge: (1) the store brand (Dominick’s) shows the strongest elasticity and highest robustness, likely because its demand is most price-driven; (2) Tropicana, the premium brand, is least robust — brand loyalty provides a confounder-like buffer.

Store demographics

We split by median household income to test whether the sensitivity varies across market segments:

# Split by median INCOME (second confounder column)
median_income = np.median(data.X[:, 1])
low_inc = data.X[:, 1] <= median_income
high_inc = ~low_inc

for label, mask in [("Low income", low_inc),
                    ("High income", high_inc)]:
    r = double_ml(
        Y=data.Y[mask], T=data.T[mask], X=data.X[mask],
        n_folds=5, model="random_forest", random_state=42,
    )
    s = compute_sensitivity_for_dml(
        theta=r.theta, se=r.se,
        n_samples=mask.sum(),
        treatment_r2=r.treatment_r2_cv,
        gamma_max=3.0, alpha=0.05,
    )
    print(f"{label:>12s}: n={mask.sum():>6d}  "
          f"θ={r.theta:>7.3f}  Γ_crit={s.gamma_crit:.2f}")

Summary

Comprehensive sensitivity summary

Diagnostic visualization

Visual diagnostics complement the numerical results above. We present three diagnostic families: overlap, sensitivity surfaces, and residual quality.

Overlap diagnostic plots

import matplotlib.pyplot as plt
import numpy as np

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# (a) Treatment residual distribution
axes[0].hist(result.T_residual, bins=80, density=True,
             alpha=0.7, color="steelblue", edgecolor="white")
axes[0].set_xlabel("Treatment Residual (T̃)")
axes[0].set_ylabel("Density")
axes[0].set_title("(a) Treatment Residual Distribution")
axes[0].axvline(0, color="red", linestyle="--", alpha=0.5)

# (b) Leverage: |T_resid| vs influence score
influence = np.abs(result.influence_scores)
axes[1].scatter(np.abs(result.T_residual), influence,
                alpha=0.1, s=5, color="steelblue")
axes[1].set_xlabel("|Treatment Residual|")
axes[1].set_ylabel("|Influence Score|")
axes[1].set_title("(b) Leverage vs. Influence")

plt.tight_layout()
plt.savefig("figures/ch04_overlap_diagnostics.pdf",
            bbox_inches="tight")

The treatment residual histogram (a) should be roughly symmetric and centered at zero, confirming the first-stage model removes confounding without introducing systematic bias. The leverage–influence scatter (b) identifies high-leverage observations that disproportionately drive the estimate.

Sensitivity surface

The sensitivity surface shows how the p-value changes across $(\Gamma, |\theta|)$ space, revealing the joint robustness of effect size and hidden bias:

fig, ax = plt.subplots(figsize=(8, 5))

gamma_grid = np.linspace(1.0, 3.5, 100)
effect_grid = np.linspace(0.5, 4.0, 100)
G, E = np.meshgrid(gamma_grid, effect_grid)

# Compute p-value surface using Rosenbaum formula
n = data.n_samples
var_T = np.var(result.T_residual)
bias_shift = (G - 1) / (G + 1) * np.sqrt(n * var_T)
z_stat = E / result.se
p_surface = 1 - stats.norm.cdf(z_stat - bias_shift)

# Contour plot
cs = ax.contourf(G, E, p_surface,
                 levels=[0, 0.01, 0.05, 0.10, 0.50, 1.0],
                 cmap="RdYlGn_r", alpha=0.8)
ax.contour(G, E, p_surface, levels=[0.05],
           colors="black", linewidths=2)
plt.colorbar(cs, label="p-value")

# Mark OJ and gasoline estimates
ax.plot(2.80, abs(result.theta), "k*", markersize=15,
        label=f"OJ (Γ_crit={2.80})")
ax.plot(1.28, abs(result_gas.theta), "rs", markersize=10,
        label=f"Gasoline (Γ_crit={1.28})")

ax.set_xlabel("Sensitivity Parameter Γ")
ax.set_ylabel("|Treatment Effect|")
ax.set_title("Sensitivity Surface: p-value vs. (Γ, |θ|)")
ax.legend(loc="upper right")

plt.savefig("figures/ch04_sensitivity_surface.pdf",
            bbox_inches="tight")

Residual Q–Q and influence diagnostics

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# (a) Q-Q plot of outcome residuals
stats.probplot(result.Y_residual, dist="norm",
               plot=axes[0])
axes[0].set_title("(a) Outcome Residual Q-Q Plot")
axes[0].get_lines()[0].set_color("steelblue")

# (b) Influence score distribution
axes[1].hist(result.influence_scores, bins=80,
             density=True, alpha=0.7,
             color="coral", edgecolor="white")
axes[1].set_xlabel("Influence Score")
axes[1].set_ylabel("Density")
axes[1].set_title("(b) Influence Score Distribution")

# Flag extreme influence (|score| > 3*IQR)
iqr = np.percentile(result.influence_scores, 75) - \
      np.percentile(result.influence_scores, 25)
extreme = np.abs(result.influence_scores) > 3 * iqr
axes[1].axvline(3 * iqr, color="red", linestyle="--",
                alpha=0.5, label=f"3×IQR ({extreme.sum()} obs)")
axes[1].axvline(-3 * iqr, color="red", linestyle="--",
                alpha=0.5)
axes[1].legend()

plt.tight_layout()
plt.savefig("figures/ch04_residual_diagnostics.pdf",
            bbox_inches="tight")

The Q–Q plot reveals departures from normality in the residual tails. Moderate departures are acceptable (the CLT handles them), but extreme heavy tails warrant larger samples or robust standard errors.

Practical recommendations

When to trust DML estimates

Based on our experience with the OJ application, we recommend:

1. First-stage $R^2$ in moderate range: outcome $R^2 > 0.1$ (confounders matter); treatment $R^2 < 0.8$ (overlap maintained).
2. Estimates in expected range: compare to literature benchmarks when available.
3. Sensitivity analysis passes: $\Gamma_{\text{crit}} > 1.5$ for policy recommendations.
4. Results stable across specifications: try different nuisance models, fold counts.

Reporting checklist

When presenting DML results in an applied paper, every item below should appear. The checklist below maps each item to its location and rationale.

DML reporting checklist

Exercises

1. Brand-specific elasticity: use OJDataLoader(brand="tropicana") to estimate elasticity for a single brand. How does it compare to the pooled estimate? What does heterogeneity across brands suggest about market segmentation?
2. Extended confounders: add EDUC, ETHNIC, and HHLARGE to the confounder set. Does the elasticity estimate change? What does this suggest about the adequacy of the original three-variable specification?
3. Sensitivity exploration: for your brand-specific estimate from Exercise 1, conduct sensitivity analysis. Is the single-brand estimate more or less robust than the pooled estimate? Why might this be?
4. Nuisance model comparison: re-run the OJ analysis using model="ridge" and model="gradient_boosting". How do the estimates and diagnostics change? When might you prefer each model?
5. E-value computation: compute E-values for (a) the pooled OJ estimate, (b) your brand-specific estimate from Exercise 1, and (c) the gasoline estimate from the contrastive section. Compare the E-value ranking with the $\Gamma_{\text{crit}}$ ranking. Do they agree? When might they diverge?
6. Gasoline with FRED controls: modify the gasoline DGP to include crude oil price as an observed confounder. How does $\hat{\theta}$ change? How does $\Gamma_{\text{crit}}$ change? This previews Chapter 7’s use of macroeconomic controls.
7. Overlap failure DGP: construct a DGP where treatment is nearly deterministic given covariates (treatment $R^2 > 0.95$). Run DML on this data. What happens to the standard errors? How does trimming help?
8. Sensitivity surface interpretation: using the sensitivity surface code from the visualization section, generate the surface for your brand-specific estimate. Identify the region where the effect would become insignificant. How does the shape of this boundary compare to the pooled estimate’s surface?

Chapter summary

This chapter demonstrated the complete DML pipeline on real observational data, from estimation through multi-layered robustness assessment.

Looking ahead

Chapter 5 extends these methods to time series data, where temporal dependence introduces new challenges for causal inference. The sensitivity framework developed here carries forward: time series DML estimates require the same robustness assessment, but with additional complications from autocorrelation and non-stationarity. The gasoline example’s confounder gap foreshadows Chapter 7, where FRED macroeconomic controls close exactly this type of omitted variable problem.