O Bouaziz

A Guilloux
email: agathe.guilloux@upmc.fr

A penalized algorithm for event-specific rate models for recurrent events

Keywords: Recurrent events process, total variation penalization, Aalen model, Cox model

We introduce a covariate-specific total variation penalty in two semiparametric models for the rate function of recurrent event process. The two models are a stratified Cox model, introduced in Prentice et al. (1981), and a stratified Aalen's additive model. We show the consistency and asymptotic normality of our penalized estimators. We demonstrate, through a simulation study, that our estimators outperform classical estimators for small to moderate sample sizes. Finally an application to the bladder tumour data of Byar (1980) is presented.

Introduction

Recurrent events arise in clinical or epidemiological studies when each subject experiences repeated events over time. Clinical examples include repetition of asthma attacks, epileptic seizures or tumour recurrences for individual patients. In this context, proportional hazards models have been largely studied in the literature to model the rate or mean functions of recurrent event data.

For instance, [START_REF] Andersen | Cox's regression model for counting processes: a large sample study[END_REF] introduce a conditional Cox model where the recurrent events process is assumed to be a Poisson process. Similar proportional hazards models and extensions are considered in [START_REF] Lawless | Some simple robust methods for the analysis of recurrent events[END_REF], [START_REF] Lin | Accelerated failure time models for counting processes[END_REF], [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF] and [START_REF] Cai | Marginal means/rates models for multiple type recurrent event data[END_REF].

To model rate functions in a recurrent events context, a different approach consists in fitting a Cox model for each event count. Along these lines, [START_REF] Prentice | On the regression analysis of multivariate failure time data[END_REF] introduce two stratified proportional hazards models with event-specific baseline hazards and regression coefficients. Gap times and conditional models are presented in their paper and a marginal event-specific model is studied in [START_REF] Wei | Regression analysis of multivariate incomplete failure time data by modeling marginal distributions[END_REF]. We refer to [START_REF] Kelly | Survival analysis for recurrent event data: an application to childhood infectious diseases[END_REF] for a complete review of existing Cox-based recurrent event models.

Additive models provide an useful alternative to proportional hazards models. For classical counting processes, the Aalen model was first introduced in [START_REF] Aalen | A model for nonparametric regression analysis of counting processes[END_REF] and is extensively studied in [START_REF] Mckeague | Asymptotic theory for weighted least squares estimators in Aalen's additive risk model[END_REF], [START_REF] Huffer | Weighted least squares estimation for aalen's additive risk model[END_REF], [START_REF] Lin | Semiparametric analysis of the additive risk model[END_REF]. It is considered in the context of recurrent events in [START_REF] Scheike | The additive nonparametric and semiparametric Aalen model as the rate function for a counting process[END_REF]. We propose in this paper to consider an event-stratified version of the Aalen model, in the manner of [START_REF] Prentice | On the regression analysis of multivariate failure time data[END_REF].

As demonstrated in what follows, event-stratified models allow more flexibility but suffer from over-parametrization when the sample size is not large enough compared to the numbers of covariates and recurrent events. To address this drawback, we construct estimators which do not vary much between two consecutive recurrent events, or equivalently with small total variations (see Section 2.4 for details). To achieve this goal, we consider minimizers of empirical risks penalized via a covariate-specific total variation penalty. As such, our algorithms are part of the class of fused lassos algorithms. The latter have been introduced and studied, for noised piecewise constant signals, by [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF], [START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF] or [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF]. Lasso estimators in the context of survival analysis with high dimensional covariates have been introduced and studied in [START_REF] Tibshirani | The lasso method for variable selection in the cox model[END_REF], [START_REF] Huang | Oracle inequalities for the lasso in the cox model[END_REF] in the Cox model and Martinussen & Scheike (2009b), [START_REF] Gaiffas | High-dimensional additive hazards models and the lasso[END_REF] in the Aalen model, among others.

The settings of recurrent events with a terminal event, and the two models studied in this paper are presented in the next section. In paragraph 2.4, we describe our novel algorithm, which involves a total variation penalization of criteria, specific to either the multiplicative or additive models. It requires preliminary details on inference in these two models, which are given in paragraph 2.2 and 2.3. Consistency and asymptotics distributions of our estimators are derived in Section 3. Simulation studies and a real data analysis are provided in Sections 4 and 5. We conclude with a discussion in Section 6.

Models and algorithm

Models

Let D denote the time of the terminal event and Ñ (t) the càdlàg process that counts the number of recurrent events occurring in the interval (0, t], with the convention Ñ (0) = 0. As no recurrent events can occur after D, the process Ñ has jumps only on (0, D]. The p-dimensional process of covariates is denoted by X and is assumed to be left continuous. The event-specific rate function of the process Ñ , denoted by ρ 0 , is defined as

E d Ñ (t) | X(t), D t, Ñ (t-) = s -1 1( Ñ (t-) = s -1) = 1(D t, Ñ (t-) = s -1)ρ 0 (t, s, X(t))dt,
(2.1) for s = 1, 2, . . . and t in A s = {t : P[ Ñ (t-) = s -1, D t] > 0} and is null outside of A s . The set A s represents the time intervals where the s-th recurrent events can occur with a positive probability and lies in Supp(D), the support of D. Apart from the stratification, this definition of the rate function can be found in [START_REF] Scheike | The additive nonparametric and semiparametric Aalen model as the rate function for a counting process[END_REF].

We consider two semiparametric models for the function ρ 0 . The first one is an event-specific multiplicative rate model introduced in [START_REF] Prentice | On the regression analysis of multivariate failure time data[END_REF]. In this model, the rate function is specified, for s = 1, 2, . . . and t in A s , by

ρ 0 (t, s, X(t)) = α 0 (t, s) exp (X(t)β 0 (s)) (2.2)
where for each event number s, β 0 (s) is an unknown p-dimensional vector of parameters and α 0 is an unknown baseline function.

Following [START_REF] Scheike | The additive nonparametric and semiparametric Aalen model as the rate function for a counting process[END_REF], and [START_REF] Zeng | A semiparametric additive rate model for recurrent events with an informative terminal event[END_REF], we also propose to consider its additive counterpart. The rate function in our event-specific additive model is then for s = 1, 2, . . . and t in A s :

ρ 0 (t, s, X(t)) = α 0 (t, s) + X(t)β 0 (s). (2.3)
When β 0 is assumed to be constant over the events, models (2.2) and (2.3) are usually referred to as stratified Cox and Aalen models (see e.g. [START_REF] Martinussen | Dynamic regression models for survival data[END_REF], page 190). To insist on the constancy of β 0 (as a function of s), these particular cases of models are hereafter designated as "constant coefficient models".

As, in practice, the individuals experience only a finite number of recurrent events, we will concentrate on the estimation of the rate function for the first B events, where B is an user-chosen integer (see the example in Section 5). Mathematically, this means that we only consider the observation of the process Ñ on the interval [0, E(B)], where E(B) is the hitting time of [B, ∞).

Equivalently, we consider that we observe the stopped process N * , defined through N * (t) = Ñ (t ∧ E(B)), for all t 0. Noticing that, for all s = 1, . . . , B and all t 0, {N * (t-) = s -1} = { Ñ (t-) = s -1} ⊂ {t E(B)}, the event-specific rate function of N * equals the one of Ñ , such that equation (2.1) holds with Ñ replaced by N * . This is equivalent to assuming that the total number of observed events is almost surely bounded by B, which is the classical framework for inference for recurrent event processes (see e.g. [START_REF] Scheike | The additive nonparametric and semiparametric Aalen model as the rate function for a counting process[END_REF], [START_REF] Dauxois | Non-parametric tests for recurrent events under competing risks[END_REF] or [START_REF] Bouaziz | Nonparametric estimation of the intensity function of a recurrent event process[END_REF]).

We consider the problem of estimating the unknown parameter β 0 = (β 0 (1), . . . , β 0 (B)), in stratified models (2.2) and (2.3) on the basis of data from n independent individuals. Introducing the censoring time C, the data consist of n independent replications

{N i (t), T i , δ i , X i (t), t τ }, i = 1, . . . , n, where N i (t) = N * i (t ∧ C i ), T i = D i ∧ C i is the minimum between D i and C i , δ i = 1(D i C i ), (X i (t), 0 t T i )
is the covariates process and τ represents the end-point of the study. In addition, we define the event-specific at-risk function Y s and the overall at-risk function Y . For each individual i, for all t in [0, τ ]:

Y s i (t) = 1(T i t, N i (t-) = s -1), Y i (t) = B s=1 Y s i (t) = 1(T i t). Let A τ s = A s ∩ [0, τ ].
The following two assumptions are mandatory to perform estimation.

Assumption 1 For all s = 1, . . . , B, and

t in A τ s , E[Y s (t)] > 0 and P[E(B) τ ] > 0.
The fist part is classical in survival analysis (see for instance [START_REF] Andersen | Statistical models based on counting processes[END_REF]). The second part implies that, for s = 1, 2, . . . , B, the sets A τ s are non-empty. Note also that the processes Y s i are almost surely null on the complementary of A s in [0, τ ].

Assumption 2 For all s = 1, . . . , B, and

t in A τ s , E dN * (t) | X(t), D ∧ C t, N * (t-) = s -1 = E dN * (t) | X(t), D t, N * (t-) = s -1 .
This assumption is classical in recurrent events context and can be found for instance in [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF]. It is the analog of the independent right censoring definition III.2.1. of [START_REF] Andersen | Statistical models based on counting processes[END_REF]. We refer the reader to the Supplementary Material for a discussion on these assumptions.

In particular, sufficient conditions are presented for these assumptions to hold.

In our framework, the unknown vector of parameters β 0 has p × B unknown coefficients to be estimated. For reasonable sizes of sample n, these models are over-parametrized in the sense that, when √ n p × B, the estimators show very poor behaviour (see Section 4 for an illustration).

On the other hand, simpler forms of models (2.2) and (2.3), in which the unknown parameter does not change with the event, β 0 (s) = β 0 , might be too poor to accurately fit the data (see also Section 4 and the discussion in [START_REF] Kelly | Survival analysis for recurrent event data: an application to childhood infectious diseases[END_REF]). In this paper, we aim at providing estimators realizing a compromise between these two situations.

Inference in the multiplicative model

As in [START_REF] Prentice | On the regression analysis of multivariate failure time data[END_REF], in the multiplicative event-specific model (2.2), an estimator βES/mult of the unknown parameter β 0 ∈ R p×B is defined as the maximizer of the stratified partial loglikelihood, or equivalently as βES/mult ∈ argmin

β∈R p×B L P L n (β) (2.4) = argmin β∈R p×B   - 1 n B s=1 n i=1 [0,τ ]    X i (t)β(s) -log   n j=1 Y s j (t) exp (X j (t)β(s))      Y s i (t)dN i (t)   .
An estimator βC/mult in the constant coefficient model is defined as βC/mult ∈ argmin

β∈R p   - 1 n n i=1 B s=1 [0,τ ]    X i (t)β -log   n j=1 Y s j (t) exp (X j (t)β)      Y s i (t)dN i (t)   .
(2.5)

The right term of equation (2.5) is the standard log-likelihood of a stratified Cox model (see for instance [START_REF] Therneau | Modeling survival data: extending the Cox model[END_REF] p.44-45 or Kalbfleisch & Prentice (2011) p.118-119).

Inference in the additive model

As noticed in Martinussen & Scheike (2009a,b) or [START_REF] Gaiffas | High-dimensional additive hazards models and the lasso[END_REF], in the usual additive hazards model, the estimator βES/add of the unknown parameter β 0 ∈ R p×B can be written as the minimizer of a (partial) least-squares criterion:

βES/add ∈ argmin

β∈R p×B L P LS n (β) = argmin β∈R p×B B s=1 β(s) ⊤ H n (s)β(s) -2h n (s)β(s) , (2.6)
where for all s ∈ {1, . . . , B}, H n (s) are p×p symmetrical positive semidefinite matrices and h n (s) are p-dimensional vectors equal to

H n (s) = 1 n n i=1 [0,τ ] Y s i (t) X i (t)-Xs (t) ⊗2 dt and h n (s) = 1 n n i=1 [0,τ ] Y s i (t) X i (t)-Xs (t) dN i (t) with Xs (t) = n i=1 X i (t)Y s i (t)/ n i=1 Y s i (t)
and the convention that 0/0 = 0.

On the other hand, an estimator βC/add in the constant coefficient model is defined as βC/add ∈ argmin

β∈R p β ⊤ H n β -2h n β , with H n = B s=1 H n (s) and h n = B s=1 h n (s). (2.7)
This formula gives an analogue of the so-called stratified Cox model to the Aalen case. Note also that the quantities H n and h n are not the same as the ones involved in a standard Aalen model with no stratification (see for instance the terms D n and d n in Martinussen & Scheike (2009b)).

In our constant coefficient models β 0 is constant and therefore its estimators are also constants, but the baseline α 0 is still stratified with respect to s.

A total-variation penalty

To overcome the possible over-parametrization of models (2.2) and (2.3), we propose to define penalized versions of criteria (2.4) and (2.6). For all β = (β(s), s = 1, . . . , B) with β(s) = (β 1 (s), . . . , β p (s)), define for all j = 1, . . . , p

β j = (β j (1), . . . , β j (B)) and tv(β j ) = B s=2 |β j (s) -β j (s -1)| = B s=2 |∆β j (s)|.
(2.8)

We now consider the minimizers of the partial log-likelihood (respectively the partial leastsquares) penalized with a covariate specific total variation. Define the penalized estimators in models (2.2) and (2.3) as:

βtv/mult ∈ argmin

β∈R p×B    L P L n (β) + λ n n p j=1 tv(β j )    and
(2.9) βtv/add ∈ argmin

β∈R p×B    L P LS n (β) + λ n n p j=1 tv(β j )    .
(2.10)

These penalized algorithms are part of the class of fused-lasso algorithms (see e.g. [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF], [START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF] for a definition) and can be rewritten as lasso algorithms (the details are given in Supplementary Material).

Asymptotic results

We successively provide the asymptotic results for the estimators βtv/add in the additive model and βtv/mult in the multiplicative model. Their proofs are postponed in Supplementary Material.

In both models, the following conditions are mandatory.

Assumption 3 1. The covariate process X(•) is of bounded variation on [0, τ ].

2. There exists a constant M such that for all t in [0, τ ], X(t) ∈ [-M, M ] p almost surely.

For all

s = 1, . . . , B, A τ s α 0 (t, s)dt < ∞.
Parts 1 and 2 of Assumption 3 are equivalent to Assumption (ii) in [START_REF] Scheike | The additive nonparametric and semiparametric Aalen model as the rate function for a counting process[END_REF]. Part 3 is a stratified form of Assumption 7.2.1, for instance, in [START_REF] Andersen | Statistical models based on counting processes[END_REF].

Define, for all s = 1, . . . , B and t in A τ s , the centered process

M s (t) = N (t) - t 0 E dN (r) | X(r), D ∧ C r, N (r-) = s -1
and the p × p matrix

H(s) := A τ s E[Y s (t)X(t) ⊤ X(t)]dt - A τ s (E[Y s (t)X(t)]) ⊗2 E[Y s (t)] dt,
which, thanks to Assumptions 1 and 2, is well defined.

Theorem 3.1 Assume that, for each s = 1, . . . , B, H(s) is non-singular and that Assumptions 1, 2 and 3 are fulfilled.

1. If λ n /n → 0 as n → ∞ then βtv/add converges to β 0 in probability.

If λ

n / √ n → λ 0 0 as n → ∞ then √ n( βtv/add -β 0 ) converges in distribution to argmin u∈R p×B Λ add (u) = argmin u∈R p×B B s=1 u(s) ⊤ H(s)u(s) -2u(s) ⊤ ξ add (s) + λ 0 p j=1 B s=2
|∆u j (s)|1(∆β j 0 (s) = 0) + sgn(∆β j 0 (s))(∆u j (s))1(∆β j 0 (s) = 0) , and for each s, ξ add (s) is a centered p-dimensional gaussian vector with covariance matrix equal to

E A τ s (X(t) -E[Y s (t)X(t)]/E[Y s (t)])Y s (t)dM s (t)

⊗2

.

We now state an analogous result in the multiplicative model. Define for all s = 1, . . . , B and for all t in A τ s ,

s (l) (s, t, β) = E[Y s (t)X(t) ⊗l exp(X(t)β(s))], l = 0, 1, 2. Introduce e(s, t, β) = s (1) (s, t, β)/s (0) (s, t, β), v(s, t, β) = s (2) (s, t, β)/s (0) (s, t, β) -e(s, t, β) ⊗2 and Σ(s, β) = A τ s v(s, t, β)E[Y s (t)dN (t)].
For any s = 1, . . . , B and for any t in A τ s , the three functions s (l) (s, t, β 0 ) are bounded due to Assumption 3 and e(s, t, β), v(s, t, β) and Σ(s, β) are finite due to Assumptions 1 and 3.

Theorem 3.2 Assume that for each s = 1, . . . , B, Σ(s, β 0 ) is non-singular and that Assumptions 1, 2 and 3 are fulfilled.

1. If λ n /n → 0 as n → ∞ then βtv/mult converges to β 0 in probability.

If λ

n / √ n → λ 0 0 as n → ∞ then √ n( βtv/mult -β 0 ) converges in distribution to argmin u∈R p×B Λ mult (u) = argmin u∈R p×B B s=1 1 2 u(s) ⊤ Σ(s, β 0 )u(s) -u(s) ⊤ ξ mult (s) + λ 0 p j=1 B s=2
|∆u j (s)|1(∆β j 0 (s) = 0) + sgn(∆β j 0 (s))(∆u j (s))1(∆β j 0 (s) = 0) , and for each s, ξ mult (s) is a centered p-dimensional gaussian vector with covariance matrix equal to

E   A τ s (X(t) -e(s, t, β 0 )) Y s (t)dM s (t) ⊗2   .
First results of Theorems 3.1 and 3.2 prove the consistency of our estimators, on the mandatory condition that λ n /n tends to 0. This ensures that they behave better than the constant coefficient estimators when β 0 is non constant. In addition, the considered penalty will induce sparsity for each covariate j = 1, . . . , p in the successive differences ∆β j (s), s = 1, . . . , B. As a consequence, the effects of a covariate on two consecutive events will often be equal. We show, in the following simulation study, that this induced sparsity ameliorates the behaviour of our estimators compared to the unconstrained ones (defined in Equations (2.4) and (2.6)).

The second results show that asymptotic normality can be achieved only if λ n / √ n tends to 0. In that case, the asymptotic variance of the limiting distribution can be estimated by means of the analog of the optional variation in this context, see Martinussen & Scheike (2006, page 150-151) for details.

However, when λ 0 = 0, the algorithm is no longer consistent in selection, in the sense that in both multiplicative and additive cases, as n tends to infinity,

P {(j, s) ∈ R p × R B , ∆( βj T V )(s) = 0} = {(j, s) ∈ R p × R B , ∆(β j 0 )(s) = 0} → 0.
Even in the case where λ n / √ n → λ 0 > 0, this probability is asymptotically stricly less than 1, see [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] for details. To enhance the sparsity in the covariate-specific successive differences (or equivalently to force βT V /mult and βT V /add to have several constant coefficients), we consider a reweighted lasso, in the manner of [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] or [START_REF] Candès | Enhancing sparsity by reweighted l 1 minimization[END_REF]. In both models the two steps (or reweighted lasso) estimators are defined as βT V /mult ∈ argmin

β∈R p×B    L P L (β) + λ n n p j=1 B s=2 1 |∆ βj T V /mult (s)| + |∆ β0 | |∆β j (s)|    (3.11)
and βT V /add =∈ argmin

β∈R p×B    L P LS (β) + λ n n p j=1 B s=2 1 |∆ βj T V /add (s)| + |∆ β0 | |∆β j (s)|    , (3.12)
where βT V /mult and βT V /add are defined in Equations (2.9) and (2.10) and in both cases:

|∆ βj 0 | = min{|∆ βj T V (s)|, j = 1, . . . , p , s = 1, . . . , B , |∆ βj T V (s)| = 0}.

Implementation and simulation studies

We compare the performance of the penalized estimators (2.9) and (2.10), the constant coefficient ones (2.5) and (2.7), and the unconstrained ones (2.4) and (2.6). To mimic the bladder tumour cancer dataset studied in Section 5, we set p = 4 and consider B = 5 recurrent events for the estimation. In the multiplicative and additive models, the sample size n varies from n = 50 = 2.5 pB to n = 1000 ≃ (pB) 2.3 .

Implementation of penalized estimators.

The minimizers of the partial log-likelihood, respectively, the partial least-squares, penalized with a covariate specific total variation of Equations (2.9) and (2.10) can be seen as lasso estimators by introducing a block matrix D of size (pB × pB) with p diagonal blocks being equal to a lower triangular matrix with nonzero elements equal to 1 and p 2 -p off-diagonal blocks being matrices of zeros.

The minimization problems of Equations (2.9) and (2.10) can then be rewritten as

βT V = Dγ T V with γT V ∈ argmin γ∈R p×B    L n (Dγ) + λ n n p j=1 B s=2 |γ j (s)|    = argmin γ∈R p×B L n (Dγ) + λ n n γ 1 , (4.13)
where L n is either L P L n or L P LS n and, in both cases, γT V = ( β1

T V (1), ∆ β1 T V (2), . . . , ∆ βp T V (B)) ⊤ .
The related R functions can be found at http://www.lsta.upmc.fr/guilloux.php?main=publications.

The regularization parameter is chosen in both multiplicative and additive cases via 5-fold cross-validation, defined by the lasso formulations of Equation (4.13). For details in the multiplicative model, see Simon et al. (2011, pages 9-10) or [START_REF] Van Houwelingen | Cross-validated cox regression on microarray gene expression data[END_REF]. Details for the additive model may be found in Martinussen & Scheike (2009b, page 608).

Simulation scheme

We draw p = 4 constant covariates from uniform distributions on [0, 2] and set the parameters values at β 1 0 = 0.25(0, 0, 1, 1, 0), β 2 0 = (1, . . . , 1), β 3 0 = b(1, 2, 3, 4, 5) and β 4 0 = (0, . . . , 0). We generated recurrent event times from the multiplicative (2.2) and additive (2.3) models with baseline defined through the Weibull distribution with shape parameter a W and scale parameter 1 (see Supplementary Materials for a more detailed description of the simulation scheme). The death and censoring times are generated from exponential distributions with parameters a D and a C respectively. We set the value of parameter a W at 1.5, and of b at 4 in the additive case and -1 in the multiplicative case. Finally, the values of a D and a C are empirically determined to obtain 14 -15% of individuals experiencing the fifth event. More results for p obs = 28 -29% were obtained and are reported in Supplementary Material.

Performance evaluation

To evaluate the performance of the different estimators, we conduct a Monte Carlo study with M = 500 experiences. The estimation accuracy is investigated for each method via a mean squared rescaled error defined as

mse = 10 3 M M m=1 βm -β 0 2 β 0 2 , (4.14)
where βm is the estimation in the sample m. We furthermore study the detection power of non-constant (respectively constant) covariate effects by computing specificities (spec) and sen-sitivities (sens) for each method. For an estimation βm , define false (fp) and true (tp) positives

such that 1 is the ideal value for both indicators of detection power. The results are presented in Tables 1-2.

Results

As expected, the constant model is biased and, in particular, for our choice of a non-constant β 0 , the MSE does not decrease with the sample size n. The comparison between the unconstrained and penalized estimators is in favour of our estimators in all cases as long as n is much smaller than (pB) 2 . For p = 4, B = 5, n = 100 and p obs = 14% (which are values close to those encountered in the bladder tumour cancer dataset studied in the next section) our penalized estimators are respectively, 2.41 in the additive model and 1.38 in the multiplicative model, times better than the unconstrained ones in terms of estimation error.

Regarding the sensitivity and specificity indexes defined in Section 4.3, the unconstrained estimators, which, by definition, has no constant coefficients, has a perfect sensitivity and a null specificity. On the opposite, the constant estimator can not detect a non-constant effect and, as a consequence, produces no false positive, nor true positive, with the consequence that its sensitivity is null and its specificity equals one.

The one-step penalized estimators of Equations (2.9) and (2.10) have almost perfect sensitivities, but their specificities decrease with the sample size. In our opinion, this is due to the cross validation, which tends to choose smaller regularization parameters for larger sample sizes. This choice is consistent with the theoretical values of λ n displayed in Theorems 3.1 and 3.2. This however leads to the classical non-consistency in selection of one-step estimators, as discussed at the end of Section 3.

To overcome this drawback, we defined in Equations (3.11) and (3.12), two-steps estimators.

They are expected to enhance the sparsity in the estimated successive differences and consequently, produce less false positives and more true negatives. Results indeed show increased specificities for comparable sensitivities, as compared to the one-step estimators. These two-steps estimators however are expected to be more biased, this phenomenon can be seen in Table 2 for example. When the specificities of the one-step and two-steps estimators are of the same order, the MSE of the latter is greater.

We repeat the simulation study for a W = 0.5 and then for a Gompertz baseline with shape parameter a G = 1.5 and 0.5 and scale parameter 1. The results are reported in Supplementary Material. As expected, we observe that a drop in the percentage of individuals experiencing the fifth event drops affects the performances of all estimators. Other conclusions are similar.

Bladder tumour data analysis

In this section we illustrate the behaviour of our estimators on the bladder tumour cancer data of [START_REF] Byar | The veterans administration study of chemoprophylaxis for recurrent stage 1 bladder tumors: comparison of placebo, pyridoxine, and topical thiotepa. Bladder Tumors and Others Topics in Urological Oncology[END_REF]. These data were obtained from a clinical trial conducted by the Veterans Administration Co-operative Urological Group. One hundred and eighteen patients were randomised to one of three treatments: placebo, pyridoxine or thiotepa. For each patient, the time of recurrence tumours were recorded until the death or censoring times. The number of recurrences ranges from 0 to 10. Two patients who were censored before the beginning of the study are removed. The dataset is therefore composed of 116 patients with 47 patients from the placebo group, 38 from the thiotepa group and 31 from the pyridoxine group. On these patients, since 13.79% experienced at least five tumour recurrences and only 6.9% patients experienced six tumour recurrences or more, we set the parameter B to 5. For interpretation purpose, the treatment variable is coded as two new binary variables, pyridoxine and thiotepa, making placebo the reference. In addition to these two treatment variables two supplementary covariates were recorded for each patient:

the number of initial tumours and the size of the largest initial tumour.

Tables 3,4 and figure 1 display the estimates obtained from the constant coefficient, unconstrained, total variation and two steps total variation estimators in the multiplicative model. The unconstrained estimator shows very strong variations and is difficult to interpret as such. On the other hand, the constant coefficient estimator gives valuable information on the impact of each covariate, but in turn cannot detect a change in variation. Our total-variation estimators reach compromise: they are not constant but easily interpretable.

For instance, a remarkable aspect of the pyrodixine treatment can be highlighted from the total variation estimation: this treatment produces a protective effect for the first three tumour recurrences but the hazard rate of further recurrences are increased by this treatment. In the same way, an increase in the effect of the initial number of tumours on recurrences is observed from the third recurrence. On the opposite, the effects of the thiotepa treatment or the size of the largest tumour are shown to be constant in the total variation model, the parameter estimates having values similar to the ones obtained in the constant model.

Our conclusions on the treatments effects are in agreement with previous studies on bladder tumours recurrences. For instance, no difference in the rate or time to tumour recurrence was found from patients using pyrodixine with patients using placebo in [START_REF] Tanaka | Pathobiology and chemoprevention of bladder cancer[END_REF] and [START_REF] Goossens | Designing the selenium and bladder cancer trial (seleblat), a phase lll randomized chemoprevention study with selenium on recurrence of bladder cancer in belgium[END_REF]. Moreover, [START_REF] Huang | Marginal regression of gaps between recurrent events[END_REF] and [START_REF] Sun | The additive hazards model for recurrent gap times[END_REF] have respec-tively studied gap time recurrences in the multiplicative and additive models. The results obtained from the former showed a small protective effect of this treatment while the latter concluded that gap times did not seem related to pyridoxine. These examples illustrate the nice features of our total-variation estimator: it provides sharper results, giving relevant informations on covariates effect with respect to the number of recurrent events experienced by a subject and it provides the ability to detect a change of variation. Further details are provided in Supplementary Material.

Discussion

In this paper, the Aalen and Cox models were studied to model the effect of covariates on the rate function. However, such models are not essential in our approach. Penalized algorithms could be easily derived for other models such as the accelerated failure time model or the semiparametric transformation model for instance.

Although we have only presented asymptotic theoretical results, the simulation studies show clear evidence that our estimators outperform standard estimators for small sample sizes. Therefore, it would be of great interest to study their finite sample properties. However, such results involve deviation inequalities for non i.i.d. and non martingale empirical processes. To our knowledge, no such results have yet been established in the context of recurrent events.

Another development of the present paper would be to establish results for the estimation of change-point locations and the number of change-points. Such results can be found for the change-point detection in the mean of a gaussian signal in Harchaoui & Lévy-Leduc (2010), for instance.

Tables Table 1. Simulation results in the multiplicative model for p obs = 14%

FIGURE CAPTIONS

  fp( βm ) = Card j ∈ {1, . . . , p} such that tv( βj m ) = 0 and tv(β j 0 ) = 0 tp( βm ) = Card j ∈ {1, . . . , p} such that tv( βj m ) = 0 and tv(β j 0 ) = 0 and false (fn) and true (tn) negatives by exchanging = and = in the above definitions (see (2.8)for the definition of tv). The specificity and sensibility of a method are defined over the M ) + fn( βm ) ,

Fig. 1 .

 1 Fig.1. Estimates for the bladder data in the multiplicative model. The crosses represent the constant estimator, the filled circles the unconstrained estimator, the circles the total variation estimator and the squares the two steps total variation estimator.

Fig. 1 .

 1 Figures

  mean square error, spec: specificity, sens: sensitivity.

	n	Unconstrained		Constant			tv		two-steps tv
		mse	spec sens	mse	spec sens	mse	spec	sens	mse	spec	sens
	50	5576.511	0	1	225.077	1	0	72.732 0.271 0.813 67.697 0.598	0.68
	100	64.231	0	1	216.658	1	0	46.484 0.226 0.844 39.562 0.583 0.709
	500	12.447	0	1	212.578	1	0	17.232	0.18	0.911 17.998 0.917 0.559
	1000	9.06	0	1	213.292	1	0	14.215 0.192	0.9	17.353 0.983 0.512
			mse:

Table 2 .

 2 Simulation results in the additive model for p obs = 14% mean square error, spec: specificity, sens: sensitivity.

	n	Unconstrained		Constant			tv		two-steps tv
		mse	spec sens	mse	spec sens	mse	spec sens	mse	spec	sens
	50	1208.174	0	1	398.849	1	0	367.377 0.312	1	480.105 0.601 0.992
	100	534.269	0	1	360.757	1	0	221.454 0.241	1	283.258 0.582	1
	500	202.669	0	1	339.446	1	0	139.481 0.154	1	171.794 0.525	1
	1000	168.751	0	1	337.813	1	0	133.39	0.103	1	157.899 0.471	1
			mse:

Table 3 .

 3 Unconstrained and constant parameters estimates for the bladder data in the multiplicative model

	s		Unconstrained			Constant
		pyridoxine thiotepa	size	number pyridoxine thiotepa size number
	1	-0.497	-0.711	-0.028	0.202	-0.037	-0.374	0.03	0.155
	2	0.466	0.013	0.044	0.014	-0.037	-0.374	0.03	0.155
	3	-0.211	0.027	0.129	0.250	-0.037	-0.374	0.03	0.155
	4	0.717	-0.095	0.064	0.274	-0.037	-0.374	0.03	0.155
	5	0.657	-0.283	0.072	0.198	-0.037	-0.374	0.03	0.155

Table 4 .

 4 Total variation and two-steps total variation parameters estimates for the bladder data in the multiplicative model

	s		TV				two-steps TV
		pyridoxine thiotepa	size	number pyridoxine thiotepa	size	number
	1	-0.080	-0.373	0.023	0.133	-0.167	-0.387	0.031	0.122
	2	-0.080	-0.373	0.023	0.133	-0.167	-0.387	0.031	0.122
	3	-0.080	-0.373	0.023	0.133	-0.167	-0.387	0.031	0.265
	4	0.241	-0.373	0.066	0.241	0.625	-0.387	0.031	0.265
	5	0.241	-0.373	0.066	0.241	0.625	-0.387	0.031	0.265

Acknowledgement

The authors are very grateful to the associate editor and the referees for their valuable comments which have helped improve the original manuscript.

Supplementary material

Supplementary material contains proofs of Theorems 3.1 and 3.2. They also include comments on the asymptotic distribution of our estimators, details of the simulation scheme used in Section 4, an extended simulation study and additional analysis on the bladder tumour data of [START_REF] Byar | The veterans administration study of chemoprophylaxis for recurrent stage 1 bladder tumors: comparison of placebo, pyridoxine, and topical thiotepa. Bladder Tumors and Others Topics in Urological Oncology[END_REF].


Aurélien Bechler

Liliane Bel

Mathieu Vrac

Conditional simulations of the extremal t process: application to fields of extreme precipitation

Keywords: Max-stable processes, extremal t processes, spatial extremes, conditional simulation, extreme precipitation

The last decade has seen max-stable processes emerge as a powerful tool for the statistical modeling of spatial extremes and there are increasing works using them in a climate framework. One recent utilization of max-stable processes in this context is for conditional simulations that provide empirical distribution of a spatial field conditioned by observed values in some locations. In this work conditional simulations are investigated for the extremal t process taking benefits of its spectral construction.

) for Brown-Resnick and Schlather models is adapted for the extremal t process with some original improvements which enlarge the possible number of conditional points. A simulation study enables to highlight the role of the different parameters of the model and to emphasize the importance of the steps of the algorithm.

An application is performed on precipitation data in the south of France where extreme precipitation events (Cevenol) may generate major floods. This shows that the model and the algorithm perform well provided the stationary assumptions are fulfilled.

Introduction

Extreme events are one of the major demonstrations of complex environmental systems. As they are inherently rare and often unexpected, society is not well prepared to face them. One important feature of extreme events is that their probability decreases rapidly with increasing intensity, while damage increases strongly as does the cost of protection against such events. Both analyses of observational data and projections from climate models suggest that the occurrence and magnitudes of some extreme events will increase in the future (IPCC (2011)). A potential consequence is an increasing number of catastrophes such as flooding, storms or crop failures. In addition to their direct effects and costs in terms of human lives, building destruction or food security, these events will also have indirect impacts such as increasing costs for building or strengthening infrastructures, in addition to higher insurance premiums, and before even mentioning the consequences of global population displacements (IPCC (2011)). A better understanding and modeling of extreme events is therefore a crucial need in many economic aspects and societal contexts for both public and private decision makers.

From a mathematical point of view, the suitable probabilistic framework for modeling the frequency and intensity of those special events is the statistical theory of extreme values. The extreme value theory is well established for univariate stationary series. There is now a growing effort to develop methods and models to deal with complex data involved in applications where space-time problems and non-stationarity are the rule rather than the exception.

Multivariate extreme distributions have been used in climate and environmental applications in order to take into account extremal dependence [START_REF] Coles | Statistical methods for multivariate extremes -an application to structural design[END_REF], De Haan and De Ronde (1998), [START_REF] Schlather | A dependence measure for multivariate and spatial extreme values: Properties and inference[END_REF], [START_REF] Fawcett | Estimating return levels from serially dependent extremes[END_REF]). For the study of extremes in a spatial context, max-stable processes are a relevant framework for modeling and inference. Introduced by De [START_REF] De Haan | A spectral representation for max-stable processes[END_REF] different parametric models have been later proposed by [START_REF] Smith | Max-stable processes and spatial extreme[END_REF], [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] or [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]. The extremal t process [START_REF] Padoan | Multivariate extreme models based on underlying skew-t and skew-normal distributions[END_REF], [START_REF] Opitz | Extremal t-processes: Elliptical domain of attraction and a spectral representation[END_REF]) is a model of that family that extends the Schlather model. It has a supplementary parameter that lets the extremal coefficient function vary in a wider range. Max-stable processes have been used in applications for modeling rainfall [START_REF] Buishand | On spatial extremes: with application to a rainfall problem[END_REF], Davison et al. (2012)), temperature [START_REF] Davison | Geostatistics of extremes[END_REF], [START_REF] Fuentes | Nonparametric spatial models for extremes: application to extreme temperature data[END_REF]), snowfall [START_REF] Gaume | Mapping extreme snowfalls in the french alps using max-stable processes[END_REF]) or snow depth [START_REF] Blanchet | Spatial modeling of extreme snow depth[END_REF]. A main objective of these works is to provide return periods of rare events for prevention purposes even in locations where no observations are available.

The aim of the present paper is to further investigate conditional simulations of extreme precipitation. Conditional simulations provide empirical distributions of a random field knowing the values observed at some given locations. Consequently, they allow us to statistically generate extreme values at locations without any recording, in agreement with surrounding observations. This capacity to generate time series at observation-free locations is of major importance for some impact models, such as hydrological models that require input data, e.g. precipitation values, correctly spatially distributed. The spatial coverage of the meteorological inputs of such models has to be as precise as possible. Furthermore, when other techniques provide individual values or indicators, conditional simulations provide distributions at every new locations which allows us to get numerous statistical quantities as e.g. quantile, dispersion indicators, etc... The conditional simulation techniques are then of major importance.

For Gaussian and Boolean fields, algorithms for conditional simulations are well established [START_REF] Lantuéjoul | Geostatistical simulation: models and algorithms[END_REF]). They have been used for representing physical phenomena by integrating heterogeneous data and assigning confidential limits to the estimates. For example, they have been used for predicting the geometry of oil reservoirs or the grade of mining blocks [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF]) and are used more and more in the environmental sciences [START_REF] Mariethoz | Truncated plurigaussian simulations to characterize aquifer heterogeneity[END_REF]). Concerning max-stable processes, [START_REF] Wang | Conditional sampling for spectrally discrete max-stable random fields[END_REF] consider max-linear random fields, a special case of max-stable processes, and provide an exact expression of the conditional distributions as well as an efficient algorithm for the simulation. In a further work, [START_REF] Dombry | Regular conditional distributions of continuous maxinfinitely divisible random fields[END_REF] derive the regular conditional distribution for the very general class of max-infinitely divisible processes which encompasses the max-stable processes. This result has been applied then to provide a methodology and an algorithm for Brown-Resnick and Schlather processes (Dombry et al. (2013)). However, for processes which are not regular, that is the exponent measure is not absolutely continuous with respect to Lebesgue measure, the conditional distribution is not explicite. [START_REF] Oesting | Conditional sampling for max-stable processes with a mixed moving maxima representation[END_REF] discuss the frameworks where [START_REF] Wang | Conditional sampling for spectrally discrete max-stable random fields[END_REF] and Dombry et al. (2013) can be applied, and, in case of non regular models, propose an algorithm based on martingale arguments.

In this paper we adapt the methodology of Dombry et al. (2013) to the extremal t process. In this case all the calculations of the algorithm can be derived explicitely, nevertheless some steps still present computational issues. In particular, one of the steps of the algorithm involves an acceptance/rejection procedure for which the rate of rejection increases very fast with the number of conditioning points. This limiting point makes the algorithm inapplicable when the number of conditioning points is greater than approximately 30. To bypass the issue, we propose some original modifications of the algorithm to replace the acceptance/rejection step by a Gibbs sampler which speeds up the calculations and enables a large number of conditioning points. Despite this improvement, some difficulties persist in the algorithm, mostly due to the calculation of high dimensional integrals that have to be approximated. Then to get a better understanding of the algorithm mechanisms, a simulation study is led to assess the influence of the parameters from the extremal t process and to evaluate the importance of the most problematic and time-consuming step of the algorithm.

We focus on the French Mediterranean basin, which is a key spot of occurrences of meteorological extremes such as heavy precipitation (e.g. [START_REF] Ducrocq | A numerical study of three catastrophic precipitating events over southern france. ii: Mesoscale triggering and stationarity factors[END_REF]; [START_REF] Lebeaupin | Sensitivity of torrential rain events to the sea surface temperature based on high-resolution numerical forecasts[END_REF]; [START_REF] Beaulant | Statistico-dynamical downscaling for mediterranean heavy precipitation[END_REF]; [START_REF] Vrac | Weather regimes designed for local precipitation modeling: Application to the mediterranean basin[END_REF] or [START_REF] Vrac | Dynamical and statistical downscaling of the french mediterranean climate: uncertainty assessment[END_REF]). Indeed, major extreme precipitation events are regularly observed in this region near the Cevennes mountains. The modeling and the forecast of this extreme precipitation -the so-called "Cevenol events" -are of major importance for hydrological studies in this complex terrain since they often trigger major floods in this region (e.g. [START_REF] Ducrocq | A numerical study of three catastrophic precipitating events over southern france. ii: Mesoscale triggering and stationarity factors[END_REF]).

The structure of this article is the following. In Section 2, a short reminder on max-stable processes is provided and the extremal t process is described. Section 3 contains theoretical considerations on conditional simulations of an extremal t process, as well as the associated algorithm and its improvement. An extensive simulation study is carried out in Section 5 to evaluate the performances of the algorithm. The use of the continuous ranking probability score and the quantile score [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]) enables us to assess the accuracy of our methodology and to compare it to other models. This evaluation allows us to quantify the importance of the model parameters and of the algorithm steps. Section 6 displays an application to maximum precipitation data in the south of France where extreme precipitation may cause major floods. Some comments and conclusions are given in section 7.

Extreme Value Theory

Max-Stable Process

We will note vectors in plain text and vectors of vectors in bold. Let us consider a stationary random field Z = {Z(x)} x∈R d with a given distribution H 0 and (Z 1 , . . . , Z n ) independent copies of Z. If there exist sequences

a n (x) ∈ R and b n (x) ∈ R + such that for all n ∈ N, max i=1,...,n Z i (x) -a n (x) b n (x) x∈R d (1)
follows the distribution H 0 , then Z = {Z(x)} x∈R d is a max-stable process [START_REF] De Haan | A spectral representation for max-stable processes[END_REF]).

As a consequence, all finite dimensional marginal distributions are max-stable, and in particular, the univariate margins Z(x) distribution belongto the Generalized Extreme Value (GEV) family:

P(Z(x) ≤ z) =        exp -1 + ξ(x) z -µ(x) σ(x) -1/ξ(x) if ξ(x) = 0, exp -exp z -µ(x) σ(x) if ξ(x) = 0, (2)
where µ(x), σ(x) and ξ(x) are respectively the location, scale, and shape parameters at site x. Furthermore, the sign of ξ(x) will determine the family distribution Z(x) belongs to among Fréchet (ξ > 0), Gumbel (ξ = 0) or Weibull (ξ < 0), and consequently the distribution tail characteristics.

Different max-stable models have been proposed and used in the past by several authors, e.g. [START_REF] Brown | Extreme values of independent stochastic processes[END_REF], [START_REF] Smith | Max-stable processes and spatial extreme[END_REF], [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] or [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]. Their margins are Fréchet distributed, i.e.

P(Z(x) ≤ z) = e -z -1 .
(3)

The extremal t process: definition and properties

The extremal t process [START_REF] Padoan | Multivariate extreme models based on underlying skew-t and skew-normal distributions[END_REF]) belongs to the max-stable process family. [START_REF] Opitz | Extremal t-processes: Elliptical domain of attraction and a spectral representation[END_REF] proposed a spectral construction of the extremal t process which leads to the general expression

Z = {Z(x)} = m -1/α α max i≥1 ξ i W i (x) , α > 0, x ∈ X ⊂ R d , with X a Borel set, (4) with m α = √ π -1 2 (α-2
)/2 Γ((α + 1)/2), ξ i a Poisson point process on (0,+∞) with intensity measure dΛ(t) = αt -(α+1) dt and W i independent copies of a stationary standard Gaussian process on R d with correlation function ρ(•). Z margins are α-Fréchet distributed, i.e. P(Z(x) ≤ z) = e -z -α .

(5)

Unlike the usual representation of max-stable processes with unit-Fréchet margins, we choose the α-Fréchet representation (as in Opitz ( 2013)) to highlight the major role of the α-parameter.

It has to be noted that for α = 1, this representation leads to a Schlather process and hence the extremal t process can be considered as a generalization of the Schlather process, with a greater flexibility.

The parameter α and the correlation function ρ(•) determine the dependence structure. Furthermore, as will be seen in the next section, α characterizes the spatial asymptotic dependence, that is the dependence in extreme values.

Extremal coefficient

For a max-stable process such that all margins are equally distributed the extremal coefficient function θ(•) [START_REF] Schlather | A dependence measure for multivariate and spatial extreme values: Properties and inference[END_REF]) is defined as

P(Z(x) < u, Z(x + h) < u) = P(Z(x) < u) θ(h) . (6)
θ(h) ∈ [1, 2] measures the strength of asymptotic dependence, θ(h) = 1 corresponds to the complete dependence case, θ(h) = 2 corresponds to the independence case. The extremal coefficient of the extremal t process depends on α and is given by [START_REF] Opitz | Extremal t-processes: Elliptical domain of attraction and a spectral representation[END_REF])

θ(h) = 2t α+1 ( 1 | ρ(h), (α + 1) -1 (1 -ρ(h) 2 )), (7)
with t ν (.|µ, Σ) the cdf (cumulative distribution function) of a Student distribution with the parameters µ, Σ and ν (respectively location, shape and degree of freedom parameters).

In particular for the extremal t process with α ≥ 0, θ ∞ , the limit of θ(h) when h tends to ∞, belongs to the interval ]1.5, 2[, while θ ∞ = 1 + √ 2/2 (≈ 1.7) for the Schlather process.

Consequently, thanks to this supplementary parameter α, the extremal t process enables us to model a greater variety of situations for asymptotic dependence. Furthermore, this flexibility of the asymptotic dependency behavior is different from the one given by the class of general Brown-Resnick processes [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]).

Conditional simulations

Conditional simulations provide simulations of random fields at locations without observation, given observations in fixed locations, according to a model. The simulated fields are required to fit the observed values at the fixed locations. In contrast to simple simulations for which simulated values follow a target distribution -i.e. Z is a realization from a distribution L -, conditional simulations take into account constraints expressed by the observations, i.e.

Z is a realization from a distribution L, Z(x i ) = z i , ∀z i observations at the locations x i , i = 1, . . . , k.

Conditional simulations are powerful tools which provide not only individual values at new locations, but probabilistic distributions at those locations, given a chosen model. Those distributions allow the calculation of many statistical quantities as for example quantiles, dispersion indicators, etc. . . . For some models, such as Gaussian or Boolean, algorithms for performing conditional simulations are quite straightforward and based on kriging techniques [START_REF] Lantuéjoul | Geostatistical simulation: models and algorithms[END_REF]). Dombry et al. (2013) proposed a general framework for conditional simulations of max-stable processes. The max-stable process Z is written as the point-wise maximum over an infinity of continuous functions

Z(x) = max i≥1 φ i (x). (8) Let x = (x 1 , x 2 , . . . , x k ) ∈ X k the k conditioning points with X ⊂ R d , z = (z 1 , z 2 , . . . , z k ) ∈ R k such that Z(x i ) = z i and s = (s 1 , s 2 , . . . , s m ) ∈ X m
the m locations the process Z is to be simulated. For each conditioning point x i , an unique function φ j almost surely exists, called extremal function and denoted φ + j such that φ + j (x i ) = z i . One extremal function may be extremal for several conditioning points. Let τ = (τ 1 , . . . , τ l ) be a partition of the set {x 1 , . . . ,

x k } with l = |τ | (≤ k) such that φ + j (x i ) = z i if x i ∈ τ j and φ + j (x i ) < z i if x i / ∈ τ j , (i = 1, . . . , k; j = 1, . . . , l). (9)
The partition τ is called a hitting scenario [START_REF] Wang | Conditional sampling for spectrally discrete max-stable random fields[END_REF]). The conditional simulation algorithm is then divided into three steps. The first one consists in drawing a partition of {x 1 , . . . , x k } such that each class of the partition gathers conditioning points that have the same extremal function. Then, for each class in the drawn partition, each extremal function φ + j is drawn according to the constraints (9) given by the conditional observations. Finally a number N sufficiently large of non extremal functions denoted φ - j are drawn such that for all conditioning points they verify

φ - j (x i ) < z i for all x i (i = 1, . . . , k; j = 1, . . . N ). ( 10
)
The simulated values on the set s are obtained by taking the point-wise maximum of {φ - j (s), φ + j (s)}. To perform the three steps of the algorithm, the conditional distribution of the partition, the {φ + } and the {φ -} are needed. Dombry et al. (2013) derived the expression of the conditional distributions for the Schlather and the Brown-Resnick process. The conditional distribution for the extremal t process is obtained thanks to a slight adaptation.

Conditional distribution for the extremal t process

Let Z be an extremal t process. It can be written as

Z(x) = max i≥1 φ (α) i (x) (11) with φ (α) i (x) = m -1/α α ξ i W i (x), i = 1, 2 . . . . (12)
Φ (α) = {φ (α) i } i≥1 is a Poisson point process on C 0 , the space of positive continuous func- tion on X ⊂ R d . The intensity measure Λ (α) x of the point process Φ (α) (x) = {φ (α) i (x)} i≥1 on (0, +∞) k is Λ (α) x (A) = +∞ 0 P[Φ (α) (x) ∈ A]ξ -(α+1) dξ, A ⊂ (0, +∞) k a Borel set. ( 13
)
The extremal t process fulfills the conditions necessary for Λ (α)

x being absolutely continuous with respect to the Lebesgue measure [START_REF] Dombry | Regular conditional distributions of continuous maxinfinitely divisible random fields[END_REF]). The density (Radon-Nikodym derivative) λ

(α) x (z) such that Λ (α) x (dz) = λ (α)
x (z)dz and the conditional density

λ (α) s|x,z (u) = λ (α) s,x (u, z) λ (α) x (z) , (s, x) ∈ X m+k , z ∈ (0, +∞) k , u ∈ R m . ( 14
)
are given by :

λ (α) x (z) = π α-1 2 w k+α-2 |Σ x | -1/2 [z T Σ -1 x z] -( k+α 2 ) Γ((α + 1)/2) -1 (15)
and

λ (α) s|x,z (u) = w k+α-2 w k+m+α-2 |Σ s,x | -1 |Σ x | [z T Σ -1 x z] (k+α)/2 [(u, z) T Σ -1 s,x (u, z)] (k+m+α)/2 , (16)
with

w d = π d/2 /Γ(d/2 + 1), | . | the determinant operator and Σ x , Σ s,x = Σ s Σ s:x Σ x:s Σ
x the covariance matrices of respectively W (x) and W (s, x). For the details of the calculations see Appendix A.

For α = 1, we note that the expressions ( 15) and ( 16) are the ones obtained for the Schlather process (Dombry et al. (2013)). They are equivalent to those obtained by Ribatet (2013a).

Conditionally to Z(x) = z, the Poisson point process Φ (α) can be decomposed into two independent Poisson processes, such that Φ

(α) = Φ (α)+ x,z ∪ Φ (α)- x,z , with Φ (α)- x,z = {φ (α) ∈ Φ (α) | ∀i = 1, . . . , k, φ (α) (x i ) < z i } = {φ (α)- j } j≥1 , Φ (α)+ x,z = {φ (α) ∈ Φ (α) | ∃i = 1, . . . , k, φ (α) (x i ) = z i } = {φ (α)+ j } j≥1 .
(17)

Algorithm

The algorithm proposed by Dombry et al. (2013) for the conditional simulation of Brown-Resnick and Schlather processes can easily be adapted for extremal t processes. It is based on the density λ (α)

x and the decomposition of Φ (α) . The following notations have been chosen to be consistent with notations of Dombry et al. (2013). Let us note x τ j = {x i | x i ∈ τ j , i = 1, . . . , k}, z τ j = {z i | Z(x τ j ) = z i , i = 1, . . . , k} and τ c j the complementary of τ j in {x 1 , . . . , x k }.

Step 1 Draw a random partition τ of the set (x 1 , . . . , x k ) conditionally to Z(x) = z, according to the distribution:

P[T = τ | Z(x) = z] ∝ |τ | j=1 λ (α) xτ j (z τ j ) {u j <z τ c j } λ (α) x c τ j |xτ j ,zτ j (u j )du j . ( 18
)
Step 2 Given a partition τ = (τ 1 , . . . , τ l ), draw l independent extremal functions (φ

(α)+ 1 , . . . , φ (α)+ l ) at s with finite dimensional distribution P[φ (α)+ j (s) ∈ dv | Z(x) = z, T = τ ] ∝ 1 {u j <z τ c j } λ (α) (s,x c τ j )|xτ j ,zτ j (v,u j )du j dv. (19) Set Z + (s j ) = max(φ (α)+ 1 (s j ), . . . , φ (α)+ l (s j )) for j ∈ 1, . . . , m.
Step 3 Independently of the two previous steps, draw a number N sufficiently large of φ (α)-at (s, x) such that φ (α)-(x) ≤ z and set

Z -(s j ) = max i≥1 {φ (α)- i (s j )} for j ∈ 1, . . . , m, (20)
Then Z(s) = max{Z + (s), Z -(s)} is distributed according to the conditional distri- bution of Z(s) given Z(x) = z.
Step 1 is the most challenging one because of the tremendous number of partition candidates. Indeed the number of possible partitions is equal to the Bell number [START_REF] Spivey | A generalized recurrence for bell numbers[END_REF]). This number grows very fast with the number of conditioning points. For instance, for k = 10, it is already equal to 115 975. Direct calculation of the probability in equation ( 18) is therefore not possible for a realistic number of conditioning points. In order to bypass this computational issue, Dombry et al. (2013) propose drawing the partition using a Gibbs sampler [START_REF] Casella | Explaining the gibbs sampler[END_REF]) which changes the class partition of one observed value at a time. Then once the transition probabilities involved simplified, the chain runs through the space of all possible partitions and converges to the target law described in equation ( 18).

Nevertheless, there are still some issues that need to be solved. The major problem is that the integrals involved in the transition distribution (see equation ( 18)) have to be calculated at each step of the Gibbs sampler. The dimension of these integrals can be as high as (k -1). The values of the integrals can be estimated by classic Monte-Carlo methods or quasi Monte-Carlo methods [START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF]) but in both cases, the estimation requires a huge number of points to get satisfying results and the Gibbs sampling can become computationally heavy.

For Step 2, it is necessary to simulate from the distribution λ (α) s|x,z . Using the same arguments as in Dombry et al. (2013), it can be shown that λ

(α) s|x,z (u), u ∈ R m is the density of a multivariate Student distribution with (k + α) degrees of freedom, with location parameter µ = Σ s:x Σ -1 x z, (21)
and scale matrix

Σ = a x (z) k + α (Σ s -Σ s:x Σ -1 x Σ x:s ) with Σ s,x = Σ s Σ s:x Σ x:s Σ x and a x (z) = z T Σ -1 x z. (22)
Therefore 19) is a truncated multivariate Student distribution and step 2 results in a simulation under constraints.

1 {u j <z τ c j } λ (α) (s,x c τ j )|xτ j ,zτ j (v,u j ) involved in the integral of (
A procedure to perform this simulation is proposed and detailed in the following section.

Step 3 is the most accessible one. It consists in simulating values for an extremal t process at every conditioning and unconditioning location and check whether they respect the constraints given by the conditioning points using an acceptance/rejection procedure. In this context, the rejection rate of the algorithm is not very high even for a large number of conditioning points. Indeed, it rarely happens that one of the φ (α)i becomes greater than the conditioning value, because in this model, the conditioning values are meant to be, by construction, extreme realizations of these functions.

To simulate Z -, especially the Poisson point process ξ, we adapted the simulation techniques described in [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] for the Schlather process by noting that

ξ = n i=1 e i -1/α n∈N , (23)
with e i ∼ E(1) a standard exponential distribution, is a Poisson point process with intensity measure αt -(α+1) dt.

By construction the sequence ξ is decreasing and for i large enough ξ i W i (.) will not participate to the maxima Z -for any x ∈ X provided that W i (x) < C for all x ∈ X .

C can be set to 5, and we choose the number N of sub-extremal functions such that the probability that the last built one contributes to the maximum is lower than 0.001. This number depends on α and increases with it.

Multivariate Student under constraints

A straightforward way to sample from a truncated multivariate Student distribution is to simulate from a multivariate Student distribution and add an acceptance/rejection procedure (Dombry et al. (2013)). Nevertheless, the rate of acceptance of this algorithm becomes very low for a number of conditioning points higher than 30. We propose using the general methodology developed by [START_REF] Geweke | Efficient simulation from the multivariate normal and student-t distributions subject to linear constraints and the evaluation of constraint probabilities[END_REF] for sampling a p-dimensional vector y from a multivariate Student distribution with p linear constraints i.e. y ∼ t ν (µ, Σ), with a ≤ Dy ≤ b, with D a full rank p × p matrix and some elements of the p-dimensional vectors a and b can be equal respectively to -∞ and +∞. This method involves two steps simulating alternatively the denominator and the numerator of the multivariate Student random vector. The denominator is a constrained χ 2 variable simulated via an acceptance/rejection procedure (on the ratio between the multivariate normal realization and the χ 2 variable) and the numerator is updated sequentially according to a Gibbs sampler. The whole procedure is repeated until convergence.

It has to be noted that the rejection rate of the χ 2 variable simulation is low because the constraint is univariate, and that the Gibbs sampler avoids using an acceptance/rejection procedure on vectors.The adaptation of this methodology to the step (2) (equation ( 19)) of the algorithm is straightforward setting for each j ∈ {1, . . . , l}, all a components to -∞, D to the (m + #(τ c j ))-dimensional identity matrix , the first m components of b to +∞ and the #(τ c j ) last ones to z τ c j . Consequently, all tools are available to perform conditional simulations of extremal t processes with this algorithm.

Comparison criteria

For assessing the quality of the results, criteria comparing distributions resulting from the conditional simulations and a value corresponding to an observation have to be used.

The accuracy of conditional simulations is assessed using first the continuous ranking probability score crps [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]). Letting F be the predictive cumulative distribution function derived from the n conditional simulations {ŷ (j) } j=1,...,n and y the test observation, the continuous ranking probability score is written

crps(F ; y) = +∞ -∞ F (t) -H(t -y) 2 dt, (24)
where H(u) is the Heavyside function which is equal to 0 for u < 0 and 1 otherwise.

The steeper the cdf is around the real value, the better. A good score is a score close to 0. It also appears that the calculation of the crps will result in a value with the same unit as the variable of interest. Generally, there is no closed-form expression for this integral. Nevertheless, when F denotes a GEV distribution, a closed-form expression exists (Friederichs and Tho-rarinsdottir (2012)).

Moreover, the conditional simulations are also evaluated based on the quantile score (qs) defined as

qs p (F ; y) = ν p (y -F -1 (p)) ( 25
)
where p is the probability of interest and ν p (u) = pu if u ≥ 0 and ν p (u) = (p -1)u otherwise. As shown by [START_REF] Friederichs | Forecast verification for extreme value distributions with an application to probabilistic peak wind prediction[END_REF], the crps is proportional to the integral of the qs p over all possible probabilities p. Then a good qs p is a qs p close to 0 as well. By choosing a high probability p, the quantile score focuses on the predictive performance in the upper tail [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]).

Then, overall scores are built from equations ( 24) and ( 25), by averaging on all points:

CRPS((F 1 , . . . , F m ); (y 1 , . . . , y m )) = 1 m m i=1 crps(F i ; y i ). QS p ((F 1 , . . . , F m ); (y 1 , . . . , y m )) = 1 m m i=1 qs p (F i ; y i ). ( 26
)
In addition to these two univariate scores, which assess the marginal qualities of the conditional simulation, we propose a multivariate score: i |i ∈ I}. This score enables us to assess the spatial coherency of the conditional simulations.

In order to compare models, skill scores can be built [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]):

S skill = S mod -S ref S opt -S ref (28)
with S mod , S ref and S opt the score of the current model, the score of a reference model and the optimal score respectively. For models that improve on the reference model, S skill belongs to [0, 1] with 1 correspond to the optimal case. For CRPS, QS p and MS d , CRPS opt = QS which corresponds to the percentages of CRPS, QS p or MS d improvement brought by the current model with respect to the reference model.

Simulation Study

In order to check the properties of the conditional simulations provided by this algorithm, simulations according to extremal t processes are drawn, and conditional simulations at non conditioning points (test sample) are compared to their real values. 100 points are sampled at random in a unit square. 20 are kept as conditioning points, and the other 80 are used as test values. For all simulations the covariance function is the exponential with scale parameter 1.

Influence of the partition

As stated in the section 3, the first and hardest step for conditional simulation is to draw a partition of the conditioning points according to the distribution given in equation ( 18). Indeed, it appears that the Gibbs sampler converges very slowly to the target law and then it is likely that we do not draw the partition from the target distribution. We investigate how the resulting partition distribution, obtained at different convergence stages of the Gibbs sampler, influences the performance of the whole algorithm. For each α ∈ {2, 4, 6}, we simulate extremal t fields and compare the performance of the algorithm when (i) the real partition is known, (ii) when it is drawn from a Gibbs sampler (10 000 iterations) according to the distribution ( 18) and (iii) when the partition is drawn randomly. 1: Importance of a good estimation of the partition: values of the CRPS and CRPSS for α = 2, 4, 6, (i) when the real partition is known, (ii) when the partition is drawn from a Gibbs sampler (10 000 iterations) according to the distribution ( 18) and (iii) when the partition is drawn randomly.

When the partition is drawn from the Gibbs sampler, the CRPSS is slightly greater than 5% for α=2 and α=4 and around 1% for α=6. However, it is greater than 30% when the real partition is used. This demonstrates the strong gain to be had by drawing close to the real one and confirms that the efforts made for this step of the algorithm are inevitable and should even be developed. Indeed, the gain in drawing the partition with the Gibbs sampler is not very high compared to using a random partition (less than 6%) and far from the 30% obtained with the real partition. We can think that drawing randomly the partition instead of performing a Gibbs sampler for time-saving reasons for instance, would not change considerably the results especially when the number of conditioning points is high.

Influence of α

We investigate how the value of the parameter α influences the quality of the conditional simulations. For each value of α, 15 simulations of an extremal t process are drawn and 20 conditioning points are sampled for each simulation. Then 50 conditional simulations are performed according to each conditioning point set. We suppose that the real partition is known, just as α and ρ(.) which are not estimated.

The conditional simulations of the extremal t process are compared to conditional simulations of a Schlather process [START_REF] Schlather | Models for stationary max-stable random fields[END_REF]), a Brown-Resnick process [START_REF] Brown | Extreme values of independent stochastic processes[END_REF]), and a Gaussian process with the same conditioning points and unconditional simulations of the extremal t process. Ordinary Kriging [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF]) is set as the reference prediction for the CRPSS calculation. The conditional simulations of the Schlather and Brown-Resnick processes are performed according to the algorithm of Dombry et al. (2013) and with SpatialExtremes (Ribatet (2013b)). For the Gaussian process conditional simulations, we apply first a Gaussian margin transformation and then perform the methodology described in [START_REF] Lantuéjoul | Geostatistical simulation: models and algorithms[END_REF] for Gaussian processes (even if it is only marginally Gaussian). For the ordinary kriging and the Gaussian simulations, the covariance function is estimated from the data.

Figure 1 describes the evolution of the different CRPSS with α. As expected since it is based on the right model, the extremal t model gives better CRPSS than do the other models. It also appears that the unconditional extremal t simulations become more and more efficient as α increases. Indeed, when α increases, the variance of the extremal t process tends to 0 [START_REF] Opitz | Extremal t-processes: Elliptical domain of attraction and a spectral representation[END_REF]) and the simulated samples show very little variability, making the information added by the conditional values less influential. It is worth noticing that Gaussian margin transformation followed by conditional simulations of Gaussian fields shows surprisingly good results and it appears to be very efficient specially when α is high. The relevance of this approach in conditional simulation framework has also been stated in [START_REF] Oesting | Conditional sampling for max-stable processes with a mixed moving maxima representation[END_REF]. q q q q q q 1 2 3 4 5 6

-30

-20

-10 0 10 20 30 40 α CRPSS q q q q q q q q q q q q q q q q q q q q q q q q t extremal Schlather Brown-Resnick Gaussian Non conditional t extremal Nevertheless, for every values of α, the CRPSS of the extremal t simulations is around 20% showing more stability in good performances than any other tested models.

Asymptotic dependence properties

We now investigate a notion related to the asymptotic dependence for the conditional simulations. As mentioned by Dombry et al. (2013), Z(.)|{Z(x) = z} is not max-stable and cannot be asymptotically dependent. However, the process integrated on all the conditional events recovers the max-stability property, and its asymptotic dependence can be characterized by the madogram [START_REF] Cooley | Variograms for spatial max-stable random fields[END_REF]). Using the madogram we compare the pairwise extremal coefficient estimates θ est to the theoretical extremal coefficient function θ th (h) at a distance h in order to check that the asymptotic dependence behavior of the conditional simulations is similar to that of the theoretical one.

For each α ∈ {2, 4, 6}, we simulate 1 000 extremal t fields at the same 20 conditional points and estimate the extremal coefficient function from the remaining 80 points with the R package SpatialExtremes (Ribatet (2013b)). Figure 2 compares the theoretical extremal function θ th (h) and the estimated one θ est (h) calculated from the madogram of the integrated conditional simulations.

For the three values of α, the conditional simulations respect the expected asymptotic dependence behavior. Indeed, the theoretical extremal coefficient function fits relatively well the pairwise extremal coefficient estimates, although it tends to overestimate them when α becomes high. The choice of the autumnal period allows us to make the temporal stationarity assumption, as it is commonly accepted that rainfall data can be considered as stationary within a given season.

In order to perform the conditional simulations presented in section 3 the margins have first to be transformed into an α-Fréchet distribution. If Z is distributed according to a GEV(µ, σ, ξ), then Y defined as

Y = 1 + ξ(Z -µ) σ 1/αξ + (30)
follows an α-Fréchet distribution [START_REF] Beirlant | Statistics of extremes[END_REF]).

To avoid the estimation of the parameters (µ, σ and ξ) of the GEV for all points and to deal with spatial non-stationarity, these parameters are modeled as functions of the spatial coordinates LAT and LON except ξ that is assumed to be constant for the whole region:

µ = exp µ 0 + µ 1 * LAT + µ 2 * LON σ = exp σ 0 + σ 1 * LAT + σ 2 * LON ξ = ξ 0 (31)
It can be observed that the exponential formulation is employed for µ and σ as σ has to be positive for model definition reasons and µ is chosen to be positive because of the nature of the observed phenomenon (precipitation).

Inference

The parameters of the model, i.e.:

1. µ 0 , µ 1 , µ 2 , σ 0 , σ 1 , σ 2 , ξ, the parameters for the margin model, 2. α, the extremal t parameter, 3. and a the scale parameter of the exponential correlation function of W , have to be estimated from the data. Nevertheless, for dimensions higher than 2, the multivariate density of the extremal t process is not tractable and usual likelihood maximization methods are not available. Alternatively, composite likelihood [START_REF] Lindsay | Composite likelihood methods[END_REF], [START_REF] Varin | An overview of composite likelihood methods[END_REF]) can be used. It consists in calculating bivariate densities for each pair of points and maximizing their weighted sum. Let us note {z m,x i } the autumnal maximum on m = 1, . . . , M years for a given x i in (x 1 , . . . , x k ). The composite log-likelihood is

LC(Ψ; z) = M m=1 k-1 i=1 k j=i+1 w ij l ij (Ψ) (32)
with l ij (Ψ) = logf (z m,x i , z m,x j ; Ψ) the bivariate log-density for points x i and x j with Ψ the vector of parameters to estimate, and w ij the corresponding weights. These weights w ij have been chosen to be decreasing (exponentially) with the distance between x i and x j . The maximization of this composite likelihood leads to consistent estimations of the target parameters [START_REF] Padoan | Likelihood-based inference for maxstable processes[END_REF], [START_REF] Xu | On the robustness of maximum composite likelihood estimate[END_REF]).

Results

The inference has been driven by the data over the whole time period. The maximization of the composite likelihood gives α = 2.2 (with θ ∞ = 1.83) which means that even for long distances there is still some asymptotic dependence. Furthermore, the estimated range of the exponential covariance function is equal to 360 km.The estimated parameters and the confidence intervals are presented in Table B.1.

The validation is performed only for the year 2007 with 100 conditional simulations. The most represented partition in the history of the Gibbs sampler is composed of 5 classes. However, there is no notable spatialization of these classes (not shown).

Figure 4 shows the individual crps scores (as defined in equation ( 24)) computed for the year 2007 taking the observed values as test observations (y i ) and the empirical distributions of the conditional simulations as the predictive cdf (F i ). In addition, Figure 5 describes the distributions of the conditional simulations at four arbitrarily chosen grid cells specified on Figure 4. It appears naturally that the closer to conditioning points the new location is, the better the simulation (e.g. points 2 and 4 are close to conditioning points whereas point 3 is far). However, this is clearly not the only important criterion. Indeed, point 1 is quite close to conditioning points yet its simulations are not as good as expected. In order to evaluate the quality of these conditional simulations, other models and algorithms are considered and compared: For simulations based on the analog model (C) (e.g. [START_REF] Barnett | Multifield analog prediction of short-term climate fluctuations using a climate state vector[END_REF]; [START_REF] Zorita | The analog method as a simple statistical downscaling technique: comparison with more complicated methods[END_REF] or [START_REF] Pascal | Ensemble reconstruction of the atmospheric column from surface pressure using analogues[END_REF]), at each non conditioning location we consider the nearest conditioning grid cell. Then we look at the history of this point to find the year with the closest observed value: the value observed this year at the non conditioning location is taken as the estimate value for the non conditioning point. For the so-called "historical" model (D), we randomly draw a value in the history of a non conditioning location and consider it as the value of the simulation. The historical model appears as the best one regarding the CRPS score. In term of MS 2 , the model analog, which respects by construction the spatial coherency, gives the best results. Model (A) shows good results in terms of QS 95 and generally appears to be the most versatile method adapted to different purposes (for each score, it gets at worst the second best value). It is also important to note that method (A) is better than method (B): there is a real gain in running a Gibbs sampler to estimate the partitions.

Generally, models (A)-(B) are the best to catch high quantile behavior (see QS 95 ). Nevertheless, as crps is the overall score of qs p integrated over all probabilities p, it indicates that for lower probabilities, model (A) is not as good as the others. It appears that for p ≥ 0.75, QS p of model (A) calculated over the whole region is the best.

Furthermore, by deepening the analysis of the results, it appears that for some points the error (in terms of crps) is fairly high, especially for model (A). Those points are located in the upper left corner of Figure 4 (i.e., the northwest of the region), corresponding to the higher region(s) of the Cevennes mountains. The climate characteristics in this area are very different from the rest of the region and it is likely that the chosen covariates (LAT and LON) are not sufficient to overcome the inhomogeneity of this region.

Then, for comparison purposes inside the stationary framework, we removed these points for the calculation of the scores. The CRPSS is now better for method (A), with a gain of 4.4% with respect to the reference model and both QSS 95 and MSS 2 remain good. It also appears that the CRPS and CRPSS have been improved for all methods.

It is also interesting to check whether the conditional simulations follow the temporal variations of the real values. To this aim, we have simulated 50 conditional fields on 50 grid cells from the SAFRAN database, on the period 2000 -2007 with the same estimated parameters. As running the Gibbs sampler for each year with 50 conditional points would be too time-consuming, the partitions are drawn randomly for each year. Figure 6 shows the temporal evolutions of the means of the conditional simulations at the four previously chosen grid cells. For points 2 and 4, the temporal evolutions of the conditional simulations are almost the same as for the test values, whereas for point 3, there is a big discrepancy in the first period and a good match in the second period. A constant significant bias is noted for point 1, which is a location in the already mentioned mountain region, and it can be suspected that it behaves in a quite different way than do the conditioning points even though it approximately follows the temporal variations. The scatterplot of the simulated values vs. the test values is represented in Figure 7. The Pearson correlation coefficient is equal to 0.71. It appears that the conditional simulations roughly follow the intrinsic variations of the data over the years. Indeed, a correlation coefficient of 0.71 is a satisfying result which indicates that the method is able to reproduce the general evolution of the observed phenomenon. In comparison, the correlation coefficient is 0.47 for the historical model and 0.33 for the analog model.

Discussion and conclusions

We have investigated conditional simulations for maxima precipitation using the extremal t model. Thanks to its parametrization, this model offers a great flexibility in modeling the asymptotic dependence. The conditional simulation algorithm designed by Dombry et al. (2013) for Brown-Resnick and Schlather models has been adapted to the extremal t process and the second step, consisting in simulating from a multivariate Student distribution under constraints, has been improved by a two-state Gibbs sampler proposed by [START_REF] Geweke | Efficient simulation from the multivariate normal and student-t distributions subject to linear constraints and the evaluation of constraint probabilities[END_REF] enabling us to handle a greater number of conditioning points.

A simulation study highlighted that the methodology of conditional simulations for the extremal t process is effective for a large range of α (∈ [1,6]). Moreover, when the value of α is too high, its estimation becomes more and more uncertain. In such cases, we would advise the use of the conditional simulations of Gaussian processes with transformed Gaussian margins [START_REF] Oesting | Conditional sampling for max-stable processes with a mixed moving maxima representation[END_REF]Schlather (2014), Lantuéjoul (2002)) which are faster than the general Brown-Resnick processes [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]) which show similar properties to those of an extremal t process for high values α.

A major aspect that needs to be improved is the way to generate the partitions of the conditioning points. Indeed, the Gibbs sampler requires too many computational resources to provide better results. Some alternative methods have to be considered in order to avoid the calculations of the multiple integrals involved in the transition probabilities.

The application performed on a dataset of rainfall in the south of France showed that the model is relevant and that the algorithm performs well, provided that the spatial and temporal assumptions of stationarity are fulfilled. The margins modelling and transformation (equations ( 33) and ( 34)) intended to make the transformed process spatially stationary. However, the sole use of the coordinates (LAT+LON) is likely not sufficient to remove all the spatial heterogeneity of the process especially at the boundaries or where some other major climatic factors have to be taken into account. For this example, a small sub-region with different climate characteristics has been identified where the spatial coordinates do not seem to be sufficient to get satisfying results. Hence, the introduction of other covariates, such as altitude or information from climate or meteorological situations, could help to improve the simulations.

However, for simulations of fields of maxima precipitation in a future climate change context, no conditioning points (i.e., really observed at locations) are available. The only information stems from Global Climate Models (GCM) or Regional Climate Models (RCM) driven by various greenhouse gas emission scenarios. Those models provide many data but at low spatial resolution (from about 300 km down to about 10 km at best, so far). Hence, the best those models can give will be mean precipitation values over relatively large regions. Those data are clearly not adapted to feeding hydrological, ecological or, more generally, environmental impact models, which require high-resolution input climate data. Hence, the use of a max-stable process to simulate data for the future and provide those simulations to environmental impact models needs to be based on a change of spatial resolution. In other words, the max-stable process will not have "station" values available as conditioning points but only large-scale climate conditioning information. Therefore, one important need is to be able to condition max-stable processes by those large-scale climate data, which corresponds to performing a "statistical downscaling" to generate high-resolution fields of extreme precipitation. This downscaling in a spatial modeling framework will enable us to statistically generate fields and times series that will be useful for many impact studies trying to anticipate the potential consequences of extreme precipitation events in a climate change context.

Appendix A. Calculation of the density and the conditional density

Proof. According to equation (13), we have Λ (α) x (A) = ξ -l dξ = 2 (l-3)/2 [z T Σ -1

x z] -( l-1 2 ) m -(l-1)/α α Γ l -1 2 (A.4)

λ (α)
x (z) can be written:

λ (α) x (z) = π -k-1 2 |Σ x | -1/2 [z T Σ -1 x z] -( k+α 2 ) Γ k + α 2 Γ((α + 1)/2) -1 = π α-1 2 w k+α-2 |Σ x | -1/2 [z T Σ -1
x z] -( k+α 2 ) Γ((α + 1)/2) -1 (A. [z T Σ -1

x z] (k+α)/2 [(u, z)

I

  MS d (ŷ 1 , . . . , ŷn ; y 1 , . . . , y m ) = -y I || d (27) with ||.|| d the euclidean norm in R d , P m (d) the ensemble of indices of length d in {1, . . . , m}, y I = {y i |i ∈ I} and ŷ(j) I = {ŷ (j)

Figure 1 :

 1 Figure 1: Evolution of the CRPSS (in %) with the α with ordinary kriging as reference.

Figure 2 :

 2 Figure 2: Extremal coefficient: comparison of θ th (h) and θ est (h) for α ∈ {2, 4, 6}.

Figure 3 :

 3 Figure 3: Study area of SAFRAN data subset.

Figure 4 :

 4 Figure 4: Map of the crps (in mm) for 2007. For further illustrations, four grid cells have been arbitrarily selected and are indicated here (from 1 to 4).

Figure 5 :

 5 Figure 5: Histograms of the conditional simulations at 4 locations, with the real value (in red).

  (A) conditional simulation of extremal t model, with Gibbs sampler, (B) conditional simulation of extremal t , with random partitions, (C) analog model, (D) historical model, (E) unconditional simulation of extremal t model.

Figure 6 :

 6 Figure 6: Comparison of the evolution of the mean of the conditional simulations (in red), the historical model (in blue), analog model (in orange) and the test values (in black).

Figure 7 :

 7 Figure 7: Correlation between the simulated values and the test values.

1y

  {m -1/α α ξz∈A} g x (z)dzξ -(α+1) dξ with for z ∈ R k , g x (z) = (2π) -(k/2) |Σ x | -1/2 exp[T Σ -1 x y ξ -(k+α+1) dξ (A.3)By noting that for l ≥ 2,

Table

  Table 1 gives the CRPSS calculated by taking the random partition model as the reference.

			Real Part.		Gibbs sampler Part.		Random Part.
	α	2	4	6	2	4	6	2	4	6
	CRPS	0.158	0.132	0.118 0.237 0.182 0.173 0.250 0.192 0.175
	CRPSS 36.8% 31.2% 32.8% 5.1% 5.2% 1.2%	0	0	0

Table 2

 2 gives the CRPS, CRPSS, QS 95 , QSS 95 , MS 2 and MSS 2 for each method with the historical model (D) as reference.

	Method	CRPS CRPSS QS 95 QSS 95	MS 2	MSS 2
	(A) Cond. Sim. Gibbs	10.05	-4.1% 1.95 34.8% 27.55 31.3%
	(B) Cond. Sim. Random 10.96 -13.1% 2.02 32.4% 28.99 27.7%
	(C) Analog	10.53	-9.1%	5.68 -90.0% 18.41 54.1%
	(D) Historical	9.65	0%	2.99	0%	40.11	0%
	(E) Uncond. Sim.	12.71 -31.6% 3.33 -11.4% 48.07 -19.9%

Table 2 :

 2 CRPS, CRPSS, QS 95 , QSS 95 , MS 2 and MSS 2 for the methods (A), (B), (C), (D) (as reference) and (E).

  Table 3 shows the new scores.

	Method	CRPS CRPSS QS 95 QSS 95	MS 2	MSS 2
	(A) Cond. Sim. Gibbs	8.85	4.4% 1.78 37.2% 30.14 21.8%
	(B) Cond. Sim. Random 9.46	-2.1%	1.80 36.5% 30.79 20.0%
	(C) Analog	10.14	-9.4%	5.55 -95.7% 18.15 52.9%
	(D) Historical	9.27	0%	2.83	0%	38.51	0%
	(E) Uncond. Sim.	11.13 -20.3% 3.09 -9.1% 47.24 -22.7%

Table 3 :

 3 CRPS, CRPSS, QS 95 , QSS 95 , MS 2 and MSS

2 for the methods (A), (B), (C), (D)(as reference) and (E) with removal of the outliers area.

  T Σ -1 s,x (u, z)] (k+m+α)/2 (A.6) with w d = π d/2 /Γ(d/2 + 1) the volume of the unit ball in d dimensions.

	Appendix B. Estimated parameters
	Parameters Estimation	CI (95%)
	µ 0	-3.0565	[ -3.0942 ; -3.0187 ]
	µ 1	-0.2227	[ -0.2233 ; -0.2221 ]
	µ 2	0.1804	[ 0.1803 ; 0.1805 ]
	σ 0	-4.3431	[ -4.3732 ; -4.3129 ]
	σ 1	-0.3665	[ -0.3675 ; -0.3654 ]
	σ 2	0.2084	[ 0.2084 ; 0.2085 ]
	ξ	0.1145	[ 0.1128 ; 0.1163 ]
	α	2.2077	[ 2.0354 ; 2.3945 ]
	a	1.2262	[ 1.1304 ; 1.3301 ]

  Table B.1: Estimated parameters (and CI) with composite likelihood maximization.