# Bayesian analysis of restricted mean survival time adjusted for covariates using pseudo-observations

Léa Orsini, Emmanuel Lesaffre, Guosheng Yin, Caroline Brard, David

Dejardin, Gwénaël Le Teuff

► To cite this version:

Léa Orsini, Emmanuel Lesaffre, Guosheng Yin, Caroline Brard, David Dejardin, et al. Bayesian analysis of restricted mean survival time adjusted for covariates using pseudo-observations. 2025. [hal-04975637v3](https://hal.science/hal-04975637v3)

HAL Id: hal-04975637

https://hal.science/hal-04975637v3

Preprint submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
BAYESIAN ANALYSIS OF RESTRICTED MEAN SURVIVAL TIME
ADJUSTED FOR COVARIATES USING PSEUDO-OBSERVATIONS

Léa Orsini

Oncostat U1018, Inserm
University Paris-Saclay
Villejuif, France
lea.orsini@gustaveroussy.fr

Emmanuel Lesaffre

I-Biostat
KU-Leuven
Leuven, Belgium
emmanuel.lesaffre@kuleuven.be

Guosheng Yin

Department of Statistics and Actuarial Science
The University of Hong Kong
Hong Kong, China
gyin@hku.hk

Caroline Brard

Ipsen Innovation, Clinical Development Organisation
Les Ulis, France
caroline.brard@ipsen.com

David Dejardin

Product Development, Data Sciences
F. Hoffmann-La Roche AG
Basel, Switzerland.
david.dejardin@roche.com

Gwénaël Le Teuff

Oncostat U1018, Inserm
University Paris-Saclay
Villejuif, France
gwenael.leteuff@gustaveroussy.fr

ABSTRACT

The difference in restricted mean survival time (RMST) is a clinically meaningful measure to quantify
treatment effect in randomized controlled trials, especially when the proportional hazards assumption
does not hold. Several frequentist methods exist to estimate RMST adjusted for covariates based
on modeling and integrating the survival function. A more natural approach may be a regression
model on RMST using pseudo-observations, which allows for a direct estimation without modeling
the survival function. Only a few Bayesian methods exist, and each requires a model of the survival
function. We developed a new Bayesian method that combines the use of pseudo-observations with
the generalized method of moments. This offers RMST estimation adjusted for covariates without
the need to model the survival function, making it more attractive than existing Bayesian methods.
A simulation study was conducted with different time-dependent treatment effects (early, delayed,
and crossing survival) and covariate effects, showing that our approach provides valid results, aligns
with existing methods, and shows improved precision after covariate adjustment. For illustration, we
applied our approach to a phase III trial in prostate cancer, providing estimates of the treatment effect
on RMST, comparable to existing methods. In addition, our approach provided the effect of other
covariates on RMST and determined the posterior probability of the difference in RMST exceeds any
given time threshold for any covariate, allowing for nuanced and interpretable results.

**Keywords** Bayesian survival analysis · Generalized method of moments · Non-proportional hazards · Pseudo-observations · Restricted mean survival time
# 1 Introduction

The difference in restricted mean survival time (dRMST) has been proposed as a clinically meaningful measure of the treatment effect of time-to-event outcomes in randomized controlled trials (RCTs), especially when the proportional hazards (PH) assumption does not hold since it does not rely on the PH assumption (Royston and Parmar, 2011; Corro Ramos, Qendri, and Al, 2024) and provides a meaningful interpretation of the results (Pak et al., 2017). Non-PH occurs when the treatment effect is time-dependent (e.g., early or delayed effect). In recent years, a delayed treatment effect has often been observed when comparing immunotherapies with chemotherapies (Anagnostou et al., 2017). Non-PH may also occur with therapies that have different mechanisms of action. For example, treatment effects may be different between biomarker-positive and biomarker-negative groups (predictive biomarker), resulting in strongly different survival patterns. Even if the hazards are proportional within each group, non-PH can occur when considering the whole population, see for example Mok et al. (2009). In these two settings, interpreting the hazard ratio as a treatment effect measure becomes tedious (Hernán, 2010).

Extensive literature has been devoted to this measure (Karrison, 1997; Royston and Parmar, 2013; Hasegawa et al., 2020). It is interpreted as the expected lifetime experienced out of $τ$ units of time (Irwin, 1949). Suppose, for example, that dRMST at $τ = 5$ years between an experimental arm and a control arm in an RCT with an overall survival endpoint is estimated to be 1 year. This means that, on average, the experimental treatment increases lifetime expectancy over the next 5 years by 1 year, compared to the control arm.

One straightforward approach for estimating RMST at a certain time $τ$ is to numerically integrate the Kaplan-Meier curve between 0 and $τ$. However, this approach does not allow for covariate adjustment, which is a major limitation because omitting important covariates may result in less precision (Karrison and Kocherginsky, 2018). One way to adjust RMST estimation for covariates is to model the survival function with a parametric or semi-parametric model and integrate it, see Karrison (1987) and Zucker (1998). A more natural approach is directly fitting a regression model on RMST through estimating equations. In this case, censoring must be handled using either the inverse probability of censoring weights (Tian, Zhao, and Wei, 2014) or pseudo-observations (Andersen, Hansen, and Klein, 2004).

RMST estimation can also be performed with a Bayesian approach, with a particular interest in rare diseases where small sample sizes are prevalent, and frequentist inference relying on asymptotic results may not be reliable (Lesaffre, Baio, and Boulanger, 2020). The Bayesian approach naturally allows the incorporation of prior information and may help to fit complex models. In addition, the Bayesian approach provides a more intuitive interpretation of the results by estimating the whole posterior distribution, giving an additional advantage compared to the frequentist RMST estimator. However, Bayesian research on RMST is limited. Recently, Zhang and Yin (2023) proposed a Bayesian nonparametric approach to derive the RMST distribution from posterior samples for right and interval-censored data by assigning a Mixture of Dirichlet Processes (MDP) prior to the cumulative distribution function. However, this approach does not allow for covariate adjustment. Chen, Basu, and Shi (2014) overcome this limitation by considering a data-adaptive stick-breaking prior based on a probit regression model. Hanada and Kojima (2024) estimate the RMST posterior distribution by fitting a parametric model on the survival function accounting for covariates and heterogeneity among clusters. All these methods require first modeling the survival function (or related functions) and then integrating them out in the two groups independently before subtracting the results, providing an indirect way of estimating the dRMST.

In this paper, we extend the frequentist method of Andersen et al. (2004) based on the pseudo-observations to the Bayesian framework, offering a direct multivariable dRMST estimation without the need to model the survival function, thereby providing an attractive alternative to existing Bayesian methods.

The rest of the paper is organized as follows. Section 2 presents the Bayesian modeling of RMST using pseudo-observations with the Generalized Method of Moments (GMM). Section 3 compares the performance of the GMM model to benchmark methods through a simulation study representing multiple scenarios of non-PH. We illustrate the methods through the analysis of a phase III RCT, the Getug-AFU 15 trial involving patients with prostate cancer, in Section 4. We conclude with some final remarks and future extensions of the proposed approach in Section 5.
# 2 Methods

## 2.1 RMST definition

Let $\tilde{T}_i$ be the unobserved event time of interest ($i = 1, ..., n$) for the i-th subject, $Z_i$ a p-dimensional baseline covariate vector, $A_i$ the treatment allocation variable of an RCT ($A_i=1$ and $0$ for experimental and control arm, respectively), and $C_i$ a right censoring random variable, independent of $\tilde{T}_i$, $Z_i$ and $A_i$. We observe $T_i = \min(\tilde{T}_i, C_i)$ and $\Delta_i = I(\tilde{T}_i \le C_i)$ the event indicator. For a pre-specified time point of interest $\tau$, the $\tau$-RMST is defined as

$$ \text{RMST}(\tau) = E(\min(\tilde{T}, \tau)) = \int_{0}^{\tau} S(t)dt. \qquad (1) $$

Under the specification of a generalized linear model, the subject-specific RMST at $\tau$ is defined by $\text{RMST}_i(\tau) = E(\min(\tilde{T}_i, \tau))$ and can be conveniently modeled as

$$ \mu_i = E(\min(\tilde{T}_i, \tau)|A_i, Z_i) = g^{-1}(\beta_0 + \delta A_i + \beta_1 Z_{i1} + \dots + \beta_p Z_{ip}), \qquad (2) $$

where $g(\cdot)$ is a monotone differentiable link function and $\beta = (\beta_0, \delta, \beta_1, ..., \beta_p)^\top$ the vector of unknown parameters. For the identity link function, the regression coefficient $\delta$ can be interpreted as the dRMST between the two arms of an RCT.

## 2.2 Regression modeling of RMST using pseudo-observations

In the frequentist framework, regression model (2) can be fitted using the pseudo-observations approach, see Andersen, Klein, and Rosthøj (2003). Following Andersen et al. (2004), the i-th pseudo-observation is computed as

$$ y_{\tau,i} = n \int_{0}^{\tau} \hat{S}(t)dt - (n-1) \int_{0}^{\tau} \hat{S}^{-i}(t)dt, \qquad (3) $$

with n the sample size, $\hat{S}(t)$ the Kaplan-Meier (KM) estimator of the survival probability, and $\hat{S}^{-i}(t)$ the KM estimator excluding the i-th subject. Given the asymptotic properties of pseudo-observations proved in Overgaard, Parner, and Pedersen (2017), we can replace the partially observed (due to censoring) $\min(\tilde{T}_i, \tau)$ by $y_{\tau,i}$ in regression model (2). Consequently, the pseudo-observations can be treated as outcome variables in a regression model to regress $\tau$-RMST on explanatory variables, using generalized estimating equations (GEE) (Liang and Zeger, 1986). This marginal approach is based on quasi-likelihood functions where only the moments are defined (McCullagh and Nelder, 1991).

The mean model is specified as:

$$ \mu_i = E(y_{\tau,i} | A_i, Z_i) = g^{-1}(\beta_0 + \delta A_i + \beta_1 Z_{i1} + \dots + \beta_p Z_{ip}), \qquad (4) $$

The vector $\beta$ is estimated by solving the score equations:

$$ U_n(\beta) = \frac{1}{n} \sum_{i=1}^{n} u_i(\beta) = \frac{1}{n} \sum_{i=1}^{n} \frac{\partial \mu_i}{\partial \beta} (y_{\tau,i} - \mu_i) = 0, \qquad (5) $$

where $\frac{\partial \mu_i}{\partial \beta}$ is the first derivative of the mean model $\mu_i$ with respect to $\beta$. In this case, the link function is the identity function.

## 2.3 Bayesian generalized method of moments

We propose to analyze pseudo-observations with a generalized linear model using the Generalized Method of Moments (GMM), defined by Hansen (1982) in the frequentist framework and derived by Yin (2009) in the Bayesian framework. The GMM is appropriate as it is a generic way of estimating parameters in statistical models, which is advantageous when the likelihood of the data is cannot be specified. The GMM does not require assuming a specific distribution for the data but is based on a number of moment conditions (first-order moment, and possibly higher-order moments). These moment conditions are combined through equations linking the model parameters and the data. If the number of equations exceeds the number of model parameters, the system of equations is over-identified, so the estimates are then found by minimizing an objective function.
To estimate RMST with the Bayesian GMM, with RMST model (2), pseudo-observations are computed as equation (3). Here, we are specifying only the first-order moment as in (4). Therefore, the system of equations (5) holds. In this context, the objective function can be written as

$$Q_n(\beta) = U_n^{\top}(\beta)\Sigma_n^{-1}(\beta)U_n(\beta), \qquad (6)$$

$$\text{where } \Sigma_n(\beta) = \frac{1}{n^2} \sum_{i=1}^{n} u_i(\beta)u_i^{\top}(\beta) - \frac{1}{n} U_n(\beta)U_n^{\top}(\beta).$$

By applying the Central Limit Theorem,

$$U_n(\beta) \xrightarrow{d} N(0, \Sigma(\beta)), \text{ as } n \to \infty \qquad (7)$$

where $\Sigma(\beta) = \lim_{n \to \infty} \Sigma_n(\beta)$, then

$$Q_n(\beta) \xrightarrow{d} \chi_{p+2}^2. \qquad (8)$$

A chi-squared test can be derived where $Q_n(\beta)$ behaves like $-2 \log L(\beta|y)$ with $L(\beta|y)$ the likelihood function (Hansen, 1982). Thus, the objective function can be seen as an approximation of the likelihood for selected moments of the data when the full likelihood of the data is unknown, (Chernozhukov and Hong, 2003). A Bayesian sampling procedure can be developed where the likelihood function is approximated by defining a pseudo-likelihood function as follows:

$$\tilde{L}(\beta|y) \propto \exp\{-\frac{1}{2} U_n^{\top}(\beta) \Sigma_n^{-1}(\beta) U_n(\beta)\}. \qquad (9)$$

Yin (2009) demonstrated the validity of the posterior distribution sampled with this pseudo-likelihood. The posterior distribution for the regression coefficient, $\delta$, represents the posterior distribution for the dRMST between the two treatment arms.

# 3 Simulation study

A simulation study of 2-arms RCTs was conducted to assess the performance of the new Bayesian approach based on pseudo-observations for estimating RMST. The purpose of this simulation study was to (a) assess the validity of the dRMST estimation using the Bayesian GMM based on pseudo-observations, (b) compare these estimates with other frequentist and Bayesian dRMST estimates (benchmark methods), and (c) evaluate the impact of covariates adjustment and model misspecification.

## 3.1 Simulation settings

We simulated individual survival data of a 2-arm RCT (ratio 1:1), for different scenarios of the treatment effect (experimental vs control) over time with PH (scenario 1) and non-PH (scenarios 2-6). Figure 1 shows the true survival curves for each scenario. In scenarios 2 and 4, we simulate an early treatment effect, in scenarios 3 and 5, we simulate a delayed treatment effect, and in scenario 6, we simulate crossing survival curves. In scenarios 4 and 5 we simulate survival data associated with prognostic variables. These two scenarios allow us to evaluate the impact of covariate adjustments and model misspecification when omitting important prognostic variables. For each scenario, the event times were simulated according to a Weibull distribution $S(t) = \exp(-(\lambda t)^{1/\sigma})$, with scale $\lambda$ and shape $\sigma$ parameters chosen to mimic different patterns of treatment effect. The corresponding true value of the RMST is

$$\int_{0}^{\tau} \exp(-(\lambda t)^{1/\sigma}) dt = \frac{\sigma}{\lambda} \gamma(\sigma, (\lambda \tau)^{1/\sigma}), \qquad (10)$$

where $\gamma(a, x) = \int_{0}^{x} t^{a-1}e^{-t}dt$ is the lower incomplete gamma function. The true value of the dRMST ($\delta$) can be calculated by subtracting the true RMST values of the two arms. The censoring times were drawn from a uniform distribution, with an additional administrative censoring at 8 years, yielding an average of 30% of censored individuals per scenario. Prognostic covariates $Z_{ij}$ were drawn from a uniform distribution (scenario 4) or normal and Bernoulli distributions (scenario 5). Their effects on time-to-event satisfy the proportional hazards assumption. Other variables $X_{ij}$, unrelated to the survival function, were also simulated from a normal, a Bernoulli, and a uniform distribution, respectively, to assess the impact of covariate adjustments in terms of model misspecification in scenarios 4 and 5. In scenario 6, crossing survival curves between the experimental and control arms were simulated by including an interaction between the treatment indicator and a binary variable $E_i$ with a predictive value. Consequently, the treatment effect depends on whether the biomarker takes $E_i = E^{-}$ or $E_i = E^{+}$. We assume that the treatment was harmful if $E_i = E^{-}$, and the treatment was protective if $E_i = E^{+}$. The different simulated scenarios are summarized in supporting
information Table 1. For simplicity, we will omit the i index in $A_i$, $Z_{ij}$, $X_{ij}$, and $E_i$ throughout the rest of this paper. Different sample sizes of simulated RCTs were considered, i.e. 50, 100, 200, and 500. According to the follow-up of the simulated RCTs (accrual period and administrative censoring), RMST was estimated at the restriction time $\tau$ set to 5 years, often used as a clinically meaningful threshold in oncology trials (Chen et al., 2014). However, when the last simulated observed time in one arm was smaller than 5 years, then $\tau$ was redefined as the minimum of the maximum observed times of each arm. The number of times that this situation occurs is reported in supporting information Table 2.

<figure><img src="image_2.png" /><figcaption>Figure 1: True survival curves for the different scenarios.</figcaption></figure>

The performance of the Bayesian GMM and the benchmark methods (See subsection 3.2 for a description of these methods) was evaluated using the following metrics calculated from 1000 replications: bias, average standard errors (ASE: the mean of the estimated standard errors of $\hat{\delta}$ across replications), empirical standard errors (ESE: the standard deviation of $\hat{\delta}$ across replications), root mean square error (RMSE = $\sqrt{ESE^2 + Bias^2}$), and the coverage rate of the intervals (either 95% confidence intervals or equal-tailed credible intervals for the frequentist and the Bayesian approaches, respectively). The convergence of the Bayesian procedures was checked with trace plots checking and Gelman-Rubin $\hat{R}$ diagnostic (Vehtari et al., 2021).

All computations were performed using the R Language for Statistical Computing (R Core Team (2021), version 4.1.2). Pseudo-observations have been computed using the R package pseudo (Pohar-Perme, Gerster, and Rodrigues, 2017). We developed a specific script to implement the Bayesian GMM using the Stan software (Carpenter et al., 2017). The Hamiltonian Monte Carlo algorithm was performed using the No-U-turn (NUTS) sampler with 3 chains of each 1000 iterations after a warm-up of 1000 iterations, yielding 3000 iterations overall. Weakly informative priors were specified using a normal distribution $N(\mu = 0, \sigma^2 = 10)$ for all parameters of the Bayesian GMM with pseudo-observations, as recommended by the Stan Development Team (2020). Regarding computational time, for one run of the Bayesian GMM with $n = 200$ patients and 3 explanatory variables included in the model, the pseudo-observations computation took less than 1 second, and the running time of one chain for the Bayesian GMM was 10 seconds with an 11th Gen Intel(R) Core(TM) i5-1135G7 processor. Using the rstan R package (Stan Development Team, 2023), chains can run in parallel, not adding additional running time when increasing the number of chains. So, for the simulation study, the Bayesian GMM was run with a parallelized code using a server HPE DL385 (2.0 GHz) with 150 virtual cores due to the high numbers of scenarios and replicates per scenario. The median running time of one chain according to different sample sizes and numbers of explanatory variables included in the model is provided in supporting information Figure 1 for scenario 5. Regarding Bayesian GMM convergence, only scenarios 4 and 5, which involve covariate adjustment with prognostic and unrelated variables, showed a small percentage of non-convergent replications (less than 1.5%, supporting information Table 3), which were subsequently removed from the analysis.
## 3.2 Benchmark methods

Four benchmark methods, providing an estimator of dRMST, have been retained for comparison: three frequentist estimators and one Bayesian. The first one (naive approach) consists of integrating the Kaplan-Meier estimator in both arms separately before subtracting the estimated RMST values of the two arms. The second consists of applying a Cox model stratified on treatment (Cox, 1972), with Breslow (1972) estimator to estimate the baseline hazard function in each arm and model the survival function conditional on covariates before averaging over the covariates, integrating it and subtracting the estimated RMST values of the two arms, see Zucker (1998) and the third benchmark method is the pseudo-observations approach from Andersen et al. (2004) that proposed GEE to model RMST. The R package *survival* was used to estimate the RMST from the Kaplan-Meier estimator (Therneau, 2024). The R package *geepack* was used to implement the GEE approach for fitting marginal generalized linear model to pseudo-observations (Højsgaard et al., 2022). The Bayesian benchmark method is the one of Zhang and Yin (2023) who proposed a nonparametric approach for estimating the RMST by independently assigning an MDP prior to the cumulative incidence function of each arm. An MDP prior was applied under an exponential base measure with a Gamma $\Gamma(0.01, 0.01)$ mixing distribution.

## 3.3 Unadjusted estimation of the RMST

In scenario 1, the simulation parameters were set to $\sigma = 0.8$ and $\lambda = \exp(-1.2 + \log(HR)A_i)$ expressing the treatment effect as a hazard ratio HR. Under the null hypothesis with proportional hazards (no treatment effect, HR = 1 thus $\delta = 0$), the Bayesian GMM produced similar estimations of the dRMST between the two arms compared to the four benchmark methods: no bias, ASE and ESE were close and decrease when sample size increased and coverage was close to the nominal 95% rate when the sample size was large enough ($n \ge 100$) (see supporting information Table 4). More generally, the performance of the different methods became similar when the sample size increased. Under an alternative hypothesis ($\sigma = 0.8$ and the treatment effect HR = 0.6, resulting in $\delta = 0.82$ when applying the formula (10)), again the results were similar to the null hypothesis setting (see supporting information Table 5).

In scenarios 2 and 3, we simulated a treatment effect with non-proportional hazards by specifying different scale and shape parameters for each treatment group. In scenario 2, we set $\sigma_0 = 1.33$, $\lambda_0 = 0.20$ for the control group and $\sigma_1 = 0.67$, and $\lambda_1 = 0.18$ for the treatment group, resulting in survival curves presenting an early treatment effect that decreases through the years. In scenario 3, we set $\sigma_0 = 0.60$, $\lambda_0 = 0.28$ for the control group and $\sigma_1 = 0.80$, $\lambda_1 = 0.18$ for the treatment group, resulting in survival curves presenting a delayed treatment effect that increases through the years. Table 1 reports the results for these 2 scenarios with sample sizes of n=50, 100, and 200 (See supporting information Table 6 and supporting information Table 7 for an additional sample size of n=500). The findings were similar to what we observed for scenario 1, i.e., the Bayesian GMM gave a valid estimate of the dRMST for any sample size with performance closer to those of the four benchmark methods as the sample size increases. For large sample sizes (n = 500), the RMSE closely approximated the ESE. Conversely, for n = 50, in scenarios 1-3, we note that the Bayesian GMM was more conservative with a higher ASE compared to the four benchmark methods, resulting in a coverage rate closer to the 95% and around 92%-93% for the other methods.

## 3.4 Adjusted estimation of the RMST

In scenarios 4 and 5, we simulated situations where RMST was associated with one prognostic variable $Z_1$ (scenario 4) and two prognostic variables $Z_1$ and $Z_2$ (scenario 5). In scenario 4, we simulated an early treatment effect with $Z_1 \sim U(0, 2)$. In scenario 5, we simulated a delayed treatment effect with $Z_1 \sim N(0, 1)$ and $Z_2 \sim Bin(0.5)$. Table 2 reports the results with and without covariates-adjustment for an RCT of n = 200. This corresponds to the first two and the first three rows of each method (Zucker (1998); Andersen et al. (2004) and Bayesian GMM) for scenario 4 and scenario 5, respectively. All methods allowing for adjustment for covariates produced slightly more precise estimates when the model was well specified compared to the misspecified univariate analysis with a relative gain in precision from unadjusted to adjusted estimations of 4% (scenario 4) and 8% (scenario 5) for the Bayesian GMM. Omitting prognostic variables also produced estimates with less precision with n = 100 and n = 500, see supporting information Table 8. For a smaller sample size (n = 50), we did not observe an increase in precision with the Bayesian GMM when adjusting for the prognostic variables, while the frequentist models of Zucker (1998) and Andersen et al. (2004) still provided smaller ASE and satisfactory coverage. This behavior was still reasonably small when adjusting for one covariate: in scenario 4, an ASE of 0.517 (94.5% coverage rate) for the unadjusted model compared to 0.532 (95.5% coverage rate) when adjusting for $Z_1$. When we adjusted for two covariates, the variance increased more: in scenario 5, ASE of 0.517 (96.2% coverage rate) for the unadjusted model compared to an ASE of 0.742 (99.1% coverage rate) when adjusting for $Z_1$ and $Z_2$. On the other hand, the ESE of the Bayesian GMM stayed comparable to the ones of the other approaches, see supporting information Table 9.
<table><caption>Table 1: Performance of the frequentist and Bayesian methods on the univariate estimation of the difference of restricted mean survival time (dRMST) between two treatment groups with non-proportional hazards treatment effect (scenario 2: early effect $\delta = 0.7302$, scenario 3: late effect $\delta = 0.5644$), with 30% of censoring and different sample sizes.</caption><thead><tr><th>Scenario</th><th>n</th><th>Methods</th><th>Bias</th><th>ASE<sup>1</sup></th><th>ESE<sup>2</sup></th><th>RMSE<sup>3</sup></th><th>Coverage</th></tr></thead><tbody><tr><td rowspan="12">2</td><td rowspan="12">50</td><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>-0.0056</td><td>0.498</td><td>0.538</td><td>0.538</td><td>92.2</td></tr><tr><td>Zucker (1998)</td><td>-0.0198</td><td>0.493</td><td>0.526</td><td>0.526</td><td>92.3</td></tr><tr><td>Andersen et al. (2004)</td><td>-0.0058</td><td>0.499</td><td>0.538</td><td>0.538</td><td>92.1</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>-0.0061</td><td>0.487</td><td>0.538</td><td>0.538</td><td>92.1</td></tr><tr><td>GMM</td><td>0.0213</td><td>0.514</td><td>0.521</td><td>0.522</td><td>93.2</td></tr><tr><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>0.0068</td><td>0.356</td><td>0.356</td><td>0.356</td><td>94.3</td></tr><tr><td>Zucker (1998)</td><td>-0.0006</td><td>0.354</td><td>0.352</td><td>0.352</td><td>94.6</td></tr><tr><td>Andersen et al. (2004)</td><td>0.0068</td><td>0.356</td><td>0.356</td><td>0.356</td><td>94.4</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>0.0067</td><td>0.352</td><td>0.356</td><td>0.356</td><td>94.4</td></tr><tr><td>GMM</td><td>0.0192</td><td>0.359</td><td>0.351</td><td>0.351</td><td>94.8</td></tr><tr><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>-0.0071</td><td>0.252</td><td>0.254</td><td>0.254</td><td>95.4</td></tr><tr><td>Zucker (1998)</td><td>-0.0107</td><td>0.251</td><td>0.252</td><td>0.253</td><td>95.3</td></tr><tr><td>Andersen et al. (2004)</td><td>-0.0070</td><td>0.252</td><td>0.254</td><td>0.254</td><td>95.4</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>-0.0072</td><td>0.251</td><td>0.254</td><td>0.254</td><td>95.3</td></tr><tr><td>GMM</td><td>-0.0020</td><td>0.252</td><td>0.252</td><td>0.252</td><td>95.4</td></tr><tr><td rowspan="12">3</td><td rowspan="12">50</td><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>-0.0212</td><td>0.461</td><td>0.491</td><td>0.492</td><td>92.7</td></tr><tr><td>Zucker (1998)</td><td>-0.0278</td><td>0.459</td><td>0.485</td><td>0.486</td><td>93.1</td></tr><tr><td>Andersen et al. (2004)</td><td>-0.0206</td><td>0.464</td><td>0.492</td><td>0.492</td><td>93.1</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>-0.0129</td><td>0.451</td><td>0.497</td><td>0.497</td><td>92.3</td></tr><tr><td>GMM</td><td>-0.0044</td><td>0.482</td><td>0.480</td><td>0.480</td><td>94.4</td></tr><tr><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>0.0058</td><td>0.330</td><td>0.340</td><td>0.340</td><td>93.5</td></tr><tr><td>Zucker (1998)</td><td>-0.0016</td><td>0.329</td><td>0.336</td><td>0.336</td><td>93.9</td></tr><tr><td>Andersen et al. (2004)</td><td>0.0055</td><td>0.331</td><td>0.339</td><td>0.339</td><td>93.7</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>0.0059</td><td>0.326</td><td>0.340</td><td>0.340</td><td>93.4</td></tr><tr><td>GMM</td><td>0.0124</td><td>0.336</td><td>0.335</td><td>0.335</td><td>94.3</td></tr><tr><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>-0.0065</td><td>0.234</td><td>0.237</td><td>0.237</td><td>93.9</td></tr><tr><td>Zucker (1998)</td><td>-0.0101</td><td>0.233</td><td>0.235</td><td>0.236</td><td>94.0</td></tr><tr><td>Andersen et al. (2004)</td><td>-0.0063</td><td>0.235</td><td>0.237</td><td>0.237</td><td>93.9</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>-0.0065</td><td>0.232</td><td>0.237</td><td>0.237</td><td>94.1</td></tr><tr><td>GMM</td><td>-0.0035</td><td>0.236</td><td>0.236</td><td>0.236</td><td>95.1</td></tr></tbody></table>

<footnote n="1">¹ASE: Average Standard Error, ²ESE: Empirical Standard Error, ³RMSE: Root Mean Square Error</footnote>
<table><caption>Table 2: Performance of the frequentist and Bayesian methods on the estimation of the difference of restricted mean survival time (dRMST) between two treatment groups, adjusted for prognostic covariates Z_j and non-prognostic covariates X_j, with non-proportional hazards treatment effect (scenario 4: early effect δ = 0.9532, scenario 5: late effect δ = 0.4911), with 30% of censoring and a sample size of 200. The gray shaded lines represent the correctly specified model.</caption><thead><tr><th>Scenario</th><th>Methods</th><th>Adjustment variables</th><th>Bias</th><th>ASE<sup>1</sup></th><th>ESE<sup>2</sup></th><th>RMSE<sup>3</sup></th><th>Coverage</th></tr></thead><tbody><tr><td rowspan="12">4</td><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>-</td><td>-0.0056</td><td>0.257</td><td>0.266</td><td>0.266</td><td>93.8</td></tr><tr><td>Zucker (1998)</td><td>-</td><td>-0.0104</td><td>0.258</td><td>0.264</td><td>0.264</td><td>93.8</td></tr><tr><td>Zucker (1998)</td><td>Z<sub>1</sub></td><td>-0.0133</td><td>0.239</td><td>0.243</td><td>0.243</td><td>93.9</td></tr><tr><td>Zucker (1998)</td><td>Z<sub>1</sub>, X<sub>1</sub></td><td>-0.0128</td><td>0.239</td><td>0.245</td><td>0.245</td><td>93.9</td></tr><tr><td>Zucker (1998)</td><td>Z<sub>1</sub>, X<sub>1</sub>, X<sub>2</sub>, X<sub>3</sub></td><td>-0.0117</td><td>0.239</td><td>0.248</td><td>0.248</td><td>93.8</td></tr><tr><td>Andersen et al. (2004)</td><td>-</td><td>-0.0056</td><td>0.258</td><td>0.266</td><td>0.266</td><td>93.8</td></tr><tr><td>Andersen et al. (2004)</td><td>Z<sub>1</sub></td><td>-0.0088</td><td>0.246</td><td>0.251</td><td>0.251</td><td>93.9</td></tr><tr><td>Andersen et al. (2004)</td><td>Z<sub>1</sub>, X<sub>1</sub></td><td>-0.0074</td><td>0.246</td><td>0.251</td><td>0.251</td><td>93.9</td></tr><tr><td>Andersen et al. (2004)</td><td>Z<sub>1</sub>, X<sub>1</sub>, X<sub>2</sub>, X<sub>3</sub></td><td>-0.0068</td><td>0.246</td><td>0.254</td><td>0.254</td><td>93.7</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>-</td><td>-0.0058</td><td>0.256</td><td>0.266</td><td>0.266</td><td>93.8</td></tr><tr><td>GMM</td><td>-</td><td>-0.0070</td><td>0.259</td><td>0.263</td><td>0.263</td><td>94.5</td></tr><tr><td>GMM</td><td>Z<sub>1</sub></td><td>-0.0033</td><td>0.250</td><td>0.249</td><td>0.249</td><td>94.6</td></tr><tr><td>GMM</td><td>Z<sub>1</sub>, X<sub>1</sub></td><td>-0.0014</td><td>0.253</td><td>0.249</td><td>0.249</td><td>94.3</td></tr><tr><td>GMM</td><td>Z<sub>1</sub>, X<sub>1</sub>, X<sub>2</sub>, X<sub>3</sub></td><td>-0.0008</td><td>0.261</td><td>0.252</td><td>0.252</td><td>95.3</td></tr><tr><td rowspan="12">5</td><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>KM estimator</td><td>-</td><td>-0.0062</td><td>0.251</td><td>0.262</td><td>0.263</td><td>93.5</td></tr><tr><td>Zucker (1998)</td><td>-</td><td>-0.0098</td><td>0.251</td><td>0.261</td><td>0.261</td><td>93.5</td></tr><tr><td>Zucker (1998)</td><td>Z<sub>1</sub></td><td>-0.0024</td><td>0.215</td><td>0.223</td><td>0.223</td><td>94.6</td></tr><tr><td>Zucker (1998)</td><td>Z<sub>1</sub>, Z<sub>2</sub></td><td>-0.0031</td><td>0.211</td><td>0.219</td><td>0.219</td><td>94.7</td></tr><tr><td>Zucker (1998)</td><td>Z<sub>1</sub>, Z<sub>2</sub>, X<sub>1</sub>, X<sub>2</sub></td><td>-0.0031</td><td>0.211</td><td>0.220</td><td>0.220</td><td>94.0</td></tr><tr><td>Andersen et al. (2004)</td><td>-</td><td>-0.0061</td><td>0.252</td><td>0.263</td><td>0.263</td><td>93.5</td></tr><tr><td>Andersen et al. (2004)</td><td>Z<sub>1</sub></td><td>-0.0025</td><td>0.227</td><td>0.235</td><td>0.235</td><td>93.8</td></tr><tr><td>Andersen et al. (2004)</td><td>Z<sub>1</sub>, Z<sub>2</sub></td><td>-0.0045</td><td>0.225</td><td>0.233</td><td>0.233</td><td>93.9</td></tr><tr><td>Andersen et al. (2004)</td><td>Z<sub>1</sub>, Z<sub>2</sub>, X<sub>1</sub>, X<sub>2</sub></td><td>-0.0039</td><td>0.225</td><td>0.234</td><td>0.234</td><td>93.8</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Zhang and Yin (2023)</td><td>-</td><td>-0.0061</td><td>0.250</td><td>0.262</td><td>0.263</td><td>93.6</td></tr><tr><td>GMM</td><td>-</td><td>-0.0029</td><td>0.253</td><td>0.261</td><td>0.261</td><td>94.2</td></tr><tr><td>GMM</td><td>Z<sub>1</sub></td><td>0.0046</td><td>0.233</td><td>0.234</td><td>0.234</td><td>94.8</td></tr><tr><td>GMM</td><td>Z<sub>1</sub>, Z<sub>2</sub></td><td>0.0021</td><td>0.232</td><td>0.232</td><td>0.232</td><td>94.7</td></tr><tr><td>GMM</td><td>Z<sub>1</sub>, Z<sub>2</sub>, X<sub>1</sub>, X<sub>2</sub></td><td>0.0032</td><td>0.239</td><td>0.234</td><td>0.234</td><td>95.2</td></tr></tbody></table>

In practical applications, distinguishing prognostic variables Z_j from non-prognostic ones X_j can be challenging. There is a risk of model misspecification due to the accidental inclusion of variables that do not influence the survival outcome. Therefore, we also assessed the performance of the methods under model misspecification by including unrelated variables in the model in addition to the prognostic ones. We extended the analysis of RMST for scenarios 4 and 5 by adding one or multiple covariates, X_i, unrelated to the survival as adjustment covariates. In both scenarios, the nuisances variables were simulated as X_1 ~ N(0, 1), X_2 ~ Bin(0.5), and X_3 ~ U(0, 2). The associated results for n = 200 correspond to the 3rd (X_1) and 4th rows (X_1, X_2, X_3) of each method for scenario 4 and the 4th row (X_1, X_2) of each method for scenario 5 in Table 2. In scenario 4, no bias nor variance inflation was observed for all the methods. Similar results were observed for n = 500 (supporting information table 8). For smaller sample sizes (n = 50, 100), the impact of misspecification is more important: for the Bayesian GMM, the more unrelated variables were included, the larger the ASE became, while the frequentist approaches had only a slight increase of ASE with lower coverage rate than the nominal 95% level (supporting information table 8). Similar results were observed with a delayed treatment effect (scenario 5): for n = 200 (and n = 500 supporting information Table 9), when adding unrelated variables, the

<footnote n="1">¹ASE: Average Standard Error, ²ESE: Empirical Standard Error, ³RMSE: Root Mean Square Error
Prognostic variables Z₁ ~ N(0, 1), Z₂ ~ Bin(0.5); Other variables X₁ ~ N(0, 1), X₂ ~ Bin(0.5), X₃ ~ U(0, 2)</footnote>
estimations of the dRMST remained unbiased with ASE and ESE comparable to the other methods and satisfactory coverage rate. When the sample size decreases (n = 50 and 100), the inflation of the ASE became more pronounced for the Bayesian GMM and resulted in large credible intervals that were conservative: 100% coverage rate for n = 50 when adjusting for X₁ and X₂ in addition to Z₁ and Z₂ (supporting information Table 9). However, the ESE of the Bayesian GMM did not inflate significantly, suggesting that the results may be strongly influenced by the wide prior we selected. Therefore, we extended the analysis in this particular setting (i.e. scenario 5, n = 50 and adjusting the analysis for Z₁, Z₂, X₁ and X₂) by considering a normal prior with smaller variance, N(0, 1), for all parameters, resulting in an ASE of 0.609 compared to 2.091 with the default prior of N(0, 10), supporting this assumption. On the contrary, we did not observe an increase in ASE when adjusting for unrelated variables for the frequentist models. Consequently, for n = 50, we observed confidence intervals narrower and resulting in a lower coverage rate than the nominal 95% level: 92.2% coverage rate for Zucker (1998), 91.9% coverage rate for Andersen et al. (2004) (supporting information Table 9).

## 3.5 Model specification with an interaction term

In scenario 6, the survival time was simulated according to a Weibull distribution with $\sigma = 0.8$ and $\lambda = \exp(-1.2 + \log(HR_1)A_i + \log(HR_2)E_i + \log(HR_3)A_i \times E_i)$, where $HR_1 = 1.7$, $HR_2 = 0.5$, and $HR_3 = 0.3$. Under these settings, the proportional hazards assumption is satisfied within each subgroup $E = E^-$ and $E = E^+$ but is clearly violated when considering the whole population with crossing survival curves (supporting information Figure 2). Firstly, we analyzed the simulated data of scenario 6 without taking into account biomarker E and its interaction with the treatment. When $n \ge 100$, Bayesian GMM gave dRMST estimates between the two treatment arms comparable to the four benchmark methods with no bias and an appropriate coverage rate (supporting information Table 10). As observed in the previous scenarios, for n = 50, the Bayesian GMM gave an ASE (0.523) slightly higher than the other methods (range from 0.494 to 0.505), leading to a coverage rate closer to the nominal 95% level: 94.1% for the Bayesian GMM and ranging from 92.0% (Zhang and Yin, 2023) to 92.7% (Zucker, 1998). Secondly, we analyzed simulated data of scenario 6 by specifying a model including biomarker E and interaction term with treatment $A \times E$. We limited this second analysis to pseudo-observation-based methods, as they are the only methods among those considered that allow for the inclusion of an interaction term in the model definition: the naive KM method and Zhang and Yin (2023) approach do not allow for covariate adjustments, and Zucker (1998) fit a Cox regression model stratified on the treatment variable to estimate the survival function conditional on the other covariates, before estimating an average survival function, averaging out the covariate parameters. Here, we consider the following model on RMST

$$ \mu_i = E(\min(\tilde{T}_i, \tau)|A_i, E_i) = g^{-1}(\beta_0 + \delta A_i + \beta_1 E_i + \beta_2 A_i \times E_i). \quad (11) $$

The regression model (11) can be expressed equivalently as

$$ \mu_i = g^{-1}(\beta_0 + \beta_1 E_i + \delta^{-} A_i \mathbb{1}_{E_i=E^{-}} + \delta^{+} A_i \mathbb{1}_{E_i=E^{+}}), \quad (12) $$

where $\mathbb{1}_{E=E^{-}}$ and $\mathbb{1}_{E=E^{+}}$ are indicator functions representing whether the biomarker for the i-th individual was $E^{-}$ or $E^{+}$, respectively. This parametrization is adequate to estimate both parameters $\delta^{-}$, which corresponds to the dRMST between the two treatment arms for patients with $E = E^{-}$, and $\delta^{+}$, which corresponds to the dRMST between the two treatment arms for patients with $E = E^{+}$. The parameter $\beta_1$ corresponds to the dRMST between the control patients with $E^{+}$ and the control patients with $E^{-}$. For n = 200, we observed that both models produced valid estimates, ASE close to the ESE, and adequate coverage rate (Table 3). We note a slightly higher bias with the Bayesian GMM than in Andersen et al. (2004) for some coefficients. Still, the corresponding ESE was always smaller with the Bayesian GMM, resulting in an RMSE of the Bayesian GMM always smaller or equal to Andersen et al. (2004). For higher sample size (n = 500), the frequentist and Bayesian methods on pseudo-observations appeared equivalent in terms of bias, variance, and coverage. For smaller sample sizes (especially n = 50), we observed a coverage slightly lower than the nominal 95% level for Andersen et al. (2004) and slightly higher for the Bayesian GMM (supporting information Table 11).
<table><caption>Table 3: Performance of the frequentist and Bayesian methods, based on pseudo-observations, on the estimation of the difference of restricted mean survival time (dRMST) between two treatment groups, adjusted on a binary covariate E with a predictive value on the treatment effect and an interaction term A × E, with non-proportional hazards treatment effect (scenario 6: crossing hazards), with 30% of censoring and a sample size of 200. The parameter $\delta^{-} = -0.9025$ and $\delta^{+} = 0.6492$ corresponds to the dRMST between the two treatment arms for patients with E = E$^{-}$ and E = E$^{+}$, respectively and $\beta_{1} = 1.0644$ corresponds to the dRMST between E$^{+}$ and E$^{-}$ in control patients.</caption><thead><tr><th>Scenario</th><th>Methods</th><th>Estimated coefficient</th><th>Bias</th><th>ASE<sup>1</sup></th><th>ESE<sup>2</sup></th><th>RMSE<sup>3</sup></th><th>Coverage</th></tr></thead><tbody><tr><td rowspan="6">6</td><td><strong>Frequentist</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Andersen et al. (2004)</td><td>β<sub>1</sub></td><td>0.0074</td><td>0.330</td><td>0.340</td><td>0.340</td><td>93.9</td></tr><tr><td>Andersen et al. (2004)</td><td>δ<sup>-</sup></td><td>0.0060</td><td>0.307</td><td>0.310</td><td>0.310</td><td>94.0</td></tr><tr><td>Andersen et al. (2004)</td><td>δ<sup>+</sup></td><td>0.0021</td><td>0.292</td><td>0.309</td><td>0.309</td><td>92.4</td></tr><tr><td><strong>Bayesian</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GMM</td><td>β<sub>1</sub></td><td>0.0212</td><td>0.339</td><td>0.333</td><td>0.334</td><td>95.2</td></tr><tr><td>GMM</td><td>δ<sup>-</sup></td><td>0.0258</td><td>0.315</td><td>0.305</td><td>0.306</td><td>94.6</td></tr><tr><td>GMM</td><td>δ<sup>+</sup></td><td>0.0022</td><td>0.301</td><td>0.306</td><td>0.306</td><td>93.8</td></tr></tbody></table>

## 4 Illustration on a real-data example

We analyzed the data from the Getug-AFU 15, a randomized phase 3 trial comparing an androgen-deprivation therapy (ADT) alone (n = 193) or with docetaxel (n = 192) in non-castrate metastatic prostate cancer (Gravis et al., 2013). The median follow-up time was 4.2 years. Our endpoint of interest was the Prostate-Specific Antigen (PSA) progression-free survival for which the PH assumption was rejected (p = 0.00022, Grambsch and Therneau (1994) test). We analyzed the complete case (n = 345 excluding 40 patients from RCT due to missing covariates) to compare the RMST estimates. 24% of the patients were censored (supporting information Figure 3 shows the Kaplan-Meier curves for the two arms). The same methods, priors, and tuning parameters used for the simulations study were applied to the Getug-AFU 15 analysis. Figure 2 reports the 5-year dRMST for the Bayesian GMM and the four benchmark methods with similar results for unadjusted dRMST estimates between the two arms. We observed an increase in precision of the dRMST estimates with narrower 95% confidence or credibility intervals) adjusted for four binary variables: the Gleason score (≥ 8 vs. < 8), European Cooperative Oncology Group performance status (1 – 2 vs. 0), concentration of alkaline phosphatases (abnormal vs. normal), and presence of bone metastases (yes vs. no) with all methods allowing for covariates adjustment. These results are consistent with the findings of the simulation study. For Bayesian GMM, the trace plots of the NUTS sampling of the 3 chains mixed well and appeared stationary, suggesting no divergence issue, supporting information Figure 4.

With the Bayesian GMM, the adjusted 5-dRMST was estimated to be 0.57 year (95% CI: 0.24 – 0.91), meaning that receiving docetaxel in addition to ADT increases the lifetime without PSA progression over the next 5 years by 0.57 year, compared to receiving ADT alone. See supporting information Table 12 for the details of the parameter estimates with their standard deviation as well as the quantiles of the posterior distributions of the four covariates and the treatment indicator. One advantage of the pseudo-observations over the other approaches is the ability to estimate the baseline RMST, which corresponds to the intercept of the model. In the control group, with all explanatory binary covariates set to their reference levels, the estimated average RMST for an individual was $\beta_{0} = 5.60$ years (SE = 0.54). All the other parameters of the model can be interpreted as 5-dRMST for the treatment indicator but also any other explanatory variables. For example, $\beta_{3}$ can be interpreted as 5-dRMST between patients with abnormal and normal concentrations of alkaline phosphatases. Moreover, the advantage of the Bayesian model we developed is to provide a tail posterior probability for any clinically meaningful threshold. For example, the estimated posterior probability of the dRMST between the two treatment arms over 5 years to be higher than 3 months: $P(\delta \ge 3 \text{ months}) = 0.973$ with the Bayesian GMM adjusted for covariates with a Monte Carlo standard error (MC-SE) of 0.003. For comparison, this posterior probability was estimated to be 0.963 (MC-SE = 0.004) for the nonparametric model of Zhang and Yin (2023). As Bayesian GMM allows to take into account covariate adjustment, tail posterior probabilities may also be estimated for each covariate included in the model (Figure 3). For example, $P(\beta_{3} \le -1 \text{ year}) = 0.907$ (MC-SE = 0.005), meaning that, on average, over the next five years, there is a high probability that patients with an abnormal concentration of alkaline phosphatases will experience PSA progression at least one year before the patients with a normal concentration of alkaline phosphatases.

<footnote n="1">¹ASE: Average Standard Error, ²ESE: Empirical Standard Error, ³RMSE: Root Mean Square Error</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Estimation of the difference of 5-RMST between the ADT plus docetaxel and ADT alone arms for the Prostate-Specific Antigen (PSA) progression-free survival endpoint from the Getug-AFU 15 trial. The horizontal lines represent the 95% confidence or credibility intervals (CI). The dashed vertical line represents no treatment effect.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Posterior distributions of parameters estimated from the Bayesian generalized method of moments applied to the Getug-AFU 15 trial for the PSA progression-free survival analysis. The intercept coefficient represents the 5-year RMST for the control group, with all binary covariates set to their reference levels. The remaining coefficients represent the 5-year dRMST between the two groups defined by each corresponding binary explanatory variable in the model.</figcaption></figure>

While the previous results reported the estimates of the dRMST for $τ = 5$ years, we conducted a sensitivity analysis to assess how varying the restricted time $τ$ affects these results. We limited this sensitivity analysis to the unadjusted estimation of dRMST. Specifically, we examined 4 alternatives for the choice of $τ$ based on data-driven definitions: (a) the 90-th percentile of the observed times, (b) the largest time for which the standard error of the survival estimate is
lower a reasonable limit, either 5% (b1) or 7.5% (b2), as outlined by Karrison (1997), (c) the minimum of the maximum observed times in each treatment group, based on the approach suggested in Uno et al. (2022). When applied to the Getug-AFU 15 trial, the alternative $\tau$ values ranged from 4.60 years (a) to 6.66 years (c), close to the original threshold of 5 years, leading to results generally consistent with those presented in Figure 2 (supporting information Table 13). The larger $\tau$ was, the larger the dRMST estimations were, with one exception for (c). This tendency was expected because we were integrating over a larger period of time with an experimental curve above the control curve until 6 years, see supporting information Figure 5. With large values of $\tau$, the proportion of subjects still at risk is small, resulting in larger confidence and credibility intervals.

# 5 Discussion

This paper presents a novel Bayesian approach for analyzing RMST with covariate adjustment using pseudo-observations in combination with the Bayesian GMM. The approach only requires the specification of a linear model but a model for the survival function is not needed, making it more straightforward than existing Bayesian methods. In addition, next to the posterior distribution of the treatment effect on RMST, posterior distributions of other covariate effects on the RMST can also be estimated, which is unique to this method.

While situations of non-PH are becoming increasingly common to analyze (Lin et al., 2020) and alternative measures of the treatment effect are needed (Uno et al., 2014), this approach provides the option to adjust the dRMST estimation for covariates and interaction terms between the treatment and covariates in a Bayesian framework. Recently, oncologists are increasingly interested in Bayesian methods for designing RCTs (Goligher, Heath, and Harhay, 2024), one attractive feature being the interesting interpretation of the results through the estimation of the whole posterior distribution of the parameters and not only point estimate with confidence bounds. As shown in the example of the Getug-AFU study, the Bayesian GMM enables the estimation of the probability of the dRMST between treatment groups (or any other covariate of interest) to exceed a desired threshold from the posterior distribution, facilitating meaningful decision-making.

Simulation results indicate that the Bayesian GMM with pseudo-observations yielded results comparable to benchmark estimators for moderate ($n = 100, 200$) and large ($n = 500$) sample sizes in scenarios where survival is not related to covariates (scenario 1, 2, 3). When survival is related to one or two covariates, similar findings were observed between Bayesian GMM and other benchmark methods (scenarios 4 and 5). Omitting important variables produced less precision, and including covariates unrelated to the survival did not affect the treatment effect estimation for large samples ($n \ge 200$). With a small sample size, the prior significantly impacted the posterior distribution. In such cases, using a wide prior may lead to excessive variability in parameter estimates, as the Bayesian algorithm is more likely to sample across a broad range of values. This issue became more pronounced as the number of parameters increased. In contrast, frequentist approaches (Zucker, 1998; Andersen et al., 2004) have confidence intervals too narrow, potentially due to unverified asymptotic properties in small samples. Unlike the frequentist approaches, the inference of Bayesian methods only relies on the MCMC algorithm or its derivative, making them not prone to under-coverage for small sample sizes. The asymptotic proprieties of pseudo-observations, proved by Overgaard et al. (2017) and mentioned in section 2.2, as well as equation (8) are provided in this article to justify the Bayesian GMM in this context but do not play a role in the inference process. Nevertheless, it is rare to adjust for several variables with such a small sample size ($n = 50$ and 30% censoring), making this situation less common in real-world applications. One limit of our simulation study is to use a fixed 30% censoring to focus more on simulating various scenarios with non-PH. However, such high censoring is rare in practical applications with a small sample size. We can extrapolate that the results will be consistent when considering lower censoring. The second limit is to have simulated only linear relationships of the covariates in the Weibull model. Caution must be taken regarding the potential misspecification of the model with continuous covariates since the functional form of the relationship between the covariate and the outcome must be specified.

When estimating dRMST, the choice of $\tau$ is important. In our simulation study, we considered the estimation of the difference of RMST at $\tau = 5$ years. The sensitivity analysis on $\tau$ from the Getug-AFU 15 trial did not show an influence on the results, however in other settings, it may impact the RCT conclusion (de Boissieu and Chevret, 2024). The ideal choice is to specify $\tau$ based on clinical considerations at the study design stage. Intuitively, one wants to choose the restricted time as the largest time point available to capture the maximum information from the Kaplan-Meier curve. However, in this application, it drastically increased the variance of the estimates due to the small number of patients still at risk. Further research is needed to generalize these findings to other settings. An approach to mitigate this arbitrary choice would be to estimate the dRMST at multiple time points. Ambrogi, Iacobelli, and Andersen (2022) developed RMST curve estimations based on pseudo-observations. Based on this work, further research will extend the Bayesian GMM based on pseudo-observations to jointly analyze RMST at multiple time points. This can be achieved by considering not only one pseudo-observation per individual but a vector of pseudo-observations computed at different times and using a non-independence working correlation matrix. More generally, the Bayesian GMM
using the pseudo-observations approach, which currently uses an independence working correlation matrix in this paper because patients are independent, can be extended to the RMST Bayesian analysis of correlated data using an appropriate working correlation structure as proposed by Le Vilain-Abraham et al. (2023) for the frequentist analysis of RMST in a cluster randomized trial using an exchangeable working correlation structure.

In conclusion, the proposed Bayesian approach for RMST analysis with covariate adjustment provides straightforward, accurate results and meaningful interpretations, but careful consideration of sample size, prior choice, and choice of the restricted time is essential.

**Additional information**

The R code to analyze RMST with the Bayesian GMM on pseudo-observations is available on the Oncostat team's GitHub https://github.com/Oncostat/pseudo_gmm

**Acknowledgments**

The authors thank C. Zhang for providing the code to implement his Bayesian nonparametric model and for his time and dedication in answering their questions.

This work was awarded the PhD Student Award of the QuanTIM Webinars organized by the SESSTIM research unit (https://sesstim.univ-amu.fr/phd-student-award).

**Conflict of interest**

The authors declare that they have no conflict of interest.

**Funding**

This study was funded by PhD grant MESRI from the doctoral School of Public Health, Paris-Saclay University.

**References**

Ambrogi, F., Iacobelli, S., and Andersen, P. K. (2022) Analyzing differences between restricted mean survival time curves using pseudo-values. *BMC Medical Research Methodology*, 22(1), 71.

Anagnostou, V., Yarchoan, M., Hansen, A. R., Wang, H., Verde, F., Sharon, E., et al. (2017). Immuno-oncology Trial Endpoints: Capturing Clinically Meaningful Activity. *Clinical Cancer Research: An Official Journal of the American Association for Cancer Research*, 23(17), 4959–4969.

Andersen, P. K., Hansen, M. G., and Klein, J. P. (2004). Regression Analysis of Restricted Mean Survival Time Based on Pseudo-Observations. *Lifetime Data Analysis*, 10(4), 335–350.

Andersen P. K., Klein J. P., and Rosthøj S. (2003) Generalised linear models for correlated pseudo-observations, with applications to multi-state models. *Biometrika*, 90(1), 15–27.

de Boissieu, P., and Chevret, S. (2024). Difference in Restricted Mean Survival Times as a Measure of Effect Size: No Assumption Does Not Mean No Rule. *Journal of Clinical Oncology*, JCO.24.00517.

Breslow, N. E. (1972). Discussion of the paper by DR Cox cited below. *Journal of the Royal Statistical Society, Series B*, 34, 187-220.

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., et al. (2017) Stan: A Probabilistic Programming Language. *Journal of Statistical Software*, 76(1), 1-32.

Chen, R., Basu, S., and Shi, Q. (2023). Restricted Mean Survival Time Estimation Using Bayesian Nonparametric Dependent Mixture Models. *arXiv (arXiv:2305.14639)*.

Chernozhukov V., and Hong, H. (2003) An MCMC approach to classical estimation. *Journal of Econometrics*, 115(2), 293-346.

Corro Ramos, I., Qendri, V., and Al, M. (2024). Beyond hazard ratios: Appropriate statistical methods for quantifying the clinical effectiveness of immune-oncology therapies – the example of the Netherlands. *BMC Medical Research Methodology*, 24(1), 260.
Cox, D. R. (1972). Regression Models and Life-Tables. *Journal of the Royal Statistical Society: Series B (Methodological)*, 34(2), 187–202.

Gelman, A., Simpson, D., and Betancourt, M. (2017) The Prior Can Often Only Be Understood in the Context of the Likelihood. *Entropy*, 19(10), 555.

Goligher, E. C., Heath, A., and Harhay, M. O. (2024). Bayesian statistics for clinical research. *The Lancet*, 404(10457), 1067–1076.

Gravis, G., Fizazi, K., Joly, F., Oudard, S., Priou, F., Esterni, B., et al. (2013). Androgen-deprivation therapy alone or with docetaxel in non-castrate metastatic prostate cancer (GETUG-AFU 15): A randomised, open-label, phase 3 trial. *The Lancet. Oncology*, 14(2), 149–158.

Grambsch, P. M., and Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals. *Biometrika*, 81(3), 515–526.

Hanada, K., and Kojima, M. (2024). Bayesian Parametric Methods for Deriving Distribution of Restricted Mean Survival Time. *ArXiv, (arXiv:2406.06071)*.

Hansen, L. P. (1982) Large Sample Properties of Generalized Method of Moments Estimators. *Econometrica*, 50(4), 1029.

Hasegawa, T., Misawa, S., Nakagawa, S., Tanaka, S., Tanase, T., Ugai, H., et al. (2020). Restricted mean survival time as a summary measure of time-to-event outcome. *Pharmaceutical Statistics*, 19(4), 436–453.

Hernán, M. A. (2010). The Hazards of Hazard Ratios. *Epidemiology (Cambridge, Mass.)*, 21(1), 13–15.

Hoffman, M. D., and Gelman, A. (2014) The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. *Journal of Machine Learning Research*, 15(1), 1351–1381.

Højsgaard, S., Halekoh, U., Yan, J., and Ekstrøm, C. T. (2022) *geepack* R package version 1.3.9, https://cran.r-project.org/web/packages/geepack/.

Irwin, J. O. (1949). The standard error of an estimate of expectation of life, with special reference to expectation of tumourless life in experiments with mice. *The Journal of Hygiene*, 47(2), 188–189.

Karrison, T. (1987). Restricted Mean Life With Adjustment for Covariates. *Journal of the American Statistical Association*, 82(400), 1169–1176.

Karrison, T. (1997). Use of Irwin's restricted mean as an index for comparing survival in different treatment groups—Interpretation and power considerations. *Controlled Clinical Trials*, 18(2), 151–167.

Karrison, T. (1997). 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.

Klein, J. P., Gerster, M., Andersen, P. K., Tarima, S., and Perme, M. P. (2008) SAS and R functions to compute pseudo-values for censored data regression. *Computer Methods and Programs in Biomedicine*, 89(3), 289–300.

Le Vilain-Abraham, F., Tavernier, E., Dantan, E., Desmée, S., and Caille, A. (2023). Restricted Mean Survival Time to Estimate an Intervention Effect in a Cluster Randomized Trial. *Statistical Methods in Medical Research*, 32(10), 2016–2032.

Lesaffre, E., Baio, G., and Boulanger, B. (2020). *Bayesian Methods in Pharmaceutical Research*. Chapman and Hall.

Liang, K.-Y. and Zeger, S. (1986) Longitudinal data analysis using generalized linear models. *Biometrika*, 73(1), 13–22.

Lin, R. S., Lin, J., Roychoudhury, S., Anderson, K. M., Hu, T., Huang, B., et al. (2020). Alternative Analysis Methods for Time to Event Endpoints Under Nonproportional Hazards: A Comparative Analysis. *Statistics in Biopharmaceutical Research*, 89(3), 187–198.

McCullagh, P., and Nelder, J. A. (1991) *Generalized Linear Models (2nd ed)*. London: Chapman & Hall.

Mok, T. S., Wu, Y.-L., Thongprasert, S., Yang, C.-H., Chu, D.-T., Saijo, N., et al. (2009). Gefitinib or Carbo-platin-Paclitaxel in Pulmonary Adenocarcinoma. *New England Journal of Medicine*, 361(10), 947–957.

Overgaard, M., Parner, E. T., and Pedersen, J. (2017). Asymptotic theory of generalized estimating equations based on jack-knife pseudo-observations. *The Annals of Statistics*, 45(5), 1988–2015.

Pak, K., Uno, H., Kim, D. H., Tian, L., Kane, R. C., Takeuchi, M., et al. (2017). Interpretability of Cancer Clinical Trial Results Using Restricted Mean Survival Time as an Alternative to the Hazard Ratio. *JAMA Oncology*, 3(12), 1692–1696.

Pohar-Perme, M., Gerster, M., and Rodrigues, K. (2017) *pseudo: computes Pseudo-Observations for Modeling R package version 1.4.3 URL: https://cran.r-project.org/web/packages/pseudo/.*
R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Royston, P., and Parmar, M. K. B. (2011). The use of restricted mean survival time to estimate the treatment effect in randomized clinical trials when the proportional hazards assumption is in doubt. *Statistics in Medicine*, 30(19), 2409–2421.

Royston, P., and Parmar, M. K. B. (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, 152.

Stan Development Team. (2020) Prior Choice Recommendations. Available online https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations (accessed on April 2023).

Stan Development Team (2023). RStan: the R interface to Stan, R package version 2.21.7, https://mc-stan.org/ (accessed on October 2022).

Therneau T. (2024). *A Package for Survival Analysis in R*. R package version 3.7-0, https://CRAN.R-project.org/package=survival. Therneau T (2024). R package version 3.7-0

Tian, L., Zhao, L., and Wei, L. J. (2014). Predicting the restricted mean event time with the subject's baseline covariates in survival analysis. *Biostatistics*, 15(2), 222–233.

Uno, H., Claggett, B., Tian, L., Inoue, E., Gallo, P., Miyata, T., et al. (2014). Moving Beyond the Hazard Ratio in Quantifying the Between-Group Difference in Survival Analysis. *Journal of Clinical Oncology*, 32(22), 2380–2385.

Uno H., Tian L., Horiguchi M., Cronin A., Battioui C., Bell J. (2022). *survRM2: Comparing Restricted Mean Survival Time*. R package version 3.7-0, https://CRAN.R-project.org/package=survRM2

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021). Rank-Normalization, Folding, and Localization: An Improved R^ for Assessing Convergence of MCMC (with Discussion). *Bayesian Analysis*, 16(2), 667–718.

Yin, G. (2009) Bayesian generalized method of moments. *Bayesian Analysis*, 4(2), 191–208.

Zhang, C., and Yin, G. (2023). Bayesian nonparametric analysis of restricted mean survival time. *Biometrics*, 79(2), 1383–1396.

Zucker, D. M. (1998). Restricted Mean Life with Covariates: Modification and Extension of a Useful Survival Analysis Method. *Journal of the American Statistical Association*, 93(442), 702–709.