Investigating the spatial risk distribution of West Nile virus disease in birds and humans in southern Ontario from 2002 to 2005

Background The West Nile virus (WNv) became a veterinary public health concern in southern Ontario in 2001 and has continued to threaten public health. Wild bird mortality has been shown to be an indicator for tracking the geographic distribution of the WNv. The purpose of this study was to investigate the latent risk distribution of WNv disease among dead birds and humans in southern Ontario and to compare the spatial risk patterns for the period 2002–2005. The relationship between the mortality fraction in birds and incidence rate in humans was also investigated. Methods Choropleth maps were created to investigate the spatial variation in bird and human WNv risk for the public health units of southern Ontario. The data were smoothed by empirical Bayesian estimation before being mapped. Isopleth risk maps for both the bird and human data were created to identify high risk areas and to investigate the potential relationship between the WNv mortality fraction in birds and incidence rates in humans. This was carried out by the geostatistical prediction method of kriging. A Poisson regression analysis was used to model regional human WNv case counts as a function of the spatial coordinates in the east and north direction and the regional bird mortality fractions. The presence of disease clustering and the location of disease clusters were investigated by the spatial scan test. Results The isopleth risk maps exhibited high risk areas that were relatively constant from year to year. There was an overlap in the bird and human high risk areas, which occurred in the central-west and south-west areas of southern Ontario. The annual WNv cause-specific mortality fractions in birds for 2002 to 2005 were 31.9, 22.0, 19.2 and 25.2 positive birds per 100 birds tested, respectively. The annual human WNv incidence rates for 2002 to 2005 were 2.21, 0.76, 0.13 and 2.10 human cases per 100,000 population, respectively. The relative risk of human WNv disease was 0.72 times lower for a public health unit that was 100 km north of another public health unit. The relative risk of human WNv disease increased by the factor 1.44 with every 10 positive birds per 100 tested. The scan statistic detected disease cluster in the bird and human data. The human clusters were not significant, when the analysis was conditioned on the bird data. Conclusion The study indicates a significant relationship between the spatial pattern of WNv risk in humans and birds.


West Nile Virus
West Nile virus (WNv) was first isolated and identified in 1937 from the blood of a resident of the West Nile district of Uganda [1,2]. Subsequently WNv caused outbreaks of human cases in Egypt, Israel, South Africa and in some parts of Europe and Asia [3]. WNv became a veterinary public health concern in North America in August of 1999, signaled by an outbreak in New York City. There are a number of theories on how the virus was able to survive and be transmitted during the spring of 2000. One is that infected mosquitoes from the 1999 New York outbreak were able to survive by hibernating through the winter in underground sewers, abandoned buildings, and bunkers [4]. Another implicates chronically infected migratory birds that may have reintroduced the virus after returning from the south the following spring [5].
The virus appeared in Canada in 2001, where the first mosquito and bird cases were recognized in Ontario in August [6]. The first human WNv cases in Canada occurred in Ontario and Quebec, in 2002 [6,7].
WNv is an emerging pathogen in Canada causing disease in animals and humans. As a Flavivirus of the Japanese encephalitis virus serogroup [8] WNv is maintained in an enzootic cycle involving viremic birds and ornithophilic mosquitoes, particularly Culex species. As a spill over effect WNv may be transmitted to humans or other deadend hosts, if the mosquitoes change their host preference [9] or if a bridge vector is involved [10].
Birds are the most important reservoir host and are able to amplify the disease since they develop high-level viremia (increased quantity of virus that replicates and circulates within the blood of the host) and remain infectious for several days [4,5]. In North America, the virus has been found in more than 150 bird species. Of these, corvids are among the most susceptible to infection and comprise an auspicious component of the mortality [11,12].
Wild bird mortality has been investigated for tracking the geographic distribution of WN virus in North America [7,12,13]. During the outbreak in New York City in 1999, a large die-off of birds, especially corvids, was associated with the outbreak in humans, both spatially and temporally [14]. Reports of the extensive die-off of birds preceded the epidemic in humans for the majority of regions [15,16]. For example, Marfin et al. [15] found that in 2000 in Northeastern United States, all 21 infected individuals had an illness onset date that was at least 15 days after the date that WNv infected birds were first collected in the individual's county of residence. From this, it has been suggested that dead birds can provide an early warning system to help predict areas of high human risk [15].
Dead bird surveillance programs may allow prevention and control methods to be intensified, before an outbreak of human cases occurs. Dead bird surveillance data are commonly used in assessing WNv risk; however different modeling approaches have been explored.
A study by Johnson et al. [17] quantified the association between clusters of dead crow sightings and onset of human WNv case in New York State. The risk in humans was positively associated with living in towns in proximity to dead crow clusters.
Theophilides et al. [18] developed the dynamic continuous-area space-time system (DYCAST) to identify and monitor high risk areas for WNv infection. The Knox test was used to assess the significance of space-time interaction in dead bird reports as an indicator of an intense WNv amplification cycle. It successfully identified areas of high risk for human WNv infection in areas where five of seven human cases resided, at least 13 days prior to the onset of illness, in New York City in 2001.
A study by Eidson et al. [19] evaluated the usefulness of dead bird surveillance in New York City in 2000 for detecting the geographic spread of the WNv and for providing an early warning system for humans. This study found that a steep increase in the number of dead crow sightings predated the onset date for the first human case and the increase in WNv positive birds by several weeks.
Watson et al. [20] assessed the spatial relationship between the locations of dead crow sightings reported early in the transmission season and the residences of WNv infected individuals in Chicago in 2002. Smoothed dead crow density values generated using kernel estimation were reclassified into high and low crow mortality areas. This study identified a spatial association between early season crow deaths and WNv infected residences of Chicago. Among humans the crude rate for WNv infection was 10.8 cases per 100,000 inside the high crow mortality areas compared to 3.2 cases per 100,000 outside.
The aim of this study was to investigate the latent risk distribution of WNv disease among birds and humans in southern Ontario and to compare the spatial risk patterns for the period 2002-2005, with a relatively simple but robust analysis. The objectives were to (1) describe the spatial variation of crude rates and smoothed risk for WNv among the public health units (PHUs) of southern Ontario from 2002 to 2005, (2) describe the geographical risk distribution and variation of WNv disease in humans and tested dead birds in southern Ontario from 2002 to 2005 in the form of risk maps, (3) investigate the potential for the tested dead bird data to be used as an indicator of human WNv risk, (4) explore the bird and human data for disease clusters.

Disease mapping
Spatial epidemiology is the description and analysis of geographic variations in disease with respect to demographic, environmental, behavioral and infectious risk factors [21]. In public health, identification and quantification of patterns in disease occurrence provide the first steps toward increased understanding and possibly, control of that particular disease [22]. In order to better understand the spatial epidemiology of the WNv, including the trends and clusters within the risk distribution, disease maps are essential.
Disease maps provide a rapid visual summary of geographic information and may identify subtle patterns in the data that are missed in tabular presentations. They are used variously for descriptive purposes, for surveillance to highlight areas at apparent high risk, to aid resource allocation [21], to identify possible disease clusters and to show changes in disease patterns over time [23].
There are various mapping techniques that can be used depending on the type of spatial data being analyzed. Dot or spot maps are used to visualize spatial point data. These maps indicate the location of case or event data on a geographical map, such as the location of WNv cases. Choropleth maps are used to explain spatial variation in regional count data, such as the number of WNv cases per PHU. These maps symbolize regional statistical data within the boundaries of geographic regions grouped into classes. Each class stands for a range of values and is represented by a logical sequence of gray (or color) tones.
When regional health information, such as the incidence rate of WNv at the public health unit level is considered as point measurement data at the regional centre (i.e. geostatistical data), then isopleth maps can visualize the spatial distribution of the latent risk. Isopleth maps show the distribution of spatially continuous phenomena by a logical sequence of gray (or color) tones that symbolize equal values. Isolines are often overlaid on top of an isopleth map to indicate threshold values.

Data sources
In February 2000, the Public Health Agency of Canada (PHAC) organized a National Steering Committee to develop a coordinated approach to respond to WNv [24]. As a result, a surveillance program was established to monitor WNv in humans, mosquitoes, birds and horses. The human and bird data utilized in this paper were obtained from this surveillance program for the years 2002 to 2005. The bird surveillance data including the number of birds tested and the numbers of WNv positive birds for each PHU were obtained from the Canadian Cooperative Wildlife Health Centre (CCWHC) WNv database. The human surveillance data, including the number of human WNv cases for each PHU were obtained from the Ontario Ministry of Health and Long Term Care (OMHLTC).
Estimates of the total population for each public health unit were calculated in 2001 for the Census of Canada and were obtained from Statistics Canada [25].

Data collection
Dead birds were submitted for WNv diagnosis by both the general public and public health personnel as part of the National WNv surveillance program. The dead birds were tested for WNv infection by the CCWHC, provincial laboratories and the PHACs National Microbiology Laboratory in Winnipeg. In 2002 submitted dead birds were tested by the reverse transcriptase polymerase chain reaction (RT-PCR) test, which detected WNv ribonucleic acids (RNA) [26]. From 2003 to 2005 testing of submitted birds was carried out as a modified two-stage test. The first test was the oral VecTest™ (Medical Analysis Systems, Camarillo, CA), which was performed on oropharyngeal swabs from birds and detects WNv antigen. The sensitivity and specificity of the oral VecTest™ for crows collected in 2002 were 83.3% and 95.8%, respectively [27]. The second test, RT-PCR was used on only the positive birds that were the first bird in the PHU or major municipality of a PHU to be found positive. This sequential testing was performed in order to minimize false test positive results that would initiate unnecessary public health interventions. A suspected first positive dead bird in a PHU was defined as a positive WNv case if it tested positive using the oral VecTest™ as well as the RT-PCT test. After confirmation of WNv activity in a particular public health unit by RT-PCR, all positive oral VecTest™ results were assumed to be positive WNv cases.
Probable or confirmed human WNv cases were reported to local and provincial health authorities by health care providers, since the WNv infection is a reportable disease of humans in Canada [28]. Blood samples of individuals suspected to have symptoms of WNv infection were sent to the OMHLTC Central Public Health Laboratory (CPHL). Testing of the blood samples involved three tests in series. The first test was the IgM enzyme-linked immunosorbent assay (ELISA), which if positive was run again to rule out false positive results. These two tests could be followed by the plaque reduction neutralization test (PRNT) to confirm the diagnosis [29]. A human blood sample was defined as a positive WNv case if the sample tested positive on both the IgM ELISA tests (probable case) or if the sample tested positive on both the IgM ELISA tests and the PRNT (confirmed case).

Study area
The occurrence of WNv disease was investigated in the 30 PHUs of southern Ontario. PHUs are official health agencies that are responsible for administering health promotion and disease prevention programs [30]. PHUs were used as components of the study area in order for prevention and control methods to be implemented by epidemiologists or health care authorities of each unit. The PHU boundaries and the number of PHUs located in southern Ontario were consistent from 2002 to 2004. The names and distribution of the PHUs of southern Ontario are shown in Figure 1. In 2005 the Muskoka-Parry Sound PHU was dissolved. However for consistency of the results, the original 30 PHUs were used for the analyses of all four years.

Statistical analyses
Crude PHU specific human incidence rates of WNv cases for the years 2002 to 2005 as well as the average annual incidence rates were calculated. The four annual incidence rates were calculated as the number of human WNv cases divided by the population totals for each PHU. The average incidence rates over the four years for each PHU were calculated by the sum of all four annual incidence rates divided by four. The four annual incidence rates and the average of the four years for each PHU were expressed as choropleth maps using the same incidence scale. This enables the reader to visualize the changes in crude incidence rates by region from year to year. The four annual mortality fractions and the average mortality fraction over the four years for birds were calculated in the same way, except, the bird WNv cases were divided by the total number of tested dead birds for each PHU. These were also expressed as choropleth maps.
There are often problems associated with mapping crude rates and fractions. For example, crude estimates of disease occurrence can be unreliable and highly variable as they are based on sample sizes or populations at risk, which can vary greatly from region to region [31]. In order to counter these problems, empirical Bayesian smoothing, also known as shrinkage estimation, was used to smooth the estimates.
Empirical Bayesian smoothing was performed to reduce the variance heterogeneity across the regional estimates [32]. The empirical Bayesian estimate is a weighted mean of the crude regional and global estimates. The weighting is determined by the variability of the estimates [33]. When the regional population is relatively large, the estimate is shrunk towards the global mean to a lesser extent and more weight is given to the regional estimate as there is more confidence in the precision of this estimate [31,32]. When the regional population is small, the shrinkage effect towards the global mean is stronger, as more weight is given to the global estimate, since there is less confidence in the precision of the regional estimate [31,32]. For empirical Bayesian smoothing, it was assumed that the bird data follow a binomial distribution and the human data a Poisson distribution.
The human and bird annual smoothed estimates, calculated by applying the empirical Bayesian smoothing method, are smoothed estimates of incidence rates and mortality fractions respectively, but for ease of reference, we refer to them collectively as smoothed risks. The average smoothed risks over the four years were calculated by the sum of all four smoothed annual risks divided by four. Choropleth maps as described above were then created with the empirically Bayesian smoothed risks for both the human and bird data. Again these choropleth maps were based on the same colour scale. Parallel boxplots of the crude rates and mortality fractions and the empirically Bayesian smoothed risks were created to show the shrinkage effect.
Choropleth maps are effective in showing the variation in regional data but are also known to have problems associated with them as outlined in the discussion. The geostatistical prediction method of kriging was used for isopleth mapping to overcome these problems. Kriging is an approach to interpolate or spatially predict regional data onto a continuous surface [34]. The kriging predictor is a weighted average that is calculated from the entire sample with weights depending on the semi-variogram [35]. The weights are constructed to give regional risk estimates more influence on the prediction the closer they are to the prediction sites and to downplay a cluster of points that contains largely redundant information [32]. The semivariogram is a graphical representation of the variation between sampling points separated by a given distance and direction [31]. The semi-variogram was estimated by the weighted least squares estimation method [35]. Due to a limited number of PHUs, the maximum likelihood estimation method (MLE) was used for any years of data that could not be estimated by the weighted least squared estimation (WLSE) method. Risk maps for both the bird and human data were created to identify high risk areas and to investigate the potential relationship between the bird mortality fractions and the human rates. This was first done by visually checking for an overlap in the high risk areas on bird and human risk maps. The isopleth maps for all years investigated were created to have included the range of risk values for each particular year so that the bird high risk areas could be compared to the human high risk areas.
Furthermore, a Poisson regression analysis with an overdispersion parameter to control for clustering (i.e. spatial dependence) was used to model regional human WNv case counts as a function of the smoothed regional bird data and the spatial coordinates in the east and north direction [36]. A normal quantile-quantile (QQ) plot of scaled deviance residuals was plotted to evaluate the fit of the model and to identify public health units as potential outliers, i.e. disease clusters. The model is formulized as follows: The presence of disease clustering and the location of disease clusters were investigated by the spatial scan test [37]. The spatial scan test is a likelihood ratio test, which uses circular scanning windows of various sizes and positions. By continuously changing the circle center and radius, the window scans the geographic area for potential localized clusters without incorporating prior assumptions about their size and location [38]. Circular search windows began with individual PHUs and expanded to include neighboring PHUs until a maximum of 50% of all dead birds investigated was reached, or -in the case of the human data -50% of the total population at risk. A Bernoulli probability model was used to calculate the likelihood for the bird data, while the Poisson probability model was used with the human data. The spatial scan statistic tested the null hypothesis that the risk of the WNv disease within the window was equal to the risk outside the window, while the alternative hypothesis stated that there was an elevated risk for WNv infection within the windows as compared to outside the window. The P-value was obtained by Monte Carlo hypothesis testing, by comparing the rank of the maximum likelihood for the observed dataset to the maximum likelihoods of 999 simulated datasets [39]. Significant disease clusters (α = 0.05) were indicated on the choropleth maps of empirical Bayesian smoothed risks.
In order to investigate if the smoothed bird data helped to explain the human disease clusters the spatial scan test was repeated again for the human data. This time the smoothed bird data were added as a covariate. If the human clusters that were significant in the original spatial scan test (no covariate added) were no longer significant, then the covariate explained the human cluster or spatial distribution of human WNv disease.
All statistical analyses were carried out within R [40], except for the disease cluster analysis with the scan statistic, where SaTScan [39] was applied.   respectively. The annual maps show no clear spatial pattern of trend or clustering, except for 2002, when a potential cluster with a spike at the Halton PHU in the centre of the study area is indicated. This clustering seems to be consistent for all four years, since it is also indicated on the average map (Figure 5e). Perhaps most interesting is the smoothed risk map for 2004 (Figure 5c), as there were almost no human WNv cases reported and the map is constant over the entire study area at the lowest risk level.

Tables
In order to generate isopleth risk maps from smoothed data via kriging, spatial dependence was modeled by exponential semi-variograms without nugget effect. The models were fitted by weighted least squares to robust empirical semi-variograms for all of the bird data and for the first two years of human data. Semi-variogram models for the 2004, 2005 and the average human data were fitted by maximum likelihood estimation. The variogram clouds identified Halton PHU as an outlier within the 2002 and 2003 human data. Halton PHU was therefore excluded from the model fitting process. There were no outliers within the bird data. Table 3 shows the results of the range and sill values of the empirical semi-variogram for both the bird and human data for all years investigated. All empirical semi-variograms for the bird and human data over the years 2002 to 2005 leveled out and reached a sill. The semi-variograms based on the bird data indicated that the range of the semi-variogram increased from year to year. This showed that as the disease spread among birds in southern Ontario, the data were correlated over longer distances from year to year. For the human data, the semi-variogram parameters varied over the four years without any tendency.  (Figure 6e) it is seen that there are two major high risk areas in the centre and the south of the study area, i.e. around the Toronto-Peel-Halton and Windsor PHUs. Furthermore the spatial downwards trend in bird risk from south to north is clearly visible from this isopleth map. Similar results can be seen on the isopleth risk maps for humans. Although these maps do not exhibit a spatial trend, again the Toronto-Peel-Halton or Windsor PHUs stand out from the maps as high risk areas for human WNv illness (Figure 7).
The annual risks for southern Ontario were estimated by the intercept parameter of the spatial regression models used for kriging. The overall bird and human smoothed risks follow the same temporal pattern: steady decrease  Figure  5e). In addition, PHUs that were found to be outliers in the Poisson regression analysis were also identified as part of the disease cluster by the scan test. The spatial scan test applied to the human data with the smoothed bird risk estimates as a covariate, resulted in no significant clusters for any of the four years or the average of the four years.

Discussion
Using isopleth maps to display regional data has many advantages. For example isopleth mapping techniques overcome the problems or disadvantages of choropleth mapping. Choropleth maps can be misleading as the uneven shape and size of the different regions produce a visual bias [32]. For example, physically large areas that may be sparsely populated tend to dominate the perception of the map and may detract from smaller, sometimes more important regions, depending on the study [31]. Regional boundaries between census tracts also have problems associated with them. The risk of disease occurrence jumps artificially at boundaries from census tract to census tract and the risk is incorrectly assumed to be constant throughout each region or administrative area [34]. Isopleth mapping eliminates these problems since the risk distribution is mapped as a spatially continuous phenomenon that does not follow administrative boundaries.
The isopleth risk maps indicated that there were certain areas in southern Ontario where populations were at higher risk for acquiring WNv disease. The bird isopleth risk maps exhibited high risk areas that were relatively constant from year to year as the majority of the same PHUs were included in the high risk areas each year. The southern region of southern Ontario consistently contained the highest risk areas for birds over the entire study period. The high risk areas for humans did not exhibit quite as strong a pattern as the bird data, but there were a few PHUs that were high risk areas for the majority of the years investigated. These included the City of Toronto PHU, Kent-Chatham PHU and the Windsor-Essex PHU. By visual inspection of the risk maps for birds and humans for each year, it was found that there was an overlap of high risk areas. The overlaps occurred in the centralwest and south-west areas of southern Ontario. Overlapping risk areas indicated that the high risk areas for birds can help to predict the occurrence of WNv disease in humans.
Human WNv disease clusters seemed to move in the south-west direction toward the most southern tip of southern Ontario. These clusters were still located within the boundaries of the bird WNv disease cluster. This movement of the human cluster may be an effect of the limited case numbers, which are furthermore not discriminated for the various types of WNv related illness.
WNv disease clusters as highlighted on the isopleth risk maps as well as identified by the spatial scan test, were meaningful in southern Ontario. These clusters cannot be regarded as chance clusters since they exist within both birds and humans for each year investigated. Neither can they be explained by perception bias (e.g. increased public awareness of WNv in urban areas compared to rural areas) because other urbanized centers do not show increased risk.
The analysis of disease clusters was an important part of this study. Kuldorff's spatial scan test used circular scanning windows to detect the potential cluster areas. A scan test that detects potential clusters based on non-circular windows may be more appropriate. For example, the flexible spatial scan statistic proposed by Tango and Takahashi [41], has the ability to detect noncircular clusters more accurately than circular clusters. Irregular shaped cluster should be identified in future research, but may require a sample size larger than 30 PHUs.
Both the scan statistic and the Poisson regression model indicated a relationship between the WNv cause-specific Parallel box plots for the raw annual human incidence rates of WNv disease per 100,000 population (raw) and the corre-   [42]. Rainfall is another variable to consider. A heavy rainfall will decrease the incidence of WNv as it flushes out sewers and drains that contained stagnant water, thus eliminating mosquito breeding sites. This study is analyzing annual data, which makes it difficult to incorporate weather events in a statistical model, e.g. heavy rainfall can prevent WNv disease for a few days but may further it when followed by a heat wave. Weather variables available at the PHU level should be included in a study that uses weekly data.
Although the tested dead bird data were useful for investigating the relationship between bird and human WNv occurrence, there were limitations or biases associated with the data. Passive surveillance is voluntary and is thus dependent on the public for finding and reporting dead birds. Surveillance can be affected by human related variables such as public awareness, public interest, media coverage and human density [28,43]. The PHUs of southern Ontario submitted birds at their discretion, therefore the analysis was highly dependent on non-homogenous sampling techniques. This may help to explain why the risk maps varied slightly from year to year. A more concise and standard sampling technique for all PHUs would be beneficial in the future. However, the surveillance data were adequate for the surveillance program's primary motivation of finding out when the WNv was active in each PHU.
The PHUs testing period (i.e. when the PHUs started and stopped testing dead birds) also manipulated the data. PHUs that did a lot of testing earlier in the season were not going to see as many cases as the PHUs who did the majority of their testing later in the year. This problem could be minimized by using tested dead bird data for the peak outbreak time as opposed to the whole year, thus increasing the stability of the risk pattern over time. The amount of testing also varied from year to year. For example, in 2005 the majority of PHUs carried testing further into the year as there were fewer cases in the beginning of the year.
In contrast, tracking dead bird sightings would have avoided delays associated with specimen collection and  testing and would have allowed rapid recognition of trends in viral activity and the potential for occasional human cases or an outbreak [44]. Dead bird sighting data for southern Ontario were not used in the present analysis as not all public health units had sighting reporting systems in place. Furthermore, the available data were considered unreliable at the public health unit level because of sporadic and inconsistent sampling from year to year. In addition, the bird sighting data shared the human related issues mentioned above.
This study found a spatial relationship between the WNv cause-specific mortality fraction in birds and human WNv incidence rates, while other studies found a temporal relationship [20,38,45]. These studies have indicated that human WNv cases (onset of illness) will usually occur 2 to 5 weeks after bird WNv cases. By combining both the spatial and temporal relationship, prevention and control strategies can be implemented in certain high risk locations for specific high risk time periods. If public health authorities and health care authorities know where the high risk areas are for birds, then control measures such as application of adulticides or larvacides can be used to prevent human cases by reducing the vector population. Another important prevention and control method is public education. This initiative provides knowledge of the disease itself, how it is transmitted and how to prevent or reduce the risk of exposure. For example, a study by Loeb et al. [45] in the Halton region discovered that individuals who practiced at least two personal protective behavior traits such as wearing long sleeves, pants, insect repellent or eliminating standing water around one's home had a 50% reduction in the risk of infection.

Conclusion
The risk pattern and locations of disease clusters among southern Ontario dead bird data were stable over time.
The human WNv risk pattern and disease clusters were slightly more variable over time. The isopleth risk maps indicated that the WNv risk patterns for birds and humans overlapped. The spatial disease cluster analysis indicated that all human clusters fell within regions contained in the bird cluster. This indicates that WNv is established in southern Ontario.
There is a relationship between observed bird and human WNv cases as shown by the scan statistic and the Poisson regression model. Smoothed bird risk estimates can be used as indicators of high risk areas for humans. By knowing where high risk areas for birds and humans are, prevention and control methods can be intensified, thus mitigating human morbidity and mortality.