Statistical models examining the effect of proactive badger culling
The most important result from the 2006 paper is presented as Table 1.
Table 1
Log-linear regression model fitted to ‘confirmed’ incidents (breakdowns) in whole trial areas (inner and outer regions combined), since the initial proactive cull, based on VETNET location data.
Parameter
|
|
Estimate
|
Overdispersion-adjusted SE
|
Overdispersion-adjusted p-value
|
Intercept
|
|
-0.362
|
0.932
|
0.70
|
Treatment
|
Proactive
|
-0.207
|
0.073
|
0.005
|
|
Survey-only
|
-
|
-
|
-
|
Triplet
|
A
|
-0.269
|
0.234
|
|
|
B
|
0.156
|
0.166
|
|
|
C
|
0.432
|
0.161
|
|
|
D
|
-0.304
|
0.190
|
|
|
E
|
-0.026
|
0.175
|
|
|
F
|
-0.159
|
0.186
|
|
|
G
|
0.482
|
0.195
|
|
|
H
|
-0.248
|
0.189
|
|
|
I
|
-0.502
|
0.197
|
|
|
J
|
-
|
-
|
-
|
Log of baseline herds
|
|
0.047
|
0.248
|
0.85
|
Log historic incidence*
|
|
1.241
|
0.213
|
< 0.001
|
* Historic incidence was calculated over a three-year period and was restricted to ‘confirmed’ incidents.
This is a reproduction of Table 8. from the supplementary information provided in the 2006 paper
Here the parameter treatment is highly significant and is interpreted as a reduction in the herd incidence of 18.7% in areas where badgers had been culled. What is a little unusual in this analysis is that the parameter Log of baseline herds was not fixed at 1. The rate is normally the count divided by some unit of exposure, in this case the exposure would be the number of baseline herds.
If exposure was the same across all triplets, then the log linear Poisson regression would be
$$log{\mu }_{x}={\beta }_{0}+{\beta }_{1}x+...$$
To control for baseline Bx as there is varying levels of exposure, this then becomes:
$$log\frac{{\mu }_{x}}{{B}_{x}}={\beta }_{0}+{\beta }_{1}x+...$$
=>
$$log{\mu }_{x}={\beta }_{0}+{\beta }_{1}x+{logB}_{x}+...$$
The usual practice in Poisson regression with a variable exposure across group would be to fix the parameters for logBx at 1, and code it as an offset variable. This then would assume that the incidence in each triplet area is directly proportional to the number of baseline herds, that is if the number of baseline herds is doubled, there is a doubling of the incidence when other parameters are held constant. However, in the 2006 paper, the parameter for logBx was not fixed at 1, but was analysed as an explanatory variable and hence estimated by the model. This could be a reasonable approach if there was some a priori reason why the number of herd breakdowns was not linearly associated with the number of herds in the baseline. All areas in each triplet had an area of approximately 100 km2. But with a relatively high variation in the number of herds in each area, there was considerable variation in herd density. Standard theory on infectious diseases does indicate that transmission is dependent on density (ie contact rate)18 so this might be a justification for not fixing the logBx parameter at 1. Models suggest that disease transmission increases as the density or numbers of cattle and badgers increases19. Hence with a fixed area for each triplet the parameter of log(baseline) or log(herd years at risk) might be expected to be above 1. However, the parameter value that is estimated is close to and not significantly different from zero (Table 1). This is not consistent with a density dependent increase in transmission. Furthermore, the interpretation for this is that the number of bTB herd incidents is independent of the number of herds in an area within a triplet. That is if the number of baseline herds is doubled, the incidence does not change. This does not seem to be credible.
Alternatively the logBx parameter could be fixed at 1 as an offset term. Doing this, and keeping the same variables in the regression and correcting for overdispersion (as reported in the 2006 paper), the result reported in Table 2 is obtained. This is model 2:
Table 2
Log-linear regression model fitted to ‘confirmed’ incidents (breakdowns) in whole trial areas (inner and outer regions combined), since the initial proactive cull, based on VETNET location data.
Parameter
|
|
Estimate
|
Overdispersion-adjusted SE
|
Overdispersion-adjusted p-value
|
Intercept
|
|
-3.408
|
0.824
|
0.003
|
Treatment
|
Proactive
|
-0.146
|
0.117
|
0.246
|
|
Survey-only
|
-
|
-
|
-
|
Triplet
|
A
|
0.328
|
0.285
|
|
|
B
|
0.214
|
0.271
|
|
|
C
|
0.244
|
0.252
|
|
|
D
|
-0.004
|
0.284
|
|
|
E
|
0.205
|
0.268
|
|
|
F
|
-0.396
|
0.283
|
|
|
G
|
0.000
|
0.251
|
|
|
H
|
-0.048
|
0.292
|
|
|
I
|
-0.277
|
0.309
|
|
|
J
|
-
|
-
|
-
|
Log of baseline herds
|
|
1*
|
|
|
Log historic incidence
|
|
0.741
|
0.263
|
0.023
|
*Fixed at 1. Codes as an offset variable |
Model 2 suggests that badger culling had no effect on confirmed incidents. However, the two models are not comparable by model selection techniques such as AIC or AICc as Model 1 is a Poisson model and model 2 is a quasipoisson model. Therefore, competing models were analysed using the generalized Poisson model which allowed direct comparison of over- or under-dispersed Poisson models. It is also possible to compare possible models with and without the explanatory variables presented in the 2006 paper. In addition, since there are 20 data points, but up to 14 variables (9 for Triplet, baseline herds, historic incidence, treatment and over-dispersion or variance parameter), the small sample size AIC (AICc) was used for model selection and comparison. Model 3 is the full model with log(Baseline) as an explanatory variable. With model 3, the most parsimonious model was without triplet as an explanatory variable. In that case, treatment had no effect on ‘confirmed’ incidents. Model averaging also indicated that there was no evidence that treatment had an effect. All models where triplet was included as an explanatory variable had a worse AICc than the null model with intercept only. This strongly indicates a degree of overfitting in models that included these parameters.
Supplementary information provided with the 2006 paper demonstrates that each triplet had a variable time at risk, ranging 2.72 years (triplet D) to 6.73 years (triplet B) (Table 3). Thus, the
Table 3
Incidence, baseline herds and incidence rates by triplet and treatment for trial areas and neighbouring areas. Figures are using the VETNET database and are for ‘confirmed’ incidents from initial cull
|
|
|
Whole trial area
|
Neighbouring area
|
Triplet
|
Treatment
|
Years since initial cull
|
Incidence
|
Baseline herds
|
Incidence rate*
|
Incidence
|
Baseline herds
|
Incidence rate*
|
A
|
Proactive
|
5.60
|
37
|
71
|
9.3%
|
24
|
60
|
7.1%
|
A
|
Survey-only
|
5.60
|
52
|
86
|
10.8%
|
20
|
68
|
5.2%
|
B
|
Proactive
|
6.73
|
87
|
152
|
8.5%
|
65
|
154
|
6.3%
|
B
|
Survey-only
|
6.73
|
61
|
132
|
6.9%
|
44
|
68
|
9.6%
|
C
|
Proactive
|
5.85
|
29
|
105
|
4.7%
|
35
|
117
|
5.1%
|
C
|
Survey-only
|
5.85
|
84
|
174
|
8.2%
|
37
|
121
|
5.2%
|
D
|
Proactive
|
2.72
|
36
|
97
|
13.7%
|
17
|
47
|
13.3%
|
D
|
Survey-only
|
2.72
|
40
|
106
|
13.9%
|
15
|
57
|
9.7%
|
E
|
Proactive
|
5.28
|
36
|
116
|
5.9%
|
26
|
94
|
5.2%
|
E
|
Survey-only
|
5.28
|
56
|
97
|
10.9%
|
31
|
72
|
8.2%
|
F
|
Proactive
|
5.13
|
15
|
138
|
2.1%
|
13
|
61
|
4.2%
|
F
|
Survey-only
|
5.13
|
61
|
191
|
6.2%
|
32
|
127
|
4.9%
|
G
|
Proactive
|
4.82
|
72
|
245
|
6.1%
|
30
|
154
|
4.0%
|
G
|
Survey-only
|
4.82
|
40
|
131
|
6.3%
|
29
|
129
|
4.7%
|
H
|
Proactive
|
4.72
|
31
|
63
|
10.4%
|
46
|
73
|
13.3%
|
H
|
Survey-only
|
4.72
|
27
|
130
|
4.4%
|
27
|
94
|
6.1%
|
I
|
Proactive
|
2.91
|
27
|
100
|
9.3%
|
19
|
70
|
9.3%
|
I
|
Survey-only
|
2.91
|
21
|
98
|
7.4%
|
9
|
60
|
5.2%
|
J
|
Proactive
|
2.88
|
34
|
114
|
10.3%
|
32
|
123
|
9.0%
|
J
|
Survey-only
|
2.88
|
36
|
123
|
10.2%
|
19
|
102
|
6.5%
|
*Calculated as number of incidents (breakdowns) per baseline herd-year at risk |
This table is a reproduction of Table 6 of the supplementary information of the 2006 paper.
effect of triplet in the 2006 paper is likely due to the variation in exposure time in addition to number of baseline herds. Time at risk was not explicitly included in the analysis of the data in the 2006 paper. There was an independent statistical auditor to oversea the the analysis of the data generated by the RBCT. In the first report of the statistical auditor it was stated that, “to some extent, the number of triplets and the years of observation are interchangeable”20. This may have motivated the investigators in the 2006 paper to use the 9 categorical variables of “Triplet” as an alternative to the single continuous variable of time at risk. Standard Poisson regression techniques would have time at risk and baseline included in the offset variables to fully control for exposure. Multiplying time at risk by the number of baseline herds achieves this, to give herd years at risk. This then can be included in the regression models as either an explanatory variable, or preferably, as an offset variable. It is not possible to use the categorical variables of Triplet as offset variables. Model 4 repeats the analyses of models 2 to 3, but in this case using herd years at risk as a replacement for baseline.
As when using baseline herds in the analysis, the only models that suggest an effect of treatment is the full model. However, this model is not supported by AICc with the most parsimonious model only including historical incidence and herd years at risk as explanatory variables. Likewise model averaging failed to support an effect of treatment on bTB herd incidents.
Linear model
The linear model where the dependent variable was incidence/herd years at risk is arguably more easily understood by the non-specialist. This was examined using model 5. Examination of the residuals indicated this was a suitable model for examination of the data. As with the Poisson models, there was insufficient statistical evidence to conclude treatment influenced ‘confirmed’ bTB herd incidents. The most parsimonious model did not include treatment as an explanatory variable.
bTB Perturbation effect hypothesis
The 2006 paper reported an increase in herd incidence in the areas neighbouring the areas where badgers had been culled leading to the proposal that badgers might transmit bTB to cattle at enhanced rates after culling. However, using the same modelling approach as for the cull areas reported above (model 6), there is insufficient evidence of such a perturbation effect. The most parsimonious model (AICc) does not include an effect of treatment.
Latent breakdowns
The analysis confirms the result of the RBCT. That is the total herd incidents (confirmed and unconfirmed) have a similar rate of occurrence in areas where badgers were culled compared to control areas. This was analysed using a (generalized) Poisson model as described above (model 7). Again, there was insufficient evidence to conclude an effect of badger culling on the incidence of bTB at the herd level.
Bayesian analysis
The model in the 2006 paper was also compared to the most parsimonious Poisson model and linear model using an alternative Bayesian analysis. Minimally informative priors were used. The same priors for the parameters were used in the model comparisons. As expected, a Bayesian analysis of the model from the 2006 paper gave almost identical parameter estimates. However, in terms of Bayes factor, the most parsimonious Poisson model was favoured over the model in the 2006 paper by a factor of 7.2x1010. There was no evidence of an effect of proactive culling on herd incidence. Indeed, the marginal likelihood indicated that models where the treatment covariate was removed completely from the model indicated better support. The posterior probability distributions of the parameters in the (overdispersed) Poisson model is given in Fig. 4.
All frequentest models, are summarized in Table 4. The R code for each model is given in the supplementary file 4.
Table 4
Statistical models analysed
|
Model
|
AIC
|
AICc
|
Number of coefficients
|
Effect of culling
reduction (CI)
|
P (for culling)
|
Post Hoc Analysis#
|
1
|
Nature, 2006
|
142.3
|
202.9
|
13
|
6.2%-29.5%
|
0.005
|
Quantile deviation significant. Some evidence of dispersion
|
2
|
log(baseline) fixed at 1*
|
NA
|
NA
|
13
|
-8.5% − 31.4%
|
NS
|
NA
|
3
3a
3b
|
Generalized Poisson model
Most Parsimonious model
Model average
|
133.2
157.7
NA
|
217.2
160.4
NA
|
14
4
NA
|
12.3% − 24.6%
NA
-8.9%-28.4%
|
< 0.001
NS
NS
|
Quantile deviation
qqplot OK.
No issues detected
NA
|
4
4a
4b
4c
4d
4e
|
Treatment retained, log(herd years at risk) fixed at 1
Most Parsimonious model
Model average
4,but log(herd years at risk) as explanatory variable, not offset
Most parsimonious model
Model average
|
156.7
166.6
NA
133.2
151.5
NA
|
217.4
168.1
NA
217.2
154.2
NA
|
13
3
NA
14
4
NA
|
0.3%% − 25%
NA
-13.3%-35.3%
12.2%-24.6%
NA
-4.7%-25.9%
|
0.045
NA
NS
< 0.001
NA
NS
|
No issues detected
Quantile deviation. qqplot OK.
NA
Quantile deviation. qqplot OK.
No issues detected
NA
|
5
5a
5b
|
Linear model
Most parsimonious model
Model average
|
-91.8
-84.3
NA
|
-31.2
-82.2
NA
|
13
3
NA
|
53.2%--16.9%
NA
41%--37%
|
NS
NA
NS
|
No issues detected
No issues detected
NA
|
6
6a
6b
6c
|
Nature, 2006, perturbation
log(baseline) replaced by log(herd years at risk). Log(herd years at risk) fixed at 1. Generalised Poisson model
Most pasimonious model
Model Average
|
135.4
150.0
149.2
|
196.0
210.7
149.9
|
13
13
2
|
-5.5%- -57.2%
9.1%- -33.0%
NA
-20.2%-13.6%
|
0.012
NS
NA
NS
|
Some evidence of dispersion
No issues detected
No issues detected
NA
|
7
7a
7b
7c
|
Nature, 2006, all breakdowns*
log(baseline) replaced by log(herd years at risk)* Triplet removed
Incidence adjusted for herd level sensitivity of test*, generalised poisson model
Most parsimonious model
Model average
|
NA
169.3
165.2
NA
|
NA
253.3
167.87
NA
|
13
14
4
NA
|
22.7%- -9.7%
16.9%- -2.3%
NA
19.7%--7.4%
|
NS
NS
NA
NS
|
NA
Quantile deviation. qqplot OK
No issues detected
NA
|
*Over or under dispersion modelled as a quasi poisson model
# Analysed by simulating residuals with the R package DHARMa
AIC and AICc are only comparable between models with the same dependent variable on the same scale. Thus linear models cannot be compared to Poisson models
Overdispersion
The 2006 paper specifically indicated that RBCT results were corrected for overdispersion. However, the model indicated by the results in Table 1 (from the 2006 paper), had no evidence of overdispersion and was a Poisson regression model (model 1). Other statistical models that were explored usually had evidence of overdispersion. This was corrected by modelling as a generalized Poisson model. It was then possible to extract AIC and AICc for model comparison. Only the full model, which included triplet and treatment suggested a significant affect of treatment, and there was some evidence that there was underdispersion in this model. However, using AICc for model selection, it was clear that the most parsimonious model only had log(herd years at risk) and log (historical three year incidence) as significant explanatory variables. Model averaging confirmed the absence of an effect of treatment (ie badger culling). AICc is now recommended for model selection rather then AIC, although the two information critera converge with large sample sizes. An important aspect is that competing models should have normally distributed residuals. Details of the theory and assumptions can be found in Burnham and Anderson11, which includes examples of selection of Poisson models using AICc. For the model from the 2006 paper, there is evidence from the qqplot that the residuals are not normally distributed (Fig. 1), and this may be because there was under-dispersion and hence model misspecification. Better specification can be achieved with the generalized Poisson model. This also confirmed that using the generalized Poisson model the residuals were normally distributed and hence using AICc for model selection is valid, at least for this set of models.
Posthoc analysis
Analysis of the residuals generated from the model in 2006 paper indicates there may be misspecification issues for model 1 (Fig. 1). The most parsimonious model, in terms of AICc, analyzed was model was model 4d. In this model, log(baseline) and Triplet were replaced by log(herd years at risk). Treatment was removed as an explanatory variable as there was no significant effect. This model had no misspecification issues (Fig. 2). The linear model 5a also had no misspecification issues (Fig. 3).
Appropriate statistical modelling
Our analysis of the original RBCT data in the 2006 paper, strongly suggests that there is insufficient evidence to conclude that badger culling reduced the rate of ‘confirmed’ herd incidents compared to control areas with no badger culling. The only statistical models that suggested an effect, was the model used in the 2006 paper and in the 2006 model where log(baseline) was replaced by log(herd years at risk) (model 4c). However, the interpretation of the results of these models, is that herd incidence of bTB was independent of number of herds studied (2006 model), or indeed the time at risk. Hence this reduces credibility from an epidemiology perspective. In addition, post hoc analysis of the statistical model used strongly suggests a degree of misspecification of the 2006 model. Although it could be argued that in terms of AIC this model produced the best fit of the count models, in terms of AICc the model was far inferior to competing models, and statistically to the better specified models.
Examination of over-dispersed models, using the generalized Poisson model also suggested the best fitting model was one where treatment had no effect. Furthermore, the linear model which analysed incidence rates had no issues with misspecification and was also unable to detect an effect of culling badgers. Finally, a Bayesian analysis also strongly suggested that there was little evidence for an effect of a proactive cull.
It seems likely that the statistical model reported in the 2006 paper fell victim to overfitting. There were 13 parameters, yet only 20 data points. As a result, the model is useful in reference only to its initial data set, which would include the specific idiosyncrasies of the data within each triplet, but it would have little predictive power.
Models with herd years at risk either as a covariate or as an offset variable are more logical, give more credible results and have the lowest AICc of competing models. A simple linear model of incidence rate as the dependent variable had little or no evidence of misspecification and failed to show an effect of culling.
The aim of badger culling is to prevent bTB herd incidents, whether they are immediately detectable or avoid detection due to false negative surveillance tests. These concerns were first raised in 2015, prior to the roll out of mass badger culling in England21 and since then policy has altered in most countries to recognise ‘all reactors’ as bTB infected for risk management purposes. Analysing the data to model total herd incidents, included undetected (ie latent) incidents, also indicated that culling had no effect on the herd incidence of bTB. In addition, when the statistical model included either baseline number of herds, or herd years at risk as a covariate rather than an offset, the parameter values were credible, indicating the number of herd bTB incidents does increase both with years at risk and number of baseline herds, which is epidemiologically credible.
The other important issue that the 2006 paper raised was the perturbation effect hypothesis where areas within and surrounding the cull area have increased herd incidence of due to higher mobility in infected badgers. However, concerns have been raised about the veracity of this conclusion22. The present analysis has also demonstrated that any such an effect is most likely to be an artefact of unsuitable analysis. Taken together, this considerably weakens the argument that badgers are responsible for the substantial proportion of the transmission of bTB between cattle herds claimed by RBCT study and subsequent iterations based upon it. It also suggests that the data, rather than the statistical model, from the RBCT is consistent with the recently published paper by Langton et al6, which failed to demonstrate a reduction in bTB incidence and prevalence associated with the culling of badgers.
It is also worth noting the data in Table 3, which is taken directly from the supplementary material of the 2006 paper. In 6 of the triplets, the incidence rate was higher in the survey only areas, whilst in 4 of the triplets the incidence was higher in the proactive area. There was also considerable heterogeneity ranging from a reduction in the proactive areas of 66% (Triplet E) to an increase of 236% (Triplet H). This, together with very marginal differences in some triplets (eg Triplet J and Triplet D), illustrates why there was likely to be no significant effect of culling in the present analysis.
Implications
The model in the 2006 paper had a credibility issue in terms of the parameter value for baseline herds, had evidence of misspecification in post hoc analysis. Furthermore, models where there was no effect of treatment also had a lower AICc. A linear model assuming normal approximations had no evidence of misspecification and failed to demonstrate an effect of culling.
This analysis has important and profound implications. Much subsequent investment in research has assumed that the coefficient value (ie -0.207) for the proactive culling parameter in the log linear Poisson model is a true reflection of the expected effect size for reduction in incidence in response to culling. There are many studies and reports which rely on this19, 23–29. The view that the RBCT analysis represented a ‘strong evidence base involving experimental studies or field data collection on bTB with appropriate detailed statistical or other quantitative analysis.’30 must now be considered unsafe. Re-examination of data from the 2006 paper, suggests this is not true and thus undermines the assumptions and results of these reports and the badger culling policy that draws from it.
This reanalysis suggests that the relationship between badgers, badger culling and bTB herd incidence is likely to be far lower than has been reported and the relationship may be minimal or extremely infrequent. Given our findings, any method of controlling the disease in badgers, such as culling or badger vaccination are equally likely to have a negligible effect in controlling bTB in cattle.