# Treatment Effect Estimation in Causal Survival Analysis: Practical Recommendations

Charlotte Voinot

Clément Berenfeld

Imke Mayer

Bernard Sebastien

Julie Josse

2025-12-08

The restricted mean survival time (RMST) difference offers an interpretable causal contrast to estimate the treatment effect for time-to-event outcomes, yet a wide range of available estimators leaves limited guidance for practice. We provide a unified review of RMST estimators for randomized trials and observational studies, establish identification and asymptotic properties, and supply new derivations where needed. Our extensive simulation study compares simple nonparametric methods (such as unweighted Kaplan-Meier estimators) alongside parametric and nonparametric implementations of the G-formula, weighting approaches, Buckley-James transformations, and augmented estimators under diverse censoring mechanisms and model specifications. Across scenarios, classical Kaplan-Meier estimators (weighted when required by the censoring process) and G-formula methods perform well in randomized settings, while in observational data G-formula estimators remain competitive; however, augmented estimators such as AIPTW-AIPCW generally offer robustness to model misspecification and a favorable bias-variance trade-off. Parametric estimators perform best under correct specification, whereas nonparametric methods avoid functional assumptions but require large sample sizes to achieve reliable performance. We offer practical recommendations for estimator choice and provide open-source R code to support reproducibility and application.

## 1 Introduction

Causal survival analysis is a rapidly evolving field that integrates causal inference (Rubin, 1974; Hernán MA, 2020) with survival analysis (Kalbfleisch and Prentice, 2002) to evaluate the impact of treatments on time-to-event outcomes, while accounting for censoring, i.e. situations where only partial information about an event's occurrence is available. In this literature, the Restricted Mean Survival Time (RMST) difference has emerged as a promising causal measure that offers an intuitive interpretation of the average survival time over a specified time horizon.
Despite its recent emergence, the field has expanded quickly, though often in a fragmented manner. As a result, researchers are faced with a broad spectrum of estimators, from weighting and regression-based approaches to more complex hybrid strategies, and must navigate numerous methodological decisions. These include, for example, the choice of how to estimate nuisance components, whether through parametric models, or nonparametric techniques, frequently without clear guidance. The situation is further complicated by the limited availability of accessible software implementations. Addressing these challenges is especially crucial in real-world contexts where randomized trials are not feasible and carefully designed observational studies, grounded in explicit causal assumptions, represent the primary source of evidence on treatment effects (Makary and Prasad, 2025). These considerations underscore the need for a structured review that consolidates existing approaches, clarifies their underlying assumptions, and evaluates estimator performance to support informed methodological decisions, particularly in finite-sample settings. Such a survey could also help to identify key methodological gaps that warrant further research.

Related efforts have already been made to clarify methodological options in survival analysis. For instance, Denz et al. (2023) focus on confounders-adjusted survival curves under proportional hazards and assess the impact of model misspecification. Yet, most existing reviews remain confined to the Cox model framework and thus overlook the broader range of estimators and causal estimands now available.

Despite its widespread use, we chose not to cover the estimation of the hazard ratio (HR), typically derived from a Cox model, because of several well-documented limitations. It relies on the proportional hazards assumption, which may not hold in practice. Moreover, the HR does not represent a causal contrast (Martinussen et al., 2020), even in randomized trials, due to inherent selection bias that can distort group comparisons unless the treatment has no effect (Martinussen and Vansteelandt, 2013; Hernán, 2010). The HR is also non-collapsible, meaning marginal effects cannot be inferred from conditional ones (Huitfeldt et al., 2019; Greenland et al., 1999). Finally, its clinical interpretation could be difficult, as it reflects an instantaneous risk rather than a cumulative or absolute survival benefit (Averbuch et al., 2025).

In this paper, after introducing the notations and the definition of the restricted mean survival time (Section 1.1), we present the statistical framework required to define and estimate RMST in both randomized controlled trials (Section 2) and observational settings (Section 3). We characterize their key statistical properties, such as consistency and asymptotic normality, and provide a unified, self-contained set of proofs. Where existing results are incomplete or unavailable, we supply original derivations; where they exist, we rederive them within a common formalism that enhances clarity and facilitates direct comparison between estimators. In Section 5, we complement this theoretical analysis with a simulation study assessing estimator performance under various conditions, including independent and conditionally independent censoring and both correctly specified and misspecified models. Finally, in Section 6, we offer practical recommendations for estimator selection, grounded in considerations such as finite-sample behavior and computational complexity.
1.1 Notations and definition of the Restricted Mean Survival Time

We set the analysis in the potential outcome framework, where a patient, described by a vector of covariates $X \in \mathbb{R}^p$, either receives a treatment ($A = 1$) or is in the control group ($A = 0$). The patient's survival time would then be $T^{(0)} \in \mathbb{R}^+$ under control and $T^{(1)} \in \mathbb{R}^+$ under treatment. In practice, both potential outcomes cannot be observed simultaneously for the same individual, since each patient receives only one treatment condition. Instead, we observe the realized survival time $T$, defined as follows:

**Assumption 1.1.** (Stable Unit Treatment Value Assumption — SUTVA) $T = AT^{(1)} + (1-A)T^{(0)}$.

Because of potential censoring, the outcome $T$ is not always fully observed. The most common form is right-censoring, which occurs when the event has not occurred by the end of follow-up (Turkson et al., 2021). In what follows, we focus on this case: for each subject, we observe an $n$-sample $(X, A, \tilde{T}, \Delta)$, where $\tilde{T} = \min(T, C)$, $C \in \mathbb{R}_+$, and $\Delta = \mathbb{I}\{T \le C\}$ (equal to 1 if the event time is observed and 0 otherwise). These observed variables are derived from the unobservable tuple $(X, A, T^{(0)}, T^{(1)}, C)$ illustrated in Table 1.

<table><caption>Table 1: Example of a survival dataset. Each observation corresponds to a patient. In practice, only $X$, $A$, $\Delta$, and $\tilde{T}$ are observed. The column Obs. Outcome represents the observed survival outcome.</caption><thead><tr><th colspan="3">Covariates</th><th>Treatment</th><th colspan="2">Potential outcome</th><th>Outcome</th><th>Censoring</th><th>Status</th><th>Obs. Outcome</th></tr><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>A</th><th>T<sup>(0)</sup></th><th>T<sup>(1)</sup></th><th>T</th><th>C</th><th>Δ</th><th>T&#773;</th></tr></thead><tbody><tr><td>1</td><td>1.5</td><td>4</td><td>1</td><td>?</td><td>200</td><td>200</td><td>?</td><td>1</td><td>200</td></tr><tr><td>5</td><td>1</td><td>2</td><td>0</td><td>100</td><td>?</td><td>100</td><td>?</td><td>1</td><td>100</td></tr><tr><td>9</td><td>0.5</td><td>3</td><td>1</td><td>?</td><td>?</td><td>?</td><td>200</td><td>0</td><td>200</td></tr></tbody></table>

Our aim is to estimate the Average Treatment Effect (ATE) defined as the difference between the Restricted Mean Survival Time (Royston and Parmar, 2013).

**Definition 1.1.** (Causal effect: Difference between Restricted Mean Survival Time) We define $\theta_{\text{RMST}} = \mathbb{E}[T^{(1)} \wedge \tau - T^{(0)} \wedge \tau]$, where $a \wedge b := \min(a, b)$ for $a, b \in \mathbb{R}$, and $\tau > 0$ is a pre-specified time horizon over which survival is evaluated.

We define the survival functions $S^{(a)}(t) := \mathbb{P}(T^{(a)} > t)$ for $a \in \{0, 1\}$, i.e., the probability that an individual will survive beyond a given time $t$ if treated or if not treated. Likewise, we let $S(t) := \mathbb{P}(T > t)$, and $S_C(t) := \mathbb{P}(C > t)$. We also let $G(t) := \mathbb{P}(C \ge t)$ be the left-limit of the survival function $S_C$. Because $T^{(a)} \wedge \tau$ are non-negative random variables, one can easily express the restricted mean survival time using the survival functions:

$$ \mathbb{E}(T^{(a)} \wedge \tau) = \int_{0}^{\infty} \mathbb{P}(T^{(a)} \wedge \tau > t)dt = \int_{0}^{\tau} S^{(a)}(t)dt. \quad (1) $$
Consequently, $\theta_{\text{RMST}}$ can be interpreted as the mean difference between the survival function of treated and control until a fixed time horizon $\tau$ (See Figure 1). A RMST difference $\theta_{\text{RMST}} = 10$ days with $\tau = 200$ means that, on average, the treatment increases the survival time by 10 days at 200 days.

Although the present work focuses on the estimation of the difference of RMST, we would like to stress that the causal effect can be assessed through other measures, such as for instance the difference between the survival functions $\theta_{\text{SP}} := S^{(1)}(\tau) - S^{(0)}(\tau)$ for some time $\tau$ (Ozenne et al., 2020).

<figure><img src="image_1.png" /><figcaption>Figure 1: Estimated survival curves on synthetic data. The $\theta_{\text{RMST}}$ at $\tau = 50$ corresponds to the yellow area between them. Curves are estimated with the Kaplan-Meier estimator, see Section 2.1.</figcaption></figure>

We provide in Table 2 a summary of the notations used throughout the paper.

<table><caption>Table 2: Summary of the notations.</caption><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>X</td><td>Covariates</td></tr><tr><td>A</td><td>Treatment indicator (A = 1 for treatment, A = 0 for control)</td></tr><tr><td>T</td><td>Survival time</td></tr><tr><td>T<sup>(a)</sup>, a ∈ {0, 1}</td><td>Potential survival time respectively with and without treatment</td></tr><tr><td>S<sup>(a)</sup>, a ∈ {0, 1}</td><td>Survival function S<sup>(a)</sup>(t) = P(T<sup>(a)</sup> > t) of the potential survival times</td></tr><tr><td>λ<sup>(a)</sup>, a ∈ {0, 1}</td><td>Hazard function λ<sup>(a)</sup>(t) = lim<sub>h→0+</sub> P(T<sup>(a)</sup> ∈ [t, t + h) | T<sup>(a)</sup> ≥ t)/h of the potential survival times</td></tr><tr><td>Λ<sup>(a)</sup>, a ∈ {0, 1}</td><td>Cumulative hazard function Λ<sup>(a)</sup>(t) := ∫<sub>0</sub><sup>t</sup> λ<sup>(a)</sup>(s) ds of the potential survival times</td></tr><tr><td>C</td><td>Censoring time</td></tr><tr><td>S<sub>C</sub></td><td>Survival function S<sub>C</sub>(t) = P(C > t) of the censoring time</td></tr><tr><td>G</td><td>Left-limit of the survival function G(t) = P(C ≥ t) of the censoring time</td></tr><tr><td>T̃</td><td>Observed time (T̃ ∧ C = min(T, C))</td></tr><tr><td>Δ</td><td>Event status I{T̃ ≤ C}</td></tr></tbody></table>
<table><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>Δ<sup>τ</sup></td><td>Event status of the restricted time Δ<sup>τ</sup> = max{Δ, I{T̃ ≥ τ}}</td></tr><tr><td>{t<sub>1</sub>, t<sub>2</sub>, ..., t<sub>K</sub>}</td><td>Discrete times</td></tr><tr><td>e(x)</td><td>Propensity score E[A|X = x]</td></tr><tr><td>μ(x, a), a ∈ {0, 1}</td><td>E[T ∧ τ | X = x, A = a]</td></tr><tr><td>S(t|x, a), a ∈ {0, 1}</td><td>Conditional survival function, P(T > t|X = x, A = a).</td></tr><tr><td>λ<sup>(a)</sup>(t|x), a ∈ {0, 1}</td><td>Conditional hazard functions of the potential survival times</td></tr><tr><td>G(t|x, a), a ∈ {0, 1}</td><td>left-limit of the conditional survival function of the censoring<br>P(C ≥ t|X = x, A = a)</td></tr><tr><td>Q<sub>S</sub>(t|x, a), a ∈ {0, 1}</td><td>E[T ∧ τ | X = x, A = a, T ∧ τ > t]</td></tr></tbody></table>

## 2 Causal survival analysis in Randomized Controlled Trials

Randomized Controlled Trials (RCTs) are the gold standard for establishing the effect of a treatment on an outcome, because treatment allocation is controlled through randomization, which ensures (asymptotically) the balance of covariates between treated and controls, and thus avoids problems of confounding between treatment groups. In a classical RCT, treatment is randomly assigned (Rubin, 1974), and everyone could, in principle, be assigned to either arm.

**Assumption 2.1. (Random Treatment Assignment)** $A \perp T^{(0)}, T^{(1)}, X$.

**Assumption 2.2. (Trial Positivity)** $0 < \mathbb{P}(A = 1) < 1$.

Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment) and 2.2 (Trial Positivity), classical causal identifiability equations can be written to express $\theta_{\text{RMST}}$ without potential outcomes.

$$ 
\begin{align}
\theta_{\text{RMST}} &= \mathbb{E}[T^{(1)} \wedge \tau - T^{(0)} \wedge \tau] \nonumber \\
&= \mathbb{E}[T^{(1)} \wedge \tau | A = 1] - \mathbb{E}[T^{(0)} \wedge \tau | A = 0] \quad (\text{Random treatment assignment}) \tag{2} \\
&= \mathbb{E}[T \wedge \tau | A = 1] - \mathbb{E}[T \wedge \tau | A = 0]. \quad (\text{SUTVA})
\end{align}
$$

However, Equation 2 still depends on $T$, which is only partially observed due to censoring. To ensure identifiability of treatment effects, standard assumptions on the censoring mechanism are required, namely independent and conditionally independent censoring.

When censoring is independent, standard survival estimators, such as the Kaplan–Meier estimator, directly accommodate censoring, allowing RMST to be obtained from the estimated survival curve via Equation 1 (Section 2.1). Under the more general condition of conditionally independent censoring, censoring-unbiased transformations of the event time can be used to recover $\mathbb{E}[T \wedge \tau | A = a]$ (Section 2.2). These transformations enable two estimation strategies: (i) direct estimation, by computing the transformed outcome and averaging it and (ii) adjusted survival-curve estimation, in which the transformation is incorporated within
a survival estimator, and RMST is then derived using Equation 1. These two families of methods form the basis of the estimators reviewed in the next sections.

## 2.1 Independent censoring: the Kaplan-Meier estimator

**Assumption 2.3. (Independent Censoring)** $C \perp T^{(0)}, T^{(1)}, X, A$.

Under Assumption 2.3, subjects censored at time $t$ are representative of all subjects who remain at risk at time $t$. Figure 2 shows the causal graph where the treatment assignment is random, the time $T$ depends on $(A, X)$, and the censoring time $C$ is independent of all variables.

<figure><img src="image_2.png" /><figcaption>Figure 2: Causal graph in RCT survival data with independent censoring.</figcaption></figure>

We assume that there is no almost-sure upper bound on the censoring time before $\tau$.

**Assumption 2.4. (Positivity of the Censoring Process)** There exists $\varepsilon > 0$ such that $G(\tau) \ge \varepsilon$.

If indeed it was the case that $\mathbb{P}(C < t) = 1$ for some $t < \tau$, then we would not be able to infer anything on the survival function on the interval $[t, \tau]$ as all observation times $\tilde{T}_i$ would be in $[0, t]$ almost surely. In practice, adjusting the threshold time $\tau$ can help satisfy this positivity assumption.

Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.3 (Independent Censoring), and 2.4 (Positivity of the Censoring Process), the potential-outcome survival $S^{(a)}(t)$ is identified as a functional of the joint law of $(\tilde{T}, \Delta, A)$. This enables several estimation strategies, the most well-known being the Kaplan-Meier, the Maximum Likelihood Estimator (MLE) of the survival functions, see for instance Kaplan and Meier (1958).

**Definition 2.1. (Kaplan-Meier estimator, (Kaplan and Meier, 1958))**

$$ \hat{S}_{\text{KM}}(t|A=a) := \prod_{t_k \le t} \left(1 - \frac{D_k(a)}{N_k(a)}\right). $$
with $D_k(a) := \sum_{i=1}^n \mathbb{I}(\tilde{T}_i = t_k, \Delta_i = 1, A = a)$ and $N_k(a) := \sum_{i=1}^n \mathbb{I}(\tilde{T}_i \ge t_k, A = a)$, which correspond respectively to the number of deaths $D_k(a)$ and of individuals at risk $N_k(a)$ at time $t_k$ in the treated group ($a = 1$) or in the control group ($a = 0$).

Furthermore, since $D_k(a)$ and $N_k(a)$ are sums of i.i.d. random variables, the Kaplan-Meier estimator inherits convenient statistical properties. In Appendix 6, we show that it exhibits exponentially fast bias decay (Proposition 6.1) and asymptotic normality with explicitly characterized variances (Proposition 6.2). Finally, the estimator in Definition 2.1 extends naturally to the continuous setting by replacing the fixed grid $\{t_k\}$ with observed event times $(\tilde{T}_i = t_k, \Delta_i = 1)$. We refer to Breslow (1974) for a study of this version and to (Aalen et al., 2008, Sec. 3.2) for a broader point-process perspective.

The associated RMST estimator is then defined as

$$ \hat{\theta}_{\text{KM}} = \int_{0}^{\tau} \hat{S}_{\text{KM}}(t|A=1) - \hat{S}_{\text{KM}}(t|A=0)dt. \quad (3) $$

This estimator inherits KM properties. It is consistent for the target RMST and its bias decays exponentially fast. In RCTs, the independence of arms delivers asymptotic normality with variance equal to the sum of the arm-specific RMST variances.

## 2.2 Conditionally independent censoring

Conditionally independent censoring relaxes the assumption of independent censoring by allowing censoring to depend on treatment or covariates — for instance, when treated patients withdraw from a study due to side effects.

**Assumption 2.5. (Conditionally Independent Censoring)** $C \perp T^{(0)}, T^{(1)} | X, A$

Under Assumption 2.5, subjects within a same stratum defined by $X = x$ and $A = a$ have equal probability of censoring at time $t$, for all $t$. In case of conditionally independent censoring, we also need to assume that all subjects have a positive probability to remain uncensored at their time-to-event.

**Assumption 2.6. (Conditional Positivity / Overlap for Censoring)** There exists $\varepsilon > 0$ such that for all $t \in [0, \tau]$, it almost surely holds $G(t|A, X) \ge \varepsilon$.

Figure 3 represents the causal graph when treatment assignment is randomized and when there is a conditionally independent censoring mechanism.

Under dependent censoring, the Kaplan-Meier estimator as defined in Definition 2.1 is not suited to estimate survival. Alternatives include plug-in G-formula estimators (Section 2.2.1) and unbiased transformations (Section 2.2.2).
<figure><img src="image_3.png" /><figcaption>Figure 3: Causal graph in RCT survival data with dependent censoring.</figcaption></figure>

### 2.2.1 The G-formula

Building on Equation 2, one can express the RMST conditioning on the covariates:

$$ \theta_{\text{RMST}} = \mathbb{E}[\mathbb{E}[T \wedge \tau | X, A = 1]] - \mathbb{E}[T \wedge \tau | X, A = 0]] = \mathbb{E}[\mu(X, 1) - \mu(X, 0)]. $$

with $\mu(x, a) := \mathbb{E}[T \wedge \tau | X = x, A = a]$.

An estimator of the RMST difference via the G-formula arises by plugging in any estimator $\hat{\mu}$ of $\mu$:

$$ \hat{\theta}_{\text{G-formula}} = \frac{1}{n} \sum_{i=1}^{n} \hat{\mu}(X_i, 1) - \hat{\mu}(X_i, 0). \quad (4) $$

The G-formula accommodates right censoring under conditionally independent censoring, without specifying a separate censoring model: the target (e.g., $S^{(a)}(t)$ or RMST$^{(a)}$) is expressed as a functional of the conditional outcome regression $\mu(\cdot)$. The same G-formula can be applied to observational settings (in Section 3.2) to address non-random treatment assignment. Different strategies can be used to estimate $\mu$, ranging from semi-parametric models such as the Cox model to fully non-parametric approaches like survival forests. For any chosen model, implementation follows one of two paradigms: T-learner approaches (i.e., two separate learners), where $\mu(\cdot, 1)$ and $\mu(\cdot, 0)$ are estimated separately; and S-learner approaches (i.e., a single learner), where $\hat{\mu}$ is obtained by fitting one model using covariates $(X, A)$.

**Cox Model.** The Cox proportional hazards model (Cox, 1972) relies on a semi-parametric model of the conditional hazard functions and specifies the function under treatment $a$, $\lambda^{(a)}(t|X)$, as

$$ \lambda^{(a)}(t|X) = \lambda_{0}^{(a)}(t) \exp(X^{\top}\beta^{(a)}), $$

where $\lambda_0^{(a)}$ is a baseline hazard function under treatment $a$ and $\beta^{(a)}$ has the same dimension as the vector of covariate $X$. The model is semi-parametric in the sense that the baseline hazard is left unspecified, while the effect of covariates enters through a multiplicative proportional hazards structure. Alternative semi-parametric models include Accelerated Failure Time (AFT) (Saikia and Barman, 2017) which directly models the survival time and does not rely on the proportional hazards (PH) assumption.
**Weibull Model.** A popular parametric model for survival is the Weibull Model, which amounts to assume that $\lambda^{(a)}(t|X) = t^{\gamma(a)} \exp(X^\top \beta)$ where $\gamma(a) > 0$. We refer to Zhang (2016) for a study of this model.

**Survival Forests.** On the nonparametric side, survival forests (Ishwaran et al., 2008) aggregate survival trees that use log-rank splitting and terminal-node Nelson-Aalen estimates to learn $S(t | X)$ under right censoring, capturing nonlinearities and high-order interactions without parametric form.

### 2.2.2 Censoring-unbiased transformations (CUT)

Censoring-unbiased transformations involve applying a transformation to $T$. Specifically, we compute a new time $T^*$ of the form

$$T^* := T^*(\tilde{T}, X, A, \Delta) = \begin{cases} \phi_0(\tilde{T} \wedge \tau, X, A) & \text{if } \Delta^\tau = 0, \\ \phi_1(\tilde{T} \wedge \tau, X, A) & \text{if } \Delta^\tau = 1. \end{cases} \qquad (5)$$

for two wisely chosen transformations $\phi_0$ and $\phi_1$, and where $\Delta^\tau := \mathbb{I}\{T \wedge \tau \le C\} = \Delta + (1 - \Delta)\mathbb{I}(\tilde{T} \ge \tau)$ indicates whether an individual is either uncensored or censored after time $\tau$. Since we only need the expectation of survival time thresholded at $\tau$, any censoring beyond $\tau$ can be treated as uncensored ($\Delta^\tau = 1$). A censoring-unbiased transformation $T^*$ is defined such that, for any $a \in \{0, 1\}$,

$$\mathbb{E}[T^*|A=a, X] = \mathbb{E}[T^{(a)} \wedge \tau|X] \quad \text{almost surely.} \qquad (6)$$

A notable advantage of this approach is that it enables the use of the full transformed dataset $(X_i, A_i, T_i^*)$ as if no censoring occurred. Moreover, by the reweighting identity

$$\mathbb{E}[\mathbb{E}[T^*|A=a, X]] = \mathbb{E}\left[\frac{\mathbb{I}\{A=a\}}{\mathbb{P}(A=a)} T^*\right], \qquad (7)$$

this equality links the arm-specific mean of $T^*$ to its full-sample inverse-probability-weighted representation. The RMST difference can be estimated by

$$\hat{\theta} = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{A_i}{\pi} - \frac{1 - A_i}{1 - \pi} \right) T_i^*, \qquad (8)$$

where $\pi = \mathbb{P}(A = 1) \in (0, 1)$ by Assumption 2.2 and where $T_i^* = T^*(\tilde{T}_i, X_i, A_i, \Delta_i)$.

The following proposition establishes the theoretical guarantees for the class of estimators constructed from bounded censoring-unbiased transformations.

**Proposition 2.1** (Consistency of CUT as RMST estimator with known $\pi$). Under Assumptions 2.1 (Random Treatment Assignment) and 2.2 (Trial Positivity), the estimator $\hat{\theta}$ derived as in Equation 8 from square integrable censoring-unbiased transformations satisfying Equation 6 is an unbiased, strongly consistent, and asymptotically normal estimator of the RMST difference.
Square integrability will be ensured any time the transformation is bounded, which will always be the case of the ones considered in this work.

It is natural in a RCT setting to assume that probability of being treated $\pi$ is known. If not, it is usual to replace $\pi$ by its empirical counterpart $\hat{\pi} = n_1/n$ where $n_a = \sum_i 1\{A = a\}$. The resulting estimator takes the form

$$ \hat{\theta} = \frac{1}{n_1} \sum_{A_i=1} T_i^* - \frac{1}{n_0} \sum_{A_i=0} T_i^* \quad (9) $$

While this estimator introduces a slight bias due to the estimation of $\pi$ (see for instance Colnet et al. (2022), Lemma 2), it remains strongly consistent and asymptotically normal under the same assumptions:

**Proposition 2.2** (Consistency of CUT as RMST estimator with estimated $\hat{\pi}$). Under Assumptions 2.1 (Random Treatment Assignment) and 2.2 (Trial Positivity), the estimator $\hat{\theta}$ derived as in Equation 9 from a square integrable censoring-unbiased transformations satisfying Equation 6 is a strongly consistent, and asymptotically normal estimator of the RMST difference.

For the later, we consider only estimator with estimated probability of being treated in Equation 9.

The two most popular transformations are Inverse-Probability-of-Censoring Weighting (Koul et al., 1981) and Buckley-James (Buckley and James, 1979), both illustrated in Figure 4 and detailed below. In the former, only non-censored observations are considered and they are

<figure><img src="image_4.png" /><figcaption>Figure 4: Illustration of Inverse-Probability-of-Censoring and Buckley-James transformations.</figcaption></figure>

weighted while in the latter, censored observations are imputed with an estimated survival
time.

## The Inverse-Probability-of-Censoring Weighted transformation

The Inverse-Probability-of-Censoring Weighted (IPCW) transformation, introduced by Koul et al. (1981) in the context of censored linear regression, involves discarding censored observations and applying weights to uncensored data (see Figure 4). More precisely, we let

$$T_{\text{IPCW}}^* = \frac{\Delta^\tau}{G(\tilde{T} \wedge \tau | X, A)} \tilde{T} \wedge \tau, \quad (10)$$

where we recall that $G(t|X, A) := \mathbb{P}(C \ge t|X, A)$ is the left limit of the conditional survival function of the censoring. This estimator assigns higher weights to uncensored subjects within a covariate group that is highly prone to censoring, thereby correcting for conditionally independent censoring and removing selection bias (Howe et al., 2016).

Using proposition 2.1 and 2.2, we show in Appendix that the IPCW transformation is a censoring-unbiased transformation in the sense of Equation 6 (Proposition 6.3).

The IPCW depends on the unknown conditional survival function of the censoring $G(\cdot|X, A)$, which thus needs to be estimated. Estimating conditional censoring or the conditional survival function can be approached similarly, as both involve estimating a time, whether for survival or censoring. Consequently, we can use parametric methods, or semi-parametric methods, such as the Cox model, or non-parametric approaches like survival forests (see Section 2.2.1). Once $\hat{G}(\cdot|A, X)$ is provided, one can construct an estimator of the difference of RMST based on Equation 9,

$$\hat{\theta}_{\text{IPCW}} = \frac{1}{n_1} \sum_{A_i=1} T_{\text{IPCW},i}^* - \frac{1}{n_0} \sum_{A_i=0} T_{\text{IPCW},i}^* \quad (11)$$

Alternatively, instead of rescaling by $n_a$, the estimator can be stabilized by dividing by the sum of the weights, yielding a Hájek-type normalization that reduces variance without affecting consistency (Hájek, 1971; Robins et al., 2000; Hernán MA, 2020).

By Proposition 2.1 (Consistency of CUT as RMST estimator with known $\pi$), Proposition 2.2 (Consistency of CUT as RMST estimator with estimated $\hat{\pi}$) and Proposition 6.3 (IPCW is a Censoring Unbiased Transformation) (in Appendix 6), we easily deduce that $\hat{\theta}_{\text{IPCW}}$ is asymptotically consistent as soon as $\hat{G}$ is.

**Corollary 2.1.** Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), if $\hat{G}$ is uniformly weakly (resp. strongly) consistent then so is $\hat{\theta}_{\text{IPCW}}$, defined in Equation 11.

This result simply comes from the fact that $\hat{\theta}_{\text{IPCW}}$ depends continuously on $\hat{G}$ and that $G$ is lower-bounded (Assumption 2.6 Conditional Positivity / Overlap for Censoring). Surprisingly, we found limited use of this estimator in the literature, beside its first introduction in Koul et al. (1981).
A more common variant is the IPCW-Kaplan-Meier (Robins and Finkelstein, 2000; Satten and Datta, 2001; Willems et al., 2018), defined as follows. In this approach, censoring weights are dynamically evaluated at each distinct event time, allowing the estimator to update risk sets and weighted failures throughout the follow-up period.

**Definition 2.2.** (IPCW-Kaplan-Meier) We let again $\hat{G}(\cdot|X, A)$ be an estimator of the (left limit of) the conditional censoring survival function and we introduce

$$D_k^{\text{IPCW}}(a) := \sum_{i=1}^{n} \frac{\Delta_i^\tau}{\hat{G}(t_k|X_i, A=a)} \mathbb{I}(\tilde{T}_i = t_k, A_i = a)$$

$$\text{and} \quad N_k^{\text{IPCW}}(a) := \sum_{i=1}^{n} \frac{1}{\hat{G}(t_k|X_i, A=a)} \mathbb{I}(\tilde{T}_i \ge t_k, A_i = a),$$

be the weight-corrected numbers of deaths and of individuals at risk at time $t_k$. The IPCW version of the KM estimator is defined as:

$$\hat{S}_{\text{IPCW}}(t|A=a) = \prod_{t_k \le t} \left( 1 - \frac{D_k^{\text{IPCW}}(a)}{N_k^{\text{IPCW}}(a)} \right).$$

Note that the quantity $\pi$ is not present in the definition of $D_k^{\text{IPCW}}(a)$ and $N_k^{\text{IPCW}}(a)$ because it would simply disappear in the ratio $D_k^{\text{IPCW}}(a)/N_k^{\text{IPCW}}(a)$. The subsequent RMST estimator then takes the form

$$\hat{\theta}_{\text{IPCW-KM}} = \int_{0}^{\tau} \hat{S}_{\text{IPCW}}(t|A=1) - \hat{S}_{\text{IPCW}}(t|A=0)dt. \quad (12)$$

Like before for the classical KM estimator, this new reweighted KM estimator enjoys good statistical properties.

We show in Proposition 6.4 (Appendix 6) that the oracle version of this estimator is strongly consistent and asymptotically normal. These properties also hold for the estimator using a consistently estimated $\hat{G}$.

**Corollary 2.2.** Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), if $\hat{G}$ is uniformly weakly (resp. strongly) consistent then so is $\hat{S}_{\text{IPCW}}(t|A=a)$.

Because the evaluation of $N_k^{\text{IPCW}}(a)$ now depends on information gathered after time $t_k$, the previous derivations of bias (Appendix 6, proposition 6.1) and asymptotic variance (Appendix 6, proposition 6.2) no longer hold. Overgaard (2025) derived asymptotic variance expressions for related Kaplan-Meier-based inverse probability-of-censoring weighted estimators using influence-function representations. While these results theoretically extend to our IPCW-Kaplan-Meier formulation, the resulting expressions remain difficult to interpret or implement in practice. Based on our own results, we can nevertheless establish that
when G is unknown, consistency still holds provided that the estimator $\hat{G}(\cdot|A, X)$ is itself consistent.

## The Buckley-James transformation

The Buckley-James (BJ) transformation, introduced by Buckley and James (1979), takes a different path by leaving all uncensored values untouched, while replacing the censored ones by an extrapolated value (see Figure 4). Formally, it is defined as follows:

$$T_{\text{BJ}}^* = \Delta^\tau(\tilde{T} \wedge \tau) + (1 - \Delta^\tau)Q_S(\tilde{T} \wedge \tau|X, A), \quad (13)$$

where, for $t \le \tau$,

$$Q_S(t|X, A) := \mathbb{E}[T \wedge \tau|X, A, T \wedge \tau > t] = \frac{1}{S(t|X, A)} \int_t^\tau S(u|X, A) du$$

and $S(t|X, A = a) := \mathbb{P}(T^{(a)} > t|X)$ is the conditional survival function.

This transformation satisfies the censoring-unbiasedness condition defined in Equation 6 (Proposition 6.5 in Appendix 6). The BJ transformation depends on a nuisance parameter ($Q_S(\cdot|X, A)$) that needs to be estimated. Estimation strategies described in Section 2.2.1 can be used. Once provided with an estimator $\hat{Q}_S(\cdot|A, X)$, a natural estimator of the RMST based on the BJ transformation and Equation 9 is

$$\hat{\theta}_{\text{BJ}} = \frac{1}{n_1} \sum_{A_i=1} T_{\text{BJ},i}^* - \frac{1}{n_0} \sum_{A_i=0} T_{\text{BJ},i}^*. \quad (14)$$

Like for the IPCW transformation, the BJ transformation yields a consistent estimate of the RMST as soon as the model is well-specified.

**Corollary 2.3.** Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), if $\hat{Q}_S$ is uniformly weakly (resp. strongly) consistent then so is $\hat{\theta}_{\text{BJ}}$ defined as in Equation 14.

The proof is an application of Propositions 2.1 (Consistency of CUT as RMST estimator with known $\pi$), 2.2 (Consistency of CUT as RMST estimator with estimated $\hat{\pi}$) and 6.5 (BJ is a censoring unbiased transformation in Appendix, Section 6), and relies on the continuity of $S \mapsto Q_S$. Additionally, the Buckley-James transformation enjoys a notable optimality property.

**Theorem 2.1** (BJ minimizes the mean squared error among CUT). For any transformation $T^*$ of the form 5, it holds

$$\mathbb{E}[(T_{\text{BJ}}^* - T \wedge \tau)^2] \le \mathbb{E}[(T^* - T \wedge \tau)^2].$$

This result, originally stated without proof in Fan and Gijbels (1994), is proved there for completeness in Appendix, Section 6.
### 2.2.3 Augmented corrections

The main disadvantage of the two previous transformations is that they heavily rely on the specification of good estimator $\hat{G}$ (for IPCW) or $\hat{S}$ (for BJ). In order to circumvent this limitation, Rubin and van der Laan (2007) proposed the following transformations, whose expression is based on theory of semi-parametric estimation developed in Van der Laan and Robins (2003),

$$T_{\text{DR}}^* = \frac{\Delta^\tau \tilde{T} \wedge \tau}{G(\tilde{T} \wedge \tau | X, A)} + \frac{(1 - \Delta^\tau) Q_S(\tilde{T} \wedge \tau | X, A)}{G(\tilde{T} \wedge \tau | X, A)} - \int_0^{\tilde{T} \wedge \tau} \frac{Q_S(t | X, A)}{G(t | X, A)^2} d\mathbb{P}_C(t | X, A), \quad (15)$$

where $d\mathbb{P}_C(t|X, A)$ is the distribution of C conditional on the covariates X and A. As this transformation depends on both conditional survival functions G and S, we will denote it $T_{\text{DR}}^*(G, S)$ to make this dependence explicit. Importantly, this transformation not only satisfies the censoring-unbiasedness condition of Equation 6, but also enjoys a double robustness property: it remains valid as long as either G or S is correctly specified. The formal statement of this result is given in Proposition 6.6 in Appendix 6. The statement and proof of this results is found in Rubin and van der Laan (2007) in the case where C and T are continuous.

## 3 Causal survival analysis in observational studies

Unlike RCT, observational data (such as from registries, electronic health records, or national healthcare systems) are collected without controlled randomized treatment allocation. In such cases, treated and control groups are likely unbalanced due to the non-randomized design, which results in the treatment effect being confounded by variables influencing both the survival outcome T and the treatment allocation A. To enable identifiability of the causal effect, additional standard assumptions are needed.

**Assumption 3.1. (Conditional Exchangeability / Unconfoundedness) It holds**
$$A \perp T^{(0)}, T^{(1)} | X$$

Under Assumption 3.1, the treatment assignment is randomly assigned conditionally on the covariates X. This assumption states that there are no unmeasured confounders as the latter would make it impossible to distinguish correlation from causality.

**Assumption 3.2. (Positivity / Overlap for Treatment) Letting $e(X) := P(A = 1 | X)$**
be the propensity score, there holds $0 < e(X) < 1$ almost surely.

The Assumption 3.2 requires adequate overlap in covariate distributions between treatment groups, meaning every observation must have a non-zero probability of being treated and not being treated. Because Assumption 2.1 (Random Treatment Assignment) no longer
holds anymore, neither does the previous identifiability Equation 2. In this new context, we can write

$$ 
\begin{align} 
\theta_{\text{RMST}} &= \mathbb{E}[T^{(1)} \wedge \tau - T^{(0)} \wedge \tau] \nonumber \\
&= \mathbb{E}[\mathbb{E}[T^{(1)} \wedge \tau | X] - \mathbb{E}[T^{(0)} \wedge \tau | X]] \nonumber \\
&= \mathbb{E}[\mathbb{E}[T^{(1)} \wedge \tau | X, A=1] - \mathbb{E}[T^{(0)} \wedge \tau | X, A=0]] \tag{16} \\
&= \mathbb{E}[\mathbb{E}[T \wedge \tau | X, A=1] - \mathbb{E}[T \wedge \tau | X, A=0]]. \tag{SUTVA}
\end{align} 
$$

A direct application of the above equation could be a direct plug-in G-formula estimator developed in Section 2.2.1, equation 4. In particular, Chen and Tsiatis (2001) prove consistency and asymptotic normality results for Cox estimators in an observational study, and they give an explicit formulation of the asymptotic variance as a function of the parameters of the Cox model. In the non-parametric literature, Foster et al. (2011) and Künzel et al. (2019) empirically study this estimator using survival forests, with the former employing it as a T-learner (referred to as Virtual Twins) and the latter as an S-learner.

In another direction, the next sections describe alternative estimation strategies that also rely on equation 16 but adapt the estimation of the conditional expectations to the specific censoring assumptions—either independent or conditionally independent censoring.

## 3.1 Independent censoring

Under Assumptions 3.1 (Conditional Exchangeability / Unconfoundedness) and 2.3 (Independent Censoring), the data exhibit confounding effect induced by the confounders. The corresponding causal graph is shown in Figure 5.

<figure><img src="image_5.png" /><figcaption>Figure 5: Causal graph in observational survival data with independent censoring.</figcaption></figure>

Under Assumption 2.3 (Independent Censoring), we saw in Section 2.1 that the Kaplan-Meier
estimator could conveniently handle censoring. Building on Equation 16, we can write

$$
\begin{align}
\theta_{\text{RMST}} &= \mathbb{E}[\mathbb{E}[T \wedge \tau | X, A = 1] - \mathbb{E}[T \wedge \tau | X, A = 0]] \nonumber \\
&= \mathbb{E}\left[\frac{\mathbb{E}[\mathbb{I}\{A = 1\}T \wedge \tau | X]}{\mathbb{E}[\mathbb{I}\{A = 1\}|X]}\right] - \mathbb{E}\left[\frac{\mathbb{E}[\mathbb{I}\{A = 0\}T \wedge \tau | X]}{\mathbb{E}[\mathbb{I}\{A = 0\}|X]}\right] \quad (\text{Positivity for treatment}) \nonumber \\
&= \mathbb{E}\left[\frac{AT \wedge \tau}{e(X)}\right] - \mathbb{E}\left[\frac{(1-A)T \wedge \tau}{1-e(X)}\right], \tag{17}
\end{align}
$$

which suggests to adapt the classical KM estimator by reweighting it by the propensity score (see equation 1). The use of propensity score in causal inference has been initially introduced by Rosenbaum and Rubin (1983) and further developed in Hirano et al. (2003). It was extended to survival analysis by Xie and Liu (2005) through the adjusted Kaplan-Meier estimator as defined below.

**Definition 3.1.** (IPTW Kaplan-Meier estimator) We let $\hat{e}(\cdot)$ be an estimator of the propensity score $e(\cdot)$. We introduce

$$
\begin{align*}
D_k^{\text{IPTW}}(a) &:= \sum_{i=1}^{n} \left( \frac{A_i}{\hat{e}(X_i)} + \frac{1-A_i}{1-\hat{e}(X_i)} \right) \mathbb{I}(\tilde{T}_i = t_k, \Delta_i = 1, A_i = a) \\
\text{and } N_k^{\text{IPTW}}(a) &:= \sum_{i=1}^{n} \left( \frac{A_i}{\hat{e}(X_i)} + \frac{1-A_i}{1-\hat{e}(X_i)} \right) \mathbb{I}(\tilde{T}_i \ge t_k, A_i = a),
\end{align*}
$$

be the numbers of deaths and of individuals at risk at time $t_k$, reweighted by the propensity score. The Inverse-Probability-of-Treatment Weighting (IPTW) version of the KM estimator is defined as:

$$
\hat{S}_{\text{IPTW}}(t|A=a) = \prod_{t_k \le t} \left( 1 - \frac{D_k^{\text{IPTW}}(a)}{N_k^{\text{IPTW}}(a)} \right). \quad (18)
$$

We let $S_{\text{IPTW}}^*(t|A=a)$ be the oracle KM-estimator defined as above with $\hat{e}(\cdot) = e(\cdot)$. Although the reweighting slightly changes the analysis, the good properties of the classical KM carry on to this setting.

In Appendix 6, Proposition 6.7 establishes that $S_{\text{IPTW}}^*(t | A = a)$ is a strongly consistent and asymptotically normal estimator of the survival function $S^{(a)}(t)$. Corollary 3.1 further shows that the same properties hold for the practical estimator provided the estimated propensity score converges uniformly. This follows directly from $S_{\text{IPTW}}^*$ is a continuous function of $e(\cdot)$, and because $e$ and $1-e$ are lower-bounded as per Assumptions 3.2 (Positivity / Overlap for Treatment).

**Corollary 3.1.** Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.3 (Independent Censoring) and 2.4 (Positivity of the Censoring Process), if $\hat{e}(\cdot)$ satisfies $\|\hat{e} - e\|_{\infty} \to 0$ almost surely (resp. in probability), then the IPTW Kaplan-Meier estimator $\hat{S}_{\text{IPTW}}(t|A=a)$ is a strongly (resp. weakly) consistent estimator of $S^{(a)}(t)$.
The resulting RMST estimator simply takes the form:

$$ \hat{\theta}_{\text{IPTW-KM}} = \int_{0}^{\tau} \hat{S}_{\text{IPTW}}(t|A=1) - \hat{S}_{\text{IPTW}}(t|A=0)dt. \quad (19) $$

which will be consistent under the same Assumptions as the previous Corollary. Note that we are not aware of any formal results concerning the bias and the asymptotic variance of the oracle estimator $S_{\text{IPTW}}^*(t|A=a)$, and we refer to Xie and Liu (2005) for heuristics concerning these questions.

## 3.2 Conditional independent censoring

Under Assumptions 3.1 (Conditional Exchangeability / Unconfoundedness) and 2.5 (Conditionally Independent Censoring), the data exhibit both confounding and conditionally independent censoring. The corresponding causal graph is shown in Figure 6. To address this, weight the IPCW and Buckley-James (BJ) transformations from Section 2.2.2 by the propensity score.

<figure><img src="image_6.png" /><figcaption>Figure 6: Causal graph in observational survival data with dependent censoring.</figcaption></figure>

### 3.2.1 IPTW-IPCW transformations

One can check that the IPCW transformation as introduced in Equation 10 is also a censoring-unbiased transformation in observational study (see Appendix in Proposition 6.8). Deriving an estimator of the difference of RMST is however slightly different in that context. In particular, Equation 7 rewrites

$$ \mathbb{E}[\mathbb{E}[T^*|X, A=1]] = \mathbb{E}\left[\frac{A}{e(X)}T^*\right], $$

which suggests defining an estimator of the RMST difference that can, similarly to IPCW estimators, be stabilized by normalizing with the sum of the weights, and is identified by combining Equation 17 with the IPCW representation in Equation 10.

$$ \hat{\theta}_{\text{IPTW-IPCW}} = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{A_i}{e(X_i)} - \frac{1 - A_i}{1 - e(X_i)} \right) T_{\text{IPCW},i}^* \quad (20) $$
This transformation now depends on two nuisance parameters, namely the conditional survival function of the censoring (through $T_{\text{IPCW}}^*$) and the propensity score. Proposition 3.1 ensures consistency of this estimator when both nuisance quantities are consistently estimated.

**Proposition 3.1** (Consistency of $\hat{\theta}_{\text{IPTW-IPCW}}$ with estimated $\hat{e}$ and $\hat{G}$). Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), and if $\hat{G}(\cdot|X, A)$ and $\hat{e}(\cdot)$ are uniformly weakly (resp. strongly) consistent estimators, then estimator 20 defined with $\hat{e}$ and $\hat{G}$ is a weakly (resp. strongly) consistent estimator of the RMST difference.

We can also use the same strategy as for the IPCW transformation and incorporate the time-dependent weights into a Kaplan-Meier estimator.

**Definition 3.2.** (IPTW-IPCW-Kaplan-Meier) We let again $\hat{G}(\cdot|X, A)$ and $\hat{e}(\cdot)$ be estimators of the conditional censoring survival function and of the propensity score. We introduce

$$D_k^{\text{IPTW-IPCW}}(a) := \sum_{i=1}^{n} \left( \frac{A_i}{\hat{e}(X_i)} + \frac{1-A_i}{1-\hat{e}(X_i)} \right) \frac{\Delta_i^\tau}{\hat{G}(t_k|X_i, A=a)} \mathbb{I}(\tilde{T}_i = t_k, A_i = a) \\ \text{and } N_k^{\text{IPTW-IPCW}}(a) := \sum_{i=1}^{n} \left( \frac{A_i}{\hat{e}(X_i)} + \frac{1-A_i}{1-\hat{e}(X_i)} \right) \frac{1}{\hat{G}(t_k|X_i, A=a)} \mathbb{I}(\tilde{T}_i \ge t_k, A_i = a),$$

be the weight-corrected numbers of deaths and of individual at risk at time $t_k$. The IPTW-IPCW version of the KM estimator is defined as:

$$\hat{S}_{\text{IPTW-IPCW}}(t|A=a) = \prod_{t_k \le t} \left( 1 - \frac{D_k^{\text{IPTW-IPCW}}(a)}{N_k^{\text{IPTW-IPCW}}(a)} \right).$$

The RMST difference estimated with IPTW-IPCW-Kaplan-Meier survival curves is then simply as

$$\hat{\theta}_{\text{IPTW-IPCW-KM}} = \int_{0}^{\tau} \hat{S}_{\text{IPTW-IPCW}}(t|A=1) - \hat{S}_{\text{IPTW-IPCW}}(t|A=0)dt. \quad (21)$$

In Appendix 6, Proposition 6.9 shows that the oracle version of this estimator is consistent and asymptotically normal. The result extends to the empirical estimator: as shown in Corollary 3.2, uniform consistency of the nuisance estimators ensures the consistency of the estimated survival curves $\hat{\theta}_{\text{IPTW-IPCW-KM}}$.

**Corollary 3.2.** Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), and for all $t \in [0, \tau]$, if the nuisance estimators $\hat{G}(\cdot|A, X)$ and $\hat{e}$ are weakly (resp. strongly) uniformly consistent, then $\hat{S}_{\text{IPTW-IPCW}}(t|A=a)$ is a weakly (resp. strongly) consistent estimator of $S^{(a)}(t)$.
We are not aware of general formulas for the asymptotic variances in this context. We mention nonetheless that Schaubel and Wei (2011) have been able to derive asymptotic laws in this framework in the particular case of Cox-models.

### 3.2.2 IPTW-BJ transformations

As with IPCW, the BJ transformation in Equation 13 is censoring-unbiased in observational contexts (Proposition 6.10, Appendix 6). The corresponding estimator of the difference in RMST is

$$ \hat{\theta}_{\text{IPTW-BJ}} = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{A_i}{e(X_i)} - \frac{1 - A_i}{1 - e(X_i)} \right) T_{\text{BJ},i}^* \quad (22) $$

This transformation depends on the conditional survival function S (through $T_{\text{BJ}}^*$) and the propensity score. Proposition 3.2 stated that consistency of the nuisance parameter estimators implies again consistency of the RMST estimator.

**Proposition 3.2** (Consistency of $\hat{\theta}_{\text{IPTW-BJ}}$ with estimated $\hat{e}$ and $\hat{S}$). Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), and if $\hat{S}(\cdot|X, A)$ and $\hat{e}(\cdot)$ are uniformly weakly (resp. strongly) consistent estimators, then $\hat{\theta}_{\text{IPTW-BJ}}$ defined with $\hat{S}$ and $\hat{e}$ is a weakly (resp. strongly) consistent estimator of the RMST.

### 3.2.3 Double augmented corrections

Building on Equations 16 and 17, it follows that

$$ \mathbb{E} \left[ \frac{\mathbb{I}\{A=a\}}{\mathbb{E}[\mathbb{I}\{A=a\} | X]} \left( (T \wedge \tau) - \mathbb{E}[T \wedge \tau | X, A=a] \right) | X \right] = 0. $$

Thus,

$$ \theta_{\text{RMST}} = \mathbb{E} \left[ \begin{aligned} & \mathbb{E}[T \wedge \tau | X, A=1] - \mathbb{E}[T \wedge \tau | X, A=0] \\ & = \mathbb{E} \left[ \mathbb{E}[T \wedge \tau | X, A=1] - \mathbb{E}[T \wedge \tau | X, A=0] \right] \end{aligned} \right. \\ & \underbrace{+ \frac{\mathbb{I}\{A=1\}}{e(X)} (T \wedge \tau - \mathbb{E}[T \wedge \tau | X, A=1]) - \frac{\mathbb{I}\{A=0\}}{1-e(X)} (T \wedge \tau - \mathbb{E}[T \wedge \tau | X, A=0])}_{\text{add mean-zero terms}}. \quad (23) $$

This representation coincides with the classical doubly robust augmented inverse probability weighting (AIPTW) functional from causal inference (Robins et al., 1994). Incorporating
the doubly robust censoring transformations from Section 2.2.3 yields a transformation that is doubly robust with respect to both censoring and treatment assignment

$$ \Delta_{QR}^* = \Delta_{QR}^*(G, S, \mu, e) := \left( \frac{A}{e(X)} - \frac{1-A}{1-e(X)} \right) (T_{DR}^*(G, S) - \mu(X, A)) + \mu(X, 1) - \mu(X, 0), $$

where we recall that $T_{DR}^*$ is defined in Section 2.2.3. This transformation depends on four nuisance parameters: G and S through $T_{DR}^*$, and now the propensity score e and the conditional response $\mu$. This transformation does not really fall into the scope of censoring-unbiased transform, but it is easy to show that $\Delta_{QR}^*$ is doubly robust for censoring and treatment assignment. That is, among the four nuisance functions involved, only one from each pair (S, G) and (e, $\mu$) needs to be correctly specified. This property (similar to Ozenne et al. (2020), Thm 1) is formally stated in Proposition 6.11 in Appendix 6.

Based on estimators $(\hat{G}, \hat{S}, \hat{\mu}, \hat{e})$ of $(G, S, \mu, e)$, one can then propose the following estimator of the RMST, coined the AIPTW-AIPCW estimator in Ozenne et al. (2020):

$$ \begin{aligned} \hat{\theta}_{\text{AIPTW-AIPCW}} &:= \frac{1}{n} \sum_{i=1}^{n} \Delta_{\text{QR},i}^*(\hat{G}, \hat{S}, \hat{\mu}, \hat{e}) \\ &= \frac{1}{n} \sum_{i=1}^{n} \left( \frac{A_i}{\hat{e}(X_i)} - \frac{1-A_i}{1-\hat{e}(X_i)} \right) (T_{\text{DR}}^*(\hat{G}, \hat{S})_i - \hat{\mu}(X_i, A_i)) + \hat{\mu}(X_i, 1) - \hat{\mu}(X_i, 0). \end{aligned} \tag{24} $$

Asymptotic properties of this estimator, such as asymptotic normality and semiparametric efficiency, are detailed in Ozenne et al. (2020).

# 4 Summary of the estimators and variance estimation

For completeness, Table 3 offers a compact synopsis of the estimators, context of application, associated nuisance models, and large-sample behavior, with a focus on misspecification sensitivity.

<table><caption>Table 3: Estimators of the RMST difference and nuisance parameters needed to compute each estimator. Empty boxes indicate that the nuisance parameter is not needed in the estimator thus misspecification has no effect.</caption><thead><tr><th>Estimator</th><th>RCT</th><th>Obs</th><th>Ind Cens</th><th>Dep Cens</th><th>Outcome model</th><th>Censoring model</th><th>Treatment model</th><th>Robustness</th></tr></thead><tbody><tr><td>Unadjusted KM</td><td>✓</td><td></td><td>✓</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>IPCW-KM</td><td>✓</td><td></td><td>✓</td><td>✓</td><td></td><td>G</td><td></td><td></td></tr><tr><td>BJ</td><td>✓</td><td></td><td>✓</td><td>✓</td><td>S</td><td></td><td></td><td></td></tr><tr><td>IPTW-KM</td><td>✓</td><td>✓</td><td>✓</td><td></td><td></td><td></td><td>e</td><td></td></tr><tr><td>IPCW-IPTW-KM</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td></td><td>G</td><td>e</td><td></td></tr></tbody></table>
<table><thead><tr><th>Estimator</th><th>RCT</th><th>Obs</th><th>Ind Cens</th><th>Dep Cens</th><th>Outcome model</th><th>Censoring model</th><th>Treatment model</th><th>Robustness</th></tr></thead><tbody><tr><td>IPTW-BJ</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>S</td><td></td><td>e</td><td></td></tr><tr><td>G-formula</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>μ</td><td></td><td></td><td></td></tr><tr><td>AIPTW-AIPCW</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>S, μ</td><td>G</td><td>e</td><td>✓ (Prp 6.11)</td></tr></tbody></table>

Some survival estimators presented in this work have known asymptotic variance expressions, such as the Kaplan-Meier estimator (via Greenwood's formula, derived in Proposition 6.2) or AIPTW-AIPCW estimator (derived in Ozenne et al. (2020)). Under standard regularity conditions ensuring asymptotic normality and assuming independence between the two survival curves, the variance of the RMST difference can be readily obtained via the delta method.

However, in practice the two estimated survival curves are typically correlated, which complicates variance estimation. While explicit variance formulas have been derived in certain settings (such as for the G-formula combined with Cox estimators (see Chen and Tsiatis (2001))) these expressions often remain complex and difficult to implement. We therefore recommend the nonparametric bootstrap (Efron and Tibshirani, 1994) as a robust and practical solution for variance estimation and confidence interval construction in general settings. The effectiveness of bootstrap methods in the context of censored data has been investigated in several studies, including Efron (1981) and Chen and Lo (1996), which provide both theoretical justification and empirical validation for its use with Kaplan-Meier estimator.

# 5 Simulations

## 5.1 Implementation

Maintained implementations for estimating RMST under right censoring remain limited. Available tools include `survRM2` (Hajime et al., 2015), `grf` (Tibshirani et al., 2017), and `RISCA` (Foucher et al., 2019). The Appendix details these packages and their functions, mapping each to the estimators studied in this paper. We implement the remaining estimators and, for a unified and pedagogical framework, provide custom implementations for all estimators, including those with existing software. The complete codebase is available on GitHub at `github.com/Sanofi-Public/causal_survival_analysis`, with function-level usage documented in the Appendix 6.

Building on these implementations, we compare the behaviors and performances of the estimators through simulations conducted across various experimental contexts. These contexts include scenarios based on RCTs and observational data, with both independent and dependent censoring (see Appendix 6 for details about the data generating process). For each setting, we generated 100 independent datasets, computed the estimator on each, and summarized the resulting distribution using boxplots.
## 5.2 RCT

We sample $n$ i.i.d. observations $(X_i, A_i, C_i, T_i^{(0)}, T_i^{(1)})_{i \in [n]}$ where covariates $X \sim \mathcal{N}(\mu, \Sigma)$, treatment is randomized with $e(X) = \Pr(A = 1 | X) = 0.5$, and the potential event times follow Cox models. Specifically, the control potential time has hazard $\lambda^{(0)}(t | X) = \lambda_0 \exp\{\beta_0^\top X\}$, and the treated potential time is obtained via a constant additive shift, $T^{(1)} = T^{(0)} + 10$. Censoring is either independent, with a constant hazard $\lambda_C(t) = \lambda_{C,0}$ (Scenario 1), or conditionally independent given $X$, with a Cox-type hazard $\lambda_C(t | X) = \lambda_{C,0} \exp\{\beta_C^\top X\}$ (Scenario 2). The truncation horizon is fixed at $\tau = 25$. For RCT scenarios 1 and 2, the ground truth value of the RMST difference at time $\tau = 25$ was computed analytically and found to be approximately 7.1.

For each scenario, we estimate the RMST difference using the methods summarized in Table 3. The methods used to estimate the nuisance components are indicated in brackets: either logistic regression or random forests for propensity scores and either Cox models or survival random forests for survival and censoring models. A naive estimator where censored observations are simply removed and the survival time is averaged for treated and controls is also provided for a naive baseline.

Figure 7 shows the distribution of the RMST difference for 100 simulations in Scenario 1 and different sample sizes: 500, 1000, 2000, 4000. The true value of $\theta_{\text{RMST}}$ is indicated by red dashed line. For sake of clarity, only a subset of estimators with particularly informative or interesting results is displayed. In this setting, and in accordance with the theory, the

<figure><img src="image_7.png" /><figcaption>Figure 7: ATEs estimation in RCT setting with independent censoring.</figcaption></figure>

simplest estimator (unadjusted KM) performs just as well as the others, and presents, as expected, an extremely small bias.

The naive estimator is biased, as expected, and the bias in the G-formula S-learner arises because the treatment effect is additive $T^{(1)} = T^{(0)} + 10$ and violates the assumption that $T$ would follow a Cox model in the variables $(X, A)$. However, $T|A = a$ is a Cox-model
for $a \in \{0, 1\}$, which explains the remarkable performance of G-formula (Cox/T-learners) achieves the lowest finite-sample variance. Some forest-based versions also yield nearly unbiased results, although their performance can vary: for instance, the T-learner remains well-behaved, while the S-learner appears more sensitive and may require larger sample sizes to achieve convergence. For completeness, we have also evaluated the other estimators (e.g. IPTW KM, IPCW, IPCW KM, BJ, IPTW-IPCW, IPTW-IPCW KM, IPTW-BJ, AIPTW-AIPCW, causal survival forest) in Figure 11 (Appendix 6). As expected, they yield no convergence gains, exhibit larger finite-sample variance, and rely on unnecessary model specification.

<figure><img src="image_8.png" /><figcaption>Figure 8: ATEs estimation in RCT setting with conditionally independent censoring.</figcaption></figure>

Figure 8 shows the results for the RCT simulation with conditionally independent censoring (Scenario 2), focusing on a subset of estimators. The naive estimator is biased, as expected. The IPCW (Cox) is slightly biased up to 2,000 observations. In contrast, both IPCW KM (Cox) and Buckley-James (Cox) are unbiased even with only 500 observations, with a smaller variance for BJ. G-formula (Cox/T-learners) performs well, even with small sample size. As in Scenario 1, G-formula (Cox/S-learner) remains biased. The forest-based estimators are biased. Notably, all estimators exhibit higher variability compared to Scenario 1 likely due to the higher censoring rate in this scenario and the estimation of censoring model.

For completeness, we have also evaluated the other estimators (e.g. IPTW KM, IPTW-IPCW, IPTW-IPCW KM, IPTW-BJ, AIPTW-AIPCW, causal survival forest) in Figure 12 (Appendix 6). Introducing a superfluous treatment model via the propensity score inflates finite-sample variance, with a marked effect at 500 observations.
## 5.3 Observational data

We simulate an observational study under both independent and conditionally independent censoring (as in Scenarios 1–2). The data consist of an i.i.d. *n*-sample $(X_i, A_i, C_i, T_i^{(0)}, T_i^{(1)})_{i \in [n]}$ generated as in Section 5.2, except that treatment allocation now follows a logistic propensity model, $\text{logit}(e(X)) = \beta_A^\top X$, so that $e(X)$ is no longer constant. The survival data-generating mechanism is unchanged; hence the target parameter $\theta_{\text{RMST}}$ coincides with that of the RCT simulations, namely 7.1 (see Section 5.2).

Figure 9 below shows the distribution of the estimators of $\theta_{\text{RMST}}$ for the observational study with independent censoring. As expected, the naive estimator is biased while IPTW KM

<figure><img src="image_9.png" /><figcaption>Figure 9: ATEs estimation in observational setting with independent censoring.</figcaption></figure>

(Log.Reg.) is unbiased. The top-performing estimator in this scenario is the G-formula (Cox/T-learners), which remains unbiased even with as few as 500 observations and achieves the lowest variance. All estimators that use forests to estimate nuisance parameters are biased across sample sizes from 500 to 8000. This setting thus highlights that one should either have an a priori knowledge of the specification of the models or large sample size.

Figure 10 below shows the distribution of the $\theta_{\text{RMST}}$ estimates for the observational study with conditionally independent censoring. The top-performing estimators in this scenario are G-formula (Cox/T-learners) and AIPCW-AIPTW (Cox & Cox & Log.Reg.), which are unbiased even with 500 observations. The G-formula (Cox/T-learners) has the lowest variance. Surprisingly, the G-formula (Cox/S-learner) performs quite competitively, showing only a slight bias despite the violation of the proportional hazards assumption. IPTW-IPCW KM (Cox & Log.Reg.) and IPTW-IPCW (Cox & Log.Reg.) exhibit different finite-sample behaviors. The Kaplan–Meier version is already unbiased with 500 observations, whereas the moment-based estimator requires around 8,000 observations to achieve unbiasedness. This pattern mirrors the results observed for IPCW and IPCW KM in the randomized trial scenario with conditionally independent censoring. All estimators that use forests to estimate nuisance parameters are biased across sample sizes from 500 to 8000. Although
<figure><img src="image_10.png" /><figcaption>Figure 10: ATEs estimation in observational setting with conditionally independent censoring.</figcaption></figure>

Causal Survival Forest (Cui et al., 2023) and AIPTW-APCW (Forest) are expected to converge, they remain extremely demanding in terms of sample size.

For completeness, we have also evaluated the other estimators (e.g. KM, IPCW KM, IPCW, BJ, IPTW KM and package versions) in Figure 14 (Appendix 6).

## 5.4 Misspecification of nuisance components

We generate an observational study with covariate interactions and conditionally independent censoring. The objective is to assess the impact of misspecifying nuisance components; specifically, we will use models that omit interactions to estimate these components. This approach enables us to evaluate the robustness properties of various estimators. In addition, in this setting forest-based methods are expected to behave better.

We generate an i.i.d. sample $(X_i, A_i, C_i, T_i^{(0)}, T_i^{(1)})_{i=1}^n$ where $X \sim \mathcal{N}(\mu, \Sigma)$, the potential event times and censoring follow Cox models with interactions, and treatment follows a logistic propensity with interactions. Specifically, the control potential time has hazard $\lambda^{(0)}(t | X) = \exp\{\beta_0^\top Z(X)\}$, the treated potential time is obtained via a constant shift, $T^{(1)} = T^{(0)} + 0.3$, censoring follows $\lambda_C(t | X) = \exp\{\beta_C^\top Z(X)\}$, and the treatment allocation satisfies $\text{logit}(e(X)) = \beta_A^\top Z(X)$. Under correct specification, the nuisance models use $(X, Z(X))$ as inputs; under misspecification, only $X$ and a strict subset of $Z(X)$ (the main-effect terms) are supplied. The true RMST difference at time 1 is estimated to be 0.21.

We first confirm consistency by estimating $\theta_{\text{RMST}}$ at 10,000 observations with correctly specified nuisance models and assess robustness under (i) single-model misspecification (Appendix, Figure 15) and (ii) paired or more misspecifications—(outcome, censoring), (outcome, treatment), (treatment, censoring) and (outcome, treatment, censoring) (Appendix, Figure 16). Given that, in the parametric setting, IPTW-IPCW KM and IPTW-IPCW
behaved similarly asymptotically in the previous simulations, we report results for IPTW-IPCW KM only.

As expected, when exactly one nuisance component is misspecified (censoring, treatment, or outcome; see Appendix, Figure 15), AIPTW-AIPCW (Cox & Cox & Log.Reg) remains unbiased. In contrast, the other estimators are biased when any of their nuisance models is misspecified; for example, if the treatment model is misspecified, IPTW-IPCW KM is biased. This finding was previously reported in Ozenne et al. (2020).

When two nuisance components are misspecified (see Appendix, figure 16), AIPTW-AIPCW (Cox & Cox & Log.Reg.) seems to converge in case where the treatment and censoring models are misspecified. Furthermore, when all nuisance models are misspecified, all estimators exhibit bias.

# 6 Conclusion

**Summary of objectives and findings** Restricted mean survival time (RMST)-based estimators provide a valuable alternative to hazard ratios, offering interpretable and robust effect estimates that are particularly well-suited for complex, real-world data. In this work, we provided a comprehensive overview of causal survival estimators targeting RMST, emphasizing both their theoretical foundations and practical performance. As well as implementation (available on a public GitHub repository), we reviewed a wide range of methods—including regression-based models, weighting approaches, or hybrid estimators—and assessed their behavior in both their parametric and non-parametric versions through extensive simulations. These scenarios covered diverse settings, including randomized trials and observational data, with varying censoring mechanisms and model specifications. Particular attention was given to robustness, sensitivity to misspecification, and estimator variability in finite samples.

**Practical recommendations for RCTs** Whether subject to independent or dependent censoring, we advise against overcomplicating estimation strategies. Classical estimators like the Kaplan-Meier estimator (with weights for dependent censoring if necessary) remain highly effective in these settings without requiring adjustment on covariates. We also recommend using the parametric G-formula. Although randomization can make certain covariates seem unnecessary, theoretical results (Karrison and Kocherginsky, 2018) and our simulations show that adjusting for prognostic covariates reduces the asymptotic variance of the RMST-difference estimator. We also advocate for the use of the Buckley-James (BJ) transformation in this context, which demonstrated excellent performance in our simulations, particularly in reducing bias under independent or dependent censoring without introducing excessive variability. Moreover, BJ offers an attractive compromise between ease of implementation and statistical efficiency.

**Practical recommendations for Observational studies** Non-robust weighted estimators such as IPTW-IPCW and IPTW-BJ were consistently marked by high variability, particularly when applied to dependent censoring. Parametric G-formula estimators yielded the best performance under correct model specification (as in Denz et al. (2023)), showing both low
finite sample bias and variance. However, their sensitivity to model misspecification remains a limitation. Under dependent censoring, we therefore recommend more robust estimators such as AIPTW-AIPCW, which maintain good performance despite their complexity. Parametric estimators perform best when correctly specified; yet with interactions or nonlinearities, nonparametric methods can dominate, as they both avoid functional assumptions and automatically learn complex interactions that parametric models would require explicit specification for. That said, parametric versions may still be preferable in small to moderate sample sizes, as they exhibited lower bias and variance in our simulations, especially for samples under 4,000 observations. In contrast, when sample sizes are sufficiently large, non-parametric estimators become viable alternatives. In such cases, we specifically recommend the AIPTW-AIPCW estimator, since the non-parametric G-formula lacks formal asymptotic guarantees and is prone to plug-in bias (Chernozhukov et al., 2018).

**Limitations** It is important to highlight that our simulation design may disadvantage nonparametric methods, as survival and censoring times were generated from exponential distributions and treatment assignments via logistic regression, naturally favoring parametric models.

**Perspectives** Future work should address variable selection in causal survival analysis. In causal inference, selection serves for first identification to adjust for confounders while avoiding instrumental variables (IVs) that inflate variance (Lunceford and Davidian, 2004) and second to efficiency—include precision variables predictive of outcome (Lunceford and Davidian, 2004; Rotnitzky and Smucler, 2019). With survival data, censoring introduces a third class: covariates tied only to censoring; whether to include them (e.g., in IPCW) is unclear, as they may reduce bias in RMST differences yet raise variance. Systematic evaluation of confounders, precision variables, IVs, and censoring-only covariates on RMST bias-variance across modeling strategies is warranted. Other estimators could be benchmarked such as pseudo-observation approaches (Andersen et al., 2017; Wang, 2018; van Hage et al., 2024) and targeted maximum likelihood estimation (TMLE) (Stitelman et al., 2012; Talbot et al., 2025).

**Acknowledgement**

The research was supported by Sanofi and INRIA (Institut National de Recherche en Informatique et en Automatique) through the CIFRE program (Convention Industrielle de Formation par la Recherche). Charlotte Voinot and Bernard Sebastien are Sanofi employees and may hold shares and/or stock options in the company. The work of CB was supported by the Deutsche Foschungsgemeinschaft (German Research Foundation) on the French-German PRCI ANR ASCAI CA 1488/4-1 "Aktive und Batch-Segmentierung, Clustering und Seriation: Grundlagen der KI". Imke Mayer and Julie Josse have nothing to disclose.

**Conflict of Interest**

The authors declare no conflicts of interest.
References

Aalen, O., Borgan, O., and Gjessing, H. (2008). *Survival and Event History Analysis: A Process Point of View*. Statistics for Biology and Health. Springer New York.

Andersen, P. K., Syriopoulou, E., and Parner, E. T. (2017). Causal inference in survival analysis using pseudo-observations. *Statistics in Medicine*, 36(17):2669–2681.

Averbuch, I., Bareket-Samish, A., Goldstein, D. A., Eizenstein, S., Markel, G., Rosenbaum, E., Limon, D., Bomze, D., Ludmir, E. B., and Meirson, T. (2025). Challenges in interpreting survival metrics in clinical trials: The utility of restricted mean survival analyses. *International Journal of Radiation Oncology-Biology*\* *Physics*, 122(2):224–233.

Breslow, N. (1974). Covariance analysis of censored survival data. *Biometrics*, 30(1):89–99.

Buckley and James (1979). Linear regression with censored data. *Biometrika*, 66(3):429–436.

Chen, K. and Lo, S.-H. (1996). On bootstrap accuracy with censored data. *Ann. Stat.*, 24(2):569–595.

Chen, P.-Y. and Tsiatis, A. A. (2001). Causal inference on the difference of the restricted mean lifetime between two groups. *Biometrics*, 57(4):1030–1038.

Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. *The Econometrics Journal*, 21(1):C1–C68.

Colnet, B., Josse, J., Varoquaux, G., and Scornet, E. (2022). Reweighting the rct for generalization: finite sample error and variable selection. *arXiv preprint arXiv:2208.07614*.

Cox, D. R. (1972). Regression models and life-tables. *Journal of the Royal Statistical Society: Series B (Methodological)*, 34(2):187–202.

Cui, Y., Kosorok, M. R., Sverdrup, E., Wager, S., and Zhu, R. (2023). Estimating heterogeneous treatment effects with right-censored data via causal survival forests. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 85(2):179–211.

Denz, R., Klaßen-Mielke, R., and Timmesfeld, N. (2023). A comparison of different methods to adjust survival curves for confounders. *Statistics in Medicine*, 42(10):1461–1479.

Efron, B. (1981). Censored data and the bootstrap. *J. Am. Stat. Assoc.*, 76(374):312–319.

Efron, B. and Tibshirani, R. J. (1994). *An introduction to the bootstrap*. Chapman and Hall/CRC.

Fan, J. and Gijbels, I. (1994). *Local polynomial modelling and its applications*. Chapman & Hall/CRC.

Foster, J. C., Taylor, J. M. G., and Ruberg, S. J. (2011). Subgroup identification from randomized clinical trial data. *Stat. Med.*, 30(24):2867–2880.
Foucher, Y., Le Borgne, F., and Chatton, A. (2019). Risca: Causal inference and prediction in cohort-based analyses.

Gill, R. (1983). Large sample behaviour of the product-limit estimator on the whole line. *The annals of statistics*, pages 49–58.

Greenland, S., Robbins, J. M., and Pearl, J. (1999). Confounding and collapsibility in causal inference. *Statistical Science*, 14:29–46.

Hájek, J. (1971). Comment on "an essay on the logical foundations of survey sampling" by d. basu. In Godambe, V. P. and Sprott, D. A., editors, *Foundations of Statistical Inference*, pages 236–239. Holt, Rinehart and Winston, New York. ISBN 0039281035.

Hajime, U., Lu, T., Miki, H., Angel, C., Chakib, B., and James, B. (2015). *survrm2: Comparing restricted mean survival time*.

Hernán, M. A. (2010). The hazards of hazard ratios. *Epidemiology*, 21(1):13–15.

Hernán MA, R. J. (2020). *Causal Inference: What If*. Boca Raton: Chapman & Hall/CRC.

Hirano, K., Imbens, G. W., and Ridder, G. (2003). Efficient estimation of average treatment effects using the estimated propensity score. *Econometrica*, 71(4):1161–1189.

Howe, C. J., Cole, S. R., Lau, B., Napravnik, S., and Eron, J. J. (2016). Selection bias due to loss to follow up in cohort studies. *Epidemiology*, 27 1:91–7.

Huitfeldt, A., Stensrud, M., and Suzuki, E. (2019). On the collapsibility of measures of effect in the counterfactual causal framework. *Emerging Themes in Epidemiology*, 16:1–12.

Ishwaran, H., Kogalur, U. B., Blackstone, E. H., and Lauer, M. S. (2008). Random survival forests. *The Annals of Applied Statistics*, 2(3):841–860.

Kalbfleisch, J. D. and Prentice, R. L. (2002). *The statistical analysis of failure time data*. John Wiley & Sons.

Kaplan, E. L. and Meier, P. (1958). Nonparametric estimation from incomplete observations. *Journal of the American Statistical Association*, 53(282):457–481.

Karrison, T. and Kocherginsky, M. (2018). Restricted mean survival time: Does covariate adjustment improve precision in randomized clinical trials? *Clinical Trials*, 15(2):178–188. PMID: 29502444.

Koul, H., Susarla, V., and Ryzin, J. V. (1981). Regression Analysis with Randomly Right-Censored Data. *The Annals of Statistics*, 9(6):1276 – 1288.

Künzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019). Metalearners for estimating heterogeneous treatment effects using machine learning. *Proceedings of the National Academy of Sciences*, 116(10):4156–4165.

Lunceford, J. K. and Davidian, M. (2004). Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study. *Statistics in Medicine*, 23(19):2937–2960.
Makary, M. A. and Prasad, V. (2025). Priorities for a new fda. *JAMA*.

Martinussen, T. and Vansteelandt, S. (2013). On collapsibility and confounding bias in cox and aalen regression models. *Lifetime Data Anal.*, 19(3):279–296.

Martinussen, T., Vansteelandt, S., and Andersen, P. (2020). Subtleties in the interpretation of hazard contrasts. *Lifetime Data Analysis*, 26.

Overgaard, M. (2025). A comparison of Kaplan-Meier-based inverse probability of censoring weighted regression methods. *Lifetime Data Anal.*, 31(4):747–783.

Ozenne, B. M. H., Scheike, T. H., Stærk, L., and Gerds, T. A. (2020). On the estimation of average treatment effects with right-censored time to event outcome and competing risks. *Biometrical Journal*, 62(3):751–763.

Robins, J. M. and Finkelstein, D. M. (2000). Correcting for noncompliance and dependent censoring in an aids clinical trial with inverse probability of censoring weighted (ipcw) log-rank tests. *Biometrics*, 56(3):779–788.

Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994). Estimation of regression coefficients when some regressors are not always observed. *Journal of the American Statistical Association*, 89(427):846–866.

Robins, J. M., Ángel Hernán, M., and Brumback, B. (2000). Marginal structural models and causal inference in epidemiology. *Epidemiology*, 11(5):550–560.

Rosenbaum, P. R. and Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. *Biometrika*, 70(1):41–55.

Rotnitzky, A. and Smucler, E. (2019). Efficient adjustment sets for population average treatment effect estimation in non-parametric causal graphical models.

Royston, P. and Parmar, M. K. (2013). Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome. *BMC medical research methodology*, 13:1–15.

Rubin, D. and van der Laan, M. J. (2007). A doubly robust censoring unbiased transformation. *The International Journal of Biostatistics*, 3(1).

Rubin, D. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. *Journal of Educational Psychology*, 66(5):688–701.

Saikia, R. and Barman, M. P. (2017). A review on accelerated failure time models. *International Journal of Statistics and Systems*, 12(2):311–322.

Satten, G. A. and Datta, S. (2001). The kaplan–meier estimator as an inverse-probability-of-censoring weighted average. *The American Statistician*, 55(3):207–210. PMID: 28845048.

Schaubel, D. E. and Wei, G. (2011). Double Inverse-Weighted Estimation of Cumulative Treatment Effects under Nonproportional Hazards and Dependent Censoring. *Biometrics*, 67(1):29–38.
Stitelman, O. M., Gruttola, V. D., and van der Laan, M. J. (2012). A general implemen-
tation of tmle for longitudinal data applied to causal inference in survival analysis. *The
International Journal of Biostatistics*, 8(1).

Talbot, D., Diop, A., Mésidor, M., Chiu, Y., Sirois, C., Spieker, A. J., Pariente, A., Noize, P., Simard, M., Luque Fernandez, M. A., Schomaker, M., Fujita, K., Gnjidic, D., and Schnitzer, M. E. (2025). Guidelines and best practices for the use of targeted maximum likelihood and machine learning when estimating causal effects of exposures on time-to-event outcomes. *Statistics in Medicine*, 44(6):e70034.

Therneau, T. M. (2001). survival: Survival analysis.

Tibshirani, J., Athey, S., Sverdrup, E., and Wager, S. (2017). grf: Generalized random forests.

Turkson, A., Ayiah-Mensah, F., and Nimoh, V. (2021). Handling censoring and censored data in survival analysis: A standalone systematic literature review. *International Journal of Mathematics and Mathematical Sciences*, 2021:1–16.

Van der Laan, M. J. and Robins, J. M. (2003). *Introduction*, pages 8–101. Springer New York, New York, NY.

van Hage, P., le Cessie, S., van Maaren, M. C., Putter, H., and van Geloven, N. (2024). Doubly robust estimation of marginal cumulative incidence curves for competing risk analysis.

Wang, J. (2018). A simple, doubly robust, efficient estimator for survival functions using pseudo observations. *Pharmaceutical Statistics*, 17(1):38–48.

Willems, S., Schat, A., van Noorden, M., and Fiocco, M. (2018). Correcting for dependent censoring in routine outcome monitoring data by applying the inverse probability censoring weighted estimator. *Statistical Methods in Medical Research*, 27(2):323–335. PMID: 26988930.

Xie, J. and Liu, C. (2005). Adjusted kaplan–meier estimator and log-rank test with inverse probability of treatment weighting for survival data. *Statistics in Medicine*, 24(20):3089–3110.

Zhang, Z. (2016). Parametric regression model for survival data: Weibull regression model as an example. *Annals of translational medicine*, 4(24).

Zhou, M. (1988). Two-sided bias bound of the kaplan-meier estimator. *Probability theory and related fields*, 79(2):165–173.
# Appendix

## A.1. Appendix A: Proofs

**Convergence notions and consistency.** Let $(\hat{\theta}_n)_{n \in \mathbb{N}}$ be a sequence of estimators for a parameter $\theta$. We recall three standard notions of convergence, commonly used to define statistical consistency:

1. **Weak consistency** means that $\hat{\theta}_n$ converges in probability to $\theta$, i.e.,

$$ \hat{\theta}_n \xrightarrow{P} \theta \quad \text{if and only if} \quad \forall \varepsilon > 0, \mathbb{P}(|\hat{\theta}_n - \theta| > \varepsilon) \to 0 \quad \text{as } n \to \infty. $$

This implies that the estimator gets arbitrarily close to the true value with high probability as the sample size increases.

2. **Strong consistency** means that $\hat{\theta}_n$ converges almost surely to $\theta$, i.e.,

$$ \hat{\theta}_n \xrightarrow{\text{a.s.}} \theta \quad \text{if and only if} \quad \mathbb{P}\left(\lim_{n \to \infty} \hat{\theta}_n = \theta\right) = 1. $$

This is a stronger form of convergence, requiring the sequence to converge pointwise on almost every sample path.

3. **Convergence in distribution** (or **weak convergence** in the sense of probability theory) refers to the convergence of the distribution of $\hat{\theta}_n$ to that of a limiting random variable $Z$, not necessarily equal to $\theta$. It is denoted:

$$ \hat{\theta}_n \xrightarrow{d} Z, $$

and defined by the condition that for every bounded continuous function $f$,

$$ \mathbb{E}[f(\hat{\theta}_n)] \to \mathbb{E}[f(Z)]. $$

This notion is typically used in the context of asymptotic distributional results, such as central limit theorems.

### A.1.1. Proofs of Section 2.1

Note that not all proofs are novel work. Consistency results on many survival estimators have been reported multiple times as explained in the main paper. For clarity we still recall them. We indicate when the proofs are not novel or when similar proofs exist elsewhere. When we indicate nothing, this means that we have not found those results in other published work.

**Proposition 6.1** (Consistency of the Kaplan-Meier estimator). (Adapted from Gill (1983), Zhou (1988))

Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.3 (Independent Censoring) and 2.4 (Positivity of the Censoring Process), and for
all $t \in [0, \tau]$, the estimator $\hat{S}_{\text{KM}}(t|A = a)$ of $S^{(a)}(t)$ is strongly consistent and admits the following bounds for its bias:

$$0 \le S^{(a)}(t) - \mathbb{E}[\hat{S}_{\text{KM}}(t|A = a)] \le O(\mathbb{P}(N_k(a) = 0)),$$

where $k$ is the greatest time $t_k$ such that $t \ge t_k$.

The bound we give, although slightly looser, still exhibits the same asymptotic regime. In particular, as soon as $S^{(a)}(t) > 0$ (and Assumption 2.4 holds), then the bias decays exponentially fast towards 0.

**Proof. (Proposition 6.1).** Consistency is a trivial consequence of the law of large number. To show that $\hat{S}_{\text{KM}}$ is unbiased, let us introduce $\mathcal{F}_k$ be the filtration generated by the set of variables

$$\{A_i, \mathbb{I}\{\tilde{T}_i = t_j\}, \mathbb{I}\{\tilde{T}_i = t_j, \Delta_i = 1\} | j \in [k], i \in [n]\}.$$

which corresponds to the known information up to time $t_k$, so that $D_k(a)$ is $\mathcal{F}_k$-measurable but $N_k(a)$ is $\mathcal{F}_{k-1}$-measurable. One can write that, for $k \ge 2$

$$\begin{aligned}
\mathbb{E}[\mathbb{I}\{\tilde{T}_i = t_k, \Delta_i = 1, A_i = a\} | \mathcal{F}_{k-1}] &= \mathbb{E}[\mathbb{I}\{\tilde{T}_i = t_k, \Delta_i = 1, A_i = a\} | \mathbb{I}\{\tilde{T}_i \ge t_k\}, A_i] \\
&= \mathbb{I}\{A_i = a\}\mathbb{E}[\mathbb{I}\{T_i = t_k, C_i \ge t_k\} | \mathbb{I}\{T_i \ge t_k, C_i \ge t_k\}, A_i] \\
&= \mathbb{I}\{A_i = a\}\mathbb{I}\{C_i \ge t_k\}\mathbb{E}[\mathbb{I}\{T_i = t_k\} | \mathbb{I}\{T_i \ge t_k\}, A_i] \\
&= \mathbb{I}\{\tilde{T}_i \ge t_k, A_i = a\} \left( 1 - \frac{S^{(a)}(t_k)}{S^{(a)}(t_{k-1})} \right),
\end{aligned}$$

where we used that $T_i(a)$ is independent from $A_i$ by Assumption 2.1 (Random Treatment Assignment). We then easily derive from this that

$$\mathbb{E}\left[\left(1 - \frac{D_k(a)}{N_k(a)}\right) \mathbb{I}\{N_k(a) > 0\} \middle| \mathcal{F}_{k-1}\right] = \frac{S^{(a)}(t_k)}{S^{(a)}(t_{k-1})} \mathbb{I}\{N_k(a) > 0\},$$

and then that

$$\mathbb{E}[\hat{S}_{\text{KM}}(t_k|A = a)|\mathcal{F}_{k-1}] = \frac{S^{(a)}(t_k)}{S^{(a)}(t_{k-1})} \hat{S}_{\text{KM}}(t_{k-1}|A = a) + O(\mathbb{I}\{N_k(a) = 0\}),$$

By induction, we easily find that

$$\mathbb{E}[\hat{S}_{\text{KM}}(t|A = a)] = \prod_{t_j \le t} \frac{S^{(a)}(t_j)}{S^{(a)}(t_{j-1})} + O\left(\sum_{t_j \le t} \mathbb{P}(N_j(a) = 0)\right) = S^{(a)}(t) + O(\mathbb{P}(N_k(a) = 0))$$

where $t_k$ is the greatest time such that $t_k \le t$.

**Proposition 6.2 (Variance of the Kaplan-Meier estimator).** (Adapted from Kaplan and Meier (1958))

Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.3 (Independent Censoring) and 2.4 (Positivity of the Censoring Process), and for
all $t \in [0, \tau]$, $\hat{S}_{\text{KM}}(t|A = a)$ is asymptotically normal and $\sqrt{n}(\hat{S}_{\text{KM}}(t|A = a) - S^{(a)}(t))$ converges in distribution towards a centered Gaussian of variance

$$V_{\text{KM}}(t|A=a) := S^{(a)}(t)^2 \sum_{t_k \le t} \frac{1-s_k(a)}{s_k(a)r_k(a)},$$

where $s_k(a) = S^{(a)}(t_k)/S^{(a)}(t_{k-1})$ and $r_k(a) = \mathbb{P}(\tilde{T} \ge t_k, A = a)$.

The proof of Proposition 6.2 is derived below. Because $D_k(a)/N_k(a)$ is a natural estimator of $1 - s_k(a)$ and, $\frac{1}{n}N_k(a)$ a natural estimator for $r_k(a)$, the asymptotic variance of the Kaplan-Meier estimator can be estimated with the so-called Greenwood formula, as already derived heuristically in Kaplan and Meier (1958):

$$\widehat{\text{Var}}(\hat{S}_{\text{KM}}(t|A=a)) := \hat{S}_{\text{KM}}(t|A=a)^2 \sum_{t_k \le t} \frac{D_k(a)}{N_k(a)(N_k(a) - D_k(a))}. \quad (25)$$

*Proof. (Proposition 6.2).* The asymptotic normality is a mere consequence of the joint asymptotic normality of $(N_k(a), D_k(a))_{t_k \le t}$ with an application of the $\delta$-method. To access the asymptotic variance, notice that, using a similar reasoning as in the previous proof:

$$\begin{aligned} \mathbb{E}[(1 - D_k(a)/N_k(a))^2 | \mathcal{F}_{k-1}] &= \mathbb{E}[1 - D_k(a)/N_k(a) | \mathcal{F}_{k-1}(a)]^2 + \frac{1}{N_k(a)^2} \text{Var}(D_k(a) | \mathcal{F}_{k-1}) \\ &= s_k^2(a) + \frac{s_k(a)(1 - s_k(a))}{N_k(a)} \mathbb{I}\{N_k(a) > 0\} + \mathcal{O}(\mathbb{I}\{N_k(a) = 0\}). \end{aligned}$$

Now we know that $N_k(a) = nr_k(a) + \sqrt{n}O_{\mathbb{P}}(1)$, with the $O_{\mathbb{P}}(1)$ having uniformly bounded moments. So that we deduce that

$$\mathbb{E}[(1 - D_k(a)/N_k(a))^2 | \mathcal{F}_{k-1}] = s_k^2(a) + \frac{s_k(a)(1 - s_k(a))}{n r_k(a)} + \frac{1}{n^{3/2}} O_{\mathbb{P}}(1),$$

where $O_{\mathbb{P}}(1)$ has again bounded moments. Using this identity, we find that

$$\begin{aligned} n \text{Var} \hat{S}_{\text{KM}}(t|A=a) &= n \left( \mathbb{E} S_{\text{KM}}(t|A=a)^2 - S^{(a)}(t)^2 \right) \\ &= n S^{(a)}(t)^2 \left( \mathbb{E} \left[ \prod_{t_k \le t} \left( 1 + \frac{1-s_k(a)}{n s_k(a) r_k(a)} + \frac{1}{n^{3/2}} O_{\mathbb{P}}(1) \right) \right] - 1 \right). \end{aligned}$$

Expending the product and using that the $O_{\mathbb{P}}(1)$'s have bounded moments, we finally deduce that

$$\mathbb{E} \left[ \prod_{t_k \le t} \left( 1 + \frac{1-s_k(a)}{n s_k(a) r_k(a)} + \frac{1}{n^{3/2}} O_{\mathbb{P}}(1) \right) \right] - 1 = \frac{1}{n} \sum_{t_k \le t} \frac{1-s_k(a)}{s_k(a)r_k(a)} + \frac{1}{n^{3/2}} O(1),$$

$$n \text{Var} \hat{S}_{\text{KM}}(t|A=a) = V_{\text{KM}}(t|A=a) + O(n^{-1/2}),$$

which is what we wanted to show. □
### A.1.2. Proofs of Section 2.2

**Proposition 6.3** (IPCW is a Censoring Unbiased Transformation). Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), the IPCW transform 10 is a censoring-unbiased transformation in the sense of Equation 6.

*Proof.* (Proposition 6.3). Assumption 2.6 allows the transformation to be well-defined. Furthermore, it holds

$$ 
\begin{aligned}
E[T_{\text{IPCW}}^{*}|A=a, X] &= E\left[\frac{\Delta^{\tau} \times \tilde{T} \wedge \tau}{G(\tilde{T} \wedge \tau|A, X)} \middle| A=a, X\right] \\
&= E\left[\frac{\Delta^{\tau} \times T^{(a)} \wedge \tau}{G(T^{(a)} \wedge \tau|A, X)} \middle| A=a, X\right] \\
&= E\left[E\left[\frac{\mathbb{I}\{T^{(a)} \wedge \tau \le C\} \times T^{(a)} \wedge \tau}{G(T^{(a)} \wedge \tau|A, X)} \middle| A, X, T^{(a)}\right] \middle| A=a, X\right] \\
&= E\left[T^{(a)} \wedge \tau|A=a, X\right] \\
&= E\left[T^{(a)} \wedge \tau|X\right].
\end{aligned}
$$

We used in the second equality that on the event $\{\Delta^{\tau} = 1, A = a\}$, it holds $\tilde{T} \wedge \tau = T \wedge \tau = T^{(a)} \wedge \tau$. We used in the fourth equality that $G(T^{(a)} \wedge \tau|A, X) = E[\mathbb{I}\{T^{(a)} \wedge \tau \le C\}|X, T^{(a)}, A]$ thanks to Assumption 2.5, and in the last one that A is independent of X and $T^{(a)}$ thanks to Assumption 2.1. $\square$

**Proposition 6.4** (Consistency of $S_{\text{IPCW}}^{*}(t|A=a)$ with known G). Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), and for all $t \in [0, \tau]$, the oracle estimator $S_{\text{IPCW}}^{*}(t|A=a)$ defined as in Definition 2.2 with $\hat{G} = G$ is a strongly consistent and asymptotically normal estimator of $S^{(a)}(t)$.

*Proof.* (Proposition 6.4). Similarly to the computations done in the proof of Proposition 6.3, it is easy to show that, for $t_k \le \tau$,

$$ 
\begin{aligned}
\mathbb{E}\left[\frac{\Delta^{\tau}}{G(t_k|X, A)}\mathbb{I}(\tilde{T} = t_k, A = a)\right] &= \mathbb{E}\left[\frac{\mathbb{I}(t_k \le C)}{G(t_k|X, A)}\mathbb{I}(T^{(a)} = t_k, A = a)\right], \\
&= \mathbb{P}(A = a)\mathbb{P}(T^{(a)} = t_k),
\end{aligned}
$$

and likewise that

$$ 
\begin{aligned}
\mathbb{E}\left[\frac{1}{G(t_k|X, A)}\mathbb{I}(\tilde{T} \ge t_k, A = a)\right] &= \mathbb{E}\left[\frac{1}{G(t_k|X, A)}\mathbb{I}(T^{(a)} \wedge C \ge t_k, A = a)\right] \\
&= \mathbb{E}\left[\frac{\mathbb{I}(T^{(a)} \ge t_k)\mathbb{I}(A = a)}{G(t_k|X, A)}\mathbb{E}[\mathbb{I}(C \ge t_k)|X, A]\right] \\
&= \mathbb{P}(A = a)\mathbb{P}(T^{(a)} \ge t_k),
\end{aligned}
$$
so that $\hat{S}_{\text{IPCW}}(t)$ converges almost surely towards the product limit

$$ \prod_{t_k \le t} \left( 1 - \frac{\mathbb{P}(T^{(a)} = t_k)}{\mathbb{P}(T^{(a)} \ge t_k)} \right) = S^{(a)}(t), $$

yielding strong consistency. Asymptotic normality is straightforward. □

**Proposition 6.5** (BJ is a censoring unbiased transformation). Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), the BJ transform 13 is a censoring-unbiased transformation in the sense of Equation 6.

*Proof. (Proposition 6.5).* There holds

$$ \begin{aligned} \mathbb{E}[T_{\text{BJ}}^*|X, A=a] &= \mathbb{E} \left[ \Delta^\tau T^{(a)} \wedge \tau + (1-\Delta^\tau) \frac{\mathbb{E}[T \wedge \tau \times \mathbb{I}\{T \wedge \tau > C\}|C, A, X]}{\mathbb{P}(T > C|C, A, X)} \bigg|_{X, A=a} \right] \\ &= \mathbb{E}[\Delta^\tau T^{(a)} \wedge \tau|X] + \underbrace{\mathbb{E} \left[ \mathbb{I}\{T \wedge \tau > C\} \frac{\mathbb{E}[T \wedge \tau \times \mathbb{I}\{T \wedge \tau > C\}|C, A, X]}{\mathbb{E}[\mathbb{I}\{T \wedge \tau \ge C\}|C, A, X]} \bigg|_{X, A=a} \right]}_{(*)}. \end{aligned} $$

Now we easily see that conditioning wrt $X$ in the second term yields

$$ \begin{aligned} (*) &= \mathbb{E}[\mathbb{E}[T \wedge \tau \times \mathbb{I}\{T \wedge \tau > C\}|C, A, X]|X, A=a] \\ &= \mathbb{E}[(1 - \Delta^\tau)T \wedge \tau|X, A=a] \\ &= \mathbb{E}[(1 - \Delta^\tau)T^{(a)} \wedge \tau|X], \end{aligned} $$

ending the proof. □

*Proof. (Theorem 2.1, BJ minimizes the mean squared error among CUT).* We let $T^* = \Delta^\tau \phi_1 + (1 - \Delta^\tau)\phi_0$ be a transformation of the form Equation 5. There holds

$$ \mathbb{E}[(T^* - T \wedge \tau)^2] = \mathbb{E}[\Delta^\tau(\phi_1 - T \wedge \tau)^2] + \mathbb{E}[(1 - \Delta^\tau)(\phi_0 - T \wedge \tau)^2]. $$

The first term is non-negative and is zero for the BJ transformation. Since $\phi_0$ is a function of $(\tilde{T}, X, A)$ and that $\tilde{T} = C$ on $\{\Delta^\tau = 0\}$, the second term can be rewritten in the following way. We let $R$ be a generic quantity that does not depend on $\phi_0$.

$$ \begin{aligned} \mathbb{E}[(1 - \Delta^\tau)(\phi_0 - T)^2] &= \mathbb{E}[\mathbb{I}\{T \wedge \tau > C\}\phi_0^2 - 2\mathbb{I}\{T \wedge \tau > C\}\phi_0 T \wedge \tau] + R \\ &= \mathbb{E}\left[\mathbb{P}(T \wedge \tau > C|C, A, X)\phi_0^2 - 2\mathbb{E}[T \wedge \tau \mathbb{I}\{T \wedge \tau > C\}|C, A, X]\phi_0\right] + R \\ &= \mathbb{E}\left[\mathbb{P}(T \wedge \tau > C|C, A, X)\left(\phi_0 - \frac{\mathbb{E}[T \wedge \tau \mathbb{I}\{T \wedge \tau > C\}|C, A, X]}{\mathbb{P}(T \wedge \tau > C|C, A, X)}\right)^2\right] + R. \end{aligned} $$

Now the first term in the right hand side is always non-negative, and is zero for the BJ transformation. □
**Proposition 6.6** (AIPCW is a censoring-unbiased transformation). We let F, R be two conditional survival functions. Under Assumptions 1.1 (SUTVA), 2.1 (Random Treatment Assignment), 2.2 (Trial Positivity), 2.5 (Conditionally Independent Censoring) and if F also satisfies Assumption 2.6 (Conditional Positivity / Overlap for Censoring), and if F(·|X, A) is absolutely continuous wrt G(·|X, A), then the transformation $T_{\text{DR}}^* = T_{\text{DR}}^*(F, R)$ is a censoring-unbiased transformation in the sense of Equation 6 whenever $F = G$ or $R = S$.

A careful examination of the proofs from D. Rubin and Laan (2007) show that the proof translates straight away to our discrete setting.

### A.1.3. Proofs of Section 3.1

**Proposition 6.7** (Consistency of $S_{\text{IPTW}}^*(t|A = a)$ with known $e$). Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.3 (Independent Censoring) and 2.4 (Positivity of the Censoring Process) The oracle IPTW Kaplan-Meier estimator $S_{\text{IPTW}}^*(t|A = a)$ with $\hat{e} = e$ is a strongly consistent and asymptotically normal estimator of $S^{(a)}(t)$.

The proof of this result simply relies again on the law of large number and the $\delta$-method.

*Proof.* (Proposition 6.7). The fact that it is strongly consistent and asymptotically normal is again a simple application of the law of large number and of the $\delta$-method. We indeed find that, for $t_k \le \tau$

$$ 
\begin{aligned}
\mathbb{E} \left[ \frac{1}{e(X_i)} \mathbb{1}\{\tilde{T}_i = t_k, \Delta_i = 1, A_i = 1\} \right] &= \mathbb{E} \left[ \frac{A_i}{e(X_i)} \mathbb{1}\{T_i = t_k, C_i \ge t_k\} \right] \\
&= \mathbb{E} \left[ \mathbb{E} \left[ \frac{A_i}{e(X_i)} \mathbb{1}\{T_i = t_k, C_i \ge t_k\} \middle| X_i \right] \right] \\
&= \mathbb{E} \left[ \mathbb{E} \left[ \frac{A_i}{e(X_i)} \middle| X_i \right] \mathbb{P}(T_i = t_k | X_i) \mathbb{P}(C_i \ge t_k) \right] \\
&= \mathbb{P}(T_i = t_k) \mathbb{P}(C_i \ge t_k),
\end{aligned}
$$

where we used that A is independent of T conditionally on X, and that C is independent of everything. Likewise, one would get that

$$ \mathbb{E} \left[ \frac{1}{e(X_i)} \mathbb{1}\{\tilde{T}_i \ge t_k, A_i = 1\} \right] = \mathbb{P}(T_i \ge t_k) \mathbb{P}(C_i \ge t_k). $$

Similar computations hold for $A = 0$, ending the proof. $\square$
### A.1.4. Proofs of Section 3.2

**Proposition 6.8 (IPTW-IPCW is a censoring unbiased transformation).** Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), the IPTW-IPCW transform 10 is a censoring-unbiased transformation in the sense of Equation 6.

*Proof.* (Proposition 6.8). On the event {$\Delta^{\tau} = 1, A = 1$}, there holds $\tilde{T} \wedge \tau = T \wedge \tau = T^{(1)} \wedge \tau$, whence we find that,

$$ 
\begin{align*}
\mathbb{E}[T_{\text{IPCW}}^*|X, A=1] &= \mathbb{E}\left[\frac{A}{e(X)} \frac{\mathbb{I}\{T^{(1)} \wedge \tau \le C\}}{G(T^{(1)} \wedge \tau|X, A)} T^{(1)} \wedge \tau \middle| X\right] \\
&= \mathbb{E}\left[\frac{A}{e(X)} \mathbb{E}\left[\frac{\mathbb{I}\{T^{(1)} \wedge \tau \le C\}}{G(T^{(1)} \wedge \tau|X, A)} \middle| X, A, T^{(1)}\right] T^{(1)} \wedge \tau \middle| X\right] \\
&= \mathbb{E}\left[\frac{A}{e(X)} T^{(1)} \wedge \tau \middle| X\right] \\
&= \mathbb{E}\left[T^{(1)} \wedge \tau \middle| X\right],
\end{align*}
$$

and the same holds on the event $A = 0$. $\square$

*Proof.* (Proposition 3.1, Consistency of $\hat{\theta}_{\text{IPTW-IPCW}}$ with estimated $\hat{e}$ and $\hat{G}$). By consistency of $\hat{G}(\cdot|X, A)$ and $\hat{e}$ and by continuity, it suffices to look at the asymptotic behavior of the oracle estimator

$$ \theta_{\text{IPTW-IPCW}}^* = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{A_i}{e(X_i)} - \frac{1-A_i}{1-e(X_i)} \right) \frac{\Delta_i^{\tau}}{G(\tilde{T}_i \wedge \tau | A_i, X_i)} \tilde{T}_i \wedge \tau. $$

The later is converging almost towards its mean, which, following similar computations as in the previous proof, write

$$ 
\begin{align*}
\mathbb{E}\left[\left(\frac{A}{e(X)} - \frac{1-A}{1-e(X)}\right) \frac{\Delta^{\tau}}{G(\tilde{T} \wedge \tau|A, X)} \tilde{T} \wedge \tau\right] &= \mathbb{E}\left[\left(\frac{A}{e(X)} - \frac{1-A}{1-e(X)}\right) T \wedge \tau\right] \\
&= \mathbb{E}\left[T^{(1)} \wedge \tau\right] - \mathbb{E}\left[T^{(0)} \wedge \tau\right].
\end{align*}
$$

$\square$

**Proposition 6.9 (Consistency of $S_{\text{IPTW-IPCW}}^*(t|A=a)$ with known $e$ and $G$).** Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), and for all $t \in [0, \tau]$, if the oracle estimator $S_{\text{IPTW-IPCW}}^*(t|A=a)$ defined as in Definition 3.2 with $\hat{G}(\cdot|A, X) = G(\cdot|A, X)$ and $\hat{e} = e$ is a strongly consistent and asymptotically normal estimator of $S^{(a)}(t)$.
Proof. (Proposition 6.9). Asymptotic normality comes from a mere application of the δ-method, while strong consistency follows from the law of large number and the follozing computations. Like for the proof of Proposition 6.4, one find, by first conditionning wrt X, A, T(a), that, for t_k ≤ τ,

$$ \mathbb{E} \left[ \left( \frac{A}{e(X)} + \frac{1-A}{1-e(X)} \right) \frac{\Delta^\tau}{G(t_k|A, X)} \mathbb{I}\{\tilde{T} = t_k, A=a\} \right] = \mathbb{P}(T^{(a)} = t_k) $$

and likewise that

$$ \mathbb{E} \left[ \left( \frac{A}{e(X)} + \frac{1-A}{1-e(X)} \right) \frac{1}{G(t_k|A, X)} \mathbb{I}\{\tilde{T} \ge t_k, A=a\} \right] = \mathbb{P}(T^{(a)} \ge t_k) $$

so that indeed $S_{\text{IPTW-IPCW}}^*(t|A=a)$ converges almost surely towards $S^{(a)}(t)$. □

**Proposition 6.10** (IPTW-BJ is a censoring-unbiased transformation). Under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), the IPTW-BJ transform 13 is a censoring-unbiased transformation in the sense of Equation 6.

Proof. (Proposition 6.10). We write

$$ \mathbb{E}[T_{\text{IPTW-BJ}}^*|X, A=1] = \mathbb{E}\left[\frac{A}{e(X)}\Delta^\tau \times \tilde{T} \wedge \tau|X\right] + \mathbb{E}\left[\frac{A}{e(X)}(1-\Delta^\tau)Q_S(\tilde{T} \wedge \tau|A, X)|X\right]. $$

On the event $\{\Delta^\tau = 1, A=1\}$, there holds $\tilde{T} \wedge \tau = T \wedge \tau = T^{(1)} \wedge \tau$, whence we find that the first term on the right hand side is equal to

$$ \begin{aligned} \mathbb{E} \left[ \frac{A}{e(X)} \Delta^{\tau} \times \tilde{T} \wedge \tau \middle| X \right] &= \mathbb{E} \left[ \frac{A}{e(X)} \Delta^{\tau} \times T^{(1)} \wedge \tau \middle| X \right] \\ &= \mathbb{E} \left[ \Delta^{\tau} \times T^{(1)} \wedge \tau \middle| X \right]. \end{aligned} $$

For the second term in the right hand side, notice that on the event $\{\Delta^\tau = 0, A=1\}$, there holds $\tilde{T} = C < T^{(1)} \wedge \tau$, so that

$$ \begin{aligned} & \mathbb{E} \left[ \frac{A}{e(X)} \mathbb{I}\{C < T^{(1)} \wedge \tau\} \frac{\mathbb{E}[T^{(1)} \wedge \tau \times \mathbb{I}\{C < T^{(1)} \wedge \tau\}|X, A, C]}{\mathbb{P}(C < T^{(1)} \wedge \tau|C, X, A)} \middle| X \right] \\ &= \mathbb{E} \left[ \frac{A}{e(X)} \mathbb{E}[T^{(1)} \wedge \tau \times \mathbb{I}\{C < T^{(1)} \wedge \tau\}|X, A, C] \middle| X \right] \\ &= \mathbb{E} \left[ T^{(1)} \wedge \tau \times \mathbb{I}\{C < T^{(1)} \wedge \tau\} \middle| X \right] \\ &= \mathbb{E} \left[ (1 - \Delta^{\tau}) T^{(1)} \wedge \tau \middle| X \right], \end{aligned} $$

and the same holds on the event $\{A=0\}$, which ends the proof. □
Proof. (Proposition 3.2, Consistency of $\hat{\theta}_{IPTW-BJ}$ with estimated $\hat{e}$ and $\hat{S}$). By consistency of $\hat{G}(\cdot|X, A)$ and $\hat{e}$ and by continuity, it suffices to look at the asymptotic behavior of the oracle estimator

$$ \theta_{IPTW-BJ}^* = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{A_i}{e(X_i)} - \frac{1-A_i}{1-e(X_i)} \right) \left( \Delta_i^\tau \times \tilde{T}_i \wedge \tau + (1-\Delta_i^\tau) Q_S(\tilde{T}_i \wedge \tau | A_i, X_i) \right). $$

The later is converging almost towards its mean, which, following similar computations as in the previous proof, is simply equal to the RMST difference. □

**Proposition 6.11** (Double/Triple robustness of $\Delta_{QR}^*$). (adapted from Ozenne et al. (2020))

Let $F, R$ be two conditional survival function, $p$ be a propensity score, and $\nu$ be a conditional response. Then, under the same assumption on $F, R$ as in Proposition 6.6, and under Assumptions 1.1 (SUTVA), 3.1 (Conditional Exchangeability / Unconfoundedness), 3.2 (Positivity / Overlap for Treatment), 2.5 (Conditionally Independent Censoring) and 2.6 (Conditional Positivity / Overlap for Censoring), the transformations $\Delta_{QR}^* = \Delta_{QR}^*(F, R, p, \nu)$ satisfies, for $a \in 0, 1$,

$$ \mathbb{E}[\Delta_{QR}^*|X] = \mathbb{E}[T^{(1)} \wedge \tau - T^{(0)} \wedge \tau|X] \quad \text{if} \quad \begin{cases} F = G & \text{or } R = S \quad \text{and} \\ p = e & \text{or } \nu = \mu. \end{cases} $$

Proof. (Proposition 6.11). We can write that

$$ \Delta_{QR}^* = \underbrace{\frac{A}{p(X)}(T_{DR}^*(F, R) - \nu(X, 1)) + \nu(X, 1)}_{(A)} - \left( \underbrace{\frac{1-A}{1-p(X)}(T_{DR}^*(F, R) - \nu(X, 0)) + \nu(X, 0)}_{(B)} \right). $$

Focusing on term (A), we easily find that

$$ \begin{aligned} \mathbb{E}[(A)|X] &= \mathbb{E}\left[\frac{A}{p(X)}(T_{DR}^*(F, R) - \nu(X, 1)) + \nu(X, 1)\middle|X\right] \\ &= \frac{e(X)}{p(X)}(\mu(X, 1) - \nu(X, 1)) + \nu(X, 1). \end{aligned} $$

Where we used that $T_{DR}^*(F, R)$ is a censoring-unbiased transform when $F = G$ or $R = S$. Now we see that if $p(X) = e(X)$, then

$$ \mathbb{E}[(A)|X] = \mu(X, 1) - \nu(X, 1) + \nu(X, 1) = \mu(X, 1), $$

and if $\nu(X, 1) = \mu(X, 1)$, then

$$ \mathbb{E}[(A)|X] = \frac{e(X)}{p(X)} \times 0 + \mu(X, 1) = \mu(X, 1). $$

Likewise, we would show that $\mathbb{E}[(B)|X] = \mu(X, 0)$ under either alternative, ending the proof. □
## A.2. Appendix B: Implementations

### A.2.1. Available packages

#### SurvRM2

The RMST difference with Unadjusted Kaplan-Meier $\hat{\theta}_{KM}$ (Equation 3) can be obtained using the function `rmst2` which requires as input the observed time-to-event, the status, the arm which corresponds to the treatment and $\tau$.

#### RISCA

The RISCA package (Foucher et al., 2019) provides several methods for estimating $\theta_{RMST}$. The RMST difference with Unadjusted Kaplan-Meier $\hat{\theta}_{KM}$ (Equation 3) can be derived using the `survfit` function from the `survival` package (Therneau, 2001) which estimates Kaplan-Meier survival curves for treated and control groups, and then the `rmst` function calculates the RMST by integrating these curves, applying the rectangle method (type="s"), which is well-suited for step functions.

The IPTW Kaplan-Meier (Equation 18) can be applied using the `ipw.survival` and `rmst` functions. The `ipw.survival` function requires user-specified weights (i.e. propensity scores). To streamline this process, we define the `RISCA_iptw` function, which combines these steps and utilizes the `estimatepropensity_score` from the `utilitary.R` file.

A single-learner version of the G-formula, as introduced in Section 2.2.1, can be implemented using the `gc.survival` function. This function requires as input the conditional survival function which should be estimated beforehand with a Cox model via the `coxph` function from the `survival` package (Therneau, 2001). Specifically, the single-learner approach applies a single Cox model incorporating both covariates and treatment, rather than separate models for each treatment arm. We provide a function `RISCA_gf` that consolidates these steps.

#### grf

The `grf` package (Tibshirani et al., 2017) enables estimation of the difference between RMST using the Causal Survival Forest approach (Cui et al., 2023), which extends the non-parametric causal forest framework to survival data. The RMST can be estimated with the `causal_survival_forest` function, requiring covariates $X$, observed event times, event status, treatment assignment, and the time horizon $\tau$ as inputs. The `average_treatment_effect` function then evaluates the treatment effect based on predictions from the fitted forest.

### A.2.2. Implementation of the estimators

The complete codebase is available on github.com/Sanofi-Public/causal_survival_analysis

#### Unadjusted Kaplan-Meier

Although Kaplan-Meier is implemented in the `survival` package (Therneau, 2001), we provide a custom implementation, `Kaplan_meier_handmade`, for completeness. The difference in
Restricted Mean Survival Time, estimated using Kaplan-Meier as in Equation 3 can then be
calculated with the RMST_1 function.

As an alternative, one can also use the `survfit` function in the `survival` package (Therneau, 2001) for Kaplan-Meier and specify the `rmean` argument equal to $\tau$ in the corresponding summary function:

## IPCW Kaplan-Meier

We first provide an adjusted.KM function which is then used in the IPCW_Kaplan_meier function to estimate the RMST difference $\hat{\theta}_{IPCW}$ as in Equation 12. The survival censoring function $G(t|X)$ is computed with the estimate_survival_function utility function from the utilitary.R file.

One could also use the `survfit` function in the `survival` package (Therneau, 2001) by adding
IPCW weights for treated and control group and specify the `rmean` argument equal to $\tau$ in
the corresponding summary function:

This alternative approach for IPCW Kaplan-Meier would also be valid for IPTW and
IPTW-IPCW Kaplan-Meier.

## Buckley-James-based estimator

The function BJ estimates $\theta_{RMST}$ by implementing the Buckley-James estimator as in Equation 14. It uses two functions available in the utilitary.R file, namely Q_t_hat and Q_Y.

## IPTW Kaplan-Meier

The function IPTW_Kaplan_meier implements the IPTW-KM estimator in Equation 19. It uses the estimatePropensity_score function from the utilitary.R.

## G-formula

We implement two versions of the G-formula: g_formula_T_learner and g_formula_S_learner.
In g_formula_T_learner, separate models estimate survival curves for treated and control
groups, whereas g_formula_S_learner uses a single model incorporating both covariates
and treatment status to estimate survival time. The latter approach is also available in the
RISCA package (Foucher et al., 2019) but is limited to Cox models.

## IPTW-IPCW Kaplan-Meier

The IPTW_IPCW_Kaplan_meier function implements the IPTW-IPCW Kaplan Meier estimator from Equation 21. It uses the utilitary functions from the utilitary.R file estimatePropensity_score and estimate_survival_function to estimate the nuisance parameters, and the function adjusted.KM which computes an adjusted Kaplan Meier estimator using the appropriate weight.

## IPTW-BJ estimator

The IPTW_BJ implements the IPTW-BJ estimator in Equation 22. It uses the utilitary functions, from the utilitary.R file, estimatePropensity_score, Q_t_hat and Q_Y to estimate the nuisance parameters.
## AIPTW-AIPCW

The AIPTW\_AIPCW function implements the AIPTW\_AIPCW estimator in Equation 24 using the utilitary function from the utilitary.R file estimatePropensity_score, Q_t_hat, Q_Y, and estimateSurvival_function to estimate the nuisance parameters.

### A.3. Appendix C: Data Generating Process and Figures

**Simulation coefficients (Section 5.2).** In the RCT setting, we set $\mu = (1, 1, -1, 1)$, $\Sigma = \text{Id}_4$, and use a Cox baseline rate $\lambda_0 = 0.01$ with $\beta_0 = (0.5, 0.5, -0.5, 0.5)$ for $\lambda^{(0)}(t | X) = \lambda_0 \exp\{\beta_0^\top X\}$. Censoring uses $\lambda_{C,0} = 0.03$: in Scenario 1, $\lambda_C(t) = \lambda_{C,0}$; in Scenario 2, $\lambda_C(t | X) = \lambda_{C,0} \exp\{\beta_C^\top X\}$ with $\beta_C = (0.7, 0.3, -0.25, -0.1)$. Treatment is randomized with $e(X) = 0.5$.

<figure><img src="image_11.png" /><figcaption>Figure 11: Results of the ATE for the simulation of a RCT with independent censoring.</figcaption></figure>

**Simulation coefficients (Section 5.3).** In the observational setting, we use a logistic propensity model $\text{logit}(e(X)) = \beta_A^\top X$ with $\beta_A = (-1, -1, -2.5, -1)$. Recall that $\text{logit}(p) = \log(p/(1-p))$. All other components of the data-generating process are identical to those in Section 5.2.
<figure><img src="image_12.png" /><figcaption>Figure 12: Estimation results of the ATE for the simulation of a RCT with dependent censoring.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: Estimation results of the ATE for the simulation of an observational study with independent censoring.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: Estimation results of the ATE for the simulation of an observational study with dependent censoring.</figcaption></figure>

**Simulation coefficients (Section 5.4).** We set $\mu = (0.1, 0.5, 0.7, 0.4)$ and $\Sigma = \text{Id}_4$. Let $Z(X) = (X_1^2, X_2^2, X_3^2, X_4^2, X_1X_2, X_1X_3, X_1X_4, X_2X_3, X_2X_4, X_3X_4)$. The control hazard uses

$\beta_0 = (0.2, 0.3, 0.1, 0.1, 0, 1, 0, 0, 0, -1)$. Censoring uses $\beta_C = (0.05, -0.01, 0.05, -0.01, 0, 1, 0, 0, -1, 0)$.

The propensity model uses $\beta_A = (0.05, -0.1, 0.5, -0.1, -1, 0, -1, 0, 0, 0)$.

Under well specification, nuisance models use inputs $(X, Z(X))$; under misspecification, they use only $X$ and the first four components of $Z(X)$, i.e., $(X_1^2, X_2^2, X_3^2, X_4^2)$.
<figure><img src="image_15.png" /><figcaption>Figure 15: Estimation results of the ATE for an observational study with dependent censoring in case of a single misspecification.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 16: ATEs estimation in observational study with conditionally independent censoring with two or more misspecifications.</figcaption></figure>