Risk factors and variations in detection of new bovine tuberculosis breakdowns via slaughterhouse surveillance in Great Britain

Slaughterhouse surveillance through post-mortem meat inspection provides an important mechanism for detecting bovine tuberculosis (bTB) infections in cattle herds in Great Britain (GB), complementary to the live animal skin test based programme. We explore patterns in the numbers of herd breakdowns detected through slaughterhouse surveillance and develop a Bayesian hierarchical regression model to assess the associations of animal-level factors with the odds of an infected animal being detected in the slaughterhouse, allowing us to highlight slaughterhouses that show atypical patterns of detection. The analyses demonstrate that the numbers and proportions of breakdowns detected in slaughterhouses increased in GB over the period of study (1998–2013). The odds of an animal being a slaughterhouse case was strongly associated with the region of the country that the animal spent most of its life, with animals living in high-frequency testing areas of England having on average 21 times the odds of detection compared to animals living in Scotland. There was also a strong effect of age, with animals slaughtered at > 60 months of age having 5.3 times the odds of detection compared to animals slaughtered between 0–18 months of age. Smaller effects were observed for cattle having spent time on farms with a history of bTB, quarter of the year that the animal was slaughtered, movement and test history. Over-and-above these risks, the odds of detection increased by a factor of 1.1 for each year of the study. After adjustment for these variables, there were additional variations in risk between slaughterhouses and breed. Our framework has been adopted into the routine annual surveillance reporting carried out by the Animal Plant Health Agency and may be used to target more detailed investigation of meat inspection practices.


Introduction
Bovine tuberculosis (bTB) is the most economically important disease of livestock currently affecting cattle in Great Britain (GB), estimated to cost the UK government >£100 million per PLOS  in GB is manifest in the consistency of the distribution of reactor animals removed from breakdowns in different risk areas [9], and in its critical role in Scotland and low-risk areas of England.
In the Republic of Ireland (ROI), separate logistic regression models were estimated for the individual animal risk of submission of suspect lesions and laboratory confirmation [10,11]. The goal of these Irish studies was to provide separate (ranked) estimates for the rates of submission of suspect lesions and laboratory confirmation per slaughterhouse respectively, after adjusting for potential confounding risk factors in the input population of each slaughterhouse. These studies estimated the impacts of individual-and herd-level effects of age, sex, season and indirect measures of exposure to bTB (e.g. risk class, number of years clear of restrictions, whether animal was homebred), along with a fixed effect representing the relative difference between each slaughterhouse. In general the individual level risk of both submission and confirmation was found to increase with age and risk of exposure, with a slightly increased risk of confirmation for female animals (which could be confounded with differences between dairy and beef management practices).
Similar research was conducted in Northern Ireland (NI) [12], again focusing on risk factors for disclosure of bTB in slaughterhouses, and for the probability of confirmation. This study looked at the effects of age, sex, season, patch incidence (e.g. risk region), past history of bTB on the disclosing herd and whether or not the animal was homebred or purchased. In general the risk of disclosure increased with age and for animals sourced from high-risk areas, with a higher risk of disclosure in females, purchased animals and animals slaughtered in winter. There was a slightly more complex relationship with herd history. We note that in ROI and NI the number of slaughterhouses (42 and 10 respectively) is much smaller than in GB (276).
In GB, an alternative approach was used to model the risk of detecting M. bovis in samples from lesions found in cattle that were consistent with being infected; quantifying the variability between slaughterhouses using a random intercept term [13]. Their study population consisted only of animals that had been sent to slaughter from Officially TB Free (OTF) herds, and individual level factors adjusted for by the model included year of birth (rather than age), year of slaughter (2003)(2004)(2005)(2006)(2007)(2008) and further indirect measures of risk of exposure to bTB (e.g. region of farm of origin of bovine, parish testing interval of farm of origin, and number of reactors in most recent previous breakdown on farm). In common with the Irish studies, the risk of confirmation increased with indirect measures of exposure and a large degree of variability in confirmation rates between slaughterhouses was reported, ranging from 21-81%, although this could reflect varying prevalence of other, confounding, intercurrent disease. In their final model the probability of confirmation was found to decrease over time with year of slaughter, but increase with date-of-birth.
In this manuscript we tackle a slightly different question and quantify the contribution that slaughterhouses make to the overall detection rate of bTB breakdowns in GB, rather than the rate of submission or confirmation of individual animals. To this end we define our cases as all animals with suspicious lesions submitted for laboratory confirmation that were either subsequently confirmed, or that triggered the disclosure of new skin test reactor animals in the originating herd. Our case definition therefore captures the overall effectiveness of slaughterhouse surveillance at the herd, rather than the individual animal level. In addition to hierarchical terms that estimate the variability between slaughterhouses, we also include an additional hierarchical term to account for possible confounding effects between breeds in terms of susceptibility [14], reaction to diagnostic tests and post-mortem confirmation [15].

Materials and methods
Bovine TB is a notifiable infectious disease in GB. Surveillance data recording the results of all testing and slaughterhouse cases are collated by the Animal and Plant Health Agency (APHA) in the Sam database and linked with cattle movement records from the Cattle Tracing System (CTS). Data downloads from the Sam system were used to construct local databases using PostgreSQL (https://www.postgresql.org/) with subsequent analysis and modelling carried out using the open-source statistical language R [16].

Number of breakdowns per herd over time
The numbers of OTFS and OTFW breakdowns identified between 1998-2013 were extracted, along with the associated disclosing test responsible for initiating each breakdown (55,653 entries). Data where first tests were missing or unknown, or 'queried' data (entered as 'VE-QSLH') were removed (39 entries), leaving 55,614 entries.

Risk factors for identifying a breakdown through routine slaughterhouse surveillance
All slaughterhouse (SLH) cases that disclosed a breakdown and submitted between 2004 and 2013 for laboratory confirmation were extracted as our study population. This period was chosen to avoid bias resulting from the well described impact of the suspension of SICCT testing during the 2001 foot-and-mouth epidemic [2]. Since one of the principal aims of the study was to quantify differences between slaughterhouse detection rates, over-and-above the background risk of detection based on individual-level risk factors, it was necessary to sample controls proportionately between slaughterhouses. This ensures that we are able to obtain unbiased estimates of the odds ratios for each slaughterhouse, relative to an average slaughterhouse (see model specification below). In addition, since there is considerable variability in the throughput of slaughterhouses, a simpler, non-stratified random sample would have resulted in some slaughterhouses having either no, or very few, controls. Since we were using all the available cases, the overall statistical power can also be improved by sampling disproportionately from the control population. Hence we extracted a stratified random sample of 0.3% of skin test negative animals with no visible lesions at slaughter from each slaughterhouse as our control population. This ensured reasonable coverage within each slaughterhouse and that the model was computationally manageable. After records with other missing data were removed, we ended up with 77,426 samples, consisting of 8,607 cases and 68,819 controls. Summaries of the raw data can be found in the S1 File, and a copy of the raw data can be found in the S1 Data file.
Risk was modelled at the individual animal level using a range of explanatory factors ( Table 1) that capture the temporal and individual variations in risk of exposure and presentation at slaughterhouses. In common with previous studies: sex, seasonal quarter, region and age are included as baseline risk factors [10][11][12][13]. Through linking the Sam and CTS databases a range of alternative indirect measures of exposure were calculated for each animal that capture different aspects of the lifetime risk-of-exposure for each animal.
The simplest of these is the number of lifetime herd-to-herd movements. This factor distinguishes the differential risk of exposure between cattle that remain in the same herd in their lifetime (no moves) or are traded (one or more moves). We aim to approximate the duration of time that an animal may have been exposed to bTB through the number of days spent in a high-risk herd. A previous history of bTB is a consistent risk factor for predicting future breakdowns [9,13,17,18]. We choose to define high-risk herds to be any herd with a history of bTB restrictions within the 5 years previous to the date of the slaughterhouse disclosure. In contrast, low-risk herds are defined as those with no reported bTB incidents in the 5 years prior to disclosure. The duration distributions are heavy-tailed, and after taking log(duration+ 1) transformations, the duration of time in high-risk herds (in particular) is strongly bimodal (Fig A in S1 File). Hence for simplicity we decided to simplify this measure into two categorical variables: the first describes the duration of stays on low-risk herds as "None" or "Some", and the second is similar but for stays on high-risk herds. Finally, the number of herd tests that an animal was present for is included as a proxy measure for both the risk of transmission within the herd, but also the likelihood that an infected animal could be identified by SICCT testing before slaughter.
To specify the model, we define the binary response variable: and assume Y i * Bernoulli(p i ). Here, where the x i are a vector of explanatory variables (shown in Table 1 Heterogeneity in slaughterhouse detection rates is modelled through the estimation of slaughterhouse-level hierarchical intercepts. This approach was taken to allow for the highly variable throughput of each SLH, with some slaughterhouses slaughtering many more animals than others. This results in a shrinkage model, where SLHs with low throughput will be shrunk towards the overall mean, thus de-emphasising the intercept term for slaughterhouses with low throughput (see e.g. [19]). We formulated the model in a Bayesian framework, which allows full posterior distributions for the parameters to be obtained, as well as full posterior predictive distributions for various measures of interest. Due to the size of the data set under study we wrote our own Markov chain Monte Carlo (MCMC) code that allowed for parallelisation of the likelihood function to improve efficiency. We also increased efficiency by using a mix of centred and non-centered re-parameterisations [20]-see S1 File for full details. For ease of reproducibility, the code has been compiled into an R package called BayesLog, the current version of which can be found at https://github.com/tjmckinley/BayesLog. A link to the data can be found next to the published manuscript online, and the R code necessary to run the models can be found in the S1 File.
To complete the Bayesian specification, we set the prior distributions as: We ran two models: both including and excluding the slaughterhouse terms. Each model was run for two chains of 10,000 iterations each, with the first 5,000 samples from each chain discarded as burn-in. Convergence was assessed using visual inspection of MCMC trace plots. The fit of both models was examined in two ways: by using binned residual plots [21], and by exploring the goodness-of-fit using ROC-curve type analyses [22]. The goodness-of-fit measures were generated from 1,000 posterior samples, chosen at random. All results are reported to 2 significant figures.

Number of breakdowns over time
We first explore the contribution of slaughterhouse surveillance to the number of bTB breakdowns stratified by the parish testing interval (PTI) as a proxy measure of risk. (Note that we exclude the rare historical PTI 3 herds in these plots. The uncommon PTI 2 areas were abolished in January 2013, with herds now classified as being in high-risk (PTI 1) or low-risk (PTI 4) regions [4]). Fig 1 compares the number of OTFS breakdowns to the number of OTFW breakdowns initiated either through slaughterhouse surveillance or SICCT testing for each year.
The absolute number of breakdowns initiated in the slaughterhouse has increased across all three risk areas over the period of the study. However, this increase must be balanced against both the overall increase in the number of breakdowns over the period and the systematic shift of herds from the PTI 2 and 4 regions into PTI 1 accompanying this growing epidemic. This can be accounted for by considering the proportion of breakdowns initiated by different tests (Fig 2).
Overall, the proportions of total breakdowns that are OTFW are relatively stable, with some oscillation (Fig 2); as previously reported by APHA [3]. However, there was a continuing systematic increase in the proportion of OTFW breakdowns being initiated through SLH surveillance (Fig 2) across all risk areas.
A key feature of these data, although not directly related to slaughterhouse surveillance, is that there were proportionally more OTFS breakdowns in PTI 4 than in PTI 1, consistent with there being a lower positive predictive value of the screening skin test in low-risk areas compared to high-risk areas [5].
Logistic regression model fits. Trace plots for both models are shown in the S1 File, and did not show any evidence for a lack-of-convergence. The binned residual plots show no systematic lack-of-fit (see Figs K and O S1 File). (We note some possible curvature in Fig O S1 File, but feel that the fit is still acceptable in this circumstance). The posterior predictive mean AUC (Area Under the ROC Curve) and 95% credible interval is 0.88 (0.88, 0.88) for the full model, and 0.85 (95% CI: 0.85-0.85) for the model without SLH effects. However, we must be careful here, since the data are unbalanced, and so for the full model we also generated a curve exploring the relationship between the positive and negative predictive values for a range of probability thresholds (analogous to the ROC curve- Fig 3). We can see that the specificity and negative predictive value (NPV) are always high, but there is more variation in the sensitivity and positive predictive value (PPV). This is to be expected when the ratio of controls to cases in the data is unbalanced in favour of controls. Nevertheless, there are probability thresholds that can find a good balance between these values (e.g. a probability threshold of p = 0.35 gives high NPV and specificity values, but also a sensitivity and PPV of around 45-55%-as shown in Figs E and F S1 File). Overall we conclude that the model fits well and is a good classifier for these data.
Risk factors. Odds ratios (ORs) for the explanatory variables obtained from the full model (including the SLH terms) are shown in Table 2. (For brevity we exclude the hierarchical terms from Table 2). Note that we have set the baseline level of each categorical variable to be the level with the lowest risk. This means that all ORs are positive, and can be directly compared on the same scale.  The strongest single risk factor is the region where animals spend most of their lives, with (as expected) animals that spend most of their lives in high-risk areas showing greater odds of being SLH cases than animals that spend most of their lives in low-risk areas. Scotland has the lowest odds, with LFT areas in England and Wales showing odds that are on average 6.6-7.2 times larger than in Scotland. Compared to Scottish animals, the odds of being a SLH case are on average 21 times higher for animals living mostly in the HFT region of England, and 12 times higher for animals living mostly in the HFT region of Wales. Qualitatively similar effects were observed for analogous risk region (patch) variables in earlier studies [10,12,13].
As expected, there is a strong increasing effect of age, such that older animals are more likely to be picked up as slaughterhouse cases than younger animals (the odds of animals >60 months old being SLH cases are on average 5.3 times higher than the odds for than animals (0, 18] months old). Again, this general increase of risk with age was also observed in earlier studies in ROI and NI [10,12], and is likely a proxy for an increased probability of infection, and also an increased duration of infection (thus increasing the time for the development of visible lesions).
There also seems to be an increase in the odds of SLH detection over time, such that for each year that passes from the beginning of the study period, the odds of being a SLH case increase by 1.1 on average. On the surface this seems to conflict with an earlier study [13], but Risk factors associated with slaughterhouse surveillance of bTB in Great Britain it is worth bearing in mind that in [13] the authors were modelling the probability of confirmation, rather than the probability of detection. Our results suggest that there are additional year-on-year changes in risk that are not being captured by the other risk factors, which is in line with the general trend observed in the descriptive analysis (shown in Figs 1 and 2). There is also some impact of whether or not an animal spends time in a herd with a history of bTB breakdowns in the previous 5 years, such that an animal that has spent some of its life in a herd with a recent history of bTB has twice the odds of being a SLH case than an animal that has spent no time in a herd with a recent history of bTB. This is likely to be because the history of bTB in a herd captures proxy risks for individual animals within those herds, such as being more likely to harbour hidden infections [9] or having localised reservoirs of infection. This is qualitatively similar to the earlier study in the ROI [10], where the risk of submission and confirmation (in general) was smaller in animals that were sourced from herds that had been clear of bTB infection for longer. In NI [12] the risk of disclosure seems to be smaller in herds that have recently cleared infection, but then rises again for herds that have been clear for longer. This may reflect differences in the underlying epidemiology, or effectiveness of testing, in NI compared to GB.
There is some marginal evidence of a seasonal effect. The lowest odds of SLH detection are found in April-June, with a slightly higher odds of SLH detection for animals slaughtered in October-December (OR of 1.2 on average), with a relatively smooth transition in the interim quarters. This seasonal variation in detection does not necessarily represent variability in the efficiencies of slaughterhouses. We think it is more likely to be driven by a changing pattern of presentation of lesioned animals as a consequence of seasonality in patterns of slaughter and testing, and note that the pattern is similar in ROI [10] but weaker in NI [12].
Interestingly, the more times that animals have been present on herds when the herd has been tested slightly increases the odds of subsequent detection as a slaughterhouse case, but the effect is not strong (OR of 1.2 on average, relative to an animals that have never been present for a test). Animals that remain on the same herd for the whole of their lives have a slightly higher odds of SLH detection (OR of 1.2 on average, relative to animals that have moved herds at least once). This contrasts with NI [12] where homebred animals (e.g. no movements) were at slightly lower risk of disclosure, but aligns with ROI [10] where the reverse was true. We note in all cases that the ORs were not strong, and in our case this pattern may be a consequence of the introduction of pre-and post-movement testing towards the end of our study period (2007 onwards). Finally, sex seems to play a negligible role, with female animals having slightly higher odds of being a SLH case than males (with an OR of 1.1 on average, but we note that the 95% credible interval contains 1). This is in line with ROI [10] (which found a negligible effect), and NI [12] (which found a stronger positive effect for females).
Over-and-above these risk factors there seems to be an effect of breed (Fig C S1 File). However these breed effects are much smaller than the corresponding SLH effects (see below), and were effectively treated as nuisance variables in these analyses. SLH variability. Due to the large number of terms, the caterpillar plots for the SLH effects (Fig B S1 File) are hard to interpret. A key point is that the relative importance of each SLH is going to be related to both the posterior mean (how much each SLH deviates from the overall mean), and the posterior variance (how much uncertainty exists), and so in addition we produce a scatterplot of the posterior variance against the posterior mean for each SLH term (Fig  4). To aid interpretation, the point size is inversely related to the variance (so larger points have smaller variances), and the colours are mapped to the mean, such that blue points show an increased odds of picking up SLH cases, and the red points show a decreased odds (relative to the overall mean). To protect confidentiality, the spatial plot is aggregated at the county level.
The key SLHs of interest are those with large (absolute) means and small variances (i.e. those that deviate away from the mean with low uncertainty). It is important to remember that these values relate to deviations away from the overall mean as defined by the regression model, and we cannot assign a causal reason for these patterns. Rather, the model can be used to identify 'atypical' patterns across the slaughterhouses, which could be subjected to further investigation. The variables used in the regression were chosen to try to reflect different sources of risk for individual animals over their lifetime, but there is clearly some residual risk due to variables that we have not measured, which is captured in these SLH terms. To assess how much of the overall risk can be attributed to these unobserved effects, we refitted the model and excluded the SLH terms. This resulted in a posterior predictive mean estimate for the AUC of 0.85 compared to the posterior mean AUC rose of 0.88 for the full model. This suggests that the main explanatory variables chosen for the model do a good job of capturing most of the risk, but that there is some additional heterogeneity between slaughterhouses.
One possibility is that these SLH effects could result from variability in the underlying disease prevalence of animals sent to particular slaughterhouses. For example, those slaughterhouses sourcing animals from high-risk areas are more likely to find 'hidden' infections than those sourcing from low-risk areas. Plotting these estimates according to the point location of each slaughterhouse suggests that there is some spatial correlation, such that on average slaughterhouses situated in high-risk areas (south-west England and Wales) find more SLH cases than those situated in low-risk areas (Fig 4), despite those areas testing more frequently. (We note that the spatial location of a slaughterhouse is not necessarily indicative of where the animals are sourced from. It is possible to send animals to slaughterhouses that are located some distance from the farm, but generally farmers tend to send animals to SLHs that are situated locally). Another possibility is that there are differences in the detection rates within slaughterhouses due to factors such as training or throughputs. Of particular interest are  -hand panel). The larger the point, the more precisely the term has been estimated, with blue points corresponding to slaughterhouses that find more cases than expected, and red points corresponding to slaughterhouses that find less cases than expected, according to the fitted model. slaughterhouses that have a lower-than-average detection rate situated in high-risk areas, and those that have a higher-than-average detection rate situated in low-risk areas.
An important interpretational consideration is that the form of the model itself defines the structure of these SLH terms, such that we would always expect deviations around the true mean. (We model these terms using a hierarchical normal distribution with a shared variance term (on the logistic scale), so by definition we would expect spread around the overall mean).
The key point is the magnitude of these differences. The posterior mean odds ratios range from between 0.11 to 10 relative to the overall mean, which represents an epidemiologically important level of variability between (a small number) of slaughterhouses. (If there were small differences we would still see a spread around the mean, but the range would be much smaller-as is seen in the breed effects, which range between 0.65 and 1.7 on the same scale).

Conclusions
We describe an increase in both the absolute number and proportion of breakdowns disclosed through slaughterhouse surveillance in Great Britain between 1998 and 2013. This increase remains even when accounting for changes in confirmation rates and changes in the proportion of herds subject to different testing frequencies.
We used a case-control study to explore the impact of individual animal-level factors on the risk of an animal being identified as a slaughterhouse case. We showed that after adjusting for other sources of risk there are residual between-slaughterhouse variabilities in detection rates. There is also an additional increased risk over time of animals being detected as slaughterhouse cases that is not explained by the available risk factors. This variability is not necessarily due to differences in practice within individual slaughterhouses, but may manifest due to systematic differences in the risk of infection or efficiency of tuberculin testing within the input population for each slaughterhouse. As such, the key use for this model is to identify 'atypical' slaughterhouses (i.e. ones that exhibit in inflated or deflated detection rates relative to the overall mean predicted by the model). Hence, although our model cannot assign causal reasons for these differences, it will act as a guide for further investigations. Currently the Animal and Plant Health Agency have adopted the model to be used as part of their ongoing bTB surveillance programme.
Supporting information S1 File. Supplementary information. Includes all code necessary to rerun analyses, plus additional plots and tables not included in the main text. (PDF) S1 Data. A data file (.rds) containing the raw data.