Vincent Y F Tan
email: vtan@mit.edu

Cédric Févotte
email: fevotte@telecom-paristech.fr

Automatic Relevance Determination in Nonnegative Matrix Factorization

Keywords: Nonnegative matrix factorization, Bayesian model selection, Automatic relevance determination

Nonnegative matrix factorization (NMF) has become a popular technique for data analysis and dimensionality reduction. However, it is often assumed that the number of latent dimensions (or components) is given. In practice, one must choose a suitable value depending on the data and/or setting. In this paper, we address this important issue by using a Bayesian approach to estimate the latent dimensionality, or equivalently, select the model order. This is achieved via automatic relevance determination (ARD), a technique that has been employed in Bayesian PCA and sparse Bayesian learning. We show via experiments on synthetic data that our technique is able to recover the correct number of components, while it is also able to recover an effective number of components from real datasets such as the MIT CBCL dataset.

I. INTRODUCTION

Nonnegative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] is a widely used technique that is employed for non-subtractive, partbased representation of nonnegative data. There are numerous diverse applications of NMF including audio signal processing [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence with application to music analysis[END_REF], image classification [START_REF] Guillamet | Analyzing non-negative matrix factorization for image classification[END_REF] and email surveillance [1]. Given a nonnegative data matrix V ∈ R F ×N + , NMF essentially involves finding two nonnegative matrices, namely W = [w 1 , . . . ,

w K ] ∈ R F ×K + and H = [h T 1 , . . . h T K ] T ∈ R K×N + such that V ≈ V ∆ = W H. (1)
K is usually chosen such that F K + KN ≪ F N , hence reducing the data dimension. NMF is also usually performed by minimizing a pre-specified cost function D(V|WH) over W and H. In practice, however, it is often hard to "guess" the effective dimension of the latent subspace K eff i.e. the effective number of columns of W (and the number of rows of H).

To ameliorate this problem, in this paper, we formulate a Bayesian approach to find the effective dimensionality K eff via automatic relevance determination [START_REF] Mackay | Probable networks and plausible predictions -a review of practical Bayesian models for supervised neural networks[END_REF]. This model selection technique was also used by Bishop to determine K eff in the context of Bayesian PCA [START_REF] Bishop | Bayesian PCA[END_REF] and by Tipping in sparse Bayesian learning [START_REF] Tipping | Sparse Bayesian Learning and the Relevance Vector Machine[END_REF]. It is important to discover K eff because we would like a parsimonious (sparse) yet accurate representation of the nonnegative data. Sparsity here refers to the total number of coefficients K(F + N ) ≪ F N required to encode the data, and not only to H as in the standard sparse linear regression setting where the dictionary W is known and fixed.

A. Main contributions

To solve this important model selection problem, we place priors, dependent on precision-like parameters β = [β 1 , . . . , β K ] ∈ R K + , on both the columns of W and the rows of H. The values of these hyperparameters, together with the values of W and H are inferred iteratively by maximizing the posterior of the parameters given the data. As a result of this optimization, a subset of the β k 's will be driven to a large upper bound [START_REF] Bishop | Bayesian PCA[END_REF], [START_REF] Mackay | Probable networks and plausible predictions -a review of practical Bayesian models for supervised neural networks[END_REF], corresponding to irrelevant components. More precisely, we seek to find W * , H * , β * by optimizing the maximum a-posteriori (MAP) criterion:

min W,H,β C MAP (W, H, β) ∆ = -log p(W, H, β|V), (2)
where, by Bayes rule, the posterior can be written as

-log p(W, H, β|V) c = -log p(V|W, H) -log p(W|β) -log p(H|β) -log p(β). (3)
We use c = to denote equality up to a constant. For some statistical models, maximizing the log-likelihood term log p(V|W, H) in (3) over (W, H) is equivalent to minimizing a specific cost function (or measure of fit) D(V|WH) over (W, H). For instance, maximum likelihood estimation of W and H in Poisson noise model corresponds to minimizing the (generalized) Kullback-Leibler (KL) divergence D KL (V|WH) = D KL (V| V), where

D KL (V| V) ∆ = f n v f n log v f n vfn -v f n + vfn . (4)
We focus on the KL-divergence as our cost function in this paper because it has been extensively used in NMF and is free of noise parameters. The model order selection technique we present in this paper naturally extends to other cost functions such as the Euclidean distance or Itakura-Saito divergence, that may also be mapped to log-likelihood functions [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence with application to music analysis[END_REF]. The Euclidean or Itakura-Saito cost functions, however, require the specification or inference of noise parameters. We show in Section IV that our method is able to recover the correct number of components for synthetic datasets, including the Swimmer dataset [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts?[END_REF], and an effective number of component for the MIT CBCL face dataset.

B. Related Work

There is generally little literature about model order selection in NMF. Variational Bayesian methods have been considered for this task [START_REF] Cemgil | Bayesian Inference for Nonnegative Matrix Factorisation Models[END_REF], [START_REF] Winther | Bayesian independent component analysis: Variational methods and non-negative decompositions[END_REF] but this type of approach is fairly computationally demanding and require successive evaluation of the evidence (or marginal likelihood) for each considered value of K. Instead, our method automatically finds the optimal model order K eff given a large initial value of K and is computationally less involving.

The methodology presented in this work bears resemblance with existing multiplicative NMF algorithms with a sparseness constraint on the whole matrix H (typically, adding a H 1 or H 0 regularization term to the cost function D(V WH)), with which is was noted that the sparseness constraint may turn off excess components [START_REF] Mørup | Approximate L0 constrained non-negative matrix and tensor factorization[END_REF]. Here, we do not impose sparseness per se on either W or H, but rather assign a relevance weight β k to each rank-1 matrix w k h k , which is learnt from the data.

II. BAYESIAN NMF MODEL

In this section, we describe our statistical model. We will describe in detail each of the terms in (3), namely the log-likelihood term log p(V|W, H) and the log-prior terms log p(W|β), log p(H|β) and log p(β).

A. Likelihood model

We assume that the log-likelihood of a single element of the matrix V, denoted p(v f n |v f n ), is given by P(v f n |v f n ), where P(x|λ) = e -λ λ x /Γ(x + 1) is the Poisson probability density function with rate λ. Thus, E[v f n |v f n ] = vfn . The log-likelihood can then be written as:

-logp(V|W, H) = D KL (V|WH) + f n v f n (1 -log v f n ) + log Γ(v f n + 1). (5)
Thus, maximizing the log-likelihood (based on the Poisson noise model) is equivalent to minimizing the KL-divergence since other terms in (5) are independent of W and H.

B. Prior model on W and H

In Bayesian PCA [START_REF] Bishop | Bayesian PCA[END_REF], each column k of W (resp. row k of H) is given a normal prior with precision parameter β k (resp. 1). Thanks to the simplicity of the statistical model (multivariate Gaussian observations with Gaussian parameter priors), H can be easily integrated out of the likelihood, and optimization can be done over p(W, β|V). This marginalization step is not easily done in NMF. However, similar in spirit to Bayesian PCA, we define independent half-normal priors over each column k of W and row k of H, and the priors are tied together through a single, common precision parameter β k .1 More precisely, we set

p(w f k |β k ) = HN (w f k |0, β -1 k ), (6) p(h kn |β k ) = HN (h kn |0, β -1 k ), (7)
where

HN (x|0, β -1 ) = 2 π β -1/2 exp - 1 2 βx 2 (8)
is the half-normal probability density function (defined for x ≥ 0) parameterized by the precision (inverse variance) β.2 From ( 6) and ( 7), we see that the minus log-priors can be written as:

-log p(W|β) c = k f 1 2 β k w 2 f k - F 2 log β k , (9)
log p(H|β)

c = k n 1 2 β k h 2 kn - N 2 log β k . ( 10
)
As a result of inference, a subset of the β k 's will be driven to a large upper bound, with the corresponding columns of W and rows of H driven to small values, by [START_REF] Mackay | Probable networks and plausible predictions -a review of practical Bayesian models for supervised neural networks[END_REF]. The effective dimensionality can be deduced from the distribution of the β k 's, from which, we have found in practice, two clusters clearly emerge : a group of values in same order of magnitude corresponding to relevant components and a group of similar values of much higher magnitude corresponding to irrelevant components. We may then define K ef f as

K eff ∆ = |{β k : β k < L k -ǫ}|, (11)
where L k is the upper bound dependent on the prior's parameters and ǫ > 0 is a user-defined small constant. We will specify a precise value for L k in terms F , N and the parameters of the prior on β k in Section III. This Bayesian approach based on ARD automatically determines the number of components in V, thus elegantly finding a solution for the odd dichotomy between finding the best model for the training data to the model and avoiding overfitting. Furthermore, the size of β k provides an approximate measure of its relevance. If β k is large (variance 1/β k is small), then this component is not very relevant for the modeling of the data. Conversely, if β k is small (variance 1/β k is large), the corresponding component contains a significant amount of energy. We term β k the relevance weight associated to component k.

We remark that it is also possible to let W be dependent of β and let H be independent of β. In this setup, we place a prior on the columns of W only. Doing this is akin to Bayesian PCA [START_REF] Bishop | Bayesian PCA[END_REF]. However, in our experiments, we found convincing evidence that tying w k and h k together with the same prior structure and precision parameter is better for recovery of the number of components. It is also possible to put exponential distributions, (dependent on the precision parameters β 1 , . . . , β K ) instead half-normals, on either W or H (or both). Indeed, the inference method described next would only need minor changes to accomodate other priors.

C. Prior model on β

Finally, each precision parameter β k is given a Gamma distribution, which is conjugate to the half-normal density function [START_REF] Mackay | Probable networks and plausible predictions -a review of practical Bayesian models for supervised neural networks[END_REF]. Thus,

p(β k |a k , b k ) = b a k k Γ(a k ) β a k -1 k exp(-β k b k ), β k ≥ 0, (12)
and the (minus) log-prior of β is given as

-log p(β) c = k b k β k -(a k -1) log β k . (13)
In our experiments, the hyperparameters a k and b k (also known as shape and scale parameters) are fixed to a prespecified common value, say a and b. Thus, all the upper bounds L k = L, some constant dependent on a, b.

D. Overall Model and Cost Function

All the terms in in the MAP criterion (3) have been defined. The dependences between the variables can be represented by a directed graphical model (Bayesian network) as shown in Fig. 1. The observations {v f n } are contained in the matrix V and we would like to infer W, H and β. The number of β k 's that is less than the upper bound L kǫ is the effective dimensionality of the data.

By combining equations ( 5), ( 9), ( 10) and ( 12), the overall MAP objective to be minimized can be written as

-log p(W, H, β|V) c = D KL (V|W, H) + 1 2 k f w 2 f k + n h 2 kn + 2b k β k -(F + N -2(a k -1)) log β k . (14)
From this expression, we see a tradeoff involving the size of the β k 's. Since F + N is typically larger than 2(a k -1), a larger β k 's will result in a smaller (more negative) last term. However, the second to last term in the sum forces some of the β k 's to remain small since it is a linear function of β k .

III. INFERENCE

In [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF], the authors proposed efficient coordinate descent algorithms with suitable step sizes which are then turned into multiplicative update rules. This simple yet effective procedure has the advantage of maintaining nonnegativity of the inferred matrices W and H. The same procedure can be applied to our model to infer W, H as well as β. We first find the gradient of the MAP criterion C MAP with respect to W, H and β k .

∇ W C MAP (W, H, β) = WH -V WH H T + Wdiag (β) , (15)
∇ H C MAP (W, H, β) = W T WH -V WH + diag (β) H, (16)
Algorithm 1 Automatic relevance determination for NMF with the KL-divergence cost Input : Nonnegative data (observation) matrix V, fixed hyperparameters a k , b k .

Output : Nonnegative matrices W and H such that V ≈ V = W H, a nonnegative vector β and K eff .

Initialize W and H to nonnegative values.

for i = 1 : n iter do H ← H W T 1 F ×N +diag(β)H • W T V WH W ← W 1F ×N H T +Wdiag(β) • V WH H T β ← F +N +2(a-1) 11×F (W•W)+(H•H)1N×1+2b
end for Compute K eff as in [START_REF] Winther | Bayesian independent component analysis: Variational methods and non-negative decompositions[END_REF].

∇ β k C MAP (W, H,β) = 1 2 f w 2 f k + 1 2 n h 2 kn + b k - 1 2 (N + F ) + a k -1 1 β k . ( 17
)
where we use the notation A B to denote entry-wise division of A by B. The multiplicative coordinate descent algorithm is simple; we simply update each parameter θ of W and H by multiplying its current value with the ratio of the positive 15) and ( 16)) with respect to θ i.e.,

[•] + to negative [•] -part of the derivative of C MAP (in (
θ ← θ [∇ θ C MAP (θ)] + [∇ θ C MAP (θ)] - . (18)
The precision parameters β k are updated by zeroing (17) i.e.,

β k ← 1 2 (F + N ) + a k -1 1 2 ( f w 2 f k + n h 2 kn ) + b k . ( 19
)
The overall algorithm is summarized in Algorithm 1. We observed in practice the monotonicity of the MAP criterion in (2) under this algorithm. We now address the issue of choosing L k = L in [START_REF] Winther | Bayesian independent component analysis: Variational methods and non-negative decompositions[END_REF]. From (19), we can easily derive an upper bound on β k , which is dependent only on a and b, since we set a k = a and b k = b:

β k ≤ F + N + 2(a -1) 2b . ( 20
)
The bound is attained when w f k = h kn = 0 for all f, n, k i.e., the columns and rows of W and H are exactly zero. Thus, a good choice of L is precisely the right-hand-side of (20). Hence, we define

L ∆ = F + N + 2(a -1) b , (21)
and the constant L will enable us to determine K eff from the set of β k 's from [START_REF] Winther | Bayesian independent component analysis: Variational methods and non-negative decompositions[END_REF].

It is worth noting that other inference algorithms (on graphs) such as variational Bayes or Gibbs sampling may be employed to find the posterior distributions of the parameters of interest. In this paper, we are only interested in point estimates of W, H and β, hence the efficient multiplicative update rules suffice.

IV. NUMERICAL RESULTS

In this section, we describe some experiments on synthetic and real data that demonstrate the efficacy of our algorithm for model order selection in NMF.

A. Synthetic Data

Our synthetic data was generated in the following fashion: W and H are 100 × 10 and 10 × 1000 matrices respectively. Each element in the first 5 columns of W (and the first 5 rows of H) are drawn uniformly from a half normal distribution HN (0, 10) while the remaining 5 columns (remaining 5 rows of H) are drawn from HN (0, 1). Thus the effective dimensionality is K eff = 5. The 100 × 1000 data matrix V was formed by multiplying W and H. The hyperparameters were set to a = b = 1, though in practice we found that the algorithm is not very sensitive to the choice of hyperparameters for the synthetic data. In Figs 2, we see that the algorithm correctly recovers the correct dimensionality of the data even though Rank(V) = K = 10 with probability 1. In Fig 3, we observe that K -K eff = 5 β k 's diverge to a large value, which means that 5 columns of W (and 5 rows of H) are irrelevant for the purpose of modeling the N = 1000 data samples. We repeated the same procedure by fixing K = 10 and varied K eff from 1 to 5. The algorithm correctly recovers K eff in each case.

B. Swimmer Dataset

We also performed some experiments on another synthetic dataset, namely the well-known swimmer dataset [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts?[END_REF], which depicts a figure with four moving parts (limbs), each able to exhibit four articulations (different positions). Thus, there are a total of N = 4 4 = 256 images each of size F = 32 × 32. Sample images are shown in Fig. 4. The objective here is to discover that the dataset contains 16 basis elements, corresponding to the 16 limb positions. Thus, K eff = 16. We set K = 50. In this experiment, we also found that we could retrieve the correct number of components, if b in the Gamma prior in ( 12) is larger than some threshold value. Keeping a fixed at 2, we plotted the regularization path in Fig. 5. The regularization path shows how K eff changes as b in the Gamma prior in (12) changes, while being insensitive to the precise value of a. The components, which are the limbs, are plotted in Fig. 6 for b = 18. We see that each of the 16 limb positions is correctly recovered, together with the torso. All the relevance weights corresponding to the relevant components are equal, by equivalence of all the "parts" represented in the dataset (β k = 2.57 for all k = 1, . . . , K eff ). This experiment further validates our technique of automatically discovering the number of latent components in the nonnegative data.

C. MIT CBCL Faces Dataset

Finally, we applied our model order selection procedure to the CBCL face dataset, containing F = 19 × 19 images, preprocessed as in [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF], i.e., the greyscale intensities of each image were linearly scaled so that the pixel mean and standard deviation were equal to 0.25, and then clipped to the range [0,1]. The training set contains N = 2429 images. Keeping a (the shape parameter in the Gamma prior in ( 12)) fixed at 2, we plot the regularization path in Fig. 7, i.e., the effective number of components K eff against the scale parameter b. This figure shows that if b is above some threshold value, the number of components K eff stabilizes at a constant value (K eff = 12). K was set to 49 initially. For instance, if we set b = 25, we observe that 12 basis images (columns of W) are recovered. These are shown in Fig. 8, together with the corresponding the relevance weights β k 's. They correspond to several parts of the faces : we roughly obtain left, right, top and bottom halves, background artifacts and eyes/eyebrows, mouth and nose parts.

Our adaptation of ARD to NMF allows us to automatically discover a latent dimensionality of this dataset and the recovered basis images make intuitive sense. However they are not as local as the one obtained with the standard NMF algorithm with K set to 49, as in the original experiment by Lee & Seung [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF]. Other prior structures (such as the exponential distribution instead of the half-normal, or even sharper priors) may lead to more local features. However the results are subjective in every case and should instead be evaluated on well-defined tasks such as classification.

V. CONCLUSIONS

We have presented a Bayesian approach that performs model order selection for NMF by borrowing ideas from ARD. By placing appropriate prior distributions on the elements on W and H, we are able to identify those components that are 'relevant' for modeling the data. The efficacy of our method was adequately demonstrated in our experiments where we are able to recover the latent dimensionality of synthetic and real data.

Indeed, our method is not limited to the use of the generalized KL-divergence; we could, in fact, adapt the updates in a very straightforward manner for the Euclidean or Itakura-Saito costs and also use other types of priors. Furthermore, our method can be extended for the purpose of nonnegative tensor decomposition (NTD). Our preliminary results show that it is possible to recover the correct number of components from a synthetically generated 3-way tensor formed from the parallel factor analysis (PARAFAC), also named canonical decomposition (CANDECOMP), model. We hope to further extend the theory and experimentation for tensor decomposition in a longer paper.

Figure 1 :

 1 Figure 1: A directed graphical model that describes our NMF statistical model. The plate notation is used here. The plate (or rectangle) is used to group random variables that repeat. The number of replicates is shown on the bottom right corner. The shaded variables {v f n } are the observations.

Figure 2 :Figure 3 :

 23 Figure 2: The W and H matrices showing that K eff = 5 columns are retained. The other columns of W and rows of H are set to 0.

Figure 4 :

 4 Figure 4: Sample images from the Swimmer dataset.

Figure 5 :

 5 Figure 5: Regularization Path for Swimmer Dataset. The effective number of components K eff is plotted against b, a parameter of the Gamma prior distribution in (12).

Figure 6 :

 6 Figure 6: Basis images obtained by performing NMF on the Swimmer dataset. The 16 limb positions of the swimmer dataset are correctly recovered if b ≥ 18. All the relevance weights {β k } are equal by symmetry.

Figure 7 :

 7 Figure 7: Regularization Path for MIT CBCL Dataset.

Figure 8 :

 8 Figure 8: Basis images (reordered) and their corresponding relevance weights {β k } obtained by performing Bayesian model order selection using ARD for NMF on the MIT CBCL dataset. The images correspond to the parts of the face including left, right, top and bottom halves, the eyes, the mouth, background and nose.

Note that this prior is not overconstraining the scales, because of the scale indeterminacy between w k and h k .

If X is a zero mean normal random variable with variance β -1 , then |X| is a half-normal random variable with precision β.

Acknowledgements

The work of V. Y. F. Tan was supported by A*STAR, Singapore, and by a MURI funded through ARO Grant W911NF-06-1-0076. C. Févotte acknowledges support from the ANR -SARAH project (Standardisation du Remastering Audio Haute-définition).


Marion Brandolini-Bunlon

Estelle Pujos-Guillot

Blandine Comte

Benoit Jaillais

Véronique Cariou

Evelyne Vigneau

Global and partial effect assessment in metabolic syndrome explored by metabolomics

Keywords: multiblock, metabolomics, clinical data, effects, partial correlation, metabolic syndrome

Introduction

In health-related case-control studies, untargeted metabolomics is routinely involved in addition to standard clinical information in order to better predict and understand a syndrome of interest. This leads to a supervised model, with several input blocks (metabolomics, clinical data) being potential predictors of the targeted output to be explained. Such a model is commonly represented with a path diagram where directed links connect each input block to the output, and input blocks may be interrelated. Within this scope, the way to determine the effect of an input block on the output response is detailed depending on whether it is a global effect or conditionally to another block.

Material and methods

A case-control study on metabolic syndrome (MetS) was designed within the Quebec Longitudinal cohort on Nutrition and Successful Aging (NuAge) [START_REF] Gaudreau | Nutrition as a determinant of successful aging: Description of the Quebec longitudinal study NuAge and results from cross-sectional pilot studies[END_REF]. Male subjects (n=121), either affected (cases) or free (controls) of MetS and stable during the 3 years of follow-up, were selected. A binary variable, y, indicates this MetS presence. A clinical data block, Clinic, included the 6 quantitative diagnostic variables of the MetS, collected at baseline. Baseline and endpoint serum samples were also analyzed using 6 different untargeted metabolomic and lipidomic methods. Operating procedures were standardized from sample preparation to data processing. Raw data were extracted and pre-processed (quality checks and signal drift correction) to yield a data matrix containing retention times, masses, and peak intensities corrected for batch effects. The analytical redundancy inside each dataset was reduced. All the variables stable over time regardless of their link with MetS were submitted to a feature selection step to select 102 metabolomic variables predictive of MetS, constituting the metabolomic data block, Metabo [START_REF] Comte | Multiplatform metabolomics for an integrative exploration of metabolic syndrome in older men[END_REF]. A logarithmic transformation was applied to this data block. Finally, only subjects at baseline without missing values (54 cases / 45 controls) were kept.

The output block was y, and the two explanatory input data blocks were Metabo and Clinic, whose variables were scaled to unit variance. Two different path diagrams (Figure 1) were investigated. In each one, one of the input block is exogenous, designated as being the block A, and surrounded by a rectangle in Fig. 1. For the path 1, A = Metabo, and for the path 2, A = Clinic. The output block, C, is the response, y, surrounded by a circle in Fig. 1, and the mediating block is denoted B.

The effect of the input block A on the output block C is determined as the cross-validated percentage of explained variance of C, obtained by a PLS regression model of A on C. This criterion is called 'global effect' of A on C ('total effect' in the terminology adopted by Naes et al. (2020) [START_REF] Naes | Sequential and orthogonalized PLS (SO-PLS) regression for path analysis: Order of blocks and relations between effects[END_REF]). When A and C have only one variable each, and cross-validation is not performed, the percentage of explained variance is equal to the squared correlation coefficient of A and C (𝑟 𝐴,𝐶 2 ). In the case where all the blocks have only one variable each, and cross-validation is not performed, the squared partial correlation coefficient between A and C, conditionally to B (𝑟 𝐴,𝐶/𝐵 2) is obtained.

Results and discussion

For both paths, global and partial effect estimations, by means of PLS2 regression models, are given in Table 1. For each model involved, the choice of the optimal number of PLS components (# PLS comp. in Table 1) and the computation of the percentages of explained variance with their standard deviation (expl.var. ± SD in Table 1) result from a 10 folds-CV, repeated 50 times. From these models, the selection of the variables with the largest VIP values led to identify metabolomic and clinical variables directly related to carbohydrate and lipid metabolisms as good predictors of MetS. Interestingly, in path 1, metabolomic variables related to other metabolic disturbances were brought to light thanks to the inspection of the partial effect; and in path 2, it highlighted clinical variables related to functional disturbances.

Conclusion

In our study, the determination of both global and partial effects together with the identification of the most important variables from the associated models highlight the redundancy vs complementarity of the clinical and metabolomic information in the MetS explanation.

Figure 1 -

 1 Figure 1 -Two different path diagrams. The 'partial effect' of the input block A on the output block C, taking into account the mediating block B (the 'direct effect' in Naes et al. (2020) [3]), is estimated as the cross-validated explained variance of the residuals of C on B (𝐄 C⟵B ) by the residuals of A on B (𝐄 A⟵B ). It represents the part of the variance in C that may be explained by A, but not explained by B (it does not go through B).In the case where all the blocks have only one variable each, and cross-validation is not performed, the squared partial correlation coefficient between A and C, conditionally to B (𝑟 𝐴,𝐶/𝐵

Table 1 -

 1 Global and partial effects estimated for both paths described in Figure1.

	effects	path 1 -expl.var. ± SD (# PLS comp.) path 2 -expl.var. ± SD (# PLS comp.)
	global (AC)	53.43 ± 1.47 % (2)	52.37 ± 0.74 % (1)
	partial (AC/B)	21.67 ± 3.83 % (2)	22.95 ± 1.85 % (2)


Moufidath Adjibade

Camille Davisse-Paturet

Jonathan Y Bernard

Karine Adel-Patient

Amandine Divaret-Chauveau

Sandrine Lioret

Marie-Aline Charles

Blandine De Lauzon-Guillain

Enrichment of infant formula with long-chain polyunsaturated fatty acids and risk of infection and allergy in the nationwide ELFE birth cohort Running title: LCPUFA-enriched formula and infection or allergy

Keywords: allergies, infant formula, infections, LCPUFA-enrichment

Background:

The new European regulations require the enrichment of formulas with docosahexaenoic acid (DHA) because of the positive effects of long-chain polyunsaturated fatty acids (LCPUFAs) on neurodevelopment and visual acuity. In this observational study, we aimed to evaluate whether the consumption of LCPUFA-enriched formula was associated with the risk of infection and allergy in early childhood.

Methods: Analyses involved data from 8,389 formula-fed infants from the ELFE birth cohort.

Formula enrichment was identified from the list of ingredients of the formula consumed at 2 months. Infections (gastro-intestinal, lower respiratory tract [LRTI], upper respiratory tract) and allergies (wheezing, itchy rash, asthma medication, food allergy) from age 2 months to 5.5 years were reported by parents during follow-up surveys. Multivariable logistic regression models were used to assess associations between the consumption of LCPUFA-enriched formula and the risk of infection and allergy.

Results: Among formula-fed infants at 2 months, 36% consumed formula enriched with DHA and arachidonic acid (ARA), and 11% consumed formula additionally enriched with eicosapentaenoic acid (EPA). Enriched formula consumption was not associated with infection or allergy, except for an association between consumption of DHA/ARA/EPAenriched formula and lower use of asthma medications. Furthermore, as compared with non-DHA/ARA/EPA-enriched formula, consumption of formula with high EPA content (≥3.2 mg/100 kcal) was related to lower risk of LRTI and lower use of asthma medications.

Conclusion:

This study suggests that consumption of DHA/ARA/EPA-enriched formula (especially those with high EPA content) is associated with a lower risk of LRTI and lower use of asthma medications.

methodological guidance; MA, CD-P, JYB, KA-P, AD-C, SL, M-AC and BL-G were involved in interpreting the results and editing the manuscript for important intellectual content. All authors read and approved the final manuscript.

Financial support:

This study is funded by an ANR grant (InfaDiet project, ANR-19-CE36-0008).

The ELFE survey is a joint project between the French Institute for Demographic Studies The funders had no role in the study design, data collection and analysis, decision to publish or preparation of the manuscript.

INTRODUCTION

Allergic diseases and infections are the leading causes of morbidity in childhood. The prevalence of infections has decreased in recent years [START_REF] Pawankar | The WAO White Book on Allergy[END_REF][START_REF] Collaborators | Estimates of the global, regional, and national morbidity, mortality, and aetiologies of lower respiratory infections in 195 countries, 1990-2016: a systematic analysis for the Global Burden of Disease Study 2016[END_REF][START_REF] Collaborators | Estimates of the global, regional, and national morbidity, mortality, and aetiologies of diarrhoea in 195 countries: a systematic analysis for the Global Burden of Disease Study 2016[END_REF] , but that of allergic diseases and asthma (both considered a major public health issue) has increased over the last decades [START_REF] Pawankar | The WAO White Book on Allergy[END_REF][START_REF] Asher | Worldwide time trends in the prevalence of symptoms of asthma, allergic rhinoconjunctivitis, and eczema in childhood: ISAAC Phases One and Three repeat multicountry cross-sectional surveys[END_REF][START_REF] Nwaru | Prevalence of common food allergies in Europe: a systematic review and meta-analysis[END_REF] .

In addition to genetic predisposition, several modifiable factors are associated with risk of developing allergic diseases, including low microbial exposure, high exposure to air pollutants or cigarette smoke and inadequate early feeding practices [START_REF] Pawankar | The WAO White Book on Allergy[END_REF][START_REF] Beasley | Risk factors for asthma: is prevention possible?[END_REF][START_REF] Halken | Prevention of allergic disease in childhood: clinical and epidemiological aspects of primary and secondary allergy prevention[END_REF] . With regard to early feeding, the protective effect of breastfeeding on infections is well established [START_REF] Victora | Breastfeeding in the 21st century: epidemiology, mechanisms, and lifelong effect[END_REF] , but its potential effect on allergic diseases is more controversial [START_REF] Victora | Breastfeeding in the 21st century: epidemiology, mechanisms, and lifelong effect[END_REF][START_REF] Horta | Short-term effects of breastfeeding: a systematic review of the benefits of breastfeeding on diarrhoea and pneumonia mortality[END_REF][START_REF] Gungor | Infant milk-feeding practices and food allergies, allergic rhinitis, atopic dermatitis, and asthma throughout the life span: a systematic review[END_REF][START_REF] Tromp | Breastfeeding and the risk of respiratory tract infections after infancy: The Generation R Study[END_REF][START_REF] Munblit | Human Milk and Allergic Diseases: An Unsolved Puzzle[END_REF] .

Breastmilk supplies all essential bioactive factors and nutrients, including long-chain polyunsaturated fatty acids (LCPUFAs), which have an important role in immune system development and regulation and in neurodevelopment, including neurotransmission, visual function and cognitive development [START_REF] Gould | The effect of maternal omega-3 (n-3) LCPUFA supplementation during pregnancy on early childhood cognitive and visual development: a systematic review and meta-analysis of randomized controlled trials[END_REF][START_REF] Janssen | Long-chain polyunsaturated fatty acids (LCPUFA) from genesis to senescence: the influence of LCPUFA on neural development, aging, and neurodegeneration[END_REF][START_REF] Richard | Evidence for the essentiality of arachidonic and docosahexaenoic acid in the postnatal maternal and infant diet for the development of the infant's immune system early in life[END_REF] . LCPUFA content in breastmilk varies according to the lactation period and maternal diet, with an average content of 0.37% of total fatty acids for docosahexaenoic acid (DHA), 0.55% of total fatty acids for arachidonic acid (ARA) and 0.10% of total fatty acids for eicosapentaenoic acid (EPA) [START_REF] Koletzko | Human Milk Lipids[END_REF][START_REF] Bobiński | Fatty acids of human milk -a review[END_REF][START_REF] Fu | An updated review of worldwide levels of docosahexaenoic and arachidonic acid in human breast milk by region[END_REF] . Several studies have examined the association between omega-3 (n-3) LCPUFAs level in breastmilk and risk of allergic diseases during infancy, but the findings remained inconsistent [START_REF] Venter | EAACI position paper: Influence of dietary fatty acids on asthma, food allergy, and atopic dermatitis[END_REF][START_REF] Waidyatillake | Association of breast milk fatty acids with allergic disease outcomes-A systematic review[END_REF] . Several randomized controlled trials (RCTs) have also assessed the effect of prenatal or postnatal n-3 LCPUFA supplementation on allergy development, but the available systematic reviews on this topic found that evidence remains insufficient to recommend supplementation with n-3 LCPUFAs during pregnancy, lactation and early life to prevent the development of allergic diseases in childhood [START_REF] Venter | EAACI position paper: Influence of dietary fatty acids on asthma, food allergy, and atopic dermatitis[END_REF][START_REF] Waidyatillake | Association of breast milk fatty acids with allergic disease outcomes-A systematic review[END_REF][START_REF] Gunaratne | Maternal prenatal and/or postnatal n-3 long chain polyunsaturated fatty acids (LCPUFA) supplementation for preventing allergies in early childhood[END_REF][START_REF] Schindler | Polyunsaturated fatty acid supplementation in infancy for the prevention of allergy[END_REF][START_REF] Vandenplas | The role of milk feeds and other dietary supplementary interventions in preventing allergic disease in infants: Fact or fiction?[END_REF] .

Exclusive breastfeeding is recommended for the first 6 months of life or at least 4 months [START_REF] Who | Feeding and nutrition of infants and young children, guidelines for the WHO European region, with emphasis on the former Soviet countries[END_REF] , but several studies have shown a high level of non-compliance with these recommendations 8, 25, 26 . For non-breastfed or partially breastfed infants, a wide range of infant formula is available, especially in France [START_REF] De Lauzon-Guillain | Use of infant formula in the ELFE study: The association with social and health-related factors[END_REF] . From scientific evidence suggesting beneficial effects of DHA on cognitive abilities [START_REF] Drover | Cognitive function in 18-month-old term infants of the DIAMOND study: a randomized, controlled clinical trial with multiple dietary levels of docosahexaenoic acid[END_REF][START_REF] Lepping | Long-chain polyunsaturated fatty acid supplementation in the first year of life affects brain function, structure, and metabolism at age nine years[END_REF][START_REF] Liao | Event-related potential differences in children supplemented with long-chain polyunsaturated fatty acids during infancy[END_REF] and to meet the nutritional needs of healthy infants, a European regulation (EU 2016/127) required that from February 2020, all infant and followon formula had to be enriched in DHA (20-50 mg/100 kcal, equivalent to 0.5-1% of total fatty acids). The European Academy of Pediatrics further suggested that DHA enrichment had to be combined with ARA enrichment 31 .

To our knowledge, only 4 studies of small sample sizes have specifically investigated the association between DHA/ARA enrichment of formula and infection or allergy in childhood [START_REF] Foiles | Formula with long-chain polyunsaturated fatty acids reduces incidence of allergy in early childhood[END_REF][START_REF] Birch | The impact of early nutrition on incidence of allergic manifestations and common respiratory illnesses in children[END_REF][START_REF] Lapillonne | Infants fed formula with added long chain polyunsaturated fatty acids have reduced incidence of respiratory illnesses and diarrhea during the first year of life[END_REF][START_REF] Pastor | Infants fed docosahexaenoic acid-and arachidonic acid-supplemented formula have decreased incidence of bronchiolitis/bronchitis the first year of life[END_REF] . These studies support a beneficial effect of this enrichment on the incidence of diarrhea and some allergic and respiratory diseases [START_REF] Foiles | Formula with long-chain polyunsaturated fatty acids reduces incidence of allergy in early childhood[END_REF][START_REF] Birch | The impact of early nutrition on incidence of allergic manifestations and common respiratory illnesses in children[END_REF][START_REF] Lapillonne | Infants fed formula with added long chain polyunsaturated fatty acids have reduced incidence of respiratory illnesses and diarrhea during the first year of life[END_REF][START_REF] Pastor | Infants fed docosahexaenoic acid-and arachidonic acid-supplemented formula have decreased incidence of bronchiolitis/bronchitis the first year of life[END_REF] . However, no study has examined the potential effect of EPA enrichment of infant formula on infection and allergy in childhood.

The purpose of this observational study was to evaluate among formula-fed infants the associations between the enrichment of infant formula with LCPUFAs (DHA, ARA and EPA) and the risk of infection and allergy in early childhood in a large French nationwide birth cohort. The infant formula consumed in the present study were those available on the market in France before the enforcement of the new European regulation on DHA enrichment of formula.

We hypothesized that (1) consumption of LCPUFA-enriched formula at 2 months of age is associated with a lower risk of infection and allergy up to 5.5 years; (2) the observed risk reductions is greater for high LCPUFA contents in the consumed formula.

MATERIALS AND METHODS

Study population

Our analyses were based on data from the ELFE (Étude longitudinale française depuis l'enfance) multidisciplinary nationwide birth cohort including children born in France in 2011 [START_REF] Charles | Cohort Profile: The French national cohort of children (ELFE): birth to 5 years[END_REF] . This study was designed to assess the impact of environmental exposures and socioeconomic and familial factors on children's development, behaviours and health. The design and protocol of the study have been described in detail elsewhere [START_REF] Vandentorren | Rationales, design and recruitment for the Elfe longitudinal study[END_REF] . Briefly, the study enrolled singletons or twins born after 33 weeks' gestation, from mothers aged 18 years or older who were able to read French, Arabic, Turkish or English and were not planning to move outside of metropolitan France in the next 3 years.

All participating mothers provided written consent for themselves and their child. Fathers also signed the consent form for their child's participation when present at inclusion or were

Data collection

Infant feeding

Infant feeding was collected at age 2 months by phone interview. For infants consuming infant formula, parents were asked to report the brand and product names of the mainly used formula. All formulas used in the ELFE study were then listed and classified as enriched in LCPUFAs (including DHA, ARA and EPA) or not according to the ingredient list and, when available, the nutrient composition. When these data were not available, we used data from another question that directly reported whether the formula was enriched or not with LCPUFAs. Parents were also asked to report the age at infant formula introduction, which was used to derive the duration of predominant breastfeeding. Infants who received both breastmilk and infant formula were classified based on the formula consumed.

The enrichment level of the LCPUFA-enriched formula identified in this study ranged from 9.0 to 20.3 mg/100 kcal for DHA and 2.1 to 3.3 mg/100 kcal for EPA. When nutrient composition was available, the level of DHA and EPA enrichment was classified as low or high according to the median calculated in the study population (10 mg/100 kcal, equivalent to 0.19% of total fatty acids for DHA; and 3.2 mg/100 kcal, equivalent to 0.06% of total fatty acids for EPA). The median value was used because only 1 infant who consumed LCPUFAenriched formula had consumed a formula with a DHA content at the minimum dose now recommended by the European regulation (20 mg/100 kcal) and only 75 (3%) infants the dose used in most RCTs (17 mg/100 kcal). Linoleic acid and alpha-linolenic acid contents were also available and modeled as continuous variables.

Infection and allergy

Data on infection and allergy were collected by phone interviews at the 2-month and, 1-, 2-, 3.5-and 5.5-year follow-ups, each referring to events that occurred since the last follow-up.

To define incident cases of infection and allergy from age 2 months to 5.5 years, only cases occurring after the 2-month follow-up were considered.

In this observational study, infections were assessed by using 3 indicators: ever gastrointestinal infection, ever lower respiratory tract infection (LRTI) and ever upper respiratory tract infection (URTI), which were defined as follows:

-Ever gastrointestinal infection: parental report of at least one hospitalization for gastroenteritis/dehydration or one emergency consultation for diarrhea/vomiting/dehydration; -Ever LRTI: parental report of at least one bronchiolitis event or one hospitalization for bronchitis, bronchiolitis or pneumopathy; -Ever URTI: parental report of at least 3 otitis events from birth, or at least 3 laryngitis events from birth or at least 3 angina events from birth, or one hospitalization URTI or one emergency consultation for otalgia.

With regard to allergies, we considered the following indicators to enhance international comparison with other birth cohorts from the LifeCycle consortium [START_REF] Jaddoe | The LifeCycle Project-EU Child Cohort Network: a federated analysis infrastructure and harmonized data of more than 250,000 children and parents[END_REF] : ever wheezing, ever itchy rash, ever food allergy and ever use of asthma medications. Children were classified as with ever wheezing or ever itchy rash from age 2 months to 5.5 years if parents reported at least one episode. Children were classified as with ever food allergy from age 2 months to 5.5 years if parents reported at least once medical advice to avoid certain foods due to an allergy. 'Ever use of asthma medications' refers to at least one parental report of use of inhaled corticosteroids or bronchodilators at home or the hospital during the first 3.5 years of life (these data were not collected at the 5.5-year follow-up).

Infant and familial characteristics

Maternal and household characteristics were collected at the maternity and 2-month interviews. Birth characteristics were also collected from the medical records. Maternal and household characteristics considered in this study were defined as: age (< 25, 25-29, 30-34, ≥ 35 years), city size (rural, urban), region of residence (Paris region, North, East, Paris Basin-East, Paris Basin-West, West, South-West, South-East, Mediterranean), education level (up to lower secondary, upper secondary, intermediate, 3-year university degree, at least 5-year university degree), employment status in early pregnancy (employed, unemployed, out of the labor force [student/housewife/retired/handicapped pension]), pre-pregnancy body mass index (<18.5, 18.5-24.9, 25-29.9, ≥30 kg/m²), maternal allergy history (including eczema, asthma, or hay fever; yes/no), smoking status (never smoker, smoker only before pregnancy, smoker only in early pregnancy, smoker throughout pregnancy), diet quality during pregnancy (measured by using the Probability of Adequate Nutrient intake based Diet quality index [PANDiet] score, which reflects nutrient-based reference guidelines adapted for pregnancy [START_REF] Bianchi | Pregnancy Requires Major Changes in the Quality of the Diet for Nutritional Adequacy: Simulations in the French and the United States Populations[END_REF] ), n-3 fatty acids supplementation during pregnancy (yes/no) and household income per consumption unit (≤€750, €751-1111, €1112-1500, €1501-1944, €1945-2500, >€2500 /month). The infant's characteristics included: child's sex, gestational age (in weeks), caesarean section (yes/no), birth weight, type of health professional for the first visit after delivery (paediatrician, other child doctor, general practitioner, other including emergency), age at introduction of infant formula (< 1 month, ≥1 month), presence of older siblings (no sibling, one sibling, at least 2 siblings) and family history of allergy (neither parent nor sibling with allergy; at least one parent or sibling with allergy including eczema, asthma, or hay fever).

Birth weight categories (small, adequate or large for gestational age) were defined according to the French Audipog reference curves [START_REF] Mamelle | Fetal growth from the AUDIPOG study. I. Establishment of reference curves[END_REF] .

Selection of the study sample

In the present study, we randomly selected one twin out of two (n=287) to avoid family clustering and excluded infants without data at the 2-month follow-up (n=1,696). Among the 11,244 infants who were not exclusively breastfed at the 2-month follow-up (eligible infants for this study), we excluded infants with insufficient detail on infant formula consumed at 2 months (n=667), those who had a medical diagnosis of cow's milk protein allergy reported at 2 months (n=193), and those with missing data for allergy or infection outcomes (n=1,995).

Finally, our study sample included 8,389 children (Figure 1).

Statistical analyses

Infants included in this study were compared to excluded eligible infants using chi-square test or Student t test as appropriate. Characteristics of the study population were presented according to LCPUFA enrichment of infant formula at the 2-month follow-up.

To provide national statistics on LCPUFA enrichment among formula-fed infants, data were weighted to account for inclusion procedure and biases related to non-consent [START_REF] Juillard | Weighting of Elfe survey data at time 0[END_REF] . Weighting also included calibration on margins from the state register's statistical data and the 2010

French National Perinatal study [START_REF] Blondel | Trends in perinatal health in France from 1995 to 2010. Results from the French National Perinatal Surveys[END_REF] on the following variables: age, region, marital status, migration status, level of education, and primiparity. This weighting was calculated for the subsample that completed the 2-month interview.

The associations between consumption of LCPUFA-enriched formula at age 2 months and the risk of infection or allergy were assessed with logistic regression models adjusted for potential confounders identified from the literature and validated with the directed acyclic graph method [START_REF] Shrier | Reducing bias through directed acyclic graphs[END_REF] . All analyses were adjusted for maternal characteristics (age, education level, city size, region of residence, employment status in early pregnancy, pre-pregnancy body mass index, smoking status, diet quality during pregnancy and n-3 PUFA supplementation during pregnancy), household income, family history of allergy, infant characteristics (sex, gestational age, caesarean section, birth weight category, type of health professional for the first visit after delivery, presence of older siblings and age at infant formula introduction) and variables related to study design (recruitment wave and maternity size). Missing covariates data were handled by using the Hot Deck method, which consists of replacing the missing value by that for respondents with the same characteristics [START_REF] Andridge | The Use of Sample Weights in Hot Deck Imputation[END_REF] .

In this study, interactions between consumption of LCPUFA-enriched formula and family or maternal history of allergy were not significant for any of the modelled outcomes. Therefore, we evaluated the associations in the whole sample.

Sensitivity analyses

To evaluate the potential presence of bias and test the robustness of our findings, several sensitivity analyses were performed. First, because infants consuming both breastmilk and infant formula may also benefit from LCPUFAs contained in breastmilk, analyses were conducted on the subsample of infants exclusively formula-fed at age 2 months (n=6,251).

Then, to limit the potential effect of the reverse causation bias, analyses were conducted on the subsample of infants who were fed the same formula during their first 2 months of life (n=4,429). Our analyses were also repeated among infants without any infection or allergic disease before age 2 months (n=6,490).

Because our previous findings showed that the consumption of partially hydrolyzed formulas at age 2 months was associated with increased risk of wheezing and food allergy [START_REF] Davisse-Paturet | Use of partially hydrolysed formula in infancy and incidence of eczema, respiratory symptoms or food allergies in toddlers from the ELFE cohort[END_REF] , we conducted a sensitivity analysis excluding all infants who consumed an infant formula with hydrolyzed proteins at age 2 months (n=364).

We also performed weighted analyses to assess the generalization of our results, using the weighting presented above.

In addition, instead of the Hot Deck method, we handled missing covariates data by using multiple imputations [START_REF] Sterne | Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls[END_REF][START_REF] Van Buuren | Multiple imputation of discrete and continuous data by fully conditional specification[END_REF] . We assumed that data were missing completely at random conditionally to the set of variables included in the analysis and generated 5 independent datasets with the fully conditional specification method (MI procedure, FCS statement, NIMPUTE option). We then calculated pooled effect estimates by using the SAS MIANALYSE procedure. Further details are available in Supplementary table 1.

All statistical analyses were conducted using SAS 9.4 (SAS Institute Inc., Cary, NC).

RESULTS

Sample characteristics

Of the 8,389 children formula-fed at 2 months of age that we included in this study, 4,282 were boys. As compared with eligible infants who were excluded from the analysis (n=2,855), those who were included were more often born with adequate weight for gestational age (80.1% vs. 78.4%; p=0.003), to an older mother (mean age 31± 4.9 vs. 29.4 ± 5.6 years; p<0.0001), with a higher education level (18.3% vs. 9.2% with a Master's degree, p<0.0001) and a higher income (1,663 ± 977 vs. 1,356 ± 1,219 euros per consumption unit; p<0.0001) and were more likely to stop smoking before pregnancy (24.7% vs. 18.5%; p<0.0001).

At age 2 months, 36.8% (CI=35.5-38.1) of formula-fed infants consumed DHA-enriched and 36.0% (CI=34.7-37.2) ARA-enriched formula. Since almost all of the DHA-enriched formulas were also enriched in ARA, we examined the effect of both DHA and ARA enrichment. Moreover, 11.3% (CI=10.4-12.2) of infants consumed a formula enriched in both DHA/ARA and EPA.

Table 1 presents the main characteristics of the study sample according to the type of formula consumed at age 2 months.

Enrichment of the formula consumed at age 2 months and risk of infection and allergy up to age 5.5 years

Consumption of DHA/ARA-enriched formula at 2 months of age was not associated with the occurrence of infection or allergy up to age 5.5 years (Tables 2 and3). As compared with non-DHA/ARA/EPA-enriched formula, the consumption of DHA/ARA/EPA-enriched formula was associated with a lower use of asthma medications up to age 3.5 years. We found no significant association between DHA/ARA/EPA enrichment and the occurrence of infections such as gastrointestinal infection, LRTI, URTI, nor with wheezing, itchy rash and food allergy up to age 5.5 years.

As compared with non-DHA/ARA-enriched formula, the consumption of formula with low DHA content (<10 mg/100 kcal) was associated with a higher risk of LRTI and wheezing, whereas an opposite trend was observed with consumption of formula with high DHA content (Tables 2 and3). We found no significant association between DHA enrichment level and the other outcomes.

As compared with non-DHA/ARA/EPA-enriched formula, the consumption of formula with high EPA content was associated with a lower risk of LRTI and a lower use of asthma medications. We found no significant association between EPA enrichment level and the occurrence of gastrointestinal infection, URTI, wheezing, itchy rash or food allergy up to age 5.5 years.

We found no significant association for formula content of linoleic acid and alpha-linolenic acid in any model (data not shown).

Sensitivity analyses

When restricting the analyses to exclusively formula-fed infants at age 2 months, the consumption of formula with high EPA content was also associated with a lower risk of LRTI, URTI and a lower use of asthma medications (Tables 2 and3).

In sensitivity analyses conducted to deal with the potential indication bias (including only infants without any change in formula up to age 2 months or without any infection or allergy up to age 2 months [Supplementary tables 2 and 3]), the negative association between EPA enrichment and asthma medications or LRTI remained similar. In addition, as compared with non-DHA/ARA-enriched formula, the consumption of formula with high DHA content (≥10 mg/100 kcal) was associated with a lower use of asthma medications (Supplementary table 3). In addition, among infants without any infection or allergy up to age 2 months, as compared with non-DHA/ARA/EPA-enriched formula, the consumption of DHA/ARA/EPAenriched formula was also associated with a lower risk of food allergy. Similar associations were observed for weighted analyses, analyses using multiple imputation to deal with missing data (Supplementary tables 4 and 5) and analyses among infants who did not consume infant formula with hydrolyzed proteins (data not shown).

DISCUSSION

In this observational study, one quarter of formula-fed infants consumed a formula enriched in DHA/ARA and just over one tenth consumed a formula additionally enriched in EPA. In the full sample, consumption of infant formula with low DHA content (<10 mg/100 kcal) was associated with a higher risk of LRTI and wheezing up to age 5.5 years, whereas an opposite trend was observed for high DHA content. The consumption of formula with high EPA content (≥3.2 mg/100 kcal) was associated with a lower risk of LRTI and use of asthma medications up to age 5.5 years. No significant association was found between LCPUFA enrichment and occurrence of gastrointestinal infection, URTI, itchy rash or food allergy up to age 5.5 years. Among exclusively formula-fed infants, the consumption of formula with high EPA content was associated with a lower risk of URTI.

Few studies have evaluated the association between the consumption of DHA/ARA-enriched formula and the risk of infection or allergy in childhood [START_REF] Foiles | Formula with long-chain polyunsaturated fatty acids reduces incidence of allergy in early childhood[END_REF][START_REF] Birch | The impact of early nutrition on incidence of allergic manifestations and common respiratory illnesses in children[END_REF][START_REF] Lapillonne | Infants fed formula with added long chain polyunsaturated fatty acids have reduced incidence of respiratory illnesses and diarrhea during the first year of life[END_REF][START_REF] Pastor | Infants fed docosahexaenoic acid-and arachidonic acid-supplemented formula have decreased incidence of bronchiolitis/bronchitis the first year of life[END_REF] , and none has evaluated the effect of additional enrichment in EPA, which makes comparing our results difficult. In addition, the outcomes definition was quite heterogeneous across the studies, which suggests caution in comparing results.

Available studies suggest that consumption of DHA/ARA-enriched formula (17 mg /100 kcal for DHA and 34 mg /100 kcal for ARA) was associated with a lower incidence or delayed onset of bronchitis/bronchiolitis [START_REF] Lapillonne | Infants fed formula with added long chain polyunsaturated fatty acids have reduced incidence of respiratory illnesses and diarrhea during the first year of life[END_REF][START_REF] Pastor | Infants fed docosahexaenoic acid-and arachidonic acid-supplemented formula have decreased incidence of bronchiolitis/bronchitis the first year of life[END_REF] . A similar trend was observed in the present study among infants who consumed formula with high DHA content (10-20 mg/100 kcal), but the association did not reach statistical significance (OR=0.88; 95% CI=0.76-1.02 for LRTI).

Moreover, an additional enrichment of formula with high EPA content was associated with a lower risk of LRTI. As infant formula with high EPA content also had high DHA content in this study (mean DHA enrichment level=14.7 mg/100 kcal), the observed findings suggest that formula with DHA content close to the recommendations would have a better effect when it is additionally enriched with EPA.

Regarding URTI, we found no significant association with the consumption of DHA/ARAenriched formula. Currently available studies showed inconsistent findings, with one RCT which highlights a protective effect of DHA/ARA enrichment on risk and delayed onset of URTI [START_REF] Birch | The impact of early nutrition on incidence of allergic manifestations and common respiratory illnesses in children[END_REF] , and another that showed no significant association with URTI incidence [START_REF] Pastor | Infants fed docosahexaenoic acid-and arachidonic acid-supplemented formula have decreased incidence of bronchiolitis/bronchitis the first year of life[END_REF] .

Although a previous observational study [START_REF] Lapillonne | Infants fed formula with added long chain polyunsaturated fatty acids have reduced incidence of respiratory illnesses and diarrhea during the first year of life[END_REF] found that DHA/ARA enrichment of infant formula was associated with a lower risk and delayed onset of diarrhea requiring medical attention, we were not able to highlight any association between this enrichment and gastrointestinal infections. This lack of association could be due to the cases definition (diarrhea requiring medical attention vs. hospitalization for gastroenteritis/dehydration or emergency consultation for diarrhea/vomiting/dehydration) and the type of outcome investigated (gastrointestinal infections [yes/no] and not age at first episode).

Finally, in agreement with the DIAMOND study [START_REF] Foiles | Formula with long-chain polyunsaturated fatty acids reduces incidence of allergy in early childhood[END_REF] , the consumption of DHA/ARA-enriched formula was not associated with the risk of itchy rash. However, a protective effect of DHA/ARA enrichment was found for allergic illness or skin allergic illnesses in children without a maternal history of allergy [START_REF] Foiles | Formula with long-chain polyunsaturated fatty acids reduces incidence of allergy in early childhood[END_REF] . One observational study also found that consumption of DHA/ARA-enriched formula was associated with a lower risk of having an increased number of eczema episodes, but not with its incidence [START_REF] Lapillonne | Infants fed formula with added long chain polyunsaturated fatty acids have reduced incidence of respiratory illnesses and diarrhea during the first year of life[END_REF] .

Although the literature focused mainly on DHA, our findings suggest that EPA could also play an important role. Indeed, DHA and EPA are both n-3 LCPUFAs and are important for growth, development and many aspects of health including cell function and inflammation 48- 52 .

The new European regulation now requires the enrichment of all infant and follow-on formula in DHA, whereas the enrichment in ARA and EPA remains optional. The European Academy of Pediatrics further suggested that DHA enrichment had to be combined with ARA enrichment [START_REF] Koletzko | Should formula for infants provide arachidonic acid along with DHA? A position paper of the European Academy of Paediatrics and the Child Health Foundation[END_REF] . Our results suggest that EPA enrichment could also be of great interest, as even with sufficient ALA in infant formula, the conversion capacity of ALA to EPA is very limited [START_REF] Swanson | Omega-3 fatty acids EPA and DHA: health benefits throughout life[END_REF] . However, further well-designed studies are needed to examine the potential benefits of such enrichment.

The ELFE cohort is a nationwide study of birth in 2011 in metropolitan France (excluding very preterm babies). The detailed data on infant diet allowed for examining precisely, under real conditions of use, the effect of LCPUFA enrichment on infection and allergy in childhood. We have to acknowledge that health data were reported by parents and not validated by medical records, which could lead to a potential measurement bias. However, the items used were derived from international ones [START_REF] Asher | Worldwide time trends in the prevalence of symptoms of asthma, allergic rhinoconjunctivitis, and eczema in childhood: ISAAC Phases One and Three repeat multicountry cross-sectional surveys[END_REF] to limit this bias. Moreover, the prospective design and the availability of repeated data on health-related outcomes in childhood limited memory bias. The very large sample and the collection of detailed socio-demographic or economic data ensured large statistical power and favoured control for potential confounders.

Furthermore, the interpretation of our findings should be done with caution as an arbitrary cut-off point (median LCPUFA content of formula) was used to classify infant formula consumed in this study as "low LCPUFA content" or "high LCPUFA content". Indeed, this methodological approach was used since very few infants consumed infant formula with the minimum dose used on the available RCTs (17 mg /100 kcal for DHA, n=75 infants) or that recommended by the European regulation (20 mg /100 kcal for DHA, n=1 infant). Further studies are needed to confirm the present results with less arbitrary cut-offs. We also performed several sensitivity analyses to examine the robustness of our results and to address selection and reverse causation bias. As findings remained stable across these sensitivity analyses, these potential biases, if any, should have a limited impact on our findings.

In conclusion, this observational study under real life conditions suggests that LCPUFAenriched formula, especially those with additional enrichment with high EPA level was associated with a lower risk of LRTI, URTI and use of asthma medications. However, URTI, findings were inconsistent across analyses. Replication in other observational studies and larger well-designed RCTs using the recommended levels are needed to confirm these results.

(

  INED) and the French National Institute of Health and Medical Research (INSERM), in partnership with the French blood transfusion service (Etablissement français du sang [EFS]), Santé publique France, the National Institute for Statistics and Economic Studies (INSEE), the Direction générale de la santé (DGS; part of the Ministry of Health and Social Affairs), the Direction générale de la prévention des risques (DGPR; Ministry for the Environment), the Direction de la recherche, des études, de l'évaluation et des statistiques (DREES, Ministry of Health and Social Affairs), the Département des études, de la prospective et des statistiques (DEPS; Ministry of Culture), and the Caisse nationale des allocations familiales (CNAF), with the support of the Ministry of Higher Education and Research and the Institut national de la jeunesse et de l'éducation populaire (INJEP). Via the RECONAI platform, it receives a government grant managed by the National Research Agency under the "Investissements d'avenir" programme (ANR-11-EQPX-0038).

  informed about their rights to oppose it. The ELFE study was approved by the Advisory Committee for the Treatment of Health Research Information (Comité Consultatif sur le Traitement de l'Information en matière de Recherche en Santé [CCTIRS]), the National Agency Regulating Data Protection (Commission National de l'Informatique et des Libertés [CNIL]) and the National Council on Statistical Information (Conseil National de l'information Statistique [CNIS]).

Acknowledgments:

We thank the scientific coordinators (B Geay, H Léridon, C Bois, JL Lanoé, X Thierry, C Zaros), IT and data managers, statisticians (M Cheminat, C Ricourt, A Candea, S de Visme), administrative and family communication staff, study technicians (C Guevel, M Zoubiri, L G L Gravier, I, Milan, R Popa) of the ELFE coordination team, as well as the families that gave their time for the study.

Conflict of interest:

Amandine Divaret-Chauveau has been invited to speak at conferences organized by laboratories such as Nutricia, Sodilac, Mead Johnson, outside of the submitted work. Karine Adel-Patient has been invited to speak at a conference organized by Stallergènes, outside of the submitted work. None of the other authors declare any conflicts of interest.


F Harlé

F Chatelain

C Gouy-Pailler

S Achard

RANK-BASED MULTIPLE CHANGE-POINT DETECTION IN MULTIVARIATE TIME SERIES

Keywords: Rank statistics, joint segmentation, dependency structure learning, Bayesian inference, MCMC methods, Gibbs sampling

In this paper, we propose a Bayesian approach for multivariate time series segmentation. A robust non-parametric test, based on rank statistics, is derived in a Bayesian framework to yield robust distribution-independent segmentations of piecewise constant multivariate time series for which mutual dependencies are unknown. By modelling rank-test p-values, a pseudo-likelihood is proposed to favour change-points detection for significant p-values. A vague prior is chosen for dependency structure between time series, and a MCMC method is applied to the resulting posterior distribution. The Gibbs sampling strategy makes the method computationally efficient. The algorithm is illustrated on simulated and real signals in two practical settings. It is demonstrated that change-points are robustly detected and localized, through implicit dependency structure learning or explicit structural prior introduction.

INTRODUCTION

Detecting change-points in multivariate data is a classical but major issue in many fields, like bioinformatics, industrial monitoring, finance, or genomics. The literature in signal processing suggests different approaches for the multiple change-point detection problem. Existing monovariate methods (as presented in [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF][START_REF] Harchaoui | Kernel-based methods for hypothesis testing: A unified view[END_REF][START_REF] Zou | Empirical likelihood ratio test for the change-point problem[END_REF]) can be extended to the multivariate case. A fused lasso latent feature model has been proposed, for example in genomics [START_REF] Nowak | A fused lasso latent feature model for analyzing multisample aCGH data[END_REF]. Other approaches introduce a hierarchical Bayesian model [START_REF] Dobigeon | Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a bayesian sampling approach[END_REF]. Our framework is the off-line joint segmentation of piecewise constant time series: the events are abrupt changes in the distribution of the observations, their number and localizations are unknown.

We consider a matrix X of K time series of N points: x j,i is the observation of sensor j at time i, and is a changepoint if a change occurs at t ∈]i, i + 1], then x j,i is the last point of a segment. The nature of these points is represented by an indicator variable r j,i , element of the matrix R:

r j,i = 1 if x j,i is a change-point, 0 otherwise, (1)
with the convention that r j,1 = r j,N = 1, ∀j ∈ {1, ..., K}. R denotes the segmentation of X by the change-points in each signal, it is the parameter to estimate. Another parameter of interest is the implicit dependency structure that express the links between the signals, assuming that if a change-point appears in a signal j, it occurs simultaneously on the signals that depend of j. Then the probabilities that a change-point is shared across several signals is represented in the vector P . Following the Bayes' rule, we get the posterior for R and P :

f (R, P |X) ∝ L * (X|R)f (R|P )f (P ). (2)
We propose to combine a non-parametric test to Bayesian inference to express this statistical model. We chose the Wilcoxon signed-rank test, that is often presented as robust to several non-normal distributions, and makes the model free of strong prior. Its power is discussed in [START_REF] Lehmann | The power of rank tests[END_REF], the test is shown to be optimum for various alternatives in hypothesis testing. This paper is organised as follows. The methodology is described in section 2. First, we explain how the Wilcoxon test is applied to define the pseudo-likelihood L * (X|R). Then the choice of priors for R and P is detailed to model the dependency structure between time series. Finally a marginalized posterior density for the parameter R is expressed, and the procedure to estimate the maximum a posteriori (MAP), with a Markov Chain Monte Carlo (MCMC) method and a Gibbs sampling strategy, is given. In section 3, we present the results of the application on synthetic and real data. The discussion about the model is in section 4.

BERNOULLI DETECTOR MODEL

The parameters to estimate are the change-points configurations at each time index i ∈ {1, ..., N }: they are modelled by the vectors of indicators defined in (1): R i = (r 1,i , ..., r K,i ) T , for signals 1 to K, where configuration R i ∈ {0, 1} K . This leads to consider the parameter R = (R 1 , ..., R N ) ∈ {0, 1} K×N , which is inferred within a Bayesian framework.

Here the data term is derived from the rank statistics, and a prior on P is chosen to infer the dependencies between signals. All observations in X are assumed to be mutually independent, but have similar statistical properties within a segment, like the median. In the following parts, we derive the joint probability distribution for R, P and other hyperparameters given the data X.

Change-point model

Each change-point probability is computed by the Wilcoxon rank sum (aka Mann-Whitney) test. In [START_REF] Lung-Yut-Fong | Homogeneity and change-point detection tests for multivariate data using rank statistics[END_REF], the authors present a statistic inspired by this test on multivariate time series to detect multiple change-points, but they are supposed to occur simultaneously on all time series, assuming a fully connected structure. In our model, the Wilcoxon test is applied on each observation x j,i , for a given segmentation in R, between the two samples defined by the segments s 1 = {x j,i -+1 , ..., x j,i }, of length n 1 , and s 2 = {x j,i+1 , ..., x j,i + }, of length n 2 , where i -and i + denotes the previous and next change positions in signal j. The Wilcoxon statistic is defined from the sum of the ranks, denoted R 1 and R 2 for the observations of s 1 and s 2 respectively, in the global segment s = (s 1 , s 2 ):

U = min(U 1 , U 2 ) with U k = n k n l + n k (n k + 1) 2 -R k , (3)
where (k, l) = (1, 2) or (2, 1). For small n 1 and n 2 , the resulting p-values are tabulated, otherwise they are computed from a normal approximation of a standardized value of U . It yields a statistical test to reject the null hypothesis H 0 that the difference between the observations of each segment is symmetrically distributed around zero; this is for instance the case when the samples of both segments are identically distributed.

The key-point of our approach is to define the observation model on the p-values statistics derived from the Wilcoxon test and considered as random variables, rather than on the samples X or the statistic U . A classical result (see [START_REF] Sackrowitz | P values as random variables-expected p values[END_REF] for instance) is that the p-value is uniformly distributed on [0, 1] under the null hypothesis H 0 , i.e. when x j,i is not a changepoint. We have to specify now a model for the alternative hypothesis H 1 , i.e. when x j,i is a change-point. In this case, the p-value can be viewed as an outlier with respect to its null distribution. Several authors proposed some distributions or approximations for some alternatives hypotheses, see for instance [START_REF] Bahadur | Simultaneous comparison of the optimum and sign tests of a normal mean[END_REF][START_REF] Becker | Small-sample accuracy of approximate distributions of functions of observed probabilities from t tests[END_REF][START_REF] Sellke | Calibration of p Values for Testing Precise Null Hypotheses[END_REF][START_REF] Allison | A mixture model approach for the analysis of microarray gene expression data[END_REF]. To be precise, the alternative hypothesis must support the largest values of the test statistic U , i.e. the smallest p-values. This leads to consider families of distributions whose densities are decreasing in p. Inspired by [START_REF] Sellke | Calibration of p Values for Testing Precise Null Hypotheses[END_REF], we assume that the p-value is distributed under H 1 according to a Beta B(γ, 1) distribution where the parameter γ is in [0, 1]. This is a specific case of a general class of distribution for the choice of the distribution of p-values under the alternative hypothesis, for which a lower bound on Bayes factor is given [START_REF] Sellke | Calibration of p Values for Testing Precise Null Hypotheses[END_REF]. It yields the following density:

f (p j,i |R) = 1 [0,1] (p j,i ) if r j,i = 0 (H 0 ), γp γ-1 j,i 1 [0,1] (p j,i ) if r j,i = 1 (H 1 ). (4)
Note that when γ = 1, the alternative distribution reduces to the uniform distribution, i.e. the null distribution. Another example of use of a Beta-Uniform model can be found in [START_REF] Allison | A mixture model approach for the analysis of microarray gene expression data[END_REF]. In this work, we propose to calibrate the alternative distribution with respect to a frequentist significance level. Then, the Beta parameter γ is defined after the acceptance level α ∈ [0, e -1 ] (for brevity reason, the notation γ(α) is omitted) such that γ is the unique solution on [0, 1] of the following equation

f (α|r = 1) = 1 ⇔ γα γ-1 = 1. ( 5
)
This calibration is discussed in more details in section 2.3. Finally, the Wilcoxon test is applied on each observation x j,i , and following (4), the data term for our model is formed by the inference function

L * (X|R) = K j=1 N i=1 f (p j,i |R) = K j=1 N i=1 γp γ-1 j,i rj,i . (6)
It is important to note that this inference function is a composite marginal likelihood based on the univariate distributions of the (p j,i ) 1≤j≤K,1≤i≤N . In particular, the dependencies between the (p j,i ) 1≤i≤N in signal j are not taken into account. As a consequence, this is not a proper likelihood as defined in [START_REF] Monahan | Proper likelihoods for bayesian analysis[END_REF] and the coverage probabilities induced by this model should differ from the real ones. However L * (X|R) depends on the data through the marginal distributions of the p-values of the Wilcoxon statistics and is calibrated up to a significance level α. The impact this level on the detection will be discussed in a future work. As a consequence, this inference function acts like a data term, and can be used as a pseudo-likelihood in the Bayesian framework. We refer to this model as the Bernoulli detector model. It differs from classical Bernoulli Gaussian models where the likelihood is directly derived from the parametric Gaussian assumption of the samples.

Prior on configurations R i

The vectors (R i ) 1≤i≤N are assumed to be a priori independent. Thus the prior distribution can be written:

f (R) = N i=1 f (R i ). (7)
Following the approach presented in [START_REF] Dobigeon | Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a bayesian sampling approach[END_REF], where vague prior on the dependency structure is chosen, the vector P = (P ) ∈E represents the dependency structure. Here P denotes the probability of having R i = , where is a K × 1 vector of zeros and ones, called a configuration, and E is the subset of {0, 1} K of all configurations for R i . This notation means that if the probability P is high, is more likely to appear in R, then changes tend to be simultaneous in all signals j such that (j) = 1. Then [START_REF] Lung-Yut-Fong | Homogeneity and change-point detection tests for multivariate data using rank statistics[END_REF] becomes

f (R|P ) = ∈E P S (R) (8)
where S (R) is the number of times that is found in R. The parameter P follows a Dirichlet distribution:

P |d ∼ D L (d) (9)
with the hyperparameter vector d = (d ) ∈E and where L is the cardinal of E. As in [START_REF] Dobigeon | Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a bayesian sampling approach[END_REF], all the d are set to the same deterministic value d ≡ d = 1, so distribution ( 9) is uniform.

Posterior distribution

The unnormalised posterior distribution (2) of the changepoints indicators and the hyperparameters expresses as follows:

f (R, P |X) ∝   K j=1 N i=1 (γp γ-1 j,i ) rj,i   ∈E P S (R)+d -1 .
(10) Based on this posterior, it is now possible to express some properties of the classical Bayesian estimates of a changepoint for a given location (j, i) under some simple hypothesis. We assume in the following proposition that there is no other change-point in the signals. We denote as R \(j,i) 0 the void configuration event such that r l,k = 0 for all (l, k) = (j, i), and 0 and 1 stand for the configurations of the column R j when r j,i = 0 and r j,i = 1 respectively. Proposition 2.1 MAP and MMSE estimators given R \(j,i) 0 , P 0 and P 1 . Under the previous hypothesis, if P 0 = P 1 and for a chosen significance level α:

• the conditional MAP estimate is 1 iff the p-value is lower than α, according to (4) and (5),

• the conditional MMSE estimate is larger than 1/2 iff the p-value is lower than α.

The proofs are directly derived from the definition of these estimators and the expression of the posterior distribution [START_REF] Becker | Small-sample accuracy of approximate distributions of functions of observed probabilities from t tests[END_REF]. These properties illustrate the influence of the significance level α chosen in (5) to calibrate the distribution (4) under H 1 for the single change-point problem. If the priors on the configurations are equivalent then the presence of a change-point is favoured when the support against the null hypothesis is significant for the level α.

Finally, based on [START_REF] Becker | Small-sample accuracy of approximate distributions of functions of observed probabilities from t tests[END_REF], one can see that the posterior of the hyperparameter P reduces to a Dirichlet distribution

P |R, d ∼ D L (S (R) + d ), (11)
in agreement with the equation ( 22) of [START_REF] Dobigeon | Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a bayesian sampling approach[END_REF]. Thus, the P can be easily integrated out in [START_REF] Becker | Small-sample accuracy of approximate distributions of functions of observed probabilities from t tests[END_REF] as nuisance parameters, yielding the following marginalized posterior

f (R|X) ∝   K j=1 N i=1 (γp γ-1 j,i ) rj,i   × ∈E Γ(S (R) + d ) Γ(N + L) . (12)

Algorithm

The problem of estimating R is solved by a MCMC method, where samples are drawn according to distribution [START_REF] Allison | A mixture model approach for the analysis of microarray gene expression data[END_REF]. At each step, all vectors R i are simulated, following a Gibbs sampling strategy. At each modification of r j,i , the segmentation changes and the p-values of the previous and the next change-point in the signal j have to be updated. However an approximation is done the reduce the number of steps in the algorithm, and only the current p-value p j,i is computed for the new segmentation. For brevity reason, the empirical validation is not presented here. The main steps of the algorithm are detailed below in Alg. 1, for M MCMC iterations.

Algorithm 1: Bernoulli detector require E = { 0 , ..., l , ..., L } ⊂ {0, 1} K , α initialize R (0) , S (R (0) ) for m ← 1 to M do initialize the index set I = {1, . . . , N } while I = ∅ do pick randomly i in I compute (p (m) j,i ) 1≤j≤K sample R (m) i
according to [START_REF] Allison | A mixture model approach for the analysis of microarray gene expression data[END_REF] remove i from I Optional step: sample P (m) from its posterior (11) return R

APPLICATION ON DATA

Simulation

At first, the impact of the noise level and the presence of outliers is studied. Two signals of N = 100 points are generated, with a change-point at t = 50, that defines two segments in each signal. In the first case, the observations on each segment k follow a non-standardized Student's t-distribution with a heavy tail, whose parameters are (ν, µ k , σ) where ν = 3.0, and σ 2 = ν/(ν -2). In the second case, the distribution is normal N (µ k , σ). The SNR is defined as SNR = 10 log (µ0-µ1) 2 σ 2

. To validate the non-parametric approach of our model, we compare it with a classical Bernoulli Gaussian model. The performances are computed on the MAP estimator in terms of precision, that is the proportion of true change-points found in the detected ones, with a tolerance of ±t points in time. The results are shown in figure 1. It appears in 1(a) that the Bernoulli detector model has much less false positives than the Bernoulli Gaussian, as expected because of robustness of the rank-test to distributions with heavy tail in the non-parametric pseudo-likelihood. The generic nature our model is pointed out by the fact that the precision values are the same in the normal case (1(b)). It is equivalent then to the Bernoulli Gaussian model, and slightly better for small SNR. Due to the restricted length of the paper, the recall values are not plotted, however, we found results in agreement to precision values: the recall values are converging to one with high SNR and the Bernoulli detector model presents better recall performance for a Student noise. In the second test we consider five signals following Normal distributions, classified in two independent groups: group 1 with signals 1 and 2 (3 segments), and group 2 with signals 3, 4 and 5 (4 segments). The change-points are the same for all signals within a group. The SNR defined previously is 14.0 for signals 1 to 4 and -6.0 for signal 5, where changepoints are not visible (figure 2). The data are processed independently (results in blue in figure 2) and jointly with a non-informative prior on dependency structure (in red in figure 2). Despite a higher computational cost, the advantage of the of joint estimation is obvious in the results: the dependency structure of both groups is perfectly learned, and all change-points are precisely found, especially in signal 5, whereas the estimation fails with independent processing.

BD, t =0 BD, t =1 BD, t =5 BG, t =0 BG, t =1 BG, t =5 (a) Student's t-distribution
BD, t =0 BD, t =1 BD, t =5 BG, t =0 BG, t =1 BG, t =5 (b) Normal distribution

Household electric power consumption data

In this section two principled applications of the proposed algorithm are illustrated using measurements of household power consumption [START_REF] Bache | UCI machine learning repository[END_REF]. This real dataset consists of four time series. One of them depicts the global electrical energy consumption in the house, while the others are devoted to the measurements of power demand by specific devices. Learning structure from data. In a first setting, the algorithm is applied with non-informative priors on the dependency structure between the four signals, namely all configurations in {0, 1} K are considered. Thanks to the Bayesian framework, the Gibbs sampling process can be swimmingly supplemented to yield a posteriori probabilities on the configurations. Resulting histograms of the links between the global household energy consumption, denoted 1, and the three sub-metering measurements (respectively denoted 2, 3 and 4), plus the link between 3 and 4, are depicted in figure 3. Unsurprisingly these figures confirm that, when a changepoint is observed in signals 2, 3 or 4, it is likely to be observed in signal 1, whereas events appear independently between sub-metering signals. Therefore this setting amounts to building a weighted graph of the dependencies between time series by focusing on small-scale events such as changepoints. Visual results on a small portion of signals are represented in figure 4(a).

Using a priori data structure to alleviate data processing.

In a second setting the introduction of informative priors on the structure of the dependencies is considered. It is indeed known that strong additive relationships exist between signal 1 and signals 2, 3 and 4 respectively. The change-point detection algorithm is thus modified by restricting E such that change-points in signals 2, 3 or 4 must occur simultaneously in signal 1. Consequences of such an approach are two-fold. First computations are made significantly faster thanks to the sharp decrease of the number of considered configurations. Second spurious or missing change-points can be filtered out with this simple approach. Visual results are represented in figure 4(b).

DISCUSSION AND CONCLUSION

An innovative change-point detection algorithm has been proposed. Its strength relies mainly on two characteristics. First the algorithm is built upon robust foundations, namely rankbased statistics. This confers the approach with a stable behaviour for various signal distributions, a worthy robustness to noise outliers and a strong efficiency in small samples problems. Second its surrounding Bayesian formulation provides the algorithm with awesome flexibility regarding the problems at hand. While classical multivariate approaches focus on detecting simultaneous change-points across time series, the presented framework offers a wide range of possible usages, from joint independent detections across time series, to simultaneous change-points detection through the introduction of priors. Two main settings have been presented to illustrate the use of the algorithm. On the one hand it has been demonstrated that the approach can be successfully applied to infer dependency structure from change-points occurrences. This constitutes an innovative way of estimating the dependencies between signals, and will be explored in future work. This algorithm therefore provides a tool to analyse multivariate time series behaviour through the observation of small-scale events in the signal. On the other hand, it has been demonstrated that priors on the dependency structure can be favourably introduced to improve change-point detections, and also reduce the expensive computational cost, due to MCMC sampling and configuration testing.

Fig. 1 .

 1 Fig. 1. Precision for several SNR, of Bernoulli detector model (BD, red) and Bernoulli Gaussian model (BG, blue), with tolerance t.

5 Fig. 2 .

 52 Fig. 2. At top, independent MAP estimation of each signal. At bottom, joint estimation with a non-informative prior on the changepoints structure. Real change-points are represented in dashed lines.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Posterior conditional probabilities Pr(i|j) to have a changepoint in signal i when there is a change-point in signal j, for noninformative prior.