Assessing the effects of the first 2 years of industry‐led badger culling in England on the incidence of bovine tuberculosis in cattle in 2013–2015

Abstract Culling badgers to control the transmission of bovine tuberculosis (TB) between this wildlife reservoir and cattle has been widely debated. Industry‐led culling began in Somerset and Gloucestershire between August and November 2013 to reduce local badger populations. Industry‐led culling is not designed to be a randomized and controlled trial of the impact of culling on cattle incidence. Nevertheless, it is important to monitor the effects of the culling and, taking the study limitations into account, perform a cautious evaluation of the impacts. A standardized method for selecting areas matched to culling areas in factors found to affect cattle TB risk has been developed to evaluate the impact of badger culling on cattle TB incidence. The association between cattle TB incidence and badger culling in the first 2 years has been assessed. Descriptive analyses without controlling for confounding showed no association between culling and TB incidence for Somerset, or for either of the buffer areas for the first 2 years since culling began. A weak association was observed in Gloucestershire for Year 1 only. Multivariable analysis adjusting for confounding factors showed that reductions in TB incidence were associated with culling in the first 2 years in both the Somerset and Gloucestershire intervention areas when compared to areas with no culling (incidence rate ratio (IRR): 0.79, 95% CI: 0.72–0.87, p < .001 and IRR: 0.42, 95% CI: 0.34–0.51, p < .001, respectively). An increase in incidence was associated with culling in the 2‐km buffer surrounding the Somerset intervention area (IRR: 1.38, 95% CI: 1.09–1.75, p = .008), but not in Gloucestershire (IRR: 0.91, 95% CI: 0.77–1.07, p = .243). As only 2 intervention areas with 2 years of data are available for analysis, and the biological cause–effect relationship behind the statistical associations is difficult to determine, it would be unwise to use these findings to develop generalizable inferences about the effectiveness of the policy at present.


Bovine tuberculosis (TB) caused by the bacterium Mycobacterium bovis
is an important problem for the British cattle industry. Existing controls include testing and slaughter of test-positive cattle, with herd test frequency determined by local incidence, along with surveillance of all slaughtered cattle in the abattoir. But despite these controls, the incidence of TB in cattle in England has shown a generally increasing trend over the last thirty years. In 2015, around 5% of cattle herds were under movement restriction due to a TB incident in England, and in the high risk area of England, this rose to around 11% (Harris et al., 2017).
European badgers (Meles meles) are a host species for M. bovis and represent a wildlife reservoir of infection. The use of culling to control the transmission of M. bovis between this wildlife reservoir and cattle has been widely debated. There is evidence from the Randomized Badger Culling Trial (RBCT) that proactive (systematic and widespread) culling of badgers for at least 4 years could reduce the incidence of confirmed TB in cattle by 23.2% (95% CI: 12.4%-32.7%) (Donnelly et al., , 2007. The RBCT was conducted in England between 1998 and. The net effect of proactive culling per year in the RBCT was initially detrimental as a 24.5% increase in incidence (95% CI: 0.6% lower to 56.0% higher) was observed in a 2-km-wide buffer around the culled area (Donnelly et al., 2007). Detrimental effects were attributed to the disruption of badger social structures (perturbation) affecting contact rates between cattle and badgers . However, an overall benefit was observed after the third year of culling and subsequent culls (Donnelly et al., 2007).
The first round of industry-led culling took place in two areas, one in west Somerset and the other in west Gloucestershire between August and November 2013 (hereafter referred to as the Somerset intervention area and Gloucestershire intervention area). Culling licenses were issued for the two areas by Natural England (an executive nondepartmental public body advising Defra on the natural environment) under the Protection of Badgers Act 1992 to enable groups of farmers and landowners to reduce local badger populations. Licenses were subject to a number of criteria. These included the application area to be at least 150 km 2 , at least 70% of the land to be accessible for culling, cattle herds to be subject to annual TB testing and "reasonable biosecurity" to be in place (Defra 2015). It was also a requirement that the culling should plan to reduce the estimated badger population by 70% and be conducted for a minimum of 4 years (Defra, 2015). The aim of the initial culls was to assess the practicalities and impact of the intervention on badger population density. Using a combination of cage trapping and controlled shooting, 1,296 badgers were culled in the Somerset intervention area in Years 1 and 2 combined and 1,198 were culled in the Gloucestershire intervention area in the 2 years combined. The target minimum number that was estimated to result in a 70% reduction in the badger population was not achieved in either area in Year 1 and was achieved in the Somerset intervention area only in Year 2, although there is uncertainty around the calculation of these minimum numbers due to difficulties in estimating the true badger population size (AHVLA 2014, Defra 2014c,d). The third year of culling took place in both areas during autumn 2015. A further area was licensed in Dorset in 2015, and seven new areas were licensed in 2016. Data for analyzing cattle TB incidence in the year following the culls in 2015 were not available for this analysis.
Our aim was to measure and assess any association between the badger culling intervention and the incidence of TB in local cattle, taking account of the methodological constraints resulting from this not being a study based on sound intervention study design principles. TB incidence in cattle herds located within areas where industry-led culling is conducted ("intervention areas") is compared with TB incidence in herds in unculled areas matched on characteristics that affect cattle TB risk ("comparison areas") in a similar, but not identical, approach to that used for analyzing the impact of culling during the RBCT (Donnelly et al., 2003(Donnelly et al., , 2007. The null hypothesis being tested is that TB incidence is the same in the intervention areas and their comparison areas in the years since badger culling began. Cattle TB incidence is also monitored among herds in 2-km "buffer areas" surrounding the intervention areas and compared to incidence among herds in similarly defined areas around unculled comparison areas to monitor for potentially adverse effects on cattle TB incidence . The definition of incidence being used for these analyses is the incidence of OTF-W (Officially Tuberculosis Free status-Withdrawn) incidents. OTF-W incidents are those where confirmatory evidence of M. bovis infection has been identified in at least one bovine animal slaughtered for disease control purposes. They also include other incidents upgraded to OTF-W for epidemiological reasons and slaughterhouse case-disclosed incidents that are confirmed through culture of M. bovis.
Here we report an unadjusted analysis of OTF-W incidence rate per 100 herd years at risk for the first 2 years following industry-led culling, alongside more detailed analysis with adjustment for other factors that are likely to be associated with a risk of TB in cattle.

| Selection of comparison areas
Areas were selected from within England which met the following inclusion criteria: within a high-risk TB area (Defra 2014a); cattle herds subject to annual TB testing; having no land closer than 2 km to any intervention area boundary (at the time of selection).
Boundary information for the intervention areas was used in an automated programming procedure using ArcGIS 10.0 software (ESRI Release 10.0. Redlands, CA, USA) to generate a population of potential comparison areas with centroids shifted 5 km from one another and rotated at 45-degree intervals. Comparison areas were selected in May 2014, around 9 months after the start of the culls in late 2013 without any consideration of cattle TB incidence data after the culls began. They were matched to intervention areas on factors known to affect cattle TB incidence in an attempt to control for confounding factors that could explain between-area variation not attributable to culling and therefore bias the analysis. These are described in more detail below.

| Ranking of comparison areas for similarity to an intervention area
Each intervention area and potential comparison area was described by six factors or attributes previously associated with TB incidence that had been extracted from the APHA Sam surveillance database (Table 1). The distribution of each factor was summarized using deciles, and the absolute difference between the decile to which the intervention area belonged and the decile to which each potential comparison area belonged was calculated. A score based on the sum of the absolute differences for each of the attributes was then used to rank potential comparison areas, with areas having the smallest summed absolute difference ranked highest. No weighting was applied to the ranking attributes.

| Selection of matched comparison areas based on rank
The highest ranked comparison areas were preferentially selected using a semi-automated process that also sequentially excluded comparison areas that covered part of a higher ranking area or one of the two intervention areas. Ten comparison areas were selected for each intervention area, in order to minimize the consequences of future censoring of comparison areas due to initiation of new intervention areas that cover land in a previously selected comparison area. The process of ranking potential comparison areas and for manually selecting highest-ranking comparison areas was conducted independently by two members of the project team, and any differences in outputs were reconciled.
Other attributes likely to be associated with cattle TB incidence were summarized for each intervention and comparison area (Table 2).

| Statistical analysis
The observed OTF-W incidence rate in cattle herds was calculated for each of the two intervention areas, in the 2 km-wide unculled buffer area around each intervention area, and in the 20 comparison areas (10 per intervention area) and their 2 km-wide buffers for the first and second years following the baseline date (the date on which the first cull started) and the periods 0-12 months, 12-24 months, and 24-36 months prior to the baseline date. In each case, the herd yearsat-risk was calculated using results from whole herd tests as the sum of the time herds in each area were unrestricted and therefore at risk of new infection (a new OTF-W incident) during the period of interest (Downs et al., 2013). Crude OTF-W incidence rate ratios (IRRs) were calculated for both the central intervention areas and buffer areas in each reporting period. 95% confidence intervals were calculated, and a probability level of p < .05 was considered to be statistically significant.
For the adjusted analysis, data for explanatory variables were merged with the OTF-W incidence event data and the herd years-at-risk denominator data. Adjusted IRRs were estimated using Poisson regression models comparing intervention with comparison areas separately for the central areas (i.e., the intervention or comparison areas) and the buffer areas (the 2-km buffer areas surrounding intervention and comparison areas), controlling for factors known to be associated with TB incidence. Initial models contained a single intervention status variable which represented the badger culling intervention in both the Somerset and Gloucestershire intervention areas combined and a variable representing the geographical area (coded as Somerset or Gloucestershire) which included both the intervention areas and their matched comparison areas. The inclusion of a variable representing the interaction between the intervention status and area was also explored in these initial models.
T A B L E 1 Distribution of attributes used to rank comparison areas across intervention areas and matched comparison areas (as mean values across 10 comparison areas) Two models were built, to estimate the effect of intervention on OTF-W incidence in the central areas (Model A) and the buffer areas (Model B) in the first 2 years of culling. Each model included two separate intervention variables, one for each of the Somerset and Gloucestershire intervention areas in order estimate the independent effects of the intervention in each geographical area. Prior to building the models, bivariable associations between the outcome and explanatory variables were assessed, and pairwise correlations were explored to investigate any colinearity between explanatory variables. Explanatory variables considered to be important potential confounders a priori were included in the initial models. These consisted of the variables used to match the comparison areas-area (Somerset or Gloucestershire), historic (3-year) incidence, median herd size, the proportion of herds that were dairy, distance to the matched intervention area and the proportion of the land in the area subjected to proactive culling during the RBCT-to account for imperfect matching, as well as estimated badger sett density (Judge, Wilson, Macarthur, Delahay, & McDonald, 2014). Additional factors were then added to the models one at a time to assess the effect they had on the estimated impact of the intervention. Akaike's Information Criterion (AIC) was used to select models with good fit and parsimonious use of explanatory variables (as recommended by Burnham & Anderson, 2002). information. To account for any overdispersion (greater-than-expected variance between areas due to clustering of cattle TB incidents), the "robust" or "sandwich" estimator of variance (Huber, 1967;White, 1982) was used to ensure that confidence intervals were adequately adjusted. Overdispersion in the final models was then assessed using the deviance and Pearson statistics.
A sensitivity analysis of the influence of individual comparison areas on the estimated associations between industry-led culling and OTF-W incidence was conducted by systematically removing one area at a time from the analysis. A further sensitivity analysis was conducted by systematically removing one explanatory variable at a time from the models. The aim of this was to identify which of the explanatory variables most influenced the estimated associations between industry-led culling and OTF-W incidence.
As part of our validation of the final models, three further sets of models were developed. The first set examined OTF-W incidence in the first 2 years since culling began, but using a reduced number of explanatory variables (area, time at risk, and historic incidence) in an approach more similar to the analyses of RBCT data (Donnelly et al., 2003. The second set examined OTF-W incidence in the year prior to the culls to explore whether there was any association T A B L E 2 Distribution of explanatory factors across intervention areas and matched comparison areas (as mean values across 10 comparison areas) between the intervention areas and incidence prior to the start of culling; and the third set estimated effects for the individual years since culling began. The results of these validation analyses are presented in the Appendices.

| Comparison of characteristics in intervention areas and matched comparison areas
The characteristics of the intervention areas and buffers and their matched comparison areas are presented in Table 1. There were 6,658 and 6,387 potential comparison areas generated through the 5-km shifts, respectively, for the Somerset and Gloucestershire intervention areas. The 10 selected comparison areas were fairly well matched to their respective intervention areas when values for the ranking factors were averaged (Table 1). The proportion of all incidents that were OTF-W was over 80% for all compari-  Table 2).
Locations of Gloucestershire Intervention Area and 10 matched comparison areas across map of smoothed cattle TB herd incidence in the high incidence area of England, 2013. Culling and comparison areas are symbolized by a solid circle that approximates actual size. OTF-W incident density was created using the spatial analyst kernel density tool within ArcGIS 10.0. (b) Locations of Somerset Intervention Area and ten matched comparison areas across map of smoothed cattle TB herd incidence in the high incidence area of England, 2013. Culling and comparison areas are symbolized by a solid circle that approximates actual size. OTF-W incident density was created using the spatial analyst kernel density tool within ArcGIS 10.0

| Analysis of the impact of the intervention on OTF-W incidence
The observed OTF-W incidence rates per 100 herd years at risk for each of the central intervention and comparison areas, and their respective buffer areas for Year 1, Year 2, and for each of the 3 years prior to the cull are presented in Table 3. Descriptive analyses without adjustment for confounding factors showed no statistically significant difference between OTF-W incidence rates in the Somerset central intervention area and comparison areas for Year 1 or Year 2. A weak difference in incidence rate between the Gloucestershire central intervention and comparison areas was observed for Year 1 (IRR: 0.64, p = .05), but not for Year 2. No differences in incidence rates were observed between the intervention buffer areas and comparison buffer areas.
In the initial models for the multivariable analysis for the first 2 years of culling which contained a single intervention variable and an area variable, the inclusion of an interaction between the intervention and area generated a p-value of <.001 for both the central and buffer areas.
This indicated that the effects of the interventions in the two central areas should be examined separately by including area-specific intervention variables. The final models comparing OTF-W incidence in the first 2 years of culling, adjusted for explanatory variables, for the central areas (Model A), and the buffer areas (Model B) are presented in Table 4.

For the central areas (Model A), the number of badgers removed
historically and the proportion of farms with more than one fragment of land in the area were included in the model alongside the a priori explanatory variables as they affected the association between the intervention and OTF-W incidence in Somerset and reduced the AIC (Table 4). A lower IRR for the intervention was observed in the Somerset and Gloucestershire central areas (both p < .001), with a considerably lower IRR in Gloucestershire (IRR: 0.42; 95% CI: 0.34, 0.51) than in Somerset (IRR: 0.79; 95% CI: 0.72, 0.87). As expected, increases in time at risk and historic incidence were significantly associated with increased OTF-W incidence. Associations with increased OTF-W incidence were also observed for the proportion of herds that were dairy, the total number of badgers removed historically, and the proportion of farms with more than one fragment of land in the area.
Increases in the distance to the intervention area and the estimated number of badger setts per 100 km 2 were both associated with a decrease in OTF-W incidence in the first 2 years of culling. The Deviance and Pearson goodness-of-fit tests for overdispersion (p = .898 and p = .904, respectively) indicated no evidence of there being significantly more variation in the data than expected.
The sensitivity analysis of OTF-W incidence in the central areas showed that with the removal of one of the Somerset comparison areas (WS02), there was no longer an association between the intervention and OTF-W incidence in Somerset (IRR = 1.00; 95% CI: 0.80, 1.26; p = .976). Further investigation of this comparison area could not identify any major differences in the data available that might explain this effect, and this area was retained in the central area model.
The sensitivity analysis of explanatory variables demonstrated that distance to intervention had a strong influence on the association T A B L E 3 OTF-W incidence rates per 100 herd years at risk and unadjusted incidence rate ratios (IRRs) for central and buffer intervention areas versus the central and buffer comparison areas, respectively, in Somerset and Gloucestershire, for each 12-month reporting period between the intervention and OTF-W incidence in the central area of Somerset (Table A1 in Appendix 1). When the distance-to-intervention variable was removed from the model, the IRR for the association between the intervention and OTF-W incidence in Somerset increased and was no longer statistically significant (p = .476). This indicated that after adjusting for the other explanatory variables, inclusion of the distance between the intervention and comparison areas variable had an influential effect on the estimated IRR.

| Analysis of the impact of the intervention on OTF-W incidence in surrounding areas
For the buffer areas (Model B), the proportion of the area classed as urban and the proportion of farms with more than one fragment of land in the area were included alongside the a priori explanatory variables as they affected the association between the intervention and OTF-W incidence in Somerset and reduced the AIC (Table 4). An increased estimated IRR was observed for the buffer area around the intervention area in Somerset (IRR = 1.38; 95% CI: 1.09, 1.75; p = .008) while no significant association between the intervention and OTF-W incidence was observed in the Gloucestershire buffer, having adjusted for explanatory variables (IRR = 0.91; 95% CI: 0.77, 1.07; p = .243).
Increases in time at risk and historic incidence were associated with an increased OTF-W incidence in the first 2 years of culling. Associations with an increase in OTF-W incidence in the intervention areas were also observed for the proportion of herds that were dairy, the proportion of land classed as urban, and the proportion of farms with more than one fragment of land in the area. An increase in the proportion Incidence rate ratio. c Intervention is industry-led badger culling, as described in the Section 1. of RBCT proactive land was associated with a decreased OTF-W incidence in the first 2 years of culling. The Deviance and Pearson goodness-of-fit tests for overdispersion (p = .762 and p = .774, respectively) indicated no evidence of there being significantly more variation in the data than expected.
The sensitivity analysis of OTF-W incidence in the buffer areas showed that removing Somerset comparison buffer areas WS02B or WS07B affected the estimated association between OTF-W incidence in the Somerset buffer and comparison area buffers (when WS02B removed: IRR = 1.64, p = .181; when WS07B removed: IRR = 1.24, p = .125). The removal of WS07B reduced the estimated IRR for intervention in the Gloucestershire buffer area to 0.81 (p = .011). However, further examination of these areas could not identify any major differences in the data available that might explain these effects and the comparison areas were retained in the final model.
The sensitivity analysis of explanatory variables demonstrated that the proportion of the herds that were classed as dairy and the proportion of farms with more than one fragment of land in the area appeared to have a strong influence on the association between the intervention and OTF-W incidence in both intervention buffer areas (Table A2 in Appendix 1). For both areas, the inclusion of these variables increased the estimated IRR. In the Somerset intervention buffer area, this resulted in a significant association between the intervention and OTF-W incidence. In the Gloucestershire intervention buffer area, adjusting for these factors removed the observed association between the intervention and OTF-W incidence. Other factors that had an influential effect on the IRR in the Somerset intervention buffer area were median herd size, and the proportion of land that was exposed to proactive culling in the RBCT.

| DISCUSSION
Industry-led badger culling in England is a component of the current multifaceted bovine TB control policy, not a scientific trial, and as such lacks randomization and other controls. Nonetheless, it is important to monitor the effects of the industry-led culling and evaluate the impacts taking study limitations into account. The external validity (i.e., the generalizability) of the findings will be limited by the way in which the study areas were selected, and the extent to which the circumstances under which the policy is operating are similar to circumstances where culling is proposed in the future.

| Suitability of selected comparison areas
In order to assess the impact of the intervention on cattle TB incidence, we first had to develop a method for selecting suitable comparison areas. For the RBCT, areas were prospectively selected in regions similar in terms of location and geographical characteristics and where the incidence of cattle TB in England was highest (Independent Scientific Group on Cattle TB 2007). RBCT areas were randomized to receive proactive or reactive culling or no culling (survey only). As randomization was not possible with the current policy, comparison areas were matched to purposively selected intervention areas on factors associated with cattle TB incidence including distance between areas and historic OTF-W incidence. Herd size was used as a criterion because it has been consistently and positively associated with TB in cattle herds Goodchild & Clifton-Hadley, 2001;Ramirez-Villaescusa, Medley, Mason, & Green, 2010). The proportion of land previously in an RBCT proactive culling area was also used to match areas, whereas land previously part of an RBCT reactive area was not. This is because no long-term effect from reactive culling on recurrence of cattle TB incidents has been detected . Reductions in cattle TB incidents in proactively culled RBCT trial areas were measurable during culling and after the cessation of culling (Donnelly et al., 2007;Jenkins, Woodroffe, & Donnelly, 2010). Further exploratory analyses are consistent with an ongoing, but diminishing (test for temporal trend p = .008), benefit of proactive culling at 55-60 months post-trial (Donnelly, Jenkins, & Woodroffe, 2011).
The method used to select comparison areas maximized the probability of identifying areas that were very similar to the intervention areas in terms of risk factors for TB, thereby reducing the impact of confounding factors in descriptive comparisons of rates. Nevertheless, the matching was not perfect and matching variables were incorporated into the regression models as a priori confounders in an attempt to account for the imperfect matching.

| Impact of the intervention on OTF-W incidence
In annually published surveillance reports describing the cull areas (APHA 2015, 2016), observed OTF-W incidence was reported for the intervention areas combined to maximize the sample size. In the analyses reported here, a difference in effects was observed between the Somerset intervention area and the Gloucestershire intervention area, also demonstrated in a statistically significant area interaction term in the modeling, so results have been reported separately for the individual areas.
When the effect of intervention in each area was analyzed with adjustment for potentially confounding factors, significantly lower IRRs for intervention were revealed in the Gloucestershire central area in Year 1 (Table A5), and in each of the central areas in Year 2 (Table A6) and when both Year 1 and Year 2 were analyzed together (Table A4).
Additional modeling to better understand interarea comparison indicated that OTF-W incidence was lower in the Gloucestershire intervention area than its comparison areas before culling began (Table S4), as was also observed in the unadjusted analysis. Similar analyses of TB incidence prior to culling were not performed in the RBCT because the randomized selection of areas in that trial was expected to have removed any biases related to differences between areas.

| Impact of the intervention on OTF-W incidence in surrounding areas
The multivariable analysis revealed an increased IRR for the intervention in the Somerset buffer area (IRR: 1.38, p = .008) and no effect in the Gloucestershire buffer area (IRR: 0.91, p = .243). The observed effect in the Somerset buffer area is consistent with the results of the RBCT, which indicated that an increase in cattle TB incidence could be expected in the buffer areas due to perturbation of the badger population Woodroffe et al., 2006). When the individual years were investigated (Tables A5 and A6 in Appendix 2), the increased IRR for intervention in the Somerset buffer area was only observed in Year 1, and a lower IRR for intervention in the Gloucestershire buffer area was observed in Year 1 and Year 2 individually. None of the differences in the buffer areas were observed when the year prior to the culls was analyzed (Table A4).

| Reliability of the models
The associations between the intervention and OTF-W incidence in Somerset that were observed in the models including all associated covariates (Table 4) were not observed in the simpler RBCT-like models which included a small number of covariates (Table A3). The modeling performed for the RBCT did not include additional explanatory factors as differences in these potential confounders between areas were assumed to have been accounted for by the randomization. Indeed, there was no substantial over dispersion despite the limited number of explanatory variables included (Donnelly et al., 2003. In this study, where randomization was not possible, simpler models that do not control for potential confounding factors may be less effective in detecting an association between industry-led badger culling and OTF-W incidence. The sensitivity analysis performed on the buffer area model identified two Somerset comparison buffer areas as being particularly influential. There was more variation in incidence among the Somerset buffer comparison areas, which might explain why some areas are more influential than others in the models. The greater variation in incidence between buffer areas is not surprising as they are geographically smaller than the central areas and contain fewer herds. The positive influence of the proportion dairy variable and the proportion of farms with more than one fragment of land on the association between culling and OTF-W incidence in both intervention buffer areas is not unexpected as both are factors that have been found to be positively correlated with incidence in previous studies (Broughan et al., 2016;Goodchild & Clifton-Hadley, 2001;Porphyre, Stevenson, & McKenzie, 2008). The sensitivity analysis demonstrates that there is considerable uncertainty around the estimates that cannot currently be explained.
Our analysis will be subject to confounding in common with all observational studies despite attempts to control for these effects. Further follow-up and inclusion of data from additional culled areas will likely make future analyses more informative (Donnelly, Bento, Goodchild, & Downs, 2015).

| OTF-W incidents versus all TB incidents
OTF-W incidents were used as the outcome measure rather than all TB incidents because the RBCT only showed an association between OTF-W incidence and culling (Independent Scientific Group on Cattle TB, 2007). This reasoning assumes that the current case definition of OTF-W herds is comparable to the definition of a confirmed incident used during the RBCT. This may not be the case due to operational and policy changes in the management of TB incidents, particularly since the introduction of the TB eradication strategy in 2014 (Defra 2014b). It is not understood why the effects of culling observed in the RBCT were only observed for OTF-W-like incidents. To investigate further, we extended our analysis of industry-led culling to include all TB incidents (Table A7 in Appendix 3). Qualitatively similar results were obtained to those obtained when modeling only OTF-W incidents, although the estimated effects were slightly weaker. This might suggest that culling has more influence on OTF-W incidents than those without postmortem evidence of infection, suggesting that there is a stronger association between badgers and the occurrence of OTF-W incidents compared with all TB incidents. The probability of M. bovis infection in a herd will be higher for OTF-W incidents than for TB incidents that have not been confirmed by postmortem evidence of infection which could possibly increase the power to detect effects (De la Rua-Domenech et al., 2006).

| Difficulties with assessing the impact of industry-led culls
In Year 1, the target of reducing the estimated precull population by at least 70% was not achieved, nor was spatial coverage of the culling homogeneous (Independent Expert Panel 2014). In Year 2, the minimum number of badgers to be removed was achieved in the Somerset intervention area, but not in the Gloucestershire area (Defra, 2014c). It has been estimated that in order to have sufficient power to be confident of observing significant differences in the incidence of OTF-W herd incidents, matched intervention and comparison areas will need to be observed for at least 3 years after culling begins, and that this increases to 4 years if only two intervention areas are licensed (Donnelly et al., 2015). As such, it was not expected that significant differences would be observed in the first 2 years of follow-up, particularly since the proportion of badgers removed was estimated to be lower than in the RBCT. While it is possible that the significant differences observed could be true differences, selection bias cannot be ruled out due to the nonrandomized selection of the cull areas.
It is also possible (although unlikely given the low p-values) that the observed effects are due to chance. Caution should always be applied in the interpretation of models which attempt to fit a relatively large number of parameters (Babyak, 2004;Green 1991), particularly when the dataset is small. However, these findings are coherent with findings from the RBCT.
With the limited data available to date, it is possible that spurious associations could have arisen. For example, the central area model for the first 2 years since culling began identified estimated badger sett density as being negatively associated with OTF-W incidence in the area. This counterintuitive result at first sight casts doubt on the model, but an increase in estimated badger sett density could increase the likelihood that the transmission of M. bovis to cattle is from badgers, and thus, culling badgers may have more of an impact on cattle TB incidence. The impact of the efficacy of the culls could not be investigated in this analysis because of the uncertainty around the estimates of badger population reduction (AHVLA 2014).
A number of potential confounding factors have been included in this analysis. However, even with the inclusion of this additional information, it is unlikely that all sources of bias have been estimated or detected. For example, we do not have all the information that may have been considered by groups of farmers and land owners selecting the culling areas. Factors that affect cattle TB risk may also change with time. Furthermore, farmers and landowners aware of the intervention may change their behaviors in ways that affect cattle TB incidence independent of any changes in Government policy (Gale, 2004).

| CONCLUSIONS
We have identified reductions in OTF-W incidence in both the Somerset and Gloucestershire intervention areas and an increase in OTF-W incidence in the Somerset buffer area in the first 2 years of industry-led culls. As this analysis has been performed on only two intervention areas with only 2 years of follow-up data, and the biological cause-effect relationship behind the statistical associations observed is difficult to determine, it would be unwise to use the findings of this analysis to develop generalizable inferences about the effectiveness of the policy at present. These findings are consistent with findings from the RBCT, although a time lag of around 4 years was observed between culling in the RBCT and measureable significant effects on cattle incidence (Donnelly et al., 2007).
Although a trial randomizing culling to different matched areas would be the most rigorous design for the evaluation of the effect of a badger culling policy, this type of design is not possible when the areas where culling is conducted are selected by stakeholders. The long-term value of monitoring information from industry-led culling will depend on the conduct of the culls, the number of areas eventually licensed and the extent to which other bovine TB control policies remain stable. Continued delivery of the intervention in these areas, and further roll out of the intervention to other areas may enable better assessments to be made of the impact of industry-led culling on TB incidence in cattle.

APPENDIX 1 Sensitivity analysis of individual explanatory variables
T A B L E A 1 Model estimates for intervention (industry-led badger culling) in the Somerset and Gloucestershire central areas in the first 2 years since culling began when individual explanatory factors are either added (in yellow) or removed (in blue) from three models: a simple RBCT-like model, the baseline model which was the RBCT-like model plus factors considered important a priori, and the final model which was the baseline model plus additional explanatory factors. Where there are notable changes in associations from the starting models, the p values are highlighted in bold font. (Final model reported in Table 4  T A B L E A 2 Model estimates for intervention (industry-led badger culling) in the Somerset and Gloucestershire buffer areas in the first 2 years since culling began when individual explanatory factors are either added (in yellow) or removed (in blue) from three models: a simple RBCT-like model, the baseline model which was the RBCT-like model plus factors considered important a priori, and the final model which was the baseline model plus additional explanatory factors. Where there are notable changes in associations from the starting models, the p values are highlighted in bold font. (Final model reported in Table 4

Replicating the RBCT models
In order to compare the current results with those of the RBCT, we assessed the effect of including just those factors that had been identified as important in the RBCT modeling when analyzing the first 2 years of culling (Table A3) T A B L E A 3 Multivariable Poisson regression models describing the effect of intervention separated by area on the number of OTF-W a incidents in the first 2 years since culling began, adjusted for a reduced set of explanatory variables as applied in analysis of the RBCT (Donnelly et al., 2003

Statistical analysis of the year prior to the culls
In the initial models for the year prior to the culls which contained a single intervention variable and an area variable, the inclusion of an interaction between the intervention and area generated a p-value of <.001 for the central areas and p = .775 for the buffer areas, so the effects were examined separately for each area.
For the central areas (Table A4,  For the buffer areas (Table A4, Model 2), the length of motorway through an area was included in the model alongside the a priori explanatory variables as it modified the coefficient for the association between the intervention and OTF-W incidence in the Gloucestershire buffer area. There appeared to be no effect of intervention on OTF-W incidence in the buffer areas, having adjusted for explanatory variables. Increases in time at risk, historic two-year incidence rate and median herd size were significantly associated with an increase in the numbers of OTF-W incidents in the year prior to the culls. An increase in the estimated number of badger setts per 100 km 2 was associated with a decrease in OTF-W incidence in the year prior to the culls. The Deviance and Pearson goodness-of-fit tests for overdispersion generated large p values (p = .592 and p = .593, respectively) indicating that there was no evidence of there being significantly more variation in the data than expected.

Statistical analysis of each of the first two individual years after culling began
In addition to examining the association between the intervention and OTF-W incidence in the 2 years of culling, we also examined the as- T A B L E A 6 Multivariable Poisson regression models describing the effect of intervention (badger culling) stratified by area on the number of OTF-W a incidents in the second year since culling began, adjusted for explanatory variables. Model 1 describes the effect in the central intervention areas compared with matched central comparison areas, and Model 2 describes the effect in the 2-km buffer area for the intervention areas compared with buffer areas for matched comparison areas Incidence rate ratio. b Intervention is industry-led badger culling, as described in the Section 1.