J.-F Aujol

Ch Dossal
email: charles.dossal@insa-toulouse.fr

H Labarrière

A Rondepierre
email: aude.rondepierre@insa-toulouse.fr

‡ Jean

FISTA restart using an automatic estimation of the growth parameter

Keywords: FISTA, restart, convex optimization, Lojasiewicz property, convergence rate, growth parameter

 show that restarting FISTA could ensure a fast convergence for this class of functions without requiring any knowledge on the growth parameter. We improve these restart schemes by providing a better asymptotical convergence rate and by requiring a lower computation cost. We present numerical results emphasizing the efficiency of this method.

Introduction

Let N ∈ N * and consider a composite optimization problem:

min x∈R N F (x). (1)
The objective function F : R N → R ∪ {+∞} belongs to the class H L of composite functions: F = f + h where f is a convex, differentiable function having a L-Lipschitz gradient and h is a convex function whose proximal operator is known.

The set of minimizers of F denoted by X * is assumed to be non empty.

For the class H L of composite functions, a classical minimization algorithm is the Forward-Backward algorithm (FB) that produces a sequence (x k ) k 1 ensuring

F (x k ) -F * = O 1
k [START_REF] Chen | Convergence rates in Forward-Backward splitting[END_REF] where F * = inf F . Under the same assumptions, FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], based on the ideas of acceleration of Nesterov [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o(1/k 2 )[END_REF], ensures a better asymptotic bound:

F (x k ) -F * = O 1 k 2 .
In many optimization problems, as for instance in statistics or in image processing, the objective function F naturally satisfies some additional geometry assumptions such as a quadratic growth condition:

∀x ∈ R N , µ 2 d(x, X * ) 2 F (x) -F * , (2)
which allows to reach better decay rates. Note that in the convex setting, the quadratic growth condition ( 2) is equivalent to a Lojasiewicz property with an exponent 1 2 [START_REF] Lojasiewicz | Une propriété topologique des sous-ensembles analytiques réels[END_REF][START_REF] Lojasiewicz | Sur la géométrie semi-et sous-analytique[END_REF][START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF]. This set of functions includes the set of strongly convex functions but it is much larger. For example, it contains functions associated to the mean square problem and the LASSO [15, Corollary 9]. Note that the uniqueness of the minimizer of F is not required.

On this set of functions, FB reaches an exponential decay rate: O e -κk where κ := µ L is the ratio between the growth parameter µ defined in [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF] and the Lipschitz constant L of ∇f [START_REF] Garrigos | Convergence of the Forward-Backward algorithm: Beyond the worst case with the help of geometry[END_REF]. It turns out that this exponential decay cannot be observed in many numerical experiments because the condition number κ = µ L 1 can be very small, especially in large dimension problems. Most of the time, up to a large accuracy, the quadratic bound of FISTA is numerically better than the theoretical exponential decay of FB. Indeed the number of iterations required to reach a precision ε is proportional to L µ for FB while it is proportional to L µ for FISTA. Assuming that F has a unique minimizer, the variation of Polyak's Heavy Ball method [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF] introduced in [START_REF] Aujol | Convergence rates of the Heavy-Ball method with Lojasiewicz property[END_REF] reaches a better decay rate:

O e -(2- √ 2) √ κk .
Unfortunately, this scheme requires an accurate a priori estimation of the growth parameter µ to get this fast decay. Moreover, in many situations, µ is small and unknown. An incorrect estimation of µ may significantly reduce the speed of the algorithm.

In this paper, we introduce a new algorithm to minimize efficiently composite functions of the class H L having additionally the quadratic growth property [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF]. More precisely this new algorithm is based on an original restart rule of FISTA ensuring:

F (x k ) -F * = O e -1 12 √ κk (3)
without any a priori knowledge of µ, and relies on an iterative estimation of µ comparing some values F (x k ). The restart rule is inspired by the one proposed by Alamo et al [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF]. Compared to [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF], we improve the decay rate and reduce the number of estimations of F (x k ) during the iterations (notice that these evaluations may heavily impact the numerical cost and calculation time). The article is organized as follows. Section 2 provides a state of the art of this optimization problem. Section 3 is then devoted to the definition and the properties of our restart algorithm. We also propose some numerical experiments and comparisons. The proofs of the results of Section 3 are postponed to Section 4.

State of the art 2.1 Framework

We introduce some notations. Given a differentiable function f , ∇f denotes the gradient of f . Given a convex lower semicontinuous function h, ∂h denotes the convex subgradient of h defined by:

∀x ∈ R N , ∂h(x) = s ∈ R N | ∀y ∈ R N , h(y) h(x) + s, y -x . (4)
Given a vector x, x and x 1 respectively denote its Euclidean norm and its 1-norm. We define the proximal operator of h : R N → R ∪ {+∞} a convex lower semicontinuous function as follows:

∀x ∈ R N , prox h (x) = argmin y∈R N h(y) + 1 2 y -x 2 . (5)
In this paper we focus on the class H L of composite functions defined as: F = f +h where f is a convex, differentiable function having a L-Lipschitz gradient and h is a convex function whose proximal operator is known. The set X * of minimizers of F is assumed to be non empty.

Definition 1 (Quadratic growth condition G 2 µ ). Let F : R N → R ∪ {+∞} be a proper lower semicontinuous convex function with a non empty set of minimizers X * . Let F * = inf F . The function F satisfies a quadratic growth condition G 2 µ for some µ > 0 if it satisfies:

∀x ∈ R N , µ 2 d(x, X * ) 2 F (x) -F * . (6)
Classically the quadratic growth condition G 2 µ can be seen as a relaxation of the strong convexity. In the convex setting, the condition G 2 µ is equivalent to a global Lojasiewicz property with an exponent 1 2 [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF][START_REF] Garrigos | Convergence of the Forward-Backward algorithm: Beyond the worst case with the help of geometry[END_REF]: Lemma 1. Let F : R N → R ∪ {+∞} be a proper lower semicontinuous convex function with a non empty set of minimizers X * . Let F * = inf F . If F satisfies a quadratic growth condition G 2 µ for some µ > 0, then F has a global Lojasiewicz property with an exponent 1 2 :

∀x ∈ R N , 2µ (F (x) -F * ) d(0, ∂F (x)) 2 . ( 7
)
To set some definitions we say that:

Definition 2.
An optimization algorithm provides a fast exponential decay on the class of functions F satisfying a growth condition G 2 µ with parameter µ and having a L-Lipschitz gradient if there exists a constant K > 0 independent of µ and L such that the sequence (x k ) k∈N provided by this algorithm satisfies

F (x k ) -F * = O e -K √ µ L k .
Definition 3. An optimization algorithm provides a low exponential decay on the class of functions F satisfying a growth condition G 2 µ with parameter µ and having a L-Lipschitz gradient if there exists a constant K > 0 independent of µ and L such that the sequence (x k ) k∈N provided by this algorithm satisfies

F (x k ) -F * = O e -K µ L k .
In most practical cases the function F is ill-conditioned i.e µ L. The difference between a fast and a low exponential decay is then significant since µ L µ L .

Literature review

Let L > 0. Let F : R n → R be a composite function of the class H L satisfying additionally a growth condition for some parameter µ > 0, namely µ-strong convexity or the quadratic growth condition G 2 µ . The Lipschitz constant L is assumed to be known.

A classical method to minimize such a function is the Forward-Backward algorithm which is an adaptation of Gradient Descent method. If F is µ-strongly convex or satisfies G 2 µ for some µ > 0, this scheme provides a low exponential decay

F (x k ) -F * = O e -µ
L k [START_REF] Garrigos | Convergence of the Forward-Backward algorithm: Beyond the worst case with the help of geometry[END_REF]. Polyak introduces in [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF] the Heavy-Ball method which ensures a fast convergence for µ-strongly convex functions that are twice differentiable. There exist many variations of this scheme and [START_REF] Aujol | Convergence rates of the Heavy-Ball method with Lojasiewicz property[END_REF] proposes a Heavy-Ball method adapted to composite functions satisfying G 2 µ for some µ > 0. This scheme provides a fast exponential decay

O e -(2- √ 2) √ µ L k
for this set of functions under a uniqueness of minimizers assumption. This method requires a prior estimate of the growth parameter µ.

Restarting schemes of FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] are efficient methods in this setting. FISTA is an inertial first order algorithm adapted from Nesterov accelerated gradient introduced in [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o(1/k 2 )[END_REF] which ensures that F (x k )-F * = O(k -2 ). Inertia generated in this scheme allows a fast convergence but also produces oscillations. Restarting FISTA is equivalent to set inertia to zero which helps to reduce oscillating behavior. A classical strategy introduced in [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o(1/k 2 )[END_REF] is to restart the algorithm at regular intervals. [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF] show that restarting Nesterov accelerated gradient every 2e L µ iterations ensures that

F (x k ) -F * = O e -1 e √ µ L k
for µ-strongly convex functions. This restart scheme and its convergence rate can be extended to FISTA and composite functions satisfying G 2 µ [START_REF] Candes | Adaptive restart for accelerated gradient schemes[END_REF][START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF]. Note that an estimate of the growth parameter µ is required. This may be restrictive in numerical experiments since the growth parameter is rarely known.

O'Donoghue and Candès propose heuristic restart rules for accelerated gradient method and FISTA in [START_REF] Candes | Adaptive restart for accelerated gradient schemes[END_REF]. These rules are based on empirical observations and they give highly effective numerical results especially for functions satisfying G 2 µ . However, no improved convergence rates was theoretically found so far.

Adaptive restart schemes take advantage of each iteration to estimate the geometry of F . This approach enables to fit the parameters of the algorithm progressively. [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF][START_REF] Gonzaga | Fine tuning nesterov's steepest descent algorithm for differentiable convex programming[END_REF][START_REF] Gonzaga | An optimal algorithm for constrained differentiable convex optimization[END_REF] propose non restart schemes based on this strategy to compute an estimate of the strong convexity parameter. [START_REF] Fercoq | Adaptive restart of accelerated gradient methods under local quadratic growth condition[END_REF] provides an adaptive restart scheme for FISTA for the set of functions satisfying G 2 µ which builds a sequence of estimates of the growth parameter µ. It requires a prior estimate µ 0 and the convergence rate of the method is given by O e

- √ 2-1 2 √ e(2- √ µ µ 0 ) √ µ L k if µ 0 µ.
As this rate significantly depends on µ µ 0 this method might be less effective if µ 0 is highly overestimated. This rate is faster than O e as long as µ 0 < 4µ. [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF] and [START_REF] Alamo | Gradient based restart FISTA[END_REF] adaptive restart schemes for FISTA under a G 2 µ assumption on F where no prior estimate of µ > 0 is required. The method proposed in [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF] provides a fast exponential decay O e -1 16

Alamo et al. introduce in

√ µ L k
and it relies on computations of F (x k ) at least at half of the iterations. There exist techniques that reduce the computational cost of these evaluations. In [START_REF] Goldstein | A field guide to forward-backward splitting with a fasta implementation[END_REF] the authors introduce FASTA which is an adaptation of FISTA optimized for composite functions satisfying H L such that f (x) = f (Ax) for some function f and some operator A. In this setting, this scheme reduces the computational cost of f and therefore of F . However in some common cases computing F is expensive despite these strategies. In these cases additional calculations may significantly slow down the method. The scheme introduced in [START_REF] Alamo | Gradient based restart FISTA[END_REF] is based on gradient information and therefore it does not require to compute F . It provides the convergence rate O e

- 1 4e(1+ √ µ+1) µ L k which
can be seen as a low exponential decay as e

- 1 4e(1+ √ µ+1) µ L k > e -1 8e µ L k for all k > 0.
O       e - √ 2-1 2 √ e(2- √ µ µ 0 ) √ µ L k      
Requires an estimate µ 0 of µ

FISTA restart by

Alamo et al.

Alamo et al. [2]

O e -1 16

√ µ L k Requires to compute F regularly Gradient based FISTA restart Alamo et al. [1] O e - 1 4e(1+ √ µ+1) µ L k - Adaptive restart Theorem 1 O e -1 12 √ µ L k -
We propose a restart scheme of accelerated gradient method and FISTA which provides a fast exponential decay in O e -1 12 √ µ L k . This algorithm does not require to estimate the growth parameter µ. Moreover, the value of F (x k ) is computed for a reduced number of iterations. As a consequence, this method has a fast convergence rate in terms of iterations and computational time as well.

Contributions

In this section we introduce a restart scheme for FISTA in order to minimize a composite function F satisfying H L and the growth condition G 2 µ for some L > 0 and µ > 0. We give a convergence rate and we describe the underlying strategy.

Adaptive FISTA restart scheme

Let us introduce some notations to simplify the writing of the algorithm. Let:

y + = prox 1 L h (y - 1 L ∇f (y)) (8)
be the vector given by a step of the Forward-Backward algorithm on y with s = 1 L . The composite gradient mapping g is then defined as follows:

g(y) = L(y -y + ). ( 9
)
Given initial condition z ∈ R N and a number of iterations n, FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] can be written as Algorithm 1.

Algorithm 1 : FISTA Require: z ∈ R N , n ∈ N y 0 = x 0 = z, k = 0 repeat k = k + 1 x k = y + k-1 y k = x k + k-1 k+2 (x k -x k-1 ) until k n return r = x k
We introduce a restart scheme which aims to accelerate the convergence of FISTA. Algorithm 2 relies on a few sequences:

• the sequence (r j ) j∈N corresponds to the iterates of the algorithm. For all j > 0, r j is the output of the jth execution of FISTA.

• the sequence (n j ) j∈N refers to the number of iterations of FISTA following the j-th restart. Namely, for all j 0 we have:

r j+1 = FISTA(r j , n j ),
• the sequence (μ j ) j 2 estimates the growth parameter µ at each restart. This estimation is built from the known convergence results of FISTA and considering a comparison of the cost function F computed at three iterates.

For all j 2, the number of iterations n j is defined according to n j-1 , μj and the predefined parameter C > 0. The algorithm verifies if a condition called doubling condition is fulfilled. If the condition holds true, then n j-1 is considered too small and n j is set to 2n j-1 . In the other case, the number of iterations is not increased. The exit condition g(r j ) ε is motivated by the following lemma which is proven in Section 4.5: Lemma 2. Let F be a function satisfying H L and G 2 µ for some L > 0 and µ > 0. Then for all x ∈ R N we have:

F (x + ) -F * 2 µ g(x) 2 . ( 10
)
The exit criteria combined to the inequality [START_REF] Aujol | Convergence rates of the Heavy-Ball method for quasi-strongly convex optimization[END_REF] enable the algorithm to ensure that the vector r j satisfies:

F (r + j ) -F * 2ε 2 µ , (11)
without computing F * . For a given initial condition r 0 ∈ R N , we propose the algorithm shown in Algorithm 2.

Algorithm 2 : Restart scheme Require: r 0 ∈ R N , j = 1 n 0 = 2C r 1 = FISTA(r 0 , n 0 ) n 1 = 2C repeat j = j + 1 r j = FISTA(r j-1 , n j-1 ) μj = min i∈N * i<j 4L (n i-1 + 1) 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) if n j-1 C L μj then n j = 2n j-1 else n j = n j-1 end if until g(r j ) ε return r = r j Theorem 1.
Let F be a function satisfying H L and G 2 µ for some L > 0 and µ > 0. Let (r j ) j∈N and (n j ) j∈N be the sequences provided by Algorithm 2 with parameters C > 4 and ε > 0. Then the number of iterations 1 + j i=0 n i required to guarantee g(r j )

ε is bounded and satisfies

j i=0 n i 4C log C 2 4 -1 L µ 2 log C 2 4 -1 + log 1 + 16 C 2 -16 2L(F (r 0 ) -F * ) ε 2 . ( 12
)
Corollary 1. Let F be a function satisfying H L and G 2 µ for some L > 0 and µ > 0. If C > 4 and ε > 0, then the sequences (r j ) j∈N and (n j ) j∈N provided by Algorithm 2 satisfy

F (r + j ) -F * = O       e - log C 2 4 -1 4C √ µ L j i=0 n i       . ( 13
)
Specifically, if C is chosen to maximize

log C 2 4 -1 4C
, namely C ≈ 6.38, then there exists K > 1 12 such that the sequences (r j ) j∈N and (n j ) j∈N provided by Algorithm 2 satisfy

F (r + j ) -F * = O       e -K √ µ L j i=0 n i       . ( 14
)
Corollary 1 states that Algorithm 2 provides asymptotically a fast exponential decay. This convergence rate is faster than any method when considering a function F satisfying H L and G 2 µ where the parameter µ cannot be estimated. In this setting, Forward-Backward algorithm provides a low exponential decay and FISTA has no improved theoretical convergence rate. The variation of Heavy-Ball method introduced in [8], the FISTA restart scheme introduced in [START_REF] Fercoq | Adaptive restart of accelerated gradient methods under local quadratic growth condition[END_REF] and fixed restart of FISTA require to estimate the growth parameter to ensure a fast exponential decay.

The restart scheme introduced by Alamo et al. in [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF] ensures the following rate:

F (r + j ) -F * = O       e -1 16
√ µ L j i=0 n i       , (15)
which is slightly slower than [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF]. The gradient based FISTA restart proposed in [START_REF] Alamo | Gradient based restart FISTA[END_REF] guarantees the theoretical rate O e

- 1 4e(1+ √ µ+1) µ L k
which can be considered as an low exponential decay. As a consequence, this decay may be significantly slower than the decay of Algorithm 2 when µ L µ L which is the case in numerous practical examples.

We recall that convergence results such as ( 14) and ( 15) are worst-case bounds. Consequently a faster theoretical rate does not necessarily guarantee a faster algorithm in general.

Structure of the algorithm

Algorithm 2 is an adaptive restart scheme based on FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. At each step j, FISTA is restarted for n j-1 iterations and it gives the vector r j :

r j = FISTA(r j-1 , n j-1 ).
The key parameter of such an algorithm is the sequence (n j ) j∈N . The strategy behind Algorithm 2 is to estimate the growth parameter µ > 0 satisfying (6) at each step j. The estimation of µ denoted μj is then used to define n j .

The estimation of µ is based on well-known convergence results on FISTA stated in Proposition 1 and proven in Section 4.4. This proposition is adapted from [2, Property 2] and introduces slightly more general claims leading to the estimate μj [START_REF] Kim | Optimized first-order methods for smooth convex minimization[END_REF]. Proposition 1. Let F be a function satisfying H L and G 2 µ for some L > 0 and µ > 0. Then the sequence (x k ) k∈N provided by Algorithm 1 satisfies

(i) ∀k ∈ N * , F (x k ) -F * 4L µ(k + 1) 2 (F (x 0 ) -F * ) , (16)
(ii) ∀k ∈ N * , F (x k ) F (x 0 ), (17)
(iii) ∀k ∈ N, F (x k ) -F * > γ (F (x 0 ) -F (x k )) =⇒ (k + 1) 2 < 4L µ 1 + 1 γ . (18)
Given the first claim of Proposition 1, the most direct strategy to estimate µ is to compare F (r j ) -F * and F (r j-1 ) -F * at each step j as we have:

∀j ∈ N * , F (r j ) -F * 4L µ(n j-1 + 1) 2 (F (r j-1 ) -F * ) , (19)
which is equivalent to:

∀j ∈ N * , µ 4L (n j-1 + 1) 2 F (r j-1 ) -F * F (r j ) -F * . ( 20
)
In many cases the optimal value F * is not known and µ cannot be estimated in this way. This issue can be avoided by rewriting [START_REF] Garrigos | Convergence of the Forward-Backward algorithm: Beyond the worst case with the help of geometry[END_REF] and by considering a third point r j+1 . For all j ∈ N * ,

F (r j ) -F * > F (r j ) -F (r j+1 ) = γ(F (r j-1 ) -F (r j )), (21)
where γ = F (r j ) -F (r j+1 ) F (r j-1 ) -F (r j ) . The third claim of Proposition 1 gives us that:

µ 4L (n j-1 + 1) 2 1 + 1 γ . (22)
Thus,

∀j ∈ N * , µ 4L (n j-1 + 1) 2 F (r j-1 ) -F (r j+1 ) F (r j ) -F (r j+1 ) . (23)
In this way, it is possible to get an estimation of µ at each restart for j 2 by comparing F (r j-1 ) -F (r j ) and F (r j-2 ) -F (r j ). This strategy is implemented in Algorithm 2 to build the sequence (μ j ) j 2 :

∀j 2, μj = min i∈N * i<j 4L (n i-1 + 1) 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) . (24)
This sequence satisfies the following property proven in Section 4.6:

Lemma 3. Let F be a function satisfying H L and G 2 µ for some L > 0 and µ > 0. Then the sequence (μ j ) j 2 provided by Algorithm 2 satisfies ∀j 2, μj μj+1 > µ.

(

) 25
It is worth noticing that this strategy could be adapted to other schemes than FISTA as long as theoretical convergence bounds are known. The estimate of the growth parameter indeed relies only on [START_REF] Chen | Convergence rates in Forward-Backward splitting[END_REF].

Once the growth parameter µ is estimated it is possible to set the number of iterations accordingly. It is known [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o(1/k 2 )[END_REF][START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF][START_REF] Candes | Adaptive restart for accelerated gradient schemes[END_REF] that the optimal number of iterations before a restart of FISTA is k * = 2e L µ . However, considering μ an upper estimation of µ, the growth condition G 2 μ is not satisfied by F . Consequently, setting k * = 2e L μ does not ensure a fast exponential decay. Since μj µ for all j 2, setting n j = 2e L μj may not be efficient. The strategy of Algorithm 2 is to reset n j by using a doubling condition that depends on the parameter C > 0. At each step the doubling condition n j C L μj is evaluated. If the condition is satisfied then the number of iterations n j is considered too small and it is doubled before the next restart. This process ensures that (n j ) j∈N is an increasing and bounded sequence which satisfies the following lemma. Its proof can be found in Section 4.7.

Lemma 4. Let F be a function satisfying H L and G 2 µ for some L > 0 and µ > 0. Then the sequence (n j ) j∈N provided by Algorithm 2 satisfies

∀j ∈ N, n j 2C L µ . ( 26
)
Remark 1. An alternative approach would be to multiply the number of iterations n j by γ > 1 instead of 2 when n j C L μj . In that case the upper bound of n j becomes

∀j ∈ N, n j γC L µ . ( 27
)
The fastest asymptotical convergence rate is obtained for the choice γ = 2 (Algorithm 2).

Numerical experiments

In this section we illustrate the convergence results stated in previous sections with numerical experiments. We compare Algorithm 2 to the following set of methods:

1. The Forward-backward algorithm, 2. FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF],

3. Empirical restart scheme of FISTA by O'Donoghue and Candès [START_REF] Candes | Adaptive restart for accelerated gradient schemes[END_REF] with the restart condition:

F (x k ) > F (x k-1 ), (28)
4. Empirical restart scheme of FISTA by O'Donoghue and Candès [START_REF] Candes | Adaptive restart for accelerated gradient schemes[END_REF] with the restart condition:

g(y k-1 ), x k -x k-1 > 0, ( 29

L 1 -regularized least squares problem

The problem considered reads as follows:

min x∈R N F (x) = 1 2 Ax -b 2 + ρ x 1 , (30)
where

A ∈ M N,N (R), b ∈ R N , ρ > 0 and N 1.
The function F is a composite function satisfying H L where L > 0 can easily be computed. The L 1 -regularization prevents F to be strongly convex but it satisfies the assumption G 2 µ for some µ > 0. There exists no explicit formula for the growth parameter µ. Experiments are carried out by setting A as a random matrix of size N × N . Results are provided for N = 1000 which ensures that the problem is ill-conditioned. The computational cost of F is negligible as optimization techniques such as FASTA [START_REF] Goldstein | A field guide to forward-backward splitting with a fasta implementation[END_REF] or speculative estimations [START_REF] Florea | A generalized accelerated composite gradient method: Uniting nesterov's fast gradient method and fista[END_REF] can be applied.

Figure 1 shows the decrease of the error according to the number of iterations. One can observe that Algorithm 2 and the restart scheme introduced by Alamo et al. in [START_REF] Alamo | Restart FISTA with global linear convergence[END_REF] provide similar precision with equal number of iterations. In this example additional computations of F are not expensive and consequently both methods perform similarly. This figure also illustrates the low exponential decay provided by Gradient based restart [START_REF] Alamo | Gradient based restart FISTA[END_REF] as this scheme is significantly slower. Note that both empirical restart schemes provide significantly fast decay.

Inpainting

Consider an image x 0 and a masking operator M . Let y = M x 0 be the damaged version of x 0 . The objective is to get an approximation of x 0 knowing y and M . This problem can be written as follows:

min x F (x) = 1 2 M x -y 2 + λ T x 1 , (31)
where T is an orthogonal transformation ensuring that T x 0 is sparse. In this example, x 0 is piecewise smooth so T is set as an orthogonal wavelet transform. In this setting, the objective function F satisfies H L and G 2 µ where L = 1 and µ > 0 cannot be computed directly. This inpainting problem [START_REF] Candes | Adaptive restart for accelerated gradient schemes[END_REF] is an example of a highly expensive function F to compute. As shown in Figure 3 and Figure 4, this high cost significantly penalizes restart schemes that rely on regular computations of F . It appears that Algorithm 2 is still efficient in this situation as its additional computations of the cost function are reduced.

Proofs

Sketch of proof

As the proof of Theorem 1 is technical, it is split into the following parts:

1. We show that there is at least one doubling step every T iterations for a well-chosen T .

(a) We suppose that there is no doubling step from j = s + 1 to j = s + T .

(b) We exhibit the geometrical decrease of (F (r j-1 )-F (r j )) j∈ s+1,s+T which represents the gain of the j-th execution of FISTA.

(c) We apply Lemma 5 to show that we can find an upper bound of g(r j-1 ) which depends on F (r j-1 ) -F (r j ) for all j ∈ s + 1, s + T .

(d) We exploit the geometrical decrease (F (r j-1 )-F (r j )) j∈ s+1,s+T to show that the exit condition is satisfied for j = s + T .

2. We use the first point to show that the number of iterations j i=0 n i is necessarily bounded by 2T n j . The conclusion of Theorem 1 comes from Lemma 4 which gives an upper bound of n j .

Proof of Theorem 1

Let C > 4 and ε > 0. We define

T = 1 + log 1+ 16 C 2 -16 2L(F (r 0 )-F * ) ε 2 log C 2 4 -1
where • is the ceiling function. We show that there is a doubling step at least every T iterations. Assume that there is no doubling step from j = s + 1 to j = s + T where s 1 which means that for a given s 1:

∀j ∈ s + 1, s + T , n j-1 > C L μj , (32)
and consequently:

∀j ∈ s + 1, s + T , n j = n s . (33)
Hence:

∀j ∈ s + 2, s + T , μj = min i∈N * i<j 4L (n i-1 + 1) 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) min i∈N * s<i<j 4L (n i-1 + 1) 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) min i∈N * s<i<j 4L n i-1 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) min i∈N * s<i<j 4L n s 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) 4L n s 2 min i∈N * s<i<j F (r i-1 ) -F (r j ) F (r i ) -F (r j ) , ( 34
)
as there is no doubling step for j s. This then implies

∀j ∈ s + 2, s + T , μj 4L n s 2 F (r j-2 ) -F (r j ) F (r j-1 ) -F (r j ) . (35)
Combining [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF] with [START_REF] Polyak | Introduction to optimization. optimization software[END_REF] and [START_REF] Roulet | Sharpness, restart, and acceleration[END_REF] we get that:

n s > C L 4L ns 2 F (r j-2 )-F (r j ) F (r j-1 )-F (r j ) = n s C 2
F (r j-1 ) -F (r j ) F (r j-2 ) -F (r j ) . ( 36
)
This leads to the following inequality

F (r j-2 ) -F (r j ) > C 2 4 (F (r j-1 ) -F (r j )), (37)
and then,

F (r j-2 ) -F (r j-1 ) > C 2 4 -1 (F (r j-1 ) -F (r j )). ( 38
)
Since C > 2 we get that

F (r j-1 ) -F (r j ) < 4 C 2 -4 (F (r j-2 ) -F (r j-1 )). (39)
We consider the case j = s + 1:

μs+1 = min i∈N * i<s+1 4L (n i-1 + 1) 2 F (r i-1 ) -F (r s+1 ) F (r i ) -F (r s+1 ) 4L (n s-1 + 1) 2 F (r s-1 ) -F (r s+1 ) F (r s ) -F (r s+1 ) 4L ( ns 2 + 1) 2 F (r s-1 ) -F (r s+1 ) F (r s ) -F (r s+1 ) 16L n s 2 F (r s-1 ) -F (r s+1 ) F (r s ) -F (r s+1 ) , (40)
as n s 2n s-1 .

By similar computations we get

F (r s ) -F (r s+1 ) < 16 C 2 -16 (F (r s-1 ) -F (r s )). ( 41
)
Since C > 4 we finally obtain the following inequalities

F (r s ) -F (r s+1 ) < 16 C 2 -16 (F (r s-1 ) -F (r s )). ( 42
) ∀j ∈ s + 2, s + T , F (r j-1 ) -F (r j ) < 4 C 2 -4 (F (r j-2 ) -F (r j-1 )). ( 43
)
We introduce Lemma 5 which links the composite gradient mapping g to the function F . This lemma is proven in Section 4.8:

Lemma 5. Let F satisfy the assumption H L for some L > 0. Then the sequence (r j ) j∈N provided by Algorithm 2 satisfies

∀j 1, 1 2L g(r j-1 ) 2 F (r j-1 ) -F (r j ). ( 44
)
From Lemma 5 and inequalities ( 42) and ( 43) we obtain the following sequence of inequalities

1 2L g(r s+T -1 ) 2 F (r s+T -1 ) -F (r s+T ) 4 C 2 -4 (F (r s+T -2 ) -F (r s+T -1 )) 4 C 2 -4 T -1 16 C 2 -16 (F (r s-1 ) -F (r s )) 4 C 2 -4 T -1 16 C 2 -16 (F (r 0 ) -F * ) 4 C 2 -4     log 1+ 16 C 2 -16 2L(F (r 0 )-F * ) ε 2 log ( C 2 4 -1 )     16 C 2 -16 (F (r 0 ) -F * ) 4 C 2 -4 log 1+ 16 C 2 -16 2L(F (r 0 )-F * ) ε 2 log ( C 2 4 -1 ) 16 C 2 -16 (F (r 0 ) -F * ) 1 1 + 16 C 2 -16 2L(F (r 0 )-F * ) ε 2 16 C 2 -16 (F (r 0 ) -F * ) ε 2 2L .
As a consequence, if there are T consecutive steps of Algorithm 2 without doubling the number of iterations, then the exit condition g(r j ) ε is eventually satisfied. This means that there is a doubling step at least every T steps and for all s 1 there exists j ∈ s + 1, s + T such that

n j-1 < C L μj . (45)
This implies that n j = 2n j-1 . As (n j ) j∈N is an increasing sequence, we get that n s+T n j = 2n j-1 2n s . And thus

n s n s+T 2 , ∀s 1. ( 46
)
Let us rewrite j as j = m + nT where 0 m < T and n 0. The increasing nature of (n j ) j∈N gives us that

j i=0 n i = m+nT i=0 n i = m i=0 n i + n-1 l=0 T i=1 n m+i+lT (47) T n m + T n l=1 n m+lT = T n l=0 n m+lT = T n l=0 n j-lT . (48)
According to equation ( 46) we have n j-T n j 2 and therefore

n j-lT 1 2 l n j , ∀l ∈ [[0, n]]. (49)
We obtain the following inequalities

j i=0 n i T n l=0 n j-lT T n l=0 1 2 l n j T ∞ l=0 1 2 l n j = 2T n j . (50)
From (50) and Lemma 4 we get that for j > 0

j i=0 n i 2T n j 4C L µ T 4C L µ   1 +     log 1 + 16 C 2 -16 2L(F (r 0 )-F * ) ε 2 log C 2 4 -1       (51) 4C log C 2 4 -1 L µ 2 log C 2 4 -1 + log 1 + 16 C 2 -16 2L(F (r 0 ) -F * ) ε 2 .
(52)

Proof of Corollary 1

Let F satisfy H L and G 2 µ for some L > 0 and µ > 0. Let (r j ) j∈N and (n j ) j∈N be the sequences provided by Algorithm 2 with C > 4 and ε > 0.

Minimizing the function

C → 4C log( C 2 4 -1)
gives us the optimal value Ĉ ≈ 6.38. This choice leads to the following rate:

F (r + j ) -F * = O e -1 12 √ µ L j i=0 n i . (62)
4.4 Proof of Proposition 1

(i) It is well known (see [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF][START_REF] Alamo | Restart FISTA with global linear convergence[END_REF]) that as F satisfies H L for some L > 0, the sequence (x k ) k∈N provided by FISTA (Algorithm 1) satisfies

∀k ∈ N * , F (x k ) -F * 2L (k + 1) 2 x 0 -x * 2 , (63)
where x * is any minimizer of F . This inequality is true for all x * ∈ X * so (63) can be rewritten

∀k ∈ N * , F (x k ) -F * 2L (k + 1) 2 d(x 0 , X * ) 2 . ( 64
)
Furthermore, F satisfies the growth condition G 2 µ so we can conclude by combining ( 6) and ( 64). (ii) We first prove the following claim. Let y ∈ R N . Then we have

∀x ∈ R N , F (y + ) + L 2 y + -x 2 F (x) + L 2 x -y 2 . ( 65
)
By definition of the proximal operator (5), y + is the unique minimizer of the function defined by

x → h(x) + ∇f (y), x -y + L 2 x -y 2 . ( 66
)
As this function is L-strongly convex we get that for all x ∈ R N ,

h(y + )+ y + -y, ∇f (y) + L 2 y + -y 2 + L 2 y + -x 2 h(x)+ x-y, ∇f (y) + L 2 x-y 2 .
(67) f has a L-Lipschitz gradient which implies that

∀(x, y) ∈ R N × R N , f (x) f (y) + x -y, ∇f (y) + L 2 x -y 2 . ( 68
)
Thus we get that

h(y + ) + f (y + ) -f (y) + L 2 y + -x 2 h(x) + x -y, ∇f (y) + L 2 x -y 2 . ( 69
)
The convexity of f gives us that f (x) f (y) + x -y, ∇f (y) and then

F (y + ) + L 2 y + -x 2 F (x) + L 2 x -y 2 . ( 70
)
By applying this inequality to y = y k = x k + k-1 k+2 (x k -x k-1 ), y + = x k+1 and x = x k for k 1 we have

F (x k+1 ) + L 2 x k+1 -x k 2 F (x k ) + L 2 k -1 k + 2 2 x k -x k-1 2
(71)

F (x k ) + L 2 x k -x k-1 2 . ( 72
)
Moreover, if we apply (65) to y = x = x 0 and y + = x 1 we get that

F (x 1 ) + L 2 x 1 -x 0 2 F (x 0 ). ( 73
)
This implies that

∀k ≥ 1, F (x k ) + x k -x k-1 2 F (x 0 ), (74)
As a consequence, we have that for all k > 0 such that (k + 1)

2 4L µ 1 + 1 γ , F (x k ) -F * γ (F (x 0 ) -F (x k )) . (76)
The contraposite of this proposition leads us to the expected conclusion.

Proof of Lemma 2

Suppose that F satisfies H L and G 2 µ for some L > 0 and µ > 0. Then Lemma 1 states that there exists c > 0 such that ∀x ∈ R N , F (x) -F * cd(0, ∂F (x)) 2 .

(77)

In particular, this assertion is true for c = 1 2µ .

Thus:

µ < 4L 1 + F (r i-1 )-F (r i ) F (r i )-F (r i+k )

(n i-1 + 1) 2 = 4L (n i-1 + 1) 2 F (r i-1 ) -F (r i+k ) F (r i ) -F (r i+k ) .

(89)

As a consequence, µ < μj . Furthermore, we have

The second claim of Proposition 1 implies that F (r j+1 ) F (r j ). As a consequence the function defined by y → F (r i-1 )-y F (r i )-y is an increasing homographic function and we get that ∀j 2, μj+1 min

i∈N * i<j 4L (n i-1 + 1) 2 F (r i-1 ) -F (r j+1 ) F (r i ) -F (r j+1 ) (92) min i∈N * i<j 4L (n i-1 + 1) 2 F (r i-1 ) -F (r j ) F (r i ) -F (r j ) (93) μj .
(94)

Proof of Lemma 4

The sequence (n j ) j∈N is defined such that for all j 2, n j = 2n j-1 if the following condition is satisfied:

n j-1 C L μj . (95)
The second claim of Proposition 1 implies that

n j-1 C L µ . (96)
This inequality ensures n j 2C L µ if j 2. For j = 0 and j = 1, n j 2C

2C L µ and we get the final conclusion.

Proof of Lemma 5

Let j 1. We can rewrite the inequality (65) for x = y ∈ R N :

L 2 y + -y 2 F (y) -F (y + ). ( 97
)
By setting y = r j-1 and as F (r + j-1 ) F (r j ) we conclude that 1 2L g(r j-1 ) 2 F (r j-1 ) -F (r j ).

(98)

Conclusions

We introduced an adaptive FISTA restart scheme for convex composite functions satisfying a quadratic growth condition around their minimizers. This method relies on an automatic estimation of the growth parameter which is generally not known in practice. The theoretical convergence rate provided is the fastest in the literature so far if the growth parameter cannot be estimated. Numerical experiments emphasize the efficiency of this method especially in the case when computing the function to minimize is expensive.

) 5 . 6 .

 56 Adaptive restart scheme introduced by Alamo et al. in [2], Gradient based restart scheme introduced by Alamo et al. in [1]. Note that methods 3, 4, 5 and 6 are respectively refered as Empirical restart scheme 1, Empirical restart scheme 2, Restart scheme 1 by Alamo et al. and Restart scheme 2 by Alamo et al. in the following figures.

Figure 1 :

 1 Figure 1: Comparison of several algorithms in terms of iterations for a L 1regularized least squares problem with N = 1000.

Figure 2 :

 2 Figure 2: Example of image inpainting: the damaged image y is on the left and an approximation of the solution of (31) is on the right.

Figure 3 :

 3 Figure 3: Comparison of several algorithms in terms of iterations for an inpainting problem.

Figure 4 :

 4 Figure 4: Comparison of several algorithms in terms of computation time for an inpainting problem.

  and thus we can conclude. (iii) By rewriting the first claim of Proposition 1 we get that∀k ∈ N * , F (x k ) -F * 4L µ(k + 1) 2 -4L (F (x 0 ) -F (x k )) .

2 F 2 F

 22 (r i-1 ) -F (r j+1 ) F (r i ) -F (r j+1 ) (r i-1 ) -F (r j+1 ) F (r i ) -F (r j+1 ).

Table 1 :

 1 Convergence rates of classical algorithms for a possibly non differentiable function F satisfying H L and G 2 µ for some L > 0 and µ > 0.

	Algorithm	References	Convergence	Limitations
			rate
	Forward-	Garrigos et al.	O e -µ L k	-
	Backward Heavy-Ball	[19] Aujol et al. [8]	O e -(2-√	2)	√ µ L k	Requires an
	variation						estimate of µ and
							uniqueness of
	Optimal FISTA	Necoara et al.	O e -1 e	√ µ L k	minimizer Requires an
	restart	[28]					estimate of µ
	Empirical FISTA	O'Donoghue and	O(k -2 )	Numerically fast
	restart	Candès [31]					but no improved
							theoretical
							convergence rate
	FISTA restart by	Fercoq and Qu
	Fercoq and Qu	[17]

Acknowledgements J-F Aujol acknowledges the support of the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement No777826. The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-PRC-CE23 Masdol and the support of FMJH Program PGMO 2019-0024 and from the support to this program from EDF-Thales-Orange.

We consider the case in which the exit condition g(r j ) ε is satisfied at first for at least 8C L µ iterations. We define the function ψ µ : R * + → 8C L µ , +∞ such that:

γ .

(53) According to Theorem 1, the number of iterations required to ensure that g(r j ) ε satisfies:

As ψ µ is a strictly decreasing function and j i=0 n i > 8C L µ we can write that:

where ψ -1 µ is the inverse function of ψ µ . By applying Lemma 2 we get that:

Elementary computations give us that:

and thus we get:

We can then conclude that

Let x ∈ R N . By definition of the proximal operator (5), x + is the unique minimizer of the function defined by

and thus x + satisfies 0 ∈ ∂h(x

As a consequence we get that

Moreover as f has a L-Lipschitz gradient we have

By combining these inequalities we conclude that

Proof of Lemma 3

The sequence (μ j ) j 2 is defined such that

On the other hand, for all i ∈ N * and k ∈ N * , we have:

where γ = F (r i )-F (r i+k ) F (r i-1 )-F (r i ) . Moreover, the third claim of Proposition 1 gives us that:


Cécile Proust-Lima

Jeremy M G Taylor

March

Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of post-treatment PSA: a joint modelling approach

Keywords: Joint latent class model, Mixed model, Predictive accuracy, Error of prediction, Posterior probability, Prostate cancer prognosis 2

PSA is a biomarker routinely and repeatedly measured on prostate cancer patients treated by radiation therapy. It was shown recently that its whole pattern over time rather than just its current level was strongly associated with prostate cancer recurrence. To more accurately guide clinical decision-making, monitoring of PSA after radiation therapy would be aided by dynamic powerful prognostic tools that incorporate the complete post-treatment PSA evolution.

In this work, we propose a dynamic prognostic tool derived from a joint latent class model, and provide a measure of variability obtained from the parameters asymptotic distribution. To validate this prognostic tool, we consider predictive accuracy measures and provide an empirical estimate of their variability. We also show how to use them in the longitudinal context to compare the dynamic prognostic tool we developed with a proportional hazard model including either baseline covariates or baseline covariates and the expected level of PSA at the time of prediction in a landmark model. Using data from three large cohorts of patients treated after the diagnosis of prostate cancer, we show that the dynamic prognostic tool based on the joint model reduces the error of prediction and offers a powerful tool for individual prediction.

Introduction

Prostate Specific Antigen (PSA) is a commonly used biomarker to monitor patients after treatment who received radiation therapy (RT) for localized prostate cancer. A rise of post-treatment PSA is highly predictive of clinical recurrence [START_REF] Sartor | Rate of psa rise predicts metastatic versus local recurrence after definitive radiotherapy[END_REF][START_REF] D'amico | Prostate specific antigen doubling time as a surrogate end point for prostate cancer specific mortality following radical prostatectomy or radiation therapy[END_REF] and definitions of biochemical recurrence have been suggested based on PSA crossing a threshold [START_REF] Roach | Defining biochemical failure following radiotherapy with or without hormonal therapy in men with clinically localized prostate cancer: recommendations of the RTOG-ASTRO Phoenix Consensus Conference[END_REF].

Recently, [START_REF] Thompson | Operating characteristics of prostate-specific antigen in men with an initial PSA level of 3.0 ng/ml or lower[END_REF] argued that the rise of PSA above a given threshold was not a satisfactory surrogate for detecting a clinical recurrence, that PSA was a continuous marker of disease progression and that its whole trajectory over time should be considered. In practice, detecting early signs of a recurrence is of major importance to assist in the patient's care and may facilitate the decision to initiate further treatment, such as salvage androgen deprivation therapy.

To more accurately guide clinical decision-making, monitoring of PSA after RT would be aided by dynamic powerful prognostic tools that incorporate the complete post-treatment PSA pattern. In this paper we refer to the pattern of PSA values for a subject as they evolve over time as the PSA trajectory. [START_REF] Tsiatis | Modeling the relationship of survival to longitudinal data measured with error. Applications to survival and CD4 counts in patients with AIDS[END_REF] stressed the importance of incorporating the complete biomarker information as a time-continuous process in order to avoid biases due to the periodically measured biomarker and measurement errors [START_REF] Prentice | Covariate measurement errors and parameter estimation in cox's failure time regression model[END_REF].

Joint modeling of repeated measures of PSA and time-to-recurrence provides such modeling in an efficient way by combining a mixed model for the change over time of the marker and a survival model that describes the associated risk of the event [START_REF] Henderson | Joint modelling of longitudinal measurements and event time data[END_REF]. In the prostate cancer context, using a shared randomeffect model, [START_REF] Pauler | Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes[END_REF] demonstrated that accounting for the trajectory of PSA improved the fit of the data compared to including only a summary measure of PSA dynamics in the Cox model. [START_REF] Yu | Joint longitudinal-survival-cure models and their application to prostate cancer[END_REF] showed that a joint model of longitudinal PSA measures and risk of recurrence could reduce the bias of the time-to-event parameters due to informative censoring, and that the posterior distribution of the probability of event could be used to monitor progression of the disease [START_REF] Taylor | Individualized predictions of disease progression following radiation therapy for prostate cancer[END_REF]. However, the numerical complexity of joint models has so far limited their application as a prognostic tool [START_REF] Pauler | Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes[END_REF]. The joint latent class model (JLCM), a different type of joint model, avoids many of the numerical complexities of the shared randomeffects model [START_REF] Lin | Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostate-specific antigen readings and prostate cancer[END_REF][START_REF] Proust-Lima | Joint modelling of multivariate longitudinal outcomes and a time-to-event: a nonlinear latent class approach[END_REF]. The JLCM assumes that the dependency between the risk of event and the trajectory of the biomarker is entirely captured by a latent class structure rather than by individual random effects. This class of models is particularly useful for heterogeneous populations, such as encountered in the study of recurrence of prostate cancer [START_REF] Sartor | Rate of psa rise predicts metastatic versus local recurrence after definitive radiotherapy[END_REF].

Using the JLCM as a computationally attractive example of a joint model, this paper builds a dynamic prognostic tool for early detection of prostate cancer recurrence and assesses its predictive ability on two large cohorts of patients treated by RT for prostate cancer. We specifically evaluate whether accounting for posttreatment PSA measures via a JLCM reduces the error of prediction compared to models with only pre-treatment prognostic factors. We also compare its predictive performance with those of a landmark (or conditional) model [START_REF] Van Houwelingen | Dynamic prediction by landmarking in event history analysis[END_REF][START_REF] Zheng | Partly conditional survival models for longitudinal data[END_REF][START_REF] Schoop | Quantifying the predictive performance of prognostic models for censored survival data with time-dependent covariates[END_REF] that can also incorporate post-treatment PSA measures.

In section 2, we describe the JLCM, the dynamic prognostic tool that is derived from this model and the computation of its standard error, as well as alternative predictive tools. Section 3 focuses on predictive accuracy measures used to compare predictive abilities of the tools. In section 4, we build the prognostic tool on a large cohort of patients, illustrate its use on individual patients, and show that the dynamic prognostic tool from the JLCM has better predictive accuracy compared to simpler prognostic tools on two independent large cohorts. Finally, we discuss the methodology and the results.

Dynamic prognostic tool from a joint model

Joint latent class model

Latent class structure

Following the model formulation of [START_REF] Lin | Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostate-specific antigen readings and prostate cancer[END_REF] and [START_REF] Proust-Lima | Joint modelling of multivariate longitudinal outcomes and a time-to-event: a nonlinear latent class approach[END_REF], we assume that the population of patients after RT can be divided into G latent classes. The latent class membership is defined by a categorical latent variable c i . The probability π ig that subject i (i = 1, ..., N ) belongs to latent class g (g = 1, ..., G) is related to covariates X pi in a multinomial logistic regression model:

π ig = P (c i = g|X pi ) = e ξ 0g +X T pi ξ 1g G l=1 e ξ 0l +X T pi ξ 1l (1)
where ξ 0g is the intercept for class g and ξ 1g is the vector of class-specific parameters associated with the vector of time-independent covariates X pi . For identifiability, ξ 01 = 0 and ξ 11 = 0. In the JLCM, the latent class structure is assumed to capture the entire dependency between the biomarker trajectory and the risk of the event so that, as shown also in the directed graph in Figure S1 of Supplementary Material (http://biostatistics.oxfordjournals.org), PSA trajectory and risk of recurrence are independent given the latent class membership.

Pattern of PSA changes

Let Y * i (t) be the value of PSA at time t for subject i, i = 1, ..., N , whose times of measurements are t ij , j = 1, ..., n i . The PSA trajectory is described by a linear mixed model [START_REF] Laird | Random-effects models for longitudinal data[END_REF] specific to class g on the logarithm scale:

Y i (t) | c i =g = ln(Y * i (t) + 0.1) = (u 0ig + β 0 X 0i ) + (u 1ig + β 1 X 1i )f 1 (t) + (u 2ig + β 2 X 2i )f 2 (t) + i (t) (2)
The parametric functions f 1 and f 2 were chosen in a preliminary analysis of 5 large cohorts of patients that described the progression of PSA after RT (Proust-Lima et al., 2008). The function f 1 represents the initial decline of PSA after the end of RT. Using a profile maximum likelihood technique for the transformation family defined by f 1 (t, η) = ((1+t) η -1), we found that η = -1.5 provided the best fit. The function f 2 represents the long-term rise in PSA after the end of radiation.

By considering the profile likelihood for the family of functions f 2 (t, ν) = t ν+1 /((t+ 1) ν ), we found that f 2 (t, 0) = t gave the best fit over the cohorts of patients.

We note that this corresponds to a long term exponential rise in PSA, which has been previously used [START_REF] Pauler | Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes[END_REF]Finkelstein, 2002, Yu et al., 2004). The vector of class-specific random-effects u ig = (u 0ig , u 1ig , u 2ig ) T follows a multivariate Gaussian distribution with mean vector µ g = (µ 0g , µ 1g , µ 2g ) T and unstructured variance-covariance matrix ω 2 g B, where ω 1 = 1 for identifiability. The vector µ g represents the mean trajectory of ln(PSA+0.1) over time in latent class g. The vectors X 0i , X 1i and X 2i are subvectors of X i , the total vector of covariates, and are associated with PSA trajectory through the regression parameters β 0 , β 1 and β 2 .

For the application, we chose those effects to be common over the classes. Finally, i (t) are independent Gaussian measurement errors with mean 0 and variance σ 2 .

Risk of recurrence

Let T * i be the time-to-recurrence and C i the censoring time. The observed event time is

T i = min(T * i , C i ). The indicator of recurrence E i equals 1 if T * i ≤ C i and 0 if C i < T * i .
We describe the risk of the event in latent class g by a proportional hazard model:

λ(t | c i = g, X ri ; ζ g , δ) = λ 0g (t; ζ g )e X ri (t)δ (3)
where X ri (t) is a vector of covariates that can be time-dependent. We assumed that the effect of covariates on the risk of recurrence was common over the latent classes, but class-specific effects could also be specified. Finally, λ 0g (t; ζ g ) is the baseline hazard in latent class g; in our application we used either a Weibull or a piecewise constant risk function. The latent class structure captures all the dependency between the biomarker evolution and the time-to-recurrence through λ 0g (t; ζ g ), so that neither the current value of the biomarker nor any other function of the random-effects appears in the survival model.

Maximum likelihood estimates

We denote by θ the vector of all the parameters. The log-likelihood of the observed data is

L(θ) = N i=1 ln G g=1 π ig f (y i | c i = g; θ)λ(T i | c i = g; θ) E i S(T i | c i = g; θ) (4)
where π ig = P (c i = g; θ) is defined in (1) and S(T i | c i = g; θ) is the survival function derived from (3). The density f (y i | c i = g; θ) of the vector of PSA measures y i in latent class g is the multivariate Normal density φ g (ỹ i ; θ) with mean and variance-covariance matrix described in [START_REF] Proust-Lima | Joint modelling of multivariate longitudinal outcomes and a time-to-event: a nonlinear latent class approach[END_REF]. For a given number of latent classes, maximum likelihood estimates are computed from (4) using a modified Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]. Convergence is assessed by stringent criteria based on the second derivatives of the log-likelihood and by using a grid of initial values. The optimal number of latent classes is determined using the Bayes Information Criterion (BIC) [START_REF] Schwartz | Estimating the dimension of a model[END_REF], as is typical in mixture models [START_REF] Hawkins | Determining the number of components in mixtures of linear models[END_REF]). An estimate of the variancecovariance matrix V (θ) of the parameters θ is given by the inverse of the Hessian matrix at the point estimate.

Posterior probability of recurrence

A posterior probability of recurrence can be easily derived from the JLCM and its parameters θ. This probability, which can be computed for a new subject using his available data at the current time, constitutes a dynamic prognostic tool of recurrence.

Dynamic prognostic tool derived from the joint model

Consider a new subject i free of recurrence at time s for whom the vector of repeated measures until s is denoted by Y (s) i

= {Y i (u), u ≤ s} and all the covariates X i included in (1), ( 2) and (3) are available. Let T * i denote the time of recurrence for subject i. Then, the posterior probability of recurrence between s and s + t for the parameter value θ can be easily computed:

P d i (s, t; θ) = P (T * i ≤ s + t | T * i ≥ s, Y (s) i , X i ; θ) = = G g=1 P (T * i ≤ s + t | T * i ≥ s, c i = g, X i ; θ)P (c i = g | T * i ≥ s, Y (s) i , X i ; θ) (5) 8 inserm-00367752,

version 1 -12 Mar 2009

The latent class structure entirely captures the dependence between the trajectory of the marker and the risk of recurrence. The predicted probability in ( 5) is the sum over the classes of the product of the class-specific conditional probability of the event that does not involve the marker measurements Y i (s)

P (T * i ≤ s + t | T * i ≥ s, c i = g, X i ; θ) = S(s|c i = g, X i , θ) -S(s + t|c i = g, X i , θ) S(s|c i = g, X i , θ) , (6)
and the posterior probability of class-membership given by

P (c i = g | T * i > s, Y (s) i , X i ; θ) = π ig f (y i | c i = g, X i ; θ)S(s|c i = g, X i , θ) G l=1 π il f (Y (s) i | c i = l, X i ; θ)S(s|c i = l, X i , θ
) .

(7)

Posterior probability of recurrence in simpler models

The standard proportional hazard model can be viewed as a specific case of the JLCM. If G = 1, the risk of recurrence is independent of the evolution of the biomarker, and the posterior probability of recurrence given in (5) reduces to:

P 0 i (s, t; θ) = P (T i ≤ s + t | T i ≥ s, X i ; θ) (8)
where X i are baseline covariates and θ the vector of parameters from the proportional hazard model given in (3) with G = 1.

When interest is in the prediction of an event after a certain time s, given the history of the event and covariates until that time, a model for the residual time distribution is required. A landmark (or partly conditional) model that does not specify the longitudinal process can be used [START_REF] Shi | Replacing time since HIV infection by marker values in predicting residual time to AIDS diagnosis[END_REF][START_REF] Zheng | Partly conditional survival models for longitudinal data[END_REF][START_REF] Van Houwelingen | Dynamic prediction by landmarking in event history analysis[END_REF][START_REF] Schoop | Quantifying the predictive performance of prognostic models for censored survival data with time-dependent covariates[END_REF]. This approach consists of fitting a survival model only to subjects still at risk at time s with covariates collected until s, specifically the repeated measures of the marker. It is common to reduce this general model and use only the value of the marker at time s. In our case, this reduces to a proportional hazard model fitted on subjects free of recurrence at time s with covariates X i and Y i (s). A different vector of parameters θ s is obtained for each time s and the predictive probability of recurrence is

P la i (s, t; θ s ) = P (T i ≤ s + t | T i ≥ s, X i ; Y i (s); θ s ). (9)
In practice, PSA is measured at discrete times, so that Y i (s) is not observed for all s. We considered two landmark models that differed in the imputation of Y i (s).

First we considered a "naive" landmark model that includes the last measure of PSA before s to approximate the value at time s. Although [START_REF] Prentice | Covariate measurement errors and parameter estimation in cox's failure time regression model[END_REF] showed that using the last measure of the marker to approximate the current marker value could induce bias, we used that approach because we wanted to provide a very simple model that included information about PSA. Second we instead extrapolated the value of the marker at the landmark point s using a two-stage approach [START_REF] Tsiatis | Modeling the relationship of survival to longitudinal data measured with error. Applications to survival and CD4 counts in patients with AIDS[END_REF].This method first estimates the mixed model for PSA evolution described in (2) with G = 1 on a training sample, and then uses the estimates to compute the empirical Bayes estimates of the random effects and to estimate the PSA level at time s for any new subject based on his PSA repeated measures before s.

We note that the landmark analysis does require a separate estimation of θ s for any time s needed for prognosis. A more complex approach proposed by Van Houwelingen (2007) defines θ s as a parametric function of s so that θ s can be 10 inserm-00367752, version 1 -12 Mar 2009 obtained at any s.

Estimate and measure of variability

For each of the prognostic tools given by ( 5), ( 8) or ( 9), the predictive probability of recurrence for patient i before time s + t given that he was free of recurrence before time s can be computed using the vector of parameters θ previously estimated on a sample, as P i (s, t) = P i (s, t; θ). Equivalently the predictive probability of being free of recurrence is denoted by S i (s, t) = 1-P i (s, t; θ). These are point estimates. To give a measure of their variability, we approximated the Bayesian posterior distribution of P i (s, t). Vectors (θ d ) d=1,...,D were drawn from the Normal approximation of the asymptotic distribution of θ, N θ, V ( θ) , so that the standard error could be estimated from the empirical standard deviation,

SE(s, t) = D -1 D d=1 (P i (s, t; θ d ) -P i (s, t)) 2
. This method does not involve any further estimation procedure. It requires only the point estimate θ and the variance V ( θ) computed once from the training sample. It avoids the need for a computationally intensive bootstrap resampling scheme and it is not model-specific in comparison to the ∆-method. However, it can only be used for parametric models. The 95% confidence bands of P i (s, t) can also be obtained as the 2.5 and 97.5 percentiles of P i (s, t; θ d ). The method can be extended to calculate standard errors of summary measures of predictive accuracy presented in Section 3 on testing samples. In that case, the calculation also involves bootstrapping the testing data.

Validation of predictive tools using predictive accuracy measures

When validating and comparing predictive rules in a survival context, there is no agreement in the literature about which measures should be preferred [START_REF] Pencina | Evaluating the added predictive ability of a new marker: From area under the ROC curve to reclassification and beyond[END_REF]. One can be interested either in the discriminative power of the predictive rule and use concordance measures derived from ROC methodology [START_REF] Heagerty | Survival model predictive accuracy and roc curves[END_REF]Zheng and Heagerty, 2007) or in the predictive accuracy of the rule [START_REF] Schemper | Predictive accuracy and explained variation in cox regression[END_REF]. In this work, we chose to focus on predictive accuracy measures that compare the actual value of predictions with the observed data. We used summary measures derived from the error of prediction (EP), err L;

X (t) = E[L(η(t) -Ŝ(t | X))],
where Ŝ(t | X) is the predictive rule regarded as fixed, η(t) the event status at time t, L is a loss function and the expectation is with respect to the joint distribution of T and X. Several estimators of EP were proposed either for time-independent rules [START_REF] Schemper | Predictive accuracy and explained variation in cox regression[END_REF]Henderson, 2000, Graf et al., 1999) or for dynamic rules [START_REF] Henderson | Identification and efficacy of longitudinal markers for survival[END_REF][START_REF] Schoop | Quantifying the predictive performance of prognostic models for censored survival data with time-dependent covariates[END_REF]. They differ in the loss function and the method used to account for censoring. In this work, we chose to focus on the estimate of absolute error of prediction proposed by [START_REF] Henderson | Identification and efficacy of longitudinal markers for survival[END_REF] that we found had the best properties in a simulation study.

For a dynamic rule Ŝ(s + t|X i (s)), the estimator of absolute error of prediction is computed at time s+t given information collected at time s and earlier:

ê rr X,s (s + t) = 1 N s Ns i=1 I(T i > s + t) 1 -Ŝ(s + t|T i > s, X i (s)) + E i I(T i ≤ s + t) 0 -Ŝ(s + t|T i > s, X i (s)) + (1 -E i )I(T i ≤ s + t) × [ 1 -Ŝ(s + t|T i > s, X i (s)) Ŝ(s + t|T i > s, X i (s)) Ŝ(T i |T i > s, X i (s)) + 0 -Ŝ(s + t|T i > s, X i (s)) 1 - Ŝ(s + t|T i > s, X i (s)) Ŝ(T i |T i > s, X i (s)) ] ( 10
)
where N s is the number of subjects still at risk at time s.

In the dynamic prognosis context, the error of prediction is a two-dimensional curve so that summary measures are useful. We used two summary measures over a [0, τ ] window for a given time s: the absolute error of prediction at horizon τ (EP (τ )) and the weighted average absolute error of prediction (WAEP) over [0, τ ],

as proposed by [START_REF] Henderson | Identification and efficacy of longitudinal markers for survival[END_REF]. This weighted average integrates the absolute error of prediction over [0, τ ] using weights that correct for the reduction in the number of observed events at longer times due to censoring. The estimator of WAEP is:

ÎW AEP X,s (s + t) = n (s) τ k=1 d (s) k ( Ĝ(s)/ Ĝ(t k )) ê rr X,s (t k ) n (s) τ k=1 d (s) k ( Ĝ(s)/ Ĝ(t k )) (11)
where

d (s)
k is the number of events at time t k among subjects still at risk at time s and Ĝ(t k ) and Ĝ(s) are the Kaplan-Meier estimates of the censoring distribution at times t k and s.

Whatever the predictive accuracy measure of interest (EP(τ ) or WAEP over [0, τ ]), a relative measure of predictive accuracy can be developed that is analogous 13 inserm-00367752, version 1 -12 Mar 2009 to R 2 in linear regression: the proportion of predictive accuracy explained by covariates. For example, the proportion of predictive accuracy at time t added when using covariates X 1 and X 2 rather than X 1 alone would be (err L;X 1 (t)err L;X 1 ,X 2 (t))/err L;X 1 (t).

Application to Prostate Cancer

Three cohort studies

We considered data from three large prospective cohorts of patients treated by external beam radiation therapy for localized prostate cancer. The three cohorts were from University of Michigan (UM) [START_REF] Taylor | Individualized predictions of disease progression following radiation therapy for prostate cancer[END_REF], William Beaumont Hospital (WBH) [START_REF] Kestin | Defining biochemical cure for prostate carcinoma patients treated with external beam radiation therapy[END_REF] and Radiation Therapy Oncology Group (RTOG9406) [START_REF] Roach | Penile bulb dose and impotence after three-dimensional conformal radiotherapy for prostate cancer on RTOG 9406: findings from a prospective, multiinstitutional, phase i/ii doseescalation study[END_REF]. Patients were included in the analysis if they had a clinical stage T1-4 and neither positive nodes or metastases, they had at least one year follow-up without clinical recurrence or salvage androgen deprivation therapy (SADT), and they had at least two repeated measures of PSA before the end of the follow-up. The endpoint of interest was the first clinical recurrence so that all the PSA measures collected after the end of RT and before this point were included unless a SADT was received, in which case measures after SADT were deleted. Clinical recurrence was defined as any of the following: distant metastases, nodal recurrence, any palpable or biopsy-detected local recurrence 3 years or later after radiation; any local recurrence within 3 years of RT if the most recent PSA was >2ng/mL; death from prostate cancer. This definition was to allow for the possibility of residual local disease up to 3 years after RT.

Three pre-treatment prognostic factors were considered: Gleason score category (2-6,7,8-10), T-stage category (1,2,3-4) and the pre-treatment level of PSA (iPSA) transformed to ln(iPSA+0.1). As the risk of recurrence should be markedly reduced after SADT, a time-dependent indicator of SADT, equal to 0 before the time of SADT and 1 after, was included in the hazard model. The three cohorts are described in Table 1.

[Table 1 about here.] A prognostic tool should be validated on external data, that is data not used for the creation of the prognostic tool. Following the validation strategy suggested by [START_REF] Altman | What do we mean by validating a prognostic model?[END_REF], we developed the prognostic tool on WBH, the largest cohort, and evaluated its predictive performances on two external samples, UM and RTOG9406. Although they are different cohorts, WBH and UM were comparable in terms of pre-treatment covariates (T-stage, Gleason and iPSA) and proportion of recurrences, whereas subjects in RTOG9406 were usually in a earlier stage of the disease with a larger proportion of subjects with T-stage 1 and Gleason below 7. The number of recurrences in RTOG9406 was also smaller (6.8% vs. 16.9% and 15.0% in WBH and UM). classes differed mainly by their long-term rise of PSA and the risk of recurrence, a higher long term increase of PSA being associated with a higher risk of recurrence.

Classes 1 and 2, which were relatively close in terms of trajectory and risk of recurrence, were different in terms of class-membership parameters.

[Figure 1 about here.]

Dynamic prediction of prostate cancer recurrence

To illustrate the use of the posterior probabilities of recurrence at time s + t given the information collected until time s, we show in Figure 2 the predicted probability of recurrence for two patients from the UM cohort. For each one, the predicted cumulative risk of recurrence was computed using the joint latent class model (denoted 5LCM), the proportional hazard model with baseline covariates and a 5-step risk function (denoted Baseline), and the two-stage landmark model with a 5-step risk function (denoted PSA(s)). Predicted probabilities using the naive landmark model with the latest PSA measure were very similar to the twostage landmark model. 95% confidence bands were computed as described in Section 2.3.3 using 2000 draws. Predictions were made up to three years ahead, which is a reasonable time horizon in this clinical setting.

[Figure 2

about here.]

For the subject on the left who recurred within the first 4 years, updated risk of recurrence was computed every 6 months from 1 to 3.5 years after end of RT. The PSA pattern for this patient is characteristic of an early recurrence with a drop of PSA the first year after the end of RT and a subsequent rise of PSA. However, levels of PSA are relatively low. The 5-LCM-based prediction that accounts for the shape of the PSA trajectory rather than the level of PSA captures the high risk of recurrence, while the PSA(s) level based prediction remains relatively low.

In Supplementary Material, Figure S2 and S3 show predictions for 3 other patients (Rec2, Rec3, Rec4) with different PSA profiles who recurred. For each of them the 5-LCM model detects higher risk of recurrence earlier after the end of RT.

For the subject on the right who did not experience any recurrence within 6 years after RT, the updated risk of recurrence was computed every year from 1 to 6 years. This subject has a PSA trajectory characteristic of a cured patient.

However, as he has relatively bad prognostic factors (T-stage=3, Gleason=8 and iPSA=62.4 ng/mL), the prediction based on baseline covariates predicts a high probability of recurrence while the two dynamic prognostic tools update the risk of recurrence according to the PSA trajectory so that, as soon as 3 years after RT, the probabilities of recurrence they provide become very low. In Supplementary Material, predictions for a second cured patient (Cens1 in Figure S4) also show how accounting for PSA repeated measures allows a better understanding of the cancer progression.

These individual predictions underline the usefulness of dynamic prognostic tools that can adapt to the PSA trajectory. In the next section, we used predictive accuracy measures to corroborate these suggestions at a population level on the UM and RTOG datasets.

Validation of the prognostic tool on UM and RTOG

We evaluated the predictive accuracy of the prognostic tool based on the JLCM four times a year from s=1 year to s=6 years after the end of RT for a time horizon of 3 years. For each s, we computed the absolute error of prediction curves and display three summary measures: the weighted averaged error of prediction (WAEP) over 3 years, and the error of prediction (EP) at 1 and 3 year horizons.

The 5-latent class model (denoted 5LCM) performances were compared to those of a proportional hazard model including baseline covariates and a 5-step risk function (denoted Baseline), a proportional hazard model with a 5-step risk function but without any covariate (denoted No covariate) and two landmark models with a 5-step risk function (denoted either PSA(s) or Last PSA depending on how the level of PSA at time s was computed). The summary measures for cohort UM and RTOG are displayed in Figure 3. The estimates and standard errors of WAEP in the 5LCM model for UM and RTOG cohorts were 0.0816 (SE=0.0090) and 0.0422 (SE=0.0068) respectively after 1-year follow-up and 0.0614 (SE=0.0095) and 0.0472 (SE=0.0074) respectively after 3-year follow-up. Table 2 gives the relative gain in WAEP for the two landmark models and the JLCM compared to the model including only baseline information.

[Figure 3 about here.]

[Table 2

about here.]

For UM, whatever the summary measure, inclusion of baseline covariates improved the predictive accuracy for prostate cancer recurrence only when using information from the first 3 years (s ≤ 3). After that point, the model without covariates gave similar predictive accuracy even though covariates effects were highly significant. Accounting for the PSA measures in addition to the baseline covariates using either a joint model or a landmark model reduced markedly the absolute error of prediction in the first 6 years with relative gain in WAEP varying from 3% to close to 20% (see Table 2). Moreover, the joint model improved the predictive accuracy more than the landmark models when using information from 1 year (11.7% vs. 3.3% of gain compared to baseline information) to 2 years (15.3% vs. 7.7% of gain) after RT while the landmark models gave a better predictive accuracy at s=6 (gain of >19.5% vs. 7.3% for JLCM). This means that in the first 2 years, the expected PSA at time of landmark was not sufficiently predictive while later, the level of PSA mainly drove the predictions for UM.

For the RTOG cohort, the baseline covariates improved the predictive accuracy during the whole follow-up. Furthermore, accounting for PSA repeated measures improved markedly the predictive accuracy by reducing the absolute error during the whole follow-up and especially in the first 2 years, as seen also in Table 2 with gain in WAEP of 18.4% or 24.9% at 1 and 2 years vs. <4.8% and <9.2% for the landmark models. In this cohort, which included earlier stages of the disease, the landmark models did not capture all the predictive value of the PSA trajectory, underlining the relevance of models like JLCM that include the whole trajectory.

For the two cohorts, the two landmark models gave similar predictiveness.

From a specific time of prediction, the curve of EP describes the change in EP over horizon of prediction (Figure 4). At 1 year after RT, including the post-treatment PSA measures in the JLCM substantially reduced the error of prediction for the UM cohort at any horizon compared to the proportional hazard model including baseline covariates (e.g. 19% improvement at 3 year horizon while including the expected PSA(s) value improved the predictive accuracy by only 5% at 3 year horizon). In contrast, when using information until 3 years after RT, the landmarking analysis EP approached or surpassed the joint model EP suggesting that after a time, the level of PSA may be sufficient for determining the risk of recurrence for UM. Conversely, for RTOG cohort, both from 1 year and 3 years after RT, the joint model reduced markedly the error of prediction at any horizon compared to the landmarking analysis.

[Figure 4 about here.]

Discussion

Although it is well known that PSA is highly predictive of prostate cancer recurrence, its use for monitoring progression of the disease is still rather limited and typically restricted to a binary summary of the PSA dynamics. Joint models offer an efficient framework to quantify the probability of recurrence utilizing the repeated measures of PSA. We have shown how a joint latent class model (JLCM) could be used to provide a dynamic prognostic tool of recurrence that can be updated for each new measurement of PSA. The methodology would be similar for a shared random effects model, except that the computations would be more burdensome [START_REF] Pauler | Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes[END_REF]Finkelstein, 2002, Yu et al., 2004). The JLCM relies on the conditional independence of the PSA repeated measures and the recurrence of prostate cancer given the latent classes. The practical advantages of this are that the log-likelihood has a closed form, and that the predictive tool can be computed analytically. The construction of the tool requires the estimation of the parameters on a single population only once. The prognostic tool can then be computed analytically for any new subject, using any information about PSA repeated measures and at any time. Moreover, to aid the user in evaluating the variability of the prediction, standard error and confidence bands are computed using an approximation of the Bayesian posterior distribution. This technique can be used for parametric rules whenever the ∆-method is not straightforward and a bootstrap is too computationally intensive. One limitation of the JLCM would be that the number of latent classes cannot be directly estimated and has to be selected according to a criterion, commonly the BIC. We note that using BIC to choose the number of classes has become standard practice in mixture modeling (e.g. [START_REF] Hawkins | Determining the number of components in mixtures of linear models[END_REF]). Furthermore we found that the predictive accuracy was not markedly impacted by the choice of the number of latent classes.

For example, the gain in predictive accuracy compared to the PHM with baseline covariates was roughly the same for 4, 5 and 6 latent classes (e.g. gain in WAEP for RTOG: 16.5%, 17.5%, 22.5%, 21.6% and 19.8% for G=2,3,4,5 and 6).

The validation of prognostic tools on different cohorts is of primary importance in the process of developing a prognostic tool, especially when using a complex statistical model. Indeed, a complex model can be fine-tuned for the dataset on which it is estimated but have a poor fit on new data. Following the [START_REF] Altman | What do we mean by validating a prognostic model?[END_REF] hierarchy of increasingly stringent validation strategies, we directly validated our model on different datasets, from other centers and other investigators. That gave us a good appreciation of whether a prognostic tool based on a relatively complex model may improve the predictive accuracy in practice.

We found that the landmark approaches gave similar predictive accuracy as a joint model for one dataset, but its predictive accuracy was lower for the other dataset, suggesting its potential lack of generalizability. For these two datasets, we found consistently that updating the risk of recurrence using the trajectory of PSA was an important refinement and that, at least in the first years, the level of PSA at time of prognosis could not capture the whole predictiveness of the PSA trajectory.

There are many choices for how to validate a model [START_REF] Pencina | Evaluating the added predictive ability of a new marker: From area under the ROC curve to reclassification and beyond[END_REF].

We chose predictive accuracy measures that focus on predictiveness rather than discrimination. Measures of predictive accuracy have been criticized because of their lack of interpretation. We showed through the application that our chosen measures give an easily interpretable assessment of the relative gain by quantifying the gain in predictiveness of a new model compared to a standard one.

To conclude, joint modeling of a marker trajectory and a clinical outcome is an attractive approach for developing powerful prognostic tools that can help clinical decision making in chronic diseases. Predictive accuracy measures offer an umbrella of criteria on which a prognostic tool can be validated and compared to other existing rules.

4. 2

 2 Estimation of the joint latent class model on WBH The three baseline covariates were included in both the survival model and the latent class membership model. For the mixed model, as recommended by Proust-Lima et al. (2008), Gleason was only included in the long-term rise part of the model while T-stage was included in both short-term and long-term parts, and iPSA was included in all three terms in equation (2). The hazard was defined by a class-specific Weibull function since the Akaike criterion was systematically better compared to a 5-step hazard function. The JLCM was fitted for different numbers of classes. The values of BIC as the number of classes varied from 2 to 6 were 13514.5, 13386.2, 13347.6, 13327.1 and 13354.5 respectively, and the associated numbers of parameters were 39, 51, 63, 75 and 87. The class-specific predicted mean trajectories and survival functions are displayed in Figure 1. The latent

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Predicted mean evolution (A) and survival function (B) in the 5 latent classes of the selected joint latent class model on WBH data (N=1268). Predictions given for a subject with T-stage=2, Gleason=7 and iPSA=10n/mL

Table 1

 1 Description of the three cohorts BM (N=1268) UM (N=503) and RTOG (N=615), categorical variable shown as number (frequency), and continuous variables as mean (standard deviation)

	Variable		BM	UM	RTOG
			(N=1268)	(N=503)	(N=615)
	Event		190 (15.0)	85 (16.9)	42 (6.8)
	T-stage	1	431 (34.0) 163 (32.4) 348 (56.6)
		2	792 (62.5) 290 (57.7) 253 (41.1)
		3,4	45 (3.5)	50 (9.9)	14 (2.3)
	Gleason	<7	902 (71.1) 276 (54.9) 421 (68.4)
		7	252 (19.9) 188 (37.4) 156 (25.4)
		>7	114 (9.0)	39 (7.7)	38 (6.2)
	Hormonal Therapy		170 (13.4)	44 (8.8)	47 (7.6)
	ln(iPSA+0.1) (ln(ng/mL))	2.16 (0.84) 2.23 (0.92) 2.00 (0.61)
	Age (years)		72.7 (6.5)	69.0 (7.1)	68.0 (7.0)
	Time-to-recurence (years)	5.03 (2.71) 3.82 (2.49) 4.61 (2.00)
	Time-to-last-contact (years) 5.91 (3.30) 6.21 (3.41) 5.92 (2.03)

Table 2

 2 Relative gain (in %) of WAEP over 3 years of forecast for the two landmark models (Last PSA and PSA(s)) and the 5-class JLCM compared to the proportional hazard model with baseline information. Gain in WAEP are computed from s=1 to s=6 years after RT, and are given for each cohort UM and RTOG

	Cohort Time Landmark models	JLCM
		s	Last PSA PSA(s)
	UM	1	3.2	3.3	11.7
		2	7.7	7.7	15.3
		3	17.6	19.5	13.7
		4	7.2	11.3	8.9
		5	5.9	9.5	8.4
		6	19.9	19.5	7.3
	RTOG 1	4.8	4.2	18.4
		2	9.2	9.0	24.9
		3	19.2	25.1	26.7
		4	17.4	20.0	24.7
		5	12.6	8.9	14.8
		6	3.1	1.5	12.8

inserm-00367752, version 1 -12 Mar 2009

This research was supported by grants CA110518 and CA21661 from the US National Cancer Institute and by a grant in aid for post-doctoral fellows from LEEM recherche France 1 inserm-00367752, version 1 -12 Mar 2009

 Zheng, Y. and Heagerty, P. J. (2007). Prospective accuracy for longitudinal markers. Biometrics 63, 332-41.

List of Figures


Bénédicte Pincemin

Serge Heiden

Franck Mazuet

on Statistical Analysis of Textual

Keywords: textometry, active corpus, corpus annotation, corpus projection, film grammar, Les Actualités françaises newsreels, TXM software

Active corpus provides the possibility to apply searching and statistical computing as if corpus were reduced to selected words, whereas full text still remains visible in context display. This is mainly implemented in paradigmatic processing, yet it may concern syntagmatic processing or text display too. Here we experiment active corpus in syntagmatic processing. A projection generates a new corpus, in which words are semantic tags that were automatically assigned in a first step to the original data. This new corpus makes it easy to explore tag sequences, with any generic textometric tool available, however sparse the original annotation may be. This methodological path was applied to film grammar analysis on 10,000 archival descriptions of news reports. 19 camera shot and angle types were ed through queries and tagged. This annotation became the lexicon of the projected corpus that was used to study shot sequences. The annotation and projection tools we have run are available as utilities in TXM open-source software and should usefully serve many research projects.

1. Definition and state of the art 1.1. What is active corpus?

We define active corpus as the subset of textual data that is taken into account in the textometric processing. It may be called background as well, as it provides the total amount of frequencies for contrastive analysis [START_REF] Pincemin | Lexicométrie sur corpus étiquetés[END_REF]. This concept is useful to manage two different realities: the full source text and its relevant content for various processing.

Active corpus plays a role at several moments of textometric analysis. (a) When the text is displayed, parts of the text that do not belong to active corpus may be shown anyway (because the whole text matters in order to achieve the reading and interpretation) but with a distinctive layout (such as a gray print), so as to know what is actually considered. (b) For syntagmatic processing dealing with the flow of the text and considering words as successive tokens (like phrase search or n-gram computing), active corpus is a means for passing over secondary items (insertions, expansions). (c) For paradigmatic processing listing words out of context and operating on types (that unify repeated occurrences), such as statistics on lexical tables, active corpus selects the focus (the items for which results are requested) and the overall set that stands for the frequency distribution (reference corpus and table margins).

Current implementations

Active corpus is commonly implemented for paradigmatic processing, through several ways to adjust the rows of the lexical tables submitted to statistical computing that can be:  qualitative filtering: stop word list, part-of-speech selection, or analytic categories as assigned by an editable dictionary, like analytic keys in IRaMuTeQ that specify if a word is active, supplementary or excluded [START_REF] Loubère | Documentation IRaMuTeQ 0.6 alpha 3[END_REF];  quantitative threshold: either the n most frequent words, or words with a frequency higher than a given threshold;  table advanced building: compiling several selection results, merging or discarding rows, like the List function in Hyperbase [START_REF] Brunet | Hyperbase, Manuel de référence[END_REF].

Active corpus is then defined at the word-type level (in a dictionary, in indexes, in table rows) and is used in statistics computing. A preselection at the word-token level can also be combined (restriction to a subcorpus). Most frequently are the selected lexical rows both the statistical reference and the vocabulary under study -the two sets are identical. However, supplementary elements make it possible for a word to be considered in the analysis without weighing in the statistical reference. Conversely, one could display in the results only a subset of the words involved in the statistical calculus.

Concerning the active corpus in text display, this is rarely available in text analysis tools. Such a feature is more suitable to TEI-based tools that encode and render complex text internal structures, like TEITOK [START_REF] Janssen | TEITOK: Text-Faithful Annotated Corpora[END_REF].

As for active corpus in syntagmatic processing, a current illustration may be provided by the possibility, in a concordance view, to consider only some word positions when sorting contexts: for instance, the 2L position meaning "the second word to the left of the target word" (see AntConc for an example of implementation [START_REF] Anthony | AntConc (Windows, MacOS, Linux)[END_REF]). Nevertheless, active corpus appears still underused in syntagmatic processing.

Active corpus in TXM software

TXM implements the usual means to refine a lexical table to fit any active corpus (paradigmatic processing): token preselection (subcorpus), qualitative filtering (CQP search engine for rows selection), quantitative thresholds, operations on table rows. Moreover, the lexical table margins can be set either to the current row selection, or to the entire corpus or subcorpus, which allows focusing results on a subset of words while keeping a larger corpus as reference for frequencies. The direct access to specificity score computation through its four parameters (PlotSpecif utility, Fig. 4 below) provides the most accurate settings for the active corpus and the pattern to be studied, for any particular case of specificity evaluation.

TXM also provides means to deal with active corpus when displaying the text: in the XML-TEI Zero import process (Heiden, 2020: §4.9.5), the textual planes setting labelled "Out of text to edit" lists XML elements whose textual data content will be displayed without being indexed (these words are neither searched nor counted). And a CSS stylesheet can be edited to distinguish the layout for these elements if necessary. An example is given with red section titles in the text displayed in Figure 1.

A recent utility (WordProperty2Word) carries out a projection of a word annotation that creates a corpus whose words are selected annotation values. In doing so, it can retain a subset of original tokens and delete others. We argue this can be a pragmatic way to implement active corpus for subsequent syntagmatic analyses, as illustrated in the following section.

JADT 2022: 16th International Conference on Statistical Analysis of Textual Data

A use case scenario

Research context and addressed question

The ANTRACT research project [START_REF] Carrive | Transdiciplinary Analysis of a Corpus of French Newsreels: The ANTRACT Project[END_REF] studies the audiovisual corpus of Les Actualités françaises, composed of about 10,000 news reports that were broadcast as weekly editions in French movie theaters between 1945 and 1969. One of the fields of investigation is an historical description of the production process and its impact on the newsreel's features. From this perspective, the question was raised to capture data about the film grammar of these newsreels. That is, what kind of camera shots and angles do the newsreels actually use, and whether recurring patterns of shots and angle types contribute to film standardization.

We considered that such information might be obtained from the newsreel descriptive forms.

In these text documents, the Sequences section is a shot-by-shot description of the film footage (see example in Figure 1 and Appendix). We must keep in mind that this is secondary data, as it relies on librarian mediation, which can vary from one librarian to another, and which cannot be totally exhaustive. Nevertheless, we thought these professional descriptions were worth mining and might give insights on film grammar.

Analysis roadmap within TXM

A TXM formatted AF-NOTICES corpus has been created from the descriptive forms [START_REF] Pincemin | Textometry on Audiovisual Corpora: Experiments with TXM software[END_REF]. Its size amounts to 2 million words, among which 38% hapax (more information in Appendix). It is a structured corpus, in which Sequences sections can be targeted. The idea is to set a typology of shot scales: Close-Up, Medium Shot, etc. For each shot type, CQL search engine queries are elaborated to match the various wording librarians happen to use to mention the shot. An automatic annotation utility applies the queries to the corpus and adds a normalized encoding for shots, available as a new analytic lexical property. Then, we project the corpus along this property, so as to derive a new corpus that retains only sequences of shot types. In this new corpus, shot pattern analysis is no longer hindered, neither by wording variations nor by lexical insertions that break up shot strings in the original text.

Annotation step

Several annotation modes are available in TXM platform [START_REF] Heiden | Annotation-based Digital Text Corpora Analysis within the TXM Platform[END_REF]. We opted here for a lexical annotation, assigning a property to a word, for two reasons. First, this annotation is indexed in the CQP model so that it benefits from all TXM usual functionalities. Secondly, a lexical annotation fits the targeted one-shot/one-unit model.

We also chose to systematize the description as every association between a term and a shot value was formally encoded through a CQL query. In addition to the consistency and clarity of the model, this procedure is compatible with corpus updates, and even with description updates (categories can be refined while getting more experienced with the data).

Our model lists 19 camera shot and angle types, plus one miscellaneous type ("divers plans"), and one off-topic type that is used to manage polysemy and mark "plan" or "vue" occurrences that do not denote a shot or an angle. For our needs, 60 queries were designed (see Appendix).

We benefited from the CQLList2WordProperties utility to automatically apply these 60 queries to the corpus and encode the new "plan" word property, so that annotating was systematic and fast.

JADT 2022: 16th International Conference on Statistical Analysis of Textual Data Figure 1 shows a summary of the resulting annotation (a developed example is provided in the Appendix). The shot that is the most frequently noted in descriptions is the Extreme Wide Shot (10PG) with 23,721 occurrences, which represents one quarter of all annotated shot mentions (93,553 annotations) in a 2-million-word corpus. For instance, one can examine how shot mentions are distributed along time (Figure 2) and note that Handheld Shots seems to have been replaced by Medium Shots. Another pattern involves the use of Zoom lenses during newsreel shootings. The appearance of this type of camera optics in the progression chart signals its adoption by Les Actualités françaises cameramen in 1964 and its increasing use among the company crews from then on.

Projection step

The projection (WordProperty2Word) of the annotated corpus along its shot-type annotation ("plan" word property) creates a new corpus whose words are the sequence of shot types (see Figure 1, right panel). Thus, shot mentions are normalized, and other words are ignored. This new corpus is 91,898-word long. It can be explored with all usual textometric tools. An illustration of the usefulness of this new representation is the study of the patterns of shots beginning a report (Fig. 3 & Appendix): we can search and count for 3-shot patterns, whether they begin the report (ii) or occur anywhere within a report (i); we can formulate a basic query matching the first shot in a report (iv), and measure specificities of shot types in this leading position (v). Specificity can be computed for patterns too (Figure 4).

We then were able to conduct a comprehensive analysis (see Appendix) of how scale evolves in the very first shots. 9,551 reports have got at least two shots mentioned in their text description. Among these, 5,582 reports (58 %) have a second shot type that is at a closer scale than the first one, which we call a funnel pattern. 2,065 other reports (22 %) use a unique and constant type of shot in their description. So all other patterns (either a second shot type wider than the first one, or a funnel pattern occurring later than in the very first shots) happen in the few remaining 1,904 reports (20 %). We interpret this prominence of the funnel shot pattern as a kind of standardization of film grammar. This approach has a dual goal. It offers the public a didactical presentation of the news based upon a simple editing where the gradual narrowing of successive shot frames enables a fluid visual exposition of any given subject. It also uses the standardization and reiteration of the ad hoc sequences to accelerate the shooting and editing process, thus providing an economical production method to the newsreel companies.

A drawback of this projection-based approach for active corpus is some loss of context, since original words are no longer present. However, one can work simultaneously with both corpora (the one before and the one after projection) and view both representations of the text (see example Fig. 1 with panels (ii) and (iv)). Textual structures (sections, paragraphs) and their metadata are preserved, and each shot-type word in the projected corpus has got a property that records the lexical form of the original tagged word in the shot mention.

Summary and prospects

Active corpus is a critical parameter in textometric analysis as it allows finetuning of textual data for searching domain and for statistical computing. On the one hand, active corpus is quite fully implemented in textometric software as regards paradigmatic processing (on lexical tables); on the other hand, it is much less advanced concerning text display and syntagmatic processing. For this last case, we devised and tested a new analysis scenario based on utilities that were recently released in TXM open-source software. First step is annotating through queries, which enriches the corpus with consistent and relevant analytic categories, even on evolving data with periodical releases. Second step is a projection that derives a new corpus focused on the chosen analytic categories. Both corpora are useful to the analysis: the original one provides the lexical and textual contexts when looking back at text content; and the projected one is efficient to address textometric questions about analytic patterns. Such a strategy is worth disseminating since it provides an effective solution for many similar research scenarios and needs. A more powerful implementation of active corpus for syntagmatic processing may still be achieved with a full integration of analytic representation and textual context into a unique corpus.

This research was supported by the ANTRACT ANR project (ANR-17-CE38-0010).

Figure 1 .

 1 Figure 1. The Les Actualités françaises corpus in TXM, with shot type annotations. From left to right: (i) annotated shot tags and their frequencies in the corpus; (ii) view of a documentary description with all annotated words highlighted; (iii) view of the shot tag assigned to each word of this text; (iv) same text in its projected version. Below: log of executed commands.

Figure 2 .

 2 Figure 2. Progression chart showing the diachronic evolution of 4 shot types. By decreasing frequency: Handheld Shot (red), Medium Shot (blue), Title & Credits (green), Zoom In & Out (yellow)

  JADT 2022: 16th International Conference on Statistical Analysis of Textual Data

Figure 3 .

 3 Figure 3. Study of the film grammar in Les Actualités françaises. From left to right: (i) frequency list of 3-shot patterns and (ii) same thing for initial 3-shot patterns only; (iii) frequency list of shots and (iv) same thing for first shot in a report only; (iv) specificity score of shot types for the first position in reports. Bottom right panel: log of executed commands.

Figure 4 .

 4 Figure 4. Specificity result +7.7 for the funnel pattern10GP 12PM 14GP (Extreme Wide Shot, Medium Shot, close-Up) as first shot sequence in a report (TXM PlotSpecif utility output).

Appendix

Full resolution figures and supplementary materials are available in the HAL-SHS archive: https://halshs.archives-ouvertes.fr/halshs-03667319