Vivien Goepp
email: vivien.goepp@parisdescartes.fr

Jean-Christophe Thalabard

Grégory Nuel

Olivier Bouaziz

Regularized Bidimensional Estimation of the Hazard Rate

Keywords: Scientific communications Spline regression [Submitted] 2, Spatial regularization [In progress] Conferences Survival Analysis, Penalized Likelihood, Piecewise Constant Hazard, Age-Period-Cohort Analysis, Adaptive Ridge Procedure Spline Regression, B-splines, Penalized Likelihood, Adaptive Ridge, Bandlinear Systems, Changepoint Detection

This thesis deals with the development of regularized methods using penalized maximum likelihood estimation. More specifically, I use a sparsity-inducing iterative method called adaptive ridge. The latter is competitive compared to other approaches, namely in terms of ease of implementation and computational cost. My work consists in the application of this method to a wide range of problems: survival analysis, spline regression, and spatial segmentation. Applications in several problematics show that the adaptive ridge's good performance in selection, great ease of implementation and low computational cost can make it a good starting point in penalization-base variable selection.

In survival analysis, data are often collected by following a cohort, in which case the events are widely spread through time and the sample is suspected to present heterogeneity. I first focus on developing a method for the inference of the incidence, which allows to detect heterogeneity with respect to the date of birth (or cohort). A closely related problem is the study of the evolution of the inference as a joint function of the age, the date of birth (cohort), and the calendar time (period). Epidemiologists have long resorted to the age-period-cohort model or its submodels. The latter assume linear effects of each variable, which is deemed too simplistic to estimate potentially important features of the incidence. In this framework, I develop a model allowing for the joint estimation of two variables' effects and of their interaction.

Spline regression is known to be a competitive method for non-parametric regression. However the estimated spline depends highly on the initial choice of knots and choosing the best knots is a computationally hard problem. I propose an approach for the estimation of the best knots jointly with the spline function. By initiating a large number of knots and successively removing the least relevant ones, my method makes a slightly restrictive hypothesis to remove much of the computational burden.

In spatial statistics, the spatial domain is often divided into "units" and data are gathered at the unit level. The spatial effect is estimated on each unit and its representation is subject to the arbitrary of the unit division, which makes its interpretation difficult. This can be resolved by regularization, which reduces the variance and increases the interpretability. I present a model for segmentation of spatial data based on the adjacency structure of the units.

i RÉSUMÉ Cette thèse porte sur l'élaboration de méthodes régularisées utilisant l'estimation par maximum de vraisemblance pénalisée. Plus précisément, j'utilise une méthode parsimonieuse itérative, appelée adaptative ridge. Cette dernière est compétitive par rapport à d'autres approches, notamment en termes de facilité de mise en oeuvre et de temps de calcul. Mon travail consiste à appliquer cette méthode à un large éventail de problèmes : l'analyse de survie, la régression par splines et la segmentation spatiale. Ces applications dans différentes problématiques montrent que la bonne performance de l'adaptive ridge en sélection, sa grande facilité de mise en oeuvre et son faible coût de calcul peuvent en faire un bon point de départ dans les méthodes de sélection de variable par pénalisation.

En analyse de la survie, les données sont souvent recueillies en suivant une cohorte, auquel cas les événements sont largement répartis dans le temps et l'échantillon peut présenter une hétérogénéité. Je me concentre d'abord sur le développement d'une méthode d'estimation de l'incidence qui permet de détecter l'hétérogénéité par rapport à la date de naissance (ou cohorte). Un problème proche est l'étude de l'évolution de l'inférence en fonction de l'âge, de la date de naissance (cohort) et de la date calendaire (period). Les épidémiologistes ont longtemps eu recours au modèle age-periodcohort ou à ses sous-modèles. Ces dernières supposent des effets linéaires de chaque variable, ce qui est jugé trop simpliste pour estimer des caractéristiques potentiellement importantes de l'incidence. Dans ce cadre, j'élabore un modèle estimant cojointement l'effet de deux variables et de leur interaction.

La régression par splines est connue pour être une méthode performante de régression non paramétrique. Cependant, la spline estimée dépend fortement du choix initial des noeuds et le choix des meilleurs noeuds est un problème difficile en pratique. Je propose une approche permettant l'estimation des meilleurs noeuds conjointement avec la fonction spline. En initiant un grand nombre de noeuds et en supprimant successivement les moins pertinents, ma méthode fait une hypothèse légèrement restrictive pour diminuer grandement le temps de calcul.

En statistiques spatiales, le domaine spatial est souvent divisé en "unités" et les données sont recueillies au niveau des unités. L'effet spatial est estimé sur chaque unité et sa représentation est soumise à l'arbitraire de la division de l'unité, ce qui rend son interprétation difficile. Ceci peut être résolu par la régularisation, ce qui réduit la variance et augmente l'interprétabilité. Je présente un modèle de segmentation des données spatiales basé sur la structure d'adjacence des unités.

iii

Plan of this thesis

This thesis is organised as follows:

Chapter 1 introduces, compares, and discusses the main statistical approaches to model selection. A great emphasis is laid on penalized likelihood methods. First, I detail the most famous methods of penalized maximum likelihood estimate: the lasso, the elastic-net, and further refinements. I develop on the use of these penalties in the linear model. Secondly, I touch on two of the methods of model selection which are not based on penalization: best subset selection and stepwise selection. In the third section, I introduce the Majorize-Minimization (MM) optimization scheme, which, applied to penalized likelihoods, yields two important iterative penalized methods: the Local Linear Approximation (also termed adaptive lasso) and the Local Quadratic Approximation. I then introduce the iteratively defined penalized method used throughout this work: the adaptive ridge. Finally, I develop on the statistical methods which enforce an a priori structure on a parameter.

Chapter 2 deals with the application of the adaptive ridge to the context of hazard estimation in survival analysis. After an introduction on the topic, and an illustration of the method to a onedimensional case, we detail how the fused adaptive ridge allows for a new method of regularized estimation of the bi-dimensional hazard rate, with detection of breakpoints. Even though this methods applies to a wide array of problems, we illustrate it through the angle of age-period-cohort analysis.

Chapter 3 deals with the same problem than Chapter 1, but through the context age-period-cohort analysis. We first introduce the topic of age-period-cohort analysis, as well as the use and drawbacks of the age-period-cohort model. Then, we develop on the "Age-Cohort-Interaction " model, which builds on the works of the previous chapter. This model can be viewed as a generalization of the age-period-cohort analysis, which does not suffer from its defects, at the added cost of computation time.

Chapter 4 deals with a different application of the adaptive ridge. In this part, we apply this method to the problem of finding the best knots to support a regression spline. This problem has long been deemed computationally intractable. We show that provided some simplifying assumptions, our new spline regression method can select the best knots as well as the regression spline in a fast fashion. Before developing on our method, we introduce the topic of spline regression and present the main tools and issues in this topic.

Chapter 5 deals with the applications of the adaptive ridge to regularization of spatially-correlated data. When the statistical problem has a spatial structure that is given in already chosen zones, the problem of regularization becomes challenging. We introduce a method for regularization along a graph, with an application to inference for medical data.

The works of Chapters 2 and 4 has lead to two preprints currently under revision. These two chapters consist of these preprints, reproduced as is and preceded by introductory talks on their respective matters. These two papers are given in Sections 2.2 ad 4.2, with respective supplementary materials given in Sections 2.3 and 4.3 respectively. These papers define their own notations, which v vi are for the most part consistent with the rest of the manuscript. Except for these sections, the present document forms a consistent manuscript, even if the Chapters were divided so as to be able to be read separately.

Notations and Definitions

Definitions

• The L q norm of a vector u is noted u q . The L 2 norm is also abbreviated u .

• E [X ] and V [X ] denote the respective expectation and variance of the random variable X .

• Tr denotes the trace.

• → d denotes the convergence in distribution.

• A T denotes the transpose of the matrix A.

• min + {u} denotes the minimum taken only over the positive values of the vector u • (x) + denotes the positive part of x, that is max (x, 0).

• diag{u i } i is the diagonal matrix whose non-zero entries are u 1 , . . . , u n .

• N (., .) denotes the (possibly multivariate) normal distribution. The two arguments are the expectancy and variance.

• #A denotes the cardinal of the set A.

• When f is a function and x an element from its domain, f (x -) and f x + denote the limits lim t →x f (t ) in the cases t < x and t > x, respectively.

Notations

• Vectors and matrices are noted in bold. When necessary, vectors are identified with column matrices.

• n is the sample size

• i ∈ {1, . . . , n} is the index of the individuals

• p is the number of covariates

• j ∈ 1, . . . , p is the index of the covariates • In the context of penalized likelihood methods, λ ∈ R is the penalty constant; in the context of survival analysis λ : R → R + is the hazard rate.

• β is the parameter to be estimated

• I refers to the identity matrix, whose dimension depend on the context.

• = -log L is the negative log-likelihood and L is the likelihood.

• I use the symbol " " when the equality serves as a definition. ix x

Abbreviations

• MM: Majorize-Minimization optimization

Conflicts in notation between chapters

We have tried to use coherent and non-conflicting notations for the mathematical objects defined in this thesis. However, for the sake of consistency with the conventions of the field, we made the choice to keep conventional notations for known quantities. The instantaneous hazard rate for instance, is noted λ(t ) (as a function of t ) as is standard in survival analysis and in the study of stochastic processes. In other parts of the manuscript, we also used the variable λ to denote the penalty constant in penalized maximum likelihood methods. These notational conflicts have been kept to ease the understanding of the manuscript. They occur between different chapters but not inside each chapter. We stress that the potential uncertainty is removed when the context is taken into consideration.

Introduction

This thesis deals with the application of penalized methods to different statistical problems. These works all have in common the use of a penalized maximum likelihood method-called the "adaptive ridge" -which performs model selection. The latter comes from a long line of historically, practically, and theoretically important regularization methods, which date back to the very beginning of computational statistics. These methods are at the intersection of different fields: statistics, optimization, and computer science.

It is therefore necessary to introduce the most important of these methods in order to (i) highlight the principle of penalized estimation methods, (ii) put the adaptive ridge in the context of the other model selection methods, and (iii) compare the adaptive ridge with the competing methods. This introduction serves this purpose.

Variable selection methods have been developed since the 1970's, with the use of the stepwise selection and the ridge regression in the context of linear regression. But their use and fame have exploded only in the 1990's and early 2000's, with the development of penalized regression method, and in the first place of the lasso, introduced in the field of statistics by Tibshirani (1996) and in signal processing by [START_REF] Shaobing | Atomic Decomposition by Basis Pursuit[END_REF] under the name "basis pursuit". These methods have been introduced for linear regression, but their principle applies to any statistical model whose likelihood is computable and practical to maximize. Penalized regression methods consist in adding a term in the negative log-likelihood (NLL) to minimize. This term enforces the estimate to be close to an a priori shape or distribution. Certain penalty terms have been found to induce the desired properties of the estimate: the resulting estimate is both better quantitatively, i.e. it has better estimation performance, and qualitatively, i.e. it infers models which are more relevant and easier to interpret, than standard estimates. These methods sparked a revolution in the field of computational statistics: the added penalty term increases the complexity of the computation by a low margin and the benefits are huge in many practical applications.

A whole array of these penalized methods have then been developed and improved upon, with penalties ever more refined and adapted to the problem at hand. However, the initial application of penalized estimation with variable selection is the linear model, and its extension to non-Gaussian errors, the generalized linear model. It seemed that in this case, the goal to find easy-and fast-tocompute methods that are performant both theoretically and practically, has been met. Two other major fields of applications of penalized methods have emerged around this period: high throughput data and wavelet analysis. The former became widely popular in the last decade due to the onset of next-generation sequencing technologies. Since many genes are studied at once, we are in a case where p n, usually by a factor ∼ 100. The problem posed by this high-dimensional setting was seldom met in usual regression settings, and led to the development of penalized regression methods for what is called high-dimensional statistics. The latter comes from the development of wavelets bases and their applications to all fields of signal processing. Wavelets are families of functions that form an orthogonal family of L 2 (R), are located in both the time and frequency domains, and are all the scaled and shifted versions of one another. These properties make it the tool of choice for representing signals (i.e. audio signals, images, videos, etc.) efficiently. Indeed, the wavelet bases allow for sparse representation of signals, with major applications to denoising, compression, and compressed sensing. In these problems, the signal sparsity is enforced using a penalized approach. This has sparked the development of penalized likelihood approaches that are specifically fit to the topology of the signal. Note that the leap forward made in that domain around the year 2010 was also enabled by the progress in convex optimization [START_REF] Boyd | Convex Optimization[END_REF]. We refer to [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF] for more details on the wavelet analysis and to [START_REF] Bach | Optimization with Sparsity-Inducing Penalties[END_REF] for more details on the application of penalized methods to wavelet analysis.

This introduction is organized as follows. The first part draws a panorama of the main penalized methods for model selection, with comparisons, explanations, and insights. The second part deals with subset selection, the main model selection method that does not used a penalized approach. The third part puts the iterative penalized methods in the theoretical framework of Majorize-Minimization (MM) optimization and compares the main iterative penalized methods amongst which is the adaptive ridge. The last section presents modifications of the penalty terms which enforces the estimate to have a certain sparsity structure, which is more general than being sparse. We arrive at the fused adaptive ridge, which is used in Chapters 2 and 3.

Regularized estimation

Linear regression. Consider the linear regression setting y = X β + ε, (1.1) where y is the n × 1-sized response variable, X = x i , j is the n × p-sized design matrix, β is the p × 1 parameter vector of linear effects, and ε is the p × 1 vector of random errors. The columns of X are called the covariates and are noted x j . We consider the case of deterministic design, that is, X is deterministic.

Linear regression is sometimes written

ŷ = β0 + X 1 β1 + • • • + X p βp ,
which includes the estimation of an intercept β 0 , adding a row of 1s in the design matrix X . Since the maximum likelihood estimate of β 0 is ȳ = 1 n n i =1 y i , we will instead consider the model (1.1) and always assume that the response variable is centered: ȳ = 0.

It is also necessary that the covariates be scaled, both for the numerical stability of the computations and to be able to compare the effects of each covariate. For instance if the unit of x j is changed from meters to millimeters, βj is multiplied by 1000. Hence the effect of that covariate will be artificially inflated and it will always be selected by a model selection method. Thus, the covariates have to be scaled for their effects to be comparable. Note that there is an equivalence between fitting the parameters on the unscaled covariates and fitting them on the scaled covariates and applying the corresponding scaling to each parameter. Therefore, without loss of generality, we will assume that n i =1

x i j = 0 and

n i =1 x 2 i j = 1, 1 ≤ j ≤ p.
The ordinary least squares (OLS) is βols arg min β y -X β 2 2 = X T X -1 X T y.

(1.2)

If not stated otherwise, we assume the two classical following assumptions to hold.

Assumption 1. y i = x i β * +ε i , where ε i are independent and identically distributed random variables of mean 0 and variance σ 2 , and β * is the real value of β that we estimate. Note that unless stated otherwise, the errors are not assumed to be normally distributed.

Assumption 2. The matrix 1 n X T X → C , where C is positive definite.

Penalized likelihood. In the forthcoming sections we will discuss penalized regression methods, in which the estimator is defined as the minimizer of the least squares residuals with an added penalty term. In the linear model, when the residuals are normally distributed, the least squares residuals y -X β 2 2 is (proportional to) the negative log-likelihood β of the model (this condition is in fact an equivalence: the negative log-likelihood writes as a sum of squares only when ε ∼ N (0, σ 2 I )). The regularization -and variable selection -methods present in this introduction have been developed in the setting of linear regression, with the case of normal residuals as an important specific case. In all generality, they can apply to any parametric model that we want to regularize. In this introduction, the notation β refers to the negative log-likelihood in the linear model or in any unspecified parametric model. In applications and illustrations however, we will mostly use the linear model to illustrate the effect of penalized estimation.

We give two examples of such models, which are generalizations of the linear model:

• Robust regression, in which the negative log-likelihood writes

β = n i =1 ρ(y i -x i β) (1.3)
where ρ is a function giving little importance to large values of the residuals. For instance, when ρ(x) = |x|, the model becomes the least absolute deviation, a robust alternative to the linear regression. When ρ is a ρ-function (as defined in [START_REF] Peter | Robust estimation of a location parameter[END_REF], the model is a robust version of the linear model. See [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF] for a thorough explanation of robust regression.

• The generalized linear model the response variable has a distribution in the exponential family, where the density of y i (in its canonical form, see [START_REF] Mccullagh | Generalized Linear Models[END_REF], Section 2.2.2) writes

f (y) = exp θ * i T (y) -A(θ * i )) + B (y) ,
where A and B are known deterministic functions and θ * i is the true unknown parameter to be estimated. In the generalized linear model, we estimate the density's parameter vector θ * by a function of the linear effect x i β, and we have E [Y ] = g -1 X β where g is a link function. Usually we take the canonical link function given by the distribution of y: g -1 = A. The negative log-likelihood (nll) writes

β = - n i =1 log f g x i β , y i , (1.4)
where f is the distribution of y. For example, if y is Poisson distributed, the canonical link function is the log function and the nll writes

- n i =1
y i x i β -exp(x i β) -log(y i !) .

Ridge regression

The OLS estimate is unbiased and has minimal variance amongst all unbiased estimates. However, when the covariates x j are highly correlated, X T X becomes close to singular, i.e. when its largest eigenvalue gets close to zero, the variance of the OLS diverges to infinity. In this case, a biased modification of the OLS estimate is necessary to control the variance of the estimate. The ridge regression, defined below, is the most famous of such modifications -owing first to its simplicity. It was introduced in the context of linear regression and we use this specific case to illustrate its principle and effect, but it extends naturally to any other model. It is obtained by adding a regularizing term to the nll, in the way defined below. [START_REF] Arthur | Ridge Regression: Biased Estimation for Nonorthogonal Problems[END_REF] introduced the ridge estimator. It is defined as the solution to

min β y -X β 2 2 + λ β 2 2 , (1.5)
where λ > 0 is a trade-off hyper-parameter to be chosen. Problem (1.5) has the explicit solution βridge = X T X + λI -1 X T y.

(1.6) Equation 1.6 can be seen as a relaxation of the matrix X T X . When the covariates are too correlated, the columns of X are close to being linearly dependent, and X T X tends to be singular. Equivalently, this means that the lowest eigenvalue of X T X tends to zero. If λ 1 is the lowest eigenvalue of X T X , the lowest eigenvalue of X T X + λI is λ 1 + λ. Equation 1.6 can be seen as a modification of the OLS that ensures the design matrix does not become singular. for t > 0. The two problems have been shown (Boyd and Vandenberghe, 2004, B.1) to have strong duality, that is, for every λ > 0, there exists a t > 0 such that the solution of (1.7) is also the solution of (1.5), and conversely. The decreasing one-to-one relation between λ and t depends on the data. The equivalent constrained optimization problem helps understand the effect of the penalty included in ridge regression. Figure 1.1 illustrates the shrinkage of β due to the projection onto the L 2 norm ball. We recall that in a normed vector space, for a chosen norm, the ball of radius r > 0 and of center c is the set of elements at a distance r or less to c, i.e. the set {x| x -c ≤ r }. In this manuscript, we will only consider centered balls, that is, balls of center 0. Balls of radius 1 are called unit balls.

Ridge regression in orthogonal design. Consider the case of orthogonal design, that is, when the design matrix is orthogonal: X T X = I . Notice that in this case, the quadratic terms y -X β 2 = y T y -2y T X β + β T β and X T y -β 2 = y T X X T y -2y T X β + β T β differ only by a constant. Then the linear regression boils down to minimizing z -β 2 , where z = X T y is the transformed data. The ridge problem (1.5) rewrites min

β z -β 2 + λ β .
This particular case is important because it is the case in which the covariates are uncorrelated, and the estimation is done component by component. Also, the last equation can be interpreted as finding the closest vector to the data z with a small L 2 norm. This result can also be seen directly through the simplifications in (1.2) and (1.6).

In this case, then, we have βols = X T y and βridge = βols 1 + λ .

Consequently, in orthogonal design, the ridge estimate simply shrinks every coordinate of the OLS towards 0 by the same factor.

Minimization of the MSE. Define the mean square error as

MSE ols = E βols -β * 2
Define λ 1 ≥ • • • ≥ λ p > 0 as the eigenvalues of X T X . Recall that the OLS is unbiased and from the equality βolsβ * = (X T X ) -1 X T ε, we obtain

MSE ols = σ 2 Tr{ X T X -1 } = σ 2 p i =1 1 λ i . (1.8)
This MSE takes problematically large values when λ 1 is small. The ridge estimate aims at reducing the MSE to a value lower than (1.8). The ridge estimate reduces the variance of the estimate, at the price of an non-zero bias. More precisely, the MSE of the ridge estimate decomposes as follows. Consider the singular value decomposition of X :

X = U DV T ,
(1.9)

where U and V are n ×p and p ×p orthogonal matrices respectively, such that U T U = V T V = I p , and D is a diagonal matrix whose entries d 1 ≥ • • • ≥ d p ≥ 0 are the singular values of X . From the equality X T X = V D 2 V T , we can see that the d 2 j s are the eigenvalues of X T X . By plugging (1.9) into (1.6), we get

MSE(λ) = E[ β -β * 2 ] = λ 2 β * T (X T X ) -2 β * + E εX (X T X + λI ) -2 X T ε = λ 2 β * T V D -4 V T β * + σ 2 p j =1 d 2 j (d 2 j + λ) 2 .
(1.10)

In the last equation, the first term is the bias of the ridge estimate, the second term is its variance. Since d MSE d λ (λ)

λ=0 + < 0,
there exists a value of λ such that the ridge estimate has a smaller MSE than the OLS.

Lasso estimation

Why select variables. In many practical situations, the statistician does not know what covariates has an influence on the response variable or does not want to make such an assumption a priori. A solution is to add all possible variables and find a regression procedure which simultaneously selects the relevant covariates and estimates their effect. This is the framework of regression with variable selection.

Lasso as a relaxation of the L 0 norm. Define the L 0 "norm" x 0 = #{i |x i = 0} as the number of nonzero elements of a vector. This quantity is not a norm, because for α = 0, αx 0 = |α| x . We still refer to it as a norm since it can be seen as the limit of the L q norm when q → 0. The L 0 penalized likelihood approach offers to minimize y -X β 2 2 + λ β 0 .

(1.11)

The minimizer of (1.11) represents a trade-off : it has both to be close to the minimum of the squared residuals and to have few non-zero coordinates. Consequently this estimate will only select the most relevant covariates in the model and estimate their effect. Unfortunately, to minimize the L 0 norm is NP-hard [START_REF] Natarajan | Sparse Approximate Solutions to Linear Systems[END_REF], which implies that there is no known algorithm that can solve an L 0 norm-constrained problem in a polynomial time complexity with respect to p. Consequently, this problem has to be solved by trying all the possibilities, which has exponential time complexity in p. This is unfeasible in practice for p 50 or larger. L q Norms q = 0.5 q = 1 q = 2 Figure 1.2: Illustration of L q unit balls for different values of q. Tibshirani (1996) introduced the lasso estimator defined by βlasso arg min

β y -X β 2 2 + λ β 1 (1.12)
where λ > 0 is a hyper-parameter which tunes the trade-off between goodness of fit and regularization. The choice of λ is not touched on in this chapter. Lasso uses the L 1 norm which can be seen as a relaxation of the L 0 norm. The L 1 norm has been used for many decades to recover sparse vectors. In signal processing, a signal y is expressed as a combination of functions forming a dictionary φ. Since φ is over-complete, it has more columns than rows, and the coefficient vector x of the decomposition of the signal in the dictionary is not uniquely determined by the problem y = φx. However one assumes that y has a sparse representation in the dictionary. By enforcing a sparsity constraint over x, a compressed representation of the signal y is obtained. Using the L 1 penalized problem to recover a signal is called basis pursuit. This method was developed by [START_REF] Shaobing | Atomic Decomposition by Basis Pursuit[END_REF]; it was further analyzed by [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF] and [START_REF] Elad | A Generalized Uncertainty Principle and Sparse Representation in Pairs of RˆN Bases[END_REF].

Using the Lagrange multipliers, the lasso estimator can be defined as the projection of the OLS onto the L 1 ball of radius t : βlasso = arg min

β y -X β 2 2 s.t. β 1 ≤ t , (1.13)
where t > 0 has a decreasing one-to-one relation with λ. Figure 1.2, which represents L q unit balls for different values of q, illustrates why the L 1 norm induces sparsity. When q > 1, the L q ball is smooth and the projection onto the ball does not set any coordinate to zero. When q ≤ 1, the L q ball has singularities at the axes and the projection onto the ball can set some coordinates to zero. Moreover for q ≥ 1, L q norms are convex, which makes computational minimization way easier. The L 1 norm naturally comes out as the only easily minimizable penalty which enables variable selection, amongst all L q norm penalties.

L 1 norm and sparsity. As previously explained, the projection onto an L 1 ball can set some coordinates to zero. This phenomenon is illustrated in Figure 1.3 with p = 2. The grey areas represent the half-cones from which a projection on the L 1 ball sets one coordinate to zero. Consequently, when βols has a coordinate of small value, its projection maps it to zero. The lasso has the effect of selecting the covariates with a non-negligible effect on y. As illustrated in the figure, as λ increases the 3: Illustration of the lasso estimation under orthogonal design. Projection onto two L 1 norm balls of radius t . As t decreases, the OLS will most probably end up in a greyed area and the lasso estimate will become sparse.

radius of the ball decreases and the grey area takes up a larger proportion of the parameter space. Consequently, the lasso will tend to set many covariates to zero as the penalty increases.

Lasso in orthogonal design: the soft-thresholding operator. In orthogonal design (see previous Section), the lasso problem reduces to the component-wise problem βlasso j = arg min β j f (β j ) where f (β j ) =1 2 ( βols j -β j ) 2 + λ|β j | and where βols = X T y is the OLS estimate in orthogonal design. Since f is the sum of two functions, who minima are in βols j and 0 respectively, its minimum is located either at 0 or at β ols j . By considering the two cases |β ols j | ≤ λ and |β ols j | ≤ λ separately, we get that the lasso solution has the explicit solution given by βlasso = S( βols ) sgn( βols ) βolsλ + , (1.14) where the function S is called the soft-thresholding function, and where the latter formula is to be considered component by component. This function is equal to zero over the interval [-λ, λ], which means that in orthogonal design, the lasso sets small value of the data to zero. We refer to Section 1.1.7 for a graphical representation of the soft-thresholding operator and a discussion of its properties.

Computation of the lasso solution. Nevertheless, since the L 1 norm is not differentiable at the minimum, the computation of the lasso estimation is not straightforward. Two algorithms for solving the lasso have gained popularity since its development.

The Least Angle Regression (LARS) algorithm has been proposed by [START_REF] Efron | Least Angle Regression[END_REF] and offers to compute the whole regularization path 1 . It initiates with no covariates in the model and sequentially adds the covariate that has the maximum correlation with the residuals. The estimate is increased linearly in the direction of equal correlation between all active covariates (hence the names least angle). Once a new covariate has the same correlation with the residuals as the current covariate, it is added to the set of active covariates and becomes the current covariate. This method is explained in more details in the Appendix.

Chapter 1. Introduction

Another approach is the coordinate descent, a simple optimization method that is particularly fit to solving the lasso. This method is explained here.

The coordinate descent was initially proposed for the lasso by Fu (1998) under the name of shooting algorithm. It was generalized to other related penalties by [START_REF] Friedman | Pathwise coordinate optimization[END_REF]. The method solves the lasso problem for a fixed λ.

The algorithm starts with an initial guess β(0) and minimizes the penalized residuals with respect to each component in a fixed cyclical order until convergence. For 1 ≤ k ≤ p, we denote by β k the k-th component of β. The minimization of (1.12) with respect to β k writes min

β k n i =1 ỹ(k) i -x i ,k β k 2 + λ|β k | (1.15)
where ỹ(k) i y ij =k x i j β j . The solution to this problem is simply the soft-thresholding operator applied to the modified output:

β * k = n i =1 ỹ(k) i x i k -λ + sgn n i =1 ỹ(k) i x i k (1.16)
Hence the coordinate descent algorithm for the lasso is simple, easy to implement and fast to compute. However, contrarily to LARS, it does not compute the solution path for all penalties; we fix a grid of penalties and compute the coordinate descent for each penalty. To speed up the computation, the initial value in the coordinate descent is taken as the minimum found for the previous penalty. This computational trick is often referred to as warm start. The coordinate descent procedure with one penalty is given in Algorithm 1.

Algorithm 1 Coordinate descent with one penalty 1: function COORDINATE-DESCENT(y, X , λ)

2:

β ← X T X -1 X y 3:

k ← 1 4:

while not converge do 5: ỹ(k) i ← y ij =k x i j βj 6:

βk

← S n i =1 ỹ(k) i x i k , λ 7: k ← k mod p + 1 8:
end while 9: end function In the lasso regression, the function to minimize is the sum of the residuals, which is convex differentiable, and the L 1 norm, which writes as the sum of univariate convex functions. The rational behind the use of coordinate descent to solve the lasso problem is the following property.

Property 1. Consider a function f : R p → R which writes f β = g β + p j =1 h j β j where g is convex differentiable and the h i , 1 ≤ i ≤ p are convex. Assume that β * is a minimum of f along all coordinates, that is: ∀i ∈ 1, . . . , p , ∀a ∈ R, f β * + ae i ≥ f β * where e j is the j t h basis vector. Then β * is a global minimum.

Proof. Let β * be a minimum of f along all coordinates. For any β:

f β -f β * = g β -g β * + p j =1 h j β j -h j β * j ≥ p j =1 ∂g ∂β j β * β j -β * j + h j β j -h j β * j
The inequality comes from the convexity and differentiability of g . Since f is maximal along all coordinates, every term in the sum is positive. This completes the proof.

This property gives an incentive to use the coordinate descent when the function to minimize has the aforementioned form. This is the case for the lasso penalized regression as well as for any lasso problem where the likelihood of the model is log-concave. This property alone does not guarantee the coordinate descent converges to the global optimum. When g is strictly convex, the coordinate descent has been proven to converge [START_REF] Tseng | Coordinate Ascent for Maximizing Nondifferentiable Concave Functions[END_REF].

Theoretical properties of the lasso. In this section we consider the asymptotic performances of the lasso. Since this selection method performs selection and estimation, we need to define asymptotic properties for both. The definitions and properties present in this section are found in [START_REF] Zou | The Adaptive Lasso and Its Oracle Properties[END_REF].

Let A n = { j : βj = 0} be the selected covariates and A * = { j : β * j = 0} the true non-zero covariates. The cardinal of A * , i.e. the number of true non-zero covariates, is noted p 0 .

Definition 1. An estimate β is said to be consistent in selection if lim n P (A n = A * ) = 1.
For simplicity, we assume without loss of generality that A * = {1, . . . , p 0 }. Let

C = C 11 C 12 C 21 C 22 , (1.17)
where C 11 is a p 0 × p 0 matrix. Let βA * (resp. β * A * ) be the first p 0 elements of the estimate β (resp. β * ). If an estimate is consistent in selection, the support A * will be known with large probability for n large enough. The question arises of the estimation quality of βA * .

Definition 2. The estimate

β is said to be v n -consistent in estimation if v n ( βA * -β * A * ) → d N (0,C 11 ).
Note that C 11 is the variance matrix of the estimate knowing the right model, i.e. knowing A * . An estimate which is both consistent in selection and n-consistent in estimation is said to have the oracle properties: it performs as well as the maximum likelihood estimate would if we knew the true covariate indices. Consequently, the oracle properties is the best we can ask for in a variable selection method, in terms of prediction accuracy. This means that when an estimate has the oracle properties, it is always more advantageous to use it instead of the OLS.

These considerations have to be moderated by the fact that (i) in practice n is not always very large (ii) the conditions about λ are related to its asymptotic speed. [START_REF] Zou | The Adaptive Lasso and Its Oracle Properties[END_REF] proved the following property.

Property 2. The lasso estimate is not necessarily consistent in selection. More precisely, if

λ n = O n , lim sup n P (A n = A ) ≤ c < 1,
where the constant c depends on the true model.

A small modification of the lasso has been proven to enjoy the oracle properties. This estimate is presented hereafter.

The adaptive lasso. [START_REF] Zou | The Adaptive Lasso and Its Oracle Properties[END_REF] introduced a slight modification of the lasso with a weighted L 1 norm:

βal arg min β y -X β 2 + λ n p j =1 w j |β j |.
(1.18)

The weights are defined as w j = 1/| βj | γ , where β is any consistent estimate of β (for instance the OLS) and where γ > 0 is a hyper-parameter to be chosen. The weights have the effect to penalize more coordinates which have a small estimate. Due to this rescaling, the adaptive lasso is proven to have the oracle properties when λ n / n → 0 and λ n n (γ-1)/2 → ∞. Moreover, the adaptive lasso estimate has the same computational cost as the lasso.

Elastic-net

The elastic-net penalty was introduced by [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]:

λ α β 2 2 + (1 -α) β 1 (1.19)
where α ∈ [0, 1] is an hyper-parameter to be fixed. The elastic-net penalty defines a norm that is a compromise between the L 1 and the L 2 norm: as α varies, the elastic-net estimate varies continuously between the lasso and the ridge estimates. For α = 0 it becomes the lasso penalty and for α = 1 it becomes the ridge penalty. Contrarily to the L q norm (1 ≤ q ≤ 2), the elastic-net penalty is non-differentiable, as illustrated in Figure 1.4. Consequently, the elastic-net always performs selection (for α = 1) like the lasso and shrinkage like the ridge (which reduces the MSE). Comparatively, the L q norm (1 < q < 2) does not perform model selection and is of little interest. As discussed hereafter, the elastic-net enjoys the desirable properties from both estimates.

The thresholding function of the elastic-net under orthogonal design is

βen i = (| βols i | -λ 1 /2) + 1 + λ 2 sgn{ βols i }. (1.20)
This thresholding function has the same behavior as the soft thresholding for small values of βols i and the same behavior as the ridge thresholding for large values of βols i . An illustration is given in Figure 1.7 with λ 1 = 2 and λ 2 = 1.

The elastic net is also more stable in selection than the lasso. More precisely, if several covariates are very correlated, the lasso will arbitrarily select one covariate and discard the others. This is proven [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] not to be the case for the elastic-net: the difference β i -β j is bounded by 1-corr{x i , x j }, so that highly correlated covariate effects have close to equal estimates. This property is known as the grouping effect: the elastic net selects important covariates in group.

Moreover, in the p n scenario, lasso does not perform well in selection because it cannot select more than n covariates. The elastic-net overcomes this obstacles and, even with small α, can select any number of covariates.

Finally, the elastic-net estimate writes also as the minimizer of

β T X T X + λ 2 I 1 + λ 2 β -2X T y + λ 1 β 1 , (1.21)
hence it is computed as a lasso estimate with a modified design matrix.

The LARS algorithm can be adapted to estimate the elastic-net. Recall that LARS computes the whole regularization path. Thus the LARS algorithm is called on a grid of values of λ 2 and the selected value is the one which minimizes the (ten-fold) cross-validation error. The computational cost of this procedure is K times more than that of a lasso fit, where K is the length of the grid of λ 2 .

Bridge regression

Bridge regression is based on the L q penalty: βbridge arg min

β β + λ p j =1 β j γ (1.22)
where γ > 0 is a parameter to be fixed. Bridge regression was first introduced by Frank and [START_REF] Ildiko | A Statistical View of Some Chemometrics Regression Tools[END_REF] as a generalization of the ridge (γ = 2) and the lasso (γ = 1) and was later studied by Fu (1998), who compared it with the lasso. For 0 < q ≤ 1 the L q penalty is singular around zero and the bridge regression performs variable selection. This property is illustrated in Figure 1.2 which represents unit balls of L q norms for different values of q. Figure 1.4: Illustration of why the elastic net induces sparsity and not the L q norm (1 < q < 2). Unit radius balls of the elastic-net norm (α = 0.2) and the L 1.2 norm.

Asymptotic properties of the bridge regression has been studied by [START_REF] Knight | Asymptotics for lasso-type estimators[END_REF] in both cases γ ≤ 1 and γ > 1.

Fu (1998) proposed an algorithm based on the Newton-Raphson procedure to compute the bridge estimate when γ > 1. However, when 0 < γ < 1, the penalty is not convex and its minimization is computationally challenging.

The bridge penalty is part of a general family of sparsity inducing penalties called non-concave penalties These are defined in the next section. The oracle properties of the bridge estimate and the algorithms available to compute it are given in this section.

Berhu

The Berhu penalty was introduced by [START_REF] Owen | A robust hybrid of lasso and ridge regression[END_REF]:

p λ (|θ|) = |θ| |θ| ≤ λ θ 2 +λ 2 2λ |θ| > λ (1.23)
It is equal to the absolute value around zero and is quadratic away from zero. Consequently it is sometimes referred to as inverse Huber function [START_REF] Peter | Robust estimation of a location parameter[END_REF]. The presence of the singularity of the penalty at zero makes it a variable selection method. The Berhu penalty replaces the absolute value penalty by a quadratic penalty for large values in order to benefit from the properties of the ridge regression: the grouping effect, and an increased accuracy in estimation. The Berhu penalty can be seen as mixing the L 1 and L 2 penalties; in this regards it is very similar to the elastic-net.

Penalized likelihood methods as a Bayesian prior

Let the variable M denote the model. Recall that in the Bayesian framework, the data are generated from a model which is also assumed to be random. The distribution of the model is an a priori information that we must provide. The goal of Bayesian inference is to infer the posterior likelihood P (M |data). From the Bayes formula: where π(M ) is the prior distribution of the model, which includes a priori information about the model distribution.

P (M |data) = P (data|M ) π(M ) P (data) , ( 1
Here, the symbol for probability must be understood in the sense of density. The prior on the data does not intervene in the minimization of the integrated negative log-likelihood. The first probability in the right-hand side of (1.24) is the likelihood (of the data). In this case, the model M is the parameter β from which the data is generated. Then, Bayesian inference is made by minimizing:

-2 log P β|data = 2 β -2 log π(β), (1.25) where we recall that (β) is the negative log-likelihood (of the data conditionally on the model). Equation 1.25 is very similar to penalized maximum likelihood estimation. Notice that Equations (1.25) and (1.5) are the same on the condition that

-2 log π β = λ β 2 .
Consequently, the ridge penalty can be seen as a Bayesian inference with a normal prior distribution on the parameter: β ∼ N 0, λ 2 . The hyper-parameter λ is the standard deviation of the centered normal prior on the parameter. More generally, the L q norm penalty corresponds to the prior

π(β) = exp -λ β q q 2 .
The lasso regression corresponds to a prior on the parameter with a Laplace distribution. These priors are represented in Figure 1.5 in the univariate case. Likewise, the elastic-net corresponds to the prior

π β ∝ exp -λ α β 2 2 + (1 -α) β 1 2 ,
which is also represented in Figure 1.5.

Non-concave penalties

Non-concave penalized estimates are defined as the minimizers of

β + p j =1 p λ β j (1.26)
where p is a non-concave function, or non-concave penalty, in the following sense.

Definition 3 (Non-concave function).

A real-valued function p λ (|θ|) defined on R is said to be nonconcave if it is (i) even on R with p λ (0) = 0 (ii) non-decreasing and concave on (0, ∞) (iii) differentiable over R * .

The differentiability condition can be relaxed to piecewise differentiability without without loss of generality. Since it is not a restrictive assumption in practice, we will assume differentiability for the sake of the presentation. We emphasize the fact that in this work, the term "non-concave" is not used to say that a function is not concave.

In the cases where p λ (|β j |) is proportional to λ, we will denote p λ (|.|) by λp(|.|) and drop the index λ. One may also want to penalize some variables more than others, because of some a priori knowledge about their importance. This is done by setting a different penalty for each coefficient. For the sake of the presentation, we will assume that all the components of β are penalized the same.

The estimates defined by (1.26), where p λ is a non-concave penalty, will be called non-concave estimates. Many important sparse estimates write as non-concave estimates, which makes (1.26) a good framework to analyze and compare sparsity-inducing penalized likelihood methods. Using the definition of non-concave penalty, we can (i) give conditions on the penalty for the corresponding estimate to be sparse, (ii) give theoretical results shared by non-concave sparse estimates, (iii) compare the behavior and performance of sparsity-inducing estimates by comparing their corresponding penalty functions.

The lasso and the bridge (0 < q < 1) are both non-concave estimates ; this follows directly from (1.12) and (1.22) respectively. In the following sections we define two important non-concave penalties: the hard-thresholding and the smoothly clipped absolute deviation (SCAD). on (0, ∞), they verify p (0 + ) > 0. Since they are even functions, we also have p (0 -) < 0. Thus by definition, they are not differentiable at 0. This property is proven [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF] to be a sufficient condition for a penalty to perform variable selection. Thus all non-concave penalties such that p λ (0 + ) > 0 perform variable selection. Figure 1.6 illustrates the shapes of several non-concave penalties. The sharper the penalty, the better it approximates the L 0 norm, which is the desired penalty. However, sharp penalties are harder to minimize, which makes the estimation procedure more computationally intensive. First, penalties like the SCAD and hard-thresholding, which are constant away from zero, are hard to minimize. The lasso is an outstanding non-concave penalty: it is minimally non-concave, since it is linear on (0, ∞), which is the least concave function. Moreover, the derivative of the L 1 penalty is constant on (-∞, 0) and (0, ∞) and thus its minimization cannot be performed using classic optimization methods. Finally, note that empirically, the sharper the spike around zero, the "stronger" the estimates will be enforced to be sparse.

Figure 1.7 represents the thresholding functions of some non-concave penalties under the case of orthogonal design: βj = f ( βols j ) = f (X T y). Asymptotically the variables with a non-zero effect on y will have an OLS estimate far away from zero and the variables with no effect on y will be around zero. Consequently, the penalty (i) Performs variable selection only if the thresholding function is identically equal to zero around a neighborhood of zero;

(ii) Estimates the non-zero effect without bias only if the thresholding function is close to the identity function (dotted lines on the figure) for large values of the data.

(iii) Is stable in selection if the thresholding function is continuous. (By stable we mean that the selected model is not highly sensitive to small variations in the data.)

The convergence properties of the non-concave penalties was studied by [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF]. Under mild regularity conditions on λ n and p λ the non-concave penalties have the oracle properties. The only restrictive condition is that lim sup n→∞ lim sup θ→0 p λn (θ) /λ n > 0.

(1.27)

The minimization of non-concave functions is a computationally difficult task. Derivative-based optimization methods are not efficient in this context, because (i) the derivative of the penalty is not always continuous and (ii) the second order derivative p λ (|θ|) takes infinitely small values. Specific optimization methods have to be used to derive numerical estimation procedures that are stable and fast to compute. These methods are introduce in the framework of MM optimization, in Section 1.3.1.

Hard-thresholding

Consider the L 0 penalty

p λ (|θ|) = λ 2 1 θ =0 ,
(1.28) also called the entropy penalty.

Note that in orthogonal design, the penalized least squares is y i -β i 2 + λ 2 1 β i =0 , and since x 0 is constant on (-∞, 0) and (0, ∞), it takes its minimal value either at y i or at 0. The minimal value is taken at 0 when λ 2 < y 2 i . Thus, the thresholding function of the L 0 penalty is the so-called hard thresholding rule (Figure 1.7): βhard = β1 |β|<λ .

(1.29)

This thresholding only clips the values below the threshold λ to zero, and leaves the other values unmodified. As discussed, this penalty makes the computation of the estimation NP-hard.

Chapter 1. Introduction

Let us now consider the penalty function

p λ (|θ|) = λ 2 -(|θ| -λ) 2 1 |θ|<λ .
(1.30) Surprisingly, under orthogonal design, this penalty also gives the hard thresholding rule (see [START_REF] Antoniadis | Wavelets in statistics: A review[END_REF][START_REF] Fan | Comments on "Wavelets in Statistics: A Review[END_REF]. Since it is a much simpler quantity to minimize, this penalty is used when we want to have an unbiased variable selection method. This penalty is also referred to as hard thresholding penalty; it is represented in Figure 1.6.

The discontinuity in the hard thresholding makes the model selection unstable. The next penalty offers to remedy this problem, while still having an unbiased estimate of the parameters.

SCAD

The smoothly clipped absolute deviation (SCAD) penalty was introduced by [START_REF] Fan | Comments on "Wavelets in Statistics: A Review[END_REF] and analyzed by [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF]. It is defined by

p λ ( β j ) = λI β j ≤ λ + λ aλ -β j + (a -1) λ I β j > λ , (1.31)
for some a > 2 and for β j > 0. This penalty is linear over [0, λ], parabolic over [λ, aλ], and constant over [aλ, ∞). Therefore it behaves like the lasso for small values of the data and like the hardthresholding penalty for large values of the data. Under orthogonal design, the SCAD has the following thresholding function (see Figure 1.7):

βj =        sgn β j β j -λ + when β j < 2λ (a -1) β j -sgn β j aλ / (a -2) when 2λ < β j ≤ aλ β j when β j > aλ (1.32)
This thresholding rule equals that of the lasso for β j < λ and equals the hard thresholding for β j > aλ.

In addition to λ, the parameter a also needs to be determined. We can use cross-validation over a two-dimensional grid, but the computational cost can be deterrent. [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF] performed simulations under orthonormal design with a prior distribution β j ∼ iid N (0, aλ) and estimated the L 2 risk E[ ββ * 2 ]. They proposed to take the value a = 3.7 and noted that the risk does not vary a lot for different values of a.

Note that the SCAD is the only non-concave penalty which is not proportional to λ. This has no influence on the minimization of (1.26) for a fixed value of λ. However, this can make it more complicated to minimize (1.26) for a sequence of values of λ. This difference can make the added computational cost of the SCAD deterrent. We refer to [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF] for more details.

Logarithmic penalty

In order to recover highly sparse estimates, non-concave penalties are required to have a "sharp" spike at zero, or in other words, to have a derivative that vanishes away from zero. We have discussed that penalties which are constant away from zero are harder to minimize. An obvious and simple choice is then to use the logarithmic function as a penalty: θ → log(|θ|). Since the logarithm equals -∞ at zero, this function does not comply as a penalty function. To remedy this issue, we cap the log function away from -∞ in the following manner: define the "log penalty"

p(|θ|) = log(|θ| + ε), (1.33)
where ε > 0 is a numerical constant bigger than computer precision, but smaller than the order of magnitude of the estimates. Typically, we chose a value of ε such that ε min{|β ols j |, 1 ≤ j ≤ p}. Surprisingly, this penalty has never been formally introduced or even used in penalized mle estimation, to the best of our knowledge. It was however used in the field of compressed sensing by [START_REF] Candès | Enhancing Sparsity by Reweighted 1 Minimization[END_REF], which also uses a small constant ε to cap the values of the penalty away from -∞.

Note that [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF] make reference to the penalty p(|θ|) = log(|θ|). More precisely, they define a numerical procedure whose corresponding non-concave penalty would be the logarithm function.

Available implementations

We gather in this section a short summary table of the aforementioned penalties with references to their available implementations in R: see Table 1.1. No known implementations of the Berhu penalty and the log-penalty have been found.

Subset selection

In this section we present the methods of variables selection which are not based on penalized estimation. Instead, the methods select the best set of variables by choosing a set of variables and fitting the model without penalty. These methods are simple to implement and are historically important, but since the set of subsets of variables has a cardinal growing exponentially with the number of variables, they become impractical in the high dimensional case p n. They have also been criticised for being unstable in selection (see [START_REF] Breiman | Heuristics of instability and stabilization in model selection[END_REF].

Best subset selection

The naive approach to model selection is to fit all possible models without penalization and to compare them with a predefined metric. The metric can be a cross-validation procedure or a model selection criterion, as the AIC [START_REF] Akaike | A New Look at the Statistical Model Identification[END_REF] or the BIC (Schwarz, 1978). Since each model is defined by a subset of the variables 1, • • • , p , this approach is called best subset selection. This method is computationally feasible only for small number of variables. The total number of subsets of 1, • • • , p is 2 p . When p = 30, there are approximately one billion models to fit. In the case of model regression, where model fitting is computationally fast, this is the limiting order of magnitude: when p 40, best subset selection if no longer viable.

Branch and Bound. [START_REF] Furnival | Regressions by Leaps and Bounds[END_REF] introduced the "branch and bound" approach (also called leaps and bound) to variable selection (see also [START_REF] Hand | Branch and bound in statistical data analysis[END_REF]. This algorithm uses a trick to reduce the computational cost of best subset selection. Consider the subsets of 1, • • • , p as a binary tree, with the full model at the root and all the possible submodels at the leaves; one moves from a vertex to its two children by either keeping or removing one variable. The naive implementation of the best subset method fits all submodels in this tree. The branch and bound does the same, but uses bounds over the model selection criterion to prune some parts of the tree that are suboptimal compared to the rest of the tree. For the sake of the presentation, assume that the model criterion is the AIC: for model m with p m parameters and with MLE estimate θm , we select the model with the lower value of AIC(m) = -2 ( θm ) + 2p m . Consider a set of models M and define the model m u which includes all variables present in any model m ∈ M and the model m i which includes only the variables present in all models m ∈ M . From the inequalities

( θm i ) ≤ ( θm ) ≤ ( θm u ) (1.34) p m i ≤ p m ≤ p m u , (1.35)
we know that for any model m ∈ M , the AIC has the following lower bound: AIC(m) ≥ -2 (θ m u ) + 2p m i . (We can derive similar lower bounds for other selection criteria.) The branch and bound explores the binary tree from the root and keeps track of the lowest AIC encountered so far. If the current vertex M has a lower bound greater than the lowest AIC in memory, then M cannot include the model with the lowest AIC. It is useless to explore M , and the branch and bound removes it from the tree search. This method is implemented in R in the package leaps [START_REF] Lumley | Leaps: Regression Subset Selection[END_REF].

The branch and bound makes the number of explored models way smaller in practice, but it is still exponential with p. For p > 40, the best subset selection becomes unviable. We then have to compare a smaller number of submodels; this is the approach developed in the following section.

Stepwise Selection

Stepwise selection consists in exploring a small subset of all the possible submodels. Consider the Hasse diagram of the submodels partially ordered by the inclusion between models (see Figure 1.8). On one extremity of the diagram is the model with no variables and on the other extremity is the model with all variables. The stepwise selection methods explores one path between the null model and the full model, using a greedy approach. We differentiate between three strategies of exploration. The forward-stepwise method -or forward selection -starts at the null model and incrementally adds variables to the model, drawing a path in the direction of the arrows in Figure 1.8. The backward-stepwise method -or backward elimination -starts at the full model and incrementally removes variables from the model, drawing a path in the opposite direction than the arrows in Figure 1.8. The forward-backward method is a combination of the two methods: at each step, we choose to either remove or add a variable.

The greedy decision rule for the removal or addition of a new variable x j is based on testing whether βj = 0. When comparing the potential variables to add (or remove), it is not necessary to fit the whole model. The new estimate is computed using the current estimate and an update of the QR decomposition of X T c X c , where X c is the current design matrix. This considerably improves the computational burden of this method. For more information on stepwise selection procedures, see Miller (2002, Section 3). These methods are implemented in the R package leaps [START_REF] Lumley | Leaps: Regression Subset Selection[END_REF].

Iterative penalized methods

MM optimization

MM optimization (see Lange, 2004, Section 6) is a family of optimization procedures. The MM algorithm generalizes the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data Via the EM Algorithm[END_REF]. The acronym MM stands for both majorization-minimization (if we want to minimize a function) and minorization-maximization (if we want to maximize it). For the sake of the presentation, we focus on the majorization-minimization version. Let f β be the function to minimize. Let β (k) be the current point. We say that a function g β|β (k) majorizes f if it satisfies:

   g β (k) |β (k) = f β (k) g β|β (k) ≥ f β , ∀β ∈ C (1.36)
where C is the convex set on which both f and g are defined. The idea of the MM algorithm is to minimize g (β|β (k) ) instead of f (β) directly.

In MM optimization, the parameter at the next step β (k+1) is defined as any vector verifying g (β (k+1) |β (k) ) < g (β (k) |β (k) ). We could set β (k+1) = arg min β g β|β (k) , but this is not required -and is not considered optimal for most choices of g and of the minimization algorithm for g (β|β (k) ). The MM optimization comes from the fact that for any g satisfying (1.36):

f β (k+1) ≤ g β (k+1) |β (k) ≤ g β (k) |β (k) = f β (k) ,
(1.37) so that the MM iterations always reduces the value of f . An optimization procedure with this descend property is praised for being stable. Notice that (1.37) does not imply convergence of the algorithm to the global minimum. Under mild regularity conditions [START_REF] Lange | Optimization. Springer Texts in Statistics[END_REF], the MM converges to a local minimum. Except when f is convex, we cannot ensure that MM converges to the global minimum.

The effectiveness of the MM optimization depends highly on the choice of the majorizing function. Its expression must be simple, so that its minimization is not too computationally intensive. Yet it must be "close" enough to f , in the sense that g (β|β (k) )f (β (k) ) does not increase too much in the neighborhood of β (k) . Each choice of a majoring function g yields a different MM algorithm. The function g is usually a local (linear or quadratic) approximation of the function f around the current point β (k) .

In the following section, we will show that many iteratively weighted penalized estimation arise as the MM algorithm of some non-concave penalty.

Remark: EM is a special case of MM. We show here that the Expectation-Maximization (EM) algorithm is a specific case of the MM optimization. The EM algorithm allows to perform inference in the presence of missing, or unobserved data, and was introduced by [START_REF] Dempster | Maximum Likelihood from Incomplete Data Via the EM Algorithm[END_REF]. We refer to McLachlan and Krishnan (2007) for an exhaustive reference to the EM.

We give a definition in a simple setting. Assume that you observe data from a random variable Y . Assume also that Y is incomplete and that it comes from a random variable X that is not observed. Define θ the parameter of the distribution of X and note f (x|θ) and g (y|θ) the respective densities of X and Y conditionally on θ. The goal is to make inference on g (y|θ). In the case of missing data, we have X = (Y , Z ), where Z is the missing data. But the EM works in a more general setting, where we assume that there is a function t that maps X to Y . The trick for the statistician is to find a complete r.v. X whose distribution (and likelihood) is simple.

We can now define the EM algorithm. It consists of iterating between two steps: the E steps defines the conditional expectancy of the missing data conditionally on the observed data (and of the parameter). The M step maximizes this density. By sequentially maximizing a marginalized density, we can infer the parameter of the missing data. Q(θ|θ (k) ) ← E log f (X |θ)|Y , θ (k) 5:

θ (k) ← arg max θ Q(θ|θ (k) ) 6: k ← k + 1 7:
end while 8: end function It is clear from the algorithm that the E-step plays the role of a local approximation of g (y|θ). The next proposition shows that the E step consists in fact of a function minorizing log g (y|θ). (k) )) minorizes log(g (y|θ (k) )) in the sense of (1.36):

Property 3. The function q(θ|θ

(k) ) = Q(θ|θ (k) ) -Q(θ (k) |θ (k) ) + log(g (y|θ
q(θ|θ (k) ) ≤ log(g (y|θ)), ∀θ q(θ k |θ (k) ) = log(g (y|θ (k) )).

Thus E-step is the local approximation using q(θ|θ (k) ) and the M-step is the maximization of q(θ|θ (k) ).

Proof. The equality case is trivial. It suffices to prove that Q(θ|θ (k) ) -log(g (y|θ)) ≤ Q(θ (k) |θ (k) )log(g (y|θ (k) )). We need to apply the Lemma 1. The detailed proof of the applicability of this lemma is somewhat technical, and is skipped here (we refer to Lange, 2004, p. 139). We have

Q(θ|θ (k) ) -log(g (y|θ)) =E log( f (X |θ))|Y = y, θ (k) -E log(g (Y |θ))|Y = y, θ (k) =E log f (X |θ) g (Y |θ) |Y = y, θ (k) ≤E log f (X |θ (k) ) g (Y |θ (k) ) |Y = y, θ (k) =Q(θ (k) |θ (k) ) -log(g (y|θ (k) ))
This completes the proof.

Lemma 1. Given two almost surely positive densities u and v, E h log(h) ≥ E h log(k) .

MM optimization applied to non-concave penalties

MM optimization can be used to minimize (1.26). Define f (β) = (β) + p j =1 p λ (|β j |) as the function to minimize, where p λ (|θ|) is a non-concave function. The optimization algorithm will depend on the choice of the dominating function g . In this section, we present two possibilities for g , which yield to different numerical procedures. Even though they aim at minimizing the same quantity, theses two algorithms have different properties in practice.

Local Quadratic Approximation

Fan and Li (2001) proposed a dominating function based on the Local Quadratic Approximation (LQA) of the penalty function. The original paper introduced this approximation as a mere numerical trick. It was formalized as an MM procedure by [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF].

The first order approximation of p λ (|θ|) around θ (k) writes

p λ (|θ|) = p λ (|θ|)sgn(θ) p λ (|θ k |) |θ k | θ (1.38)
Taking the antiderivative of the both terms in the previous equation yields the following approximation of the penalty (see Figure 1.9b):

p λ (|θ|) q λ (θ|θ (k) ) p λ (|θ (k) |) + (θ 2 -θ (k) 2 ) p λ (|θ (k) |) 2|θ (k) | (1.39)
We will establish that q λ (θ|θ k ) majorizes p λ (|θ|). We prove the property for θ ∈ (0, ∞) and since the same reasoning applies to θ ∈ (-∞, 0) (and to θ = 0 by continuity), the property will be proven over R. For θ ≥ 0, we have

q λ (θ|θ (k) ) -p λ (|θ|) = θ p λ (|θ (k) |) |θ (k) | -sgn(θ)p λ (|θ|) (1.40) = θ p λ (|θ (k) |) |θ (k) | - p λ (|θ|) θ .
(1.41)

Since θ → p λ (|θ|)/θ is non-increasing over (0, ∞), the right hand-size term of the last equation is strictly positive over (0, θ (k) ) and strictly negative over (θ (k) , ∞). Consequently q λ (θ|θ (k) )p λ (|θ|) is decreasing over (0, θ (k) ) and increasing over (θ (k) , ∞). With q λ (θ (k) |θ (k) )-p λ (|θ (k) |) = 0, this completes the proof. The function p(|θ|) and its LQA q(θ|θ (k) ) around θ (k) = 0.5 are represented in Figure 1.9a. When minimizing q(θ|θ (k) ), we can neglect the constant terms. We derive the following dominating function for the penalized nll:

g (β|β (k) ) = (β) + p j =1 p λ (|β (k) j |) 2|β (k) j | β 2 j .
(1.42)

The minimization of g (β|β (k) ) is a weighted ridge problem whose solution is explicit (from Equation 1.6). The MM algorithm iterates over (1.43) until convergence, where Ω λ (β (k) ) is a diagonal matrix with diagonal entries p λ (|β (k) j |)/2|β (k) j | . A notable special case is the linear model, for which (1.43) simplifies to:

β (k+1) = β (k) -∇ 2 (β) (k) + Ω λ (β (k) ) -1 ∇ (β (k) ) + Ω λ (β (k) )β (k) ,
β (k+1) = (X T X + Ω λ (β (k) )) -1 X T y,
(1.44)

The denominator in (1.39) is source of numerical instabilities when β j gets close to zero. To overcome this problem, we can either (i) Set β j to zero if its absolute value becomes smaller than a fixed threshold, and continue the iterations without the j th component This solution is easy to implement, but once a component of β is set to zero, it will never be selected again. Thereby, the iterative estimate has same disadvantage as stepwise backward elimination: once a variable is removed from the model, it is never considered again. This restriction can be very limiting in practice.

(ii) Bound the denominator away from zero in (1.42) by fixing a small ε > 0 and solving

β (k) = arg min β (β) + p j =1 p λ (|β (k) j |) 2(|β (k) j | + ε) β 2 j , (1.45)
as suggested by [START_REF] Hunter | Variable selection using MM algorithms[END_REF]. This perturbed version is not an MM algorithm any more, and the descent property is not verified. [START_REF] Hunter | Variable selection using MM algorithms[END_REF] showed that (1.45) is the dominating function of a perturbed version of the negative nll, and the LQA estimator converges towards the arg max of this quantity. Under a mild regularity condition [START_REF] Hunter | Variable selection using MM algorithms[END_REF] on f (β), the limit estimator obtained by (1.45) is close to that of (1.44) when ε tends to zero.

Local Linear Approximation

Zou and Li (2008) introduced a dominating function based on the local linear approximation of the penalty function. A first order approximation of p λ (|θ|) around the current point θ (k) gives the Local Linear Approximation (LLA):

p λ (|θ|) l λ (θ|θ (k) ) p λ (|θ (k) |) + (|θ| -|θ (k) |)p λ (|θ (k) |) (1.46)
Lets us prove that this approximating function majorizes p λ (|θ|). Since both terms of (1.46) are even with respect to θ, one can assume θ > 0 in the following proof. The local linear approximation l (θ|θ (k) ) of the bridge penalty p(|θ|) = |θ| 0.5 around the current point θ (k) = 0.5 is represented in Figure 1.9a. From

l λ (θ|θ (k) ) -p λ (|θ|) = p λ (|θ (k) |) -p λ (|θ|), (1.47)
and the concavity of p λ , it follows that l λ (θ|θ (k) )-p λ (|θ|) is minimal at θ (k) , where it equals zero. This completes the proof. Derive the following dominating function for the penalized nll:

g (β|β (k) ) = (β) + p j =1 p λ (|β (k) j |)|β j |. (1.48)
Minimizing this quantity is a weighted lasso problem. When the model at stake is the linear model or the generalized linear model, g (β|β (k) ) can be minimized efficiently using the LARS algorithm.

The LLA iteration step is more computationally intensive than that of the LQA, which is explicit. However the fact that each iteration step of the LLA is a lasso problem has deep consequences. The LLA converges towards a sparse model using a sequence of sparse models. Consequently, the LLA need not be very close to converge to provide a sparse model. Comparatively, the LQA's converging sequence is only asymptotically sparse, and many iterations might be necessary to obtain a good approximation of the limiting sparse model. In practice, much fewer iteration steps are carried out for the LLA than for the LQA.

The MM converges to a local minimum, so the choice of the initial guess β (0) is important. For both LLA and LQA, the standard choice is to set β (0) as the unpenalized maximum likelihood. When the model is not overparametrized, this initial guess is considered to be not too far from the global minimum. [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF] go as far as introducing the one-step estimator, which is obtained by performing only one step of the MM algorithm:

β (1) = arg min β y -X β 2 + p j =1 p λ (|β (0) j |)|β j |, (1.49)
where the initial guess is the OLS:

β (0) = (X T X ) -1 X T y.
This estimator is a weighted lasso estimator, similarly to [START_REF] Zou | The Adaptive Lasso and Its Oracle Properties[END_REF]'s adaptive lasso. In fact, for 0 < q < 1, the one-step estimate with the L q penalty is the adaptive lasso with parameter γ = 1q > 0.

Adaptive Ridge

Definition of the adaptive ridge. The Adaptive Ridge is an iterative penalized regression method.

It shares some similarities with the adaptive lasso, with the notable difference that each iteration is less computationally expensive. As its name suggests, it consists in iterating over a weighted ridge penalty, with the value of the weights adapted at each iteration. The adaptive ridge method was used by [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF] as a numerical trick for approaching an L 0 norm penalty. It was analyzed in a more general setting by [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF] as an iterative procedure to numerically approach any L q penalty, for q ≥ 0. This method iterates over a weighted L 2 penalty problem while adaptively changing the weights at each iteration. It is defined as the following procedure:

(i) Set w (0) = 1 (ii) β(k) = arg min β β + λ p j =1 w (k-1) j β 2 j (1.50) (iii) w (k) j = (|β (k) j | γ + ε γ ) (q-2)/γ (1.51)
(iv) Iterate between the last two steps until convergence.

Equations (1.50) and (1.51) are justified empirically by the following reasoning. Assume β (k) is reasonably close to the limit estimate β (∞) . Then β (k+1) will be not too far from β (k) and the quantity can be well approximated by

w (k-1) j β 2 j = β 2 j (|β (k) j | γ + ε γ ) (2-q)/γ β 2 j (|β j | γ + ε γ ) (2-q)/γ .
This function is differentiable and approximates well the L q penalty. It is represented in Figure 1.10 for different values of q. The solid lines represent the penalties and the dotted lines represent their approximations. The parameter q is usually chosen a priori ; it is rarely cross-validated in q = 0.7 q = 1 q = 0 q = 0. practice. Namely, for q = 1, the Adaptive Ridge gives another algorithm for computing the lasso. The estimating function is continuous even in the case where q = 0 and the penalty is not continuous.

The approximation is very close to the L q penalty on the sets (-∞, ε) and (ε, ∞), thus choosing a very small value for ε will increase the quality of the approximation. However ε has to be large enough to avoid the numerical instabilities in (1.51). Using a theoretic argument, [START_REF] Hunter | Variable selection using MM algorithms[END_REF] proposed a value that depends on the smallest valued OLS component. Using a simulation study, [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF] proposed the value 10 -8 .

The parameter γ > 0 determines how well the approximating function tries to fit to the penalty around zero. Simulations carried by [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF] highlight that setting a large value of γ is however not necessary in practice and show that γ = 2 seems a reasonable choice.

The L 0 adaptive ridge. The adaptive ridge is numerically efficient to approximate any L q penalty, q ≥ 0. Setting q > 0 yields a procedure to numerically estimate the bridge estimator. Since, to our knowledge, the L 0 norm has never been efficiently approximated in the contexted of penalized maximum likelihood estimation, the case q = 0 is of particular interest. We refer to this case as the "L 0 adaptive ridge". Note that this term does not signify that the penalty used is the L 0 norm. The reweighing step of the L 0 adaptive ridge with γ = 2 writes:

w (k) j = 1 β (k) j 2 + δ 2
In orthogonal design, the limit estimate (k → ∞) of the L 0 adaptive ridge is:

f (β ols j ) =      0 |β ols j | ≤ 2 λ 1 2 β ols j + sgn(β ols j ) β ols j 2 -4λ |β ols j | > 2 λ (1.52)
which is represented in Figure 1.7. This thresholding function is equal to zero on the interval [-2 λ, 2 λ], which proves that the limit estimate of the L 0 adaptive ridge is sparse. The adaptive ridge iteration (with γ = 2) has already been used by [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF] in the context of compressed sensing. In compressed sensing, the signal y can be recovered completely from its encoding β since the matrix X has more columns than rows. The author iteratively solves min

β p j =1 w j β 2 j s.t. y = X β,
instead of (1.50), but the weights w i are adaptive as in (1.51).

In their implementation, the authors adaptively reduce the value of ε as the absolute error between successive values of β(k) decreases. To ensure numerical accuracy, the value of ε is capped at 10 -8 . The authors compare the square log penalty to [START_REF] Candès | Enhancing Sparsity by Reweighted 1 Minimization[END_REF]'s log penalty. The square log penalty yields higher recovery rates of y with very sparse reconstructed signal β.

The L 1 adaptive ridge. The adaptive ridge could be used to approximate the L 1 norm, by setting q = 1 in (1.50). The L 1 adaptive ridge bypasses the numerical difficulty posed by the minimization of the L 1 norm. However, different tools have been developed to solve the lasso problem efficiently (see Section 1.1.2). Consequently, this version of the adaptive ridge is not particularly superior to the lasso as far as computational speed is concerned. We note that it was already introduce by [START_REF] Vogel | Iterative Methods for Total Variation Denoising[END_REF], who define the iterative procedure as a "lagged diffusivity fixed points iteration". The authors developed this method in the framework of image denoising and the minimization of the penalized likelihood (Equation 1.50) was conducted using the Newton-Raphson method.

Relation to similar procedures

Weighted L 1 regression. [START_REF] Candès | Enhancing Sparsity by Reweighted 1 Minimization[END_REF] introduced a weighted L 1 regression problem in the context of compressed sensing. The algorithm uses a weighted L 1 penalty (see Algorithm 3). The constant ε > 0 ensures that the procedure stays stable when a coefficient is set to zero. It also helps against numerical instability issues arising when dividing by small numbers. [START_REF] Candès | Enhancing Sparsity by Reweighted 1 Minimization[END_REF] offer to adapt the value of ε at each iteration, setting it to the smallest absolute value of the coefficients of the current estimate (but no smaller than 10 -3 ).

The formula for the weights is very close to that of the non-iterative procedure adaptive lasso (see Equation 1.18) with γ = 1. In compressed sensing, the matrix X has more columns than rows and the problem min y -X β is ill-posed. Consequently there is no OLS estimate to plug in the formula of the weights. This justifies the use of an adaptive procedure where the OLS estimate in Equation 1.18 is replaced by the current estimate.

Algorithm 3 Candes's weighted L 1 minimization

1: function ITERATIVE L 1 (y, X , λ, ε) 2: w (0) ← 1 3: k ← 0 4:
while not converge do 5:

β (k) ← arg min β β + λ p j =1 w (k) j |β j | 6: w (k) j ← 1 |β (k) j |+ε 7:
end while 8: end function Consider now the non-concave penalized problem min

β β + λ p j =1 log |β j | + ε .
(1.53)

Numerical performance

This procedure is strikingly easy to implement and fast to compute. Indeed, the ridge problem is easier to solve than the lasso problem and the reweighing step has complexity O p . Comparatively, the lasso estimate is more computationally intensive. This makes the adaptive ridge quicker to compute than [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF]'s one step estimates an [START_REF] Candès | Enhancing Sparsity by Reweighted 1 Minimization[END_REF]'s adaptive procedure for the log penalty. The downside of the adaptive ridge is that the k-th step estimate is not sparse: only the limit estimate β(∞) is sparse. In practice, the adaptive ridge procedure must be stopped when it is estimated to have converged, and the estimate's coordinates that are close to zero must be set at zero.

Structured variable selection

The methods introduced previously are useful in cases with the prior assumption that β * is sparse, i.e., many of its components are zero. In some applications, the true parameter is assumed to have a simple structure other than having many components to zero. This specific sparsity structure can vary greatly and depends on the relationship between parameters. For instance if β * is assumed piecewise constant, this prior is expressed as the sparsity of (β j -β j -1 ). In general, many prior structures of the parameters are expressed as the sparsity of a transformation of the data. This section presents different possible transformations available to ensure different sparsity structures.

Group lasso

In some practical situations in statistical regression, the variables belong to a group. An important example of such situation is the regression against factor variables that are represented as dummy variables. Another instance of grouped variable occurs when we know a priori that some variables have a common effect on the response variable. In these cases, it makes little sense to perform variable selection that would select some variables inside a group and not the others. Variable selection techniques have been developed to take into account the presence of groups, and to select groups of variables together. [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] introduced the group lasso, building on ideas by [START_REF] Bakin | Adaptive Regression and Model Selection in Data Mining Problems[END_REF] and [START_REF] Lin | Component Selection and Smoothing in Smoothing Spline Analysis of Variance Models[END_REF]. Consider a partition G = (G k ) 1≤k≤K of the set of indices of variables 1, • • • , p . Assume that the variables are partitioned into K groups which are given by G . For k ∈ {1, • • • , K }, p k denotes the number of covariates in the k-th group (we have k p k = p) and β k denotes the vector containing Figure 1.12: Unit ball of the group lasso norm with three variables and with G = {{1, 2} , {3}}. We represent the 3d ball (a) and its restriction to the variables {x 1 , x 2 } (b) and {x 1 , x 3 } (c).

only the variables in G k . Without loss of generality, we can write

β = (β T 1 , • • • , β T K ) T .
The group lasso is defined as the mimimizer of

(β) + λ J k=1 w k β k q ,
(1.54) with positive weights w p j > 0 . This group norm, also called L 1 /L q norm or mixed norm [START_REF] Bach | Optimization with Sparsity-Inducing Penalties[END_REF] is the norm obtained by a L 1 averaging of the L q norm defined on each group. It behaves as an L q norm between variables of each group and as an L 1 norm between groups. Consequently, for q > 1, the L 1 /L q norm performs selection between groups and performs shrinkage (without selection) within each group. It is easy to see that the L p /L q penalty is a norm for every p ≥ 1, q ≥ 1. Notice that when the partition of 1, • • • , p is the trivial partition composed of singletons, (1.54) becomes the lasso penalty and when there is only one group, (1.54) becomes the ridge penalty. The group lasso can therefore be seen as a compromise between the ridge and the lasso penalty.

There are several possible choices for q. Some authors [START_REF] Berwin | Simultaneous Variable Selection[END_REF], [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF]) consider the L 1 /L ∞ mixed norm, mainly because optimization with this penalty is numerically efficient (see [START_REF] Mairal | Convex and Network Flow Optimization for Structured Sparsity[END_REF][START_REF] Zhao | Grouped and Hierarchical Model Selection through Composite Absolute Penalties[END_REF]. However the preferred choice in the literature is the Euclidean norm, q = 2, which shrinks equally in all directions.

We now illustrate the group penalty on a simple example. Consider the case of p = 4 variables divided into K = 2 groups: G 1 = {1, 2} and G 2 = {3, 4}. The group lasso penalty writes

β 2 1 + β 2 2 + β 2 3 + β 2 4 .
This mixed norm is a L 1 norm with respect to the groups p 1 and p 2 , and is a L 2 norm inside each group. Consequently, the L 1 norm will enforce selection of (β p 1 , β p 2 ) and will not enforce selection of variables inside each group. The variables 1, 2 can be selected together but not aloneand the same goes for 3, 4. The effect of the variables 1, 2 inside the group p 1 is given by the L 2 norm. The selection effect of the mixed norm is illustrated in Figure 1.12, which represents the unit ball in the case where the p = 3 variables are partitioned by G = {{1, 2} , {3}}. This ball is a double cone with singularities at its top and bottom and at its central circle. As is the case for any penalty norm, these singularities give the sparsity structure encoded by the mixed norm: either β 1 and β 2 are jointly set to zero (in which case x 3 is estimated with shrinkage), or only β 3 is set to zero (in which case x 1 and x 2 are estimated with ridge shrinkage).

The weights w p j are included to adjust for the relative importance of the groups. If one group has twice as many components as another, it will be regularized 2 times less -assuming all variables have equal importance. Consequently, the choice w p j = p j is often used to renormalize the importance of each group. [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] introduced a group LARS algorithm, a group version of the LARS algorithm which successively adds groups of variable to the active set. Contrarily to the LARS, which can be easily modified to compute the lasso solution [START_REF] Efron | Least Angle Regression[END_REF], the group LARS generally differs from the solution of (1.54). However, since the L 1 /L 2 penalty is a norm, the penalized likelihood (1.54) is convex and the group lasso solution can be computed in reasonable time using the coordinate descent algorithm (or any convex optimization method). [START_REF] Roth | The Group-Lasso for Generalized Linear Models: Uniqueness of Solutions and Efficient Algorithms[END_REF] presented a comp utationally efficient algorithm to compute the group lasso.

The group lasso can be naturally extended to any mixed norm. One could for instance consider the L p /L q mixed norm, where p < 1 and q > 1, such that the norm between groups induces sparsity and the norm inside does not induce sparsity. [START_REF] Berwin | Simultaneous Variable Selection[END_REF] introduced a version of the group lasso using the L 1 /L ∞ norm.

We emphasize that the group lasso is to be used when clear a priori knowledge is known about the variable group. Theoretically, we could set the partition as a Bayesian variable and try to find the best partition. We could explore the set of partitions of {1, • • • , d } by following the vertices of the Hasse diagram defined by the partial order on partitions (see Figure 1.8). We could explore this set using a discrete MCMC algorithm. However, in practice, this Hasse diagram is highly connected so exploring this space with an MCMC approach is computationally unfeasible.

Overlapping groups

The previous section considers a bagging of the covariates based on a partitioning of 1, • • • , p . This group sparsity structure corresponds mainly to the case of factor covariates represented using dummy variables. In this section we consider generalizations to subsets of 1, • • • , p which do not form a partition. It turns out that the sparsity obtained from such overlapping groups are of great use to include interesting a priori structure between the variables. This extension of the group lasso was introduced by [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF] and [START_REF] Jacob | Group Lasso with Overlap and Graph Lasso[END_REF].

We keep the same notation as in the previous section, except that here

G = (G 1 , • • • , G K
) is an element of the power set of 1, • • • , p . Without loss of generality, we can assume that each variable is present in at least one group G k . We consider the same penalty as before:

(β) + λ K k=1 β k q k , (1.55)
where q k > 1 is the order of the penalty norm inside G k . As already discussed in the previous section, setting q k = 2 for all k will give a ridge penalty inside each group. First, note that the mixed norm (1.55) remains a norm when G is no longer a partition. Consequently the penalty is convex and the group lasso has the same computational cost when the groups have overlapping.

We now develop on the type of sparsity structure that overlapping groups can encode. Overlapping groups allow to encode a priori information of the type "if one variable is removed, other variables have to be selected". For the sake of the illustration, consider the following example, with 4 covariates x 1 , • • • , x 4 . The a priori knowledge here is that the variable 1 and 3 are closely related, as well as 1 and 2. Define G = ({1, 2} , {2, 3}). The group lasso estimate is the minimizer of

(β) + λ β 2 1 + β 2 2 + β 2 2 + β 2 3 .
This penalty sets all the variables inside one of the groups to zero. This penalty is represented in Figure 1.13. As illustrated, the singularities enforce that either β 1 = β 3 = 0 or β 2 = β 3 = 0. Consequently the estimated model will either only select the variable 3 or only the variable 2.

There are a large number of possibilities to use overlapping groups to include prior knowledge on the link between variables. However, defining groups that encode the desired sparsity structure is a difficult task, and is somewhat counter-intuitive. Indeed if a variable is selected in the model, all the variables which have one group in common with this variable will be kept in the model. Consequently the selected variables will be formed by the complementary of a union of the groups G k . It is often easier to express the relationship between variables through a measure of their proximity and we want to use the groups defined using this proximity measure. The relevant penalty corresponds to selecting variables whose support is formed by a union of the groups G k . [START_REF] Jacob | Group Lasso with Overlap and Graph Lasso[END_REF] introduced this group norm. For instance, with the group G = {{1, 3} , {2, 3}}, this penalty removes either β 2 or β 3 from the model, and does not select β 3 . This method is of great use in signal processing, where it is easy to specify the prior structure of an object using a union of basic structures. For examples of applications in signal processing and image processing, see [START_REF] Jenatton | Structured Variable Selection with Sparsity-Inducing Norms[END_REF][START_REF] Huang | Learning with structured sparsity[END_REF].

Sparse Group lasso. Friedman et al. (2010a) introduced an extension of the group lasso which enforces the estimate to be both sparse and sparse with respect to its group partition. We use the notations of Section 1.4.1, i.e. G = (G k ) 1≤k≤K denotes a partition of 1, • • • , p . The penalty of the sparse group lasso is obtained by adding an L 1 penalty to (1.54) and choosing q = 2. The sparse group lasso estimate minimizes

(β) + λ 1 β 1 + λ 2 K k=1 p k β k 2 .
Note that this penalty is a mere particular case of the overlapping group penalty. Indeed, the sparse group penalty rewrites as (1.55) with the overlapping group G = G ∪ {1} , • • • , p . In its generalized version [START_REF] Simon | A Sparse-Group Lasso[END_REF], the sparse group lasso uses a weighted mean between the lasso and the group lasso.

Hierarchical structured sparsity

A particular case of the overlapping group lasso is of great interest in many applications. Suppose that we dispose of a hierarchy between the variables, in the sense that one variable can be more important or more primary than another one. We may want to include this a priori knowledge as a constraint on β. If x 1 is "above" x 2 in this hierarchy, we say that x 1 precedes x 2 , and it is noted

x 1 ≺ x 2
The hierarchy corresponds to a directed acyclic graph (DAG) between the variables, with the property that the parent variables precede their children. In the following, we assume that the hierarchy includes all the variables. We identify the variables to the vertices of the DAG.

We can build a grouping of the variables such that each variable will only be considered for removal from the model if all of its children have already been removed. We define G as follows: it contains all the singletons of child-less nodes, then all the groups of size 2 composed of the child-less nodes and their respective parents, then all the groups of size 3 composed of the child-less vertices and their grand-parents, etc until all the variables are included. Note that the groups are nested: if a group contains a variable, it overlaps with all groups containing this variable. There is a one-to-one correspondence between a DAG and its nested group G .

We define the hierarchical group lasso in the same way as before: The following argument illustrates why the overlapping group penalty encodes the desired hierarchical structure. From the way G was constructed, we know that any variable x k appears in (1.56) through an L 1 norm of β k and an L 2 norm of all the groups containing x k and its dominating variables. Consequently, x k will always be removed from the model before any of its ancestors. We provide a simple example to illustrate hierarchical sparsity structures. Consider the three variables

(β) + λ K k=1 β k 2 . (1.56) 1 2 3 (a) Hierarchy
x 1 ≺ x 2 ≺ x 3 . The group for this hierarchy is G = {{1, 2, 3} , {2, 3} , {3}} (see Figure 1.14a). The corresponding penalty is β 2 1 + β 2 2 + β 2 3 + β 2 2 + β 2 3 + β 2 3 .
The unit ball is represented in Figure 1.14b. The two "spikes" in the unit ball encode the sparsity β 2 = 0. But since they are located in the plane β 3 = 0, the sparsity induced by these singularities necessarily induce that β 3 = 0. The shape of the unit ball correctly encodes the sparsity structure of Figure 1.14a. Finally, Figure 1.14c represents the projection of the unit ball on the variables β 1 and β 3 . It illustrate that the third variable will be removed from the model before the first variable.

Fused lasso and total variation

In this section, the parameter vector is assumed to have a natural ordering. This is the case in signal processing, when the response variable y i = y(t i ) is indexed by time and each parametrization reflects this time ordering. [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] was among the first to introduces a variable selection method on a transformation of the data. In their method, called fused lasso, one assumes that the vector is sparse and that it deviates from zero only on intervals over which it is constant. The fused lasso minimizes

(β) + λ 1 p j =1 |β j | + λ 2 p j =2 |β j -β j -1 |, (1.57)
where λ 1 > 0 and λ 2 > 0 are trade-off constants. When taking λ 2 = 0, this estimate becomes the lasso estimate. When taking λ 1 = 0, the penalty becomes the total variation of the parameter. The total variation, defined by TV(x) = n j =2 |x jx j -1 |, quantifies how much a sequence varies. Thus, the total variation penalty enforces the parameter to be piecewise constant. In the linear model with orthogonal design X T X = I , the total variation regularized problem writes arg min

β y -β 2 2 + λ 2 p j =2 |β j -β j -1 |, (1.58)
which is referred to as total variation denoising in the signal processing community (Rudin et al., 1992a). It is used to recover noisy signals under the assumption that the original signal is piecewise constant.

-4

-2 0 2 4 -4 -2 0 2 4 β 1 β 2 Figure 1
.15: Admissible solution set for the fused lasso penalty (dark grey) and for the lasso and total variation penalties (light grey).

The joint use of lasso and total variation penalty enforces the signal to be both sparse and piecewise constant. The admissible solution set for the fused lasso penalty is illustrated in Figure 1.15 (dark grey). In light grey, the diagonal band represents the admissible set of

|β 2 -β 1 | ≤ c 2 and the rotated square represents the L 1 ball |β 1 | + |β 2 | ≤ c 1 , with c 1 = 3 and c 2 = 1.
The fused lasso inherits the following properties from the lasso:

(i) It is as computationally efficient, since [START_REF] Friedman | Pathwise coordinate optimization[END_REF] showed that a small modification of the coordinate descent is viable to solve the fused lasso.

(ii) It enforces the same level of sparsity on the parameter than the lasso. Recall that if p > n, the lasso selects at most n parameters (see Rosset et al., 2004, Appendix A). For the fused lasso, this property remains true when replacing n by the number of constant plateaux:

p j =1 1 β j =β j -1 (with β * 0 = 0).
Two-dimensional fused lasso. Consider now the case where the parameter β is parametrized by two indices: β = (β j , j ), where 1 ≤ j ≤ q and , 1 ≤ j ≤ q , with q × q = p. In some applications, this parametrization appears naturally in the model, which then has a bi-dimensional structure. Then the fused regularization introduced in the previous section generalizes naturally to the two-dimensional case. Notice that in (1.58), we penalize over the difference operator ∆β j β j -β j -1 of the parameter. This operator generalizes in two dimensions to the penalty

q j =1 q j =1 |β j , j -β j -1, j | + |β j , j -β j , j -1 |. (1.59)
which defines the two-dimensional total variation penalty As for the one-dimensional case, we can add regularization over β j directly, so as to form the bi-dimensional equivalent of (1.57), called bidimensional fused lasso. Note that there is not a unique difference operator in two dimensions. This two-dimensional setting arises in image processing: each parameter corresponds to the grey level of one pixel in a numerically-stored image. See [START_REF] Friedman | Pathwise coordinate optimization[END_REF] for an illustration of the bi-dimensional fused lasso, applied to image denoising. Many statistical models dealing with spatial data also falls within the scope of this bi-dimensional parameter structure.

Minimization of (1.59) does not enter in the framework of the coordinate descent, as presented in Property 1. However, the optimization method can be slightly modified so as to solve (1.59), as was shown by [START_REF] Friedman | Pathwise coordinate optimization[END_REF]. Optimization methods for the bi-dimensional fused regularization is not given here in all generality. Instead, we present the computational approach used for the fused adaptive ridge later in this section.

Many extensions of the bi-dimensional fused lasso are possible. First, the L 1 penalty can be generalized to any penalty presented in this Chapter. In particular, this thesis focuses on the adaptive ridge applied on bi-dimensional regularization, which we will call fused adaptive ridge. Second, (1.59) makes use of the one-dimensional finite difference operator of the first order. To enforce β to be close to linear instead of constant, we would use the penalty

4β j , j -β j -1, j -β j +1, j -β j , j +1 , β j , j -1 or 8β j , j -β j -1, j -β j +1, j -β j , j +1 , β j , j -1 -β j -1, j -1 -β j -1, j +1 -β j +1, j +1 + β j +1, j -1 .
The two serve the same role, the difference between them being out of the scope of this work. Extensions to any order are available, although the first order given in (1.59) is the most used in many applications.

Fused Adaptive Ridge

The work presented in this thesis makes extensive use of the (L 0 ) adaptive ridge, and mostly of the fused adaptive ridge, used in Chapters 2 and 3. We recall the method here and discuss some computational considerations that will be developed further in the corresponding Chapters.

Generalizing (1.59) to the adaptive ridge penalty is straightforward. The weighted penalty of the adaptive ridge becomes a weighted difference along each dimension. Thus, we use two arrays of weights: v and w .

Algorithm 4 Fused Adaptive Ridge

1: function FUSED ADAPTIVE RIDGE(λ) 2: w (0) ← 1, v (0) ← 1 3: k ← 1 4:
while not converge do 5:

β (k) ← arg min (β) + λ q j =2 q j =2 v j +1, j +1 (β j , j -β j -1, j ) 2 + w j +1, j +1 (β j , j -β j , j -1 ) 2 6: v (k) j +1, j +1 ← (β j , j -β j -1, j ) 2 + ε 2 -1 7: w (k) j +1, j +1 ← (β j , j -β j , j -1 ) 2 + ε 2 -1 8: k ← k + 1 9:
end while 10: end function

The method is given in Algorithm 4. Notice that line 6 and 7 are done in complexity O(q ×q ). The computational bottleneck is the minimization of the penalized likelihood in line 5. This is typically done using the Newton-Raphson method. This method requires a number of successive inversions of the Hessian (i.e. the second order derivative) of the penalized likelihood. This matrix equals the Hessian of (β) to which is added the second order derivative of the penalty term (given in line 6). Note by H (β) this function, whose values are (q q ) × (q q ) matrices. Note that the form of the fused penalty induces that H is a (symmetric) band matrix. More precisely, assuming that β is ordered β = (β 1,1 , . . . , β 1,q , β 2,1 , . . . , β q,1 , . . . , β q,q ) it has 5 non-zeros diagonals: the main three diagonals, and two diagonals located at distance ±q of the main diagonal. Inverting H can be done cleverly by making use of its symmetry and bandedness, using QR decomposition for banded matrices. We have found no way of taking advantage of the fact that most diagonals are zero to speed up this inversion; hence the QR decomposition method inverts H in O(q × q × q ), which is usually of the same order as p 3/2 . This remark is only relevant is ∇ 2 (β) is also banded, with a band no bigger than q . This turns out to be the case in our applications. More details of computational considerations are given in the corresponding chapters.

Conclusion

In this introducing chapter I have developed on the many and diverse methods of penalized likelihood estimation that perform selection of variables.

In a first part, I explain in details the simplest and mostly used methods of penalized likelihood: the lasso and its extensions: elastic-net and bridge regression. These methods' have historically been developed from the generalization of the norm penalty: the lasso penalty (L 1 norm), which is often introduced as a relaxation of the L 0 norm, was considered an extension of the ridge penalty (L 2 norm) to enable both shrinkage and sparsity. Most other methods are built as extensions of the L 1 norm, as is the case for the elastic-net and the bridge. A breakthrough was made when [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF] introduced the idea of "non-concave" penalty, generalizing the previous methods to a wide class of sparsity inducing penalties. These non-concave penalties benefit from satisfying asymptotic properties of consistence in selection and asymptotic normality in estimation, with the lowest possible variance. These penalties are hard to minimize, making the computation of their corresponding estimate a difficult problem. [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF] and [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF] introduced two different iterative procedures to minimize these penalties: the LQA and the LLA.

The next Section is dedicated to briefly introducing the other important approach to variable selection: stepwise selection. Instead of enforcing a prior on the distribution of the estimated parameter, like the penalized methods do, the stepwise selection methods (including forward-backward selection) performs variable selection by successively fitting the model with different subset of variables.

In Section 1.3, I resume to penalized methods and I develop on the iterative methods for solving penalized likelihood problems. I also introduce the iterative penalized method that is central to this thesis: the adaptive ridge. I discuss its strengths and particularities and I develop on a particular case of interest: the L 0 adaptive ridge.

Finally in the last section I introduce an important application of sparsity-inducing estimation: structured sparsity. In many applications, the parameter is assumed to have a certain characteristic. Two examples have a historical importance and have played a large role in the development of structured sparsity: (i) in image processing, images are assumed to be mainly composed of large regions of equal brightness and (ii) in many models, the variables are assumed to belong to groups, and we want to select the variables by groups. Structured sparsity focuses on finding ways to apply the penalized likelihood variable selection methods to a specific function of the parameter.

The work of this thesis makes an extensive use of structure sparsity in different applications. We will consistently use the (L 0 ) adaptive ridge with a structured sparsity adapted to the problem at hand.

Why the adaptive ridge? Many efficient penalized likelihood methods have been successfully brought to light and the class of non-concave penalties seem to enjoy both optimal theoretical properties and numerical procedures for fast computation. We see the adaptive ridge as a new and interesting approach to iterative model selection methods. Its simplicity of implementation makes it easily applicable to new statistical problems. Indeed, the adaptive ridge procedure is a weighted ridge problem; for which readily available numerical solvers can be used.

The adaptive ridge was introduced by [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF] and [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF] with an empirical justification for its formulation. In Chapter 1.3.3, we make an important connection between the adaptive ridge and non-concave penalty. The adaptive ridge is the numerical procedure obtained when applying [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF]'s LQA to a particular penalty, which is a relaxation of the logarithmic function. This provides the pertinent framework for a study of the properties of the adaptive ridge. It can be shown that although the relaxation used to generalize the non-concave penalty to the adaptive ridge procedure makes its computation easier, the theoretical properties of the latter do not apply for the former. We lay the ground for a study of these properties.

Appendix: Least angle regression

The least angle regression (LARS) is an iterative method for penalized regression introduced by [START_REF] Efron | Least Angle Regression[END_REF]. A modified version of the LARS, called lasso-modified LARS or abusively LARS, allows to compute the solution path of the lasso. Indeed, the lasso estimate, as a function of the parameter λ, is piecewise linear. As long as the active variables are the same, the estimates are a linear function of λ. Each time a variable is included in or removed from the model, the slope of the other variables' effects shift. The LARS -and the LARS-lasso -make use of this property: the estimates are only computed at the values of λ were a variable is added to or removed from the model. The entire solution path is then drawn by connecting the dots. In this respect, the LARS-lasso compares favorably with the coordinate descent, which computes the lasso estimate only on a grid of values of λ.

The LARS-lasso is initiated with no "active" variables in the model and it iteratively adds variables one by one in the model. The current estimate of y is noted ŷ = X β and is initiated at zero. The LARS-lasso operates by selecting the variable with the greatest correlation with the current residuals ĉ = X T (y -ŷ) and adding is to the active set of variables. The current estimate is then updated in the direction u such that the current residuals have equal correlation with the newly added variable than with the other active variables, that is: ŷ ← ŷ + αu. The value of α is chosen such that the correlation between the current residuals and the new variable becomes equal to the correlation between the current residuals and any other non-active variable. The first non-active variable to satisfy this condition is in turn added as an active variable. Then, the current estimate ŷ is updated in the direction with equal correlation between the two last variables to be activated. By iteration, the current estimate follows the path of equal correlation -i.e. of equal angle -between all covariates, until they are all added to the model. Hence the name "least angle regression".

The LARS-lasso is the same as the LARS algorithm, except for the following small modification. When moving the current estimate along the least angle direction, if an active variable j makes its estimate cross zero: (i) update the current estimate at the corresponding value of α, (i) remove j from the set of active variables, and (iii) carry on the LARS procedure. The LARS-lasso was proven by [START_REF] Efron | Least Angle Regression[END_REF] to make the algorithm compute the lasso solution path, under the condition that no two variables have equal correlation with the current estimate. This condition is always met in practice. The proof is technical and there does not seem to be an intuitive explanation to the connexion between least-angle regression and L 1 norm.

Chapter 2

Regularized estimation of the hazard rate

In this section we introduce a statistical method for hazard rate estimation with automatic detection of breakpoints. This methods uses the adaptive ridge penalization to effectively produce a piecewise constant estimate of the hazard rate, where the cuts are chosen from the data. This method is general and could be applied to a wide range of statistical problems, like density estimation. We only present it in a more specific setting, where instead of estimating a density, we estimate a typical quantity of interest in survival analysis: the hazard rate, which can be seen as a conditional density. This choice is motivated by two reasons. First, this problem was motivated by the study of heterogeneity in survival analysis. Secondly, we present another model in Chapter 3 which answers the problems of age-period-cohort analysis and is built on the method presented in this Chapter. Thus, we choose to present our method in the specific context of survival analysis, and more specifically age-periodcohort analysis.

More precisely, we focus on the estimation of the hazard rate λ(t , u) taken as a function of two variables t and u, were t is the time dependent variable and u is another variable. In particular, we present our method in the context of age-period-cohort analysis. We will show that we can consider only the case of two variables: one is time-dependent (either the age or the period) and the other, the cohort (i.e. date of birth), is not time-dependent. For the sake of simplicity, our method is illustrated with u as the cohort variable.

This chapter is organized as follows. In a first part, we give a general introduction to time-toevent data and the statistical study of censored variables. We go on to introduce a simple model for the estimation of the hazard rate in the case of right-censored data. This model is the piecewise constant hazard model, with automatic detection of the breakpoints. This first model was developed by Bouaziz and Nuel (2017). We present it because it provides a good illustration of application of the fused adaptive ridge to hazard estimation. In fact, the model developed in this Chapter is built as an extension of the work from Bouaziz and Nuel (2017). We then present the setting of age-periodcohort analysis, and more specifically, the estimation of the hazard rate in presence of another continuous variable. Finally, Sections 2.2 and 2.3 present our contribution.

Introduction

Survival Analysis

Survival analysis is the statistical study of duration variables, sometimes called time-to-event variables. This duration variable is the time between an origin and an event of interest. The origin is the moment where the individual starts being at risk of having the event. (In the example where death is the event of interest, the origin is the birth of the individual.) In order to collect time-to-event data we have to wait until the event of interest has occurred. In cases where the duration time is long (as in demographics and epidemiology), studies are stopped before all events have occurred.

Thus, not all the duration times are observed. We can only observe the event of interest (for instance death) if it happened before all events that would make its observation impossible (for instance that the individual exit the epidemiological study). This problem can occur for a large proportion of the sample, for example with demographic studies where the event of interest is death. Removing the corrupted individual from the study and carrying the statistical inference with the remaining sample is a mistake because it would introduce a bias in the observed sample.

This problem is called censoring. When the variable cannot be observed after the occurrence of another incompatible event, it is called right-censoring. In this section, right-censoring is statistically defined and the assumptions under which inference on the event of interest can be recovered are discussed.

Time-to-event data. Let T * ≥ 0 be the time-to-event random variable of interest. In the case of right-censoring, T * is not observed. Instead we observe (T i , ∆ i ), an i.i.d sample of the observed time T and the right-censoring time C :

T min T * ,C (2.1) ∆ 1 T * ≤C = 1 T =T * .
(2.2)

The term right censoring comes from the fact it is the higher values of T * that are censored by Cand hence the right part of the distribution. To observe ∆ i is to know whether we observe T * i or C i .

Terminology. We provide and recall a few terms used in time-to-event data. The event of interest is the event which marks the end point of the duration T * at stake. It is typically death of the onset of a disease. A duration T i is said to be observed if ∆ i = 1, otherwise, it is said to be censored (and

∆ i = 0).
An individual is said to be at risk at a time t if T i > t , that is, if he is at risk of having the event of interest.

Survival function.

The survival function is defined as

S (t ) P (T > t ) = 1 -F (t -) , (2.3)
where F is the cumulative probability distribution function of the variable of interest. The survival function at time t gives the probability of survival at time t -i.e. the probability that the individual has not had the event of interest at time t .

Hazard Rate. The hazard rate (or instantaneous hazard) rate is defined as the infinitesimal probability of the event to occur now, conditionally on the fact that it has not occurred yet:

λ * (t ) lim δt →0 P (t ≤ T * < t + δt | T * ≥ t ) δt . (2.4)
This function is important in survival analysis for (i) theoretical reasons as is illustrated in the next paragraph (ii) practical reasons because it has an important meaning in demographics and epidemiology.

Many simple statistical models are formulated using the hazard rate rather than the survival function. Denote by f the density of T * ; then simple calculations yield λ (t ) = f (t ) /S (t ). This equation shows that there is a one-to-one correspondence between a distribution and its hazard rate function (if it exists).

Assumptions on the censoring variable. Define the following assumptions.

Assumption 3 (Independent censoring). The random variables C and T * are independent

Assumption 4 (Non-informative censoring). The censoring variable C does not depend on the model parameters

Inference under right-censoring. Define the crude hazard rate as:

λ (t ) lim δt →0 P (t ≤ T < t + δt | T ≥ t ,C ≥ t ) δt . (2.5)
We recall the following property (Fleming and Harrington, 2011, p. 25).

Property 4. Under Assumption 3, λ(t ) = λ * (t ).
This property is essential for inference with censored data. Indeed, the crude hazard depends on T and C and not on T * . Consequently, under Assumption 3, inference on T * is straightforward and we can perform inference on T * (through its hazard rate) without directly observing it.

This assumption is not equivalent to Property 4. It is stronger than necessary and we can construct examples of dependent censoring where λ = λ * (Fleming and Harrington, 2011, Exercice 1.8). The precise assumption that is equivalent to Property 4 is

lim δt →0 1 δt P t ≤ T * < t + δt |T * ≥ s = lim δt →0 1 δt P t ≤ T * < t + δt |T * ≥ t ,C ≥ t , (2.6)
for all t such that P (T * > u) > 0 (see Fleming and Harrington, 2011, , Theorem 1.3.1). This equality signifies that C does not convey any information supplementary on T * . However, for the sake of simplicity, we assume independent censoring and not the latter, weaker assumption. This is a simplification often made in the literature. Moreover, wherever necessary, it is implicitly assumed that T * and C are absolutely continuous with respect to the Lebesgue measure.

Remark on censoring and competing risk. As explained previously, the censoring is due to the individual leaving the study. We make remark here to note that not all events preventing the event of interest from being observed are censoring events. Take the example where the event of interest is the birth of the first child. If the individual leaves the study before his first child is born, this individual is censored. But the death of the individual is not a censoring variable, because after dying, the individual is no longer at risk. Censoring requires that the individual still be at risk of having the event of interest. The question "is the individual still at risk ?" allows to understand whether an alternative event is a censoring event or a competing risk. This simple question is often overlooked in applications.

The piecewise constant hazard estimation

The piecewise constant hazard model. In this section we introduce a simple model upon which this chapter's work is built. Let c = c 1 , . . . , c J be a vector of positive-valued cuts sorted in increasing order: c 1 < • • • < c J . In the piecewise constant hazard model the hazard is constant on the intervals defined by c:

λ (t ) = J j =1 exp(α j )1 c j -1 ≤t <c j , (2.7)
where α j are the values of the log-hazard rate over each interval and where c 0 = 0 and c J +1 = ∞ by convention. The logarithmic transformation ensures that the parameter α can take any real value, and thus, no constraints need to be added to its estimation. The parameter of the model is α ∈ R J . The interest of the piecewise constant hazard model comes from its simplicity. Indeed the hazard rate is directly interpretable and in many applications (e.g. medicine or epidemiology), the practitioner is interested in a simple hazard function in order to draw conclusions. On the other hand, with sufficiently enough cuts, this model is very flexible. One of the aim of this chapter is to choose the number and values of the cuts from the data. In this manuscript we use the acronym "PCH" to denote the piecewise constant hazard function.

In this section we will discuss the estimation in the PCH model with an automatic detection of the knots. This gives a good introduction to the contribution of this chapter, which is estimation of bi-dimensional hazard rate with automatic selection of the knots. This work was developed by Bouaziz and Nuel (2017) who also introduced a procedure for inferring α.

Maximum likelihood estimation. Under Assumptions 3 and 4, the likelihood writes

L(α) = n i =1 f (T i ) ∆ i S (T i ) 1-∆ i .
(2.8)

This product has n terms: the information provided by individual i is f

(T i ) if the event is observed (i.e. ∆ i = 1) and S(T i ) = P (T * > T i ) if the event is censored (i.e. ∆ i = 0). The NLL = -log L writes (α) = n i =1 T i 0 λ (t ) d t -∆ i log (λ (T i )) .
(2.9)

Combining (2.7) into (2.9), we get

(α) = n i =1 T i 0 λ(t )d t - n i =1 ∆ i log(λ(T i )) = n i =1 J j =1 e α j T i 0 1 c j -1 ≤t <c j d t - n i =1 J j =1 α j ∆ i 1 c j -1 ≤T i <c j . Define R i , j T i 0 1 c j -1 ≤t <c j d t = 1 T i ≥c j -1 (max(T i , c j ) -c j -1 )
(2.10)

and

O i , j ∆ i 1 c j -1 ≤T i <c j .
(2.11)

The former quantity is equal to the time spent by individual i in the interval [c j -1 , c j ). It is called time at risk. The latter quantity is equal to the number of observed events in the interval [c j -1 , c j ) -note that the censored events are not taken into account here. Then (2.12) where R j = n i =1 R i , j is the total time at risk and O j = n i =1 O i , j is the total number of events, both in the interval [c j -1 , c j ).

(α) = J j =1 exp(α j )R j -α j O j ,
A few remarks need to be made at this point. The two quantities (O j ) j and (R j ) j are exhaustive statistics. They gather all the information from the data (T i , ∆ i ) i with respect to the PCH model. They have a crucial role in any survival model with a piecewise constant hazard. When there is no ambiguity, they will be referred to as the "exhaustive statistics".

Moreover, (α) writes as a very simple function of the exhaustive statistics. In fact Equation 2.12 is also the likelihood of the Poisson regression model whose response variable is O j ∼ P (µ j ) with mean µ j = log(α j )R j . For this reason, the PCH model is often called "Poisson model", even though it does not assume that the O j s be Poisson distributed or that the R j s be fixed quantities. Notice that the MLE is explicit here, since the minimum of (α) is

α mle j = log O j R j .
Fused penalization using the adaptive ridge. The former considerations yield an explicit maximum likelihood estimate of the hazard in the PCH model. We want to regularize over the successive differences of α j . Define the penalized negative log-likelihood corresponding to the fused adaptive ridge regularization:

pen (α, w ) (α) + λ J -1 j =1 w j (α j +1 -α j ) 2 , (2.13)
where λ > 0 is a trade-off constant and (w j ) 1≤ j ≤J -1 is the vector of weights used in the adaptive ridge.

The penalized estimate is obtained by iterating between (2.13) and an update of the weights w , as in Algorithm 4. The principle of the estimation follows the iterative procedure of the adaptive ridge, which is given in Section 1.3.3.

Proportional hazard model with piecewise constant hazard

The previous model generalizes easily to the setting of right censored data with covariate, that is, the data consists of (T i , ∆ i , X i ) 1≤i ≤n , where x i is a p ×1 vector of covariates. We assume that the covariate are not time-dependent. Consider now the PCH model with proportional effect of the covariates, that is:

λ(t |x i ) = λ 0 (t ) exp(x T i β), (2.14)
where λ 0 (t ) is the piecewise constant baseline hazard (see Equation 2.7) and β = (β j ) 1≤ j ≤p is the plength parameter of the covariates' effect. Define θ = (α, β) the (J + p)-length parameter vector of this model: α represents the values of the hazard rate across time and β represents the proportional effect due to the covariates. Equation 2.14 is called proportional hazard model. In this model, the parameter β has an important interpretation. Assume for the sake of the illustration that the j -th variable is binary. Then under this model, two individuals whose covariates are equal except the j -th have proportional hazard rates: In this model, log β j can be interpreted as the multiplicative effect of the j -th covariate on the baseline hazard rate λ 0 (t ). We choose to model a proportional effect of the covariates because it offers a simple model with highly interpretable parameter β. Unlike Cox's model [START_REF] Cox | Regression Models and Life-Tables[END_REF] -which is semiparametric -the baseline hazard is specified here, and the model is fully parametric.

λ(t |x i , j = 1) = λ(t |x i , j = 0) exp(β j ).

Cohort data

In demographics and epidemiology, an important problem is the study of the evolution of the hazard rate of death (or disease onset) with respect to time. But there are different time scales that we can consider.

• the cohort refers to the date of the time origin.

• the period refers to the calendar time (i.e., the date).

• The age refers to the time spent since the time origin.

When the event of interest is the death from a specific disease, the time origin is the onset of the disease and the age variable is the time since the disease onset.

When the event of interest is death or the onset of a disease, the time origin is the birth of each individual. Then, the age variable is the age of the individual and the cohort variable is the date of birth of the individual.

The two time-dependent variables are period and age. In all cases, we have the linear relation:

period = cohort + age. (2.15)
This section introduces inference of the hazard rate, when one considers more than one of these three variables.

The Lexis Diagram. Consider Figure 2.1a, called the Lexis diagram. It represents the plane with the age in the y-axis and the period in the x-axis. The life of 2 individuals is represented by the black solid lines. These lines are oblique at a 45-degree angle: at the increment of each year, the age of every person increments by one year. In the Lexis diagram, the time origin is the birth of the individual and the lines represent the time spent at risk. The event of interest is represented by the dot at the end of each line. There is another variable that is hidden in these two variables, the date of birth of each individual, or cohort, which is the x-intercept of the lines. Figure 2.1b represents the age-cohort plane: the only time-dependent variable is age, on the yaxis. The cohort is on the x-axis, so that the life of the individuals are vertical lines, growing upwards as time passes. The same two individuals represented in Figure 2.1a are also represented here.

The age-cohort plane may be easier to represent the evolution of the time at risk, because they are represented by vertical lines. But the Age-period plane is a more natural presentation of the problem, because the period appears directly in most time-to-event data bases. Indeed, cohort studies have a starting date and an end date, and consequently the times for the event are comprised between two limiting periods. Between these two periods, events are registered for individuals with any age, and any date of birth (i.e. cohort). Consequently, cohort study data are represented as events comprised in a rectangle in the age-period plane.

Inference of the hazard rate. We present here the different approaches taken to infer the hazard rate over the Lexis diagram.

Age-period-cohort analysis can be interested in estimating the hazard as a bivariate function, defined in the following way:

λ (t |u) = lim δt →0 P (t ≤ T < t + δt |T ≥ t ,U = u) δt , (2.16)
where U is the cohort variable. An important approach is that of the factor models. The factor models assume a discretization of the variables: suppose the age and cohort variables are discretized into J and K intervals respectively. Then the period variable is naturally divided into J +K intervals (see Figure 2.1b, where J +K = 10 and we count 20 period intervals). Consider the discretized hazard rate (λ j ,k ) j ,k obtained by discretizing (2.16). The factors model assume that the age, cohort, and period variables have an additive effect on log(λ j ,k ): log(λ j ,k ) = α j + β k + γ j +k , (2.17)

where the parameter vectors α, β k , and γ j +k are the effects of each variable. Equation (2.17) is called the age-period-cohort model. Because of the linear dependency between the three variables (Equation 2.15), this model is not identifiable [START_REF] Clayton | Models for Temporal Variation in Cancer Rates. II: Age-period-cohort models[END_REF] and arbitrary constraints need to be added. We also define three submodels of (2.17): the age-period (log λ j ,k = α j + γ j +k ), agecohort (log λ j ,k = α j + β k ), and period-cohort (log λ j ,k = β k + γ j +k ) models, which are identifiable. The factors models denote the age-period-cohort model and its submodels. They play an important role in epidemiology (see Section 3.1) but the additive effect that they assume on the variables can be too restrictive.

In this chapter, we focus on developing a regularized estimation of the discretized hazard. In the next chapter, we develop an extension of the factor models that uses regularization to infer the effect of two variables and the interaction between the two effects. Both chapters rely on the application of the adaptive ridge to the discretized hazard rate.

Regularized estimation of the hazard rate

In this section, I include the preprint entitled "Regularized Bidimensional Estimation of the Hazard Rate". This work was conducted with my PhD supervisors and with the help of Jean-Christophe Thalabard.

Introduction

In epidemiological or demographic studies, with variable age at onset, a typical quantity of interest is the incidence or the hazard rate of a disease (for example the cancer incidence). In these studies, individuals are recruited and followed-up during a long period of time, usually from birth. The data are then reported either in the form of registers, which contain the number of observed cases and the number of individuals at risk to contract the disease, or in the form of the observed time for each individual. These types of studies are of great interest for the statistician, especially when the event of interest will tend to occur at late ages, such as in cancer studies. However, these data are usually highly heterogeneous in terms of dates of birth and with respect to the calendar time. In such cases, it is therefore very important to take into account the variability of the age, the cohort (date of birth) and the period (the calendar time) in the hazard rate estimation. This is usually done using age-period-cohort estimation methods (see Yang and Land, 2013, and citations therein).

In age-period-cohort analysis, the effects of age, period and cohort are fit as factor variables in a regression model where the output is the logarithm of the hazard rate. However, this induces an identifiability problem due to the relationship: period = age + cohort. There have been several solutions proposed to this problem. [START_REF] Osmond | Age, Period and Cohort Models Applied to Cancer Mortality Rates[END_REF] proposed to compute each submodel (age-cohort, age-period, and period-cohort) and use a weighting procedure to combine the three models. Different constraints have also been proposed to make the age-period-cohort model identifiable. However, as noticed by Heuer (1997, p 162), the obtained estimates highly depend on the choice of the constraints. [START_REF] Holford | The Estimation of Age, Period and Cohort Effects for Vital Rates[END_REF] proposed to directly estimate the linear trends of each effect. This procedure leads to results that are difficult to interpret. See Carstensen (2007) for a detailed discussion of the identifiability problem of the age-period-cohort model. More recently, Kuang et al. (2008) proposed to estimate the second order derivatives of the three effects. This model is implemented in the package apc [START_REF] Nielsen | Apc: An R Package for Age-Period-Cohort Analysis[END_REF]. Finally, Carstensen (2007) proposed to first fit one submodel (say age-cohort) and then to fit the period effect over the residuals of the first model. This model is implemented in the R package Epi [START_REF] Carstensen | Epi: A Package for Statistical Analysis in Epidemiology[END_REF], [START_REF] Plummer | Lexis: An R Class for Epidemiological Studies with Long-Term Follow-Up[END_REF].

All these approaches can be viewed as parametric models, where the parameters are the age, period, and cohort vector parameters. As such they are also restrictive because they do not allow for interactions between the three effects, that is they assume that one effect does not depend on the other effect's value. A different approach consists in considering the hazard rate as a function of age and either period or cohort and to estimate this bi-dimensional function in a non-parametric setting. No specific structure of the hazard rate is assumed. However, for moderate sample sizes, non-parametric approaches are prone to overparametrization. As a consequence, regularized methods have been proposed in order to avoid overfitting in this non-parametric context. A kernel-type estimator was proposed by [START_REF] Beran | Nonparametric Regression with Randomly Censored Survival Data[END_REF] and [START_REF] Mckeague | Identifying Nonlinear Covariate Effects in Semimartingale Regression Models[END_REF] where the cumulative hazard is smoothed using a kernel function. See [START_REF] Keiding | Statistical inference in the Lexis diagram[END_REF] for a thorough discussion of methods for hazard inference in age-period-cohort analysis. More recently, [START_REF] Currie | Smoothing Age-Period-Cohort Models with P -splines: A Mixed Model Approach[END_REF] proposed a spline estimation procedure to infer the hazard rate as a function of two variables. The authors use a generalized linear model using B-splines and overfitting is dealt with using a penalization over the differences of adjacent splines' coefficients.

In this article, we propose a new non-parametric method for bi-dimensional hazard rate estimation. As the previous non-parametric approaches, this model considers the estimation of the hazard rate with respect to two variables, i.e. either age-cohort, age-period, or period-cohort, without assuming any specific structure on the hazard rate. Inference is made in two dimensions, but through the linear relationship period = age + cohort, the hazard rate can be represented as a function of any two of the three variables. Finally, in order to take into account the issue of overfitting, we use the L 0 penalization procedure introduced by [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF], [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF], and Bouaziz and Nuel (2017). This penalty offers a segmentation of the hazard rate into constant areas. It makes use of an approximation of the L 0 norm which is computationally tractable. The novelty of this method lies in the parsimonious representation of the bi-dimensional hazard rate into segmented areas. In particular, the method can efficiently exhibit cohort, age or period effects, that is, specific changes of the hazard rate due to the date of birth, the age or the calendar time. Our approach also allows L 2 norm penalization, which will induce a smoothed estimate of the hazard in a similar way as the aforementioned non-parametric methods.

Our model is introduced in Section 1. The regularization method is then presented in Section 2. In Section 3, the penalty term selection problem is discussed. Finally, the performance of our model is assessed through a simulation study in Section 4 and illustrations on the SEER cancer dataset is provided in Section 5.

Modeling strategy

In the age-period-cohort setting, the date of birth (the cohort) U of each individual is available and the variable of interest is a time-to-event variable of this individual denoted T . The data are subject to right-censoring and they are represented as tabulated data over the J cohort intervals and the K age intervals

[c 0 , c 1 ), [c 1 , c 2 ), . . . , [c J-1 , c J ) and [d 0 , d 1 ), [d 1 , d 2 ), . . . , [d K-1 , d K ) respectively, with the convention c 0 = d 0 = 0 and c J = d K = ∞.
On a sample of n individuals, the available data can then be rewritten in terms of the exhaustive statistics O = (O 1,1 , . . . , O J,K ), R = (R 1,1 , . . . , R J,K ), where for j = 1, . . . , J, k = 1, . . . , K, O j,k represents the number of observed events that occurred in the j-th cohort interval [c j-1 , c j ) and k-th age interval [d k-1 , d k ) and R j,k represents the total times individuals were at risk in this j-th cohort and k-th age interval. In the case of register data, the discretization (c j ) , (d k ) is imposed by the data and the available data is directly R and O, which are often called the cases and person-years, respectively. See for instance Carstensen (2007) for an example of such data. The aim is to use the available data to provide an estimator of the hazard rate, defined in the age-cohort setting as: 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 in the situation where λ(t, u) is assumed to be piecewise constant. That is, we assume that

λ(t|u) = lim dt→0 1 dt P(t < T < t + dt|T > t, U = u), period age
λ(t|u) = J j=1 K k=1 λ j,k 1 [c j-1 ,c j )×[d k-1 ,d k ) (t, u),
and inference is made over the J × K dimension parameter λ = (λ 1,1 , . . . , λ J,K ). Note that the hazard can be equivalently defined as a function of age and period or as a function of period and cohort where the period is defined as the calendar time, that is: period = cohort + age. For illustration, the change of coordinates between the age-period and age-cohort diagrams is represented in Figure 1. In our models, the hazard will be considered as a function of solely age and cohort since the influence of any of the two elements of age, period or cohort can be retrieved using this reparametrization. Following Aalen et al. (2008, p. 224) the negative log-likelihood takes the form

n (λ) = J j=1 K k=1 {λ j,k R j,k -O j,k log (λ j,k )}. (1)
The authors also noticed that this log-likelihood is equivalent to a log-likelihood arising from a Poisson model. However, note that no distribution assumptions are made on the data and in particular the O j,k are not assumed to be Poisson distributed (see Carstensen, 2007, for a discussion on the "Poisson" model). Minimizing n yields an explicit maximum likelihood estimator λ mle j,k = O j,k /R j,k . However, for moderate sample sizes this estimator is overfitted, especially in places of the age-cohort plane where few events are recorded.

To remedy this problem we propose in the following to penalize the differences between adjacent values of the hazard in the log-likelihood.

For computation convenience, we first reparametrize the model: η j,k = log λ j,k , for 1 ≤ j ≤ J and 1 ≤ k ≤ K. The estimate is obtained by minimizing the penalized function

κ n (η, v, w) = n (η) + κ 2 J-1 j=1 K k=1 v j,k (η j+1,k -η j,k ) 2 + κ 2 J j=1 K-1 k=1 w j,k (η j,k+1 -η j,k ) 2 , (2)
where n (η) was defined in (1), κ is a penalty constant used as a tuning parameter, and v = (v 1,1 , . . . , v J-1,K ), w = (w 1,1 , . . . , w J,K-1 ) are constant positive weights of respective dimensions (J -1)K and J(K -1). Note that the case κ = 0 corresponds to the maximum likelihood estimation and the case κ = ∞ corresponds to a hazard uniformly constant over the age and cohort intervals. The parameter κ needs to be chosen in an appropriate way in order to obtain a compromise between these two extreme situations. This model does not attempt to estimate the age, period and cohort effect as parameter vectors. Instead, it performs a regularized estimation of λ that has no age-period-cohorttype structure. Two choices for the weights v and w can be made: one will lead to a smooth hazard rate and the other to a segmented hazard rate. This will be discussed in the next section. The choice of the optimal value for κ is addressed in Section 3.

Minimization of κ n is performed using the Newton-Raphson algorithm (see Algorithm 1). Let U κ n (η, v, w) = ∂ κ n /∂η be the gradient vector of the negative log-likelihood and I κ n (η, v, w) = ∂U κ n (η, v, w)/∂η T be its Hessian matrix. For 1 ≤ j, j ≤ J and 1 ≤ k, k ≤ K, we have

∂ n ∂η j,k (η) = exp (η j,k ) R j,k -O j,k , ∂ 2 n (η) ∂η j ,k ∂η j,k = 1 j=j ,k=k exp (η j,k ) R j,k , and
∂ κ n ∂η j,k (η) = ∂ n (η) ∂η j,k + κ [-v j,k (η j+1,k -η j,k ) + v j-1,k (η j,k -η j-1,k )] + κ [-w j,k (η j,k+1 -η j,k ) + w j,k-1 (η j,k -η j,k-1 )] , ∂ 2 κ n (η) ∂η j ,k ∂η j,k = ∂ 2 n ∂η j ,k ∂η j,k (η) + κ [1 j=j ,k=k (v j ,k + v j -1,k + w j ,k + w j ,k -1 ) -v j ,k 1 j=j +1,k=k -v j -1,k 1 j=j -1,k=k -w j ,k 1 j=j ,k=k +1 -w j ,k -1 1 j=j ,k=k -1 ] .
As a consequence, the Hessian matrix can be written

I κ n (η, v, w) = ∂ 2 n (η) ∂η∂η T + κB(η),
where B(η) is a band matrix of bandwidth equal to min(J, K) -1. Thus the Hessian matrix has the same structure as B(η) and the calculation of I κ n (η, v, w) -1 U κ n (η, v, w) has a O(min(J, K)JK) complexity instead of O(J 3 K 3 ). Fast inversion of the Hessian matrix is done using Cholesky decomposition as implemented in Rcpp in the package bandsolve 1 .

Choice of the regularization parameters v and w

In this section, two different expressions of the weights v and w are proposed which correspond to two different types of regularization of the hazard rate. The first one yields a smooth estimate. The second one uses an iterated adaptation of the weights to approximate an L 0 norm penalization of the first order differences.

Algorithm 1 Newton-Raphson Procedure with constant weights

1: function NEWTON-RAPHSON(O, R, κ, v, w) 2: η ← 0 3:
while not converge do 4:

η new ← η -I κ n (η, v, w) -1 U κ n (η, v, w) 5: η ← η new 6:
end while 7:

return η 8: end function

L 2 Norm Regularization

A ridge-type penalization is performed when setting v = w = 1. In this case the penalization corresponds to the square of the first-order differences of δ. In the penalized estimation model, this choice of weights yields a globally smooth estimator of the hazard rate. Note that our penalized maximum likelihood model will yield similar results as the spline method of [START_REF] Ogata | Likelihood Analysis of Spatial in Homogeneity for Marked Point Patterns[END_REF] presented in Section 1. In our method the penalization is performed over the first order differences of the parameter while in the spline method it is performed over the second order differences. This means that for arbitrarily large values of the penalty constant, the regularized hazard will be a constant function instead of a linear function. This model will be referred to as L 2 regularized estimation or smooth estimation.

Finally, one notes that Equation 2 allows for some flexibility in the regularization. Indeed, manually setting the weights v and w will allow to tune the importance of the regularization between different regions of the plane and between the two variables.

Approximate L 0 Norm Regularization

Following the work from [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF], [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF], and Bouaziz and Nuel (2017) an adaptive ridge procedure is performed when the weights are updated at each iteration of the Newton-Raphson algorithm. At the m-th iteration of the Newton-Raphson algorithm the weights are computed from the following formulas:

       v (m) j,k = η (m) j+1,k -η (m) j,k 2 + ε 2 v -1 , w (m) j,k = η (m) j,k -η (m) j,k-1 2 + ε 2 w -1
, where ε v and ε w are constants negligible compared to 1 (in practice one typically chooses

ε v = ε w = 10 -5
). We iterate between minimizing κ n for fixed weights and reevaluat-cohort age ing the weights such that at the m-th step, v

j = 1 j = 2 j = 3 k = 1 k = 2 k = 3 cohort age j = 1 j = 2 j = 3 k = 1 k = 2 k = 3 cohort age j = 1 j = 2 j = 3 k = 1 k = 2 k = 3
j+1,k -η (m) j,k ) 2 η (m) j+1,k -η (m) j,k 0 and w (m) j,k (η (m) j,k+1 -η (m) j,k ) 2 η (m) j,k+1 -η (m) (m) j,k (η (m)
j,k 0 , where . 0 denotes the L 0 norm -i.e. u 0 = 0 if u = 0 and u 0 = 1 otherwise. In other words, this adaptive ridge procedure approximates the L 0 norm regularization over the differences of η j,k and yields a segmentation of η j,k into piecewise constant areas. As with other classical penalized methods (e.g. LASSO, ridge) and as pointed out in [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF], the adaptive ridge penalization scheme induces a shrinkage bias. Therefore, after segmentation of the η i,j s, the hazard rate is estimated on each constant area using the unpenalized maximum likelihood estimator. More precisely, at convergence of the adaptive ridge algorithm, v j,k (η j+1,kη j,k ) 2 will be approximately equal to 0 if |η j+1,kη j,k | is smaller than ε v and approximately equal to 1 if |η j+1,kη j,k | is greater than ε v -and similarly for w j,k (η j,k+1η j,k ) 2 . Then one creates the graph whose vertices are the JK discretization cells and whose edges are the connexions between adjacent cells that have differences close to 0. Each connected component of this graph is a different area over which the hazard has been estimated to be constant. The extraction of connected components from the graph is done using the package igraph [START_REF] Csardi | The igraph Software Package for Complex Network Research[END_REF]. The log-hazard η (r) of the r-th constant area is such that r) . Finally, the values of η (r) are not estimated using the results of the adaptive ridge algorithm, but by unpenalized maximum likelihood estimation: η (r) = log O (r) /R (r) where O (r) is the number of events in the r-th constant area and R (r) is the time at risk in the r-th constant area.

∀[c j-1 , c j )×[d k-1 , d k ) ∈ r, η j,k = η (
This estimation method will be called L 0 regularized estimation or segmented estima-tion. This method is illustrated through the toy-example of Figure 2 and the adaptive ridge procedure is summarized in Algorithm 2. In practice, the stopping criterion for the adaptive ridge algorithm is when the absolute difference between successive values of the weighted differences is smaller than a predefined value -we use 10 -8 in our implementation.

Algorithm 2 Adaptive Ridge Procedure

1: function ADAPTIVE-RIDGE(O, R, κ) 2: η ← 0 3: v ← 1 4:
w ← 1 5:

while not converge do 6:

η new ← NEWTON-RAPHSON(O, R, κ, v, w) 7: v new j,k ← η new j+1,k -η new j,k 2 + ε 2 v -1 8: w new j,k ← η new j,k -η new j,k-1 2 + ε 2 w -1 9:
η ← η new 10:

end while 11:

Compute (O new , R new ) for selected (η, v new , w new ) 12: η new ← log (O new /R new ) 13:
return η new 14: end function 3 Choice of the penalty constant κ

In practice, the hazard rate needs to be estimated for a set of penalty constants and the choice of κ is determined as the penalty that provides the best compromise between model fit and reduced variability of the hazard rate estimate. For the L 0 regularization model, different values of the penalty constant lead to different segmentations of the η j,k . As a consequence, the problem of choosing the optimal penalty constant can be rephrased as the problem of choosing the optimal model among a set of models M 1 , . . . , M M , where each of these models corresponds to a different segmentation of the η j,k and M is the maximum number of different models. In this section we propose different methods to select the optimal model. Comparison of the efficiency of the different methods will be analyzed in Section 4 on simulated data.

We recall that R and O are the exhaustive statistics and η is the parameter to be estimated in our two models. Bayesian criteria attempt to maximize the posterior probability P(M m |R, O) ∝ P(R, O|M m )π(M m ), where P(R, O|M m ) is the integrated likelihood and π (M m ) is the prior distribution on the model. This problem is equivalent to minimizing -2 log P(M m |R, O). By integration

P(R, O|M m ) = η P(R, O|M m , η)π(η)dη,
where P (R, O|M m , η) is the likelihood and π(η) is the prior distribution of the parameter, which is taken constant in the following. Thus Bayesian criteria are defined as

-2 log (P(M m |R, O) = 2 n ( η m ) + q m log n -2 log π(M m ) + O P (1),
where q m is the dimension of the model M m i.e., the number of constant areas selected by the adaptive ridge algorithm.

The BIC (Schwarz, 1978) corresponds to the Bayesian criterion obtained when one neglects the term π(M m ), which is equivalent to having a uniform prior on the model:

BIC(m) = 2 n ( η m ) + q m log n. (3)
As explained by Żak-Szatkowska and Bogdan ( 2011), a uniform prior on the model is equivalent to a binomial prior on the model dimension B(JK, 1/2). When the true model's dimension is much smaller than the maximum possible dimension JK, the BIC tends to give too much importance to models of dimensions around JK/2, which will result in underpenalized estimators. To this effect, [START_REF] Chen | Extended Bayesian information criteria for model selection with large model spaces[END_REF] have developed an extended Bayesian information criterion called EBIC 0 (or EBIC for short). One can write

π(M m ) = P(M m |M m ∈ M [qm] )P(M m ∈ M [qm] ) where M [qm]
is the set of models of dimension q m . The EBIC 0 criterion is defined by setting

P(M m |M m ∈ M [qm] ) = 1/ JK qm and P(M m ∈ M [qm] ) = 1. Thus π(M m ) = JK q m and EBIC 0 (m) = 2 n ( η m ) + q m log n + 2 log JK q m . (4)
Note that the EBIC 0 assigns the same a priori probability to all models of same dimension. Therefore, when the true model's dimension is not close to JK/2 the EBIC 0 will be able to select this model more easily. Namely, when the true model's dimension is very small the EBIC 0 will tend to choose very sparse models.

The last criterion that will be used is the Akaike Information Criterion [START_REF] Akaike | Information Theory and an Extension of the Maximum Likelihood Principle[END_REF], or AIC, defined as AIC(m) = 2 n ( η m ) + 2q m . This criterion is known for performing better than the BIC in terms of mean squared error, however the BIC will tend to select sparser models than the AIC.

Note that Bayesian criteria and the AIC can only be used for the L 0 regularized estimation only, since the L 2 model does not perform a model selection. An alternative to performing model selection is to use the K-fold cross validation. With this method, the data are split at random into L parts. The estimated parameter obtained when the l-th part is left out is noted η -l (κ) and the cross-validated score is defined as

CV(κ) = L l=1 κ,l n ( η -l ),
where κ,l n is the negative log-likelihood evaluated on the l-th part of the data. The optimal penalty constant is obtained by minimizing CV(κ) with respect to κ. The L-fold cross validation method can be used for both the L 0 regularized estimation and the L 2 regularized estimation. However, this method is numerically time consuming as the estimator has to be computed L times while Bayesian criteria or the AIC provide direct methods to perform model selection from the original estimator. In the simulation studies and data analysis, we set L = 10.

Simulation study 4.1 Simulation designs

In this section, our segmented estimation method with L 0 norm is compared with the AGE-COHORT model and with the smoothed hazard estimate with the L 2 norm. The different criteria for model selection are also compared with each other. We present two simulation designs. In the first one, the true hazard rate is generated from a smooth age-cohort model which includes an interaction term on a small region of the age-cohort plane. In the second case, the true hazard rate is a piecewise constant function with four heterogeneous areas. The two true hazards are displayed in Figure 3, both in greyscale and in perspective plot.

The simulation design is as follows. We set J = 10 equally spaced age intervals and K = 10 equally spaced cohort intervals. The age intervals are defined as [0, 10), . . ., [90,100] and the cohorts intervals are defined as [1900, 1910) ,. . ., [1990, 2000]. In order to simulate a dataset, the cohorts are first sampled on K = 10 cohort group intervals of 10 years length ranging from 1900 to 2000. Censoring is then simulated as a uniform distribution over the age interval [75,100] for all cohorts such that all observed events are comprised in the age interval [0, 100]. Since in practice one does not know the appropriate discretization in advance, a different discretization was used for the estimation procedure : the age and cohort intervals were defined as 5-year length intervals instead of 10 for the true hazard. As a result, a total of 20 × 20 parameters need to be estimated.

For each of the two designs, we simulated data of sample sizes 100, 400, 1000, 4000, and 10000. For each sample size, the simulation and estimation were replicated 500 times.

Smooth true hazard

The smooth true hazard (Figures 3a and3b) is generated using the age-cohort model log λ j,k = µ + α j + β k with an intercept µ = log(10 -2 ). The age effect vector α and cohort effect vector β are arithmetic sequences such that α 2 = 0, α J = 2.5, β 2 = 0, and β K = 0.3. An interaction term is added to the hazard. It corresponds to a bump in the hazard located in the neighbourhood of the region of the age-cohort plane (45,1945). The bump is defined as 10 times the Gaussian density function with mean (1945, 45) and with a diagonal variance-covariance matrix with diagonal equal to (50, 50). This true hazard displays a sharp increase for high values of the age, which implies that few events will be recorded in this region. On average, 91 % of the events are observed in this simulation design.

Piecewise constant true hazard The piecewise constant true hazard (Figures 3c and3d) has four constant areas over the age-cohort square [0, 100] × [1900,2000]. On average, 71 % of the events are observed in this simulation design.

Performance of the estimation methods in terms of MSE

Our two estimation methods (L 0 and L 2 norm) are compared in terms of the Mean Squared Errors in each simulation scenario. The different selection methods for the penalty (AIC, BIC, EBIC and cross-validation) are also compared. The results are presented in Table 1.

On the overall, the EBIC and cross-validated criteria outperform the AIC and the BIC for the two simulations scenarios. This is particularly true for small sizes where the AIC and the BIC behave very poorly. As expected, the L 2 norm estimator is the most performant of all estimators in the smooth true hazard scenario (Table 1a) and the L 0 method performs better in the piecewise constant hazard scenario (Table 1b) than in the smooth true hazard scenario. The L 2 norm estimator is also the most performant of all estimators in the piecewise constant hazard scenario except for very large sample sizes (n = 10000) where the BIC, EBIC and cross-validated criterion provide slightly better performances. Finally, in both scenarios, the EBIC always outperforms the AIC, the BIC and the cross-validated criterion. Different censoring rates were also studied which showed a degradation of the performances of the overall estimators as the percentage of censoring increased. The performance in terms of number of selected areas was also investigated. It showed that the EBIC and CV criterion perform better at selecting sparse models with few areas, while the AIC and BIC tend to overestimate the true number of areas. Indeed, for sample size 4000, the 80% inter-quantile range of the selected number of areas is [START_REF]4 Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ+α j +β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j[END_REF]5] for the EBIC and [1,5] for the CV, whereas it is [START_REF]4 Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ+α j +β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j[END_REF]13] and [36,72] for the BIC and AIC respectively. These experiments are not reported here.

In conclusion, the simulation experiments suggest to use the EBIC among all different criteria for the L 0 norm estimator as it provides the best tradeoff between computation time and estimation performance.

Perspective plots of the estimation methods

In this section the performance of our adaptive ridge (L 0 norm) and ridge (L 2 norm) estimates is assessed visually by comparison of the true hazard. The standard age-cohort model [START_REF] Holford | The Estimation of Age, Period and Cohort Effects for Vital Rates[END_REF] has also been implemented. This model assumes that the hazard has the following expression:

log λ j,k = µ + α j + β k ,
where µ is the intercept, α is the age effect and β is the cohort effect. It should be noted that this model does not allow for interactions between age and cohort effects. Perspective plots of the median hazard estimations over 500 replications are presented in Figures 4 and5 for the smooth and piecewise constant true hazard respectively. For the L 0 regularized estimate, the penalty constant is chosen using the EBIC. In Figure 4, it is seen that the age-cohort model is not able to estimate the central bump in the hazard. On the contrary, the smoothed estimate accurately recovers the shape of the true hazard except for the high values of age where few events are observed. Interestingly, it is seen that our segmentation method provides similar results as the smoothing technique even though the true hazard is not piecewise constant.

The results in Figure 5 yield similar conclusions. The age-cohort model behaves very poorly due to its constrained structure while the ridge and adaptive estimates provide satisfactory results. In particular the shape of the true hazard is correctly captured by the adaptive ridge on the majority of replicated samples.

Real data application

Our method is applied to data of survival times after diagnosis of breast cancer. The dataset is provided by the Surveillance, Epidemiology, and End Results (SEER) Program from the US National Cancer Institute (NCI). SEER collects medical data of cancers (including stage of cancer at diagnosis and the type of tumor) and follow-up data of patients in the form of a registry. Around 28 percent of the US population is covered by the program. The registry started in February 1973 and the available current dataset includes follow-up data until January 2015. We refer to the website https://seer.cancer.gov/ for information about the SEER Program and its publicly available cancer data.

In this study the duration of interest T is the time from breast cancer diagnosis to death in years, the variable U is the date of diagnosis (in years) and the period is the calendar time (in years). Patients continuously entered the study between 1973 and 2015 and rightcensoring occurred for patients that were still alive at the end of follow-up or for those that were lost to follow-up.

The breast cancer data was extracted using the package SEERaBomb. For the sake of comparison, the subsample of malignant, non-bilateral breast tumor cancers was extracted from the dataset, such that the data comprises 1, 265, 277 women with 60 percent of censored individuals. Times from diagnosis to last day of follow-up vary between 0 and 41 years, and the dates of cancer diagnosis U i vary between 1973 and 2015. Death from another cause than cancer is available in the dataset and is accounted for as right-censoring.

The implementation of our adaptive ridge method aims at two goals. Firstly we aim at simultaneously detecting a cohort effect and an age effect, that is the evolution of the mortality with respect to the time elapsed since cancer diagnosis (age effect) and with respect to the date of diagnosis (cohort effect). Secondly, our method will provide estimation of the hazard rates on the resulting heterogeneous areas. The method is first applied on the whole sample of 1265277 individuals. In order to take into account the fact that mortality from cancer highly depends on the cancer stage, we also perform a stratified analysis with respect to the stage of cancer at diagnosis. For this purpose, we use the cancer stage classification provided by the SEER data: we keep the patients with cancer stages 1, 2, and 3 at the time of diagnosis. This classification closely follows that of the American Joint Comitee on Cancer (AJCC), 3 rd Edition; the details are given at page 86 of the manual entitled Comparative Staging Guide for Cancer, available at https://seer.cancer.gov. The main difference between the two classifications is that the SEER Program classifies the cases where lymph node status cannot be assessed as if there was no regional lymph node metastasis. The L 0 estimates for the whole sample and for each cancer stage are displayed in Figure 6. We see that the different stages of cancer at diagnosis have a great impact on the survival times. For Stage 1 cancers, the mortality is low between 0 and 4-5 years after diagnosis, and steadily increases afterwards. The date of diagnosis seems to have no impact on the mortality of Stage 1 cancers. On the other hand, Stage 2 cancers exhibit a strong effect of the date of diagnosis: around 1995 -1997, the mortality significantly decreases. This can correspond to an improvement of the treatment of breast cancer around that period in the United States. Finally, Stage 3 cancers display a very high hazard rate across all dates of diagnosis. This seems to indicate that the evolution in treatments of breast cancer had a significant impact on the survival times after diagnosis, but almost exclusively when cancers were diagnosed at Stage 2. Two additional analyses of the hazard rate with stratification with respect to age at diagnosis and estrogen receptor status were performed in the Supplementary Materials. The results suggest that the shift in mortality around year 1996 could correspond to the introduction of hormone-blocking therapy [START_REF] Fisher | Tamoxifen in treatment of intraductal breast cancer: National Surgical Adjuvant Breast and Bowel Project B-24 randomised controlled trial[END_REF].

Conclusion

In this article, we have introduced a new estimation method to deal with age-period-cohort analysis. This model assumes no specific structure of the effects of age and cohort and the hazard rate is directly estimated without estimating the effects. In order to take into account possible overfitting issues, a penalty is used on the likelihood to enforce similar consecutive values of the hazard to be equal. Two different types of penalty terms were introduced. One leads to a ridge type regularization while the other leads to a L 0 regularization. Different selection methods of the penalty parameter were also introduced. To our knowledge, a segmented estimation model of this kind has never been introduced in this context.

Using simulated data, it has been shown that the cross validated ridge estimator and the EBIC 0 adaptive ridge estimator perform the best in terms of mean squared error. The cross validation criterion was shown to provide the best fit of the hazard rate, but its very high computationally cost makes it non-competitive. In this context, this modified BIC criterion comes out as a powerful tool to select the best bias-variance tradeoff.

The method was successfully applied to data of survival after breast cancer provided by the SEER program. The segmented estimate of the hazard rate displays important information about the shift in mortality after being diagnosed of breast cancer in the United States in the mid-1990s.

Our method could be directly extended to a different discretization of the age-periodcohort plane, such as 1 × 1 × 1-year triangles that are represented in dark gray in Figure 1 (see Section 3 of Carstensen, 2007, for an example of this discretization). Another extension would be to consider other types of penalizations. Instead of estimating a piecewise constant hazard, one could estimate a piecewise linear hazard by penalizing over second order differences of the hazard.

Application to the evolution of breast cancer mortality

In this section, I include the supplementary material to the paper "Regularized Bidimensional Estimation of the Hazard Rate". The mortality of breast cancer is known to greatly vary on whether the cancer is pre or postmenopausal (Consensus, 1985). Consequently, a thorough analysis of the mortality from breast cancer would require to stratify with respect to the menopausal status at diagnosis. Since this covariate is not present in the data, we decided to stratify the sample with respect to the age of the patient at diagnosis, which is a proxy of menopausal status. Most women are known to have their menopause between 45 and 55 years old [START_REF] Hill | The demography of menopause[END_REF][START_REF] Henderson | Predictors of the timing of natural menopause in the Multiethnic Cohort Study[END_REF][START_REF] Gold | The timing of the age at which natural menopause occurs[END_REF], with 25th, 50th, 75th percentiles ranging from years 47-49, 50-51, 52-54, respectively, according to countries and surveys [START_REF] Mishra | Early menarche, nulliparity and the risk for premature and early natural menopause[END_REF]. Consequently, based on the available information in SEER, for each cancer stage, the patients were divided into three classes of age at diagnosis: (., 45], (45, 55], and (55, .) as a proxy for pre-menopausal, peri-menopausal and post-menopausal ages, respectively. The resulting estimated hazards are represented in Figure 1.

Supplementary

The stage 1 cancer patients younger than 45 and the stage 3 cancer patients older than 55 display the same mortality across all dates of diagnosis, i.e. with no cohort effect.

Moreover, the mortality of stage 1 cancer patients aged 45 and older at diagnosis has a slight cohort effect corresponding to a progressive decrease in the mortality across all survival times [START_REF] Peto | UK and USA breast cancer deaths down 25% in year 2000 at ages 20-69 years[END_REF]. This could suggest a trend of slow and steady improvement of the treatment of breast cancer in the United States over the period 1887 -2005.

Finally, we observe a clear decrease of the mortality for stage 2 cancers for all three age classes. This shift is located at the year 1995 for middle-aged patients and around the years 1997 -1998 for patients younger than 45 and older than 55. The same drop in mortality is observed for stage 3 cancers with patients younger than 45 at diagnosis, around year 1995. This could correspond to the introduction of improvements in the treatments of breast cancer in the United States (Consensus, 1985). Among the three main medical innovations, which can be considered in this period, the improvement of the surgical procedures for the loco-regional control of the disease and the assessment 1 of the beneficial effect of hormone-receptor therapies could be reflected in the observed survival in stages 1-2, whereas the later emergence during this period of new classs of chemotherapeutic agents like taxoids [START_REF] Rowinsky | Taxol: The first of the taxanes, an important new class of antitumor agents[END_REF][START_REF] Crown | Docetaxel and paclitaxel in the treatment of breast cancer: A review of clinical experience[END_REF] or herceptin-based therapies targeted on new class of tumor markers [START_REF] Pegram | HER-2/neu as a predictive marker of response to breast cancer therapy[END_REF][START_REF] Emens | Trastuzumab in breast cancer[END_REF]) would be related with the changes in survival observed in stage 3. In the next section, we will use a stratified analysis to understand the effect of hormone-receptor therapies on the mortality shift in the mid-1990s.

2 Application to Breast Cancer Mortality: Stratification with respect to the Estrogen Receptor Status

The cohort effect highlighted in the previous section could correspond to the introduction of Selective Estrogen Receptor Modulator (SERM) treatments and in particular the use of Tamoxifen as a treatment for breast cancer, showing improved survival in women with estrogen receptor positive tumor, initially in post-menopausal women [START_REF] Fisher | A randomized clinical trial evaluating tamoxifen in the treatment of patients with node-negative breast cancer who have estrogen-receptor-positive tumors[END_REF], later in both post -and premenopausal women (Early Breast Cancer Trialists' Collaborative Group, 1988;[START_REF] Fisher | Tamoxifen for Prevention of Breast Cancer: Report of the National Surgical Adjuvant Breast and Bowel Project P-1 Study[END_REF][START_REF] Pritchard | Endocrinology and hormone therapy in breast cancer: Endocrine therapy in premenopausal women[END_REF][START_REF] Cochrane | Tamoxifen for early breast cancer[END_REF]. Indeed, Tamoxifen was gradually used in the early years of 1990's [START_REF] Gail | Weighing the Risks and Benefits of Tamoxifen Treatment for Preventing Breast Cancer[END_REF][START_REF] Harlan | Adjuvant therapy for breast cancer: Practice patterns of community physicians[END_REF][START_REF] Mariotto | Dissemination of adjuvant multiagent chemotherapy and tamoxifen for breast cancer in the United States using estrogen receptor information: 1975-1999[END_REF] to decrease the mortality of breast cancer patients. This treatment is only efficient on estrogen receptor-sensitive cancers. To validate our hypothesis, we conducted the estimation of mortality separately for patients with estrogen receptor sensitive and non-sensitive cancers. Since stage 2 cancers displayed a strong cohort effect across all ages at diagnosis, we only kept stage 2 cancers in this study. The estimated mortality is given in Figure 2. Note that the spikes in the mortality are an artifact of the segmentation procedure when the sample sizes tend to be too small in some regions of the age-cohort plane and are not to be taken into account in the interpretation of the mortality.

There is a clear difference in the evolution of mortality with respect to time at diagnosis between sensitive and non-sensitive estrogen cancers. For estrogen sensitive cases, the mortality displays the same sudden decrease around years 1997 -1998 as in Figure 1, across all age classes. In particular for individuals aged 55 or more at the time of diagnosis, the mortality has gradually decreased for estrogen sensitive patients, whereas it did not evolve with time for estrogen non-sensitive patients. On the other hand, the mortality for non-estrogen sensitive cancers displays almost no cohort effect for all ages at diagnosis [START_REF] Knight | Estrogen receptor as an independent prognostic factor for early recurrence in breast cancer[END_REF].

The same analysis was run with stratification with respect to progesterone receptor status, with very similar morality estimates (results not shown here). Further analyses could be carried out to better understand the effect of the introduction of hormone-blocking therapies on mortality. However, the segmentation of the hazard rate, even with this simple stratified analysis, highlighted that the adoption of SERM therapies in the United States is a potential reason for the sharp decrease of mortality in the middle of the 1990s [START_REF] Peto | UK and USA breast cancer deaths down 25% in year 2000 at ages 20-69 years[END_REF].

Introduction to age-period-cohort analysis

In the study of diseases with age at onset, the incidence -or hazard rate -is an important quantity to estimate. In these studies, data are either present in the form of registers or collected from a cohort that is being followed. In both cases, the individuals are usually very heterogeneous in terms of date of birth -also called cohort. Equivalently, the event of interest is observed at widely different calendar times -also called period. It is therefore necessary to adjust for this heterogeneity. Namely, many models have been developed to infer the hazard rate as a function of age, cohort, and period, called age, cohort, and period effects. This is the goal of age-period-cohort analysis.

The simplest models in age-period-cohort analysis are the age-cohort and age-period models [START_REF] Clayton | Models for Temporal Variation in Cancer Rates. II: Age-period-cohort models[END_REF]. These models assume that each variable has an additive effect on the logarithm of the hazard. In the age-cohort model,

log λ j ,k = α j + β k , (3.1)
where α j is the age effect at age j , β k is the cohort effect at calendar time c, and λ is the hazard rate. In this model, the hazard rate is assumed to be the product of the age effect and the cohort effect. Since there are fewer parameters than values of the hazard, this model yields a regularized estimation of λ, and thus of the age and cohort effect. But this model is very restrictive because it assumes that the age and cohort have an additive effect only, that is, the log-hazard difference between two ages is the same for all cohorts and that the log-hazard difference between two cohorts is the same for all ages. The full model in age-period-cohort analysis is called the age-period-cohort model [START_REF] Clayton | Models for Temporal Variation in Cancer Rates. II: Age-period-cohort models[END_REF]. It jointly estimates the three additive effects:

log λ j ,k = α j + β k + γ j +k , (3.2)
where γ l is the period effect at calendar time l = j + k. Because of the linear dependency (age + cohort = period) between the variables, these effects are only identifiable up to a linear trend. This makes the estimated effects hard to interpret and many attempts have been made to overcome this difficulty. Carstensen (2007) offered to choose a cohort of reference and to infer the relative hazard with respect to this cohort. Kuang et al. (2008) introduced a reparametrization based on the second order derivatives of the effects which allows estimation of the three effects when the index set i , j is a trapezoid. Both the two-factor models and the age-period-cohort model are useful to infer the influence of the corresponding variables over the hazard rate. Indeed the age, cohort, and period effects are important information to understand the evolution of the incidence. However they can represent a simplistic view of the evolution of the hazard along time. We study the age-cohort model where an interaction term is added. This term is introduced in order to capture non-linear relations between age and cohort effects. The interaction between age and cohort can be interpreted as regions of the age-cohort plane where the hazard is either unexpectedly high or unexpectedly low. Thus, the interaction between effects is an important tool to detect peculiar features of the event of interest.

In this chapter, we introduce a generalization of the age-cohort model which includes the estimation of the interaction between the age and cohort effects. This model is called age-cohortinteraction model, or "ACI" model. For the sake of simplicity, we only consider the age-cohort model in this work, but the period-cohort and age-period models can be generalized in the same way to estimate interactions between effects. We introduce a penalized likelihood approach to jointly infer the age and cohort effects as well as the interaction effect. For the model to be interpretable, we enforce a regularization of the interaction term. We offer a choice between two types of regularization: one provides a smooth estimate of the interaction and one provides a segmented estimate of the interaction. We illustrate these two regularizations and their respective advantages.

The rest of this chapter is constructed as follows. Section 3.2 introduces the age-cohort-interaction model. Section 3.3 develops on a procedure for the estimation of the model. In Section 3.4, we explain how to choose the tuning parameter of the regularization. Finally, our method is illustrated on simulated data in Section 3.5.

The age-cohort-interaction model

We assume that the time of interest is subject to independent right-censoring (Fleming and Harrington, 2011, p. 27), such that the available data is of the form (U i , T i , ∆ i ), i = 1, . . . , n, where T i is the censored time-to-event of individual i and ∆ i equals 0 if individual i is censored and 1 otherwise. The hazard rate as a function of the cohort u and age t is defined by

λ(t |u) = lim δt →0 1 δt P(t < T < t + δt |T > t ,U = u),
where U is the random variable of the cohort. The cohort variable is discretized into J intervals [c 0 , c 1 ), . . . , [c J -1 , c J ) and the age variable is discretized into

K intervals [d 0 , d 1 ), . . . , [d K -1 , d K ).
The data is rewritten using the exhaustive statistics O = O 1,1 , . . . ,O J ,K and R = R 1,1 , . . . , R J ,K , where O j ,k and R j ,k are respectively the number of observed events and total time at risk in the j -th cohort interval and in the k-th age interval. We assume that λ

(t |u) is constant over each set [c j -1 , c j ) × [d k-1 , d k )
and infer the discretized hazard rate λ 1,1 , . . . , λ J ,K . The negative log-likelihood takes a simple expression:

n (λ) = j k λ j ,k R j ,k -O j ,k log λ j ,k . (3.3)
We introduce the age-cohort-interaction ("ACI") model:

log λ j ,k = µ + α j + β k + δ j ,k , (3.4)
where α j is the cohort effect for the j -th cohort interval, β k is the age effect for the k-th age interval, and δ j ,k is the interaction term for the j , k -th rectangle. The parameter µ plays the role of the intercept (that is λ 1,1 = µ), and we impose α 1 = β 1 = δ 1,k = δ j ,1 = 0 so that the model has J K freely varying parameters and is identifiable. The parameter of the model is θ = (µ, α, β, δ), where α = (α j ) 2≤ j ≤J , β = (β k ) 2≤k≤K , and δ = (δ j ,k ) 2≤ j ≤J ,2≤k≤K . Note that the interaction term δ between the age and cohort is different from the period effect γ in the age-period-cohort model (3.2). The period effect can be seen as an interaction effect between age and cohort, which is a function of the period only. Contrarily to the age-period-cohort model, our model can estimate an interaction between age and cohorts which is any function of these two variables. Consequently, the ACI model can be seen as an extension of the age-cohort model which is also more general than the age-period-cohort model.

The negative log-likelihood takes the form

n (θ) = J j =1 K k=1 exp µ + α j + β k + δ j ,k R k j -µ + α j + β k + δ j ,k O k j . (3.5)
For the parameter δ to correctly estimate the interaction between age and cohort -that is, the part of the hazard which diverges from the age-cohort model -one needs to impose a constraint on δ.

We want the δs to be zero-valued everywhere except for values of j , k for which the hazard is too different from the age-cohort model. To this end, the parameters are estimated using the penalized negative log-likelihood:

n (θ) - κ 2 j ,k { δ j ,k -δ j -1,k 0 + δ j ,k -δ j ,k-1 0 },
where x 0 is the L 0 norm, equal to 0 if x = 0 and to 1 otherwise and κ > 0 is a trade-off parameter. This penalty is proportional to the number of non-zero differences of δ j ,k . Consequently, the penalized estimation enforces the differences of δ to be equal to zeros except in a small amount of times.

Since δ 1,k = δ j ,1 = 0, this penalty also enforces the interaction term to be equal to zero except where necessary.

This penalty term ensures that δ's second order total variation is small. Thus the parameters µ, α, and β estimate the age-cohort model that is the closest to the data, and δ estimates the divergence of the data from the age-cohort model. The penalty constant tunes the trade-off between goodness-of-fit and regularization. When κ = 0, the ACI model is subject to no regularization and the corresponding estimate is the maximum likelihood estimate. When pen → ∞, the interaction term is set to be uniformly equal to zero and the ACI model is the age-cohort model. In this sense, the ACI model is an extension of the age-cohort model which allows for interaction between age and cohort effects. This penalty is however not continuous, which makes the maximization of the penalized loglikelihood intractable in practice. The next section deals with a numerical method for approximating the latter penalized likelihood.

The estimating procedure

Following the works from [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF] and [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF], we use an iterative algorithm to approximate the L 0 norm penalty. The estimate is the minimizer of:

κ n (θ) n (θ) + κ 2 j ,k v j ,k δ j ,k -δ j -1,k 2 + w j ,k δ j ,k -δ j ,k-1 2 , ( 3.6)
where v j ,k and w j ,k are positive weights. The procedure iterates between minimizing (3.6) with fixed weights and updating the weights, using the formula

       v (m) j ,k = δ (m) j +1,k -δ (m) j ,k 2 + ε 2 v -1 , 1 ≤ j ≤ J -1, 2 ≤ k ≤ K , w (m) j ,k = δ (m) j ,k -δ (m) j ,k-1 2 + ε 2 w -1 , 2 ≤ j ≤ J , 1 ≤ k ≤ K -1, (3.7)
where (m) is the iteration step, δ (m) j ,k is the parameter from the last iteration and ε v and ε w are positive constants negligible compared to 1 (in practice we choose ε v = ε w = 10 -6 [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF]).

At convergence, v j ,k δ j +1,k -δ j ,k 2 are very close to 1 if the adjacent values δ j +1,k and δ j ,k have been estimated to have different values and to 0 if they have been estimated to have the same valueand similarly for w j ,k δ j ,k+1 -δ j ,k 2 . We typically use a threshold of 10 -8 , so that values smaller than 10 -8 are set to 0 and values larger than 1-10 -8 are set to 1. Then one creates the graph whose vertices are the J K discretization rectancles and whose edges are the connexion between adjacent rectangles that have a difference close to 0. As with other classical penalized methods (e.g. lasso, ridge) and as pointed out in [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF], the adaptive ridge penalization scheme induces a shrinkage bias. Therefore, the unpenalized maximum likelihood estimate is used to infer the value of the hazard over the constant areas estimated using the adaptive ridge procedure. The values of δ (r ) are not estimated using the results of the adaptive ridge algorithm but by unpenalized maximum likelihood estimation: δ (r ) = O (r ) /R (r ) where O (r ) is the number of events in the r -th constant area and R (r ) is the time at risk in the r -th constant area.

The estimating procedure is summarized in Algorithm 5. See also (Bouaziz and Nuel, 2017) and [START_REF] Goepp | Regularized Bidimensional Estimation of the Hazard Rate[END_REF] for implementations of the adaptive ridge procedure in similar contexts. The penalized likelihood in Equation (3.6) is minimized using the Newton-Raphson algorithm. Expressions of first and second order derivatives of κ n are given in Appendix A. The Hessian matrix has J K rows and columns, so its inversion is computationally intensive. A fast inversion method for its inversion is detailed in Appendix B. while not converge do 6:

θ new ← NEWTON-RAPHSON(O, R, κ, v , w ) 7: v new j ,k ← δ new j +1,k -δ new j ,k 2 + ε 2 v -1 8: w new j ,k ← δ new j ,k -δ new j ,k-1 2 + ε 2 w -1 9: θ ← θ new 10:
end while 11:

Compute (O new , R new ) selected (δ, v new , w new ) 12: θ ← O new /R new 13:
return θ 14: end function

Choice of the penalty constant

In practice, the hazard rate needs to be estimated for a set of penalty constants and the optimal κ is chosen as the penalty that provides the best compromise between model fit and reduced variability of the hazard rate estimate. In the piecewise constant estimation, different values of the penalty constant yield to different segmentations of δ j ,k . As a consequence, to choose the optimal penalty constant is to choose the optimal model, where each model corresponds to a different segmentation of the δ j ,k s.

The computational burden of cross validation dissuades us from using it for the ACI model. We offer to use a criterion for model selection. The AIC [START_REF] Akaike | A New Look at the Statistical Model Identification[END_REF], the BIC (Schwarz, 1978), and, following a similar work in [START_REF] Goepp | Regularized Bidimensional Estimation of the Hazard Rate[END_REF], a modified version of the BIC called "EBIC" (for "Extended BIC", see [START_REF] Chen | EXTENDED BIC FOR SMALL-n-LARGE-P SPARSE GLM[END_REF]). The latter is defined as a BIC criterion with a specific prior on the model. It is defined as

EBIC(κ) = BIC(κ) + 2 log J K m(κ) , (3.8)
where m(κ) is the model dimension. Notice that the term on the right hand side of the latter equation takes small values when the model dimension is close to 1 or to J K , the maximal dimension. Thus the prior of the EBIC gives more weight to models of small and large dimensions, compared to J K /2. Comparison of the efficiency of the different methods will be analyzed in Section 3.5 on simulated data.

Simulation results

Simulation setting

The simulation setting is as follows. The data is simulated using a "true" (conditional) hazard rate λ(t |u), given as piecewise constance function: (3.9) where (λ j ,k ) is the discrete hazard rate. We set J = 20 equally spaced age intervals [0, 5), . . . , [95,100] and K = 20 equally spaced cohort intervals [1900, 1905), . . . , [1995, 2000] 1900,2000]. For inference, we use the same discretization into 5-year age and cohort intervals. Consequently, there are J × K = 400 parameters in the full model. We use three simulation designs, with three different hazard rates; see Figure 3.1.

λ(t |u) = J j =1 K k=1 λ j ,k 1 [c j -1 ,c j )×[d k-1 ,d k ) (t , u),
For each design, we simulate data of sample sizes 400, 100, 4000, and 10000. For each sample size, the simulation and estimation were replicated 100 times. In all simulation designs, the log-hazard rate is the sum of and age effect, a cohort effect, and an interaction effect:

log λ * j ,k = µ * + α * j + β * k + δ * j ,k ,
where µ * = log(10 -2 ) and (α * j ) and (β * k ) are sequences in arithmetic progression ranging from α * 2 = 0 to α * J = 1.4 and from β * 2 = 0 to β * K = 0.3 respectively. We defined three simulation designs, with different choices for the interaction effect δ * j ,k .

Simulation Design 1. In this simulation design δ * j ,k is the sampling of the mixture of two Gaussian densities:

δ * j ,k = 15 f µ 1 ,σ 1 (c j , d k ) + 3 f µ 2 ,σ 2 (c j , d k ),
where f µ,σ (t , u) denotes the bivariate Gaussian density at age t and cohort u. We set µ 1 = (40, 1940) T , µ 2 = (50, 1965), σ 1 = 100I 2 , and σ 2 = 10I 2 . The resulting hazard is represented in Figure 3.1a: the general (increasing) trend is given by the age and cohort effects, and the two "bumps" are given by the interaction effect. The first bump is spread out and has a large amplitude; the second bump however is smaller and more spiky. The latter can be difficult to estimate. Contrarily to Design 2, the interaction term has a low amplitude compared to the combined (additive) effect of age and cohort. Consequently, the interaction effect is purposefully hard to estimate in this design, as it is hard to separate it from the age and cohort effects.

Simulation Design 2. In this simulation design, the interaction term is a piecewise function whose support is a circle at the center of the age-cohort plane:

δ * j ,k = 0.01 * 1(( j , k) ∈ C (1950, 50, 500)).
where C m k , m j , R is the circle of center m k , m j and of radius R.

The interaction term is a very simple piecewise constant function. The resulting hazard is represented in Figure 3.1b. In this design, it is required that the ACI model infers correctly the support of δ * j ,k .

Simulation Design 3. This simulation design is piecewise constant with two levels: 1955,45,200) .

δ * j ,k = 0.015 * 1 ( j , k) ∈ C (1950, 50, 700) + 0.015 * 1 ( j , k) ∈ C (
The resulting hazard is represented in Figure 3.1c. The two terms in this equations are functions whose supports are circles and one circle is inside the other. Consequently the interaction effect has a large support and a large amplitude (compared to the age and cohort effects). Hence it is easy to estimate the presence of the interaction term. However it can be harder to estimate correctly the presence and shape of the higher level.

Predictive performance

We represent the estimation performance of the ACI model compared to two other estimates: the AC model and the MLE. Note that in all simulation designs, the interaction term takes values close to zero and consequently the true hazard is purposefully close to a age-cohort model. This simulation setting is aimed at quantifying how well the interaction term can be inferred by the ACI model. is very dissimilar from an additive effect between age and cohort. This is beyond the scope of this study.

In this section, we compare the predictive performance of the ACI model with the AIC, BIC, and EBIC criteria, to the AC model and to the MLE. We simulate the data a repeated number of L = 100 times. For each repetition, we compute the ACI, AC, and MLE estimates. The MLE estimate is explicit: λj,k = O j ,k /R j ,k . The AC estimate is not explicit but is easy to compute -we provide a Newton-Raphson-based estimation in the package hazreg. For each estimate, we then compute the mean squared error

MSE = 1 L L l =1 λl -λ * 2 ,
where λl is the estimated hazard for the l -th sample. This procedure is completed for different sample sizes: 400, 1000, 4000, and 10000. A sample size of 10000 is considered small for applications in epidemiology. With J = K = 20, this sample size is 25 times larger than the number of parameters.

The results are gathered in Table 3.1. The ACI model performs always (almost) as well as the AC model, which performs way better than the MLE. This is coherent with the fact that in three simulation designs, the true hazard is close (in L 2 norm) to the age-cohort model. For larger sample sizes, the ACI model outperforms the AC model (results not shown here). No criteria stands out as better for the task of model selection. These remarks hold likewise for all three simulation designs.

Consequently, the ACI model is shown to perform as well as the AC model, even when the true hazard is close to that of the AC model. In addition, the ACI infers the interaction effect. In the next section, we will qualify, how well the interaction term is estimated.

Perspective plots

In this section, we represent the hazard rate estimated by the ACI model. We also represent the interaction effect. We simulate datasets of sample size 40000 a repeated number of 100 times. We represent the pointwise medians of the estimates over the 100 replications. The estimates are represented in Figures 3.2,3.3,and 3.4 for the three simulation designs, respectively. Figure 3.2 represents the results for Simulation Design 1. We compare the hazard estimated with the ACI model (Figure 3.2b) and with the MLE (Figure 3.2a). The MLE infers accurately the two "bumps" forming the interaction effect, as well as the general trend forming the age and cohort effects. But the regularizing effect of the ACI model improves the quality of the estimation, especially for large values of the age, where there are fewer data points and thus a lower signal-to-noise ratio. We also represent the decomposition of the estimated hazard (Figure 3.2b) in the age and cohort effects and the interaction effect: Figure 3.2c represents exp(µ + α j + β k ) j ,k and Figure 3.2d represents exp(δ j ,k ) j ,k . The interaction effect is accurately estimated, which leads to think the ACI model is useful to detect the presence of interactions between the age and cohort effects. (We stress that the functions represented here are medians over 100 repetitions. The estimation of interaction effect with such precision as that of Figure 3.2 requires a greater sample size than 40000). The estimated age and cohort effects presented in Figure 3.2c are very close to their true values. Figures 3.2e and 3.2f represent the age (α 2 , • • • , α J ) and cohort (β 2 , • • • , β K ) effects respectively The estimates of each of the 100 replications are in light grey, their median is in red and the true value is in blue. The age effect is well estimated: the estimate is close to the true value and the replicated estimate highlight a symmetrical variance thereof. The cohort effect displays the same behavior, but with a decreased quality of estimation: the empirical variance is large and the median estimate presents two bumps, whereas the true cohort effect has none. This is explained by the small values taken by β * compared to that of µ * , α * , and δ * . We used the AIC selection criterion in this simulation. The same simulations were done with the BIC and EBIC criterion (results not shown here) and yield similar, albeit visually slightly worst, results. This remark holds for . sented in Figures 3.3b and 3.3c. The AC model fails to estimate the shape of the hazard. In particular, the interaction term cannot -by definition -be estimated by the AC model, and the estimated age and cohort effects are skewed to fit the interaction term. Figure 3.3 represents the interaction term. The estimated interaction is not visually close to a piecewise constant function. Its circular support is however well estimated.

Figure 3.4 represents the same results as Figure 3.3 for Simulation Design 3. The remarks are the same: the ACI model infers both the age-cohort trend and the interaction term satisfyingly. The added difficulty of this simulation design is the presence of two levels instead of one. Figure 3.4d shows that in this case, the interaction effect is still well estimated.

Conclusion

We have introduced a new model to estimate the interaction between effects in the age-periodcohort setting. This model was developed to generalize the age-cohort model and is called agecohort-interaction model. (The age-period and period-cohort models can be generalized in the same fashion.)

This model uses penalized maximum likelihood estimation. Indeed, a fused adaptive ridge method is used to regularize the interaction term. This regularization forces the interaction term to be piecewise constant. We propose different model selection criteria for the choice of the penalty constant: the AIC, the BIC, and the EBIC. We compare the criteria using simulations. The AIC is underpenalizing and consequently, the estimated interaction effect is piecewise constant with many different values (simulation results not shown here). The BIC and EBIC tend to select similar penalties and consequently to provide the same estimate, the EBIC being slightly more penalizing than the BIC (results not shown here).

The estimation performance of our model is evaluated on simulations and compared with that of the age-cohort model and the MLE. With sample size 400, 1000, 4000, and 10000, our method performs as well as the age-cohort model and outperform the MLE. For higher sample sizes, our model slightly outperforms the age-cohort model (results not shown here). Finally, perspective plots of the interaction term illustrate that the age-cohort-interaction model performs well at estimating the presence and support of interaction between age and cohort effects. The shape of the interaction is inferred correctly in the case where it is piecewise constant.

Perspectives include using a bootstrapping of the data to provide empirical confidence intervals of the estimate. Averaging the collection of estimates provides an estimate of the interaction term that is not piecewise constant. I have performed exploratory analysis in this direction, with promising results. Moreover, I have started to apply this method to the study of breast cancer incidence among Norwegian women. I plan to apply the ACI model to the data provided by the NOWAC cohort [START_REF] Lund | Cohort Profile: The Norwegian Women and Cancer Study-NOWAC-Kvinner og kreft[END_REF]. Finally, the ACI model could be used to test for the presence of a (non-linear) interaction effect. Indeed, if we had an asymptotic distribution for the estimated interaction effect, one could build a test with null hypothesis "∀ j , k, δ j ,k = 0". This would require to establish asymptotic properties of the adaptive ridge. We refer to the conclusion of this manuscript for perspectives on this matter.

∂ ∂α

j (θ) = K k=1 exp µ + α j + β k + δ j k R k j -O k j , ∂ ∂β k (θ) = J j =1 exp µ + α j + β k + δ j k R k j -O k j ,
and

∂ ∂δ j k (θ) = exp µ + α j + β k + δ j k R k j -O k j ,
with 1 ≤ j ≤ J and 1 ≤ k ≤ K . Its second order derivative is given by

∂ 2 ∂µ 2 (θ) = J j =1 K k=1 exp µ + α j + β k + δ j k R k j , ∂ 2 ∂α j ∂µ (θ) = K k=1 exp µ + α j + β k + δ j k R k j , ∂ 2 ∂β k ∂µ (θ) = J j =1 exp µ + α j + β k + δ j k R k j , ∂ 2 ∂δ j k µ (θ) = exp µ + α j + β k + δ j k R k j , ∂ 2 ∂α j ∂α j (θ) = 1 j = j K k=1 exp µ + α j + β k + δ j k R k j , ∂ 2 ∂β k ∂α j (θ) = exp µ + α j + β k + δ j k R k j , ∂ 2 ∂δ j k ∂α j (θ) = 1 j = j exp µ + α j + β k + δ j k R k j , ∂ 2 ∂δ j k ∂β k (θ) = 1 k =k exp µ + α j + β k + δ j k R k j , ∂ 2 ∂β k ∂β k (θ) = 1 k =k J j =1 exp µ + α j + β k + δ j k R k j , ∂ 2 ∂δ j k ∂δ j k (θ) = 1 j = j ,k =k exp µ + α j + β k + δ j k R k j , with 1 ≤ j , j ≤ J and 1 ≤ k , k ≤ K .
Likelihood with piecewise constant regions. Let Z l define the l -th region, for 1 ≤ l ≤ L. Let c be the constrained negative log-likelihood function, depending on the parameter vector θ = µ, α, β, δ (1) , . . . , δ (L) of size 1

+ J -1 + K -1 + L. Then λ j ,k = exp µ + α j + β k + L l =1 1 ( j ,k)∈Z l δ (l ) and c (θ) = J j =1 K k=1 λ j ,k R j ,k -log λ j ,k O j ,k . (3.10) Then ∂ c ∂δ (l ) (θ) = J j =1 K k=1 1 (j,k)∈Z l λ j ,k R j ,k -O j ,k ∂ 2 c ∂µ∂δ (l ) (θ) = J j =1 K k=1 1 (j,k)∈Z l λ j ,k R j ,k ∂ 2 c ∂α j ∂δ (l ) (θ) = K k=1 1 (j ,k)∈Z l λ j ,k R j ,k ∂ 2 c ∂β k ∂δ (l ) (θ) = K k=1 1 (j ,k)∈Z l λ j ,k R j ,k ∂ c ∂δ (l ) ∂δ (l ) (θ) = 1 l =l J j =1 K k=1 1 (j,k)∈Z l λ j ,k R j ,k
Penalized Likelihood. The penalized likelihood's derivatives with respect to µ, α, and β are unchanged. For 2 ≤ j ≤ J and 2 ≤ k ≤ K we have: (3.11) with by convention δ J +1,k = δ j ,K +1 = 0 and v J ,k = w j ,K = 0. For 2 ≤ j ≤ J , 2 ≤ k ≤ K , 2 ≤ j ≤ J and 2 ≤ k ≤ K we have:

∂ κ ∂δ j k (θ) = ∂ ∂δ j k (θ)+κ -v j ,k δ j +1,k -δ j ,k + v j -1,k δ j ,k -δ j -1,k +κ -w j ,k δ j ,k +1 -δ j ,k + w j ,k -1 δ j ,k -δ j ,k -1 ,
∂ 2 κ ∂δ j ,k ∂δ j ,k (θ) = ∂ 2 ∂δ j ,k ∂δ j ,k (θ) + κ 1 j = j ,k =k v j ,k + v j -1,k + w j ,k + w j ,k -1 -v j ,k 1 j = j +1,k =k -v j -1,k 1 j = j -1,k =k -w j ,k 1 j = j ,k =k +1 -w j ,k -1 1 j = j ,k =k -1 (3.12)

Appendix B: Fast Inversion of the Hessian matrix

The Newton-Raphson iteration requires the computation of I κ n -1 U κ n . This linear problem is the computational bottleneck of the algorithm. The Hessian matrix has the form

I κ n = A B B T D where A = ∂ 2 κ n ∂(µ, α, β)∂(µ, α, β) T , B = ∂ 2 κ n ∂(µ, α, β)∂δ T , and D = ∂ 2 κ n ∂δ∂δ T .
The Hessian matrix is inverted using the Schur block matrix inversion formula [START_REF] Zhang | The Schur Complement and Its Applications[END_REF]: .13) where

I κ n -1 = S -1 -S -1 B D -1 -D -1 C S -1 D -1 + D -1 C S -1 B D -1 , where S = A -B D -1 C is the Schur complement of the matrix D. Consequently (I κ n ) -1 U κ n = S -1 U 1 -S -1 B D -1 U 2 -D -1 B T S -1 U 1 + D -1 U 2 + B T S -1 B D -1 U 2 , ( 3
U 1 = ∂ κ n /∂ µ, α, β and U 2 = ∂ κ n /∂δ.
Since D is a symmetric band diagonal matrix, D -1 B T is computed using Cholesky factorization and the computation complexity of I κ n -1 U κ n using Equation (3.13) reduces from O (J K ) 3 to O max (J , K ) 2 J K .

Chapter 4

Spline regression with automatic knot selection

The problem of finding a function f that links a (potentially multivariate) explanatory variable x to a univariate response variable y is central in statistics. In all generality the link between x and y is not affine and the linear model is often deemed too restrictive. Often, the link between x and y is supposed complex and is unknown, and we do not want to make restrictive assumptions on f . This is the framework of non-parametric regression, where the unknown function f is assumed to belong to a family of functions that cannot be configured by R p (or a subset thereof). Instead, we infer the function from the data without additional assumptions on f . We count four important tools for non-parametric regression: Gaussian process regression, kernel regression, spline regression, and regression trees. Among those, spline regression has been praised for being simple and efficient in many practical cases. This domain has grown in popularity since the development of computational statistics in the years 1980's and has gained mainstream appeal in the years 2000's. Strictly speaking, spline regression does not fall inside the category of nonparametric regression, because splines are families of parametric functions. They are still termed non-parametric because (i) they are over-parametrized and each parameter has no interpretability inside the model and (ii) splines with a sufficiently large number of parameters can approximate any function inside some infinite-dimensional function space with arbitrary precision (the details are beyond the scope of this work, see [START_REF] Curry | On Pòlya Frequency Functions IV: The Fundamental Spline Functions and Their Limits[END_REF]. Some refer to spline regression as being semiparametric, which is also technically incorrect.

Spline regression is simple in practice: we set the number of knots, their position, and the order of the spline (knots are positions on the x-axis where the spline displays a discontinuity in its derivative and a spline's order accounts for its regularity). Conditionally on this choice, the regression is a parametric method and is easy to compute. However, with sufficiently large order and number of knots, the spline approximates any function inside a very large (i.e. infinite-dimensional) function set. This is the simplifying trick which makes spline regression both simple and powerful.

In practice, performing regression with one choice of knots (and order) is not enough: we need to select the best number and position of knots and order a posteriori. This is a complicated problem. The order of the spline is always chosen by the statistician. With sufficiently enough knots, a spline of sufficiently high order can approximate any spine of lower order. Thus the order is chosen as the lowest value such that the fit is deemed satisfactory (rarely more than 5 in practice). The knots are however harder to choose. The regression will change drastically if there are too many knots or if they are not well placed. This is the hidden cost of spline regression. Two main solutions have been brought forward. The first avoids the problem by voluntarily over-parametrizing the splines and taking overfitting into account with penalization. The second tackles the problem of optimizing with respect to the knots, either through a Bayesian framework or through a well-chosen optimization scheme.

In this chapter, we introduce a novel approach to automatic spline selection. It consists in jointly selecting the number and position of knots as well as the fitted spline. The method relies on the adaptive ridge and we iterate between penalized spline regressions to find the best unpenalized regression spline. We iteratively remove the least relevant knots from a large initial collection of knots. This approach takes a simplifying approach to the problem, but the computational burden is vastly reduced and simulations show that the resulting fit has good prediction performance. This chapter is organized as follows. In the first part, we introduce splines and spline regression through its different methods, with emphasis on interpretation and computational cost. We elaborate on the different unpenalized and penalized spline regression methods. The paper "Spline Regression with Automatic Knot Selection" is given in the second part.

We refer to [START_REF] Green | Nonparametric Regression and Generalized Linear Models[END_REF], Härdle (1990), and Wahba (1990) for reviews on spline regression. sense of interpolation. It turns out that polynomials of degree smaller of equal than a constant k ≥ 0 are a poor choice of function family for interpolation. For example, polynomial interpolation is globally sensitive to a local perturbation of the data, while spline interpolation will only change locally if the data is perturbed locally. See [START_REF] De Boor | A Practical Guide to Splines[END_REF], Chapter 2) for a detailed discussion of the drawback of polynomial bases for function approximation.

In the following, we will use splines as approximating functions f for the regression (4.1). We stress that since few knots are used (q ≤ 40 typically), and with an order k typically smaller than 5, the spline linear space has a low dimension. This makes the regression problem a parametric estimation problem, in a parameter space of small dimension. We first have to define a well-suited base of the spline linear space.

The truncated power basis

The truncated power basis is a basis for the spline linear space. For fixed q, t and k, the j -th truncated power basis spline if defined as

T j ,k,t (x) = x j if 0 ≤ j ≤ q -1 (x -t j -k ) q-1 + if k + 1 ≤ j ≤ k + q. (4.2)
We recall that the notation (x) + denotes the positive part max(x, 0); hence the naming of truncated power basis. The indices t , k, and t will be dropped from the notation T j ,k,t when they can be inferred from the context. We will show that these functions form a basis of the spline space. It is easy to see that it has k + q elements who are linearly independent. It remains to show that each element of the basis is a spline. For j ∈ {1, • • • , k -1}, T j ,k,t is a spline of order k because is it a polynomial of degree k -1. For j ∈ {k, • • • , k + q}, T j ,k,t is polynomial on both [a, t j ] and [t j , b] and it is continuously differentiable on [a, b] up to the order k -1. This completes the proof.

Numerical stability issues. An example of truncated power basis is given in Figure 4.1, for t = (t 1 , t 2 , t 3 ), = (0.5, 0.1, 0.15) and k = 3. The truncated power basis is then composed of the k + q = 6 elements x 0 = 1, x 1 = x, x 2 , (xt 1 ) 2 + , (xt 2 ) 2 + , and (xt 3 ) 2 + . The figure also illustrates the major problem with the truncated power basis: it yields poorly conditioned problems. We develop on this issue here. When the number of knots increases, two consecutive knots t j and t j +1 can be so close that the truncated basis functions (xt j ) q + and (x - .2: Illustration of the numerical instability of the decomposition in truncated power basis, with a spline of order 2 supported by 3 knots (tick line) and its decomposition with 2-digit machine precision coefficients (dotted line).

t j +1 )
q + become very highly correlated. As is shown after, the fitting of splines to a data set (x i , y i ) is done by linear regression. The design matrix of this linear regression is the n ×k +q matrix composed of the basis functions T j ,k,t (x i ) evaluated at the data points. When the columns of this matrix are too close to being correlated, numerical precision issues occur. This is due to the subtraction of numbers that are potentially large and have almost opposite values. In the truncated power basis, the base functions span over a large number of knots, so the error made on one coefficient has an effect on whole fitted function (more precisely, on the right side of the corresponding knot). This issue makes the representation of spline in the truncated power basis prone to numerical instabilities, and thus unfit to computation. As an illustrating example, consider the example adapted from (De Boor, 1978, page 85) of a spline of order q = 2 defined on [a, b] = [0, 1] and supported over the knots (t 1 , t 2 , t 3 , t 4 ) = (0.1, 0.5h/2, 0.5 + h/2, 0.8), with values (0, 1.1, -1.1, 0) at the knots, and where h is a small constant. This spline allows the following decomposition in the truncated power basis:

f (x) = α 1 (x -t 1 ) + + α 2 (x -t 2 ) + + α 3 (x -t 3 ) + + α 4 (x -t 4 ) + , where              α 1 = 1.1 t 2 -t 1 α 2 = -2.2 h + 1.1 t 2 -t 1 α 3 = 2.2 h + 1.1 t 4 -t 3 α 4 = -1.1 t 4 -t 3 .
When h gets small, we have α 2 -α 3 , and both get large in absolute value. With an illustrative machine precision of 2 significant numbers and h = 0.05, we have α 1 2.9, α 2 -47, α 3 48, and α 4 -4. Figure 4.2 illustrates the loss of quality suffered by this decomposition. This makes the truncated power basis unfit in practice. In the next section, we introduce the spline basis which solves all the pitfalls that the truncated power basis has.

B-spline basis

The B-spline basis is a reparametrization of the truncated power basis. It is defined by recurrence, for fixed order k ≤ 1 and fixed knots t , as where by convention t 0 = a and t q+1 = b. We can easily verify by recurrence that this formula defines indeed splines. There are k + q B-splines of order k, which form a base. The suffixes of the splines are dropped when they can be inferred from context. Consider first the B-spline basis of order 1: they are the characteristic functions of the inter-knot intervals [t j , t j +1 ). Compare now these functions to the truncated power basis of order 1: T j ,1,t = I t j ≤ x . Here, the B-splines are equal to the consecutive difference of two truncated power splines. This allows them to have a small support and to have have the same order of magnitude. Likewise, the B-splines of order 2 are triangular functions such that B j ,2 (t j -1 ) = 1, B j ,2 (x) = 0 for x ∈ [t j -2 , t j ). We can show that the B-splines have the following properties for any k:

B j ,1,t = I t j ≤ x < t j +1 B j ,k,t (x) = x -t j t j +k-1 -t j B j ,k-1,t (x) + t j +k -x t j +k -t j +1 B j +1,k,t (x), for k > 1, ( 4
• B j ,k (x) has exactly the support [t j -k , t j ) ∩ [a, b].
• At any point x ∈ [a, b], there are exactly k B-splines with non-zero values.

• The B-spline base is normalized, i.e. at any point x ∈ [a, b], k+q j =1 B j ,k (x) = 1. This makes the representation of splines in this base better fitted to computations. More precisely • The decomposition coefficients of a spline are of the same order of magnitude as the spline.

• These coefficients represent the local value of the spline, and have no global effect on the spline.

This makes B-splines optimally fitted to spline representation. For instance, in the example given in the last section (see Figure 4.2), B-spline decomposition recovers perfectly the initial spline.

For the sake of completeness, we give an illustration of B-spline bases of different orders and with equally spaces knots in Figure 4.3.

Computation of B-splines.

The recurrence formulas defining B-splines are directly used to compute B-splines [START_REF] De Boor | A Practical Guide to Splines[END_REF] gave an algorithm to compute B-splines of any order based on Bspline of lower order (recall that splines of order 1 are piecewise constant). This method is efficient for computing B-splines and most B-spline implementations use FORTRAN routines implementing de Boors' algorithm (de Boor, 1977).

Regression using splines

Spline regression as linear regression. We want to solve the regression problem (4.1) using a spline to estimate f . Given an order and the sequence of knots, we define f (x) = k+q j =1 α j B j (x), (4.4) where α = (α j ) is the vector parameter which gives the decomposition of f in B-splines. We chose the L 2 norm as our loss function, i.e. we select the estimated function which minimizes the sum of square

n i =1 (y i -f (x i )) 2 .
Notice now that from this equation and (4.4), the criterion to minimize writes

y -X α 2 ,
where the n ×(k + q) design matrix X has (i , j )-th entry B j (x i ). This is the framework of linear regression, whose estimate is explicit: α = (X T X ) -1 X T y.

From the aforementioned properties of B-splines, the design matrix has the following form. For j ∈ 0, • • • , q let n j be the number of points x i in the interval [t j , t j +1 ) (we have j n j = n). Then X has non-zero values on blocks of size n j × k ordered one above the other, each block being shifted by one entry compared to the previous block. An easy calculation shows that X T X is banded, of bandwidth k: its ( j , j )-th entry is zero if

| j -j | > k.
Numerical efficiency. As noted by (Hastie et al., 2001, Chap. 5, Appendix), the numerical cost is vastly reduced in this setting. Its inversion is computed efficiently in the following steps: (i) compute the Cholesky decomposition X T X + λΩ = LL T , with complexity O(nk 2 ), (ii) solve L x = X T y for x using forward substitution, in complexity O(n), and (iii) solve L T x = x for x using back substitution, in the same complexity. Overall, smoothing spline regression is computed in O(nk 2 ) time complexity. In practice, k is small (often set to 4) and so fitting B-splines is fast.

Remark. In the following, we will discuss about the properties of the underlying function f . Indeed, the order of the spline to choose and the number and location of knots all depend on the shape of f , as does the quality of the regression. In this chapter, we say that a function is regular, or smooth, if is has few non-zero high-order derivatives. This is a local notion, which means that a function can be smoother in one place of the segment [a, b] and less smooth in another. Since no assumption is made on f , we don't give a proper mathematical definition of regularity; it is used to discuss the difficulties of spline regression.

Choice of the knots. Spline regression is fast and easy to implement and it provides a flexible enough setting to approximate a large spectrum of functions. But the quality of the fit depends on the choice of the knots. With splines of order 1 and 2, the knots ought to be placed at the cut-points in the data, which can be selected by visual inspecting (although this method is not data-driven and highly arbitrary). For splines of higher order, the selection criterion is more difficult to choose.

It turns out that equally-spaced knots is not always a good choice, because the function f may have more variability in one region than in another. One may choose equally spaced knot, which seems to be a "neutral" setting, and which works well when the true function has enough regularity. However, a spline with too many knots will (i) be prone to "wiggliness" and (ii) be overfitted, because as the number of knots gets close to the number of data points, the fit becomes an interpolation. qq q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -150 -100

-50 0 50 0 20 40 60

(a) Helmet data, equidistant knots qq q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -150 -100

-50 0 50 0 20 40 60

(b) Helmet data, quantile knots q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 1.0 1.5 2.0 600 700 800 900 1000 1100

(c) Titanium data, equidistant knots q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 A spline with more knots has more flexible fitting to the data, i.e. it can fit to functions with high variability (in the sense that it has higher order derivatives with important values). Since B-splines are locally adaptive, a region with more knots will accommodate high variability of the function in this region. Thus, we would want to choose the knots in an adaptive, data-driven way that puts more knots when it is necessary to fit the high variability, and less when the function is smoother. This is a difficult task.

A possible solution is to place the knots at the quantiles of (x i ), allowing for higher adaptability where there are many data points. But this is efficient only if highly-varying regions have more data points, which is not always verified. We illustrate this problem on two different datasets (Figure 4.4). The "helmet" data (see de Souza and Heckman, 2013) set is helmet crash test data, representing acceleration against time (n = 132). It is repeatedly used to illustrate spline regression methods and compare them (Silverman, 1985;Eilers and Marx, 1996). The "titanium" data represents a heat property of titanium as a function of its temperature (n = 49) (used for instance by [START_REF] Dierckx | Curve and Surface Fitting with Splines[END_REF][START_REF] De | Least Squares Cubic Spline Approximation, II -Variable Knots[END_REF][START_REF] Jupp | Approximation to Data by Splines with Free Knots[END_REF]. Since this data contains close to no noise, we add a Gaussian noise of standard deviation 0.05 to the y values in order to make the inference reasonably challenging. Compare now Figures 4.4b and 4.4b, which display spline regressions of order 4 with 9 knots placed at the quantiles, for two different data sets. In the first one, the x i s are not spread uniformly, and the sharp shift at x ∼ 17 is where the x i s are the most dense. The knots set at the quantiles allow for many knots in this region, and the fit is of good quality comparatively to equidistant knots (Figure 4.4a). On the other one, the x i s are uniformly spread, and there is a high variability in one region. This lead a poor fit (at least as poor as with equidistant knots, see Figure 4.4c).

The problem of finding the "best" position of the knots is a difficult one. In the next section, we present a few solutions present in the literature.

Penalized approaches

In this part, we will see some methods that circumvent the difficulty of choosing the knot's position. Instead, they set a (too) large number of knots and deal with overfitting through penalized estimation. They all differ in their choice of knot number and placement and type of penalization, but essentially follow this simple principle.

Natural cubic splines. We first introduce a family of spline called natural cubic spline.

Spline of order 2 are continuous, splines of order 3 are continuous and have continuous first order derivatives, splines of order 4 are continuous and have continuous first and second order derivatives. For order ≤ 3, the splines' shifts at the knots are visible. For order 4 and larger, it becomes almost impossible to tell where the knots are: the splines are seamless. This is why in practice splines of order 4 are chosen when f is assumed "smooth" without any other assumption. They are called cubic splines -they are made of piecewise cubic polynomials.

Around the boundary of [a, b], there are fewer data points then in the center of [a, b]. Consequently, the regression has more variance close to the margins a and b. To remedy this issue, one solution is to add constraints at the boundary. For cubic splines, we can add the constraints that the spline be linear on (-∞, t 1 ] and [t q , ∞), i.e. that the second and third order derivatives are zero left of t 1 and right of t q . This adds 4 constraints on the spline, which can now be parametrized by k + q -4 = q freely varying parameters. Definition 4. A natural cubic spline is a spline of order 4 linear outside of its boundary knots.

Smoothing splines

Consider the problem of finding the function f which solves arg min

f n i =1 (y i -f (x i )) 2 + λ R f (t ) 2 d t (4.5)
where f is any function such that this quantity is defined and λ > 0 is a smoothing parameter. Consequently, f belongs to the Sobolev space of functions whose derivatives of order 0, 1, and 2 are squareintegrable. Equation 4.5 is minimization problem in infinite dimension. Fortunately, the solution to (4.5) is proven to be a natural cubic spline with n knots placed at the x i s. These splines are called smoothing splines. Consider a basis (N j (x)) 1≤ j ≤n for the space of cubic splines defined on the knots (x i ) and write f (x) = α T N (x) (we can define such a basis with properties similar to B-splines. We do not give the details here for the sake of the presentation and refer to Hastie et al., 2001, Section 5.2.2). Then problem (4.5) rewrites

min α y -X α 2 + λα T Ωα (4.6)
where X is the design matrix of the values of N j (x i ) and

Ω i , j = R N i (t )N j (t )d t (4.7)
is the n × n symmetric matrix enforcing the smoothing of the spline. Consequently, the smoothing spline is not over-parametrized, because its estimated vector is the solution of a ridge regularized regression: α = (X T X + λΩ) -1 X T y.

Since Ω is positive definite, it acts as a regularization of the estimate, in a way similar to ridge regression (if Ω were diagonal, the previous equation would be a weighted ridge regression). Note that Ω acts as a weighting matrix and it depends on x and λ, but not on y. Hence the term smoothing spline:

the projection matrix X (X T X +λΩ) -1 X T y that maps y to ŷ is obtained by a smoothing independent from y. It comes that X T X + λΩ is symmetric definite positive and 4-banded. Since smoothing splines are of order 4, we know that X T X is also 4-banded. Consequently smoothing splines have the same computational cost as B-splines.

O'Sullivan Penalized Splines

O'Sullivan (1986) introduced a penalized regression spline that it closely related to smoothing splines. Indeed, the method (called O'Sullivan splines for clarity) uses cubic B-splines (although any order could be used without modifications to the method) with k +4 knots located anywhere on [a, b]. The estimate is the cubic spline

f (x) = k+4 j =1 α j B j ,4 (x)
Consider the n×(k+4) design matrix X with (i , j )-th entry B j (x i ). O'Sullivan penalized splines use the same penalization as smoothing splines, but in the context of B-splines (instead of natural splines). Likewise (4.5), the estimate is defined as arg min

α y -X α + λ x min x min k+4 j =1 α j B j ,4 (t ) 2 d t , (4.8)
and likewise smoothing splines, the estimate reduces to â = (X T X + λΩ) -1 X y, (

where the penalty matrix Ω is defined by (4.7), where N j is replaced by B j ,4 . The rationale behind (4.8) is similar to smoothing splines. Using many knots creates wiggliness in the spline. For a spline of order 4, the penalization of its integrated squared second order derivative forces it to have a small second order part and to have small third order derivative differences around each knot. Consequently, smoothing splines are a particular case of O'Sullivan splines, with k = n knots placed at the data points x i . Whereas smoothing splines perform a trade-off between interpolation (for λ → 0) and linear regression (for λ → ∞), this method performs a trade-off between unpenalized B-splines (for λ → 0) and linear fit (λ → ∞). It enjoys the same computational efficiency as smoothing splines, since the penalty function is also 4-banded. It is in fact more efficient since we usually define way fewer knots then there are data points (n ≤ k). This makes O'Sullivan splines a tool of choice for penalized spline regression. In the next section, a further simplification of this method is presented. Eilers and Marx (1996) introduced a penalized spline regression method called P-splines (for penalized splines). Its simplicity made it the most widely used penalized spline regression model, having been used in many applications and extended to more general settings [START_REF] Eilers | Twenty Years of P-splines[END_REF][START_REF] Currie | Smoothing and Forecasting Mortality Rates[END_REF][START_REF] Wand | On Semiparametric Regression with O'Sullivan Penalised Splines[END_REF][START_REF] Jiang | Age-period-cohort models using smoothing splines: A generalized additive model approach: B[END_REF].

P-splines

Consider the same framework as the previous section: the regression spline is a B-spline of any order k. P-splines are defined with knots equidistant over [a, b].

Recall that ∆α j α j -α j -1 denotes the difference operator and ∆ p ∆ p-1 • ∆, p ≥ 2 denotes the p-th composition of the difference operator with itself. P-splines are defined as the solution of the penalized problem arg min (4.10) This penalty has an intuitive explanation, at least with the order 2. In this paragraph, consider q = 2, i.e. piecewise linear functions. Recall that f = j α j B j ,2 (x) and that over [t j -1 , t j ), B j 0 ,2 (x) has slope -1, B j 0 -1,2 (x) has slope 1 and all the other B-splines are equal to zero. Then at t j 0 , the slope shift of f (x) is αj 0 -αj 0 -1 . When (∆ 2 αj 0 ) 2 is small, αj 0 αj 0 -1 , and the slope shift of f (x) at t j 0 -1 vanishes. This enforces the spline to have small variations in slope at each knot, and throughout [a, b]. This result generalizes to other orders, and in particular cubic splines will have small differences in their third order derivatives.

α y -X α 2 + λ q+k j =k (∆ k-2 α j ) 2 .
The penalized problem (4.10) gives the explicit estimate α

= (X T X + λD T D) -1 X T y,
where D is the q ×(q +k) matrix of the operator ∆ k-2 . The entries of D are the signed binomial coefficients. Consequently, D T D is much simpler to compute than the penalty matrix Ω used in smoothing splines or O'Sullivan splines, which depend on the knots. The difference in computational cost and complexity increases for higher orders of the spline. We will see that the penalty used in P-splines is however quite close to that of O'Sullivan splines. For comparison we also assume k = 4 in this paragraph. From (De Boor, 1978, p. 115)'s formula of spline derivatives, we have the identity

(t j -t j -1 )B j ,k (x) = k B j ,k-1 (x) -B j +1,k-1 (x) .
Using this formula, the O'Sullivan spline penalty rewrites

x n x 1 j α j B j ,4 (x) 2 d x ∝ x n x 1 j ∆ 2 α j B j ,2 (x) 2 d x (4.11) = x n x 1 j (∆ 2 α j ) 2 B j ,2 (x) 2 d x + x n x 1 j ∆ 2 α j ∆ 2 α j -1 B j ,2 (x)B j -1,2 (x)d x (4.12) = c 1 j (∆ 2 α j ) 2 + c 2 j ∆ 2 α j ∆ 2 α j -1 , (4.13)
where the indices are defined implicitly by the order of the spline in each sum. The constants c 1 and c 2 are integrals of second order splines over [x j , x j +1 ); they are proportional to 2/3 and 1/6 respectively. These results follow from simple algebra and are not detailed further (see Eilers and Marx, 1996, Section 3).

From this result we see that where the P-Spline penalty for the j -th knot is proportional to (∆ 2 α j ) 2 , that of O'Sullivan is proportional to 2/3(∆ 2 α j ) 2 +1/6∆ 2 α j ∆ 2 α j -1 . These two penalties are quite similar for many values of α (and even more so as ∆ 3 α gets small). This gives an understanding of the close relation between the two methods. We don't investigate this topic further (see [START_REF] Wand | On Semiparametric Regression with O'Sullivan Penalised Splines[END_REF]. P-splines are praised for their simplicity. In comparison, O'Sullivan splines are somewhat numerically less simple, although [START_REF] Wand | On Semiparametric Regression with O'Sullivan Penalised Splines[END_REF] provides explicit formula for the penalty matrix in (4.9). O'Sullivan splines have been generalized to any order, so both methods stand equal in terms of generality. P-splines can be deemed too restrictive for their choice of setting equal knots. But in fact the penalization that it enforces will compensate for knots that are set too close (simulation work by Ruppert, 2002, show that above a certain threshold, the number of knots has little impact on the performance of P-splines).

Spline regression with automatic knot selection

In this section, I include the preprint "Spline regression with automatic knot selection".

Introduction

After witnessing great developments in the past decades (see Wahba, 1990;[START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Ruppert | Semiparametric Regression During 20032007[END_REF][START_REF] Wood | Generalized Additive Models: An Introduction with R, 2 edn[END_REF] spline regression has become a tool of choice for semiparametric regression: it is restrictive enough to benet from the simplicity of parametric estimation and general enough to accurately approximate a large range of smooth functions. For a while, the knots over which the splines are built was arbitrary chosen by the data analyst. Yet the number of knots has an important inuence in the resulting t:

with not enough knots the regression is undertted and with too many knots it is overtted.

Choosing the knots' position is also important since uniformly distributed knots can lead to overtting in an area where there are few points and undertting in an area where there are many points.

During the 1990's, penalization methods have been developed to overcome this diculty.

The idea is to set (too) many knots and to control overtting by penalizing over the spline parameters. In smoothing splines (see Hastie et al., 2001, Section 5), knots are set at each data point and the spline's wiggliness is controlled by penalizing over the integrated squared second order derivative {f (t)} 2 dt. The smoothing spline estimate has a closedform expression and is computationally ecient. O'Sullivan (1986) generalized smoothing splines to an arbitrary choice of knots, allowing to set fewer knots than the sample size.

Later, Eilers and Marx (1996); [START_REF] Marx | Direct Generalized Additive Modeling with Penalized Likelihood[END_REF] introduced P-splines, based on a penalty over the nite order dierences of the parameters. This penalization is closely related to that of O'Sullivan (see Eilers and Marx, 1996, Section 3) and is a generalization of Whittaker (1922)'s graduation method, which can be seen as a P-spline of order 0 with knots placed at data points. See [START_REF] Eilers | Twenty Years of P-splines[END_REF] and citations therein for a review of P-splines and see [START_REF] Wand | On Semiparametric Regression with O'Sullivan Penalised Splines[END_REF] and Eilers et al. (2015, Appendix A) for a comparison between O'Sullivan's splines and P-splines. Smoothing splines are implemented in the R packages gam [START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Hastie | Gam: Generalized Additive Models[END_REF] and mgcv [START_REF] Wood | Generalized Additive Models: An Introduction with R, 2 edn[END_REF] and Psplines are implemented in the R package pspline.

These regularized spline regression methods are simple and computationally fast. However, a spline with fewer knots is easier to interpret, which in many cases is a desired goal. Indeed a knot corresponds to a discontinuity in the spline's kth derivative (where k is the spline's order) and thus to a shift in the trend of the estimated function. As pointed out by [START_REF] Wand | A Comparison of Regression Spline Smoothing Procedures[END_REF], setting a very large number of knots and exploring the set of splines dened on any subset of the knots is not tractable in practice. [START_REF] Friedman | Multivariate Adaptive Regression Splines[END_REF] has developed a multivariate variable selection technique called MARS (Multivariate Adaptive Regression Splines). It uses a recursive partitioning of the domain and sequentially selects the knots with a forward/backwards step size procedure (see [START_REF] Friedman | Flexible Parsimonious Smoothing and Additive Modeling[END_REF] and (Hastie et al., 2001, Section 9.4) for details). [START_REF] Luo | Hybrid Adaptive Splines[END_REF] have later developed a closely related approach called Hybrid Adaptive Splines which uses a forward stepwise procedure and ts penalized splines instead of using a backward procedure. Other paths have been taken to solve this computationally intensive problem. Namely, Jamrozik et al.

(2010) have oered to estimate the best location of knots using a dierential evolution algorithm. However, their approach was limited to a number of knots varying between 4 and 7 and to splines of order 1.

A dierent approach is to jointly estimate the spline coecients and the position of the knots, as well as their number. The splines estimated by allowing the knots' position to vary are called free-knot splines. [START_REF] Jupp | Approximation to Data by Splines with Free Knots[END_REF] was among the rst to introduce a method for free knot spline tting. The squared residuals are regarded as a function of both the knot positions and the spline coecients and the procedure uses a reparametrization of the knots' location and the Levenberg-Marquardt optimization algorithm [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF] to avoid staying in local minima. [START_REF] Lindstrom | Penalized Estimation of Free-Knot Splines[END_REF] improved on this method by adding a penalty over the dispersion of the knots. Secondly, Bayesian methods have also been widely used to select the knots' number and position, which are viewed as priors. In this context, [START_REF] Denison | Automatic Bayesian curve tting[END_REF] developed a reversible jump Monte Carlo-based method (RJMCMC, see [START_REF] Green | Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination[END_REF]. For each knot set, unpenalized least square regression is used to t the spline and the RJMCMC algorithm is used to roam the space of possible knot number and positions.

The best knots are selected using a model averaging approach. Biller (2000) extended this approach to the case where the coecients are also estimated in a Bayesian framework; [START_REF] Dimatteo | Bayesian Curve-Fitting with Free-Knot Splines[END_REF] and [START_REF] Holmes | Bayesian Regression with Multivariate Linear Splines[END_REF] presented developments around the same lines, but limited to 3 rd order splines. [START_REF] Leitenstorfer | Knot selection by boosting techniques[END_REF] introduced a method of knot selection using families of radial functions (in practice Gaussian densities) instead of splines. It uses a boosting algorithm to iteratively select the best tting function among a predened number of candidate functions, to estimate its parameter, and to add it to the current t. Finally, penalized estimation has been used previously by [START_REF] Osborne | Knot Selection for Regression Splines via the LASSO[END_REF] to perform knot selection. This work uses power series basis splines and selection is performed using the Lasso (Tibshirani, 1996). This penalty is only used for knot selection and an unpenalized t is used to infer the coecients selected by the Lasso.

Stepwise selection has been used for a long time for knot selection (see [START_REF] Smith | Curve tting and modeling with splines using statistical variable selection techniques[END_REF][START_REF] Stone | Polynomial Splines and their Tensor Products in Extended Linear Modeling[END_REF]. A formulation of the procedure for one-dimensional splines is available in [START_REF] Wand | A Comparison of Regression Spline Smoothing Procedures[END_REF]. Not many of these methods for knot selection have a currently available implementation in R. The method developed by DiMatteo et al. ( 2001) was implemented in C with an R wrapper by [START_REF] Wallstrom | An Implementation of Bayesian Adaptive Regression Splines (BARS) in C with S and R Wrappers[END_REF] and is publicly available 1 . The package freeknotsplines implements [START_REF] Spiriti | Knot selection for least-squares and penalized splines[END_REF]'s method for knot selection.

In this article we introduce a new computationally ecient method called A-splines, for adaptive splines to automatically select the number and position of the knots from the data. Contrarily to most free-knot spline methods, we set a high number of evenly distributed initial knots and use a penalized likelihood approach to gradually remove the least relevant knots. It is based on a regularization method with an approximate L 0 norm penalty. Therefore, it benets from the simplicity and computational speed of penalizedbased methods and is orders of magnitude faster to run than comparable knot-selection methods.

Section 2 gives a short summary of B-splines and B-spline regression. Section 3 introduces our spline regression method. In Section 4, our method is extended to the generalized linear model. Section 5 deals with the choice of the bias-variance tradeo parameter. Section 6 compares the prediction performance of our model to comparable methods through a simulation study. Section 7 gives some details about the fast implementation of the tting algorithm. Finally, A-spline is illustrated on several real datasets in Section 8. Note that with 3 knots, there are 4 splines in the basis of degree 0 and 7 splines in the basis of degree 3.

B-spline Regression

B-spline Basis

In this section we recall the denition and some basic properties of splines and B-splines.

Throughout this work, let t 1 , . . . , t k be the ordered knots included in a real interval [a, b].

A spline of order p ≥ 1 is a piecewise polynomial function of degree p -1 such that its derivatives up to order p -2 are continuous at every knot t 1 , . . . , t k . In the following, we will refer to the degree q ≥ 0 of a spline and not to its order q + 1. The set of splines of degree q ≥ 0 over the knots t = (t 1 , . . . , t k ) is a vector space of dimension q + k + 1.

A possible choice of spline basis is the truncated power basis: {x 0 , . . . , x q , (xt 1 ) q + , . . . , (xt k ) q + }, where (u) + = max (u, 0). The rst q + 1 functions of the basis are polynomials and the other k functions are truncated polynomials of degree q. Decomposing a spline into the truncated power basis brings out powers of large numbers, which lead to rounding errors and numerical inaccuracies (De Boor, 1978, p. 85).

Schoenberg introduced [START_REF] Schoenberg | Contributions to the problem of approximation of equidistant data by analytic functions. Part B. On the problem of osculatory interpolation. A second class of analytic approximation formulae[END_REF][START_REF] Curry | On Pòlya Frequency Functions IV: The Fundamental Spline Functions and Their Limits[END_REF] a spline basis called B-splines for Basic-splines. This spline basis provides more stable computations of spline regression (see [START_REF] De Boor | On calculating with B-splines[END_REF]. A B-spline is a spline which is non-zero over [x k , x k+q+1 ] for some k. For i = 1, . . . , q + k + 1, the i-th B-spline of degree q is noted B i,q (x) and is dened by

B i,q (x) = x -t i t i+q -t i B i,q-1 (x) + t i+q+1 -x t i+q+1 -t i+1 B i+1,q-1 (x) if q > 0 and B i,0 (x) = 1 t i ≤x<t i+1 .
Important properties of a B-spline are: (i) the B-spline is nonzero over an interval spanning q + 2 knots; (ii) at a point, only q + 1 B-splines are non-zero;

(iii) B i,q (x) ∈ [0, 1]. An illustration of B-spline bases of degree 0 to 3 is given in Figure 1. In practice, B-splines can be computed using the function bSpline from the R package splines2 [START_REF] Wang | Splines2: Regression Spline Functions and Classes[END_REF].

B-Spline Regression

Let (x i , y i ) ∈ R × R be the univariate data and consider the non-parametric regression setting

y i = f (x i ) + ε i , 1 ≤ i ≤ n, (1)
where f is a smooth function and where ε i are i.i.d. Gaussian errors. The function f is estimated by a spline over an interval [a, b] containing all x i s. Fitting the data consists in minimizing the sum of squares

SS (a, t) = n i=1 y i - q+k+1 j=1 a j B j,q (x i ) 2 , (2)
where a = (a 1 , . . . , a q+k+1 ) is the B-spline coecients, and t are kept as parameters to highlight that the tting procedure depends on the choice of the knots. This is the framework of ordinary least squares regression with design matrix B = [B j,q (x i )] i,j and parameter a:

SS (a, t) = y -Ba 2 2 .

(

) 3
3

Automatic Selection of Knots

Model selection

When there are many knots, spline regression is prone to overtting. In the extreme case, when there as as many parameters as data points, the tted spline interpolates the data.

In this paper, we propose to estimate the spline which makes the best tradeo between model dimension (i.e. number of knots) and goodness of t. To this eect, we choose a high number of equally spaced initial knots and penalize over the number of knots. When a B-spline is dened over the knots t 1 , . . . , t k and is such that ∆ q+1 a j * = 0 for some j * , it can be reparametrized as a B-spline over the knots t 1 , . . . , t j * -1 , t j * +1 , . . . t k . Consequently, we penalize over the number of non-zero (q + 1)-order dierences

λ 2 k j=q+2 ∆ q+1 a j 0 , (4)
where . 0 is the L 0 norm, i.e. x 0 = 0 if x = 0 and x 0 = 1 otherwise, and where the parameter λ > 0 tunes the tradeo between goodness of t and regularity. This penalty allows to remove a knot t j * that is not relevant for the regression, to merge the adjacent intervals [t j * -1 , t j * ) and [t j * , t j * +1 ), and to continue the tting procedure with a spline dened over the remaining knots. When λ → 0, the tted function is a B-spline with all knots t 1 , . . . , t k and when λ → ∞, the tted function is a polynomial of degree q.

However, the penalty in Equation ( 4) is non dierentiable and the estimation is therefore computationally non-tractable. An approximation method for the L 0 norm is introduced in the next section to overcome this diculty.

Adaptive ridge

Following the work from [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF] and [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF], we approximate the L 0 norm by using an iterative procedure called Adaptive Ridge. The new objective function is the weighted penalized sum of squares

WPSS (a, λ) = y -Ba 2 2 + λ 2 q+k+1 j=q+2 w j ∆ q+1 a j 2 , (5)
where ∆a j = a ja j-1 is the rst order dierence operator, ∆ i a j = ∆ i-1 ∆a j , and w j are positive weights. The penalty is close to the L 0 norm penalty when the weights are iteratively computed from the previous values of the parameter a following the formula:

w j = ∆ q+1 a j 2 + ε 2 -1 ,
where ε > 0 is a small constant. Indeed the function x → x 2 / (x 2 + ε 2 ) approximates the function x → x 0 if ε is suciently small. In practice, one typically sets ε = 10 -5 [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF]. At convergence, (∆ q+1 a j ) 2 w j ∆ q+1 a j 0 gives a measure of how relevant the j-th knot is. One chooses a threshold of 10 -2 and selects the knots with a weighted dierences higher than 0.99, which we note t sel j . The number of selected knots will be noted k λ , such that the number of parameters of the selected spline is q + k λ + 1.

(
Since the selected knots are present in breakpoints of the curve, one then ts unpenalized Bsplines over the knots t sel , as explained in Section 2.2. Consequently, this method provides a regression model that is both regularizing and simple, in the sense that the model dimension is small.

We note that [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF] give a more general formula for the weights that allows to approximate any L p norm, for p > 0. In particular, the L 1 norm could be chosen, which induces both shrinkage and selection of the coecient. This penalty was already used in this context by Tibshirani (1996), which used in the Lasso for knot selection.

WPSS (a, λ) of Equation ( 5) easily rewrites

y -Ba 2 2 + λD T W Da, (6)
where W = diag (w) and D is the matrix representation of the dierence operator ∆ q+1 .

The minimization of WPSS is explicit:

â = B T B + λD T W D -1 B T y. (7)
The adaptive ridge procedure is detailed in Algorithm 1.

The penalty term is conveniently written with the circulating matrix D. However, for computational eciency, D is never computed and instead we implement a fast computation algorithm for the penalty term. More details about the implementation are given in Section 7. a new ← arg min a WPSS (a, λ)

5:

w j ← ∆ q+1 a new j 2 + ε 2 -1 6:
a ← a new 7:

end while 8:

Compute t sel using (∆ q+1 a) 2 w 9:

â ← arg min a SS a, t sel

where the dierence order p is a parameter to be chosen. Thus, the implementation of A-splines can be seen as a weighted P-splines tting. The philosophies of A-splines and P-splines are however very dierent. P-splines avoid choosing the best knots by penalizing over the dierences of the coecients. Instead, we directly choose the best knots for spline regression.

Generalized Linear Model

Spline regression has also been used to t values in the general linear model setting, like in Eilers and Marx (1996); [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. In this section, we extend A-spline regression to the generalized linear model. The goal is to estimate µ = E [y|x] = g -1 (Ba), where g is the canonical link function and the variance of y is a function V of µ: Var [y] = V (µ). Like the linear model, µ can be estimated using spline regression. The generalized linear model is tted using the Iteratively Reweighted Least Squares (IRLS) algorithm (McCullagh and Nelder, 1989, Section 2.5). With weighted penalization, the IRLS iteration writes:

â(k+1) = B T Ω (k) B + λD T W D -1 B T Ω (k) Bâ (k) + y -µ (k) (9)
where k is the step index and Ω (k) is the diagonal matrix with entries

ω (k) i,i = 1 V µ (k) i g µ (k) i 2 , with µ (k) i = g -1 B i â(k) .
In practice, the estimation procedure in Algorithm 1 remains the same, except that WPSS is minimized by the Newton-Raphson procedure given in Equation ( 9).

Choice of the Penalty Constant

In this section, one selects the penalty that performs the best trade-o between goodness of t and regularity. A rst criterion is the AIC, which was used by Eilers and Marx (1996) in a similar context:

AIC(λ) = SS (â λ ) + 2 (q + k λ + 1) .

(10)

A dierent criterion is the Bayesian Information Criterion (BIC) (see Schwarz, 1978):

BIC (λ) = SS (â λ ) + (q + k λ + 1) log n. The BIC is the Bayesian criterion obtained when one chooses a uniform prior on the model: π(M λ ) = 1. As explained by ak-Szatkowska and Bogdan (2011), a uniform prior on the model is equivalent to a binomial prior on the model dimension. Therefore, the BIC tends to give too much importance to models of dimensions around q+k+1 2 . Since the adaptive knot selection is performed with a large number of initial knots, this will result in underpenalized estimators. To this eect, [START_REF] Chen | Extended Bayesian information criteria for model selection with large model spaces[END_REF] have developed an extended Bayesian information criterion called EBIC 0 . The EBIC 0 criterion is dened by choosing:

π(M λ ) = q + k + 1 q + k λ + 1 -1 and EBIC 0 (λ) = SS (â λ ) + (q + k λ + 1) log n + 2 log q + k + 1 q + k λ + 1 . ( 12
)
The EBIC 0 assigns the same a priori probability to all models of same dimension.

Therefore the EBIC 0 will tend to choose sparse models even with a high number of initial knots. These criteria's selection performances are compared in the next section through a simulation study.

6 Simulation Study

Comparing the Selection Criteria for A-spline

A simulation study has been conducted to compare the performances of the three criteria.

Data are simulated as follows. The x i are taken uniformly over [0, 1] and y i are simulated using Equation ( 1), where f is a known function and ε i ∼ N (0, σ 2 i ). We use four dierent functions: the Bump function

f 1 (x) = 0.4 x + 2 exp -{16 (x -0.5)} 2 , the Logit function f 2 (x) = 1 1 + exp {-20 (x -0.5)} ,
the Sine function f 3 (x) = 0.5 sin (6πx) + 0.5, q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0.5 0.0 0.5 1.0 0.00 0.25 0.50 0.75 1.00

(a) Logit Function q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.5 1.0 0.00 0.25 0.50 0.75 1.00 (b) Sine Function q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.4 0.8

1.2 0.00 0.25 0.50 0.75 1.00 (c) Bump Function q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0

f 4 (x) = x (1 -x) sin 2π 1 + 2 -3/5
x + 2 -3/5 + 0.5.

These functions were used by [START_REF] Wand | A Comparison of Regression Spline Smoothing Procedures[END_REF] and Ruppert (2002) in similar contexts for benchmarking the eciency of spline regression. The functions f 1 to f 4 have been rescaled

in order to vary in [0, 1], so that all simulation cases have similar signal-to-noise ratios. We choose homoscedastic errors σ i = 0.15 for the functions Logit and Sine and heteroscedastic errors for the Bump and SpaHet functions:

σ i = 0.3x i + 0.2 √ x i 2
, so that the variance increases from 0 when x = 0 to 0.25 when x = 1. Data are simulated with sample sizes 50, 100, 200, and 400. Illustration of the functions and of the simulated data are given in Figure 2. For each example 500 datasets were simulated. A-splines are tted and we compare the Mean Squared Error (MSE) of the estimated function for the three criteria:

f -f 2 2 = 1 0 f (x) -f (x) 2 dx.
The median MSEs are displayed in Table 1 for each value of the sample size. For all functions and for all criteria, the MSE decreases with the sample size, as is expected. The comparison between the criteria brings the same conclusions for all four functions: the BIC and EBIC 0 always perform better than the AIC. Moreover, note that the EBIC 0 always outperforms the BIC for the sample size 100, and performs almost as well for the sample size 200. In conclusion, the BIC and EBIC 0 are to be preferred over the AIC and overall, the EBIC 0 seems a better choice than the BIC. In the remaining of this work, A-splines will be used with the EBIC 0 criterion.

Comparing A-spline with knot selection methods

In this section, the performance of A-spline is compared to that of other knot selection methods. As mentioned in the introduction, only two knot selection methods are currently available in R. [START_REF] Dimatteo | Bayesian Curve-Fitting with Free-Knot Splines[END_REF]'s method, implemented by [START_REF] Wallstrom | An Implementation of Bayesian Adaptive Regression Splines (BARS) in C with S and R Wrappers[END_REF], is essentially an RJMCMC over the space of knot number and position. [START_REF] Spiriti | Knot selection for least-squares and penalized splines[END_REF]'s method, implemented in the package freeknotsplines initiates the knots uniformly and uses an optimization algorithm to search for the best location of each knot between its two neighboring knots. The method performs four passes through the list of ordered knots:

two passes in increasing order and two in decreasing order. Two choices are provided for the optimization algorithm: a blind search or a genetic algorithm (see [START_REF] Haupt | Practical Genetic Algorithms[END_REF]). In the genetic algorithm, the crossover step consists in keeping the leftmost knots of one parent and the rightmost knots of the other parent. The mutation step consists in choosing one knot at random and sampling it uniformly between its neighboring knots.

The selection function is set to either the residual sum of squares or the general cross validation (GCV). With both algorithms, the number of knots is chosen using an adjusted GCV criterion. Of the two optimization algorithms, we will use the genetic algorithm in our simulation setting.

We use the same simulation setting as in the previous section. Figure 3 represents the MSE distribution over 500 replications for every sample size and for every function.

We see that the three methods display close to equal performances. Across all sample sizes and all functions, the BARS seems to perform better than the genetic optimization q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Bump SpaHet Logit Sine the algorithm does not depend on the sample size.

Adaptive Spline Regression with Several Penalties

The penalty constant λ tunes the tradeo between goodness of t and regularity. To choose the optimal λ, regression is performed for a sequence of penalties λ = (λ ) , 1 ≤ ≤ L and a criterion is used to determine which regression model to select. Computing the procedure for a series of values of λ signicantly increases the computing time. Note that a small variation of λ yields a small variation of âλ = arg min a WPSS (a, λ). Consequently, âλ is a good initial point for the minimization of WPSS(a, λ +1 ). Making use of this hot start signicantly speeds up the minimization of WPSS (a, λ +1 ) and thus decreases the computation time of the adaptive ridge procedure. This implementation of the adaptive ridge is introduced in Rippe et al. ( 2012) and [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF] and a similar idea is used in the implementation of the LASSO in the package glmnet [START_REF] Friedman | Regularization Paths for Generalized Linear Models via Coordinate Descent[END_REF].

Fast Computation of the Weighted Penalty

The matrix inversion in Equations ( 7) and ( 9) is the computational bottleneck of the adaptive ridge procedure. The matrix D T W D is symmetric and q-banded, and as noticed by [START_REF] Wand | On Semiparametric Regression with O'Sullivan Penalised Splines[END_REF], so is B T B. Consequently, the inversion is done using Cholesky decomposition and back-substitution, as implemented in the package bandsolve 3 .

This reduces the temporal complexity from O (k + q + 1) 3 to O ((k + q + 1) (q + 2)). For example, if k = 50 and q = 3, the computation time will be reduced by a factor 500. It is important to note that the matrices W and D are not stored in memory: only the vector w and the rst row of D are used. This leads to improvements in spatial complexity, the details of which are not given here.

Fast Computation of the Weighted Design Matrix

In the setting of generalized linear regression, the matrix product B T ΩB in Equation ( 9) is computed at each iteration of the Newton-Raphson procedure. Since the design matrix has 3 github.com/monneret/bandsolve n rows, this operation makes the generalized linear regression computationally expensive for large datasets. Fortunately B is sparse: it has q + 1 non-zero elements in each row. Due to this structure, the product B T ΩB only has (q + k + 1) (q + 1) non-zero entries. Each entry takes O n k operations to compute on average. Thus the matrix product can be computed with a O ((q + k + 1) (q + 1) n/k) temporal complexity, compared to the O (q + k + 1) 2 n complexity of the naive implementation. For instance, even with q = 3 and k = 50, this implementation is faster by a factor ∼ 700.

Real Data Applications

Our method is illustrated with several real data applications. All the data sets used in this section and the code used to produce the gures are available in Supplementary Materials.

We rst present a dataset of simulated motorcycle accidents used to crash-test helmets.

The data consists of 132 observations of helmet acceleration (in units of g) measured along time after impact (in milliseconds). These data have being used as illustration of spline regression by Silverman (1985) and Eilers and Marx (1996) and are available in Härdle (1990). A history of the data set and its use in spline regression is given in de Souza and Heckman (2013). This dataset represents a good test for non parametric regression since the variance of the errors varies a great deal and there are several breakdown moments in the data. A-spline regression and BARS regression with order q = 3 are performed (Figures 4a and Figure 4b respectively). In Figure 4a, the solid lines represent the estimated t and the dashed lines represent the decomposition of the t onto the B-spline family. The two estimations are almost equal; the A-spline t seems visually smoother than BARS, especially around the points of high variation. Many non-parametric regression methods have been criticised for tting a function which increases slightly between 15 and 20 ms, as is the case for BARS, and (to a lesser extent) for A-spline. The rational is that since the helmet is subject to no exterior force before the impact (∼ 15 ms), the real signal is equal to zero before that time, and the tted function should not display variations in that region. However our work aims to nd a function that can simply represent a set of points, regardless of whether they respect some constraint about the real function. We could imagine extensions of spline regression that would incorporate such a priori knowledge on q q q q qq q qq q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -150 -100 (a) A-splines q q q q qq q qq q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -100 - q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0. A-splines of order 0 are tted (solid line) as well as the mean values tted using the PELT changepoint detection method (dashed line).

-
q q q q q q q q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0.75 -0.50 -0.25 0.00 400 500 600 700 A-splines of order 1 (solid line) and Multivariate Adaptive Regression Splines (dashed lines) are tted. the signal; this is beyond the scope of this article. A-spline regression has selected only 5 knots as relevant, and thus the tted function is a linear combination of 5 + 3 + 1 = 9

splines.

The second illustrative example uses a dataset of array Comparative Genomic Hybridization (aCGH) proles for 57 bladder tumor samples (see [START_REF] Stransky | Regional Copy NumberIndependent Deregulation of Transcription in Cancer[END_REF] references and access to the data). This dataset was used by Bleakley and Vert (2011) in the similar context of changepoint detection. The data represent the log-ratio of DNA quantities along 2215 probes. For the illustration, the 500 rst observations of individual 1's aCGH prole are used. We t a spline of order 0, i.e. a piecewise constant function.

Indeed, A-splines of order 0 perform a regression with changepoint detection of the data, which is a desired goal for these data. The tted spline is represented in solid line in Fig- ure 5. The estimated function performs a satisfying estimation of the changepoints and of the mean values over each interval. Our regression method estimated 9 changepoints, each corresponding to a shift in the mean value of the signal. Our method is compared to a popular changepoint detection algorithm (dashed line of Figure 5) called PELT [START_REF] Killick | Optimal detection of changepoints with a linear computational cost[END_REF]. We used the R package changepoint.np [START_REF] Haynes | Changepoint.np: Methods for Nonparametric Changepoint Detection[END_REF] . This method detects 8 changepoints, all of which correspond to a changepoint detected by the A-spline regression.

The third example is based on the LIDAR data [START_REF] Sigrist | Air Monitoring by Spectroscopic Techniques[END_REF][START_REF] Holst | Locally Weighted Least Squares Kernel Regression and Statistical Evaluation of LIDAR Measurements[END_REF], which is used by [START_REF] Ruppert | Semiparametric Regression[END_REF] to illustrate regression methods. The data come from a light detection and ranging (LIDAR) experiment. It consists of 221 observations of log-ratio of measured light intensities between two sources, as a function of the distance travelled by the light before being reected (in meters). The data are available in the R package SemiPar and are represented in Figure 6.

The y-variable looks to be linear on three intervals: it is almost constant before x = 550 m, decrease steeply between x = 550 m and x = 600 m, after which the slope increases. To highlight these shifts in slope, splines of order 1 (i.e. piecewise linear functions) are chosen to t the data. The A-spline t displays two slope changes, at x = 567 m and x = 607 m.

These moments visually correspond to the two biggest shifts in slope. We also t Friedman (1991)'s MARS procedure (in dashed line, Figure 6) and compare it to A-splines. We use an implementation of the procedure in the R package earth. This method also selects two breakpoints of the slope, at x = 558 and x = 612, which are very close to the breakpoints detected by A-splines.

The last example uses the data of the registered number of disasters in British coal mines per year between the years 1850 and 1962 [START_REF] Diggle | Equivalence of Smoothing Parameter Selectors in Density and Intensity Estimation[END_REF]. The number of coal disasters in each year is assumed to be Poisson distributed and the mean of the distribution is tted using a Poisson regression. The data are tted using A-spline regression of order 3.

The tted curve μ = g -1 (B â) is given in Figure 7. The 3 selected knots are represented by vertical dashed lines. The regression is now compared with P-splines (in dashed lines), which yields a similar estimation although less regularized.

Conclusion

In this paper we introduce a method called A-spline (for adaptive spline) performing spline regression which automatically selects the number and position of the knots. For that purpose, we set a large number of initial knots and use an iterative penalized likelihood approach (the adaptive ridge) to sequentially remove the unnecessary knots. The model achieving the best bias-variance tradeo is selected using the EBIC 0 .

Our method yields sparse models which are more interpretable than classical penalized spline regressions (e.g. P-splines). Other methods in the literature oer to jointly t the spline and select the knot's position. Through simulations, we highlight that our method performs almost as well as these methods, while being way less computationally intensive (by at least a factor 100).

Since knots correspond to a shift in the trend of the spline, their position is an important information about the data. To illustrate this, we apply A-splines to several real data sets.

When using A-spline with low order splines (e.g. 0 or 1), the approach allows performing changepoint detection. Indeed, A-splines of order 0 t a piecewise constant function to the data and hence detect changepoint in terms of mean. A-spline of order 1 ts a piecewise linear continuous function (i.e. a continuous broken line) that detects changepoints in terms of slope. For splines of degree 3, the knots indicate shifts in the third order derivative of the underlying function.

A fast implementation of A-spline is provided in R. Thanks to this, the computation of A-spline is very fast (∼ 1 sec for n ∼ 10000 and k ∼ 1000 on a standard laptop), even when tting generalized linear models with large sample sizes.

Our work can be naturally generalized to multivariate data using multidimensional Bsplines. Moreover, we limited our work to using B-splines for the sake of simplicity. But a variety of other splines can be used instead. For example M-splines, which are a basis of non-negative splines, could be used for tting non-negative functions (e.g. densities) and I-splines, which are a basis of monotonous splines, would yield a sparse isotonic regression model. Finally, our method can be used for non-parametric transformation of variables.

In particular, splines of order 0 could provide an automatic categorization of continuous covariates in regression models. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.001 0.100 50 100 200 400

(a) Logit q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.001 0.010 0.100 50 100 200 400

(b) Sine q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.001 0.010 0.100 50 100 200 400

(c) Bump q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.001 the knots selected by A-spline are mostly located at shifts when they exist in the signal: for the Bump function, the knots corresponds to places in the data with a lot variation, and for the Logit function, two knots are placed at the inection point.

In order to verify that A-spline's gain in interpretability is not detrimental to its predictive performance, we compare the Mean Square Errors between the two methods. Figure 2 shows the MSE for A-splines (solid lines) and P-splines (dotted lines) for every sample size and every function. It shows that for sample size 50, P-splines performs on average better than A-splines. For greater sample size however, A-splines performs almost as well as P-splines. These results are the same for every function. In conclusion, P-splines has better predictive performance for small sample sizes, but for data sets of size 200 and above,

A-splines and P-splines turn out to have comparable predictive performance.

Chapter 5

Segmentation of spatial data Crosses are water wells, dots are cholera outbursts.

Introduction

Spatial statistics plays a prominent role in epidemiology. The founding of epidemiology is said to occur with studies of spatial data of contamination. In 1854, the physician John Snow studies the geographical distribution of cholera outbreaks during the cholera epidemic of 1854 in the London neighborhood of Soho. He represented each outbreaks on a city map and over imposed the locations of water pumps in this neighborhood (see Figure 5.1). It enabled to identified a public water pump as the origin of the epidemic. This is the first record of a geographical analysis being used to epidemiology, and John Snow is considered to be the founding father of epidemiology.

Nowadays, the study of disease with respect to the geographical location of the patients has become widespread, and it plays a major role in epidemiology. Spatial (or geographical) epidemiology is based on the comprehension that there are geographical factors that are connected to a disease or to health indicators. It is interested in the detection and study of (i) "disease clusters" which are limited regions with unusually high values of a health indicator and (ii) "geographic correlation studies", which are the study of the effect of geographical location on the risk of disease onset. The data for these studies can be of various types, leading to different statistical tools to answer the epidemiological questions at hand. We refer to [START_REF] Andrew B Lawson | Handbook of Spatial Epidemiology[END_REF]; [START_REF] Lawson | Statistical Methods in Spatial Epidemiology. Wiley Series in Probability and Statistics[END_REF] for a review of the field and to [START_REF] Elliott | Spatial Epidemiology: Current Approaches and Future Challenges[END_REF]; [START_REF] Rezaeian | Geographical epidemiology, spatial analysis and geographical information systems: A multidisciplinary glossary[END_REF] for a summary of the tools and stakes thereof.

A regular problem in spatial statistics is the regularization of the spatial data. Most regularization methods perform a spatial smoothing of the data. These methods offer some advantages: (a) reduction of spatial covariance, (b) higher interpretability of the resulting map, and (c) interpolation. One of the most popular approaches to spatial smoothing is kriging, or Gaussian process regression (Cressie, 1993, Section 3). It assumes that the spatial data are the realization of a Gaussian process, with a specific covariance structure. The Gaussian process can be computed over the whole map using the covariance process at the observed points and the values of the process at the observed points. Since the Gaussian process is a Gaussian random variable at any point, the process's mean and variance can be estimated, and the whole distribution of the process is estimated. Computationally, this estimate requires the inversion of the covariance matrix, whose dimension is the sample size. This method is thus easy to compute and to interpret, and the assumption of an underlying Gaussian process is reasonable in many cases in absence of other information. The choice of covariance kernel tunes the estimated process: the smoother the covariance kernel, the smoother the estimated process. Other techniques of spatial smoothing have been developed; a popular approach is kernel smoothing [START_REF] Wand | Kernel Smoothing. Monographs on Statistics and Applied Probability[END_REF].

Spatial smoothing has many advantages. In some cases, one may want to obtain a segmented es-timation of the spatial distribution instead of a smoothed estimation, in cases where the underlying spatial effect is assumed to be discontinuous by nature. In many public health indicators, the spatial information is a proxy for a very located information, that is discontinuous by nature. In cities, the discretization into small neighborhoods is often very significant to tackle the spatial inequalities. When adjusted for professional category, income, wealth, and other quality-of-life variables, the neighborhood in which we live codes for the sole quality of life of each neighborhood provides, which is often a discontinuous variable. The same principle applies to demographic and epidemiologic studies where the quality of life is an important explanatory variable. This variable is of discrete nature and is indirectly present in the spatial distribution of the individuals. As an example, demographic studies on longevity focus namely on finding specific geographic areas where the longevity is unexpectedly high [START_REF] Michel Poulain | Identification of a geographic area characterized by extreme longevity in the Sardinia island: The AKEA study[END_REF]) -called blue zones. These areas are discrete by nature, and such studies use a spatial division into (administrative) areas to identify blue zones.

Our work is considering the case where the spatial distribution is indexed by a finite number of geographical units, called "areas" here (sometimes also called geographical sites). These areas represent a partition of the spatial domain. They often correspond to an administrative division of a territory, for instance census tracts, municipalities, or counties. Many studies present this type of data. This is due to the way data are collected by the administration (as is the case for demography) or to the imperative of privacy (as is the case in epidemiology). The areas are defined as polygons on R 2 . This work is built on solely using the planar graph of the adjacency structure of the areas as the information representing the spatial distribution. Doing so can seem to throw away information about the geographical proximity between areas, namely because two areas can be geographically close on average (i.e. for the Euclidean distance) and not be adjacent. However when the division into areas is regular enough, as is the case for most administrative unit types, the adjacency between areas is a good measure of proximity (see Section 5.2.1). Examples of studies based on the adjacency of areas are popular in epidemiology (van de [START_REF] Van De Kassteele | Estimating the prevalence of 26 health-related indicators at neighbourhood level in the Netherlands using structured additive regression[END_REF] and demographics [START_REF] Michel Poulain | Identification of a geographic area characterized by extreme longevity in the Sardinia island: The AKEA study[END_REF]. This field of spatial statistics has been widely studied, including modelizations using Gaussian distributions and Markov random fields (Cressie, 1993, see Sections 6.3 and 6.4). We refer to [START_REF] Cressie | Spatial Data Analysis of Regional Counts[END_REF] and [START_REF] Cressie | Spatial Modeling of Regional Variables[END_REF] for applications of such models to real-world data.

In this context, we introduce a new model for estimation of the spatial effect with segmentation based only on the adjacency of the areas. To the best of our knowledge, this problem has not been tackled in this context. We use a regularization method by penalized likelihood to perform segmentation over the graph. This approach can be seen as a generalization of the fused lasso-type regularization used in image denoising [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF], where the adjacency structure of pixels is a lattice. As in previous chapters, we use the Adaptive Ridge to penalize over the differences between adjacent areas. Our method yields a segmented estimate of the spatial effect in a computationally efficient way, because the complexity depends on the number of areas and not on the number of individuals. Particular attention has been brought to speeding up the computation, which is tractable up to a large number of areas. Our work is organized as follows. Section 5.2 introduces the model. In Section 5.3, we illustrate on simulated data that out method performs well to recover areas of constant spatial effect. Moreover, our method is shown to perform well also in the case where the true (unknown) effect has a smooth variation over space. Illustration with a real-data example of obesity prevalence in the Netherlands is provided in Section 5.4. The implementation of our method in R is publicly available (github.com/goepp/graphseg).

A model for spatial segmentation

Using graphical data for spatial segmentation

Our work is built on the following assumption: the areas provide a spatial discretization adapted to the sample distribution. This implies that the areas are not too small with respect to the overall spatial distribution, and not too large either. When the areas are too large (or equivalently when there are too few areas), most of the sample is located in few areas, and a lot of information is lost.

When the areas are too small (or equivalently when there are too many areas), many areas have no individuals and thus carry no information on the outcome variable. This case is not problematic in our work, since our objective is to find groups of adjacent areas with similar effects. However the total number of areas is a limiting factor for computation efficiency. Moreover, we assume that the areas are roughly equal in size, so that the discretization into geographical units is roughly uniform over space. In applications where there are several connected components (like islands), one can consider each component as a separate problem an perform segmentation independently from each other. If one wishes to introduce a notion of proximity between the connected components, one can connect every pair of components by artificially adding an edge between their closest areas (using the Euclidean distance). This case is not considered in the remainder of this work. In the forthcoming, we assume that the areas form a planar graph with only one connected component. Under these conditions, we can see that the distance between two areas is well approximated by its adjacency distance, that is the minimal number of other areas one needs to cross to join the two areas. This is illustrated with the following example. We use the départements of France as the geographical areas (data obtained at github.com/gregoiredavid/france-geojson, collected from the French national statistical institute (INSEE) and the French national geographical institute (IGN)). We exclude the overseas départements as well as the two départements forming the island of Corsica. The 94 départements are represented in Figure 5.2a, numbered using their codes. Since départements number 7 and 8 are removed, the numbering of the départements number 9 and more are offset by 2. We make the assumption that the distance given by the shortest path on the adjacency gives a good approximation of the geographic (i.e. Euclidean) distance between the areas -in practice, we will take the distances between the centroids of the areas. We compare the similarity matrices of these two distances, whose (i , j )-th entry is the distance between the areas i and j . There is a strong visual similarity between the two similarity matrices: see Figures 5.2b and 5.2c. We could test whether the two distances are similar in a sense to be defined. This is beyond the scope of this study.

Segmentation on a graph

Let n be the sample size and y i , i = 0, . . . , n be the response variable for the n individuals. Let D ∈ R 2 be the connected subset of R 2 representing the spatial domain of interest. Let A = (A j ) 1≤ j ≤p denote the areas and p denote the number of areas. The areas form a partition of D: ∀ j , A j ⊂ D, ∀ j , j , A j ∩ A j = , and ∪ j A j = D. Each individual is assumed to be in only one area:

n i =1 1 i ∈A j = 1. Let G = (V , E ) be
the adjacency graph of (A j ) j : each vertex v ∈ V corresponds to an area A j and there exists and edge between two edges (v, v ) if their corresponding areas have a border in common (areas connected by only one point are not said to be adjacent). The areas (A j ) are identified with their nodes V . The effect of each area over the response variable is noted γ = γ 1 , . . . , γ p . Our method applies to segmenting any area-based signal γ. When this signal is itself an estimate over the individuals in each area, our model includes its estimated variance σ2 = ( σ2 ) as a measure of the relative importance to give to each area.

We assume that γj ∼ N θ j , σ2 j for each area j . In this context of meta-analysis, γj are considered as observations and θ j as parameters. The negative log-likelihood simply writes:

(θ) = -log p j =1 1 2π σj exp - γj -θ j 2 2 σ2 j , ( 5.1)
In order to enforce the values of γ which are both adjacent and close in value to have close estimated θ, we use a penalized likelihood estimated defined as the minimizer of:

pen (θ, v ) = (θ) + κ 2 j ,k p(θ j -θ k ) (5.2)
with respect to θ. where the sum is taken over all adjacent nodes ( j , k), κ > 0 is the penalty constant, and p is a zero inducing penalty function (for instance a non-concave penalty, see Section 1.1.7).

Since we penalize over the differences of values between adjacent edges, many adjacent areas will be estimated to have the same value. Consequently, the penalized estimate is piecewise constant over few connected components of the graph. In this regard, this penalization procedure over a graph is related to the graphical lasso [START_REF] Friedman | Sparse Inverse Covariance Estimation with the Graphical Lasso[END_REF], which penalizes over the non-zero entries of the sensitivity matrix.

The adaptive ridge algorithm on a graph

To ensure segmentation of θ, we must take a zero-inducing penalty. We choose the L 0 Adaptive Ridge penalty, an iterative method developed by [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF] and [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF]. We define a matrix of weights v = (v j ,k ) 1≤ j ,k≤p between adjacent areas (the areas that are not adjacent have a weight v j ,k set to zero by convention). The L 0 adaptive ridge procedure is defined iteratively. The current estimate is defined as the minimizer of

pen (θ) = p j =1 (θ j -γj ) 2 2σ 2 j + κ j ,k v j ,k θ j -θ k 2 .
(5.

3)

The weights are initialized at the beginning of the iteration and are adapted at each step following the formula: .4) where ε > 0 is a small numerical constant. We iterate between estimating θ using (5.3) with fixed weights and adapting the weights using the previously estimated values of θ. At convergence, the weighted differences δ j ,k v j ,k (θ j -θ q ) 2 is almost equal to either zero (if the two values are set equal) or one (if the two values are set different). Using a cutoff threshold, the weighed differences are rounded to zero or one, which defines a set of connected components of the graph, that we call "regions". These connected component are the segmentation of the areas into groups of constant spatial effect. The adaptive ridge procedure performs a model selection amongst all the possible division of the graph into connected components. This effect θ is then estimated over each connected component by unpenalized estimation The weighted differences δ is used to diagnose the convergence of the algorithm. The whole procedure is given in Algorithm 6.

v j ,k = 1 θ j -θ k 2 + ε 2 . ( 5
δ old j ,k ← v j ,k (θ j -θ k ) 2 5:
while not converge do 6:

θ ← arg min θ pen (θ, v old ) 7: return θ 17: end function The last point is the numerical minimization of Equation 5.3 (Line 6 in Algorithm 6). We first reformulate (5.2) using the following property. Let L be the Laplacian matrix of the weighted adjacency graph, that is, L = D -V , where D is the degree matrix of the graph and V = (v j ,k ) j ,k is the matrix of weights between adjacent areas. Then the following identity (see [START_REF] Mohar | Some applications of Laplace eigenvalues of graphs[END_REF], Proposition 2.2) holds:

v j ,k ← θ j -θ k 2 + ε 2 -1 8: δ j ,k ← v j ,k (θ j -θ k ) 2 9: if δ -δ old <
j ,k v j ,k (θ j -θ k ) 2 = θ T Lθ
and thus (5.2) rewrites pen (θ) = W T (θ -Γ) 2 + κθ T Lθ. (5.5) Where

W = ( 1 2σ 2 j ) j
is a vector of weights of the areas (not to be mixed up with the weights between areas, noted v j ,k ). The minimizer of this quantity is explicit:

arg min θ pen (θ) = (W W T + κL) -1 W W T Γ. (5.6)
Consequently, the minimization of pen (θ) in Line 6 of Algorithm 6 is explicit. This significantly improves the computational speed of our method. Note that the computation of (5.6) is the bottleneck of our method, the complexity thereof depends on the dimension of L, that is, on p. Hence the limiting factor of our algorithm is the number of areas.

Remark. Our method is closely related to the problem of dividing an image into connected clusters of equal brightness. (This is not to be mistaken for the clustering of images into unconnected clusters of equal brightness, called image thresholding.) Consider a graph that is a lattice (or a square grid graph): each area has 4 neighboring areas (except those on the "border" and in the "corners", which have 3 and 2 neighboring areas, respectively). Since an image is a grid of pixels, the lattice is the adjacency structure of the image, where each area is a pixel. Then penalizing over the adjacent areas (see Equation 5.3) is like penalizing over the differences of adjacent pixels. This is the approach taken by penalty-based denoising methods, such as total variation (Rudin et al., 1992b). The image processing community usually uses a penalization of the L 1 norm of the differences of adjacent pixel values. But this choice is mainly motivated by the development of methods for solving this problem numerically, like the coordinate descent [START_REF] Friedman | Pathwise coordinate optimization[END_REF] or the primal-dual algorithm [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]. We could instead apply another penalty, like that of the L 0 adaptive ridge. Consequently, the present work generalizes directly to segmentation of images (into a piecewise constant image). This is however beyond the scope of this chapter.

Choice of the trade-off constant The penalized likelihood approach makes use of a penalty constant: κ. In this section, we tackle the choice of this tradeoff parameter. In practice, the choice of κ is done a posteriori: we run the penalized estimation for a series of values of κ, and then select the best value. We denote κ = (κ (1) , • • • , κ (L) ), the sequence of values of the penalty constant, in increasing order. As previously explained, the adaptive ridge performs model selection: to each κ (l ) corresponds one model. Several values of the penalty constant can induce the same model and according to the procedure of the adaptive ridge, if two models are the same, their estimate are the same. Consequently the criteria used for model selection depend only on the model dimension.

We use the Bayesian Information Criterion (BIC, see Schwarz, 1978):

B IC (m) = 2 ( θm ) + m log n, (5.7) where m is the model dimension (i.e., the number of estimated regions) and θm denotes here its corresponding estimate.

Remark. This criterion is shown to perform well in simulations and on real data. However, the BIC does not ensure that the selected models are of small enough dimensions, especially when p gets too large. This is explained here. The BIC is the criterion that maximizes the posterior probability of the model conditionally on the data, in the Bayesian framework. The BIC implicitly sets a uniform prior over all the models. Now consider what this prior implies on the dimension of the model: there are very few models of small dimension (close to 1) and of maximal dimension (close to p) and there are many models of medium dimension (close to p/2). This implies that the a priori probability of models of small dimensions is very small and consequently, the BIC will very probably not select very sparse models. When p gets large, the BIC fails to select sparse models (in practice, with ∼ 20 parameters), which is a desired goal. We note that other well-known model selection criteria are less penalizing than the BIC, so are of no use here. One can add a corrective prior term to the BIC which would assign equal probability to all model dimensions. This requires to enumerate the number of models for each dimension. That is: for k ∈ {1, • • • , p}, we want to compute the number of ways of dividing the planar graph G in k connected components. There seems to be no easy approach of this problem, which, interestingly, is not mentioned in the literature on planar graphs. Finding an exact formula seems out of reach. I believe there are good hopes to develop an algorithm that computes this enumeration. However (i) finding a polynomial-time algorithm seems much harder a problem (ii) even in that case, such an algorithm could take a long time to run on a large graph (p 1000). The study of this problem has been relegated as potential future work.

Visual inspection.

As an alternative to the BIC, we propose to use visual inspection to choose the best-fitting penalty. This does not provide an automated and non-arbitrary way of selecting the best model and as such, is not a satisfying method. Instead, we advise the practitioner to use visual inspection as a verification that the BIC selects a "coherent" model. Concretely, our method is applied on a sequence of penalty and one can consider that the resulting estimate is a sequence of segmentations, ranging from the least regularized to the most regularized. In this case, the practitioner can inspect if some segmented regions correspond to a sensible division or if they are isolated areas which

Simulation

We illustrate the model on several simulation settings. We use the spatial structure of iris (for Ilots Regroupés pour l'Information Statistique, i.e. "grouped islets for statistical information" ), a geographical division of the French territory used by the national French statistical institute. There are 16100 iris areas dividing the territory into roughly uniform areas geometrically and demographically. The iris data can be openly accessed from the French government's website 1 ; more information about iris can be found from the INSEE's website 2 . We choose iris because they are a good example of uniform division into administrative areas, as is the case in many problems in spatial statistics. Moreover, they are numerous enough to illustrate the numerical efficiency of our approach and the iris are connected (in the sense of graph connectivity). We limit the domain to the city of Paris, which comprises 987 iris of roughly equal size (see Figure 5.3). For the sake of the visual representation, we removed the 5 disproportionally large iris that constitute two forests east and west of Paris. We simulate a spatial distribution γ of different types, leading to different simulation settings. In both simulation settings, the variance term σ2 is set to 1 for each area.

First simulation setting. We defined the observed spatial effect to be piecewise constant over a small number of regions. These regions are the arrondissements of Paris. We merge the arrondissements 1 and 2 together as well as the 3rd and 4th arrondissements, so that there are a total of 18 regions of comparable size, noted (R k ) 1≤k≤18 in the increasing order of the arrondissement numbering. The aim of this simulation setting is to quantify how well our method can recover these regions.

The spatial term is defined as γ j = α j + ε j , (5.8) where α j = 18 k=1 1 A j ∈R k η k is the effect over A j , η = (η k ) 1≤k≤K is the vector of values taken by this function, and ε = (ε j ) 1≤ j ≤p is a vector of iid Gaussian errors:

ε j ∼ N (0, 0.25).

(5.9)

We set η = (4,2,4,[START_REF]Level set of the square log penalty 2 j =1 log(β 2 j + δ 2 )[END_REF][START_REF]4 Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ+α j +β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j[END_REF]7,[START_REF]4 Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ+α j +β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j[END_REF]5,[START_REF]4 Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ+α j +β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j[END_REF]1,[START_REF]Level set of the square log penalty 2 j =1 log(β 2 j + δ 2 )[END_REF]8,1,5,2,4,8,[START_REF]Level set of the square log penalty 2 j =1 log(β 2 j + δ 2 )[END_REF]. This piecewise constant "function" is illustrated in Figure 5.4a and the observed effect is given in Figure 5.4b. The resulting estimate obtained using our method is given in Figure 5.4c. Our method accurately estimates most regions R j . The limits between the regions are accurately estimated, except for a total number of 7 areas. However, the BIC tends to underestimate the total number of regions: it selects 12 regions but there are 18 regions.

To illustrate that the segmentation method performs well if the correct number of regions is correctly selected. We represent in Figure 5.4d the segmented estimate with a penalty constant κ chosen such that there are 18 selected regions. This estimate is close to the original piecewise constant signal; the only difference in terms of segmentation being that the arrondissements 9 and 18 are fused together, as well as 1/2 and 3/4. Since this information is not known a priori, this estimate cannot be used to assess the method's performance in selection. It is however used here to highlight that the correct segmentation can appear somewhere on the regularization path. This indicates that with additional work done on finding a model selection method more fitted to this problem, our method could show an even good performance in selecting regions of constant effect. We represent the regularization path of the estimation procedure in Figure 5.5. Each line represent the estimate value for on area across different values of κ, or rather for different values of number of estimated parameters, which is equivalent (the greater κ the fewer regions are estimated). The vertical line represents the model chosen by the BIC. The path illustrates the functioning of the adaptive ridge: for high values of κ, only one region is estimated, and as κ decreases, there are sudden shifts in the number of estimated regions. The breaking point where an additional region is estimated corresponds to a splitting of a colored line in two. The trajectories illustrate that the segmentations are encapsulated models: once two regions have been fused together (as κ increases), the region is never split again for higher values of κ. Consequently, the dendrogram provided that the regularization path forms is a hierarchical clustering of all the areas. Second simulation setting. In the first simulation setting, γ is constant over a set of regions and our method satisfyingly estimates these regions. In order to test the performance of our model it is necessary to evaluate it in a setting where the underlying signal is not piecewise constant. What it means for a signal over a graph to be "smoothed" is defined here. Signal processing on graphs have been developed using the shortest-path distance. We refer to [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF] for a review on this topic. The Laplacian of the graph plays a major role in the study of signals on graphs. This work naturally extends the results in signal processing of signals with discrete support. For example when the graph is a two-dimensional lattice, the Laplacian of the graph is the same as the Laplacian operator for discrete two-dimensional signals. We can apply filters to graph-based signals as for the signals indexed by integers (see Kalofolias, 2016, and references therein).

In this simulation, we define the spatial signal as the smoothing of γ (0) obtained by heat filtering [START_REF] Zhang | Graph spectral image smoothing using the heat kernel[END_REF]. This writes simply γ = exp(-sL)γ (0) , (5.10) where the exponential is the matrix exponential, L is the Laplacian of the graph, and s ∈ (0, ∞) is the smoothing parameter. We take s = 2 here. The resulting spatial effect is represented in Figure 5.6a.

The estimated segmentation is represented in Figure 5.6b. The BIC has selected 13 regions and their estimated effect θ has similar values. In fact, the "segmented" estimate (Figure 5.6b) displays a gradient of values and is close to the underlying smooth effect γ. This highlights that our method is robust to models where the true spatial effect is not in fact piecewise constant.

Real data application: Overweight prevalence in the Netherlands

We apply our method to the analysis of overweight prevalence indicators in the Netherlands. The data comprises 387, 195 persons in the Netherlands. The individuals are spread over a domain spanning 415 Dutch municipalities and divided into 11, 432 areas, which are neighborhoods. We have limited our studies to the region of Utrecht, comprising p = 2955 areas. The response variable is binary: it equals 1 is an individual is overweight and 0 if he is not. There are 11 explanatory variables (outside of the area): 6 qualitative variables (sex, ethnicity, marital status, household type, household income source, and home ownership) and 5 quantitative variables (age, household size, household capital, household income, and neighborhood urbanization). An additive binomial regression model is used to link η = logit(E y|x ) (the prevalence of being overweight) to the covariates and to the area. The model is a generalized additive model [START_REF] Wood | Generalized Additive Models: An Introduction with R, 2 edn[END_REF]:

η = X (1) T β (1) + p j =1 f j (x (2) j )β (2) j + p j =1 1 i ∈ j γ j ,
where ( f j ) are (known) spline functions, X (1) is the (binary-coded) design matrix for the 6 qualitative variables, X (2) is the design matrix for the spline regression for the 5 quantitative variables, and β (1) and β (2) are their corresponding variables. Finally, γ is the log-odds-ratio of the area location. Conditionally on the choice of the spline knots, the model writes

η = X β + p j =1
1 i ∈ j γ j , (5.11) where X = (X (1) , X (2) ) is the design matrix for all covariates except the area location. This model was fitted and I have used γ and its estimated variance σ2 . From the theory on logistic regression, we know that

σ2 j Var γj = n i =1 1 i ∈ j pi 1 -pi -1
, (5.12) where pi is the estimate probability that individual i is overweight. The computation of σ was made by people with access to the design matrix. From (5.12) we see that it is not possible to recover personal information from the estimated variance, as long as there are enough individuals in each area. Thus, ( γ, σ) can be communicated publicly without concerns of data privacy.

Remark. Define X , the extended design matrix obtained by column-concatenation of X (1) and the matrix with (i , j )-th entry 1 i ∈ j . Equation 5.11 is a linear model with design matrix X and parameter β = (β, γ). We could define a regularized estimate of γ by introducing a penalty term over γ in this linear model. We take another two-fold approach. We first estimate β from (5.11), considering γ as a nuisance parameter. We obtain an estimated γ corresponding to the spatial effects from the latter step. Then, we perform a regularization of γ using the previously estimated spatial effects γ, its estimated variance σ2 , and the adjacency graph. This has some advantages. The first one is regarding the privacy of the data. The second step of spatial segmentation can be performed without having to know the potentially sensitive data present in X . This methodology allows to apply the segmentation method to many studies where access to privacy-sensitive data is impossible, as is the case here. Results. The vector γ is represented in Figure 5.7a. The segmented estimate θ is represented in Figure 5.7b. We identify 40 regions of constant overweight prevalence using the BIC. The values with a higher risk (with respect to the overall risk) of being overweight are coded in red, the ones with the lower risk of being overweight are coded in blue. Our method highlights three disconnected regions with a specifically low overweight risk: two in light blue, in the central left and central right part of the map and in the center of the map one very large light blue zone, which includes three smaller regions with a very low overweight risk. The two light blue regions correspond to urban areas around the cities of Gouda (west-most region) and Arnhem (east-most region). The big blue region consists of a densely populated area spanning from Amsterdam to Utrecht. The three dark blue regions correspond almost exactly to the cities of Amsterdam (north-west), Hilversum (centernorth) and Utrecht (center-south). These three urban regions represent sites where the quality of life is better, since individuals are less likely to be overweight, even when adjusted for a series of important financial, social, and demographic variables. The red regions span across the west, east and (mostly) south of the Amsterdam-Utrecht axis. We identify two major dark red regions: one south of the city of Gouda (west on the map) and one in the south-east, the two being rural areas. There are also a few notable informations that can be seen from Figure 5.7b. There is a (light) red area in the eastern suburbs of Amsterdam, which seems like an important message for public health policy. Moreover there are few areas that have a different overweight prevalence than their neighboring areas. Further investigation by public health specialist need be carried about these specific neighborhoods to determine if there is indeed an abnormally increase risk of overweight or if this phenomenon is an artifact due to statistical uncertainty.

To conclude, the segmented estimate obtained by our method (Figure 5.7b) is not so different from the unregularized estimate (Figure 5.7b) and has some advantages:

• It defines few zones with different risks.

• It removes the arbitrariness of having to define what regions are part of a blue zone or a red zone.

• It provides a more strikingly visual message which is easier to communicate with.

Conclusion

Many public health studies now require precise geographic information to enforce policies and take action. This faces two issues: that of the privacy issues that go with sensitive (medical) personal information and that of finding models that make sense of the sometimes very geographically precise data. These two problems are solved by using a post-treatment of the spatial effect that estimates a segmentation into constant regions. We show that using a Bayesian framework, this approach only requires data aggregated at a unit of spatial resolution, called "area" (for instance, neighborhoods).

We develop an approach to perform this segmentation using the adjacency structure of the areas.

We show through simulations that (i) it performs well to recover the initial segmentation in the case where the underlying spatial signal is piecewise constant (ii) it is robust to other models for instance smoothly varying spatial signals. We found no record of another statistical method performing segmentation over an adjacency graph. Consequently, we had no reference method to compare to. We illustrated our method on the motivating example for this work: an application to the study of prevalence of overweight individuals in the Amsterdam-Utrecht area. We have used data from a pre-existing epidemiological study (van de [START_REF] Van De Kassteele | Estimating the prevalence of 26 health-related indicators at neighbourhood level in the Netherlands using structured additive regression[END_REF] and perform segmentation of the parameter coding for the spatial effect. Note that the personal data of individuals was not accessible in this case, and our method does not require to rerun the logistic regression of the original study. Instead, we perform segmentation of the estimated parameters.

The computation time is very sensible to the number of areas, which is the limiting factor. The method is computation-intensive due to important number of possible edges in a planar graph, which makes the adjacency matrix a computational burden. We highlight that our method scales relatively well with the number of areas in practical applications, up to p 1000. For that number of areas, the typical computation time is 7 minutes on a laptop equipped with an Intel Core™i3 CPU.

An R implementation of the method is available in the package graphseg3 . It makes use of sparse matrices and warm start to boost up computations. Some work is still needed for this method to be fully validated. Future works include a quantitative analysis of the performance in selection and estimation. An extension of the BIC to take into account the specifics of the model seems to be of interest. The EBIC (for Extended BIC Chen and Chen, 2008) is an extension of the BIC that includes a Bayesian prior on the model to give equal importance between the models of different dimension. This does not apply directly to our graphical model. We would require to compute the number of ways to divide a (planar) graph into k connected components, for varying values of k. Work in that direction has been undertaken, and would significantly add to the performance of our method.

Chapter 6

Conclusion The adaptive ridge and its applications

This thesis provides several applications of the adaptive ridge in age-period-cohort analysis, spline regression, and spatial epidemiology. When I started this thesis, the tentative title of my work was "Latent classes in survival analysis and heterogeneity of response to cancer treatment". The initial project was to study different approaches to tackle heterogeneity in survival analysis. My work was not initially meant to be focused on the adaptive ridge, but it evolved in this direction.

My work began with the study of the possibilities to extend the work of Bouaziz and Nuel (2017) to more general frameworks in survival analysis. Age-period-cohort analysis was our first choice to generalize the PCH model to bidimensional hazard. Age-period-cohort analysis is concerned with inference of follow-up data, where in most cases, regularization is needed. The factor models implicitly perform regularization of the log-hazard as they infer the effect of the age, period, and cohort variables. But they are too restrictive to infer the interaction between the effects of these variables. I have applied the fused adaptive ridge method to two models: the direct estimation of the log hazard (Chapter 2) and age-cohort-interaction model (Chapter 3). The former does not infer the effects of age, period, and cohort but instead performs segmentation of the log-hazard into regions of equal value. The latter provides an extension of the age-cohort model (or whatever two-variable model) and performs regularization of the interaction term.

The work on spline regression, in Chapter 3, was motivated by a discussion with G. Nuel. We thought that using a sparsity-inducing penalty on the coefficient of a spline would enforce nonrelevant knots to be selected out. At that time, I knew little about P-splines, which performs a ridge penalty over the difference of the coefficients. Only when implementing the adaptive ridge on this problem did I understand that it makes an iterated use of the same penalized likelihood estimation as the P-splines. The resulting estimate is, to my knowledge, the fastest procedure to perform selection of the knots.

Finally, I have had the occasion to apply the adaptive ridge on a widely different topic: spatial epidemiology (Chapter 5). This work was motivated by the application to the study of overweight prevalence (see Section 5.4). When discussing with the author of van de [START_REF] Van De Kassteele | Estimating the prevalence of 26 health-related indicators at neighbourhood level in the Netherlands using structured additive regression[END_REF], I realized that segmentation of the spatial effect would enable to detect a simpler and more interpretable spatial effect of the prevalence of being overweight. The data used in the generalized additive model in van de Kassteele et al. ( 2017) is sensitive and I could not have access to the design matrix. Thus I developed a method that used only the estimate for each area (i.e. spatial unit) and an estimate of its variance. This method performs a fused adaptive ridge over a graph: the graph of adjacency of the areas. This work can be seen as a generalization of the fused adaptive ridge over a grid, which is used in Chapters 2 and 3. I believe the method I developed, albeit somewhat simple, can be of great interest for the epidemiologist: a segmented spatial signal gives a clear message, which is important to take decisions of public health policy. 149 Chapter 6. Conclusion Why the adaptive ridge? The adaptive ridge is just one in a wide range possible model selection methods. The adaptive ridge was introduced under different names in different contexts, but has only been introduced as an iterative penalized maximum likelihood methods in [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF]. Contrarily to the adaptive lasso, it requires many iteration to provide a sparse estimate. It is however easier to implement, since its uses the L 2 norm penalty at each iteration, which in many cases has an explicit solution. I chose to use this method because it seemed like a good trade-off between computational cost and ease of implementation. However, no proof has yet been given of the asymptotic properties of the adaptive ridge. Thus, an important perspective of the present thesis is the theoretical study of the adaptive ridge. The next section presents the framework for this analysis.

Perspective: the adaptive ridge and its properties

In this Thesis I have used the adaptive ridge extensively. This method was first introduced as a numerical trick by [START_REF] Rippe | Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty[END_REF] and an investigation into its numerical performance was provided by [START_REF] Frommlet | An Adaptive Ridge Procedure for L0 Regularization[END_REF]. In this section, we lay the ground for a statistical study of this method. In a first part, we demonstrate that the L 0 adaptive ridge can be equivalently defined as the LQA of a penalized estimate, as defined by [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF]. The corresponding penalty is a function that has interesting properties. It can be seen as a relaxation of the log penalty, which was already introduced in the LLA one-step estimator by [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF]. This framework gives good insight into the strength of the adaptive ridge: it gives up the non-concave property in exchange for a gain in computational cost. Indeed, the oracle properties of non-concave penalties [START_REF] Fan | Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[END_REF] do not apply for the adaptive ridge because its function has an inflexion point.

Adaptive ridge as the LQA of a penalized estimate

Define the "square-log" penalty p(|θ|) = log(θ 2 + δ 2 ), (6.1) where δ 2 1 is a small numerical constant. This function is not a non-concave penalty (see Section 1.1.7): it is convex on (-δ, δ). However, this penalty converges pointwise towards the penalty log(θ 2 ) = 2 log(|θ|). Since δ is set to a very small value in practice, Equation 6.1 can be seen as an heuristic approximation of the log penalty.

The level set of the square log penalty in two dimensions is represented in figure 6.1.

Figure 1 . 1 :

 11 Figure 1.1: Visualization of the ridge estimate as the projection of the OLS onto an L 2 norm ball. The ellipses are level curves of the quadratic form y -X β 22 . The projection onto the circle of radius t = 1 has the effect of shrinking the estimate's coordinates together.

  Figure1.3: Illustration of the lasso estimation under orthogonal design. Projection onto two L 1 norm balls of radius t . As t decreases, the OLS will most probably end up in a greyed area and the lasso estimate will become sparse.

Figure 1 . 5 :

 15 Figure1.5: Bayesian priors on the parameter with the Ridge, Lasso, Elastic-net (α = 0.5), and Bridge (q = 0.5) penalties.

Figure 1 . 6 :

 16 Figure 1.6: Several non-concave penalties with different values of λ.

Figure 1 . 7 :

 17 Figure 1.7: Thresholding functions in orthogonal design.

Figure 1 . 8 :

 18 Figure 1.8: All possible submodels with three variables, forming a Hasse diagram of the power-set of {1, 2, 3}.

  Figure 1.9: Local linear approximation (a, in dashed line) and local quadratic approximation (b, in dashed line) of the bridge penalty (solid line) p(|θ|) = |θ| 0.5 around the current point θ (k) = 0.5.

Figure 1

 1 Figure1.11: Approximations of L q penalties by the Adaptive Ridge for different values of γ, and with q = 0 and ε = 10 -4 .

Figure 1

 1 Figure1.13: Unit ball for the L 1 /L 2 mixed norm with overlapping groups {{1, 3} , {2, 3}}.

  Figure 1.14: (a) Hierarchical structure of the overlapping group G = {{1, 2, 3} , {2, 3} , {3}}, (b) Unit ball of the corresponding mixed norm, (c) Projection of the unit ball on the plane x 2 = 0.

  Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.1.2 The piecewise constant hazard estimation . . . . . . . . . . . . . . . . . . . . . . 40 2.1.2.1 Proportional hazard model with piecewise constant hazard . . . . . . 41 2.1.3 Cohort data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.2 Regularized estimation of the hazard rate . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.3 Application to the evolution of breast cancer mortality . . . . . . . . . . . . . . . . . . 66

  Figure 2.1: Lexis Diagrams in the Age-Period plane (a) and in the Age-Cohort plane (b).

  Figure 1: Diagrams representing the lives of individuals: in the age-period plane (a) -called Lexis diagram -and in the age-cohort plane (b). Solid lines represent lives of individuals until occurrence of the event of interest. The same age, cohort, and period intervals are displayed in light gray. The intersection of two intervals forms a parallelogram and the intersection of three intervals forms a triangle.

Figure 2 :

 2 Figure2: Representation of the method used to select the constant areas for the adaptive ridge procedure. In this example, J = K = 3. In Panel (a), the circles represent the values of the differences v j,k (η j+1,kη j,k ) 2 and w j,k (η j,k+1η j,k ) 2 : empty circles correspond to the value 0 and filled circles correspond to the value 1. Panel (b) represents the graph that is generated from these values. Adjacent nodes whose difference is null are connected by a vertice. Panel (c) represents the last step, where the connected components of the graph are extracted. Each connected component corresponds to one constant area. The numbering is arbitrary.

Figure 3 :

 3 Figure 3: True hazard of the two simulation designs: smooth hazard in heatmap (a) and perspective plot (b) and piecewise constant hazard in heatmap (c) and perspective plot (d).

Figure 4 :Figure 5 :

 45 Figure 4: Smooth true hazard and corresponding estimates. The sample size is 4000 and the hazard estimates are medians taken over 500 simulations. The estimations are performed in the age-cohort plane and with different methods. Panel (a) represents the true hazard used to generate the data, Panel (b) represents the hazard estimated using the age-cohort model, Panel (c) represents the smoothed estimate, and Panel (d) represents the segmented estimate with the EBIC criterion.

Figure 6 :

 6 Figure 6: Estimated hazard of death after diagnosis of breast cancer for different stages of cancer. The estimate is obtained with the L 0 regularization. The upper right corner of every graph corresponds to the region where no data are available. Note that the scales are different between panels.

Figure 1 :

 1 Figure 1: Estimated hazard of death since diagnosis of breast cancer for different cancer stages and for different ages at diagnosis. The estimate is obtained with the L 0 regularization. The upper right corner of every graph corresponds to the region where no data are available. All graphs share the same scale.

Figure 2 :

 2 Figure 2: Estimated hazard of death since diagnosis of breast cancer for Stage 2 cancers. The estimation is carried separately for three classes of age at diagnosis: (., 45], (45, 55], and (55, .) and for sensitive and non-sensitive estrogen receptor cancers. Inference is made with the L 0 regularization. All graphs share the same scale.

Algorithm 5

 5 Adaptive Ridge Procedure 1: function ADAPTIVE-RIDGE(O, R, κ)

Figure 3 . 1 :

 31 Figure 3.1: True hazard (λ * j ,k ) for simulation design.

Figure 3 .

 3 Figure 3.2: Estimates of the ACI model in Simulation Design 1. (a) Maximum likelihood estimate (b) Estimated hazard (log λ j ,k ) in the ACI model (c) Age and cohort effects (µ + α j + β k ) in the ACI model (d) Interaction effect (δ j ,k ) from the ACI model (e) Age effects (α j ) (f) Cohort effects (β k ).

Figure 3 . 3 :

 33 Figure 3.3: Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ + α j + β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j ,k ) in the ACI model.

Figure 3 . 4 :

 34 Figure 3.4: Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ + α j + β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j ,k ) in the ACI model.

Figure 4 . 1 :

 41 Figure 4.1: Representation of the truncated power basis (T j ,k,t (x)) 1≤ j ≤k+q for k = 3 and t = (0.5, 1, 1.5).

  Figure 4.2: Illustration of the numerical instability of the decomposition in truncated power basis, with a spline of order 2 supported by 3 knots (tick line) and its decomposition with 2-digit machine precision coefficients (dotted line).

Figure 4 . 3 :

 43 Figure 4.3: B-spline bases of order 1 to 4 with three equally spaced knots t = c(0.25, 0.5, 0.75).

Figure 4 . 4 :

 44 Figure 4.4: Spline regression with 9 knots placed at equi-distance (c, d) and at the quantiles (a, b), for two different data sets.

Figure 1 :

 1 Figure 1: Bases of B-spline of degree 0 to 3 (Panels a to d) with 3 knots: (0.25, 0.5, 0.75).

  Splines tting It is interesting to note that A-splines are closely related to P-splines(Eilers and Marx, 1996), whose objective function writes PSS (a, λ) = SS(a) + λ 2 k+q+1 j=p+1 (∆ p a j ) 2 ,

  the posterior probability P(M λ |data) ∝ P(data|M λ )π(M λ ), where P(data|M λ ) is the integrated likelihood and π (M λ ) is the prior distribution on the model M λ . This problem is equivalent to minimizing -2 log P(M λ |data). By integration P(M λ |data) = a P(data|M λ , a)π(a)da, where P (data|M λ , a) is the likelihood and π(a) is the prior distribution of the parameter, which is taken constant in the following. Thus Bayesian criteria are dened as -2 log P (M λ |data) = SS(â λ ) + (q + k λ + 1) log n -2 log π(M λ ) + O P (1).

Figure 2 :

 2 Figure 2: Simulated data using functions Logit (a), Sine (b), Bump (c) and SpaHet (d), in solid line. Each dataset has size 200. The errors are chosen homoscedastic (σ = 0.15) for (a) and (b) and heteroscedastic (σ i = 0.3x i + 0.2 √ x i 2 ) for (c) and (d).

Figure 3 :

 3 Figure 3: Mean squared errors of A-spline (solid line), BARS (dashed line), and Spiriti's genetic optimization (dotted line) estimates for dierent sample sizes: 50, 100, 200, and 400. The simulations are performed with the Bump, Logit, Sine, and SpaHet functions and repeated 500 times.

Figure 4 :

 4 Figure 4: Motorcycle helmet crash data: helmet acceleration (unit of g) as a function of time (in ms). A-spline (a) regression and BARS (b) regression are tted. In Figure (a), the bold lines represent the estimates and grey lines represent the decomposition of the estimates onto the B-spline bases.

Figure 5 :

 5 Figure 5: aCGH data of bladder cancer: probes 1 through 500. A-splines of order 0 are

Figure 6 :

 6 Figure 6: LIDAR data: log-ratio of light intensity as a function of the travelled distance.

Figure 7 :

 7 Figure 7: Yearly number of coal accidents in Britain (grey bars) with P-splines regression (dashed curve) A-spline regression (solid curve). The three knots selected by A-splines are represented by vertical lines.

Figure 2 :

 2 Figure 2: Mean squared errors of A-spline (solid line) and P-spline (dashed line) estimates for dierent sample sizes: 50, 100, 200, and 400. The simulations are performed with the Bump, Logit, Sine, and SpaHet functions and repeated 500 times.

Figure 5 . 1 :

 51 Figure 5.1: Map of cholera outbreaks in the neighborhood of Soho (London) by John Snow (1854). Crosses are water wells, dots are cholera outbursts.

Figure 5

 5 Figure 5.2: (a) Division of continental France into 94 départements of roughly uniform size. (b) Euclidean distances between the centroids of each départements, using the latitude and longitude as coordinates. (c) Shortest path distances between départements, using the adjacency graph.

Algorithm 6

 6 Adaptive ridge procedure over a graph 1: function σ2 )

Figure 5 . 3 :

 53 Figure 5.3: Iris of the city of Paris used for simulation

Figure 5 . 4 :

 54 Figure 5.4: Result of simulation setting 1: piecewise constant effect α (a), noisy piecewise constant effect with additive Gaussian noise γ (b), segmented estimate using the L 0 adaptive ridge penalization θ with the BIC (c), and with 18 regions (d).

Figure 5 . 5 :

 55 Figure 5.5: Regularization path in simulation setting 1: trajectories of the estimated (θ j ) j for different values of κ. For a better visualization, the x axis gives the number of estimated parameters instead of κ. The vertical line represents the estimate selected by the BIC. Each color represents a different area.

Figure 5 . 6 :

 56 Figure 5.6: (a) Smooth spatial signal γ used in simulation setting 1 (b) and segmented estimate θ using the BIC.

Figure 5 .

 5 Figure 5.7: Unregularized (a) and segmented (b) values of overweight prevalence score in the Netherlands, adapted for 11 ethnic, financial, and social covariates.

Figure 6 . 1 :

 61 Figure 6.1: Level set of the square log penalty 2 j =1 log(β 2 j + δ 2 ) with δ 2 = 10 -2 .

Table 3 .

 3 1: Mean squared errors of the AC model, the ACI model, and the maximum likelihood estimate. For each sample size, the mean squared error was computed over 100 repetitions. The smallest mean squared error in each row is highlighted in bold.

  .[START_REF]4 Estimates in Simulation Design 2 with the ACI model: (a) True hazard, (b) AC model (log λ j ,k = µ+α j +β k ), (c) Estimated hazard in the ACI model, and (d) Interaction effect (δ j[END_REF]

			Order: 1				Order: 2
	1.00
	0.75
	0.50
	0.25
	0.00
			Order: 3				Order: 4
	1.00
	0.75
	0.50
	0.25
	0.00
	0.00	0.25	0.50	0.75	1.00 0.00	0.25	0.50	0.75	1.00

Table 1 :

 1 Mean squared errors of adaptive spline regression for dierent selection criteria and for dierent sample sizes. Dierent datasets are simulated using four dierent functions: the Bump function (a), the Logit Function (b), the Sine function (c) and the SpaHet function (d). The smallest value of each row is highlighted in bold. and the SpaHet for spatially heterogeneous function

	Sample size	AIC	BIC	EBIC	Sample size	AIC	BIC	EBIC
	50	0.02220	0.02	0.02418	50	0.02239	0.02001 0.02459
	100	0.00754	0.00324	0.00248	100	0.00755	0.00486	0.00458
	200	0.00285	0.00136	0.00127	200	0.00316	0.00231 0.00247
	400	0.00131	0.00071 0.00072	400	0.00156	0.00132 0.00141
		(a) Logit Function			(b) Sine Function
	Sample size	AIC	BIC	EBIC	Sample size	AIC	BIC	EBIC
	50	0.02000	0.01801 0.02211	50	0.02082	0.01784 0.02138
	100	0.00735	0.00627	0.00479	100	0.00727	0.00509	0.00371
	200	0.00354	0.00234	0.00217	200	0.00333	0.00194	0.00161
	400	0.00177	0.00106	0.001	400	0.00170	0.00081	8e -04
		(c) Bump Function			(d) SpaHet Function

  Using graphical data for spatial segmentation . . . . . . . . . . . . . . . . . . . . 5.2.2 Segmentation on a graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3 The adaptive ridge algorithm on a graph . . . . . . . . . . . . . . . . . . . . . . .

Contents 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 A model for spatial segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2.1

  10 -8 then

	10:	break
	11:	end if
	12:	θ old ← θ
	13:	v old ← v
	14:	δ old ← δ
	15:	end while
	16:

The regularization path of a regularized model is the path of solutions β(λ) represented as a function of λ. This path illustrates the progressive effect of λ on the fitted parameters and is often obtained over a grid of values of λ

http://github.com/Monneret/bandsolve

github.com/goepp/graphseg

Acknowledgement The authors are thankful to the National Cancer Institute for providing U.S. mortality data on cancer.

Conflict of Interest

The authors have declared no conflict of interest.

Appendix A: Expressions of the score and Hessian matrix

Unpenalized likelihood in the ACI model. We give the expressions of the (unpenalized) negative log-likelihood given in Equation 3.5. Its first order derivative is given by

Introduction to spline regression

Spline regression

We consider the problem of regressing a real-valued response variable y = (y i ) 1≤i ≤n onto a realvalued explanatory variable x = (x i ) 1≤i ≤n , where n is the sample size. We assume that there exists an unknown function f such that

where ε i are realizations of a centered random variable ε, not necessarily independent. The values of the x i s are assumed to be inside an interval [a, b]. Without loss of generality, we order the x variable, so that x 1 ≤ • • • ≤ x n . In this setting, the x variable is considered deterministic and the aim of regression is to find an approximation f of the f that links x to y. This problem being ill-posed, we must assume that f satisfies some additional property. This property is most often a level of smoothness, for instance, an order of differentiability.

Definition of splines

Splines are a parametric family of functions that are smooth and general enough to approximate a wide range of functions. Splines are defined using a series of knots, which are base points for the spline located on the x-axis. Let q ≥ 1 denote the number of knots. Let t = (t j ) define the knots in increasing order, such that a ≤ t 1 < • • • < t q ≤ b. Furthermore, we denote by k ≥ 1 the order of the spline. For j ∈ 1, • • • , q -1 , the splines of order k are the functions defined on [a, b] which are piecewise polynomials of degree k -1 on each interval [t j , t j +1 ) and whose k -2-th derivatives are defined and continuous at each knot t j (we use the convention that a function's derivative of a negative order is the function itself). For illustration, splines of order 1 are piecewise constant functions with possible discontinuities at each t j ; we easily verify that these functions are polynomials of degree 0 on each interval [t j , t j +1 ) and that they are possibly discontinuous at the knots. Splines of order 2 are broken lines, i.e. piecewise linear functions which are continuous but possibly non-differentiable at each knot t j . We can easily verify that these are the functions that are polynomials of degree 1 on each knot interval and are continuous everywhere. The set of splines of order k defined over the same knots are a linear space of dimension k + q (this property is easy to verify and we omit the proof). Moreover, splines are slightly more general functions than polynomials, which allows them to approximate a wider range of functions, in the method, which in turn performs slightly better than the A-spline. However these dierences are very small compared to the variation of the MSE inside each method, which remains important even for high sample sizes. Moreover, with sample size 800, the A-spline performs almost to as well as the genetic optimization method. Furthermore, we compare A-spline's predictive performance with P-splines through a similar simulation study, as shown in the Supplementary Materials. It results that P-splines have lower MSE on average, but the dierence between the two become negligible for large data sets (n ∼ 400). It is worth noticing that the comparison also illustrates how A-splines provide models that are easier to interpret.

The three methods compared here have very dierent computational costs. Simulations were carried on an Intel Core i3 CPU running at 2.3 GHz. For a sample size of 400, the median running time over 500 repetitions is 0.28 seconds for A-spline, while it is 2.66 seconds and 123.02 seconds for the BARS and genetic optimization methods respectively. This illustrates the computational eciency of A-spline compared to state-of-the-art knotselection methods.

In conclusion, the BARS and Spiriti's genetic optimization knot selection methods perform slightly better than A-spline, but this comes at the cost of inconveniently intensive computation.

Practical Implementation

In this section, the implementation of A-splines is explained in details. Particular attention has been brought to the computation of matrix products. Consequently, tting A-splines is almost instantaneous: on a standard laptop, it takes 1.3 seconds with 200 initial knots and 5000 data points. In the next three sections, several bottlenecks in the computation of A-splines are addressed. Matrix products computations are accelerated using an Rcpp [START_REF] Eddelbuettel | Seamless R and C++ Integration with Rcpp[END_REF] implementation. An R implementation of the A-spline estimation procedure is publicly available in the package aspline

Comparison of A-spline with P-spline

The following pages consist in the supplementary materials to the previous section. In this section, I provide additional simulations aimed at comparing the A-splines with the P-splines, both qualitatively and quantitatively. In this Appendix, the performance of A-splines is compared to a penalized spline regression method. For the sake of simplicity, we limit our study to the comparison with P-splines (Eilers and Marx, 1996), which is one of the most successful penalized spline regression methods. The goal of A-splines is t a model that is sparser and consequently more interpretable then that of P-splines. Therefore, A-splines are not required to have a better predictive performance than P-splines. This simulation study aims to illustrate that the loss in performance is acceptable when using A-splines instead of P-splines. We use the same simulation setting as in Section 6.2 of the article. We use the EBIC 0 criterion to select the penalty. Around zero, the square log penalty is a convex function and does not induce sparsity. For greater values of (β 1 , β 2 ), the square-log penalty is concentrated around the axes. This function has no singularity, so it enforces the penalized estimate to have coordinates set to almost zero, but not equal to zero.

Appendix to the article

We will show in this section that the Local Quadratic Approximation (see Section 1.3.2.1) of the square-log penalty is the L 0 adaptive ridge. Since the square-log penalty is not a non-concave function, the LQA of the square-log penalty is not proven to have the oracle properties that non-concave penalties enjoy. But writing the L 0 adaptive lasso as the MM Optimization procedure with the squarelog penalty allows to prove that the procedure is stable and to give a criterion to ensure that the procedure converges.

From the definition (1.42), the local quadratic approximating function of the square-log penalty writes

Hence the dominating function of the penalized nll writes

where

3), we can neglect this term in the minimization of g (β (k) j ). In the minimization step of the MM optimization, we minimize g (β|β (k) ). The criteria for the MM procedure is g k) ). In practice, we minimize g with the Newton-Raphson algorithm (or another second-order derivative-based method). These methods are unstable when the second-order derivative of (β) is too close to zero. It is therefore important to verify that the minimization step has decrease the function g . Equation 6.3 gives the following criteria for a candidate β (k+1) :

From the theory on MM optimization, if this criteria is met at each step, the procedure is stable and converges to a local minima. This framework motivates the study of the theoretical properties of the adaptive ridge.

List of Figures

1.1 Visualization of the ridge estimate as the projection of the OLS onto an L 2 norm ball. The ellipses are level curves of the quadratic form y -X β 2 2 . The projection onto the circle of radius t = 1 has the effect of shrinking the estimate's coordinates together. . . . 1.2 Illustration of L q unit balls for different values of q. . . . . . . . . . . . . . . . . . . . . . . 1.3 Illustration of the lasso estimation under orthogonal design. Projection onto two L 1 norm balls of radius t . As t decreases, the OLS will most probably end up in a greyed area and the lasso estimate will become sparse. . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Illustration of why the elastic net induces sparsity and not the L q norm (1 < q < 2). Unit radius balls of the elastic-net norm (α = 0.2) and the L 1.2 norm. . . . . . . . . . . . . . . . 1.5 Bayesian priors on the parameter with the Ridge, Lasso, Elastic-net (α = 0.5), and Bridge (q = 0.