# Sensitivity analysis for high quantiles of ochratoxin A exposure distribution

Isabelle Albert, Jean-Pierre Gauchi

► To cite this version:

Isabelle Albert, Jean-Pierre Gauchi. Sensitivity analysis for high quantiles of ochratoxin A exposure distribution. International Journal of Food Microbiology, 2002, 75 (1-2), pp.143-155. 10.1016/S0168-1605(01)00747-4. hal-01263591

HAL Id: hal-01263591

https://hal.science/hal-01263591v1

Submitted on 31 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Sensitivity analysis for high quantiles of ochratoxin A exposure distribution

I. Albert*, J.-P. Gauchi

Food Risk Analysis Group, Biometrics Unit, National Institute for Agricultural Research (INRA),
Domaine de Vilvert, 78352 Jouy-en-Josas cedex, France

Received 28 March 2001; received in revised form 29 October 2001; accepted 17 November 2001

## Abstract

Using available data from a consumption survey and contamination data on ochratoxin A (OA) in food, a sensitivity analysis (SA) for high quantiles (95th and 99th quantiles) of OA exposure distribution was carried out, obtained by a Monte Carlo simulation in French children. Exposure assessment for food contaminants is important to control the risk of foodborne diseases. Risk assessors are interested in high quantiles of contaminant exposure distributions. As these exposure distributions are generally very asymmetrical, it is difficult to obtain relevant and stable high quantiles in such a context. Determining OA exposure distribution is complex because it is based on the sum of elementary exposure distributions (eight foodstuffs are analysed here), and each one of these is the product of a consumption distribution and a contamination distribution. The SA enables us to quantify the influences of the parameter variability of the consumption and contamination probability density functions (pdf) which have been fitted to the data, our simulation model inputs, on the 95th and 99th quantiles of the output exposure distribution. After some preliminary trials, we have postulated a quadratic polynomial regression model for the quantiles of OA exposure distribution in view of undertaking this SA. This regression model comprises 32 main factors, their 496 two-factor interactions and their 32 quadratic terms. The 32 factors are the parameters of the fitted pdf: 16 parameters of Gamma distributions relative to the eight consumed foods and 16 parameters of Gamma distributions relative to the eight food OA contaminations. For an optimal parameter estimation of such a large model, we used an experimental design approach depending on a resolution-V fractional factorial design of 6561 experiments. The factor ranges are established by a preliminary study of bootstrap sampling. From the bootstrap samples, the factor ranges are obtained taking into account the correlation between the two parameters of the fitted Gamma pdf. A full exposure distribution is simulated for each of the 6561 experiments. The consumption dependencies are taken into account by the Iman and Conover method. On the basis of this analysis, validated and useful models for each desired quantile are obtained showing a major influence of the parameters of “Cereals” (consumption and contamination) and slightly less so for parameter of “Pork” consumption in the sensitivity of the quantiles. © 2002 Elsevier Science B.V. All rights reserved.

**Keywords:** Sensitivity analysis; Exposure assessment; Monte Carlo; High quantiles; Regression; Design of experiments; Food safety

<footnote>* Corresponding author. Tel.: +33-1-3465-2843; fax: +33-1-3465-2217.</footnote>

<footnote>E-mail addresses: Isabelle.Albert@jouy.inra.fr (I. Albert), Jean-Pierre.Gauchi@jouy.inra.fr (J.-P. Gauchi).</footnote>
# 1. Introduction

Ochratoxin A (OA) is a mycotoxin produced by fungi of the genera *Aspergillus* and *Penicillium*. This mycotoxin may be a contaminant of grain stored in poor conditions, especially in Europe, and through the food chain also contaminates other foodstuffs, especially pork and poultry meat. Several toxicological studies (e.g. Pfohl-Leszkowick, 1999; JECFA, 1995; CCFAC, 1999) have shown this mycotoxin is an acute kidney toxin in humans. For these reasons, this mycotoxin remains under the close surveillance of the «Conseil Supérieur d'Hygiène Publique» in France and thus evaluating human exposure to this mycotoxin is an important public health question.

A consumption survey (ASPCC, 1999)—was conducted from June 1993 to June 1994 on 1161 French individuals (children, women and men) to estimate risk assessment exposure to OA in food. During 7 days, the participants were questioned on their consumption of several types of food. Our sensitivity analysis (SA) was restricted to children (individuals younger than 15) and to eight types of food, given in Table 1, keeping in mind that certain cover several foods, consumed by children. Table 2 gives descriptive statistics of these data for child consumers and non-consumers, and for men and women descriptive statistics are given in Verger (1999) and Gauchi (in press).

From consumption and contamination data, several exposure statistics have been computed by a non-parametric method—the NP-NP method—and a parametric method—the MP-P method (Gauchi, in press). In this paper, we use the MP-P method to

<table><caption>Table 2
Consumption in grams over an average day for all children (French children) and only those who consume a particular type of food</caption><thead><tr><th>Food</th><th>n</th><th>Mean</th><th>Standard deviation</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td colspan="6"><em>All the children</em></td></tr><tr><td>Cereals</td><td>232</td><td>165.36</td><td>89.87</td><td>26.14</td><td>619</td></tr><tr><td>Raisins</td><td>232</td><td>0.25</td><td>1.19</td><td>0</td><td>11.71</td></tr><tr><td>Other dried fruit</td><td>232</td><td>0.41</td><td>2</td><td>0</td><td>17.14</td></tr><tr><td>Pork</td><td>232</td><td>39.12</td><td>25.98</td><td>0</td><td>140</td></tr><tr><td>Poultry</td><td>232</td><td>16.38</td><td>16.78</td><td>0</td><td>95.71</td></tr><tr><td>Fruit juices</td><td>232</td><td>81.67</td><td>149.72</td><td>0</td><td>1497.14</td></tr><tr><td>Wines</td><td>232</td><td>0.53</td><td>2.47</td><td>0</td><td>23.86</td></tr><tr><td>Coffee</td><td>232</td><td>0.34</td><td>1.47</td><td>0</td><td>12.57</td></tr><tr><td colspan="6"><em>The children who consume:</em></td></tr><tr><td>Cereals</td><td>232</td><td>195.36</td><td>89.87</td><td>26.14</td><td>619</td></tr><tr><td>Raisins</td><td>17</td><td>3.39</td><td>3</td><td>0.43</td><td>11.71</td></tr><tr><td>Other dried fruit</td><td>15</td><td>6.29</td><td>5.12</td><td>1.29</td><td>17.14</td></tr><tr><td>Pork</td><td>219</td><td>41.44</td><td>24.87</td><td>0.57</td><td>140</td></tr><tr><td>Poultry</td><td>167</td><td>22.76</td><td>15.68</td><td>1.43</td><td>95.71</td></tr><tr><td>Fruit juices</td><td>157</td><td>120.69</td><td>168.68</td><td>1.43</td><td>1497.14</td></tr><tr><td>Wines</td><td>25</td><td>4.94</td><td>6</td><td>0.43</td><td>23.86</td></tr><tr><td>Coffee</td><td>24</td><td>3.33</td><td>3.36</td><td>0.14</td><td>12.57</td></tr></tbody></table>

<table><caption>Table 1
Definition of the eight types of food analysed</caption><tr><td>Cereals:</td><td>rice, wheat, corn starch, bread, pasta, flour, crispbread, breakfast cereals, crackers, biscuits, pastries, croissants and the like, popcorn, cornflakes, cooked sweet corn</td></tr><tr><td>Raisins:</td><td>dried raisins</td></tr><tr><td>Other dried fruits:</td><td>apricots, bananas, dates, figs, prunes</td></tr><tr><td>Pork:</td><td>kidneys, cutlets, fillet, roast, ribs, loin, black pudding, ham, pâté, salami, rillettes, andouillettes, bacon, sausages, foie gras</td></tr><tr><td>Poultry:</td><td>chicken liver, chicken, duck, goose</td></tr><tr><td>Fruit juices:</td><td>fruit juices</td></tr><tr><td>Wines:</td><td>several different wines</td></tr><tr><td>Coffee:</td><td>dry coffee (correction by a dilution coefficient of 17)</td></tr></table>

evaluate the high quantiles of the OA exposure distribution in the SA. Indeed, the parametric method (the mixed parametric/parametric method, named the MP-P method in the text) seems to be very suitable for estimating relevant and stable high quantiles. This MP-P method, based on a Monte Carlo simulation, depends on the random samplings from fitted pdf (mixed parametric pdf for consumption and parametric pdf for contamination). The 95th and 99th quantiles of the exposure distribution are the model outputs (the responses) of our SA. The inputs are the parameters of the pdf fitted to the data. The SA enables us to quantify the influences of these parameters (and their interactions) on the responses via a quadratic polynomial model postulated for modelling the variation of these responses. For an optimal estimation of the parameters of the polynomial model, we used an experimental design approach depending on a resolution-V fractional factorial design of 6561 experiments. The factors of this design are the fitted pdf parameters. To simulate a meaningful range for each factor, we determine it by a parametric bootstrap approach (Efron and Tibshirani, 1993) carried out on the consumption and contamination samples. In
order to consider the correlation between the pdf parameters, three experimental domains of variation were designed inside a 95% concentration ellipse.

This type of SA should be distinguished from marginal analysis: what is the effect of infinitely small changes (perturbations)? The effects in the SA are measured by the parameter vector $\beta$ of the postulated quadratic polynomial regression model.

In Section 2, we present how we obtain the OA exposure distribution and the methodology of our SA. In Section 3, we give the results, and we conclude in Section 4 by a discussion.

## 2. Methods

The contamination data (µg kg⁻¹ OA in food) come from elements independent of the consumption survey: they originate from SCOOP task 3.2.7. (see detailed reference). Table 3 gives descriptive statistics of OA contamination in the eight types of food. For poultry contamination data, the only information available was 31 values between 0 and 0.03, and 62 values between 0.03 and 0.18: the mean of 0.075 (Table 3) is calculated by the weighted arithmetic mean of the two means calculated from the two uniform distributions on [0; 0.03] and [0.03; 0.018]. On the other hand, the standard deviation of the unknown data is estimated by the Bienaymé–Tchebysheff inequality by supposing that the range [0; 0.18] encompasses 100% of the data. This inequality indicates that 75% of the data are contained in the interval [μ±2σ], which, with an estimation of 0.075 for μ,

<table><caption>Table 3
Descriptive statistics of the contamination data (µg kg⁻¹) for each food</caption><thead><tr><th>Food</th><th>n</th><th>Mean</th><th>Standard deviation</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>Cereals</td><td>183</td><td>0.65</td><td>0.98</td><td>0.20</td><td>7.20</td></tr><tr><td>Raisins</td><td>13</td><td>0.71</td><td>1.16</td><td>0.10</td><td>4.30</td></tr><tr><td>Other dried fruit</td><td>33</td><td>0.28</td><td>0.29</td><td>0.10</td><td>1.60</td></tr><tr><td>Pork</td><td>1011</td><td>0.30</td><td>0.29</td><td>0.10</td><td>6.10</td></tr><tr><td>Poultry</td><td>93</td><td>0.075</td><td>0.0342</td><td>0.00</td><td>0.18</td></tr><tr><td>Fruit juices</td><td>19</td><td>0.25</td><td>0.78</td><td>0.02</td><td>3.45</td></tr><tr><td>Wines</td><td>104</td><td>0.16</td><td>0.26</td><td>0.00</td><td>1.64</td></tr><tr><td>Coffee</td><td>155</td><td>1.12</td><td>1.70</td><td>0.02</td><td>3.10</td></tr></tbody></table>

leads here to an estimation of 0.0342 for $\sigma$. In the following, poultry contamination data are simulated by a Gamma distribution with mean = 0.075 and standard deviation = 0.0342 corresponding to a shape parameter = 4.81 and a scale parameter = 64.12. Histograms of consumption and contamination data are given in Figs. 1, 2 and 3.

Several methods and tools are needed to carry out this SA. First, we present the MP–P method used for estimating an exposure distribution and the quantiles that we are interested in. We refer the reader to previous work (Gauchi, in press) for more detailed information about this method. Second, for the SA, we expose the four-step methodology used to attain our objectives: (i) postulation of a polynomial regression model; (ii) construction of the appropriate experimental design; (iii) determination of the factor ranges; (iv) validation of the model.

### 2.1. Construction of the exposure distribution

#### 2.1.1. Definition of exposure

Determining OA exposure distribution is rather complicated because it is based on the sum of eight elementary exposure distributions (one per type of food analysed), and each one of these is the product of a consumption distribution and a contamination distribution. Let us give a rigorous definition of the exposure in our context. Consider the consumption $C_j$ of one food $j$, divided by the individual weight, i.e. a normalized consumption, and the contamination $T_j$ in OA of this food. Then, for $p$ foods consumed per day, we define the normalized global exposure $E$ as the sum of the $p$ normalized products consumption×contamination: $E = \sum_{j=1}^{p} C_j T_j$, expressed in nanograms per kg of body weight per day, i.e. in ng×bw⁻¹×day⁻¹ units.

#### 2.1.2. The MP–P method

The MP–P method, described in detail previously (Gauchi, in press), is essentially parametric in the sense that a mixed pdf is adjusted on each food consumption and a parametric pdf is adjusted on each food contamination. For consumption, apart from cereals, the typical histogram that can be plotted for each foodstuff is shown in Fig. 4, where two zones appear corresponding to the two sub-classes: non-consumers (sub-class of proportion $h$, represented by
<figure><img src="image_2.png" /><figcaption>Fig. 1. Consumption histograms (relative frequencies) for children studied (French children).</figcaption></figure>

the tall bar on the left) and consumers (sub-class of proportion (1-h), represented by the asymmetrical part). In this situation, we propose for each consumed foodstuff j to fit a mixed distribution defined by a continuous uniform distribution for the non-consum-

ers and a Gamma distribution for the consumers. An example of such a mixed pdf has been given by Vose (1996). The Gamma distribution has been validated by means of the Anderson–Darling statistic (D’Agostino and Stephens, 1986) among several asymmetrical
<figure><img src="image_3.png" /><figcaption>Fig. 2. Normalized consumption histograms (relative frequencies) and fitted Gamma pdf for children consumers (French children).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Contamination histograms (relative frequencies) and the fitted Gamma pdf.</figcaption></figure>

distributions. In the same way, this statistic has been
used to decide Gamma distributions for contamina-
tion data.

The Gamma pdf for a continuous random variable
X defined in [0, +∞>]
[is: f(x, r, λ) = [1/Γ(r)]λ⁻¹
(λ⁻¹(x - θ))r⁻¹e⁻¹e⁻⁻¹], where r, λ, θ are the shape,
<figure><img src="image_5.png" /><figcaption>Fig. 4. Typical consumption histogram.</figcaption></figure>

scale and threshold parameters, respectively, and $\Gamma(r)$ is the usual Euler's integral. The mean and variance of the Gamma distribution are connected to $r$ and $\lambda$ in the following way:

$$E(X) = r\lambda \quad \text{and} \quad V(X) = r\lambda^2. \tag{1}$$

The estimated parameters of the corresponding pdf are given in Table 4. The fitted Gamma distributions are presented in the histograms Figs. 2 and 3.

An exposure value is calculated as follows:

$$\hat{E}_i^{(S)} = \sum_{j=1}^{p} \tilde{c}_{i,j} \tilde{t}_{i,j}, \tag{2}$$

where:

- $\tilde{c}_{i,j}$ is a random normalized consumption for the foodstuff $j$, drawn from the corresponding

<table><caption>Table 4
Parameters of the Gamma pdf fitted to the normalized consumption data (only to those who consume) and parameters of the Gamma pdf fitted to the contamination data</caption><thead><tr><th rowspan="2">Food</th><th colspan="2">Parameters of the fitted Gamma pdf</th></tr><tr><th>Consumption</th><th>Contamination<br>(θ=0)</th></tr></thead><tbody><tr><td>Cereals</td><td>r̂ = 2.594; λ̂ = 0.480; θ̂ = 1.610</td><td>r̂ = 0.439; λ̂ = 0.674</td></tr><tr><td>Raisins</td><td>r̂ = 0.701; λ̂ = 5.484; θ̂ = 0.010</td><td>r̂ = 0.375; λ̂ = 0.530</td></tr><tr><td>Other<br>dried fruit</td><td>r̂ = 0.582; λ̂ = 2.286; θ̂ = 0.043</td><td>r̂ = 0.998; λ̂ = 3.503</td></tr><tr><td>Pork</td><td>r̂ = 1.714; λ̂ = 1.127; θ̂ = 0.041</td><td>r̂ = 1.085; λ̂ = 3.591</td></tr><tr><td>Poultry</td><td>r̂ = 1.813; λ̂ = 2.499; θ̂ = 0.080</td><td>r̂ = 4.810; λ̂ = 64.12</td></tr><tr><td>Fruit juices</td><td>r̂ = 0.939; λ̂ = 0.209; θ̂ = 0.083</td><td>r̂ = 0.108; λ̂ = 0.423</td></tr><tr><td>Wines</td><td>r̂ = 0.764; λ̂ = 5.594; θ̂ = 0.092</td><td>r̂ = 0.413; λ̂ = 2.404</td></tr><tr><td>Coffee</td><td>r̂ = 0.784; λ̂ = 8.961; θ̂ = 0.033</td><td>r̂ = 0.431; λ̂ = 0.386</td></tr></tbody></table>

cumulative density function (cdf) of the above-mentioned mixed distribution, whose density parameters are given in Table 4;

- $\tilde{t}_{i,j}$ is a random contamination for the foodstuff $j$, drawn from the fitted Gamma cdf, whose density parameters are given in Table 4.

Then, with $N=10,000$ outputs of form (2), we obtain an exposure simulation set $S$. The eight vectors of the $N$ consumption random deviates have been arranged by means of the Iman and Conover (1982) method to be correlated with the aim of respecting the consumption dependencies. The $N$ outputs lead to calculate various usual statistics and notably the 95th and 99th quantiles. Fig. 5 shows the simulated output exposure histogram.

### 2.1.3. Calculation of quantiles

The $\alpha$th theoretical exposure quantile is here referred to as $Q_\alpha$ and defined by: $F_E(Q_\alpha) = \alpha$, where $F_E$ is the theoretical cdf of the exposure of the population studied. The $\alpha$th empirical exposure quantile $\hat{Q}_\alpha^{(S)}$ is defined by:

$$\#[\hat{E}_i^{(S)} \le \hat{Q}_\alpha^{(S)}]/N = \alpha, \tag{3}$$

where the notation $\#[x_i \le K]$ means the number of $x_i$ less than or equal to $K$. For example, if $N=10,000$, the estimate of the 0.95th quantile is the 9500th largest value of all the $\hat{E}_i^{(S)}$. If $\alpha N$, is not an integer, the

<figure><img src="image_6.png" /><figcaption>Fig. 5. Exposure output histogram (relative frequencies) with the MP-P method.</figcaption></figure>
following procedure can be used. Assuming $\alpha \le 0.5$, let $k = [(N+1)\alpha]$, the largest integer $\le (N+1)\alpha$. Then, we define the empirical $\alpha$th and $(1-\alpha)$th quantiles by the $k$th largest and $(N+1-k)$th largest values of all the $\hat{E}_i^{(S)}$, respectively. These quantiles are known as the empirical quantiles.

## 2.2. Methodology of the SA

The SA used is based on two main joint tools: a quadratic polynomial model and a simulation design which is optimal relative to this model. Our aim is to determine 95th and 99th quantile models for clearly expressing the influence of the factors chosen and, secondly, to have at our disposal a simulation tool to easily evaluate the amplitude of a quantile variation when new inputs are given.

### 2.2.1. Polynomial model

After some preliminary trials, we assumed a full quadratic polynomial regression model for the responses (outputs) studied. This regression model is based on 32 main factors, their 496 two-factor interactions and their 32 quadratic terms. The 32 factors are the pdf parameters: 16 parameters of the Gamma distributions relative to the eight consumed foods (each of these distributions is described by three parameters; the influence of the threshold parameter is not studied because it merely demarcates the lower position of the distribution of the consumers to consider the non-consumers), and 16 parameters of the two parameter Gamma distributions relative to the eight food OA contaminations. The general form of the model is:

$$y = \beta_0 + \sum_{j=1}^{32} \beta_j z_j + \sum_{j=1}^{31} \sum_{k=j+1}^{32} \beta_{jk} z_j z_k + \sum_{j=1}^{32} \beta_{jj} z_j^2 + \varepsilon, \quad (4)$$

where the $\beta$ are unknown parameters, the factors $z$ are the consumption and contamination pdf parameters, and $\varepsilon$ is an error term for which no particular probability distribution is assumed and we only suppose its expectation is zero. Model (4) is a second degree (quadratic)

polynomial in $z$. In the following, model (4) is successfully applied to the 95th and 99th empirical quantiles, $\hat{Q}_{0.95}$ and $\hat{Q}_{0.99}$ respectively, defined by Eq. (3).

### 2.2.2. Experimental design

In order to estimate optimally the $\beta$ parameters of model (4), an appropriate simulation design is chosen. As the three-factor interactions and over are assumed negligible, a full factorial design of $3^{32}$ experiments was not necessary, especially as it would have taken an extremely long time. Instead of this unfeasible design, we built a resolution-V fractional factorial design (Box et al., 1978) of 6561 experiments by means of the Factex procedure of the SAS software (SAS/QC® Institute). A resolution-V design allows us to estimate independently the parameters of the main effects of the factors (the $z_j$), the parameters of the 496 two-factor interactions (the $z_j z_k$), and the parameters of the 32 quadratic effects (the $z_j^2$) assuming that three-factor and over interactions are negligible. The factor levels were coded as -1, 0 and +1.

### 2.2.3. Factor ranges

The $z$ factors, i.e. the parameters of the consumption and contamination pdf, are the inputs of the SA. The aim is to make the inputs vary in a reasonable range but large enough to observe how the response outputs vary. Because we had only one data set for consumption and contamination, no hypothesis of the factor variations could be made, so no supplementary surveys or contamination analyses were available. Thus, to simulate a meaningful range for each factor, we decided to determine it by a parametric bootstrap approach (Efron and Tibshirani, 1993). We considered the correlation between the shape and scale parameter of the pdf. For each consumption and contamination sample (16 in all), 10,000 bootstrap samples were obtained. A parametric bootstrap sample is obtained by sampling $n$ times $[n$ is the size of the data sample; see Table 2 (consumers) and Table 3 for the values of $n]$ in the distributions fitted to the data. Table 4 gives the parameters of those distributions. From the mean and the standard deviation of each bootstrap sample, we calculated the scale and shape parameters from relation (1). As the bootstrap means were biased due to the sparseness of the data (see histograms Figs. 1–3), all the bootstrap scale and shape parameters were re-centred on the fitted param-
eter values given in Table 4. As the scatter plot of the 10,000 parameter pairs (scale and shape) showed an elliptical appearance for each food consumption and food contamination, we determined a variation domain by calculating a 95% concentration ellipse, one in each zone. In each ellipse, three squared variation domains were designed for the shape and scale parameters: the first one in the lower part of the ellipse (zone I), the second one in the middle part (zone II) and the third one in the upper part (zone III). Thus, in this approach, three polynomial models were determined for the 95th and 99th quantiles. Fig. 6 gives an example of the three zones for the consumption of cereals. For example, Tables 5, 6 and 7 indicate the three factor levels −1, 0 and +1 for the lower, middle and upper zones, respectively.

### 2.2.4. Model validation

To check whether the estimated regression model was a valid approximation, we proceeded as follows:

#### 2.2.4.1. Calculation of the adjusted squared multiple correlation coefficient.

$$R_{adj}^2 = 1 - (1 - R^2) \frac{(n - 1)}{(n - p)},$$

where R is the usual multiple correlation coefficient, n is the number of observations (here n = 6561) and p is the number of parameters in the model (in Eq. (4), p = 1+32+496+32 = 561).

#### 2.2.4.2. Examination of the residuals. We considered the DFFITSᵢ statistic defined as:

$$DFFITS_i = \frac{(\hat{y}_i - \hat{y}_{(i)})}{s_{(i)} \sqrt{h_i}},$$

where ŷᵢ is the prediction for the ith observation and ŷᵢ(i) is the ith prediction but without the ith observation, sᵢ(i) is the estimated standard deviation after deleting the ith

<figure><img src="image_7.png" /><figcaption>Fig. 6. Scatter plot of the shape and scale parameters for cereals obtained by bootstrap showing the three experimental domains.</figcaption></figure>
<table><caption>Table 5
Levels of the inputs of the sensitivity analysis for the lower part (zone I)</caption><thead><tr><th rowspan="2">Food</th><th colspan="6">Factor range</th></tr><tr><th colspan="3">Consumption</th><th colspan="3">Contamination</th></tr><tr><th></th><th>-1</th><th>0</th><th>1</th><th>-1</th><th>0</th><th>1</th></tr></thead><tbody><tr><td>Cereals</td><td>r 2.148<br>λ 0.395</td><td>2.215<br>0.404</td><td>2.282<br>0.412</td><td>0.257<br>0.306</td><td>0.292<br>0.359</td><td>0.327<br>0.412</td></tr><tr><td>Raisins</td><td>r 0.0524<br>λ 0.492</td><td>0.165<br>1.574</td><td>0.278<br>2.656</td><td>0.2822<br>0.3055</td><td>0.375<br>0.5305</td><td>0.468<br>0.755</td></tr><tr><td>Other<br>dried fruit</td><td>r 0.1691<br>λ 0.500</td><td>0.328<br>1.116</td><td>0.487<br>1.732</td><td>0.402<br>1.189</td><td>0.582<br>1.685</td><td>0.762<br>2.180</td></tr><tr><td>Pork</td><td>r 1.381<br>λ 0.855</td><td>1.455<br>0.930</td><td>1.529<br>1.005</td><td>0.547<br>1.722</td><td>0.577<br>1.808</td><td>0.607<br>1.894</td></tr><tr><td>Poultry</td><td>r 1.470<br>λ 1.859</td><td>1.535<br>2.050</td><td>1.600<br>2.241</td><td>3.705<br>46.784</td><td>3.891<br>51.17</td><td>4.077<br>55.56</td></tr><tr><td>Fruit juices</td><td>r 0.5618<br>λ 0.1051</td><td>0.6466<br>0.1254</td><td>0.7314<br>0.1457</td><td>0.0693<br>0.2486</td><td>0.1078<br>0.4233</td><td>0.1464<br>0.598</td></tr><tr><td>Wines</td><td>r 0.323<br>λ 0.777</td><td>0.485<br>2.525</td><td>0.647<br>4.273</td><td>0.294<br>1.416</td><td>0.3161<br>1.753</td><td>0.338<br>2.090</td></tr><tr><td>Coffee</td><td>r 0.342<br>λ 4.122</td><td>0.508<br>5.888</td><td>0.673<br>7.654</td><td>0.271<br>0.1796</td><td>0.2918<br>0.1928</td><td>0.3125<br>0.2059</td></tr></tbody></table>

observation, and $h_i$ is the ith diagonal of the projection matrix for the predictor space (see Belsley et al., 1980 for more details). As the distribution of the model residuals shows an approximately normal distribution, a possible threshold is the size-adjusted cut-off recom-

<table><caption>Table 6
Levels of the inputs of the sensitivity analysis for the middle part (zone II)</caption><thead><tr><th rowspan="2">Food</th><th colspan="6">Factor range</th></tr><tr><th colspan="3">Consumption</th><th colspan="3">Contamination</th></tr><tr><th></th><th>-1</th><th>0</th><th>1</th><th>-1</th><th>0</th><th>1</th></tr></thead><tbody><tr><td>Cereals</td><td>r 2.3690<br>λ 0.4688</td><td>2.5944<br>0.4796</td><td>2.8198<br>0.4904</td><td>0.3987<br>0.6210</td><td>0.4395<br>0.6740</td><td>0.4802<br>0.7260</td></tr><tr><td>Raisins</td><td>r 0.5122<br>λ 3.6110</td><td>0.7012<br>5.4840</td><td>0.8902<br>7.3570</td><td>0.2822<br>0.3055</td><td>0.3750<br>0.5305</td><td>0.4680<br>0.7550</td></tr><tr><td>Other<br>dried fruit</td><td>r 0.4042<br>λ 1.5492</td><td>0.5818<br>2.2862</td><td>0.7594<br>3.0230</td><td>0.8160<br>3.0400</td><td>0.9980<br>3.5020</td><td>1.1790<br>3.9640</td></tr><tr><td>Pork</td><td>r 1.6320<br>λ 0.5635</td><td>1.7140<br>1.1269</td><td>1.7950<br>1.6903</td><td>1.0240<br>3.5477</td><td>1.0850<br>3.5907</td><td>1.1450<br>3.6340</td></tr><tr><td>Poultry</td><td>r 1.7590<br>λ 2.2890</td><td>1.8310<br>2.4950</td><td>1.9030<br>2.7020</td><td>4.6150<br>59.0400</td><td>4.8100<br>64.1200</td><td>5.0000<br>69.2000</td></tr><tr><td>Fruit juices</td><td>r 0.8450<br>λ 0.1862</td><td>0.9388<br>0.2087</td><td>1.0330<br>0.2312</td><td>0.0693<br>0.2486</td><td>0.1078<br>0.4233</td><td>0.1464<br>0.5980</td></tr><tr><td>Wines</td><td>r 0.5895<br>λ 3.5740</td><td>0.7645<br>5.5940</td><td>0.9390<br>7.6140</td><td>0.3657<br>1.7800</td><td>0.4127<br>2.4040</td><td>0.4596<br>3.0300</td></tr><tr><td>Coffee</td><td>r 0.5862<br>λ 7.0300</td><td>0.7842<br>8.9610</td><td>0.9820<br>10.8920</td><td>0.3861<br>0.3421</td><td>0.4310<br>0.3860</td><td>0.4760<br>0.4300</td></tr></tbody></table>

<table><caption>Table 7
Levels of the inputs of the sensitivity analysis for the upper part (zone III)</caption><thead><tr><th rowspan="2">Food</th><th colspan="6">Factor range</th></tr><tr><th colspan="3">Consumption</th><th colspan="3">Contamination</th></tr><tr><th></th><th>-1</th><th>0</th><th>1</th><th>-1</th><th>0</th><th>1</th></tr></thead><tbody><tr><td>Cereals</td><td>r 2.8838<br>λ 0.544</td><td>2.9738<br>0.5552</td><td>3.0637<br>0.566</td><td>0.552<br>0.9363</td><td>0.587<br>0.9893</td><td>0.622<br>1.0425</td></tr><tr><td>Raisins</td><td>r 1.124<br>λ 8.312</td><td>1.237<br>9.394</td><td>1.350<br>10.476</td><td>1.305<br>3.095</td><td>1.404<br>3.339</td><td>1.502<br>3.583</td></tr><tr><td>Other<br>dried fruit</td><td>r 0.676<br>λ 2.84</td><td>0.835<br>3.456</td><td>0.994<br>4.072</td><td>1.233<br>4.825</td><td>1.413<br>5.3202</td><td>1.593<br>5.815</td></tr><tr><td>Pork</td><td>r 1.899<br>λ 1.2488</td><td>1.973<br>1.3241</td><td>2.047<br>1.399</td><td>1.5625<br>5.287</td><td>1.5925<br>5.373</td><td>1.6225<br>5.459</td></tr><tr><td>Poultry</td><td>r 2.062<br>λ 2.749</td><td>2.127<br>2.940</td><td>2.192<br>3.131</td><td>5.543<br>72.68</td><td>5.729<br>77.07</td><td>5.915<br>81.44</td></tr><tr><td>Fruit juices</td><td>r 1.146<br>λ 0.272</td><td>1.231<br>0.292</td><td>1.316<br>0.312</td><td>0.238<br>0.997</td><td>0.2689<br>1.141</td><td>0.299<br>1.285</td></tr><tr><td>Wines</td><td>r 0.882<br>λ 6.915</td><td>1.044<br>8.663</td><td>1.206<br>10.41</td><td>0.487<br>2.718</td><td>0.509<br>3.055</td><td>0.531<br>3.392</td></tr><tr><td>Coffee</td><td>r 0.895<br>λ 10.26</td><td>1.060<br>12.03</td><td>1.225<br>13.796</td><td>0.549<br>0.565</td><td>0.570<br>0.578</td><td>0.591<br>0.591</td></tr></tbody></table>

mended by Belsley et al. (1980) equal to $2\sqrt{p/n}$. Therefore, a value of DFFITS$_i$ above this cut-off indicates an influential residual $e_i$. The DFFITS$_i$ was calculated by SAS software (Proc REG).

### 2.2.4.3. Test experiments.
A final test was performed to evaluate the model validation: the comparison between the predicted outputs for K, typically 10,000, new input combinations randomly sampled in the three squared experimental domains and the corresponding K simulated outputs. Then, we check if the K residuals are Gaussian and if their variation ranges are reasonable. We calculate the percentages $\delta e_i$:

$$ \delta e_i = \frac{e_i}{y_i} \times 100 = \frac{\hat{y}_i - y_i}{y_i} \times 100, \quad (5) $$

where $y_i$ is the observation and $\hat{y}_i$ the value predicted by the model to see the percentages of discrepancy between predicted and observed values for the test experiments.

## 3. Results

We fitted the model (4) to the 6561 experiments for the two responses $\hat{Q}_{0.95}$ and $\hat{Q}_{0.99}$ in the three zones as described in Section 2.2.3, after the factors $z_j$ had been
<table><caption>Table 8
$R_{\text{adj}}^2$ statistic for main models (95th quantile and zone II)</caption><thead><tr><th>Model</th><th>Number of parameters</th><th>R<sup>2</sup><sub>adj</sub></th></tr></thead><tbody><tr><td>Model 1: y = β<sub>0</sub> + Σ<sub>j=1</sub><sup>32</sup> β<sub>jzj</sub> + Σ<sub>j=1</sub><sup>31</sup> Σ<sub>k=j+1</sub><sup>32</sup> β<sub>jk</sub>z<sub>j</sub>z<sub>k</sub> + Σ<sub>j=1</sub><sup>32</sup> β<sub>jj</sub>Z<sup>2</sup><sub>j</sub> + ε</td><td>560</td><td>0.9715</td></tr><tr><td>Model 2: y = β<sub>0</sub> + Σ<sub>j=1</sub><sup>32</sup> β<sub>jzj</sub> + Σ<sub>j=1</sub><sup>32</sup> β<sub>jj</sub>Z<sup>2</sup><sub>j</sub> + ε</td><td>64</td><td>0.9682</td></tr><tr><td>Model 3: y = β<sub>0</sub> + Σ<sub>j=1</sub><sup>32</sup> β<sub>jzj</sub> + Σ<sub>j=1</sub><sup>31</sup> Σ<sub>k=j+1</sub><sup>32</sup> β<sub>jk</sub>z<sub>j</sub>z<sub>k</sub> + ε</td><td>528</td><td>0.9695</td></tr><tr><td>Model 4: y = β<sub>0</sub> + Σ<sub>j=1</sub><sup>32</sup> β<sub>jzj</sub> + ε</td><td>32</td><td>0.9664</td></tr></tbody></table>

appropriately coded (use of orthogonal polynomials and scaling). Table 8 gives the values of the $R_{\text{adj}}^2$ of the four main models in the middle part for the 95th quantile. According to these $R_{\text{adj}}^2$, the model with only the main effects was chosen in the first step (model 4 in Table 8). Indeed, the model has just 32 terms and its value of $R_{\text{adj}}^2$ is still very high ($R_{\text{adj}}^2 = 0.97$). The results being equivalent for the 99th quantile and the $R_{\text{adj}}^2$ values being slightly smaller (about 3% less), the model with only the 32 main effects was also chosen for the 99th quantile. For both quantiles, the results were similar in the other two zones.

In the second step, with a view to parsimony, for each zone, variables were selected considering their explained sum of squares. The variables with the higher sum of squares were kept in the final model. Table 9 gives the final models with their parameter estimations for the 95th and 99th quantiles and for the

three zones. The model for the 95th quantile includes only four terms (the main effects) in the lower and upper zones and five terms in the middle zone. The results are similar for the 99th quantile with a slightly smaller percentage of explained variability for each model. The results show a major influence of the parameters of the distributions on cereals (contamination and consumption). In the middle part, a slight influence of the pork consumption scale parameter appears. We noted that both the orthogonality of the simulation designs and the above-mentioned coding enable us to compare rigorously the influences of all the food distribution parameters. We obtained the marginal effect of each parameter and all of the parameters could be compared with each other because they have the same scale. On the other hand, the constant terms of the models of Table 9 represent the quantile means inside the cuboidal zones: for the 95th quantile they are 28.1, 20.1 and 16.7, and for the

<caption>Table 9</caption>

Polynomial models of the 95th and 99th quantiles in the three zones of the 95% concentration ellipse with $z_1$ = cereals contamination scale parameter, $z_2$ = cereals contamination shape parameter, $z_3$ = cereals consumption shape parameter, $z_4$ = cereals consumption scale parameter, and $z_5$ = pork consumption scale parameter

<table><thead><tr><th>95% concentration elliptical zone</th><th>Percentile</th><th>Final model</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="2">I (Lower part)</td><td>95th</td><td>y&#770; = 28.1-8.2z<sub>1</sub>+5.4z<sub>2</sub> + 1.2z<sub>3</sub>-z<sub>4</sub></td><td>0.97</td></tr><tr><td>99th</td><td>y&#770; = 61.5-18.2z<sub>1</sub>+7.3z<sub>2</sub> + 2.7z<sub>3</sub>-2.2z<sub>4</sub></td><td>0.96</td></tr><tr><td rowspan="2">II (Middle part)</td><td>95th</td><td>y&#770; = 20.1-3.1z<sub>1</sub>+2.6z<sub>2</sub> + 2.6z<sub>3</sub>-0.7z<sub>4</sub>-0.5z<sub>5</sub></td><td>0.97</td></tr><tr><td>99th</td><td>y&#770; = 39.1-6z<sub>1</sub>+3.6z<sub>2</sub> + 4.9z<sub>3</sub>-1.5z<sub>4</sub>-0.5z<sub>5</sub></td><td>0.94</td></tr><tr><td rowspan="2">III (Upper part)</td><td>95th</td><td>y&#770; = 16.7-1.7z<sub>1</sub>+1.3z<sub>2</sub> + 0.7z<sub>3</sub>-0.5z<sub>4</sub></td><td>0.94</td></tr><tr><td>99th</td><td>y&#770; = 30.4-3.2z<sub>1</sub>+1.9z<sub>2</sub> + 1.3z<sub>3</sub>-z<sub>4</sub></td><td>0.89</td></tr></tbody></table>
99th quantile they are 61.5, 39.1 and 30.4 according to the upper, middle and lower zone, respectively. We can compare these values to the value given in Gauchi (in press): for the 95th and 99th quantiles, the following 95% confidence intervals were calculated [14.07; 29.56] and [28.08; 76.92], respectively. We noted that the quantile means obtained with this SA lie inside the 95% confidence intervals. Thus, we assume that the SA outputs are reasonable values.

The distributions of the residuals of the different models did not indicate outliers. Only a few DFFITS measures slightly exceeded the threshold as defined in Section 2.2.4.2. Ten thousand new combinations of inputs have been randomly sampled inside the three experimental domains. In the three zones, the distributions of the residuals of each model are approximately Gaussian and have moderate variation ranges. For example, in the middle part for the 95th quantile, the estimated standard deviation of the model residuals is estimated at 2.92 (for the 10,000 test experiments). It is higher than the estimated standard deviation (0.37) of the residuals on the 6561 experiments from which the model has been chosen, but their variation range is reasonable because 95% of these residuals have percentage, as defined in Eq. (5), between ±28%. The discrepancy between the predicted and observed values for the test experiments is acceptable. Fig. 7 represents the histogram of the 10,000 residuals for the 95th quantile in the middle zone.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Histogram of the 10,000 residuals for the 95th quantile for the 10,000 new observations in the middle part.</figcaption></figure>

## 4. Discussion

On the basis of the observation, well known to risk analysts, that to obtain a good estimation of the high quantiles of exposure is a difficult operation when it is based on data which are both sparse and characterised by very asymmetrical distributions, we undertook a sensitivity analysis of the 95th and 99th quantiles of exposure to OA. We conducted an analysis of their sensitivity to probability density parameters adjusted to consumption and contamination data. That is, we measured the influence of the shapes of the consumption and contamination distributions of the studied foods to highlight those foods that possibly require particular investigation in order to better estimate the high exposure quantiles.

One of the original ideas of this work was to generate a realistic variability of data available by a parametric bootstrap method. This enabled us to construct variation intervals of the density parameters concerned, intervals presumed to be representative of the child population. Another original idea was to take into account the natural correlation of the parameters of the Gamma laws using concentration ellipses. However, we have not explored all the variation domains of the inputs: we studied the sensitivity of the 95th and 99th quantiles in three subdomains of cuboidal symmetry distributed judiciously in the variation domain (lower, middle and upper zone). The advantage of these subdomains is that they make it possible to set up orthogonal experimental designs and to take into account a large variation of the quantiles (see the constant terms of the models in the three zones). On these subdomains, it appears very clearly that a small number of parameters have a great influence on estimating these quantiles compared with the very large number of parameters envisaged at the beginning. These include mainly the parameters of the laws of consumption and contamination of the category cereals. Cereals are known by epidemiologists and toxicologists to play a major role in health problems concerning OA in Europe.

It is difficult to take into account all the volume of the variation domain of the inputs (for example, because of the non-convexity of this domain). The necessary tools are derived from the theory of optimum experimental designs. This could possibly improve our work. However, it would be appropriate beforehand to
guarantee the pertinence of our domain of variation (obtained by bootstrap sampling) ideally with data from further surveys.

Finally, it is difficult to compare our results with other results because, to our knowledge, no other research group has conducted a study on contamination by OA. The other estimations of the OA exposure quantiles found elsewhere are very empirical as indicated previously in Gauchi (in press) and the sensitivity analysis presented here is the first SA on OA exposure quantiles.

## Acknowledgements

We thank P. Verger and J.C. Leblanc (INRA, Scientific Directorate for Human Nutrition and Food Safety, Paris) for their important preliminary data handling and providing the consumption survey data and the contamination raw data.

## References

ASPCC, 1999. Observatoire des Consommations Alimentaires: Etude des distributions statistiques des consommations de huit groupes d'aliments susceptibles de contenir de l'Ochratoxine A. Internal Technical Report of CREDOC, for INRA, based on data from ASPCC survey, April 1999.

Belsley, D.A., Kuh, E., Welsch, R.E., 1980. Regression Diagnostics John Wiley & Sons, New York.

Box, G.E.P., Hunter, W.G., Hunter, J.S., 1978. Statistics For Experimenters.

Codex Committee on Food Additives and Contaminants (CCFAC), 1999. Position paper on ochratoxin A. Codex Committee on Food Additives and Contaminants, Thirty-first Session, The Hague, The Netherlands, 22–26 March 1999.

D'Agostino, R.B., Stephens, M.A., 1986. Goodness-of-Fit Techniques. Marcel Dekker, New York.

Efron, B., Tibshirani, R.J., 1993. An Introduction to the Bootstrap Chapman & Hall, London.

Gauchi, J.P., 2002. Quantitative assessment of exposure to the mycotoxin ochratoxin A in food. Risk Anal., in press.

Iman, R.L., Conover, W.J., 1982. A distribution-free approach to inducing rank correlation among input variables. Commun. Stat., Simul. Comput. 11, 311–334.

JECFA (Joint FAO/WHO Expert Committee on Food Additives), 1995. Ochratoxin A Report: TRS 859-JECFA 44/35 Tox monograph: FAS 35-JECFA 44/363. World Health Organization, Geneva.

Pfohl-Leszkowick, A., 1999. Les Mycotoxines dans L'alimentation: Évaluation et Gestion du Risque. Tec&Doc Editor, Paris, France Chap. 3.

SCOOP task 3.2.7. Report of the European Communaut y, French contribution, in press. Origin of data: Direction Générale de l'Alimentation (Ministère de l'Agriculture et de la Pêche), Direction Gènèrale de la Concurrence, de la Consommation et de la Répression des Fraudes (Ministère de l'Economie, des Finances et de l'Industrie).

Verger, P., 1999. Exposure assessment to certain contaminants, working group «Contaminants». Internal Technical Report, March 1999. INRA, Paris.

Vose, D., 1996. Quantitative Risk Analysis—A Guide to Monte Carlo Simulation Modelling Wiley, Chichester, UK.