Paper deep dive
Observing the unobserved confounding through its effects: toward randomized trial-like estimates from real-world survival data
Vasiliki Stoumpou, Dimitris Bertsimas, Samuel Singer, Georgios Antonios Margonis
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 95%
Last extracted: 4/15/2026, 1:33:32 AM
Summary
The paper introduces a three-step framework to mitigate unobserved confounding in observational survival data by inferring a latent prognostic factor (U) from restricted mean survival time (RMST) discrepancies among similar patients, balancing this factor alongside observed covariates, and applying multivariable survival analysis to estimate hazard ratios.
Entities (5)
Relation Signals (3)
Restricted Mean Survival Time (RMST) → usedtoinfer → Latent Prognostic Factor (U)
confidence 98% · We infer a latent prognostic factor (U) from restricted mean survival time (RMST) discrepancies.
Latent Prognostic Factor (U) → addresses → Unobserved Confounding
confidence 95% · The framework aims to address unobserved confounding in observational survival data by inferring a latent prognostic factor.
Prognostic Matching → balances → Latent Prognostic Factor (U)
confidence 95% · We balance U with observed baseline covariates using prognostic matching.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Background: Randomized controlled trials (RCTs) are costly, time-consuming, and often infeasible, while treatment-effect estimation from observational data is limited by unobserved confounding. Methods: We developed a three-step framework to address unobserved confounding in observational survival data. First, we infer a latent prognostic factor (U) from restricted mean survival time (RMST) discrepancies between patients with similar observed factors, the same treatment, and divergent outcomes, leveraging the idea that the aggregate effect of unmeasured factors can be inferred even if individual factors cannot. Second, we balance U with observed baseline covariates using prognostic matching, entropy balancing, or inverse probability of treatment weighting. Third, we apply multivariable survival analysis to estimate hazard ratios (HRs). We evaluated the framework in three observational cohorts with RCT benchmarks, two RCT cohorts, and six multicenter observational cohorts. Results: In three observational cohorts (nine comparisons), balancing U improved agreement with trial HRs in all cases; in the strongest settings, it reduced absolute log-HR error by approximately ten-fold versus using observed covariates alone (mean reduction 0.344; p=0.001). In two RCT cohorts, U was balanced across arms (most SMDs <0.1) and adjustment had minimal impact on log-HRs (mean absolute change 0.08). Across six multicenter cohorts, balancing U within centers reduced cross-center dispersion in chemotherapy log-HR estimates (mean reduction 0.147; p=0.016); when populations were directly balanced across centers to account for case-mix differences, cross-center survival differences were narrowed in 75%-100% of comparisons. Conclusions: Inferring and balancing a latent prognostic signal may reduce unobserved confounding and improve treatment-effect estimation from real-world data.
Tags
Links
- Source: https://arxiv.org/abs/2604.12137v1
- Canonical: https://arxiv.org/abs/2604.12137v1
Trouble viewing inline? Open PDF directly →
Full Text
93,737 characters extracted from source content.
Expand or collapse full text
Observing the unobserved confounding through its effects: toward randomized trial-like estimates from real-world survival data Vasiliki Stoumpou 1 , Dimitris Bertsimas 1,2 , Samuel Singer 3 , Georgios Antonios Margonis 2,3,4 1 Operations Research Center, Massachusetts Institute of Technology, Cambridge, MA, USA 2 Sloan School of Management, Massachusetts Institute of Technology, Cambridge, MA, USA 3 Department of Surgery, Memorial Sloan Kettering Cancer Center, New York, NY, USA 4 Charité-Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität zu Berlin, Berlin, Germany margonig@mskcc.org Abstract Background: Randomized controlled trials (RCTs) are costly, time-consuming and often infeasible, whereas treatment-effect estimation from observational data is limited by unobserved confounding. Methods: We developed a three-step framework to address unobserved confounding in observational survival data. First, we infer a latent prognostic factor ( ̃ U) from restricted mean survival time (RMST) discrepancies between patients with similar observed prognostic factors, the same treatment assignment, and divergent outcomes, based on the premise that the aggregate effect of unmeasured prognostic factors can be inferred even when individual factors cannot be observed. Second, we balance ̃ U together with observed baseline prognostic factors between treatment groups using prognostic matching, entropy balancing, or inverse probability of treatment weighting. Third, we apply multivariable survival analysis to estimate hazard ratios (HRs). We evaluated the framework in three settings: three observational cohorts with RCT benchmark HRs, two RCT cohorts, and six multicenter observational cohorts. Results: In three observational cohorts (nine cohort-method comparisons), balancing ̃ U improved agreement with benchmark trial HRs in all comparisons; in the strongest settings, it reduced absolute log-HR error by approximately ten-fold relative to adjustment using observed covariates alone (mean reduction, 0.344; p=0.001). In the two RCT cohorts, ̃ U was balanced across treatment arms (most SMDs <0.1) and adjusting for it had little effect on estimated treatment log-HRs (mean absolute change ≃ 0.08). Across six multicenter cohorts, balancing ̃ U reduced cross-center dispersion in chemotherapy log-HR estimates for overall- and recurrence-free survival (mean reduction 0.147; p = 0.016) when balancing was performed within centers to address residual confounding, and narrowed cross-center survival differences in 75% to 100% of informative comparisons when populations were directly balanced across centers to account for differences in patient case mix. Conclusions: Inferring and balancing a latent prognostic signal may reduce unobserved confounding in observational survival analyses and support more reliable treatment-effect estimation from real-world data. 1. Introduction Observational real-world data are increasingly used to eval- uate treatment effects, particularly when randomized controlled trials (RCTs) are infeasible, underpowered, outdated, or not fully representative of contemporary practice. However, such analyses remain vulnerable to bias, including immortal time bias, selection bias, and confounding by indication. Although recent methodological advances have improved our ability to address some of these problems, and standard regression, match- ing, and weighting approaches can mitigate confounding from measured covariates, unobserved confounding remains funda- mentally difficult to resolve. By definition, it is not possible to adjust for variables that are not measured; one cannot regress or match on what is not observed. As a result, even after careful adjustment for measured co- variates, treatment groups may still differ in important unmea- sured aspects of prognosis, potentially leading observational studies to produce biased or even misleading treatment-effect estimates. This challenge is especially relevant in surgical on- cology, where decisions regarding surgery and perioperative treatment often depend on nuanced aspects of disease sever- ity, patient fitness, and clinician judgment that are not fully captured in routine datasets. In practice, this problem is of- ten addressed not by solving it, but by assuming that unob- served confounding is negligible for the clinical question at hand. Yet that assumption is frequently implausible, as illus- trated by settings in which treatment-effect estimates from ob- servational data were later contradicted by randomized evidence. 1, 2 These discrepancies highlight the need for methods that can ac- count, at least indirectly, for latent prognostic differences that remain outside the recorded data. arXiv:2604.12137v1 [stat.AP] 13 Apr 2026 Research in context Evidence before this study Unobserved confounding arising from unknown, unmea- sured, or unrecorded prognostic factors remains a fundamen- tal challenge in clinical outcomes research. When such fac- tors are imbalanced across treatment groups, treatment effect estimates can be biased even after adjustment for observed covariates. Existing approaches generally do not directly ad- dress this problem; instead, they rely on the assumption that unobserved confounding is limited or negligible in the clini- cal setting under study. As a result, the validity of many ob- servational analyses depends heavily on design restrictions and unverifiable assumptions. Added value of this study Existing methods for treatment- effect estimation in real-world data, including inverse prob- ability of treatment weighting, entropy balancing, and prog- nostic matching, adjust only for observed covariates and therefore remain fundamentally limited by unobserved con- founding.The added value of this study is that it ex- tends these standard frameworks by incorporating a patient- level latent quantity, ̃ U, that indirectly captures the aggre- gate effects of unobserved prognostic factors. Thus, rather than only balancing what is measured, the proposed ap- proach seeks to also balance part of what is unmeasured. Empirically, this added value was supported by improved agreement with benchmark randomized trials, by preserva- tion of balance of ̃ U across randomized trial arms, and by findings in multicenter colorectal liver metastasis cohorts showing that ̃ U captured both hidden factors contributing to treatment-selection bias within cohorts and hidden dif- ferences in case mix across centers. These are sources of heterogeneity that conventional observed-covariate methods cannot directly address. Implications of all the available evidence This framework provides a practical strategy to strengthen causal inference from real-world data in the presence of likely unobserved confounding. By augmenting standard adjustment methods with a patient-level estimate of latent prognostic burden, it may help observational analyses better approximate trial- like treatment-effect estimates. In doing so, it could be par- ticularly useful for revisiting unresolved clinical debates in settings where randomized evidence is unavailable, limited, outdated, or difficult to generalize to contemporary practice. In this study, we propose a framework to address this prob- lem by inferring and balancing a latent prognostic signal that reflects the aggregate effect of otherwise unobserved prognos- tic factors. Rather than attempting to recover each missing factor individually, the framework aims to capture their com- bined prognostic influence before estimating treatment effects as hazard ratios (HRs). We evaluate the framework in three complementary settings. First, in three observational surgical oncology cohorts, we test whether it moves treatment-effect es- timates closer to corresponding RCT benchmarks. Second, in two oncology RCT datasets, where confounding should be ab- sent by design, we test whether it leaves the trial-reported HRs largely unchanged. Third, in multicenter colorectal liver metas- tasis (CRLM) data, we test whether the framework reduces two distinct forms of between-center heterogeneity: hidden bias in who receives adjuvant chemotherapy within a center, reflected in variation in chemotherapy HR estimates across centers, and hidden differences in patient case mix across center popula- tions, reflected in residual cross-center survival differences after balancing recorded covariates. 2. Methods We first provide an overview of the proposed methodology and then describe each component in detail. Our approach aims to transform an observational cohort to a population that closely mimics the conditions of an RCT, namely balance in both ob- served and unobserved covariates. The ultimate goal is to en- able reliable treatment effect estimation within a survival anal- ysis framework, which can approximate the results we would expect if the data had come from an actual RCT evaluating the same therapy. The first step addresses unobserved prognostic factors and constitutes the main contribution of our methodology. We esti- mate an unobserved feature ̃ U, which serves as a latent prognos- tic factor capturing outcome-related differences not explained by the observed covariates. In other words, this variable is con- structed to capture variation in the outcome that cannot be ade- quately justified by the observed features alone. The second step incorporates this latent factor together with the observed covariates in the balancing procedure, using stan- dard balancing methods such as prognostic matching followed by sample weighting, 3 entropy balancing, 4 or inverse propen- sity weighting. 5 By balancing both observed covariates and the inferred latent prognostic factor across treatment arms, this step aims to address both observed and unobserved confounding. The third step involves multivariable analysis, where we fit a Cox proportional hazards model, which serves as a final cor- rection for any residual observed confounding and remaining imbalances. The coefficient of the treatment variable in the re- sulting model is then utilized to estimate the treatment Hazard Ratio, which is used to evaluate the treatment effect. The full pipeline is depicted in Figure 1. The following Sections elaborate on each of these three steps. 2.1. Step 1: Constructing the latent factor U This step introduces our core methodological contribution: estimating a hidden prognostic factor that can later be used to correct for unobserved confounding. Our goal is to define a scalar quantity ̃ U i that captures the portion of patient i’s outcome that cannot be explained by ob- served covariates. One way to approach the construction of ̃ U i is to utilize patient i’s nearest neighbors, namely patients with similar observed characteristics, treatment, and clinical presen- tations. Such patients are expected to have similar outcomes to patient i; however, in practice, individuals with nearly identical observed profiles may experience very different survival trajec- tories (e.g., recurrence vs. non-recurrence), suggesting the in- fluence of unobserved patient-level factors. 2 Step 1 Estimating latent factor 푈 " capturing unobserved prognostic influences Step 2 Balancing 푋,푈 " between treated and untreated cohorts ≈ w Step 3 Multivariable analysis + w Figure 1: Full pipeline of the proposed framework. Intuitively, we quantify how “unexpected” a patient’s out- come is relative to similar patients within the same treatment arm. When two such patients exhibit different outcomes, the discrepancy primarily reflects latent factors; as these patients are similar in their observed features and share the same treat- ment assignment, the difference in outcomes is mostly driven by unobserved influences such as biological differences, undoc- umented comorbidities, or other unmeasured sources of hetero- geneity. Quantifying this discrepancy between the neighbors therefore provides a way to capture the influence of these un- measured factors on the outcome: a small discrepancy suggests that unobserved factors have limited impact, whereas a larger discrepancy indicates a stronger latent influence. Importantly, we do not attempt to recover any specific latent variable, as this would be infeasible for factors that are not directly observed or even known. Instead, ̃ U should be interpreted as a summary measure of the combined influence of unobserved factors af- fecting the outcome, which is expected to be balanced in ran- domized controlled trials. There are two important nuances in the construction of ̃ U i , without which differences in outcomes between neighboring patients cannot be safely attributed to unobserved factors. The first is that focusing on neighbors within the same treat- ment arm is crucial. Under this restriction, we can reasonably assume that the interaction between treatment and the observed covariates is approximately the same for nearby patients in the same treatment arm, and that the baseline effect of observed covariates is similar for nearby patients in the untreated arm. Consequently, differences in outcomes between such neighbors are not attributed to treatment-related heterogeneity, but to hid- den, underlying factors that affect their outcomes in different ways. The second is that ̃ U i must satisfy a directional sign con- vention. By forcing the comparison to occur against neighbors with the opposite survival trajectory, we ensure that the sign of the residual reflects latent influence with respect to the ob- served outcome. Under this convention, negative values cor- respond to latent adverse factors associated with worse-than- expected survival, whereas positive values correspond to latent favorable factors associated with better-than-expected survival, conditional on the observed covariates and treatment status. For example, among patients who do not receive chemotherapy de- spite having similar observed disease characteristics and appar- ent indication for treatment, negative values may reflect unmea- sured frailty or comorbidity burden that both discouraged treat- ment and contributed to earlier death, whereas positive values may reflect preserved physiologic reserve or other favorable un- measured characteristics in patients who remained untreated for non-clinical reasons, such as treatment refusal. The farther ̃ U i lies from zero, the stronger the inferred influence of these unob- served factors. Without this directional selection of neighbor- hoods, a patient with an unexpectedly early event could have a value of ̃ U i similar to that of a patient with unexpectedly long survival, as contributions from neighbors with mixed outcome trajectories could cancel out. These two elements are essential: selecting neighbors purely based on proximity in the covariate space would remove the directional interpretation of the residual, while selecting neigh- bors from a different treatment arm would reintroduce treatment effects into the comparison. In either case, the interpretation of the residual as reflecting latent factors would be undermined. 2.1.1. Directional neighbor sets (within arm) We now formalize the neighborhood-based construction in- troduced above. Let X i denote the observed covariates of patient i, T i their observed follow-up time, and E i ∈ 0, 1 the event indicator (1 = event observed; 0 = censored). Let V i denote un- observed factors (e.g., frailty, latent tumor aggressiveness, un- documented comorbidities) that may influence both treatment assignment and outcomes. 3 In addition to the observed covariates, we use a prognostic score s i summarizing the baseline risk of patient i. The score s i is obtained either from an externally validated prognostic model or from a predictive model trained on the observed covariates, using cross-fitting on the untreated cohort to avoid overfitting. Throughout the methodology, s i is treated as an additional base- line feature: it may be used when constructing neighborhoods for the latent factor U, and later as a balancing feature when aligning treated and untreated cohorts. For each patient i and a fixed neighborhood size k, we define a directional neighbor set N i,k restricted to the same treatment arm (t j = t i ) and selected among patients whose observed sur- vival trajectory is opposite in the relevant sense. • If patient i experienced the event (E i = 1) at time T i , candidate neighbors are those who survived longer (i.e., had a later event time or were censored after T i ). Oper- ationally, we select neighbors among patients satisfying T j > T i . • If patient i was censored (E i = 0) at time T i , candidate neighbors are those who experienced the event earlier, i.e., patients satisfying E j = 1 and T j < T i . Among these candidates, we select the k nearest neighbors in the X-space according to a distance metric d(·,·) (Euclidean distance after standardizing the covariates). Formally, for event patients (E i = 1), N i,k =arg min S⊆ j: t j =t i , T j >T i |S|=k X j∈S d(X i , X j ),if E i = 1,(1) and analogously, for censored patients (E i = 0), N i,k =arg min S⊆ j: t j =t i , E j =1, T j <T i |S|=k X j∈S d(X i , X j ),if E i = 0.(2) This construction ensures that comparisons are both (i) within treatment and (i) anchored on patients whose observed survival trajectories represent the opposite behavior relative to patient i, conditional on observed follow-up information. 2.1.2. Residual-type latent factor We then define the raw latent factor as the difference be- tween patient i’s outcome signal and the average outcome sig- nal of its directional neighbors: U i = Y i − 1 |N i,k | X j∈N i,k Y j .(3) This definition admits a natural structural interpretation. Let the outcome signal of interest be written as Y i = f t (X i , V i ) + ε i ,(4) where t indexes the treatment arm (t = 0 untreated, t = 1 treated), f t (·) captures the joint contribution of observed and unobserved covariates under treatment t, and ε i represents id- iosyncratic noise. This functional form is intentionally general and does not impose structural assumptions on the relationship between observed covariates, latent factors, and outcomes. Substituting (4) into (3) yields U i = f t (X i , V i )− 1 |N i,k | X j∈N i,k f t (X j , V j ) + ε i − ̄ε N i,k = f t (X i , V i )− 1 |N i,k | X j∈N i,k f t (X i , V j ) + 1 |N i,k | X j∈N i,k f t (X i , V j )− f t (X j , V j ) + ε i − ̄ε N i,k .(5) Because neighbors are selected to be close in the covariate space (X j ≈ X i ), the second bracketed term, which represents differences due to observed features, is expected to be small. Under this local approximation, U i can therefore be inter- preted as a contrast in the unobserved component: U i ≈ f t (X i , V i )− E[ f t (X i , V )|N i,k ],(6) where the expectation is taken over the directional neighbor set. Thus, U i captures how much the realized prognosis of pa- tient i deviates from that of patients with similar observed char- acteristics receiving the same treatment, which is attributed to unobserved factors. Using a neighborhood rather than a sin- gle closest neighbor is important, as averaging across neigh- bors smooths variation and reduces the influence of extreme outcomes that could otherwise distort the estimation of U i . 2.1.3. Choice of Y i in survival settings: pseudo-RMST In survival settings, the quantity Y i in Equation (4) is not directly observed as a scalar due to censoring. A naive option would be to use the event indicator E i , but this discards time information and is highly sensitive to censoring. Instead, we use the restricted mean survival time (RMST) as the outcome signal because it summarizes the survival curve up to a clinically meaningful horizon τ and remains robust to cen- soring when estimated appropriately. We compute individual- level pseudo-observations for RMST using a jackknife construc- tion. 6 Let b μ denote the RMST estimated from the full sample of size n, and b μ (−i) the estimate obtained after leaving patient i out. The jackknife pseudo-observation is Y i = n b μ − (n− 1) b μ (−i) .(7) This construction transforms censored survival data into a scalar outcome contribution reflecting how patient i perturbs the estimated survival curve, providing a smoother and more infor- mative signal than the raw follow-up time and event indicator alone. With the directional neighbor construction above, the sign of U i has a natural interpretation: 4 • For event patients (E i = 1), neighbors are selected among those who survive longer; hence Y i is typically smaller than the neighbor average, yielding U i < 0, which corre- sponds to higher latent risk. • For censored patients, neighbors are selected among those who experience an earlier event; hence Y i tends to exceed the neighbor average, yielding U i > 0, corresponding to a protective latent profile. Thus the sign of U reflects latent risk direction, while its magnitude captures how atypical the outcome signal is relative to nearby patients. 2.1.4. Percentile normalization across arms To make U comparable across treatment arms while pre- serving its directional interpretation, we apply a signed normal- ization within each arm, producing ̃ U ∈ [−1, 1]. Because the raw values of U may be skewed and contain extreme observa- tions, the normalization is performed separately for positive and negative values within each arm using winsorization that limits the influence of outliers while preserving the sign interpreta- tion. The key comparability assumption is that the directional meaning of ̃ U is consistent across treatment arms: more nega- tive values correspond to stronger latent risk, and more positive values correspond to stronger latent protection relative to pa- tients with similar observed characteristics. The signed within- arm normalization places | ̃ U| on a comparable, unitless scale across treatment arms, allowing ̃ U to be interpreted as a relative deviation from expected prognosis within arm. Under random treatment assignment, latent prognostic factors are expected to be similarly distributed across arms; therefore, after this nor- malization, the distribution of ̃ U should also be approximately balanced across treatment arms in RCT or RCT-like cohorts. The resulting normalized quantity ̃ U is subsequently used as an additional covariate in the balancing step between treatment arms, summarizing the latent prognostic component of each pa- tient. 2.2. Balancing the Treatment Arms After constructing ̃ U, the next step is to balance the treated and untreated cohorts. Specifically, we construct a weighted pseudo-cohort in which the distribution of baseline covariates X and the latent feature ̃ U is comparable between treated and un- treated patients; this can enable reliable estimation of the treat- ment effect. To assess the framework’s robustness across dif- ferent balancing procedures, we consider three complementary strategies: prognostic matching followed by sample weighting, entropy balancing, and inverse propensity weighting. 2.2.1. Prognostic matching followed by weighting Prognostic matching 3 reduces observed confounding by (i) stratifying patients into buckets based on a predicted baseline risk score (prognostic score), and (i) performing one-to-one matching within each bucket using standardized baseline co- variates. This design encourages similarity both in baseline co- variate profiles and in baseline prognosis. Although ̃ U could in principle be treated as an additional baseline covariate and included in the stratification, we avoid doing so during the prognostic-matching step because ̃ U is con- structed to be strongly associated with the outcome, and would distort the interpretation of the prognostic buckets. Further- more, prognostic matching primarily enforces similarity in the prognostic score rather than strict similarity in the individual covariates X. Treating ̃ U as an additional feature when calcu- lating distances for the matching step would therefore not guar- antee adequate balance in the effect of the unobserved factors. Instead, we extend the weighting approach of prognostic matching. After performing prognostic matching, we compute a scalar weight w that is applied to a pre-specified subgroup of the treated cohort in order to align the mean of ̃ U between treated and untreated patients. More specifically, let T = i : t i = 1 and U = i : t i = 0 denote the treated and untreated sets, respectively. We denote by G = i : E i = 0 the set of patients without an observed event and by B = i : E i = 1 the set of patients with an observed event. In our implementa- tion, we apply this scalar weight to the subgroupT∩G (treated patients without an observed event), which yields a simple and interpretable reweighted treated cohort. We note that one could alternatively assign weights to other subgroups. We use a single scalar weight to avoid instability and overfitting that may arise from highly flexible, individual-specific weights. Concretely, we calculate w such that the mean of ̃ U among the untreated patients equals the mean of ̃ U in the reweighted treated cohort: 1 |U| X i∈U ̃ U i = P i∈T∩B ̃ U i + w P i∈T∩G ̃ U i |T ∩B| + w|T ∩G| .(8) Solving (8) for w yields a reweighted treated cohort whose mean ̃ U matches that of the untreated cohort. To avoid extreme weights and to ensure positivity, we constrain w∈ [0.5, 20]; in practice, extreme values are rare. 2.2.2. Entropy balancing As a second balancing approach, we use entropy balancing to reweight control patients so that selected moments of certain covariates match those of the treated cohort. Specifically, we compute nonnegative weights for control patients that approx- imately match the first moments of the treated group for the baseline covariates, the prognostic score, and the latent factor ̃ U. In addition, we match the second moments of the latter two quantities to better account for differences in their distributions across treatment arms, as the prognostic score summarizes ob- served features and ̃ U summarizes latent prognostic factors. Let Z i denote the vector of features to be balanced. In our implementation, Z i = X i , s i , ̃ U i , s 2 i , ̃ U 2 i , where X i are standardized baseline covariates and s i is the prog- nostic score. Entropy balancing assigns weights w i to controls (t i = 0) by 5 solving the convex optimization problem min w i i∈U X i∈U w i log w i (9) s.t. 1 |T| X i∈T Z i = 1 |T| X i∈U w i Z i , X i∈U w i =|T|, w i ≥ 0, ∀i∈U. We set weights for treated patients to one (w i = 1 for i ∈ T ) and only reweight the controls. 2.2.3. Inverse propensity weighting (IPTW) The third balancing method is inverse propensity weight- ing, aimed at estimating weights for the untreated cohort so that it resembles the treated cohort. We fit a logistic regression propensity model on the full cohort to estimate p i = P(t i = 1 | X i , s i , ̃ U i ). We use logistic regression (rather than more flexible models) to reduce overfitting and improve stability. We then apply the weights by assigning weight 1 to treated patients and reweighting controls according to w i = 1,t i = 1, p i 1− p i , t i = 0. (10) This weighting scheme upweights control patients who look more likely to have been treated, making the weighted control cohort comparable to the treated cohort, similarly to the entropy balancing method. 2.3. Multivariable analysis After balancing, we estimate the treatment effect using a multivariable Cox proportional hazards model fit on the weighted samples, with the weights determined by the chosen balancing strategy. The Cox model includes the treatment indicator and all baseline observed covariates (and prognostic score when ap- plicable). We do not include ̃ U as a covariate in this final re- gression. The influence of ̃ U is already incorporated through the balancing weights; including ̃ U directly in the Cox model would introduce instability due to the strong association be- tween ̃ U and the outcome. The coefficient of the treatment indi- cator from the weighted Cox model is reported as our adjusted estimate of treatment effect. Taken together, these three steps produce balanced cohorts with respect to both observed covariates and the latent prognos- tic feature ̃ U, enabling more reliable estimation of treatment effects in observational cohorts. 2.4. Statistical analysis of validation experiments To evaluate our framework and its robustness across dif- ferent settings, we experimented with multiple hyperparame- ter configurations; all tuning ranges were defined a priori, and for each balancing method the baseline and ̃ U-augmented ver- sions were run over identical hyperparameter grids. We then summarized performance across the full set of pre-specified pa- rameter combinations (rather than selecting the best-performing setting), which reduces the risk that apparent gains of the aug- mented approach are attributable to post hoc hyperparameter selection. For descriptive visualization, we report average HRs and standard errors across parameter combinations (Tables 2–6). Formal inference was based on paired comparisons between the corresponding baseline and ̃ U-augmented analyses within matched settings (same dataset, same balancing method, and same hyperparameter configuration). For Validation 1 (observational cohorts with benchmark RCT HRs), the performance metric was absolute distance to the bench- mark on the log-HR scale. For each matched run c, we com- puted d base c = log HR base c − log HR RCT , d aug c = log HR aug c − log HR RCT , and defined the paired improvement ∆ c = d base c − d aug c , so that ∆ c > 0 indicates that adding ̃ U moved the estimate closer to the benchmark RCT HR. To avoid treating hyperparameter combinations as independent observations, we first summarized ∆ c within each dataset-by-balancing-method cell (mean across parameter combinations), then tested whether these cell-level summaries were systematically positive using an exact sign test (primary analysis). Ties (cell summaries exactly equal to zero) were excluded from the sign-test denominator. A Wilcoxon signed-rank test on the same cell-level summaries was used as a sensitivity analysis. As an effect-size summary, we report the median and mean cell-level improvement in benchmark dis- tance on the log-HR scale (Table 2). For Validation 2 (RCT cohorts), the objective was to ver- ify that adding ̃ U did not materially distort treatment-effect es- timates in settings where confounding should be absent by de- sign. We therefore used equivalence testing. The augmentation- induced shift was defined as ∆ shift = log HR aug − log HR RCT , and evaluated using two one-sided tests (TOST) against a pre- specified equivalence margin δ, testing whether−δ < ∆ shift < δ. The equivalence margin was specified a priori as a clinically negligible relative change in HR (e.g. δ = log(1.10)). Equiva- lence was concluded only if both one-sided tests were signif- icant. We additionally report the observed log-HR shifts as effect-size summaries (Table 4). For Validation 3 (cross-center consistency in CRLM), we consider two complementary analyses corresponding to distinct estimands: (i) the consistency of treatment-effect estimates across centers, and (i) differences in survival outcomes driven by vari- ation in patient case mix. (i) Treatment-effect consistency across centers.. Between-center heterogeneity in treatment-effect estimates was quantified using 6 the mean absolute pairwise difference of center-specific treat- ment effects on the log-HR scale. For each matched setting c (same outcome, same balancing method, and same hyperpa- rameter configuration), we computed D base c and D aug c , where D c = 1 C 2 X i< j log(HR i,c )− log(HR j,c ) , and defined ∆D c = D base c − D aug c , so that positive values indicate reduced between-center variabil- ity after incorporating ̃ U. This analysis targets within-center confounding: if latent confounders are shared across centers and captured by ̃ U, adjustment should yield more comparable HR estimates. As in Validation 1, we summarized ∆D c within each outcome- by-balancing-method cell (mean across parameter combinations) and tested whether the cell-level summaries were systemati- cally positive using an exact sign test (primary analysis), with a Wilcoxon signed-rank test as a sensitivity analysis. Ties were excluded from the sign-test denominator. As an effect-size sum- mary, we report the median and mean cell-level reduction in pairwise dispersion on the log-HR scale (Table 6). (i) Cross-center survival comparisons.. For the CRLM cen- ters, we performed an additional validation experiment to as- sess the role of ̃ U beyond X when comparing adjusted 5-year survival across centers. In contrast to the HR analysis, this setting does not target confounding in treated-versus-untreated comparisons, but rather differences in patient selection (case mix) across centers. Accordingly, instead of balancing treated and untreated arms within each center, we use the same balanc- ing methods to construct pairwise cross-center comparisons by reweighting patients from a source center to match the covari- ate and latent-factor distribution of a target center. This allows survival differences to be interpreted conditional on comparable patient populations across centers. For each pair of centers (A, B), suppose without loss of gen- erality that center A has lower crude 5-year survival than center B. We define D raw AB = S raw B (5)− S raw A (5), D base AB = S base B (5)− S base A (5), D aug AB = S aug B (5)− S aug A (5), where S raw A (5), S base A (5), and S aug A (5) denote the crude, X-adjusted, and (X, ̃ U)-adjusted 5-year survival estimates for center A, re- spectively. We focused on the subset of center pairs for which the or- dering was preserved after adjustment on X (i.e., S X A (5) < S X B (5)) and for which adjustment using only observed covariates widened the survival gap, namely D X AB > D raw AB . This pattern is informative: if balancing on observed charac- teristics makes the patient populations more comparable, yet the survival disadvantage of center A becomes larger, it sug- gests the presence of important unobserved prognostic differ- ences between centers. Since ̃ U is designed to capture such latent prognostic vari- ation, we would expect that adjusting for ̃ U in addition to X reduces this enlarged gap, that is, D aug AB < D base AB . For this validation experiment, we compute the percentage of selected center pairs for which adjustment using (X, ̃ U) re- duces the X-adjusted survival gap. A higher percentage indi- cates that ̃ U captures previously unobserved prognostic differ- ences between centers. We then perform an exact binomial test against the null hypothesis of success probability 0.5. We note that this analysis targets a diagnostic subset of cen- ter pairs for which adjustment on X reveals potential unob- served prognostic differences; for other center pairs, the ad- dition of ̃ U may either increase or decrease the survival gap depending on the relative contributions of observed and latent prognostic factors. All analyses were performed in Python (version 3.8) using standard scientific computing libraries, including NumPy, pan- das, scikit-learn, lifelines, and SciPy. Statistical tests were con- ducted using the scipy.stats module. 3. Results In this Section, we present the results of our experiments evaluating how incorporating the latent factor ̃ U during the bal- ancing step influences downstream estimates, including both treatment effects (hazard ratios) and cross-center survival com- parisons. We conducted three groups of experiments designed to assess performance under complementary settings regarding the types of datasets employed. To estimate ̃ U, several hyperparameters must be specified. We explored multiple design choices and report average results across all configurations. These include: (i) the number of near- est neighbors k used to compute ̃ U i for anchor patient i, (i) whether the prognostic score s i is included in the distance met- ric for neighbor selection, and (i) which prognostic score is utilized when applicable. Each balancing method also involves its own design param- eters. For prognostic matching, we vary the number of risk strata used during stratification. For entropy balancing, we ex- amine whether matching only first moments or both first and second moments of ̃ U i and the prognostic score s i affects per- formance. For IPTW, we vary clipping thresholds to mitigate extreme weights. An extended sensitivity analysis for all hy- perparameters is provided in Appendix Section 1. 3.1. Real-World Observational Data In the first set of validations, we applied our framework to observational oncology cohorts. The objective was to deter- mine whether the derived cohorts produce treatment-effect esti- mates that align with ground-truth HRs reported in RCTs eval- uating the same therapies. 7 For each dataset, we estimated the HR of the treatment within the constructed cohort and compared it with literature-established HRs from prior RCTs. The HR was computed as the exponen- tiated coefficient of the treatment indicator from a Cox propor- tional hazards model. We evaluated our approach on three observational oncology datasets with tumor recurrence as the event of interest. Table 1 summarizes cohort size and treatment/outcome distributions. For each dataset, we used a clinically established feature set X i , a prognostic score s i , and the latent feature U i to define the profile of patient i. We examined multiple prognostic scores: an externally validated score when available, and an internally trained model fitted on untreated patients only, predicting event probability within a fixed time horizon using cross-fitting. The Memorial Sloan Kettering Cancer Center (MSKCC) gastrointestinal stromal tumor (GIST) dataset comprises patients with resected GIST treated with or without adjuvant imatinib; details of this cohort have been reported previously. 3 According to the benchmark RCT, adjuvant imatinib was associated with an HR of 0.35 for recurrence-free survival. 7 Our feature set included the most commonly used prognostic factors in GIST risk models, namely tumor size, tumor site, and mitotic index. 8 The external prognostic model used in our analysis 9 incorpo- rates these same variables. The MSKCC retroperitoneal sarcoma (RPS) dataset includes patients with retroperitoneal or pelvic sarcoma treated with surgery alone or surgery plus perioperative radiotherapy; details of this cohort have been reported previously. 10 According to the per- protocol analysis of the benchmark RCT, pre-operative radio- therapy was associated with an HR of 0.75 for abdominal recurrence- free survival. 2 Our feature set included the most commonly used prognostic factors in RPS local recurrence risk models, namely age, tumor size, grade, histology, and margin status. The external prognostic model used in our analysis similarly incorporates these same variables. 8, 10 The CRLM multi-institutional dataset comprised patients with resected colorectal liver metastases treated with or with- out adjuvant chemotherapy after surgery; these patients were selected to match the inclusion criteria of the most recent RCT evaluating adjuvant chemotherapy in CRLM. Details of this co- hort have been reported previously. 11 According to that bench- mark RCT, adjuvant FOLFOX chemotherapy was associated with an HR of 0.67 for disease-free survival. 12 The external prognostic model used in our analysis incorporates age, sex, primary tumor T stage, lymph node status (N stage), KRAS mutational status, tumor sidedness, CEA level, disease-free in- terval, tumor size, bilobar involvement, and number of metas- tases. 13 To quantify agreement with the benchmark RCT estimates, we measure the absolute log-HR error between the estimated hazard ratio and the RCT reference value. For dataset d, balanc- ing method m, variant v (X-adjustment or (X, ̃ U)-adjustment), and hyperparameter configuration h, the error is defined as Error d,m,v,h = log( c HR d,m,v,h )− log(HR RCT d ) .(11) The plotted value corresponds to the mean error across all hyperparameter configurations: Error d,m,v = 1 H H X h=1 Error d,m,v,h .(12) Figure 2 presents these average log-HR errors across all configurations for simple covariate adjustment as a baseline, compared with prognostic matching, entropy balancing, and inverse propensity weighting. Across all three observational datasets and across all balancing methods, incorporating the la- tent factor ̃ U consistently reduces the estimation error. The raw HR values are presented in Figure 3. The only dataset that exhibits instability is RPS. This is expected, as the RPS cohort is relatively small and contains very few events in the treated group, making reliable HR estimation more diffi- cult. Despite this limitation, all methods incorporating ̃ U im- prove upon adjustment using X alone. Among the evaluated ap- proaches, IPTW generally produces the most stable estimates. Entropy balancing occasionally overshoots the benchmark HR, likely due to extreme weights produced by the balancing con- straints. To quantify whether incorporating ̃ U systematically improves agreement with the benchmark RCT hazard ratios, we evaluated the directional improvement in absolute log-HR error between models using X alone and those using (X, ̃ U). Statistical testing procedures are described in Section 2.4. Across all datasets and balancing methods, incorporating ̃ U reduced the error relative to the RCT benchmark, resulting in consistent directional im- provements (9/9 comparisons). As summarized in Table 2, this improvement was statistically significant under both the sign test and the Wilcoxon signed-rank test. The differences for all pairs of datasets and balancing methods can be found in the Appendix Table 2.1. 3.2. Randomized Controlled Trials In the second set of validations, we applied our methodol- ogy to RCT datasets. Because treatment assignment in RCTs is randomized, both observed and unobserved confounders are expected to be balanced between study arms. Therefore, incor- porating ̃ U in the balancing methods should not materially alter the estimated HR. In contrast to the observational validation, where estimates are compared against benchmark RCT results, the purpose of this validation is to verify that incorporating the latent factor ̃ U does not materially change treatment-effect es- timates when treatment assignment is already randomized. We evaluated two RCT datasets: STRASS and BRT. The STRASS RCT 2 compared preoperative radiotherapy plus surgery versus surgery alone in patients with RPS, with ab- dominal recurrence-free survival as the primary trial endpoint. The reported HR in the intention-to-treat analysis was 1.00, in- dicating no treatment benefit. The prognostic features used in our analysis were age, tumor size, tumor grade and histological subtype. The BRT RCT included patients with soft tissue sarcomas of the extremity or superficial trunk treated with either adjuvant brachytherapy (BRT) or no further therapy after complete re- section, with local recurrence as the primary trial endpoint. 14 8 Table 1: Cohort characteristics of the observational datasets used in the study. DatasetNSurgery + adjuvant therapySurgery aloneEvents (Surgery + adjuvant therapy)Events (Surgery alone) GIST480105 (21.9%)375 (78.1%)34 (32.4%)52 (13.9%) RPS20129 (14.4%)172 (85.6%)3 (10.3%)49 (28.5%) CRLM778140 (18.0%)638 (82.0%)83 (59.3%)445 (69.8%) Table 2: Validation approach 1: directional improvement in closeness to the benchmark RCT hazard ratio after adding ̃ U. StatisticValue Directional improvements (∆ > 0)9 / 9 Sign test p-value (one-sided)0.00195 Wilcoxon signed-rank p-value (one-sided)0.00195 Median cell-level ∆0.2782 Mean cell-level ∆0.3443 The prognostic features used in our analysis were tumor histol- ogy, margin status, tumor grade, tumor size, and tumor depth. Summary statistics are provided in Table 3. As expected, incorporating ̃ U in the balancing step does not materially change the HR estimates relative to those obtained using X alone, as we observe in Figures 4 and 5. In smaller samples, such as BRT, simple covariate adjustment may yield deviations from the reported HR due to sampling variability; however, the inclusion of ̃ U reduces the error relative to the benchmark HR in two of the three balancing methods. This supports the interpretation that ̃ U corrects latent imbalance only when such imbalance exists. The two-sided statistical test described in Section 2.4 was applied using an equivalence margin of log(1.10), correspond- ing to a maximum allowable change of approximately 10% in the hazard ratio after incorporating ̃ U. Specifically, we evalu- ated the shift ∆ shift = log(HR aug) − log(HR RCT ), and assessed whether the 95% confidence interval of ∆ shift lies entirely within the interval [− log(1.10), log(1.10)]. Using this criterion, equivalence was achieved in four of the six dataset–method combinations. In STRASS, one configuration narrowly exceeded the predefined equivalence margin, while in BRT the most de- viating estimate remained close to the corresponding X-only adjustment. All corresponding values are reported in the Ap- pendix Table 2.2. An additional validation examines whether the estimated ̃ U differs between treatment arms in the RCT datasets. Because treatment assignment is randomized, both observed covariates and latent prognostic factors are expected to be balanced be- tween arms. To assess this, we computed standardized mean differences (SMDs) for both the observed covariates and the es- timated ̃ U values across treatment groups (Figure 6). Across both datasets, the majority of SMDs for ̃ U are below or close to the commonly used balance threshold of 0.1 and are compa- rable to, or even smaller than, those observed for the baseline covariates. Together, these results confirm that the latent factor ̃ U be- haves as expected in randomized settings, remaining balanced across treatment arms and not introducing systematic shifts in treatment-effect estimates. 3.3. Cross-Center Consistency Analysis The third set of validations examined six observational CRLM cohorts from six academic centers; details of these cohorts have been reported previously by our group. 15 These cohorts span the USA, Europe, and Asia and are therefore expected to differ in patient selection for surgery, eligibility for adjuvant chemother- apy, and treatment protocols (e.g., chemotherapy regimen). Summary statistics are shown in Table 5. The heterogene- ity in both treatment allocation and outcomes across centers is evident. The first component of the cross-center validation exam- ines whether balancing on (X, ̃ U) between adjuvant treatment arms within each center reduces between-center differences in treatment-effect estimates. If a meaningful component of the observed heterogeneity is driven by differences in the distri- bution of observed and unobserved confounders across cen- ters, then this balancing should reduce inter-center variability in HR estimates. However, even if unmeasured prognostic fac- tors are successfully adjusted for, we would not expect identical HRs for adjuvant chemotherapy across all cohorts, because both treatment exposure (e.g., chemotherapy regimen and treatment duration) and the characteristics of the treated populations may differ between centers. Across all balancing methods and both outcomes (recur- rence and overall survival), incorporating ̃ U reduces the mean absolute pairwise deviation of HR estimates across centers, as shown in Figure 7. For a given balancing method m, variant v, and outcome e, this deviation is computed as D m,v,e = 1 C 2 X i< j log(HR i,m,v,e )− log(HR j,m,v,e ) ,(13) where C denotes the number of centers and HR i,m,v,e is the esti- mated hazard ratio for center i under method m, variant v, and outcome e. In most settings, the reduction in D m,v,e after in- corporating ̃ U is statistically meaningful, indicating that adjust- ing for ̃ U within each center reduces residual confounding and leads to more consistent treatment-effect estimates across cen- ters. The first component of Validation 3 (Table 6) indicates that adjusting for both X and ̃ U reduces the average absolute pair- wise differences in HR estimates across centers in all six evalu- ated cases. This reduction is statistically significant under both 9 MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 Absolute log-HR error GIST MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 1.0 Absolute log-HR error RPS MatchingEntropy Balancing IPTW 0.0 0.1 0.2 0.3 Absolute log-HR error CRLM Standard (X)Augmented (X + U)Multivariable Analysis Figure 2: Comparison of treatment-effect estimation accuracy across observational datasets. Bars represent the mean absolute log-HR error relative to the benchmark RCT hazard ratio, averaged across hyperparameter configurations. Results are shown for matching, entropy balancing, and inverse propensity weighting using either baseline covariates X or the augmented representation (X, ̃ U). The dashed line denotes multivariable Cox adjustment using baseline covariates only. Error bars represent the standard error across configurations. the sign test and the Wilcoxon signed-rank test. The differences in all pairs of datasets and balancing methods are reported in the Appendix Table 2.3. The second component of the cross-center validation, in contrast, examines differences in patient case mix across cen- ters. In this setting, we construct pairwise comparisons by bal- ancing center populations to each other through reweighting. We performed the statistical analysis described in Section 2.4, focusing on center pairs where adjustment using X alone in- creased disagreement in adjusted 5-year survival relative to the crude estimates. Such cases suggest the presence of unobserved prognostic differences in patient populations across centers. We therefore examined whether incorporating ̃ U in the adjustment could reduce these discrepancies. This is precisely what we observe in Table 7. Across both outcomes (overall survival and recurrence), and in both the full cohort and the surgery alone subset, adding ̃ U reduced the sur- vival gap in the majority of center pairs for which adjustment on X alone had increased disagreement. Moreover, across sev- eral adjustment methods, these improvements occur at rates that are statistically significant under the corresponding binomial tests, supporting the interpretation that ̃ U captures the unob- served prognostic component of cross-center differences in pa- tient case mix, which in turn contributes to cross-center survival variation. 4. Discussion In this study, we developed and evaluated a framework to address unobserved confounding in observational survival anal- yses by inferring and balancing a latent prognostic factor, ̃ U. We examined the framework across three complementary vali- dation settings: comparison of observational HR estimates against benchmark randomized trial effects, application to randomized datasets as a negative-control setting, and multicenter colorec- tal liver metastasis analyses evaluating whether balancing ̃ U reduces between-center heterogeneity in treatment effects and patient case mix. Taken together, these findings suggest that ̃ U captures clinically meaningful latent prognostic information not contained in the recorded covariates alone. The first validation strategy, comparison with benchmark randomized trials, addressed the central question of whether balancing ̃ U improves treatment-effect estimation in confounded observational data. Across the observational cohorts studied here, this was generally the case. Importantly, improved agree- ment with the benchmark trial should not be interpreted as im- 10 UnivariableMatchingEntropy Balancing IPTW 0.0 0.5 1.0 1.5 2.0 2.5 Hazard Ratio GIST UnivariableMatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Hazard Ratio RPS UnivariableMatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 Hazard Ratio CRLM UnivariableStandard (X)Augmented (X + U)MultivariableGround Truth Figure 3: Comparison of treatment-effect estimates across observational datasets. Bars represent the mean hazard ratio (HR), averaged across hyperparameter con- figurations. Results are shown for matching, entropy balancing, and inverse propensity weighting using either baseline covariates X or the augmented representation (X, ̃ U). The dotted line denotes multivariable Cox adjustment using baseline covariates only, and the dashed green line indicates the benchmark RCT hazard ratio. Error bars represent the standard error across configurations. MatchingEntropy Balancing IPTW 0.00 0.02 0.04 0.06 0.08 0.10 Absolute log-HR error STRASS MatchingEntropy Balancing IPTW 0.00 0.05 0.10 0.15 0.20 Absolute log-HR error BRT Standard (X)Augmented (X + U)Multivariable Analysis Figure 4: Comparison of treatment-effect estimation accuracy across RCT datasets. Bars represent the mean absolute log-HR error relative to the trial-reported hazard ratio, averaged across hyperparameter configurations. Results are shown for matching, entropy balancing, and inverse propensity weighting using either baseline covariates X or the augmented representation (X, ̃ U). The dotted line denotes multivariable Cox adjustment using baseline covariates only. Error bars denote the standard error across configurations. plying that observational estimates ought to become identical to the randomized estimate. Even with better control of con- founding, hazard ratios may still differ because they depend not only on confounding, but also on differences in eligibility cri- teria, treatment definitions, follow-up duration, and other fea- tures distinguishing real-world cohorts from trial populations. For example, in RPS, the benchmark STRASS trial used ab- dominal recurrence-free survival, a trial-specific endpoint that 11 Table 3: Cohort characteristics of the RCT datasets. TreatmentEvents DatasetNSurgery + periop. therapy Surgery alone Surgery + periop. therapy Surgery alone BRT16275 (46.3%)87 (53.7%)12 (16.0%)26 (29.9%) STRASS 265132 (49.8%)133 (50.2%)59 (44.7%)61 (45.9%) MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 1.0 Hazard Ratio STRASS MatchingEntropy Balancing IPTW 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Hazard Ratio BRT Standard (X)Augmented (X + U)MultivariableGround Truth Figure 5: Comparison of treatment-effect estimates across RCT datasets. Bars represent the mean hazard ratio (HR), averaged across hyperparameter configurations, for matching, entropy balancing, and inverse propensity weighting using either baseline covariates X or the augmented representation (X, ̃ U). Error bars denote the standard error across configurations. The dotted line indicates the multivariable estimate, and the dashed green line indicates the trial-reported ground-truth HR. Table 4: Validation approach 2: equivalence testing for RCT datasets. A “win” indicates that the augmentation with ̃ U passes the equivalence test, meaning the 95% confidence interval of the shift ∆ shift = log(HR aug )− log(HR RCT ) lies entirely within the pre-specified equivalence margin (± log(1.10)). MetricValue Equivalence wins4 / 6 Mean|∆ shift |0.081 cannot be directly recreated in real-world datasets, where the available outcome is local recurrence, a similar but not identi- cal endpoint. The relevant finding is therefore not exact numer- ical convergence, but a systematic shift toward more plausible estimates when an external benchmark is available. The second validation strategy, application to randomized datasets, provides an important complementary check. If ̃ U reflects latent prognosis rather than arbitrary statistical noise, it should be balanced across treatment arms under randomiza- tion and should not materially distort treatment-effect estimates when added to the analysis. This is what we observed. In that sense, the randomized analyses function as a negative control: they support the interpretation of ̃ U as a prognostic construct while arguing against the possibility that its apparent utility in observational cohorts simply reflects overfitting. The multicenter CRLM analyses extend the evaluation be- yond comparison with an external randomized benchmark and instead test whether ̃ U helps explain two distinct forms of between- center heterogeneity. The first concerns hidden bias in who re- ceives adjuvant chemotherapy within a center. If treatment allo- cation remains influenced by latent prognosis after adjustment for recorded covariates, then estimated chemotherapy effects may differ across centers partly for confounding-related reasons rather than because the true treatment effect is fundamentally different. The observed reduction in between-center dispersion of chemotherapy hazard-ratio estimates after balancing ̃ U sug- gests that part of the latent structure underlying chemotherapy selection was shared across cohorts and was captured by ̃ U. This is biologically and clinically plausible, because many of the factors that influence receipt of chemotherapy in routine practice—such as frailty, limited tolerance for systemic treat- ment, occult disease burden, and other incompletely measured aspects of patient fitness—are not specific to any one center but recur across cohorts. Even so, we did not expect identical haz- ard ratios across centers, because treatment itself was not fully standardized across cohorts, with differences in chemotherapy regimens, biologic agents, treatment sequencing, overall patient selection, and follow-up. The second multicenter question concerns hidden differences in patient case mix across centers. Here, the aim was not to re- duce confounding in treated-versus-untreated comparisons within a cohort, but to determine whether latent differences in progno- sis between center populations contributed to cross-center sur- vival differences beyond those explained by recorded covari- ates. In several center pairs, balancing on measured covari- ates alone widened the survival gap, indicating that observed case mix had partly masked underlying outcome differences. Subsequent incorporation of ̃ U narrowed these residual gaps in many settings, supporting the interpretation that ̃ U captures the 12 Baseline Covariates Latent U Configurations 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 Absolute SMD (STRASS) Baseline Covariates Latent U Configurations 0.00 0.05 0.10 0.15 0.20 0.25 Absolute SMD (BRT) Balance of Baseline Covariates and Latent Factor U in RCT Cohorts Figure 6: Standardized mean differences (SMDs) between treatment arms for observed covariates and estimated latent factors ̃ U in the RCT datasets. The SMD values for ̃ U are generally below or close to the commonly used balance threshold of 0.1 and are comparable to those of the observed covariates. Table 5: Cohort characteristics of CRLM multi-center datasets. DatasetEventNSurgery + chemotherapySurgery aloneEvents (Surgery + chemotherapy)Events (Surgery alone) Graz and Vienna UniversitiesDeath12888 (68.8%)40 (31.2%)47 (53.4%)12 (30.0%) Graz and Vienna UniversitiesRecurrence12888 (68.8%)40 (31.2%)68 (77.3%)27 (67.5%) Johns HopkinsDeath461296 (64.2%)165 (35.8%)160 (54.0%)116 (70.3%) Johns HopkinsRecurrence461296 (64.2%)165 (35.8%)183 (61.8%)78 (47.3%) Kumamoto UniversityDeath9046 (51.1%)44 (48.9%)23 (50.0%)26 (59.1%) Kumamoto UniversityRecurrence9046 (51.1%)44 (48.9%)39 (84.8%)38 (86.4%) University of BergenDeath22086 (39.1%)134 (60.9%)30 (34.9%)108 (80.6%) University of BergenRecurrence22086 (39.1%)134 (60.9%)48 (55.8%)112 (83.6%) Yokohama UniversityDeath17664 (36.4%)112 (63.6%)27 (42.2%)40 (35.7%) Yokohama UniversityRecurrence17664 (36.4%)112 (63.6%)51 (79.7%)81 (72.3%) Cleveland ClinicDeath338218 (64.5%)120 (35.5%)81 (37.2%)52 (43.3%) Cleveland ClinicRecurrence338218 (64.5%)120 (35.5%)156 (71.6%)70 (58.3%) MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS Standard (X)Augmented (X + U) Figure 7: Between-center heterogeneity in treatment-effect estimates for the CRLM multi-center cohort. Bars represent the mean pairwise absolute difference in log-hazard ratios (log-HR) across centers, averaged over all center pairs within each method and configuration. Results are shown for matching, entropy balancing, and inverse propensity weighting using either baseline covariates X or the augmented representation (X, ̃ U). Error bars denote the standard error across pairwise comparisons. 13 Table 6: Validation approach 3, first component: reduction in mean pairwise absolute deviation after adjusting for X, ̃ U for 6 CRLM centers. StatisticValue Cells with improvement (∆ > 0)6 / 6 Sign test p-value (one-sided)0.0156 Wilcoxon signed-rank p-value (one-sided)0.0156 Median cell-level ∆0.1343 Mean cell-level ∆0.1467 unobserved prognostic component of cross-center differences in patient case mix. At the same time, these differences were not expected to disappear entirely. After balancing on recorded covariates and ̃ U, any remaining survival differences are more plausibly attributable to residual case-mix differences not cap- tured by the framework or to true center-level differences in care including differences in surgical quality, perioperative manage- ment, surveillance, and follow-up. An additional strength of the framework is that its perfor- mance was not restricted to a single analytic specification. We evaluated it across multiple balancing strategies and prespeci- fied hyperparameter ranges, and summarized performance across parameter combinations rather than selecting a single favorable run. This reduces the likelihood that the apparent benefit of ̃ U-augmented analyses was driven by post hoc tuning. The framework also remained informative across different prognos- tic model inputs, including externally developed prognostic mod- els, suggesting that its usefulness does not depend on a single modeling pipeline. This is particularly important because, in medical cohorts and especially in surgical oncology cohorts, prognostic models often achieve only modest discrimination. Retroperitoneal sarcoma and colorectal liver metastases are good examples of settings in which few published models achieve a c-index above 0.7. 16, 17, 18 Demonstrating that the framework can still perform well in such settings is important for its broader applicability. At the same time, the scope of the framework should be stated clearly. It was designed to address residual confounding arising from latent prognostic heterogeneity, not to resolve all sources of bias in observational treatment com- parisons. In particular, it was not explicitly designed to correct immortal time bias which may be relevant in comparisons such as surgery alone versus surgery plus adjuvant therapy. Such bias can arise because patients who go on to receive adjuvant treatment must, by definition, survive event-free long enough after surgery to initiate therapy, thereby creating a potential ar- tificial advantage if outcomes are compared from the date of surgery. Accordingly, our findings should not be interpreted as showing that ̃ U augmentation corrects immortal time bias per se. Rather, they suggest that the framework can remain in- formative in real-world settings where multiple sources of bias coexist. In applications where immortal time bias is likely to be important, the framework should ideally be combined with methods developed specifically for that problem, such as land- mark or time-dependent analyses. The disease-specific patterns observed here also help clarify the settings in which the framework is most likely to add value. In GIST, conventional multivariable adjustment already shifted estimates toward the benchmark trial effect, consistent with the strong prognostic and treatment-selection relevance of routinely recorded variables such as tumor size, mitotic index, and tumor site. In contrast, in RPS, latent-factor augmentation appeared more informative, likely because prognostic models are weaker and radiotherapy allocation may depend more heavily on un- recorded clinical judgment. Indeed, both the published anal- ysis of the RPS dataset used in our experiments and our own findings showed that multivariable regression produced little change relative to the univariable estimate, suggesting that the commonly recorded prognostic factors were not the principal drivers of perioperative radiotherapy use. 10 Rather, selection for RT may have been influenced more strongly by unmeasured factors reflecting radiation oncologists’ clinical judgment, in- cluding subtler aspects of disease presentation not captured in the available covariates. Taken together, these findings suggest that the framework may be especially valuable in clinical set- tings where substantial residual prognostic heterogeneity per- sists despite adjustment for standard recorded covariates. The implications of this work extend beyond estimation of average treatment effects. If observational cohorts can be made more comparable not only with respect to measured covariates but also with respect to latent prognostic imbalance, they may provide a stronger basis for downstream analyses aimed at iden- tifying subgroups with differential treatment benefit. This is particularly important because real-world datasets are typically much larger than RCTs, which are costly and time-intensive, and may better reflect contemporary practice patterns and pa- tient selection. However, the intention of the proposed frame- work is not to replace randomized trials—either for estimating average treatment effects or for discovering subgroup-specific effects—but to supplement them in settings where trials are in- feasible, underpowered for subgroup analyses, or not fully rep- resentative of routine care. If confirmed in independent external datasets, this approach may support more credible treatment- effect estimation and more reliable downstream precision-oncology applications. Declaration of interests The authors declare no competing interests. Data sharing Patient data used in this study were obtained from partic- ipating institutions under data use agreements and cannot be publicly shared due to privacy restrictions. Acknowledgments The authors would like to thank Dr. Nefeli Bampatsikou for the excellent technical support and Ms. Eftychia Ioustine Margonis for her support. This research was supported by the NCI Cancer Center Support Grant P30 CA008748 (G.A.M.) and the Abdul Latif Jameel Clinic for Machine Learning in Health (D.B.). 14 Table 7: Validation 3, second component: Cross-center validation of the latent prognostic factor ̃ U. Among center pairs where adjustment for observed covariates X increased disagreement in adjusted 5-year survival, we evaluated whether adding ̃ U reduced the gap. Results are shown for the full cohort and for the surgery alone subset. The adjustment happens in the form of weights calculated to balance the pairs of centers with different balancing methods. These weights are then used in the Kaplan-Meier estimation. SubsetEventMethodPairs Evaluated Pairs Widened by X U Fixed (%) Binomial p Full cohortMortalityMatching3019 (63.3%)16/19 (84.2%)0.0022 Full cohortMortalityIPTW3015 (50.0%)13/15 (86.7%)0.0037 Full cohortMortalityEntropy Balancing1610 (62.5%)9/10 (90.0%)0.0107 Full cohortRecurrence Matching3010 (33.3%)8/10 (80.0%)0.0547 Full cohortRecurrence IPTW3018 (60.0%)15/18 (83.3%)0.0038 Full cohortRecurrence Entropy Balancing169 (56.2%)7/9 (77.8%)0.0898 Surgery alone MortalityMatching3012 (40.0%)9/12 (75.0%)0.0730 Surgery alone MortalityIPTW3014 (46.7%)14/14 (100.0%)0.0001 Surgery alone MortalityEntropy Balancing64 (66.7%)4/4 (100.0%)0.0625 Surgery alone Recurrence Matching3014 (46.7%)11/14 (78.6%)0.0287 Surgery alone Recurrence IPTW3017 (56.7%)14/17 (82.4%)0.0064 Surgery alone Recurrence Entropy Balancing65 (83.3%)5/5 (100.0%)0.0312 References [1] D. P. Nussbaum, C. N. Rushing, W. O. Lane, D. M. Car- dona, D. G. Kirsch, B. L. Peterson, D. G. Blazer, Preoper- ative or postoperative radiotherapy versus surgery alone for retroperitoneal sarcoma: a case-control, propensity score-matched analysis of a nationwide clinical oncology database, The Lancet Oncology 17 (7) (2016) 966–975. [2] S. Bonvalot, A. Gronchi, C. Le Péchoux, C. J. Swal- low, D. Strauss, P. Meeus, F. Van Coevorden, S. Stoldt, E. Stoeckle, P. Rutkowski, et al., Preoperative radiother- apy plus surgery versus surgery alone for patients with primary retroperitoneal sarcoma (eortc-62092: Strass): a multicentre, open-label, randomised, phase 3 trial, The Lancet Oncology 21 (10) (2020) 1366–1377. [3] D. Bertsimas, A. G. Koulouras, G. A. Margonis, The road to precision medicine, npj Digital Medicine 7 (1) (2024) 307. [4] J. Hainmueller, Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies, Political analysis 20 (1) (2012) 25–46. [5] P. R. Rosenbaum, Model-based direct adjustment, Journal of the American statistical Association 82 (398) (1987) 387–394. [6] P. K. Andersen, M. G. Hansen, J. P. Klein, Regression analysis of restricted mean survival time based on pseudo- observations, Lifetime data analysis 10 (4) (2004) 335– 350. [7] R. P. DeMatteo, K. V. Ballman, C. R. Antonescu, R. G. Maki, P. W. Pisters, G. D. Demetri, M. E. Black- stein, C. D. Blanke, M. von Mehren, M. F. Bren- nan, et al., Adjuvant imatinib mesylate after resection of localised, primary gastrointestinal stromal tumour: a randomised, double-blind, placebo-controlled trial, The lancet 373 (9669) (2009) 1097–1104. [8] D. Bertsimas, G. A. Margonis, S. Tang, A. Koulouras, C. R. Antonescu, M. F. Brennan, J. Martin-Broto, P. Rutkowski, G. Stasinos, J. Wang, et al., An inter- pretable ai model for recurrence prediction after surgery in gastrointestinal stromal tumour: an observational co- hort study, EClinicalMedicine 64 (2023). [9] J. S. Gold, M. Gönen, A. Gutiérrez, J. M. Broto, X. García-del Muro, T. C. Smyrk, R. G. Maki, S. Singer, M. F. Brennan, C. R. Antonescu, et al., Development and validation of a prognostic nomogram for recurrence- free survival after complete surgical resection of localised primary gastrointestinal stromal tumour: a retrospective analysis, The lancet oncology 10 (11) (2009) 1045–1052. [10] K. J. Kelly, S. S. Yoon, D. Kuk, L.-X. Qin, K. Duk- leska, K. K. Chang, Y.-L. Chen, T. F. Delaney, M. F. Brennan, S. Singer, Comparison of perioperative radiation therapy and surgery versus surgery alone in 204 patients with primary retroperitoneal sarcoma: a retrospective 2- institution study, Annals of surgery 262 (1) (2015) 156– 162. [11] D. Bertsimas, A. Koulouras, H. Nagata, C. Gao, J. Mizu- sawa, Y. Kanemitsu, G. A. Margonis, The road to clinical trial emulation, Research Square (2025) rs–3. [12] Y. Kanemitsu, Y. Shimizu, J. Mizusawa, Y. Inaba, T. Ham- aguchi, D. Shida, M. Ohue, K. Komori, A. Shiomi, M. Sh- iozawa, et al., Hepatectomy followed by mfolfox6 ver- sus hepatectomy alone for liver-only metastatic colorectal cancer (jcog0603): a phase i or i randomized controlled trial, Journal of Clinical Oncology 39 (34) (2021) 3789– 3799. 15 [13] Y. Kawaguchi, S. Kopetz, H. Tran Cao, E. Panettieri, M. De Bellis, Y. Nishioka, H. Hwang, X. Wang, C.-W. Tzeng, Y. Chun, et al., Contour prognostic model for pre- dicting survival after resection of colorectal liver metas- tases: development and multicentre validation study us- ing largest diameter and number of metastases with ras mutation status, British Journal of Surgery 108 (8) (2021) 968–975. [14] P. Pisters, L. B. Harrison, D. Leung, J. M. Woodruff, E. S. Casper, M. F. Brennan, Long-term results of a prospec- tive randomized trial of adjuvant brachytherapy in soft tis- sue sarcoma., Journal of Clinical Oncology 14 (3) (1996) 859–868. [15] P. B. Olthof, S. Buettner, N. Andreatos, J. Wang, I. M. Løes, D. Wagner, K. Sasaki, A. Macher-Beer, C. Kam- phues, I. Pozios, et al., Kras alterations in colorectal liver metastases: shifting to exon, codon, and point mutations, British Journal of Surgery 109 (9) (2022) 804–807. [16] S. Kokkinakis, I. A. Ziogas, J. D. Llaque Salazar, D. P. Moris, G. Tsoulfas, Clinical prediction models for prog- nosis of colorectal liver metastases: A comprehensive re- view of regression-based and machine learning models, Cancers 16 (9) (2024) 1645. [17] M. C. Tan, M. F. Brennan, D. Kuk, N. P. Agaram, C. R. Antonescu, L.-X. Qin, N. Moraco, A. M. Crago, S. Singer, Histology-based classification predicts pattern of recur- rence and improves risk stratification in primary retroperi- toneal sarcoma, Annals of surgery 263 (3) (2016) 593– 600. [18] K. Sasaki, J. Gagnière, A. Dupré, V. Ardiles, J. M. O’Connor, J. Wang, A. Moro, D. Morioka, S. Buettner, L. Gau, et al., Performance of two prognostic scores that incorporate genetic information to predict long-term outcomes following resection of colorectal cancer liver metastases: An external validation of the md anderson and jhh-msk scores, Journal of Hepato-Biliary-Pancreatic Sci- ences 28 (7) (2021) 581–592. 16 1. Sensitivity Analysis All sensitivity plots report the mean absolute log-HR error relative to the benchmark hazard ratio, averaged across hyperparam- eter configurations, with standard errors shown as error bars. 45678910 Number of bins 0.3 0.4 0.5 0.6 0.7 0.8 Absolute log-HR error GIST 45678910 Number of bins 0.4 0.6 0.8 Absolute log-HR error RPS 45678910 Number of bins 0.1 0.2 0.3 0.4 Absolute log-HR error CRLM 45678910 Number of bins 0.10 0.15 0.20 0.25 Absolute log-HR error BRT 45678910 Number of bins 0.04 0.06 0.08 0.10 0.12 0.14 Absolute log-HR error STRASS Sensitivity: Matching Bins Standard (X)Augmented (X + U) Figure 1.1: Sensitivity of treatment-effect estimation to the number of bins used in prognostic matching for the observational and RCT cohorts. The plotted values represent the mean absolute log-HR error relative to the benchmark hazard ratio, averaged across hyperparameter configurations. Results are shown with and without the latent factor ̃ U. Error bars denote the standard error across configurations. 45678910 Number of bins 0.3 0.4 0.5 0.6 0.7 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS 45678910 Number of bins 0.3 0.4 0.5 0.6 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS CRLM centers: sensitivity to matching bins Standard (X)Augmented (X + U) Figure 1.2: Sensitivity of treatment-effect estimation to the number of bins used in prognostic matching for the 6 CRLM centers. The plotted values represent the mean pairwise absolute difference of log-HR between centers, averaged across hyperparameter configurations. Results are shown with and without the latent factor ̃ U. Error bars denote the standard error across configurations. 2. Permutation and Ablation Experiments 2.1. Permutation Experiments To examine whether the improvement observed after incorporating the latent factor ̃ U could arise from mechanical artifacts rather than genuine latent structure, we conducted a series of permutation experiments. 17 0.010.020.030.040.05 IPTW clipping threshold 0.3 0.4 0.5 0.6 Absolute log-HR error GIST 0.010.020.030.040.05 IPTW clipping threshold 0.0 0.2 0.4 0.6 0.8 1.0 Absolute log-HR error RPS 0.010.020.030.040.05 IPTW clipping threshold 0.05 0.10 0.15 0.20 Absolute log-HR error CRLM 0.010.020.030.040.05 IPTW clipping threshold 0.10 0.12 0.14 0.16 0.18 0.20 Absolute log-HR error BRT 0.010.020.030.040.05 IPTW clipping threshold 0.070 0.075 0.080 0.085 0.090 0.095 Absolute log-HR error STRASS Sensitivity: IPTW Clipping Threshold Standard (X)Augmented (X + U) Figure 1.3: Sensitivity of inverse propensity weighting (IPTW) to the clipping threshold applied to the estimated propensity scores for the observational and RCT cohorts. The plotted values represent the mean absolute log-HR error relative to the benchmark hazard ratio across hyperparameter configurations. Results are shown for models using baseline covariates X and the augmented representation (X, ̃ U). Error bars denote the standard error across configurations. 0.010.020.030.040.05 IPTW clipping threshold 0.2 0.3 0.4 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS 0.010.020.030.040.05 IPTW clipping threshold 0.30 0.35 0.40 0.45 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS CRLM centers: sensitivity to IPTW clipping Standard (X)Augmented (X + U) Figure 1.4: Sensitivity of inverse propensity weighting (IPTW) to the clipping threshold applied to the estimated propensity scores for the 6 CRLM centers. The plotted values represent the mean pairwise absolute difference of log-HR between centers, averaged across hyperparameter configurations. Results are shown for models using baseline covariates X and the augmented representation (X, ̃ U). Error bars denote the standard error across configurations. 18 12 Entropy moments 0.0 0.2 0.4 0.6 Absolute log-HR error GIST 12 Entropy moments 0.0 0.2 0.4 0.6 0.8 1.0 Absolute log-HR error RPS 12 Entropy moments 0.00 0.05 0.10 0.15 0.20 0.25 Absolute log-HR error CRLM 12 Entropy moments 0.0 0.1 0.2 0.3 Absolute log-HR error BRT 12 Entropy moments 0.00 0.02 0.04 0.06 0.08 0.10 Absolute log-HR error STRASS Sensitivity: Entropy Moments Standard (X)Augmented (X + U) Figure 1.5: Sensitivity of entropy balancing to the number of balanced moments for the observational and RCT datasets. Bars represent the mean absolute log- HR error relative to the benchmark hazard ratio, averaged across hyperparameter configurations. Results compare models using baseline covariates X with those augmented by the latent factor ̃ U. Error bars denote the standard error across configurations. 12 Entropy moments 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS 12 Entropy moments 0.0 0.2 0.4 0.6 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS CRLM centers: sensitivity to entropy moments Standard (X)Augmented (X + U) Figure 1.6: Sensitivity of entropy balancing to the number of balanced moments for the 6 CRLM centers. Bars represent the mean pairwise absolute difference of log-HR between centers, averaged across hyperparameter configurations. Results compare models using baseline covariates X with those augmented by the latent factor ̃ U. Error bars denote the standard error across configurations. 19 46810 Number of neighbors (K) 0.2 0.3 0.4 0.5 0.6 Absolute log-HR error GIST 46810 Number of neighbors (K) 0.1 0.2 0.3 0.4 0.5 0.6 Absolute log-HR error RPS 46810 Number of neighbors (K) 0.025 0.050 0.075 0.100 0.125 0.150 0.175 Absolute log-HR error CRLM 46810 Number of neighbors (K) 0.05 0.10 0.15 0.20 0.25 0.30 Absolute log-HR error BRT 46810 Number of neighbors (K) 0.04 0.06 0.08 0.10 0.12 Absolute log-HR error STRASS Sensitivity to Neighbor Count (K) MatchingEntropy BalancingIPTW Figure 1.7: Sensitivity of the framework to the number of nearest neighbors (K) used when constructing the latent prognostic factor ̃ U for the observational and RCT datasets. Lines represent the mean absolute log-HR error relative to the benchmark hazard ratio across configurations for each balancing method. Error bars denote the standard error across configurations. 345678910 Number of neighbors (K) 0.2 0.3 0.4 0.5 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS 345678910 Number of neighbors (K) 0.3 0.4 0.5 0.6 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS CRLM centers: sensitivity to neighbor count (K) MatchingEntropy BalancingIPTW Figure 1.8: Sensitivity of the framework to the number of nearest neighbors (K) used when constructing the latent prognostic factor ̃ U for the 6 CRLM centers. Lines represent the the mean pairwise absolute difference of log-HR between centers, averaged across hyperparameter configurations for each balancing method. Error bars denote the standard error across configurations. 20 Not includedIncluded Prognostic score in distance 0.0 0.1 0.2 0.3 0.4 Absolute log-HR error GIST Not includedIncluded Prognostic score in distance 0.0 0.1 0.2 0.3 0.4 0.5 Absolute log-HR error RPS Not includedIncluded Prognostic score in distance 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Absolute log-HR error CRLM Not includedIncluded Prognostic score in distance 0.00 0.05 0.10 0.15 0.20 0.25 Absolute log-HR error BRT Not includedIncluded Prognostic score in distance 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Absolute log-HR error STRASS Impact of Adding Prognostic Score in Distance MatchingEntropy BalancingIPTW Figure 1.9: Impact of incorporating the prognostic score in the distance metric used to construct ̃ U for the observational and RCT datasets. Bars represent the mean absolute log-HR error relative to the benchmark hazard ratio across hyperparameter configurations for each balancing method. Error bars denote the standard error across configurations. Not includedIncluded Prognostic score in distance 0.0 0.1 0.2 0.3 0.4 0.5 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS Not includedIncluded Prognostic score in distance 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS CRLM centers: impact of prognostic score in distance MatchingEntropy BalancingIPTW Figure 1.10: Impact of incorporating the prognostic score in the distance metric used to construct ̃ U for the 6 CRLM centers. Bars represent the mean pairwise absolute difference of log-HR between centers, averaged across hyperparameter configurations for each balancing method. Error bars denote the standard error across configurations. 21 Cox (CF)EXT (2y)EXT (5y)RF (CF)XGB (CF) Prognostic score model 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Absolute log-HR error GIST Cox (CF)MSKRF (CF)XGB (CF) Prognostic score model 0.0 0.2 0.4 0.6 0.8 1.0 Absolute log-HR error RPS Cox (CF)OS EXTRF (CF)XGB (CF) Prognostic score model 0.000 0.025 0.050 0.075 0.100 0.125 0.150 Absolute log-HR error CRLM Cox (CF)MSKRF (CF)XGB (CF) Prognostic score model 0.0 0.1 0.2 0.3 0.4 Absolute log-HR error BRT Cox (CF)EXTRF (CF)XGB (CF) Prognostic score model 0.000 0.025 0.050 0.075 0.100 0.125 0.150 Absolute log-HR error STRASS Impact of Prognostic Score Model MatchingEntropy BalancingIPTW Figure 1.11: Sensitivity of results to the choice of prognostic score model used to estimate the baseline risk component for the observational and RCT datasets. Bars represent the mean absolute log-HR error relative to the benchmark hazard ratio across hyperparameter configurations for each balancing method. Error bars denote the standard error across configurations. Cox (CF)OS EXTRF (CF)XGB (CF) Prognostic score model 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS Cox (CF)OS EXTRF (CF)XGB (CF) Prognostic score model 0.0 0.2 0.4 0.6 0.8 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS CRLM centers: impact of prognostic score model MatchingEntropy BalancingIPTW Figure 1.12: Sensitivity of results to the choice of prognostic score model used to estimate the baseline risk component for the 6 CRLM centers. Bars represent the mean pairwise absolute difference of log-HR between centers, averaged across hyperparameter configurations for each balancing method. Error bars denote the standard error across configurations. 22 MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Absolute log-HR error GIST MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 1.0 Absolute log-HR error RPS MatchingEntropy Balancing IPTW 0.0 0.1 0.2 0.3 Absolute log-HR error CRLM Permutation: U Global Permutation: U Within-Arm Permutation: Y Within-Arm Standard (X) Ablation: Signed Outcome Ablation: Signed RMST Ablation: Random Neighbors Augmented (X + U) Multivariable Analysis Figure 2.13: Permutation and ablation test for the observational datasets. The latent factor ̃ U is randomly permuted across patients before balancing. Bars represent the mean absolute log-HR error relative to the benchmark hazard ratio. Performance degradation relative to the non-permuted case indicates that the predictive signal of ̃ U contributes to improved treatment-effect estimation. Error bars denote the standard error across configurations. MatchingEntropy Balancing IPTW 0.00 0.05 0.10 0.15 0.20 0.25 Absolute log-HR error STRASS MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 Absolute log-HR error BRT Permutation: U Global Permutation: U Within-Arm Permutation: Y Within-Arm Standard (X) Ablation: Signed Outcome Ablation: Signed RMST Ablation: Random Neighbors Augmented (X + U) Multivariable Analysis Figure 2.14: Permutation and ablation test for the RCT datasets. The latent factor ̃ U is randomly permuted across patients before balancing. Because treatment assignment is randomized, permuting ̃ U does not materially change treatment-effect estimates, providing a negative control for the method. Error bars denote the standard error across configurations. 23 MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on OS MatchingEntropy Balancing IPTW 0.0 0.2 0.4 0.6 0.8 Mean pairwise absolute difference of log-HR Effect of adjuvant chemotherapy on RFS Permutation: U Global Permutation: U Within-Arm Permutation: Y Within-Arm Standard (X) Ablation: Signed Outcome Ablation: Signed RMST Ablation: Random Neighbors Augmented (X + U) Figure 2.15: Permutation test for the CRLM multi-center analysis. The latent factor ̃ U is randomly permuted across patients prior to balancing. Bars represent the mean absolute pairwise HR difference across centers. The increase in cross-center variability after permutation indicates that the estimated ̃ U captures meaningful latent prognostic structure shared across cohorts. Error bars denote the standard error across center pairs. Permutation of Y (RMST) within treatment arms. In the first experiment, we permuted the pseudo-RMST outcome values Y within each treatment arm prior to computing ̃ U, while keeping the original (true) outcomes when estimating the final hazard ratios (HRs). This manipulation breaks the individual-level association between X and Y within each arm, and consequently between Y and ̃ U during the construction of ̃ U. Importantly, however, it preserves the marginal distribution of Y | T within each arm. In other words, the arm-level survival distribution remains intact, but the patient-specific mapping X i ↔ Y i is destroyed. If the method continued to produce improved HR estimates under this manipulation, it would suggest that the observed cor- rection arises merely from injecting outcome-derived quantities into the balancing step (a mechanical artifact) rather than from ̃ U capturing meaningful latent structure tied to the observed covariates. Instead, the experiment fails: the improvement disappears. This indicates that the specific pairing between X i and Y i is necessary for ̃ U to encode meaningful information. The latent factor is not functioning as a generic outcome-based adjustment; it requires the genuine individual-level relationship between covariates and survival to operate properly. Permutation of ̃ U within treatment arms. In the second experiment, we permuted ̃ U within each treatment arm after its construction. This preserves the distribution of ̃ U within each arm but breaks the individual-level mappings ̃ U i ↔ X i and ̃ U i ↔ Y i . A deterioration in results relative to the original method indicates that individual-level alignment matters: it is not sufficient for ̃ U to encode only arm-level imbalance. Rather, the patient-specific latent information contained in ̃ U i contributes materially to the balancing step. The observed degradation supports the interpretation that ̃ U carries individualized latent risk information rather than merely reflecting group averages. Global permutation of ̃ U. In the third experiment, we permuted ̃ U across the entire dataset, ignoring treatment arms. This destroys not only the individual-level relationships but also the arm-level distribution of ̃ U. Under this manipulation, ̃ U becomes unrelated to treatment assignment and loses any confounding-relevant structure. As expected, this permutation fails entirely, serving primarily as a completeness check. It confirms that the improvement observed in the main analysis depends on the structural relationship between ̃ U and treatment allocation. All permutation experiments fail to reproduce the original improvement. This strongly suggests that the method does not operate through arbitrary or mechanistic adjustment. In particular, the Y -permutation experiment directly addresses the potential Y → U → T → Y feedback concern: because permuting Y breaks the individual-level relationship while preserving arm-level survival distributions, the failure of the method under this manipulation demonstrates that incorporating Y in the construction of ̃ U does not automatically induce artificial correction. The improvement arises only when the genuine structure linking X, Y , and T is preserved. 2.2. Ablation Experiments To further understand which components of the ̃ U construction drive the results, we conducted several ablation studies. ̃ U = signed outcome. In this simplified variant, we replaced ̃ U with the sign of the outcome only, testing whether directional information alone (risk versus protection) is sufficient to drive the correction. 24 ̃ U = signed RMST. Here, we used the continuous signed RMST without neighbor-based differencing. This preserves outcome magnitude but removes the comparative (neighbor-based) component of the construction. Random neighbors (direction preserved). Finally, we computed ̃ U using randomly selected neighbors rather than nearest neighbors in X-space, while preserving directionality (i.e., selecting frailer neighbors for anchors with no events and vice versa). In this case, the covariates X no longer influence neighbor selection. The sign-based experiments produce unstable and often problematic results. In particular, in GIST and BRT, these ablations substantially distort the HR estimates. They are also inconsistent in the CRLM multi-center analysis. This indicates that directional information alone is insufficient; the magnitude and structure of the neighbor comparison are essential. Interestingly, the random-neighbor experiment yields results comparable to the full method in some datasets (e.g., CRLM), but performs substantially worse in others, particularly GIST and BRT. We interpret this heterogeneity as follows. In datasets where the observed covariates X have limited prognostic strength, the signal available for constructing meaning- ful neighborhoods is inherently weak. As a result, both nearest-neighbor and random direction-consistent constructions rely on similarly noisy information, leading to comparable performance. In contrast, in datasets such as GIST, where the covariates are strongly prognostic, accurately identifying nearest neighbors in X-space becomes essential. In these settings, replacing structured neighborhoods with random ones leads to a clear degradation in performance. Taken together, the permutation and ablation studies demonstrate that: • The method does not operate via mechanical outcome injection. • Individual-level alignment between X, Y , and U is necessary. • The neighbor-based construction contributes meaningfully beyond simple outcome transformations. • The effectiveness of nearest-neighbor selection depends on the prognostic strength of the observed covariates. These findings reinforce the interpretation of U as a data-derived proxy for latent, pre-treatment summary of unobserved con- founders, rather than an artifact of post-treatment outcome manipulation. 2.3. Statistical Results Table 2.1: Comparison of balancing methods across datasets. ∆ represents directional improvement in closeness to the benchmark RCT hazard ratio after adding ̃ U. Mean difference is reported on the log hazard ratio scale. Dataset MethodMean ∆ Std. Error GIST Matching0.3340.012 Entropy Balancing0.2560.022 IPTW0.3260.026 RPS Matching0.0840.008 Entropy Balancing0.5920.031 IPTW0.9030.007 CRLM Matching0.2780.006 Entropy Balancing0.1620.007 IPTW0.1620.004 Table 2.2: Equivalence testing of treatment effect differences across datasets. Mean difference is calculated as ∆ shift = log(HR aug )− log(HR RCT ). Dataset MethodMean Diff. (log HR) Std. Error95% CI BRT Matching0.0750.0085[0.058, 0.092] Entropy Balancing0.1830.0221[0.139, 0.227] IPTW0.0350.0163[0.002, 0.067] STRASS Matching-0.0910.0058[-0.1025, -0.796] Entropy Balancing-0.0380.0095[-0.057, -0.019] IPTW-0.0660.0079[-0.081, -0.050] 25 Table 2.3: Comparison of balancing methods across datasets. ∆ represents the reduction in mean pairwise absolute deviation after adjusting for X, ̃ U for 6 CRLM centers. Mean difference is reported on the log hazard ratio scale. EventMethodMean ∆ Std. Error Mortality Matching0.1290.0087 Entropy Balancing0.0980.0327 IPTW0.2240.0029 Recurrence Matching0.1390.0047 Entropy Balancing0.1690.0360 IPTW0.1210.0035 26