Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
65,103 characters · 20 sections · 33 citation commands
Synthetic Survival Control: Extending Synthetic Controls for “When-If” Decision
Understanding and modeling the time until an event of interest happens is a fundamental problem across various domains, ranging from patient survival in healthcare han-petal-2024, equipment failure in engineering fi15050153, criminal recidivism in criminal justice han2025fairnessalgorithmsracialdisparities, to customer churn gao2024causalcustomerchurnanalysis and dynamic pricing in insurance pmlr-v202-choi23c. Unlike classical machine learning algorithms focusing on binary or average outcomes, survival analysis explicitly models how risk evolves over time generally through a survival function, which represents the probability of remaining event-free up to a certain time point, or a hazard (risk) function, which describes the instantaneous probability of the event occurring at a given moment, conditional on not having experienced the event yet. Moreover, survival analysis accounts for censoring, where the event of interest is not observed within the study window or is preempted by a censoring event. As Samuel Beckett writes in Waiting for Godot, “Nothing happens, nobody comes, nobody goes, it's awful”–yet even the wait is informative. The absence of an event still revises our beliefs about underlying risk JMLR:v25:23-0888. This time-to-event perspective enables us to understand not just whether an event happens but also when and more specifically, how the probability of occurrence evolves as time passes.
In applied settings, survival analysis plays more than the aforementioned descriptive role through estimation of survival distributions; it can also be used inferentially to assess how interventions affect the timing and risk of events. For example, do novel cancer therapies prolong patient survival? Our work is motivated precisely by such a real-world problem in oncology: evaluating the effectiveness of new therapies for patients with relapsed or refractory T-cell lymphoma (TCL), a rare and aggressive form of blood cancer with historically poor prognosis han-petal-2024. Typically, cancer patients first receive an initial, or first-line, treatment. The disease is considered refractory if this initial treatment fails to work at all and called relapsed if it initially succeeds but the cancer later returns. In either scenario, patients then require an alternative, known as second-line, treatment. Recently, novel single agent (SA) have been adopted as choices for second-line therapies in the United States, while many other countries continue to rely primarily on conventional chemotherapy (CC). This staggered adoption across countries provides a natural quasi-experimental setting for investigating the causal effects of novel therapies on patient survival. Such investigations require constructing counterfactual survival functions, which are estimates of how patient survival trajectories would have evolved under alternative treatment strategies.
Inferring such counterfactual survival functions presents additional challenges beyond those faced in classical statistical analyses. First, survival analysis inherently deals with sparse data: individual-level right-censored data typically provides only a single observed time point—the minimum of the actual time-to-event and time-to-censoring, yet the goal is to recover full survival or hazard trajectories under different interventions. In the context of relapsed or refractory TCL, these challenges are further exacerbated. As a rare disease, TCL has limited participation from clinical centers, resulting in small sample sizes and magnified statistical uncertainty. Moreover, its aggressive progression makes it difficult to distinguish survival patterns between treatment groups over short follow-up periods, intensifying the need for more effective therapies and reliable estimation methods. Adding to these difficulties, the available data are observational and collected retrospectively across different countries and institutions. Treatment assignment is therefore non-random and occurs by country-level drug availability, making it particularly challenging to isolate causal effects from spurious associations.
To formalize this problem, it is helpful to situate it within a panel data setup. In classical panel data applications, the outcomes associated with a unit (e.g., country or institution) are repeated measurements of a scalar variable (a single real number) across multiple periods. Methods such as Synthetic Control abadie2010synthetic and Synthetic Interventions agarwal2024syntheticinterventions operate in such {\em data-rich} scenarios to leverage the structure among units for counterfactual estimations in the presence of confounding. However, our setup is {\em data-sparse} by the nature of time-to-event data: each individual often contributes only one observed time-to-event, or none if censored. Moreover, our target is not a scalar but an entire hazard function over time. This leads to the central question of our work: {\em Can methods like Synthetic Control, originally designed for scalar outcomes in data-rich scenarios, be generalized to handle functional, time-evolving survival processes under censoring in data-sparse setups?}
In this work, our main contribution is to address this question by bridging the gap between the literature on Synthetic Control and Survival Analysis. We propose a panel data setting for survival analysis where multiple units are present and each unit experiences one of the possible treatments in each period. Within a unit and period, we observe multiple observations of, potentially censored, time-to-event. The potential outcomes associated with units in periods under a treatment can be characterized through “potential survival functions”, which can be factorized into (latent) characteristics associated with units, periods, and treatments. This factorized representation is a natural extension of classical survival models such as Cox, where the log-hazard decomposes into additive contributions from baseline risk, covariates, and treatment. Generalizing this decomposition to higher dimensions suggests a low-rank or factorized structure, with latent unit and period factors playing the role of unobserved confounders. Crucially, we allow for treatment assignment to arbitrarily depend on these (latent) factors of, thereby permitting (unobserved) confounding.
Within this causal framework, we consider a causal estimand analogous to that in the traditional framework of synthetic control (see Section (ref)). For such a causal estimand, we establish identification under minimal conditions (see Theorem (ref)). To estimate the causal estimand, we propose the method of Synthetic Survival Control (SSC) (see Section (ref)). We establish finite sample consistency guarantee (and hence inference) property for such an estimator (see Theorem (ref)).
To demonstrate the soundness of the framework, we first validate it on synthetic data. We apply it to a real-world scenario in the context of evaluating outcomes associated with TCL cancer treatment. Specifically, through collaboration with leading cancer research institutes, we obtain clinical data across multiple countries (units) over different periods. For each patient, typically, there are multiple stages of treatment. We consider two stages (called first-line and second-line). Across all countries, patients receive the treatment CC (called control) in stage one and two, with the exception of the United States where in the second stage, patients receive SA (called intervention). In each stage, the outcome of interest is the time to treatment failure, defined as either the initiation of a new line of therapy or death. Therefore, this naturally forms a panel setup for survival analysis with two periods (first stage and second stage of treatment). As it happens, the dataset also contains patients in the US who receive CC in the second stage as well. We withhold such data and use SSC to do counterfactual estimation for patient outcomes in the US under CC in the second stage. We find that SSC faithfully recovers such counterfactual hazard function and thus validates the method. This is particularly significant given all the data challenges present in such observational clinical data.
Together, our work illustrates that counterfactual estimation on survival functions is feasible even with sparse, censored data when leveraging group-level quasi-experiments and survival-specific estimands. By integrating time-to-event modeling with causal inference, our work contributes to a growing literature on counterfactual survival analysis. In summary, this work offers a practical framework for survival-based counterfactual estimation, bridging methodological gaps and addressing real-world questions.
Section (ref) reviews related work on causal inference for survival outcomes and synthetic control methods. Section (ref) introduces our panel framework with low-rank structure, including identification, SSC estimator, and theoretical guarantees. Section (ref) describes the TCL dataset and and presents our empirical findings on treatment efficacy. Section (ref) discusses implications, limitations, and future direction.
Synthetic Control, Panel Methods and Matrix Completion. Synthetic control methods reconstruct counterfactuals for a single treated unit by optimally weighting a pool of control units, under the assumption that a linear combination of controls can approximate the treated unit's pre‑treatment path. abadie2003economic first applied this idea in economics to assess the costs of conflict in the Basque Country, and abadie2010synthetic formalized it in a statistical framework for comparative case studies. Subsequent work has extended and unified synthetic control with other panel‑data approaches: doudchenko2016balancing proposed a unified framework that integrates elements from synthetic controls, difference-in-differences (DiD), and regression-based approaches to address limitations inherent in each method when applied separately, and xu2017generalized proposed a generalized synthetic control via interactive fixed effects models. More recently, arkhangelsky2021synthetic blended synthetic control with difference‑in‑differences in the synthetic DiD estimator, improving robustness to heterogeneous trends. agarwal2024syntheticinterventions extended the synthetic control estimator to a synthetic intervention method which addresses the complementary question of what could have happened had a control unit received the treatment. Extending scalar outcomes to distributions, feitelberg2025distributionalmatrixcompletionnearest developed a latent-factor nearest-neighbor method for imputing missing empirical distributions in a sparsely observed matrix by leveraging Wasserstein geometry, but do not model temporal structure, censoring, or interventions.
Survival Analysis. Classical survival analysis focuses on modeling hazard and survival functions under censoring, with foundational works including Kaplan-Meier's non-parametric estimator kaplan1958nonparametric and Cox's proportional hazards model cox1972regression. These works typically reply on assumptions like random assignment of treatment, non-informative censoring, no noncompliance and so on. Along with the development of causal inference, these strong assumptions are gradually relaxed. robins1992correcting developed rank‑preserving structural failure time models to correct for noncompliance; ipcw proposed the inverse probability of censoring weighted (IPCW) method to address the missing data problem due to dependent censoring, and marginal structural models in robins2000marginal use inverse‑probability weighting to estimate time‑varying treatment effects.
Causal Survival Analysis. Several recent studies have explored how to adapt causal inference methods to time-to-event outcomes. One line of work focuses on adapting DiD to the survival setting wuHazardLinearProbability2022, wooldridgeSimpleApproachesNonlinear2023, deanerCausalDurationAnalysis2024. cui2023 developed doubly robust and tree-based survival estimators. Other works take the perspective of matrix completion in past literature. For example, gao2024causalcustomerchurnanalysis proposed a tensor block hazard model using 1-bit tensor completion for binary churn outcomes and intervention clustering but their method was focused on modeling the underlying propensity parameters that lead to binary events. This is different from our approach of adapting synthetic control's weighting philosophy to hazard trajectories under a panel framework with low-rank structure.
In this section, we develop a causal framework for estimating counterfactual survival outcomes using panel data, with a motivating application introduced in Section (ref), to evaluate the impact of novel cancer therapies on patient survival.
We begin by describing the causal survival analysis problem in its most fundamental form. Our TCL study is naturally panel-like: multiple individuals are observed across countries and lines of therapy. Assume there are $N$ units (e.g. countries) and $P$ periods (e.g. $p \in \{0,1\}$ first-line vs. second-line therapy).
\paragraph{Data Generating Process and DAG.}
To formalize the data generating process, we describe it in a directed acyclic graph (DAG) shown in Figure (ref)(a). For each unit $n\in [N]$ and period $p \in \{0, 1\}$, treatment assignment $D_{p, n} \in \{0, 1 \} $ acts on the event time $\tau_{p, n} \in \mathbb{R}^+$ (through a hazard process defined below). Unit- and period-level covariates $X_{p, n}$ influence both treatment and outcome. However, we do not observe $\tau_{p, n}$ directly, as it can be masked by the censoring time $C_{p, n} \in \mathbb{R}^+$ if $C_{p, n} < \tau_{p, n}$; we observe time $T_{p, n} = \min\{\tau_{p, n}, C_{p, n}\}$ and the event indicator $\Delta_{p, n} = \mathbb{I}\{\tau_{p, n} \leq C_{p, n}\}$.
\paragraph{Potential Outcomes.}
Following the potential outcomes framework of nayman1923 and rubin1974 , we define the potential hazard function \[ h^{(d)}_{p, n}(t) := \lim_{\Delta t \to 0} \frac{\mathbb{P}(t \leq \tau^{(d)}_{p, n} < t+\Delta t \mid \tau^{(d)}_{p, n} \geq t)}{\Delta t}, t \in [0, \widetilde{\tau}] \] which represents instantaneous risk of event occurrence under treatment $d$ at time $t$ for unit $n$ in period $p$ since time-to-event tracking, as opposed to the typical calendar time. Because late-time survival data are sparse and unstable, we restrict attention to a finite evaluation horizon $\widetilde{\tau} > 0$.
As such, $h^{(d)}_{p, n}(t)$ induces the potential survival function $S^{(d)}_{p, n}(t) = \exp(-\int_0^t h^{(d)}_{p, n}(s) ds)$ and density $f^{(d)}_{p, n}(t) = h^{(d)}_{p, n}(t) S^{(d)}_{p, n}(t)$ of the potential event time $\tau^{(d)}_{p, n}$.
In our TCL example, $d=0$ denotes conventional chemotherapy (CC), while $d=1$ corresponds to novel single-agent (SA) therapy. Each patient can only receive one type of therapy in a treatment line and thus, we can never observe all possible potential outcomes – this is the fundamental challenge of causal inference.
Given the model encoded in Figure (ref), we define our causal estimand formally. For a specific unit $n$, our target is $\theta_n(\cdot) = h^{(0)}_{1, n}(\cdot)$, which captures the potential hazard trajectory in the post-period $p=1$ under control ($d = 0$). Returning to our motivating example, the estimand corresponds to the potential hazard trajectory in the second-line therapy under the control treatment CC.
For identification and estimation, we have the following assumptions.
Assumption (ref) ensures that factual hazards match counterfactual potential hazards under the assigned treatment, and there are no cross-individual spillovers.
Assumption (ref) allows us to treat censored subjects as if they were removed at random from the risk set, so the remaining risk set stays representative \footnote{Note that we assume this mainly for readability and clarity. A natural extension of our current framework using techniques like inverse probability censoring weighting (IPCW)ipcw should work if this assumption is relaxed.}.
In the setting when we have a rich collection of covariates, a natural modeling assumption is that all pertinent covariates (and hence confounders) are observed. We first discuss how observed data can be utilized for estimating potential outcomes of interest in a simple single-period setting. This allows us to better understand the data generating process and separate the fundamental difficulties of causal survival analysis from those specific to the panel structure.
For identification, we also need the following positivity assumption.
We now present the identification result.
See Proof in Appendix (ref).
\paragraph{Data.} In period \(p\in\{0,1\}\), we observe survival data for each unit (country in our example) \(n\in[N]\) consisting of \(K_{p, n}\) independent observations (patients in our example). For observation \(i=1,\dots,K_{p, n}\), they contribute the assigned treatment $D_{p,n,i}$, features $X_{p, n, i}$ and censored time‐to‐event observation \((T_{p,n,i}, \Delta_{p,n,i})\).
\paragraph{Time-to-event Modeling.} Based on Proposition (ref), we boil down the problem of estimating the potential hazard function to estimating the observational conditional hazard $h_{p, n}(t\mid D_{p, n}=0, X_{p, n})$ and survival functions $S_{p, n}(t\mid D_{p, n}=0,X_{p, n})$ . This can be done using standard survival models with certain modeling structures. For example, the Cox model encodes the additional structure assumption that covariates contribute to the hazard function multiplicatively, i.e.,
where $\widetilde h(t)$ is a baseline hazard.
In fact, Equation (ref) implies that the log-hazard decomposes additively into a baseline term $\alpha(t)=\log \widetilde h(t)$, covariate contribution $\beta^\top X$, and treatment effect $\gamma d$. More generally, this suggests representing the transformed hazard surface by an additive (or low-rank) decomposition: $ \alpha(t) + \beta(X) + \gamma(d). $ Such a representation naturally generalizes to higher-dimensional arrays of transformed hazards indexed by units, periods, treatments, and time, motivating the low-rank structure formulation developed in Section (ref) later.
\paragraph{Handling Censoring.} Note that Assumptions (ref) and (ref) ensure that $h_{p, n}(t\mid D_{p, n}=d,X_{p, n})$ and $S_{p, n}(t\mid D_{p, n}=d,X_{p, n})$ are estimable from right-censored data. For example, if the event-time model has parameter $\eta$, a probability density function $f(\cdot)$ and a survival function $S(\cdot)$, and the censoring model has parameter $\psi$, a probability density function $\widetilde f(\cdot)$ and a survival function $\widetilde S(\cdot)$, then under Assumption (ref), we have: for unit $n$, period $p$ and $dt \to 0_+$,
Hence, likelihood-based survival methods like the Cox model can continue working with the observed-data's log-likelihood
Note that the $\psi$-terms and $\eta$-terms nicely decompose. As such, maximizing over $\eta$-terms only is equivalent to maximizing $\ell(\eta, \psi)$ with observed data.
\paragraph{Estimation Procedure.} We now formally state the estimation procedure for $\theta_n(t) = h^{(0)}_{1, n}(t)$ for all $t$ in some fixed time horizon $[0, \widetilde{\tau}]$ for a specific unit $n$. We introduce the usage of the Cox proportional hazards model and present its consistency result in Section (ref). Note that we are utilizing known results in literature for Cox proportional hazards model. In its place, other methods from literature can be utilized in its place and prior results from literature may provide similar guarantees.
(1) Step 1: Fit a Cox proportional hazards model \[ h_{1, n}(t \mid D_{1, n},X_{1, n}) = \widetilde h(t)\exp(\beta^\top X_{1, n} + \gamma D_{1, n}), \] obtaining $(\widehat\beta, \widehat {\widetilde \Lambda}(t) := \int_{0}^t \widehat{\widetilde h}(s)ds)$ via the Breslow estimator breslow1975.
(2) Step 2: Predict potential conditional hazard and survival for each individual by setting the treatment to be the control. Given $n$ and indices of observed data $i = 1,\dots, K_{1, n}$, regardless of their observed treatment,
(3) Step 3: Plug-in. \[ \widehat h^{(0)}_{1, n}(t) = \frac{ \sum_{i} \widehat h_{1, n, i}^{(0)}(t) \widehat S_{1, n, i}^{(0)}(t) }{ \sum_{i} \widehat S_{1, n, i}^{(0)}(t) }. \]
Our estimator for the marginal potential hazard $h_{1,n}^{(0)}(t)$ builds on the Cox proportional hazards model introduced in Section (ref). Under standard regularity assumptions and independent censoring, the Cox partial-likelihood estimator $\widehat{\beta}$ is asymptotically consistent and asymptotically normal, and the Breslow estimator of the baseline cumulative hazard $\widehat{\widetilde \Lambda}(t)$ converges uniformly to $\widetilde \Lambda(t)$ on compact time intervals cox1972regression, Tsiatis1981. Consequently, \[ \widehat{S}_{1, n}(t \mid D=0, X_{1,n,i}) = \exp\!\left( -\widehat{\widetilde \Lambda}(t)\exp(\widehat{\beta}^\top X_{1,n,i}) \right) \xrightarrow{p} S_{1, n}(t \mid D=0,X_{1,n,i}) \] uniformly in $t \in [0, \widetilde{\tau}]$. Thus, our estimator of the marginal counterfactual hazard, \[ \widehat{h}^{(0)}_{1,n}(t) = \frac{ \sum_i \widehat{h}_{1, n}(t\mid D=0,X_{1,n,i}) \widehat{S}_{1, n}(t\mid D=0,X_{1,n,i}) }{ \sum_i \widehat{S}_{1, n}(t\mid D=0,X_{1, n,i}) }, \] inherits consistency by the continuous mapping theorem.
In practice, only a subset of the covariates may be observed. The unobserved components confound both treatment assignment and outcomes, making direct identification of causal survival functions challenging, if not infeasible. Moreover, Assumption (ref) may be violated due to limited availability of certain treatment. These motivate need for additional structure for our setting with unobserved confounding.
We notice that persistent latent characteristics of units (e.g., country-level differences in healthcare infrastructure, patient mix, or access to medical technology) influence both treatment assignment and outcomes. Similarly, there are persistent latent characteristics of period (e.g., second-line patients have specific biological, disease-specific conditions that are persistent across countries). These characteristics or factors can be considered to be stable over the duration of the study. This suggests {\em factorization} of characteristics as $X_{p, n} = (V_n, U_p)$ where $V_n$ represents unit characteristics or factors and $U_p$ represent period characteristics or factors. We shall consider $V_n$ and $U_p$ to be latent or unobserved. As we discussed near the end of the manuscript, this framework has potential to extend to additional observed covariates and is left for future work.
Because each unit can only adopt one treatment per period, at most one of the potential hazard functions is observed. Estimating counterfactual survival trajectories is therefore equivalent to imputing missing hazards across the unit-period-treatment tuples, which is illustrated in Figure (ref).
\paragraph{Factorization of Survival Function.} The distribution of time-to-event $\tau_{n,p}$, i.e. the survival function, depends on $X_{p,n}$. We imposed restriction that $X_{n, p}$ separates or factorizes into unit specific component $V_n$ and period specific component $U_p$. More precisely, we impose the following {\em low-rank} factorization structure on survival function that factorizes into components that unit specific and period specific. We state the formal assumption first and then provide its justification.
To understand why Assumption (ref) follows naturally once we impose structural separation of $X_{n, p}$ into unit specific component $V_n$ and period specific component $U_p$, let us recall the Cox model. In such a model, we have that logarithm of hazard function having form
where $\widetilde h(\cdot)$ is a baseline hazard, $\beta$ being model parameter, $\gamma$ capturing treatment effect. Writing $X_{p,n} = [V_n U_p]$ as a vector with corresponding element-wise split of $\beta = [\beta_1 \beta_2]$, we have
\footnotetext{This is similar to two-way fixed effect in econometrics literature.} Denote $a_n = \exp(\beta_1^\top V_n)$, $b_p = \exp(\beta_2^\top V_p)$, and using the fact that survival function is exponential of negative integral of hazard function, we obtain (for $d = 0$)
where $\tilde b_p(t) = - \int_{0}^t b_p \widetilde h(s)ds$.
We consider the following elementary fact (see Appendix (ref) for a proof):
This suggests that
By writing $v_n = [a_n^k: 0\leq k\leq r-1] \in \mathbb{R}^r$ and $u_p(t) = [\frac{1}{k!}\tilde b_p^k(t): 0\leq k \leq r-1]$ with $u_p: \mathbb{R}_{\geq 0} \to \mathbb{R}^r$, we obtain
That is, Assumption (ref) is (approximately) satisfied. More generally, for a number of parametric models, it can be checked that $S^{(0)}_{p,n}(t)$ is a smooth nonlinear transformation of a bilinear form similar to (ref). For such settings, like Lemma (ref), this would result into approximate low-rank factorization like (ref) (also see Proposition 1, abadie2025causalinferenceframeworkdata). This provides natural justification for Assumption (ref).
\paragraph{Selection on Latent Factors.} Analogous to the unconfoundedness assumption in classical causal inference, we require that once latent factors are conditioned upon, treatment assignment is independent of potential hazards:
Together, Assumptions (ref) and (ref) provide the structural basis for imputing missing entries of $\mathcal{S}$ and thereby estimating counterfactual hazard trajectories.
\paragraph{Causal Estimand.} Building upon the definition in Section (ref), our target causal estimand is $ \theta_{n}(\cdot) = S^{(0)}_{1, n}(\cdot), $ for a given treated unit $n \in \mathcal{I}^{(1)}$. Thus it captures the potential transformed hazard trajectory for an treated unit $n \in \mathcal{I}^{(1)}$ under control during the post-treatment period, conditional on latent factors. Returning to our motivating example, the estimand corresponds to the true transformed hazard trajectory from second-line treatment onward had the US continued using CC instead of adopting novel SA.
We now turn to the identification question: under what assumptions can we recover the unobserved counterfactual trajectories $\theta_{n}(\cdot)$ for treated units $n\in \mathcal{I}^{(1)}$?
Assumption (ref) states that within the space of latent unit factors, the treated unit $n$ can be represented as a weighted average of control during the post-intervention period and such weights are period- and time-invariant.
It is a direct relaxation of the “convex span” condition in synthetic control and utilized in prior works such as synthetic interventions. It ensures that counterfactual outcomes for experimental units can be approximated as weighted averages of control group outcomes. In the context of our motivating example, Assumption (ref) requires USA's latent factors to be expressed as a weighted sum of the remaining countries' latent factors in the post-intervention period.
It is worth noting that Assumption (ref) implies that all unit-specific latent factors \(\{v_n\}_{n=1}^N \subset \mathbb{R}^r\) lie in an \(r\)-dimensional subspace for some fixed \(r\). Let \(\mathcal{V}_i = \mathrm{span}\{v_1, \ldots, v_i\}\) denote the subspace generated by the first \(i\) vectors. Then the sequence \((\dim \mathcal{V}_i)_{i=1}^N\) is nondecreasing and bounded above by \(r\). Hence, there exist at most \(r\) indices \(i_1, \ldots, i_r \in [N]\) such that $ \dim(\mathcal{V}_{i_j}) = \dim(\mathcal{V}_{i_j-1}) + 1,$ i.e., vectors \(v_{i_j}\) that introduce new linearly independent directions. Equivalently, all but at most \(r\) of the vectors \(\{v_n\}_{n=1}^N\) lie in the span of their predecessors. Thus, with $N$ large enough, Assumption (ref) should be satisfied for almost all (with exception of $r$ units in the worst case).
We now state the identification result.
See proof in Appendix (ref). This means the counterfactual post-treatment transformed hazard for treated unit $n$ under control is identified by a linear combination of post-treatment transformed hazards from control units, using weights determined by pre-treatment similarity in latent space.
\paragraph{Data.} In period \(p\in\{0,1\}\), we observe survival data for each unit \(n\in[N]\) consisting of \(K_{p, n}\) independent observations. For observation \(i=1,\dots,K_{p, n}\), they contribute the assigned treatment $D_{p,n,i}$ and censored time‐to‐event observation \((T_{p,n,i}, \Delta_{p,n,i})\).
\paragraph{Estimation Procedure.} Consider a treated unit $n \in \mathcal{I}^{(1)}$ of interest. Let $\mathcal I^{(0)}$ be $N_0$ control donor units and $\mathcal I^{(1)}$ the treated unit(s).
(1) Time-to-event estimation. For every donor $m\in\mathcal I^{(0)}$ and both periods $p\in\{0,1\}$, and for the target treated unit $n$ in the pre–period $p=0$, fit a Kaplan–Meier estimator to obtain \[ \widehat{S}^0_{p,l}(t),\quad t\in[0,\widetilde\tau],\qquad l\in\{n\}\cup\mathcal I^{(0)}. \] This inherently handles censoring in patient-level data $\{T_{p,n,i}, \Delta_{p,n,i}\}$ to estimate the underlying hazard functions for each unit-period combination. The accuracy of these estimates depend on the number of observations $K_{p,n}$ within each unit-period. See Section (ref) for more details of time-to-event estimation with Kaplan-Meier empirically.
(2) Common grid subsampling. Choose a finite grid $\mathcal T=\{t_1,\ldots,t_{T_0}\}\subset[0,\widetilde\tau]$ and subsample from the estimated survival functions to form vectors for all $l\in\{n\}\cup\mathcal I^{(0)}$, \[ \widehat{S}_{p,l} := \big(\widehat{S}_{p,l}(t_1),\ldots,\widehat{S}_{p,l}(t_{T_0})\big)^\top . \]
(3) Weights from the pre–period. Stack all control units’ pre–period columns into the $T_0\times N_0$ matrix $\widehat{S}_{0,\mathcal I^{(0)}}:=[\widehat{S}_{0,m}]_{m\in\mathcal I^{(0)}}$ and the treated unit’s pre vector $\widehat{S}_{0,n}(\cdot)$. We wish to utilize the Principal Component Regression (PCR) algorithm (cf. agarwal2024syntheticinterventions) to learn weights. Specifically, compute the top-$r_0$ singular value decomposition of $\widehat{S}_{0,\mathcal I^{(0)}}$ (with $r_0$ chosen by a gap rule, elbow, or cross-validation as explained in agarwal2024syntheticinterventions): \[ \widehat{S}_{0,\mathcal I^{(0)}} = \sum_{i=1}^{r_0} \widehat{s}_i \widehat{u}_i \widehat{v}_i^\top, \qquad \widehat{V}_0 := [\widehat{v}_1,\ldots,\widehat{v}_{r_0}] \in \mathbb{R}^{N_0 \times r_0}. \]
We estimate weights $\widehat{w}$ by projecting onto the rank-$r_0$ space and imposing the constraint: \[ \widehat{w} \in \arg\min_{w\in\mathbb R^{N_0}} \big\| \widehat{S}_{0,n} - \widehat{S}_{0,\mathcal I^{(0)}} w \big\|_2^2 \quad\text{s.t.}\quad (I - \widehat V_0\widehat V^{\top}_0)w = 0, \] which has a closed form solution $\widehat{w} = \left( \sum_{i=1}^{r_0} 1/{\widehat{s}_{i}} \widehat{v}_{i} \widehat{u}_{i}^{\top} \right) \widehat{S}^0_{0,n}.$
{(4) Post--period extrapolation and inverse transform.} Form the counterfactual transformed trajectory for the experimental unit $n$ in the post period: \[ \widehat{S}^{(0)}_{1,n}(t) = \sum_{m\in\mathcal I^{(0)}} \widehat{w}_{n,m} \widehat{S}^0_{1,m}(t), \qquad t\in[0,\widetilde\tau]. \]
In the context of TCL example, the estimator learns a set of weights that minimizes pre-intervention discrepancy between the US and a weighted combination of countries that maintained CC throughout the study period after pre-processesing $\widehat{S}_{0,\mathcal I^{(0)}}$ by obtaining its rank-$r_0$ approximation. Then this weighted average forms our synthetic USA. Accordingly, the estimator applies the learned weights to the post-intervention outcomes of the control countries to estimate the counterfactual hazards trajectory of the US under continued CC after drug offering shift.
Note that in this work, we introduce the use of the Kaplan-Meier estimator for the first step time-to-event estimation. The subsequent finite-sample consistency guarantee is also focused on this approach. Other time-to-event estimators can be applied as well and their finite-sample consistency property analyses are left for future works.
We analyze the consistency of our SSC estimator.
First, we state additional assumptions for such analysis. Let \( \mathcal{S}_{0, \mathcal{I}^{(0)}} \in \mathbb{R}^{T_0 \times N_0} \) be the pre-intervention matrix for units remaining under control, where $N_0=|\mathcal{I}^{(0)}|$ is the number of control units and $T_0$ is the number of evaluation timestamps in $ \mathcal{T} \subseteq [0, \widetilde{\tau}]$. We also note that the survival functions $S_{p,n}^{(0)}(\cdot)$ are naturally bounded between 0 and 1.
We now state the main result characterizing the consistency of our estimator.
See proof in Appendix (ref).
To begin with, we consider a simulation setting where we generate observations as per the data generating model introduced in this work. As argued, the SSC method should be able to provide a consistent estimator for such a setting, and we can verify it given the prior knowledge of the data generating process. As we do so, we shall aim to validate a few aspects. First, SSC method is able to provide consistent estimation across a number of popular parametric models from survival analysis literature. Since the SSC method is agnostic to the underlying data generating mechanism, this would suggest robustness of SSC across different types of data characteristics. Second, understanding the role of different aspects of the data on method performance. While theoretical results provide performance guarantees, such data-driven analysis helps set guidelines for using SSC in empirical studies.
\paragraph{Simulation Setup.}
In each simulation, we consider $N=20$ units in the same panel setup introduced in Section (ref), with one unit designated as the treated unit that switches from the control treatment $D=0$ in the pre-period $p=0$ to the intervention treatment $D=1$ in the post-period $p=1$, and all other units remaining under control ($D=0$) in both periods. For each period--unit pair $(p, n)$ we generate $K_{p,n}\in\{100,300,700\}$ independent samples of right–censored time–to–event data \[ (T_{p,n,i},\Delta_{p,n,i}) = \big(\min\{\tau_{p,n,i},C_{p,n,i}\}, \mathbf{1}\{\tau_{p,n,i}\le C_{p,n,i}\}\big), \] where censoring times $C_{p,n,i}\sim\mathrm{Exp}(\nu)$ are independent of $\tau_{p,n,i}$ conditional on $\big(D_{p,n},U_p,V_n\big)$, satisfying the independent censoring Assumption (ref). We set the censoring rate to be below $10\%$ by choosing $\nu$ so that the empirical censoring proportion is approximately $10\%$ across designs. Latent factors $(V_n,U_p)$ are drawn i.i.d.\ from $\mathcal N(0,I_r)$ with $r=4$ and then fixed throughout the simulation.
For the event-time models, we consider two canonical survival specifications.
1.\ Cox proportional hazards model. For each unit $n$ and period $p\in\{0,1\}$, we generate potential event times under control from \[ h^{(0)}_{p,n}(t) = \widetilde h(t)\exp\!\big(\beta_1^\top V_n + \beta_2^\top U_p\big), \qquad \widetilde h(t)\equiv \lambda > 0, \] with $\lambda$ fixed (we take $\lambda = 0.05$ in all experiments) and the coefficient vectors $\beta_1,\beta_2\in\mathbb R^r$ are drawn as $\beta_1,\beta_2 \overset{\text{i.i.d.}}{\sim} \mathcal N(0,I_r)$. The implied survival function under control is thus \[ S^{(0)}_{p,n}(t) = \exp\big[-\lambda t \exp\big(\beta_1^\top V_n + \beta_2^\top U_p\big)\big], \] and potential event times $\tau^{(0)}_{p,n}$ are sampled via inverse transform: \[ \tau^{(0)}_{p,n} = \frac{-\log W}{\lambda \exp(\beta_1^\top V_n + \beta_2^\top U_p)}, \qquad W\sim\mathrm{Unif}(0,1). \]
2.\ Aalen’s additive hazards model. For each unit $n$ and period $p\in\{0,1\}$, we generate potential event times under control from \[ h^{(0)}_{p,n}(t) = \widetilde h'(t)+\beta_1'^\top V_n+\beta_2'^\top U_p, \qquad \widetilde h'(t)\equiv \lambda_0 > 0. \] with the coefficient vectors $\beta_1',\beta_2'\in\mathbb R^r$ are drawn as $\beta_1',\beta_2' \overset{\text{i.i.d.}}{\sim} \mathcal N(0,I_r)$. To guarantee nonnegative hazards, we choose $ \lambda_0 = 0.05 - \min_{p\in\{0,1\}, n\in[N]} \big(\beta_1'^\top V_n + \beta_2'^\top U_p\big) $. Thus, the potential event times $\tau^{(0)}_{p,n}$ are sampled via inverse transform: \[ \tau^{(0)}_{p,n} = \frac{-\log W}{\lambda_0 + \beta_1'^\top V_n + \beta_2'^\top U_p}, \qquad W\sim\mathrm{Unif}(0,1). \]
\paragraph{Evaluation metrics.} For the treated unit $n$, let $\theta_n(t)=S^{(0)}_{1,n}(t)$ denote the true post–period control survival trajectory based on the true data generating process and $\widehat{\theta}_n(t)$ its SSC estimate, both evaluated on a common time grid $\mathcal T=\{t_1,\ldots,t_{T_0}\}\subset[0,\widetilde{\tau}]$. Specifically, the evaluation horizon is set to be $ \widetilde{\tau} :=\mathrm{quantile}_{0.90}\!\left(\{T_{p,n,i}\}\right), $ which trims extreme right tails where empirical support is sparse, stabilizing the downstream survival estimates. And the grid is then chosen as $T_0=100$ equally spaced timestamps $\mathcal T=\mathrm{linspace}(0,\widetilde{\tau},T_0)$. We quantify estimation accuracy using the sup–norm error \[ \|\widehat{\theta}_n-\theta_n\|_{\infty} :=\max_{t\in\mathcal T}\big|\widehat{\theta}_n(t)-\theta_n(t)\big|, \] which directly captures the worst–case deviation between the estimated and true trajectories.
\paragraph{Evaluation results.} Figure (ref) compares, in one example simulation with $K=700$ observations for each period-unit pair, the post-treatment control survival trajectory of the treated unit from the true data generating proces, SSC estimator and Confounder-Aware Parametric Estimator, which observes the true latent factors and fits the correctly specified parametric survival model. In both the Cox and Aalen settings, the SSC estimator closely matches the ground-truth trajectory. The agreement is particularly notable because SSC does not rely on any knowledge of the underlying survival model or the latent confounders. Its estimates also closely track those of the Confounder-Aware Parametric Estimator. This demonstrates that SSC can recover the correct counterfactual trajectory even when the data are generated from structurally distinct models and fully latent confounding variables that are agnostic to the SSC estimator itself.
We repeat the simulation process for 20 times for $K = \{100, 300, 700\}$ each respectively and Figure (ref) reports the distribution of sup-norm errors for SSC as the number of observations per unit-period increases. Both DGPs exhibit a clear and rapid reduction in error as $K$ grows: both the median error and variability shrink substantially from $K=100$ to $K=700$.
Table (ref) summarizes the mean and standard deviation of SSC’s sup-norm error across repeated simulations. In both DGPs, SSC achieves diminishing average errors as $K$ increases. The Confounder-Aware Parametric Estimator consistently attains lower error, as expected, since it is given access to the true confounders and the correct parametric form, but the performance gap narrows quickly with increasing $K$. This convergence in accuracy highlights SSC’s ability to fully exploit the latent low-rank structure without requiring model knowledge or observing the confounders themselves.
Overall, the synthetic experiments demonstrate that SSC is accurate, robust across model classes, and exhibits the expected improvement with additional data.
To complement the synthetic validation, we now demonstrate the use of SSC in our motivating clinical example for evaluating the second-line therapies for relapsed and refractor TCL. This application fits precisely into our SSC framework as all countries prescribe CC in both lines of treatment except the USA, as the treated unit, switches to prescribe the novel SA in second-line. Our multi-country dataset naturally supplies a donor pool of control units so that we can use our SSC estimation procedure to leverage the control units to construct a synthetic control for the USA under CC in the second-line and compare it against the factual USA receiving the novel SA treatment. This section therefore applies SSC to quantify the superiority of treatment efficacy of the novel SA in TCL and to assess estimator fidelity using a held-out USA-CC cohort.
We use a retrospective clinical dataset on 925 patients with relapsed and refractory TCL from 13 institutions across 10 countries. The data includes time-to-event information and treatment assignments (SA, as intervention, or CC, as control). Follow-up occurred from the date of cancer diagnosis until death, loss to follow-up, or end of study.
The analysis defines the first-line treatment period as pre-period ($p=0$) and the second-line treatment as post-period ($p=1$). The treatment assignment $D_{p,n}$ indicates whether a patient received SA ($D_{p,n}=1$) or CC ($D_{p,n}=0$) in period $p$. For each patient, the observed time $T_{p,n,i}$ is the minimum of the time-to-event ($\tau_{p,n,i}$) and time-to-censoring ($C_{p,n,i}$). The censoring indicator is $\Delta_{p,n,i}$ = 1 if $\tau_{p,n,i} \leq C_{p,n,i}$, and 0 otherwise. Since the target outcome is failure of treatment, in the pre-period, the event is death or start of second-line treatment, and censoring is loss-to-followup; in the post-period, the event is death or start of third-line treatment, and censoring is loss-to-followup.
Our empirical estimand is the counterfactual second-line survival trajectory for the USA had it continued prescribing CC rather than adopting SA. Applying the SSC procedure, we estimate the USA-SA cohort's counterfactual survival under CC in second-line and report the pre- and post-period survival curves in Figure (ref). The factual USA trajectory under SA (denoted in orange line) exhibits consistently higher survival than its synthetic control (denoted in blue dash line), suggesting superiority in efficacy of SA relative to CC.
To validate the estimator, we exploit the fact that we observe a cohort of USA patients who received CC in the second-line treatment. A necessary step is to establish that the held-out USA-CC and USA-SA cohorts do not differ systematically at baseline. We compare the two cohorts along a comprehensive set of prognostic, clinical, demographic, and laboratory variables, including International Prognostic Index (IPI), Prognostic Index for T-cell lymphoma (PIT), disease characteristics, and standard laboratory markers. Continuous covariates are assessed using two-sample $t$-tests and binary covariates using chi-squared tests. In addition, we compute standardized mean differences (SMDs), interpreting SMD $<0.1$ as negligible imbalance and SMD $<0.2$ as acceptable. Across all covariates, hypothesis tests yielded $p$-values $>0.1$, and all SMDs were below $0.2$ (with the majority below $0.1$), indicating that the two cohorts are well balanced at baseline. This comparability enables the held-out USA-CC cohort to serve as a real-world benchmark for assessing estimator fidelity. We also note that the USA-CC cohort is completely hidden from our SSC estimator, thus no information leakage. Figure (ref) shows that SSC closely matches the true USA--CC trajectory in the second-line, supporting the credibility of the synthetic control estimation.
To further assess robustness, we perform $500$ bootstrap resamples of the donor pool and recompute the synthetic USA-CC counterfactual for each resample. The resulting 95% bootstrap confidence interval is shown in Figure (ref) alongside the ground truth USA-CC trajectory. The synthetic estimates remain stable across resamples and the confidence band closely overlaps the held-out USA-CC curve throughout the post-period. This bootstrap analysis highlights the robustness of SSC to sampling variability and provides additional empirical support for the accuracy of our estimator in a real clinical setting.
This paper develops a causal framework for survival analysis in panel settings and applies it to evaluate impact of novel TCL therapies. Our method extends synthetic control techniques to the hazard domain, enabling counterfactual estimation with observational panel data.
Empirically, we find that access to novel SA therapies in the US is associated with a significant reduction in post-treatment hazard rates compared to a synthetic control constructed from countries that continue to rely on CC. This suggests that SA may provide a meaningful clinical benefit in prolonging survival among patients with aggressive TCL subtypes. Our framework also enables rigorous validation: by comparing synthetic USA under CC to actual USA-CC patient outcomes, we show that the synthetic control reproduces key survival trends in held-out data. This supports the internal validity of the estimated counterfactual survival.
Several limitations merit discussion. First, we assume non-informative censoring, which may not hold in all clinical settings. Extensions using IPCW ipcw or joint modeling should help relax this. Second, we do not directly leverage the observed covariates in the unobserved confounding panel setup. Future work could explore extensions to time-dependent treatment regimes and the incorporation of high-dimensional covariates. We expect our framework to be extended to applications in other domains where timing matters.