Gerda Claeskens
email: gerda.claeskens@kuleuven.be

Emilie Devijver
email: emilie.devijver@univ-grenoble-alpes.fr

Irène Gijbels
email: irene.gijbels@kuleuven.be

Nonlinear mixed effects modeling and warping for functional data using B-splines

Keywords: B-spline approximation, Nonlinear functional modeling, Phase and amplitude variation, Random effects, Warping function

In functional data the interest is to find a global mean pattern, but also to capture the individual curve differences in phase and amplitude. This can be done conveniently by building in random effects on two levels: in the warping functions to account for individual phase variations; and in the linear structure to deal with individual amplitude variations. Via an appropriate choice of the warping function and B-spline approximations, estimation in the nonlinear mixed effects functional model is feasible, and does not require any prior knowledge on landmarks for the functional data. Sufficient and necessary conditions for identifiability of the flexible model are provided. A theoretical study is conducted: we establish asymptotic normality and consistency of the estimators of the registration and amplitude models, convergence of the iterative process, and consistency of the final estimator provided by the iterative process. The finite-sample performance of the proposed estimation procedure is investigated in a simulation study, which includes comparisons with existing methods. The added value of the developed method is further illustrated via the analysis of a real data example.

Introduction

Functional data are encountered in many fields, a multitude of examples can be found in the books by [START_REF] Ferraty | Nonparametric Functional Data Analysis[END_REF]; Ramsay andSilverman (2002, 2005). When analyzing functional data it is of particular interest to provide answers to the questions: (i) is there a common main (mean) functional pattern to be distinguished?; (ii) can we quantify the significant individual fluctuations with respect to such a mean pattern? While the common functional mean is capturing main features such as peaks and valleys, differences between individual curves are often exposed via differences in phase and in amplitude of the main features. In Figure 1 the Pinch Force data are depicted. These data were collected as part of an experiment to investigate the force (measured in Newton) exerted by thumb and forefinger when pinching a 6 cm width force meter. See [START_REF] Ramsay | A functional data analysis of the pinch force of human fingers[END_REF]. Data on 20 recordings of such force measurements, recorded every 2 milliseconds during a time period of 0.3 seconds, are presented in Figure 1. There seems to be a clear maximum for each curve, but the position and the size of this maximum differs considerably from curve to curve. See further Section 5. Often there is no prior information available regarding the number of important features, and where, in which region, they occur. A flexible method should thus not rely on such information, and be able to various elements of the modeling framework, and provide the identifiability of the model. All the proofs of the results stated in this section can be found in Appendix A.

A functional mixed model with warping function

We consider the following functional mixed model. For i = 1, . . . , n, and for t ∈ [0, 1], we define the process Y i (t) = µ w -1 (t; θ i ) + U i w -1 (t; θ i ) + ε i w -1 (t; θ i ) ,

(1) with µ the unknown common mean and where U i denotes the unknown random effect on the amplitude for the observation i. The flexible warping function w : [0, 1] → [0, 1] is strictly increasing and depends on a random variable θ i ∼ N r (θ 0 , Σ θ ), that describes the individual phase variability. Details about the warping function are provided in Section 2.3. We rather use the discretization of model ( 1) with time points (t i,j ) for j = 1, . . . , T i ; i = 1, . . . , n, where T i denotes the number of fixed (non-random) time points for the individual i: Y i (t i,j ) = µ w -1 (t i,j ; θ i ) + U i w -1 (t i,j ; θ i ) + ε i,j .

(2)

We assume that for all i, the error vectors ε i = (ε i,1 , . . . , ε i,Ti ) ⊤ with ε i ∼ N Ti (0 Ti , σ 2 ε I Ti ) are i.i.d., meaning that the error terms are independent of t i,j and of the warping effects θ i .

An analogous model was used in [START_REF] Rakêt | A nonlinear mixed-effects model for simultaneous smoothing and registration of functional data[END_REF], where the warping function stands only in the common mean and not in the individual effect, and in [START_REF] Gervini | Warped functional analysis of variance[END_REF], where a group level is added. We argue that from this general formulation, many things have to be defined to allow the estimation of this model. The specificities of our study and its novelty will be described in the next two subsections, through the decomposition of the signals onto a B-splines basis and the warping function.

B-spline basis decomposition

The warping function w, the unknown mean function µ and the individual random effect amplitude functions U i are modeled in a flexible fashion via B-splines. In this paper, we make the following assumption.

Assumption A. We assume that the functions µ, (U i ) i=1,...,n and w belong to the space spanned by the considered spline basis.

Assumption A ensures to have unbiased estimators for the curves, and avoids having to theoretically deal with a modeling bias. When using a spline basis in practice, the curves are well approximated when utilizing a finite (maybe large) number of knots.

For the mean function µ, we define a sequence of K µ interior knots 0 = κ µ 0 < κ µ 1 < . . . κ µ Kµ < κ µ Kµ+1 = 1. In addition, we put p µ + 1 boundary knots at 0 as well as at 1, and denote κ µ -pµ = . . . = κ µ -1 = κ µ 0 and κ µ Kµ+1 = κ µ Kµ+2 = . . . = κ µ Kµ+pµ+1 . We denote by κ µ = {κ µ -pµ , . . . , κ µ Kµ+pµ+1 } the set of all knots involved in estimation of µ. The B-spline basis functions of degree p µ are defined by induction as

B µ l,1 (t; κ µ ) = 1 if κ µ l ≤ t ≤ κ µ l+1 ; 0 otherwise; B µ l,pµ+1 (t; κ µ ) = t -κ µ l κ µ l+pµ -κ µ l B µ l,pµ (t; κ µ ) + κ µ l+pµ+1 -t κ µ l+pµ+1 -κ µ l+1
B µ l+1,pµ (t; κ µ );

for l = -p µ , . . . , K µ . With the use of the additional (boundary) knots, this gives precisely m µ = K µ + p µ + 1 basis functions. The function µ is decomposed in the B-spline basis, with coefficient vector α µ = (α µ -pµ , . . . , α µ Kµ ) ⊤ , µ(t) = Kµ l=-pµ α µ l B µ l,pµ+1 (t; κ µ ).

(3)

Note that if µ(•) does not belong to the space spanned by the basis functions, then the equality in (3) should be replaced by an approximation. The induced modeling bias can be controlled by taking a large number of knots.

Similarly, for i = 1, . . . , n each individual random function U i is decomposed in a basis of B-splines of degree p Ui . Denote the B-spline basis for U i by (B U i,-p U i ,p U i +1 , . . . , B U i,K U i ,p U i +1 ) with knots sequence κ Ui , resulting in m Ui = K Ui + p Ui + 1 basis functions, and consider

U i (t) = K U i l=-p U i α U i,l B U i,l,p U i +1 (t; κ Ui ),
where

α U i = (α U i,-p U i , . . . , α U i,K U i
) ⊤ . For this random vector α U i of B-spline coefficients, attaining values in R m U i ×1 , we assume that, for all i = 1, . . . , n,

α U i =    α U i,-p U i . . . α U i,K U i    ∼ N m U i 0 m U i , Σ Ui with Σ Ui =      σ 2 U,1 0 . . . 0 0 σ 2 U,2 . . . 0 0 0 . . . 0 0 . . . 0 σ 2 U,m U i     
, the covariance matrix for which we assume a diagonal structure, and which needs to be estimated. Further we denote α

∼ U = (α U 1 ) ⊤ , . . . , (α U n ) ⊤ ⊤ , a random vector taking values in R n i=1 m U i ×1 .

The warping function

A flexible way to model the warping function is as follows. For every t ∈ [0, 1], we define

w -1 (t; θ i ) = t 0 exp h -1 (u; θ i ) du 1 0 exp {h -1 (u; θ i )} du , (4)
with h -1 as indicated below. Note that w -1 (and hence w) is by construction an increasing function.

The advantage of using the exponential function is that it warrants the positivity of the function. A non-random version of this warping function was introduced in [START_REF] Ramsay | Functional Data Analysis[END_REF] and used in [START_REF] Hadjipantelis | Analysis of spike train data: A multivariate mixed effects model for phase and amplitude[END_REF]; [START_REF] Wagner | Nonparametric registration to low-dimensional function spaces[END_REF]. There are many other choices of warping functions that could be made (see for example [START_REF] Marron | Functional data analysis of amplitude and phase variation[END_REF]). In short, the warping function w -1 (or w) in (4) satisfies the following necessary conditions: increasing, and from [0, 1] to [0, 1]. To ensure identifiability, the function h -1 will be decomposed using a basis of centralized B-splines, i.e.

h -1 (u; θ i ) = K h l=-p h θ i,l Bh l,p h +1 (u; κ h ), (5)
where ( Bh l,p h +1 ) l=-p h ,...,K h satisfy

1 0 Bh l,p h +1 (u; κ h )du = 0.
The vector of random effects θ i = (θ i,-p h , . . . , θ i,K h ) ⊤ describes the individual phase variability, for which we assume a linear mixed effects model

θ i = θ 0 + E i + εi , (6)
with

E i ∼ N r (0 r , Σ E ) and εi ∼ N r (0 r , σ 2 εI r ) independent. Then θ i ∼ N r (θ 0 , Σ θ ), with Σ θ = Σ E + σ 2 εI r .
To ensure identifiability, we assume that σ 2 ε is known. In [START_REF] Gervini | Warped functional analysis of variance[END_REF] the parameter θ 0 is considered to be a Jupp transform of the landmarks of the mean function µ. In contrast, we avoid the use of landmarks, and θ 0 is a parameter to be estimated.

By construction, the warping function is injective, as proved in Lemma 1.

Lemma 1. The warping function w defined via (4) and ( 5) is injective with respect to the second parameter: for every t ∈ [0, 1],

t = w(w -1 (t; θ 1 ); θ 2 ) ⇒ θ 1 = θ 2 .
Further, we assume that the α U i s and θ i s, the random effects describing the individual phase and amplitude variability, are independent of each other.

The model in matrix form

For further analysis it will be useful to introduce some matrix notation. The matrix B µ i of dimension T i × m µ contains (j, l)th element B µ l,pµ+1 (t ij ; κ µ ), and

B U i is the matrix of dimension T i × m Ui with (j, l)th element B U i,l,p U i +1 (t ij ; κ Ui ). Further, (B µ i ) θi = ([(B µ i ) θi ] j,l ) j=1,...,Ti; l=-pµ,...,Kµ , with [(B µ i ) θi ] j,l = B µ l,pµ+1 {w -1 (t i,j ; θ i ); κ µ }. Define [(B µ i ) θi, θi ] j,l = B µ l,pµ+1 [w -1 {w(t i,j ; θ i ); θi }; κ µ ] for j = 1, . . . , T i and l = -p µ , . . . , K µ . Similarly, we define (B U i ) θi and [(B U i ) θi, θi ] j,l . The model (2) in matrix representation is Y ∼ = B ∼ µ α µ + B ∼ U α ∼ U + ε ∼ (7)
with

Y ∼ = (Y 1 (t 1,1 ), . . . , Y 1 (t 1,T1 )) ⊤ , . . . , (Y n (t n,1 ), . . . , Y n (t n,Tn )) ⊤ ⊤ ∈ R n i=1 Ti ×1 , ε ∼ = (ε 1 ) ⊤ , . . . , (ε n ) ⊤ ⊤ ∼ N n i=1 Ti (0 n i=1 Ti , σ 2 ε I n i=1 Ti ), α µ ∈ R mµ×1 , B ∼ µ = (B µ 1 ) θ1 ; . . . ; (B µ n ) θn ∈ R n i=1 Ti×mµ , α ∼ U ∼ N n i=1 m U i (0 n i=1 m U i , Σ ∼ U ), Σ ∼ U =      Σ U1 0 0 0 0 Σ U2 0 0 0 0 . . . 0 0 0 0 Σ Un      , B ∼ U = (B U 1 ) θ1 ; . . . ; (B U n ) θn ∈ R n i=1 Ti×mµ , θ ∼ = (θ 1 ) ⊤ , . . . , (θ n ) ⊤ ⊤ ∼ N r×n ( θ ∼ 0 , I r ⊗ Σ θ ), θ ∼ 0 ∈ R (r×n)×1 ,
where A ⊗ B denotes the Kronecker product of two matrices A and B.

In summary, the unknown parts in the model consist of

α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ .

Identifiability of the model

In this section, we provide sufficient and necessary conditions to ensure the identifiability of model ( 7).

First, the joint model ( 7) is identifiable if and only if at least one (approximate) individual model ( 2) is identifiable. We thus focus on a fixed i, and on the set of parameters (α µ , σ 2 ε , Σ Ui , θ 0 , Σ θ ) which consists of the subparts (α µ , σ 2 ε , Σ Ui ) and (θ 0 , Σ θ ), where the latter is linked to the warping modeling part, and the former with the other parts. We start by investigating identifiability in each part.

Identifiability of the warped process

Take any i ∈ {1, . . . , n}. For j = 1, . . . , T i , let X i be the warped process:

X i,j = Y i {w(t i,j ; θ i )} = µ(t i,j ) + U i (t i,j ) + ε i,j = Kµ l=-pµ α µ l B µ l,pµ+1 (t i,j ; κ µ ) + K U i l=-p U i α U i,l B U i,l,p U i +1 (t i,j ; κ Ui ) + ε i,j . (8)
Since ε i ∼ N Ti (0 Ti , σ 2 ε I Ti ), we obtain

X i |α U i ∼ N Ti (B µ i α µ + B U i α U i , σ 2 ε I Ti ).
First, remark that if we know (σ 2 ε , Σ Ui ), or if we know σ 2 ε , or if we know Σ Ui , model ( 8) is identifiable. In the following theorem, we give sufficient and necessary conditions for model (8) to be identifiable when both variance parameters are unknown. Since, for given (θ 0 , Σ θ ), and due to the use of B-spline approximations, the warped process leads to a linear mixed effects model, we can use general results on the identifiability of such models, as obtained by [START_REF] Wang | Identifiability of linear mixed effects models[END_REF]. Theorem 1 is an adaptation of Corollary 4.2 in [START_REF] Wang | Identifiability of linear mixed effects models[END_REF] to the current setting.

Theorem 1. Let i∈{1, . . . , n} be given. Model (8) is not identifiable if and only if the three conditions are fulfilled: 8) is identifiable if at least one of the three conditions in Theorem 1 is not satisfied.

1. (B U i ) ⊤ B U i ̸ = 0 m U i ; 2. H U i = B U i {(B U i ) ⊤ B U i } -1 (B U i ) ⊤ = I Ti ; 3. (B U i ) ⊤ B U i is diagonal. Consequently, model (

Identifiability of the warping function

Here, we assume that we know the parameters of the mixed effects model (α µ , σ 2 ε , Σ Ui ), and we want to prove the identifiability of the warping process

Y i (t) = X i {w -1 (t; θ i )}, (9)
involving the parameters (θ 0 , Σ θ ). Sufficient and necessary conditions for identifiability of this part are established in Theorem 2.

Theorem 2. Let i ∈ {1, . . . , n} be given. Let θ i ∼ N r (θ 0 , Σ θ ) and θi ∼ N r ( θ0 , Σ θ ) be used to define two warping functions w -1 (.; θ i ) and w -1 (.; θi ), and let X i and Xi be the corresponding warped functions, such that

Y i (t) = X i {w -1 (t; θ i )} = Xi {w -1 (t; θi )}.
Then model (9) is identifiable if and only if

B µ i = E θi, θi (B µ i ) θi, θi ; (B U i ) ⊤ Σ Ui B U i = Var θi, θi (B µ i ) θi, θi α µ + E θi, θi (B U i ) θi, θi ⊤ Σ Ui (B U i ) θi, θi .

Identifiability of the global model

We proved that, when knowing the warping parameters, the functional linear mixed effects model is identifiable, and that when knowing the functional linear mixed effects model, the warping parameters are identifiable. Then, by iterating between these two identifiable steps until convergence, we have a procedure which is identifiable and leads to the estimation of all the parameters of the model defined in (2). Note that the identifiability conditions are essentially conditions on the englobing B-spline basis structure.

Estimation procedure and asymptotic properties

Recall that the unknown parameters of model ( 2) are (α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ ). Model (2) is a nonlinear functional mixed effects model due to the composition by the warping function, which is an essential ingredient to describe the individual phase variability. First, we analyse each part of the model, that is, the warping parameters and the linear mixed effect model, by providing an estimator and theoretical guarantees. Then, we propose an iterative estimation procedure, where in a first step we fix the warping parameters (θ 0 , Σ θ ) and estimate the functional parameters (α µ , σ 2 ε , Σ ∼ U ); and next, we start from these estimated parameters, and estimate the warping parameters. Further, we obtain the convergence of the method and the consistency of the global estimator.

We have access to a sample (Y i (t i,j )) j=1,...,Ti; i=1,...,n of n curves, the ith curve being evaluated in T i points. Given are the knot sequences in the B-splines approximations, the degree of the B-splines, and the dimension of the warping parameters r = K h + p h + 1.

Proofs of the results in this section are provided in Appendix B.

Parameters of the warping function

Estimators for the parameters of the warping function

Suppose we know the functional parameters (α µ , Σ ∼ U , σ 2 ε ), and the predictors α U i for all i = 1, . . . , n. The goal is to estimate (θ 0 , Σ θ ).

We construct pseudo-observations by minimizing the following empirical L 2 criterion:

θTi i = argmin θi∈R r   Ti-1 j=1 Y i {w(t i,j ; θi )} -µ(t i,j ) -U i (t i,j ) 2 (t i,j+1 -t i,j )   . ( 10
)
Note that the criterion which is minimized tends to the L 2 distance between Y i • w(.; θi ) and µ

+ U i , if T i → +∞.
However, as we want to consider the warping parameter as a random effect, we fit a mixed effects model as defined in Eq. ( 6) on the pseudo-observations θT1 1 , . . . , θTn n , that is,

θTi i = θ 0 + E i + εi , (11)
with E i ∼ N r (0 r , Σ E ), εi ∼ N r (0 r , σ 2 εI r ) for all i = 1, . . . , n. The random variables E i and εi are independent. As we assume that σ 2 ε is known for identifiability reasons, we use the empirical mean of { θT1 1 , . . . , θTn n } to estimate θ 0 , and the empirical covariance to estimate Σ θ = Σ E + σ 2 εI r . The prediction of E i is easy to get because σ 2 ε is known. We consider the following estimators:

θ0 = 1 n n i=1 θTi i ; Σθ = 1 n -1 n i=1 ( θTi i -θ0 )( θTi i -θ0 ) ⊤ ; Êi = Σθ -σ 2 εI r
Σθ -1 θTi i .

Asymptotic normality of θTi

i First, we focus on the distribution of θTi i conditional on θ i . To do so, we rely on the theory of nonlinear least squares estimators developed in [START_REF] Jennrich | Asymptotic properties of non-linear least squares estimators[END_REF], which uses the weighted tail product defined as follows.

Definition 1. Let p be a nonnegative integer and (t j ) j=1,...,p be fixed time points. Let x = (x p ) and y = (y p ) be two sequences of real numbers and let

(x, y) π p = 1 p p-1 j=1
x j y j (t j+1 -t j ).

If (x, y) π p converges to a real number when p → +∞, its limit (x, y) π is called the weighted tail product of x and y.

Let g and h be two sequence valued functions on Θ.

If (g(α), h(β)) π p → (g(α), h(β)) π when p → +∞ uniformly for all α and β in Θ, we define [g, h] : (α, β) ∈ Θ × Θ → (g(α), h(β)) π .
This function is called the weighted tail cross product of g and h.

Then, we define the r × r -matrix a i ( θi ) as follows.

Definition 2. For l = 1, . . . , r, we denote by

∂ l (µ + U i ) = ∂[(µ(w -1 (t i,j ; θi )) + U i (w -1 (t i,j ; θi ))) j=1,...,Ti ] ∂[ θi ] l .
the partial derivative of the aligned signal. We define

a i,Ti ( θi ) = (∂ l (µ + U i ), ∂ l ′ (µ + U i ))
π Ti l=1,...,r;l ′ =1,...,r the matrix with coefficients the weighted tail product between two partial derivatives, and a i ( θi ) its limit when T i → +∞.

Assumption B.

For all i = 1, . . . , n, the r × r-matrix a i (θ i ) is non-singular.

Theorem 3. Fix i ∈ {1, . . . , n}, and α ∈ (0, 1). Let ( θTi i ) Ti be a sequence of weighted least squares estimators of θ i . We assume that the model is identifiable and satisfies Assumption A. With probability 1 -α, conditional on θ i , θTi i is a strongly consistent estimator of θ i (convergence a.s.). If we assume that the model satisfies Assumption B, conditional on θ i ,

T 1/2 i ( θTi i -θ i ) d -----→ Ti→+∞ N r (0, σ 2 ε a -1 i (θ i ));
and a i,Ti ( θTi i ) is a strongly consistent estimator of a i (θ i ).

Let f θT i i (.|θ i ) be the conditional distribution function. Denote by φ the Gaussian density function.

Theorem 3 implies that for all η ∈ R,

f θT i i (η|θ i ) → Ti→+∞ φ(η; θ i , σ 2 ε a -1 i (θ i )).
By the dominated convergence theorem, we get the asymptotic marginal distribution, for all η:

m θT i i (η) → Ti→+∞ m θ∞ i (η) = φ(η; θ i , σ 2 ε a -1 i (θ i ))φ(θ i ; θ 0 , Σ θ )dθ i = R r φ(θ i ; c, C)φ(θ 0 ; η, σ 2 ε a -1 i (θ i ) + Σ θ )dθ i
where θ∞ i = lim Ti→∞ θTi i , and

C = (σ -2 ε a i (θ i ) + (Σ θ ) -1 ) -1 ; c = σ -2 ε a i (θ i ) + (Σ θ ) -1 -1 σ -2 ε a i (θ i ) θ∞ i + (Σ θ ) -1 θ 0 .
The last line comes from a computation with Gaussian densities, see Lemma 4 in Appendix C. We discuss two cases where the limiting distribution m θ∞ i is computed explicitly:

• the case when the noise ε tends to disappear, which makes the theory easier, but also requires a strong assumption for the limit to hold;

• the case when the eigenvalues of the matrix a i (θ i ) are bounded. Under this weak assumption the limiting distribution will be more complicated (see below).

We next discuss these two cases in more detail.

In the first case we assume that the noise tends to disappear, when the number of points in the time grid increases.

Assumption C. We assume that σ ε → min Ti→∞ 0.
It is important to remark the following. If, however, there is a non-negligible noise, the method will warp the observed noise curve on some global mean, and the warping parameter will depend on this noise, whereas the true warping parameter does not, as it would be based on the denoised data.

Theorem 4. Fix i ∈ {1, . . . , n}. Let ( θTi i ) Ti be a sequence of weighted least squares estimators of θ i , and θ∞ i = lim Ti→∞ θTi i . We assume that the model is identifiable and satisfies Assumptions A, B and C.

Then, for all η ∈ R,

lim σ 2 ε →0 m θ∞ i (η) = φ(η; θ 0 , Σ θ ).
We now turn to the second case. With weaker assumptions, the limiting distribution m θ∞ i (η) is more complicate to describe. Denote by E α (θ 0 , Σ θ ) the following ellipsoid:

E ρ (θ 0 , Σ θ ) = x ∈ R r | (θ 0 -x) t (Σ θ ) -1 (θ 0 -x) ≤ ρ . Assumption D. For all θ i , a -1 i (θ i ) -a -1 i (θ 0 ) is positive definite.
Assumption E. The eigenvalues of a i are bounded: there exist λ m , λ M such that, for all θ i ,

λ m ≤ min eigen(a i (θ i )) ≤ max eigen(a i (θ i )) ≤ λ M .
The following theorem establishes that in this second setting, the limiting distribution m θ∞ i (η) is close to a Gaussian distribution.

Theorem 5. Fix i ∈ {1, . . . , n}. Let ( θTi i ) Ti be a sequence of weighted least squares estimators of θ i , and

θ∞ i = lim Ti→∞ θTi i .
We assume that the model is identifiable and satisfies Assumptions A, B, D and E. Let ρ > 0, and A a positive definite matrix. Then,

m θ∞ i (η) = (1 + O(ρ))φ(η; θ 0 , σ 2 ε a -1 i (θ 0 ) + Σ θ ) if η ∈ E ρ (θ 0 , A).
In summary, we get that under the two settings, the distribution of θ∞ i is close to a Gaussian distribution.

Asymptotic normality for θ0 and Σθ

We consider the linear mixed effect model given in Eq. ( 11). We assume that σ 2 ε is known for identifiability reasons. Remark that Σ θ = Σ E + σ 2 ε I Ti . We are now interested in the estimators θ0 and Σθ .

Theorem 6. Fix i ∈ {1, . . . , n}. Let ( θTi i ) Ti be a sequence of weighted least squares estimators of θ i . We assume that the model is identifiable and satisfies Assumptions A, B and C.

Let

b n,T = n i=1 T -1 i . Then, b -1/2 n,T θ0 -θ 0 d -→ n→+∞,min Ti→+∞ N r (0, Σ θ ); Σθ d -→ n→+∞,min Ti→+∞ W(Σ θ , n -1),
where W(Σ, p) denotes the Wishart distribution with scale matrix Σ and p degrees of freedom.

Note that this implies that T i has to go to infinity faster than n goes to infinity, i.e. n = o(min T i ). Indeed,

1 n i=1 1 Ti ≥ min T i n .

Functional parameters

Suppose we know the warping parameters (θ i ) i=1,...,n . Then, we warp the observations (Y i (t i,j )) j=1,...,Ti; i=1,...,n onto the estimated warped curves X i,j = Y i {w(t i,j ; θ i }, and we fit a functional linear mixed model on (X i ) i=1,...,n as defined in Eq. ( 8) using maximum likelihood estimation, which leads to estimators ( αµ , ΣU , (σ 2 ε )) and predictors ( αU i ) i=1,...,n . Following the ideas described in (Pinheiro, 1994, Chapter 3), we need the following assumption: Assumption F. Existence and positive definiteness of I, which is the limit of minus the expected Hessian matrix of the log-likelihood function based on the model given in Eq. (8).

Then, we get the asymptotic normality of the estimator.

Theorem 7. Let ( αµ , ΣU , (σ 2 ε )) be a sequence of maximum likelihood estimator of the functional linear mixed model, computed over the observations (Y i (t i,j )) j=1,...,Ti; i=1,...,n . We assume that the model is identifiable and satisfies Assumption F. Then,

√ n   αµ -α µ σU -σ U σε -σ ε   d -→ n→∞ N (0, I -1 ).

Global model and iterative estimation procedure

We propose to directly estimate the nonlinear model. Working with the L 2 -distance, we want to fit the model of which the coefficients minimize

Y - Kµ l=-pµ α µ l B µ l,pµ+1 {w -1 (.; θ); κ µ } - K U i l=-p U α U l B U l,p U i +1 {w -1 (.; θ); κ Ui } 2 .
Using the steps described previously, we propose an iterative process that approximates the following minimizer:

argmin α µ ,α U ,(θi)i=1,...,n    1 n n i=1 Ti-1 j=1   Y i (t i,j ) - Kµ l=-pµ α µ l B µ l,pµ+1 {w -1 (t i,j ; θ i ); κ µ } - K U i l=-p U i α U i,l B U i,l,p U i +1 {w -1 (t i,j ; θ i ); κ Ui }   2 (t i,j+1 -t i,j )     
.

Algorithm 1 presents the steps in the iterative procedure. Further details are provided regarded the initialization, the convergence criterion, the theoretical convergence and the consistency of the resulting estimator.

Algorithm 1 WarpMix

Initialization: Computation of the deepest function μ(0) = µ deep . for i = 1, . . . , n do Approximation of θ (0) i by θ (0) i = argmin θi∈R r   Ti-1 j=1 Y i {w(t i,j ; θ i )} -μ(0) (t i,j ) 2 (t i,j+1 -t i,j )   .
Fit a linear mixed model on the pseudo-observations, θ

(0) i = θ (0) 0 + E (0)
i + εi ; and deduce θ(0) 0 , ( Σθ ) (0) and θ(0

) i = θ(0) 0 + Ê(0) i . for ite = 1, . . . until convergence do
Warp the observed curves according to w -1

θ(ite-1) i ; Estimate {( αµ ) (ite) , ( αU i ) (ite) , ( ΣUi ) (ite) , (σ 2 ε ) (ite) } with the R package nlme; Approximate (θ (ite) 1 , . . . , θ (ite) n ) by computing, for every i = 1, . . . , n, θ (ite) i = argmin θi∈R r   Ti-1 j=1 {Y i {w(t i,j ; θ i )} -µ(t i,j ) -U i (t i,j )} 2 (t i,j+1 -t i,j )   .
Fit a linear mixed model on these observations: θ

(ite) i = θ(ite) 0 + Ê(ite) i + εi ; and define θ(ite) i = θ(ite) 0 + Ê(ite) i .

Details about the initialization

First, we initialize the mean function µ. There exist several ways to define a central curve in functional data analysis. Here we use band depth for functional data as introduced in Sun and Genton (2011), and compare every observed curve with the deepest function μ(0) = µ deep . We then deduce α(0) µ , the projection of the function μ(0) onto B µ the B-spline basis considered. In the initialization step, we do not consider individual amplitude effects, i.e. ( αU i ) (0) = 0 m U i for all i = 1, . . . , n.

Convergence of the algorithm

We define

C n = 1 n n i=1 Ti-1 j=1 Y i (t i,j ) -μ{w -1 (t i,j ; θi )} -Ûi {w -1 (t i,j ; θi )} 2 (t i,j+1 -t i,j ).
The iterations are stopped when C n < 10 -4 during five successive iteration steps.

Note first of all that the various iterations in Algorithm 1 involve three operations Ψ 1 , Ψ 2 and Ψ 3 , and that the update function to go from one iteration to the next is composed of three parts

Ψ = Ψ 3 • Ψ 2 • Ψ 1 : R mµ+nm U +1+nr → R mµ+nm U +1+nr . (12)
Herein Ψ and its components are defined as follows.

Ψ : ((α µ ) (ite) , (α U i ) (ite) i=1,...,n , (σ 2 ε ) (ite) , (θ i ) (ite) i=1,...,n ) → ((α µ ) (ite+1) , (α U i ) (ite+1) i=1,...,n , (σ 2 ε ) (ite+1) , (θ i ) (ite+1) i=1,...,n ) Ψ 1 : ((α µ ) (ite) , (α U i ) (ite) i=1,...,n , (σ 2 ε ) (ite) , (θ i ) (ite) i=1,...,n ) → ((α µ ) (ite) , (α U i ) (ite) i=1,...,n , (σ 2 ε ) (ite) , (θ i ) (ite+1) i=1,...,n ) Ψ 2 : ((α µ ) (ite) , (α U i ) (ite) i=1,...,n , (σ 2 ε ) (ite) , (θ i ) (ite+1) i=1,...,n ) → ((α µ ) (ite) , (α U i ) (ite) i=1,...,n , (σ 2 ε ) (ite) , (θ (ite+1) 0 + E (ite+1) i ) i=1,...,n ) Ψ 3 : ((α µ ) (ite) , (α U i ) (ite) i=1,...,n , (σ 2 ε ) (ite) , (θ (ite+1) 0 + E (ite+1) i ) i=1,...,n ) → ((α µ ) (ite+1) , (α U i ) (ite+1) i=1,...,n , (σ 2 ε ) (ite+1) , (θ i ) (ite+1) i=1,...,n ).
In Ψ 1 the vector θ i is updated. This is used as input for Ψ 2 where observations are denoised, through the linear model defined in (6). Then, this is used as input for Ψ 3 , where α µ , (α U i ) i=1,...,n , σ 2 ε are updated. In Theorem 8 we prove that the algorithm is converging. A condition under which this holds is that Ψ 1 is a contraction mapping, as stated in the following assumption.

Assumption G. There exists k Ψ1 < 1 such that, for (x, y)

∈ (R mµ+nm U +1+nr ) 2 , ∥Ψ 1 (x) -Ψ 1 (y)∥ 2 ≤ k Ψ1 ∥x -y∥ 2 .
Under Assumption G we show the convergence of the algorithm, seen as iterations of Ψ. We denote by (( αµ

) (∞) , (σ ε ) (∞) , ( Σ ∼ U ) (∞) , θ(∞) 0 , ( Σθ ) (∞)
) the estimator obtained at the end of the algorithm.

Theorem 8. Fix n and T. Suppose (Y 1 , . . . , Y n ) is a sequence of iid random variables satisfying the functional nonlinear mixed model (1) observed on fixed time points: for i = 1, . . . , n, for j = 1, . . . , T i ,

[Y i ] j = Y i (t i,j
). Moreover, suppose that the model is identifiable and satisfies Assumption G.

Then, (( αµ

) (∞) , (σ ε ) (∞) , ( Σ ∼ U ) (∞) , θ(∞) 0 , ( Σθ ) (∞)
) exists and is unique, and the algorithm converges to this solution with a geometric rate with respect to the Euclidean distance.

This theorem gives the pointwise convergence of the algorithm. The randomness has not been taken into account. We rather focus on the iterations of steps. Theorem 8 relies on Assumption G, which appears as rather technical. To get some insights into this assumption, we investigate it in a specific setting in Example 1.

Example 1. We focus on µ, do not consider U i , and restrict the family of warping functions to translations: w -1 (t; θ i ) = θ i + t. The global mean is supposed to be a linear function µ(t) = α + βt. Let t i,1 = 0 and t i,Ti = 1 Finally, we set θ 0 = 0.

Fix i. Recalling (10), in this case we are looking for

θ i = argmin θi    1 T i Ti-1 j=1 (α -α + (β -β)t i,j -β θi ) 2 (t i,j+1 -t i,j )    .
This is a polynomial function of degree 2 in θi with nonnegative coefficient of the quadratic term: there exists a unique minimizer:

θi = α - α β + β - β β Ti-1 j=1 t i,j (t i,j+1 -t i,j ).
We know that the Lipschitz constant is bounded by the norm of the differential. Here, the function we consider is (α, β) → θ, so we compute the differential, evaluated in (α, β):

∥D α, β Ψ 1 ∥ 2 2 = - 1 β 2 +   -α + α β2 + -β β2 Ti-1 j=1 t i,j (t i,j+1 -t i,j )   2 ; ∥D α,β Ψ 1 ∥ 2 2 = 1 β 2   1 + { Ti-1 j=1 t i,j (t i,j+1 -t i,j )} 2   .
These expressions reveal that for β small, the problem is more complicated (as one could expect). Note that in this special case Assumption G in fact leads to an assumption on β.

Example 1 also shows that, in some particular settings, Assumption G might be translated into a condition on µ and w.

Consistency of the estimators

To conclude, we provide the statistical consistency of the full procedure. This has the following meaning. When the sample size and the number of time points are going to infinity, the parameters estimated by the iterative process are converging almost-surely to the true parameter. Finally, the consistency is deduced for the common mean, seen as a functional parameter.

Theorem 9. Suppose (Y 1 , . . . , Y n ) is a sequence of iid random variables satisfying the functional nonlinear mixed model (1) observed on fixed time points:

for i = 1, . . . , n, for j = 1, . . . , T i , [Y i ] j = Y i (t i,j ).
Suppose that the model is identifiable, and Assumptions A, B,C, F and G hold. Then,

(( αµ ) (∞) , (σ ε ) (∞) , ( Σ ∼ U ) (∞) , θ(∞) 0 , ( Σθ ) (∞) ) a.s. -→ n→∞ min T i →∞ (α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ ).
As a consequence, we get that, from a functional viewpoint, for µ ∈ span(B µ ), if we denote μ = ( αµ

) (∞) B µ , ∥µ -μ∥ L2[0,1] a.s.
--→ n→∞ min T i →∞ 0.

Simulation study

We investigate the finite-sample performance of the proposed estimation method, and we compare it with four state-of-the-art methods, described below. An R package, called warpMix, has been developed for the proposed method and is available at https://cran.r-project.org/web/packages/warpMix/ index.html.

Description of the simulation settings

Warping functions

The warping process is the same in most of the settings (with exception of Model M2), and defined via (4) and (5). Three different interior knots {0.2, 0.5, 0.7} are used for a basis of cubic splines for h -1 . So in this setting κ h = 3, p h = 3 and r = κ h + p h + 1 = 7.

The random variables (θ i ) i=1,...,n are distributed according to N r (0 r , Σ θ ), with Σ θ = Σ E + σ 2 εI, where σ 2 ε = 10 -3 , and Σ E a diagonal matrix with elements {2, 0.8, 0.4, 0.3, 0.4, 0.8, 2}.

Figure 2 depicts a sample of size 100 of the warping function, and the empirical covariance matrix of (w -1 (.; θ i )) i=1,...,n computed on this sample. This highlights the differences between the correlation in (θ i ) i=1,...,n and that in (w -1 (.; θ i )) i=1,...,n . Note that due to the random warping structure, there is a variability induced on the whole time period.

Elements of the functional model

The elements determining the functional model are the function µ and the individual random effects U i . For the mean function µ we consider four different functions. For t ∈ [0, 1], µ 1 (t) = {sin(3πt) + 3πt}/4, µ 2 (t) = exp -(t-0.25) 2 /0.04 + exp -(t-0.75) 2 /0.02 , µ 3 (t) = cos(2πt + π/2), µ 4 (t) = sin(6πt)/(6πt).

These functions are plotted in Figure 3. The modeling framework in Section 2 assumes that the functions µ and U i are well approximated using a B-spline basis. This is in practice not always the case, for example when a too limited number of knots is considered in the B-spline bases. In the simulation study we present results on the B-spline approximations of the µ-functions, denoted by μ1 , μ2 , μ3 and μ4 respectively. To illustrate the impact of modeling bias, we provide for the fourth function simulation results for its B-spline approximation μ4 as well as for the function µ 4 itself. We refer to model ( 2 to refer to the highdimensional sample setting.

The estimators of µ and U are computed using quadratic splines (p µ = p Ui = 2), with interior knots at {0.12, 0.24, 0.36, 0.48, 0.60, 0.72, 0.84}, so that K µ = K Ui = 7, and m µ = m Ui = 7 + 2 + 1 = 10. The individual effects U i in the functional model have a centered multivariate normal distribution with diagonal isotropic covariance matrix Σ U = 0.1I 10 , except for models M 4 and M HD 4 which are harder to fit, where we use Σ U = 0.05I 10 . The variance of ε in the functional linear model equals σ 2 ε = 0.02. In the numerical study, we simulate 100 times from each setting, and report the evaluation criteria based on these 100 simulated samples.

Variability

To generate the data, we first construct a sample of the process (X 1 , . . . , X n ), defined in ( 8) and then un-warp them via (9) and the warping function described in Section 4.1.1. To understand the variability induced by each modeling aspect, we plot in Figure 4, a warped sample and the un-warped sample for model M 1 .

The signal-to-noise ratio expresses the ratio of the variability caused by the signal [µ{w -1 (t i,j ; θ i )} + U i {w -1 (t i,j ; θ i )}] and that due to the noise ε i,j

SNR(t i,j ) = Var µ{w -1 (t i,j ; θ i )} + U i {w -1 (t i,j ; θ i )} Var(ε i,j ) .
To compute the numerator we use the conditional variance formula, for a random variable V seen as a function of two random variables U and θ,

Var (θ,U ) (V ) = E θ {Var U (V |θ)} + Var θ {E U (V |θ)}.
For each given time point t i,j we compute this SNR function 50 times to get 50 values for SNR at each time point. To compute the function once, we proceed as follows. For a fixed θ, we compute the empirical conditional variance Var U (Z|θ) and the conditional expectation E U (Z|θ) over a sample of size 100. By varying θ 60 times, we compute the global variance. This whole process is then repeated 50 times. In Figure 5 the resulting approximations for the SNR functions for models M 1 and M 2 are plotted. For higher values of SNR we expect the estimation problem to be somewhat easier. Some caution regarding this interpretation is needed though. In our functional mixed effects model there are several sources of variability in the signal part (the individual effect related to U i and the warping effect due to θ i ).

The SNR-criterion does not distinguish between these variabilities, and just considers the global signal variability against the error variability. Note from the SNR plots in Figure 5 that the estimation task can be harder in some time-regions. At the endpoints of the interval [0, 1] the SNR-values for the different models are equal, since the warping is not effecting these parts, and the only effect is coming from the covariance matrix Σ U , the noise variance σ 2 ε , and their relative contribution.

Comparison with existing methods and performance criteria

To illustrate the numerical performance of the proposed method, we compare with four methods available in the literature. Since our nonlinear functional mixed effects model is closely related to that of Gervini and Carter (2014) with major differences as indicated in Section 1, we include a comparison with this method. Some procedure parameters have to be chosen in the method of [START_REF] Gervini | Warped functional analysis of variance[END_REF]: we took p = q = 1 for the number of components in the Karhunen-Loève decompositions; λ = 1 for the regularization parameter; and τ 0 = {0.3, 0.6} as the set of average landmarks. Their estimation method involves a Monte Carlo approximation part, for which we considered 100 iterations; and an EM algorithm part in which we also considered at most 100 iterations. Convergence was said to be reached when the difference in norm between estimated parameters in two consecutive iteration steps was less than 10 -2 . We also would like to mention that in our simulation study we use a rewritten Matlab version of the original Table 1: Evaluation criteria for the estimation tasks. [START_REF] Gervini | Warped functional analysis of variance[END_REF], since the latter was no longer running properly. The use of the Matlab code can make computations a bit slower.

quantity of interest evaluation criterion µ ∆ µ = T -1 j=1 [μ(t j ) -µ(t j )] 2 (t j+1 -t j ) w ∆ w = T -1 j=1 [w -1 (t j ; θ0 ) -t j ] 2 (t j+1 -t j ) Σ U ∆ U = Tr( ΣU -Σ U ) Σ θ ∆ θ = Tr( Σθ -Σ θ ) σ 2 ε ∆ ε = |σ 2 ε -σ 2 ε | code used in
The elastic square-root slope is a promising framework, so we include a comparison with the method developed in [START_REF] Tucker | Generative models for functional data using phase and amplitude separation[END_REF]. We use the default settings: no elasticity, Karcher mean, do not smooth the data and at most 20 iterations. We use the code available in the R package fdasrvf.

Bayesian methods are also of interest, and we choose to compare with [START_REF] Cheng | Bayesian registration of functions and curves[END_REF], also available in the R package fdasrvf. Also here we considered the default settings: 150000 iterations and a uniform prior distribution.

Finally, we compare the performances with that of the algorithm of [START_REF] Sangalli | k-mean alignment for curve clustering[END_REF], available in the R package fdakma, that allows for clustering misaligned data. We assume that there is no cluster, consider affine alignment and compute the similarity through the cosine of the angle between the two function.

Since the available inference in those studies does not fully match our modeling inference, we can only report on the comparison related to estimating µ.

To evaluate the estimation performance for the various components of the target (α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ ) we need some criteria. Note that the modeling framework involves two unknown functions, namely the overall mean function µ and the warping function w, unknown matrices Σ U , Σ θ , as well as the unknown variance σ 2 ε . For each sample we obtain estimates μ, ŵ, ΣU , Σθ and σ2 ε . Since in our simulation setting we have the same observational time points for each individual curve, i.e. t i,j = t j , and j = 1, . . . , T i , with T i = T , we use the criteria in Table 1 to evaluate the estimation performance in each sample. Herein Tr(A) denotes the trace of a matrix A.

For each simulated sample we calculate the estimates, and the corresponding evaluation criteria of Table 1. To report on the bias of an estimator, we compute the empirical mean of a criterion over the 100 simulations. To report on the variance of an estimator, we proceed as follows. For example, when estimating the function µ, we calculate in each point t j the empirical mean over all 100 estimated values of µ(t j ) and denoting this by μ(t j ). For each simulated sample we then calculate ∆µ = T -1 j=1 [μ(t j )μ(t j )] 2 (t j+1 -t j ). The empirical variance of the estimator for µ is then computed by taking the average over the 100 obtained ∆µ values. In a similar way we obtain ∆w , ∆U , ∆θ and ∆ε .

A final remark is that for ∆ θ and ∆θ , we use medians rather than means across all simulations as a measure of central position, since sometimes estimation of some components of Σ θ resulted in large outlying values. However, even in the latter cases the quality of the estimated warping function was still (n = 200 and T = 150, bottom) using the proposed method. 0.0020 0.0019 0.0007 0.0009 very good, as will be seen from the reported results.

Simulation results for the proposed method

Models M 1 and M HD 1 . Figure 6 depicts the simulation results for estimating μ1 in M 1 and M HD 1 . In the left panels we depict, for each time point t j , the boxplots of the obtained estimated values for μ1 (t j ), whereas in the right panels we use a functional boxplot, as developed in [START_REF] Sun | Functional boxplots[END_REF]. The true curve μ1 is in all plots presented as the solid (red) curve. The black solid curve in the centre of the functional boxplots indicates the deepest function among all estimated mean functions.

The quality of estimating μ1 is quite good for the proposed method. Passing from low dimension to high dimension (from the top row to the bottom row plots), we see that the results improve for larger values of n and T . Note that the largest variability occurs in the region where there was also most variability noticed in the SNR plot for model M 1 in Figure 5. Table 2 further summarizes the simulation results for models M 1 and M HD 1 . The results on estimation of µ are in correspondence with what was observed from Figure 6. Note that in estimation of Σ θ there are quite some extreme estimation results. However, the resulting estimation of the warping function w is still good, as can be noticed from the last row in Table 2. In Figure 7 we present boxplots of the estimation results for the components of Σ θ for models M 1 (left) and M HD 1

(right), with the true component values indicated as red horizontal lines. Outliers have been excluded for plotting the boxplots for clarity of presentation. To complement these boxplots, we summarize in Table 3 the maximum (across simulations) of the estimated values for each of the seven components of Σ θ . Note that the most extreme values occur for the first coefficient. For larger n and T there are less extreme estimates.

Next, we focus on estimating the individual curve amplitude variability, which is captured by the estimation of the ten diagonal components of Σ U . In Figure 8 we provide boxplots of the estimation results. The horizontal red line presents the true value 0.10 for all diagonal components. As can be seen the estimation results tend to be better for larger value of n and T , as expected.

Models M

w CDG 2 and M w GC 2

. To study the finite-sample performance of the proposed estimation method when there is a misspecification with respect to the warping function, we consider the model with μ2 and simulate data under two different warping schemes:

• scheme w CDG : the warping scheme of Section 2.3;

• scheme w GC : the warping scheme of [START_REF] Gervini | Warped functional analysis of variance[END_REF]. In the scheme w GC , θ i are generated via a linear mixed effects model, and then mapped into the set of landmarks using a Jupp transform. This is followed by interpolation by cubic splines to get to the corresponding parameters. Simulations were carried out from the two models, referred to as models M wCDG 2 and M wGC 2

. Table 4 summarizes the simulation results for all elements in the functional mixed effects model. Overall conclusions remain as above. Note that also under the misspecified warping scheme w GC the proposed method continues to perform very well. q q q q q q q q q q 1 2 3 4 5 6 7 8 9 10 0.00 0.05 0.10 0.15 0.20 Boxplots for Σ U q q q q q q q q q q q q q q q q q 1 2 3 4 5 5 presents the simulation results for the low-and high-dimensional settings for model M 3 . Also in these settings the method performs well. For the fourth model we include simulation results (in the low-dimensional sample setting) when simulating from the unprojected function µ 4 , for which the B-spline approximation induces a modeling bias. As can be seen from columns 2-5 in Table 6 there is only a little loss in performance when modeling bias is present. . The simulation results are summarized in Table 7. Note that the proposed method often has low/lowest bias, but at the price of having a larger estimation variance. On model M 1 , the method fdakma performs the best, with a very low variance, but it has a comparable performance (in terms of bias) to the proposed method for M HD 1 . On models M wCDG 2 and M wGC 2

, our method has particularly good results in mean, but with a larger variance. As GC's method is very slow (i.e high computational cost) and does not provide very good results whereas the modelling is close to the proposed one, we restrict further comparisons, for Models M 3 , M HD 3 , M 4 and M HD 4 , to the other three methods. On all models M 3 , M HD 4 , M 4 and M HD 4 , fdakma performs the best among the competitive methods, followed by the Bayesian warping method (method (3)), but both are less good than the proposed method in terms of bias. Finally, we see that the elastic square-root slope method (method (4) in the table) does not perform well on those simulated datasets.

Real data analysis

We analyze the Pinch Force dataset, available in the R package fda. These data were described and analyzed in [START_REF] Ramsay | A functional data analysis of the pinch force of human fingers[END_REF]. The data consist of measurements, at every second millisecond, on the exerted force (in Newton) during a period of 0.3 seconds. The resulting measurements consist of 20 curves recorded on 151 points in [0, 0.3]. See Figure 1. For convenience the data were rescaled to the domain [0, 1].

We analyzed these data, using B-splines of degree 2 for µ and U (i.e. p µ = p Ui = 2), with interior knots {0.25, 0.50, 0.75}, resulting in six B-spline basis functions. For the function determining the warping in (5) we use B-splines of degree 3 (i.e. p h = 3) and the same set of interior knots {0.25, 0.50, 0.75}, leading to seven B-spline basis functions for h -1 .

From the analysis with the proposed method, we get the estimated individual warping functions as in the left panel of Figure 9, and the warped (aligned) functions X i,j for each individual (right panel). The estimated covariance matrix ΣU (respectively Σθ ) is presented in the left (respectively right) panel of Figure 10. From this, we observe that there is more time variability induced by the coefficient of the second B-spline basis function in the decomposition of h -1 (., θ), whereas there is more amplitude variability caused by the coefficient associated to the third basis function in the decomposition of U , see also Figure 9 (right panel).

Conclusion and further discussion

this paper we considered a nonlinear mixed effect model for functional data. We apply a B-splines approximation on three different levels: on the inverse of the warping function describing the individual phase variability; on the global mean function and on the individual amplitude random effects. Random effects enter to model the individual amplitude as well as the phase variability. The main advantage of the proposed method is that it avoids the (costly) choice of landmarks, and that we can provide important theoretical support for the procedure: (i) convergence of the iterative algorithm to the target function(s); (ii) consistency and asymptotic normality of the estimators.

In this paper we considered the discrete T i time points to be fixed (non-random). However the methodology could be generalized fairly easily to random time points. Typically one would then need to assume that the distribution of the random time points is regular enough (meaning that there are no empty regions in the observed pattern of discretized time points). This would require, for example, an adaptation of the criteria used in Sections 3.2 and 3.3.2. An analysis of the assumptions, particularly modeling assumptions of the noise, with the aim to see the robustness of the method, would be of interest.We postpone this analysis to an experimental work.

A Proof of results in Section 2

A.1 Proof of Lemma 1

Let θ 1 and θ 2 such that t = w(w -1 (t; θ 1 ); θ 2 ). Then it follows that,

w -1 (t; θ 2 ) = t 0 exp(h -1 (u; θ 2 ))du 1 0 exp(h -1 (u; θ 2 ))du = t 0 exp(h -1 (u; θ 1 ))du 1 0 exp(h -1 (u; θ 1 ))du = w -1 (t; θ 1 ) ⇔ t 0 exp(h -1 (u; θ 2 ))du t 0 exp(h -1 (u; θ 1 ))du = 1 0 exp(h -1 (u; θ 2 ))du 1 0 exp(h -1 (u; θ 1 ))du def = δ(θ 2 , θ 1 ) ⇔ t 0 δ(θ 2 , θ 1 ) exp(h -1 (u; θ 1 )) -exp(h -1 (u; θ 2 )) du = 0.
As this equation is true whatever the value of t ∈ [0, 1], we have that the integrand is equal to 0 for every u ∈ [0, 1] except for possibly a countable number of points. As B-splines are continuous, it is equal to 0 for every u ∈

[0, 1]. It holds that, log(δ(θ 2 , θ 1 )) = h -1 (u; θ 2 ) -h -1 (u; θ 1 ) = K h l=-p h (θ 2 l -θ 1 l ) Bh l (u; κ h ).
As the left hand side does not depend on u, so the right hand side should be equal to 0 for all u. As (B h l ) l is a B-spline basis, it induces that θ 2 = θ 1 .

A.2 Proof of Theorem 1

First, suppose that

H U i = I Ti . Let Σ ε = σ 2 ε I Ti , and choose 0 < σ2 ε < σ 2 ε . Then, H U i (Σ ε -Σε ) = Σ ε -Σε . Since σ2 ε -σ 2 ε < 0, and (B U i ) ⊤ B U i ̸ = 0 m U i , ΣU = Σ U + (σ 2 ε -σ2 ε ){(B U i ) ⊤ B U i } -1 is semi positive, definite
, and is diagonal. Thus, we have found two parameters (Σ ε , Σ U ) and ( Σε , ΣU ) which define the same model. Hence we do not have identifiability.

Suppose now that model ( 8) is not identifiable. Then, according to Theorem 4.1 in [START_REF] Wang | Identifiability of linear mixed effects models[END_REF], for all (Σ ε , Σ U ), there exists ( Σε , ΣU )

̸ = (Σ ε , Σ U ) such that (B U i ) ⊤ Σ ε B U i ̸ = (B U i ) ⊤ Σε B U i ; H U i (Σ ε -Σε ) = Σ ε -Σε ; ΣU = Σ U + {(B U i ) ⊤ B U i } -1 (Σ ε -Σε )B U i {(B U i ) ⊤ B U i } -1 = Σ U + (σ 2 ε -σ2 ε ){(B U i ) ⊤ B U i } -1 . As Σ ε ̸ = Σε , H U i = I Ti . Since ΣU ̸ = Σ U , one gets (B U i ) ⊤ (Σ ε -Σε )B U i = (σ 2 ε -σ 2 ε )(B U i ) ⊤ B U i ̸ = 0 m U i , which implies (B U i ) ⊤ B U i ̸ = 0 m U i .
Moreover, as ΣU is diagonal, (B U i ) ⊤ B U i must be diagonal.

A.3 Proof of Theorem 2

We have

X i (t) = Xi w -1 w(t; θ i ); θi = µ(t) + U i (t) + ε i (t).
As X i (t) is Gaussian distributed, the identifiability of the warping process is equivalent to the identifiability of the mean and of the variance of X i . We thus investigate the expectation and the variance of X i (t).

Expectation:

E[X i (t)] = E Xi {w -1 {w(t; θ i ); θi }} = E θi, θi E[ Xi {w -1 {w(t; θ i ); θi }}|θ i , θi ] ; µ(t) = E θi, θi µ{w -1 {w(t; θ i ); θi }} ;
which is equivalent to, by projecting onto the B-spline basis, for all l = -p µ , . . . , K µ , and for all j = 1, . . . , T i , B µ l,pµ+1 (t i,j ; κ µ ) = E θi, θi B µ l,pµ+1 {w -1 {w(t i,j ; θ i ); θi }; κ µ } or, in a matrix representation, B µ i = E θi, θi (B µ i ) θi, θi .

Variance.

By the formula of the conditional variance,

Var Xi {w -1 {w(t; θ i ); θi }} = E θi, θi Var Xi {w -1 {w(t; θ i ); θi }} θ i , θi + Var θi, θi E Xi {w -1 {w(t; θ i ); θi }} θ i , θi . (13)
Now,

Var Xi {w -1 {w(t i,j ; θ i ); θi }} j=1,...,Ti = E θi, θi {(B U i ) θi, θi } ⊤ Σ Ui (B U i ) θi, θi + E θi, θi Var ε i {w -1 {w(t i,j ; θ i ); θi }} j=1,...,Ti + Var θi, θi (B µ i ) θi, θi α µ , and Var [{X i (t i,j )} j=1,...,Ti ] = (B U i ) ⊤ Σ Ui B U i + Var[{ε i (t i,j )} j=1,...,Ti ].
As Var[{ε i (t i,j )} j=1,...,Ti ] = σ 2 ε I Ti , and

Var ε i {w -1 {w(t i,j ; θ i ); θi }} j=1,...,Ti = σ 2 ε I Ti , (13) becomes
(B U i ) ⊤ Σ Ui B U i = Var θi, θi (B i ) θi, θi α µ + E θi, θi {(B U i ) θi, θi } ⊤ Σ Ui (B U i ) θi, θi .
B Proof of results in Section 3

B.1 Proof of Theorem 3

Before passing to the proof of Theorem 3 we study the weighted tail cross products, that are appearing in this context. This is done in Lemmas 2 and 3.

Lemma 2. The weighted tail cross product of µ(w -1 (t i,j ; .)) + U i (w -1 (t i,j ; .)) with itself (with respect to T i ) exists and Q θi ( θi ) =|µ(w -1 (t i,j ; θ i )) + U i (w -1 (t i,j ; θ i ))

-µ(w -1 (t i,j ; θi )) -U i (w -1 (t i,j ; θi ))| 2 has a unique minimum at θi = θ i .

Proof. We assume that µ and U i can be decomposed onto a spline basis. As every function is defined on [0, 1], µ • w -1 and U i • w -1 are L 2 . Then the tail cross product of µ(w -1 (t i,j ; θ i )) + U i (w -1 (t i,j ; θ i ))

with µ(w -1 (t i,j ; θi )) + U i (w -1 (t i,j ; θi )) exists because it converges uniformly for θ i ∈ R r and θi ∈ R r .

The function Q has a unique minimum at θ i : if Q( θi ) = 0, for every t ∈ [0, 1],

(µ + U i )(w -1 (t; θi )) = (µ + U i )(w -1 (t; θ i )) ⇒ t = w(w -1 (t; θ i ); θi ).

Then, as the warping function is injective (see Lemma 1), θi = θ i .

Lemma 3. For l = 1, . . . , r and l ′ = 1, . . . , r, the derivatives ∂ µ(w -1 (t i,j ; θ i )) + U i (w -1 (t i,j ; .)) ∂[θ i ] l and ∂ 2 µ(w -1 (t i,j ; θ i )) + U i (w -1 (t i,j ; .))

∂[θ i ] l ∂[θ i ] l ′
exist and are continuous on R r and all weighted tail cross products in between µ(w -1 (t i,j ; θ i )) + U i (w -1 (t i,j ; .))

and its first and second derivatives exist.

Proof. Each function is defined on [0, 1], and is either decomposed onto a spline basis, or composed with the exponential function, then there is no problem to intervert derivation and integrals. As explained before, the weighted tail cross products exist because those functions are L 2 ([0, 1]).

Equipped with Lemmas 2 and 3 we can adapt results from [START_REF] Jennrich | Asymptotic properties of non-linear least squares estimators[END_REF] to prove Theorem 3.

Proof of Theorem 3. We consider the following compact set:

E α (θ 0 , Σ θ ) = x ∈ R r | (θ 0 -x) t (Σ θ ) -1 (θ 0 -x) ≤ χ 2 r (1 -α)
where χ 2 r (1 -α) denotes the 1 -α quantile of the χ 2 -distribution with r degrees of freedom. With probability 1 -α, θ i and θTi i belongs to E α (θ 0 , Σ θ ). (Jennrich, 1969, Theorem 6) is used to get the strong consistency and (Jennrich, 1969, Theorem 7) is used to get the asymptotic normality of θTi i .

B.2 Proof of Theorem 4

Proof. We use the dominated convergence theorem. As the variance is lower bounded by Σ θ , the mean is upper bounded by a fixed constant, so we can construct a dominating function. Then, φ(θ i ; c, C)φ(η; θ 0 , σ 2 ε a -1 i (θ i ) + Σ θ ) dθ i = R r φ(θ i ; θ 0 , Σ θ )φ(η; θ 0 , Σ θ )dθ i = φ(η; θ 0 , Σ θ ).

B.3 Proof of Theorem 5

We consider η ∈ E ρ (θ 0 , A), with A positive definite.

m θ∞ i (η) = R r φ θ i ; η, σ 2 ε a -1 i (θ i ) φ θ i ; θ 0 , Σ θ dθ i = R r φ θ i ; η, σ 2 ε a -1 i (θ 0 ) φ θ i ; θ 0 , Σ θ × det(a -1 i (θ 0 )a i (θ i )) × exp - 1 2σ 2 ε (θ i -η) T a -1 i (θ i ) -a -1 i (θ 0 ) (θ i -η) dθ i = φ(η; θ 0 , σ 2 ε a -1 i (θ 0 ) + Σ θ ) R r φ(θ i ; c 0 , C 0 ) × det(a -1 i (θ 0 )a i (θ i )) × exp - 1 2σ 2 ε (θ i -η) T a -1 i (θ i ) -a -1 i (θ 0 ) (θ i -η) dθ i
where the last equality comes from Lemma 4, with C 0 = (σ -2 ε a i (θ 0 ) + (Σ θ ) -1 ) -1 ; c 0 = C 0 σ -2 ε a i (θ 0 )η + (Σ θ ) -1 θ 0 = θ 0 + C 0 σ -2 ε a i (θ 0 )(η -θ 0 ).

Our goal is now to prove that

I = R r φ(θ i ; c 0 , C 0 ) det(a -1 i (θ 0 )a i (θ i )) × exp - 1 2σ 2 ε (θ i -η) T a -1 i (θ i ) -a -1 i (θ 0 ) (θ i -η) dθ i is close to 1.
We then divide this integral into two parts: let ρ > 0, if θ i ∈ E ρ(θ 0 , A), a -1 i (θ i ) = a -1 i (θ 0 ) + D θ0 a -1 i (θ i -θ 0 ) + O(ρ), ⇒ (θ i -η) T (a -1 i (θ i ) -a -1 i (θ 0 ))(θ i -η) = O(ρ + ρ), and det(a -1 i (θ 0 )a i (θ i )) = 1 + O(ρ). If θ i ∈ E ρ(θ 0 , A) c , by Assumption E we get that det(a -1 i (θ 0 )a i (θ i )) = O(1)

⇒ exp - 1 2σ 2 ε (θ i -η) T a -1 i (θ i ) -a -1 i (θ 0 ) (θ i -η) ≤ 1.
This leads to

I = (1 + O(ρ))
Eρ(θ0,A) φ(θ i ; c 0 , C 0 )dθ i + O(1)

Eρ(θ0,A) c φ(θ i ; c 0 , C 0 )dθ i = 1 + O(ρ).

B.4 Proof of Theorem 6

From Theorem 4,

T i ( θTi i -θ 0 ) d -→
Ti→+∞ N r (0, Σ θ ).

For i = 1, . . . , n, let Z Ti i ∼ N (θ 0 , 1 Ti Σ θ ), and ∆ Ti i = θTi i -Z Ti i . We know that ∆ Ti For i = 1, . . . , n, let Z Ti i ∼ N (θ 0 , 1 Ti Σ θ ), and ∆Ti i = √ T i ( θTi i θTi i -Z Ti i ). We know that √ T i Z Ti i ∼ N (θ 0 , Σ θ ) and then n i=1 T i Z Ti i (Z Ti i ) T ∼ W(Σ θ , n -1). Moreover, we know that ∆Ti Thus, we get

1 n -1 n i=1 T i θTi i ( θTi i ) T = 1 n -1 n i=1 T i Z Ti i (Z Ti i ) T + 1 n n i=1 T i ∆ Ti i (∆ Ti i ) T d -→
n→+∞,min Ti→+∞ W(Σ θ , n -1).

Figure 1 :

 1 Figure 1: Pinch Force dataset.

Figure 2 :Figure 3 :

 23 Figure 2: A sample of size 100 of warping functions (left), and the empirical covariances of these functions (right).

  ) with mean function μk (k = 1, 2, 3, 4) as model M k , and with mean function µ 4 as model M 4 .We consider a low-dimensional setting in which n = 100 and T i = 70, as well as a high-dimensional setting in which n = 200 and T i = 150. We use M

Figure 4 :

 4 Figure 4: A sample of the warped process X (left) and the un-warped process Y (right) for M 1 .

Figure 5 :

 5 Figure 5: SNR functions for M 1 (left) and M 2 (right).

Figure 6 :

 6 Figure 6: Comparison between μ1 and μ1 for model M 1 (n = 100 and T = 70, top) and for model M HD 1

Figure 7 :

 7 Figure 7: Simulation results for the components of ( Σθ i,i ) 1≤i≤7 for the proposed procedure for model M 1 (left) and M HD 1 (right). Boxplot without the outliers. True coefficient values: red horizontal lines.

Figure 8 :

 8 Figure 8: Simulation results for the components of ( ΣU i,i ) 1≤i≤10 for models M 1 (left) and M HD 1(right) for the proposed procedure. The horizontal red line present the true value.

Figure 9 :Figure 10 :

 910 Figure 9: Estimated warping functions (left) and warped individual curves (right).

  i ; c, C)φ(η; θ 0 , σ 2 ε a -1 i (θ i ) + Σ θ )dθ i

1

  Ti i , we use Lindeberg Central limit Theorem, recalled in Theorem A of Appendix C. Let S 2 n = Σ θ n i=1 Ti . We assume that (14) holds. Then,

Table 2 :

 2 Simulation results for the proposed procedure for models M 1 and M

	HD 1 .

Table 3 :

 3 Simulation results for Σ θ using the proposed procedure for models M 1 and M

	HD 1 .

Table 4 :

 4 Simulation results for the proposed procedure for models M

	wCDG 2	wGC , M 2

HD 3 , and Models M 4 , M 4 and M HD 3 . Table

Table 5 :

 5 Simulation results for the proposed procedure for models M 3 and M

	HD
	3

Table 6 :

 6 Simulation results for the proposed procedure for models M 4 , M 4 and M We compare the four methods introduced in Section 4.2 with the proposed one on models M 1 , M

	HD 4 .

Table 7 :

 7 Simulation results for µ(•) for the proposed and competitive methods. Method (3)=Bayesian warping; Method (4)= elastic square-root slope.

		Simulation results for models M 1 , M HD 1 , M wCDG 2	wGC and M 2
	Method	M 1 bias variance	M HD 1 bias variance	M bias variance wCDG 2	wGC M 2 bias variance
	proposed 0.0102 0.0090 0.0055 0.0065 0.0230 0.0236 0.0053 0.0052
	GC	0.0454 0.0043 0.0588 0.0048 0.0703 0.0134 0.0280 0.0014
	(3)	0.0435 0.0009 0.0201 0.0001 0.0493 0.0003 0.0423 0.0001
	(4)	0.1062 0.0011 0.0745 0.0009 0.1277 0.0014 0.1245 0.0011
	fdakma	0.0075 7.10 -6 0.0069 3.10 -6 0.0281 2.10 -5 0.0283 1.10 -5
		Simulation results for models M 3 , M	HD 3 , M 4 and M HD 4
	Method	M 3 bias variance	M HD 3 bias variance	M 4 bias variance	HD M 4 bias variance
	proposed 0.0183 0.0153 0.0126 0.0105 0.0088 0.0044 0.0057 0.0018
	(3)	0.0794 0.0024 0.0491 0.0002 0.0280 0.0003 0.0212 0.0001
	(4)	0.1851 0.0009 0.1951 0.0007 0.1555 0.0016 0.1754 0.0026
	fdakma	0.0277 4.10 -5 0.0272 2.10 -5 0.0113 5.10 -6 0.0104 2.10 -6

The first and third author gratefully acknowledge support from the C1-project C16/20/002 of the KU Leuven Research Fund. The resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation -Flanders (FWO) and the Flemish Government. Part of this work was accomplished when the second author was a Postdoctoral Researcher at the KU Leuven.

Acknowledgements

The authors thank Professor Alois Kneip for helpful discussions, and Professor Daniel Gervini for help with his computer codes. E. Devijver also sincerely thanks R. Molinier for fruitful discussions.

B.5 Proof of Theorem 7

We check that our model satisfies the assumptions given in (Pinheiro, 1994, Chapter 3):

The U i are independent and follow a N (0, Σ Ui ) distribution, ε i follows a N Ti (0, σ 2 ε ) distribution and the U i are independent of ε i

The matrix B µ i is of full rank, as it is a functional basis n ≥ m µ + 1 + 1

The concatenated matrix [B µ i , B U i ] has rank greater than m µ if we don't take the same basis for µ and U i

The matrices I Ti and

So we get the asymptotic normality of the estimator.

B.6 Proof of Theorem 8

Recall the different operator parts of the iterative algorithm in (12). The Banach fixed point theorem, recalled in Theorem B of Appendix C, is used to prove that there is a unique fixed point, and that the algorithm converges. To use this theorem, we work in R mµ+nm U +1+nr with the Euclidean distance. It is a non-empty complete metric space. The mapping we consider is Ψ, as defined in (12). We need to prove that Ψ is a contraction mapping.

Denote by k f the Lipschitz constant for the function f . We want to find k Ψ such that, for (x, y)

As Ψ 2 and Ψ 3 are linear, for (x, y)

The proof relies on Lemma 5 applied to Ψ 3 and Ψ 2 , defined via the mixed effect models. The statement of Lemma 5 and its proof can be found in Appendix C. Then, using Assumption G, there exists k Ψ < 1 such that Ψ is k Ψ -Lipschitz. Banach fixed point theorem concludes.

B.7 Proof of Theorem 9

First, we prove that for a fixed iteration ite, the several computations we are doing to the true parameters (α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ ) are keeping it fixed under Assumption C.

Fix the iteration number (ite) and consider (α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ ). We know (θ 0 , Σ θ ) and we predict θ i with the BLUP. As n → ∞ and min i T i → ∞, predictions are good. Then, we estimate (

by Theorem 7 we get (α µ , σ 2 ε , Σ ∼ U ) (strong consistency). Then, we approximate θ i for all i = 1, . . . , n. By Theorem 4, those pseudo-observations are close to the true random variables, with the good distribution function. Finally, we estimate a linear mixed model on those observations: by Theorem 6, we get θ 0 and Σ θ (consistency).

Then, under identifiability, Assumptions A, C and G, (α µ , σ 2 ε , Σ ∼ U , θ 0 , Σ θ ) is a fixed point of Ψ. By Theorem 8, there exists only one fixed point: then

This convergence is almost surely, as the convergence in each step is almost surely.

C Additional useful results and tools.

Theorem A (Lindeberg Central Limit Theorem). Let (Ω, F, P) be a probability space, and X k : Ω → R, k ∈ N, be independent random variables defined on that space. Assume that the expected values

If this sequence of independent random variables X k satisfies Lindeberg's condition: for all ε > 0,

where 1 is the indicator function, then the central limit theorem holds, i.e. the random variables

converge in distribution to a standard normal random variable as n → ∞.

Theorem B (Banach fixed point theorem). Let (X, d) be a non-empty complete metric space with a contraction mapping T : X → X with Lipschitz constant q ∈ [0, 1). Then, T admits a unique fixed-point x ∞ in X. Furthermore, x ∞ can be found as follows: starts with an arbitrary element x 0 in X and define a sequence {x n } by

Lemma 4 [START_REF] Ahrendt | The multivariate gaussian probability distribution[END_REF]). The following holds, for all x ∈ R p , for A, B positive definite matrices of size p × p and (a, b) ∈ (R p ) 2 :

Lemma 5. Let π 1 and π 2 be two orthogonal projections, and denote |||.||| the operator norm. Then,

Proof. We first prove (15). Let π be an orthogonal projection. As it is a projection, its norm is larger than 1. As it is an orthogonal projection, we can use Pythagorean theorem to prove that its norm is smaller than 1.

We now prove (16). If E π1 ∩ E π2 ̸ = {0}, let x ∈ E π1 ∩ E π2 . Then, π 1 • π 2 (x) = x so |||π 1 • π 2 ||| ≥ 1. If E π1 ∩ E π2 = {0}, assume that |||π 1 • π 2 ||| = 1: there exists x ̸ = 0 such that ∥π 1 • π 2 (x)∥ = ∥x∥. But as π 1 and π 2 are projections, it means that x ∈ E π1 ∩ E π2 : contradiction.


Geo-photographers

Danièle MÉAUX

In the field of contemporary photography, as in other sectors of cultural life), a real geographical tropism is undeniable. Many use the expression of "spatial turn" 1 when referring to this rise of interest among artists and researchers from varying disciplines in the different types of spatiality in which man evolves. The exhibition, New Topographics: Photographs of a Man-Altered Landscape, presented in 1975 at the George-Eastman House in Rochester, New York, partly initiated a new vision of territories in the field of photography, where the imprint of man was now apparent 2 . The DATAR 3 Photographic Mission (1983-1989) confirmed the trend between 1983 and 1989 and had an important impact on the way landscapes were to be viewed, no longer presenting them as external objects of contemplation but as spaces of human life. Geographical tropism is now a developing theme in many photographic practices, including in particular those where photographers use protocols, more or less directly inspired by those used by geographers. Those protocols are often more revealing than the actual photographic subjects. Photography is an art of space, but not all landscapes relate to a geographical approach. Nevertheless, using the protocols confers upon the works a form of scientific empiricism inspired by the field practices of certain geographers.

'Repeat' Photography

The DATAR Photographic Mission was a commission that brought together town and country planners and photographers in pursuit of a triple aim: to observe the evolution of French landscapes, to re-create a culture of landscape and finally to promote photography in the artistic domain. 4 The photographers, selected on the basis of their work and projects, were given freedom in pursuing the aims of the project. _________________________________ 1 See Edward W. Soja, Postmodern Geographies. The Reassertion of Space in Critical Social Theory (London: Verso Book, 1989); Lévy J., Le Tournant géographique. Penser l'espace pour lire le monde (Paris : Belin, 1999). 2 William Jenkins, New Topographics. Photographs of a Man-altered landscape (Rochester: International Museum of Photography at George Eastman House, 1975). This exhibition catalogue was reproduced in New Topographics (Göttingen: Center for Creative Photography / University of Arizona / George-Eastman House / Steidl, 2010). In his introduction to the catalogue, William Jenkins insisted on the importance of protocols which allowed a topographical approach.

3 DATAR : Délégation à l'Aménagement du Territoire et à l'Action Régionale. (Delegation for Territorial Planning and Regional Action) 4 See on the subject the analyses of Jordi Ballesta, in Le Projet photographique comme experience et document photographiques, Doctoral Thesis defended on 10th December, 2011, under the direction of Yannis Tsiomis, Director of Studies at the EHESS, 37-87.

Created in 1991 on the initiative of the Ministry of the Environment, its goal was to "constitute a collection of photographic seriesshots repeated in an identical way at regular intervalsto allow an analysis of the mechanisms of spatial change and the actors involved, in order to orientate the future development of landscape areas". 5 The protocol imposed was one of "repeat" photography. This constraining feature was to guarantee that the results could be exploited. All the photographs taken are precisely dated. The Landscape Photographic Observatory set up a form of monitoring activity around a certain number of sites. 6 The aim was to observe the transformations of sites and premises, generated by the interaction of multiple parameters, of which certain cannot be controlled. There is the work of the town and country planners but also the individual initiatives of inhabitants, the effects of climate, the growth or retreat of vegetation … For the Observatory, landscape is viewed as a complex and dynamic reality, with changes subject to different rhythms. The changes observed concern for example the colours of facades, the development of industrial or commercial zones, the growth of urban sprawl, phenomena of abandon of agricultural land, deforestation or reforestation …. This type of observation and monitoring is now considered essential: the period of modern architecture has moved on and it is now accepted that projects can only be conceived by taking into account the complexity of the present.

The photographs, in the manner in which they are repeated, take on the function of providing information and knowledge. They are designed as documents which will subsequently be subject to analysis. The adoption of a strict protocol places the operator in the role of an observer of the obtained results. It restricts his freedom and limits subjectivity. A certain neutrality is required to enable comparisons to be made. The photographers who are solicited to carry forward the project are nevertheless renowned as authors or artists.

Itineraries and Transects

In all, nineteen "itineraries" 7 have been officially designated by the Landscape Photographic Observatory. The use of this term emphasises the importance of the field experience. The sites must be surveyed; the visits carried out with the technical help of maps and maybe even a compass. In Pilat Itineraries, Sophie Ristelhueber accompanies each image with a section of the 1:50,000 scale map indicating her position and orientation.

__________________________________________________________

Each "observation point" is numbered. A text specifies each time the characteristics of the site and the theme adopted. 8 All these elements contrive to illustrate the geographical nature of the approach.

Thierry Girard is a photographer who constantly associates itinerary with his photographs (whether it relates to a commission or not). He notes: "My work has now for some time consisted of the notion of a restricted 9 circuit […]". Many of his outings are in fact regulated by specific protocols. For "From one sea to another" 10 a straight line drawn on the map determined his route from Nice to Ouessant. The colour photograph on the book cover shows sand and water trails; the image is cut by a prominent diagonal slash in the form of winding and irregular roads on a map with the red colour usually used to depict national arterial roads.

In a similar fashion, Dennis Adams and Laurent Malone followed a rectilinear itinerary from Manhattan to the John F. Kennedy Airport; each time one of the two associates took a photograph, the other did the same in a diametrically opposite direction, without taking the time to compose the image, to regulate the aperture or to adjust the focus. The 243 pairs of photographs they produced were brought together in a book 11the two shots taken from the same view-point were placed opposite each other, one on each side of the page, thus more or less materialising the trajectory they followed.

Adopting a pre-arranged itinerary differentiates the journey form a simple promenade or stroll, resembling more an exploratory circuit, a controlled progression interspersed with the taking of test samples. The protocol restricting the freedom of the practitioner appears as the condition of a methodical collecting process, allowing a form of observation of places and their appropriation by man.

This method is similar to the transect process used by geographers; a virtual line is drawn in order to study the incidence of such and such a phenomenon. Along this axis, samples are taken or observations maderelating perhaps to the vegetation, the geological composition of the soilwhich enable the region to be mapped. The method could also be applied to urban zones in order to study different types of housing, the land values or the population density. The geographer's transects however are not drawn up in a haphazard form. In forest zones they are perpendicular to contour lines. To be representative 12 they need to cross all the geological formations. To establish a transect requires preliminary study before going into the field. When a transect is drawn up for an urban zone its route is also deliberately delineated. According to Bernard Jouret, who carried out a study of the Brussels conglomeration, the transect should be drawn radially to allow pertinent observation. 13 Transects of photographers are not always based on a scientific analysis of data; they have an arbitrary and playful dimension. The diagonal line followed by Thierry Girard contrasts with the usual route, more often than not based on a north/south axis. His trajectory allows him to pass through less populated areas, poor and often disinherited regions. Following a restricted route also leads to a focusing on ordinary sites which are not usually taken into consideration.

Typologies

Elements, be they natural or man-made, can be combined on the basis of determined criteria to establish a category. Deciding on classifications of objects based on defined parameters appears as the means to attain a certain form of generality or abstraction, both in geography and observation sciences in general. Establishing a typology refers back to the world of science (even if rather bizarre categories exist such as those defined by Borges). Discerning a "category" is « "[…] the methodical process, the path (odos), the simplest way to allow the recognition and appreciation of an object obscured by an indistinct crowd of qualitative differences or particularities" 14 noted Pierre-Henry Frangne. He continued: "classifying allows one to leave behind variety (necessarily cumulative and inchoative) to take on variation which places the difference within the identity and the identity within the difference". 15 Hence the constitution of a series of images, corresponding to a given category, legitimises comparison. It enables the observation of discernable variations and potentially incites a questioning of their causes. The visible leads to one wondering about the invisible.

From the 1960's Bernd and Hilla Becher promoted the conception of typologies, driven by a powerful attraction for categories and systems which tended to place them with the conceptual artists. They were also drawn by the aptitude of series to contribute to an analysis of reality. _________________________________ 12 « A transect is a rectilinear itinerary of prospection and/or sampling covering a maximum diversity of topographical, geomorphological and vegetal situations » See René Delpech et alii, Vocabulaire : Typologie des stations forestières (Paris: Institute of Forestry Development, 1986). 13 Bernard Jouret, « The transect method applied to urban analysis. Example of Brussels », Revue de géographie de Lyon, volume 47, n° 1 (1972): 77-96. 14 Pierre-Henry Frangne, « Inventaires et inventions », in Classifications (Rennes: Rennes University Press, 1996), 29. 15 Ibid., p. 32.

The observable differences between the listed industrial buildings were highlighted in pictures featuring several views; it's up to the spectator to wonder if they were used for markedly different purposes, as properties in natural settings or for aesthetic considerations. At the Düsseldorf Kunstakademie (Art Academy) the couple's teaching was based on a few simple principles: they asked their students to choose an architectural or social subject and to adopt a uniform and neutral style so as to produce a large quantity of images to extract a typological perspective. Erik Verhagen noted: "All the students went through this initiatory stage … making it easier for them to subsequently rid themselves of the constraint". 16 While Becher's students were later able to assert their own personal styles, there is a strong tendency amongst current photographers to produce typological work (to the point perhaps of establishing a trend).

These productions often involve landscapes and geography since they highlight the way in which man appropriates and organizes the space in which he lives. They emphasise the fact that the planning and lay out of sites is a collective task of societies since human activity leaves its mark an imprint on the ground and transforms it. Jean-Marc Besse wrote: "The landscape is an expression of human effort, always fragile and to be recommenced, to inhabit the world". 17 Laurent Guéneau recently produced a series entitled "Dominante verte" 20 (Dominant green), consisting of photographs from workers' allotments in the Saint-Etienne region 21 . Through the colours and the subtle lighting the spectator tends to focus on the varied lay-outs of the plants, the structures and tones of the garden sheds, the rudimentary systems set up to collect the rainwater, the common functioning, but with all its diversity, and the integration of the gardens in the environment.

"Living in a place […], transforming a site into an inhabited environment, is generally […] a question of forming habits, of living there in a systematic way […]. What Jackson calls the vernacular landscape corresponds in fact to this manner of gradually inhabiting the world by establishing habits" 22 . The series of Laurent Guéneau seeks to highlight this imprint of lifestyles on the site.

Surveys

A survey based on the use of a questionnaire can be used in geography 23 . During the presentation of the results the record of people's answers is often associated with a photographic representation of the sites concerned. In this way sites are addressed through a relational and cultural dimension. The aim is to know how they are appropriated by the inhabitants. It is more or less this relation between man and his environment that Augustin Berque refers to as "mediance" 24 .

Sophie Calle, a proponent of "relational art" 25 , according to the definition of Nicolas Bourriaud, takes up this survey protocol to produce "Memories of East Berlin" 26 ; a title used both for the installation presented at the Berlin gallery of Matthias Arndt, at the Strasbourg Museum of Modern and Contemporary Art in 1999, and then for a book published by Actes Sud. Following reunification, the Berlin Senate created a commission to decide on the approach to be adopted towards the structures symbolizing Soviet power, located in the East Berlin sector. Some were covered up, others dismantled. For ten of these constructions Calle used the following procedure: she photographed the site where the monument had once stood and questioned the inhabitants on their feelings towards the site and their recollections. In the installation, as in the book, she juxtaposes the coloured photograph of the now empty site, the listed statements and an old image of the site, smaller in black and white, showing the monument which had been removed. _______________________ 20 « Dominante verte » was exposed at the Arles International Photographic festival in 2010, at the Montmajour Abbey, then at the French National Library from September to November, 2010. 21 Saint-Étienne was one of the first cities to introduce workers allotment gardens on a large scale in 1895. 22 Jean-Marc Besse, Le Goût du monde. Exercices de paysage, op. cit., p. 131 23 Lydie Goeldmer-Gianella & Anne-Lise Humain-Lamoure, « Survey by questionnaire in geography of the environment », L'Espace géographique n° 4, December 2010. 24 Augustin Berque, Écoumène. Introduction à l'étude des milieux humains (Paris: Belin, 2000), 128. 25 Nicolas Bourriaud, Esthétique relationnelle (Dijon: Les Presses du réel, 1998). 26 Sophie Calle, Souvenirs de Berlin-Est (Arles: Actes Sud, 1999).

The collected phrases reveal an assortment of feelings, from satisfaction to regret, from indifference to an impression of humiliation.

But what emerges in general is the powerful relationship which unites the inhabitants with statues and symbols that accompanied their daily lives for so long. On the subject of the bust of Lenin which stood in front of the Russian Embassy and which is now concealed in a type of box, we can read: "I don't miss it, but I would prefer it still to be there. Even now when I pass in front of the façade I keep looking for it" 27 . On the gigantic statue of Lenin which stood in the United Nations Square and is now gone: "I think it had its place here, it was part of the architectural design of the district. It was a prominent landmark, occupying the front of the stage. Everything around had to measure up to it.

The traffic system in place meant you had to drive around it, whatever your direction. I no longer know if Lenin was put up before the surrounding buildings, or vice versa, but the square with its houses in an S shape, almost in the form of waves, created an ensemble with the monument. Today it no longer works" 28 .

A survey allows sites to be considered in terms of the meaning inhabitants attribute to them, the living memories they have. Memories differ from history 29 ; in the sense that they invest the site with a relational, even an identity dimension, it interests geographers today who are committed to studying symbolic spatiality 30 .

A Geo-photograph ?

The way in which certain photographers now contemplate space (in its complexity, its plurality and dynamics, in its relational and cultural dimension), takes inspirationin a direct or indirect mannerfrom geographical research, recent reflections about the landscape, the place and the region, with an increasingly developed awareness of the environment. A certain number of contemporary photographers, adopting protocols inspired by geography, now use large format view cameras. This cumbersome apparatus corresponds both to the highly concerted character of their approach and to a willingness to record the appearance of a site in the most exhaustive way possible. It is however more especially the adoption of protocols (repeat shooting, itineraries and transects, typologies and surveys) which enable practitioners to come close to a relation with space resembling that of geographers. Through the intermediary of these protocols a whole epistemological dimension is called upon. ________________________________________ 27 Ibid., p. 14-15. 28 Ibid., p. 20. 29 Pierre Nora, Lieux de mémoire 1 (Paris: Gallimard, 1997), 24. Nora distinguishes between history which is « always a problematic and incomplete reconstruction of what no longer exists » and memory which is « life still led by living groups and consequently […] in constant evolution ». 30 Michel Bédard, « Une typologie du haut-lieu, ou la quadrature d'un géosymbole », Cahiers de géographie du Quebec, volume 46, n° 127 (2002): 69.

It is not necessary to dwell on the objectivity with which the medium is normally credited. Every image stems from a procedure which is partly automated, in such a way that Gary Winogrand was able to write "I take photographs to find out what something looks like when it is photographed" 31 , considering the shot taking to be an operation controlled by a device and making possible a particular form of observation of the world. Scientists from the 19 th century immediately understood the service the medium could render to science. This related to its precision but especially to its automaticity. To the extent that it is not entirely filtered by the choice of the operator, the photograph can retain traces which are not considered pertinent at the time by the latter, but are susceptible to become instructive thereafter. They constitute, in short by provision, an archive of reality which can be exploited at a later stage. 32 Combining a defined protocol with the photographic procedure reinforces and extends the inherent empirical dimension of the medium: the photographer seems to be setting in motion an experiment with strict parameters, and thereafter to view the result. Once the constraint is accepted, any deviation is prohibited. If the photographer assumes more or less the role of initiator of the protocol, he then loses the initiative and in the end seems only able to enjoy the fruits of the sampling conditioned by his "programme". The fact that the adopted protocols are more or less inspired by geographical practices epitomises a marked attraction for the discipline.

The images are not conceived in isolation, they are taken in correspondence with other images or with a text. Above all they are taken in terms of a global approach, hence their worth. More than the photograph itself, it is obviously the projectincluding the shooting and the restrictive approachwhich demonstrate the geographical model. The protocol places the photograph (already subject itself to empiricism) in a context of the uncertainties of research. The project as a whole appears in the light of the dynamic aspects of a previously designed experiment carried out in the field, a venture aimed at the acquisition of observations and knowledge. _______________________ 31 Gary Winogrand G., « Interview », in Peninah R. Petruck, The Camera Viewed. Writings on Twentieth Century Photography, volume 2 (Syracuse: Dutton, 1979), 126. 32 "When a zoologist makes a drawing he only represents what he has noticed of his model, and consequently, the line traced by his pencil just interprets the more or less complete idea he has formed of the object to be reproduced, and it is very rare that the detail thus obtained clearly shows characteristics that the author would not have taken into account. Hence, when in the progress of science, one of his successors introduces, in his quest for zoological solutions, elements which the first drawer would not have used, it is again rare that he finds them faithfully represented in the other's drawings. In order to identify the presence or absence of these structural characteristics, he cannot make do with the examination of figures which have already been published, but is obliged to observe objects in nature once more. But with photography he can act differently, since a photographic image represents not only what the author himself has seen and seeks to show, but also what is actually visible in the reproduced object. Another naturalist can therefore focus on aspects the other has not noticed, and make real discoveries thanks to the image, as he would have done when observing the objects in nature." In « Rapport sur un ouvrage inédit intitulé: Photographie zoologique, par MM. Rousseau et Devéria », Comptes rendus hebdomadaires des scéances de l'Académie des sciences, t. XXXXVI, 1853, 991-993, in La Photographie en France. Textes et Controverses : une anthologie 1816-1871 (Paris: Macula, 1989), 78.

Artialisation / articities

The inclusion of photographs in a project, governed by a protocol, tends to resemble sample taking, listings undertaken according to defined criteria and meant for further analysis. They would appear therefore to oppose the notion of "artialisation" 33 which some consider to be the very prerequisite for "landscape". In reverse of a pictorial tradition, an image is not valued for its aesthetic potential recognized in its own right; it does not claim to be the equivalent of a "state of mind" 34 . Making use of a protocol is a way for the practitioner to place himself in another field, in fact more precisely to bestow much more importance on the "field" itself. The transcription of reality is highlightedwithout being validated by any objectivity of the views. Does this imply however that these images do not belong to the field of art ? In this regard the Observatory of Landscape Photography called on professionals recognised as authors or artists, and not technicians. While this choice is likely to serve the promotional needs of regions, it is also based on a recognition of the capacity of artist to choose an instructive and embracing focus which render visible the relations existing between things.

Certain planners today have become aware of what they can learn from artists. As for these artists, they have a tendency to assume certain restrictive procedures; accepting commissions which may serve to renew their approach, moving more towards decentering, far removed from usual practices or simplistic effects. 35 This is roughly what Thierry Girard says, talking about his responsibility for an itinerary in the Vosges Natural Park:

[…] the constraints and methodological rigour, inherent in the notion of « observatory » have a tendency to restrict the desire for aesthetic experimentation […], which in the end obliges the authors to remain very modest and to accept the restrictions of a particular commission, rather than subjecting the commission to their artistic ego. […] I have to admit that to begin with, I found the constraints a little too strict and restrictive, but I soon accepted the rules to a point of even finding a form of revitalisation. 36 In fact photographers who work on typologies, surveys or transects, are often doing this on their own initiative. While adopting a protocol distances the images from any artialisation, by following a more or less scientific model, it in no way diminishes the artistic dimensions of the production. Quite the opposite: the operator comes to resemble a designer, moving away from the notion of testimony (to which photography has for a long time been confined).

______________________________ 33 Alain Roger, Court traité du paysage (Paris: Gallimard, 1997). 34 Henri Frédéric Amiel, Fragments of a personal diary [1880] (Paris: Stock, 1931), 36. 35 See, with regard to the subject, the statements of Gabriele Basilico who notes how much the DATAR photographic commission modified his way of working, imposing more slowness and patience. 36 Thierry Girard, Un Voyage en Saintonge (Saintes: CCR de l'Abbaye aux Dames, 1995).

Referring to an exhibition of photographs during the International Photographic Festival of Arles in 2012, Pascal Beausse remarked on a tendency to redefine "artistic activity as a laboratory of knowledge" 37 . Today photographs reveal a capacity to propose complex and original visibilities, critical representations which even interest the human sciences and in particular geography. Results which can be exploited are of course not always the case. Certain projects carried out in this field serve more to promote the photographer or the accreditation of regions than in-depth studies of spatial phenomena. Nevertheless the protocols planned here are the emblem of the capacity of photography to take part in geographical ventures.

These protocols have also an intrinsic importance: the spectator is appreciative of the simplicity, the scientific nature of the approach which convey a certain recoil from subjectivity. The protocol combined with the shooting constitute, in a manner of speaking, an "attitude", a "form" 39 , where the patience and regularity, the empiricism and humility are capable of arousing emotion. It contributes to what could be called a "geographical aestheticism". _____________________ 37 Pascal Beausse, « Documents pour une information alternative », in Les Rencontres de la Photographie d'Arles: une école française (Arles: Actes Sud, 2012), 274. 39 « When attitudes become form: live in your head » is the title of an exhibition organised by Harald Szeeman in1969 at the Berne Kunsthalle.

______________________________8

  Sophie Ristelhueber, Pilat. Itinéraires (Pilat Natural Park / Ministry of the Environment, 1994). The following themes were planned: « development of wastelands and urbanization » (p. 34), « development of afforestation and sites (restoration and promotion)» (p. 64), « development of infrastructure equipment» (p. 48), « development of agriculturewastelandafforestation » (p. 44). 9 Thierry, Girard, « La route de Tôkaidô » in Une route, un fleuve. Itinéraires récents au Japon et en Europe Centrale, catalogue of the exhibition organised at the Marine Gallery, Nice from the 5th of December 1998 to the 24th of January 1999, Published by the City of Nice, 1999, work unpaginated. 10 Thierry Girard, D'une mer l'autre (Paris: Marval, 2002). 11 Dennis Adams, Laurent Malone, JFK (Paris: LMX, 2002).

  Further on, when referring to the work of John B. Jackson, he pointed out that "Landscape in reality is a reflection of nature and society, an integration of natural elements and human projects, a synthetic reality […]".18 Many series produced over the last few years have questioned this hybrid reality. They don't always possess the austere rigour of the typologies elaborated by the Becher couple, but the inventory of a category (of sites or lay-outs) does however allow an analysis of how social activity has come to organise space.NaoyamaHatakayama recently published a work entitled Terrils (Slag heaps) 19 . Page after page is full of these artificial mountains, usually reflecting extraction activity which has ceased, and which now in differing ways fit in to the surroundings. The varied shapes are undoubtedly the result of different operational methods, and various solutions have been found to transform the sites once the extraction over. The series highlights certain visibilities, and the causes deserve being examined. ____________________________ 16 Erik Verhagen, « Une objectivité en trompe-l'oeil. Bernd et Hilla Becher et leurs élèves », in Objectivités (Paris: Museum of Modern Art of the City of Paris / ARC, 2006), 63. 17 Jean-Marc Besse, Le Goût du monde. Exercices de paysage (Arles/Paris : Actes Sud / ENSP, 2009), 41. 18 Ibid., p. 44. 19 Naoya ma Hatakayama, Terrils (Light Motiv, 2011).

« Presentation of the Landscape Photographic Observatory », Séquences / Paysages n° 2 (2000): 5.

Sophie Ristelhueber works in the Pilat Natural Park, Raymond Depardon on an itinerary from the « Causses » region to the Mediterranean, Dominique Auerbacher in the area of Northern France crossed by the motorway and the HST, Gérard