Diagnostics

compute_positivity_metrics

compute_positivity_metrics(df: DataFrame, ps_col: str = PS_COL, treatment_col: str = TREATMENT_COL, eps: float = 0.01, common_support_threshold: float = 0.05) -> dict

Compute positivity/overlap diagnostics for the propensity score.

Parameters:

Name	Type	Description	Default
df	DataFrame	Input DataFrame with treatment and propensity score columns.	required
ps_col	str	Name of the propensity score column.	PS_COL
treatment_col	str	Name of the treatment status column (1 treated, 0 control).	TREATMENT_COL
eps	float	Propensity scores outside [eps, 1 - eps] are counted as extreme.	0.01
common_support_threshold	float	Quantile threshold passed to get_common_support_range.	0.05

Returns:

Type	Description
dict	dict with sample sizes, shares of extreme propensity scores (overall and
dict	per arm), the trimmed common support range (intersection of the arms'
dict	quantile-trimmed PS ranges — observations can fall outside it even under
dict	perfect overlap when common_support_threshold > 0) and shares outside it,
dict	the KS test comparing the arms' propensity score distributions, and a
dict	flag_extreme_ps bool (True if any propensity scores fall outside
dict	[eps, 1 - eps] — evidence of a practical positivity problem, not proof of
dict	a formal violation).

Source code in CausalEstimate/diagnostics/positivity.py

def compute_positivity_metrics(
    df: pd.DataFrame,
    ps_col: str = PS_COL,
    treatment_col: str = TREATMENT_COL,
    eps: float = 0.01,
    common_support_threshold: float = 0.05,
) -> dict:
    """
    Compute positivity/overlap diagnostics for the propensity score.

    Args:
        df: Input DataFrame with treatment and propensity score columns.
        ps_col: Name of the propensity score column.
        treatment_col: Name of the treatment status column (1 treated, 0 control).
        eps: Propensity scores outside [eps, 1 - eps] are counted as extreme.
        common_support_threshold: Quantile threshold passed to
            get_common_support_range.

    Returns:
        dict with sample sizes, shares of extreme propensity scores (overall and
        per arm), the trimmed common support range (intersection of the arms'
        quantile-trimmed PS ranges — observations can fall outside it even under
        perfect overlap when common_support_threshold > 0) and shares outside it,
        the KS test comparing the arms' propensity score distributions, and a
        flag_extreme_ps bool (True if any propensity scores fall outside
        [eps, 1 - eps] — evidence of a practical positivity problem, not proof of
        a formal violation).
    """
    if not 0 < eps < 0.5:
        raise ValueError(f"eps must be in (0, 0.5), got {eps}.")
    n_treated, n_control = validate_ps_and_treatment(df, ps_col, treatment_col)

    treated_ps = get_treated_ps(df, treatment_col, ps_col)
    control_ps = get_untreated_ps(df, treatment_col, ps_col)
    ps = df[ps_col]
    support_low, support_high = get_common_support_range(
        df, treatment_col, ps_col, common_support_threshold
    )
    ks_stats = compute_propensity_score_stats(df, ps_col, treatment_col)

    def prop_extreme(s: pd.Series) -> float:
        return float(((s < eps) | (s > 1 - eps)).mean())

    def prop_outside(s: pd.Series) -> float:
        return float(((s < support_low) | (s > support_high)).mean())

    prop_ps_extreme = prop_extreme(ps)
    return {
        "n_total": int(len(df)),
        "n_treated": n_treated,
        "n_control": n_control,
        "prop_ps_below_eps": float((ps < eps).mean()),
        "prop_ps_above_1_minus_eps": float((ps > 1 - eps).mean()),
        "prop_ps_extreme": prop_ps_extreme,
        "prop_treated_ps_extreme": prop_extreme(treated_ps),
        "prop_control_ps_extreme": prop_extreme(control_ps),
        "common_support_low": float(support_low),
        "common_support_high": float(support_high),
        "prop_outside_support": prop_outside(ps),
        "prop_treated_outside_support": prop_outside(treated_ps),
        "prop_control_outside_support": prop_outside(control_ps),
        "ks_statistic": float(ks_stats["ks_statistic"]),
        "ks_p_value": float(ks_stats["p_value"]),
        "flag_extreme_ps": bool(prop_ps_extreme > 0),
    }

compute_weight_diagnostics

compute_weight_diagnostics(df: DataFrame, ps_col: str = PS_COL, treatment_col: str = TREATMENT_COL, weight_type: Literal['ATE', 'ATT'] = 'ATE', clip_percentile: float = 1) -> dict

Compute IPW weight diagnostics: effective sample size and weight summaries.

Weights are computed with compute_ipw_weights; clip_percentile=1 (default) means raw, unclipped weights, matching the estimators' default. Propensity scores of exactly 0 or 1 raise ValueError: their IPW weights are undefined and would only reflect the numerical stabilizer.

Parameters:

Name	Type	Description	Default
df	DataFrame	Input DataFrame with treatment and propensity score columns.	required
ps_col	str	Name of the propensity score column.	PS_COL
treatment_col	str	Name of the treatment status column (1 treated, 0 control).	TREATMENT_COL
weight_type	Literal['ATE', 'ATT']	"ATE" or "ATT".	'ATE'
clip_percentile	float	Upper-tail clipping passed to compute_ipw_weights.	1

Returns:

Type	Description
dict	dict with per-arm ESS, ESS as a fraction of arm size, and per-arm weight
dict	summaries (max, mean, 95th and 99th percentile). Pooled equivalents
dict	(ess_total, ess_fraction_total, max_weight, mean_weight, weight_q95,
dict	weight_q99) are included only for weight_type="ATE": ATT weights put the
dict	two arms on different scales (treated weights are identically 1), so
dict	pooled Kish ESS and pooled quantiles would be misleading there.

Source code in CausalEstimate/diagnostics/weights.py

def compute_weight_diagnostics(
    df: pd.DataFrame,
    ps_col: str = PS_COL,
    treatment_col: str = TREATMENT_COL,
    weight_type: Literal["ATE", "ATT"] = "ATE",
    clip_percentile: float = 1,
) -> dict:
    """
    Compute IPW weight diagnostics: effective sample size and weight summaries.

    Weights are computed with compute_ipw_weights; clip_percentile=1 (default)
    means raw, unclipped weights, matching the estimators' default.
    Propensity scores of exactly 0 or 1 raise ValueError: their IPW weights are
    undefined and would only reflect the numerical stabilizer.

    Args:
        df: Input DataFrame with treatment and propensity score columns.
        ps_col: Name of the propensity score column.
        treatment_col: Name of the treatment status column (1 treated, 0 control).
        weight_type: "ATE" or "ATT".
        clip_percentile: Upper-tail clipping passed to compute_ipw_weights.

    Returns:
        dict with per-arm ESS, ESS as a fraction of arm size, and per-arm weight
        summaries (max, mean, 95th and 99th percentile). Pooled equivalents
        (ess_total, ess_fraction_total, max_weight, mean_weight, weight_q95,
        weight_q99) are included only for weight_type="ATE": ATT weights put the
        two arms on different scales (treated weights are identically 1), so
        pooled Kish ESS and pooled quantiles would be misleading there.
    """
    n_treated, n_control = validate_ps_and_treatment(df, ps_col, treatment_col)
    A = df[treatment_col].to_numpy()
    ps = df[ps_col].to_numpy()

    W = compute_ipw_weights(
        A, ps, weight_type=weight_type, clip_percentile=clip_percentile
    )
    W_treated = W[A == 1]
    W_control = W[A == 0]
    ess_treated = compute_ess(W_treated)
    ess_control = compute_ess(W_control)

    result = {
        "n_total": int(len(W)),
        "n_treated": n_treated,
        "n_control": n_control,
        "ess_treated": ess_treated,
        "ess_control": ess_control,
        "ess_fraction_treated": min(ess_treated / n_treated, 1.0),
        "ess_fraction_control": min(ess_control / n_control, 1.0),
        **_weight_summary(W_treated, "treated"),
        **_weight_summary(W_control, "control"),
    }
    if weight_type == "ATE":
        ess_total = compute_ess(W)
        result.update(
            {
                "ess_total": ess_total,
                "ess_fraction_total": min(ess_total / len(W), 1.0),
                "max_weight": float(W.max()),
                "mean_weight": float(W.mean()),
                "weight_q95": float(np.percentile(W, 95)),
                "weight_q99": float(np.percentile(W, 99)),
            }
        )
    return result

compute_ess

compute_ess(weights: ndarray) -> float

Effective sample size of a weighted sample (Kish): (sum w)^2 / sum(w^2).

Equals n for uniform weights and approaches 1 as a single weight dominates. Defined for nonnegative weights; nonfinite, negative, or all-zero weights raise ValueError.

Source code in CausalEstimate/diagnostics/weights.py

def compute_ess(weights: np.ndarray) -> float:
    """
    Effective sample size of a weighted sample (Kish): (sum w)^2 / sum(w^2).

    Equals n for uniform weights and approaches 1 as a single weight dominates.
    Defined for nonnegative weights; nonfinite, negative, or all-zero weights
    raise ValueError.
    """
    w = np.asarray(weights, dtype=float)
    if w.size == 0:
        raise ValueError("weights must be non-empty.")
    if not np.all(np.isfinite(w)):
        raise ValueError("weights must be finite (no NaN or inf).")
    if np.any(w < 0):
        raise ValueError("weights must be nonnegative.")
    if w.max() == 0:
        raise ValueError("weights must not be all zero.")
    w = w / w.max()  # ESS is scale-invariant; normalizing prevents overflow
    return float(w.sum() ** 2 / (w**2).sum())

compute_ipw_weights

compute_ipw_weights(A: ndarray, ps: ndarray, weight_type: Literal['ATE', 'ATT'] = 'ATE', clip_percentile: float = 1, eps: float = 1e-09) -> np.ndarray

Inverse Propensity Score weights with optional per-arm clipping.

Formulas: - ATE: w = A/ps + (1-A)/(1-ps) - ATT: w = A + (1-A) * ps/(1-ps)

Note these are the IPW estimator's weights, which differ from the TMLE targeting weights of compute_arm_weights: the ATT variants are normalised differently (1 vs 1/p_treated on the treated arm), and here the two arms are summed into a single vector rather than kept separate.

Parameters:

Name	Type	Description	Default
A	ndarray	Binary treatment assignment vector (1 for treated, 0 for control).	required
ps	ndarray	Propensity score vector.	required
weight_type	Literal['ATE', 'ATT']	"ATE" or "ATT".	'ATE'
clip_percentile	float	Upper quantile at which to clip, as a PROPORTION in (0, 1]. 1.0 (default) applies no clipping; 0.99 clips the top 1%.	1
eps	float	Small constant added to denominators.	1e-09

Returns:

Type	Description
ndarray	An array of IPW weights.

Raises:

Type	Description
ValueError	on an invalid weight_type, mismatched shapes, or a clip_percentile outside (0, 1].

Source code in CausalEstimate/estimators/functional/utils.py

def compute_ipw_weights(
    A: np.ndarray,
    ps: np.ndarray,
    weight_type: Literal["ATE", "ATT"] = "ATE",
    clip_percentile: float = 1,
    eps: float = 1e-9,
) -> np.ndarray:
    """
    Inverse Propensity Score weights with optional per-arm clipping.

    Formulas:
    - ATE: w = A/ps + (1-A)/(1-ps)
    - ATT: w = A + (1-A) * ps/(1-ps)

    Note these are the IPW estimator's weights, which differ from the TMLE
    targeting weights of `compute_arm_weights`: the ATT variants are
    normalised differently (1 vs 1/p_treated on the treated arm), and here the
    two arms are summed into a single vector rather than kept separate.

    Args:
        A: Binary treatment assignment vector (1 for treated, 0 for control).
        ps: Propensity score vector.
        weight_type: "ATE" or "ATT".
        clip_percentile: Upper quantile at which to clip, as a PROPORTION in
            (0, 1]. 1.0 (default) applies no clipping; 0.99 clips the top 1%.
        eps: Small constant added to denominators.

    Returns:
        An array of IPW weights.

    Raises:
        ValueError: on an invalid `weight_type`, mismatched shapes, or a
            `clip_percentile` outside (0, 1].
    """
    if weight_type not in ("ATE", "ATT"):
        raise ValueError("weight_type must be 'ATE' or 'ATT'")
    if not (0 < clip_percentile <= 1.0):
        raise ValueError("clip_percentile must be in the interval (0, 1.0].")
    if A.shape != ps.shape:
        raise ValueError("A and ps must have the same shape.")
    check_ps_not_exact_zero_one(ps)

    if weight_type == "ATE":
        weights = A / (ps + eps) + (1 - A) / (1 - ps + eps)
    else:  # ATT: treated units have weight 1.
        weights = A + (1 - A) * ps / (1 - ps + eps)

    if clip_percentile < 1.0:
        q = clip_percentile * 100
        treated_mask = A == 1
        for mask in (treated_mask, ~treated_mask):
            if np.any(mask):
                threshold = np.percentile(weights[mask], q)
                weights[mask] = np.clip(weights[mask], a_min=None, a_max=threshold)

    return weights

compute_evalue

compute_evalue(estimate: float, ci_lower: Optional[float] = None, ci_upper: Optional[float] = None, scale: Literal['RR', 'RD'] = 'RR', baseline_risk: Optional[float] = None) -> dict

E-value for unmeasured confounding (VanderWeele & Ding, 2017).

The E-value is the minimum strength of association, on the risk-ratio scale, that an unmeasured confounder would need with both treatment and outcome (conditional on measured covariates) to fully explain away the observed effect. For RR >= 1, E = RR + sqrt(RR * (RR - 1)); protective effects use 1 / RR. The CI E-value applies the same formula to the confidence bound closest to the null and is 1 when the interval contains the null.

Parameters:

Name	Type	Description	Default
estimate	float	Effect estimate on the given scale.	required
ci_lower	Optional[float]	Lower confidence bound; pass together with ci_upper.	None
ci_upper	Optional[float]	Upper confidence bound.	None
scale	Literal['RR', 'RD']	"RR" for risk ratios (RR, RRT effect types), or "RD" for risk differences (ATE, ARR, ATT, ATC on a binary outcome). RD values are converted to risk ratios via (baseline_risk + RD) / baseline_risk.	'RR'
baseline_risk	Optional[float]	Outcome risk in the untreated (e.g. the estimator's effect_0). Required when scale="RD".	None

Returns:

Type	Description
dict	dict with keys "evalue" and "evalue_ci" (None if no CI was given).

Source code in CausalEstimate/diagnostics/sensitivity.py

def compute_evalue(
    estimate: float,
    ci_lower: Optional[float] = None,
    ci_upper: Optional[float] = None,
    scale: Literal["RR", "RD"] = "RR",
    baseline_risk: Optional[float] = None,
) -> dict:
    """
    E-value for unmeasured confounding (VanderWeele & Ding, 2017).

    The E-value is the minimum strength of association, on the risk-ratio scale,
    that an unmeasured confounder would need with both treatment and outcome
    (conditional on measured covariates) to fully explain away the observed
    effect. For RR >= 1, E = RR + sqrt(RR * (RR - 1)); protective effects use
    1 / RR. The CI E-value applies the same formula to the confidence bound
    closest to the null and is 1 when the interval contains the null.

    Args:
        estimate: Effect estimate on the given scale.
        ci_lower: Lower confidence bound; pass together with ci_upper.
        ci_upper: Upper confidence bound.
        scale: "RR" for risk ratios (RR, RRT effect types), or "RD" for risk
            differences (ATE, ARR, ATT, ATC on a binary outcome). RD values are
            converted to risk ratios via (baseline_risk + RD) / baseline_risk.
        baseline_risk: Outcome risk in the untreated (e.g. the estimator's
            `effect_0`). Required when scale="RD".

    Returns:
        dict with keys "evalue" and "evalue_ci" (None if no CI was given).
    """
    if scale not in {"RR", "RD"}:
        raise ValueError("scale must be 'RR' or 'RD'.")
    if (ci_lower is None) != (ci_upper is None):
        raise ValueError("ci_lower and ci_upper must be given together.")
    if scale == "RD":
        if baseline_risk is None:
            raise ValueError("baseline_risk is required when scale='RD'.")
        if not 0 < baseline_risk <= 1:
            raise ValueError("baseline_risk must be in (0, 1].")
    if ci_lower is not None and ci_lower > ci_upper:
        raise ValueError("ci_lower must not exceed ci_upper.")

    rr = _to_rr(estimate, scale, baseline_risk)
    result = {"evalue": _evalue_from_rr(rr), "evalue_ci": None}
    if ci_lower is None:
        return result

    lo = _to_rr(ci_lower, scale, baseline_risk)
    hi = _to_rr(ci_upper, scale, baseline_risk)
    if lo <= 1 <= hi:
        result["evalue_ci"] = 1.0
    else:
        result["evalue_ci"] = _evalue_from_rr(lo if rr > 1 else hi)
    return result