Principal Stratification and Survivor Average Causal Effect

Principal stratification

classifies patients not by what was observed but by their joint potential outcomes on an intermediate, post-baseline event under every treatment arm being compared. For a binary intermediate event S (survive vs. die by the assessment time, or persist vs. discontinue) and two arms, each patient belongs to exactly one of four latent principal strata: always-survivors (would survive under both arms), protected (would survive only under the arm they are assigned, e.g. treatment), harmed (would survive only under control), and never-survivors (would die under both). Because principal strata are defined from potential outcomes under both arms simultaneously, they are, by construction, unaffected by treatment assignment -- a genuine baseline characteristic, even though S itself is a post-baseline variable. This is what makes stratifying on them a legitimate causal operation where stratifying on the observed value of S is not. The survivor average causal effect (SACE) -- the term coined by Hayden, Pauler & Schoenfeld (2005) and developed within the principal-stratification framework by Rubin (2006) -- names the treatment effect on a downstream outcome Y (quality of life, functional status, a biomarker, a patient-reported outcome), restricted to the always-survivor stratum: it is the canonical instance, E[Y(1) − Y(0) | always-survivor]. ICH E9(R1) (2019) formally lists principal stratum as the fifth intercurrent-event (ICE) strategy alongside treatment policy, while-on-treatment, hypothetical, and composite, for exactly this situation: an ICE (death, discontinuation) that makes the outcome of interest undefined, not merely unobserved, for some patients.

Why the naive survivor-only comparison is biased

Restricting the analysis to observed survivors and comparing mean Y between arms conditions on S = "survived," a variable that is itself caused by treatment. If treatment affects survival, the two observed-survivor groups are different populations -- a coincidence of comparable baseline covariates does not make them comparable, because survival is a collider on the causal path from unmeasured severity to both the treatment-survival relationship and Y. Concretely: a drug that reduces mortality does so by keeping alive some patients who were sicker at baseline than the patients who would have survived on the comparator anyway (the "protected" stratum). The treated survivor pool is therefore diluted with these sicker "rescued" patients, and a naive treated-survivors-vs-control-survivors comparison of Y is not a causal contrast at all: the control side is the always-survivor control mean, E[Y(0) | always-survivor], but the treated side is a mixture of the always-survivor treated mean and the protected stratum's treated mean, π·E[Y(1) | always-survivor] + (1−π)·E[Y(1) | protected], where π (the always-survivor share of observed treated survivors) is itself unidentified without further assumptions. Because the protected stratum's control-arm outcome Y(0) is undefined -- they die under control -- there is no common population underlying the naive comparison, so it is not an average of two well-defined causal effects; it is a difference between a stratum-specific mean and an unidentified mixture. It can flip sign entirely from the SACE (the drug "looks worse" on quality of life among survivors purely because it saved sicker patients) -- this is truncation-by-death bias, and it is structurally identical to the collider/selection bias that drives immortal-time and healthy-survivor artifacts elsewhere in the catalog, but here there is no fix by better covariate adjustment: Y(1) and Y(0) are not defined for a patient who dies under that arm, so no amount of regression adjustment recovers them. The problem is not missing data; it is missing potential outcomes.

Identification is the hard part

Principal strata are latent -- no one observes both S(1) and S(0) for the same patient -- so SACE is not point-identified without further assumptions. The near-universal starting assumption is monotonicity: treatment does not decrease the probability of the favorable-direction intermediate outcome for any unit -- formally S(1) ≥ S(0) for survival, i.e. no patient who would survive on control dies on treatment -- which rules out the "harmed" stratum and lets the always-survivor proportion be identified as the smaller of the two arm-specific survival rates. Even under monotonicity, however, the treated-survivor pool remains a mixture of always-survivors and the unidentified "protected" stratum, so estimation still requires one of: (1) nonparametric bounds (Zhang & Rubin 2003) that trim the most extreme observed outcomes from the mixed arm to bracket the always-survivor mean without further assumptions -- wide but assumption-light; (2) a sensitivity-parameter model (Chiba & VanderWeele 2011) that posits a single interpretable parameter for how the protected stratum's outcome differs from the always-survivor mean and reports SACE as a tipping-point function of that parameter; or (3) a Bayesian mixture / latent-class model -- following the joint S-and-Y latent-class modeling template that Hirano, Imbens, Rubin & Zhou (2000) developed for principal stratification in a noncompliance/encouragement-design setting, later adapted directly to truncation-by-death outcomes (e.g. Tong et al. 2023) -- that jointly models S and Y with baseline covariates, estimating posterior strata-membership probabilities and outcome distributions directly -- narrower intervals than bounds, at the cost of parametric assumptions about the mixture components that are themselves untestable. All three approaches identify SACE from arm-level survival rates and arm-level outcome summaries that are exchangeable across arms by design -- guaranteed by randomization in an RCT, but requiring a separate identification layer in observational RWE data (consistency, positivity, and no-unmeasured-confounding for the treatment assignment itself, on top of monotonicity for the intermediate event): arm-level survival rates and outcome means must first be adjusted for baseline confounders (e.g., via propensity weighting or standardization) before being fed into the bounds/sensitivity/mixture machinery below -- see Egleston et al. (2009) for the observational-data version of this identification problem, and the worked example's Estimation step for how the adjustment layer and the SACE layer compose.

Pros, cons, and trade-offs

(specific and comparative, naming the alternatives).

• vs the hypothetical ICE strategy: Both target an effect "as if" the ICE had not disrupted the outcome, but they ask different questions and require different machinery. The hypothetical strategy (IPCW, g-methods) asks what Y would have been had the patient not discontinued or died -- a coherent question for a "soft" ICE like discontinuation, where the counterfactual undisrupted outcome plausibly exists. For death, that counterfactual quality-of-life value does not exist even in principle; you cannot hypothesize a dead patient's PROMIS score. Principal stratification sidesteps this by restricting the population, not imputing the outcome: it asks about the effect among patients for whom Y(1) and Y(0) both exist, which is philosophically cleaner for death but answers a narrower, sub-population question rather than the full-population "no ICE" counterfactual. Use principal stratum, not hypothetical, whenever the ICE makes the outcome structurally undefined rather than merely unmeasured.
• vs the composite strategy: Composite (fold death into a "failure," e.g., death-or-poor-function) is trivially identified from observed data and answers a genuine net-clinical-benefit question, but it collapses two distinct mechanisms -- a survival effect and a functional effect -- into one number. A drug that saves severely impaired patients can score worse on a composite than a drug that saves no one, because the composite treats every saved-but-impaired survivor as equivalent to every death. Principal stratification keeps the survival result and the always-survivor functional result as two separate, cleanly interpretable numbers. Prefer composite when the true decision genuinely trades survival against function on one scale (e.g., a single QALY-style summary is the deliverable); prefer principal stratum when the audience needs to see the mortality effect and the functional effect separately, which is most RWE and most regulatory submissions.
• vs a naive complete-case (survivor-only) analysis: The naive analysis is trivial to compute and is exactly the biased estimator described above. Principal stratification costs identifiability (bounds, a sensitivity parameter, or a parametric mixture model) in exchange for a well-defined causal target. Never report a survivor-only comparison as a treatment effect without at least characterizing how differential the survival rates were and whether SACE bounds meaningfully differ from the naive number, as they do whenever the two arms' survival rates diverge and any subgroup's outcome differs systematically by severity.
• vs treatment-policy with worst-outcome imputation: Assigning a fixed floor value ("if dead, function = 0" or a worst-rank score) to deceased patients sidesteps identifiability the way the composite strategy does, and is common in oncology PRO analyses (e.g., PFS-adjusted quality-of-life scoring). It is defensible when the audience truly treats death as the worst possible functional outcome, but it silently redefines the estimand's scale and is highly sensitive to the arbitrary floor value chosen. Principal stratification avoids inventing a number for an outcome that does not exist, at the cost of describing a narrower population.

When NOT to use -- and when it is actively misleading

• Mortality/discontinuation barely differs between arms. If the two arms' survival rates are close to each other (not necessarily close to 100%), the estimated protected-stratum share, |S_high − S_low|, is small, so under a defensible monotonicity assumption nearly all observed survivors are always-survivors and the naive survivor-only estimate converges with the SACE bounds -- report the simpler complete-case result with a note on how similar the survival rates were and the estimated protected fraction. This does not mean the always-survivor stratum is close to the whole population: if both arms survive at only 20%, the always-survivor stratum is still roughly 20% of all patients, just a large share of the much smaller surviving subgroup. And without a defensible monotonicity assumption, equal marginal survival rates can mask equal-sized, offsetting protected and harmed strata within each arm, in which case survivor-only bias need not shrink at all -- close arm-level survival rates only license the "close enough" simplification once monotonicity itself is justified.
• Monotonicity is implausible. In comparative-safety or heterogeneous-effect settings, some patients may genuinely be harmed by the "beneficial" arm even while the average effect favors it (e.g., a drug that reduces overall mortality but a rare fatal idiosyncratic reaction in a subgroup). Sharp Zhang-Rubin bounds without monotonicity are frequently uninformative (they can span the entire outcome range). Monotonicity is a within-patient statement about two potential outcomes and is not directly testable from observed data: non-crossing arm-level or subgroup-level survival curves are consistent with monotonicity but do not verify it (a crossing curve rules monotonicity out; a non-crossing curve does not rule it in, since offsetting individual-level harm and benefit within a subgroup can still produce a non-crossing average curve). State monotonicity as an explicit, scientifically justified assumption -- grounded in mechanism of action and prior evidence, not curve-checking alone -- and pair it with a sensitivity analysis that allows a nonzero harmed stratum, rather than treating it as a default.
• The policy question is about the full population, including those who die. A payer or regulator asking "what happens to all patients assigned this drug, mortality included" needs the survival curve and, if a single summary is required, a composite or QALY-style estimand -- not SACE, which by design describes only the always-survivor subgroup and is silent about the patients it excludes.
• Reporting SACE without also reporting the survival difference. SACE deliberately discards the mortality result to isolate the functional effect; presenting the SACE point estimate or bounds alone, without the underlying survival curves and the always-survivor stratum's estimated size, hides exactly the information a reader needs to judge whether "the drug that saves more people who function well" is actually the better drug overall.
• Principal-strata membership is well predicted by rich baseline covariates but the analysis ignores them. Reporting only unconditional Zhang-Rubin bounds when covariate-sharpened bounds (Long & Hudgens 2013) built from baseline severity could narrow them substantially wastes information routinely available in claims/EHR data (comorbidity burden, prior utilization, frailty indices).

RWE-specific extension: principal strata beyond literal death

A structurally analogous but not identical problem arises whenever an outcome is only measured conditional on a post-baseline behavioral event other than death -- most commonly persistence/tolerability: pain scores, PROMIS function, or a lab value collected only among patients who remain on therapy long enough for the assessment. The key difference from literal death: a discontinued patient's pain score or function still exists in principle -- it is merely unobserved, not structurally undefined the way a dead patient's PROMIS score is. That makes the persistence case closer to a missing-data problem with a hypothetical-strategy remedy (IPCW, multiple imputation targeting the off-treatment value) than to a true truncation-by-death problem. Choosing the principal-stratum framing anyway -- defining an "always-persister" (or "always-tolerator") stratum of patients who would tolerate and stay on treatment regardless of which of the two drugs they were assigned, and estimating the effect only within it -- is a distinct scientific choice, not a mechanical transfer of the SACE machinery: it answers "among patients who can tolerate either drug, which works better while they are on it" rather than "what would this patient's outcome have been had they stayed on therapy." It requires its own persistence-measurement definition (fill-gap rules), its own monotonicity-in-discontinuation-risk justification (often less plausible than monotonicity in mortality, since adverse-event profiles differ by mechanism rather than uniformly by severity), and its own exchangeability assumptions -- borrowing the Zhang-Rubin/sensitivity-parameter/Bayesian-mixture machinery is legitimate, but borrowing the justification ("the outcome does not exist, so restrict the population") from the death case is not. This still answers a question RWE stakeholders ask constantly and rarely formalize correctly: "among patients who can tolerate either drug, which works better?" -- as opposed to the naive, badly confounded comparison of outcomes among whoever happened to still be filling prescriptions at the assessment date in each arm.

Data-source operational depth

(claims vs EHR vs registry vs linked).

• Claims: Vital status requires a linked death source -- state vital records, the CMS Master Beneficiary Summary File / denominator file, or the National Death Index; claims-inferred death (inpatient discharge disposition "expired," a hospice episode with no subsequent claims) is noisy and systematically misses out-of-network deaths. Failure mode: treating disenrollment as death (or vice versa) corrupts strata membership at the source -- Medicare Advantage plans in particular have historically lagged or incomplete death-date capture without a denominator-file linkage. PRO/functional outcomes are rarely in claims; where they exist (IRF-PAI, OASIS functional items, some Medicare Advantage encounter supplemental data), construct validity is weaker than a dedicated instrument, and a claims-only study may need a utilization-based proxy (SNF admission, DME order) with its own measurement-error caveats.
• EHR: Mortality capture is incomplete outside the health system -- a patient who dies elsewhere looks like loss to follow-up or an ongoing "alive, just not seen recently" patient, which contaminates the "alive" stratum with unknown deaths and generally biases the observed survival-rate estimate upward (deaths misclassified as censored/alive inflate the apparent survival rate, they do not deflate it); the size of that inflation, and therefore of any resulting strata-membership error, depends on linkage/capture completeness and can differ by arm if that completeness itself differs by severity or by health-system retention pattern. Link to state vital records or SSA where possible. PRO instruments (PROMIS, EQ-5D) are more often structured in EHR, but non-response among the frailest survivors creates a second truncation-by-selection problem layered on top of truncation-by-death -- survivors who are too sick to complete the questionnaire are missing not at random, which the Bayesian mixture approach can be extended to model jointly (a three-state S-Y-response model) but nonparametric bounds cannot handle without an additional missing-data assumption.
• Registry: Frequently the strongest source -- adjudicated death dates and structured, prespecified functional or PRO assessments -- but still subject to survivor non-response, and typically requires linkage to claims for complete longitudinal drug-exposure history to define the persistence-based principal strata described above.
• Linked claims-EHR-vital records: The ideal substrate (reliable mortality + EHR severity for principal-score covariates + claims completeness for exposure), at the cost of reconciling death dates across sources before fixing the strata definition, and of linkage selection potentially differing by the same severity that drives strata membership.

Worked claims-linked example (oncology)

Question: does a newer targeted therapy improve 6-month patient-reported physical function (PROMIS Physical Function T-score) versus chemotherapy, in a linked claims-EHR-registry oncology cohort, given the targeted therapy also improves 6-month survival? (1) Eligibility: adults with incident metastatic disease initiating first-line therapy, continuous enrollment, no prior use of either regimen class. (2) The ICE: death before the 6-month PROMIS assessment -- PROMIS is structurally undefined for a deceased patient, not merely missing. (3) Monotonicity justification: prior randomized evidence for this drug class shows no crossing of the survival curves in any prespecified subgroup, together with a mechanism of action with no known mortality-increasing subpopulation; this supports "the targeted therapy does not increase near-term mortality risk relative to chemotherapy for any patient" as a scientifically defensible assumption -- non-crossing curves are consistent with monotonicity, not proof of it (monotonicity is a within-patient statement that observed curves cannot verify), so the analysis plan pre-specifies a sensitivity check allowing a small harmed stratum. (4) Strata: always-survivors = min(6-month survival under targeted therapy, 6-month survival under chemotherapy), by monotonicity equal to the chemotherapy arm's survival rate; protected = the survival-rate difference; never-survivors = one minus the targeted-therapy survival rate. (5) Estimation: because this is an observational claims-EHR-registry cohort rather than a randomized trial, first fit an exposure model (e.g., propensity score or entropy-balancing weights on baseline ECOG status, comorbidity burden, disease extent, and prior treatment history) so the two arms' survival rates and outcome distributions are exchangeable; compute Zhang-Rubin trimming bounds on the weighted targeted-therapy arm's observed PROMIS distribution (assumption-light beyond the exchangeability weighting), then a Chiba-VanderWeele sensitivity-parameter tipping-point curve incorporating the same baseline covariates as a principal score to narrow the bounds; report both. (6) Reporting: present the 6-month survival curves by arm and the SACE bounds/tipping-point result side by side -- never the functional result alone -- so a reader can judge the mortality benefit and the functional benefit among likely survivors as two distinct, honestly-scoped findings.

Interpreting the output

In a Drug A (targeted therapy) versus Drug B (chemotherapy) study (200 patients per arm, 6-month follow-up) -- the same arm proportions as the 10-patient-per-arm dataset in the worked example below, scaled up 20-fold so every observed QoL score recurs 20 times at each arm and rank -- 6-month survival is 80% (160/200) on Drug A versus 60% (120/200) on Drug B. Mean PROMIS physical-function score among observed survivors is 72 on Drug A (160 survivors) versus 65 on Drug B (120 survivors) -- a naive survivor-only difference of +7.0 points favoring Drug A.

(1) Formal interpretation. Under monotonicity (Drug A's survival benefit means it does not increase death risk for any patient), the always-survivor stratum comprises 60% of each arm (120 of every 200 patients) -- exactly the patients who would survive under Drug B, all of whom are therefore always-survivors. The 160 observed Drug A survivors are a mixture: 120 always-survivors plus 40 patients "protected" by Drug A (they would have died on Drug B) whose identity in the data is unknown. Because every one of the 8 distinct Drug A QoL scores and 6 distinct Drug B QoL scores in the worked example below recurs exactly 20 times at this scale, trimming the 40/160 = 25% most extreme PROMIS scores from the Drug A survivor distribution in each direction reproduces, point for point, the same trimmed means as the 10-patient-per-arm calculation: SACE bounds of [+3.5, +11.7] points -- a genuinely different quantity from the naive +7.0, which averages the always-survivor contrast with an unidentified mixture component. See the worked example below for the row-level arithmetic that produces these bounds. The true SACE could be anywhere in that interval without further assumptions; a sensitivity-parameter or principal-score model is needed to narrow it to a point estimate.

(2) Practical interpretation. A clinician asking "if my patient is going to survive six months either way, will Drug A give better function than chemotherapy?" needs the SACE bounds (a genuine functional advantage of at least +3.5 points, arguably more), not the naive +7.0, which partly reflects the fact that Drug A's survivor pool includes sicker rescued patients who pull the average down, not up -- so the naive number is not even a conservative estimate; its relationship to the true SACE is unknown without the strata calculation. Reporting the 20-percentage-point survival advantage alongside the SACE bounds, rather than the functional result alone, gives the complete, honestly-scoped picture: more patients alive at 6 months, and among the patients who would be alive regardless of arm, meaningfully better physical function on Drug A.

Decision diagram

flowchart TD
  Pop[Eligible patient at baseline] --> Y1{Would survive<br/>under Treatment?}
  Y1 -->|Yes| Y0a{Would survive<br/>under Control?}
  Y1 -->|No| Y0b{Would survive<br/>under Control?}
  Y0a -->|Yes| LL[Always-survivor LL<br/>outcome defined under both arms<br/>= SACE target stratum]
  Y0a -->|No| LD[Protected-by-treatment LD<br/>outcome defined only under Treatment<br/>unidentified without extra assumptions]
  Y0b -->|Yes| DL[Harmed-by-treatment DL<br/>ruled out under monotonicity]
  Y0b -->|No| DD[Never-survivor DD<br/>outcome undefined under either arm]
  style LL fill:#e6ffe6
  style DD fill:#ffe6e6
  style DL fill:#ffe6e6,stroke-dasharray: 5 5

flowchart TD
  Q[Outcome truncated by death<br/>or by non-persistence] --> M{Monotonicity plausible?<br/>treatment cannot increase<br/>death/discontinuation risk}
  M -->|No| Warn[Sharp bounds likely uninformative;<br/>report survival curve + composite estimand instead]
  M -->|Yes| B{Baseline covariates predict<br/>strata membership well?}
  B -->|Weak or none| Bounds[Nonparametric Zhang-Rubin<br/>trimming bounds]
  B -->|Strong| Score[Principal-score adjusted<br/>bounds or mixture model]
  Bounds --> Wide{Bounds too wide<br/>to be decision-useful?}
  Wide -->|Yes| Sens[Sensitivity-parameter<br/>tipping-point model]
  Wide -->|No| Report[Report SACE bounds<br/>alongside the survival result]
  Sens --> Report
  Score --> Bayes[Bayesian mixture / latent-class model]
  Bayes --> Report
style M fill:#e6f3ff
style B fill:#e6f3ff

Worked example

Scenario

A researcher compares two heart-failure drugs (Drug A, the newer therapy, and Drug B, standard of care) using a claims-linked cohort of 10 patients per arm. The outcome is a 6-month quality-of-life (QoL) functional score from 0-100, only measured in patients who survive to 6 months. Drug A has better 6-month survival than Drug B. The researcher wants to know whether Drug A also gives better QoL among patients who would survive regardless of which drug they took -- not just among whoever happens to still be alive in each arm.

Dataset

6-month status and QoL functional score for 10 patients per arm. QoL is blank (undefined, not missing) for patients who died before the 6-month assessment.

Steps

Result

SACE lies in the bound [+3.5, +11.7] QoL points among the always-survivor stratum, under monotonicity -- a genuinely different quantity from the naive survivor-only difference of +7.0 points, which mixes the always-survivor contrast with an unidentified 2-patient "rescued" subgroup of unknown QoL. Even the lower bound (+3.5) is positive, so the functional advantage among always-survivors is directionally robust here, but the exact magnitude is not identified without a sensitivity-parameter or principal-score model.

Strata Bound Table

How the 8 observed Drug A survivors decompose into the always-survivor and protected strata, and the resulting SACE bounds versus the naive comparison.

Trade-offs

Runnable example

import numpy as np
import pandas as pd


def _weighted_trim_mean(values, weights, trim_frac, drop):
    """Mean of the weighted distribution after dropping the `trim_frac` share of total WEIGHT from the
    low or high tail (drop='low' or 'high'). Weight-based (not row-count-based) trimming lets trim_frac
    be a non-integer share of the sample -- required once rows carry unequal IPTW weights, and more
    honest than rounding to the nearest observation even when every weight equals 1."""
    values = np.asarray(values, dtype=float)
    weights = np.asarray(weights, dtype=float)
    order = np.argsort(values)
    values, weights = values[order], weights[order]
    cum_weight = np.cumsum(weights)
    total_weight = cum_weight[-1]
    if drop == "low":
        keep = cum_weight > trim_frac * total_weight
    else:
        keep = cum_weight <= (1 - trim_frac) * total_weight
    if not keep.any():
        return np.nan
    return np.average(values[keep], weights=weights[keep])


def sace_bounds_monotonicity(df, arm_col="arm", higher_surv_arm="A", lower_surv_arm="B",
                              survive_col="survived", outcome_col="outcome", weight_col=None):
    """Zhang & Rubin (2003) nonparametric trimming bounds for SACE under monotonicity.
    Assumes higher_surv_arm has survival(A) >= survival(B); raises otherwise (check direction first).

    weight_col: None for an RCT (all rows weight 1). For observational data, the name of a column
    holding pre-fit IPTW weights (see the implementation description above) -- this function does not
    fit the exposure model itself, it only consumes already-exchangeable weighted summaries."""
    if not set(df[survive_col].dropna().unique()) <= {0, 1}:
        raise ValueError(f"{survive_col} must be binary (0/1).")

    w = pd.Series(1.0, index=df.index) if weight_col is None else df[weight_col].astype(float)

    hi_mask = df[arm_col] == higher_surv_arm
    lo_mask = df[arm_col] == lower_surv_arm
    s_hi = np.average(df.loc[hi_mask, survive_col], weights=w.loc[hi_mask])
    s_lo = np.average(df.loc[lo_mask, survive_col], weights=w.loc[lo_mask])
    if s_hi < s_lo:
        raise ValueError("Monotonicity direction violated: higher_surv_arm must have survival >= lower_surv_arm.")
    if s_hi <= 0:
        raise ValueError("No survivors observed in higher_surv_arm; cannot estimate SACE bounds.")

    pi_always = s_lo                              # P(always-survivor) = P(survive lower-survival arm)
    pi_protected = s_hi - s_lo                     # unidentified subgroup within the higher-survival arm's survivors
    frac_protected = pi_protected / s_hi

    hi_surv_mask = hi_mask & (df[survive_col] == 1)
    hi_outcome = df.loc[hi_surv_mask, outcome_col].to_numpy()
    hi_weight = w.loc[hi_surv_mask].to_numpy()
    if hi_outcome.size == 0:
        raise ValueError("No observed survivors with a non-missing outcome in higher_surv_arm.")

    upper_mean_hi = _weighted_trim_mean(hi_outcome, hi_weight, frac_protected, drop="low")   # protected assumed lowest-scoring -> trim bottom
    lower_mean_hi = _weighted_trim_mean(hi_outcome, hi_weight, frac_protected, drop="high")  # protected assumed highest-scoring -> trim top

    lo_surv_mask = lo_mask & (df[survive_col] == 1)
    lo_mean = np.average(df.loc[lo_surv_mask, outcome_col], weights=w.loc[lo_surv_mask])
    naive_mean = np.average(hi_outcome, weights=hi_weight)

    return dict(
        survival_higher_arm=s_hi, survival_lower_arm=s_lo,
        pi_always_survivor=pi_always, pi_protected=pi_protected, frac_protected=frac_protected,
        sace_lower_bound=lower_mean_hi - lo_mean,
        sace_upper_bound=upper_mean_hi - lo_mean,
        naive_survivor_only_diff=naive_mean - lo_mean,
    )


def sace_sensitivity_tipping_point(df, arm_col="arm", higher_surv_arm="A", lower_surv_arm="B",
                                    survive_col="survived", outcome_col="outcome",
                                    weight_col=None, delta_grid=None):
    """Chiba & VanderWeele (2011)-style sensitivity parameter: assume the unidentified 'protected'
    stratum's mean outcome differs from the always-survivor mean (within the higher-survival arm) by
    a fixed delta. Solve the mixture equation for the implied always-survivor mean and SACE at each
    delta on a grid, producing a tipping-point curve rather than a single bound. weight_col as above."""
    if delta_grid is None:
        delta_grid = np.linspace(-20, 20, 41)
    w = pd.Series(1.0, index=df.index) if weight_col is None else df[weight_col].astype(float)

    hi_mask = df[arm_col] == higher_surv_arm
    lo_mask = df[arm_col] == lower_surv_arm
    s_hi = np.average(df.loc[hi_mask, survive_col], weights=w.loc[hi_mask])
    s_lo = np.average(df.loc[lo_mask, survive_col], weights=w.loc[lo_mask])
    if s_hi <= 0:
        raise ValueError("No survivors observed in higher_surv_arm; cannot run the sensitivity analysis.")
    w_always = s_lo / s_hi                              # fraction of higher-arm survivors who are always-survivors

    hi_surv_mask = hi_mask & (df[survive_col] == 1)
    lo_surv_mask = lo_mask & (df[survive_col] == 1)
    hi_mean = np.average(df.loc[hi_surv_mask, outcome_col], weights=w.loc[hi_surv_mask])
    lo_mean = np.average(df.loc[lo_surv_mask, outcome_col], weights=w.loc[lo_surv_mask])

    rows = []
    for delta in delta_grid:
        # hi_mean = w_always * always_survivor_mean + (1 - w_always) * (always_survivor_mean + delta)
        always_survivor_mean = hi_mean - (1 - w_always) * delta
        rows.append((delta, always_survivor_mean, always_survivor_mean - lo_mean))
    return pd.DataFrame(rows, columns=["delta_protected_minus_always", "implied_always_survivor_mean", "implied_sace"])


bounds = sace_bounds_monotonicity(df)   # pass weight_col="iptw" for observational claims/EHR data
print(f"SACE bounds under monotonicity: [{bounds['sace_lower_bound']:+.2f}, {bounds['sace_upper_bound']:+.2f}]")
print(f"Naive survivor-only difference (biased): {bounds['naive_survivor_only_diff']:+.2f}")

tp = sace_sensitivity_tipping_point(df)
print(tp[tp["delta_protected_minus_always"].isin([-10, 0, 10])])   # spot-check a few sensitivity values

weighted_trim_mean <- function(values, weights, trim_frac, drop = c("low", "high")) {
  # Mean of the weighted distribution after dropping the trim_frac share of total WEIGHT from the low or
  # high tail. Weight-based trimming lets trim_frac be non-integer once rows carry unequal IPTW weights,
  # and is more honest than rounding to the nearest observation even when every weight equals 1.
  drop <- match.arg(drop)
  ord <- order(values)
  values <- values[ord]; weights <- weights[ord]
  cum_weight <- cumsum(weights)
  total_weight <- cum_weight[length(cum_weight)]
  keep <- if (drop == "low") cum_weight > trim_frac * total_weight else cum_weight <= (1 - trim_frac) * total_weight
  if (!any(keep)) return(NA_real_)
  weighted.mean(values[keep], weights[keep])
}

sace_bounds_monotonicity <- function(df, higher_surv_arm = "A", lower_surv_arm = "B", weight_col = NULL) {
  # weight_col: NULL for an RCT (all rows weight 1). For observational data, the name of a column holding
  # pre-fit IPTW weights (see the implementation description above); this function does not fit the
  # exposure model itself, it only consumes already-exchangeable weighted summaries.
  if (!all(df$survived %in% c(0, 1))) stop("survived must be binary (0/1).")
  w <- if (is.null(weight_col)) rep(1, nrow(df)) else df[[weight_col]]

  hi <- df$arm == higher_surv_arm
  lo <- df$arm == lower_surv_arm
  s_hi <- weighted.mean(df$survived[hi], w[hi])
  s_lo <- weighted.mean(df$survived[lo], w[lo])
  if (s_hi < s_lo) stop("Monotonicity direction violated: higher_surv_arm must have survival >= lower_surv_arm.")
  if (s_hi <= 0) stop("No survivors observed in higher_surv_arm; cannot estimate SACE bounds.")

  pi_always <- s_lo
  pi_protected <- s_hi - s_lo
  frac_protected <- pi_protected / s_hi

  hi_surv <- hi & df$survived == 1
  hi_outcome <- df$outcome[hi_surv]; hi_weight <- w[hi_surv]
  if (length(hi_outcome) == 0) stop("No observed survivors with a non-missing outcome in higher_surv_arm.")

  upper_mean_hi <- weighted_trim_mean(hi_outcome, hi_weight, frac_protected, "low")   # protected assumed lowest-scoring -> trim bottom
  lower_mean_hi <- weighted_trim_mean(hi_outcome, hi_weight, frac_protected, "high")  # protected assumed highest-scoring -> trim top

  lo_surv <- lo & df$survived == 1
  lo_mean <- weighted.mean(df$outcome[lo_surv], w[lo_surv])
  naive_mean <- weighted.mean(hi_outcome, hi_weight)

  list(survival_higher_arm = s_hi, survival_lower_arm = s_lo,
       pi_always_survivor = pi_always, pi_protected = pi_protected, frac_protected = frac_protected,
       sace_lower_bound = lower_mean_hi - lo_mean,
       sace_upper_bound = upper_mean_hi - lo_mean,
       naive_survivor_only_diff = naive_mean - lo_mean)
}

sace_sensitivity_tipping_point <- function(df, higher_surv_arm = "A", lower_surv_arm = "B",
                                            weight_col = NULL, delta_grid = seq(-20, 20, by = 1)) {
  w <- if (is.null(weight_col)) rep(1, nrow(df)) else df[[weight_col]]
  hi <- df$arm == higher_surv_arm
  lo <- df$arm == lower_surv_arm
  s_hi <- weighted.mean(df$survived[hi], w[hi])
  s_lo <- weighted.mean(df$survived[lo], w[lo])
  if (s_hi <= 0) stop("No survivors observed in higher_surv_arm; cannot run the sensitivity analysis.")
  w_always <- s_lo / s_hi

  hi_surv <- hi & df$survived == 1
  lo_surv <- lo & df$survived == 1
  hi_mean <- weighted.mean(df$outcome[hi_surv], w[hi_surv])
  lo_mean <- weighted.mean(df$outcome[lo_surv], w[lo_surv])
  always_survivor_mean <- hi_mean - (1 - w_always) * delta_grid
  data.frame(delta_protected_minus_always = delta_grid,
             implied_always_survivor_mean = always_survivor_mean,
             implied_sace = always_survivor_mean - lo_mean)
}

b <- sace_bounds_monotonicity(df)   # pass weight_col = "iptw" for observational claims/EHR data
cat(sprintf("SACE bounds under monotonicity: [%+.2f, %+.2f]\n", b$sace_lower_bound, b$sace_upper_bound))
cat(sprintf("Naive survivor-only difference (biased): %+.2f\n", b$naive_survivor_only_diff))

tp <- sace_sensitivity_tipping_point(df)
print(subset(tp, delta_protected_minus_always %in% c(-10, 0, 10)))

/* (1) Weighted arm-level survival rates and identified stratum proportions
       (monotonicity: arm A = higher survival). */
proc means data=work.cohort noprint;
  class arm;
  var survived;
  weight weight;
  output out=surv(where=(_TYPE_=1)) mean=surv_rate;
run;
data _null_;
  set surv; by arm;
  if arm = 'A' then call symputx('s_hi', surv_rate);
  if arm = 'B' then call symputx('s_lo', surv_rate);
run;
%macro check_monotonicity;
  %if %sysevalf(&s_hi < &s_lo) %then %do;
    %put ERROR: Monotonicity direction violated -- arm A must have survival >= arm B.;
    %abort cancel;
  %end;
%mend;
%check_monotonicity;
%let pi_protected = %sysevalf(&s_hi - &s_lo);
%let frac_protected = %sysevalf(&pi_protected / &s_hi);

/* (2) Zhang-Rubin trimming bounds: trim the protected-fraction share of cumulative WEIGHT (not raw row
       count) from each tail of the higher-survival arm's observed outcomes. */
proc sort data=work.cohort(where=(arm='A' and survived=1)) out=hi_surv;
  by outcome;
run;
data hi_surv_cum;
  set hi_surv; retain cum_weight 0; cum_weight + weight;
run;
proc sql noprint;
  select sum(weight) into :total_weight from hi_surv;
quit;

proc means data=hi_surv_cum noprint;
  where cum_weight > &frac_protected * &total_weight;         /* drop bottom-weight tail -> upper bound */
  var outcome; weight weight; output out=upper_stat mean=upper_mean_hi;
run;
proc means data=hi_surv_cum noprint;
  where cum_weight <= (1 - &frac_protected) * &total_weight;  /* drop top-weight tail -> lower bound */
  var outcome; weight weight; output out=lower_stat mean=lower_mean_hi;
run;
proc means data=hi_surv noprint;                                /* FULL (untrimmed) arm-A survivor mean --
                                                                    needed for the naive diff and the
                                                                    sensitivity block below, computed
                                                                    separately from the trimmed bounds. */
  var outcome; weight weight; output out=hi_full_stat mean=hi_mean_full;
run;
proc means data=work.cohort noprint;
  where arm = 'B' and survived = 1;
  var outcome; weight weight; output out=lo_stat mean=lo_mean;
run;

data sace_bounds;
  merge upper_stat(keep=upper_mean_hi) lower_stat(keep=lower_mean_hi)
        hi_full_stat(keep=hi_mean_full) lo_stat(keep=lo_mean);
  sace_upper_bound = upper_mean_hi - lo_mean;
  sace_lower_bound = lower_mean_hi - lo_mean;
  naive_survivor_only_diff = hi_mean_full - lo_mean;
run;
proc print data=sace_bounds; run;

/* (3) Sensitivity-parameter tipping-point curve: implied always-survivor mean and SACE across a delta
       grid, using the FULL (untrimmed) weighted arm-A survivor mean, not a trimmed bound. */
data tipping_point;
  set hi_full_stat(keep=hi_mean_full);
  if _n_ = 1 then set lo_stat(keep=lo_mean);
  w_always = &s_lo / &s_hi;
  do delta = -20 to 20 by 1;
    implied_always_survivor_mean = hi_mean_full - (1 - w_always) * delta;
    implied_sace = implied_always_survivor_mean - lo_mean;
    output;
  end;
  keep delta implied_always_survivor_mean implied_sace;
run;
proc print data=tipping_point; where delta in (-10, 0, 10); run;

Citations

1. [1]Frangakis CE, Rubin DB. Principal stratification in causal inference. Biometrics. 2002;58(1):21-29.
2. [2]Hayden D, Pauler DK, Schoenfeld D. An estimator for treatment comparisons among survivors in randomized trials. Biometrics. 2005;61(1):305-310.
3. [3]Rubin DB. Causal inference through potential outcomes and principal stratification: Application to studies with "censoring" due to death. Statistical Science. 2006;21(3):299-309.
4. [4]Chiba Y, VanderWeele TJ. A simple method for principal strata effects when the outcome has been truncated due to death. American Journal of Epidemiology. 2011;173(7):745-751.
5. [5]Bornkamp B, Rufibach K, Lin J, et al. Principal stratum strategy: Potential role in drug development. Pharmaceutical Statistics. 2021;20(4):737-751.

1. [6]Zhang JL, Rubin DB. Estimation of causal effects via principal stratification when some outcomes are truncated by "death". Journal of Educational and Behavioral Statistics. 2003;28(4):353-368.
2. [7]Long DM, Hudgens MG. Sharpening bounds on principal effects with covariates. Biometrics. 2013;69(4):812-819.
3. [8]Hirano K, Imbens GW, Rubin DB, Zhou XH. Assessing the effect of an influenza vaccine in an encouragement design. Biostatistics. 2000;1(1):69-88.
4. [9]Tong G, Li F, Chen X, Hirani SP, Newman SP, Wang W, Harhay MO. A Bayesian approach for estimating the survivor average causal effect when outcomes are truncated by death in cluster-randomized trials. American Journal of Epidemiology. 2023;192(6):1006-1015.

1. [10]Egleston BL, Scharfstein DO, MacKenzie E. On estimation of the survivor average causal effect in observational studies when important confounders are missing due to death. Biometrics. 2009;65(2):497-504.
2. [11]International Council for Harmonisation (ICH). ICH E9(R1): Addendum on Estimands and Sensitivity Analysis in Clinical Trials to the Guideline on Statistical Principles for Clinical Trials. Step 5 version, 2019.