Human case data were taken from the weekly U.S. Geological Survey West Nile maps for the 2003 transmission season based on county-level data provided by ArboNet through voluntary reporting by state and local health officials to CDC (4). The case numbers comprise reports of mild West Nile fever as well as the more severe West Nile meningitis or encephalitis. Crude county-specific incidence rates were calculated by using the Census 2000 county population totals.

We created a human risk map for WNV based on the crude human incidence early in the transmission season, on August 13, 2003. A disease map that displays observed human incidence will show not only spatial variation in risk but also random variation resulting from low case numbers relative to the base populations. Removing random noise permits improved estimates of disease risk (2). We have approached this procedure by finding the estimates of expected incidence from a conditional autoregressive model (5,6). The model helps remove random variation based on the premise that contiguous regions tend to have similar disease risks, when compared to regions that are far apart. We applied the conditional autoregressive model to calculate expected WNV incidence rates (Appendix).

The first step was to identify the adjacent neighbors of each county by using a geographic information system (GIS, ArcView 3.2, ESRI, Redlands, CA). A data file that included the number of cases, total population, and number and names of neighboring counties for each county was then generated in SAS (SAS Institute Inc., Cary, NC). The file was imported into WinBUGS v1.4 (Imperial College, St. Mary's, UK; and Medical Research Council, Cambridge, UK). This software implements a simulation process to estimate model parameters, including improved estimates of WNV incidence rates. These estimates were then brought back into GIS to display the human WNV risk map.

To verify the method's potential use as an early warning system for human risk, we calculated the validity of the model-estimated risk map versus the raw incidence map from August 13 for predicting the case distribution for October 1, 2003. The two time points represent an ≈14-fold increase in total cases, from 399 to 5,685. For each of the three disease maps, counties were grouped into high- and low-risk classes on the basis of WNV incidence. High risk was defined as human incidence >1 case per 1 million population for the August 13 maps and 1 case per 100,000 for the October 1 map, findings that reflect the change in risk over time. The sensitivity of the method for predicting risk was calculated as the proportion of high-risk counties on October 1 that was correctly identified as such by the model-estimated August 13 risk map. Similarly, specificity was defined as the proportion of low-risk counties on October 1 that was correctly identified as such by the modeled risk map. The sensitivity and specificity values were compared to those obtained when the raw August 13 incidence map was used to predict risk on October 1. Measure of agreement between risk classes of the August 13 map and the October 1 map was assessed by the k statistic, which accounts for the degree of overlap expected by chance alone; k has a range of 0 to 1; values of <0.4 represent poor agreement (7).

We assessed the quantitative predictive ability of the nonhuman surveillance systems by fitting a regression model to the rate of WNV human cases for counties with the final USGS maps for the 2002 season (4). The model includes covariates for the number of virus-positive tissue samples from dead and diseased wild birds and virus-positive mosquito pools, both provided by state health officials at the county level (Appendix). Each covariate was considered together and separately to determine its contribution for predicting WNV incidence. The model was fitted by using GENMOD in SAS (SAS Institute). The contribution of nonhuman surveillance systems to variability in human risk was determined by calculating the proportion of the deviance explained (R²).

Because the human West Nile virus (WNV) case number is low relative to the base population, it was assumed to have a Poisson distribution. Under general conditions, the Poisson provides a good description for the distribution of the numerator for an incidence rate (15). However, our model also allowed for the estimation of "extra-Poisson" variation in case it is also needed to provide an accurate description of these data.

The log linear model used for spatial smoothing assumed that the number of disease cases in the i-th county, n_i, has a mean, P_i2_i, where P_i, is the denominator for the rates and 2_i = exp [a₀ + b_i + h_i} where a₀ is the intercept, b_i is the spatially correlated random variation with mean 0 and variance , and h_i the unstructured extra-Poisson variability with mean 0 and variance . In addition, we assume that both the spatial and the unstructured variability have Gaussian distributions, which are independent in the latter case. On the other hand, the mean for the spatial component, conditional on the means for the contiguous neighbors, is

where r_i is the number of neighbors for region i. The adjacent neighbors for each county were determined by using a geographic information system (GIS, ArcView 3.2, ESRI, Redlands, CA). Thus, the overall log linear model for the number of cases in the i county that incorporates both spatial correlation and unstructured variability is log n_{i = log} P_{i +} b_{i +} h_{i +} a_0.

The population size for county i (P_i) was determined from the Census 2000 data.

Markov Chain Monte Carlo (MCMC) simulation methods were used to find Bayesian estimates of the model parameters as implemented in WinBUGS v1.4 (Imperial College and Medical Research Council) (16,17). Gamma prior distribution parameters were assumed for the variances of the Gaussian distributions, and a plot of the history of the simulation was used to determine the number of iterations required for the process to equilibrate. The approach provides improved estimates of county-specific rates that have been spatially smoothed.

In the MCMC method, parameters estimated from each step are used in turn to determine values for the next step; therefore, a good set of initial values is essential before gleaning the values that will be used in the estimation. To accomplish robust parameter estimates, an arbitrary set of values was chosen, and the number of successive steps taken to stabilize the simulations was noted, which is known as the burn-in. The burn-in period was determined through the use of two chains and the modified Gelman-Rubin convergence statistic. This statistic indicates the point at which the process stabilized by describing how well the chains overlap. Final estimates were obtained by using 1,000 iterations as the burn-in period, and the next 9,000 were used as the sample for deriving the Bayes estimates of the smoothed WNV incidence rates.

The quantitative predictive ability of the nonhuman surveillance systems was assessed by once again fitting a log-linear model to the rate of WNV human cases. For this analysis, we instead used a maximum likelihood approach, in which we assumed a Poisson distribution for the number of cases, allowing for extra-Poisson variation by estimating the scale factor. In this model, log n_i = log P_i + β_AM_i + β_MA_i + a_o where P_i is the population offset, A_i is avian mortality attributable to WNV, and M_i is the number of virus-positive mosquito pools. The model was implemented by using GENMOD in SAS (SAS Institute Inc., Cary, NC). Only counties that submitted both mosquito and bird samples were included in the analysis (N = 382) (Figure A1).