Refutation of DML IRM inference

This notebook explain refutation tests that implemented in Causalis DGP took from benchmarking notebook

Inference

As we see our estimate is accurate and CI bounds include ground-truth ATE. In real life we can't compare estimation with truth so we need check robustness of it: run some tests on assumptions and answer questions about research design

Overlap

What “overlap/positivity” means

Binary treatment T \in \{0,1\}: for all confounder values $x$ in your target population,

often strengthened to strong positivity: there exists an \varepsilon > 0 such that

Why it matters

• Identification: Overlap + unconfoundedness are the two pillars that identify causal effects from observational data. Without overlap, the effect is not identified — you must extrapolate or model-specify what never occurs.

• Estimation stability: IPW/DR estimators use weights

  $$w_1 = \frac{D}{e(X)}, \qquad w_0 = \frac{1 - D}{1 - e(X)}.$$

  If $e(X)$ is near 0 or 1, these weights explode, causing huge variance and fragile estimates.

• Target population: With trimming or restriction, you may change who the effect describes — e.g., ATE on the region of common support, not on the full population.

It's summary for overlap diagnostics

edge_mass

edge_0.01_below, edge_0.01_above, edge_0.02_below, edge_0.02_above are shares of units whose propensity is below or above the percents

To keep in mind: DML IRM is clipping out the interval [0.02 and 0.98]

Huge shares are dangerous for estimation in terms of weights exploding

Flags in Causalis:

For $ε=0.01$: YELLOW if either side 0.02 (2%), RED if 0.05 (5%).

For $ε=0.02$: YELLOW if either side 0.05 (5%), RED if 0.10 (10%).

ks - Kolmogorov–Smirnov statistic

Here KS is the two-sample Kolmogorov–Smirnov statistic comparing the distributions of the propensities for treated vs control:

$D=\max_t |\hat F_A(t)-\hat F_B(t)|$

Interpretation:

• (D=0): identical distributions (perfect overlap).
• (D=1): complete separation (no overlap).
• Your value KS = 0.5116 means there exists a threshold (t) such that the share of treated with $(m\le t)$ differs from the share of controls with $(m\le t)$ by ~51 percentage points. That’s why it’s flagged RED (your thresholds mark RED when (D>0.35)): treatment assignment is highly predictable from covariates ⇒ poor overlap / strong confounding risk.

AUC

Probability definition (most intuitive)

where $(s^+)$ is a score from a random positive and $(s^-)$ from a random negative. So AUC is the fraction of all $(n_1 n_0)$ positive–negative pairs that are correctly ordered by the score (ties get half-credit).

Rank / Mann–Whitney formulation

1. Rank all scores together (ascending). If there are ties, assign average ranks within each tied block.

2. Let $(R_1)$ be the sum of ranks for the positives.

3. Compute the Mann–Whitney (U) statistic for positives:

   $$U_1 = R_1 - \frac{n_1(n_1+1)}{2}.$$

4. Convert to AUC by normalizing:

   $$\boxed{\text{AUC} = \frac{U_1}{n_1 n_0} = \frac{R_1 - \frac{n_1(n_1+1)}{2}}{n_1 n_0}}$$

   This is exactly what your function returns (with stable sorting and tie-averaged ranks).

ROC-integral view (equivalent)

If (\text{TPR}(t)) and (\text{FPR}(t)) trace the ROC curve as the threshold $(t)$ moves,

i.e., the geometric area under the ROC.

Properties you should remember

• Range: (0 \le \text{AUC} \le 1); 0.5 = random ranking; 1 = perfect separation.
• Symmetry: (\text{AUC}(s,y) = 1 - \text{AUC}(s,1-y)).
• Monotone invariance: Any strictly increasing transform $(f)$ leaves AUC unchanged (only ranks matter).
• Ties: Averaged ranks ⇒ adds the $(\tfrac12\Pr(s^+=s^-))$ term automatically.

In the propensity/overlap context

• A higher AUC means treatment (D) is more predictable from covariates (bad for overlap/positivity).
• For good overlap you actually want AUC close to 0.5.

ESS_treated_ratio

Weights used

For ATE-style IPW, the treated-arm weights are

so on the treated subset ({i:D_i=1}) the weights are simply $(1/m_i)$.

Effective sample size (ESS)

Given the treated-arm weights (w_1,\ldots,w_{n_1}) (only for $(D=1)$),

This is exactly what _ess(w) computes.

• If all treated weights are equal, ESS $(= n_1)$ (full efficiency).
• If a few weights dominate, ESS drops (information concentrated in few units).

The reported metric

This lies in $(0,1]$. Near 1 ⇒ well-behaved weights; near 0 ⇒ severe instability.

Why it reflects overlap

When propensities $(m_i)$ approach 0 for treated units, weights $(1/m_i)$ explode → large CV → low ESS_treated_ratio. Hence this metric is a direct, quantitative read on how much usable information remains in the treated group after IPW.

tails_w1_q99/med

Interpretation

• It’s a tail-heaviness index for treated weights: how large the 99th-percentile weight is relative to a typical (median) weight.
• Scale-invariant: if you re-scale weights (e.g., Hájek normalization), both numerator and denominator scale equally, so the ratio is unchanged.
• Bigger $(\Rightarrow)$ heavier right tail $(\Rightarrow)$ more variance inflation for IPW (since variance depends on large $(w_i^2)$). It typically coincides with a low ESS(_\text{treated ratio}).

Edge cases & thresholds

• If (\text{median}(W_1)=0) or undefined, the ratio is not meaningful (your code returns “NA” in that case; with positive treated weights this is rare).
• Defaults: YELLOW if any of ({q95/med,q99/med,q999/med,\max/med}) exceeds 10; RED if any exceed 100. tails_w1_q99/med is one of these checks, focusing specifically on the 99th percentile.

Quick example

If \mathrm{median}(W_1) = 1.2 and Q_{0.99}(W_1) = 46.8, then

indicating heavy tails and a likely unstable ATE IPW.

ATT_identity_relerr

With estimated propensities $(m_i=\hat m(X_i))$ and (D_i\in{0,1}):

• Left-hand side (controls odds sum): $$\text{LHS} = \sum_{i=1}^n (1-D_i),\frac{m_i}{1-m_i}.$$
• Right-hand side (treated count): $$\text{RHS} = \sum_{i=1}^n D_i = n_1.$$

If $(\hat m\approx m)$ and overlap is ok, LHS $(\approx)$ RHS.

You report the relative error:

(when (n_1>0) otherwise it’s set to $(\infty)$).

How to read the number

• Small (\texttt{relerr}) (e.g., (\le 5%)) ⇒ propensities are reasonably calibrated (especially on the control side) and ATT weights won’t be wildly off in total mass.
• Large (\texttt{relerr}) ⇒ possible miscalibration of $(\hat m)$ (e.g., over/underestimation for controls), poor overlap (many controls with $(m_i\to 1)$ inflating $(m_i/(1-m_i)))$, or clipping/trimming effects.

Your default flags (same as in the code):

• GREEN if (\texttt{relerr} \le 0.05)
• YELLOW if (0.05 < \texttt{relerr} \le 0.10)
• RED if (> 0.10)

Quick intuition

The term $(m/(1-m))$ is the odds of treatment. Summing that over controls should reconstruct the treated count. If it doesn’t, either the odds are off (propensity miscalibration) or the data lack support where you need it—both are red flags for ATT-IPW stability.

clip_m_total

look at edge_mass

calib_ECE, calib_slope, calib_intercept

calib_ECE = 0.018 (GREEN)

Math: with 10 equal-width bins,

(weighted average gap between observed rate (\bar y_k) and mean prediction (\bar p_k) per bin). Result: ~1.8% average miscalibration → overall probabilities track outcomes well. Note the biggest bin error is in 0.5–0.6 (abs_error ≈ 0.162) but it’s tiny (95/10,000), so ECE stays low.

calib_slope (β) = 0.889 (GREEN)

Math (logistic recalibration):

Interpretation: (\beta<1) ⇒ predictions are a bit over-confident (too extreme); the optimal calibration slightly flattens them toward 0.5.

calib_intercept (α) = −0.107 (GREEN)

Math: same model as above; (\alpha) is a vertical shift on the log-odds scale. Interpretation: Negative (\alpha) nudges probabilities downward overall (your model is, on average, a bit high), consistent with bins like 0.5–0.6 where (\bar p_k &gt; \bar y_k).

All three fall well within your GREEN thresholds, so calibration looks solid despite minor mid-range overprediction.

Score

We need this score refutation tests for:

• Catch overfitting/leakage: The out-of-sample moment check verifies that the AIPW score averages to ~0 on held-out folds using fold-specific θ and nuisances. If this fails, your effect can be an artifact of leakage or overfit learners rather than a real signal.
• Verify Neyman orthogonality in practice: The Gateaux-derivative tests (orthogonality_derivatives) check that small, targeted perturbations to the nuisances (g₀, g₁, m) don’t move the score mean. Large |t| values flag miscalibration (e.g., biased propensity or outcome models) that breaks the orthogonality protection DML relies on.
• Assess finite-sample stability: The influence diagnostics reveal heavy tails (p99/median, kurtosis) and top-influential points. Spiky ψ implies high variance and sensitivity—often due to near-0/1 propensities, poor overlap, or outliers.
• ATTE-specific risks: For ATT/ATTE, only g₀ and m matter in the score. The added overlap metrics and trim curves show how reliant your estimate is on scarce, high-m controls—common failure mode for ATT.

psi_p99_over_med

• Let $\psi_i$ be the per-unit influence value (EIF score) for your estimator. We look at magnitudes $a_i \equiv |\psi_i|$.

• Define the 99th percentile and the median of these magnitudes:

  $$q_{0.99} \equiv \operatorname{Quantile}_{0.99}(a_1,\dots,a_n),\qquad m \equiv \operatorname{median}(a_1,\dots,a_n).$$

• The metric is the scale-free tail ratio:

  $$\boxed{ \texttt{psi\_p99\_over\_med} = \frac{q_{0.99}}{m} }$$

Why this works (brief):

• Uses $|\psi_i|$ to ignore sign (only tail size matters).
• Dividing by the median makes it scale-invariant and robust to a few large values.
• Large values $\big(\gg 1\big)$ mean a small fraction of observations dominate uncertainty (heavy tails → unstable SE).

Quick read:

• $\approx 1!-!5$: tails tame/stable
• $\gtrsim 10$: caution (heavy tails)
• $\gtrsim 20$: likely unstable; check overlap, trim/clamp propensities, or robustify learners.

psi_kurtosis

• Let $\psi_i$ be the per-unit influence values and define centered residuals

  $$\tilde\psi_i \equiv \psi_i - \bar\psi,\qquad \bar\psi \equiv \frac{1}{n}\sum_{i=1}^n \psi_i.$$

• Sample variance (with Bessel correction):

  $$s^2 \equiv \frac{1}{n-1}\sum_{i=1}^n \tilde\psi_i^2.$$

• Sample 4th central moment:

  $$\hat\mu_4 \equiv \frac{1}{n}\sum_{i=1}^n \tilde\psi_i^4.$$

• The reported metric (raw kurtosis, not excess):

  $$\boxed{ \texttt{psi\_kurtosis} = \frac{\hat{\mu}_4}{s^4} }$$

Interpretation (quick):

• Normal reference $\approx 3$ (excess kurtosis $=0$).
• Much larger $\Rightarrow$ heavier tails / more extreme $\psi_i$ outliers.
• Rules of thumb used in the diagnostics: $\ge 10$ = caution, $\ge 30$ = severe.

$max_|t|_g1$, $max_|t|_g0$, $max_|t|_m$

We work with a basis of functions

Let (m_i^\tau \equiv \mathrm&#123;clip&#125;(m_i,\tau,1-\tau)) be the clipped propensity (guards against division by zero).

ATE case

For each basis function $(b)$, form a sample mean (Gateaux derivative estimator) and its standard error, then compute a t-statistic; finally take the maximum absolute value across bases.

$(g_1)$ direction

$(g_0)$ direction

$(m)$ direction

Interpretation: under Neyman orthogonality, each derivative mean (\widehat d_&#123;\bullet,b&#125;) should be approximately zero, so all (|t_&#123;\bullet,b&#125;|) should be small. Large (\max_&#123;|t|&#125;) values flag miscalibration of the corresponding nuisance.

ATTE / ATT case

Let (p_1 = \mathbb&#123;E&#125;[D]) and define the odds $(o_i = m_i^\tau / (1 - m_i^\tau))$.

• The $(g_1)$ derivative is identically zero:

  $$\Rightarrow\quad \max_{|t|_{g_1}} = 0.$$

• $(g_0)$ direction

  $$\widehat d_{g_0,b} = \frac{1}{n}\sum_i h_b(X_i)\frac{(1-D_i)o_i - D_i}{p_1}, \qquad t_{g_0,b} = \frac{\widehat d_{g_0,b}}{\mathrm{se}(\widehat d_{g_0,b})}, \qquad \max_{|t|_{g_0}} = \max_b |t_{g_0,b}|.$$

• $(m)$ direction

  $$\widehat d_{m,b} = -\frac{1}{n}\sum_i h_b(X_i) \frac{(1-D_i)(Y_i - g_{0,i})} {p_1(1 - m_i^\tau)^2}, \qquad \max_{|t|_{m}} = \max_b |t_{m,b}|.$$

Rule of thumb: (\max_&#123;|t|&#125; \lesssim 2) is “okay”; larger values indicate orthogonality breakdown — fix by recalibrating that nuisance, changing learners, features, or trimming.

oos_tstat_fold, oos_tstat_strict

Here’s the math behind the two OOS (out-of-sample) moment t-stats used in the diagnostics. Assume K-fold cross-fitting with held-out index sets $(I_k)$ (size $n_k$) and complements $(R_k)$.

Step 1 — Leave-fold-out (\hat\theta_&#123;-k&#125;)

For the moment condition (\mathbb&#123;E&#125;[\psi_a(W)\,\theta + \psi_b(W)] = 0), the leave-fold-out estimate used on fold $(k)$ is

Step 2 — Held-out scores on fold $(k)$

Define the fold-specific held-out score for $i\in I_k$:

Compute per-fold mean and variance:

OOS t-stat diagnostics

(\texttt&#123;oos\_tstat\_fold&#125;)

A fold-aggregated, variance-weighted t-statistic:

Intuition: averages fold means and scales by a fold-pooled standard error.

(\texttt&#123;oos\_tstat\_strict&#125;)

A “strict” t-stat using every held-out observation directly:

Intuition: computes a single overall mean and standard error across all held-out scores (often slightly more conservative).

Interpretation

Under a valid design and correct cross-fitting (so that \mathbb&#123;E&#125;[\psi]=0 out-of-sample), both statistics are approximately standard normal:

Values near $0$ indicate that the moment condition holds out of sample. Large $|t|$ suggests overfitting, leakage, or nuisance miscalibration.

SUTVA

Those assumptions are statistically untestable. We need design of research for them

Uncofoundedness

balance\_max\_smd

For each covariate $(X_j)$, the (weighted) standardized mean difference is

Group means and variances are computed under the IPW weights implied by your estimand:

• ATE: w_&#123;1i&#125; = \tfrac&#123;D_i&#125;&#123;\hat m_i&#125;, w_&#123;0i&#125; = \tfrac&#123;1-D_i&#125;&#123;1-\hat m_i&#125;
• ATTE: w_&#123;1i&#125; = D_i, w_&#123;0i&#125; = (1-D_i)\tfrac&#123;\hat m_i&#125;&#123;1-\hat m_i&#125;

(If normalize=True, each weight vector is divided by its mean.)

Weighted means and variances:

Special cases in the code:

• If both variances are $\approx 0$ and |\mu_&#123;1j&#125;-\mu_&#123;0j&#125;| \approx 0 ⇒ \mathrm&#123;SMD&#125;_j = 0
• If both variances are $\approx 0$ but means differ ⇒ \mathrm&#123;SMD&#125;_j = \infty
• If denominator is $\approx 0$ otherwise ⇒ \mathrm&#123;SMD&#125;_j = \text&#123;NaN&#125;

Then

implemented as a nanmax over the vector of \mathrm&#123;SMD&#125;_j. NaNs are ignored; if any feature produced $\infty$, the max is $\infty$.

balance\_frac\_violations

Let the SMD threshold be $\tau$ (default $0.10$). Define the set of finite SMDs:

Then the fraction of violations is

So it’s the share of covariates whose weighted SMD exceeds the threshold, computed only over finite SMDs (NaN / Inf are excluded from the denominator).

Quick interpretation

• Smaller is better. A common rule of thumb is \mathrm&#123;SMD&#125; \le 0.10.
• balance_max_smd tells you the worst residual imbalance across covariates.
• balance_frac_violations tells you how many covariates (as a fraction) still exceed the chosen threshold.

Sensitivity analysis

1) sensitivity_analysis: bias-aware CI

Goal. Start from your estimator $(\hat\theta)$ with sampling standard error $(se)$. Allow a controlled amount of worst-case hidden confounding through three knobs $(cf_y, cf_d, \rho)$. Inflate the uncertainty by an additive “max bias”.

Step A — Sampling part

• Point estimate $\hat\theta$, standard error $se$, and $z_\alpha$ for level $(1-\alpha)$.
• Usual sampling CI:
• $$[\,\hat\theta - z_\alpha\,se,\ \hat\theta + z_\alpha\,se\,].$$

Step B — Confounding geometry

The code pulls sensitivity elements from the fitted IRM:

• $\sigma^2$: the asymptotic variance of the estimator’s EIF (so that se = \sqrt&#123;\sigma^2&#125; in the module’s normalization).

• $m_\alpha(i) \ge 0$: per-unit weight for the outcome channel (how outcome-model misspecification moves the EIF).

• $r(i)$ (“riesz_rep”): per-unit weight for the treatment channel (how propensity-model misspecification moves the EIF).

We turn the user’s sensitivity knobs into a quadratic budget for adversarial confounding:

• $cf_y \ge 0$: strength of unobserved outcome disturbance
• $cf_d \ge 0$: strength of unobserved treatment disturbance
• $\rho \in [-1,1]$: their correlation

This $\nu^2$ is a dimensionless bias multiplier — how sensitive the EIF is to those perturbations.

Step C — Max bias and intervals

Two equivalent forms appear in the code:

Then the module reports:

• Confounding bounds for $\theta$:

  $$[\,\hat\theta - \text{max\_bias},\; \hat\theta + \text{max\_bias}\,].$$

• Bias-aware CI (sampling + confounding, worst-case additive):

  $$\Big[\,\hat\theta - (\text{max\_bias} + z_\alpha\,se),\; \hat\theta + (\text{max\_bias} + z_\alpha\,se)\,\Big].$$

(So you’re adding sampling error and the adversarial bias linearly for a conservative envelope.)

Notes & edge handling

• Numeric PSD clamping ensures \text&#123;base&#125;_i \ge 0; $\rho$ is clipped to $[-1,1]$.

• If $cf_y = cf_d = 0 \Rightarrow \nu^2 = 0 \Rightarrow$ bias-aware CI collapses to the sampling CI.

• Internally, a delta-method IF for \text&#123;max\_bias&#125; is

  $$\psi_{\text{max}}(i) = \frac{\sigma^2\,\psi_{\nu^2}(i) + \nu^2\,\psi_{\sigma^2}(i)} {2\,\text{max\_bias}},$$

  matching \text&#123;max\_bias&#125; = \sqrt&#123;\sigma^2\nu^2&#125; (used for coherent summaries).

2) sensitivity_benchmark: calibrating $(cf_y, cf_d, \rho)$ from omitted covariates

Goal. Pick a set $Z$ of candidate “omitted” covariates (the benchmarking_set). Refit a short IRM that excludes $Z$ and compare it to the long (original) model. Use how well $Z$ explains residual variation to derive plausible $(cf_y, cf_d, \rho)$.

Step A — Long vs short estimates

• Long: \hat\theta_&#123;\text&#123;long&#125;&#125; (original model).
• Short: \hat\theta_&#123;\text&#123;short&#125;&#125; (drop $Z$, same learners/hyperparams).
• Report \Delta = \hat\theta_&#123;\text&#123;long&#125;&#125; - \hat\theta_&#123;\text&#123;short&#125;&#125;.

Step B — Residuals from the long model

Let $(g_1, g_0, \hat m)$ be the outcome and propensity learners:

These are the EIF’s outcome and treatment residual components.

Step C — How much of each residual does $Z$ explain?

Regress $r_y$ on $Z$ and $r_d$ on $Z$ (unweighted OLS; ATT case uses ATT weights):

• Obtain $R^2_y$ and $R^2_d$.

• Convert to signal-to-noise ratios (the “strength” of confounding channels):

  $$cf_y = \frac{R^2_y}{1 - R^2_y}, \qquad cf_d = \frac{R^2_d}{1 - R^2_d}.$$

  (These are the same $R^2 / (1 - R^2)$ maps used in modern partial-$R^2$ robustness frameworks.)

Compute the correlation between the fitted pieces from those two regressions:

weighted for ATT when applicable, then clipped to $[-1,1]$.

Outputs

A one-row DataFrame (indexed by the treatment name) with

You can pass $(cf_y, cf_d, \rho)$ straight into sensitivity_analysis to get the associated bias-aware interval. Intuitively, this calibrates how strong hidden stuff would need to be by using a concrete, observed proxy $Z$.

How to read them together

1. Use sensitivity_benchmark with a plausible omitted set $Z$ to derive $(cf_y, cf_d, \rho)$ and observe the actual estimate shift $\Delta$.

2. Plug those $(cf_y, cf_d, \rho)$ into sensitivity_analysis to get:

   $$\text{max\_bias} = \sqrt{\nu^2}\,se, \qquad \text{Bias-aware CI} = \hat\theta \pm (\text{max\_bias} + z_\alpha\,se).$$

Small $cf$ values (or $\rho \approx 0$) ⇒ tiny $\nu^2$ ⇒ bias-aware CI near the sampling CI. Large $cf$ values and $|\rho|\approx 1$ widen it, reflecting stronger plausible hidden confounding.