Statistical Analysis of Chemical Transformation Kinetics Using

ARTICLE pubs.acs.org/est

Statistical Analysis of Chemical Transformation Kinetics Using Markov-Chain Monte Carlo Methods Linus G€orlitz,† Zhenglei Gao,‡ and Walter Schmitt*,‡ † ‡

Bayer Technology Services, Leverkusen, Germany Bayer CropScience, Monheim, Germany ABSTRACT: For the risk assessment of chemicals intentionally released into the environment, as, e.g., pesticides, it is indispensable to investigate their environmental fate. Main characteristics in this context are transformation rates and partitioning behavior. In most cases the relevant parameters are not directly measurable but are determined indirectly from experimentally determined concentrations in various environmental compartments. Usually this is done by fitting mathematical models, which are usually nonlinear, to the observed data and such deriving estimates of the parameter values. Statistical analysis is then used to judge the uncertainty of the estimates. Of particular interest in this context is the question whether degradation rates are significantly different from zero. Standard procedure is to use nonlinear least-squares methods to fit the models and to estimate the standard errors of the estimated parameters from Fisher’s Information matrix and estimated level of measurement noise. This, however, frequently leads to counterintuitive results as the estimated probability distributions of the parameters based on local linearization of the optimized models are often too wide or at least differ significantly in shape from the real distribution. In this paper we identify the shortcoming of this procedure and propose a statistically valid approach based on Markov-Chain Monte Carlo sampling that is appropriate to determine the real probability distribution of model parameters. The effectiveness of this method is demonstrated on three data sets. Although it is generally applicable to different problems where model parameters are to be inferred, in the present case for simplicity we restrict the discussion to the evaluation of metabolic degradation of chemicals in soil. It is shown that the method is successfully applicable to problems of different complexity. We applied it to kinetic data from compounds with one and five metabolites. Additionally, using simulated data, it is shown that the MCMC method estimates the real probability distributions of parameters well and much better than the standard optimization approach.

’ INTRODUCTION For the characterization of the environmental fate of pesticides and other chemicals frequently compound properties have to be determined indirectly from experimental data, because they cannot be measured directly. Examples for this are transformation rates of the substances and their metabolites in various environmental compartments. Such rates can be inferred from respective transformation studies in which a substance is applied to the system, e.g., soil, water, etc., in which its fate shall be determined. The concentrations of the compound and eventually also of its degradation products are then measured after different incubation times. For inferring the compound properties from the resulting concentration vs time series’ then usually a mathematical model that contains the properties of interest as parameters is fitted to the observed data. The parameter values leading to the best fit are considered the actual values of the compound properties. Usually fitting of the model is done using nonlinear optimization methods. The standard approach for this is nonlinear least-squares (LS) optimization. Due to the nature of the optimization problem the inferred parameter values cannot be seen as absolute values but are the most probable values in a r 2011 American Chemical Society

distribution of values. The width of this distribution is determined by different factors as e.g., the variation of the observed data and the mathematical structure of the model itself. The latter might, for example, impose strong correlation between parameters which prevents the identifiability of a certain parameter value. For a sound assessment of substance properties inferred by model optimization it is therefore indispensable not only to determine the most likely value but also to assess its uncertainty. In case of degradation rates it is, for example, an important question whether a degradation is observed at all or if the nonzero degradation rate determined is so uncertain that with high probability the rate could also be zero. The question of uncertainty is of particular interest if the parameter values are determined to serve as substance properties used in a regulatory environment, e.g., the approval process for pesticides. For the case of kinetic evaluation of degradation or dissipation data and determination of respective kinetic data, Received: December 16, 2010 Accepted: April 13, 2011 Revised: April 11, 2011 Published: April 28, 2011 4429

dx.doi.org/10.1021/es104218h | Environ. Sci. Technol. 2011, 45, 4429–4437

Environmental Science & Technology i.e., degradation and formation rates, for use in pesticide risk assessments in Europe a guidance is given by the report of the FOCUS Group on Kinetic Evaluations.1 Apart from recommendations for kinetic models to be used in such assessments this guidance also requests to show that degradation rates are nonzero. It prescribes the use of either a p-value of a t test or a quantile of the normal approximation of the probability distribution over the respective degradation rate around the most probable parameter vector. There is no guidance given on which methods to use for inferring the parameter values. However, it is common practice in the field to use nonlinear least-squares methods with Gauss or LevenbergMarquardt algorithms.2 Applying the recommended statistical evaluations of the parameter distributions these methods inherently use a symmetric approximation around the parameter value giving the best model approximation to the observed data. If the applied models are linear this approximation will globally coincide with the distribution,3 and all decisions based on it are correct. Models for chemical reaction kinetics, the transformation models described here, are, however, highly nonlinear. In this case the aforementioned approximation is only locally valid around the most probable parameter values. For this reason it may be that conclusions drawn from such derived distributions are unrealistic and misleading. This is particularly the case for the determination of the probability of a parameter to be nonzero. For this the marginal distribution may be more relevant than the central distribution for which the approximation is valid. We propose the Markov-Chain Monte Carlo (MCMC) methodology for determination of probability distributions of model parameters. This computationally intensive approach consists of an algorithm generating samples from a distribution without making explicit assumptions on the type or parametric family of the distribution. Posthoc analysis can then be performed using computational statistics. It can thus be used to perform statistical analysis without any approximations and is thus capable of generating correct results in any situation.

’ METHODS Data Sets. For demonstration purposes in the following we use different experimental and artificial data sets. All experimental data result from standard aerobic degradation studies with soil following the respective OECD guidance.4 In these studies radiolabeled parent compound is applied to samples of natural soil and incubated under defined and controlled moisture and temperature conditions. At different time points the soil samples are extracted, and the radioactive solutes in the extraction solution are quantified and analyzed for their fraction of parent compound and metabolites formed. The samples analyzed at different time points are independent in that they are incubated in separate testing vessels. Moreover for each time point two samples were prepared. In all cases the full material balance is checked including nonextractable residues in soil, and the results of analysis are given in percentage of applied radioactivity. One data set comprises only data on the applied substance and one metabolite, while the second set contains data on five metabolites in addition to the applied compound. Because the data here just serve as more or less arbitrarily chosen realistic examples which show certain purely data-related traits needed for the intended demonstration we refrain from describing the compounds behind the data and their properties in detail. The artificial data were simulated using the same equations as also used in the

ARTICLE

Figure 1. Schematic representation of the type of models used. This model is constructed for situations where one parent substance can be either metabolized or degraded. The active metabolite can only be degraded. This is the simplest model used for degradation data analysis. Much more complex models with several metabolites and complex metabolization pathways must often be used.

model for their subsequent evaluation. Simulations were run with mean values for the model parameters, and an error was subsequently introduced by adding random values to each data point drawn from normal distributions with adequate standard deviation. Kinetic Reaction Models. In the following we assume that fj(θ,ti) is the concentration value at time ti for species j in a transformation model (e.g., like the one shown in Figure 1) and θ contains all parameters of the model (i.e., degradation rates, formation fractions, and initial concentrations). The measured data set is denoted by D, and p always denotes a probability distribution. Conceptual transformation models consist of boxes and links. Each box represents either the parent drug or any of its active metabolites. A link between two boxes constitutes that the species in the starting box of the link is metabolized into the target species. According to the FOCUS guidance1 the kinetic model assigned to a link can be chosen to be single or multiple first-order, or biphasic. Any species is furthermore degraded into inactive substances which are typically represented by one sink. In the case of the model shown in Figure 1 and linear degradation kinetics the concentration fj({kp,km,M0,c},ti) of species j at time ti is given by f1 ðfkp , km , M0 , cg, tj Þ ¼ M0 expkp tj f2 ðfkp , km , M0 , cg, tj Þ ¼ cM0 kp

expkm tj expkp tj kp km

with degradation rates kp (parent) and km (metabolite), initial parent concentration M0, and formation fraction c. For a description of the more complex metabolic pathway underlying the second, more complex experimental data set the above scheme has respectively been extended by adding metabolites and reactions. The structure of the metabolic transformation model was 4430

dx.doi.org/10.1021/es104218h |Environ. Sci. Technol. 2011, 45, 4429–4437

Environmental Science & Technology

ARTICLE

chosen to appropriately consider the chemical and kinetic properties of the metabolic pathway. This process of model development is not described here in detail because it is beyond the scope of the present work. In this complex case the model equations cannot be given as analytical equations but were implented as the respective differential equations of the reactions. Bayesian Modeling. There are two views on parameter identification: optimization-based approaches and the Bayesian framework. Both approaches rely on the same statistical modeling principles, but optimization-based parameter identification finds the parameter vector θ which best explains the measured data. Only this parameter vector is later used to perform any analysis with the transformation model. It ignores the fact that almost never a single parameter vector is clearly identified. Often only a parameter region where the true value will lie can be outlined. Uncertainties of parameter estimates can be estimated using the Fisher information matrix which yields a local approximation of the “real” distribution.5 A second disadvantage of this approach is that prior knowledge on the parameters (e.g., upper and lower bounds, covariate restrictions) is hard to include in this approach. A conceptually different approach to dealing with uncertain parameters in models which allows a theoretically sound way to combine prior information with the data is to use the Bayesian framework. Here, any unknown parameter is modeled as a random variable and no single parameter vector but the probability distribution over all parameter vectors is sought. General information on Bayesian methodology, its applications, and many examples can be found in refs 6 and 7. Bayes’ theorem, given by

trivial lower and upper bounds). The uncertainties {σj} which represent the measurement error for species j are also unknown a priori. In the Bayesian framework they are also modeled as random variables and assigned a prior distribution. It has to be noted that in contrast to classical approaches no single parameter vector is selected to contain all information in the model, but the whole posterior distribution of the parameter vector given the data is estimated. Consequently, the model prediction ^f (ti) at time ti cannot be calculated using only one parameter vector. It is given by the expectation value of the model if all parameter vectors are used ^f ðti Þ ¼ Epðθ, fσ gjDÞ ½fj ðti , θÞ j j Z ¼ fj ðθ, ti Þ 3 pðθ, fσ j gjDÞdðθ, ðfσ j gÞÞ θ

with weight given by the posterior probability of the used parameter vector θ. According to the FOCUS guideline, to show that the degradation rate for a substance is nonzero, it suffices to show that the 5%- or 10%-quantile of the posterior distribution for the degradation rate (for mother substance and metabolites, respectively) is larger than zero. In the Bayesian case this has to be done by estimating the quantile of the marginal distribution of the degradation rate θi. The 5%-quantile for a rate θi, for example, is defined by Z xi pðθi jDÞdθi e 0:05 ð4Þ inf xi j ¥

Z pðθi jDÞ ¼

provides a way to combine prior information on the parameter θ and closeness to the measurements D = {(xij,ti)}i=1,...,N given by the likelihood function p(D|θ). Here xij denotes the measurement of species j at time point ti. It is exactly the negative of the logarithm of the likelihood function which is used in optimization approaches. Thus, Bayesian as well as optimization-based approaches will identify the same “optimal” parameter vector. The prior distribution can be very complex but is often assumed to be flat within the admissible parameter range and zero elsewhere. The optimal formulation of the likelihood function for transformation models is explained in detail in ref 8. In this work likelihood function and prior distributions are given as N Y j 2 2 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 exp0:5 3 ððxi fj ðθ, ti ÞÞ Þ=σj i ¼ 1 2 3 π 3 σ 2j Q Q ¼ pðθi Þ 3 pðσ j Þ i j 8 1 > < if ai e θi e bi 1 ¼ bi ai pðσj Þ 2 > σ : 0 else j

pðDjθ, fσj gÞ ¼ pðθ, fσ j gÞ pðθi Þ

ð2Þ with unequal variances for different species and ai, bi denoting the lower and upper limit of θi, respectively. Thus no prior information on the value of the model parameters is assumed (apart from

ð3Þ

~ , fσ j gÞjθ ~ ¼ θi g fðθ i

~ fσj gÞjDÞdðθ, ~ fσj gÞ pððθ,

ð5Þ

In contrast to optimization-based approaches, where the quantiles are estimated from a distribution where all other parameters are fixed at their optimal values, this approach judges the information for the degradation rates by taking all possible parameter values into account. The explicit dependence of predictions and marginals on {σj} will now be dropped from notation whenever no confusion can arise. Prediction and Quantile Estimation. In general, high-dimensional integrals like eq 3 and eq 4 are analytically intractable, and thus estimators for these quantities are required. The most widespread approach for numerical intergration of high dimensional functions is to use Monte Carlo simulations. To compute aforementioned integrals, a set of K independent parameter vectors Θ = {θk}k∈{1,...,K} sampled from the posterior distribution shown in eq 1 is required, with θk = {θk,l}l denoting a parameter vector. Then, because eq 3 can be written as an expectation, its value can be estimated by K X ^f ðti Þ ¼ 1 f ðθk , ti Þ K i¼1

ð6Þ

Estimation of the quantile eq 4 of θi can be achieved using the same set of points by only using the values of this degradation rate in the sampled parameter vector. These values are used to estimate the marginal distribution of θi which is used to estimate the quantile eq 4. Both estimates are derived from Θ and can therefore only be reliably estimated within confidence interval. This is especially relevant for the quantile estimate of the degradation rate as the 4431


Environmental Science & Technology confidence interval is the range of values being compatible with the data. Situations where the quantile estimate of the degradation rate is nonzero but the confidence interval contains the zero indicate nonsufficient proof for a nonzero rate. These confidence intervals can be calculated using bootstrapping procedures (cf. chapter 7 in ref 3). Markov-Chain Monte Carlo Method. To perform the calculations outlined in the previous section and to estimate the quantile of the distribution in eq 4 or to calculate model predictions eq 3 it thus suffices to have an algorithm which generates random samples from the posterior distribution in eq 1. If this distribution were to depend only on one parameter several approaches exist to generate such samples (e.g., using the cumulative distribution). In higher dimensions such approaches do not work. In these situations Markov-Chain Monte Carlo (MCMC) methods are employed. Under certain technical assumptions they are guaranteed to generate samples from any distribution. MCMC can be largely divided into two approaches: Metropolis-Hastings type algorithms9,10 and the Gibbs sampler.11,12 The Metropolis-Hastings algorithm is applicable to general probability densities. It generates a sequence of points which are distributed according to a pre-specified probability density. Its ~) which is used to main ingredient is a proposal density Q(θ,θ generate new points given the current vector. The algorithm prescribes when to accept and when to reject this new point. In detail it proceeds as follows: ~). (1) propose a new point θ ~ by random sampling from Q(θ,θ (2) generate R∈[0;1] uniformly distributed. If Re

~ ~ pðθjDÞ 3 Q ðθ, θÞ ~ θÞ pðθjDÞ 3 Q ðθ,

then θtþ1 = θ ~ else θtþ1 = θ goto (1). By construction the (tþ1)-th point depends on the t-th point. Estimating the quantiles of the marginal for degradation rates requires detailed knowledge of this distribution close to zero. Therefore, we use an asymmetric, mixed proposal density which samples the region of zero degradation rates more closely. This leads to increased confidence at small degradation rate values with small sample sizes. Let Ξ denote the set of indices of model parameters representing degradation rates, then the used asymmetric proposal density is given as ( ~Þ ¼ Q ðθi , θ i

β 3 δ0 þ ð1 βÞ 3 N ðθi , σ2 Þ i ∈ Ξ iˇΞ N ðθi , σ 2 Þ

with δ0 denoting the Dirac distribution localized at zero and N (θi,σ2) representing the Gaussian distribution with mean θi and standard deviation σ. But it has no influence on the asymptotic behavior of the algorithm. Often successive parameter vectors in the chain exhibit significant autocorrelation across more than one point. This phenomenon can be relieved using an adaption of the classical Metropolis-Hastings algorithm called “Delayed Rejection”.13 This is relevant, since in order to perform quantile estimation of the marginal distribution and the prediction given by eq 3 statistically independent or at least uncorrelated samples lead to reduced estimation error. Subsampling of the chain such that consecutive samples do not show any statistically significant autocorrelation ensures this. Decreased autocorrelation length thus requires fewer points to be generated by the Metropolis-Hastings algorithm.

ARTICLE

We use a 1-stage delayed rejection Metropolis-Hastings algorithm with equal proposal density at both stages. The chain generated by an MCMC procedure can be roughly divided into two phases. During the first phase called “burn-in” the points are not yet sampled from the correct distribution and should therefore not be used for estimating any properties of the model. After this phase the chain is said to be converged. Many statistical approaches exist to judge whether the stable phase is already reached1416 but usually formal convergence criteria accept convergence too quickly. Therefore, for all presented results the burn-in period was judged by visual inspection of the traces of all parameters.

’ RESULTS The methodology described in the previous section will be applied to three data sets: data where optimization-based approaches counterintuitively indicate zero degradation rate for the metabolite, an artificial data set generated to demonstrate the ability to identify situations with a degradation rate of zero, and one measurement data set of a six species model. The metabolization processes are always modeled by first-order kinetics. The parameters to be estimated are degradation rates for all substances, all formation fractions, and the initial concentration of the parent substance. The degradation rate for a substance C will be denoted by Kc. All shown results were generated using the proposal density (eq 6) with β = 0.05 and σ = 0.1. These values were chosen based on prior experiments. It has to be noted that the results are uninfluenced by the choice of parameters, they only change the speed of convergence of the Markov-Chain. As the initial concentration of the parent substance takes values on a significantly different scale than the degradation rates and the formation fractions, all values were scaled to [0;1] prior to proposing a new point and rescaled afterward. This guarantees that all parameters are sampled on a comparably detailed grid and no bias is introduced. The 5%-quantile estimates are shown for all degradation rates included in the model. To judge the reliability of this estimate the 90% confidence interval (CI) of the quantile estimate are indicated via the 5% and the 95% value. If one of the estimated quantiles is zero the respective value is indicated by a red line and a green line if the estimate is larger than zero. Impact of Local Approximation. Globally the probability distribution over parameters might be heavily skewed and consequently the approximation and subsequent decisions can be highly misleading even as far as falsely not rejecting the hypothesis that a degradation rate is zero. The latter behavior is regularly perceived in the analysis of transformation experiments. The approximation error induced by symmetric approximations for the probability distributions is demonstrated by the example shown in Figure 2. Here, simulated data of a transformation model with two species are shown. The right figure displays the symmetrical distribution implicitly used for estimation of standard errors by standard nonlinear least-squares methods in comparison to the “real’’ (posterior) distribution. The distribution is heavily skewed toward larger values of the degradation rate. The symmetric approximation overestimates the probabilities of the values close to or equal to zero. For assessing if the degradation rate is significantly non-zero the symmetric distribution is further transformed into a Student’s t-distribution as the standard deviation has to be estimated from the data. This is additionally increasing the problems associated with the symmetrical approximation as this distribution shows larger tail areas 4432



ARTICLE

Figure 2. Left: Concentration measurements for parent substance and solution found by optimization based approach. Middle: Measured data and optimal model approximation for metabolite data. Right: Analytical error distribution and the symmetric Gaussian approximation for this model.

Figure 3. Results of the MCMC analysis. Left: The parent data are well estimated with tight confidence intervals. Right: The metabolite data are well approximated, too.

than the displayed Gaussian. The classical analysis returns a p-value of 0.13 for the degradation rate of the metabolite thus requiring additional experiments although the data clearly indicate a positive degradation rate for the metabolite. This problem can be overcome if the test used to judge whether the data sufficiently indicate a degradation rate larger than zero is based on the correct distribution. This can be achieved by using the Bayesian model formulation combined with the MCMC method which guarantees to generate samples from the posterior distribution and estimating the quantiles of the marginal eq 4. Consequently, the final decision is based on the correct distribution. The results of the MCMC approach will be discussed in detail in the following section. Two-Species Model with Nonzero Degradation Rate. The first data set is already shown in the left and middle images in Figure 2. It is modeled by the previously mentioned two compartment model. Visual inspection of the measurements indicates nonzero degradation rates for the parent substance and the metabolite. The MCMC method was run to generate 500 000 points with the first 50 000 being classified as burn-in period. The remaining points were subsampled such that consecutive points did not show any significant correlation. The results are shown in Figure 3. The MCMC prediction curves (solid red line) were estimated by randomly sampling 500 parameter vectors and averaging the model predictions. Furthermore,

the 5% and 95% quantiles of the prediction accuracy are reported in the plots (dashed black line). The prediction curves nicely represent the parent and metabolite data indicating that the chosen model agrees well with the measurements. In contrast to the classical approach the analysis of the marginals for both substances indicates that the data contain sufficient proof for the degradation rates to be positive (cf. Figure 4). The estimated confidence intervals of the quantile for both substances are larger than zero. This result agrees with visual inspection of the data. Two-Species Model with Zero Degradation Rate. To demonstrate the ability of the proposed MCMC approach to correctly identify situations where one species does not have a nonzero degradation rate artificial data were generated. The twocompartment model shown in Figure 1 was used. The model parameters were set to the following: • degradation rate for the parent substance 5.8 103 L/d • degradation rate for the metabolite 0 L/d • initial concentration of parent substance 194 g/L • formation fraction 22.4% Setting these parameters the model was simulated for 365 days, and the predicted time vs concentration profile for parent and metabolite were sampled at 13 time points. Measurement noise was simulated by adding scale-dependent Gaussian noise to the data. The MCMC method was applied to the two substance model with these data. It again generated 500 000 points which were 4433



ARTICLE

Figure 4. Analysis of the marginal distributions shows that the data indicate that both rates (left: degradation rate for parent; right: degradation rate for metabolite) are positive.

Figure 5. Top: MCMC predictions for parent substance (left) and metabolite (right). Bottom: Marginal distributions for the degradation rate of the parent substance (kp, left image) and for the degradation rate of metabolite (km, right image). Vertical lines indicate 5%, 50%, and 95% values of the confidence region for the 5% quantile of the respective parameter. The model correctly identifies that km is zero and that kp is significantly larger than zero.

subsampled to show no significant autocorrelation. The MCMC prediction curves and the marginals are shown in Figure 5

Again, the Bayesian prediction for parent and metabolite nicely coincide with the simulated data. The confidence regions 4434


Environmental Science & Technology for the model prediction for parent and metabolite substance seem plausible. The measurements for the parent substance show a significant decline. Analysis of the marginal of the degradation rate of the parent substance clearly identifies a nonzero degradation rate. Analyzing the marginal of the metabolite indicates that the data do not suggest a nonzero degradation rate. Even the 95%-limit of the confidence interval is exactly zero reflecting that the model is highly confident that this quantile is zero. The method thus correctly identifies situations where the degradation rate is zero. Six-Species Model. The previous two examples relied on the simplest transformation model consisting of one drug and one metabolite. Although problems are already encountered in this type of models the distortion problem is more severe if the number of substances is increased and the metabolization pattern

Figure 6. Structure of the complex model.

ARTICLE

is more complex. Therefore, the method is applied to a very complex example. It consists of one parent substance plus five metabolites with complex transformation pathways. A schematic representation of this substance is show in Figure 6. The soil concentrations for all 6 substances were measured at 12 time points over the course of 120 d with 2 measurements per time point. These measurements were like measurements from different time points (i.e., uncorrelated measurement noise was assumed). During the generation of new parameter vectors in the Metropolis-Hastings algorithm it was explicitly taken care of that the formation fractions of any substance always sum to one. The MCMC chain was run to provide 750 000 points (after burnin). This is more than in the previous example as more parameters have to be estimated for this model. Again, these were subsampled to give an uncorrelated sample. The achieved predictions are shown in Figure 7. Analysis of the measurement data suggests that only metabolite B3 does not clearly support a nonzero degradation rate. The degradation rate for metabolite B3 seems to be small, and therefore the data will have to contain highly reliable support for a nonzero degradation rate of B3. All other substances show an obvious decrease in concentration. The Bayesian predictions were calculated according to eq 3, i.e., by averaging the model predictions belonging to different parametrizations and can therefore be different from model predictions based on only one, optimized parameter vector. The resulting marginal distributions for all six substances as well as the estimated 5% quantiles are shown in Figure 8. Analysis of the marginal distributions indicates that the degradation rates for all substances but metabolite B3 are nonzero. This nicely agrees with the visual perception of the measurement data. For the metabolite B3 the method indicates that there is not sufficient evidence in the measurement data that this metabolite has a nonzero degradation rate. The lower end of the confidence interval of the estimated quantile is zero. All results are in good agreement with the data and human perception.

Figure 7. Timeconcentration curves for the complex model Figure 6 with 90% confidence intervals. The data measured for metabolite B3 seem to have a small degradation rate. 4435



ARTICLE

Figure 8. Marginals for all degradation rates in the complex example shown in Figure 6. Rate for metabolite B3 is identified to be borderline.

In all three examples the decisions based on the MCMC analysis agree with visual analysis of the data. This indicates that founding the decision whether a degradation rate is equal to zero on the true posterior distribution of the respective parameter leads to decisions which better reflect the situation represented in the data than the classical LS approach which relies on assumed distributions.

’ DISCUSSION Using different typical examples of data on chemical transformation in soil it was shown that Markov-Chain Monte Carlo methods are able to infer the probability distributions of parameters of kinetic models used to describe those data. In particular the distributions of degradation rates of the applied chemicals themselves as also of their metabolites can be determined, and thus the uncertainty of their estimated values, used for environmental risk assessments, can be estimated. An important application in risk assessments is to judge whether the observed data sufficiently support nonzero degradation rates, as it is for example requested by the guidance given by the FOCUS working group on kinetic evaluation. For this purpose usually a t test is performed based on standard errors of the parameters derived from nonlinear least-squares regression. We have shown that performing this approach, borrowed from asymptotic theory of maximum likelihood estimation, can lead to incorrect results as it implicitly assumes that the locally valid symmetric approximation to the probability distribution of parameters is also globally valid. This assumption is obviously not always correct. As has been shown the distributions can be significantly asymmetric leading to erroneous results for assessments if the symmetric distributions are used. Therefore, we have proposed the use of MCMC methods which do neither perform any approximation to the probability distribution on the parameters nor base their (global) decision on any local asymptotic property of the distribution. They furthermore do not have

to deal with problems of local convergence of the nonlinear optimization algorithm. A special setup of the MCMC method and a problem-specific proposal density to generate samples from the real posterior distribution was proposed. It was shown how to estimate quantiles and model predictions using these samples. The capability of this approach was demonstrated using artificial and real data sets. Concerning the assessment of probabilities of degradation rates being equal to zero it could be demonstrated that the method correctly identifies situations with zero and nonzero degradation rates by applying it to artificial data. On real data it performs in close agreement with human perception of the data sets. Disadvantages of MCMC methods are that they are more time-consuming than regression methods and that in some cases it is hard to estimate whether the chain has already reached its stationary phase or whether it is still in “burning-in”. Although these problems might occasionally induce problems, the MCMC approach is capable of eliminating statistical shortcomings of the classically used method to characterize non-zero degradation rates and yield more reliable and plausible characterizations in such cases.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ REFERENCES (1) FOCUS. Guidance Document on Estimating Persistence and Degradation Kinetics from Environmental Fate Studies on Pesticides in EU Registration; Report of the FOCUSWork Group on Degradation Kinetics, EC Document Reference Sanco/10058/2005 version 2.0, 2006; 434 pp. (2) Nocedal, J.; Wright, S. Numerical Optimization; Springer: New York, 2000. (3) Hastie, T.; Tibshirani, R.; Friedman, J. The Elements of Statistical Learning, 1st ed.; Springer: New York, 2001. 4436



ARTICLE

(4) OECD. Aerobic and Anaerobic Transformation in Soil; OECD Guideline for Testing of Chemicals 307. 2002; adopted 24 April 2002. (5) Shao, J. Mathematical Statistics; Springer: New York, 1998. (6) Bolstad, W. Introduction to Bayesian Statistics, 2nd ed.; Wiley Interscience: 2007. (7) Gelman, A.; Carlin, J.; Stern, H.; Rubin, D. Bayesian Data Analysis, 2nd ed.; Chapman & Hall, CRC: 2003. (8) Gao, Z.; Schmitt, W. Improving Uncertainty Analysis in Kinetic Evaluations Using Iteratively Reweighted Least Squares. In Preparation. (9) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087–1092. (10) Hastings, W. Monte Carlo Sampling Methods using Markov Chains and their Applications. Biometrika 1970, 57, 97–109. (11) Geman, S.; Geman, D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Trans. Pattern Anal. Machine Intell. 1984, 6, 721–741. (12) Gelfand, A.; Smith, A. Sampling Based Approaches to Calculating Marginal Densities. J. Am. Stat. Assoc. 1990, 85, 298–409. (13) Mira, A. On Metropolis-Hastings algorithms with delayed rejection. Metron - In. J. Stat. 2001, 59, 231–241. (14) Gelman, A.; Rubin, D. Inference from Iterative Simulation Using Multiple Sequences (with discussion). Stat. Sci. 1992, 7, 457–511. (15) Raftery, A.; Lewis, S. In Bayesian Statistics 4; Bernardo, J., Berger, J., Dawid, A., Smith, A., Eds.; Oxford University Press: 1992; pp 763773. (16) Heidelberger, P.; Welch, P. Simulation Run Length Control in the Presence of an Initial Transient. Oper. Res. 1983, 31, 1109–1144.

4437


Statistical Analysis of Chemical Transformation Kinetics Using

Recommend Documents