Introduction
Bovine tuberculosis (bTB) is prevalent in New Zealand (NZ) and has been the focus of eradication efforts [Reference Livingstone1–Reference Nugent4]. These efforts have been made ongoing for several decades and involve disease testing of livestock (cattle and farmed deer), culling of positive reactors, livestock movement controls, and control of wildlife hosts, especially of brushtail possums (Trichosurus vulpecula). Brushtail possums are an important wildlife host of bTB [Reference Coleman5–7]. Experimental control measures targeting possums have proven effective in decreasing the incidence of bTB in cattle [Reference Caley8]. In response to the bTB challenge, the broad bTB control approach was revised, and a new strategy was adopted in 2004, which was again revised in 2011 and a new plan was adopted in 2016 [Reference Livingstone1–Reference Nugent4]. The current national aims are to achieve bTB-free livestock by 2026 and bTB-free brushtail possums by 2040 and aim for NZ to be biologically free of bTB by 2055 [7].
Making predictions and evaluating them is a fundamental part of science and can be highly useful in the management of wildlife and their diseases. Predictions generated by bTB models in possums have largely remained untested [Reference Caley9]. Predictions from an earlier study indicated that bTB eradication in NZ could be achievable within 30 years from 2009 [Reference Livingstone2]. In Britain, predictions regarding breakdowns in bTB control demonstrated high success rates [Reference Karolemeas10]. The validation of predictions in wildlife management studies is encouraged as part of both causal and statistical inference. Many examples have been highlighted in the literature [Reference Hone, Drake and Krebs11, Reference Hone and Krebs12]. This study distinguishes between causal inference, which focuses on using logic and evaluating the strength of evidence to determine whether a proposed cause has the observed effect, and statistical inference, which focuses on parameter estimation and inferring population parameter values from a random sample.
The aim of this study is to evaluate the bias and precision of predictions of the number of livestock herds infected with bTB in New Zealand, as part of bTB eradication efforts. These predictions are generated by statistical modelling, which uses trends across years in bTB-infected herds, and mechanistic (causal) modelling, which uses the costs associated with control efforts. These predictions are compared with observed numbers of bTB-infected herds to evaluate bias and precision. The implications for bTB eradication are then described through future predictions.
Methods
The framework of analyses is analogous to a factorial design with three factors, each having two levels, resulting in a total of eight (23) combinations. The factors are statistical versus mechanistic inference, power versus exponential models, and no updating versus updating of predictions. Statistical modelling focuses on estimating relationships and parameters without specifying any mechanism(s). Trends refer to patterns in diseased herds over time, and since time is not a cause [Reference Harre13], trends do not reveal or infer mechanisms or causes. Mechanistic modelling focuses on evidence of a mechanism(s) generating a pattern. A mechanism is an important aspect of evidence supporting a cause-and-effect relationship [Reference Hone, Drake and Krebs11, Reference Hone and Krebs12, Reference Hill14, Reference Howick, Glasziou and Aronsen15]. The comparison of predictions not using and those using updating of predictions examines if updating of predictions influences their bias and precision.
Data consisted of the number of livestock herds (cattle and farmed deer) with at least one positive reactor to a bTB test at 30 June each year. It also included the annual costs associated with disease control and vector control published in annual reports of the sequence of agencies responsible for coordinating bTB eradication, which include the Animal Health Board, TB Free NZ, and OSPRI [7]. Disease control costs include expenses related to testing livestock herds, such as bTB skin tests. Vector control costs include expenses incurred in controlling brushtail possums. Total annual costs were calculated by adding of disease and vector costs, and the annual costs were summed starting from and including 1995 to estimate cumulative costs. Other annual costs, as detailed in the annual reports, were not included in the present study. Recent data from 2005 onwards were used, corresponding to the adoption of the new bTB eradication strategy in 2004 extending up to the present year. This period also encompasses a later review in 2011, during which a new objective was developed and adopted [7].
Statistical modelling of trends
Two models for analysing trends were evaluated: a power equation (H = aY b) and an exponential equation (H = ae bY), where H denotes the number of bTB herds and Y denotes year. Additionally, a linear regression which assumed a constant annual change in the number of positive herds, was evaluated, and preliminary analysis showed it was highly biased; hence, it was not examined further. Both the power and exponential models are based on the assumption of diminishing returns. The exponential model assumes a constant proportional change annually. In all analyses involving the power and exponential equations, logarithms with base e were used, either as loge-loge or loge–linear regression. Furthermore, for each model, the mean squared error (MSE = bias2 + variance) [Reference Snedecor and Cochran16], was also calculated. Statistical modelling allows a comparison of model predictions with management aims; however, it does not make any inferences about causes.
Mechanistic (causal) modelling
The relationship between the annual number of bTB herds (H) and the cumulative costs of eradication (C) was examined, assuming a curved relationship, as reported previously [Reference Hone17], and was derived from a cattle–possum bTB model [Reference Barlow18]. The power model was analysed using loge-loge regression, corresponding to an inverse relationship (H = aC b when b < 0) on an arithmetic scale and a linear relationship on a loge-loge scale. This power relationship is expected to occur as disease cases, like pest abundance in general, tend to be inversely related to effort (cost) [Reference Hone17, Reference Mendes19]. Additionally, considering the possibility of exponential growth or decline in disease cases [Reference Mendes19, Reference Eales20], an exponential model H = ae bC, was also assessed. In the case of exponential decline, b < 0. The costs of eradication were represented by the cumulative costs (NZ$ millions) of possum and livestock disease control, starting from and including 1995. The cumulative costs were used as they provided an overall representation of the combined effects of both current and past controls, without an attempt to separate the effects of each. These costs were not discounted over the years, reflecting a higher value placed on future costs than the case where discounting is applied. In reviews on climate change in the UK [Reference Stern21] and Australia [Reference Garnaut22], very low discount rates were used, with the latter using discount rates of 1.35% and 2.65%. Previous studies have used discount rates of 7.5%, 8% [Reference Livingstone2], and 5% [Reference Mendes19, Reference Gormley23] when discounting future costs.
Validation of predictions
In this study, validation refers to the process of comparing observed and predicted values of a parameter. The analysis involved the estimation of a regression relationship, hereafter referred to as the calibration regression, using data from 2005 to 2011 inclusive and then using the out-of-sample [Reference Lipton24] data from 2012 to 2021 to evaluate the accuracy of predictions. The predictions were of two types, namely, those derived from the original calibration regression, which projected further into the future progressively, hereafter referred to as no updating. This approach is similar to modelling trends in foot-and-mouth disease (FMD) in the UK in 2001 [25]. The second set of predictions were annually updated by re-estimating the calibration regression each year, resulting in one-year-ahead predictions, hereafter referred to as updated, which is similar to the method used for one-year-ahead predictions of mallard (Anas platyrhynchos) abundance in parts of North America [Reference Nichols, Kendall and Boomer26]. The terms prediction, projection, and forecasting are used interchangeably. The focus is on evaluating quantitative predictions of continuous data, such as abundance, rather than presence/absence data.
The bias and precision of predictions were evaluated using two measures: strength of association (which examines the relationship between predicted and observed data) and difference (which determines whether the observed difference between predicted and expected values includes 0.0, if there is no bias on average). In the association analysis, predicted data were placed on the y axis and observed data were placed on the x axis [Reference Pauwels, Guyot and Walker27].
Causal criteria of strength of association (coefficient of determination R2) and consistency [Reference Hill14] (using its inverse measure of mean squared error (MSE)) [Reference Snedecor and Cochran16] can be combined graphically (Figure 1). MSE is the vertical (y) axis, and R2 is the horizontal (x) axis (Figure 1). A perfect association in a calibration regression has an R2 of 1.0 (or 100%), and perfect predictions have no bias and a zero variance, resulting in an MSE of 0. This ideal scenario is a point located at the bottom right corner of the graph at (1,0) (Figure 1). The worst association and predictions have an R2 of 0 and a high MSE, positioning them towards the upper left corner of the graph. Multi-model analysis can measure association using Akaike weights which range, like R2, between 0 and 1.0 [Reference Anderson28], thus can be an alternative measure used on the x axis.
Results
Statistical modelling of trends: Calibration regression
The curved relationship between bTB herds (H) and years (Y) was strongly supported. The power model on a logeH vs logeY scale revealed a significant negative linear relationship (F1,5 = 137.00, P < 0.0001, R2 = 0.97) from 2005 to 2011 inclusive. The estimated slope was −328.259 (95%CI -400.350 to −256.168) and intercept was 2501.300 (95%CI of 1953.056 to 3049.544). After back-transformation, the relationship exhibited a curved, concave-up pattern (Figure 2a). The exponential model on a logeH vs Y scale also demonstrated a highly significant linear relationship (F1,5 = 137.08, P < 0.0005, R2 = 0.97), with a slope of −0.1635 (95%CI -0.1994 to −0.1276) and an intercept of 333.1882 (95%CI 261.1155 to 405.2608).
Statistical modelling of trends: No updating
In the association analyses, the predicted number of bTB-infected herds of the power and exponential models was significantly and positively related (Figure 2b) to the observed number of bTB-infected herds (Table 1). The 95%CI of the estimated slopes included 1.0 and those of the intercepts included 0.0 (Table 1), indicating that both parameters had no bias. However, in the difference analysis, the power model showed a mean bias of 12.0 herds (+/− 4.0 SE), with a 95%CI of 3.0 to 21.0, and a paired t value of 3.0 (df = 9, P < 0.025), indicating that the power model was significantly biased as it underestimated the number of bTB herds (Table 2). The exponential model exhibited a similar bias (mean = 12.3 herds, +/− 4.0 SE), with a 95%CI of 3.3 to 21.2 (paired t = 3.1, df = 9, P < 0.025) (Table 2). The mean squared errors (MSEs) of the power and exponential models were 301 and 308, respectively (Table 3, Figure 1).
Note: Analyses use statistical inference of trends in bTB herds over the years and mechanistic inference of bTB herds relative to costs. For each form of inference, there are two models (power and exponential).
Abbreviations: LCL, lower limit of 95%CI; UCL, upper limit of 95%CI.
a negatively biased.
b positively biased.
Note: Shown are the numerical values of the mean differences (+/− SE) between observed and predicted numbers of bTB-infected herds in New Zealand from 2012 to 2021 inclusive.
* P < 0.05 in a paired t test. All significant means are underestimates.
Statistical modelling of trends: With annual updating
In the association analysis, one-year-ahead updated predictions of the power model showed a significant and positive relationship with the observed number of bTB herds, with unbiased slope and intercept (Table 1). In the difference analysis, the mean bias (observed – predicted) was 5.8 herds (+/− 4.6 SE), and it was not different from 0.0 (paired t = 1.3, df = 9, P > 0.2) (Table 2). The MSE was 244 (Table 3, Figure 1). However, for the one-year-ahead updated predictions of the exponential model, there was no significant relationship with the observed number of bTB herds (Table 1). In the difference analysis, the exponential model provided biased predictions (Table 2). The MSE of the exponential model predictions was 486 (Table 3, Figure 1).
Statistical modelling of trends: Future predictions
There was a significant negative linear relationship (F1,15 = 138.10, P < 0.0001, R2 = 0.90) between number of bTB herds (H) and years (Y) on a logeH –logeY scale from 2005 to 2021 inclusive. The power model had a slope of −249.513 (95%CI -294.768 to −204.258) and an intercept of 1902.433 (95%CI 1558.161 to 2246.706). After back-transformation, the relationship displayed a clear curved, concave-up pattern (Figure 3). The exponential model on a logeH vs Y scale also demonstrated a similar relationship (F1,15 = 137.65, P < 0.0001, R2 = 0.90), with an estimated slope of −0.124 (95%CI -0.147 to −0.101) and an intercept of 253.766 (95%CI 208.443 to 299.088). Both fitted regressions predicted that the number of positive bTB herds was 24 in 2022, which coincided with the observed number. They also projected a decrease to 3 in the year 2040 and ultimately the eradication of bTB by 2055 (Table 4).
Note: The selected years correspond to aims of the national programme of bTB eradication with either or both of livestock and brushtail possums ‘free’ of bTB. The independently observed number of bTB herds in 2022 was 24 [7].
Abbreviation: NA, not applicable.
Mechanistic modelling: Calibration regression
The curved relationship between bTB herds (H) and cumulative costs (C) was strongly supported. The power model showed a significant negative linear relationship (F1,5 = 107.29, P = 0.00014, R2 = 0.96) between bTB herds (H) and cumulative costs (C) on a logeH vs logeC scale. After back-transformation, the relationship with predictions was a clear curved, concave-up pattern (Figure 4a). The exponential model showed a significant negative linear relationship (F1,5 = 136.230, P < 0.00001, R2 = 0.97) on a logeH vs C scale. After back-transformation, this relationship also displayed a clear curved, concave-up pattern (Figure 4a).
Mechanistic modelling: With no updating
In the association analysis, the power model showed a significant positive linear relationship between predicted and observed values (Table 1, Figure 4b). However, the slope of the linear regression (0.501, 95%CI 0.191 to 0.812) was less than 1.0, indicating a negative bias, and the intercept (28.733, 95%CI 11.987 to 45.479) was greater than 0.0, indicating a positive bias (Table 1). In the difference analysis, the mean bias of predictions was −3.7 (+/− 4.1 SE), which was not significantly different from 0.0 (paired t = −0.9, df = 9, P > 0.20) (Table 2). On the other hand, in the association analysis, the exponential model showed that predicted values were significantly related to the observed values of herds with bTB (Table 1, Figure 4b). The estimated slope was 0.656 (95%CI 0.253 to 1.059), which was not different to 1.0 and hence not biased, and the intercept on the y axis was 6.985 (95%CI -14.760 to 28.730), which was not different to 0 and was unbiased (Table 1). However, in the difference analysis, the mean difference between observed and predicted bTB herds was 10.3 herds (+/− 3.9 SE), indicating bias (paired t = 2.6, df = 9, P < 0.05), with predictions tending to be, on average, underestimates (Table 2). The mean squared errors of the power and exponential models were 181 and 261, respectively (Table 3, Figure 1).
Mechanistic modelling: With annual updating
In the association analysis, the power model generated one-year-ahead predictions of the number of livestock herds with bTB, and these predictions were positively associated with the observed number of bTB herds (Table 1, Figure 5). However, the fitted linear regression had an estimated slope of 0.551 (95%CI 0.126 to 0.976), indicating a negative bias, and an estimated intercept of 23.830 (95%CI 0.882 to 46.777), suggesting a positive bias (Table 1). In the difference analysis, the mean bias (observed – predicted) was −1.2 herds (+/− 4.5 SE), which was not different from 0.0 (paired t = −0.3, df = 9, P > 0.5) (Table 2). The mean squared error of this model was 203 (Table 3, Figure 1). Similarly, in the association analysis, the exponential model provided one-year-ahead predictions that were positively associated with observed herds (Table 1, Figure 5). However, the estimated slope was 0.524 (95%CI 0.098 to 0.950), indicating a negative bias, and the intercept was 19.428 (95%CI -3.547 to 42.404), indicating no bias (Table 1). In the difference analysis, the mean bias was 4.5 herds (+/− 4.6 SE), which was not different from 0.0 (paired t = 1.0, df = 9, P > 0.20) (Table 2). The mean squared error was 233 (Table 3, Figure 1). The mechanistic power model with updating shows a temporal trend of damped oscillations in bias (Figure 6). This suggests that differences between the observed number of bTB herds and the predicted number vary over the years, but the amplitude of the differences decreases, approaching zero (Figure 6).
Mechanistic modelling: Future predictions
Predictions were also made for the future number of herds with bTB using data from 2005 to 2021 inclusive. The power model showed a high level of significance (F1,15 = 163.72, P < 0.0001, R2 = 0.92). The fitted equation is as follows:
The equation predicted that the number of bTB herds (H) would be 31 in 2022, 24 in 2026, and 7 in 2055 (Table 4), assuming annual costs remain at the current level of NZ$59million. The actual number of bTB herds in 2022 was 24. According to the fitted power equation, there is a marginal benefit of a reduction of 6 bTB herds when the cumulative cost increases from NZ$1,600 million to NZ$1,800 million and a reduction of 2 bTB herds when the cumulative cost increases from NZ$2,800 million to NZ$3,000 million.
The exponential model was highly significant (F1,15 = 159.11, P < 0.0001, R2 = 0.91). The fitted is as follows:
The exponential equation predicted the number of bTB herds (H) would be 26 in 2022, 17 in 2026, and 1 in 2055 (Table 4). The actual number of bTB herds in 2022 was 24. According to the fitted exponential equation, there is a marginal benefit of a reduction of 8 bTB herds when the cumulative cost increases from NZ$1,600 million to NZ$1,800 million and a reduction of 1 bTB herd when the cumulative cost increases from NZ$2,800 million to NZ$3,000 million.
The power model estimates that reducing the number of bTB-infected herds to 1.0 would require a cumulative cost of NZ$10,437 million. On the other hand, the exponential model estimates that achieving the same goal would require a cumulative cost of NZ$3,382 million. In 2021, the cumulative cost was NZ$1,587 million.
Discussion
The present study has demonstrated that the number of bTB-infected herds in New Zealand (NZ) can be predicted with low bias. Mechanistic models that imply causal relationships received more support than statistical models, and these models described diminishing returns on investments. Specifically, the power model had less biased predictions than the exponential model. Both measures of association and difference yielded some similar findings, but there were also some differences in the results. Future predictions suggest that bTB eradication may take longer than 2055 if the predictions remain unbiased and current strategies and costs of control are maintained.
The control and eradication efforts for bTB in NZ have shown impressive progress [Reference Livingstone1, Reference Livingstone2, Reference Nugent4]. Previous predictions regarding the feasibility of bTB eradication [Reference Livingstone2] are broadly supported by the findings of the present study, although some differences exist. The current predictions imply that NZ’s livestock will likely be provisionally free of bTB by 2026, which is consistent with the current situation. The statistical models predict biological freedom from bTB by 2055, aligning with the national aims. In contrast, the mechanistic models do not predict biological freedom by 2055 (Table 4), assuming that funding for disease control and possum control remains at current levels. It is worth noting that as the number of bTB-infected herds reaches very low levels, predicting the number may be become difficult because of the very low numbers. Additionally, if strategies change, for example, depopulation of bTB-infected herds and using sentinel animals to detect infection in livestock and wildlife, it may further complicate predictions.
The study does not formally evaluate the aim of achieving bTB-free possums in NZ by 2040 because it is challenging to distinguish the separate effects of concurrent livestock disease control and possum control. Other studies have evaluated various aspects of surveillance and control of possums as bTB hosts [Reference Gormley23]. It is noteworthy that Australia eradicated bTB in livestock [Reference Bradshaw29, Reference More, Radunz and Glanville30], though Australia did not have a widespread wildlife bTB host, in contrast to possums in NZ. Of the cumulative costs incurred since 2004, 75% (NZ$841.2 m) was spent on possum control and related administration and research, and 25% (NZ$287.4 m) was spent on disease testing and control in livestock. Advancements in disease surveillance, including efforts in wildlife [Reference Gormley23, Reference Anderson31, Reference Anderson32], increase the probability of achieving bTB eradication.
Given the substantial costs associated with bTB eradication, it is important to demonstrate benefits [Reference McInerney, Howe and Schepers33]. In this study, the clear reduction in bTB-affected herds over the years serves as a demonstrated benefit. It is important to note that the cumulative costs reported here may not be directly comparable to the previously reported results [Reference Livingstone2]. The curved trends observed in bTB herds in this study (Figure 3) are similar to the curved trends observed in countries with smallpox during the final stages of its eradication [Reference Fenner34, Figure 10.4]. A previous study [Reference Hone17] and the current study clearly demonstrate diminishing returns on investments. This pattern of diminishing returns holds implications for other countries engaging in bTB eradication, such as Ireland, which has shown a trend in bTB cases with a possible asymptote above zero, rather than a decline towards zero [Reference More35].
A comparison between statistical and mechanistic modelling of climate change revealed that mechanistic modelling yielded less biased results than statistical modelling [Reference Hausfather36, Reference Supran, Rahmstorf and Oreskes37]. Similarly, the results reported here show that mechanistic models produced fewer biased predictions than statistical models (Table 2). The damped oscillations in the bias of the updated mechanistic power model are encouraging and implies increased confidence in the predictions. These damped oscillations are a form of convergence, which has been reported as the correlation between observed and predicted values increasing with larger datasets, suggested as a metric for causal inference [Reference Sugihara38]. However, correlation alone is not a basis for causal inference. Occurrence of convergence through updating of predictions highlights the importance of annual updating to assess on-going progress in bTB eradication efforts. Annual updating is already established in the management of mallards and their harvest in North America [Reference Nichols, Kendall and Boomer26]. The utility of validating predictions both for statistical and causal inference has been emphasised in wildlife studies and management [Reference Hone, Drake and Krebs11, Reference Hone and Krebs12]. The present study demonstrates that this can be extended to livestock production and disease control.
The study has several limitations. The use of cumulative costs for bTB control may introduce limitations. Firstly, how the same amount of money is spent may differ between years, thus generating extra variation in the x variable in the mechanistic models. This variation could reduce the estimated slope of the linear regression between predicted and observed numbers [Reference Snedecor and Cochran16] in the association analyses and lead to biased predictions. Secondly, in some analyses, the x variables are not independent as they use cumulative costs. This may increase the risk of a type 1 error (concluding there is a significant relationship when it actually does not exist); however, due to the lack of multiple sites, there was no alternative. Thirdly, this study is observational, rather than experimental, which may lead to weaker causal inferences. Fourthly, the use of bTB numbers from one site across multiple years could provide more precise data than having the same sample size from different sites in a single year. The higher precision might increase the risk of a type 1 error.
In summary, the present study concludes that mechanistic (causal) models provide less biased predictions of the number of bTB-infected herds in New Zealand than statistical models. The results indicate that New Zealand is progressing towards bTB eradication, although it may take longer than the target year of 2055.
Data availability statement
Data are available from the author on request.
Acknowledgements
The author acknowledges the support of the Institute for Applied Ecology, Jason Weber, and the University of Canberra.
Author contribution
Conceptualization: J.H.; Data curation: J.H.; Formal analysis: J.H.; Methodology: J.H.; Writing – original draft: J.H.; Writing – review & editing: J.H.
Financial support
The study had no external funding.
Competing interest
The author declares none.
Ethical standard
The study did not require animal ethics committee approval.