Campylobacter QMRA: A Bayesian Estimation of Prevalence and Concentration in Retail Foods Under Clustering and Heavy Censoring

A Bayesian statistical temporal‐prevalence‐concentration model (TPCM) was built to assess the prevalence and concentration of pathogenic campylobacter species in batches of fresh chicken and turkey meat at retail. The data set was collected from Finnish grocery stores in all the seasons of the year. Observations at low concentration levels are often censored due to the limit of determination of the microbiological methods. This model utilized the potential of Bayesian methods to borrow strength from related samples in order to perform under heavy censoring. In this extreme case the majority of the observed batch‐specific concentrations was below the limit of determination. The hierarchical structure was included in the model in order to take into account the within‐batch and between‐batch variability, which may have a significant impact on the sample outcome depending on the sampling plan. Temporal changes in the prevalence of campylobacter were modeled using a Markovian time series. The proposed model is adaptable for other pathogens if the same type of data set is available. The computation of the model was performed using OpenBUGS software.


INTRODUCTION
Quantitative microbial risk assessment (QMRA) models are used to model food pathways and describe changes in the prevalence and concentration of pathogenic microorganisms over the various stages of the food chain. (1) Microbiological data in QMRA models consist of the detection and enumeration of pathogens in samples from food production. At the retail level, such samples are often clustered in food batches or due to some other factors. Assessment of the ratio of within-batch and between-batch variability is of central importance in risk assessment because such variability may have a significant impact on the sample outcome and the statistical inference thereof. For example, samples collected from the same batch may all be positive if the batch is contaminated. Only a small number of samples per batch are often taken, which leaves room for uncertainty in both the between-batch prevalence and the within-batch prevalence of the food batches. In 2012, a hierarchical Bayesian model was published by Commeau et al. (2) in order to represent the withinbatch variability and the between-batch variability in foods (pork breasts and cold smoked salmon) contaminated by Listeria monocytogenes (bacteria). In that model, the probability of having a positive detection result depended on the number of colony forming units (cfu) in the sample, which was modeled using a log-normal or Poisson distribution. Three levels of aggregated data (batch, unit, and test portion) could be used in the model by Commeau et al. A log-normal model is widely applied to model this type of enumeration data due to the usually good fit. However, other probability distributions have also been examined in order to model the concentrations of pathogens. (3) Another feature of microbiological data is that the data set may contain small positive concentrations that are below the limit of determination of the microbiological method. In statistics, such observations are called left-censored. Ignoring these observations or setting them to fixed values may have a significant impact on the concentration estimation. (4) Bayesian methods for low concentration estimation have been considered in some previous microbiological studies. (5)(6)(7) These models have been based on a hierarchical log-normal model that takes into account the low concentrations by treating them as censored data. Bayesian inference is found to be an appropriate method for this type of microbiological modeling because it provides a way to take into account the uncertainty related to core model parameters as well as censored data in a coherent way. (8) Campylobacter is one of the major zoonotic agents causing human cases, as it is in the whole European Union. Poultry meat has been considered one of the most important campylobacter sources worldwide, and in Finland as well as in the whole European Union its consumption has increased rapidly during the last two decades. Furthermore, the seasonal occurrence of campylobacter varies strongly in Finnish poultry flocks, while the concentration has been found relatively low. (9) However, the occurrence and amount of campylobacter at the retail level has been less explored.
This article presents a Bayesian model that was developed to assess both the prevalence and the concentration of pathogenic campylobacter species (C. jejuni and C. coli) in chicken meat and turkey meat on the basis of retail-level data (detection and enumeration data). The model was built to perform with clustered data containing only a small number of samples per batch, and when a large proportion of the observed concentrations is below the limit of determination. In particular, the data set includes several contaminated food batches with only one observation from each of the batches and those observations are all below the limit of determina-tion. The Bayesian hierarchical structure borrows strength across food batches, enabling overall estimation as well as batch-specific inference. The model consists of two submodels separately treating detection and enumeration data. A Markovian time series was added to our model in order to take into account and assess the seasonal variation in the estimation of prevalence. In addition, a hierarchical gamma model was compared to the widely used log-normal model for the concentration estimation. Besides the retail level, this model can be used in other stages of the production chain if the same type of microbiological data set is available. The results of the model will be exploited in further risk assessment to predict foodborne infections.

Sample Data
From November 2012 to December 2013 and during the summer of 2014, samples of packed fresh Finnish poultry meat were collected from grocery stores in the Helsinki area. The data collection was carried out throughout the year in order to take into account the seasonal variation. These samples were grouped according to their retail batch numbers. The number of samples taken from a batch was small compared to the batch size because a single retail batch may consist of thousands of meat units. These meat units originate from one or several animal flocks. Hence, samples taken from the same food batch may be correlated due to their common origin. Although the samples were collected from the Helsinki area, the whole country was effectively represented because these companies cover over 90% of the total production volume in Finland and deliver their products all over the country. The detection and quantitative determination of campylobacters were performed using the accredited modified NMKL 119:2007 and ISO 10272-2:2006 methods at the national reference laboratory of Evira. All the samples were examined for campylobacters in the laboratory before their expiration date in order to mimic the true time of consumption after purchase. A 50-g portion of meat was homogenized in 50 mL of buffered peptone water using a stomacher, and 25 mL of the liquid (corresponding to 25 g of meat) was used for enrichment. A 10-fold dilution series was prepared from the liquid for quantitative determination. The total size of the units varied from 250 g to 500 g. The sensitivity and the specificity of the 3 0 microbial detection method were assumed to be 100%. The structure of the data is presented in Fig. 1. The total number of studied batches per meat type was 226 (chicken) and 185 (turkey), and 1 to 16 samples per batch were taken (see Fig. 2 and Appendix A). The total number of contaminated samples was 108, of which 76 were chicken meat and 32 turkey meat. The contaminated samples of chicken meat originated from 31 batches and turkey meat from 17 batches. The number of campylobacters in the majority of these samples was below the limit of determination of the microbiological method, which was 0.5 cfu/g (see Table I). In turkey meat, only two samples yielded a quantitative result. Sixteen out of 31 contaminated chicken meat batches and 16 out of 17 contaminated turkey meat batches only contained censored observations. The maximum observed concentration was 0.5 cfu/g in turkey meat and 38 cfu/g in chicken meat.

Prevalence Model
In the model, the detection and enumeration data were treated separately. The detection data were used to estimate the percentage of retail foods with campylobacter contamination (pf). The unstructured use of the detection data in the modeling of the prevalence may lead to biased results if the samples are clustered in batches. In addition, only a proportion of all batches and a relatively small number of samples per batch were examined, which leaves room for uncertainty in both the between-batch prevalence and the within-batch prevalence. The number of contaminated samples (X i ) among the samples taken from one food batch (N i ) was modeled as a binomial distribution: where B is the total number of studied batches per meat type. The probability of sampling a contaminated meat unit depends on both the between-batch prevalence (pb) and the within-batch prevalence (pw). A conventional uninformative uniform(0,1) distribution was used as a prior for the parameter pw. A binary latent variable I i takes the value of 1 if at least one of the units in a batch i is contaminated, and 0 otherwise. Variable I i describes the "true contamination status" of a batch, with P(I i = 1|pb) = pb.
According to a previous study (10) and the register data, (11) the incidence of campylobacter in both humans and chickens in Finland is significantly higher in late summer than in other seasons. The number of samples collected in different months varied much (see Appendix A). A Markovian time series was added to the model in order to take into account the seasonal variation in campylobacter prevalence in food batches: where the parameter e m represents the changes between consecutive months. In this model the parameter e m is normally distributed with mean = 0 and precision = τ e = 1/σ e 2 . Conventional uninformative distributions were given as a prior for the parameter τ e (gamma(0.001,0.001)) as well as for the parameter pb 1 (Logit(pb 1 )ßNorm(0,0.001)). The monthly proportion of meat units with campylobacter contamination (pf m ) is a product of the between-batch prevalence (pb m ) and the within-batch prevalence (pw). Possible seasonality in within-batch prevalence was not modeled due to insufficient data. The detailed structure of the prevalence model is presented in Fig. 3.

Concentration Estimation Under
Heavy Censoring

Hierarchical Log-Normal Model
A hierarchical log-normal model was used to model the enumeration data (c j,k ) derived from contaminated samples (k) in contaminated batches (j). The hierarchical structure was used in order to evaluate and take into account the clustering of the samples. The following model was first separately applied to the data sets of chicken meat and turkey meat: where J is the number of detected contaminated batches and t j denotes the number of detected contaminated samples in a batch j. The parameter μ j describes the mean concentration of contaminated units in batch j and the precision parameter τ w denotes the inverse of the within-batch variance (τ w = 1/σ w 2 ). Parameter μ 0 is the mean concentration of contaminated units in all contaminated batches. The precision parameter τ b denotes the inverse of the between-batch variance (τ b = 1/σ b 2 ). A con-ventional uninformative distribution was given as a prior for the parameter μ 0 (Norm(0,0.001)) as well as for the parameter τ w (gamma(0.001,0.001)). Instead of placing a prior distribution for the precision parameter τ b , an uninformative prior distribution (uniform(0,100)) was assigned for the standard deviation parameter σ b . The effect of this choice is considered in Section 3.4. The data set included a very high percentage of observations below the limit of determination. The minimum concentration for positive samples was assumed to be 1 cfu/ws, where the parameter ws denotes the weight of one test portion (25 g). This assumption was based on the fact that a positive test result from the detection test means that at least one bacterium was present. A similar assumption was also made in some previous studies. (6,12) Hence, all observations were included in the full likelihood as a product of probabilities: Fig. 3. A directed acyclic graph of the whole model, which was developed on the basis of retail data. In this graph, the log-normal model is presented for concentration estimation. A temporal structure was not included in concentration estimation.
where the indicator vector δ j,k is 1 if the observation is censored and 0 otherwise. The whole data set of measured concentrations ({c}) thus provides information on the concentration model parameters. Finally, the predicted concentration for contaminated meats (c rep ) was obtained as a posterior predictive distribution, based on the posterior distribution P(μ 0 ,τ b ,τ w |{c}) and the following conditional distributions: In other words: where the parameter μ rep indicates the predicted mean concentration in a random batch. A log-normal distribution is known as a heavy-tailed distribution. Hence, there is a small positive probability of extremely high concentrations, even though the typical concentration is low. However, the assumption of the lower limit is beneficial because it prevents the unrealistically small values for censored data. The censoring limits of the data affect the estimation of model parameters, particularly if the information on the variance given by the data is poor, as in this study.
A better way to estimate can be achieved by borrowing strength from other groups via a hierarchical structure. This can be performed by setting the same variance parameters σ w 2 and σ b 2 for both meat types. However, this assumption should be reevaluated if more data become available. The whole prevalenceconcentration model is graphically presented in

Hierarchical Gamma Model
A hierarchical gamma model was used as an alternative approach to model concentrations in contaminated foods: where J is the total number of contaminated batches per meat type and t j is the total number of contaminated meat units in a batch j. The parameters s j and r j are the batch-specific shape and rate (inverse scale parameter r j = 1/λ j ) parameters. A gamma prior distribution was used for these parameters, s j ßgamma(α s ,β s ) and r j ßgamma(α r ,β r ). Separate parameters were given for each batch in order to allow flexibility in batch-related distributions. Uninformative priors were assigned to the hyperparameters α s ,β s ,α r ,β r ßExp(0.001). The intervalcensored observations (below the limit of determination) were included in the gamma model in a similar manner to the log-normal model (see Section 2.3.1).
The predicted concentration for a contaminated meat unit (c rep ) was derived based on the posterior distribution P(α s ,β s ,α r ,β r |{c}) and the conditional distributions: s rep ßgamma(α s ,β s ), r rep ßgamma(α s ,β s ), c rep ßgamma(s rep ,r rep ).
In other words: The benefit of the gamma distribution is that it can present a wide selection of distribution shapes and has directly the same range as the concentration measurements. The gamma distribution is also capable of presenting high concentrations due to its flexibility. Unlike the normal distribution, where distinct parameters represent the mean and variance of the distribution, the mean and variance of the gamma distribution depend on both parameters (shape and scale). Hence, the variance components in different levels of the hierarchical model cannot be compared in a similar way as in the hierarchical normal model for log concentrations. The hierarchical gamma model could also be used to borrow strength across data sets by adding a level of hierarchy to the model. However, this was not done due to the small number of groups (only two meat types).

MCMC Computation
The computation of the model was performed using OpenBUGS software, (13,14) with 100,000 iterations after the burn-in phase of 10,000 iterations. The convergence of the MCMC (Markov chain Monte Carlo) chains was achieved by 10,000 iterations. The OpenBUGS code is provided in Appendix C.

Prevalence in Retail Meats
The average monthly percentage of batches with contamination was estimated to be 13.9% (95% CI: 9.3-19.6%) in chicken meat batches and 12.3% (95% CI: 6.3-22.0%) in turkey meat batches. The seasonal variation in the between-batch prevalence (pb) is illustrated in Fig. 4. The same type of seasonality can be seen in both meat types.
The within-batch prevalence (pw) was estimated to be 59.5% (95% CI: 49.8-68.8%) in contaminated chicken meat batches and 28.6% (95% CI: 18.3-40.0%) in contaminated turkey meat batches. The prevalence of campylobacter in meat units at retail is a product of the between-batch prevalence and the within-batch prevalence. For comparison, the posterior distribution for prevalence in foods was first obtained from the full model and then by ignoring either the temporal effect or clustering of the samples. The sample prevalence was considerably higher than the estimated population prevalence (see Fig. 5) in meat units. This is because more samples were collected during the summer than in other seasons of the year (see Appendix A). The temporal pattern also decreased the uncertainty in the overall prevalence in meat units, resulting in a narrower credible interval for prevalence in meat units. The same type of effect can be seen in both meat types. The posterior mean was slightly greater and the credible interval slightly narrower compared to the full model if the clustering of the samples was ignored.

Concentration in Retail Meats
The between-batch variability in the concentration of campylobacter was estimated to be greater than the within-batch variability for both meat types (see Table II). This indicates that both variations have an important role in the estimation and prediction. The total variance of the log concentration was

Note:
The figures are presented separately for chicken meat (c) and turkey meat (t). The results obtained by using the same variance parameters for chicken meat and turkey meat batches are given in the two bottom rows (c&t). The lower total variation in turkey meat is explained by the fact that all observations from turkey meat lie in the interval of 0-0.5 cfu/g. In addition, the results were calculated by using the same variance parameters for chicken meat and turkey meat. In this case, the total variance was 0.35 (95% CI: 0.24-0.52). The same variance structure may be necessary to assume if the information given by one data set is very poor.
The posterior predictive distribution for the concentration based on a log-normal model is presented in Table III. This is a conditional distribution for random contamination at retail, given the full data set and the priors. The predicted concentration was first obtained with the hierarchical model and then without the batch variability. In addition, the prediction was obtained by using the same variance parameters for chicken meat and turkey meat. The main difference in the predictive distributions from the different models is the length of the tail (see Table III). However, small differences can also be seen in the posterior means and medians. The nonhierarchical model leads to a predictive distribution in which the tail is shorter compared to the hierarchical model, which was better supported by the data than the Note: The results were calculated with the full model (hierarchical) and then by omitting the structure that describes the clustering of the samples (nonhierarchical). The results were also calculated by using the same variance parameters for chicken meat and turkey meat (hierarchical comb.).
nonhierarchical model (see Section 3.5). In addition, the posterior predictive distribution appeared to be slightly sensitive to the prior choice, if the nonhierarchical model was used for the limited and heavily censored turkey data set. The results (Table III) were obtained using a gamma(0.01,0.01) prior for the variance parameter and a normal(0,0.01) prior for the mean parameter in the nonhierarchical model. The same variance structure led to a wider credible interval for turkey because the data set from turkey meat alone was extremely homogeneous. For the same reason, the upper limit of the credible interval becomes smaller for chicken meat. It may be beneficial to borrow strength over similar data sets if some of them are very small with a very limited range of concentrations.
The gamma model led to a right-skewed predictive distribution of the campylobacter concentration, which has a very long tail compared to the log-normal model for chicken meat. Thus, the posterior mean obtained with the gamma model is also much greater than with the log-normal model (see Table IV). A quite similar predictive distribution of concentration was obtained for turkey meat with the gamma and the log-normal model. The fit of these models is compared in Section 3.5. The posterior predictive distribution obtained with the nonhierarchical gamma model was slightly sensitive to the prior choice similar to the nonhierarchical log-normal model.

Sensitivity Analysis
The choice of the prior distribution for the between-batch variance parameter (σ b 2 ) is usually a key problem in hierarchical models of this type. An inverse-gamma distribution was not chosen as the prior for the parameter σ b 2 = 1/τ b in the hierarchical log-normal model because inference might have become sensitive to the values of parameters set for this prior distribution. (15) As a more robust choice, an uninformative uniform(0,100) prior distribution was assigned for the standard deviation parameter σ b . The sensitivity of this prior choice was studied by using an uninformative half-normal (0,100) and a weakly informative half-normal (0,1) distribution as alternatives for the standard deviation parameter. The posterior distributions for the predicted concentration in both chicken meat and turkey were almost identical with all these priors. This indicates that the results are not sensitive to various reasonable choices of priors. These priors are convenient when data are sufficient for the estimation. If not, more informative priors could be used but these may be difficult to obtain from expert opinions if there is virtually no previous knowledge of the kinds of samples studied.
In this model, a conventional uninformative gamma(0.001,0.001) distribution was used as a prior for the precision parameter τ m = 1/σ m 2 that describes the inverse variance of parameter e m in the Markovian time series (see Section 2.2). The effect of this prior choice was studied by setting different values of parameters to the gamma distribution (gamma(0.01,0.01) and gamma(0.1,0.1)) and also by using a uniform(0,100) prior distribution for the parameter σ m . No significant difference in the posterior distributions of monthly between-batch prevalences was seen using different priors. However, the convergence of the MCMC chains was slower if the uniform distribution was used.

Model Comparison
The deviance information criterion (DIC) is a commonly used metric when comparing different Bayesian models with the same data set. DIC is the sum of the posterior mean of the deviance (Dbar) and the number of effective parameters in the model (pD). It is possible to obtain negative values of DIC, but only the difference in the values of DIC between different models is meaningful. The model with the smallest DIC is best supported by the data. As a rough basic rule, (16) the models within 1-2 units of the best model are almost equally supported by the data as the model with smallest DIC. The models within 3-7 units of the best model are significantly less supported than the best model. If the difference is more than seven units, these models have no support compared to the best model. The different concentration models proposed in this article were compared by using the DIC tool, which is built into OpenBUGS software. DIC values obtained from different models are presented in Table V. The hierarchical log-normal model has a much smaller DIC value compared to the nonhierarchical log-normal model for both data sets (chicken and turkey). The hierarchical log-normal model was also significantly better supported by both data sets than the hierarchical gamma model. The nonhierarchical gamma model has a smaller DIC value than the nonhierarchical log-normal model for the turkey data. The difference is more than two units but less than seven units. The nonhierarchical gamma model was also slightly better supported by the turkey data than the hierarchical gamma model. This may be due to the extremely homogeneous turkey data. However, DIC values obtained with the nonhierarchical models for the turkey data set may be sensitive to the prior choice.

DISCUSSION
The already known clustering of campylobacter contamination in slaughter batches propagates down to the batches at retail. Omitting the clustering of the samples in the estimation would lead to biased estimates and an inaccurate level of uncertainty. The extreme amount of censoring in the enumeration data presented a large proportion of observations below the limit of determination of the microbiological method, which is not an entirely unusual situation in quantitative microbiological studies, but can make statistical inference challenging. Careful statistical analysis of these observations is crucial in producing a predictive posterior distribution for the concentration. For instance, if censored observations had been ignored in this study, the predicted mean concentration for contaminated chicken meat would have been much higher (6.4 cfu/g vs. 1.9 cfu/g) and the 95% credible interval much wider (0.1-35.7 cfu/g vs. 0.0-12.2 cfu/g) compared to the results obtained with full data. In turkey meat, nearly all observations were below the limit of determination. Hence, the whole analysis would have been challenging without censored data. The proposed model makes it possible to take into account all observations, even when the majority of the contaminated batches only included concentrations that were below the limit of determination. The assumption behind this method is that all positive concentration observations follow the same probability distribution, i.e., provide information about the same model parameters. The Bayesian methods provide a consistent framework for this type of modeling. This is because the uncertainty related to core model parameters can be handled at once and exploit several data sets, for example, by borrowing strength across meat batches or even meat types when the information given by one data set is insufficient. Inference can be performed even with a very small number of observations per batch. The proposed hierarchical gamma model offers a competing method to the hierarchical log-normal model to predict campylobacter concentrations. However, the results indicate that the hierarchical log-normal model was better supported by the data than the hierarchical gamma model.
Previous studies have not considered the seasonal variation in the prevalence estimation. It was beneficial to include a time series model for the prevalence estimation because temporal variation is known to exist. In the future, the model could be expanded by also including the temporal variation in the concentration estimation.
Bayesian methods are also well suited to making predictions that can easily be produced once the posterior distribution of parameters has been obtained. The resulting predictive distribution presents the total epistemic and aleatoric uncertainty in the unknown quantities based on all available information. The results of this model can be exploited in the wider context of campylobacter QMRA. The posterior predictive distribution for the pathogen concentration can be used as an input for the doseresponse model in making predictions of the risk of foodborne campylobacter infections. In addition, the resulting posterior distribution for the percentage of meat units with campylobacter contamination can be exploited when the number of human campylobacteriosis cases due to the consumption of poultry meat is estimated at the population level.
The proposed model is easily applicable to other microbiological data sets. Other types of censored data can be included in a similar manner. The model does not require strongly subjective expert opinions, but aims to be data-driven and applicable under clustering and heavy censoring. The conventional uninformative prior distributions were sufficient for the purpose and used for all unknown model parameters. The minimum value for the number of observations that is required to use uninformative prior distributions may depend not only on the number of samples but also on the number of batches they belong to as well as the proportion of samples that is censored. An optimal sampling plan is an issue that needs to be further studied.