Lin et al. (2023) E-value unmeasured confounder sensitivity

lin_e_value.py

"""Lin, Ni & Lu (2023) §4.2 E-value sensitivity for a matched RMST difference.

Propensity-score matching balances *measured* covariates. An unmeasured
confounder ``U`` can still bias the matched ATT. The E-value answers a
single question:

    "How strong would U's association with treatment *and* U's association
     with the outcome each have to be, on a risk-ratio / mean-ratio scale,
     for that confounder to fully explain away the observed effect?"

A large E-value means only an implausibly strong unmeasured confounder
could account for the result; a small one means a modest confounder
suffices. It is a single number, comparable across studies, that sits
alongside the point estimate and its CI.

Ding & VanderWeele (2016) introduced E-values for risk ratios. Lin, Ni &
Lu (2023, §4.2, https://doi.org/10.1515/jci-2022-0035) adapt the
bounding-factor logic to a marginal RMST *difference* under a matched
design. All the identifying algebra runs on two scalars: the arm-marginal
RMSTs at the truncation horizon ``tau``. That is what this file uses.

Formulas we implement
---------------------
Let ``r1 = E[min(T, tau) | A = 1]`` and ``r0 = E[min(T, tau) | A = 0]`` be
the arm-marginal RMSTs on the matched cohort. The observed ATT is
``delta = r1 - r0``. Set the two sensitivity parameters equal,
``RR_AU = MR_UZ = e`` for some ``e >= 1`` — the standard "equal both
sides" convention — and Lin et al.'s eq (1) reduces to the bounding
factor

    BF(e) = e**2 / (2*e - 1).

Eq (2): if the observed effect is positive, the bias-adjusted *lower*
bound on the true ATT is

    L(e) = 0.5 * (1 + 1/BF) * r1  -  0.5 * (1 + BF) * r0.

Eq (3) is the symmetric *upper* bound for a negative observed effect:

    U(e) = 0.5 * (1 + BF) * r1  -  0.5 * (1 + 1/BF) * r0.

Both reduce to ``delta`` at ``e = 1``. As ``e`` grows, bias eats into
the effect (L decreases, U increases). The E-value is the smallest ``e``
at which the adjusted bound just reaches the null.

1:k matching
------------
Lin et al. §2 use **1:1** pair matching and build the CI from Murray &
Cole's closed-form paired-KM covariance. We use **1:k matching with
replacement** and bootstraps *matched sets* to build the CI — matched
sets are the resampling unit, and each treated unit carries its k
controls with it. That correctly propagates the induced 1:k dependence
into the SE on ``delta``.

The §4.2 identification result depends only on the arm-marginal RMSTs
``r1`` and ``r0``; it is silent on *how* those marginals or the CI were
computed. The bounding-factor / adjusted-bound algebra therefore
transfers to the 1:k setting unchanged. The genuinely new piece is that
the CI we plug into the CI-based E-value is a bootstrap-over-matched-sets
limit rather than Lin's pair-based limit; for the CI adjustment we take
the standard VanderWeele–Ding stance and hold the bootstrap SE fixed
while shifting the CI by the same amount that bias adjusts the point.
"""

from __future__ import annotations


def _bounding_factor(e: float) -> float:
    """Lin et al. eq (1) with the ``RR_AU = MR_UZ = e`` convention."""
    return e * e / (2.0 * e - 1.0)


def _adjusted_att_bound(rmst_1: float, rmst_0: float, e: float, sign: int) -> float:
    """Lin et al. eq (2) for ``sign=+1`` (lower bound) or eq (3) for ``sign=-1``."""
    bf = _bounding_factor(e)
    if sign > 0:
        return 0.5 * (1.0 + 1.0 / bf) * rmst_1 - 0.5 * (1.0 + bf) * rmst_0
    return 0.5 * (1.0 + bf) * rmst_1 - 0.5 * (1.0 + 1.0 / bf) * rmst_0


def _solve_e(rmst_1: float, rmst_0: float, target: float, sign: int) -> float:
    """Smallest ``e >= 1`` driving the adjusted bound to ``target``; 1.0 if already reached."""
    def f(e: float) -> float:
        return _adjusted_att_bound(rmst_1, rmst_0, e, sign) - target

    if sign * f(1.0) <= 0.0:
        return 1.0
    lo, hi = 1.0, 2.0
    while sign * f(hi) > 0.0:
        hi *= 2.0
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        if sign * f(mid) > 0.0:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


def lin_evalue(
    rmst_1: float, rmst_0: float, ci_bound: float | None = None
) -> dict[str, float]:
    """Lin et al. (2023, §4.2) E-value for a matched RMST difference.

    Parameters
    ----------
    rmst_1, rmst_0
        Arm-marginal RMSTs at ``tau`` on the matched cohort (treated, control).
    ci_bound
        Optional bootstrap CI limit closest to the null on
        ``delta = rmst_1 - rmst_0`` — the one-sided *lower* limit if
        ``delta > 0``, the *upper* limit if ``delta < 0``. If given, an
        E-value for that CI limit is returned alongside the point E-value.

    Returns
    -------
    dict
        ``{"e_value_point": ...}`` and, if ``ci_bound`` is provided,
        ``{"e_value_ci": ...}``. A returned value of exactly 1.0 means no
        confounding is needed to reach the null — the observed CI already
        crosses zero at ``e = 1``.
    """
    delta = rmst_1 - rmst_0
    if delta == 0.0:
        out = {"e_value_point": 1.0}
        if ci_bound is not None:
            out["e_value_ci"] = 1.0
        return out
    sign = 1 if delta > 0 else -1
    out = {"e_value_point": _solve_e(rmst_1, rmst_0, 0.0, sign)}
    if ci_bound is not None:
        # Shift approximation: bias-adjusted CI limit = ci_bound + (bound(e) - delta).
        # Setting the adjusted CI limit to zero gives target = delta - ci_bound.
        out["e_value_ci"] = _solve_e(rmst_1, rmst_0, delta - ci_bound, sign)
    return out