Modelling unit-to-unit variation¶
The basic model for an acute exposure assessment assumes that the concentration of the substance displays the variation of residues between units in the marketplace. In general, both monitoring data and controlled field trial data are obtained using composite samples. As a result some of the unit-to-unit variation is averaged out. The model for unit variability aims to adjust the composite sample mean such that sampled concentrations represent the originally unit-to-unit variation of the units in the composite sample.
MCRA offers three distributions to sample from:
the beta distribution,
and the bernoulli distribution.
The beta distribution simulates values for a unit in the composite sample. It requires knowledge of the number of units in a composite sample and of the variability between units.
The lognormal distribution simulates values for a new unit in the batch. It requires only knowledge of the variability between units.
The bernoulli distribution is considered as a limiting case of the beta distribution when knowledge of the variability between units is lacking and only the number of units in the composite sample is known. For the beta and lognormal distribution, estimates of unit variability are either realistic (no censoring at the value of the monitoring residue) or conservative (unit values are left-censored at the value of the monitoring residue). For the lognormal distribution sampled concentrations have no upper limit. Whereas for the beta distribution, sampled concentration values for a unit are never higher than the monitoring residue times the number of units in the composite sample.
Variability between units is specified using a variability factor v (defined as 97.5th percentile divided by mean) or a coefficient of variation cV (standard deviation divided by mean). Following FAO/WHO recommendations, the default variability factor v = 1 for small crops (unit weight < 25 g). For large crops (unit weight ≥ 25 g) v = 5. For foods which are processed in large batches, e.g. juicing, marmalade/jam, sauce/puree, bulking/blending the variability factor v = 1 is proposed.
Estimation of intake values using the concept of unit variability¶
A composite sample for food k is composed of nuk units with nominal (whole food/RAC) unit weight wuk. The weight of a composite sample is wmk=nuk⋅wuk with mean residue value cmk.
For each iteration i in the MC-simulation, obtain for each food k a simulated intake xik, and a simulated composite sample concentration cmik.
Calculate the number of unit intakes nuxik in xik (round upwards) and set weights wikl equal to unit weight wuk, except for the last partial intake, which has weight wikl=xik−(nuxik−1)wuk .
For the beta or bernoulli distribution: draw nuxik simulated values bcikl from a beta or bernoulli distribution. Calculate concentration values as cikl=bcikl⋅cmik,max=bcikl⋅cmik⋅nuk=svfikl⋅cmik, where nuk is the number of units in a composite sample of food k, and svfikl is the stochastic variability factor for this simulated unit, i.e. the ratio between simulated concentration cikl and the simulated composite sample concentration cmik. Sum to obtain the simulated concentration in the consumed portion:
For the lognormal distribution: draw nuxik simulated logconcentration values lcikl from a normal distribution with (optional) a biased mean μ=ln(cmik) or (default) unbiased mean μ=ln(cmik)−1/2σ2 and standard deviation σ. Calculate concentration values as
where svfikl is the stochastic variability factor for this simulated unit, i.e. the ratio between simulated concentration cikl and the simulated composite sample concentration cmik. Back transform and sum to obtain the simulated concentration in the consumed portion:
For cumulative exposure assessments, a sensitivity analysis may be performed by specifying a full correlation between concentrations from different substances on the same unit. As a result, high (or low) concentrations from different substances occur together on the same unit. In MCRA, for each unit the random sequence is repeatedly used to generate concentration values for all substances.
Beta distribution¶
Under the beta model simulated unit values are drawn from a bounded distribution on the interval (0, cmax) with cmax=nuk⋅cmk. The standard beta distribution is defined on the interval (0, 1) and is usually characterised by two parameters a and b, with a>0, b>0 (see e.g. Mood et al. (1974)). Alternatively, it can be parameterised by the mean
and the variance
or, as applied in MCRA, by the mean μ and the squared coefficient of variation
For the simulated unit values in each iteration of the program we require an expected value cmk. This scales down to a mean value μ=cmk/cmax=1/nuk in the (standard) beta distribution. From this value for μ and an externally specified value for cvk the parameters a and b of the beta distribution are calculated as:
and
From the second formula it can be seen that cvk should not be larger than √nuk−1 in order to avoid negative values for b. When the unit variability is specified by a variability factor
instead of a coefficient of variation cvk then MCRA applies a bisection algorithm to find a such that the cumulative probability
for b=a(nuk−1).
Sampled values from the beta distribution are rescaled by multiplication with cmmax to unit concentrations cijk on the interval (0,cmmax).
Lognormal distribution¶
The lognormal distribution is characterised by μ and σ, which are the mean and standard deviation of the log-transformed concentrations. The unit log-concentrations are drawn from a normal distribution with mean μ=ln(cmik)−1/2σ2. The coefficient of variation cv is turned into the standard deviation σ on the log-transformed scale with:
The variability factor is defined as the 97.5th percentile of the concentration in the individual measurements divided by the corresponding mean concentration seen in the composite sample. A variability factor v is converted into the standard deviation σ as follows:
with μ and σ representing the mean and standard deviation of the log-transformed concentrations. So
Solving for σ gives:
with roots for σ according to:
The smallest positive root is taken as an estimate for σ.
Bernoulli distribution¶
The bernoulli model is a limiting case of the beta model, which can be used if no information on unit variability is available, but only the number of units in a composite sample is known (see van der Voet et al. 2001). As a worst case approach we may take the coefficient of variation cv as large as possible. When cv is equal to the maximum possible value √nuk−1, the (unstandardised) beta distribution simplifies to a bernoulli distribution with probability
or
for the value 0 and probability
or
for the value cmax=nuk⋅cmk.
In MCRA values 0 are actually replaced by cmk, to keep all values on the conservative side. For example, with nuk = 5, there will be 80% probability at cijk=cmk and 20% probability at cijk=cmax. When the number of units nuk in the composite sample is missing, the nominal unit weight wuk is used to calculate the parameter for unit variability.