Empirical Bayes estimation of farm prevalence adjusting for multistage sampling and uncertainty in test performance: a Brucella cross-sectional serostudy in southern Kazakhstan

Summary

Estimation of farm prevalence is common in veterinary research. Typically, not all animals within the farm are sampled, and imperfect tests are used. Often, assumptions about herd sizes and sampling proportions are made, which may be invalid in smallholder settings. We propose an alternative method for estimating farm prevalence in the context of Brucella seroprevalence estimation in an endemic region of Kazakhstan. We collected 210 milk samples from Otar district, with a population of about 1000 cattle and 16 000 small ruminants, and tested them using an indirect ELISA. Individual-level prevalence and 95% confidence intervals were estimated using Taylor series linearization. A model was developed to estimate the smallholding prevalence, taking into account variable sampling proportions and uncertainty in the test accuracy. We estimate that 73% of households that we sampled had at least one Brucella-seropositive animal (95% credible interval 68–82). We estimate that 58% (95% confidence interval 40–76) of lactating small ruminants and 14% (95% confidence interval 1–28) of lactating cows were seropositive. Our results suggest that brucellosis is highly endemic in the area and conflict with those of the official brucellosis-testing programme, which found that in 2013 0% of cows and 1·7% of small ruminants were seropositive.

INTRODUCTION

Brucellosis is a bacterial zoonosis that reduces reproductive performance and milk production in cattle (predominantly Brucella abortus), sheep and goats (predominantly B. melitensis) and other livestock species [1]. Both B. abortus and B. melitensis infections cause a range of syndromes in humans that include fever and joint pain and range from mild to debilitating. Brucella sp. can be transmitted to humans via contact with aborted fetuses, parturition fluids or via consumption of unpasteurized milk and dairy products. [2]

Similarly to other Central Asian countries, Kazakhstan has a high human incidence of brucellosis. According to reports from the Ministry of Health, nine cases were reported per 100 000 humans in 2012. The reported incidence is higher in the south of the country, for example in Almaty Oblast (region) where incidence increased between 2007 and 2010, from 19·5 to 30 cases/100 000 [3, 4]. Despite widespread testing of livestock, the percentage of positive tests has been consistently very low in all oblasts according to data provided by the Kazakh Ministry of Agriculture in Astana (generally <1% of cattle and small ruminants). Considering the high human incidence, and the likelihood of under-reporting of human cases, this suggests there may be either a systematic bias in selection of animals for testing, an inaccurate testing regimen or inaccurate reporting. There are very few alternative reliable data sources on the current prevalence and distribution of the disease in livestock or humans in Kazakhstan [5–8].

Control of livestock diseases, particularly zoonoses, is important in Kazakhstan because a large proportion of the population relies on small-scale agriculture for subsistence, particularly in rural areas [9]. Brucellosis has been identified as a priority by the Kazakh Ministry of Agriculture and the current control strategy is in the process of being revised.

Accurate estimation of herd-level prevalence (or seroprevalence) is essential to the planning and implementation of cost-effective disease control programmes. Theoretical aspects of defining ‘herd status’ based on testing of individual animals were reviewed by Christensen & Gardner [10].

In prevalence surveys of livestock diseases, often not all animals in each household are sampled, the number and proportion of animals sampled varies per household, the test(s) used is (are) not perfect and the sensitivity and specificity of the test(s) are uncertain. It is not uncommon to ignore these potential biases when calculating herd-level prevalence, and if the sampling fraction within each herd is high, and the test has a high sensitivity and specificity, this may be justified [11, 12]. However, if this is not the case, it is good practice to adjust ‘apparent’ prevalence to generate ‘true prevalence’ estimates that account for potential misclassification. One approach for herd-prevalence estimation has been to estimate the herd-level sensitivity and specificity from individual-level sensitivity and specificity values, and then to calculate a so-called ‘true prevalence’ taking into account the likely numbers of false-positive and false-negative herds based on the values of herd-level sensitivity and specificity [13–15]. However, in order to calculate a single herd-level sensitivity and specificity it is necessary to make the implicit assumption that the herd-level sensitivity and specificity are the same for each and every herd/household, which is unlikely to be the case when there is variability in herd/household size and the number of animals sampled in each herd/household. More critically, a somewhat arbitrary cut-off has to be chosen for ‘design prevalence’, a threshold below which it is assumed the disease cannot be present. The trade-offs between sample size, ‘design prevalence’ and herd-level sensitivity and specificity have been explored in detail elsewhere [16]. Simulation has been used to account for uncertainty in the true herd-level sensitivity and specificity, but it has still been assumed that there are a single set of true values that are applied to every herd [13, 14].

These problems are exacerbated by very small herd sizes, which often occur in smallholder settings where brucellosis and other livestock diseases are often most prevalent. For example, when testing one cow out of three in a household, if the cow tests negative, it is nonsensical to state that the household is ‘negative’, or that the prevalence is ‘0·05 after adjusting for herd-level sensitivity and specificity’. However, it does make sense to say that the house has a given probability of being ‘negative’ (meaning that all three cows are negative). If this approach was applied to each household/farm in the study one-by-one, household-level prevalence could then be estimated. This general approach has been used widely in risk assessment. An extension of this probabilistic approach is Bayesian estimation of prevalence, which has been reviewed by Branscum et al. [17], who also proposed a method for estimation of herd-level prevalence implemented in WinBUGS, which was adopted by Verdugo et al. [18], Pruvot et al. [19] and others; however, the method relies on the Binomial approximation for sampling of animals, which is not suitable for small herds. Suess et al. [20] simulated the true status of each animal within each herd, removing the need to assume a certain herd-level sensitivity and specificity; however, in their model they assumed that the same number of animals was sampled within each herd. In both cases, prior distributions for prevalence and uncertainty distributions for sensitivity and specificity were generated using the opinions of experts.

Here, we present a method of calculating exact probabilities of positive household status for each household, one-by-one, based on a discrete (rather than continuous) probability distribution of the number of true positives in the household, that is generated separately for each household. Then in a second step we estimate household-level prevalence. We make no assumptions about the sizes of the herds or the numbers of animals sampled per herd. The prior distribution of within-household prevalence for each herd is a discrete distribution that gives the probability of each possible number of positives in that herd, given the number of animals in the herd. This prior distribution is generated from the within-household prevalence distribution of the other herds in the survey, and sensitivity and specificity uncertainty distributions are generated from the data used by the manufacturer to validate the test. The self-contained program runs in R [21] and the only input required by the user is the input of the survey data and the sensitivity and specificity validation data (or point values of sensitivity and specificity, if preferred).

The aims of this study were: (1) to estimate the prevalence of brucellosis antibodies in milk, from cattle and small ruminants in a typical rural village in one region in the south of Kazakhstan; and (2) to develop a method for estimating true herd-level prevalence taking account of a range of sampling fractions used on each smallholding and uncertainty in the sensitivity and specificity of the test used.

METHODS

Study area

We selected Otar Selskiy Okrug (district) for our study site as it is a typical rural Kazakh district, is conveniently located near to the laboratory and the necessary permissions from the local veterinarians, regional veterinary office and Ministry of Agriculture were granted. Otar is also in southern Kazakhstan where the human incidence of brucellosis is high. It has a population of 10 759 humans, 1054 cattle and 16 050 small ruminants (sheep and goats), many of which are kept by smallholders. These animals are kept at the household during the night and share grazing around the villages during the day. In Otar Selskiy Okrug there are 1525 households with livestock and 1300 of these are in the main village, Otar. There is only one other large village, Matybulak (with 110 households with livestock), which is located ~3 km from Otar village (in Almaty Oblast) and animals from both villages share grazing areas. Throughout Kazakhstan, brucellosis vaccination is prohibited, and a national test-and-slaughter programme is being carried out, that involves twice-yearly testing of all sheep, goats and cattle.

Original study design

We planned to conduct systematic random sampling of households. We estimated that we needed to sample 250 cows and 360 small ruminants in order to estimate individual-level prevalence within 1·7% absolute error with 95% confidence, with a hypothesized prevalence of 2%, a design effect of 1·4 and a finite population size of 500 lactating cows and 8000 lactating small ruminants. We planned to select the required number of animals needed to detect disease (assuming a minimum of 10% seroprevalence) with 95% confidence in each household according to the number of livestock present.

Amendments to study design

Due to practical limitations, we had to resort to convenience sampling of households. The households were selected either by the local veterinarians because they were already planning to visit them for routine brucellosis blood testing, or they were relatives or friends of the research team, or they were neighbours or friends of these people. In addition we sampled four large farms in the area, one of which was the research farm belonging to the laboratory and the remaining three were contacts of the research team. These animals grazed on the steppe during the day and were brought into an enclosure at night.

Due to limitations in accessing households, we sampled as many as possible of the cattle, sheep and goats on each household or farm. We tried to avoid any obvious bias in selection of animals for sampling, but random sampling was not possible.

Data collection

The identification number or description, species, breed and age of each animal was recorded when available. We completed an interview (in Russian or Kazakh) with each owner, using a pre-designed form including questions on the number of animals of each type on the household and the gender and age of members of the household who regularly milked the livestock or assisted with parturition.

Sample collection, processing and testing

Ten millilitres of milk were collected from each quarter into a single plain polyethylene tube, after cleaning and drying the teats. The samples were placed immediately into a cool box and placed in a refrigerator within a few hours. The samples were left to stand to allow the lactoserum to separate from the fat layer, or they were centrifuged, and the lactoserum was pipetted into Eppendorf tubes. The samples were then frozen for up to 5 months before de-frosting at room temperature.

The milk samples were tested using an indirect ELISA for brucellosis antibodies (ID Screen^® Brucellosis Milk Indirect; ID.vet, France) according to the manufacturer's instructions.

Individual-level seroprevalence

The individual-level percentages of test-positive milk samples were estimated using the ‘survey’ package [22] in R [21], which produces both a point estimate that is weighted according to the sampling fraction, and a confidence interval that is adjusted for clustering (in this case within a smallholding or farm) by adjusting the standard error using Taylor series linearization [23]. The sampling fraction was calculated as the number sampled/the number of lactating animals in the household at the time of the visit. There were seven households for which the number of lactating small ruminants in the household was not recorded, and five for which the number of lactating cows was not recorded. The mean of the available sampling fractions was used for these households. Uncertainty distributions for the sensitivity and specificity of the iELISA were generated based on data used for validation by the manufacturer (see Table 1).

Table 1. Sensitivity and specificity values assumed in this analysis

CFT, Complement fixation test; RBT, Rose Bengal test; PAT, plate agglutination test.

The point estimates and confidence intervals were adjusted for the sensitivity and specificity using the uncertainty distribution, in R, using the following formula:

$${\rm TP} = {\rm} \left( {{\rm AP} + {\rm Sp} - 1} \right)/\left( {{\rm Se} + {\rm Sp} - 1} \right),$$

where TP = true prevalence, AP = apparent prevalence, Se = sensitivity and Sp = specificity. [24].

The median values of the resulting uncertainty distributions for point estimate, lower and upper confidence intervals are presented.

Household-level prevalence

The household-level prevalence of brucellosis (the proportion of households with at least one animal with antibodies against Brucella spp. in the milk) was estimated taking into account the imperfect sensitivity and specificity of the test and uncertainty arising from (1) estimation of the true sensitivity and specificity from a sample of ‘truly infected’ and ‘truly non-infected’ animals and (2) sampling of only a proportion of each household (the proportion being different in each household). A model was constructed in R (R Core Team, 2015) as follows. (The complete model is available as Supplementary material.)

Step 1. The probability that each given household was negative (P _n ) (i.e. that there were no lactating animals with antibodies in the milk) was calculated for each household individually as follows:

$$P_n = O_n/\left( {1{\rm} + O_n} \right),$$

where O _n  = the (posterior) odds that the household was negative, and was calculated according to Bayes theorem as follows:

$$O_n = {\rm} \left( {{\rm Prio}{\rm r}_n \times {\rm Likelihoo}{\rm d}_n} \right)/\sum \left( {{\rm Prio}{\rm r}^i \times {\rm Likelihoo}{\rm d}^{ijk}} \right),$$

where Prior _n  = the prior probability of that there were zero positives on the farm (this was a discrete probability for each iteration of the model); Likelihood _n  = the likelihood of obtaining the laboratory results, given that there were zero positives on the farm; Prior ⁱ  = the prior probability that there were i positives on the farm (where i is a vector from 1 to the total number of lactating animals on the farm) (these were discrete probability values for each iteration of the model); and Likelihood ^i,j,k  = the likelihood of obtaining the laboratory results, given that there were i disease-positives in the household, j disease-positives among the tested animals and k false positives. For each household, each possible permutation of the number of true positive lactating animals in the household ^† (i) and the sample (j), and the number of false positives (k) that could result in the laboratory results was generated programmatically and the likelihood was calculated for each permutation, as shown in the worked example in Figure 1.

Fig. 1. Worked example showing how Prior _n, Likelihood _n, Prior ⁱ and Likelihood ^i,j,k were calculated. In this example household, there were three lactating animals, two of which were sampled, and there was one positive result. For demonstration purposes, a sensitivity of 0·99 and a specificity of 0·99 were used; however, in the final model, an uncertainty distribution was used.

The prior distribution of the number of true positives on the farm was also calculated for each household separately as follows. First, the frequency distribution of within-household prevalence values from all households in this study was multiplied by the number of lactating animals in the particular household, and second, the resulting numbers of positives were rounded to whole numbers, to create a discrete probability distribution of each possible number of true positives on the farm.

Step 2. For each run of the model, each household was simulated to be positive or negative by drawing one random sample from a binomial distribution with probability of success of (1 – P _n ).

Step 3. Steps 1–2 were repeated 1000 times to create an uncertainty distribution, where the 2·5th and 97·5th percentiles give a 95% credible interval and the 50th percentile gives the most likely household-level prevalence.

Comparison of field-study results with official seroprevalence data

In order to compare the field-study results with the official data, we conducted the following analysis. The number of official tests conducted in 2013 and the number of seropositives were obtained from the local veterinary office. It was assumed that the results were obtained by combining the plate agglutination test (PAT), complement fixation test (CFT) and Rose Bengal test (RBT) in series. Sensitivity and specificity estimates for each test were obtained based on a published meta-analysis [25] where available, or literature review (Table 1). Combined sensitivity and specificity were estimated as shown in Table 1 (which assumes the tests are independent from one another).

Based on the combined sensitivity and specificity, adjusted seroprevalence estimates were obtained using the method by Reiczigel et al. [26] implemented online in ‘Estimated true prevalence and predictive values from survey testing’ (http://epitools.ausvet.com.au/content.php?page=TruePrevalence). This assumes that the results were obtained by simple random sampling.

The difference between the adjusted seroprevalence values obtained based on the field-study and official results, and an exact 95% confidence interval (CI) was estimated using the following formula:

$$\eqalign{{\rm 95}\% \,{\rm CI} &= {\rm RD} - \surd \left[ {{(\,p_{1} - l_1)}^2 + {\rm} {\left( {u_2 - p_2} \right)}^2} \right] \cr &\quad\ \hbox{to RD } + {\rm} \surd \left[ {{\left( {\,p_2 - l_2} \right)}^2 + {\rm} {\left( {u_1 - p_1} \right)}^2} \right],} $$

where RD = the risk difference, l ₁ to u ₁ is the 95% CI of the first proportion, p ₁ and l ₂ to u ₂ is the 95% CI of the second proportion, p ₂ [27].

Ethical standards

Ethical approval for the study was granted by the Royal Veterinary College ethics committee (URN 2011 1097). Written consent was obtained from one member of each household.

RESULTS

Descriptive analysis

Farms and households

In total, three farms and 31 smallholdings in Otar Selskiy Okrug were visited. Interviews were conducted on all three farms, 22 of the smallholdings from which milk samples were taken and a further three smallholdings which had no lactating animals at the time of the visit.

The farms were all mixed species and had between one and three lactating cows (median 2) and 150–318 small ruminants (median 284). Of the 25 smallholdings for which data were available, seven kept only small ruminants (range 12–40, median 30), three kept only cattle (range 1–3, median 2) and 15 kept both small ruminants (range 7–85, median 20, data missing for two farms) and cattle (range 1–16, median 2).

A range of 1–5 individuals on each smallholding had contact with the livestock via milking or delivering newborn animals (median 2, missing data for 7/31 smallholdings). Of these individuals (n = 48), 21 were female and 27 were male, and their ages ranged from 13 to 75 years (median 47, missing data for 5/48 individuals).

Milk samples

We collected 210 milk samples from all three farms and 28 of the smallholdings. The milk samples came from 43 cows and 167 small ruminants (129 sheep, 23 goats, 15 not specified). The cows ranged in age between 3 and 15 years (median 5, data missing for 15/43 cows); the small ruminants ranged in age between 1 and 7 years (median 4, data missing for 86/167 small ruminants). The breeds of the livestock were frequently unknown or not recorded, but cattle breeds included mixed, Kazakh, Zerno-pestreesa and Alatau; sheep breeds included mixed, Kazakh, Merino and Yedilbai; and goat breeds included mixed and Angora.

After accounting for the multistage sampling that resulted in variable sampling fractions with households, it was estimated that 13·6% (95% CI −1·1 to 28·3) of lactating cattle and 57·9% (95% CI 40·2–75·6) of lactating small ruminants were seropositive (Table 2). The percentage of official serological tests reported positive in Otar Selskiy Okrug was 53% lower (95% CI 38–73) for small ruminants, suggesting the results were incompatible with one another. For cattle our results were 13·2% lower than the official serological results but the 95% CI included zero, suggesting that the results were not incompatible at the 95% level, although our sample size was possibly not large enough to detect a statistically significant difference.

Table 2. Seroprevalence of brucellosis according to this field study of milk samples and official blood testing in Otar Selskiy Okrug

CI, Confidence interval.

The apparent percentage of households with at least one seropositive lactating animal was 64% (14/22 households for which sufficient data were available) and the adjusted percentage was 73 (95% credible interval 68–82). (Note, this is a sample estimate not a population estimate.)

DISCUSSION

The results strongly suggest that brucellosis is endemic in livestock in Otar Selskiy Okrug and the seroprevalence in small ruminants is much higher than suggested by official results. (A larger sample of cows would be needed to estimate the seroprevalence in cows more precisely.) It is unlikely that this difference is due to systematic bias in selection of animals for official testing, due to the large numbers officially tested, and the size of the difference between our results and official results. It seems unlikely that selection of animals in our study was heavily biased towards animals with brucellosis, as owners would be less likely to allow testing if disease was suspected, and many of the households were known to have been included in the official testing programme. It also seems unlikely that such a large difference would be observed due to testing milk as opposed to blood samples.

Finally, we used an indirect ELISA, whereas the official testing involves the PAT, RBT and CFT. We used a commerical indirect ELISA to detect Brucella antibodies because there were good laboratory facilities available to perform this test, and the laboratory staff had experience with ELISAs. The agglutination tests (RBT and PAT) and CFT also measure antibodies; however, sensitivity and specificity values vary between the different tests, and there can be slight differences in the sensitivities to different classes of antibody (IgM, IgG or IgA). There is no gold standard test for detecting antibodies to Brucella [28]; however, all of the tests have been validated by previous authors using samples from known infected and non-infected animals [25]. We adjusted the official prevalence values according to published sensitivities and specificities, and the ELISA results were adjusted according to sensitivity and specificity data provided by the manufacturer of the ELISA. These values were based on cows and B. abortus only, and the cows were from different contexts to Kazakhstan. Sensitivity and specificity values can vary in different populations. Data specific to Kazakhstan or to small ruminants generally are lacking, and it is possible that there is a higher incidence of cross-reactions in our study; however, it seems unlikely that this would explain the high proportion of positives that we found.

About 33 cases of (culture-positive) brucellosis were reported per 100 000 people in 2009 in Kordai rayon (a larger administrative division than Selskiy Okrug – data not available at Selskiy Okrug level), which would suggest that a higher seroprevalence in livestock would be expected than is officially reported, based on data from other endemic countries.

Possible explanations for the difference in results include poor sensitivity of tests used in official laboratories, or false reporting. The stringent requirement for all three official tests to be positive for a sample to be classified as positive means that a deficient sensitivity in any one of the tests could result in a very low combined sensitivity. The consequences of culling a positive animal may be very severe for households in Otar Selskiy Okrug, despite the compensation that is given, as the owners are heavily reliant on a few animals, and this may help to explain a bias towards classifying animals as negative in doubtful situations.

There were several additional limitations to this study. Despite planning to do random sampling, we had to resort to convenience sampling. However, obvious biases were avoided where possible, and we sampled typical households in a village setting, as far as we could tell. In addition, we were sampling within a small geographical area with a population of ~ 16 000 small ruminants, most of which share grazing pastures. Although there could be some bias in the sampling, it seems reasonable to conclude that brucellosis is highly endemic in the area, given the high proportion of positives that we obtained.

We conclude that brucellosis is highly endemic in Otar Selskiy Okrug and due to the structure of the livestock system, a large proportion of people in the area are likely to be at risk of exposure to Brucella. The official census data shows that there are >1500 households with livestock in Otar Selskiy Okrug, an area with a human population of ~ 11 000. Our interview data suggested that in most households more than one person had regular contact with the household livestock, and the ages of those in contact ranged from 13 to 75 years.

The methods we developed for estimating the herd-level prevalence could be applied to any multi-stage prevalence study in order to account for variable sampling proportions and imperfect tests, simultaneously. The method also incorporates uncertainty arising from estimation of sensitivity and specificity, based on a small validation study. Furthermore, this approach could readily be tailored to various study designs, could be extended to use in analytical studies aimed at identifying risk factors, or could also be used as a ‘module’ in a quantitative risk assessment involving heterogeneous farms, for example.

Brucellosis control has historically been extremely challenging in Central Asia, and there is very little precedent for the success of test-and-slaughter in endemic areas [5, 29]. A wide-scale review of the national test-and-slaughter scheme in Kazakhstan is needed, bearing in mind that Kazakhstan is a large country and our study was only conducted in a small geographical region.