Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 https://doi.org/10.1186/s12917-020-02312-8 METHODOLOGY ARTICLE Open Access Enhancing the monitoring of fallen stock at different hierarchical administrative levels: an illustration on dairy cattle from regions with distinct husbandry, demographical and climate traits Amanda Fernández-Fontelo1,2*, Pedro Puig2, German Caceres3, Luis Romero3, Crawford Revie4,5, Javier Sanchez4, Fernanda C. Dorea6 and Ana Alba-Casals7,8 Abstract Background: The automated collection of non-specific data from livestock, combined with techniques for data mining and time series analyses, facilitates the development of animal health syndromic surveillance (AHSyS). An example of AHSyS approach relates to the monitoring of bovine fallen stock. In order to enhance part of the machinery of a complete syndromic surveillance system, the present work developed a novel approach for modelling in near real time multiple mortality patterns at different hierarchical administrative levels. To illustrate its functionality, this system was applied to mortality data in dairy cattle collected across two Spanish regions with distinct demographical, husbandry, and climate conditions. Results: The process analyzed the patterns of weekly counts of fallen dairy cattle at different hierarchical administrative levels across two regions between Jan-2006 and Dec-2013 and predicted their respective expected counts between Jan-2014 and Jun- 2015. By comparing predicted to observed data, those counts of fallen dairy cattle that exceeded the upper limits of a conventional 95% predicted interval were identified as mortality peaks. This work proposes a dynamic system that combines hierarchical time series and autoregressive integrated moving average models (ARIMA). These ARIMA models also include trend and seasonality for describing profiles of weekly mortality and detecting aberrations at the region, province, and county levels (spatial aggregations). Software that fitted the model parameters was built using the R statistical packages. Conclusions: The work builds a novel tool to monitor fallen stock data for different geographical aggregations and can serve as a means of generating early warning signals of a health problem. This approach can be adapted to other types of animal health data that share similar hierarchical structures. Keywords: Fallen stock, Dairy cattle, Syndromic surveillance, Hierarchical time series, ARIMA models, Spain * Correspondence: amanda.fernandez-fontelo@hu-berlin.de; amanda@mat.uab.cat 1 Chair of Statistics, School of Business and Economics, Humboldt Universität zu Berlin, Berlin, Germany 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Barcelona, Spain Full list of author information is available at the end of the article © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 2 of 13 Background detecting local and moderate variations was limited. The Animal health surveillance is essential for planning and problem was that the selection of an appropriate and implementing efficient prevention and control measures parsimonious ARIMA model accounting for seasonality that are important for animal and public health [1, 2]. and trend for specific time series involved the testing Current advances in data mining techniques and spatial- and checking of many combinations of values for each temporal analysis facilitate the development of novel ani- parameter. This process was not automatic and model- mal health syndromic surveillance (AHSyS) approaches. ling many time series from different subpopulations was These methods allow us to obtain information regarding very tedious. Previously, Alba et al. plotted the time the health status of animal populations to be extracted series of smaller subpopulations using HTS structures. from diverse automated data sources of a non-specific However, these HTS did not quantify temporal patterns, nature in near real-time [3–5]. AHSyS may enhance forecast mortality based on previous observations, or de- traditional animal surveillance systems by identifying tect aberrations over time at different administrative subpopulations at high risk, assessing the impact of pre- levels. The present study solves this problem by creating vious prevention or control measures, supporting claims an automated process to analyze patterns of multiple of freedom from disease, and serving as an early warning fallen cattle subpopulations with lower numbers and system [6]. This AHSyS approach involves the continu- high variability. This modelling process aimed to analyze ous analysis of mortality data registered at farm level. and predict fallen stock patterns at both low and high The potential of mortality data for this purpose has been geographical levels. Moreover, in this work, we used as demonstrated in previous studies conducted in a num- an indicator of fallen dairy cattle the “counts of fallen ber of European countries [7–12]. In Catalonia (North- bovines”, which is more accurate than the “counts of Eastern Spain), patterns of fallen bovine were modelled, carcass disposal visits” or the “total of kg collected” used for the main cattle production types between 2006 and before by Alba et al. 2015 [7]. 2013, using time series analyses [7, 8]. In those studies, To illustrate the functionality of our method in differ- bovine fallen stock data weekly aggregated at the re- ent contexts, it was put in motion by using the retro- gional level were fitted using autoregressive integrated spective fallen stock data of dairy cattle collected and moving average (ARIMA) models. These models between 2006 and 2015 in two regions of Spain with dif- included trend and seasonal patterns as covariates. ferent demographical structure, husbandry systems, and According to these ARIMA models and assuming regu- climate conditions. larity in the bovine population patterns over time, both the number of visits to carcass disposal and the total Results weight collected in these visits (in kg) were predicted at These results describe the basic traits of two Spanish re- n-week ahead. These models provided information to as- gions with distinct dairy cattle populations (R1 and R2) sess the impacts of different events occurred over time and the monitoring of counts of fallen bovines collected at the regional level. However, to identify high-risk sub- by week at three administrative levels (i.e., region, prov- populations and/or to serve as an early warning system, ince, and county) over a 10-year period. the number of fallen bovines must be assessed at more The region R1 had registered 3469 dairy farms (with spatially discrete administrative aggregations. With this an annual median of 2958 dairy farms and 117,572 heads aim, fallen bovine data aggregated at province and of cattle), while region R2 had 984 farms registered (with county levels were plotted using hierarchical time series an annual median of 797 farms and 124,953 heads of structures (HTS) according to the methodology pro- cattle). This system covered 77 and 81% of all these dairy posed by Hyndman et al. [13, 14]. These plots indicated farms of R1 and R2, respectively. Figure 1 shows the ad- that patterns in fallen bovine data varied substantially ministrative levels monitored in this study, i.e. across R1 among different subpopulations. This implied that mul- seven of its 17 counties, and across R2 two of its four tiple time series would have to be modelled at more re- provinces, and in seven of its 26 counties. fined spatial scales in order to get well-fitting models for Over this period of study in these two regions, the different spatial aggregations. number of dairy cattle farms decreased and the sizes of This work is a methodological extension to the model- the herds increased; but the total number of bovines in ling approach introduced in the previous papers of Alba these populations did not change substantially. R1 had et al. [7, 8]. In [7, 8] the indicator of fallen stock was 3.4 times more dairy farms than R2 and the herd size only modelled at the region level. Although the model- was 2.8 times smaller. This difference was also evidenced ling of these data aggregated at higher levels allowed the consistently at smaller administrative levels, such as quantification of the total impact (or lack thereof) on a provinces and counties. The proportion of herds in ex- population in the event of a potential animal threat, tensive production also varied, being around 0% in R2 using only this level of analysis, the likelihood of and 25% in R1. Regarding the climate, this was quite Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 3 of 13 Fig. 1 Spanish regions, provinces, and counties within which the levels of fallen dairy cattle were monitored. The figure is own-created similar across counties of R1, but differed among the seasonal effects, was fitted based on the training dataset counties of R2. Moreover, in some counties of R2 the collected between Jan-2006 and Dec-2013. high proportion of herds housed all-year-round com- The main traits of the two dairy cattle subpopulations bined with a continental climate with warm summers and the monitored fallen stock data are summarized in could cause heat stress problems and an increase of Table 1. mortality. Initially, each administrative level was de- In most of these administrative levels of both regions scribed through the number of dairy farms, the herd R1 and R2, the fallen dairy cattle presented annual and size, the type of husbandry, the climate, and the number biannual seasonality, with an increasing trend over time. of fallen bovines collected in total and by week. This Interestingly, this pattern was consistent for all the preliminary exploration indicated that most of these sub- counties of R2, but differed among counties of R1. Both populations differed substantially in relation to demo- regions (i.e., all the series in R2 and six series of R1) ex- graphical and climate traits as shown in Table 1. hibited an increased trend and seasonal pattern that rose Plots of the HTS of weekly counts of fallen bovines at during January and February, coinciding with the coldest region, province and county levels are shown in S2 and season (Figs. 2 and 3). During July and August, in the S3 Figs. warmest period of the year, there was also an increase of The respective plots of HTS for each region indicated mortality, but this rise was only very evident in R2 and that the fallen stock patterns varied among subpopula- especially in some counties (see Table 1). The weekly tions. For this reason, it was important to model these counts of fallen bovines collected were much lower in patterns at each administrative level to obtain robust in- R1 than in R2 (i.e., a median of 221 heads in R1, versus formation, detect and investigate abnormal mortality 308 in R2). However, the number of farms involved in events and support the decision-making. each peak was higher in R1 than in R2 (a median of 145 For each administrative aggregation the most appro- farms in R1, versus 38 in R2). Moreover, it should be priate ARIMA model, which could include trend and noted that in R1 the trend and seasonal patterns among Table 1 Description of dairy cattle populations under monitoring with their respective patterns of mortality and mortality aberrations detected at each level Zones Demographical and climate traits of populations Carcasses collected Patterns of fallen bovines by week (Jan 2006- Dec 2013) Mortality Aberrations of study monitored (Jan 2014-Jun 2015) Nr farms Farm size Climate Total nr By week Median ARIMA Trend Seasonal effect (Warm) (% Seasonality effect (Cold) Nr of Nr of affected farms (Median) Median (min-max) (p, d, q) of increase and period) (% of increase and period) peaks Median (min-max) (min-max) R1(=P1) 2681 74 (1–561) Oceanic coastal 109,744 221 (151–326) 1,0,1 + 0.086 8% Jun-Aug 22% Oct-Feb 4 145 (105–184) C1 343 77 (1–369) Oceanic coastal 16,114 32 (13–61) 0,1,1 + 0.007 14% May-Aug 26% Oct-Feb 6 35 (30–44) C2 302 71 (1–407) Oceanic coastal 15,440 31 (10–58) 0,1,1 _ 9% May-Aug 17% Oct-Feb 4 28 (13–37) C3 400 63 (1–561) Oceanic coastal 15,145 30 (14–58) 1,0,1 + 0.009 16% May- Aug 16% Oct -Feb 6 29 (3–35) C4 425 75 (1–430) Oceanic coastal 17,205 34 (12–68) 0,1,1 + 0.022 8% May-Aug 11% Oct - Feb 2 48 (44–51) C5 143 84 (2–240) Oceanic coastal 5924 12 (3–25) 1,0,1 + 0.004 _ 15% Oct-Jan 4 15 (13–18) Fernández-Fontelo et al. BMC Veterinary Research C6 314 82 (1–487) Oceanic coastal 14,850 29 (13–62) 0,1,1 + 0.125 22% May-Aug 35% Oct-Jan 1 30 C7 71 88 (5–218) Oceanic coastal 4350 8 (0–20) 1,0,1 −0.006 _ 15% Oct-Mar 2 10 (9–11) R2 799 198 (1–3639) Diverse climates: 153,520 308 (144–502) 0,1,1 + 0.229 61% May-Aug 85% Oct-Jan 1 38 Mediterranean and Alpine (2020) 16:110 P2 212 220 (1–3639) Mediterranean 49,557 104 (40–197) 0,1,1 + 0.110 96% May-Aug 138% Oct-Jan 1 49 continental C8 22 297 (3–3369) Mediterranean 9331 17(2–52) 0,1,1 + 0.037 145% May- Set 184% Oct- Jan 7 5 (2–9) continental warm summers C9 21 526 (14–2005) Mediterranean 10,309 19 (3–55) 2,1,2 + 0.031 145% May- Set 164% Oct- Jan 4 7 (1–7) continental warm summers C10 98 192 (1–1556) Mediterranean 18,107 36 (9–75) 0,1,1 + 0.040 83% May-Aug 125% Oct- Jan 1 28 continental subhumid C11 25 206 (17–1403) Mediterranean 9055 17(4–58) 4,1,2 + 0.012 100% May-Set 134% Oct-Jan 1 5 continental warm summers P3 308 191 (6–1933) Mediterranean 56,274 106 (54–200) 1,0,1 + 0.079 59% May-Aug 78% Oct- Jan 7 5 (2–9) Continental subhumid /costal C12 41 197 (6–559) Mediterranean 7701 15 (2–37) 1,0,1 + 0.009 72% May-Aug 111% Oct-Jan 1 11 continental subhumid C13 61 231 (7–905) Mediterranean costal 14,265 28 (11–73) 0,1,1 + 0.021 67% May- Aug 72% Oct- Jan 3 9 (8–15) C14 54 228 (6–1933) Mediterranean costal 14,778 29 (10–70) 3,1,1 + 0.011 41% May-Aug 72%Set-Jan 2 9 (7–10) Zones of study: R Region, C County, P Province ARIMA(pdq) where p = order of autocorrelation, d = differentiation, q = order moving average Page 4 of 13 Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 5 of 13 Fig. 2 Time series plots of fallen bovines in the training data set collected weekly together with mortality peaks (A) detected in the testing dataset for R1 (P1) counties were more disparate and irregular, mainly in a province level would be also detected at a county level counties with lower counts such as C5 and C7. Finally, (see Figs. 2 and 3). based on these models, the fallen bovines expected by In both regions there was a marked increase in counts week were forecasted for each subpopulation during Jan- of mortality during Jan-2015, being identified as unusual 2014 and Jun-2015. Those observations that exceeded events in all administrative aggregations except for one the predicted upper confidence limits were signaled as county (C14 in R2). mortality aberrations. Considering a conventional upper 95% confidence Discussion limit, six mortality aberrations would be signaled at a Despite the growing interest of using syndromic infor- province level (four in R1 and two in R2) affecting a total mation for AHSyS [6], its implementation as a part of an of 553 farms (466 in R1 and 87 in R2). Two of the four early warning system has to be validated in the field yet. peaks detected in R1 at a province level were also de- In many European countries, bovine fallen stock data tected at a county level, while in R2 all peaks detected at have been gathered consistently over several years by Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 6 of 13 Fig. 3 Time series plots of fallen bovines in the training dataset collected by week together with mortality peaks (A) detected in the testing dataset for R2 Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 7 of 13 carcass disposal services, covering a significant propor- aggregation, forecast the counts of fallen dairy cattle that tion of active farms. The potential use of bovine fallen would be expected at n-ahead period and identify those stock data for AHSyS was described in previous studies counts that exceeded the upper 95% confidence limit, as [7–12]. In [9, 11], some authors also demonstrated the presumed mortality peaks. This system can be applied to usefulness of these methods for assessing the impact of different dairy cattle population structures (such as those past events, such as the incursion of bluetongue or the seen in R1 or R2). heat waves in 2003 and 2006 in the French cattle popu- Classical time series models such ARIMA had the ad- lation. However, many technical challenges must still be vantage that they can deal with trend and seasonal com- faced when putting in place the monitoring of bovine ponents as well as temporal correlation structures, fallen stock as part of an alert system or as a support which in the case of fallen bovine data were especially tool for decision making [7, 10]. marked (see S1 and S2 Tables). The trend and seasonal One of the main challenges relates to the heterogen- components could be modelled using trigonometric eity of the animal populations under surveillance. In the functions [15], simple exponential smoothing [16, 17] or context of animal health surveillance, the target popula- Holt-Winters [18, 19]. The temporal autocorrelation tion usually comprises different subpopulations with dis- could be fitted using autoregressive integrated moving parate managements, structures and risks. Moreover, the average (ARIMA) models [20, 21]. In this work ARIMA official animal health services of different administrative modelling, including possible trends and seasonal effects levels coordinate which prevention and control measures as covariates, demonstrated to best fit data related to should be applied in each population. The modelling ap- weekly counts of fallen bovines without requiring exces- proach proposed in this study allows for an assessment sive computation time and using standard software. of bovine fallen stock baselines and the detection of Moreover, if the response variable would depart from mortality peaks in different populations. The results of the normality assumption (e.g. time series with counts this study indicate that the baselines of fallen bovines less than 10 in average), these models could be replaced vary among different administrative aggregations. Conse- by Integer-valued Autoregressive (INAR) models [12, quently, an effective monitoring system based on fallen 22]. Other more sophisticated methodologies of time stock data should be able to provide information and de- series have been proposed. For instance, Le Strat and tect alerts for each population in near real time. Carrat [23] used Hidden Markov Models (HMM) to In this work we combined the expertise from veterin- monitor epidemiological data by segmenting series into ary epidemiologists and statisticians to lead to a robust epidemic and non-epidemic unobserved phases. These tool that constitutes an important step towards its im- unobserved phases could be modelled by a two-state plementation. Its functionality has been illustrated across homogeneous Markov chain of order 1 (epidemic and two regions of Spain from populations with different non-epidemic); higher orders could also be considered. husbandry systems and environmental conditions. Ac- As an alternative, Martinez-Beneito [24] proposed a cordingly, this study demonstrated the capability of this Markov switching model in which the response variable approach to analyze, model and predict, in an automatic did not depend only on the hidden states, but also on way, the number of fallen dairy cattle by week for differ- the lagged observable variable. However, according to ent subpopulations. [25], these methodologies both require extensive compu- The system used as an indicator of mortality the num- tation time and are therefore not suited to monitoring ber of carcasses collected each week, assuming that the epidemiological data in near real time. Furthermore, in studied populations remained relatively constant over general, Bayesian methods were computationally inten- time. To support this assumption, the behavior of the sive and strongly dependent on prior information, which dairy population was explored year-by-year for both re- in our case did not exist [26]. Parametric regression gions. This exploration indicated that, although the num- techniques, including a variety of models such as simple ber of dairy cattle farms decreased throughout the study, regression, could be explored. These include naive re- the number of dairy cattle heads slightly increased. gression that estimated the response variable at time by Initially, the weekly counts of fallen bovines were de- the mean of responses at times t-1, t and t + 1 to adjust scribed and plotted for all administrative levels using for possible seasonal effects [27] or extended regression HTS. This technique was used to explore, compare fallen that incorporated trigonometric functions with linear dairy cattle patterns among different administrative levels, trend accounting for both trends and seasonal effects and select those series that would be immediately mod- [28]. Additionally, when the response variable departed elled using a time series methodology such ARIMA mod- from normality (i.e. low counts with many zeros), some elling, including trend and seasonality as covariates. authors [29–31] have proposed Poisson regression A computer programming routine was developed to fit models with logarithmic link to model possible trends an ARIMA model for each selected administrative and seasonal patterns. Hierarchical generalized linear Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 8 of 13 mixed models [32] are another methodology for fitting also important to underline that there were mortality non-normal correlated data. Another approach of semi- peaks that were only detected at a county level, and were parametric regression could be considered; these com- thus likely linked to the occurrence of more local events. bine parametric models for representing series data and These results provide evidence that the analysis of bo- non-parametric models for including possible trends and vine fallen stock data throughout several hierarchical ad- seasonal effects, such as smoothing methods [33] and ministrative levels may bring relevant information to generalized additive models [34]. Although these regres- assess the evolution and the impact of different events sion methods take into account trends and seasonal occurred and/or detect earlier the start of a health threat components, the major disadvantage of using these tech- that can spread to other areas. The assessment of mor- niques for fallen stock data is that the models did not tality at lower aggregation level may put in evidence the consider the temporal correlation structure. occurrence of a local health problem (e.g. a food-borne Similar drawbacks exist for methods related to statistical infection locally transmitted), but also may sign the start control process, which are mainly based on control charts of a contagious disease outbreak. In this last case, if mor- (cumulative sum and exponentially weighted moving aver- tality counts are uniquely studied at higher aggregation age) [35, 36], temporal scan statistics [37] or intervened level, some problems may be unnoticed and get worse due times [38]. As in the case for regression techniques, these to the lack of implementation of prevention measures on models seemed to be simpler in terms of interpretation due time. Contrarily, the analyses conducted at higher ag- and computation, but they did not adequately handle the gregation level provide information to determine more autocorrelation structures of our data. easily the domain and global impact of a specific event on Other possible alternatives might be those based on the study population and support the policy decision- models which take into account spatial information (i.e. making to coordinate actions of prevention and/or control cumulative sum charts, temporal scan statistics or spatial at this level if it is necessary. Our outcomes (see Table 1 regression) [39–43] or methods of multivariate outbreak and Figs. 2 and 3) are consistent with the rationale that detection that dealt with several correlated series related the analysis of this kind of data, aggregated at different to the same process [44, 45]. In contrast to previous hierarchical levels, enhances the detection of unusual techniques, spatial and multivariate models allow for the aberrations, both globally and locally. incorporation of correlation structures in the data, but there was the inconvenience that they would likely re- Limitations of the study quire the inclusion of new data, such as point locations Despite the benefits of our approach for AHSyS use, it is of farms [39–43], or the application of dimension reduc- also important to discuss some constraints. First, the use tion methods [46]. of counts of fallen cattle as a proxy measure of mortality without taking into account the farm size could cause an Comparison of baselines among different populations over-expression of larger farms, masking unusual mor- obtained from the models tality events occurred in small farms. Besides, the base- Based on the selected ARIMA models, we observed that lines of fallen bovines might vary according to other in region R2 all counties exhibited an increasing linear factors, such as age, breed or sex. To enhance the accur- trend, and an annual or biannual seasonality. This region acy of the system and identify unusual events of mortal- (and its provinces and counties) presented a more ity for different populations, it would be important to homogeneous behavior in mortality data than region R1. include these covariates in the future. Unfortunately, This method could provide relevant information to these covariates could not be incorporated within our point out or discard different health problems in these present study due to the lack of accuracy of the available populations by comparing the number of mortality peaks records relating to these fields. detected among subpopulations, the time and season of Second, it is important to recall that ARIMA models detection, and the number of farms involved in each can only be used when populations and subpopulations peak. For example, some mortality peaks detected at a show regular average counts. On this matter, it was first regional level were also detected at a province and necessary to disaggregate and describe the series using county levels (marked with dotted lines in Figs. 2 and 3). HTS, to select those that could adequately be fitted In particular, in R1, there were two peaks in weeks 466 using ARIMA models. Thus, the current results only and 470 that were detected at a region and county (C1, covered a part of the subpopulations at province and C2, and C3) levels. On the other hand, in R2, there was county levels. A future innovation of the system to a mortality peak in week 469 that was identified at a re- model the mortality in these subpopulations with very gion, province (P3), and county (C14) level. This sug- low counts would be to use integer-valued autoregres- gested that the cattle populations of these hierarchical sive models (INAR) [22], Hermite integer-valued autore- series were affected by a common cause. However, it is gressive models (HINAR) [12], Bayesian approaches for Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 9 of 13 seasonal count time series [20] or one of the other tech- production. Region R2 was characterized by a complex niques reviewed in [25]. and diverse orography and a combination of a continen- Third, the use of historical data without any prior ex- tal, coastal Mediterranean and alpine climates. These re- ploration or filtering as a training set presented important gions were monitored between 1st of Jan-2006 and 31th limitations. For instance, if the population under study of Jun-2015. had suffered any unusual health problem during this During this period in both regions, and according to period, the problem could be misclassified as a regular the European Union regulations [47], fallen animals were event and be included as a part of the basal pattern of the removed from farms using specialist disposal services. series, leading to a loss of sensitivity for similar alerts. According to the Royal Decree 728/2007 (Annex III) Fourth, to illustrate how our system forecasts and de- [48], all the cattle farmers were obliged to use national tects unusual aberrations, we used a conventional upper fallen stock companies to collect and dispose their fallen confidence limit of 95% a year and a half ahead. How- stock and communicate the cattle casualties to the offi- ever, considering the unspecific nature of fallen stock cial animal health services. data and the heterogeneity among subpopulations, it is plausible to think that the causes that underlie each Types and sources of data mortality pattern and produced unusual aberrations Two types of data from dairy cattle populations were could differ for each subpopulation. This implies that used: fallen stock and demographic data. Fallen bovines the optimal upper limits may have to be set according to data were provided by two national insurance compan- the context of each subpopulation. Unfortunately, cur- ies: Entidad de Seguros Agrarios (ENESA) and Agrupa- rently, the available information to determine these ción Española de Entidades Aseguradoras de los Seguros limits is scarce and unspecific, and any validation of this Agrarios Combinados S.A. (AGROSEGURO). For every system in the field would require significant additional visit where carcass disposal occurred, the following in- effort. Accordingly, we believe that the validation of this formation was recorded: the identification code of the method as an early warning system would require more farm, date of pick-up, number of carcasses, and number robust and updated data as well as the establishment of of kilograms collected. On the other hand, demographic a more effective hierarchical communication between data yearly updated, were provided by the Sub- the field and each administrative level. All these aspects Directorate General for Animal Health and Hygiene and have led us to plan for a follow-on study to conduct an Traceability of the Ministry of Agriculture, Fisheries and accurate validation using a multidisciplinary approach Food of Spain (MAPA). This data set included: the iden- and more recent data. tification code of the farm, type of production, number of dairy cattle per holding, and administrative aggrega- Conclusions tions (i.e. region, province, and county). Both sets of This work proposes a novel methodology to enhance the data, fallen bovine stock and bovine population, were monitoring of bovine fallen stock data at different ad- merged in a final dataset using the identification code of ministrative levels and detect unusual mortality events. the farm as a primary key. The system may provide essential information to identify According to the legal basis on the protection of personal spatiotemporal populations at high risk and allocate data [49], the non-public staff and scientists involved in the more effectively resources destined to control and pre- analyses of these data signed a confidentiality agreement to vent potential health problems. Moreover, the method- set out the terms and conditions to limit the collection, ology proposed in this work can be applied to other handling, and disclosure of non-public information from populations and adapted to other AHSyS initiatives that farmers. require the follow-up of many populations with dispar- ate husbandry systems. Statistical analysis Methods We built a routine that consisted of a sequence of ana- The populations under surveillance were dairy cattle in lytical methods designed to: two Spanish regions: R1 (Principado de Asturias) and R2 (Catalonia). The region R1 included a large number of (1) describe and select longitudinal fallen stock data farms with small herd sizes, mostly in extensive hus- aggregated at different levels (i.e., region, province, bandry systems. This region R1 had green landscapes and county); and was characterized by coastal oceanic climate with (2) fit, for each subpopulation, appropriate time series moderate to high precipitations and moderate tempera- models; tures. In contrast, the region R2 had fewer farms but (3) forecast the number of fallen bovine stock at n- with larger herd sizes, most of which were in intensive ahead period; and Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 10 of 13 (4) detect and register aberrations in mortality counts noise. These ARIMA models can also be used even if over time, considering the conventional upper the series are non-stationary (i.e., when they present confidence limit of 95%. positive, negative, or quadratic trends, and/or annual and biannual seasonality, among other). In these cases, Hierarchical time series structures some initial differentiation can be applied one (d = 1) or Firstly, for exploring and comparing baseline patterns more times (d > 1) to make the series stationary. Con- from different sub-populations, fallen stock data aggre- ceptually, when series are differentiating, the count at gated at each of the three administrative levels were de- time t-1 is subtracting from the count at time t. The scribed and represented using HTS [13, 14]. HTS general ARIMA model is defined by the parameters p provide a rapid means to visualize, aggregate, and disag- (autoregressive part), d (differentiation), and q (moving gregate series from different subpopulations. The HTS average part). were designed by combining the information on the Let the random variable Xt be an ARMA (p, q) model fallen dairy cattle at the most disaggregated level of the such that: time series (in our case, counties) and the hierarchical organization. This dictates how the information, at the Xt ¼ α þ ρ1 Xt−1 þ ρ2 Xt−2 þ … þ ρp Xt−p þ Zt county level, had to be aggregated into provinces and re- þ θ1 Zt−1 þ θ2 Zt−2 þ … þ θq Zt−q ð1Þ gions. An explanatory illustration has been included in where Xt is a stationary series at time t, α is the intercept the supplemental material (S1 Fig) to explain the meth- of the model, ρ1, ρ2, …, ρp are the coefficients of the odology used to build the HTS structures. The ‘base’ autoregressive term, θ1, θ2, …, θq are the coefficients of and ‘hts’ packages in the software R were used to sup- the moving average term, and Zt, Zt − 1, …, Zt − q are the port their creation [14]. error terms of the model which are normally distributed. In case of lack of stationarity, the series Xt can be differ- ARIMA modelling including trend and seasonal terms entiated (d > 0) before estimating the parameters in Eq.1 The patterns and forecasting of weekly fallen bovines at to avoid possible trend and seasonal effects, leading to different administrative levels were analyzed using an ARIMA (p, d, q) model. Additionally, trends and sea- ARIMA models that included trend and seasonal com- sonal influences can be included in the ARIMA (p, d, q) ponents as covariates. It is important to remark that model as covariates. In particular, we accounted for ARIMA models are preferable for regular time series of trend and seasonality effects in the ARIMA (p, d, q) normally distributed observations. Therefore, once the model by using the following equation: data were disaggregated into regions, provinces, and counties through HTS, each series that showed counts     2πt 2πt greater than or equal to 10 in average and normally dis- Yt ¼ γ0 þ γ1 t þ γ2 sin þ γ3 cos   52   52 tributed [50] were considered for ARIMA modelling. 2πt 2πt Since we were dealing with counts, most of the time þ γ4 sin þ γ5 cos þ Xt ; ð2Þ 26 26 series were Poisson distributed [12]. According to [51], for those series with an average higher than 10 counts, where Yt is the observed series at time t (which can be the Poisson distribution is well approximated by the non-stationary), Xt is the stationary ARIMA (p, d, q) Normal distribution; thus our selected series could be model expressed in Eq. 1., γ1is the trend coefficient, γ2 properly modelled through a classical parametric ap- and γ3 are the coefficients of the annual seasonality ag- proach such ARIMA. These series also coincided with gregated at a weekly level, and γ4 and γ5 the coefficients those administrative levels in which the dairy popula- of the biannual seasonality aggregated at a weekly level. tions were larger within these two regions of study. Here the trigonometric part corresponds to the first and The ARIMA model, an extension of the ARMA (Auto- second-order Fourier terms commonly used in the regressive and Moving Average) model, has been widely analysis of time series [54]. used in the classical time series analysis and has also These models were used to monitor the number of been applied in many contexts related to veterinary and weekly fallen bovines for each subpopulation at every se- public health [7, 50, 52, 53]. The autoregressive part of lected administrative level and identify abnormal peaks of an ARIMA model (AR) indicates that the temporal re- mortality by comparing the predictions to the observations. sponse variable (i.e., weekly fallen bovine counts) is The number of weekly fallen bovines recorded be- regressed by its lags. While the moving average part tween Jan-2006 and Jun-2015 at the regional, provincial, (MA) indicates that the model error is regressed by its and county levels were divided into training and testing lags (i.e., at time t, the lags of a random variable Zt are data sets. The data collected between Jan-2006 and Dec- essentially temporal delays such as Zt − 1, Zt − 2 …). The 2013 were used for estimating the model parameters for latter means that the error does not behave like white each population, while the data collected between Jan- Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 11 of 13 2014 and Jun-2015 were used for testing the models and codes of all involved farms. Additionally, the list also identifying unusual increases of mortality. provided information on the number of carcasses re- To estimate the most appropriate values of parameters corded during the peak and over the preceding two p, d and q for the ARIMA models (Eq. 1), and also the weeks in those farms. linear trend and seasonal coefficients (Eq. 2) we built a The programming routine was built using the base programming routine. This routine explored a range of package of R [28]. values between 0 and 5 for parameters p and q. For in- stance, if p = 5, the observed count at week t was Supplementary information regressed by the counts at weeks t-1, t-2, t-3, t-4 and t-5; Supplementary information accompanies this paper at https://doi.org/10. 1186/s12917-020-02312-8. while if q = 5, the error at week t was regressed by the errors at weeks t-1, t-2, t-3, t-4 and t-5. Here, the par- Additional file 1: S1 Fig Illustration of hierarchical time series structures (hts). ameter d was constrained to a binary indicator, taking 0 Additional file 2: S2 and S3 Figs show the hierarchical time series by (not differentiating) or 1 (differentiating once and re- region, provinces and county in both regions R1 and R2. These Figures moving the linear trend). Note that in those series in are computed by using the hts package included in R. which the trend was not linear, but quadratic, the series Additional file 3: S1 Table. Estimates and standard errors of the selected ARIMA (p,d,q) models for region R1, its province and counties was differentiated, and a coefficient of the linear trend between 2006 and 2013. S2 Table. Estimates and standard errors of the (Eq. 2) was also included in the candidate model. The selected ARIMA(p,d,q) models for region R2, its provinces and counties inclusion of this coefficient was because the trend was between 2006 and 2013. not linear, and hence it was not entirely removed by dif- ferentiating only once. In some scenarios, this problem Abbreviations ACF: Auto-Correlation Function; AGROSEGURO: Agrupación Española de could have been addressed by considering values of d Entidades Aseguradoras de los Seguros Agrarios Combinados S.A; AIC: Akaike greater than 1. Information Criterion; AHSyS: Animal Health Syndromic Surveillance; In order to select the most suitable ARIMA model for AR: Autoregressive; ARIMA: Autoregressive Integrated and Moving Average; ARMA: Autoregressive and Moving Average; BIC: Bayesian Information each time series among all candidates, three criteria were Criterion; C: County; d: Order of differentiation; ENESA: Entidad de Seguros considered. The first one was based on the Bayesian In- Agrarios; HINAR: Hermite Integer-Valued Autoregressive; HMM: Hidden formation Criterion (BIC) proposed by Schwarz [55]. Markov Models; HTS: Hierarchical time series structures; INAR: Integer-valued Autoregressive; MA: Moving Average; MAPA: Ministry of Agriculture, Fisheries The model with the smallest BIC was preferred over the and Food; p: Order of autoregressive part; P: province; PACF: Partial Auto- whole set of candidates. The BIC was used because it Correlation Function; q: Order of moving average part; R: Region typically provides more parsimonious models (in terms of the number of parameters) than those proposed by Acknowledgements The authors of this work want to thank the technical staff of ENESA, the Akaike Information Criterion (AIC). In addition, the AGROSEGURO and REGA for their valuable contribution and also for BIC is also more robust than the AIC when coping with supplying the data. We also thank the technical staff of Los Servicios heterogeneity in samples [56]. The second and third cri- Oficiales del MAPA, the researchers of CReSA-IRTA and, a special thanks to Dr. Joaquim Segales (CReSA-IRTA), Dr. Jordi Casal (CReSA-IRTA) and Dr. teria consisted of assessing the statistical significance of Andres Perez (University of Minnesota) for their support during the elabor- the model parameters given a significance level (i.e., 5%), ation of this work. and checking the behavior of the model residuals Authors’ contributions through the Auto-Correlation Function (ACF) and the A.F.F and A.A.C. designed the system, analyzed and interpreted data. A.F.F. Partial Auto-Correlation Function (PACF), respectively. built the programming routine. A.A.C. designed the maps. P.P. was a major The best ARIMA model was the one for which the third contributor in the performance of statistical modeling techniques. G.C. and L. R participated in the design of the study, the data organization, collection criterion was not only completely satisfied but also dem- and interpretation of result. C.R, F.C.D and J.S. contributed to the discussion onstrated good results for both the first and the second and to the preparation of the manuscript. All authors read and approved the criteria. Full details can be found in [20, 21, 54]. final manuscript. Subsequently, each selected model was used to predict Funding the expected number of carcasses that would be col- The Spanish Ministry of Agriculture, Food and Environment have supported lected weekly between Jan-2014 and Jun-2015 in each this work as a minor contract, whose primary role was as data providers. population from every selected administrative aggrega- Availability of data and materials tion. Those observed counts of fallen dairy cattle that The authors (AFF & AA) were required to analyze these data. Public access to exceeded the 95% confidence limits predicted by the se- these data is only possible through written permission of the Spanish lected ARIMA models were identified as mortality peaks. Ministry of Agriculture and Fisheries and Food The access of programming routines and analyses generated during this study are possible through the Reasonably, any relevant increase in mortality could sign contact of the corresponding author of this paper. an unusual event in the population, and thus require fur- ther investigation at field level to determine the specific Ethics approval and consent to participate This work has been written in agreement with the data providers: causes. In case a presumed mortality peak would be de- AGROSEGURO SA and the Spanish Ministry of Agriculture, Food, and tected, this routine provided a list with the identification Environment. The data have been analyzed aggregately, and guaranteeing Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 12 of 13 the digital rights, their privacy and security under the signed agreement syndromic surveillance. Prev Vet Med. 2012;105:244–52. https://doi.org/10. -following the Organic Law 3/2018, of December 5, on the Protection of 1016/j.prevetmed.2011.12.015. Personal Data- and the Council Regulation (EC) 1069/2009 that sets the 12. Fernández-Fontelo A, Fontdecaba S, Alba A, Puig P. Integer-valued AR health rules as regards animal by-products and non-human-consumption de- processes with Hermite innovations and time-varying parameters: an rived products. This research did not involve any specific clinical study using application to bovine fallen stock surveillance. Stat Model. 2017;17(3):1–24. animal experimentation. 13. Hyndman RJ, Ahmed RA, Athanasopoulos G, Shang HL. Optimal combination forecasts for hierarchical time series. Comput Stat Data Anal. 2011;55:25792589. https://doi.org/10.1016/j.csda.2011.03.006. Consent for publication 14. Hyndman RJ, Ahmed RA, Shang HL, Wang E. hts: Hierarchical and grouped Not Applicable. time series. R package version 4.00, URL: http://CRAN.R- project.org/ package=hts; 2014. Competing interests 15. Serfling R. Methods for current statistical analysis of excess pneumonia- The authors declare that they have no competing interests. influenza deaths. Publ Hlth Rep. 1963;78:494–506. 16. Healy MJR. A simple method for monitoring routine statistics. Statistician. Author details 1983;32:347–9. 1 Chair of Statistics, School of Business and Economics, Humboldt Universität 17. Ngo L, Tager IB, Hadley D. Application of exponential smoothing for zu Berlin, Berlin, Germany. 2Departament de Matemàtiques, Universitat nosocomial infection surveillance. Am J Epidem. 1996;143:637–47. Autònoma de Barcelona, Cerdanyola del Vallès, Barcelona, Spain. 18. Holt CC. Forecasting Seasonals and trends by exponentially weighted 3 Subdirección General de Sanidad e Higiene Animal y Trazabilidad. Ministerio moving averages, Memorandum ONR 52/1957. Pittsburgh: Carnegie de Agricultura y Pesca, Alimentación (MAPA), Madrid, Spain. 4Centre for Institute of Technology; 1957. Veterinary Epidemiological Research, AVC, University Prince Edward Island 19. Winters PR. Forecasting sales by exponentially weighted moving averages. (UPEI), Charlottetown, Canada. 5Department of Computer and Information Management Sci. 1960;6:324–42. Sciences, University of Strathclyde, Glasgow, Scotland, UK. 6Department of 20. Chatfield CH. The analysis of time series: an introduction. 6th ed. Boca Disease Control and Epidemiology, National Veterinary Institute (SVA), Raton: Chapman and Hall/CRC; 2004. Uppsala, Sweden. 7Department of Veterinary Population Medicine, College of 21. Cowpertwait PSP, Metcalfe AV. Introductory time series with R. Heidelberg, Veterinary Medicine, University of Minnesota, St. Paul, USA. 8Centre de London, New York: Springer Dordrecht; 2009. Recerca en Sanitat Animal (CReSA), Institut de Recerca i Tecnologia 22. Moriña D, Puig P, Rios J, Vilella A, Trilla A. A statistical model for hospital Agroalimentàries (IRTA), Cerdanyola del Vallàs, Barcelona, Spain. admissions caused by seasonal diseases. Stat Med. 2011;30(26):3125–36. https://doi.org/10.1002/sim.4336. Received: 16 October 2018 Accepted: 11 March 2020 23. Le Strat Y, Carrat F. Monitoring epidemiologic surveillance data using hidden Markov models. Stat Med. 1999;18:3463–78. 24. Martinez-Beneito MA, Conesa D, Lopez-Quilez A, Lopez-Maside A. Bayesian References Markov switching models for the early detection of influenza epidemics. 1. Hoinville LJ, Alban L, Drewe JA, Gibbens JC, Gustafson L, Hasler B, et al. Statt Med. 2008;27:4455–68. Proposed terms and concepts for describing and evaluating animal-health 25. Unkel S, Farrington CP, Garthwaite PH. Statistical methods for the surveillance systems. Prev Vet Med. 2013;112(1):1–12. https://doi.org/10. prospective detection of infectious disease outbreaks: a review. J R Stat Soc. 1016/j.prevetmed.2013.06.006. 2012;175(1):48–82. 2. Kuiken T, Leighton FA, Fouchier RAM, LeDuc JW, Peiris JSM, Schudel A, et al. 26. Jornet-Sanz M, Ana Corberan-Vallet A, Santonja F, Villanueva R. A Bayesian Public health. Pathogen surveillance in animals. Science. 2005;309(5741): stochastic SIRS model with a vaccination strategy for the analysis of 1680–1. https://doi.org/10.1126/science.1113310. respiratory syncytial virus. SORT. 2017;41(1):159176. 3. Dorea FC, Sanchez J, Revie CW. Veterinary syndromic surveillance: current 27. Stroup DF, Williamson GD, Herndon JL, Karon JM. Detection of aberrations in the initiatives and potential for development. Prev Vet Med. 2011;101(1):1–7. occurrence of notifiable diseases surveillance data. Statist Med. 1989;8:323–9. https://doi.org/10.1016/j.prevetmed.2011.05.004. 28. R Core Team. R: A language and environment for statistical computing. 4. Dorea FC, McEwen BJ, McNab WB, Revie CW, Sanchez J. Syndromic Vienna: R Foundation for Statistical Computing; 2015. URL https://www.R- surveillance using veterinary laboratory data: data pre-processing and project.org/. algorithm performance evaluation. J R Soc Interface. 2013;10(83):20130114. 29. Parker RA. Analysis of surveillance data with Poisson regression: a case https://doi.org/10.1098/rsif.2013.0114. study. Statist Med. 1989;8:285–94. 5. Dupuy C, Bronne A, Watson E, Wuyckhuise-Sjouke L, Reist M, Fouillet A, 30. Jackson MJ, Baer A, Painter I, Duchin J. A simulation study comparing et al. Inventory of veterinary syndromic surveillance initiatives in Europe aberration detection algorithms for syndromic surveillance. BMC Med (triple-S project): current situation and perspectives. Prev Vet Med. 2013;111: Informat Decsn Makng. 2007;7:6. 220–9. https://doi.org/10.1016/j.prevetmed.2013.06.005. 31. Farrington CP, Andrews NJ, Beale AD, Catchpole MA. A statistical algorithm 6. Dorea FC, Vial F. Animal health syndromic surveillance: a systematic for the early detection of outbreaks of infectious disease. J R Statist Soc. literature review of the progress in the last 5 years (2011-2016). VMRR. 2013; 1996;159:547–63. 7:157170. 32. Molenberghs G, Verbeke G, Demétrio CGB. Hierarchical models with normal 7. Alba A, Dorea FC, Arinero L, Sanchez J, Cordon R, Puig P, et al. Exploring the and conjugate random effects: a review. SORT. 2017;41(2):191–254. surveillance potential of mortality data: nine years of bovine fallen stock 33. Stern L, Lightfoot D. Automated outbreak detection: a quantitative data collected in Catalonia. PLoS One. 2015a;10(4):e0122547. https://doi.org/ retrospective analysis. Epidem Infectn. 1999;122:103–10. 10.1371/journal.pone.0122547. 34. Wieland SC, Brownstein JS, Berger B, Mandl KD. Automated real-time 8. Alba A, Fernández-Fontelo A, Revie CW, Dorea F, Sanchez J, Romero L, et al. constant specificity surveillance for disease outbreaks. BMC Med. Informat. Development of new strategies to model bovine fallen stock data from DesnMakng. 2007;7:15. large and small subpopulations for syndromic surveillance use. Epidemiol et 35. Shewart, W.A. (1931) Nostr and Reinhold. Roberts, S. W. (1959) Control chart sante anim. 2015b;67:67–76. tests based on geometric moving averages. Technometrics, 1, 239–250. 9. Morignat E, Perrin JB, Gay E, Vinard JL, Calavas D, Henaux V. Assessment of 36. Page ES. Continuous inspection schemes. Biometrika. 1954;41:100–15. the impact of the 2003 and 2006 heat waves on cattle mortality in France. 37. Naus J, Wallenstein S. Temporal surveillance using scan statistics. Statist PLoS One. 2014;9(3):e93176. https://doi.org/10.1371/journal.pone.0093176. Med. 2006;25:311–24. 10. Perrin JB, Ducrot C, Vinard JL, Morignat E, Gauffier A, Calavas D, Hendrikx P. 38. Sego LH, Woodall WH, Reynolds MR. A comparison of surveillance methods Using the National Cattle Register to estimate the excess mortality during for small incidence rates. Statist Med. 2008;27:1225–47. an epidemic: application to an outbreak of bluetongue serotype 8. 39. Rogerson PA. Surveillance systems for monitoring the development of Epidemics. 2010;2(4):207–14. spatial patterns. Statist Med. 1997;26:2081–93. 11. Perrin JB, Ducrot C, Vinard JL, Morignat E, Calavas D, Hendrikx P. Assessment 40. Rogerson PA. Monitoring point patterns for the development of spacetime of the utility of routinely collected cattle census and disposal data for clusters. J R Statist Soc. 2001;164:87–96. Fernández-Fontelo et al. BMC Veterinary Research (2020) 16:110 Page 13 of 13 41. Wallenstein S, Gould M, Kleinman M. Use of the scan statistic to detect time-space clustering. Am J Epidem. 1989;130:1057–64. 42. Lawson AB. Spatial and spatio-temporal disease analysis. In: Lawson AB, Kleinman K, editors. Spatial and Syndromic surveillance for public health. Chichester: Wiley; 2005. p. 55–76. 43. Clark AB, Lawson AB. Surveillance of individual level disease maps. Statist Meth Med Res. 2006;15:353–62. 44. McCabe GJ, Greenhalgh D, Gettingby G, Holmes E, Cowden J. Prediction of infectious diseases: an exception reporting system. J Med Informat Technol. 2003;5:67–74. 45. Kulldorff M. A spatial scan statistic. Communs Statist Theor Meth. 1997;26: 1481–96. 46. Jolliffe IT. Principal component analysis. 2nd ed. New York: Springer; 2002. 47. European Union Commission, et al. Regulation (EC) no. 1774/2002 of the European Parliament and of the Council of 3 October 2002 laying down health rules concerning animal by-products not intended for human consumption. Off J Eur Communities L. 2002;273:1–95. 48. Royal Decree 728 / 2007, Of 13 Of June, By Which Is Sets And Regulates The Register General Of Movements Of Won And The Register General Of Identification Individual Of Animals. (Real Decreto 728/2007, de 13 de junio, por el que se establece y regula el Registro general de movimientos de ganado y el Registro general de identificación individual de animales). 49. Royal Decree 1720/2007, of December 21th, which approves the regulation implementing Organic Law 15/1999 of December 13th on the protection of personal data in any type of clinical research. (Real Decreto por el que se aprueba el Reglamento de desarrollo de la Ley Orgánica 15/1999, de protección de datos de carácter personal, BOE 17 of January 19th, 2008). 50. Lee HS, Her M, Levine M, Moore GE. Time series analysis of human and bovine brucellosis in South Korea from 2005 to 2010. Prev Vet Med. 2013; 110(2):190–7. https://doi.org/10.1016/j.prevetmed.2012.12.003. 51. Rosner B. Fundamentals of biostatistics. 8th ed. Boston: Cengage Learning; 2016. 52. Neumann E, Hall W, Stevenson M, Morris R. Descriptive and temporal analysis of post-mortem lesions recorded in slaughtered pigs in New Zealand from 2000 to 2010. N Z Vet J. 2014;62(3):110–6. 53. Chamba Pardo FO, Alba-Casals A, Nerem J, Morrison RB, Puig P, Torremorell M. Influenza herd-level prevalence and seasonality in breed-to-wean pig farms in the Midwestern United States. Front Vet Sci. 2017;4:167. 54. Brockwell PJ, Davids RA. Introduction to time series and forecasting. 2nd ed. New York: Springer-Verlag New York, Inc; 2002. 55. Schwarz G. Estimating the dimension of a model. AnnStatist. 1978;6(2): 461464. 56. Brewer MJ, Butler A, Cooksley SL. The relative performance of AIC, AICC and BIC in the presence of unobserved heterogeneity. MEE. 2016;7(6):679–92. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.