Reliability of environmental sampling culture results using the negative binomial intraclass correlation coefficient

The Intraclass Correlation Coefficient (ICC) is commonly used to estimate the similarity between quantitative measures obtained from different sources. Overdispersed data is traditionally transformed so that linear mixed model (LMM) based ICC can be estimated. A common transformation used is the natural logarithm. The reliability of environmental sampling of fecal slurry on freestall pens has been estimated for Mycobacterium avium subsp. paratuberculosis using the natural logarithm transformed culture results. Recently, the negative binomial ICC was defined based on a generalized linear mixed model for negative binomial distributed data. The current study reports on the negative binomial ICC estimate which includes fixed effects using culture results of environmental samples. Simulations using a wide variety of inputs and negative binomial distribution parameters (r; p) showed better performance of the new negative binomial ICC compared to the ICC based on LMM even when negative binomial data was logarithm, and square root transformed. A second comparison that targeted a wider range of ICC values showed that the mean of estimated ICC closely approximated the true ICC.


Introduction
In the simple case of estimating the correlation among 2 factors with a set of quantitative observations, an investigator may elect to utilize the Spearman Rank correlation coefficient or Pearson's correlation coefficient assuming the observations are independent. The measure of agreement κ can be estimated for correlation between binary observations. For more complex data structures that may include either crossed or nested factors of a latent character, the investigator may utilize the Intraclass Correlation Coefficient (ICC). The ICC is related to unexplained variance at the subject level. More specifically, the ICC is defined as the ratio of the covariance of measurements from the factor of interest to the marginal variance of the observations. Ranging between 0 and 1, an ICC close to 1 indicates that the difference in observations due to the factor of interest are ignorable. Hence, using variance estimates attributable to each of a study's factors, the ICC can be used as a measure of similarity in observations between subjects due to a particular factor. A direct application of the ICC is a measure of the correlation between subjects in a reliability and repeatability gauge study (Aly et al. 2009;Kittawornrat et al. 2012).
Investigators analyze and obtain variance estimates for normally distributed data using linear mixed models (LMM) or non-normally distributed data using generalized linear mixed models (GLMM). Health science researchers more commonly work with count data and while the ICC for the LMM has been extended to the Poisson case (Carrasco and Jover 2005), its equivalence for count data with overdispersion was only recently described (Carrasco 2010). Until the ICC for negative binomial distributed data was developed, researchers transformed such data using different transformations to make their data normally distributed in order to use LMM and their ICC.
An example of count data that may commonly be overdispersed is bacterial culture results. Culture results are commonly reported as colony forming units per specimen mass or culture medium tube. Another example is parasite counts which are commonly reported as parasitic stage count per gram of specimen. Given the nature of such infectious agents, they can exist in very large numbers within their hosts, at the same time not all potential hosts in a population are infected. In fact, more hosts tend to be uninfected leading to the inequality of the mean and variance of the data, hence overdispersion. In the current study, we report on a reliability analysis for environmental sampling to quantify Mycobacterium avium subspecies paratuberculosis (MAP) on California free-stall dairies (Aly et al. 2009). A previous study with these data was unique in that it involved the use of nested and crossed factors and used the natural logarithm to attain normally distributed data for a LMM analysis and ICC estimation. Such transformations may normalize the data provided the number of replicates was large and the variance components were small (Solomon and Taylor 1999). Both sample size and magnitude of variance conditions may be difficult to attain with negative binomial distributed data especially when replicates are limited due to cost or subject use limitations such as in health sciences research. The performance of the negative binomial ICC has not been compared to LMM ICC using previously described data transformations in multilevel models with crossed and nested random effects.
Hence, the objectives of this study were to specify a negative binomial mixed model and estimate and contrast the performance of the resulting ICC to that based on estimates from linear mixed models of several data transformations. In addition to the reliability study on environmental sampling to quantify MAP in dairy pens, a wide variety of negative binomial distributed data was simulated to contrast estimator performance.

ICC for the negative binomial mixed model
For the purpose of deriving the ICC that estimates the similarity of samples collected by two veterinarians on the same day from the same pen. Here y ijkl denotes the observed value of the j th pen in i th dairy, the k th day, and the l th veterinarian (i = 1, 2,..m; j = 1, 2,..n m ; k = 1, 2,..s; l = 1, 2,..t); the total number of observations is N = st∑n m . We assume that the conditional distribution of y ijkl given the random effects a, b, c, d (dairy, pen, day and veterinarian, respectively) is distributed negative binomial with the number of successes needed r and probability of success p ijkl , or NB(r; p ijkl ). In this distribution, r is fixed for all y ijkl . Furthermore, it is assumed that μ ijkl ¼ e βþa i þb ijk þc k þd i , where μ ijkl is the conditional expectation of y ijkl given p ijkl . Recall p ijkl ¼ r μ ijkl þr (Casella and Berger 2002), where a i (i = 1, 2,…,m) are independent and distributed as N(0; σ 2 a ), b i (i = 1,2,…n m ) are independent and distributed as N(0; σ 2 b ), c i (i = 1,2,..s) are independent and distributed as N(0; σ 2 c ), and d i (i = 1, 2,.. t) are independent and distributed as N(0; σ 2 d ). Hence the conditional expectation and the conditional probability thereby the conditional distribution of y ijkl is NBðr, r= ðe βþa i þb ij þc k þd l þ rÞ Á which is a special case of the GLMM. The ICC for the similarity in Herrold's egg yolk medium (HEYM) culture results for MAP in samples collected by 2 different collectors (l 1 and l 2 ) will be derived as an example. Given the model assumptions and study design, y ijkl 1 and y ijkl 2 are conditionally independent if l 1 ≠ l 2 , so the conditional expectation of their product is the product of their conditional expectations. Therefore, According to the expectation of the log-normal distribution, we have: Similarly, The covariance between two measurements that are generated by different veterinarians but otherwise are identical in all factors is the difference between the expectation of their product and the product of their expectations, that is: On the other hand, according to the variance of the log-normal distribution, Hence, the expectation of conditional variance of the observations can be expressed as: Therefore, the variance of the observations is: It follows then that the ICC for the negative binomial mixed model for the similarity between samples collected by two different veterinarians on the same day and from the same pen is: When the variance of data is much larger than its expectation, the negative binomial distribution is often used as an alternative to the Poisson distribution. The random effects follow the normal distribution and the link function is the logarithm. Based on this formula, the ICC is no longer just based on the random effects, but also the fixed effect intercept and the number of successes. Thus, the negative binomial mixed model may be more reasonable than the LMM or the Poisson GLMM when count data are overdispersed.

Simulations
Simulations were used to compare the performance of the new negative binomial ICC to that estimated after traditionally transforming count data to normalize it using transformation such as the logarithm, square root, square, or their inverse values (Carrasco and Jover 2005). To compare the performance of the ICC estimator derived for the negative binomial GLMM to the ICC used in traditional methods such as LMM of normalized count data, 16 scenarios were generated. Fixed estimates of input parameters were used in each of the 16 scenarios and their respective true ICC as summarized in Table 1. The scenarios included 2 different estimates of r (r = 1 and r = 2), numbers of successes. In addition, zero and non-zero intercepts (β =0 and β =2), 2 different between-dairy variance estimates (0.5 and 1), and 2 different between-veterinarian variance estimates (0.1 and 0.5) were assumed. The justification behind the use of fixed estimates for the betweenpen and between-day variances is that by equation (1), these variances influence the ICC in the same way, as the between-dairy variance; therefore it is reasonable to vary only one of them. Based on the study by Aly et al. (2009) there were 4 factor levels: dairy, pen, veterinarian and day. Pens were nested within dairy. In dairy i , i = 1, 2, 3, 4 has j pens; where for i = 1, j = 1,…, 8; for i = 2, j = 1,…, 11; for i = 3, j = 1,…, 7; and for i = 4, j = 1,…, 4. Pens were cross-classified by veterinarian l; l = 1, 2; and day k; k = 1, 2, 3. Data were generated under the assumption of negative binomial GLMM with log link using all four factors a, b, c, d included as random effects. Each sample dataset consisted of 180 observations. Each simulation followed the following procedure: 1. Randomly generate normal random effects a i , b ij , c k , d l (i = 1, 2,.. n m ; k = 1, 2.. s; l = 1, 2,.. t) with respective scenario's variances 2. Sum the intercept and random effects as conditional expectation μ ijkl ¼ e βþa i þb ij þc k þd l , β is estimated intercept from field data 3. Randomly generate negative binomial variable Y ijklÑ B(r, μ ijkl ) r is number of successes 4. Estimate model parameters: intercept β, number of successes r and random effects σ 2 a þσ 2 b þσ 2 c þσ 2 d

Calculate the ICC
One hundred simulated data sets were generated under each scenario. For each simulated data set, the ICC was estimated using four different methods: 1) the negative binomial GLMM, 2) LMM of raw data (untransformed); 3) LMM of square root transformed data; and 4) LMM of logarithm transformed data where taking logarithm of zero was avoided by replacing zeros with 0.5. For LMM, restricted maximum-likelihood estimation (REML) was used, while maximum-likelihood (ML) estimation was used for the GLMM. Relative bias, variance of the ICC, and mean square error (MSE) of the ICC estimate were calculated to evaluate the performance of the ICC. The relative bias was calculated as the difference between the mean of estimated ICC and it's true value, variance was calculated by unbiased estimation based on the simulation, and MSE was calculated as the sum of squared bias and variance.
A second simulation explored the performance of the ICC estimate over a wider range. The mean estimated ICC was computed using 400 simulations per combination of number of successes (r = 5 and r = 30) and variance estimates for dairy and veterinarian (0 to 1 in increments of 0.2).

Field data analysis
Finally, field data used in the report by Aly et al. (2009) were analyzed using the negative binomial GLMM. Briefly, environmental samples were collected every other day on 3 different occasions from 4 California dairies between November 2006 and June 2007. Samples were cultured using bacterium-specific medium using standard microbiological procedures as reported by Aly et al. (2009). Confidence intervals for model parameters were obtained based on parameter estimates from the field data and using parametric bootstrap similar to that described in Table 1 (Efron and Tibshirani 1993). The resulting negative binomial based ICC was contrasted to that estimated from transformed data and reported previously by Aly et al. (2009). The R package lme4 was used for LMM analysis, and the package glmmADMB for GLMM analysis. All packages were loaded in the R 2.15.1 environment.

Results
Results of the first simulation targeted a range of ICC values based on 16 combinations of input parameters (r, β, variances of dairy, pen, veterinarian and day) and are presented in Table 2. The relative bias in the ICC, variance and MSE were compared for the ICC estimates based on the negative binomial model to those based on the LMM of raw (untransformed) and transformed data. The negative binomial model ICC had the least absolute relative bias in 5 of the 16 scenarios (3, 4, 5, 6 and 8) that were characterized by small variance estimates for veterinarian (0.1). In comparison, the ICC based on LMM of raw data had the most number of scenarios with the least absolute relative bias (9, 10, 11, 12, 14, 15 and 16) characterized by large variance estimates for veterinarian (0.5). In terms of variance, the negative binomial ICC had the most number of scenarios with the least variance (1 to 5,7,8,11,12,15). Similarly for MSE, the ICC based on the negative binomial model had the least MSE in 11 of the 16 scenarios (1 to 8, 11, 12, 15 and 16).
The second simulation performed to investigate the effect of larger number of successes (r = 5 and r = 30) and a wider range of variance estimates for dairy and veterinarian that also include zero. Figure 1 showed that the mean of the estimated ICC and the true ICC were similar as estimates of variance due to veterinarian ranged from 0.1 to 0.3 even as variance due to dairy increased to 1. However, as depicted in the diverging planes, the difference between the estimated and true ICCs increased towards extreme variance estimates. Both behaviors were consistent in a higher number of successes (r = 30). Figure 1 depicts the differences between the true ICC and the mean of the respective estimated ICC based on simulations.
Results of the negative binomial GLMM are summarized in Table 3. The negative binomial ICC was estimated to be 0.5207 (95% CI = 0.4033, 0.6091) compared

Discussion
The current study updates an earlier report on the reliability of environmental sampling to quantify MAP in freestall dairy pens utilizing the negative binomial ICC for count data. A unique character of the negative binomial ICC is the inclusion of the fixed effect intercept estimate unlike the ICC based on LMM which is based soley on variance components. Fixed effects are similarly included in the formula for the poisson ICC however the negative binomial ICC also includes r, the distribution parameter for number of successes. A performance comparison of the ICC estimates showed that the negative binomial ICC was more suitable for count data that is overdispersed given the smaller MSE and variance estimate than the ICC from LMM. Relative bias tended to the least in more scenarios (7 out of 16) with LMM compared to the GLMM based ICC. The lower relative bias with LMM may be explained by the use of REML estimation. The choice of MLE for GLMM was justified by that REML for GLMM has not been well defined, unlike for LMM (Jiang 2007). Nevertheless, the ICC for the   random effect for pen j, where for i = 1, j = 1,…,8, for i = 2, j = 1,…,11, for i = 3, j = 1,…,7 and for i = 4, j = 1,…,4. c random effect for day k of sample collection, where k = 1,2,3. d random effect for collector l, where l = 1, 2 and day k; k = 1, 2, 3.