Environ. Sci. Technol. 2010, 44, 1720–1727
Particle and Microorganism Enumeration Data: Enabling Quantitative Rigor and Judicious Interpretation M O N I C A B . E M E L K O , * ,† PHILIP J. SCHMIDT,† AND PARK M. REILLY‡ Department of Civil and Environmental Engineering and Department of Chemical Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
Received August 5, 2009. Revised manuscript received December 21, 2009. Accepted January 6, 2010.
Many of the methods routinely used to quantify microscopic discrete particles and microorganisms are based on enumeration, yet these methods are often known to yield highly variable results. This variability arises from sampling error and variations in analytical recovery (i.e., losses during sample processing and errors in counting), and leads to considerable uncertainty in particle concentration or log10-reduction estimates. Conventional statistical analysis techniques based on the t-distribution are often inappropriate, however, because the data must be corrected for mean analytical recovery and may not be normally distributed with equal variance. Furthermore, these statistical approaches do not include subjective knowledge about the stochastic processes involved in enumeration. Here we develop two probabilistic models to account for the random errors in enumeration data, with emphasis on sampling error assumptions, nonconstant analytical recovery, and discussion of counting errors. These models are implemented using Bayes’ theorem to yield posterior distributions (by numerical integration or Gibbs sampling) that completely quantify the uncertainty in particle concentration or log10-reduction given the experimental data and parameters that describe variability in analytical recovery. The presented approach can easily be implemented to correctly and rigorously analyze single or replicate (bio)particle enumeration data.
Introduction There are many applications in which the rigorous quantification of discrete, uniform particles is of interest; applications concerning the abundance of microorganisms (particularly pathogens and associated indicators) or colloidal particles used as surrogates for microorganisms are numerous. Common environmental examples of such applications include (bio)colloid transport in groundwater (1) and the monitoring and treatment of pathogens in drinking and recreational water (2, 3). Other examples include detection and quantification of microbial contamination in food and food-processing equipment (4), detection and disinfection of biological pathogens such as Bacillus anthracis (5), and * Corresponding author phone: (519) 888-4567 x. 32208; fax: (519) 888-4349; e-mail:
[email protected]. † Civil & Environmental Engineering. ‡ Chemical Engineering. 1720
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 44, NO. 5, 2010
numerous health applications such as the enumeration of parasites in fecal samples (6) or cryopreserved cells in autologous transplantation (7). In many of these applications, quantitative data are obtained by enumeration; the abundance of a specific type of particle or microorganism is estimated by counting the number present in a sample of specified volume. Enumeration data, however, are often highly variable due to sampling and analytical errors, and this variability contributes to uncertain concentration (or log10-reduction) estimates. Moreover, the number of particles per unit volume may be a biased estimate of concentration if a fraction of the particles in the sample are not observed (i.e., recovered) due to various types of losses during sample processing and counting (8). For this reason, analytical recovery (the capacity of the analyst to successfully enumerate each particle in a sample using a specific enumeration method) is an important characteristic of enumeration methods (9). Enumeration-based methods include microscopic assays, flow cytometry, and plating; quantitative PCR and multiple tube fermentation methods, however, are also susceptible to many of the same types of errors. Measurement of the concentration of Cryptosporidium oocysts in water is an example in which enumeration methods are routinely used and analytical recovery is low and commonly quantified. Enumeration of oocysts in surface water used as a source for large drinking water systems was required in the United States (Long Term 2 Enhanced Surface Water Treatment Rule) to determine the level of treatment that is required to adequately reduce risk to consumers (2). Enumeration is also common in evaluations of environmental occurrence (10-12) and removal and/or disinfection during treatment (13, 14). Cryptosporidium oocyst (and Giardia cyst) data are typically obtained using USEPA Method 1623 (15), and accredited laboratories must meet or exceed specified quality control criteria (16). The method, which includes concentration, sample purification by immunomagnetic separation, and immunofluorescence assay, is prone to substantial losses. Laboratory and method validation includes a recovery experiment in which replicate samples containing approximately known numbers of oocysts are enumerated to determine the mean recovery and relative standard deviation (24-100% and c · 10-F)dc
(14)
Results and Discussion Model Implementation. Through the use of Bayes’ theorem, it is possible to obtain a posterior distribution that describes what the source concentration is likely to be given the experimental data and parameters for nonconstant analytical recovery. Integration (by numerical approximation) and Gibbs sampling have each been presented above as methods VOL. 44, NO. 5, 2010 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1723
FIGURE 1. Gibbs sampling algorithms: (a) beta-Poisson model, (b) negative binomial model. to obtain the posterior distribution. Each of these methods has been successfully implemented using the Visual Basic Editor in Microsoft Excel or other similar programming languages, but implementation is not a trivial task nonetheless. Brief discussion of the algorithms follows, and more specific details are provided in the Supporting Information. In the integration method, with improper uniform prior, the posterior distribution for concentration is proportional to the likelihood function. An algorithm is required to numerically approximate the integral of the likelihood function for all values c > 0. This is relatively easily implemented by estimating the sum of the area under the likelihood function for a large number of narrow segments (i.e., the summation of the likelihood at the midpoint of each segment multiplied by segment width). It is sometimes challenging to maintain efficiency while choosing a sufficiently high upper concentration limit and sufficiently narrow segment width to ensure an acceptably accurate approximation of the integral. The posterior probability density for any value of concentration can be calculated by dividing its likelihood by the integral of the likelihood function, while the cumulative density function must be calculated using numerical integration. The results were output as a matrix of concentration values, probability densities, and cumulative densities. Summary statistics such as the mean, mode, and variance, and a 95% credible interval (an interval within which the analyst can be 95% confident that the true value of concentration lies given the Bayesian model) can be calculated. Gibbs sampling can be implemented using the algorithms presented in Figure 1. The recommended estimates that are used to start either Markov chain will ensure rapid convergence, so only a small number of burn-in iterations need to be discarded. Convergence and mixing of the Markov chain is discussed in the Supporting Information. Random number generators for gamma, beta, and Poisson distributions can be developed based on published algorithms. For these relatively simple models, large numbers of iterations (e.g., one million) can usually be completed in at most a few minutes. The concentration value generated in each iteration can be stored, and the mean, mode, variance, and 95% credible interval can be calculated. The mode must be estimated from a histogram, which is sensitive to the number of bins and the number of Gibbs sampling iterations (i.e., the estimated mode may vary between runs of the Gibbs sampling algorithm). Consequently, accurate estimation of the mode requires many more iterations than may be required for estimation of the mean or variance. The results were output as a matrix of concentration values, relative frequencies, and cumulative relative frequencies from a histogram with a large number of bins. 1724
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 44, NO. 5, 2010
Either approach to calculate posterior concentration distributions can usually be completed in seconds to a few minutes. Both methods have been observed to yield very comparable summary statistics and distributions (results not shown). Calculation of the posterior log10-reduction distribution using eq 14 has been observed to give poor results if density function values are tabulated for too few concentration values in either posterior concentration distribution. The credible intervals shown herein are the narrowest credible intervals with specified probability. Because the posterior distributions are often skewed, an interval that excludes 2.5% of the distribution in either tail will be wider than some other 95% credible intervals. For a unimodal distribution, the narrowest credible interval of specified probability will occur where the probability density is equal at either limit. This is relatively easily calculated using the numerical integration algorithm, but is more difficult using Gibbs sampling because the probability density for a specific concentration value can not be precisely calculated with small numbers of iterations. Statistical Analysis of Enumeration Data Using Probabilistic Models. The posterior distributions for concentration and log10-reduction, which can be obtained by integration or Gibbs sampling, provide more complete quantitative information than a point estimate or confidence interval. The distributions can be used (1) to calculate summary statistics such as the mean, mode, and variance, (2) to determine credible intervals (a convenient summary of uncertainty in a point estimate), (3) to evaluate the statistical evidence against specific hypotheses, and (4) as complete knowledge of uncertainty in more complex models (e.g., quantitative microbial risk assessment). While not strictly consistent with a hypothesis test, the weight of evidence against a null hypothesis can be obtained from the cumulative posterior density functions. For example, the weight of evidence against an alternative hypothesis that the concentration is less than 10 microorganisms/L would be the posterior probability that the concentration is greater than 10 microorganisms/L; if this value (which is analogous to a P-value) is small, then the concentration is significantly less than 10 microorganisms/L. Similarly, evidence that the difference between two concentrations is more than 10 microorganisms/L would be a small posterior probability that the log10-reduction is less than 1.0. A set of example posterior distributions follows, which are based on a simulated experiment to quantify log10reduction. The enumeration data shown in Table S2 of the Supporting Information were simulated using the betaPoisson model with initial and final concentrations of 50 and 0.5 microorganisms/L, respectively, and recovery parameters a ) 287.08 and b ) 94.76 (µ ) 75.18%, σ ) 2.208%). The posterior distributions (based on numerical integration
FIGURE 2. Sample posterior distributions for concentration and log10-reduction obtained by numerical integration of the beta-Poisson model using enumeration data from Table 2. (a) Posterior distribution for initial concentration and 95% credible interval (CI). [mean ) 50.91, 95% CI ) (46.78, 55.11)]. (b) Posterior distribution for final concentration and 95% CI. [mean ) 0.539, 95% CI ) (0.423, 0.659)]. (c) Posterior distribution for log10-reduction and 95% CI. [mean ) 1.978, 95% CI ) (1.877, 2.081)]. of the beta-Poisson model) for initial and final concentration as well as log10-reduction are shown in Figure 2. As expected, the 95% credible intervals contain the true values of the initial and final concentrations and log10-reduction. Credible intervals for experimental data, however, may not contain the true values with specified probability due to departures from the assumptions of the models (e.g., validity of the Poisson distribution for random sampling error, validity of the beta or gamma distributions for nonconstant analytical recovery, accuracy of the recovery parameters). As with any model (including a conventional t test or confidence interval), the results may be wholly incorrect if critical assumptions are violated. Comparison to Conventional Confidence Intervals and Hypothesis Tests. The results obtained using the statistical approaches presented in this paper (e.g., estimation of credible intervals) cannot be compared with those obtained using conventional techniques (e.g., confidence intervals and hypothesis tests for a single mean that use the t-distribution) because the assumptions are very different. When constructing confidence intervals that use Student’s t-distribu-
tion, it is necessary to assume that the data are unbiased, independent replicates that are normally distributed with mean µ and variance σ2. The resulting (1 - R) × 100% confidence limits, xj ( tR/2 · s/n, are a function of the data j sample standard deviation s, and only (sample mean x, number of replicates n) and make no further assumptions about the distribution that the data follow. In contrast, credible intervals assume that the data follow a distribution defined by a probabilistic model (e.g., eqs 1 and 2) as well as a prior distribution for the parameter. If the probabilistic model and recovery parameters are assumed to be true, then the resulting credible intervals more correctly quantify uncertainty in the source concentration than the corresponding confidence intervals. The first step in conventional analysis is to ensure that the data are unbiased estimates of source concentration. The enumeration data must be converted to concentration estimates and corrected for mean analytical recovery using the equation cˆ ) x/Vµp, where x is the observed count, V is the sample volume, and µp is the mean analytical recovery. The mean analytical recovery of the enumeration method is unknown, but it can be estimated from the results of a recovery experiment (9). The data must also have uniform variance (i.e., measurement error arising from sample collection and processing), which can be assumed if all samples have equal volume and analytical recovery distributions. Determining the suitability of the normal assumption for a small set of data, or evaluating the possible error due to non-normality in computed confidence intervals and P-values, is difficult. A simple indicator of non-normality is the P-value (calculated from the t-distribution) for the alternative hypothesis that the mean of the concentration estimates is greater than zero; the normal approximation is obviously not valid if the probability of a negative concentration is not very small. It is conceptually possible to transform data so that they meet the required assumptions of normality and constant variance, but this is often impractical with the limited amounts of replicate data that are available in typical practice. The simulated data in Table S2 are used below to illustrate the differences between conventional and Bayesian statistical analysis. In these simulations, the recovery varied according to a beta distribution with a ) 287.08 and b ) 94.76, so each concentration estimate must be adjusted using µp ) 75.18%. Figure 3a shows that the count distribution from which the final concentration data were generated is approximately normally distributed; the same is true for the initial concentration data (results not shown). Accordingly, the approximate 95% confidence intervals are (40.67, 60.95) and (0.271, 0.793) for the initial and final concentration respectively. In contrast, the associated 95% credible intervals shown in Figure 2 are (46.78, 55.11) and (0.423, 0.659); however, these intervals cannot be directly compared to confidence intervals because the assumptions and interpretation are different. The stated level of confidence (or P-value) associated with conventional analysis of enumeration data will only be an estimate if the data are only approximately normally distributed. Non-normal count distributions are prevalent when the numbers of observations (i.e., counts) are typically low (e.g., Figure 3b). A major limitation of conventional statistics such as confidence intervals is that the results depend only on the data and do not incorporate knowledge of the method by which the data were collected (e.g., knowledge that random sampling error follows a Poisson distribution or about the variability in analytical recovery observed in a recovery experiment). In contrast, this information is included in the posterior distributions presented herein, so these represent more complete knowledge about uncertainty in the source concentration. Furthermore, there are a number of practical VOL. 44, NO. 5, 2010 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1725
herein. This information is available free of charge via the Internet at http://pubs.acs.org/.
Literature Cited
FIGURE 3. Normal approximation of counts obtained from the beta-Poisson model. (a) (a ) 287.08, b ) 94.76, c ) 0.5, V ) 50). (b) (a ) 287.08, b ) 94.76, c ) 0.5, V ) 10). situations where the Bayesian methods presented herein can easily be implemented, but direct calculation of confidence intervals or hypothesis testing based on the t-distribution is not possible. These include quantification of uncertainty (1) when the data are not approximately normally distributed, (2) when samples have unequal variance (e.g., different sample volumes), (3) when the sample variance is zero, and (4) when there is only one datum. Posterior distributions are also useful in risk analysis applications in which the uncertainty in important parameters needs to be quantified (e.g., in demonstrating effective treatment of potable water). Further research on the models presented herein is warranted to determine the sensitivity of the models to various assumptions (e.g., nonrandom distribution of particles in the source, accuracy of recovery parameters) and to make recommendations for improved experimental design (e.g., the relative value of replication and increased sample size).
Acknowledgments We thank Dr. Mary Thompson of the Department of Statistics and Actuarial Science at the University of Waterloo for her assistance and the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canadian Water Network for financial support. The comments of three anonymous reviewers contributed to the significant improvement of this manuscript.
Supporting Information Available Derivation of the negative binomial model, a table summarizing the model components, discussion of how to incorporate sample-specific recovery parameters in the presented analyses, discussion of numerical integration, example Gibbs sampling code, discussion of convergence and mixing of the Gibbs sampling algorithm, and the enumeration data used in the example analyses presented 1726
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 44, NO. 5, 2010
(1) Wang, M.; Ford, R. M.; Harvey, R. W. Coupled effect of chemotaxis and growth on microbial distributions in organicamended aquifer sediments: Observations from laboratory and field studies. Environ. Sci. Technol. 2008, 42, 3556–3562. (2) USEPA. National Primary Drinking Water Regulations: Long Term 2 Enhanced Surface Water Treatment Rule; Final Rule. Fed. Regist. 2006, 71, 653-786. (3) USEPA. Water Quality Standards for Coastal and Great Lakes Recreation Waters; Final Rule. Fed. Regist. 2004, 69, 6721767243. (4) APHA. Compendium of Methods for the Microbiological Examination of Foods, 4th ed.; APHA: Washington, DC, 2001. (5) Edmonds, J. M.; Collett, P. J.; Valdes, E. R.; Skowronski, E. W.; Pellar, G. J.; Emanuel, P. A. Surface sampling of spores in drydeposition aerosols. Appl. Environ. Microbiol. 2009, 75, 39–44. (6) Teixeira, C. F.; Neuhauss, E.; Ben, R.; Romanzini, J.; GraeffTeixeira, C. Detection of Schistosoma mansoni eggs in feces through their interaction with paramagnetic beads in a magnetic field. PLoS Negl. Trop. Dis. 2007, 1, e73. (7) Sartor, M.; Antonenas, V.; Garvin, F.; Webb, M.; Bradstock, K. F. Recovery of viable CD34+ cells from cryopreserved hemopoietec progenitor cell products. Bone Marrow Transpl. 2005, 36, 199– 204. (8) Emelko, M. B.; Schmidt, P. J.; Roberson, J. A. Quantification of uncertainty in microbial data - reporting and regulatory implications. J. Am. Water Works Assoc. 2008, 100, 94–104. (9) Schmidt, P. J.; Emelko, M. B.; Reilly, P. M. Quantification of analytical recovery in particle and microorganism enumeration methods. Environ. Sci. Technol. 2010 (doi:10.1021/es902237f). (10) Rose, J. B.; Gerba, C. P.; Jakubowski, W. Survey of potable water supplies for Cryptosporidium and Giardia. Environ. Sci. Technol. 1991, 25, 1393–1400. (11) Crainiceanu, C. M.; Stedinger, J. R.; Rupert, D.; Behr, C. T. Modeling the U.S. national distribution of waterborne pathogen concentrations with application to Cryptosporidium parvum. Water Resour. Res. 2003, 39, 1235–1249. (12) Petterson, S. R.; Signor, R. S.; Ashbolt, N. J. Incorporating method recovery uncertainties in stochastic estimates of raw water protozoan concentrations for QMRA. J. Water Health 2007, 5, 51–65. (13) Emelko, M. B.; Huck, P. M. Microspheres as surrogates for Cryptosporidium filtration. J. Am. Water Works Assoc. 2004, 96, 94–105. (14) Assavasilavasukul, P.; Lau, B. L. T.; Harrington, G. W.; Hoffman, R. M.; Borchardt, M. A. Effect of pathogen concentrations on removal of Cryptosporidium and Giardia by conventional drinking water treatment. Water Res. 2008, 42, 2678–2690. (15) USEPA. Method 1623: Cryptosporidium and Giardia in Water by Filtration/IMS/FA; EPA 815-R-05-002; U.S. Environmental Protection Agency, Office of Water: Washington, DC, 2005. (16) Clancy, J. L.; Connell, K.; McCuin, R. M. Implementing PBMS improvements to USEPA’s Cryptosporidium and Giardia methods. J. Am. Water Works Assoc. 2003, 95, 80–93. (17) Gronewold, A. D.; Borsuk, M.; Wolpert, R. L.; Reckhow, K. H. An assessment of fecal indicator bacteria-based water quality standards. Environ. Sci. Technol. 2008, 42, 4676–4682. (18) Nahrstedt, A.; Gimbel, R. A statistical method for determining the reliability of the analytical results in the detection of Cryptosporidium and Giardia in water. J. Water Supply Res. Technol. AQUA 1996, 45, 101–111. (19) Fisher, R. A. The negative binomial distribution. Ann. Eugenic. 1941, 11, 182–187. (20) Gosset, W. S. On the error of counting with a haemacytometer. Biometrika 1907, 5, 351–360. (21) Eisenhart, C.; Wilson, P. W. Statistical methods and control in bacteriology. Bacteriol. Rev. 1943, 7, 57–137. (22) Haas, C. N.; Heller, B. Test of the validity of the Poisson assumption for analysis of Most-Probable-Number results. Appl. Environ. Microbiol. 1988, 54, 2996–3002. (23) El-Shaarawi, A. H.; Esterby, S. R.; Dutka, B. J. Bacterial density in water determined by Poisson or negative binomial distributions. Appl. Environ. Microbiol. 1981, 41, 107–116. (24) Christian, R. R.; Pipes, W. O. Frequency distribution of coliforms in water distribution systems. Appl. Environ. Microbiol. 1983, 45, 603–609. (25) Haas, C. N.; Heller, B. Statistics of enumerating total coliforms in water samples by membrane filter procedures. Water Res. 1986, 20, 525–530.
(26) Fisher, R. A.; Thornton, H. G.; MacKenzie, W. A. The accuracy of the plating method of estimating the density of bacterial populations. With particular reference to the use of Thornton’s agar medium with soil samples. Ann. Appl. Biol. 1922, 9, 325–359. (27) Parkhurst, D. F.; Stern, D. A. Determining average concentrations of Cryptosporidium and other pathogens in water. Environ. Sci. Technol. 1998, 32, 3424–3429. (28) Teunis, P. F. M.; Evers, E. G.; Slob, W. Analysis of variable fractions resulting from microbial counts. Quant. Microbiol. 1999, 1, 63–88. (29) USEPA. Occurrence and Exposure Assessment for the Final Long Term 2 Enhanced Surface Water Treatment Rule; EPA 815-R06-002; U.S. Environmental Protection Agency, Office of Water: Washington, DC, 2005.
(30) Margolin, B. H.; Kaplan, N.; Zeiger, E. Statistical analysis of the Ames Salmonella/microsome test. Proc. Natl. Acad. Sci. U.S.A. 1981, 78, 3779–3783. (31) Box, G. E. P.; Tiao, G. C. Bayesian Inference in Statistical Analysis; Addison-Wesley: Reading, MA, 1973. (32) Gilks, W. R.; Richardson, S.; Spiegelhalter, D. J. Markov Chain Monte Carlo in Practice; Chapman & Hall: London, 1996. (33) Lee, P. M. Bayesian Statistics, An Introduction, 2nd ed.; John Wiley: Toronto, ON, 1997.
ES902382A
VOL. 44, NO. 5, 2010 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1727
