M El Rhabi
email: elrhabi@gmail.com

H Fenniri
email: hassan.fenniri@univ-reims.fr

A Ghazdali
email: a.ghazdali@gmail.com

A Hakim

A Keziou
email: amor.keziou@univ-reims.fr

Blind noisy mixture separation for independent/dependent sources through a regularized criterion on copulas

Keywords: Blind source separation, Noisy instantaneous mixtures, Copulas, Total variation, Mutual information, Kullbak-Leibler divergence . . . . . . . . . . . . 10 dependent . . . . . . . . . . .

Blind noisy mixture separation for independent/dependent sources through a regularized criterion on copulas

Introduction

The blind source separation (BSS) problem is a fundamental issue in applications of many different fields such as signal and image processing, biological and medical data analysis, communications, ... etc. The BSS aims to recover unknown source signals from a set of observations which are unknown mixtures of source signals. In order to separate the mixtures, different assumptions on the sources have to be made. In the literature, the most common assumptions are statistical independence of the source components and the condition that at most one of the components is gaussian. Under these assumptions, the BSS problem is linked to the well known problem of Independent Component Analysis (ICA), see for instance [START_REF] Comon | Independent component analysis, a new concept?[END_REF]. Using the above conditions, particularly, for free-noisy mixtures, many methods of BSS have been proposed in the literature, see e.g. [START_REF] Mansour | A direct solution for blind separation of sources[END_REF][START_REF] Taleb | Entropy optimization[END_REF][START_REF] Hyvärinen | A fast fixed-point algorithm for independent component analysis[END_REF][START_REF] Pham | Blind separation of instantaneous mixture of sources based on order statistics[END_REF]. Some of these procedures use second or higher order statistics [START_REF] Pesquet | Cumulant-based independence measures for linear mixtures[END_REF],

maximize likelihood [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF], maximize nongaussianity [START_REF] Novey | Ica by maximization of nongaussianity using complex functions[END_REF], minimize the mutual information [START_REF] Pham | Mutual information approach to blind separation of stationary sources[END_REF] or ϕ-divergences [START_REF] Keziou | Séparations aveugle de sources par minimisation des α-divergences[END_REF]. An interesting overview on the subject can be found in [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. Recently, under

the assumption that the source components are independent, based on estimation of mutual information and total variation (TV) regularization, [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] provided a new BSS algorithm to separate noisy mixtures of instantaneous or covolutive mixtures of independent source components.

Further, it has been shown in [START_REF] Keziou | New blind source separation method of independent/dependent sources[END_REF] that, based on copula models, without the assumption of the independence of the source components, we can still identify both mixing matrix and sources uniquely (up to scale and permutation indeterminacies) of (free-noisy) mixtures of both independent and dependent source components. Motivated by various applications, we investigate, in the present paper, models of noisy linear instantaneous mixtures of independent/dependent sources, for which we propose, based on the previews works [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] and [START_REF] Keziou | New blind source separation method of independent/dependent sources[END_REF], a new BSS procedure.

We successfully combined the two methods of [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] and [START_REF] Keziou | New blind source separation method of independent/dependent sources[END_REF] to obtain a new BSS algorithm able to separate noisy mixtures of independent/dependent source components. Our methodology consists in denoising the observed signals through the minimization of their total variation, and then minimizing Kullback-Leibler divergence, between copula densities, penalized by the total variation of the estimated source signals. The outline of this paper is organized as follows.

In Section 2, we state the BSS problem in a noisy mixture context. In Section 3, we give a bref introduction on copulas and some of their main properties. In Section 4, we describe our approach. Section 5 presents how to implement the proposed approach using both numerical and statistical techniques. Section 6 gives some numerical results, illustrating the efficiency and robustness of the proposed method. Conclusions are drawn in Section 7.

The BSS model

BSS model in the noise-free case

BSS problem can be modeled as follows. Denoting A[•] the (unknown) mixing operator, the relationship between the observed and source signals can be written as

x(t) := A[s(t)] ∈ R p , (1)
where s(t) ∈ R q is the unknown vector of source signals to be estimated, and x(t) represents the observed vector signal at time t ∈ [0, T ]. The goal of BSS, is therefore to estimate the unknown sources s(t) from the observed mixtures x(t). The estimation is performed with no prior information about either the sources or the mixing operator A[•]. Specific restrictions are made on both the mixing model and the source signals in order to limit the generality. We will restrict our self to the case where the number of source component and the number of observed mixture ones are equal (p = q), and we assume, in the present paper, that the mixtures are linear and instantaneous, so that the mixing operator A can be considered as a p × p matrix.

In this case, supposing in addition that A is invertible, the candidate estimates of the sources will be obviously of the form

y(t) := B x(t) ∈ R p , (2)
where B ∈ R p×p represents an appropriate demixing matrix. In other words, the problem is to obtain an estimate, denote it B, "closing" to the ideal solution B = A -1 , by the use of only the observation x(t), which leads to accurate estimation of the source s(t):

s(t) := B x(t) s(t). (3)

BSS model for noisy mixtures

In the present paper, we will focus on the BSS problem for noisy linear instantaneous mixtures

x(t) := A s(t) + n(t) ∈ R p , t ∈ [0, T ], (4)
where x(t) ∈ R p is the vector of noisy observations, s(t) ∈ R p is the unknown vector of sources to be estimated, n(t) ∈ R p is the (unknown) noise, and A ∈ R p×p is the unknown mixing matrix. We assume that the source signal s(t) and the noise n(t) are independent. The mixing model (4) can also be written as

x(t) = x(t) + n(t), (5)
with

x(t) := A s(t) (6)
is the inaccessible noise-free mixed vector signals. The aim here, is to estimate the sources s(t) using only the observed noisy signal x(t). The presence of noise n(t) within the mixing model (4), as well as the possible dependency structure of the source components, complicate significantly the BSS problem. Unlike the above case of free-noise mixtures in Subsection 2.1, the separating system of the form

y(t) := B x(t) (7)
does not provide a valid estimation for the source signals s(•). In fact, the obtained signals according to the above system, can be written as

y(t) := B x(t) = B x(t) + B n(t) =: y(t) + n(t), (8)
where y(t) ∈ R p is the "noisy estimate" of s(t). That is the noisy estimated source y(t) is the sum of

y(t) := B x(t) (9)
the "noisy-free" estimate of the sources, which is inaccessible, and the unknown "noise"

n(t) := B n(t). (10)
Ideally, we would like to retrieve y(t) by denoising y(t), but it is rather difficult since the noise n(t) is unknown. During last years, several algorithms have been proposed to tackle the noisy BSS problem; In [START_REF] Sahmoudi | Robust approach for blind source separation in non-gaussian noise environments[END_REF], the authors propose a two-step approach by combining the fraction allower order statistic for the mixing estimation and minimum entropy criterion for noise-free source component estimation. In [START_REF] Belouchrani | Robust whitening procedure in blind source separation context[END_REF], a whitening procedure is proposed to reduce the noise effect. At our knowledge, the noisy BSS problem for possibly dependent sources has not been considered in the literature. In [START_REF] Keziou | New blind source separation method of independent/dependent sources[END_REF], the authors proposed a new criterion to successfully separate noise-free mixtures of both independent and dependent sources, based on estimation of Kullback-Leibler divergence between copula densities. [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] provides a BSS procedure for noisy mixtures using TV-regularization of mutual information, which applies exclusively under the assumption that the source components are independent. Based on these two last works, the present paper introduces a new procedure for estimating sources in the context of noisy mixtures of possibly dependent source components. It proceeds on two stages : (i) denoising the observed signal x(t) before demixing; (ii) a simultaneous BSS-denoisy procedure via minimizing a TV-regularized measure between copula densities.

Bref recall on Copula

When modeling multivariate distributions, one has to take into account the effects of the marginal distributions as well as the dependence between them. This can be achieved by using the copula approach, which allows to deal with the margins and the dependency structure separately. It was introduced by [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF] as a function which couples a joint distribution function with its univariate margins. Precisely, it can be presented as follows. Consider any random vector

Y := (Y 1 , . . . , Y p ) ∈ R p , p ≥ 1,
with joint distribution function (d.f.)

F Y (•) : y ∈ R p → F Y (y) := F Y (y 1 , . . . , y p ) := P(Y 1 ≤ y 1 , . . . , Y p ≤ y p ),
and continuous marginal d.f.'s

F Y j (•) : y j ∈ R → F Y j (y j ) := P(Y j ≤ y j ), ∀j = 1, . . . , p.
The characterization theorem of Sklar [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF] shows that there exists a unique p-variate function

C Y (•) : [0, 1] p → [0, 1], such that, F Y (y) = C Y (F Y 1 (y 1 ), . . . , F Yp (y p )), ∀y := (y 1 , . . . , y p ) ∈ R p .
The function C Y (•) is called a copula and it is in itself a joint d.f. on [0, 1] p with uniform margins. We have

∀u := (u 1 , . . . , u p ) ∈ [0, 1] p , C Y (u) = P F Y 1 (Y 1 ) ≤ u 1 , . . . , F Yp (Y p ) ≤ u p .
Conversely, for any marginal d.f.'s F 1 (•), . . . , F p (•), and any copula function C(•), the function

C(F 1 (•), . . . , F p (•)
) is a multivariate d.f. on R p . On the other hand, since the marginal d.f.'s F Y j (•), j = 1, . . . , p, are assumed to be continuous, then the random variables

F Y 1 (Y 1 ), . . . , F Yp (Y p )
are uniformly distributed on the interval [0, 1]. Moreover, we have that the components

Y 1 , . . . , Y p are statistically independent if and only if (iff) the copula C Y (•) of the random vector Y := (Y 1 , . . . , Y p ) ∈ R p writes C Y (u) = p j=1 u j =: C (u), ∀u ∈ [0, 1] p ,
which will be denote by C (•); It is called the copula of independence. Define, when it exists, the copula density

c Y (u) := ∂ p C Y (u) ∂u 1 • • • ∂u p , ∀u ∈ [0, 1] p .
Hence, the copula density of independence, denote it by c (•), is the function taking the value 1 on [0, 1] p and zero otherwise, namely,

c (u) := 1 [0,1] p (u), ∀u ∈ [0, 1] p . (11)
Let f Y (•), if it exists, be the probability density of the random vector

Y := (Y 1 , . . . , Y p ) ,
and, respectively, f Y 1 (•), . . . , f Yp (•), the marginal probability densities of the random variables Y 1 , . . . , Y p . Then, a straightforward computation shows that, for all vector y := (y 1 , . . . , y p ) ∈ R p , we have the relation

f Y (y) = p j=1 f Y j (y j ) c Y F Y 1 (y 1 ), . . . , F Yp (y p ) , (12)
As previously highlighted, copulas play an important role in the construction of multivariate d.f.'s. Therefore, several investigations have been carried out concerning the construction of different families of copulas and their properties. For more details on modeling theory as well as surveys of the commonly semiparametric used copulas, we can refer to the monographs by [START_REF] Nelsen | An introduction to copulas, 2nd Edition[END_REF][START_REF] Joe | Multivariate models and dependence concepts[END_REF].

A combined TV-Copula approach for BSS

Our TV-Copula approach for BSS proceeds in two steps. (i) uses the TV regularization technique for denoising the observed signals. (ii) separates the mixtures by minimizing penalized Kullback-Leibler divergence between copulas with a TV-regularization term. In the following, we describe these two steps.

Setp 1 : Denoising the observed signals

As in [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF], we would obtain the denoisy mixture signal

x(t) := (x 1 (t), . . . , x p (t)) , t ∈ [0, T ], (13)
from the (noisy) observed one x(t), see [START_REF] Hyvärinen | A fast fixed-point algorithm for independent component analysis[END_REF], by means of the following variational minimization problem

x i (•) = arg min w i ∈Xc 1 T T 0 1 2 (w i (t) -x i (t)) 2 dt + λ 1 T T 0 φ (|w i (t)|) dt , i = 1, . . . , p, (14)
where λ > 0 is the parameter of penalization, φ(•) : R + → R + is a well chosen function, w i (t) is the first derivative of the signal w i (t) with respect to time t, and X c is some appropriate space of signals. In practice, in order to reduce the noise, the function φ(•) is chosen to encourage smoothing in regions where the variations of the signal are weak, that is |w i | ≈ 0, and to preserve discontinuities where |w i | is strong. The total variation case corresponds to the choice φ(x) = x which will be adopted here. One can make other choices like φ(x) = √ 1 + x 2 . A natural space when dealing with the above continuous variational problem [START_REF] Sahmoudi | Robust approach for blind source separation in non-gaussian noise environments[END_REF], with φ(x) = x,

is X c = BV ([0, T ]
), the space of all real valued functions on the interval [0, T ] with bounded variation; see e.g. [START_REF] Evans | Measure theory and fine properties of function[END_REF] and [START_REF] Cohen | Harmonic analysis of the space bv[END_REF]. Here T is the observation time of the signals. The mathematical treatment of the above problem, in its continuous form, is beyond the scope of this paper. We will consider only the corresponding discrete version which will be described in Subsection 5.1 below.

Step 2 : BSS via minimizing KL-divergence between copulas and denoising

The aim of the following step is to construct an estimated source signal

s(t) := B x(t)
from the denoised signal x(t) := (x 1 (t), . . . , x p (t)) obtained as solution of [START_REF] Sahmoudi | Robust approach for blind source separation in non-gaussian noise environments[END_REF]. Our approach consists in minimizing, with respect to B on the demixing matrices space, an estimate (to be defined below) of some criterion of the form

B → J (B) := J sep (B) + J reg (y), (15)
where

y(t) := B x(t) ∈ R p ,
and x(t) = (x 1 (t), . . . , x p (t)) ∈ R p is obtained from ( 14); J sep (•) is some separating criterion, while J reg (y) is a regularization term devoted to denoising. For the denoising phase, for both cases, of independent or dependent source components, we propose to use as regularization term the following one

J reg (y) := γ 1 2T T 0 y(t) -y(t) 2 dt + µ 1 T T 0 |∇y(t)| dt, γ > 0, µ > 0, (16)
where y = B x; the real numbers γ and µ are regularization parameters to be suitably chosen by the user. In the above display, • is used to denote the Euclidian norm on R p , and

|∇y(t)| = p i=1 |y i (t)|
, for all t. The first term J sep (•), in [START_REF] Belouchrani | Robust whitening procedure in blind source separation context[END_REF], is a measure of dependence, or similarity between copula densities, and it depends obviously on the fact that the source components are independent or dependent. Hence, to describe the criterion J sep (•) hereafter, we will consider separately, the case where the source components are independent and the case where they are dependent. First, we will consider the following stochastic modeling of the analog signals s(t), x(t) and y(t) := Bx(t), t ∈ [0, T ]. We assume that the above last random processes are stationary, so that the corresponding sampled versions, with certain time period, say T e , s(n) := s(nT e ), x(n) := x(nT e ) and y(n

) := y(nT e ) = B x(n), n = 1, . . . , N, (17)
can be considered as copies (realizations) of random vectors in R p , which will be denoted S, X and Y := BX,

respectively. It has been shown, in [START_REF] Keziou | New blind source separation method of independent/dependent sources[END_REF], for the noisy-free mixtures case, that if we dispose of some prior information about the copula density of the random source vector S, we can detect, under some assumptions, both the mixing matrix and the sources uniquely (up to scale and permutation indeterminacies) for both cases of independent or dependent source components.

Denote by

F Y (•) : y ∈ R p → F Y (y) := F Y (y 1 , . . . , y p ) := P(Y 1 ≤ y 1 , . . . , Y p ≤ y p ), (19)
the joint distribution function of the random vector Y := (Y 1 , . . . , Y p ) ∈ R p , with continuous marginal distribution functions

F Y i (•) : y i ∈ R → F Y i (y i ) := P(Y i ≤ y i ), ∀i = 1, . . . , p.
Let, if it exists, f Y (•) be the probability density of the random vector Y , and, respectively,

f Y 1 , . . . ,
I(Y ) := R p log f Y (y) p i=1 f Y i (y i ) f Y (y) dy = E log f Y (Y ) p i=1 f Y i (Y i ) (20)
which equals the Kullback-Leibler divergence, denote it K(•, •), between the joint density f Y (•) and the product density p i=1 f Y i (•) of the margins f Y i (•), i = 1 . . . , p, that is

I(Y ) = K f Y , p i=1 f Y i . (21)
In ( 20), E(•) is the mathematical expectation. Note also that the above criterion I(Y ) is nonnegative and achieves its minimum value zero iff

f Y (•) = p i=1 f Y i (•), i.e.
, iff the components of the vector Y are statistically independent. We give now, the form of the separating term J sep (•) of the criterion (15).

The case of independent source components

Assume that the source components are independent. Using the relation [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF], and applying the change variable formula for multiple integrals, we can show that I(Y ) can be written, in terms of copula densities, as follows The separating term will be chosen then to be

I(Y ) = [0,1] p log (c Y (u)) c Y (u) du (22) = E log c Y F Y 1 (Y 1 ), . . . , F Yp (Y p ) =: -H(c Y ) (23) = [0,1] p log c Y (u) 1 c Y (u) du (24)
= [0,1] p log c Y (u) c (u) c Y (u) du (25) = E log c Y (F Y 1 (Y 1 ), . . . , F Yp (Y p )) c (F Y 1 (Y 1 ), . . . , F Yp (Y p )) = K (c Y , c ) , (26)
J sep := J ind sep (B) := KL (c Y , c ) = E log c Y F Y 1 (Y 1 ), . . . , F Yp (Y p ) . ( 27
)
We have that the function B → J ind sep (B) is nonnegative and attains its minimum value zero iff B = A -1 (up to scale and permutation indeterminacies of rows).

The case of dependent source components

In the case where the source components are dependent, we assume that we dispose of some prior information about the copula density of the random source vector S. Note that this is possible for many practical problems, it can be done, from learning samples of S, by a model selection procedure, see e.g. [START_REF] Chen | Estimation and model selection of semiparametric copula-based multivariate dynamic models under copula misspecification[END_REF], among semiparametric copula density models {c θ (•); θ ∈ Θ ⊂ R d }, typically indexed by a multivariate parameter θ; see e.g. [START_REF] Joe | Multivariate models and dependence concepts[END_REF] and [START_REF] Nelsen | An introduction to copulas, 2nd Edition[END_REF] for many examples of such models. The parameter θ can be estimated using maximum semiparametric likelihood, see e.g. [START_REF] Tsukahara | Semiparametric estimation in copula models[END_REF]. Denote then by θ, the obtained value of θ and c θ (•) the copula density modeling the dependency structure of the source components. Obviously, since the source components are assumed to be dependent, c θ (•) should be different from the density copula of independence c (•). Hence, we naturally replace in [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], c by c θ , and propose the following separating criterion J sep (•)

J sep (B) := J dep sep (B) := [0,1] p log c Y (u) c θ (u) c Y (u) du := E log c Y (F Y 1 (Y 1 ), . . . , F Yp (Y p )) c θ (F Y 1 (Y 1 ), . . . , F Yp (Y p )) , (28)
which equals K(c Y , c θ ), the Kullback-Leibler divergence between the copula densities c Y (•) and c θ (•). Moreover, we can show that the function B → J dep sep (B) is nonnegative, and attains its minimum value iff B = A -1 (up to scale and permutation indeterminacies of rows), provided that the copula density c θ (•) of the random vector source S satisfies the following assumption : for any regular matrix M , if the copula density of the random vector M S equals c θ (•), then M = D P , for some diagonal matrix D and permutation matrix P .

Remark 1. The criterion (29) proposed in the above section supposes the knowledge of the copula density model c θ (•) (with known parameter θ), or the availability of learning samples of the random vector source S by which we can obtain the copula density c θ (•) modeling the dependency structure of the source components, as described in Subsection 4.2.2. When only the model is known and the parameter is unknown, we can adapt the above criterion as follows.

Denote by {c θ (•); θ ∈ Θ} the density copula model of the dependency structure of the sources.

The criterion (29) can be in this case replaced by

J sep (B) := inf θ∈Θ [0,1] p log c Y (u) c θ (u) c Y (u) du := inf θ∈Θ E log c Y (F Y 1 (Y 1 ), . . . , F Yp (Y p )) c θ (F Y 1 (Y 1 ), . . . , F Yp (Y p )) . ( 29
)
In the more general case, when both the model and the parameter are unknown, the criterion may be adapted in a similar way as in our previous work [START_REF] Keziou | New blind source separation method of independent/dependent sources[END_REF] Section 5.3.

Algorithm 1 The denoising step.

1: Data: X the observed noised signal 2: Result: X the denoised signal 3: Initialization: q = 0, X (0) = X. Given > 0, λ > 0 and ρ > 0 suitably chosen 4: Do: q = q + 1 5:

for j = 1, . . . , p 6:

X (q) j = λG (X (q-1) j + 2λρ(λ(X (q-1) j ) * -X j ) 7:
end

8: Until ||X (q) -X (q-1) || < 9: X = X (q) .
where J ind sep (B) is given by, see (27),

J ind sep (B) := KL (c Y , c ) = E log c Y F Y 1 (Y 1 ), . . . , F Yp (Y p )
and

J reg (y) := γ 1 2T T 0 y(t) -y(t) 2 dt + µ 1 T T 0 |∇y(t)| dt, γ > 0, µ > 0,
with y(t) = B x(t) and y = B x(t). Using the stochastic modeling [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF] and the relation (27), we propose to approximate the criterion (35) by

B → J ind (B) := J ind sep (B) + J reg,d (y), (36)
where

J reg,d (y) := γ 2N N i=1 y(i) -y(i) 2 + µ N N i=1 |∇y(i)| (37)
is the discrete version of J reg (y), and B → J ind sep (B) is the statistical estimate, of the separating criterion B → J ind sep (B), which we define as follows

B → J ind sep (B) := 1 N N i=1 log c Y ( F Y 1 (y 1 (i)), . . . , F Yp (y p (i))) , (38)
where

c Y ( F Y 1 (y 1 (i)), . . . , F Yp (y p (i))) := 1 N H 1 • • • H p N =1 p j=1 k F Y j (y j (i)) -F Y j (y j ( )) H j ,
is the kernel estimate of the copula density c Y (•), and F Y j (•), ∀j = 1, . . . , p, is the estimate of the marginal distribution function F Y j (•) of the random variable Y j , defined for any real number r ∈ R, by

F Y j (r) := 1 N N =1 K r -y j ( ) h j ,
where K(•) is the primitive of a kernel k(•), a symmetric centered probability density. In our forthcoming simulation study, we will take as kernel k(•) a standard Gaussian probability density. A more appropriate choice of the kernel k(•), for estimating the copula density, can be done according to [START_REF] Omelka | Improved kernel estimation of copulas: weak convergence and goodness-of-fit testing[END_REF], which copes with the boundary effect. The bandwidth parameters H 1 , . . . , H p and h 1 , . . . , h p will be chosen according to Silverman's rule of thumb, see [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], i.e., for all j = 1, . . . , p, we take

H j = 4 p + 2 1 p+4 N -1
p+4 Σ j , and h j = 4 3

1 5 N -1
5 σ j , where Σ j and σ j are, respectively, the empirical standard deviation of the data F Y j (y j (1)), . . . , F Y j (y j (N )) and y j (1), . . . , y j (N ). The source vector signal s(i), i = 1, . . . , N , will be then estimated by

s(i) = B x(i), i = 1, . . . , N,
where

B := arg inf B J ind (B),
which can be computed using a gradient descent type algorithm. In fact, straightforward computation shows that, the gradient in B of the estimated criterion B → J ind (B), can be written as

d J ind (B) dB = 1 N N i=1 d dB c Y (u(i)) c Y (u(i)) + γ N N i=1 (y(i) -y(i))(x(i) -x(i)) + µ N N i=1 ∇y(i) |∇y(i)| * x(i) (39)
where,

d dB := ∂ ∂B l,j
, l, j = 1, . . . , p, u(i) := ( F Y 1 (y 1 (i)), . . . , F Yp (y p (i))) and,

∂ c Y ( F Y 1 (y(i)), . . . , F Yp (y p (i))) ∂B l,j = 1 N H 1 • • • H p N m=1 p j=1,j =i k F Y j (y j (m)) -F Y j (y j (i)) H j × k F Y j (y j (m)) -F Y j (y j (i)) H j 1 H j ∂( F Y j (y j (m)) -F Y j (y j (i))) ∂B l,j , (40)
with

∂( F Y j (y j (m)) ∂B l,j = 1 N h j N i=1 k y j (i) -y j (m) h j (x j (i) -x j (m)). ( 41
)
We can then derive the following Algorithm 2.

Algorithm 2 BSS algorithm for independent source components. Data: x the observed signal

Result: s the estimated source signal

Initialization: Compute x = x -D x from Algorithm 1, B (0) = I p , y (0) = B (0) x. Given ε > 0, ν > 0.
Do: Update B and y:

B (q+1) = B (q) -ν d J ind (B) dB y (q+1) = B (q+1) x.
Until ||B (q+1) -B (q) || < s = y (q+1) .

The case of dependent source components

Recall that, in the case where the source components are dependent, the criterion function [START_REF] Belouchrani | Robust whitening procedure in blind source separation context[END_REF] is defined by

B → J dep (B) := J dep sep (B) + J reg (y), (42)
where J dep (B) is defined by, see (29),

J dep sep (B) := [0,1] p log c Y (u) c θ (u) c Y (u) du := E log c Y (F Y 1 (Y 1 ), . . . , F Yp (Y p )) c θ (F Y 1 (Y 1 ), . . . , F Yp (Y p )) .
and

J reg (y) := γ 1 2T T 0 y(t) -y(t) 2 dt + µ 1 T T 0 |∇y(t)| dt, γ > 0, µ > 0.
Using the stochastic modeling [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF] and the relation (29), we propose to approximate the criterion (42) by

B → J dep (B) := J dep sep (B) + J reg,d (y), (43)
where

J reg,d (y) := γ 2N N i=1 y(i) -y(i) 2 + µ N N i=1 |∇y(i)| (44)
is the discrete version of J reg (y), and B → J dep sep (B) is the statistical estimate, of the criterion B → J dep sep (B), which we define as follows

B → J dep sep (B) := 1 N N i=1 log c Y ( F Y 1 (y(i)), . . . , F Yp (y p (i))) c θ ( F Y 1 (y(i)), . . . , F Yp (y p (i))) . ( 45
)
The source vector signal s(i), i = 1, . . . , N , will be then estimated by

s(i) = B x(i), i = 1, . . . , N,
where

B := arg inf B J dep (B),
which can be computed using a gradient descent type algorithm. In fact, straightforward computation shows that, the gradient in B of the estimated criterion B → J dep (B), can be written as

d J dep (B) dB = 1 N N i=1 d dB c Y (u(i)) c Y (u(i)) c θ (u(i)) d dB c θ (u(i)) + γ N N i=1 (y(i) -y(i))(x(i) -x(i)) + µ N N i=1 ∇y(i) |∇y(i)| * x(i) , (46)
where u(i) := ( F Y 1 (y(i)), . . . , F Yp (y p (i))); the gradients d dB c Y (u(i)) and d dB c θ (u(i)) can be explicitly computed in a similar way as in Subsection 5.2.1 above. We obtain then the following Algorithm 3. Algorithm 3 BSS algorithm for dependent source components.

Data: x the observed signal Result: s the estimated source signal

Initialization: Compute x = x -D x from Algorithm 1, B (0) = I p , y (0) = B (0) x. Given ε > 0, ν > 0.
Do: Update B and y:

B (q+1) = B (q) -ν d J dep (B) dB y (q+1) = B (q+1) x. Until ||B (q+1) -B (q) || < s = y (q+1) .
We dealt with noisy instantaneous mixtures of five kinds of sample sources: i.i.d. random vector sources (with uniform marginals in Fig. 5, and binary phase-shift keying (BPSK) marginals in Fig. 6), with dependent components generated from Ali-Mikhail-Haq (AMH) copula model where θ = 0.6. i.i.d. vector sources (with uniform marginals in Fig. 7, and binary phaseshift keying (BPSK) marginals in Fig. 8) where the components are dependent, generated from Clayton copula with θ = 1. i.i.d vector sources (with uniform marginals) with dependent components, generated from Fairlie-Gumbel-Morgenstern (FGM) copula model where θ = 0.8 (Fig. 9). It can be seen from the results, in Fig. 5 to Fig. 9, that the proposed method is able to separate, with good performance (with level equivalent to that of independent source case, around 27 dB). Fig. 10 shows the convergence of the criterion J Dep (•) to 0 when the separation is achieved for the last example. We compare, for the following two examples of dependent mixtures, the present approach (Algorithm 2) and the MI-TV of [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF]. We consider i.i.d. vector sources (with uniform marginals in Fig. 14, and BPSK marginals in Fig. 15) where the component are dependent, generated from FGM copula (with θ = 0.8). We can see, clearly, that separating dependent sources -with classical criterion devoted to the independent casefails, and that our criterion is able to deal with the dependent case with good performance.

Conclusion

We have presented a new BSS approach in presence of noise, by minimizing a new regularized criterion. The approach is able to separate instantaneous linear mixtures of both independent and dependent source components. It proceeds in two steps: a preprocessing which consists in reducing the noise on the observations, and a second one minimizing a separation criterion based on Kullback-Leibler divergence between copula densities regularized through total variation. In Section 7, the accuracy and the consistency of the obtained algorithms are illustrated by simulation, for 2 × 2 mixture-source. It should be mentioned, however, that our proposed algorithms based on copula densities, rather than the classical ones based on probability densities, are more time consuming, since we estimate both copulas density of the vector and the marginal distribution function of each component. The present approach can be extended to deal with convolutive mixtures of independent/dependent sources, that will be addressed in future communications.

  model in the noise-free case . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

  where c Y (•) is the copula density of Y , c (•) := 1 [0,1] p (•) is the copula density of independence, and E(•) the mathematical expectation operator. The above equations mean that the MI of a random vector Y can be seen as the opposite of the Shannon entropy H(•) of the copula density c Y (•) of Y , or as the Kullback-Leibler divergence between the copula density c Y (•) of Y and the copula density c (•) of independence. Moreover, I(Y ) = K (c Y , c ) is nonnegative and achieves its minimum value zero iff c Y (u) = c (u), ∀u ∈ [0, 1] p , namely, iff the components of the vector Y are independent.

Figure 1 :

 1 Figure 1: Average output SNRs versus iteration number : Uniform independent sources.

Figure 2 :

 2 Figure 2: Average output SNRs versus iteration number : ASK independent sources.

Figure 3 :

 3 Figure 3: The criterion value vs iterations : uniform independent sources.

Figure 4 :

 4 Figure 4: The criterion value vs iterations : Bpsk independent sources.

Figure 5 :

 5 Figure 5: Average output SNRs versus iteration number : Uniform dependent sources from AMH-copula.

Figure 6 :

 6 Figure 6: Average output SNRs versus iteration number : Bpsk dependent sources from AMH-copula.

Figure 7 :

 7 Figure 7: Average output SNRs versus iteration number : Uniform dependent sources from Clayton-copula.

Figure 8 :

 8 Figure 8: Average output SNRs versus iteration number : Bpsk dependent sources from Clayton-copula.

Figure 9 :

 9 Figure 9: Average output SNRs versus iteration number : Uniform dependent sources from FGM-copula.

Figure 10 :

 10 Figure 10: The criterion value vs iterations : Uniform dependent sources from FGM-copula.

Figure 11 :

 11 Figure 11: Average output SNRs versus iteration number: Uniform independent sources.

Figure 12 :

 12 Figure 12: Average output SNRs versus iteration number: Bpsk independent sources.

Figure 13 :

 13 Figure 13: Average output SNRs versus iteration number: Uniform dependent sources from FGM-copula.

Figure 14 :

 14 Figure 14: Average output SNRs versus iteration number: Bpsk dependent sources from FGM-copula.

Figure 15 :

 15 Figure 15: Average output SNRs versus iteration number: Uniform dependent sources from Clayton-copula.

Figure 16 :

 16 Figure 16: Average output SNRs versus iteration number: Bpsk dependent sources from AMH-copula.

f

  Yp , the marginal probability densities of the components Y 1 , . . . , Y

p . The mutual information (MI) of Y is defined by

Discretization and Statistical estimation

In this section, we describe how to make operational the above methodology. Precisely, using both numerical and statistical techniques, we will describe hereafter the discrete versions of the regularization terms [START_REF] Sahmoudi | Robust approach for blind source separation in non-gaussian noise environments[END_REF] and [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF], as well as statistical estimates of the separation criterions ( 27) and ( 29). In all the sequel, analog signals are sampled with a certain period T e . We need to define additional notation. To each scalar signal, say u(t), t ∈ [0, T ], one associate a vector

) for all k = 1, . . . , N. The vector space X d := R N will be equipped with the euclidian inner product Here, the backward derivatives (•) * is the discrete adjoint operator of -(•) . That is, for all U, V ∈ X d , we have

All these definitions can be extended in a natural way to elements of the cartesian product space X p = R p × R N . In the sequel, elements of the vector space X d = R N will be denoted by double capital symbols, while the elements of the vector space X p d are denoted by bold double capital symbols, which can be considered as p × N dimension matrices.

Denoising the discrete observed signal

In this section, we show how to compute the denoisy observations

from the observed noisy ones

This will be obtained by solving the discrete version

of the optimization problem [START_REF] Sahmoudi | Robust approach for blind source separation in non-gaussian noise environments[END_REF]. We start with the following proposition which gives a simple characterization of the solution; see [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] for a proof.

Proposition 5.1. The discrete version (32) of the optimization problem ( 14) has a unique solution given by

where λG is the euclidian orthogonal projection operator on the convex set λG with

Thus, to compute X, we are led to compute the projection operator λG on the convex set λG. For the denoising step, we propose then the algorithm 1 below; the reader may refer to [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] for more details, particularly, for a proof on convergence of this algorithm.

In all the sequel, the filtered observed signal X will be considered as the denoised version of the observed signal X.

BSS by minimizing an estimate of KL-divergence, between copulas, penalized by TV regularization

The case of independent source components

Recall that, when the source components are independent, the criterion function ( 15) is defined by

Simulation results

In this section, we present simulation results illustrating the performance of our proposed method. For simplicity, we will deal with only the case of two mixtures/two sources. We will consider many examples of independent/dependent source component signals. The results will be compared with those of [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] under the same conditions. Recall that the separating algorithm proposed in [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] is based on minimizing a kernel type estimate of mutual information between probability densities, penalized by the same TV-regularization term. In all examples the number of samples is N = 2000. The two sources are mixed with the matrix A := [1 0.5; 0.5 1]. A centered gaussian noise with standard deviation 0.01 was added to the normalized mixtures, so that the signal-to-noise ratio equals -20 dB. The gradient descent parameter is taken to be ν = 0.1. In the denoising step, see Algorithm 1, we take ρ = 0.1, λ = 0.01, and = 0.001.

In the second step, see Algorithm 2 and Algorithm 3, we chose γ = 0.001 and µ = 0.01. All simulations are repeated 100 times, and the accuracy of source estimation is evaluated through the following criterion, called again the signal-to-noise-ratio (SNR), defined by

(47)

Examples where the source components are independent

We consider two noisy mixed signals of two kinds of sample sources: uniform i.i.d with independent components (Fig. 1); i.i.d sources with independent components drawn from the 4-ASK (Amplitude Shift Keying) alphabets (-3,-1,1,3 with the same weights 0.25) at which was added a centered gaussian random variable with variance equals 0.25 (Fig. 2). We observe from Fig. 1 and Fig. 2, that our proposed method (Algorithm 2) gives accurate estimation of sources in the standard case of independent component sources, with performance around 27 dB in term of SNR. Fig. 3 and Fig. 4 shows the criterion values versus iterations. Observe that the criterion values J Ind (•) converges to 0 when the separation is achieved. In Fig. 11 and Fig. 12, we compare Algorithm 2 with the MI-TV presented in [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF], to separate the same instantaneous mixtures with independent sources. We can see, from figures 11 and 12, that the two approches are equivalent, in this standard case, as expected.

Examples where the source components are dependent

In this subsection, we show the capability of the proposed algorithm using BSS-copula-TV (Algorithm 3) to successfully separate two noisy mixed source signals with dependent components.


Bachirou O Taddé

Hélène Jacqmin-Gadda

Jean-Franois Dartigues

Daniel Commenges

Cécile Proust-Lima

Dynamic modeling of multivariate latent processes and their temporal relationships: Application to Alzheimer's disease

Keywords: . causality, mixed models, difference equations, latent process, longitudinal data. 1

Alzheimer's disease gradually affects several components including the cerebral dimension with brain atrophies, the cognitive dimension with a decline in various functions and the functional dimension with impairment in the daily living activities. Understanding how such dimensions interconnect is crucial for AD research. However it requires to simultaneously capture the dynamic and multidimensional aspects, and to explore temporal relationships between dimensions. We propose an original dynamic model that accounts for all these features. The model defines dimensions as latent processes and combines a multivariate linear mixed model and a system of difference equations to model trajectories and temporal relationships between latent processes in finely discrete time. Parameters are estimated in the maximum likelihood framework enjoying a closed form for the likelihood. We demonstrate in a simulation study that this dynamic model in discrete time benefits the same causal interpretation of temporal relationships as mechanistic models defined in continuous time. The model is then applied to the data of the Alzheimer's Disease Neuroimaging Initiative. Three longitudinal dimensions (cerebral anatomy, cognitive ability and functional autonomy) are analyzed and their temporal structure is contrasted between different clinical stages of Alzheimer's disease.

Introduction

Dementia is a general syndrome characterized by a long term and gradual decrease in the ability to think and remember with consequences on the person's daily functioning. It represents a pressing public health problem and major research priority. Alzheimer disease (AD) is the most common form of dementia, representing 60% to 80% of the cases [START_REF] Reitz | Alzheimer's disease: epidemiology, diagnostic criteria, risk factors and biomarkers[END_REF]. AD gradually affects multiple components long before any clinical diagnosis including brain atrophies, cognitive decline in various functions (memory, language, orientation in space and time, etc.) and loss of autonomy in the daily living activities. There is currently no cure for AD although many novel compounds are under development. Natural history of AD and its progression are still misunderstood; yet their understanding constitutes a crucial step for the assessment of compounds, and more generally AD research. Theoretical schemes have been proposed to understand the disease taken as a whole, highlighting the expected dynamic and multidimensional aspects [START_REF] Jack | Tracking pathophysiological processes in Alzheimer's disease: an updated hypothetical model of dynamic biomarkers[END_REF]). Jack's scheme, which is now central in AD research, hypothesizes a sequence of alterations: mainly the accumulation of proteins in the brain (amyloid-β and tau proteins), the atrophy of brain regions (e.g., hippocampus) and finally clinical manifestations with cognitive and functional declines on which the diagnosis of dementia currently relies. However, because of its complexity, such theoretical scheme is very difficult to translate into a statistical model as it necessitates to combine the dynamic and multidimensional aspects, and to properly explore the temporal relationships between dimensions.

The dynamic aspect of AD has been apprehended mostly using the linear mixed model theory with models that analyzed one dimension according to time and covariates [START_REF] Hall | A change point model for estimating the onset of cognitive decline in preclinical Alzheimer's disease[END_REF][START_REF] Amieva | Prodromal Alzheimer's disease: successive emergence of the clinical symptoms[END_REF][START_REF] Donohue | Estimating long-term multivariate progression from short-term data[END_REF]. One difficulty is that the dimension of interest is usually not directly observed but measured by a series of outcomes; for instance, cognition is usually measured by a battery of neuropsychological tests. [START_REF] Proust-Lima | A nonlinear model with latent process for cognitive evolution using multivariate longitudinal data[END_REF], 2013) considered latent process models which distinguish a structural model for analyzing the unobserved quantity over time and equations of observations for linking the unobserved quantity to the observed outcomes.

To account for the multidimensional aspect of AD and understand how dimensions are inter-related, many contributions explored pre-determined relationships by examining change over time of one biomarker (e.g., cognitive measure) according to another observed biomarker (e.g., hippocampal initial volume or change) and assumed the latter was ob-served without measurement error (e.g. [START_REF] Landau | Associations between cognitive, functional, and fdg-pet measures of decline in AD and MCI[END_REF]; [START_REF] Han | Beta amyloid, tau, neuroimaging, and cognition: sequence modeling of biomarkers for alzheimer's disease[END_REF]). This approach quantifies temporal relationships but it relies on a specific a priori determined sequence and does not consider all biomarkers as stochastic processes, which limits the interpretation and may induce biases. Some authors also proposed to use bivariate mixed models to dynamically model two dimensions and account for their correlation through correlated random effects (e.g., [START_REF] Mungas | Longitudinal volumetric MRI change and rate of cognitive decline[END_REF][START_REF] Robitaille | Multivariate longitudinal modeling of cognitive aging: Associations among change and variation in processing speed and visuospatial ability[END_REF]). However such models remain descriptive; they do not rely on asymmetric temporal relationships between processes and thus do no allow causal interpretations.

Temporal asymmetric relationships between processes have usually been apprehended with Dynamic Bayesian Networks (DBN) [START_REF] Song | Time-varying dynamic bayesian networks[END_REF]. This approach extends the concept of Directed Acyclic Graphs (DAG) [START_REF] Greenland | Causal Analysis in the Health Sciences[END_REF] to the dynamic case by modelling constant or time-varying temporal relationships between successive states of a network of processes. However, firstly time is usually too grossly discretized so that spurious temporal associations might appear as showed by [START_REF] Aalen | Can we believe the dags? a comment on the relationship between causal dags and mechanisms[END_REF]. Second, the trajectories of processes are not modelled over time. Finally, DBN quantify the association between successive levels of processes while a dynamic view of causality would seek local dependence structures linking the network of processes to its infinitesimal change over time [START_REF] Aalen | What can statistics contribute to a causal understanding[END_REF].

Local dependence structures can be naturally investigated with mechanistic models which are dynamic models relating a system of processes over time using differential equations. Proposed in HIV studies [START_REF] Prague | Dynamic models for estimating the effect of HAART on CD4 in observational studies: application to the Aquitaine Cohort and the Swiss HIV Cohort Study[END_REF] or cancer studies [START_REF] Desmée | Nonlinear joint models for individual dynamic prediction of risk of death using Hamiltonian Monte Carlo: application to metastatic prostate cancer[END_REF], they allowed retrieving causal associations between disease components. Yet mechanistic models are numerically very demanding so that their application to complex diseases such as AD is compromised. In addition they require precise biological knowledge which lacks for AD.

The objective of this work is to propose a statistical model that simultaneously describes the dynamics of multiple dimensions involved in Alzheimer's disease and assess their temporal relationships similarly as in a mechanistic model but with much less numerical complexity. We consider a system of latent processes, each one representing a latent dimension possibly observed through one or several longitudinal markers. In contrast with mechanistic models, we rely on a discrete-time framework and define a set of difference equations to model the change over time of the system according to the previous state of the system. In addition, we account for individual differences through subject-specific covariates and random effects. As discretization can distort the causal interpretations of temporal relationships compared to a model in continuous time, we specifically evaluate the impact of discretization on the temporal influence structure in a simulation study. The methodology is applied to the Alzheimer's Disease Neuroimaging Initiative database to explore the temporal structure between cerebral, cognitive and functional dimensions at different stages of AD.

Motivating Data

Data used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). The ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). For up-to-date information, see www.adni-info.org. We focused in this work on the ADNI 1 phase which is a longitudinal multisite observational study that included approximately 800 individuals aged between 55-90 enrolled in three different stages of progression to AD: normal aging (CN, N ≈ 200), mild cognitive impairment (MCI, N ≈ 400), and diagnosed Alzheimer's Disease (dAD, N ≈ 200). The individuals were followed-up every 6 months up to 3 years for CN and MCI groups and up to 2 years for the dAD group. At each visit, a lot of information was collected, including for the present work, volumes of brain regions measured by MRI, a battery of 19 cognitive tests and a functional assessment with the Functional Activities Questionnaire(FAQ) [START_REF] Pfeffer | Measurement of functional activities in older adults in the community[END_REF]. Further details on the ADNI 1 study can be found in [START_REF] Mueller | The Alzheimer's disease Neuroimaging Initiative[END_REF].

Methodology

3.1 Structural model for the system of latent processes

Consider D latent processes Λ i (t) t≥0 (with Λ i (t) = (Λ d i (t)) d=1,.
..,D ) representing a system of D dimensions (e.g. cerebral anatomy, cognitive ability and functional autonomy dimensions for Alzheimer's disease) for individual i with i = 1, . . . , N . We assume Λ i (t) t≥0 is defined at discrete times t = j × δ with j = 0, 1, . . . , J, and δ a constant discretization step. Let us denote ∆Λ i (t + δ) = Λ i (t + δ) -Λ i (t) the change of the system between two successive times and ∆Λ i (t+δ) δ the rate of change of the system. The model for the trajectories of the processes is split in two parts: (i) the level of the processes at baseline Λ i (0) is modelled with a multivariate linear mixed model, and (ii) the rate of change of the system over time ∆Λ i (t+δ) δ is modelled by difference equations combined with a multivariate linear mixed model:

Λ i (0) = X 0 i β + u i ∆Λ i (t+δ) δ = X i (t + δ)γ + Z i (t + δ)v i + A i,δ (t)Λ i (t), ∀t > 0, (1)
where X 0 i is the D × p 0 -matrix of covariates associated with the p 0 -vector of fixed effects β, and u i is the D-vector of individual random intercepts u d i in the initial system Λ i (0). The (D × p)-matrix X i and the (D × q)-matrix Z i include time-dependent covariates associated with the p-vector of fixed effects γ and the q-vector (q

= D d=1 q d ) of individual random effects v i = (v d i ) d=1,...,D , respectively. A i,δ (t) is the D × D-matrix of transition intensities.
For each dimension d, the (q d + 1) vector of individual random effects (u d i , v d i ) is assumed to have a multivariate normal distribution with variance-covariance matrix

1 B d uv B d uv B d v
where B d uv and B d v are unstructured. The variances of the u d i s are constrained to 1 so that the D processes have the same magnitude at baseline. Random effects are assumed independent between dimensions so that the entire (D + q)-vector of individual random effects w i = (u i , v i ) has a multivariate normal distribution,

w i ∼ N 0 0 , B = I D B uv B uv B v .
with I D the D×D identity matrix, B v the D-block diagonal matrix with d th block B d v , and B uv the D ×q matrix with

d th row O d-1 l=1 q l , B d uv , O D l=d+1 q l
where O x is the x-row vector of zeros. In the estimation process, B is replaced by its Cholesky decomposition: B = LL , where L is a (D + q) × (D + q) lower triangular matrix. Independence between random effects across dimensions is handled by fixing corresponding parameters to zero in L, and unit-variances of random intercepts are handled by fixing the corresponding parameters to 1 in L. In addition to the variances of random intercepts fixed at one, and without loss of generality, the matrix of covariates X 0 i excludes intercepts for each process so that the D processes are standardized at baseline. These two constraints are compensated by parameters in the measurement models.

The temporal influences between processes are modelled through the D × D-matrix of time-dependent transition intensities A i,δ (t):

A i,δ (t) =          a i,
        
This matrix captures the directed temporal influences between latent processes at time t and subsequent rates of change of latent processes between times t and t + δ. Specifically, coefficient a i,dd (t) quantifies the temporal effect of process d at time t on process d. Each effect can be modelled according to time/covariates through a linear regression a i,dd (t) = R i (t)α dd where R i (t) is a r-vector of time-dependent covariates associated with the r-vector of regression coefficients α dd = (α m dd ) m=0,(r-1) . When the discretization step is not too large, the temporal influences intend to have the same causal interpretations as those of a model in continuous time (see Simulation Study 2 for the demonstration).

Measurement Models of the longitudinal markers

Consider K (K ≥ D) continuous longitudinal markers Y ij = (Y ijk ) k=1,...,K that have been measured for subject i at (n i + 1) successive discrete times t ij = j × δ with j ∈ τ i = {j 0i , . . . , (j 0i + n i )}. Following [START_REF] Proust-Lima | A nonlinear model with latent process for cognitive evolution using multivariate longitudinal data[END_REF], 2013), we assume that each latent process Λ d i is the underlying common factor of K d markers (K = D d=1 K d ) and we note K d the set of marker subscripts associated with latent process d. We assume that a marker measures only one latent process.

The link between a marker and its underlying latent process is defined by a markerspecific measurement model. If marker k is Gaussian, the measurement model is a linear equation:

Y ijk -η 0k η 1k = Ỹijk = Λ d i (t ij ) + ˜ ijk , ∀k ∈ K d & ∀j ∈ τ i (2)
where the vector of transformation parameters η k = (η 0k , η 1k ) is used to get the standardized form Ỹijk of the marker and ˜ ijk are independent Gaussian errors with variance σ 2 k . In the more general case of a continuous marker (possibly non Gaussian), one may consider a nonlinear equation of observation:

H k (Y ijk ; η k ) = Ỹijk = Λ d i (t ij ) + ˜ ijk , ∀k ∈ K d & ∀j ∈ τ i (3)
where the link transformation H k comes from a family of monotonic increasing and continuous functions parameterized with η k . Again Ỹijk is the transformed marker and ˜ ijk are independent Gaussian errors with variance σ 2 k . Following Proust-Lima et al. ( 2013), the link transformation H k can be defined from a basis of I-splines (which are integrated M-splines [START_REF] Ramsay | Monotone regression splines in action[END_REF]) in association with positive coefficients, thus providing an increasing bijective flexible transformation. We used here a quadratic I-splines basis with p k internal knots, (I m ) m=1,p k +3 , so that

H k (Y ijk ; η k ) = Ỹijk = η 0k + p k +3 m=1 η 2 mk I m (Y ijk ),
with (η mk ) m=0,p k +3 the vector of parameters of the transformation.

In the following, we denote Σ = diag((σ 2 k ) k=1,...,K ) the diagonal variance matrix of the vector of errors ˜ ij = (˜ ijk ) k=1,...,K and η = (η k ) {1,...,K} , the total vector of transformation parameters for the K markers. The vector of transformed markers Ỹij is mapped to the system of latent processes Λ i (t ij ) through a K × D matrix P with element (k, d) equal to 1 if marker k measures latent process d:

Ỹij = P Λ i (t ij ) + ij (4)
In practice the observation process may include intermittent missing observations for a subset of markers or for all the markers at a given occasion j ∈ τ i , so that K * ij ≤ K markers are actually observed at occasion j for subject i. We assume that observations are missing at random and note Ỹ * ij the transformations of the actual K * ij -vector of observed markers Y * ij at occasion j. The model linking the system of processes Λ i (t ij ) to the transformed markers Ỹ * ij can be easily adapted to the presence of intermittent missing data by considering a K * ij × K observation matrix M ij (for j ∈ τ i ) where element (k * , k) equals 1 if marker k is the k * th observed marker at occasion j and 0 if not for k = 1, ..., K and k * = 1, ..., K * ij :

Ỹ * ij = M ij P Λ i (t ij ) + * ij (5)
with * ij the vector of independent Gaussian errors with covariance matrix M ij ΣM ij .

Estimation

The parameters are estimated in the maximum likelihood framework.

Distribution of the Latent Processes and transformed observations

We first derive the marginal distributions of the latent processes and of the transformed observations. By recurrence, the structural model ( 1) can be rewritten:

Λ i (t ij ) =                  X 0 i β + u i if j = 0 j-1 l=0 Ãi,δ (t il )(X 0 i β + u i ) if j > 0 +δ j s=1 j-1 l=s Ãi,δ (t il )(X i (t is )γ + Z i (t is )v i ) (6)
where t ij = j × δ for j ∈ {0, ..., J} and Ãi,δ

(t ij ) = I D + δA i,δ (t ij ).
By introducing Ψ i,δ (t 0 , j, s) for t 0 ≥ 0 and s ≤ j so that

Ψ i,δ (t 0 , j, s) =      I D , if s = j j-1 l=s Ãi,δ (t 0 + t il ) if s < j (7)
Equation ( 6) can be rewritten

Λ i (t ij ) = Ψ i,δ (0, j, 0)(X 0 i β +u i )+ δ j s=1 Ψ i,δ (0, j, s) (X i (t is )γ + Z i (t is )v i ) 1 {j>0} (8)
As a result, the vector Λ i (t ij ) has a multivariate normal distribution with expectation µ Λ ij and variance covariance matrix V Λ ijj = var (Λ i (t ij )) and the vector

Λ i = Λ i (t ij ) j∈τ i has a multivariate normal distribution with expectation µ Λ i = µ Λ ij j∈τ i and variance-covariance matrix V Λ i = V Λ ijj (j,j )∈τ 2 i where µ Λ ij = E[Λ i (t ij )] = Ψ i,δ (0, j, 0)X 0 i β + δ j s=1 Ψ i,δ (0, j, s)X i (t is )γ 1 {j>0} (9)
and

V Λ ijj = cov(Λ i (t ij ); Λ i (t ij )) = Ψ i,δ (0, j, 0)Ψ i,δ (0, j , 0) +   Ψ i,δ (0, j, 0)B uv δ j s =1 Ψ i,δ (0, j , s )Z i (t is )   1 {j >0} + δ j s=1 Ψ i,δ (0, j, s)Z i (t is ) B uv Ψ i,δ (0, j , 0) 1 {j>0} +   δ j s=1 Ψ i,δ (0, j, s)Z i (t is ) B v δ j s =1 Ψ i,δ (0, j , s )Z i (t is )   1 {min(j,j )>0}
(10) It can be easily deduced that the vector of incomplete and transformed data

Ỹ * ij at occasion j is multivariate Gaussian with expectation µ Ỹ * ij = M ij P µ Λ ij and variance- covariance matrix V Ỹ * ij = M ij P V Λ ijj P + Σ M ij ,

and the total vector of incomplete and transformed data Ỹ

* i = Ỹ * ij j∈τ i is multivariate Gaussian with expectation µ Ỹ * i = µ Ỹ * ij j∈τ i and variance-covariance matrix V Ỹ * i , a block matrix with M ij P V Λ ijj P M ij + M ij ΣM ij 1 {j=j } the (j, j ) block.

Likelihood

As the N subjects of the sample are independent, the log-likelihood of the model is

L(Y * ; θ) = N i=1 log(L i (Y * i ; θ)) with L i (Y * i ; θ) the individual contribution to the likeli- hood. Here, θ = (β , γ , vec(B) , (α dd ) d,d ∈{1,...,D} 2 , (σ k ) k∈{1,.
..,K} , η ) is the whole vector of parameters. Using the Jacobian of the link functions H k (k = 1, ..., K), the individual contribution is:

L i (Y * i ; θ) = φ i ( Ỹ * i ; µ Ỹ * i , V Ỹ * i ) j∈τ i K * ij l=1 J κ(l) H κ(l) ( Y * ijκ(l) ; η κ(l) ) (11)
where φ i (.; µ, V ) denotes the density function of a multivariate Gaussian vector with expectation µ and variance-covariance matrix V , and J κ(l) (H κ(l) ( Ỹ * ijκ(l) ; η κ(l) )) denotes the Jacobian of the link function H κ(l) used to transform Y * ijκ(l) , the l th observed marker at oc-casion j for subject i. For instance, with a linear link function,

J κ(l) H κ(l) ( Ỹ * ijκ(l) ; η κ(l) ) = 1 η 1κ(l) .

Optimisation Algorithm and Implementation

The maximum likelihood estimates are obtained using an extended Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] because of its robustness and good convergence rate. At each iteration p, if necessary, the Hessian matrix H (p) is diagonal-inflated to obtain a positive definite matrix H * (p) which is used to update the parameters

θ (p+1) = θ (p) - ν[H * (p) ] -1 U (θ (p) ), with U (θ (p)
) the gradient at iteration p and ν the improvement control parameter. Convergence is reached when

||θ (p+1) -θ (p) || 2 < 10 -3 , |L(Y * ; θ (p+1) ) - L(Y * ; θ (p) )| < 10 -3 and U (θ (p) ) [H (p) ] -1 U (θ (p) )
npara < 10 -2 , with n para the total number of parameters. The latter criterion is by far the most stringent one and specifically targets maximum search. The variances of the estimators are obtained from the inverse of H (p) .

Given the possibly high number of parameters, we first estimate the parameters for each process taken separately, then we start the maximization of the likelihood of the multivariate model from these simple estimates, setting initial values of the inter-dimension parameters to zero. The model estimation is implemented in R (program available on request); it combines R and C ++ languages, and includes parallel computations.

Marginal and subject-specific predictions

The goodness-of-fit of the model can be assessed by comparing predictions with observations of the markers in their transformed scales. From notations defined in section 3.3.1, marginal and conditional distributions of the markers are:

Ỹi ∼ N P µ Λ i , P V Λ i P + Σ i . ( 12
)
Ỹi Λ i ∼ N (P Λ i , Σ i ) (13)
where Σ i is the block-diagonal matrix constituted of n i Σ blocks.

The marginal (Y (M ) i

) and subject-specific (Y

(SS) i
) predictions in the transformed scales are then respectively obtained by taking the expectations of the marginal and conditional distributions of the transformed markers at the parameter estimates θ and at the predicted latent processes Λi of Λ i given the observations

Ỹ * i : Λi = E Λ i | Ỹ * i = µ Λ i +C Λ i Ỹ * i V -1 Ỹ * i Ỹ * i -µ Ỹ * i , where C Λ i Ỹ * i = V Λ ijj P M ij (j,j )∈τ i 2
is the covariance matrix between Λ i and Ỹ * i . Using these individual predictions, one can graphically compare either the marginal predictions Y (M ) or subject-specific predictions Y (SS) averaged within time intervals to the observations averaged within the same time intervals. Marginal and subject-specific predictions in the natural scale of the markers can also be derived from the marginal and conditional distributions by using a Monte-Carlo approximation [START_REF] Proust-Lima | Analysis of multivariate mixed longitudinal data: a flexible latent process approach[END_REF].

Simulations

We performed two series of simulations to evaluate the estimation program and the impact of time discretization on the interpretation of A(t) matrix.

Simulation Study 1: Validation of the estimation program

Design

To evaluate the estimation program, we generated a system of two Gaussian processes ((Λ 1 (t)) t≥0 and (Λ 2 (t)) t≥0 ), each one measured by one longitudinal marker (Y 1 and Y 2 ) and two covariates, one binary C 1 and one continuous C 2 . We considered two scenarios: a covariate-specific temporal influence structure (Scenario 1) and a time-dependent temporal influence structure (Scenario 2). In Scenario 1, we assumed linear trajectories over time for the system of latent processes (random intercepts and simple effects of both covariates in the sub-models for the initial level and the change over time) and a temporal influence structure

A(t) different in each level of C 1 :            Λ 1 i (0) = β 1 0 + β 1 1 C 1,i + β 1 2 C 2,i + u 1 i Λ 2 i (0) = β 2 0 + β 2 1 C 1,i + β 2 2 C 2,i + u 2 i ∆Λ 1 i (t+δ) δ = γ 1 0 + γ 1 1 C 1,i + γ 1 2 C 2,i + v 1 i + (α 0 i,11 + α 1 i,11 C 1,i )Λ 1 i (t) + (α 0 i,12 + α 1 i,12 C 1,i )Λ 2 i (t) ∆Λ 2 i (t+δ) δ = γ 2 0 + γ 2 1 C 1,i + γ 2 2 C 2,i + v 2 i + (α 0 i,21 + α 1 i,21 C 1,i )Λ 1 i (t) + (α 0 i,22 + α 1 i,22 C 1,i )Λ 2 i (t) Y ijk -η 0k η 1k = Λ k ij + ˜ ijk , k = 1, 2 (14) where u i = (u 1 i , u 2 i ) , v i = (v 1 i , v 2 i
) and (u i , v i ) ∼ N (0, LL ) with L such that the random effects are independent between dimensions, and ˜ ijk ∼ N (0, σ 2 k ), ∀k ∈ {1, 2}.

In scenario 2, we assumed linear trajectories for the latent processes adjusted for C 2 and a transition matrix that evolved with time:

           Λ 1 i (0) = β 1 0 + β 1 1 C 2,i + u 1 i Λ 2 i (0) = β 2 0 + β 2 1 C 2,i + u 2 i ∆Λ 1 i (t+δ) δ = γ 1 0 + v 1 i + a 11 (t)Λ 1 i (t) + a 12 (t)Λ 2 i (t) ∆Λ 2 i (t+δ) δ = γ 2 0 + v 2 i + a 21 (t)Λ 1 i (t) + a 22 (t)Λ 2 i (t) Y ijk -η 0k η 1k = Λ k ij + ˜ ijk , k = 1, 2, ( 15
)
where

u i = (u 1 i , u 2 i ) , v i = (v 1 i , v 2 i ) and (u i , v i ) ∼ N (0, LL
) with L such that the random effects are independent between dimensions, and ˜ ijk ∼ N (0, σ 2 k ), ∀k ∈ {1, 2}. Each element of the transition matrix a kk (t) is defined from a basis of quadratic B-splines with one internal knot at the median (S m ) {m=1,3} so that a kk (t

) = α 0 kk + α 1 kk S 1 (t) + α 2 kk S 2 (t) + α 3 kk S 3 (t), ∀k, k ∈ {1, 2} 2 .
The design of the simulations and the parameters were chosen to mimic the ADNI data. Dimensions 1 and 2 were cerebral anatomy and cognitive ability, respectively measured by a specific composite score. Covariate C 1 corresponded to the two groups CN (for normal elderly, in reference) and MCI (for Mild Cognitive Impairment). It was generated according to a Bernoulli distribution with probability P (C 1 = 1) = 0.64. Covariate C 2 corresponded to the standardized age at baseline and was generated according to a standard Gaussian distribution. We considered a discretization step δ of 6 months. Markers observations were generated every 6 months up to 3 years. Thus a subject had 7 repeated measures at occasions j ∈ {0, 1, . . . , 6}. We also considered a design in which scheduled visits could be missed completely at random with a probability of 0.15, and when a visit was not missed, a marker could be missing with a probability of 0.07. For each design and scenario, we generated 100 samples of 512 subjects.

Results

Table 4 and Table 1 provide the results of the simulations for scenarios 1 and 2 respectively. In both settings, all the parameters were correctly estimated with satisfying coverage rates in the absence of missing data (left part) and in the presence of missing data (right part).

Design

To formally assess whether the interpretation under the discretized time was the same as the one obtained under continuous time, we assessed the type-I error rate associated with each non diagonal element of the transition matrix A under three discretization steps δ = 1/3, 1/2 and 1 when data were actually generated under continuous time (approximated by a step δ=0.001). We considered for this a system of three latent processes ((Λ 1 (t)) t≥0 , (Λ 2 (t)) t≥0 , and (Λ 3 (t)) t≥0 ), each one measured by one Gaussian repeated marker (Y 1 , Y 2 , Y 3 ). The trajectories were linear, with no adjustment for covariates and a transition matrix A constant over time:

                     Λ 1 i (0) = β 1 0 + u 1 i Λ 2 i (0) = β 2 0 + u 2 i Λ 3 i (0) = β 3 0 + u 3 i ∆Λ 1 i (t+δ) δ = γ 1 0 + v 1 i + a 11 Λ 1 i (t) + a 12 Λ 2 i (t) + a 13 Λ 3 i (t) ∆Λ 2 i (t+δ) δ = γ 2 0 + v 2 i + a 21 Λ 1 i (t) + a 22 Λ 2 i (t) + a 23 Λ 3 i (t) ∆Λ 3 i (t+δ) δ = γ 3 0 + v 3 i + a 31 Λ 1 i (t) + a 32 Λ 2 i (t) + a 33 Λ 3 i (t) Y ijk -η 0k η 1k = Λ k ij + ˜ ijk , k = 1, 2, 3 (16)
The simulation model mimicked the ADNI data with cerebral anatomy, cognitive ability and functional autonomy as dimensions, respectively measured by a specific composite score. We estimated the model on the ADNI data with δ=1 (6 months) and transformed the parameters to the scale δ=0.001 to generate the data in almost continuous time. We provide in the Appendix A: the formulas that we used to relate model components defined under δ=1 and δ=0.001. Elements of the transition matrix were set one by one to 0 in δ = 0.001 scale to evaluate the associated type-I error rate. Latent processes were generated with solver from dsolve package [START_REF] Soetaert | Solving differential equations in r: Package desolve[END_REF] and observations were derived every 6 months up to 3 years. For each design, we considered 200 samples of 512 subjects.

Results

Table 2 displays the type-I error rates (in percentage) associated with each non diagonal element of the transition matrix A when estimated with discretization steps δ=1/3, δ=1/2 and δ=1. All the type-I error rates are close to the nominal 5% rate (expected rates in the 95% interval [1.98 , 8.02] with 200 replicates). This shows that with such discretization steps, the causal interpretation expected in continuous time is not altered.

Table 2: type-I error rates (in %) associated with each non diagonal element of the transition matrix A when the transition matrix is generated approximately in continuous time (δ=0.001) and estimated with discretization steps: δ=1/3, δ=1/2, δ=1 (expected 95% interval [1.98 , 8.02] for the nominal type-I error of 5%).

Parameter δ = 1/3 δ = 1/2 δ = 1 a 12 5.5 4.5 6.5 a 13

3.5 5.5 5.0 a 21 7.0 7.5 8.0 a 23 5.0 4.5 5.0 a 31 4.0 3.5 5.5 a 32 7.0 7.5 5.5

Application

The application aimed to describe the decline over time of cerebral anatomy, cognitive ability and functional autonomy in three clinical stages of AD (normal aging (CN), Mild Cognitive Impairment (MCI) and diagnosed with Alzheimer's Disease (dAD)) and to quantify the temporal influences between these dimensions by assessing especially whether the relationships differed according to the clinical stage.

The sample

We considered the nineteen available cognitive tests as outcomes of global cognitive ability [START_REF] Crane | Development and assessment of a composite score for memory in the alzheimer's disease neuroimaging initiative (ADNI)[END_REF], the volume of hippocampus relative to total intracranial volume and the cortical thickness of nine regions (Freesurfer version 4.4.0 for longitudinal data) as outcomes for cerebral anatomy following regions of interest identified by [START_REF] Dickerson | The cortical signature of Alzheimer's disease: regionally specific cortical thinning relates to symptom severity in very mild to mild ad dementia and is detectable in asymptomatic amyloidpositive individuals[END_REF], and the sumscore of the FAQ composed of 30 items [START_REF] Pfeffer | Measurement of functional activities in older adults in the community[END_REF] to assess functional autonomy. For numerical simplicity, we considered a unique marker for each dimension by summarizing the cerebral and cognitive markers into two composite scores [START_REF] Proust-Lima | Are latent variable models preferable to composite score approaches when assessing risk factors of change? Evaluation of type-I error and statistical power in longitudinal cognitive studies[END_REF]

The dynamic multivariate model

The dynamic model applied to ADNI is summarized in Figure 1. The three latent processes (cerebral anatomy (Λ A ), cognitive ability (Λ C ), functional autonomy (Λ F )) were repeatedly measured by the three composite scores Y A , Y C and Y F . We assumed that changes of the processes were linear in time; thus, the processes had quadratic trajectories over time, adjusted for gender, education, APOE genotype, clinical stage, and age at baseline. In order to account for the correlation between individual repeated measures, we included random intercepts on the initial level of the processes and random intercepts and slopes on the change of the processes over time. Note that random effects on initial levels and on changes are correlated within each dimension but not between dimensions. The transition matrix that captured the interrelations between processes was constant over time and adjusted for clinical stage. To correct the possible departure from normality of the composite scores, we transformed them using a link function approximated by quadratic I-splines with 2 internal knots chosen at the terciles for cerebral and cognitive composite scores and by quadratic I-splines with one internal knot chosen at the median for the functional composite score.

In the main analysis, we took a 6 months discretization step. The strategy of analysis consisted in finding first the best adjustment for each process taken separately (assuming independency between processes ie., a dd {d =d } = 0 in Figure 1) with a significance threshold for covariate effects at 25%. Then, the whole multivariate model including all elements of the transition matrix was estimated. In secondary analyses, we re-estimated the final model by considering discretization steps of 3 months and of 2 months in order to evaluate whether the interpretations varied with smaller steps.

Results

Latent process-specific trajectories

Estimates of fixed effects, Cholesky's decomposition parameters (for the random effects covariance matrix) and measurement model parameters are provided in Tables 5, 6 and7. In summary, older age, male gender, APOE 4 carrying, and clinical stages MCI and dAD were associated with lower cerebral anatomy and lower cognitive ability levels at baseline. Higher education was associated with lower cerebral anatomy level but associated with higher cognitive ability level at baseline. Only clinical stages MCI and dAD were associated with lower functional autonomy level at baseline. APOE 4 carrying was associated with steeper declines in cerebral anatomy and cognitive ability. Clinical stage dAD was associated with steeper cerebral anatomy decline and both clinical stages MCI and dAD were associated with steeper declines in cognitive ability and functional autonomy. Higher education was associated with smaller cognitive ability decline. Figure 4 depicts the expected trajectories of each dimension according to stage for two profiles of individuals (women noncarrier of APOE 4 and with lower educational level; men carrier of APOE 4 and with higher educational level).

Temporal influence structure between processes

Estimates of the transition intensity matrix are given in Table 3 andsummarized in Figure 2. Figure 2 shows that the temporal influences between cerebral anatomy (Λ A ), cognitive ability (Λ C ), and functional autonomy (Λ F ) slightly evolve from healthy stage to Alzheimer's disease stage. At all stages, cerebral dimension affects significantly the dynamics of both cognitive and functional dimensions. Moreover, the cognitive and functional dimensions are interrelated, each one affecting the change over time of the other. However from MCI stage, the cognitive dimension becomes central by affecting significantly both the dynamics of cerebral and functional dimensions. When considering a discretization step of 3 or 2 months instead of 6 months, the Akaike Information criterion was slightly improved (AIC=1386.68 with 6 months, AIC= 1370.45 with 3 months and AIC=1367.21 with 2 months) but the results regarding the temporal influence relationships between dimensions remained the same (see Table 8). mean transformed observations at each visit (every 6 months). As shown in Figure 3, the mean predictions are very close to the mean observations showing the good fit of the model to the data. These trajectories also illustrate the differences in trajectories according to clinical stage with each dimension more and more impaired and decline steeper and steeper with disease progression.

Discussion

We proposed an original dynamic model which aims to simultaneously describe multivariate processes over time and retrieve temporal relationships between the processes involved. Our model aims to be a dynamic causal model in the spirit of mechanistic models. We relied on discrete time with difference equations rather than continuous time with differential equations to largely reduce the numerical problems, notably with a closed-form likelihood.

We fundamentally explored in this work temporal influences between dynamic processes and the causal interpretation of these influences has to be done with caution. We relied on the dynamic approach to causality notably developed by [START_REF] Aalen | What can statistics contribute to a causal understanding[END_REF], [START_REF] Arjas | Causal reasoning from longitudinal data[END_REF], [START_REF] Didelez | Graphical models for marked point processes based on local independence[END_REF], [START_REF] Commenges | A General Dynamical Statistical Model with Causal Interpretation[END_REF]. There is also a large literature in statistics on the potential outcomes formalism introduced by [START_REF] Rubin | Causal inference using potential outcomes[END_REF]; it is however not adapted for modeling dynamic relationships between non-directly observed factors, which is the aim in this paper so we do not discuss this approach further.

Following the dynamic approach to causality, we chose to analyze the changes of processes (as seeking local dependence structures) and defined causal relationships in a structural model, that is at a latent process level rather than at an observation level. One major simplifying assumption in our approach was the time discretization. Fundamentally, causal relationships are to be explored at an infinitesimal level and thus, a causal model is to be defined in continuous time [START_REF] Commenges | A General Dynamical Statistical Model with Causal Interpretation[END_REF][START_REF] Aalen | Can we believe the dags? a comment on the relationship between causal dags and mechanisms[END_REF]. To assess whether we could claim the same type of interpretation as mechanistic models, we thus assessed in a simulation study the impact of the discretization of the underlying time-continuous causal structure; we found that the type-I error rates of the temporal influence parameters were not altered by the use of a model in discrete time, provided the discretization step remained small in regards with the dynamics of the disease under study. Indeed in our simulations and application on Alzheimer's disease, our discretization step was between 2 and 6 months while the disease progresses over decades [START_REF] Amieva | Prodromal Alzheimer's disease: successive emergence of the clinical symptoms[END_REF][START_REF] Jack | Tracking pathophysiological processes in Alzheimer's disease: an updated hypothetical model of dynamic biomarkers[END_REF]. In addition to the assessment that the absence of causal associations could be retrieved in the discretized time, we also provided in supplementary material analytic approximate relationships between transition matrices defined in two discretization steps or in continuous time and discrete time. These relationships quantify the associations between temporal influence structures; they remain however approximate as they were found under the assumption that the model structure did not differ substantially between the two discretization steps. This assumption deserves more attention in future work.

We acknowledge that despite our intention to provide a dynamic causal model, causal interpretation still has too be made cautiously as in any causal model since it is always subject to a correct specification of the system and of the statistical model.

The estimation we proposed relies on the assumption that missing data, both intermittent and monotone, are missing at random. This assumption is realistic in ADNI 1 which had a very short follow-up (3 years max). However, in applications with a longer followup, and notably in population-based cohorts, such assumption might become restrictive and it will be probably necessary to account for the occurrence of informative clinical events, such as dementia and death in our context. This will be possible by extending our model to a joint analysis of times to event but probably at the price of additional numerical complexities. In addition, in studies on chronic diseases with a long follow-up, it is reasonable to assume that the temporal influence structure may evolve with time. Our model can handle this as shown in the simulations where the temporal influence structure was modelled according to time using regression splines. We did not consider however time-dependent temporal influence structures in the application due to the short follow up. In addition, our objective was to contrast the influence structures between clinical stages, and as such we rather included stage-specific intensities of transition.

The model was primarily motivated by the study of the multiple alterations involved in the dementia process in the elderly. Indeed, dementia is characterized by very long and progressive alterations in different cerebral and clinical dimensions. However, although there is a global agreement on the relevant dimensions in AD, their relationships are still poorly understood and a confusion persists between alterations due to normal aging and alterations due to pathological aging leading to dementia. Thanks to the ADNI data which included subjects at different clinical stages of dementia, we were able to exhibit the global structure of dependence between cerebral, cognitive and functional domains in normal aging characterized as expected by the cerebral structure which explains a part of the cognitive and functional changes, and the functional and cognitive changes which have bidirectional relationships. We also found out that the main evolution in this structure due to the pathological process towards dementia (in subjects with mild cognitive impairment or diagnosed with AD) was on cognitive functioning with some effect of cognitive functioning on change in cerebral structure in addition to its effect on change in functional structure. Such backward influence of cognition on anatomic structure may appear counterintuitive from a biological point of view. Yet, some studies showed that interventions made on cognitive functioning (e.g., literacy at adulthood) could induce an improvement in cerebral structure [START_REF] Carreiras | An anatomical signature for literacy[END_REF][START_REF] Boyke | Training-induced brain structure changes in the elderly[END_REF]. This application which gives a first insight on the possibilities of this dynamic model, could be now refined by targeting specific brain regions (separating for instance cortical thicknesses from hippocampal volume) and specific cognitive functions (separating for instance episodic memory from executive functioning).

As a conclusion, we proposed here a new methodology that may help identify temporal structures in multivariate longitudinal data. Although motivated and applied in dementia context, this approach has some potential to address unsolved questions in many other chronic diseases where multiple processes are in play, and or time-dependent exposures are to be assessed.

Appendix A: Relationships between model components under different discretization steps

For computational reasons, the latent processes in the structural model are assumed to evolve in discrete time with a constant step δ which may vary depending on the application. In reality the processes evolve in continuous time and causal inference should be done in continuous time. We thus established the relationship between transition matrices obtained under different discretization steps: δ * and a smaller discretization step δ with δ * = ρ × δ and integer ρ > 1.

With a discretization step δ = δ * ρ , the second line of equation ( 1) can be rewritten:

Λ i (t + δ) = δ (X i (t + δ)γ + Z i (t + δ)v i ) + (I D + δA i,δ (t)) Λ i (t), ∀t > 0, ρ > 1 (17)
By recurrence, the network level at t + ρδ, that also correspond to t + δ * , can be expressed as a function of Λ i (t + δ * ) as follows:

Λ i (t + δ * ) = Ψ i,δ (t, ρ, 0)Λ i (t) + δ ρ s=1 Ψ i,δ (t, ρ, s) (X i (t + δs)γ + Z i (t + δs)v i ) (18)
with function Ψ i,δ defined in the relation (7).

Considering directly a discretization step of δ * , the second line of equation (1) can also be rewritten

Λ i (t+δ * ) = δ * (X i (t + δ * )γ * + Z i (t + δ * )v * i )+(I D + δ * A i,δ * (t)) Λ i (t), ∀t > 0 and δ * > 0. ( 19
)
Relationships between transition matrices under different steps Considering that the second parts of equations ( 18) and ( 19) are close enough although the model specifications are different, we obtain:

I D + δ * A i,δ * (t) ≈ Ψ i,δ (t, ρ, 0) = ρ-1 l=0 I D + δ * ρ A i,δ (t + l δ * ρ ) (20)
Equation ( 20) can be rewritten to highlight the relationship between the causal structure A i,δ * and A i,δ :

A i,δ * (t) = 1 δ * ρ-1 l=0 I D + δ * ρ A i,δ (t + l δ * ρ ) -I D (21)
When assuming that the causal structure A i,δ is constant in each interval [t, t + δ * ], ∀t ∈ τ , the relationship (21) becomes :

A i,δ * (t) = 1 δ * I D + δ * ρ A i,δ (t) ρ -I D (22)
Or equivalently,

A i,δ (t) = ρ (δ * ) (ρ+1) 1 δ * I D + A i,δ * (t) 1 ρ -I D (23)
Note that the limits of the right part of relation ( 22), when ρ tends to infinity provide the relationship between A i,δ * (t) and its continuous time counterpart. Assuming that lim ρ→∞ A i,δ (t) = A i (t) is finite, the relationship (22), become:

A i,δ * (t) ≈ 1 δ * (exp (δ * A i (t)) -I D ) (24)
Relationships between trend parameters under different steps To relate the transition matrices under different steps, we considered second parts of equations ( 18) and ( 19) were close enough. We identified:

           X i (t + δ * )γ * ≈ 1 ρ ρ s=1 Ψ i,δ (t, ρ, s)X i (t + δ * ρ s)γ Z i (t + δ * )v * i ≈ 1 ρ ρ s=1 Ψ i,δ (t, ρ, s)Z i (t + δ * ρ s)v i (25)
In the particular case where X i and Z i reduce to intercepts,

           γ * ≈ 1 ρ ρ s=1 Ψ i,δ (t, ρ, s)γ v * i ≈ 1 ρ ρ s=1 Ψ i,δ (t, ρ, s)v i (26)
Or equivalently,

           γ ≈ ρ ρ s=1 Ψ i,δ (t, ρ, s) -1 γ * v i ≈ ρ ρ s=1 Ψ i,δ (t, ρ, s) -1 v * i ( 27
)
Appendix B: Additional results in the simulations studies

Figure 1 :

 1 Figure 1: Graph of the dynamic causal model considered on ADNI data with three dimensions (labelled Λ A for cerebral anatomy, Λ C for cognitive ability and Λ F for functional autonomy), each one repeatedly measured by one marker Y A , Y C and Y F .

Figure 2 :

 2 Figure 2: Temporal relationships estimated between cerebral anatomy (Λ A ), cognitive ability (Λ C ), and functional autonomy (Λ F ) at healthy (CN), MCI and dAD stages. Arrows represent effects of one dimension on the change of another dimension between two discretized times. Numbers indicate intensity effects, stars indicate the significance of the effects. Arrows in solid line indicate effects that exist at all stages; arrows in dashed line indicate effects that exist only for some stages.

Figure 3 :

 3 Figure 3: Means of subject-specific predictions at each visit (with crosses) along with the corresponding means of transformed observed scores (plain lines) and their 95% confidence interval (dashed lines). Columns refer to the group (controls CN, Subjects with Mild Cognitive Impairments MCI and subjects diagnosed with Alzheimer's disease dAD). Rows refer to the dimensions (cerebral anatomy, cognitive ability and functional autonomy).

Figure 4 :

 4 Figure 4: Predicted trajectories of cerebral, cognitive and functional dimensions at healthy (CN), MCI, and dAD clinical stages of Alzheimer's disease for a woman non-carrier of APOE 4 allele with lower educational level (plain lines) and for a man carrier of APOE 4 allele with higher educational level (dashed lines). "n", "m" and "d" indicate trajectories at healthy, MCI and dAD stages, respectively.

Table 1 :

 1 Results of the simulations for scenario 2 (100 replicates of samples of size 512).

			without missing values (conv † =100%)	with missing values (conv † =98%)
		θ	θ ESE ‡ ASE ‡	CR(%)	θ ESE ‡ ASE ‡	CR(%)
	β	-1.661 -1.655 0.104 0.112	96.0	-1.653 0.111 0.113	98.0
	β	-1.803 -1.831 0.146 0.128	92.0	-1.834 0.150 0.133	90.0
	γ	0.278	0.281 0.043 0.046	96.0	0.285 0.045 0.050	98.0
	γ	-0.152 -0.154 0.023 0.023	96.0	-0.156 0.024 0.026	96.0
	γ	0.117	0.121 0.080 0.080	94.0	0.111 0.093 0.086	91.0
	γ	-0.143 -0.145 0.040 0.041	96.0	-0.140 0.046 0.045	95.0
	L(3,1)	0.041	0.040 0.019 0.022	97.0	0.039 0.020 0.024	98.9
	L(4,2)	0.012	0.026 0.021 0.039	97.0	0.029 0.025 0.044	97.0
	L(3,3)	0.198	0.197 0.015 0.014	94.0	0.196 0.016 0.014	92.0
	L(4,4)	0.373	0.372 0.027 0.027	95.0	0.371 0.029 0.029	93.0
	α 1	0.112	0.109 0.040 0.044	96.0	0.107 0.046 0.048	97.0
	α 2	0.003	0.010 0.047 0.050	95.0	0.015 0.052 0.056	96.0
	α 3	-0.033 -0.050 0.070 0.070	97.0	-0.056 0.070 0.079	98.0
	α 4	-0.132 -0.117 0.095 0.093	95.0	-0.108 0.099 0.104	97.0
	α 1	0.168	0.170 0.039 0.037	93.0	0.172 0.044 0.041	92.0
	α 2	0.079	0.067 0.043 0.043	92.0	0.067 0.049 0.048	92.0
	α 3	0.143	0.160 0.062 0.057	92.0	0.161 0.065 0.064	92.0
	α 4	0.070	0.050 0.080 0.073	90.0	0.048 0.088 0.081	89.0
	α 1	0.238	0.225 0.071 0.068	95.0	0.220 0.076 0.073	92.0
	α 2	0.138	0.159 0.092 0.090	94.0	0.162 0.097 0.097	93.0
	α 3	0.359	0.331 0.128 0.126	94.0	0.330 0.134 0.135	95.0
	α 4	0.014	0.048 0.170 0.165	94.0	0.045 0.179 0.178	95.0
	α 1	0.108	0.108 0.077 0.074	91.0	0.102 0.085 0.081	89.0
	α 2	0.037	0.024 0.082 0.084	97.0	0.025 0.091 0.088	95.0
	α 3	-0.067 -0.046 0.105 0.107	90.0	-0.049 0.116 0.113	91.0
	α 4	-0.066 -0.093 0.133 0.133	96.0	-0.085 0.154 0.142	94.0
	η	0.387	0.391 0.015 0.015	93.0	0.391 0.015 0.016	96.0
	η	0.713	0.713 0.035 0.033	94.0	0.709 0.035 0.035	95.0
	η	3.769	3.725 0.192 0.197	95.0	3.719 0.199 0.199	93.0
	η	2.580	2.561 0.092 0.090	95.0	2.563 0.094 0.092	95.0
	σ	2.489	2.511 0.138 0.120	92.0	2.515 0.144 0.124	91.0
	σ	1.412	1.413 0.063 0.060	93.0	1.415 0.063 0.064	95.0

  . Covariates were the age at entry (centered around 75.4 and indicated in decades), gender, educational level (low level if ≤ 12 years vs high level if >12 years), the APOE genotype ( 4 carrier vs 4 non-carrier ) and the 3 clinical stages (CN, MCI, dAD). We selected in the sample all the subjects who had no missing data for the covariates and had at least one measure for each dimension during the follow-up. The main analysis included 656 subjects (82% of the initial sample).

	The sample consisted of 190
	(29%) subjects at healthy stage (CN), 322 (49%) subjects at mild cognitive impairment
	stage (MCI), and 144 (22%) subjects diagnosed with Alzheimer's disease (dAD); 43%
	were females, 83% had a high educational level and 49% carried APOE 4 allele. The
	mean age at entry was 75.4 years old (sd=6.6). The mean number of visits was 5, 6 and
	4 for CN, MCI and dAD subjects, respectively.

Table 4 :

 4 Results of the simulations (100 replicates of samples with 512 subjects) considering two latent processes, each one repeatedly measured by a marker, with linear trajectories and constant covariate-specific causal structure. ASE is the asymptotic standard error, ESE is the empirical standard error and CR is the coverage rate of the 95% confidence interval.

			without missing values (conv † =100%)	with missing values (conv † =99%)
		θ	θ ESE ‡ ASE ‡	CR(%)	θ ESE ‡ ASE ‡	CR(%)
	β 1 1	-0.268 -0.267 0.068 0.073	97.0	-0.269 0.070 0.075	97.9
	β 1 2	-1.695 -1.698 0.112 0.112	92.0	-1.705 0.113 0.116	92.9
	β 2 1	0.057	0.055 0.075 0.078	96.0	0.060 0.079 0.082	95.9
	β 2 2	-1.749 -1.777 0.139 0.127	94.0	-1.775 0.148 0.137	94.9
	γ 1 0	0.042	0.044 0.018 0.019	97.0	0.044 0.020 0.021	97.9
	γ 1 1	-0.033 -0.032 0.020 0.019	94.0	-0.032 0.020 0.020	94.9
	γ 1 2	-0.242 -0.232 0.048 0.052	97.0	-0.225 0.054 0.063	98.9
	γ 2 0	-0.097 -0.098 0.041 0.040	95.0	-0.098 0.047 0.045	92.9
	γ 2 1	-0.014 -0.016 0.023 0.020	93.0	-0.015 0.023 0.022	95.9
	γ 2 2	-0.066 -0.062 0.063 0.056	92.0	-0.061 0.068 0.062	94.9
	L(3,1)	0.188	0.181 0.027 0.031	97.0	0.177 0.033 0.038	97.9
	L(4,2)	0.067	0.066 0.021 0.022	96.0	0.069 0.027 0.026	96.9
	L(3,3)	0.149	0.145 0.013 0.014	95.0	0.143 0.014 0.016	93.9
	L(4,4)	0.250	0.247 0.019 0.019	96.0	0.247 0.022 0.022	93.9
	a 0 11	-0.230 -0.223 0.031 0.035	98.0	-0.219 0.037 0.043	96.9
	a 1 11	0.099	0.100 0.017 0.017	95.0	0.100 0.017 0.019	96.9
	a 0 12	0.118	0.120 0.024 0.025	92.0	0.121 0.027 0.027	93.9
	a 1 12	-0.040 -0.044 0.024 0.023	93.0	-0.045 0.026 0.027	95.9
	a 0 21	0.095	0.095 0.022 0.023	95.0	0.094 0.024 0.025	95.9
	a 1 21	0.043	0.046 0.026 0.027	97.0	0.046 0.027 0.029	98.9
	a 0 22	-0.399 -0.400 0.043 0.042	97.0	-0.401 0.054 0.050	90.9
	a 1 22	0.319	0.325 0.034 0.033	96.0	0.321 0.040 0.037	95.9
	η 01	3.793	3.794 0.188 0.195	94.0	3.795 0.195 0.199	93.9
	η 02	2.601	2.608 0.135 0.127	93.0	2.607 0.148 0.136	93.9
	η 11	1.597	1.594 0.026 0.027	97.0	1.593 0.027 0.028	96.9
	η 12	1.226	1.225 0.026 0.025	93.0	1.225 0.026 0.027	96.9
	σ 1	0.397	0.399 0.015 0.015	95.0	0.399 0.015 0.016	96.9
	σ 2	0.672	0.672 0.031 0.029	96.0	0.673 0.032 0.033	95.9

Table 7 :

 7 Parameter estimates of the measurement model. Parameters (σ 1 , σ 2 , σ 3 ) represent the standard deviation of the measurement errors for the three composite scores. Other parameters given by series of six (η 10 to η 15 ; η 20 to η 25 and η 30 to η 34 ) correspond to the parameters of the quadratic I-splines link function used to transform the anatomic composite score, the cognitive composite score and the functional composite score, respectively.

	score	Parameter Estimate	SE	p-value
		η 10	-5.611	0.249 <0.001
		η 11	0.527	0.115 <0.001
		η 12	0.263	0.097 0.006
	cerebral anatomy score	η 13	0.846	0.046 <0.001
		η 14	0.671	0.074 <0.001
		η 15	0.558	0.208 0.007
		σ 1	0.236	0.008 <0.001
		η 20	-6.429	0.310 <0.001
		η 21	-0.122	0.285 0.668
		η 22	0.579	0.107 <0.001
	cognitive score	η 23	1.436	0.042 <0.001
		η 24	0.461	0.101 <0.001
		η 25	0.884	0.136 <0.001
		σ 2	0.394	0.015 <0.001
		η 30	-7.332	0.299 <0.001
		η 31	0.953	0.082 <0.001
		η 32	0.272	0.135 0.045
	functional score	η 33	0.528	0.064 <0.001
		η 34	0.631	0.041 <0.001
		σ 3	0.581	0.025 <0.001

Table 8 :

 8 Estimates of the transition matrix intensities on ADNI data when considering three discretization steps: 6 months (AIC=1386.68), 3 months (AIC=1370.45) and 2 months (AIC=1367.21).

		δ = 6 months	δ = 3 months	δ = 2 months
		Estimate	SE Estimate	SE Estimate	SE
	α 0	-0.054 0.014	-0.054 0.015	-0.054 0.015
	α 1	0.008 0.009	0.008 0.009	0.008 0.009
	α 2	0.001 0.013	0.001 0.013	0.001 0.013
	α 0	0.010 0.008	0.011 0.008	0.011 0.008
	α 1	0.022 0.010	0.022 0.010	0.022 0.010
	α 2	0.025 0.014	0.025 0.015	0.025 0.015
	α 0	0.010 0.017	0.008 0.017	0.008 0.017
	α 1	-0.012 0.017	-0.011 0.018	-0.010 0.018
	α 2	-0.011 0.019	-0.010 0.020	-0.009 0.020
	α 0	0.077 0.023	0.083 0.026	0.084 0.027
	α 1	0.018 0.026	0.019 0.029	0.019 0.029
	α 2	-0.009 0.034	-0.010 0.037	-0.011 0.039
	α 0	-0.396 0.074	-0.423 0.088	-0.432 0.093
	α 1	0.081 0.030	0.088 0.033	0.090 0.034
	α 2	0.026 0.042	0.029 0.045	0.029 0.047
	α 0	0.146 0.054	0.156 0.062	0.159 0.064
	α 1	-0.068 0.054	-0.073 0.060	-0.074 0.062
	α 2	-0.078 0.059	-0.080 0.066	-0.079 0.068
	α 0	0.077 0.035	0.103 0.048	0.109 0.050
	α 1	0.025 0.043	0.017 0.056	0.016 0.058
	α 2	-0.001 0.052	-0.010 0.068	-0.012 0.071
	α 0	0.132 0.042	0.177 0.062	0.185 0.063
	α 1	0.018 0.048	0.000 0.066	-0.004 0.066
	α 2	0.012 0.061	-0.006 0.079	-0.009 0.082
	α 0	-0.892 0.108	-1.195 0.207	-1.237 0.219
	α 1	0.569 0.078	0.770 0.142	0.798 0.149
	α 2	0.478 0.082	0.667 0.141	0.692 0.148

† Rate of convergence, (15% missing occasions, 7% missing outcomes), ‡ ASE is the asymptotic standard error, ESE is the empirical standard error and CR is the coverage rate of the 95% confidence interval.

† Rate of convergence (15% missing occasions, 7% missing outcomes) ‡ ASE is the asymptotic standard error, ESE is the empirical standard error and CR is the coverage rate of the 95% confidence interval.

Acknowledgements

This work benefited from the support of the project SMALA (ANR-15-CE37-0002) of the French National Research Agency (ANR). Computer time was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) at the Université de Bordeaux and the Université de Pau et des Pays de l'Adour. Data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012).

Goodness-of-fit of the model

We assessed the goodness-of-fit of the model by comparing the subject-specific predicted trajectories (conditionally to the random effects) with the observed trajectories in the transformed scale of the markers. Specifically, we computed the mean predictions and

Appendix C: Additional results in ADNI application