 Research
 Open Access
 Published:
Forecasting United States heartworm Dirofilaria immitis prevalence in dogs
Parasites & Vectorsvolume 9, Article number: 540 (2016)
Abstract
Background
This paper forecasts next year’s canine heartworm prevalence in the United States from 16 climate, geographic and societal factors. The forecast’s construction and an assessment of its performance are described.
Methods
The forecast is based on a spatialtemporal conditional autoregressive model fitted to over 31 million antigen heartworm tests conducted in the 48 contiguous United States during 2011–2015. The forecast uses countylevel data on 16 predictive factors, including temperature, precipitation, median household income, local forest and surface water coverage, and presence/absence of eight mosquito species. Nonstatic factors are extrapolated into the forthcoming year with various statistical methods. The fitted model and factor extrapolations are used to estimate next year’s regional prevalence.
Results
The correlation between the observed and modelestimated countybycounty heartworm prevalence for the 5year period 2011–2015 is 0.727, demonstrating reasonable model accuracy. The correlation between 2015 observed and forecasted countybycounty heartworm prevalence is 0.940, demonstrating significant skill and showing that heartworm prevalence can be forecasted reasonably accurately.
Conclusions
The forecast presented herein can a priori alert veterinarians to areas expected to see higher than normal heartworm activity. The proposed methods may prove useful for forecasting other diseases.
Background
Heartworm disease, caused by the mosquitoborne filarial nematode Dirofilaria immitis, is arguably the most medically important parasitic infection of domestic dogs in the United States (US), affecting at least 115,000 dogs in 2015. Beyond the US, heartworm disease is a global veterinary healthcare problem, with D. immitis affecting dogs in many parts of South America, Europe, Asia, and Australia [1, 2]. Infection is associated with lifethreatening complications and significant financial burden, costing millions in veterinary care annually for disease treatment [3–7]. Although less common and less studied, heartworm disease is also a health concern for other mammals such as domestic cats, domestic ferrets, and some wildlife species [8]. Clinical signs of heartworm disease in domestic dogs include exercise intolerance, coughing, dyspnea, cachexia, anorexia, epistaxis and ascites. Dogs with a high burden of adult heartworms can suffer from pulmonary arterial occlusion and inflammation, leading to pulmonary hypertension and potentially right heart failure. Cats and ferrets may experience similar signs but acute death may occur, even with very low worm burdens. Humans can also be infected with D. immitis, but infections are rare, with fewer than 100 cases reported in the US over the last 60 years [9]. Human infection is most commonly asymptomatic, with people considered deadend hosts for the parasite. While rare, human D. immitis infection is highly problematic in that it most often manifests as “coin lesions” in the lungs that may be mistaken for a neoplasm on chest radiographs; surgical excision is necessary to differentiate the two entities [10].
Heartworm disease in dogs is most commonly diagnosed through the detection of circulating D. immitis antigen in the blood [3, 11]. The prevalence of heartworm infection in the US varies considerably by geographical region. Two nationwide surveillance studies of D. immitis infection seroprevalence (henceforth prevalence) in domestic dogs found the highest prevalence in the Southeast and the lowest in the Northeast [11]. For unknown reasons, a noted decrease in the prevalence of D. immitis occurred between the 2001–2007 and 2010–2012 in these studies. Importantly, regardless of time period and even within areas where heartworm infection is considered common, there can be considerable local variation, with prevalence reaching as high as 13 % [12, 13].
Numerous factors are purported to be associated with regional and local variations in D. immitis prevalence in domestic dogs. Highly effective commercially available anthelmintics (e.g. macrocyclic lactones (ML) [3], including the avermectins (ivermectin, selamectin) and the milbemycins (moxidectin, milbemycin oxime) can be administered monthly to prevent the development of immature stages into adult worms. Yearround preventive use is recommended throughout the US, yet the majority of dogs only receive seasonal treatment [14]. Even within highly endemic regions, anthelminthic use varies based on client compliance, knowledge, or dog owner’s demographics. In addition, resistance of D. immitis to ML has been recently documented and is a growing concern in the Gulf States, but the current extent of resistant phenotypes remains unknown [7, 15]. Dirofilaria immitis can be transmitted by over 70 species of mosquitoes, although certain species (e.g. Aedes trivittatus, Aedes sierrensis and Culex quinquefasciatus) are considered more important vectors [16]. Because the density of mosquitoes and community composition of competent vector species is influenced greatly by habitat use and climate, these factors should be considered when investigating factors influencing heartworm disease. In support of this, a previous study found that temperature, median household income, population density, precipitation, elevation, relative humidity, forestation coverage, and surface water coverage all significantly influence D. immitis prevalence in dogs [16].
Clearly, it would be advantageous to accurately forecast D. immitis prevalence on a local scale, providing an a priori alert to veterinarians in problem areas where immediate remediation measures could be taken. Annual forecasts of emergent infection will also inform veterinary and public health officials to shifting areas of infection, particularly in temperate regions of the US where D. immitis is generally absent, rare, or prevalence is highly influenced by annual variation in biotic or abiotic factors.
Methods
Data structure
The data studied here contain 31,345,244 heartworm antigen test results from dogs in the conterminous United States from 2011 to 2015, and various climate, geographic and socioeconomic factors purported to influence heartworm prevalence. The raw tests were obtained from the Antech and IDEXX laboratories [17, 18]. Over all 5 years in the study, 384,905 of the tests were positive (1.23 %). The test data contain the county/parish of the testing clinic and the month when the tests were conducted; however, no measure of uncertainty is given with the individual test results.
The test data were aggregated into the number of positive and negative tests for each year in each conterminous United States county/parish. Table 1 lists 16 explanatory factors that are purportedly related to dog heartworm prevalence, as well as their time period of record and geographic scale of collection. These 16 factors include the climatic variables of annual temperature, precipitation and relative humidity, the geographic variables of county elevation, forestation coverage and surface water coverage, the socioeconomic variables of county population density and median household income, and the presence of eight mosquito vectors. For more details on these factors, see [19].
Figure 1 displays countylevel raw heartworm prevalences obtained by dividing the number of positive tests by the number of tests over all 5 years in the study. The raw prevalences exhibit a large degree of spatial correlation in that neighboring counties tend to report similar prevalences. Significant temporal dependence is also present in the data: the current prevalence is similar to past prevalence. Therefore, this data set requires a statistical model with both spatial and temporal dependence.
Figure 2 provides a spatially smoothed prevalence map, using a headbanging smoothing procedure, based on the Fig. 1 prevalences. In the headbanging smoothing procedure, 45 triples were employed. The smoothing was also weighted proportionally to the number of tests in each county over the 5year period. This prevents the map from signaling a high/low prevalence that is more likely attributed to a small sample (one positive out of three tests has the same prevalence as one hundred positive in 300 tests, though the latter is more indicative of infection risk). Details on headbanging smoothing and its uses in disease mapping are contained in [16]. Figure 2 serves as a contemporary depiction of the “baseline” heartworm risk for dogs in the United States.
Statistical modeling
The model and methods used to statistically analyze the heartworm tests are now described. The goal here is to assess the significance of the 16 factors and accurately estimate regional heartworm prevalence.
Let Y _{ s }(t) and n _{ s }(t) denote the number of positive and total tests conducted in county s during year t, respectively, for counties s ∈ {1, …, S} and years t ∈ {1, …, T}. Methodologies for modeling spatial and temporal dependence have received much recent attention in the statistics literature [20–27]. Among many choices, Bayesian hierarchical models have been prominent due to their flexibility. In a Bayesian hierarchical model, spatial and temporal dependence is modelled in a hierarchy via a series of random effect terms with prescribed structures; see [20, 21] for a modern review of spatiotemporal models. Typically, when modeling the spatial or spatialtemporal dependent count data via parametric models, a Poisson distribution is preferred [21–24]. The following hierarchical regression model is used here:
where log(⋅) denotes natural logarithm, X _{ s }(t) = (X _{ s,1}(t), …, X _{ s,16}(t)) ' is a vector of covariate information for county s at time t ('denotes matrix transpose), β = (β _{0}, …, β _{ p }) ' is a vector of regression coefficients, p _{ s }(t) denotes the heartworm prevalence of county s at time t, the symbol ~ means has the distributional type, and ∣ indicates given quantities. Equation (1) indicates that Y _{ s }(t) has a Poisson distribution with mean n _{ s }(t)p _{ s }(t). In addition, it is assumed that the positive test counts (i.e. Y _{ s }(t)) are conditionally independent of each other given the number of tests, factor information, and random effects. This does not imply that Y _{ s }(t) is independent across varying space s or time t.
To relate prevalence to the factors and build spatial and temporal dependence, the model in (2) is proposed, as is common in Poisson regressions [21–24]. There are many ways to induce spatial and temporal dependence from the random effects {ξ _{ s }(t)}. One natural and popular choice is the conditional autoregressive (CAR) structure
where ξ _{ t } = (ξ _{1}(t), …, ξ _{ S }(t)) ' and ϕ _{ t } = (ϕ _{1}(t), …, ϕ _{ S }(t)) ' are random vectors. Equation (4) indicates that the spatial effects (i.e. ϕ _{ t }, for t = 1,…,T) are independent and identically distributed random vectors that follow a conditional autoregressive (CAR) model [25], which is a popular choice for modeling spatial dependence [26].
More specifically, let ϕ = (ϕ _{1}, …, ϕ _{ S }) ' denote a random vector which follows a CAR model. There are several varieties of CAR models. Typically, the CAR model is specified via a series of univariate conditional distributions. Spatial dependence is induced through a neighboring system involving geographically adjacent counties. Our version of the CAR model, taken from [26], uses
Here, ϕ _{− k } = (ϕ _{1}, …, ϕ _{ k − 1}, ϕ _{ k + 1}, …, ϕ _{ S }) ' is a vector that contains county effects for all counties except the kth one and N(μ, σ ^{2}) denotes a normally distributed quantity with mean μ and variance σ ^{2}. In addition, W = {w _{ k,i }} is an S × S dimensional matrix that describes the neighborhood structure of all counties. Specifically, the entries of W are either zero or unity; w _{ k,i } = 1 if and only if the ith and kth counties share some common border.
The parameter τ ^{2} in (5) is a scaling variance parameter. In fact, it can be seen that the conditional variance of ϕ _{ k } given its neighbor’s random effects is inversely proportional to the number of counties bordering county k. Hence, counties with more neighboring counties tend to have a smaller variance, which is reasonable since data from the bordering counties helps predict the prevalence in the said county.
In Eq. (5), ρ ∈ [0, 1] is an autocorrelation parameter that governs correlation between bordering counties. Notice that the conditional expectation of ϕ _{ k } is the weighted arithmetic average of the neighboring random effects, multiplied by ρ. When ρ = 0, the conditional expectations of ϕ _{ s } are zero and all random effects are independent of each other; antipodally, ρ close to unity indicates strong spatial dependence between bordering counties.
In (3), timedependence is modeled through a temporal autoregressive model of order one (AR(1)), which is a time series staple [28]. Here, it describes prevalence for a fixed county across different years. The parameter φ is the temporal correlation between consecutive years and lies within (−1, 1). This ensures a causal and stationarity solution to the time series model [28], which is needed in estimation.
From (5), it is possible to explicitly identify the joint distribution of ϕ, which is multivariate normal:
where W is the previously mentioned neighborhood matrix and D = {d _{ i,j }} is an S × S diagonal matrix whose i th diagonal element is the number of neighboring counties for county i.
Bayesian techniques are used to estimate the model parameters, which are β, φ, ρ, and τ ^{2}. Thus, to complete the Bayesian hierarchical model, the following prior distributions for these parameters are introduced:
Prior distributions for φ and ρ are taken as uninformative in that all admissible possibilities are equally likely. Priors for the regression coefficients β _{0}, …, β _{16} are taken as diffuse so that inferences for these parameters are based primarily on the data. The prior for τ ^{− 2} is chosen as a conjugate prior (the posterior and prior distributions are from the same distributional family) for ease of computation. The random effects and model parameters are estimated based on posterior samples from a Markov chain Monte Carlo (MCMC) simulation. The MCMC simulation for our model uses a combination of Gibbs and MetropolisHastings steps. In the implementation of the algorithm, the test data for nonreporting counties was viewed as being latent and was subsequently sampled along with the model parameters. To run our MCMC simulation and assess significance of model parameters, a program was developed and implemented in R and C++.
Results
Model assessment
The spatiotemporal Poisson regression model in (1) has 16 explanatory factors, all of which may not have predictive power. To assess this issue, a full model with all 16 factors was first fitted. Credible intervals, Bayesian analogs to confidence intervals in frequentist statistics, were then created for the parameters of interest. Table 1 summarizes our full model findings, showing 16 regression coefficients estimates (posterior median) and their 95 % highest posterior density (HPD) intervals; for further details about credible and HPD intervals, see [29, 30].
Table 2 implies that not all factors are significant, e.g. 95 % HPD intervals of annual precipitation, elevation, percentage surface water coverage, and all mosquito presence factors except A. albopictus contain zeroes. To develop a parsimonious model with only significant factors, all explanatory variables whose 95 % HPD intervals contain zeroes were removed and the model was refitted. This leaves a “reduced model” with the six explanatory factors: annual temperature, annual relative humidity, percentage forest coverage, population density, median household income and A. albopictus absence/presence. Parameter estimates (posterior median) and 95 % HPD intervals for the regression parameters for the reduced model are shown in Table 3. The estimates (posterior median) of the other model parameters are φ = 0.914, ρ = 0.998, and τ ^{2} = 0.802.
Most of the significant factors have an intuitive interpretation. For example, the positive regression coefficient for the temperature and relative humidity factors implies that heartworm is more prevalent in warmer and humid locations. On the other hand, as is seen by negative regression coefficients, heartworm prevalence decreases with increasing population densities and median household incomes. Given the presence of relative humidity in the model, it is not overly surprising that precipitation drops out of the model fit. Finally, the negative regression coefficient on A. albopictus presence is not a contradiction: in the presence of all other factors (which include space and time prevalence histories), presence of this mosquito is associated with lessened heartworm prevalence. It is worthwhile to note that in a separate analysis (results not shown) a model with only A. albopictus was fitted, and the accompanying regression coefficient was nonnegative.
To assess the overall performance of our reduced model, Fig. 3 graphically portrays our fitted model by plotting the average (over all 5 years) of modelestimated prevalence for each county after smoothing (standard Kriging with default parameters were used here). The modelestimated prevalence in Fig. 3 compares well to the headbanging smoothed baseline in Fig. 2. In fact, the correlation between the Figs 2 and 3 graphics is 0.727 (only counties reporting at least one test during the 5 year study period were used in this calculation). Clarifying further, our correlation between the two observation sets {A _{ s }} _{ s = 1} ^{S} and {B _{ s }} _{ s = 1} ^{S} is
where
are the samplesize weighted averages of {A _{ s }} _{ s = 1} ^{S} and {B _{ s }} _{ s = 1} ^{S} , and n _{ s } is the number of tests conducted in county s. Since the correlation here is between smoothed and modelestimated prevalence (these are non sample size dependent quantities), the weights were taken as n _{ s }≡1 (and not the countybycounty sample sizes). The 0.727 correlation achieved indicates that the regression model has explained most of the data structure.
The fitted model has a number of uses. In the next section, it is used to construct annual heartworm prevalence forecasts. The model could also be used to estimate how climate change could alter heartworm disease risk.
Forecasting
This section shows how to use our model to forecast next year’s regional heartworm prevalence. For this, all six significant explanatory factors and the spatialtemporal effects will need to be forecasted for the forthcoming year. To see how our forecast performs, the 2015 test and factor data was removed from the analysis, and the proposed sixfactor model was refitted using data from 2011 to 2014 only. Our forecasts simply “plug in” 2015 forecasted factors for A. albopictus, annual temperature, annual relative humidity, percent forest coverage, population density, median household income, and a randomly generated random effect component into our model; for further details see [29, 30].
Two of the six factors are relatively stable over time: county forestation and the presence of A. albopictus. For these two factors, the most recent observations are used as 2015’s forecasted values.
To forecast annual temperature, historical temperature records were collected from 1895 to 2014 for each county and modeled as an autoregressive model of order one. The AR(1) model for an annual temperature series {F _{ t }} (previously denoted by {X _{ s,1}(t)} in Section 3 for county s) obeys the difference equation
where {ω _{ t }} is zero mean white noise; for further time series forecasting information, see [28]. The AR(1) model can be fitted to the temperature observations using practically any statistical software package. Let \( \widehat{\delta} \) and \( \widehat{\gamma} \) denote estimates of δ and γ, respectively. A prediction of the annual temperature at year t + 1 from temperatures from year 1 to year t is
In our forecast, \( {\widehat{F}}_{t+1} \) is used as next year’s forecasted temperature factor. Figures 4 and 5 compare forecasted and observed annual temperatures for 2015. The correlation between these two figures in (6) is r = 0.996, which is very good (n _{ s }≡1 here).
A simple linear regression was used to forecast next year’s relative humidity (previously denoted by {X _{ s,3}(t)}) and median household income (previously denoted by {X _{ s,8}(t)}) in each county. Historical relative humidities from 2006 to 2014 and median household incomes from 1997 to 2014 were used to fit a regression model of form
for each county. Here, {I _{ t }} denotes the relative humidity (previously denoted by {X _{ s,3}(t)}) or median household income (previously denoted by {X _{ s,8}(t)}), {η _{ t }} is zeromean random noise. Least squares estimators of α and κ, denoted by α and \( \widehat{\kappa} \), respectively, were computed from the data at each county. The forecasted value for year t + 1 is simply
Forecasting the county population density for next year requires the county areas and their recent population counts. The US Census provides good county population estimates for 2010, but not in years since 2010. Estimated state populations were obtained for each state between 1969 and 2014. A simple linear regression was fitted to these data for each state and 2015 state populations were forecasted. This forecasted state population was then partitioned to the counties within the state at a proportion that agrees with 2010 Census proportions.
To forecast the next year’s spatial and temporal random effects, formula (3) is applied. Since the ϕ _{ t } s are independent and identically distributed over various years, based on the values of τ ^{2} and ρ (available from the posterior samples), for the current time t, ϕ _{ t + 1} can be generated from the multivariate normal distribution N(0, τ ^{2}(D − ρ W)^{− 1}). Then ξ _{ t + 1} is calculated by ξ _{ t + 1} = φξ _{ t } + ϕ _{ t + 1}. This process is repeated for each value of (ρ, τ ^{2}, φ), which are available from the posterior sample, thus yielding predictive posterior samples of the next year’s random effect ξ _{ t + 1}. See [29, 30] for additional detail on obtaining predictive posterior samples.
Figures 6 and 7 compare observed and forecasted heartworm prevalence during 2015. One can discern where heartworm is forecasted to be higher/lower than normal by comparing Figs. 2 and 7. The correlation between the Figs. 6 and 7 prevalence, as measured in (6), is 0.940. Hence, the model is accurately forecasting in locations that report more tests. Finally, Fig. 8 presents our heartworm forecast for 2016. When 2016 concludes, we will be able to compare this forecast to 2016 test results.
Discussion
Generally, the management of emerging infectious diseases is approached reactively, with efforts focused on managing outbreaks after onset. The ability to reliably forecast transmission risk, particularly for diseases influenced by dynamic factors such as climate, could shift our paradigm from reaction to prevention. This is particularly true for vectorborne diseases, as specific environmental needs for vector survival are well documented [31]. One approach to infectious disease modeling is to use these factors to predict transmission and model the data in both space and time. This has been used successfully to estimate the incidence of malaria during eradication campaigns in Namibia and cutaneous leishmaniasis in highrisk areas of Columbia [32, 33].
Although preventable, heartworm disease is a relatively common and serious vectorborne disease of domestic dogs. Annual disease incidence, as reported by IDEXX and Antech, averages greater than 100,000 new cases annually. Annual data likely represent the true annual incidence of heartworm infection in domestic dogs: when diagnosed with heartworm, most dogs are either treated, or in some cases euthanized, due to poor outcome or financial constraints [34]. While fulminant infection with D. immitis may be due to lack of owner compliance in use of preventatives, it also may be due to misunderstanding the disease risk; mosquito vectors are known to be dynamic in their range and survival under changing climatic conditions.
To enhance veterinary client education and illuminate the benefits of preventatives, factors associated with D. immitis transmission [16, 35] were identified and used to develop a spatialtemporal conditional autoregressive forecast model of heartworm prevalence. A comparison of observed versus forecasted heartworm prevalence was made in 2015 and was quite accurate. This may be attributed, in part, to the fact that many of the factors influencing heartworm prevalence do not change significantly from year to year (e.g. forest coverage, population density, household income). While temperature and humidity change annually and are important disease risk factors [36], these factors are still reasonably predictable; however, environmental or climate catastrophes (e.g. regional climate shifts, flooding, hurricanes) could impact heartworm incidence. Finally, mosquito populations can fluctuate greatly from year to year as they depend on numerous local landuse and environmental factors. Some competent vectors of D. immitis are still expanding their range in the US [37].
Several of the mosquito presence/absence factors were not included in our final model; this may be because our currently available data are only presence/absence, whereas mosquito abundance, a purportedly more powerful factor, varies annually at a local level. More accurate mosquito counts would likely yield more accurate forecasts. In addition, human activities such as treatment abatement programs may impact mosquito abundances. Since the introduction of West Nile virus into the US, localities have developed or expanded mosquito control programs, including reducing breeding habitats and application of pesticides. With increased concern over the Chikungunya and Zika viruses, it is possible that increased mosquito control may be initiated in the coming year(s). If such programs are initiated, mosquito abundance counts should take these programs into account.
Conclusion
In conclusion, our 2016 heartworm disease forecast (Fig. 8) has some noteworthy implications for veterinary practitioners, including an increased prevalence in northern California, eastern Montana, and central New Mexico. A relatively small increase in risk is predicted in some areas where heartworm is likely under appreciated, such as parts of the Dakotas and Nebraska. Importantly, our data indicate that all regions of the lower 48 United States have some risk for heartworm infection. Our maps and forecasts provide veterinarians with evidencebased recommendations for use of preventive in nonendemic regions of the US and support the recommendation of yearround use of preventive in high risk areas. Ultimately, we believe that these methods can be used to forecast multiple vectorborne diseases with veterinary and human health impacts, including Lyme disease, ehrlichiosis and anaplasmosis. Currently, the Companion Animal Parasite Council (CAPC, www.capcvet.org) provides monthly updates of heartworm prevalence on a county by county scale. Through a combination of realtime updates and forecasting efforts, we hope to see fewer cases of heartworm disease in dogs and cats in the future.
Abbreviations
 AR(1):

First order temporal autoregression
 CAPC:

Companion Animal Parasite Council
 CAR:

Conditional autoregressive model
 HPD:

Highest posterior density
 MCMC:

Markov chain Monte Carlo
 ML:

Macrocyclic lactones
References
 1.
Morchon R, Carreton E, GonzálezMiguel J, MelladoHernández I. Heartworm disease (Dirofilaria immitis) and their vectors in Europe  new distribution trends. Front Physiol. 2012;3:196.
 2.
Simon F, Morchon R, GonzalezMiguel J, MarcosAtxutegi C, SilesLucas M. What is new about animal and human dirofilariosis? Trends Parasitol. 2009;25:404–9.
 3.
Bowman DD, Atkins CE. Heartworm biology, treatment, and control. Vet Clin North Am Small Anim Pract. 2009;39:1127–58.
 4.
Bowman DD, Georgi JR. Georgis’ Parasitology for Veterinarians. 9th ed. Saint Louis: Saunders/Elsevier; 2009.
 5.
Colby KN, Levy JK, Dunn KF, Michaud RI. Diagnostic, treatment, and prevention protocols for canine heartworm infection in animal sheltering agencies. Vet Parasitol. 2011;176:333–41.
 6.
Maxwell E, Ryan K, Reynolds C, Pariaut R. Outcome of a heartworm treatment protocol in dogs presenting to Louisiana State University from 2008 to 2011: 50 cases. Vet Parasitol. 2014;206:71–7.
 7.
Wolstenholme AJ, Evans CC, Jimenez PD, Moorhead AR. The emergence of macrocyclic lactone resistance in the canine heartworm, Dirofilaria immitis. Parasitol. 2015;142:1249–59.
 8.
Venco L, Marchesotti F, Manzocchi S. Feline heartworm disease: A ‘Rubik’scubelike’ diagnostic and therapeutic challenge. J Vet Cardiol. 2015;17 Suppl 1:S190–201.
 9.
Lee AC, Montgomery SP, Theis JH, Blagburn BL, Eberhard ML. Public health issues concerning the widespread distribution of canine heartworm disease. Trends Parasitol. 2010;26:168–73.
 10.
Orihel TC, Eberhard ML. Zoonotic filariasis. Clin Microbiol Rev. 1998;11:366–81.
 11.
Little SE, Beall MJ, Bowman DD, Chandrashekar R, Stamaris J. Canine infection with Dirofilaria immitis, Borrelia burgdorferi, Anaplasma spp., and Ehrlichia spp. in the United States, 2010–2012. Parasit Vectors. 2014;7:257.
 12.
Levy JK, Lappin MR, Glaser AL, Birkenheuer AJ, Anderson TC, Edinboro CH. Prevalence of infectious diseases in cats and dogs rescued following Hurricane Katrina. J Am Vet Med Assoc. 2011;238:311–7.
 13.
Yabsley MJ, DresdenOsborne C, Pirkle EA, Kirven JM, Noblet GP. Filarial worm infections in shelter dogs and cats from northwestern South Carolina, USA. Comp Parisotol. 2004;71:154–7.
 14.
Gates MC, Nolan TJ. Factors influencing heartworm, flea, and tick preventative use in patients presenting to a veterinary teaching hospital. Prev Vet Med. 2010;93:193–200.
 15.
Bourguinat C, Lee AC, Lizundia R, Blagburn BL, Liotta JL, Kraus MS, et al. Macrocyclic lactone resistance in Dirofilaria immitis: Failure of heartworm preventives and investigation of genetic markers for resistance. Vet Parasitol. 2015;210:167–78.
 16.
Wang D, Bowman DD, Brown HE, Harrington LE, Kaufman PE, McKay T, et al. Factors influencing US Canine heartworm (Dirofilaria immitis) prevalence. Parasit Vectors. 2014;7:264–82.
 17.
Antech Diagnostics, Incorporated: http://www.antechdiagnostics.com/Main/Home.aspx
 18.
IDEXX Laboratories, Incorporated: http://www.idexx.com/view/xhtml/en_us/corporate/home.jsf
 19.
Stich RW, Blagburn BL, Bowman DD, Carpenter C, Cortinas MR, Ewing SA, et al. Quantitative factors proposed to influence the prevalence of canine tickborne agents in the United States. Parasit Vectors. 2014;7:417–24.
 20.
Cressie N, Wikle CK. Statistics for Spatiotemporal data. New York: John Wiley & Sons; 2011.
 21.
LópezQuılez A, Munoz F. Review of spatiotemporal models for disease mapping. In: Final report for the EUROHEIS 2 project. 2009. http://www.uv.es/~famarmu/doc/Euroheis2report.pdf.
 22.
MartínezBeneito MA, LópezQuilez A, BotellaRocamora P. An autoregressive approach to spatiotemporal disease mapping. Stat Med. 2008;27:2874–89.
 23.
Xia H, Carlin BP. Spatiotemporal models with errors in covariates: mapping Ohio lung cancer mortality. Stat Med. 1998;17:2025–43.
 24.
Nobre AA, Schmidt AM, Lopes HF. Spatiotemporal models for mapping the incidence of malaria in Pará. Environmetrics. 2005;16:291–304.
 25.
Besag J. Spatial interaction and the statistical analysis of lattice systems. J R Stat Soc Series B. 1974;36:192–236.
 26.
Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data. 2nd ed. Boca Raton: Chapman and Hall/CRC; 2014.
 27.
Brockwell PJ, Davis RA. Introduction to time series and forecasting. 2nd ed. New York: Springer; 2002.
 28.
Harrison J, West M. Bayesian forecasting & dynamic models. New York: Springer; 1999.
 29.
Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. Bayesian data analysis. 3rd ed. Boca Raton: Chapman & Hall/CRC; 2014.
 30.
Nene V, Wortman JR, Lawson D, Haas B, Kodira C, Tu ZJ, et al. Genome sequence of Aedes aegypti, a major arbovirus vector. Science. 2007;316(5832):1718–23.
 31.
Parham PE, Waldock J, Christophides GK, Hemming D, Agusto F, Evans KJ, et al. Climate, environmental and socioeconomic change: weighing up the balance in vectorborne disease transmission. Philos Trans R Soc Lond B Biol Sci. 2015;370(1665). doi:10.1098/rstb.2013.0551.
 32.
Alegana VA, Wright JA, Nahzat SM, Butt W, Sediqi AW, Habib N, et al. Modelling the incidence of Plasmodium vivax and Plasmodium falciparum malaria in Afghanistan 2006–2009. PLoS One. 2014;9(7), e102304.
 33.
PerezFlorez M, Ocampo CB, ValderramaArdila C, Alexander N. Spatial modeling of cutaneous leishmaniasis in the Andean region of Colombia. Mem Inst Oswaldo Cruz. 2016. doi:10.1590/00742760160074.
 34.
Polak KC, SmithBlackmore M. Animal shelters: managing heartworms in resourcescarce environments. Vet Parasitol. 2014;206(1–2):78–82.
 35.
Brown HE, Harrington LC, Kaufman PE, McKay T, Bowman DD, Nelson CT, et al. Key factors influencing canine heartworm, Dirofilaria immitis, in the United States. Parasit Vectors. 2012;5:245.
 36.
Beugnet F, ChalvetMonfray K, Loukos H. FleaTickRisk: a meteorological model developed to monitor and predict the activity and density of three tick species and the cat flea in Europe. Geospat Health. 2009;4(1):97–113.
 37.
Ogden NH, Milka R, Caminade C, Gachon P. Recent and projected future climatic suitability of North America for the Asian tiger mosquito Aedes albopictus. Parasit Vectors. 2014;7:532.
Acknowledgements
The IDEXX and Antech corporations are thanked for their data contributions. The Companion Animal Parasite Council financially supported this work. Two referees are thanked for comments that improved this paper.
Funding
YL and RBL were partially supported by Grant DMS 1407480 from the National Science Foundation. CM was partially supported by Grant R01 AI121351 from the National Institutes of Health.
Availability of data and material
The data used in this study are available from www.capcvet.org.
Authors’ contributions
All authors made substantial contributions to the manuscript. YL, CSM, and RBL performed the statistical analysis and contributed to the writeup. DDB, SN, and MJY wrote parts of the manuscript. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable (retrospective study using data that are publically available).
Author information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.
About this article
Received
Accepted
Published
DOI
Keywords
 Autoregression
 CAR Model
 Headbanging
 Heartworm
 Kriging
 Prevalence
 Spatiotemporal correlation
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate. Please note that comments may be removed without notice if they are flagged by another user or do not comply with our community guidelines.