A comparison of least squares regression and geographically weighted regression modeling of West Nile virus risk based on environmental parameters

Introduction

West Nile Virus (WNV) is a vector-borne disease that was first detected in the United States in 1999 (Nash et al., 2001). Within a few years the virus had spread across the North American continent (Hayes et al., 2005). WNV has had important environmental and human impacts, including a decline in numerous bird species (CDC) and increased morbidity and mortality among humans. This has also resulted in increased economic burdens due to initial acute health care needs of infected individuals and subsequent long-term costs associates with infection, estimated at approximately $56 million per year between 1999 and 2012 (Barrett, 2014). Because that study indicated how difficult predicting and planning for WNV outbreaks was, we became interested in developing a spatially explicit model using environmental factors in an attempt to improve WNV risk predictions.

There are two important considerations that should typically be examined when developing spatially explicit environmental disease risk models (Miller, 2012). The first should be an examination of potential spatial autocorrelation (the degree to which a set of spatial features and their associated data values tend to be clustered together in space). This involves accounting for whether environmental factors and the corresponding disease rates in geographically proximate areas are more or less clustered together than they are in geographically distant areas. Second, data non-stationarity (changing means, variances and covariances in data across space) should be investigated and controlled when necessary (Fotheringham, 2009a; Miller, 2012). Geographically weighted regression (GWR) can be used for these two considerations and can often produce improved models that enable better spatial inference and prediction. Recent studies have applied GWR modeling to drug-resistant tuberculosis versus risk factors (Liu et al., 2011); environmental factors versus typhoid fever (Dewan et al., 2013); local climate and population distribution versus hand, foot, and mouth disease (Hu et al., 2012); and environmental factors and tick-borne disease (Atkinson et al., 2012; Atkinson et al., 2014; Wimberly, Baer & Yabsley, 2008; Wimberly et al., 2008), all showing that predictor variables varied spatially across large geographic regions, implying that the results for such studies may be improved using GWR.

The spatially explicit model that is discussed in this paper uses GWR to account for spatial heterogeneity for two reasons: (a) WNV disease risk observed across space may be related to similar environmental variables that increase vector habitat suitability and (b) environmental variables that influence WNV risk are not typically uniformly distributed across geographic space. Although many epidemiological models of WNV risk have been developed, it appears that there has been little research to explicitly examine techniques that account for spatial heterogeneity. Most models assume that the impact of various environmental factors are constant across the study region, which is unrealistic as larger areas display substantial variations in distribution of environmental, socio-economic, and demographic conditions (Goovaerts, 2008).

Due to the unavailability of reliable and complete data, developing models of WNV risk pose additional challenges. Human case data is lacking due to issues of under-reporting and limited surveillance. Our alternative strategy was to assess WNV infected dead bird counts as a surrogate measure of human risk because “infection rates” in dead birds can be more precise because of the genetic markers tested in dead birds may be more reliable than case data and/or surveillance data. Additionally, others have also used mosquito habitat suitability as a surrogate for estimating WNV risk for human infection (Cooke, Grala & Wallis, 2006). For our study, we followed a similar approach and used a model of mosquito habitat suitability condition as a predictor of the spatial distributions of infected birds, which in turn can be used to estimate WNV disease risk among human populations. Further, because the environmental variables considered in this study are known to vary across space, we account for spatial autocorrelation and non-stationarity using GWR following the approach of (DeGroote et al., 2008) in order to improve the predictability of a model.

Statistical Considerations

Miller (2012) suggests that a ‘global’ model is the one that assumes that the parameters (commonly mean and variance) of some process are constant across geographic space (commonly mean and variance), typically referred to as the spatial stationarity of a process. Miller suggests that in the case these parameters vary across geographic space (spatial heterogeneity), then such models may lead to inaccurate predictions and subsequent problems for decision-making. In an ecological context, spatial heterogeneity usually results from the interaction of various environmental processes that operate at different scales (Legendre, 1993). Fotheringham (2009b) used local statistics for linking the concepts of spatial autocorrelation and heterogeneity that are deemed important when developing spatial models. Local statistics disaggregate a global mean value into locally computed values for each spatial unit. It is based on a conceptualization of Tobler’s first law in Geography (Tobler, 1970) that specifies that “everything is related to everything else, but near things are more related than distant things.” Spatial autocorrelation is a commonly used measure of the degree of spatial heterogeneity.

GWR is a local regression method that can be used for diagnosing spatial heterogeneity between dependent and explanatory variables over space (Fotheringham, Brunsdon & Charlton, 2003). It is performed within local windows centered on the nodes of a regular grid. Each observation within the local window is weighted based on its proximity to the center of that window. This approach has several advantages: it avoids abrupt changes in the local statistics computed for adjacent windows, helps visualize spatial variability within the geographic entity, and allows analysis of regionally aggregated data (Goovaerts, 2008). A model’s predictive ability, particularly in ecological modeling, is influenced not only by the strength of relationships between the species and its environment, but whether the model recognizes if the relationships are operating at multiple spatial scales. GWR provides a framework for exploring scale-dependent effects. It tests the effect on a model’s predictive ability by systematically increasing the size of the local window (Miller, 2012).

GWR can be used for mapping the spatial distribution of a model’s coefficient values in order to identify potential missing variables or to suggest other underlying factors associated with the observed non-stationarity (Miller, 2012). GWR is also useful for exploratory data analysis and visualization; for example Kupfer & Farris (2007) used a ‘leave-one-out’ (jackknifing) methodology to compare residuals from GWR and ordinary least squares regression. They found that GWR often had more accurate predictions for sites that were difficult to predict (where both models had overall higher residuals), which is why we used a GWR framework for explicitly modeling the spatial relationships between WNV and its environmental risk factors.

Materials & Methods

Results

LSR modeling identified land surface temperature (VIF = 1.046), stream density (VIF = 1.177), and road density (VIF = 1.143) as statistically significant (p < 0.05) variables related to WNV risk: (1)${\displaystyle \mathrm{WNV\; risk}=-75.87+595.60\left(RD\right)+1.89\left(LST\right)-146.89\left(SD\right)}$ Where:

WNV risk = average infected dead bird count

RD = road density

LST = land surface temperature

SD = stream density.

The histogram of the LSR model’s residuals approximates that of a normal curve, with a non-significant (0.134, p > 0.05) Jarque–Bera statistic (Jarque & Bera, 1980), and the Moran’s I Index Z-score (1.23) all imply that the model is unbiased and significantly different than random.

However, the Koenker statistic (0.000007*, p < 0.05) confirmed non-stationarity in the LSR model indicating that there is not a consistent relationship between the explanatory variables and WNV risk across the study area. Further, the presence of mild heteroskedasticity was noted in the LSR model. We conclude that the LSR model is stable but non-stationary, suggesting that proceeding with GWR model was warranted.

The GWR model in this study was implemented using the following algorithm: (2)${\displaystyle \mathrm{WNV}{\mathrm{risk}}_{\mathrm{(i)}}={\mathrm{\beta }}_{\mathrm{i0}}+{\mathrm{\beta }}_{\mathrm{(i1)}}{RD}_{\mathrm{(i)}}+{\mathrm{\beta }}_{\mathrm{(i2)}}{LST}_{\mathrm{(i)}}-{\mathrm{\beta }}_{\mathrm{(i3)}}{SD}_{\mathrm{(i)}}+{\epsilon }_{\mathrm{(i)}}}$ where β coefficients are county (i) specific, and RD is road density, LST is land surface temperature, and SD is stream density.

Comparing the fit of the global LSR model (assumes homogeneity of variables across space) and local GWR model (makes no assumption of homogeneity), we found that the global LSR adjusted R² was 0.61 (R² was 0.66, P < 0.05, Fig. 1) with analysis run on all 58 counties. The local GWR adjusted R² was 0.71 (R² is 0.75, p < 0.05, Fig. 2) with a bandwidth of 54, which suggests that there has been some improvement by using a local modeling approach. Our preferred measure of model fit, AIC, gave a value of 567.7 for the global model and 551.4 for the local model. The difference of 16.3 is relatively strong evidence of an improvement in the model fit to the data. Further, the problem of hetroskedasticity that was noted in the OLS model was not observed in the GWR model.

We also tested the results using different bandwidth parameter. Several iterations were run but it was observed that although a smaller band-width criterion gave an improved combination of AICc and adjusted R² values, it also compromised the model diagnostics by introducing local collinearity and thus instability in the model. Addressing local collinearity by removing the Counties having condition number greater than 30 affected the model’s overall results. Thus, it is better to have a larger band-width rather than violating model assumptions and to avoid the unstable prediction (Charlton & Fotheringham, no date; Nakaya, 2014).

Mapping the values of the standardized residual across California (Fig. 3A) provides a representation of: (a) areas with unusually high or low residuals and (b) whether the residuals were spatially autocorrelated. Counties with excessively large positive residuals would under-predict WNV risk, and counties with excessively large negative residuals would over-predict WNV risk. The spatial autocorrelation of GWR residuals for our model resulted in a Moran’s I value of −0.11 (p = 0.18), implying little evidence of any autocorrelation in them.

Local coefficient estimates for significant factors were mapped using quantile classification method. Figure 3B shows the variation in the model’s coefficient estimates for the land surface temperature (LST) variable. The map for the local coefficients reveals that the influence of this variable in the model varies considerably over California, with a strong north-south direction. The range of the local coefficient is from 1.26 for the northernmost counties to 3.06 for the southernmost counties—evidence that points to heterogeneity in the model structure within California. The global coefficient and all the local coefficients for this variable are positive—there is agreement between the two models on the direction of the influence of this variable. Figure 3C shows a similar distribution in north-south direction of positive road density coefficient. Figure 3D reveals the opposite for stream density coefficients, with larger values in the north and smaller values in south. Contrary to our initial thoughts, stream density demonstrated a negative relation to disease risk. This may reflect that flowing water is normally not suitable for larval development of the various species of mosquitos that commonly transmit WNV or that rasterizing the stream database into stream density introduces a component that is not yet fully understood.

Our best ordinary least squares model, the global LSR model (Eq. 1) produced an adjusted R² of 0.61 (p < 0.05) with a corresponding corrected AIC of 567.70. Utilizing the same environmental variables, our best local GWR model (Eq. 2) produced an adjusted R² of 0.71 (p < 0.05) with a corresponding corrected AIC of 551.4. A 16 point decrease in the AIC and approximately 16% improvement in the model performance suggest that incorporating spatial data improves the predictive ability of WNV risk.

Conclusions

One of the frequent technical issues in modeling disease risk is to incorporate local rather than global associations in these models (Foley, Charlton & Fotheringham, 2009). In spatial regression models, a global model can be used to examine the relationship between disease risk and potential explanatory factors which are based on the assumption that the relationship is a stationary spatial process (Miller, 2012). For a small and homogenous region of interest, it is reasonable to assume that the explanatory factors would not change significantly across the region, and the relationship between WNV risk and the potential factors would also be unchanged. However, important variables such as topography, climate, and population distribution change greatly when it comes to a large region like California with an area of over 163,000 square miles. California is geographically diverse and is equally varied in its range of climates with several climatic sub-regions recognized. It would be unexpected to find that the spatial stationarity assumption holds in such a large area.

The distinct north–south pattern revealed in our study could be attributable to typical latitudinal expressions of temperature and precipitation, especially since California has a north-south length of 1,350 km. This environmental pattern is also a likely contributor to the distribution of different mosquito species in the United Stated, especially notable in its manifestation in California. A recent report (CDC, 2013) shows that while Cx. Tarsalis is distributed throughout California, Cx. pipiens is a more important mosquito vector in northern California, while Cx. quinquefasciatus is more important in southern California. While WNV can be found in a wide variety of ecosystems, the north–south pattern of infected birds detected in this study may be expressed more noticeably in California due to the north–south differences in mosquito species distributions as observed in the Centers for Disease Control report.

Our results concur that understanding WNV risks is improved when considering spatial heterogeneity of the variables that affect the risk (Beck et al., 1994). Besides improving prediction accuracy, spatial heterogeneity can also provide insights into the underlying ecological processes controlling the distributions of vector populations and zoonotic pathogens (Wimberly, Baer & Yabsley, 2008) because GWR models consider spatial heterogeneity by separating the large heterogeneous region into smaller, more homogeneous local regions. Fotheringham (Fotheringham, 2009a), stated that an advantage of using GWR is that it accounts for much of the spatial autocorrelation in the residuals that is usually found in global modeling. Further, it is possible that a variable that is insignificant at the global level might be important locally.

There are several limitations of this study. First, it is assumed that factors suitable for mosquito habitat increase the likelihood of WNV spread in human populations. On the surface this seems to be reasonably apparent; however, we do not have specific evidence that this is true. Second, it is also assumed that the probability of human infection is higher in counties with multiple confirmed WNV bird cases, another reasonable conjecture with several references in the literature, but without direct confirmation. Potential problems with this assumption include varying human population density (e.g., two areas with the same number of infected dead birds reported but one area’s human population density is substantially different than the other), variations in level of public concern (as reports of infected dead birds increase, more people begin looking for dead birds), and resource availability might bias the reporting of dead birds (wealthy areas devote disproportionate resources to the issue). Thus, proper surveillance methods that take into consideration these limitations while collecting infected dead bird data will contribute to more meaningful results. Third, our approach assumes that people are infected within the county of their residence, ignoring the possibility of contracting an infection while traveling outside the county limits. Lastly, road density could also be correlated with dead bird surveillance effort and might be a potential bias for reporting dead birds. We recognize that if these assumptions do not hold, modeling WNV risks based on infected dead birds may yield biased results. However, if the assumptions do hold, the local modeling approaches should improve predictions of WNV risks.

The research described in this paper suggests that a spatially explicit local model using GWR approaches to adjust for spatial autocorrelation and non-stationarity can yield improved predictions compared to ordinary LSR modeling of WNV risk. A spatially explicit modeling technique may be useful in policy-making and decision-making depending on the granularity and resolution of available data. Identifying the spatial variations in relationships by estimating local regression parameters allows the spatial distribution and interaction of predictor variables to be explored. Analyzing local variations in relationships provides those concerned with public health policy the ability to target resources and to achieve improved outcomes through location-specific activities (Comber, Brunsdon & Radburn, 2011) because spatial heterogeneity can improve predictions by capturing geographic shifts in the ecological drivers (Wimberly, Baer & Yabsley, 2008). While environmental data used in this research were of fine resolution, WNV disease human incidence data and infected dead bird data that is used is available only at a coarser county scale. We had to assume that aggregating the environmental data up to the county adequately represented the environmental conditions presented in the county, but we knew that data aggregation was likely to introduce some uncertainty into the model. The dead bird model applied in this study can be used for better understanding of WNV risk and the techniques used could be replicated at finer spatial scales thus leading to better intervention efforts.

In summary, WNV, a globally emerging infectious disease, was found to be heterogeneously related to environmental factors at the county level throughout California during the time that our data were collected. Our findings may assist those conducting risk assessments for WNV transmission in local areas by helping local public health entities allocate resources and improve preparedness for an outbreak according to region-specific conditions.