Amélie Crépet

Jessica Tressou

February

Bayesian nonparametric model for clustering individual co-exposure to pesticides found in the French diet

Keywords: Dirichlet process, Bayesian nonparametric modeling, multivariate Normal mixtures, clustering, multivariate exposure, food risk analysis

Bayesian nonparametric model for clustering individual

co-exposure to pesticides found in the French diet.

Amélie Crépet, Jessica Tressou

Introduction

Each food product may contain several pesticide residues, consequently meals ingested daily may include a wide range of pesticides. All consumers are thus exposed to complex cocktails of pesticides 1 whose combined effects on health are still unknown. This work proposes an original methodology to respond to the following question "what are the cocktails of pesticides to which the French population is simultaneously and most heavily exposed?" Cocktails of pesticides were selected based on their joint probability of occurring at high doses in the French diet, according to the following process. The French population exposure to P different pesticides found in the diet, called the co-exposure, was estimated considering both the food residue level patterns generated from national food monitoring administrations in charge of food control and the dietary habits of n surveyed individuals from the INCA 2 French consumption survey [START_REF] Afssa | Report of the Individual and the National Study on Food Consumption[END_REF]. A nonparametric mixture model was developed to cluster the population co-exposures in order to define groups of individuals with similar patterns of exposures to the P pesticides. For the biggest clusters made of the highly exposed individuals, the correlations between exposures to the P pesticides were studied to characterize the relevant cocktails.

To define homogeneous groups of individuals, the population co-exposure to the P pesticides was modeled with a Bayesian nonparametric model relying on the use of the Dirichlet process, [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF][START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF]. This approach consists in building a mixture model, in which the number of mixture components is potentially unlimited, and is itself a random variable that is part of the overall model. In our pesticide study, a multivariate Normal distribution was chosen as the kernel density of the mixture. Thus, the correlations of the P pesticides were modeled and individuals were clustered according to both their co-exposure levels and their co-exposure correlations. The mixing distribution was modeled with a Dirichlet process (DP) which is the most popular choice of prior distribution for infinite mixture models in a Bayesian nonparametric context, [START_REF] Lo | On a class of Bayesian nonparametric estimates: I. Density Estimates[END_REF].

Indeed, the DP can be viewed as a probabilistic measure on the space of probability measures which has the required properties defined by [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF], to be used as a prior distribution.

Models combining both mixture model and Dirichlet process are called Dirichlet Process Mixture (DPM) models. This approach provides a way of putting a distribution on all possible partitions of the data and then, through the corresponding posterior distribution or classification likelihood, selecting the most likely cluster formations. To include the uncertainty of individual exposure to each pesticide, we modified the base model using a hierarchical DP approach, similar to the one proposed by [START_REF] Teh | Hierarchical Dirichlet Processes[END_REF]. In the case of a real application to a large sample size, like the one under study, scanning all possibilities is computationally unfeasible and thus there is a need to resort to simulation methods such as Markov Chain Monte Carlo techniques. Popular choices for such simulation technique include Gibbs samplers based either on the Pólya urn scheme which is closely related to the Chinese restaurant process (CRP, [START_REF] Blackwell | Ferguson distributions via Polya urn schemes[END_REF][START_REF] Pitman | Some developments of the Blackwell-MacQueen Urn scheme[END_REF], or on the stick-breaking (SB) representation of the DP [START_REF] Sethuraman | A constructive definition of Dirichlet prior[END_REF]. The latter was retained to account for the complexity of our model's hierarchy within an effective algorithm, [START_REF] Ishwaran | Gibbs Sampling Methods for Stick-Breaking Priors[END_REF]. The stick-breaking priors can be simply constructed using a sequence of independent Beta random variables and the SB algorithm does not require individuals to be reassigned one by one like the CRP, which becomes computationally complex with a large sample size. Still, to reduce this computational burden further, the random-block sampling proposed in the technical report of the Department of Informations Systems, Business Statistics and Operations Management, Hong Kong University of Science and Technology [START_REF] Cabrera | Random Block Sampling for high dimensional clustering (from the Bayesian point of view)[END_REF], for the CRP was applied to the SB algorithm to reduce computational time. This procedure consists in subsampling d variates among the available P dimensions at each iteration.

The following section describes the data on residue levels, quantities consumed and the coexposure estimate. In Section 3, the Dirichlet Process Mixture model applied to the modeling of pesticide co-exposure is outlined. Then, the chosen prior distributions and the model in its hierarchical form are presented. For reproducibility purposes, the SB and the random-block SB algorithms are detailed. Finally, in section 4 the models are applied to a set of simulated data and to the French population's co-exposure to pesticides.

Co-exposure to pesticides

Pesticide food exposure was estimated with the individual food consumption data from the French dietary survey and with residue levels obtained from French pesticide residue monitoring programs.

It is necessary to identify and take into account all foodstuffs in which significant residues might occur, as well as all pesticides that may be present in the food. Therefore, a first step consists in identifying the food / pesticide combinations to consider for the exposure assessment.

Data on pesticide residues levels

The data source on pesticide residues in food and drinking water corresponds to the annual monitoring programmes implemented in 2006 by the French administrations (Ministry of Economy, Ministry of Agriculture, Ministry of Health). These surveys provide sample distributions of residues for up to 300 pesticides measured in about 150 raw agricultural commodities (RACs). The number of samples varies from about 10 for minor commodities up to 480 for staples (apple, lettuce, etc.).

Most of analytical results are left-censored, i. e. the residue level stands between 0 and the limit of reporting (LOR) from the laboratory but cannot be quantified or detected. For one pesticide, when more than 90% of results were censored in each food commodity then, it was considered of no interest to take it into account for the co-exposure calculation. Some pesticides have however been included in the study, when the determined residue levels were of the same order as the corresponding LOR. In such cases, it was considered that the pesticide may really be present but could not be determined due to analytical limitations. Thus, 79 pesticides were selected for the analysis.

Residues of the selected pesticides were analyzed in 120 RACs and in drinking water consumed by the INCA2 population. A total of 306, 899 analytical results corresponding to 8, 364 combinations of pesticide/commodity were used in this work.

Food consumption data

Consumption data were provided by the second "Individual and National Study on Food Consumption", INCA2 survey, carried out by the French Food Safety Agency, AFSSA (2009). The study was conducted in three fieldwork phases between late 2005 and April 2007 in order to cover seasonal variations. Two independent population groups were included in the study: 2, 624 adults aged 18-79 years and 1, 455 children aged 3-17 years. Participants were selected using a three-stage random probability design stratified by region of residence, size of urban area and population category (adults or children). Each subject was asked to complete a seven-day food diary as well as other questionnaires on anthropometric and socio-economic factors. Food were subsequently categorized into 1, 305 "as consumed" food items (INCA2 classification). In order to match the consumption data to the residues ones, the food items defined in the INCA2 survey were broken down into 181 RACs. To do this, 763 standardized recipes were used, which were defined by the French Food Safety Agency taking account of industrial processes, home cooking habits and edible portions for the INCA2 survey [START_REF] Afssa | Report of the Individual and the National Study on Food Consumption[END_REF].

In the INCA2 survey, sampling weights are provided for each surveyed individual representing its frequency in the entire French population. Based on the sampling weights provided for each individual, two samples of adults and children were built out of the original samples by carrying out random trials with replacement. Only normal-reporters, i.e. individuals whose energy needs are covered by the declared consumptions, were considered for this study. Therefore, two normalreporter samples of 1, 898 adults and 1, 439 children were used for the analysis.

Dietary co-exposure assessment

To estimate acute exposure, i.e the exposure during 24-hour day, one day of consumption was randomly selected for each individual from the 7 days recorded in the INCA2 survey. The individual daily consumption of a commodity denoted c ia corresponds to the sum of all the quantities of commodity a consumed by the individual i during the selected day. For each commodity a treated with the pesticide p, the daily consumption c ia was multiplied by one residue level q pa and adjusted by the body weight w i of the consumer i. The intakes calculated for each commodity were summed to obtain a total daily exposure in milligrams of the pesticide p per kilogram of body weight of the consumer per day (mg/kg bw/d). This process was performed for m = 1, . . . , M values randomly selected from the contamination distribution of each pesticide/commodity combination to account for the residue level uncertainty. The final data set consisted of a series of M possible daily exposures x pim = Ap a=1 (c ia × q pam )/w i to each pesticide p = 1, . . . , P , for each individual i = 1, . . . , n.

In order to deal with quantitative values, each censored datum was uniformly selected between 0 and its censoring value (LOR). Similarly, for each pesticide/commodity combination, random contamination values were uniformly selected among the different residues levels, with respect to the probability of lying in the interval between two consecutive observed residues levels. Scaling problems between pesticide exposure levels were ruled out by using a log scale, centering around the mean and rescaling by the standard deviation for each pesticide and across individuals.

Then, two datasets were created, one considering the 95 th percentile of the distribution of the M exposures to the pesticide p of each individual i (one high exposure per individual), the other one considering the entire distribution empirically described by the M exposure values. High percentiles such as the 95 th , 99 th or the maximum of the exposure distribution are usually employed to study the worst case of exposure to a chemical. In the context of co-exposure assessment, this scenario could be quite unrealistic in the sens that along a day the probability of being jointly exposed to high values for all pesticides is very low. Random co-exposures are more likely to occur in current life. In this way, using the distribution of exposure for each pesticide permits to attribute to each individual random values of exposure integrating the uncertainty of residue levels of each pesticide. Therefore, computations were carried out, in the first case with a co-exposure matrix of size n × P and in the second case with an array of size n × P × M .

Sensitivity of the exposure estimate to the value of M was tested in comparing the distributions of the exposure obtained with M = 100 and M = 1000 using a Wilcoxon test. For the alternative hypothesis "two sided", the p-value ranges between 0.80 and 0.96 depending on the pesticide and for the alternative hypothesis "greater", the p-value ranges between 0.40 and 0.55. So no statistical difference between both distributions were observed, and therefore computations were performed with M = 100.

Methodology

Bayesian nonparametric model-based clustering

Dirichlet Process Mixture model

A common approach to assign data to clusters is to construct a model in which data are generated from a mixture of probability distributions. In this way, the co-exposures of the n individuals to the P pesticides arise from a distribution composed of different sub-distributions, namely the mixture components. Therefore, the groups of individuals with similar patterns of pesticides co-exposure are identified as the ones sharing the same sub-distributions. Let the observed co-exposures x = (x 1 , . . . , x i , . . . , x n ) with x i a P dimensional vector x i = (x i1 , . . . , x ip , . . . , x iP ), be distributed with a probability density

f (x i ) = Θ k(x i |θ)G(dθ) (1)
where k(.|θ) is the known density of the mixture components called the kernel density, with parameter θ ∈ Θ and G the unknown mixing distribution. Under a nonparametric perspective, the unknown distribution G is one of an infinite-dimensional function space. Equation (1) can be broken down by introducing the latent variables θ i that will be shared among individuals with similar patterns of pesticides co-exposure, as follows

x i |θ i ∼ k(dx|θ i ) (2) θ i |G ∼ G(dθ) G ∼ P (G).
In a Bayesian approach, the challenge is to place an appropriate prior P (G) on the distribution G. A collection of distribution functions called random probability measures (RPMs) can be assigned to the density G [START_REF] Walker | Bayesian nonparametric inference for random distributions and related functions[END_REF][START_REF] Muller | Nonparametric Bayesian data analysis[END_REF]. [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF]

G(A 1 ), . . . , G(A k )) is drawn from a Dirichlet distribution (G(A 1 ), . . . , G(A k )) ∼ Dir(γH(A 1 ), . . . , γH(A k )).

Partition Models

A clustering of n objects can be represented by a partition, denoted as p [START_REF] Quintana | Bayesian Clustering and Product Partition Models[END_REF]. The partition p separates the n vectors x i into n(p) groups of individuals. The partition of size n(p) ∈ {1, . . . , n}, can be represented as p = {C 1 , . . . , C n(p) } where C j denotes the j th cluster for j = 1, . . . , n(p). Equation ( 1) can be expressed conditionally on the partition p, as a classification likelihood [START_REF] Lau | Model based clustering and weighted Chinese restaurant processes[END_REF])

f (x|p) = n(p) j=1 k(x i , i ∈ C j ),
where

k(x i , i ∈ C j ) is given by k(x i , i ∈ C j ) = Θ i∈C j k(x i |θ)G(dθ).
In the context of this classification likelihood, the partition p is the parameter for which a prior/posterior analysis is required. The prior distribution of p is induced by the distribution P (G)

and is proportional to

n(p) j=1 g(C j )
, where g (known as the cohesion) is a function of the cluster only, e.g. its size. By Bayes theorem, the posterior distribution of p is the following product over partition components

π(p|x) ∝ n(p) j=1 g(C j )k(x i , i ∈ C j ). ( 3
)
Choosing a Dirichlet process DP (γ, H) as the prior distribution P (G) is equivalent to considering g(C j ) = γΓ(e j ), where e j is the size of cluster C j , see [START_REF] Lau | Model based clustering and weighted Chinese restaurant processes[END_REF] for other settings.

An estimate of the optimal partition is the one that maximizes the posterior distribution (3), which is approximated in this paper with a Gibbs sampler as described in Section 3.3. The number of clusters in the optimal partition represents the number of sub-populations of individuals with similar pesticide co-exposure patterns.

Stick-breaking representation of the Dirichlet process

If G ∼ DP (γ, H) then it can be represented as an infinite mixture of point masses

G = ∞ k=1 β k δ φ k
where φ k are random variables sampled from H, δ φ k refers to a point mass concentrated at atom φ k and β k are the "stick-breaking weights" depending on γ. Drawing θ i from G in equation 2 means that θ i is equal to one of φ k with the associated probability β k .

One way to deal with such infinite mixtures is to use truncation (see [START_REF] Walker | Sampling the Dirichlet mixture model with slices[END_REF] for a slice sampler based alternative). [START_REF] Ishwaran | Gibbs Sampling Methods for Stick-Breaking Priors[END_REF] have shown that when truncating the sum to a reasonable N (N < ∞), the quality of approximation of G is good. The random weights β k can therefore be built from auxiliary weights β * k ∼ Beta(1, γ) through the stick-breaking procedure given by

β 1 = β * 1 , β k = β * k k-1 l=1 (1 -β * l ) for k = 2, . . . , N -1, and β N = 1 - N -1 k=1 β k .
In our study, clusters of individuals with similar patterns of co-exposure are identified as the ones sharing the same atoms φ k .

Specific models

Multivariate Normal mixture model A P dimensional multivariate Normal distribution N P (µ, τ -1 ) with mean vector µ ∈ R P and random covariance matrix τ -1 ∈ R P ×P is assigned to the kernel density k.

The distribution H is chosen to be a Wishart-Normal (µ, τ ) ∼ W N (α, Ψ, m, t) distribution due to its conjugate properties with the multivariate Normal distribution. Writing (µ, τ ) ∼ W N (α, Ψ, m, t) means that a Wishart distribution is used for the symmetric and positive definite precision matrix τ as τ ∼ W (α, Ψ), where α is a scalar degree of freedom and Ψ a P × P scale matrix; and conditionally to τ , the random vector µ is assigned a P -dimensional Normal distribution µ|τ ∼ N P m, (tτ ) -1 . Hence, we have

H(dµ, dτ ) = 2 -αP/2 |Ψ| α/2 (Γ P (α/2)) -1 × |τ | (α-P -1)/2 exp - 1 2 T r (Ψτ ) × (2π) -P/2 |tτ | 1/2 exp - t 2 (µ -m) τ (µ -m) dµdτ
where Γ P is the multivariate Gamma function Γ P (α/2) = π P (P -1)/4 d r=1 Γ( α+1-r

2

) and T r(A) is the trace of the matrix A.

The marginal density k(x i , i ∈ C j ) is obtained as

k(x i , i ∈ C j ) = i∈C j k(x i |µ, τ )H(dµ, dτ ).
With respect to our distribution choice, the marginal density is written as

k(x i , i ∈ C j ) = Γ P (α * j /2) Γ P (α/2) t P/2 π P e j /2 t * j P/2 |Ψ| α/2 Ψ * j α * j /2 ,
where ( * j ) indicate the updated values of the parameters of the Wishart-Normal denoted W N (α * j , Ψ * j , m * j , t * j ), and equal to

α * j = α + e j , m * j = tm + e j x j t * j , t * = t + e j , Ψ * j = Ψ + S j + e j t t * j (m -x j )(m -x j )
, where e j is the number of observations classified in the cluster C j , x j = 1 e j i∈C j

x i is the mean of the cluster C j and S j = i∈C j (x i -x j ) (x i -x j ) is the corresponding covariance matrix.

The optimum number of mixture components n(p) and the posterior distribution of the component mixture parameters are obtained by maximizing the following criterion (similar to Eq. 14 in [START_REF] Lau | Model based clustering and weighted Chinese restaurant processes[END_REF]) corresponding up to a constant to the posterior empirical log likelihood (i.e.

the log of equation ( 3))

Q(p) = n(p) × ln(γ) + n(p) j=1 ln Γ(e j ) + n(p) j=1 ln k(x i |i ∈ C j ). ( 4
θ im |G i ∼ G i G i ∼ DP (α i , G 0 ) G 0 ∼ DP (γ, H).
With such a model, the co-exposure to the P pesticides of each individual i is composed of several sub-distributions identified as the ones sharing the same θ im . Therefore, the form of the individual exposures to each pesticide given by the M values is considered into the clustering process. The algorithm of the base model described in Equation ( 2) is presented below in three steps1 .

Algorithm

Considering G ∼ DP (γ, H), we use the stick-breaking representation of the Dirichlet process G = N k=1 β k δ φ k (.), where φ k denotes the hidden parameters of the multivariate Gaussian distribution (µ k , τ k ), β = (β 1 , . . . , β N ) are the stick-breaking weights, and N is the maximum number of atoms of G. A vector K = (K i ) is introduced to store the assignment of each data point x i = (x ip , p = 1, . . . , P ) for i = 1, . . . , n to an atom (φ k ) k=1,...,N so that K i is an integer from 1 to N . Only N * of the N available atoms are distinct values whose set is denoted by K * .

1. Sampling (φ|K, β, X) : for those k in K * , sample φ k with respect to the "updated" base measure H * k (a Wishart-Normal with parameters α * k , Ψ * k , m * k , t * k obtained from the posterior distribution given by {x i , K i = k}) and for the remaining (N -N * ) atoms, get φ k from the base measure H (the prior Wishart-Normal(α, Ψ, m, t)).

2. Sampling (K|β, φ, X) : for k = 1, . . . , N, and i = 1, . . . , n

Pr(K i = k) ∝ β k × k(x i |φ k ). ( 6
)
3. Sampling (β|φ, K, X) : based on the β * k ∼ Beta(1 + e k , γ + N l=k+1 e l ) with e k corresponding to #{x i , K i = k} ≤ n for k = 1, . . . , N, then

β 1 = β * 1 , β k = β * k k-1 l=1 (1 -β * l ), for k = 2, . . . , N -1, β N = 1 - N -1 l=1 β l 3.3.1.

b Hierarchical model

The algorithm for the stick-breaking representation of the hierarchical Dirichlet process presented in equation ( 5) requires the sampling of additional intermediate weights π = (π ik ) and the definition of a matrix K = (K im ) describing the assignment of each data point x im = (x pim , p = 1, . . . , P ) to one of the N atoms, for i = 1, . . . , n and m = 1, . . . , M . Steps 2 and 3 are replaced with steps 2' and 3' below.

2'. Sampling (K|π, β, φ, X) for k = 1, . . . , N , i = 1, . . . , n and m = 1, . . . , M

Pr(K im = k) ∝ π ik × k(x im |φ k )
3'. Sampling (π|β, φ, K, X), independently on the fixed i's and based on

π * ik ∼ Beta α 0 β i + e ik , α i 1 - k l=1 β l + N l=k+1 e il ,
with e ik corresponding to #{x im , K im = k} ≤ M for i = 1, . . . , n and k = 1, . . . , N , then

π i1 = π * i1 , π ik = π * ik k-1 l=1 (1 -π * il ), for k = 2, . . . , N -1, π iN = 1 - N -1 l=1 π il . Note that π * iN ∼ Beta(α i β N + e iN , 0
) is a properly defined Beta distribution. Finally sampling (β|π, φ, K, X) is performed exactly as described in the original step 3.

Learning about the parameter γ

Considering γ as a random parameter with prior distribution Γ(a γ , b γ ) leads to an additional last step as proposed in [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF].

4. Sampling (γ|φ, K, π, β, X) based on an auxiliary variable γ * ∼ Beta(γ + 1, n), according to the following mixture distribution

γ ∼ w γ * × Γ(a γ + k, b γ -ln γ * ) + (1 -w γ * ) × Γ(a γ + k -1, b γ -ln γ * ), (7)
with weights w γ * defined by

w γ * 1-w γ * = aγ +k-1 bγ -ln γ * = c γ * that is w γ * = c γ * 1+c γ * .

Starting values of hyperparameters

Starting values of the hyperparameters were taken to be equal to α = P, Ψ = 0 P ×P , m = 0 P , t = 1, as proposed in [START_REF] Cabrera | Random Block Sampling for high dimensional clustering (from the Bayesian point of view)[END_REF] in order to use vague prior distributions. Different settings were tested on the simulated datasets for the starting value related to the parameter γ: no prior distribution but γ is set to 1, a prior Gamma distribution Γ(a γ , b γ ) with hyperparameters equal to:

(a γ , b γ ) = (2, 4) for informative prior and equal to (a γ , b γ ) = (1, 1) and (a γ , b γ ) = (0.01, 0.01) for more vague prior distributions, yielding the posterior given in ( 7). In the case of the hierarchical model, the weights α i were set to 1 for each i = 1, . . . , n.

Algorithm convergence checking

The convergence of the Gibbs sampler to the optimal partition was visually checked by plotting the Q-criterion and the number of clusters n(p) against the number of realized iterations. The number of iterations necessary to reach the optimal partition depends on the size of the dataset. Therefore, the algorithm is stopped after checking that the Q-criterion stabilizes over a minimum of 20, 000 iterations. The optimal partition is determined over all the iterations.

Random-block Gibbs Stick-breaking

Cabrera et al. ( 2009) introduced a novel procedure called the random-block Gibbs weighted Chinese restaurant process algorithm to reduce the heavy computational burden in estimating the optimal partition induced by the Gibbs sampler and the high dimensionality of the data. We propose to apply this method to the SB algorithm. The principle is to randomly reduce the dimension of the data by selecting a number d (d < P ) of pesticides from the original number P at each Gibbs cycle. Therefore, given the sequence of random integers v d = {l 1 , . . . , l d }, a subset of observations

x i = (x il 1 , . . . , x il d
) is used instead of the x i = (x i1 , . . . , x iP ) for the i = 1, . . . , n individuals. This procedure will be referred to as RB-SB below.

Application

Simulated datasets

To investigate the quality of the clustering estimates under various settings, a simulation study was conducted for the stick-breaking algorithms applied to two datasets based on the settings proposed by [START_REF] Cabrera | Random Block Sampling for high dimensional clustering (from the Bayesian point of view)[END_REF]. The datasets were built from a three component mixture of P = 5 dimensional multivariate Normal distributions denoted N P (µ k , Σ k ) k=1,...,3 . The parameters (µ k , Σ k ) are detailed in Table 1.

Dataset for the base model

A sample of 1, 000 values (x i ) was built from f (x) = 3 k=1 β k N P (µ k , Σ k ) with (β 1 , β 2 , β 3 ) = (0.3, 0.3, 0.4).
After 30, 000 iterations, and with the parameter γ set to 1, the optimal partition was obtained at the 34 th iteration with the Q-criterion of equation ( 4) being equal to -3, 687, see Fig. 1(a). The optimal partition was composed of 3 clusters corresponding to the 3 components of the mixture dataset, see Fig. 1(b). Performing the RB-SB algorithm with the dimension reduced to d = 2, the maximum Q-criterion is reached at the 650 th iteration, also with an optimal partition matching the generating one (see Figure 1(b)). When using different prior distributions for the parameter

k µ k Σ k 2
1.0 0.5 0.2 0.1 0.1 2 0.5 2.0 0.5 0.2 0.1 1 4 0.2 0.5 1.0 0.5 0.2 5 0.1 0.2 0.5 1.0 0.5 6 0.1 0.1 0.2 0.5 3.0 -2 1.0 0.5 0.2 0.1 0.1 -2 0.5 2.0 0.5 0.2 0.1 2 -4 0.2 0.5 3.0 0.5 0.2 -5 0.1 0.2 0.5 2.0 0.5 -6 0.1 0.1 0.2 0.5 1.0 -5 3.0 0.5 0.2 0.1 0.1 5 0.5 1.0 0.5 0.2 0.1 3 -7 0.2 0.5 2.0 0.5 0.2 7 0.1 0.2 0.5 2.0 0.5 -9 0.1 0.1 0.2 0.5 3.0 γ, the same maximum value of Q = -3, 687 is reached, as with γ being set to 1. Figure 2 shows the 3 different prior distributions attributed to the parameter γ and their corresponding posterior distributions. 50% of posterior values of γ are below 0.27 < 1, showing that it is possible, though perhaps not crucial here, to learn about γ.

Dataset for the hierarchical model

To reproduce the co-exposure data structure including uncertainty of exposure, 240 individuals were generated from the f (x) = 3 k=1 β k N P (µ k , Σ k ) with (β 1 , β 2 , β 3 ) = (0.33, 0.17, 0.5). For each individual i, M = 100 values were sampled resulting in a total sample of 24, 000 observations. The convergence of the SB and the RB-SB algorithms to the optimal partition was very slow. After 200, 000 iterations, the number of clusters which maximized the Q-criterion was 11. The size of the clusters ranged from 269 to 16, 089 observations. The biggest clusters had means and covariances similar to those used to generate the dataset. This suggested that if the number of iterations was increased, three main clusters would become apparent and the optimal partition would eventually converge to that generating the data. Note: To read parallel coordinates consider that the x-axis is the dimension j, the y-axis is the value of x j , a line connects x 1 to x 2 , x 2 to x 3 and so on. The color indicates the 3 clusters resulting from the SB or RB-SB algorithms, which both match the generating ones described in Table 1 For the children sample, the optimal partition was obtained after 98, 362 iterations and is composed of 16 clusters. As with the adults, 14 clusters accounting for only 45 individuals were discarded to focus on the 2 main ones. The first cluster is composed of 743 children who are highly exposed to a large number of pesticides, and the second one consists of 651 children. As with the adults sample, the correlation matrix of the most heavily exposed sub-population was analyzed to determine cocktails of pesticides. There are 39 pesticides with correlations over 0.95, divided axis 1, which represents 68% of the total variance is mainly determined by the 34 pesticides selected with our base model. Indeed, among the 24 pesticides which mainly contribute to axis 1, 23 were also selected with the base model. The second axis only represents 6% of the total variance. The coordinates of the 79 pesticides are positive on axis 1, while the individuals are found on both sides of the axis. The results obtained with the PCA lead to the same conclusion that there are two groups of individuals, one highly exposed and one only marginally exposed to all the pesticides.

Hierarchical model for distributions of individual co-exposures

The hierarchical model was developped to deal with distributions of individual co-exposure to the P pesticides in order to account for uncertainty of residue levels. Because of the algorithm's long computational time, the hierarchical model was only applied to the 34 pesticides selected with the base model for the adult population. From 30, 000 iterations, the optimal partition was found after 25, 853 iterations and was composed of 6 clusters. From these clusters, three were composed of large sets of observations according to the clustering obtained with the base model (see previous subsection). Moreover, two of these clusters were composed of highly exposed individuals to a large proportion of the 34 pesticides. The distribution of co-exposure of each individual was mostly found to consist of 2 or 3 component mixtures. The analysis of the correlation matrix of the two main clusters showed that the correlations between pesticides were very low, ranging between 0.2 and 0.45, compared to those obtained with the base model. These low correlations could be due to the high uncertainty associated with the exposure to each pesticide. Indeed, for a random set of contamination selected for the P pesticides, an individual can be exposed to a low level for one pesticide and a high level for another, leading to low correlations. The inclusion of the exposure uncertainty is more realistic in terms of exposure assessment but implies some difficulties in the definition of the cocktails of pesticides as we can place ourselves at the different hierarchy levels to define the cocktails, which results in too many different cocktails.

Conclusion

This paper presents a Bayesian nonparametric model based on Dirichlet process mixtures, applied to cluster the co-exposures of the French population to various pesticides in order to define cocktails of pesticides which are relevant to study in terms of human health effects. The hierarchical model applied in this work is original in the field of food risk analysis although it has been applied in other fields, mainly to deal with functional data [START_REF] Rodríguez | Bayesian nonparametric functional data analysis through density estimation[END_REF][START_REF] Teh | Hierarchical Dirichlet Processes[END_REF]. Compared to the PCA or other clustering methods, this type of Bayesian nonparametric model has several advantages, for example it is not based on the linearity assumption. Also, it can comprise an infinite number of mixture components which is particularly suited for the case of high data dimensionality.

Moreover, the number of clusters is automatically determined through the estimation process, no parametric assumption with regard to the form of the co-exposure distribution is required and the structure of the data set may be introduced through a specific hierarchy to account for exposure uncertainty. As a result, the model applied to the 95 th percentile of the French population coexposure to pesticides clustered the adult population into two sub-populations: individuals highly exposed to a large proportion of pesticides and individuals slightly exposed to a large proportion of pesticides. Thus, cocktails defined from individuals highly exposed are numerous and composed of various pesticides, for example for the adults sample 34 pesticides shared into 20 cocktails have been selected. Nevertheless, this approach relies on the major assumption that the non-detect values are uniformly distributed between 0 and the LOR. This assumption may not be realistic for pesticides residues, for which a real 0 can occur when the pesticide is not used on a crop, (EFSA, 2010). Then, refinement on the way to deal with non-detects is necessary in order to define more realistic cocktails.

)

  Hierarchical model to account for the uncertainty of the exposure An additional Dirichlet process is used to account for the uncertainty of the exposure when considering the set of data x im = {x pim , p = 1, . . . , P } for each individual i = 1, . . . , n and the contamination value m = 1, . . . , M x im |θ im ∼ k(.|θ im )(5)

Figure 1 :

 1 Figure 1: The Stick-Breaking (SB) and the Random-Block Stick-Breaking (RB-SB) algorithms for the base model applied to a simulated dataset (N = 30 atoms and 30, 000 iterations).Note: To read parallel coordinates consider that the x-axis is the dimension j, the y-axis is the value of x j , a line connects x 1 to x 2 , x 2 to x 3 and so on. The color indicates the 3 clusters resulting from the SB or RB-SB algorithms, which both match the generating ones described in Table1

Figure 2 :

 2 Figure 2: Densities of the Gamma(a γ , b γ ) prior (solid line) and the posterior (dashed line) distributions of the parameter γ. Red line: (a γ , b γ ) = (2, 4), blue line: (a γ , b γ ) = (1, 1), green line: (a γ , b γ ) = (0.01, 0.01), prior densities are plotted in logarithmic scale.

  stated the properties of this class of measures and introduced the Dirichlet Process (DP) as one of the RPMs. The DP is defined by two parameters, a scaling parameter γ and a base probability measure

H. The probability distribution G is drawn from a DP, denoted G ∼ DP (γ, H), if and only if for any partition (A 1 , . . . , A k ) of Ω, the vector of random probabilities (

Table 1 :

 1 Parameters of the multivariate Normal distributions N P (µ k , |Sigma k ) k=1,...,3 used to generate the simulated datasets

The code is available on the webpage http://www.paris.inra.fr/metarisk/members/tressou_jessica/ publications

Acknowledgment

This work is a part of a project funded by the French National Research Agency (ANR) and the French Agency for Environmental and Occupational Health Safety (AFSSET). The authors thank Fanny Héraud for her relevant comments.

into 13 different cocktails. Of these 39 pesticides, 28 are similar to those obtained with the adult population.

The analyses of the pesticide exposures of the main clusters show that the individuals are either highly exposed to all pesticides or have a low exposure to all pesticides. A preliminary conclusion is that this clustering is mostly related to the consumption behavior of individuals rather than the contamination levels.

-

Comparison with PCA

To compare the results obtained with the base model to those from a classical method, a principal component analysis (PCA) for the adult population was realized and described in this section. The