Olivier Bouaziz

Fast approximations of pseudo-observations in the context of right-censoring and interval-censoring

Keywords: Pseudo-observations, Restricted Mean Survival Time, Von Mises expansions, Jackknife, Interval-censoring

In the context of right-censored and interval-censored data we develop asymptotic formulas to compute pseudo-observations for the survival function and the Restricted Mean Survival Time (RMST). Those formulas are based on the original estimators and do not involve computation of the jackknife estimators. For right-censored data, Von Mises expansions of the Kaplan-Meier estimator are used to derive the pseudo-observations. For intervalcensored data, a general class of parametric models for the survival function is studied. An asymptotic representation of the pseudo-observations is derived involving the Hessian matrix and the score vector. Theoretical results that justify the use of pseudo-observations in regression are also derived. The formula is illustrated on the piecewise-constant-hazard model for the RMST. The proposed approximations are extremely accurate, even for small sample sizes, as illustrated on Monte-Carlo simulations and real data. We also study the gain in terms of computation time, as compared to the original jackknife method, which can be substantial for large dataset.

Introduction

In order to study censored data in time to event analysis it is common to model the hazard rate. This allows to correctly take into account censoring in the estimation procedure and provides hazard ratio estimates in the framework of proportional hazard models. However, in some contexts other quantities, that have a more direct interpretation related to the studied problem, might be of interest. One example is the Restricted Mean Survival Time (RMST) which is defined as the average survival time up to a fixed point. In that case, it is common (see [START_REF] David | Restricted mean life with covariates: modification and extension of a useful survival analysis method[END_REF], [START_REF] Chen | Causal inference on the difference of the restricted mean lifetime between two groups[END_REF], [START_REF] Zhang | Estimating differences in restricted mean lifetime using observational data subject to dependent censoring[END_REF]) to first model the hazard rate using for instance a Cox model, to derive a survival estimator from the estimated hazard rate and to obtain an estimator of the RMST by integrating out this function. This procedure results in a cumbersome computation where it might be difficult to disentangle the effect of each covariate on the RMST. This is a serious drawback for medical applications and there is a need for more direct approaches. There are several other contexts that are concerned by the difficulty of direct modelling of the quantity of interest. This is typically the case for cumulative incidence functions in a competing risk setting or transition probabilities in a multi-state framework.

Pseudo-observations have been developed in the seminal work of [START_REF] Kragh Andersen | Generalised linear models for correlated pseudo-observations, with applications to multi-state models[END_REF] to answer this problem. Those pseudo-observations are constructed using the jackknife method from an estimator of the survival function. Theoretical results in [START_REF] Jacobsen | A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations[END_REF] show that the pseudo-observations can then be used as response variables in a regression model for the quantity of interest, such as the conditional RMST, the cumulative incidence functions in a competing risk setting or the transition probabilities in a multi-state framework. This offers the possibility to directly model the quantity of interest and it is often performed by use of a generalised linear model.

Another more recent areas of development involving pseudo-observations concerns the study of machine learning methods for time to event analysis. In this context, the problematic is similar: one aims at deriving a complex model, based for instance on neural networks, for quantities of interest such as the survival function (see [START_REF] Zhao | Deep neural networks for survival analysis using pseudo values[END_REF]), the cumulative incidence function (see [START_REF] Michael C Sachs | Ensemble prediction of time to event outcomes with competing risks: A case study of surgical complications in crohn's disease[END_REF], [START_REF] Gonzalez | A deep cnn approach for predicting cumulative incidence based on pseudo-observations[END_REF]) or the RMST (see for instance [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF]). The use of pseudo-observations is then appealing since, once the pseudo-observations are obtained, it is possible to directly use any standard machine learning algorithm by considering those pseudo-observations as (non-censored) response variables.

Methods based on pseudo-observations are also attractive for interval-censored data. With those data, it is challenging to build a regression model based on semi-parametric methods, for quantities of interest such as the RMST. This is due to the lack of informations induced by interval-censoring. As a matter of fact, even in a nonparametric setting it may be problematic to perform estimation of the survival function. In this context, one usually relies on the Turnbull estimator or the convex minorant method which were introduced in [START_REF] Turnbull | The empirical distribution function with arbitrarily grouped, censored and truncated data[END_REF] and [START_REF] Groeneboom | Information bounds and nonparametric maximum likelihood estimation[END_REF], respectively. In [START_REF] Groeneboom | Information bounds and nonparametric maximum likelihood estimation[END_REF] it has been proved that these estimators achieve the slow rate of convergence of order n 1/3 and their distribution is not Gaussian and cannot be explicitly computed. In a regression context, for the estimation of the hazard rate, the Cox model with nonparametric baseline was studied in [START_REF] Huang | Efficient estimation for the proportional hazards model with" case 2" interval censoring[END_REF] but again, the baseline survival function has the n 1/3 slow rate of convergence and the asymptotic distribution of this estimator could not be derived. As a result, it is common to rely on fully parametric models for modelling quantities such as the survival function or the hazard rate in a regression context. In [START_REF] Lindsey | A study of interval censoring in parametric regression models[END_REF] and [START_REF] Sun | The statistical analysis of interval-censored failure time data[END_REF] a Cox model was studied using parametric baselines such as Weibull or piecewise constant. The methods used to perform estimation are based on maximum likelihood theory where the parametric estimators are derived by maximising the likelihood of the observed data. This allows to recover the classical √ n rate of convergence of the parametric estimators. However, the derivation of the estimators is not explicit, even in the absence of covariates, and rely on a maximisation algorithm such as the Newton-Raphson procedure. In [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF] a different approach was proposed based on the EM algorithm by considering the true event times as unobserved variables. This method has the advantage that direct estimators can be computed in the E-step of the algorithm when no covariates are present, which results in a stable and robust estimation procedure. All the aforementioned methods consider estimation of the survival function or of the hazard rate through proportional hazard assumptions, but they are not suited for direct modelling of the RMST, in a regression context. However, this can be achieved by using the pseudo-observations approach. In [START_REF] Sabathe | Regression analysis in an illness-death model with interval-censored data: A pseudo-value approach[END_REF], an illness-death model was considered, and conditional transition probabilities or RMST were computed based on this approach. In order to compute the pseudo-observations, the cumulative transition intensities were estimated using either a penalised spline approach or assuming a Weibull distribution. Similarly, in [START_REF] Nygård Johansen | Regression models using parametric pseudo-observations[END_REF] pseudo-observations were computed using a spline approach in order to estimate parameters related to the cumulative incidence function in a competing risk setting.

The key concept about pseudo-observations is that they are built based on the unconditional jackknife estimator of the quantity of interest. While applying the jackknife is straightforward in practice, a limitation of this method comes from the computation burden of calculating the initial estimator n times, where n is the sample size. There exists some R functions designed to improve the computation time, such as the jackknife function in the prodlim package, or the pseudo independent function in the eventglm package, which rely on a C++ implementation but the gain is limited as the initial estimator still needs to be implemented n times. The computational burden is particularly important for interval-censored data where there is no direct calculation of the estimators, even in the absence of covariates. In this paper, we develop approximated formulas for pseudo-observations where the jackknife technique does not need to be implemented. In our formulas, the pseudo-observations can be directly computed based on the initial estimator. In the case of right-censored data, we provide formulas based on Von Mises expansion of the Kaplan-Meier estimator. In the case of interval-censored data, we derive general formulas for parametric models that only involve the original estimator, the score function and the Hessian of the density. Those formulas are approximations of the original jackknife procedure in the sense that they are equal to the original pseudo-observations up to a remainder term that tends towards 0 as n tends to infinity.

However, they turn out to have a very high precision even for moderate sample sizes. Since they only involve the original estimator, the score vector and the Hessian matrix in a parametric context, they are extremely fast to compute thus resulting in a drastic reduction of time.

In the next section, we present a brief summary on the pseudo-values approach. In Section 3 we develop asymptotic formulas for computing pseudo-observations of the survival function and the RMST in the context of right-censored data. The case of interval censored data is studied in Section 4. We first discuss the context of nonparametric estimation of the survival function in Section 4.1. Then the asymptotic pseudo-observations formulas are developed for general parametric models in Section 4.2. In Section 4.3, theoretical validation of pseudo-observations for parametric models are provided: those results show that the conditional expectation of pseudo-observations approximate the conditional expectation of the response variable of interest. Simulations studies for modelling the conditional RMST in the context of right-censored or interval-censored data are conducted in Section 5 where precision and computation time of the approximate formulas are evaluated. Finally, two real data are analysed using the proposed methodology in Section 6.

Backgrounds on pseudo-regression estimation methods

Let T * 1 , . . . , T * n be n independent and identically distributed (i.i.d.) time to event variables of interest, let θ be a parameter of the form

θ = E[h(T * i )],
where h is a known function. Then introduce Z 1 , . . . , Z n n i.i.d. covariates and define the conditional expectation θ

(l) = E[h(T * l ) | Z l ].
We further assume there exists an invertible link function g such that g(θ (l) ) = Z l β, where β is a vector of regression parameters of interest. Instead of observing the T * i 's one usually observes a sample X 1 , . . . , X n of i.i.d. variables, from which an estimator θ is constructed. The l th pseudo-observation is then given by:

θ(l) = n θ -(n -1) θ(-l) , (1)
where θ(-l) is the jackknife estimator of θ, that is the estimator θ computed on the sample where the l th observation has been removed. It has been suggested (see [START_REF] Kragh Andersen | Generalised linear models for correlated pseudo-observations, with applications to multi-state models[END_REF]) to estimate β based on the estimating equation

U (β) = n l=1 θ(l) V -1 l ( θ(l) -θ (l) ) = 0,
where θ(l) denotes the derivative with respect to β of θ (l) = g -1 (Z T l β) and V l is a weight matrix. As a result, the estimator β verifies the equation U ( β) = 0 and it has been suggested to use a sandwich estimator to estimate the variance of β (see for instance [START_REF] Kragh Andersen | Generalised linear models for correlated pseudo-observations, with applications to multi-state models[END_REF] for more details).

In the context of right-censored data, where the θ parameter is the survival function evaluated at some time point and θ is its Kaplan-Meier estimator, it has been proved in [START_REF] Graw | On pseudo-values for regression analysis in competing risks models[END_REF] and [START_REF] Jacobsen | A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations[END_REF] that the resulting estimating function has a mean asymptotically equal to zero. More specifically, one observes T = min(T * , C), ∆ = I(T * ≤ C), where C is a right-censoring variable and we define X = (T, ∆). Let θ = S(t) = P(T * > t), for t ∈ [0, τ ]. We further assume:

(i) C ⊥ ⊥ (T * , Z) (ii) P(T ≥ τ ) > 0.
We set X i = (T i , ∆ i ), i = 1, . . . , n be n i.i.d replications of (T, ∆). The authors have proved that:

θ(l) = θ + ψ(X l ) + o P (1), ( 2
)
where ψ is a first order influence function that verifies E( ψ(X l ) | Z l ) = θ (l) -θ. On the other hand, [START_REF] Jacobsen | A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations[END_REF] have shown that the sandwich estimator used to estimate the variance of β is asymptotically biased. However, the authors have concluded that the difference between their corrected variance estimator and the usual sandwich estimator is of minor importance and as a consequence it is customary to use the sandwich estimator for pseudo-regression.

Once the pseudo-observations have been computed, implementation of the estimating equation along with the sandwich variance estimator can be easily performed from the geese function in the geepack R package.

In this article, we will present approximate formulas for computing pseudo-observations in the context of right-censoring in Section 3 and in the context of interval-censoring in Section 4. Instead of directly using Equation ( 1) to obtain the pseudo-observations, we will present approximated formulas that only involve the estimator θ computed on the whole sample. In both Sections 3 and 4, we will first focus our attention on the problem of modelling θ (l) = S(t | Z l ), the conditional survival function evaluated at time t given the covariate Z l . The pseudo-observations can be computed using an estimator of θ = S(t) the unconditional survival function. A standard link function is g(•) = log(-log(•)) which gives rise to the Cox model. More complex functions can be chosen for g, such as neural-networks (see for instance [START_REF] Zhao | Dnnsurv: Deep neural networks for survival analysis using pseudo values[END_REF], [START_REF] Zhao | Deep neural networks for survival analysis using pseudo values[END_REF], [START_REF] Feng | Bdnnsurv: Bayesian deep neural networks for survival analysis using pseudo values[END_REF]), which will provide performant prediction methods for the conditional survival function. Based on those results we will also consider the problem of modelling

θ (l) = E(T * ∧ τ | Z l ) = τ 0 S(t | Z l )dt, (3)
for some τ > 0. This allows to estimate the RMST in a regression context by considering for instance the identity function for g or again more complex link functions such as neural-networks (see for instance [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF]). In the context of right-censored data only, θ will be estimated based on the Kaplan-Meier estimator and in the context of interval-censored data, it will be estimated based on parametric models.

Approximate pseudo-observations for right-censored data

In this section, we use the same notations as in Section 2 for right-censored data. We denote by Ŝ the Kaplan-Meier estimator of S and we define for l = 1, . . . , n, the l th jackknife estimator Ŝ(-l) of Ŝ as the estimator constructed when omitting the l th observation X l = (T l , ∆ l ). Introduce Ĥ(t) = i I(T i ≥ t)/n and the observed counting process

N i (t) = I(T i ≤ t, ∆ i = 1), where I(•) is the indicator function. Let Λ(t) = 1 n n i=1 t 0 dN i (u) Ĥ(u)
be the standard Nelson-Aalen estimator (see [START_REF] Andersen | Statistical models based on counting processes[END_REF]) of the cumulative hazard function and define the martingale residuals

Ml (t) = N l (t) - t 0 I(T l ≥ u)d Λ(u).
Proposition 1. Under Assumptions (i) and (ii) in Section 2 the following results hold:

Ŝ(l) (t) = n Ŝ(t) -(n -1) Ŝ(-l) (t) = Ŝ(t) -Ŝ(t) t 0 d Ml (u) Ĥ(u) + o P (1)
,

and τ 0 Ŝ(l) (t)dt = τ 0 Ŝ(t)dt - τ 0 τ u Ŝ(t)dt d Ml (u) Ĥ(u) + o P (1)•
From those formulas it is clear that the pseudo-observations can be approximated from quantities computed on the original sample. In other words pseudo-observations can be computed without performing the jackknife procedure. This results in a drastic reduction of the computation time of those pseudo-observations as illustrated in the simulation section. We will also see that this approximation is very accurate even for moderate sample sizes. Besides, the main interest for using this formula is for large sample sizes, in particular in a machine learning context where computing pseudo-observations is the first step of the procedure before applying algorithms such as neural networks. In those contexts, the order of the sample size is often in millions. Even though Proposition 1 is a direct consequence of the results from [START_REF] Graw | On pseudo-values for regression analysis in competing risks models[END_REF] and [START_REF] Jacobsen | A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations[END_REF], a separate proof is provided in the Appendix section.

It should be noted that a different approach for obtaining a fast approximation of pseudoobservations has been considered, based on the infinitesimal jackknife method (see [START_REF] Louis | The infinitesimal jackknife[END_REF]). It has been implemented in the survival package through the pseudo function. Two versions of the infinitesimal jackknife have been implemented, depending on the method used to estimate the survival function. The pseudo function takes as input a survival function that has either been estimated using the Kaplan-Meier estimator or using the Breslow estimator (the exponential of minus the Nelson-Aalen estimator). Both give very similar formulas as the one proposed in this paper. In particular, all three formulas are asymptotically equivalent. Details about the infinitesimal jackknife and its implementation in the survival package are given in the Supporting Information.

Approximate pseudo-observations for interval-censored data

In this section, we suppose that instead of directly observing T * we observe a random interval [L, R], L ≥ 0 and L ≤ R, which almost surely contains the event time: P(T * ∈ [L, R]) = 1. The right end interval is allowed to take the infinite value such that:

• if 0 < L < R < ∞ the data are strictly interval-censored, • if 0 = L < R < ∞ the data are left-censored, • if 0 < L < R = ∞ the data are right-censored, • if 0 < L = R < ∞ the data are exactly observed.
Using the notations of Section 2 the data then consist of i.i.d. replications X i = (L i , R i ), i = 1, . . . , n. This situation is often called interval-censoring case 2 (see [START_REF] Groeneboom | Information bounds and nonparametric maximum likelihood estimation[END_REF]) when exact observations are not allowed and mixed interval censoring (see [START_REF] Yu | On consistency of the self-consistent estimator of survival functions with interval-censored data[END_REF]) or partly interval censoring (see [START_REF] Huang | Asymptotic properties of nonparametric estimation based on partly intervalcensored data[END_REF]) otherwise. In order to derive consistent estimators of the survival function under interval censoring one will usually assume independent censoring in the following way (see for instance [START_REF] Zhang | Regression analysis of interval-censored failure time data with linear transformation models[END_REF]): P(T * ≤ t | L = l, R = r) = P(T * ≤ t | l ≤ T * ≤ r). This supposes that the variables (L, R) do not convey additional information on the law of T * apart from assuming T * to be bracketed by L and R.

Comments on the nonparametric case

It seems appealing to use the same methodology for interval-censored data as in Section 3. A natural nonparametric estimator in the context of interval-censored data is the Turnbull estimator which can be seen as an EM estimator and is consequently rather slow to compute. The gain for avoiding computing n times the Turnbull estimator would therefore be highly significant.

However, in [START_REF] Groeneboom | Information bounds and nonparametric maximum likelihood estimation[END_REF] and [START_REF] Wellner | Interval censoring, case 2: alternative hypotheses[END_REF] it has been showed that this nonparametric maximum likelihood estimator converges at the n 1/3 or (n log(n)) 1/3 rates. Therefore it will not be possible to derive a relation of the following type:

θ = θ + 1 n n i=1 ψ(X i ) + o P (n -1/2 ), with ψ verifying E( ψ(X l ) | Z l ) = θ (l)
-θ such as derived in [START_REF] Jacobsen | A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations[END_REF] and [START_REF] Graw | On pseudo-values for regression analysis in competing risks models[END_REF]. Because if that would be the case, the convergence rate of θ would be of the order n 1/2 due to the central limit theorem. However, this result is crucial to derive Equation ( 2) and assess the validity of the procedure. An alternative could be to use results from [START_REF] Huang | Asymptotic properties of nonparametric estimation based on partly intervalcensored data[END_REF] where it is further assumed that n 1 /n tends to a positive constant as n tends to infinity, where n 1 is the number of exact observations. Under this assumption, the authors retrieved a n 1/2 rate of convergence for the nonparametric maximum likelihood estimator which converges toward a centred gaussian process. However, the covariance function of this process is not explicit and can only be determined as the solution of two integrals. This is caused by the construction of the nonparametric estimator that has no closed form but verifies a self-consistency equation. The asymptotic distribution of the nonparametric estimator was derived using results from infinite dimensional M-estimators in [START_REF] Van Der | Efficiency. of infinite dimensional m-estimators[END_REF]. The same properties in M-estimators could be used here to derive approximated formulas for the nonparametric survival estimator. However, a careful examination of the proofs in [START_REF] Huang | Asymptotic properties of nonparametric estimation based on partly intervalcensored data[END_REF] shows that such formulas would lead again to implicit expressions of the pseudo-observations in the same form as the asymptotic limit of the nonparametric survival estimator. Since it does not seem possible to approach those expressions in a straightforward manner we will not pursue this idea. We will focus instead in the next section in modelling the survival function using parametric models.

Parametric modelling of the survival function

We now assume that the common density function of T * 1 , . . . , T * n depends on α 0 ∈ Θ ⊂ R d , the true model parameter of dimension d. We will denote by f * (t; α 0 ), λ(t; α 0 ), Λ(t; α 0 ) and S(t; α 0 ) = exp(-Λ(t; α 0 )) the true density, hazard, cumulative hazard and survival functions of T * , respectively. Instead of directly observing the variables of interest, one usually observes a sample of i.i.d. variables X 1 , . . . , X n which are assumed to have a common density f (t; α 0 ). We will use the notations ∇ log f (t; α 0 ) and ∇ 2 log f (t; α 0 ) to represent the score vector and the Hessian matrix of this log-density where the derivatives are taken with respect to the model parameter α and are evaluated at α = α 0 . The same notations will be used for f * (t; α 0 ). It is important to emphasise the distinction between the notation f * , which represents the density of the true data, and the notation f which represents the density of the observed data. As an illustration, in the general context of mixed interval-censored data, X i = (L i , R i ) with 0 ≤ L i < R i ≤ ∞ and we have (see [START_REF] Sun | The statistical analysis of interval-censored failure time data[END_REF] or [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF])

f (X i ; α) = (S(L i ; α) -S(R i ; α))I(L i = R i ) + (λ(L i ; α)S(L i ; α)) I(L i = R i ),
with the slight abuse of notation S(R i ; α) = 0 if R i = ∞. In the following, we will consider maximum likelihood estimation for the parameter α 0 based on the observed variables X 1 , . . . , X n .

The results derived in this section are not limited to the case of interval-censored data and can be applied to any parametric framework for incomplete data.

The maximum likelihood estimator α of α 0 maximises with respect to α the log-likelihood i log f (X i ; α) and, subject to regularity conditions, verifies the following equality:

√ n(α -α 0 ) = 1 √ n - 1 n n i=1 ∇ 2 log f (X i ; α) -1 n i=1 ∇ log f (X i ; α 0 ),
where α lies between α and α 0 . Let I = -E ∇ 2 log f (X; α 0 ) be the Fisher information and consider the jackknife version α(-l) of the maximum likelihood estimator. It is then straightforward to write:

√ n(α -α 0 ) = 1 √ n I -1 n i=1 ∇ log f (X i ; α 0 ) + ε n , √ n -1(α (-l) -α 0 ) = 1 √ n -1 I -1 n i =l ∇ log f (X i ; α 0 ) + ε (-l) n ,
where

ε n = 1 √ n   - 1 n n i=1 ∇ 2 log f (X i ; α) -1 -I -1   n i=1 ∇ log f (X i ; α 0 ),
and ε

(-l) n
is the jackknife version of ε n . As a result, the l th pseudo-observation of α verifies the relation:

nα -(n -1)α (-l) = α 0 + I -1 ∇ log f (X l ; α 0 ) + √ n ε n - √ n -1 ε (-l) n . (4)
In the Appendix section, it is proved that the term

√ n ε n - √ n -1 ε (-l) n
tends towards 0 in probability as n tends to infinity. This entails that asymptotically, the pseudo-observation of α only depends on the true parameter, the Fisher information and the score vector, this latter quantity being only evaluated at the observation l. Since α is a consistent estimator of α 0 and Î = -n i=1 ∇ 2 log f (X i ; α)/n is a consistent estimator of I, a natural asymptotic approximation for the pseudo-observation of α is simply:

α + Î-1 ∇ log f (X l ; α).
While this result is interesting on its own, more work needs to be done in order to derive the pseudo-observations of S(t; α). The following proposition is derived based on this latter expression of the approximate pseudo-observation for α. The notation • is used to denote the transpose of a vector or a matrix. Proposition 2. Under standard regularity conditions for maximum likelihood theory, the following relations hold:

nS(t; α) -(n -1)S(t; α(-l) ) = S(t; α) -S(t; α)∇Λ(t; α) Î-1 ∇ log f (X l ; α) + o P (1),
and for τ > 0,

n τ 0 S(t; α)dt -(n -1) τ 0 S(t; α(-l) )dt = τ 0 S(t; α)dt - τ 0 S(t; α)∇Λ(t; α) dt Î-1 ∇ log f (X l ; α) + o P (1).
The same results also hold when replacing α, Î by α 0 , I in the right-hand side of the equations.

The proof of this proposition can be found in the Appendix section. As in Section 3, the main interest in this result lies in the fact that the approximated version of the pseudo-observation only depends on the parameter estimator α and not on its jackknife version. This means that pseudo-observations in parametric models can be obtained without actually computing the n jackknife estimators. Only the estimator of α 0 , along the Hessian matrix, the gradient of Λ and of the log-density are needed. This is particularly interesting in the context of intervalcensored data since parametric estimators cannot be derived explicitly and numeric methods must be implemented. Two different strategies exist for those types of data: either a direct maximisation of the likelihood can be performed using the Newton-Raphson algorithm (see [START_REF] Lindsey | A study of interval censoring in parametric regression models[END_REF] and [START_REF] Sun | The statistical analysis of interval-censored failure time data[END_REF] for instance) or the complete likelihood (based on the unobserved true times) can be used through the EM algorithm in order to maximise the likelihood (see [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF]). But in either case the method is iterative. Also, it should be noted that the Newton-Raphson algorithm requires to compute the score vector and Hessian matrix. Therefore the computational cost for implementing the approximated pseudo-observations is similar to the cost of simply computing the pointwise estimate from the Newton-Raphson algorithm.

Those approximated formulas are general and work for any parametric model. As an illustration, the piecewise-constant hazard (pch) model will be used in the simulation section. This model assumes that the hazard function verifies λ(t; α)

= K k=1 α k I k (t) where I k (t) = I(c k-1 < t ≤ c k ), c 0 = 0 < c 1 < • • • < c K = +∞ represent K + 1
cuts and I(•) denotes the indicator function. We do not specify precisely the regularity conditions for maximum likelihood theory to hold. However, two important assumptions are first to assume the model identifiable and second to impose the Fisher information to be positive definite in a neighbourhood of the true parameter. For the pch model in the context of interval-censored data, two necessary conditions for those regularity assumptions to hold are:

P(R < +∞, [L, R] ∩ (c k-1 , c k ] = ∅) > 0, ∀k = 1, . . . , K, P(L > c k-1 ) > 0, ∀k = 1, . . . , K. (5)
The first assumption is quite natural: in order to estimate α k , the probability that an interval intersects [c k-1 , c k ] should be positive. The second assumption is necessary for the existence of a maximum of the likelihood function. It should be noted that those conditions are also valid when exact observations L = R are allowed. Exact expressions of the score vector and Hessian matrix for the pch model along with the derivation of condition (5) are detailed in Section 9.4 of the Appendix. Details on the implementation of Proposition 2 for the pch model are given in Section 9.3 of the Appendix. Precision and computational cost of the approximation for the RMST are evaluated and compared to the actual jackknife version of the pseudo-observations in the simulation section. In particular, it is seen that the approximation is much faster than the jackknife method and is very accurate even for small sample sizes.

Theoretical validation of pseudo-observations for parametric models

In this section we want to investigate if the approximated formula derived in Proposition 2 provides valid observations for performing pseudo-regression, similarly to the Kaplan-Meier estimator in the context of right-censored data (see Equation ( 2)). In other words, if we set

ϕ(X l ; α 0 ) = τ 0 S(t; α 0 )dt - τ 0 S(t; α 0 )∇Λ(t; α 0 ) dt I -1 ∇ log f (X l ; α 0 ), (6)
we want to investigate under which conditions we may have

E(ϕ(X l ; α 0 ) | Z l ) = E(T * l ∧ τ | Z l ). (7)
It is easily seen that this equality will generally not hold by considering the simple scenario of exact observations. In that case, X i = T * i and f = f * is simply the density of the true variable T * . If we further assume for instance that T * follows a Weibull distribution, with shape parameter a > 0 and scale parameter b > 0 (the true parameters are noted a 0 , b 0 ) such that:

f (X; α) = f (T * ; α) = a b T * b a-1 exp - T * b a ,
with α = (a b) , then ∇ log f (X l ; α) will depend on log(T * l ), (T * l ) a-1 and (T * ) a , when a = 1. As a result, E(ϕ(X l ; α 0 ) | Z l ) will be a function of

E(log(T * l ) | Z l ), E((T * l ) a 0 -1 | Z l ) and E((T * l ) a 0 | Z l ) when a 0 = 1. When a 0 = 1 (the exponential model), then E(ϕ(X l ; α 0 ) | Z l ) will be a function of E(T * l | Z l )
, but it will still not verify Equation ( 7) unless τ = ∞. Performing the same calculation for other distributions, we can similarly conclude that Equation [START_REF] Michael C Sachs | Ensemble prediction of time to event outcomes with competing risks: A case study of surgical complications in crohn's disease[END_REF] will not hold in general.

Nevertheless, even though Equality ( 7) is not verified in most cases, it is still possible to prove that the formula given by Proposition 2 provides a good approximation. The key is to assume that there exists a value α z of the parameter α such that the conditional distribution of T * given Z = z follows a distribution with density f * (t; α z ). In that case, E(T * l ∧ τ | Z l = z) = τ 0 S(t; α z )dt and by expanding S(t; α z ) around S(t; α 0 ) from a Taylor development, we can prove that

E(T * l ∧ τ | Z l = z) is equal to E(ϕ(X l ; α 0 ) | Z l = z)
up to two remainder terms. Those remainder terms measure the distance between α 0 and α z , and between the inverse of the Fisher information I -1 and the quantity (-E(∇ 2 log f (X; α) | Z = z)) -1 for any α that is on the real line between α 0 and α z . Define Θ z = {α ∈ Θ : α -α 0 ≤ α 0 -α z } which represents the set of parameters that are on the real line between α 0 and α z . In the next proposition, we denote for k = 1, . . . , K, by α k , α 0,k , α z,k the kth component of α, α 0 and α z , respectively. Proposition 3. Assume there exists α z such that the conditional distribution of T * given Z = z follows a distribution with density f * (t; α z ). Assume also there exists

M z < +∞ such that ∀k, k ∈ {1, . . . , d}, ∀α ∈ Θ z , ∂ 2 ∂α k ∂α k τ 0 S(t; α)dt ≤ M z .
Then,

E(T * l ∧ τ | Z l = z) = E(ϕ(X l ; α 0 ) | Z l = z) + R 1,z + R 2,z ,
where

R 1,z ≤ M z 2 d k=1 (α z,k -α 0,k ) 2 , R 2,z ≤ max α∈Θz τ 0 (∇S(t; α 0 )) dt I -1 α,z -I -1 E(∇ log(f (X; α 0 )) | Z = z)) , with
I α,z = -E(∇ 2 log(f (X; α)) | Z = z).
The proposition makes the strong assumption that the conditional distribution of T * given Z = z follows a distribution with density f * (t; α z ). While this will not be true in general, it seems reasonable to assume that, if the chosen parametric distribution for T * is rich enough, there will exist a value of the parameter α such that the parametric distribution is not too far from the conditional distribution of T * given Z = z. This advocates for flexible parametric models such as the pch model or a spline approach such as proposed in [START_REF] Nygård Johansen | Regression models using parametric pseudo-observations[END_REF]. However, this needs to be imposed for all possible values of Z, which again seems reasonable if discrete covariates are considered and the number of those covariates is not too large.

Besides, it is difficult to evaluate, from the proposition, how large the remainder terms are. Previous experiments on the Weibull distribution in [START_REF] Sabathe | Regression analysis in an illness-death model with interval-censored data: A pseudo-value approach[END_REF] using the jackknife approach suggest that the approximation is quite accurate in practice. When using the pch model, we can establish a different type of theoretical result. In the context or right-censored data, we show in the next proposition that if the number of cuts in the pch model tends to infinity, then we exactly retrieve Equality [START_REF] Michael C Sachs | Ensemble prediction of time to event outcomes with competing risks: A case study of surgical complications in crohn's disease[END_REF]. We were not able to prove this result in the context of interval-censored data, we conjecture however that this result still holds in this case.

c 0 = 0 < c 1 < • • • < c K = +∞
and if we assume standard regularity conditions for maximum likelihood theory then the function ϕ(X l ; α 0 ) defined in Equation (6) converges, as K tends to infinity and

max |c k+1 -c k | tends to 0, towards a function ϕ ∞ (X l ; α 0 ) that verifies E(ϕ ∞ (X l ; α 0 ) | Z l ) = E(T * l ∧ τ | Z l ).
This result is interesting as it shows that it is theoretically possible for Equality [START_REF] Michael C Sachs | Ensemble prediction of time to event outcomes with competing risks: A case study of surgical complications in crohn's disease[END_REF] to hold true when using the pch model. Of course, in practice one has to choose a finite number of cuts. However, there are some strategies to choose the number of cuts from the data. In particular, in [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF] the authors have developed a penalised method based on the adaptive-ridge to choose the number of cuts in an efficient way. In the simulation study (Section 5.2), we show that the approximation formula or the original jackknife method provide very similar and very performant results in pseudo-regression, when modelling the distribution of T * with the pch model. The proofs of Propositions 3 and 4 are deferred to the Appendix section.

Simulation studies for the Restricted Mean Survival Time

We study two different simulation scenarios for the RMST: one with right-censored data and another one with interval-censored data. In the first scenario, the approximate pseudo observations are based on the Kaplan-Meier estimator (using Proposition 1) while in the second scenario they are based on the pch model (using Proposition 2). In both settings, the performance of the estimators derived from the approximated formulas and the ones obtained from the standard jackknife method is compared based on 500 replications. Implementation of the generalised estimation equation is performed through the geese function in the geepack R package. Computation times of the estimators were evaluated on 100 replications, from 10 different samples with 10 replications on each sample using the microbenchmark R package.

Right-censored data

The simulation setting is based on the one in [START_REF] Wang | Modeling restricted mean survival time under general censoring mechanisms[END_REF]. We assume that

T * i = β 0 Z i + ε i , i = 1, . . . , n,
with β0 = (5.5, 0.25, 0.25) ,

Z i = (1, Z i,1 , Z i,2
) , Z i,1 and Z i,2 are Bernoulli variables with parameter 0.5 and ε i ∼ U[-σ, σ], with σ = 3. Under this model it can easily be seen that

E(T * i ∧ τ | Z i ) = β 00 + β 01 Z i,1 (1 -Z i,2 ) + β 10 Z i,2 (1 -Z i,1 ) + β 11 Z i,1 Z i,2 , (8)
where β 0 = (β 00 , β 01 , β 10 , β 11 ) can be determined computationally using Monte-Carlo samples with size 10 million. We further set τ = 6 which corresponds to the 54.2% quantile of T * and to the value β 0 = (4.98, 0.14, 0.14, 0.27) . Right-censored data were simulated from an exponential distribution with parameter 0.07 yielding 33% of censoring on average. The results are presented in Table 1.

It is seen that the approximated formula gives similar results as compared to the standard jackknife method for n = 100. For larger sample sizes, the results are almost identical. We also compared the difference between the two estimators of β 0 by looking at the standard deviation for all four components and taking the maximum: the maximum value over all four components is equal to 7.06 × 10 -3 , 5.90 × 10 -4 , 2.16 × 10 -4 , 6.88 × 10 -6 for n = 100, n = 500, n = 1, 000, n = 10, 000 respectively. This shows that there is very little variations between the estimator computed from the jackknife and the one computed from the approximated formula. In terms of computation times, there is a clear advantage for the approximated formula which goes 14.3, 27.5, 25.1 and 18.7 times faster for n = 100, n = 500, n = 1, 000 and n = 10, 000 respectively. Clearly the computation time for the original jackknife method is not a linear function of the sample size and the gain for using the approximated method is considerable for large sample sizes. It should be noted that the computation time was evaluated for the pseudo-regression procedure, but it does not include the computation of the initial survival estimator, it only takes into account the computation of the pseudo-observations along with the implementation of the generalised estimating equations. Finally, the infinitesimal jackknife implemented in the pseudo function of the survival package was also briefly compared to our approach. It seems that the pseudo function provides very similar results but with a faster computational time. However, since there is no available information on the implementation of the pseudo function, it is not possible to clearly evaluate which of the two formulas (the one obtained with the infinitesimal jackknife and the one obtained from the Von Mises formula) has the smallest computational complexity. 8) based on pseudoregression with the Kaplan-Meier estimator on 33% of right-censored data. In the pseudoregression, the true jackknife is compared to the approximated pseudo-estimates.

Jackknife

Interval-censored data

For interval censored data the survival function is estimated from the pch model, as detailed in Section 4.2. Using this model, estimation of the model parameter α 0 is performed using the EM algorithm, as presented in [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF]. An alternative method could be to directly maximise the observed likelihood but this would result in implementing the Newton-Raphson algorithm for each jackknife sample with inversion of a Hessian matrix of full rank which, in turn, would result in unstable results. In the EM algorithm, the M-step is explicit and as a result the computation of the jackknife methods is always stable. We refer the reader to [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF] for more details on the two methods. The approximated method is implemented from the result in Proposition 2 and details on the computation of the score vector and Hessian matrix are detailed in Section 9.3 of the Appendix. We assume Model [START_REF] Gonzalez | A deep cnn approach for predicting cumulative incidence based on pseudo-observations[END_REF] with the same values of σ and τ . Then, in order to simulate intervalcensored data, a total of K = 5 visits were simulated such that V 1 ∼ U[0, 6] and

V k = V k-1 + U [0, 2]
, for k = 2, . . . , K. The observations for which T * i < V 1 correspond to left-censored observations with L i = 0 and R i = V 1 , the observations for which T * i > V K correspond to right-censored observations with L i = V K and R i = ∞, and the observations for which V k-1 < T * i < V k (k = 2, . . . , K) correspond to strictly interval-censored observations with L i = V k-1 and R i = V k . This resulted in 14.6% of left-censored data, 52.07% of interval-censored data and 33.33% of right-censored data. For interval-censored data, the average length of the intervals was approximately equal to 1.34. The pch model with cuts equal to 4, 5, 6, 7 was used for the computation of the survival estimator. The pseudo-observations were generated based on the standard jackknife and on the approximated formulas and the results for the RMST model are presented in Table 2.

Again the results between the jackknife and the approximate formula are almost identical while there is a huge gain in terms of computational time for the approximated formula. The approximated formula is 107, 198 and 310 times faster than the jackknife method for n = 200, n = 500 and n = 1, 000 respectively. It should be noted that the cuts must be carefully chosen in the pch model. In particular, the regularity conditions of Equation ( 5) must be satisfied. If there are only few values of L i and R i that intersect a cut [c k-1 , c k ] or if the proportion of L i 's such that L i > c k-1 is too low then the pseudo-values can be incorrect (both for the jackknife method or using our approximated formula) which will in turn result in a poor performance of the parameters estimation. On the other hand, if the regularity conditions hold, the choice of the cuts will only have a minor impact on the performance of the estimator of β 0 and will lead to similar results. A supplementary simulation scenario for interval-censored data with τ equal to infinity is also presented in the Supporting Information.

Illustrative real data examples 6.1 The Cardiovascular Health Study (CHS)

In this data example, we mimic the analysis of the Cardiovascular Health Study (CHS) as it was performed in [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF]. This study was initiated in 1987 to determine the risk factors for development and progression of cardiovascular disease (CVD) in older adults. The event of interest was time to CVD. In [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF], the author considers a subsample of 5, 380 individuals of whom 65.2% had CVD during the study period and the others were right-censored. The aim of the study was to estimate the conditional RMST with 29 covariates and τ = 5 years.

The methodology proposed in [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF] uses pseudo-observations and implements a deep neural network directly on the pseudo-observations of the RMST, that is the g link function presented in Section 2 is a neural network. Moreover, a training dataset including 75% of the observations and a test set based on the remaining 25% of the data are built in order to evaluate the prediction performance of the method. This split of the data between training and test sets is repeated 10 times. At each repetition, the pseudo-observations must be entirely computed but only on the training datasets. This results in computing the pseudo-observations for the RMST for 10 samples of size 4, 035. We computed those pseudo observations from the jackknife method and 8) based on pseudoregression with 14.6% of left-censored data, 52.07% of interval-censored data and 33.33% of right-censored data. The piecewise constant hasard model with cuts equal to 4, 5, 6, 7 was used for the estimation of the survival function in the computation of the pseudo-observations. In the pseudo-regression, the true jackknife is compared to the approximated pseudo-estimates.

the approximated formula. The former was computed in approximately 21.3 seconds while the latter took 1.9 seconds. Therefore, our approximated formula is more than 11 times faster than the original jackknife method. Since building a neural network is computationally expensive and needs to be implemented for all the training samples, this reduction in computation time is a major advantage for our approximated formula. Of note, the results of the analysis implemented with the approximated formula are identical to the original analysis (based on the jackknife method) and are therefore omitted.

The Signal Tandmobiel R data

In this section, our aim is to analyse the Signal Tandmobiel R data using three different models: the standard Cox model, a logistic model for the conditional survival function and the conditional RMST model presented in Equation ( 3). This dataset is part of the bayesSurv R package. Those data were collected from a longitudinal dental survey of 4, 468 school children born in 1989, who were annually examined by a dentist. The time scale is age in years. The dataset is composed of 0.68% of left-censored data, 61.69% of strictly interval-censored data and 37.63% of right-censored data. Our aim is to study the emergence of the tooth number 14 which is a permanent first premolar. The covariates used for the analysis are: gender (binary variable equal to 1 for boys, 0 for girls) and the number of decayed or missing deciduous first molars due to caries among teeth 54, 64, 74, 84 of the dataset. This covariate is thus discrete taking values between 0 and 4. These data were previously studied by [START_REF] Lesaffre | An overview of methods for interval-censored data with an emphasis on applications in dentistry[END_REF] using the Accelerated Failure Time model (AFT). In our analysis, the survival function is estimated from the pch model using the whole dataset and the pseudo-observations are then computed from the approximated formula in Proposition 2. There are 126 individuals with missing covariates and the generalised estimating equation used to implement our models is therefore applied to this reduced dataset composed of 4, 342 pupils.

In the pch model, the number of cuts and locations were chosen using the adaptive-ridge algorithm developed in [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF]. This led to the selection of the four cuts 7.6, 8.4, 9 and 10. Since the maximum likelihood estimator has converged this entails that the regularity conditions of Equation ( 5) are satisfied. We can also easily check them empirically: in particular there are 3% of strictly interval censored observations whose left intervals fell before 7.6, 40% of left intervals that fell after 10 and the percentage of strictly interval censored observations that intersect each other is high (values not shown). The corresponding estimated hazard and survival functions are displayed in Figure 1. We observe a low estimated hazard value (equal to 6•10 -4 ) from age 0 until age 7.6 due the low percentage of left intervals that fell before 7.6. This yields a very flat decay of the survival function on this time period, then the decay increases drastically for the four other time periods [7.6, 8.4], [8.4, 9], [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF][START_REF] Turnbull | The empirical distribution function with arbitrarily grouped, censored and truncated data[END_REF] and [10, ∞). For illustration, we estimate from the survival function that approximately 83.39% of the teeth will emerge between age 7.6 and 12.

In order to implement a Cox model from the pseudo-observations, we need to set a grid of time points on which the baseline hazard rate is estimated (see [START_REF] Kragh Andersen | Generalised linear models for correlated pseudo-observations, with applications to multi-state models[END_REF]). We choose the time points t 1 = 8, t 2 = 9, t 3 = 10, t 4 = 11, t 5 = 12 and we use the link function g(x) = log(-log(x)) in the pseudo-regression which leads to the Cox model:

log(-log(S(t m | Z i )) = log Λ 0 (t m ) + Z i β, m = 1, . . . , 5,
where Λ 0 (t m ) represents the cumulative baseline hazard function at time t m . The results are presented in Table 3. The second column provides the cumulative baseline hazard at different time points and the hazard ratios of the two covariates. The last column displays the Wald tests which are all extremely significant. We clearly see that the hazard for the emergence of the tooth is an increasing function of time. Also, under the proportional hazards assumption, boys have an increased risk as compared to girls with an hazard ratio equal to 1.47 and the hazard ratio for one supplementary decayed or missing deciduous first molar due to caries equals 1.13.

For the logistic regression model, we study the probability for the emergence of the tooth before a fixed time point. We conduct two separate analysis, one for the time point 9 and another one for the time point 12 which correspond to the 11% and 84% estimated quantiles of T * , respectively, thus corresponding to early and late emergence of the tooth. There is no guarantees that the pseudo-observations are in the interval [0, 1] and we therefore set the negative values to 0 and the values greater than 1 to 1. The model is the following:

logit(1 -S(t | Z i )) = γ + Z i β,
where t is either equal to 9 or 12, γ is the intercept and β is a two dimensional vector representing the effect of the covariates and logit is the classical logistic function (logit(x) = log(x/(1 -x)). It should be noted that the model can be directly implemented with the glm or geese functions by using one minus the pseudo-observations as input. The results are presented in Table 4. We first observe that the effect of the two covariates are highly significant but they differ depending on the time endpoint: the odds ratio for the number of decayed or missing deciduous first molars decreases from exp(0.2808) ≈ 1.32 before age 9 to exp(0.1080) ≈ 1.11 before age 12 while the odds ratio for the gender effect increases from exp(0.2978) ≈ 1.35 before age 9 to exp(0.5284) ≈ 1.70 before age 12. This strongly suggests that the number of decayed or missing deciduous first molars is mostly responsible for the early teeth emergence while gender (with boys having a higher risk) is mostly responsible for the late tooth emergence. It is also interesting to compare the probabilities of the tooth emergence before age 9 between girls with no decayed or missing deciduous first molars (5.49%), boys with no decayed or missing deciduous first molars first molars (7.26%) and boys with 4 decayed or missing deciduous first molars (19.39%). We can similarly compare the probabilities of the tooth emergence before age 12 between girls with no decayed or missing deciduous first molars (70.63%), boys with no decayed or missing deciduous first molars (80.31%) and boys with 4 decayed or missing deciduous first molars (86.27%).

Finally, two RMST analysis were conducted with τ = 9 and τ = 12. The estimated regression parameters in the RMST model along with their Wald test are presented in Table 5. For τ = 9 we observe a weak effect of the covariates with an intercept that is almost equal to τ , highlighting that most emergences of the tooth will occur after 9 years of age. As a matter of fact, gender is not significant and the number of decayed or missing deciduous first molars is highly significant but with a weak effect. The number of decayed or missing deciduous first molars will accelerate the emergence of the tooth with 1 decayed molar (respectively 4 decayed molars) yielding a reduction of 0.0097 years (respectively 0.0390 years) for the emergence of the tooth. For τ = 12 the effect of gender is now highly significant, meaning that gender only plays a role for late emergence of the tooth (a finding that was also observed with the logistic regression model). The emergence of the tooth for boys arrives on average 0.3336 years earlier than for girls. The number of decayed or missing deciduous first molars is also highly significant with 1 decayed molar (respectively 4 decayed molars) yielding a reduction of 0.1303 years (respectively 0.5211 years) for the emergence of the tooth.

We also tried to repeat the procedure using different cut values in the pch model and as already observed in the simulation study, this led to very similar results, for all three models. The results from [START_REF] Lesaffre | An overview of methods for interval-censored data with an emphasis on applications in dentistry[END_REF] obtained using the AFT models were similar to our findings except that the authors did not provide statistical tests for the effects of the covariates and it was not possible from their method to detect that gender had mainly a role for late emergence of the tooth.

Finally, based on the approximated formulas developed in this paper, the whole procedure (computation of the pseudo-observations and implementation of the generalised estimating equations) took about 1.78 seconds for the RMST analysis (the computation times are similar for the other two models). The method was not implemented using the classical jackknife method but according to the simulation study it would have taken 29.8 minutes to obtain the pseudoobservations, since in the simulation study the time for the jackknife procedure was evaluated at 6.7 minutes for n = 1, 000 (see Sections 5.2). Also, the results would have been identical, thus highlighting the relevance of the proposed approach in practical situations.

Conclusion

In this paper, we presented asymptotic formulas for computing pseudo-observations for time to event data. In the context of right-censored data, those formulas are based on the Kaplan-Meier estimator of the survival function. When dealing with interval-censoring our formulas were developed for a general class of parametric models. Pseudo-regression is an appealing tool when the goal is to directly model a complex quantity of interest, such as the RMST, cumulative incidence functions in a competing risk setting or transition probabilities for multistate models. Our formulas were precisely developed for the RMST but they could be easily extended for those other quantities of interest. While the pseudo-values approach is originally based on the jackknife procedure, our formulas only involve quantities computed on the initial sample. This results in a drastic reduction of the computational time, which is an interesting feature when dealing with large dataset or when the data are interval-censored, since in that case, the estimators are computationally intensive.

There has been an increasing interest of the pseudo-values approach in the machine learning community. After having computed the pseudo-observations, standard machine learning models can be applied to those new observations, by simply ignoring the censoring. In particular, this methodology has been applied for estimating the survival function in [START_REF] Zhao | Dnnsurv: Deep neural networks for survival analysis using pseudo values[END_REF], [START_REF] Zhao | Deep neural networks for survival analysis using pseudo values[END_REF], [START_REF] Feng | Bdnnsurv: Bayesian deep neural networks for survival analysis using pseudo values[END_REF] or the RMST in [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF] based on neural networks that were directly applied on the pseudo-observations. Therefore, our formulas are particularly interesting in those settings where the dataset can be extremely large and the algorithm usually relies on a cross-validation procedure. Using our approximated formulas results in a significant gain in terms of computation time as illustrated on the real data analysis. Also, the approximations made by our formulas are extremely precise, even for moderate sample sizes, as shown in the simulation study. Surprisingly, we also saw that our formulas are more robust than the original jackknife method which sometimes fails due to some rare extreme values. For all these reasons, we advocate the use of our asymptotic formulas in practical situations.

As a reviewer pointed out, there exists an alternative to our fast formula in the case of rightcensored data, in the pseudo function of the survival package. It is based on the infinitesimal jackknife and can be applied to the Kaplan-Meier estimator or to the Breslow estimator (defined as the exponential of minus the Nelson-Aalen estimator). A theoretical comparison of those approximations with the Von-Mises approach developed in this paper has been made in the Supporting Information. It is shown that all three approximations are very similar for small sample sizes and asymptotically equivalent. It would now be interesting to extend the infinitesimal jackknife to the parametric setting, in order to apply it to interval-censored data, as there already exists some approaches based on the infinitesimal jackknife to compute delta-beta residuals, which are very closely related to pseudo-observations.

Finally, new theoretical results for parametric pseudo-regression were also developed in this paper. Those results indicate that parametric pseudo-regression is only valid up to two extra terms but simulation studies suggest that they are reasonably small in practice and do not significantly impact the performance of the final estimator, as long as the number of covariates is not too large. They also suggest to use flexible parametric models such as the pch model. As an extension it would be interesting to study the theoretical validity of splines models for pseudo-regression such as the method developed in [START_REF] Nygård Johansen | Regression models using parametric pseudo-observations[END_REF]. This is let to future research.

Software

The asymptotic formulas developed in this paper for the pseudo-values of the survival function and the RMST can be implemented using the GitHub package FastPseudo available at https://github.com/obouaziz/FastPseudo. The package can deal with both right-censored or interval-censored data. In the latter case, the formulas are implemented for the pch model.

g(H 1 , H) := t 0 dH 1 (u)/H(u) and Λ(t) = g( Ĥ1 , Ĥ). We have the following Von-Mises expansion:

Λ(-l) (t) = Λ(t) + g ( Ĥ1 , Ĥ) ( Ĥ(-l) 1 -Ĥ1 , Ĥ(-l) -Ĥ) + o P (n -1 ),
where g is the Hadamard derivative of g, which is equal to [see [START_REF] Richard | Lectures on survival analysis[END_REF][START_REF] Aad | Weak convergence[END_REF]:

g (H 1 ,H) (h 1 , h) = t 0 dh 1 H - t 0 h 2 dH 1 H 2 •
The o P (n -1 ) term above comes from the expressions: Ĥ(-l)

1 (t) -Ĥ1 (t) = 1 n(n -1) n i=1 I(T i ≤ t, ∆ i = 1) - I(T l ≤ t, ∆ l = 1) n -1 and Ĥ(-l) (t) -Ĥ(t) = 1 n(n -1) n i=1 I(T i ≥ t) - I(T l ≥ t) n -1 ,
which entail as a consequence that Ĥ(-l)

1
(t) -Ĥ1 (t) and Ĥ(-l) (t) -Ĥ(t) are O P (n -1 ). Moreover, using those expressions we have 2 • Gathering all the different parts, we obtain Ŝ(l) (t) = Ŝ(t) + Ŝ(t)

g ( Ĥ1 , Ĥ) ( Ĥ(-l) 1 -Ĥ1 , Ĥ(-l) -Ĥ) = 1 n -1 t 0 d Ĥ1 (u) Ĥ(u) - 1 n -1 I(T l ≤ t, ∆ l = 1) Ĥ(T l ) - 1 n -1 t 0 d Ĥ1 (u) Ĥ(u) + 1 n -1 t 0 I(T l ≥ u)d Ĥ1 (u) ( Ĥ(u)) 2 , = - 1 n -1 I(T l ≤ t, ∆ l = 1) Ĥ(T l ) + 1 n -1 t 0 I(T l ≥ u)d Ĥ1 (u) ( Ĥ(u))
T l ∧t 0 d Ĥ1 (u) ( Ĥ(u)) 2 - I(T l ≤ t, ∆ l = 1) Ĥ(T l ) + o P (1) = Ŝ(t) + Ŝ(t) t 0 I(u ≤ T l ) Ĥ(u) d Λ(u) - t 0 dN l (u) Ĥ(u) + o P (1) = Ŝ(t) - t 0 d Ml (u) Ĥ(u) + o P (1)•
The approximation for the RMST is then obtained by directly integrating the previous equation as it actually follows that the convergence holds in the Skorohod space D[0, τ ] [see [START_REF] Richard | Lectures on survival analysis[END_REF][START_REF] Aad | Weak convergence[END_REF] and therefore the convergence holds uniformly with respect to t ∈ [0, τ ].

Proof of Proposition 2

Starting with Equation (4) we will first prove that

√ n ε n - √ n -1 ε (-l) n
tends to 0 in probability as n tends to infinity. Set

Ĩn = - 1 n n i=1 ∇ 2 log f (X i ; α), Ĩ(-l) n = - 1 n -1 i =l ∇ 2 log f (X i ; α),
where α lies between α and α 0 . We have:

√ n ε n - √ n -1 ε (-l) n = Ĩ-1 n -I -1 n i=1 ∇ log f (X i ; α 0 ) - Ĩ(-l) n -1 -I -1 n i =l ∇ log f (X i ; α 0 ) = Ĩ(-l) n -1 -I -1 ∇ log f (X l ; α 0 ) + Ĩ-1 n -Ĩ(-l) n -1 n i=1 ∇ log f (X i ; α 0 ).
Clearly for I positive definite, Ĩ(-l) n -1

-I -1 ∇ log f (X l ; α 0 ) tends to 0 in probability.

Since n i=1 ∇ log f (X i ; α 0 )/n tends to E(∇ log f (X; α 0 )) = 0 in probability, we just need to prove that Ĩ-1

n - Ĩ(-l) n -1 = O P (1/n) to conclude the proof. Write: Ĩ-1 n -Ĩ(-l) n -1 = Ĩ-1 n ( Ĩ(-l) n -Ĩn ) Ĩ(-l) n -1
.

From the law of large numbers, Ĩ-1 n and Ĩ(-l) n -1

tend towards I -1 in probability and

Ĩ(-l) n -Ĩn = - 1 n(n -1) i =l ∇ 2 log f (X i ; α) + 1 n ∇ 2 log f (X l ; α) = O P (1/n) . This proves that nα -(n -1)α (-l) = α 0 + I -1 ∇ log f (X l ; α 0 ) + o P (1). (9)
Using the consistency of Λ(t; α) towards Λ(t; α 0 ) from standard maximum likelihood theory, we now write a Taylor expansion for the cumulative hazard function around α 0 :

Λ(t; α) = Λ(t; α 0 ) + ( α -α 0 ) ∇Λ(t; α 0 ) + 1 2 (α -α 0 ) ∇ 2 Λ(t; α)(α -α 0 ), (10)
where α lies between α and α 0 . We also write a Taylor expansion for the function x → exp(-x) around 0:

exp(-(Λ(t; α) -Λ(t; α 0 ))) = 1 -Λ(t; α) -Λ(t; α 0 ) + 1 2 Λ(t; α) -Λ(t; α 0 ) 2 - 1 6 e ξn Λ(t; α) -Λ(t; α 0 ) 3 ,
with ξ n tends to 0 in probability as n tends to infinity. This can be rewritten as:

S(t; α) = S(t; α 0 ) + S(t; α 0 ) -Λ(t; α) -Λ(t; α 0 ) + 1 2 Λ(t; α) -Λ(t; α 0 ) 2 + o P (1/n),
using the fact that √ n(Λ(t; α) -Λ(t; α 0 )) converges in distribution towards a centred gaussian variable with finite variance from standard results on maximum likelihood theory and the deltamethod. As a result,

nS(t; α) -(n -1)S(t; α(-l) ) = S(t; α 0 ) + A n,1 + A n,2 + o P (1), (11)
where

A n,1 = -S(t; α 0 ) n Λ(t; α) -Λ(t; α 0 ) -(n -1) Λ(t; α(-l) ) -Λ(t; α 0 ) , A n,2 = S(t; α 0 ) n Λ(t; α) -Λ(t; α 0 ) 2 -(n -1) Λ(t; α(-l) ) -Λ(t; α 0 ) 2 1 2 •
Under the mixed-case of interval-censored and exact data, we can directly write the loglikelihood as the sum between the log-likelihood of strictly interval-censored observations and the log-likelihood of exact observations. For the latter part see [START_REF] Aalen | Survival and Event History Analysis[END_REF]. Recall that X i = (L i , R i ) and f (X i ; α) denotes the density of the observations with parameter α evaluated at X i . For strictly interval-censored data (L i = R i ), the log-likelihood (α) can be written as (see [START_REF] Sun | The statistical analysis of interval-censored failure time data[END_REF] or [START_REF] Bouaziz | Regression modelling of intervalcensored data based on the adaptive-ridge procedure[END_REF])

(α) = n i=1 log f (X i ; α) = n i=1 -(1 -∆ i )Λ(L i ; α) + ∆ i log 1 -exp Λ(L i ; α) -Λ(R i ; α) -Λ(L i ; α) ,
where we used the notation ∆ i = I(R i < +∞) to denote uncensored observations. The k th component of the score vector is equal to:

∂ (α) ∂α k = n i=1 ∂ log f (X i ; α) ∂α k = n i=1 -(c k ∧ L i -c k-1 )I(c k-1 ≤ L i ) + ∆ i (c k ∧ R i -c k-1 )I(c k-1 ≤ R i ) -(c k ∧ L i -c k-1 )I(c k-1 ≤ L i ) 1 -exp Λ(L i ; α) -Λ(R i ; α) × exp Λ(L i ; α) -Λ(R i ; α) • (12)
The k × k component of the Hessian matrix is equal to:

∂ 2 (α) ∂α k ∂α k = - n i=1 ∆ i (c k ∧ R i -c k-1 )I(c k-1 ≤ R i ) -(c k ∧ L i -c k-1 )I(c k-1 ≤ L i ) 1 -exp Λ(L i ; α) -Λ(R i ; α) × (c k ∧ R i -c k -1 )I(c k -1 ≤ R i ) -(c k ∧ L i -c k -1 )I(c k -1 ≤ L i ) × exp Λ(L i ; α) -Λ(R i ; α) + (c k ∧ R i -c k-1 )I(c k-1 ≤ R i ) -(c k ∧ L i -c k-1 )I(c k-1 ≤ L i ) 1 -exp Λ(L i ; α) -Λ(R i ; α) 2 × (c k ∧ R i -c k -1 )I(c k -1 ≤ R i ) -(c k ∧ L i -c k -1 )I(c k -1 ≤ L i ) × exp 2 Λ(L i ; α) -Λ(R i ; α) • (13)
The Fisher information is equal to the expectation of minus the Hessian matrix divided by n.

Looking at its expression, we directly see that a necessary condition for the Fisher information to be positive definite is to assume that

P(∆ = 1, [L, R] ∩ (c k-1 , c k ] = ∅) > 0, ∀k = 1, . . . , K.
Another important condition for the model to be identifiable is to assume that E α 0 [f (X; α)] has a unique maximum with respect to α, equal to α 0 , where the notation E α 0 means the expectation is taken with respect to the true parameter α 0 . However, it is clear from Equation ( 12) that E α 0 [∂f (X; α)/∂α] cannot vanish if P(L > c k-1 ) = 0. Therefore, a second necessary condition for the model to be identifiable is to assume

P(L > c k-1 ) > 0, ∀k = 1, . . . , K.
Those two conditions have opposite effects on the estimation method if they are violated. In case the first one is not valid for a given k then it will not be possible to compute the corresponding estimator αk from the Newton-Raphson algorithm (since the Hessian will not be invertible) while using the EM algorithm (which does not involve the Score vector nor the Hessian matrix), the estimator αk will become smaller at each iteration step until eventually reaching the value 0. This situation can be numerically resolved in the latter case, by simply setting the iterated estimate αk to 0 when it reaches a value below a fixed threshold. However this situation is problematic for the computation of the pseudo-values. This can be easily seen by recalling that pseudo-values should average to the initial estimator. In Proposition 2 this simply follows from the fact that n l=1 ∇ log f (X l ; α) = 0 from regularity conditions for maximum likelihood estimation. However, the k th component of the score vector will never vanish if the first condition is not valid, leading to incorrect pseudo-values.

On the other hand, if the second condition is not valid for a given k, the algorithm will attempt to minimise the term exp(-Λ(R i ; α)) from Equation ( 12) and as a consequence the corresponding estimator αk will become larger at each iteration step of the EM algorithm, diverging to infinity.

Finally, note that if the log-likelihood only include exact observations L = R, the conditions then translate to P(c k-1 < L < c k ) > 0, ∀k = 1, . . . , K.

Implementation of the pseudo-observations for the survival function and the RMST in the pch model

In this section we provide the precise expression of the terms involved in Proposition 2 for the pch model. We have S(t; α) = exp - while the expression of the gradient of the density ∇ log f (X l ; α) is given by the term between brackets in Equation ( 12) and Î is equal to minus the Hessian matrix (see Equation ( 13)) divided by n.

For the integrated version we need to precise how to compute the integral between 0 and τ of S(t; α) and the integral between 0 and τ of S(t; α)∇Λ(t; α). We first notice that S(t; α) = exp - Set A 1 = 0 and for l ≥ 2, define

A l = - l-1 k=1 α k (c k -c k-1 ) + c l-1 α l .
and αz is on the real line between α 0 and α z . Then, since α z maximises with respect to α the expected log-likelihood E(log f (X l ; α) | Z l = z)) we have that E(∇ log f (X l ; α z ) | Z l = z)) = 0. Then, from a Taylor expansion around α 0 we obtain:

0 = E(∇ log f (X l ; α 0 ) | Z l = z)) + (α z -α 0 ) E(∇ 2 log f (X l ; αz ) | Z l = z)),
where αz is on the real line between α 0 and α z . As a result, we have:

α z -α 0 = I -1 E(∇ log f (X l ; α 0 ) | Z l = z) + (-E(∇ 2 log f (X l ; αz ) | Z l = z))) -1 -I -1 E(∇ log f (X l ; α 0 ) | Z l = z).
Gathering all parts, we have proved that

Proof of Proposition 4

Let X i = (T i , ∆ i ) and assume the pch model for λ. In this model, the cumulative hazard function evaluated at T i , is equal to

Λ(T i ; α) = K k=1 α k (c k ∧ T i -c k-1 )I(c k-1 ≤ T i ).
However, since we want to prove the result when the mesh of the partition 0 = c 0 < c 1 < • • • < c K = +∞ tends to zero, we can write without loss of generality that there exists a δ > 0, such that for any partition whose mesh is less than δ we have for all T i , i = 1, . . . , n,

Λ(T i ; α) = K k=1 α k (c k -c k-1 )I(c k ≤ T i ).
Therefore, for a partition 0 = c 0 < c

1 < • • • < c K = +∞ such that max k |c k -c k-1 | < δ, log f (X i ; α) = ∆ i K k=1 log(α k )I k (T i ) - K k=1 α k (c k -c k-1 )I(c k ≤ T i ), ∂ ∂α k log f (X i ; α) = ∆ i α k I k (T i ) -(c k -c k-1 )I(c k ≤ T i ), ∂ 2 ∂α 2 k log f (X i ; α) = - ∆ i α 2 k I k (T i ), ∂ 2 ∂α k ∂α k log f (X i ; α) = 0, for k = k .
We therefore have that the Fisher information I is a diagonal matrix whose kth element, k = 1, . . . , K, is equal to E(∆ i I k (T i ))/(α 0 k ) 2 , where α 0 k is the kth component of α 0 . Also under standard maximum likelihood regularity conditions, the true parameter α 0 verifies that E ∂ ∂α k log f (X i ; α 0 k ) = 0, which is equivalent to

α 0 k = E(∆ i I k (T i )) (c k -c k-1 )P(T i ≥ c k )
•

Now, we have that E(∆ l I k (T l ))/(c k -c k-1 ) tends to H 1 (c k-1 ) as (c k -c k-1 ) tends to 0, where H 1 (t) = P(T ≤ t, ∆ = 1). Therefore, as the limit (c k -c k-1 ) goes to 0, α 0 k tends to λ(c k ), the true hazard rate evaluated at c k . In other words, the parametric hazard λ(t; α 0 ) tends to the true hazard function λ(t) as (c k -c k-1 ) goes to 0.

Then, ∇Λ(t; α 0 ) is a vector whose kth component is equal to (c k -c k-1 )I(c k ≤ t). As a consequence,

∇Λ(t; α 0 ) I -1 ∇ log f (X l ; α 0 ) = K k=1 I(c k ≤ t) (α 0 k ) 2 (c k -c k-1 ) E(∆ l I k (T l )) ∆ l α 0 k I k (T l ) -(c k -c k-1 )I(c k ≤ T l ) .
In the first term of this equation, there can only be one interval of the form [c k-1 , c k ] that contains T l and therefore the sum over k is not zero for only one value of k. As the limit of (c k -c k-1 ) goes to 0, T l gets closer to c k and c k-1 and α 0 k tends to λ(T l ). More precisely, lim

c k -c k-1 →0 K k=1 I(c k ≤ t) (α 0 k ) 2 (c k -c k-1 ) E(∆ l I k (T l )) ∆ l α 0 k I k (T l ) = ∆ l I(T l ≤ t)λ(T l ) H 1 (T l ) •
On the other hand, the second term of the equation is simply a Riemann integral. We have: We now want to compute the conditional expectation with respect to Z l of this quantity. Let H Z 1 (t) = P(T ≤ t, ∆ = 1 | Z), H Z (t) = P(T ≥ t | Z) and write

E ∆ l I(T l ≤ t) H(T l ) - T l ∧t 0 dH 1 (u) (H(u)) 2 | Z l = t 0 dH Z l 1 (u) H(u) - t 0 H Z l (u)dH 1 (u) (H(u)) 2 = t 0 (H Z l 1 (u)) H Z l (u) - (H 1 (u)) H(u) H Z l (u) H(u) du = t 0 (λ Z l (u) -λ(u)) H Z l (u) H(u) du,
where λ Z l represents the conditional hazard function. Setting q Z l (u) = (λ Z l (u)-λ(u))H Z l (u)/H(u) as in [START_REF] Jacobsen | A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations[END_REF] we obtain from their proofs in Section A.2, Equation [START_REF] Groeneboom | Information bounds and nonparametric maximum likelihood estimation[END_REF], that -S Z l /S is a primitive of q Z l , where S Z l is the conditional survival function of T * . As a consequence, t 0 (λ Z l (u) -λ(u)) dH Z l (u) H(u) = 1 -S Z l (t) S(t) •

Proposition 4 .

 4 Under the context of right-censored data, if T * follows the pch model with cuts

Figure 1 :

 1 Figure 1: Distribution of time to emergence of the tooth number 14. On the left: estimated hazard function. On the right: estimated survival function. Those estimates were obtained from the pch model with cuts equal to 7.6, 8.4, 9 and 10.

α

  k (t ∧ c k -c k-1 )I(c k-1 ≤ t) ∂Λ(t; α) ∂α k = (c k ∧ t -c k-1 )I(c k-1 ≤ t),

α

  k (c k -c k-1 )I(c k ≤ t) exp -K k=1 α k (t -c k-1 )I(c k-1 ≤ t ≤ c k ) , c k -c k-1 )I(c k ≤ t) exp -α l (t -c l-1 ) dt I(τ > c l-1 ).

E 0 S(t; α 0 )dt + τ 0 ( 0 (

 000 (T * l ∧ τ | Z l = z) = τ ∇S(t; α 0 )) dtI -1 E(∇ log f (X l ; α 0 ) | Z l = z) + R 1,z + R 2,z , with R 2,z = τ ∇S(t; α 0 )) dt (-E(∇ 2 log f (X l ; αz ) | Z l = z))) -1 -I -1 E(∇ log f (X l ; α 0 ) | Z l = z).

Finally 0 S

 0 , writing S(t; α 0 ) = exp(-Λ(t; α 0 )), we directly obtain τ 0 ∇S(t; α 0 )dt = -τ (t; α 0 )∇Λ(t; α 0 )dt, which concludes the proof.

lim c k -c k- 1

 1 →0 K k=1 I(c k ≤ t) (α 0 k ) 2 (c k -c k-1 ) 2 E(∆ l I k (T l )) I(c k ≤ T l ) = T l ∧t 0 (λ(u)) 2 H 1 (u) du•Noticing that λ(t) = H 1 (t)/H(t), with H(t) = P(T ≥ t), we therefore have proved limc k -c k-1 →0 ∇Λ(t; α 0 ) I -1 ∇ log f (X l ; α 0 ) = ∆ l I(T l ≤ t) H(T l ) -T l ∧t 0 dH 1 (u) (H(u)) 2 •

Table 1 :

 1 Simulation results for the estimation of β in the RMST model (

	Approximated formula

Table 2 :

 2 Simulation results for the estimation of β in the RMST model (

	Approximated formula

Table 3 :

 3 Cox Model for the time to emergence of the tooth 14 with the covariates gender and number of decayed or missing deciduous first molars due to caries among teeth 54, 64, 74, 84. The time has been discretised at t = 8, 9, 10, 11, 12 in the pseudo-regression approach. se represents the standard estimate of the regression parameter. Nb of decayed molars 0.2808 0.0260 < 10 -15 0.1080 0.0194 2.7 × 10 -8

			time = 9			time = 12
	Covariates	effect	se	p-value effect	se	p-value
	Intercept	-2.8458 0.0761 < 10 -15 0.8777 0.0501	< 10 -15
	Gender (1 = boy)	0.2978 0.0819	0.0003	0.5284 0.0599	< 10 -15

Table 4 :

 4 Two logistic models for the probability that the emergence of the tooth 14 occurred before time 9 and 12 with respect to the covariates gender and number of decayed or missing deciduous first molars due to caries among teeth 54, 64, 74, 84. se represents the standard estimate of the regression parameter.

			τ = 9			τ = 12
	Covariates	effect	se	p-value	effect	se	p-value
	Intercept	8.9851	0.0047	< 10 -15	10.8755 0.0306 < 10 -15
	Gender (1 = boy)	-0.0097 0.0066	0.1422	-0.3336 0.0361 < 10 -15
	Nb of decayed molars -0.0180 0.0024 8.1379 × 10 -14 -0.1303 0.0120 < 10 -15

Table 5 :

 5 Restricted Mean Survival Time Model for the time to emergence of the tooth 14 with the covariates gender and number of decayed or missing deciduous first molars due to caries among teeth 54, 64, 74, 84. Two values of τ are analysed in Equation[START_REF] Zhang | Estimating differences in restricted mean lifetime using observational data subject to dependent censoring[END_REF]. se represents the standard estimate of the regression parameter.

Acknowledgments

We thank the reviewers for their very constructive criticisms and comments that have helped improve the paper. We also thank Per Kragh Andersen and Terry Therneau for their very interesting discussions about the connection between the infinitesimal jackknife and the Von Mises approximation for right-censored data.

Appendix

Proof of Proposition 1

Introduce H 1 (t) = P(T ≤ t, ∆ = 1), H 0 (t) = P(T ≤ t, ∆ = 0), H(t) = P(T ≥ t) and their empirical counterparts, Ĥ1 (t) = i I(T i ≤ t, ∆ i = 1)/n, Ĥ0 (t) = i I(T i ≤ t, ∆ i = 0)/n, Ĥ(t) = i I(T i ≥ t)/n. Let ψ(A)(s, t] = s<u≤t (1 + dA(u)), such that ψ(Λ)(0, t] = S(t), where Λ(t) is the cumulative hazard function and ψ( Λ)(0, t] = Ŝ(t). We have the following Von-Mises expansion [see [START_REF] Richard | Lectures on survival analysis[END_REF][START_REF] Aad | Weak convergence[END_REF]:

We now derive a Von-Mises expansion for Λ(-l) (t) -Λ(t). The cumulative hazard function and its estimator can be defined as functions of H, H 1 and of Ĥ1 , Ĥ respectively where Λ(t) = We start with the A n,2 term. From Equation [START_REF] Turnbull | The empirical distribution function with arbitrarily grouped, censored and truncated data[END_REF] we have:

Using the consistency of α -α 0 and the asymptotic normality of √ n(α -α 0 ) from standard maximum likelihood theory, each of the last two terms in the above equation tends to 0 in probability as n tends to infinity. Therefore

where the last two lines were derived from Equation ( 9) and R n , R n both tend to 0 in probability. The consistency of α and α(-l) shows that A n,2 = o P [START_REF] David | Restricted mean life with covariates: modification and extension of a useful survival analysis method[END_REF]. We now study the term A n,1 . From Equation [START_REF] Turnbull | The empirical distribution function with arbitrarily grouped, censored and truncated data[END_REF],

Using similar arguments as before, the last two lines of this Equation tend to 0 in probability from Equation [START_REF] Zhao | Deep neural networks for predicting restricted mean survival times[END_REF] and from the consistency of α and α(-l) . Finally, using again Equation ( 9)

This equality combined with Equation [START_REF] Groeneboom | Information bounds and nonparametric maximum likelihood estimation[END_REF] give

The final result of Proposition 2 is obtained by simply replacing each quantity by its consistent estimator. Integrating the equation in Proposition 2 directly yields the approximation for the RMST. By careful examination of the remainder term, we directly see that its integral over [0, τ ] is also o P (1).

Log-likelihood, score vector and Hessian matrix in the piecewise constant hazard model

In this section, we study the parametric piecewise constant hazard model defined as follows:

For the first term we now have:

For the second term we have:

From the previous calculation on the first term, we easily see that on the one hand

On the other hand, we have:

where the last equation was obtained using integration by parts. Gathering all elements allows to implement the second equation in Proposition 2.

Proof of Proposition 3

First, from a Taylor expansion of τ 0 S(t; α z )dt around τ 0 S(t; α 0 )dt, we obtain:


Florian Naudet

Eiko Fried

Lisa Cosgrove

Erick Turner

Alain Braillon

Ioana Cristea

Psychedelic drugs: more emphasis on safety issues Florian Naudet, University of Rennes, Rennes, France. Eiko I. Fried, Leiden University, Leiden, The Netherlands. Lisa Cosgrove, University of Massachusetts Boston, Boston, Massachusetts, USA. Erick Turner, Oregon Health & Science University, Portland, Oregon, USA. Alain Braillon, Amiens, France. Ioana A. Cristea, University of Padua, Padua, Italy.

Psychedelic drugs have safety issues that are linked with effects such as hallucinations, paranoia and mood changes. In our view, your editorially independent articles on these drugs -sponsored by the biotechnology firm Atai Life Sciences -should have given more attention to such concerns to avoid serving an industry-friendly narrative that could do more harm than good (see www.nature.com/collections/edcbibiadj).

For example, the article on psilocybin for pain relief (Nature 609, S100-S102; 2022) notes that researchers are trying to avoid the mistakes that contributed to the opioid crisis. But there are strong parallels: prescription of psychedelics and of opioids is rising as a result of aggressive marketing, whereas the evidence for effectiveness comes mostly from small, short-term trials involving highly selected patients in secure settings, without adequate comparators. And attention to safety issues such as tolerance and drug dependence is crucial, particularly over the long term.

The US Food and Drug Administration has designated many psychedelics as a 'breakthrough therapy', in effect lowering the bar for approval. Esketamine, for example, was controversially approved for treatment-resistant depression without adequate monitoring to assess its long-term harm-to-benefit ratio (E. H. Turner Lancet Psychiatry 6, 977-979; 2019).

Competing Interests