Does Size Matter to Models? Exploring the Effect of Herd Size on Outputs of a Herd-Level Disease Spread Simulator

Introduction

In countries that are historically free of significant livestock diseases such as foot-and-mouth disease (FMD), the outputs of disease spread models are a useful proxy for field information on disease behavior. This information may be used by the competent authority to compare the impacts of alternative disease control policy decisions (1–8). Traditionally, the policies of FMD-free countries such as the United Kingdom, the United States, Australia, and New Zealand rely on stamping out methods to eradicate outbreaks of FMD. This involves depopulation and thorough cleaning and disinfection of detected infected premises (IPs), tracing and biocontainment of all contacts, active surveillance to detect all clusters of infection, and intensive movement restrictions to limit disease spread. Areas of ongoing research include comparing the impacts of policies that would allow animals to be vaccinated, with policies that would cull animals on all affected farms. Furthermore, if a vaccination policy is considered, the impact of vaccinating cattle only compared with vaccinating all susceptible species (2) is of interest as there are seldom sufficient human resources and vaccine doses to target every animal.

These comparisons must consider the spread of disease in different regions, the effectiveness of a variety of control options, and the economic impact of the outbreak. This includes the cost of control measures and loss of trade due to restrictions implemented by the international community. The complexity of these decisions has driven the development of ever more complex disease spread simulation models that can now incorporate detailed information on within- and between-herd spread of disease. Some models, for example, use the count of animals on each farm to estimate infectivity according to latent periods, within-herd contact rates, and incubation periods, specific to the species and numbers of each species present in each herd on each farm. Although there is still much debate over the best modeling approach (2, 9–12), a key requirement of any spatially enabled disease spread simulator is national (or district/state/county)-level data of farm locations with susceptible animal populations. The model also requires data on the contact patterns of susceptible individuals and disease-specific information for each species represented. Few countries have publicly available and accurate animal population information at the farm level, which can be used in these models; however, the impact of inaccuracies in data on model outputs and decisions is seldom discussed. This is an increasingly critical point as these modeling activities generally make use of centrally held datasets, the accuracy of which is rarely scrutinized (13, 14), while modeling becomes both progressively more complex and more highly valued by decision makers.

The objective of this study was to test the null hypothesis that uncertainty around farm-level animal population sizes is not important when interpreting the outputs of within-herd spread FMD models. While no dataset can be expected to have an exact representation of herd size at every point in time, our study is concerned with examining the performance of an FMD model which explicitly models within herd spread using a heterogeneous herd dataset with census-based herd sizes, compared with simplified herd datasets where herd size is estimated according to herd type. Three different herd datasets are used in simulations that cover six geographic areas in New Zealand, under three different disease control strategies. Each of the geographic areas have large populations of foot-and-mouth susceptible livestock in different densities. Impacts on outbreak size and duration are assessed, and the potential implications for decision makers and competent authorities of using inaccurate data discussed.

Materials and Methods

To design an experiment that provides information on the impact of herd-level animal counts on disease modeling, a disease spread simulator that utilizes the susceptible population within a herd or farm was required. The Australian Animal Disease (AADIS) model (15), which was developed in 2015 for use by the Australian Federal Government for disease response preparedness, was used for this study. AADIS is a hybrid model of livestock disease spread and control which is designed to support emergency animal disease planning. The disease simulator uses both population-based and individual-based modeling techniques. AADIS uses the herd populations to model within-herd spread with a deterministic model, and between-herd spread with a spatially explicit stochastic agent-based model (ABM). Our models use passive first IP detection which comprises two stochastic processes: detection and reporting. Detection is defined as inspecting stock (on a farm, at a saleyard, or in an abattoir), noticing clinical signs and consulting a veterinarian. An infected herd is only a candidate for detection if it meets the minimum clinical prevalence level configured for the herd type. Reporting is defined as a veterinarian suspecting FMD, sending samples to a lab, FMD being confirmed and the Chief Veterinary Officer being notified. The detection and reporting probabilities are defined per herd type and per premises type.

The AADIS model allows for representation of a situation where any type of herd may be present on any type of farm. AADIS allows disease to spread more effectively between herds on the same farm than between herds on separate farms. The ability to parameterize the spread of disease with respect to herd type and region captures the heterogeneous nature of seasonal management practices and contact patterns. For example, a beef herd on a non-commercial farm (referred to as a lifestyle block in New Zealand) will be subject to different management and marketing practices affecting the spread of disease when compared with a beef herd on a pastoral farm. AADIS also takes species and herd size into account when estimating herd susceptibility and infectivity reflecting Tildesley et al. (15) observation that a non-linear relationship between herd size and herd infectivity/susceptibility best described data from the 2001 UK FMD outbreak (15).

Model Parameters

Herd location and herd size data were obtained from AgriBase, a commercial database of farm properties and animal populations maintained by AsureQuality, a New Zealand state owned entity (16). All AgriBase farms with no animals susceptible to FMD were removed, leaving 115,618 herds on 76,487 farms in the model (Table 1). Farms were categorized into four primary types (pastoral farming, dairy cattle farming, lifestyle farming, or pig farming) based on what farmers reported in AgriBase as their main production activity. Ten secondary herd types were created based on other livestock species present on the farm: a large and small herd type for each of deer, sheep, pigs, dairy cattle, and beef cattle, which are the farmed species in NZ susceptible to FMD (Table 1). Large or small herd type was allocated based on the farm type in AgriBase (which is specified by the farmer) to make allowance for management practices and then further divided based on the size of the herd. This made it possible for some of the effects of hobby or “lifestyle” farming to be represented by assigning a “small” herd type to those herds of any kind (or size) present on hobby farms and then allocating cut points as shown in Table 2. Cut points were chosen based on experience with farming practices in New Zealand.

Table 1. Counts of secondary herd types and primary farm types used to parameterize the Australian Animal Disease model.

Table 2. Descriptive data for each of the herd types used to parameterize Australian Animal Disease.

To explore the importance of herd size, three different herd size model parameterizations were derived from the AgriBase data. The first used the actual herd sizes reported in AgriBase and represented the real or “gold standard” dataset; the second assigned each herd a size equal to the median herd size for each of the 10 herd types. The third assigned each herd a size that was sampled from a beta-pert distribution generated from each of the 10 herd types. The herd type descriptive summaries used to generate the median and beta-pert distributions are shown in Table 2. Beta-pert distributions were selected to represent the herd size distributions based on testing of the gold standard herd dataset in a quantitative risk analysis software which identified this as being the best fit for the data (17). The minimum, maximum and mode (most likely) values were then selected to describe the best fitting beta-pert distribution for each of the ten herd types. Each simulation run sampled a different value for each herd from the constructed beta-pert distributions.

The disease-specific parameterization of the AADIS model was derived from the New Zealand Standard Model (NZSM) of FMD spread, which is represented in Interspread Plus, and models an outbreak of FMD serotype O pan PanAsia (18–20).

Outbreak Seeding

Using a random seed design across the whole of New Zealand in the AADIS model results in a large degree of heterogeneity in outbreak size with an insufficient number of large outbreaks to allow comparison of the different herd size scenarios equally. Furthermore, the population densities of farms and susceptible animals are known to influence the spread of highly infectious diseases (21) as well as the efficacy of vaccination strategies for FMD (22, 23). To address these known effects, as well as to make the simulations more representative of an economically severe outbreak in New Zealand, the territorial local authorities (TLAs) that were most likely to have FMD introduced were identified and then further ranked based on where an introduced outbreak would be most likely to spread.

The greatest risks of FMD introduction to New Zealand have been reviewed and published elsewhere (24). Based on the FAO FMD contingency plans manual (25), the greatest risk for New Zealand appears to be through the feeding of FMD-infected material to non-commercially kept pigs. As there were no data available for imported materials, we used the density of small pig herds as a proxy for the risk of introduction. The likelihood of spread was based on cattle and pig population density (Table 3).

Table 3. Description of the six New Zealand territorial local authorities (TLAs) used as disease index herd for hypothetical foot-and-mouth (FMD) outbreaks in the study.

Six TLAs were chosen to provide sufficient areas to give examples in both the North and South Islands but to limit the number of TLAs so that results are still intuitively comparable (Figure 1). The goal of the study is not to predict the distribution of outbreaks sizes but to allow the comparison of the different data quality scenarios.

Figure 1. Six New Zealand territorial local authorities selected for locations of simulation models to examine the effects of variation in herd size in the within-herd spread disease model of foot-and-mouth disease, Australian Animal Disease model.

Control Strategies

The control of disease within AADIS is a part of the ABM. Measures include movement restrictions, surveillance and tracing, IP operations, resource management, and vaccination. The emergent behavior of the ABM is the spatiotemporal spread of disease across the population and the subsequent activities to control and eradicate the disease. The disease spread pathways and control measures can be thought of as components of the ABM environment. Each component of the AADIS ABM environment operates independently (26). Three control strategies were modeled. The first was “stamping out” of identified infected farms (each FMD-susceptible herd on each of these farms is culled). The second applied vaccination to all susceptible species (with no restrictions on the number of doses). The third applied vaccination to cattle only.

Research is ongoing to identify alternative methods to controlling and eradicating the FMD virus, rather than automatic culling of sometimes healthy animals. The benefits of augmenting stamping out with vaccination for disease-free countries have been explored, and the strategy of vaccinating cattle only postulated as an effective alternative to vaccinating all susceptible species (2, 4, 27–30). Here we examine the effect of the accuracy of herd-level population information on the selection between two vaccination strategies, namely vaccinate cattle only, and vaccinate all susceptible animals (note that culling of animals on IPs is still employed in these strategies). When considering the use of vaccination (vs. stamping out only), decision makers must take into account the current World Organization for Animal Health (OIE) regulations, which restrict international trade for a country for an additional time period if it employs vaccination compared with if it employs a stamp out strategy (31).

Model Simulations

One thousand simulations were performed for each of the nine model parameterizations (three control strategies and three herd populations) in each of the six selected TLAs, giving 54,000 simulations in total. The study structure is represented in Table 4. Each outbreak simulation was seeded into a small pig herd (small pig herds were identified as the most likely entry point in NZ). The same seed lists were reused in each of the 54 different model parameterizations. Simulated outbreaks, which reached 365 days’ duration (number of simulation days) were terminated.

Table 4. Tabular representation of a study designed to test the null hypothesis that uncertainty around herd size is not important when interpreting the results of a within-herd spread FMD model.

Statistical Analysis

Those outbreaks that failed to propagate were analyzed with a logistic regression model that included data type, TLA, and control strategy as explanatory variables and failure to propagate as the outcome variable. This analysis was performed to test the hypothesis that failure to propagate was independent of control type but was associated with herd data type and TLA.

For outbreaks that were eliminated within 365 days, outbreak duration and count of IPs on the last day of the outbreak were used as outcome variables as these are both important to decision makers choosing between control strategies. Further analysis was performed between the two vaccination scenarios using count of vaccinated animals (a proxy for vaccine doses) as the outcome variable with the same explanatory variables. We included interactions between each of these terms and a three-way interaction between all explanatory variables in all models.

Australian Animal Disease model outputs were described and analyzed using the R statistical computing language, Cox proportional hazard (CPH) models were fitted using the R survival package (32–35). The data were right censored because not all outbreaks had been eradicated within 365 days when the simulations were terminated. Although the study design is balanced, the data are not, because many simulations generated outbreaks that were not detected, and were removed from further analysis (Table 5). Therefore, CPH models were run with all orderings of predictor variables to ensure that the explanatory ability of the variables was assessed conditionally on other variables in the model.

Table 5. Number of simulations that generate outbreaks that are not detected (not detected), outbreaks that last 365 days without being eradicated (right censored), and number of simulations where the outbreak is detected and eradicated within 365 days (detected and eradicated), by region and by dataset.

In order to simplify how we determined the relative contribution of each predictor variable to the response variable, we augmented comparison of p-values (which are all highly statistically significant), with a comparison of the deviance values that they index. Our reasoning is that under the null hypothesis of no term effect, the deviance follows a chi-square distribution with set degrees of freedom. Ordinarily we compare these values with the index distribution to obtain a p-value. Because of the large numbers of simulations and the strength of the effects, all p-values are very small (Table 6). Our goal was to make a statement about the relative contribution of each of the experimental variables upon the response variable, but the uniformly very small p-values are difficult to interpret in this light. Therefore, we augmented our consideration of the predictor variables by interpreting the size of the deviance values relative to their expectation under the null hypothesis of no effect, which is the same as the number of degrees of freedom. So, the relative importance of model covariates was determined by averaging the deviance values from the output then dividing by the degrees of freedom for each term to estimate the relative deviance.

Table 6. Analysis of variance (ANOVA) table for the cox proportional hazards (CPH) models with infected premises (IPs) and duration as outcome variables.

To assess the effect of herd type on choice of control strategy, the strategy that had the lowest operational cost for each seed herd was recorded. This process was repeated for each of the herd data types. This allowed the percentage agreement and Fleiss and Cohen’s Kappa statistics on the lowest cost operational option to be calculated between the real, median, and beta-pert data sets (36, 37). When one of the herd type data sets did not have a completed simulation for that seed, the seed was deleted from the dataset used for comparison. This left 4,292 of a possible 6,000 data lines to compare between scenarios.

The operational cost of each strategy selected using the real herd data was compared with the operational cost that would have been incurred had an alternate strategy (based on the suggestion of the alternate herd data set) been followed. The cost of the alternate strategy was based on the cost generated by the real herd dataset. These amounts were examined as ratios rather than absolute amounts as we wish to demonstrate the value of this information rather than to predict outbreak costs which will change over time. As an example of how this comparison was made, consider a model iteration where the lowest fixed cost management strategy was to stamp out according to the real herd size dataset. This model iteration would be found in cell 2A in Table 4. Conversely, when using the median herd size data-set the lowest operational cost corresponds with the strategy to vaccinate all species in cell 3B in Table 4. To calculate the ratio between costs using the real dataset and median dataset, the real cost of the strategy to vaccinate all species (cell 2B) was divided by the cost of stamping out (cell 2 A).

Results

Descriptive statistics for the number of simulated outbreaks that reached 365 days, the number of outbreaks that ended before they were detected, and the number of outbreaks that were detected and controlled are shown in Table 5 for each of the herd size scenarios. When the subset of simulations that ended prior to spreading were analyzed in a logistic regression model, the explanatory variables representing data type and TLA were significant, but control strategy was not significant (p > 0.05).

The distribution of outbreak duration and counts of IPs between TLAs, control options (only stamping out, culling IPs and vaccinating all susceptible animals, and culling IPs and vaccinating cattle only), and data sets (beta-pert modeled herd size, gold standard herd size, and median herd size) were compared and the results displayed in Figure 2 (duration) and Figure 3 (count of IPs). The CPH models demonstrated that for outbreak duration and number of IPs, all three explanatory variables are significantly associated with the outcomes. In addition, for each of the models a three-way interaction term was statistically significant, indicating that for each region both the control strategy and the data type are significantly associated with the outcome variable. This is evidenced by the relative deviance values shown in Table 6.

Figure 2. Distribution of duration of outbreak for models using each of three herd size data sets across six New Zealand territorial local authorities for three control strategies.

Figure 3. Distribution of the log10 of final count of infected premises (IPs) for models using each of three herd size data sets across six New Zealand territorial local authorities for three control strategies.

When assessing agreement between the three herd datasets on lowest operational cost strategy, 493 instances (12% of a total of 4,292 observations where each of the three herd data sets could be compared) produced the same recommendation. In 1,681 instances (39%), the simulations based on the beta-pert dataset identified the same lowest operational cost strategy as simulations based on the real dataset. Similarly, in 1,801 instances (42%), the median and real datasets identified the same strategy as being the most cost effective. A Fleiss Kappa statistic was 0.04 (p < 0.0001) indicating very low levels of agreement. Agreement on the most cost effective control strategy between the median and real dataset and the real and beta-pert datasets had very low Kappa statistics of 0.06 (p < 0.005) and 0.02 (p < 0.005), respectively. The size of the median ratio between the lowest cost operational strategy and the strategy indicated by the alternate (median or beta-pert) dataset was 7 (5th percentile 1.2; 25th percentile 2.5; 75th percentile 30.5; and 95th percentile 613). This distribution is shown in Figure 4 at the log scale and stratified by data type.

Figure 4. Ratio between the operational costs generated by a comparison between the strategy indicated as being the lowest by real herd data-set and the costs of the strategy that would have been chosen should an alternate herd data set have been used. The x-axis is displayed at the log scale. The median ratio is 7 (not at log scale) and the median of the cost ratio for those decisions using the beta-pert data-set is 7.5 and the median ratio for the decisions made using the median data-set is 6.64.

Discussion

Simulation studies have been performed in the UK to identify the effect of modeled farm location information on the performance of disease spread models (38, 39). In these studies, simulations using modeled farm locations do not significantly differ from those for which farm locations are drawn from real data. To our knowledge, there are no published studies in which model results using simulated animal counts are compared with results of models using real animal counts.

The accuracy of herd population data used in modeling can significantly influence the preparations for responding to disease outbreaks. In our study, those model runs that use the gold standard herd sizes result in significantly different numbers of IPs and outbreak length when compared with those runs that used the modeled or median herd sizes. The size of the operational costs incurred based on decisions made using the median and beta-pert herd sizes was seven times the cost of the decision made using the real herd data. This result suggests that accurate population datasets should be a priority to ensure effective decisions on the best available information in order to minimize the impact of disease outbreaks. The lowest operational cost was variable among the three datasets for the same seed farm (low values of Fleiss and Cohens Kappa). Although the seed herds are the same for each instance, AADIS introduces stochasticity into each run—so true agreement may be greater than what our calculations show. Another possible source of bias is that large outbreaks that run more than 365 days without finishing were removed from the comparative data analysis.

The number of model simulations that ran for unexpectedly long durations (resulting in right censoring) is strongly correlated with data type. Those simulations run with the median herd size are over represented in all TLAs, but those TLAs that have higher density of farms with susceptible species are most affected. Similarly, a beta-pert estimation of herd size results in an over representation of outbreaks that burn out before spreading further and remain undetected (Table 5). It is possible that as the model parameterization of AADIS for New Zealand has not been fully tested and had as much time invested in it as the NZSM, it could cause artifacts in these results (18). However, as the objective of this study was to examine the impact of the accuracy of herd-level populations on a disease spread model, we focus on these findings, comparing results between models to gather information on the importance of regionally representative herd size information, rather than making recommendations on specific disease control policy for FMD free countries. Particularly important in the context of our hypothesis is the large relative deviance for data types identified in the CPH models. The relative deviance for data types is approximately five times larger than the relative deviance for the next most influential variable in the model (TLA) when considering both duration and final count of IPs.

An underestimation of outbreak size or duration in a particular region based on incorrect or estimated herd and farm population information could be as damaging to response decision making efforts as an over estimation. Take the example of the duration in days in the Whakatane region (Figure 2; Table A1 in Appendix). Here, the beta-pert data would suggest that there is little benefit in augmenting stamping out with vaccination (median of 26 days duration regardless of strategy). The median herd size data suggest that stamping out will result in longer outbreaks when using the stamping out and stamping out augmented with vaccinating all species (a median of about 330 days) when compared with stamping out augmented with vaccinating cattle only (median of 215 days). Running the model with the gold standard dataset (a more accurate reflection of regional population heterogeneity), results in a more complex picture where stamping out results in longer median outbreaks (332 days) compared with either of the vaccination augmented strategies (152 days when vaccinating all species and 138 days when vaccinating cattle only).

Actual farm population information may not be available for a variety of reasons which include resource limitation and legal restrictions (14). Where actual census data are not available, then a compromise may be to use modeled population data that are conditional on region as well as herd type rather than as a single function of herd type across the country, as has been done in this study. The optimal size for the regions that would best strike a balance between representing regional heterogeneity and best use of resource is not known. It might be argued that collection of representative samples from each locality to generate the conditional populations might require so much effort that the collection and use of actual data (which have other essential uses apart from disease spread modeling) would be a better use of resources. Our models indicate that population detail is more important in some areas of New Zealand than others and this is supported by previous work on farm animal populations in New Zealand (40). The current study does not address the impact of subtle animal count biases on disease spread as it compares only the “gold standard,” beta-pert distributed and median herd sizes for each of 10 herd types. The beta-pert representation of herd size was used based on the finding that the beta-pert distribution best fit the herd size distributions in the gold standard data. The herd type median was used as the final data set to include as InterSpread Plus, the current disease spread simulator used as the NZSM for FMD spread uses a probability of disease transmission based on the median herd size to drive disease spread in the model (18). While we have reasons for choosing both comparative datasets in our study, it is likely that inaccuracy in herd size data will, in practice not be distributed uniformly among all herds in the dataset and that some herd types will be more affected than others. We hope that our study will serve as a starting point for future, more nuanced studies which will explore the effect of sector-specific data inaccuracies.

It is important to note that no resource constraints were applied in this set of model simulations. Our results reflect that (as expected) in the absence of constraints, it is preferable to vaccinate all susceptible species than to vaccinate cattle only. Note that in Table A1 in Appendix, there appears to be a paradoxical effect of vaccinating all species in the median dataset when compared with vaccinating cattle only. This is explained by the fact that model runs that exceeded 365 days were removed from the dataset. As shown in Table 5, there are more of these model runs in the parameterizations that use the median herd size dataset which leads to larger numbers of IPs. Further investigation of the response of the model to vaccine dose and human resource limitations will make interesting future work as will further investigation of impacts not limited to operational costs of an outbreak.

Our model findings are aligned with other published research that indicates that the value of vaccination is associated with the start location of an epidemic. Furthermore, while susceptible animal density does affect outbreak size, it does not alone predict infectiousness or infectivity of a herd (22, 41). The presence of multiple species on a farm, the size of the holding and the distance to the closest infective property were risk factors identified in the 2001 UK epidemic (42). The authors point out that the accuracy of herd-level population and location information would be relevant to two (if not three) of these risk factors. In addition, should an actual outbreak occur, the strain type of the virus and its specific epidemiology would be hugely influential on the effectiveness of any chosen control strategy (20).

Our study indicates that when using a disease spread simulator that explicitly represents the spread of disease within farms the quality and origin of the data used to represent herd size has significant impacts on the model results. We recommend that specific attention needs to be focused on national-level animal population datasets that results in their alignment and more efficient utilization.

Author Contributions

MA parameterized and ran the AADIS model. RB provided AADIS support and advice and made necessary code changes for AADIS to be able to dynamically define herd sizes via a distribution. AR, TH, and MA performed statistical analysis of the results. MG, TC, and MA designed the study. All authors contributed to writing the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank Robert Sanson and AsureQuality for the use of the AgriBase dataset. Special thanks to Graeme Garner for his encouragement and help with parameterizing AADIS prior to this study. Further thanks are due to the Ministry for Primary Industries in New Zealand and to the Center of Excellence for Biosecurity Risk Analysis at the University of Melbourne for funding.

References

Appendix

Table A1. Descriptive five number summaries (minimum, 25th percentile, median, mean, 75th percentile and maximum) for each of 2 outcome variables (count of infected premises and duration in days) and 3 explanatory variables for the models described in Table 4.