Adjusting Second Moment Bias in Eigenspace Using Bayesian Empirical Estimators, Dirichlet Tessellations and Worldview I Data for Predicting Culex Quinquefasciatus Habitats in Trinidad Models for West Nile Vectors

Temporally weighted regression models with a spatial autoregressive component may estimate nonlinearities in spatiotemporal-sampled data of Culex quinquefasciatus, a major vector of West Nile Virus (WNV) which can help implement control strategies by determining optimal predictors associated to prolific habitats. The design of this kind of mixed model can specifically incorporate spatial autocorrelation whilst including the influence of other aspatial predictor variables. Currently, the lack of an estimation theory that allows for het-eroscedasticity and corresponding joint hypothesis testing in the presence of spatial dependence in georefer-enced Cx. quinquefasciatus habitat data is a serious shortcoming in WNV research. In this paper we used spatially lagged and simultaneous autoregressive models based on multiple predictor variables of immature Cx. quinquefasciatus and Worldview 1 (WV-1) data to help implant a remote habitat-based surveillance system in Trinidad. Initially, we used Geomatica Ortho Engine ® v. 10.2 for extracting a Digital Elevation Model (DEM) from the WV-1 raw imagery. Results of the DEM analyses indicated a statistically significant inverse linear relationship between total sampled Cx. quinquefasciatus data and elevation (m) (R 2 = −0.439; p < 0.0001), with a standard deviation of 10.41. Additional field-sampled information was derived using data from an orthogonal grid-matrix constructed in an ArcInfo 9.3 ® and overlaid onto the WV-1 data. A unique identifier was placed in the centroid of each grid cell. Univariate statistics and Poisson regression models were then generated using the georeferenced covariates in SAS/GIS ®. Coefficient estimates were also used to define expectations for prior distributions in a Bayesian estimation matrix using Markov Chain Monte Carlo (MCMC) specifications. A spatial residual trend analyses was then performed using autocorrelation indices which linked tabular data in SAS PROCLMIXED ® with the egg-raft count data in ArcInfo ®. The estimation matrix identified prolific habitats based on the covariate distance to the nearest house. An Ordinary kriged-based interpolator was then constructed in Geostatistical Analyst Extension of ArcGIS 9.3 ® based on the adjusted Bayesian estimates. For total Cx. quinquefasciatus egg-raft count, first order trend was fitted to the semivariogram at a partial sill of 5.931 km, nugget of 6.374 km, lag size of 7.184 km, and a range of 31.02 km using 12 lags. We assessed the performance accuracy of the interpolation procedures based on the magnitude and distribution of errors between observed and model-predicted values using Voroni tessella-tions. These residuals divided the space between the individual georeferenced Cx. quinquefasciatus habitats …


Introduction
West Nile virus (WNV) is among the Culex-borne encephalitides in the Flavivirus genus.The virus cycles between birds and mosquitoes, while horses, humans, and a number of other vertebrates are considered incidental hosts.Since WNV first appeared in the Western Hemisphere in New York in 1999, it has spread rapidly across the North American continent, causing large numbers of human cases with neurologic disease and death and even greater amounts of milder disease characterized principally by fever and rash [1].West Nile Virus has been spreading southward into the Caribbean Basin and Latin America as well, where its public health impact remains poorly understood [2].The virus was first detected in 2001, in Jamaica and the Cayman Islands.Thereafter, serological evidence of WNV transmission has been accumulating.Cross-reactive WNV antibodies in humans have been detected from Mexico, the Bahamas, and Cuba, in horses from Guadalupe, Mexico, Central America, Puerto Rico, Colombia, and also in resident birds from the Dominican Republic and Venezuela.In November 2006, 4 serologically confirmed human WNV encephalitis cases were reported in Argentina.West Nile Virus may be following the same patterns of southward dissemination as other arboviruses into the Caribbean and Central and South America via migratory birds.
There have been attempts to develop and implement WNV surveillance in the Caribbean but, unfortunately present systems in the region are unprepared to track WNV spread.For example, in Guadeloupe, extensive surveillance for WNV infection demonstrated a high seroprevalance in equines in 2003, despite no human or equine case being reported.A protocol for epidemiological surveillance for WNV has been developed through a collaborative electronic working tool posted on the Caribbean Animal Health Network (CaribVET) website (www.caribvet.net).CaribVET is a collaboration of veterinary services, diagnostic laboratories, research institutes, universities, and regional/international organization to improve animal health in the Caribbean through exchange of information and data management.Improving collaborations and information dissemination for understanding WNV epidemiology is vital; however, WNV surveillance and diagnostics in the Caribbean requires tools for regional data analyses.Unfortunately, to date there have been virtually no significant findings in any established WNV surveillance programs in tropical America.In september 2007, WNV infection detected by polymerase chain reaction (PCR) in post-mortem brain tissue taken from an encephalitic horse and by virus isolation from a dead bird confirmed WNV transmission in southwest Puerto Rico but, since then there have been no isolates reported from other areas in Central and South America or the Caribbean.
The explanation for the paucity in WNV detection in tropical America may be the lack of active WNV surveillance.One would hypothesize that the likelihood of identifying WNV cases would be maximized in the Caribbean Basin, given its location in the flight paths of infected migratory birds.However, due to cross-immunity of WNV with other tropically endemic Flaviruses such as St. Louis encephalitis, it is possible that WNV has not been able to establish persistent enzootic foci in the region, so that reported serological evidence may be caused by frequent but unsuccessful virus reintroductions.Regardless, presently a WNV mosquito habitat-based intervention surveillance system using a vigorous quantitative assessment of field and remote-sampled data is warranted for countries of the Caribbean.
A main challenge in devising an effective West Nile surveillance program in regions of the Caribbean; however, is the lack of sufficient empirical knowledge on the spatiotemporal patterns of Culex quinquefasciatus larval habitats and their underlying contribution to the adult population.The main tropical vector for WNV in birds and people is Cx.quinquefasciatus in the Southeastern United States and the Caribbean [3].Previous research has revealed that Cx. quinquefasciatus abundance will vary significantly in different areas due to hydrologic regime and the proximity to habitats and blood-meal hosts [4][5][6][7].These factors are not evenly distributed spatiotemporally, so it is logical to expect a variation in abundance and distribution of virus-positive mosquitoes within different ecosystems in the Caribbean.
Recently, the use of remote sensing (RS) and GIS for the quantitative prediction of geographical distributions of agent, reservoir, and vector species has been advocated to augment traditional WNV mapping methods [4,5,8].For example, the United States Geological Survey (USGS) and Center of Disease Control (CDC) have employed GIS and RS to prepare interpretive maps showing WNV activity in many regions of North America (http://www.cdc.gov/ncidod/dvbid/westnile/resources/wnvguidelines1999.pdf).Common GIS/RS technology for WNV data analyses include: 1) time series overlay analyses of thematic data and spatial intersection, 2) buffer generation and neighborhood analysis, 3) vector-based grid generation and network analysis; and, 4) surface modeling [9].Once an area has been imaged, associations between risk variables and environmentalsampled covariates can be quantified using the spatial data analysis capabilities of GIS for implementing a WNV mosquito surveillance system [9][10][11].
Quantification of vector-host interactions incorporated on Land Use Land Cover (LULC) maps generated in GIS time series overlay analyses suggest that Culex habitat productivity vary greatly over short distances, depending on the space and time at which field measurements are made [5,[12][13][14][15].For example, Jacob et al. (2009a) generated a land use land cover (LULC) classification in Erdas Imagine based on Landsat-7 ETM+ data acquired July 2003 and Landsat-5 TM data acquired July 1991 for examining two ornithophillic mosquito vector species Culex pipiens and Culex restuans in residential areas of Urbana-Champaign, Illinois [4].A maximum likelihood unsupervised classification and LULC detection was performed using a cross-tabulation detection method and a grid-based algorithm stratified by drainage for detection of environmental predictors associated to field-sampled egg-raft abundance count data.The resulting LULC change matrix revealed that between 1991 and 2003 there was a total of 12.1% LULC change in the study site.Of 20,166 egg-rafts sampled, 76.7% was in land cover classified as maintained urban followed by 14.8% in maintained non-urban.The well-drained stratum harbored 56.4% of the total sampled data.Additionally, a remote stratification of the urban land cover site based on QuickBird visible and near infrared (NIR) data revealed that 83.3% of the egg-raft distribution was in LULC classified as residential high canopy coverage.In another study, Jacob et al. (2009b) constructed multiple geospatial models based on LULC and meteorological variables to predict adult abundance of the same Culex species in Des Planes mosquito abatement district in northern Cook County, Illinois [5].An ecological database of the study site was constructed in ArcGIS 9.2 ® also using multitemporal QuickBird visible and NIR data overlaid with total adult Culex abundance and WNV infection rates from 15 georeferenced gravid traps sampled weekly from May 2002 to October 2005.A regression analyses of LULC parameters revealed that adult Culex abundance in the georeferenced gravid traps was positively associated with temperature and negatively associated with precipitation.In 2002, the overall accuracy of the model was 73.3%, 65.1% in 2003, 47.2% in 2004 and 60.0% in 2005.Temperature and precipitation explained > 80.1% of the variation in each of the models.Estimates from the multivariate analysis also revealed that land cover that was classified as forest LULC along with middle range built environment (< 40.1%) predicted heightened infected vector activity.Because abiotic factors such as LULC are measurable in GIS using submeter resolution data, these pre-epidemic abiotic signatures may be quantified for predicting arboviral epidemic transmission in isolated areas of the Caribbean.These epidemic signatures may be tracked, measured, quantified, and presented visually to help forecast human arboviral epidemics.While remote models of the observed spread of WNV have been used to evaluate potential vector and host mechanisms of dispersal and to define the suitability of the various landscape features and meteorological variables for mosquito transmission of WNV in many portions of the US, there is no literature using LULC and high resolution satellite data for evaluating ecologically sampled Cx. quinquefasciatus habitats in any region of the Caribbean.
Recently, digital elevation models (DEMs) have become popular with the use of GIS for surface modeling flood and swamp water mosquito abundance in many different geographic regions [16][17][18].A digital elevation model (DEM) is a digital file consisting of terrain elevations for ground positions at regularly spaced horizontal intervals.It is a digital representation of ground surface topography which is usually represented as a raster or as a triangular irregular network.DEMs include a number of production strategies for determining terrain-related parameters associated to prolific Cx. quinquefasciatus habitats such as manual profiling from photogrammetric stereo models, digitizing of contours, digitizing topographic map contour plates, converting hypsographic and hydrographic tagged vector files and also for performing autocorrelation via automated photogrammetric systems.To obtain DEMs from satellite images for determining vector amplification and host population dynamics for mapping the spatiotemporal distribution of WNV from high resolution data two methods are possible: along-track stereoscopy from the same orbit using fore and aft images and across-track stereoscopy from two adjacent orbits [16].The four possibilities of stereo image processing for DEM generation are Fore-Aft, Fore-Nadir, Aft-Nadir and Fore-Aft-Nadir images as a pair [19].The simultaneous acquisition of along-track stereo data has a strong advantage in terms of radiometric variation versus the multi-date acquisition of across-track stereo data which has been traditionally employed in past research for identifying terrain-related explanatory covariates of floodwater mosquito habitats [20].The across-track approach has been applied frequently since 1980, first with Landsat from two adjacent orbits, then with SPOT using across-track steering capabilities and finally with IRS-I C/D.But these satellite data were limited due to spatial resolution.Fortunately, along-track stereoscopy has recently gained renewed popularity using higher resolution satellite data systems including JERS-i's optical sensor(OPS), German Modular Opto-electronic, Multi-Spectral Stereo-Scanner (MOMS), ASTER, IKON-OS, QuickBird, Orthoview, SPOT-5, FORMOSAT II, Cantosat and the latest addition of Digitalglobe's Worldview (WV)-1.These data can define georeferenced terrain covariates associated with prolific Cx. quinquefasciatus habitats.For example, Jacob et al. (2010) used multispectral QuickBird data and DEM-based GIS methods to evaluate stream flow direction for determining satellite-derived characteristics of drainage basin physiographic predictors associated with Culex erraticus, a mosquito vector of eastern equine encephalitis virus (EEEV) and Northern Cardinals, an avian host associated with WNV and EEEV in Tuskegee, Alabama [16].Eastern equine encephalitis virus is a virulent arbovirus maintained in an enzootic cycle between local avian fauna and ornithophilic mosquitoes in coastal regions of the eastern United States [1,14,16].Euclidian distanceto-nearest hydrological body was calculated as the distance in geo-space from a grid cell to a stream grid cell defined by a Stream Raster Grid.A Terrain Analysis using a DEM (i.e., TauDEM) was then performed in ArcGIS 9.2 ® to retrieve multiple geomorphological parameters.Additionally, a three-dimensional model of the study area was constructed based on the DEM using ArcScene extension of ArcGIS ® .A Pearson correlation analyses computed the ratio of covariance between the sampled variables to the product of their standard deviations.Substituting estimates of the covariances and variances based on the Cx.erraticus samples provided the sample correlation coefficient which revealed significant inverse linear relationships between the total field-sampled adult count data and elevation (R 2 = −.4711;p < .0001)with a standard deviation (SD) of 11.16.The analyses of the Northern Cardinal data and elevation also revealed a similar statistical relationship (R 2 = −.5831;p < .0001)with a SD of 11.42.Even though the usefulness of selected static catchment predictor covariates (e.g., stream, density, and slope) in assessing Culex habitats has been verified in previous research, characteristics of drainage networks and basin physiographic parameters have not been applied in modeling WNV mosquito abundance and distribution in any country of the Caribbean.
Predictive interpolation algorithms have also not been used for forecasting productive Cx. quinquefasciatus habitats based on field and remote-sampled data in any country of the Caribbean.One of the most powerful and potentially valuable GIS/RS vector mosquito data modelling mechanisms is that of automated prediction.Generally, given relevant spatiotemporal field and re-mote-sampled explanatory predictor variables with appropriate annotation, a WNV mosquito model can be constructed for making accurate forecasts of one or more dependent "target" quantities of interest.These quantities of interest may be continuous real values (a scenario often referred to as regression" or "interpolation"), in vector mosquito data analyses or, it can be in the format of discrete labels ("aquatic habitat classification" or "habitat distribution pattern recognition").Both these cases are referred to as "supervised learning" in the GIS/RS vector mosquito data learning vernacular.Predictive interpolation algorithms can also map urban variation in vectors based on differentials in environmental-sampled parameters which can identify and validate models linking human ecological factors, while controlling for potential spatial autocorrelation in model residuals [21][22][23][24][25]. Spatial autocorrelation, the correlation among values of a single variable strictly attributable to their relatively close positions on a two-dimensional surface is frequently encountered in the analyses of spatiotemporal georeferenced vector mosquito data [4,8,21,26,27].Autocorrelation violates the ordinary least squares (OLS) assumption that the error terms are uncorrelated.While it does not bias the OLS coefficient estimates in a predictive vector mosquito habitat distribution model, the standard errors tend to be underestimated and the t-scores overestimated when the autocorrelations of the errors at low lags are positive [16].
Kriged-based interpolation may analyze spatiotemporal-sampled WNV mosquito data in the presence of fine-scale autocorrelation.Kriging is mathematically related to regression analyses.Both theories derive a best linear unbiased estimator based on assumptions of covariance and make use of the Gaussian-Markov theorem to prove independence amongst georeferenced vector mosquito data exploratory parameters [26].Kriging vector mosquito data is commonly made by interpolation of a single realization of a randomized study area based on field and remote observations in multivariate environmental-sampled datasets.Under the regionalized assumption, local habitat variation is not generally unstructured but, it is spatially dependent at some scale as georeferenced points within a given habitat distance apart depend on one another in a statistical sense (16) Any variable is stationary when its distribution is invariant under translation, that is, the stationary region is the location to make estimates (26).Stationary variables provide a preliminary model decision for uncertainty about the sample value at a specific sampled habitat location z(x) (16).For any space lag (h) increment in a particular direction, the distribution of z(x 1 ), z(x 2 ) ... z(x n ) habitat based observations is the same as that of z(x 1 + h), z(x 2 + h)...z(x n + h), where z(x) is regarded as the sum of three terms: (1) mean, drift or trend, (2) random dependent variable plus (3) error (16) This, however, requires the first two moments (i.e., mean and variance) to be constant.This spatial homogeneous region in the distribution of the georferenced Cx. quinquefasciatus habitats then will make random functions homogeneous and self-repeating in space.The variance of the random spatial variable will equal covariance at zero lag habitat distance (σ 2 = C (0), while the covariance between z(x) and z(x + h) becomes σ x,x+h = σ 2 f(h) = C (h) (16).This means that the stationary covariance only depends on the sampled habitat distance and direction between two sampled habitats and not on their georeferenced locations, or the second order stationary assumption.Hence, the expected difference between sampled values at two georferenced Cx. quinquefasciatus habitat locations separated by distance h is zero.If E (μ(x)) = m then z(x) = m + e(x) and E (z(x) -z(x + h)) = 0 (23).Covariance at h lag habitat distance then will be E{(z(x) -m)(z(x + h)-m)} = E(z(x) x z(x + h)) -m 2 = Cov(h), being that z(x) the value at x sampled habitat site, z(x + h) the value at x sampled habitat site plus h lag habitat distance is based on m representing the overall mean and E() symbolizing the expected value (16).Although the covariance does not exist theoretically in a predictive vector mosquito habitat distribution model, the variance of the difference between two sampled habitat locations will be assumed to be present, as it does not depend on the georeferenced location and Var (z(x + h) -z(x)) = E({z(x + h) -z(x)} = 2γ(h), being E(z(x) -z(x + h) = 0.That is, differences between sampled habitat sites are merely a displacement function between them (i.e, the intrinsic hypothesis).The γ(h) (i.e., dissimilarity average function) will be the key tool for the structural interpretation of the predicted Cx.Quinquefasciatsus habitats as well as for uncertainty estimation parameters.Mathematically, the variogram is defined as γ(h) = 0.5xVar(z(x + h) -z(x) (26).Since, the mean of z(x + h) -z(x) is zero for intrinsic variables then the variance is just the mean square difference.Consequently, the variogram for a predictive vector habitat distribution model is as follows: γ(h) = 0.5xE{(z(x + h)z(x)) 2 }(16).These predictive interpolation vector mosquito habitat distribution models can measure unknown local mean habitat values as having a local linear or quadratic trend by combining regression of the dependent variable (e.g., larval/pupal count data) based on auxiliary variables, such as terrain-related parameters (i.e., slope), RS imagery, and thematic maps, with kriged residuals where auxiliary predictors are used directly to solve the weights.For instance, in our previous example of research conducted in Tuskegee, Alabama for examining vector-host activities of Culex erraticus, a krigedbased interpolator was also constructed in Geo-statistical Analyst Extension of ArcGIS 9.2 ® for predicting prolific habitats based on adult abundance count data [16].For total adult Cx. erraticus count, a first-order trend ordinary kriging process was fitted to the semivariogram at a partial sill of 5.764 km, nugget of 6.114 km, lag size of 7.472 km, with a range of 32.62 km using 12 lags, suggesting that this species forage several kilometers beyond their habitats.The values that the semivariogram model attains at the range (i.e., the value on the y-axis) is called the sill, while the nugget or nugget variance, C 0 , occurs when a curve is fitted to the set of points that lie within the range and sill and the model used has a non-zero intercept [28].
Furthermore, to improve interpolation accuracy in a kriged-based algorithm, predictive mean standard error distributions can be quantified for attaining asymptotically optimal predictors.Traditional kriging is motivated by an expected squared prediction error from a stochastic model [28].Presently many algorthmitic paradigms of computational geometry have been applied efficiently in Geomatics and GIS using Dirichlet/Voronoi/Thiessen tessellation techniques for assessing interpolation models for quantifying the maximum error variance and the extended prediction variance [28].A Voronoi diagram is a special kind of decomposition of a metric space determined by distances to a specified discrete set of objects in the space, (e.g., dataset of Culex habitat ground reference coordinates).Generalization towards dynamic operations and kinematic motion have been investigated for about a decade using Dirichlet neighborhood analyses.Voroni residuals are amenable to assessment of goodness-of-fit tests in a predictive vector mosquito habitat distribution models as they are approximately Gamma distributed and tend to be far less skewed than Pearson's residuals and other residual diagnostic tests [16].Such tessellations have been used in grid generation in volumes and surfaces, data compression, image segmentation and edge detection, optimal allocation of resources, optimal quadrature rules, optimal representation, quantization, finite difference tests and finite volume schemes.Unfortunately, there is no published literature using interpolation algorithms or Dirichlet tessellations for quantifying spatial error and predicting habitat locations using georeferenced spatiotemporal-sampled predictor variables in WNV mosquito research for any country of the Caribbean.
In this research, we generated multiple predictive geospatial models using field and remote-sampled explanatory variables of Cx. quinquefasciatus habitats in Trinidad.Elucidation of variability in mosquito productivity requires spatial accounts of oviposition processes [29].Indeed, this is important because of the recent identification of sero-positive horses and domestic birds in Trinidad [2].Presently there is no remote habitat-based intervention for early detection of WNV activity in any regions of the Caribbean or Central and South America.A remote habitat-based intervention may provide early detection of arboviral activity by providing useful information for projection of the potential course of WNV transmission and also other arthropod-born viral diseases (e.g., St Louis encephalitis, Venezuelan equine encephalitis) later in the season while unraveling the environmental factors associated with the maintenance and establishment of enzootic cycles in Trinidad and in other countries of the neotropics.
For mosquito larval control to be effective, however, it may be crucial to quantify individual sampled habitat variation to maximize control efforts.For example, the importance of variation in mosquito production among breeding sites in design of control programs is now broadly accepted in suppressing domestic containerbreeding Aedes aegypti, a mosquito vector of Dengue in Trinidad.This research involved quantifying population densities of A. aegypti in four towns using standard house-to-house inspections of all water-holding containers to determine whether persistently positive containers and premises existed over a three-month period in the wet season, from May to July 2002.From a total of 1,503 houses inspected, 15% were positive with 3% persistently positive over the three month period and classified as 'key premises'.From a total of 24,439 containers inspected, 1.3% were positive for A. aegypti larvae and pupae.Of 17 container types inspected, only water drums (54%), buckets (22%), tubs and basins (8%), water tanks (5%), brick holes (4%), and tires (2%) were significant larval habitats (P < 0.001) producers.The role that key premises play in the introduction and reinfestation of A. aegypti-free communities was determined.The results suggested that A. aegypti control programs could be more cost effective and sustainable by targeting efforts on key premises and key containers to control mosquito densities and dengue transmission while reducing manpower needs and insecticide use.
In this research, we constructed an autocovariate error matrix in SAS/GIS ® for spatially targeting estimates generated from individual sampled Cx. quinquefasciatus habitat data.In our Bayesian model we did not assume that the georeferenced predictor variables were conditionally independent, as is commonly the case in hierarchical generalized linear model (HGLM) construction practices for predicting WNV mosquito habitats.In our Bayesian network we assumed that the variables were conditionally independent.In practice, this assumption may not hold for WNV mosquito habitat models and may give rise to incorrect inferences.Violation of the assumption of statistical independence will give rise to incorrect model specifications [28].Thus, we proposed the creation of a hidden node for modeling the dependency in the spatiotemporal-sampled parameters.In order to determine the conditional probability matrices for the hidden node, we used a gradient descent method.The objective function to be minimized was the squared-error between the measured and computed values of the instantiated nodes.Both forward and backward propagation were used to compute the node probabilities for quantifying the dependence in the data.We also restricted our attention to the simultaneous autoregressive Gaussian spatial process and the autoregressive Gaussian response model (i.e., the spatial lag model), for quantifying latent autocorrelation error components in the Bayesian uncertainty estimates.We assumed that observed parameter error estimation patterns in the response variable decomposed into three statistically independent components: 1) a systematic spatial trend component that could be specified by a parsimonious set of exogenous variables; 2) a stochastic signal that reflects either an underlying spatial process and/or a set of missing exogenous factors with an inherent spatial pattern; and, 3) independent white-noise disturbances.A specific subset of eigenvectors was then used to determine local and regional variation in the field and remote-sampled covariates by capturing dependencies among disturbances generated from the Bayesian probabilistic model.Eigenvectors are a special set of vectors associated with a linear system of equations (i.e., a matrix equation) [4,5,8,10].
We constructed our remote habitat-based intervention by adoption of a landscape approach to elucidate mechanisms underlying Cx.quinquefasciatus productivity.Our landscape perspective was based on multiple georeferenced explanatory variables of individual sampled habitats to determine their role in WNV disease transmission in Trinidad.Arguably, our framework evolved from the traditional concept of regression models generated from LULC covariates for identification of hot zones (i.e., cluster), but our remote habitat-based intervention also focused on critical spatial autoregressive modifications of the georeferenced predictors.We developed our models to emphasize the role of oviposition foraging in Cx. quinquefasciatus in Trinidad.Habitat-based interventions should emphasize the link between foraging behaviors of egg-laying mosquitoes and the availability of breeding sites in evaluation of environmental management programs [29].The basic assumption in our models was that heterogeneity in immature productivity was substantially large among Cx.quinuefaciatus habitats and that the number of productive habitats was relatively fewer than that of unproductive ones in focal areas.Untargeted interventions are inefficient because habitats are randomly selected without consideration of the distribution of production (29).Spatially targeting prolific breeding habitats using a remote habitat-based intervention may reduce Cx.quinuefaciatus populations in Trinidad and help implement a cost-effective surveillance mechanism.Therefore, our objectives of this research were to: 1) construct a DEM to identify terrain-related explanatory covariates associated with spatiotemporal-sampled mosquito data; 2) generate regression models to determine statistical significance associated with the sampled covariates, 3) quantify latent autocorrelation components in model residual variance estimates generated from an eigenfunction decomposition spatial filter algorithm; and, 4) develop spatial interpolation models of potential prolific habitats site based on field-sampled count data.Although the discussion is centered on WNV vectors specifically in Trinidad, the framework and derived guidelines in this research may be applicable to integrated control programs for other mosquito species and insect-borne diseases in other regions.

Study Site
The island of Trinidad is the most southerly part of the Lesser Antilles located at approximately 10° north latitude and 15 km from the eastern coast of Venezuela.The geology, flora and fauna are typical of the adjacent South American mainland [31,32].Trinidad is approximately 80 km from north to south and is on average 60 km wide, roughly rectangular in shape with large promontories on the northwest and southwest of the country.The climate of Trinidad is tropical with a "rainy" season between May and November and a "dry" season from December to May.Details about climate, physical features, vegetation, and population size have been previously described [33,34].

Mosquito Egg Sampling
The Cx. quinquefasciatus habitat population was surveyed by counting the number of egg-rafts deposited in oviposition traps.Field sampling was conducted from July 2008 to July 2009.One hundred and fifty two temporary, permanent, and semi-permanent mosquito larval habitats in the study site were mapped and classified, and recalibrated using a CSI-Wireless differential corrected global positioning systems (DGPS) Max receiver which yields a positional error of .179m (+/− 392 m) [35].Water bodies were inspected for Cx.quinquefasciatus eggrafts.Oviposition traps were located in shady areas around residences and municipal buildings, as well as the edge of parks and woodlots.

Habitat Base Mapping
Initially, Landsat images were collected using SLC-off mode (http://edc.usgs.gov/#Find_Data/Products_and_Data_Available/ETM).We used the USGS Global Visualization Viewer (GloVis) system (http://glovis.usgs.gov/) to search for images for performing a land cover classification.Path 233/row 53 covered the Trinidad study area.Unfortunately, the images acquired from GloVis, were all cloudy, due to the climate conditions of the island.Cloudy image is an issue for remote classification and change detection for identification of mosquito habitats [24].We did find three Landsat 7 passes of path 233/row 53 which were less than 10% cloud free.Theses collection periods were April 30, May 16, and June 1, 2005.These data; however, were not able to discriminate the land cover due to poor resolution (i.e., 30m resolution) associated with the overlaid georeferenced Cx. quinquefasciatus habitats coordinates.
We acquired WV-1 data (www.digitalglobe.com)for the study site on July 2010.Operating at an altitude of 496 kilometers, WV-1 has an average revisit time of 1.7 days and is capable of collecting up to 750,000 km 2 per day of half-meter resolution imagery.The satellite is equipped with state-of-the-art geolocation accuracy capabilities and exhibits stunning agility with rapid targeting and efficient in-track stereo collection.Various collection sizes of surveying are available: frame mode, route survey (along coastline, roads and linear objects), areal survey (60 × 60 km), and stereoscopic areas on a single pass as well (www.digitalglobe.com).
In this research, the WV-1 imagery was radiometrically and sensor corrected but was not projected to a plane using a map projection or datum.The sensor correction blended the pixels from all the detectors into synthetic arrays to form the image of the study site.The geolocation accuracy specification of the data was 4.0 to 5.5 m CE 90 at nadir excluding terrain and off-nadir effects.The imagery was delivered at full-swath cut into 14 km lengths.The WV-1 imagery was classified using the Iterative Self-Organizing Data Analysis Technique (ISODATA) unsupervised routine in ERDAS Imagine v.8.7™ (Atlanta, USA) which is commonly used to classify land covers associated with disease foci [18,26,29].The images were co-registered manually using ground control points (GCPs) and by applying a first order polynomial algorithm with a nearest neighbor re-sampling method.The Universal Transverse Mercator (UTM) Zone 20P datum WGS-84 projection was used for all of the spatial datasets.
Base maps were prepared from the field-sampled Cx.
quinquefasciatus habitat data and WV-1 imagery for the study site using the DGPS ground coordinates.We generated GIS data layers in the vector shapefile format, ArcInfo ® Coverage format and raster GRID format for use in the GIS software package ArcGIS v 9.3 ® .We constructed a discrete tessellation of the region using a digitized grid-based algorithm (Figure 1).For remote identification of vector mosquito habitats the first step is often to construct a discrete tessellation of the region [36].A unique identifier was placed in each grid cell.Each Cx. quinquefasciatus habitat/polygon was assigned a unique identifier.Field attribute tables were then associated to the polygons.

Environmental Data Analyses
Base maps were also generated using the spatiotemporalsampled Cx. quinquefasciatus habitat data which was stored in a Vector Control Management System ® (VCMS) (Clarke Mosquito Control Products, Inc. 159 N. Garden Avenue, Roselle, IL 60172) database.The VCMS supported the mobile field data acquisition in the study site through a PocketPC™.All two-way, remote synchronization of data, geocoding and spatial display were processed using the embedded GIS Interface Kit™ that was built using ESRI's MapObjects™ 2 technology.The VCMS database plotted and updated the DGPS ground coordinates of the Cx.quinquefasciatus habitats and supported exporting data in spatial format in which any combination of the sampled habitats and explanatory covariates were described in ESRI shapefile format for use in GIS.The database displayed this information onto a user-defined field base map.

Spatial Hydrological Model
The latest version of PCI Geomatics Orthoengine ® software was used to generate a DEM using the WV-1 data and the georeferenced Cx. quinquefasciatus habitat covariates (Figure 2).OrthoEngine offers an industryleading variety of control sources including, manual entry, geocoded imagery, geocoded vectors, chip database, digitizing tablet, or a text file (www.pcigeomatics.com).We used the DGPS coordinates to orient the images to our map coordinate system.We used PCI Geomatics OrthoEngine to extract the DEM.PCI software supports automatic GCP/tie point (TP) collection, using Toutins rigorous models, Rational Polynomial Coefficient (RPC) model construction, automatic DEM generation, orthorectification and automatic mosciaking (www.pcigeomatics.com).
WV-1 stereo pairs were supplied by the full scene of the Trinidad study site using Basic 1B level data de-signed for the creation of the DEM.A set of WV-1 Basic 1B stereo images were provided by Digitalglobe.The set of stereo data was in Basic 1B OR 2 A product format.Due to the 2 gigabyte limit of the TIFF data, WV-1 data was distributed in tile format.Each tile came with its own RPC file.Because our OR 2 A format was map projected, the effects of any high frequency movements had already been removed from the scene resulting in a good fit for the RPCs.
In order to leverage the WV-1 image for attaining topographic indices from the Cx.quinquefasciatus data, a geometric model was required.We generated an RPC model.The RPC model has proven to be the most popular geometric model for high resolution images (www.digitalglobe.com).Since bias or errors may still have been present in the RPCs, our results were preprocessed with a polynomial adjustment using several accurate GCPs (Table 1).A zero order RPC adjustment is adequate to obtain an RPC model accuracy within 1m (www.digitalglobe.com).
To generate a DEM for determining the terrain-related Cx. quinquefasciatus habitat parameters in the Trinidad study site, a pair of quasi-epipolar images were generated from the fore and aft stereo images to retain elevation parallax in the X-direction only.An automated image match procedure was then used to produce the DEM through a comparison of the respective gray values of these images.To find the corresponding pixels in the left and right quasi-epipolar images, we used a hierarchical sub-pixel mean normalized cross-correction matching method.The actual matching method employed generated correlation coefficients between 0 and 1 for each matched pixel, where 0 represented a total mismatch and 1 represented a perfect match.A second order surface was then fitted around the maximum correlation coefficients.The difference in the georeferenced Cx. quinquefasciatus habitats between the images gave the disparity  or parrallex arising from terrain relief which was then converted to absolute elevation values above the WGS-84 absolute ellipsoid using a 3-d space interaction solution.Digital elevation model statistics were then generated employing the Pearson's correlation to find if an association existed between the georeferenced Cx.Quinuefasciatus habitats and any of the sampled covariates (Table 2).

Regression Analyses
Variable selection for the multiple regression models was carried out by a combination of automatic (stepwise) procedures and goodness-of-fit criteria, and by selecting cluster covariate measurements that explained WNV prevalence.Routinely, the Analysis of Variance (ANOVA) are used to provide calculations about spatio-temporal sampled WNV habitat data regarding levels of variability within a regression model while forming a basis for tests of significance for the parameter estimates [18].Analysis of variance is a collection of statistical models and their associated procedures in which the observed variance is partitioned into components due to different sources of variation.
In this research, the regression line concept we used was: ( ) ( ) ( ) , where the first term was the total variation in the response y (egg-raft count of Cx. quinquefasciatus habitats ), the second term was the variation in mean response based on the sampled parameters and the third term was the residual value in the model estimates.Squaring each of these terms and adding over all of the sampled observations gave the equation  mate of σ² for determining whether or not the null hypothesis is true [25].The null hypothesis for the ANO-VA analyses was based on the average value of the dependent variable (i.e., egg-raft count of Cx.Quinquefasciatus).The ANOVA calculations for multiple regression are nearly identical to the calculations for simple linear regression, except that the degrees of freedom are adjusted to reflect the number of explanatory variables included in the model [25].
For the spatiotemporal-sampled explanatory variables, (p) the DFM was equal to p and the error degrees of freedom (DFE), which was equal to (n -p -1) and the total degrees of freedom (DFT) which was equal to (n -1), the sum of DFM and DFE.The corresponding ANOVA table generated from the field and remotesampled Cx. quinquefasciatus habitat parameters is shown below: In the multiple regression analyses, the test statistic MSM/MSE had an F (p, n -p -1) distribution.In this research, the null hypothesis stated that , and the alternative hypothesis simply stated is that at least one of the sampled Cx. quinquefasciatus habitat parameters j 0, j 1, 2, p β ≠

=
 .Large values of the test statistic can provide evidence against the null hypothesis [25].The F test did not indicate which of the parameters j β ≠ was not equal to zero, only that at least one of them was linearly related to the response variable (i.e., total egg-raft count of Cx.Quinquefasciatus habitat parameters).The ratio SSM/SST = R 2 (i.e., squared multiple correlation coefficient) was the proportion of the variation in the response variable that was explained by the field and remote-sampled Cx. quinquefasciatus habitat data.The square root of R² (i.e., the multiple correlation coefficient) was the correlation between the sampled observations y i and the fitted values ˆi y .The relationship between WNV prevalence and each individual potential predictor variable sampled in the Trinidad study site was also investigated by single variable regression analysis using PROC MIXED.Since parasite prevalence data are binomial fractions, a regression model was used, as is standard practice for the analysis of the spatiotemporal-sampled WNV mosquito data [4].Poisson regression analyses were also constructed in PROC MIXED to determine the relationship between Cx. quinquefasciatus habitat egg-raft count data and the sampled habitat characteristics.These regression models were generalized linear models with the logarithm as the (canonical) link function, and a Poisson distribution function.
The regression analyses assumed independent counts (i.e., n i ), taken at sampled habitats I = 1, 2… n.The fieldsampled count data were described by a set of explanatory variables denoted by matrix X i , a 1 × p vector of covariate measurements for a sampled habitat i.The expected value of these data was given by μ i (X i ) = n i (X i ) exp (X i β), where β was the vector of non-redundant parameters and the Poisson rates parameter was given by λ i (X i ) = μ i (X i )/n i (X i ) [37].The rates parameter λ i (X i ) was both the mean and the variance of the Poisson distribution for each sampled Cx. quinquefasciatus habitat location i.The dependent variable was total field-sampled egg-raft count data.The Poisson regression model assumed that the data was equally dispersed, meaning that the conditional variance equaled the condition mean.The procedure used maximum likelihood estimation to find the regression coefficients.The data was log-transformed before analyses to normalize the distribution and minimize standard error.The predictor variables used in the regression analyses are listed in Table 3.

Bayesian Analyses
In the Bayes formulation, the specification of the Cx.quinquefasciatus habitat model was completed by assigning priors to all unknown parameters.We used our dataset of spatiotemporal-sampled Cx. quinquefasciatus observations X = [x 1 , , x n ] where each x i for i = 1  n was assumed to be distributed according to some distribution p(x i |θ).In this research, θ was a parameter that was unknown and had to be inferred from the data.Our Bayesian procedure began by assuming that θ was distributed according to some prior distribution p (θ|α), where the parameter α was a hyperparameter.The joint of the Cx.quinuefasciatus data was then generated using: and ( ) ( ) were conditionally independent of the hyperparameter.Bayesian inference then determined the posterior distribution of the parameter ( ) For the fixed regression parameters, a suitable choice was the diffuse prior, i.e., p (γ) ∝ const, but a weakly informative Gaussian prior was also possible.A secondorder Gaussian random walk prior was used to allow enough flexibility while penalizing abrupt changes in the function [17].The prior was expressed as: τ was the variance, with diffuse priors 1 f ∝ const, 2 f ∝ const.For additional spatial effects in the sampled Cx. quinquefasciatus data, smoothness priors were assigned whose joint distribution, s, was given by: ( ) Again, assuming diffuses priors, the variance that controlled the degree of smoothness in the model was 1 s  , 11 s , and 2 f τ .The unstructured spatial heterogeneity term, u i, was assumed to follow an exchangeable Gaussian prior with zero mean and variance, , .., ~0, similar prior was then assigned to the heterogeneity term for the sampled predictor variables, i.e., ( ) Finally, for the spatial components, v i , a Markov random field (MRF) prior was assigned.The conditional distribution of v i , given an adjacent covariate, v j , was a univariate normal distribution with mean equal the average v j coefficient values of v i 's neighboring sampled site and variance equal to 2 v τ divided by the number of total adjacent sampled Cx. quinquefasciatus predictors.This lead to a joint density of the form: where i ~ j denoted a sampled site i adjacent to j, and where the parameter estimates v i and v j in the adjacent sampled sites were similar.The degree of uncertainity in the spatiotemporal-sampled data was determined by the unknown precision parameter By writing for a well defined design matrix Z, a vector of regression parameters β, with all different priors, was expressed in a general Gaussian form using the expression:  with an appropriate penalty matrix K j .The model structure was dependent on the sampled explanatory variables and smoothness of the function.In most cases, K j is rank deficient and, hence, the prior for β j is improper [18].For the variances 2 j τ inverse Gamma priors IG (a j , b j ) was assumed, with hyperparameters a j , b j chosen such that this prior was weakly informative.
The Bayesian framework in this research was defined by conditional probabilities constructed from environmental-sampled Cx. quinquefasciatus habitat data.The observation nodes in the Bayesian estimation model were denoted by a vector ( ) , and the set of states of the observation node x j generated from the sampled data that was represented by { } . The conditional probability in this research was that the jth observation node x j was l, 1 j i Y ≤ ≤ , which was based on the condition that the states of hidden nodes were ( ) , generated by ( ) , and let

{ }
, ω α β = be the set of all the spatiotemporal-sampled Cx. quinquefasciatus habitat parameters in the study site.Then the joint probability that the states of observation nodes were ( ) and the states of hidden nodes were ( ) , , , The marginal probability that the states of observation nodes were x was generated using , and we used the notation for the summation over all states of hidden nodes.We assumed the sampled Cx. quinquefasciatus habitat predictor variables were independently and identically taken from the true distribution ( ) o p x .In Bayesian learning, the prior distribution ϕω on the parameter ω is set [18].As such, the posterior distribution ( ) n p X ω was computed from the spatiotemporal-sampled Cx. quinquefasciatus habitat dataset and the prior by which was generated using the expression ( ) ( ) which was used as a criterion by which the model was selected and the hyperparameters in the prior were optimized.We let [ ] n x E  be the expectation over all the sampled Cx. quinquefasciatus parameters.The Bayesian stochastic complexity had the following asymptotic form using where λ and m were the rational number and the natural number respectively, and the two were determined by the singularities of sampled predictor variables.In regular models, 2λ is equal to the number of parameters and m = 1, while in non-identifiable models, 2λ is not larger than the number of parameters and m ≥ 1 [38].However, Bayesian frameworks using field and remote-sampled vector mosquito data requires integration over the posterior distribution, which typically cannot be performed analytically [4].
In this research, we let { } , n n X Z be the sampled Cx. quinuqefasciatus parameter estimates corresponding to the hidden error variables in the equation . The variational Bayesian framework approximated the Bayesian posterior ( ) p Z X ω of the hidden variables and the Cx.quinuqefasciatus parameters using the variational posterior ( ) where ( ) were posteriors based on the inconspicuous error variables and the sampled data respectively.The variational posterior ( ) , l og , , , , which was then further defined by the equation and the true Bayesian posterior was ( ) ( ) where r C and Q C were the normalization constants.It is important to note that these equations gave only necessary conditions for the functional [ ] F q to be minimized in the Cx.quinuqefasciatus habitat model.The variational posteriors were computed by an iterative algorithm.We defined the variational stochastic complex- We then generated variational posterior for the estimation matrix.We assumed that the prior distribution and covariates estimates were provided by which were Dirichlet distributions with hyperparameters generated using 0 φ > 0 and 0 ξ > 0. The Dirichlet distribution [i.e., Dir(α)] is a family of continuous multivariate probability distributions parameterized by the vector α of positive reals which can generate the multivariate generalization of the beta distribution and conjugate prior of the categorical distribution and multinomial distribution in Bayesian statistics for predictive vector mosquito habitat models [4].The Dirichlet distribution is the multinomial extension to the beta distribution for a binomial process which can also be used in quantifying probabilities in predictive vector mosquito habitat probability models [39].
We then let ( ) n δ be 1 when n = 0 and 1 otherwise, and then defined the sampled parameter uncertainty estimates using ( ) : In this research we also used was the state of the jth observation node and ( ) was the state of the kth hidden node.The variational posterior distribution of parameters { } , a b ω = was given by using the equation where each of the spatiotemporal-sampled Cx.Quinqefasciatus habitat covariates were generated using the equation denoted the sum over ( ) In this research we evaluated the statistical efficiency of the MCMC sequence in WinBUGS ® using the following steps: Step1: Sample S (m) from f(S|X,Z (m-1) ) Step2: Sample P (m) from f(P|X,S (m) ,Z (m-1) ) Step3: Sample Z (m) from f(Z|X,P (m) ,S (m) ,Φ (m-1) ) where X was the field and remote-sampled Cx.Quinquefasciatus habitat data and Z was randomly assigned a starting value Z (1) using a uniform prior distribution.
Step 1 was performed by drawing from an inverse Wishart distribution.In Bayesian statistics inverse Wishart distribution is used as the conjugate prior for the covariance matrix of a multivariate normal distribution [9].The probability density function of the inverse Wishart was: A probability density function or density of a continuous random variable is a function that describes the relative likelihood for this random variable to occur at a given point [25].
The marginal and conditional distributions from the inverse Wishart-distributed matrix was quantified using . We then partitioned the matrices for determining if ψ was conformable with each other using: where A ij and ij ψ were Pi x Pj matrices.We then determined if: 1) 11 A was independent of A ⋅ , when  ,A ) , where ( , ) pxq MN ⋅ ⋅ was a matrix normal distribution generated from the spatiotemporal-sampled Cx. quinquefasciatus habitat parameters; 4) ( ) We used the Conjugate distribution to make inference about a covariance matrix  whose prior ( ) are independent p-variate Gaussian variables drawn from a ( ) where A =XX T is n times the sample covariance matrix [3].Because the prior and posterior distributions are the same family, the inverse Wishart distribution was the conjugate to the multivariate Gaussian generated from the georeferenced Cx. quinquefasciatus habitat parameters.
Due to its conjugacy to the multivariate Gaussian it was possible to "integrate out" the Gaussian's parameter  using: This task was simplified by assuming that the spatiotemporal-sampled Cx. quinquefasciatus habitat data analyses had covariance matrices with common eigenvectors.If covariance's differ among sources, the inverse Wishart draws often produce invalid, especially for sources that are small components of the mixture [9].In this research we developed a different approach, noting that a covariance matrix S can be decomposed into a vector of standard deviations, V, and a correlation matrix, R using: S = diag(V)Rdiag(V) where diag(V) was a matrix with diagonal elements V and zeros elsewhere.This decomposition permitted the independent sampling of V and R. We simulated the standard deviations, V, from an inverse gamma distribution: , where n k was the number of individual habitats assigned to cluster k by Z (m-1) (i.e., an estimate of the Cx.quinquefasciatus habitat sample size), s 2 k,l was the sample variance of element l in cluster k as assigned by Z (m-1) , and α and β were constants, both set to the non-informative prior value of 1. Quasi-likelihood techniques in a logistic regression equation and Bayesian prior distributions can enumerate intra-cluster correlations in urban vector mosquito aquatic habitats [4].
We simulated the elements of the correlation matrices, R, from a hyperbolic-tangent transformed distribution using : was the current estimate of the correlation between the i th and j th elements (i ≠ j) in a cluster k given Z (m-1) , and n k was the same.After sampling values for both V and R, we reassembled the covariance matrix, S k , for each Cx.quinquefasciatus habitat cluster, thus, completing Step 1.
Step 2 required drawing values for the elemental means, P. The field and remote-sampled data X had an approximate multivariate normal distribution as such Step 2 was performed using the vector of mean concentrations for cluster k.The multivariate normal distribution using the Cx.quinquefasciatus habitat parameters was given by the sample means calculated from X, generated from the cluster assignments, Z (m-1) ) and the covariance S k from Step 1.If cluster k was empty at m-1 (that was, no individual sampled habitat parameters were assigned to k by Z (m-1) ), then the grand mean and covariance matrix of X were used as the sample mean and covariance matrix of k.
Step 3 involved drawing new cluster assignments using each individual sampled Cx. quinquefasciatus habitat in the Trinidad study site.To do so, it was necessary to calculate Pr(z i = k) for every i, k combination (z i was the i th element of Z).Again we assumed multivariate normal distributions where Z (m) was simulated from:  and where the ( ) f ⋅ terms on the right-hand side were multivariate normal likelihoods.Thus, the likelihood that the i th element of X was present in the Cx.Quinquefasciatus habitat population k, sampled in the study site normalized by the sum of likelihoods was quantified for all potential sources influencing presence of immatures.Finally, the sample Φ (m) from f(Z|X, P (m) , S (m) , Z (m) ) was used to generate asymptotically efficient estimates from the spatiotemporal-sampled data.

Eigenvector Mapping
We conducted a residual analyses using a misspecification perspective for the error estimation models generated in SAS/GIS ® assuming that the basic regression model the misspecification term was Eγ .Quantification of topographic patterns generated from the distribution of the georeferenced Cx. quinquefasciatus habitat covariates was required to describe independent key dimensions and underlying spatial processes for defining a pattern in the misspecification term.This was accomplished by expanding the matrix term: ( ) , as an infinite power series, which was feasible under the assumption that the underlying spatial process in the data was stationary; note that in this research 0 0 V ρ gave the identity matrix I .The simultaneous autoregressive error model was then rewritten as y Vy X VX ρ β ρ β ε − = − + .Substituting these transformation gave: The misspecification term remained uncorrelated with the exogenous variable, X , as the standard OLS assumption of the disturbances, ε , were uncorrelated with the predictor variables generated from the model.The spatial lag model was expressed as: ( ) Substituting the transformation gave: The misspecification term The correlation, or lack thereof, between the exogenous variables and the misspecification terms of both habitat models were used to design spatial proxy variables, so that the properties of either model could be satisfied.We considered two different projection matrices.Different sets of eigenvectors were established using a spatial regression model.This model feature enabled us to also use the eigenvector spatial filtering approach for predictions of the endogenous variable y .
We selected an evaluation criterion using the standardized Moran's statistic.The standardized Moran's coefficient (MC) for regression residuals * ε was defined as

(
)( ) ( ) , where N was the number of sampled habitats and their respective parameter estimates indexed by i and j respectively; X was the spatiotemporal-sampled egg-raft count data; X was the mean of X; and w ij was a matrix of spatial weights.The expected value of Moran's I under hypothesis of no spatial autocorrelation was: ( ) . Its variance equaled: ) ( ) In this research, negative (positive) values indicated negative (positive) spatial autocorrelation.Values ranged from -1 (indicating perfect dispersion) to +1 (perfect correlation).A zero values indicated a random spatial pattern.For statistical hypothesis testing predictive WNV mosquito habitat modeling, MC values can be transformed to z-scores in which values greater than 1.96 or smaller than −1.96 indicate spatial autocorrelation that is significant at the 5% level [8,29].
Geary's Contiguity Ratio (GC) another measure of spatial autocorrelation was also generated.GC is inversely related to the Moran's statistic but it is not identical.Moran's indices are measure of global spatial autocorrelation, while Geary's C is more sensitive to local spatial autocorrelation in vector mosquito data analyses [16].In this research, the sampled Cx. quinquefasciatus habitat data was quantified by GC which was defined as

  
where N was the number of the sampled habitats and explanatory parameters indexed by i and j respectively; X was the egg-raft count data; X was the mean of X; w ij was a matrix of spatial weights; and W was the sum of all w ij .The value of GC lies between 0 and 2. 1 means no spatial autocorrelation while 0 refers to a chaotic random distribution.Smaller (larger) than 1 means positive (negative) spatial autocorrelation [8].
In this research, distance between sampled habitats was defined in terms of an n-by-n geographic weights matrix, C, whose c ij values were; 1 if the sampled Cx. quinquefasciatus egg-raft count locations i and j in the Trinidad study site were deemed nearby, and 0 otherwise.Adjusting this matrix by dividing each row entry by its row sum gave C1, where 1 was an n-by-1 vector of one's which converted this matrix to matrix W. The resulting autoregressive model specification with no sampled covariates present (i.e., the pure spatial autoregression specification) took on the following form: , where μ was the scalar conditional mean of Y and ε was an n-by-1 error vec- tor residuals of the sampled data which in this research was represented as statistically independent and identically distributed (iid) normally random variates.The spatial covariance matrix for equation (2.1), using the sampled covariates, was where E ( ⋅ ) indicated the calculus of expectations, I was the n-by-n identity matrix denoting the matrix transpose operation, and 2 σ was the error variance.Two different autoregressive parameters appeared in the spatial covariance matrix.When autocorrelation is present in residual data, a more explicit representation of both effects leads to a more accurate interpretation of empirical results [38].Consequently, we used an SAR model specification ( )( ) where the diagonal matrix of autoregressive parameters, diag ρ contained two sampled habitat parameters: ρ + displaying positive spatial dependency, and - ρ for those sampled habitat pairs displaying negative spatial dependency based on the aggregation of the Cx.Quinquefasciatus egg-raft count data.For example, by letting 2 σ = 1 and employing a 2-by-2 regular square tessellation, for the vector enabled positing a positive relationship between the sampled egg-raft count data and the explanatory covariates, y 1 and y 2 , a negative relationship between covariates, y 3 and y 4 , and, no relationship between covariates y 1 and y 3 and between y 2 and y 4 .This covariance yielded: ( ) where I + was a binary 0-1 indicator variable which denoted those parameter estimates which displayed positive spatial dependency, and I -was a binary 0-1 indicator variable denoting those sampled covariates displaying negative spatial dependency, where I + + I -= 1. Expressing the preceding 2-by-2 example in terms of an equation yielded: 1 2 3 4 y 1 0 0 0 1 0 0 0 y 0 1 0 0 0 1 0 0 μ ρ y 0 0 1 0 0 0 0 0 y 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 1 0 0 ρ ρ 0 0 1 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 ρ 0 0 1 0 0 0 0 1 if either ρ 0 + = then I + was = 0 and I -= I or if ρ 0 − = then I -= 0 and I + = I.
In this research, this spatial structuring was achieved by constructing a linear combination of a subset of the eigenvectors of a modified geographic weights matrix, using (I -11 ' / n )C(I -11 ' / n ) that appeared in the numerator of the MC.Spatial filtering furnishes an approach that incorporates spatial structuring into an intercept term that accounts for latent autocorrelation [38].A subset of eigenvectors was selected with a stepwise regression procedure.Because (I -11 ' / n )C(I -11 ' / n ) = ' EΛE , where E was an n-by-n matrix of eigenvectors and Λ was an n-by-n diagonal matrix of the corresponding eigenvalues, the resulting egg-raft count model specification was given by: , where μ was the scalar mean of Y, Ek was an n-by-k matrix containing the subset of k < < n eigenvectors selected with a stepwise regression technique and β was a k-by-1 vector of regression coefficients.
A number of the eigenvectors were extracted from the matrix (I -11 ' / n )C(I -11 ' / n ) which were affiliated with geographic patterns of the sampled Cx. quinquefas- ciatus habitat covariates portraying a negligible degree of spatial autocorrelation.Consequently, only k of the n eigenvectors were of interest for generating a candidate set for a stepwise regression procedure.One way of constructing a candidate set is to include an eigenvector in the set only when the vector's corresponding MC represents a minimum degree of spatial autocorrelation [38].In this research, we used a spatial statistical threshold of |MC/MC max | > 0.25 where MC max was the maximum MC for a given surface partitioning, as then each candidate eigenvector represented a level of autocorrelation tending to account for the redundant information in the dataset.As redundant information spatial autocorrelation represents pseudo-replicated data linking it to missing values estimation and interpolation, as well as to notions of effective sample size and degrees of freedom [38].

Spatial Analyses
In this research, spatial linear prediction of the prolific Cx. quinquefasciatus habitat and its spatiotemporal sampled covariates was performed using an Ordinary kriged-based interpolator in ArcGIS 9.3 ® using Geostatistical Analyst.The Semivariogram Properties dialog box has several models to choose from.We set the Kriging method to Ordinary.A default value for lag size was initially set to the default output cell size.For Major range, Partial sill, and Nugget, a default value was calculated internally.The length of the longer axis to reach the sill is called the major range, and the length of the shorter axis to reach the sill is called the minor range; the angle of rotation of the line formed the major range [28].The optional output variance of prediction raster contained the kriging variance at each output raster cell.Assuming the kriging errors are normally distributed in a predictive vector mosquito habitat distribution model there is a 95.5 percent probability that the actual z-value at the cell is the predicted raster value, plus or minus two times the square root of the value in the prediction raster [16].Low values within the output variance of prediction raster indicated a high degree of confidence in the predicted estimates.High values in interpolation based algorithms may indicate a need for more data points (26).
Geostatistical techniques were used to interpolate the value Z(x 0 ) of a prolific Cx. quinquefacistus habitat site, Z(x) at an unobserved sampling site, x 0 , where z i was Z(x i ), i = 1  n based on the dimensions of the sampled habitat locations [i.e., (x 1 ,  , x n )].The Ordinary interpolator also computed the best linear unbiased estimator, Ž(x o ) of Z(x 0 ), for the Cx.quinquefacistus data based on a stochastic model quantified by the variogram γ(x, y), by expectation μ(x) = E[Z(x)], and by the covariance function c(x,y) of the random field.In this research, the kriging estimator was given by a linear combination of the sampled field and remote parameters using ; whereby, z i = Z(x i ) was the weights w i (x o ), and i = 1  n was the variance used to minimize any biased condition.The dependent variable was the spatiotemporal-sampled Cx. quinquefasciatus egg-raft count data, which were transformed to fulfill the diagnostic normality test prior to performing the kriging.The kriging weights of the algorithm were used to fulfill the unbiasedness condition in the spatial interpolation of the ecological-dependent variables using: which was given by the equation system: The Ordinary kriging was given by: , The semivariogram generated described the dependence, between the Cx.quinquefasciatus parameters as a function of the distance between the sample sites.This allowed for abundance estimations at any point in the Trinidad study site.The value of prevalence, Z, at the coordinate (x 0 , y 0 ) was estimated from the n nearest habitat sampling values using Zobs(x 1 , y 1 ), Zobs(x 2 , y 2 ).
Zobs(x n , y n ) and the linear formula which generated: The a i was found by the Lagrange multiplier λ and solving the system using: ( ) ( ) , 1  under the constraint. 1


, where h i,j denoted the distance between any two sampled Cx. quinquefasciatus habitat locations in the Trinidad study site at (x i , y i ) and (x j , y j ), and h j,0 was the distance between the two habitats (x 0 , y 0 ).In the field of calculus of variations in mathematics, the method of Lagrange multipliers can be used to solve certain infinite-dimensional constrained optimization problems [26].In this research the semivariance was defined as γ (h).The magnitude of the semivariance was dependent on the distance between sampling sites.Semivariance of the deviance residuals of the Cx.quinquefasciatus data was calculated and a variogram was constructed to determine if there was evidence of residual spatial correlation in the data.The plot of the semivariances as a function of distance from a point is referred to as a semivariogram [26].
In this research, parameters of a fitted mathematical function (i.e., the semiovariogram model) included generating a range, a nugget and a sill.Since our data was spatiotemporal, the square of the difference between expected values at habitat points had to be added: 2γ(x,y) = C(x,x) + C(y,y) -2C(x,y) -(E(Z(x)) -E(Z(y)) 2 .If the covariance function of a stationary process exists, it is related to variogram by 2γ(x,y) = C(x,x) + C(y,y) -2C(x,y) [26].We checked for spatial dependence also in the sampled Cx. quinquefasciatus habitat data.We then determined the best fit variogram equation based on the number of lags (i.e., bins) and the lag distance (the minimum value of h) used in the calculation of experimental variogram equation.Lag size is the width (distance) of the bins into which these vectors are grouped (www.esri.com).In this research the lag size was dependent upon the minimum and maximum distances between sampled Cx. quinquefasciatus habitats in the study site.In an effort to uncover our variogram's structure, similar lags were grouped together (i.e., pairs of habitat that points aligned in roughly the same direction and roughly the same distance from each other) into the bins.
The semivariogram nugget coefficient allowed an interpretation of the required scale for defining spatial variability of the Cx.quinquefasciatus habitat data for providing minimized unbiased prediction error for residuals on the logit scale.Standard error in the models was calculated.Variations of kriging can produce noise--free predictions (www.esri.com).The kriging interpolator produced a map that was smooth and free of "jumps" at the sampled Cx. quinquefasciatus habitat locations.The algorithms incorporated in Geostatistical Analyst provided exact filtered value predictions at the sampled habitat locations.This prevented discontinuities in the habitat predictions and their associated standard errors.A simple quantitative measure of the interpolation was determined by using a root mean square error (RMSE) for the models.
The kriging procedure returned the observed values at the sampled Cx. quinquefasciatus habitat locations and their interpolated values using the intensity and shape of the variograms.Using a neighborhood and/or distance search radius provided the standard errors of the interpolated values.Prediction errors have often been used to optimize sampling design by identifying areas where sampling effort should be increased or decreased [26].In this research, prediction errors were generated as a function of the variogram models.A simple quantitative measure of the uncertainty in the interpolation model was determined by generating the RMSE values for the models.
Because we were working in two-dimensional space, the semivariogram and covariance functions changed not only with distance but with direction (i.e., anisotropy).We quantified sampled habitats points and the vector that separated them.This vector had a distance on the x-coordinate as well as the y-coordinate.Alternatively, the vector had a distance and an angle in polar coordinates.Anisotropy was then described for the semivariogram.The isotropic model was the same in all directions; whereas, the anisotropic model reached the sill more rapidly in some directions than others.
In Geostatistical Analyst, an outline of the range was given over the empirical semivariogram surface.Having created a covariance function with which to form the kriging weights matrix, the next task was to tessellate the region into patches.Voronoi maps were used to explore the distribution of the sampled Cx. quinquefasciatus habitat data for global and local outliers and to quantify the covariation among the spatiotemporal-sampled data.Voronoi diagram is a special kind of decomposition of a metric space determined by distances to a specified discrete set of objects in the space [28].Voronoi polygons were also used for interpolating neighborhoods and patches.
A Voronoi diagram was constructed with the sample points as the centers of the polygons using the Weighted Voronoi Diagram Extension in ArcGIS 9.3 ® .Using the Generate tab, we generated a weighted Voronoi diagram from the geo-sampled point features.The Graphical User Interface (GUI) has two tabs: Generate and Update (www.esri.com).The Voronoi Diagram was generated whereby a sampled Cx. quinquefasciatus habitat was associated with p i (spatiotemporal-sampled covariate); whereby, P = {p 1 ,  , p n }, where 2 ≤ n ≤ ∞ and x i ≠ x j for i ≠ j, i, j ∈ I n .The region was given by V(p i ) = {x: || xx i || ≤ || xx j || for j ≠ i, i ∈ I n } was the ordinary Voronoi polygon associated with p i or the Voronoi polygon of p i and the set given by V = {V(p i ),.., V(p n )} was the planar ordinary Voronoi diagram generated by P or simply Voronoi diagram of P. We also defined a planar ordinary Voronoi diagram with half planes as follows where we let P = {p i ,  , p n } ⊆ R 2 , where 2 ≤ n ≤ ∞ and x i ≠ x j for i ≠ j, i, j ∈ I n .We called the region which was the ordinary Voronoi polygon associated with p i and set V(P) = {V(p 1 ), , V(p n )} the planar ordinary Voronoi diagram generated by P.
A raster image showing normal Euclidean distance and adjusted Euclidean distance was generated, as well as a Voronoi polygon shapefile.The Cx. quinquefasciatus habitat point attributes were then transferred to Voronoi polygons automatically by appending the spatial attributes of one layer to another.Once a distance raster is generated, the user can use new points to update the distance raster and create updated Voronoi polygons [28].
The Voronoi tessellation produced a set of polygons V i with area f i (i=1,..., n).The (unknown) global mean z D was estimated by a weighted mean of the spatiotemporal-sampled Cx. quinquefasciatus habitat values using: . The extension error of each polygon was calculated by using a discrete version of . If the elementary error terms are uncorrelated, errors can be combined to an estimate of the global error: where the summation is over all polygons V i [28].By convention we assumed a normal probability distribution for the global error and achieved a 90% confidence limits for D Z , 1.65* 1.65* . The global estimation error decreased with an increasing number of samples although some local deviation from this tendency did occur due to large polygons.As the number of samples increased the global estimation error converged to zero.

Results
The ANOVA used in data analyses tested the null hypothesis that the sampled parameters of the immature Cx. quinquefasciatus population in the Trinidad study site means were equal (i.e., H 0 : µ 1 = µ 2 =  = µ a) , by comparing two estimates of variance.If the null hypothesis is false, then Mean Square Between (MSB) estimates generated from a predictive vector mosquito habitat distribution model is something larger than σ² [25].In this research the MSE was an estimate of variance for determining whether or not the null hypothesis was true, while the second estimate MSB was based on the variance of the sample means.We tested the null hypothesis as follows: if the null hypothesis was true, then MSE and MSB would be about the same since they are both estimates of the same quantity (i.e., σ 2 ); however, if the null hypothesis was false then MSB would have been expected to be larger than MSE, since MSB is estimating a quantity larger than σ².The ANOVA model encompassed all possible sources of variation in sampled Cx. quinquefasciatus habitat data which allowed further testing of our research hypotheses.In this research the significance test of spatiotemporal-sampled data involved the ratio of MSB to MSE: F = MSB/MSE.The F-statistic was used to calculate the p-value to determine deviations from normality.If the kurtosis is greater than 0, then the F tends to be too small and we cannot reject the null hypothesis even though it is incorrect; the opposite is the case when the kurtosis is less than 3 [25].The skewness of the distribution in the Cx.quinquefasciatus data did not have a sizable effect on the F statistic.If the n per cell is fairly large, then deviations from normality do not matter much at all because of the central limit theorem, which implicitly states the sampling distribution of the mean approximates the normal distribution, regardless of the distribution of the variable in the population [25].The P value was the estimated probability of rejecting the null hypothesis (H0).Traditionally for WNV mosquito models, if p < .05 the null hypothesis is rejected [18].If the null hypothesis is true in the estimation of a model then the F ratio should be approximately one, since MSB and MSE should be about the same [25].In order to conduct a significance test using the prolific spatiotemporal-sampled Cx. quinquefasciatus habitats predictor variables, it was thus necessary to know the sampling distribution of F.
Additionally, from the sampling distribution, generated from the sampled t parameters, the probability of obtaining an F as large or larger than the one was also calculated.When there are only two means to compare, the t-test and the F-test are equivalent; the relation between ANOVA and t is given by F = t 2 [25].If the probability values were lower than the significance level, then the null hypothesis was rejected.Significant differences by ANOVA were noted for mean numbers of Cx.Quinquefasciatus captured throughout the sampling frame (F = 42.2,degrees of freedom [df] = 1).
The Poisson regression modeled the spatiotemporalsampled Cx. quinquefasciatus habitat count data and contingency tables.Poisson regression assumed the response variable Y had a Poisson distribution and assumed the logarithm of its expected value can be modeled by a linear combination of the spatiotemporalsampled Cx. quinquesciatus parameters.A Poisson regression model is sometimes known as a log-linear model, especially when used to model contingency tables [25].Our model with a single independent variable x, took on the form: ( ) ( ) If Y i are independent observations with corresponding values x i of the sampled predictor variable, then a and b can be estimated by maximum likelihood if the number of distinct x values is at least 2 [25].In this research the maximum-likelihood estimates lacked a closed-form expression and had to be quantified using numerical methods.The probability surface for maximum-likelihood Poisson regression was convex, making Newton-Raphson or other gradient-based methods appropriate estimation techniques for data analyses of the spatiotemporal-sampled Cx. quinquefas- ciatus habitat data.Linearized iteration scheme converges quickly and accurately using inversion equations solved through the Newton-Raphson method as each iteration, requires solving the finite difference equations and the linear adjoint equations only once, respectively [25].
OrthoEngine, PCI Geomatics' Geocorrection and Extraction tools generated a DEM of the Trinidad study site.OrthoEngine included a number of automation tools to improve the efficiency of georeferencing the sampled Cx. quinquefasciatus data including Automatic GCP collection from Chip Databases, and Automatic TP measurement.We stitched the image files.Stitching operation merged all the raster tiles and recalculated the image RPC from individual RPCs.Stitching automatically revealed the final RPCs as a binary segment in an output PIX file.Coefficients were reviewed using 'Project Report Option' in 'Reports' procedure step.We then selected geometric Model under the image information panel.The final step was the 'Schedule Ortho Generation.'We proceeded to the 'Orthogeneration' processing step and selected the files to be processed.We selected an appropriate DEM file and set other processing options before generating the final orthorectified image.This allowed for the precise orthorectification of the WV-I imagery.The WV-1 over the Trinidad study site had look angles of approximately south at 30 degrees-off nadir and vertical.For simple corrections, polynomial, thin plate spline and rational functions were used.This allowed using image correlation from the stereo images.Once the images were geometrically corrected, they were mosaicked into a seamless image database, to form a conventional image map product and a GIS base layer for derivation of vector hypsographic and hydrographic covariates associated to the sampled Cx. quinquefascia- tus habitat data.The range of the elevation in the DEM had a minimum value of 0 m with a maximum value of 425 m.
The spatial autocorrelation analysis results are re-ported in Table 4. Eigenvectors were extracted from the matrix (I−11′/n) C (I−11′/n) using the ecological-sampled predictor variables.Estimation results for these models appear in Table 5. Positive and negative spatial autocorrelation spatial filter component pseudo-R 2 values are reported.Positive and negative spatial autocorrelation spatial filter component pseudo-R 2 values are reported.These values did not exactly sum for the complete spatial filter; however, they are very close to their corresponding totals, suggesting that any induced multicollinearity was quite small.To summarize Table 5 results, the Trinidad study site had a random effects intercept with a near-zero mean and only trace spatial autocorrelation.The random effects term accounted for two-thirds or more of the variation in the spatiotemporal-sampled Cx. quinquefasciatus habitat count data.In addition, the random effects term resulted in the loss of a counterbalancing spatial autocorrelation component in the spatial filters.The spatial autocorrelation components revealed 17% redundant information in the Cx.quinquefasciatus habitat datasets.
Table 6 lists the improvements of fit for all model specifications and random error in the spatial analyses.The unadjusted model compared the univariate model to a model containing only the intercept term.Improvement of fit was also calculated for the first-order interaction models to determine whether including significant interactions improved fit compared to the full main effects model.
Table 7 presents the results of the regression for the interactions model.These results provided information for estimates of the prior distribution of main effect coefficients in the Bayesian analysis.The values for parameter estimates and standard were used as mean values and standard errors to parameterized prior expected values for quantifying the field and remote-sampled Cx. quin- quefasciatus habitat explanatory variables.The prior expected mean value for the error term was assumed to be zero ('0'), with a standard deviation of 0.01.Initial values for the MCMC chains were generated.
The difference in the deviances between a simple model and the more complex model provided the improvement χ 2 values listed in Table 8.We examined all interaction between the sampled predictor variables and found that an interaction model did not improve the fit; therefore, no interaction terms were included in the final model.We could not examine the improvement of fit between a saturated model and the full effects model, as the number of the sampled parameters that needed to be estimated exceeded the maximum number that could estimate.
To derive the improvement of fit values listed in Table 9, the posterior mean deviance values were obtained with Deviance Information Criterion (DIC) spatial analytical tools.Factors that did not improve fit were omitted from the final model.DIC tools were used to obtain mean posterior deviance values to construct improvement of fit tables and DIC statistics for identification of the best fitting model.In this research, this deviance was defined   Median parameter values, as well as the 95% credibility intervals (2.5 percentile and 97.5 percentile values), are listed in Table 10.As the sampling sites increased based on the covariate distance from the nearest house, the median log-egg-raft count of count data changed.The adjusted model assumed independence among the sampled predictor variables of Cx. quinquefasciatus aquatic habitat egg-raft count data fit better that the model that adjusted for correlation within the study site based on the RMSE (Figure 3).
The kriged-based model resulted in two images, a surface of estimates and a surface of estimated variances.
The latter image was used to identify problems with the fit of the model to the sampled data by revealing relative differences in the model fit across the Trinidad study site.A measure of the degree of spatial dependence between samples was also generated (i.e., semivariance).The magnitude of the semivariance between any two sampled Cx. quinquefasciatus habitats was dependent on the distance between the georeferenced habitats.Thus, given two habitat locations x and (x + h), in the study site, a measure of one half of the mean square difference (i.e.,  A kriged map of deviance residuals was then calculated, which was added to the predicted values on the logit scale before transforming the result back to proportions.Spatial dependence displayed by these plots was modeled using the constructed semivariogram.The maps developed from the kriged uncertainty residuals allowed deviation from the model and incorporated the field sampled egg raft count values.The final maps were greatly improved in sensitivity and did not deviate too severely from the field-sampled Cx. quinquefasciatus habitat data.The exponential model fitted to the semivariogram at a partial sill 0 and nugget of 5.8 km, lag size 6.5 km number of lags 12 with a range of 73.6 km which formed the basis of the model to map predictions of residuals in an 73.6 km radius around each sampled Cx. quinquefasciatus habitat (Figure 4).
Voronoi polygons were then modeled.An approximation, based on the midpoints of the legs and center of the triangle were used to evaluate the area covered by each Voronoi polygon.A new sample point was calculated using the importance metric (Gamma) of each pixel in the Voronoi polygons generated according to ( ) where (u, v) were the coordinates of a WV-I pixel in the polygon, Gamma(u, v) was the importance metric of that pixel, and C(p) was the coordinates of the new polygon center (i.e. the sampled Cx. quinuqefasciatus habitat).To do this calculation on each polygon we used a scan conversion.We ran a scan conversion and examined each pixel in the frame buffer for each polygon.The Voronoi tessellation provided a spatial trend analyses of the error in the model which revealed that all coefficients were within normal statisticcal limitations (Figure 5).

Discussion
The ANOVA tested the significant differences between means in the Cx.quinquefasciatus habitat data.In repeated measures ANOVA; however, containing factors such as spatiotemporal-sampled Cx. quinquefasciatus with more than two levels, additional special assumptions may be required for quantification of sampled predictor variables; specifically, the compound symmetry assumption and the assumption of sphericity.The compound symmetry assumption requires that the variances derived from the data (i.e., pooled within-group) and  covariances of the different repeated measures be homogeneous (25).This is a sufficient condition for the univariate F test for repeated measures to be valid (i.e., for the reported F values to actually follow the F distribution); however, it is not a necessary condition in spatio-temporal Cx. quinquefasciatus habitat data analyses.When the compound symmetry or sphericity assumptions have been violated, the univariate ANOVA table will give erroneous results (23).A remedy may be to stress the necessity of independent hypotheses testing in spatiotemporal Cx. quinquefasciatus data analyses.A general algorithm implemented can attempt to generate, for each effect, a set of independent (orthogonal) contrasts.In repeated measures ANOVA, these contrasts specify a set of hypotheses about differences between the levels of the repeated measures factor.However, if these differences are correlated across subjects, then the resulting contrasts are no longer independent.In fact, in most instances where a repeated measures ANOVA is used, we would probably suspect that the changes across levels are correlated across the sampled Cx. quinquefas- ciatus habitat parameters.When this happens, the compound symmetry and sphericity assumptions have been violated, and independent contrasts cannot be computed.
Thus, the problem of compound symmetry and sphericity in Cx. quinquefasciatus habitat data analyses pertains to the fact that multiple contrasts involved in testinsg repeated measures effects with more than two levels are not independent of each other.However, they do not need to be independent of each other if we use multivariate criteria to simultaneously test the statistical significance of the two or more repeated measures contrasts.This insight is the reason why Multivariate Analysis of Variance (MANOVA) methods are increasingly applied to test the significance of univariate repeated measures factors with more than two levels.Multivariate analysis of variance is a generalized form of univariate ANOVA which is used in cases where there are two or more dependent variables (23).Analogous to ANOVA, MAN-OVA is based on the product of model variance matrix and error variance matrix inverse; thus, invariance considerations generated from spatiotemporal-sampled Cx. quiquefasciatus habitat predictors will imply the statistic is a suitable measure of magnitude of a singular value decomposition from this matrix product.Thus, instead of a univariate F value, we would obtain a multivariate F value (Wilks' lambda) based on a comparison of the error variance/covariance matrix and the effect variance/covariance matrix.The problem of compound symmetry and sphericity pertains to the fact that multiple contrasts involved in testing repeated measures of Cx. quinquefasciatus habitat effects with more than two levels are not independent of each other but, fortunately a multivariate criteria can simultaneously test the statistical significance of the two or more repeated measures contrasts.We thus endorse a MANOVA approach for spatiotemporal data analyses of Cx. quinquefasciatus habitat parameters as it simply bypasses the assumption of compound symmetry and sphericity altogether.
The spatial hydrological model generated from the environmental-sampled Cx. quinquefasciatus covariates revealed that elevation was an important variable in the DEM model.Anderson et al. (2004) captured signifycantly more Culex quinquefasciatus in a variety of traps (bird-baited, MMX mosquito magnet, CO 2 -baited CDC trap) in the canopy (7.6 m) than at near ground level (1.5 m) [40].In our research local differences among WV-1 DEM grid cells were also analyzed to calculate slope and other land surface parameters using relative and absolute accuracy estimates for determining the statistical significance of terrain-related parameters associated to the Cx.quinquefasciatus habitats.In this research the slope gradient was defined as the maximum rate of change in altitude (i.e., tan (Θ), aspect (ψ) as the compass direction of this maximum rate of change).Slope is defined by a plane tangent to a topographic surface, as modeled by the DEM at a point, whereby, it is classified as a vector; as such it has a quantity (gradient) and a direction (aspect).More analytically, the slope gradient at a sampled Cx. quinquefasciatus habitat in the Trinidad study site was the first derivative of elevation (Z) with respect of the slope (S), where S was in the aspect direction (ψ).We calculated the first derivative of a function at a sampled habitat as the angular coefficient of the tangent to the function of a sampled habitat.The resolution of our WV-1 DEM quantified the georeferenced explanatory covariates at each 0.5m pixel (i.e., absolute accuracy).At the same time the relative accuracy of the sampled data generated a geometrically correct reference frame for validating the covariate elevation associated with the sampled Cx. quinquefasciatus habitats.WV-1 DEMs can allow for hydrologic modeling, view-shed determination, slope/aspect analyses, and 3-d surface visualization of georeferenced terrain-related parameters associated to spatiotemporal-sampled Cx. quinquefasciatus habitat data.
In the regression model, PROC MIXED quantified correlations among error measurements using random effects estimates and random regression coefficients.PROC MIXED used three options for the method of esimation including: 1) Maximum Likelihood (ML); 2) Restricted or Residual maximum likelihood (REML), which was the default method; and, 3) Minimum Variance Quadratic Unbiased Estimation (MIVQUE).In this research, the ML was the regular maximum likelihood method as the Cx.quinquefasciatus parameter estimates maximized the likelihood function, while the REML was a variant of maximum likelihood.The sample mean was the maximum likelihood estimator of the sampled habitat population mean, based on the egg-raft count data.The sample variance was also a close approximation to the maximum likelihood estimator of the immature sampled variance.REML estimated the variance parameters using a marginal likelihood which produced unbiased estimates.REML estimators generated estimates from the spatiotemporal-sampled parameters by maximizing only that part of the likelihood function that was invariant to the fixed effects part of the linear model.The variance of a quadratic function of the random variables in the linear model was then minimized to obtain locally best unbiased estimators of the variance components.The matrices of these quadratics were computed from relationship matrices and from prior estimates of variances used in the model equations.
In the Bayesian matrix we were able to instantiate the nodes directly which was equivalent to supplying, a measurement for a radius based on the distribution of the georeferenced Cx. quinquefasciatus habitats.We were able to quantify an uncertainty estimate in our likelihood information by supplying a probability distribution over a set of radii generated from the sampled data.We used a backward propagation to compute the node probabilities.The error gradients were then propagated in any direction throughout the Bayesian framework.This type of probability propagation allowed choosing between the sampled predictor variables using the computed probability of each sampled covariate value based on one or more of the nodes.
Our Bayesian approach allowed flexible model fitting and estimation and geomapping of all "high risk" Cx. quinquefasciatus habitats based on sampled egg-raft count data in the study site.The MCMC algorithm for the Cx.quinquefasciatus egg-raft count model produced a sequence of parameter vectors that represented random draws from the posterior distribution.The DIC comprised two goodness-of-fit measures and the posterior distribution of the deviance, which was the number of effective sampled habitat covariates for measuring complexity in our Cx.quinquefasciatus model.Habitats with high egg-raft counts were compared with the results of a Monte Carlo simulation, which established the probabilities and occurrences of these habitats.Our results indicated that likelihood weights influenced the resulting posterior distributions of the field and remote-sampled Cx. quinquefasciatus egg-raft count parameters, which, in turn influenced the spatial trends in the variance uncertainty estimates for model prediction.Distance from the nearest house was a significant predictor associated with prolific habitats in the Trinidad study site.
Environmental factors commonly associated with Cx. quinquefasciatus habitats include the distance from the nearest human habitation.The location of WNV cases at the place of residence can be based on the fact that Culex mosquitoes are crepuscular feeders, and are most likely to come in contact with people while they are at home during the evening hours.Because egg-laying is a spatial process dependent upon the location of the focal habitat relative to sources of gravid mosquitoes, habitats closer to human inhabitations tend to receive more eggs and thus are more productive if conspecific competition is negligible [29].Furthermore, housing age can also play a big role in the explanation of the focal areas.This is expected since older houses have poorer drainage systems and open foundations which may provide larval and adult habitats.For example, houses built in the post World War II era in the study site may suffer from a general inattention.These houses may contain catch basins, rich with organic material which are breeding grounds for Culex mosquitoes.Most of the Culex species require organic-rich water for larval development and are commonly found in high numbers in artificial containers, unattended pools, retention ponds, storm drains, catch basins, sewage systems and treatment plants [41,42].Furthermore, combined sewage and street runoff discharges into natural waterways are considered the main sources of urban stream pollution, which are associated to WNV as Culex mosquitoes proliferate in eutrophic hydrological networks.The identification of social factors that characterize these focal areas in the Trinidad study site may also provide insight into human risk and can help to target control and prevention strategies.Future analyses should explicitly test whether enzootic and bridge transmission are related to human population factors includeing household income, population age, age of housing, housing density, and population density.Since it can be assumed from this analysis that the WNV is transmitted near a person's residence in the Trinidad study site, extensive engineering to improve the runoff of water in residential areas should be considered.
A framework for assessing autocorrelation in residual error from the spatial analyses was constructed in SAS/ GIS ® .We generated a semiparametric filtering model using proxy variables constructed from the spatiotemporal-sampled georeferenced explanatory variables.In our Cx.quinquefasciatus habitat model, the misspecification terms was replaced by the proxy vari-ables for implementing a conditional covariance matrix employing the autoregressive specification.An unknown misspecification term can be approximated by a set of spatial proxy variables [38].After modeling the mis-specification terms, the remaining residuals ε become white noise.This result implied that the estimated re-gression parameters β were unbiased for the basic regression model, * y Xβ ε = + , where * ε incorpo-rated the misspecification term and the white-noise dis-turbances.This allowed us to calibrate the autoregressive models.All of these eigenvectors exhibited weak posi-tively autocorrelated spatial patterns in the habitat pa-rameter estimates.Local weather patterns, mosquito con-trol programs and movement of hosts can have spatial expressions, which cause immature Culex mosquetoes to aggregate in geographic space rendering positive autocorrelation in sampled predictor variables [8,26].
In this research, we used an Ordinary interpolation in ArcGIS ® for identifying initial spatial structures in the sampled Cx. quinquefasciatus habitats and for quantifying the variance and autocorrelation in the data.The estimation method determined the best linear unbiased estimator.It was linear since the estimated values were weighted linear combinations of the spatiotemporalsampled data and it was unbiased because the mean of the error was zero.For determining the optimal predicttors, a semivariogram was constructed which expressed the spatial variation in the sampled data.In this research the variogram [i.e., 2γ(x,y)] was a function describing the degree of spatial dependence between the sampled Cx. quinquefasciatus habitats [i.e., Z(x)].This was defined as the expected squared increment of the values between the sampled habitat locations x and y.Our semivariogram was nonnegative since it was the expectation of a square.The covariance function was related to variogram by 2γ(x,y) = C(x,x) + C(y,y) -2C(x,y).In this research, the γ(x,y) = E(|Z(x) − Z(y)| 2 ) = γ(y,x) was a symmetric function, consequently, γ s (h) = γ s (h) was an even function.A function is a semivariogram if, and only if, it is a conditionally negative definite function, i.e. for all weights 1 , , N w w  subject to 1 0 , thus, it holds: , 0   [28].The function corresponded to the fact that the variance [i.e., var(X) of ( ) generated from the sampled Cx. quinquefasciatus data was given by the negative of this double sum, which in this research was nonnegative.All spatial dependence was quantified in the model.If a stationary random field has no spatial dependence (i.e.C(h) = 0 if 0 h ≠ the semivariogram is the constant var(Z(x)) everywhere except at the origin, where it is zero [28].
Theoretically, at zero separation distance (lag = 0), the semivariogram value is zero in a predictive vector mos-quito habitat distribution model; however, at an infinitesimally small sampled habitat distances, the semivariogram often exhibits a nugget effect (i.e., measurement errors or spatial sources of variation at distances smaller than the sampling interval), which is some value greater than zero [16].Additionally, since the random field was stationary and ergodic, the corresponded to the variance generated by the Cx.quin- quefasciatus habitat data.We defined a practical range and defined the distance at which 95% of the sill was reached for the asymptotic variogram model.The semivariogram calculated estimates of the surface at the different habitat locations.The increase in the semivariogram values with increasing lags diminished with distance and leveled off, at around 70 km (i.e., the range or active lag distance).This was the approximate distance at which spatial autocorrelation between the Cx.Quinque- fasciatus habitat data point pairs ceased or became much more variable.In the future sampling variograms should be constructed using spatiotemporal-sampled Cx.Quin- quefasciatus habitat data.The sampling variogram, unlike the semivariogram, shows where a significant degree of spatial dependence in the sample space or sampling unit dissipates into randomness when the variance terms of a temporally or in-situ ordered set are plotted against the variance of the set and the lower limits of its 99% and 95% confidence ranges [28].Sampling variograms may quantify extraneous measurement variances before spatial dependence is verified in a predictive Cx. quinquefasciatus habitat distribution model.In this research interpolation accuracy was also measured by the natural logarithm of the mean squared interpolation error and Voroni polygons which revealed all uncertainty effects while quantifying several covariate interaction terms in Euclidean space.The Voroni cells generated were a d-dimensional polytope (i.e., a convex hull of a finite set of georeferenced points).A set of hyperplanes were generated using the polygonal shapes around each sampled Cx. quinuqefasciatus habitat in the Trinidad study site.If S contains only two points, a and b, then the set of all points equidistant from a and b is a hyperplane-an affine subspace of codimension 1 [28].The set of the polytopes covered the whole space and was the Voronoi tessellation corresponding to the georeferenced Cx. quinquefasciatus habitats.The convex Voronoi polygons were visualized as space-filling polygons constructed around the dataset of the sampled covariates (i.e.,Voronoi centers), such that each polygon contained all points closer to its Voronoi center than to the center of any other polygon.The relationship of the Voronoi centers to edges of polygons were used to test whether any convex tessellation existed in the polygons.This test amounted to finding Voronoi centers that best fit the given tessellation.In this research, Voronoi centers were found by solving two systems of linear equations.These equations represented 1) conditions on the slope of polygon edges relative to the slope of lines through Voronoi centers, and 2) conditions on the distance from edges to Voronoi centers.Least squares and constrained least-squares solutions were used to solve the two systems.In order to quantitatively characterize the error in the interpolation model, we examined normalized cell size distribution related to the homogenous and inhomogeneous fixation densities (i.e., Cx. quinquefasciatus habitat egg-raft count), which provided additional insight as to how the convex tessellations varied between the Voronoi polygons.A goodness-of-fit statistic was derived using the Voronoi tessellations which determined the error propagation in the multivariate model.The Voroni residuals reflected sampling errors and analytical errors based on the sampled Cx. quinquefasciatus habitat egg-raft count density.
In the future new network nodes using Dirichlet tessellations in spatiotemporal-sampled Cx. quinque- fasciatus habitats should be analyzed.For example, Okabe and Satoh, (2008) have developed a software package, SANET, which can be installed as an ArcGIS add-in for "Spatial Analysis on a Network" for determining new network nodes [43].The SANET software includes a network-based Voronoi tool, very similar to some of the network analysis tools provided in ArcGIS itself and in network-oriented packages such as TransCAD and Cube.The point set (i.e., off-network sites) can be linked to the nearest network edge, and then non-directed network Voronoi regions can be generated from georferenced Cx. quinquefasciatus habitat locations using shortest path distances.Separate inward and outward-directed Voronoi regions can then be computed and each habitat can be derived using unweighted calculations.
For future research, the risk of WNV transmission in the Trinidad study site should also be investigated using wild bird roosting behavior and seasonal migratory pattern data.Quantifying WNV amplification at a local scale, with explicit attention to distances measurements from mosquito and hatch-year bird sites and human habitation areas may affect timing of infection in birds, mosquitoes, and humans.Seropositive resident and migratory birds have been documented along the Atlantic and Mississippi flyways into the Caribbean [2].While transient birds may suggest a means by which WNV has rapidly spread across North and Central America, resident breeding birds may be more critical for zoonotic transmission in Trinidad [44][45][46][47].In North America, high seroprevalence is associated with avian species that breed year-round and live close to humans (e.g.American robin, House sparrow, Mourning dove, Northern cardinal, Northern mockingbird), though other wild birds may be important [47][48][49].In contrast to the urban avian species commonly associated as WNV hosts in North America, avian WNV hosts that occur in Trinidad are often non-urban, wild bird species [2].Caribbean birds that have tested seropositive for WNV and are resident to Trinidad or breed seasonally in Trinidad include: Bananaquits (Coereba flaveola), Palm-tanagers (Phaenicophilus palmarum), and Grossquit species [2].Since the significance of avian species in WNV transmission may relate to their host competency, natural exposure to mosquito vectors and spatiotemporal availability, an avian GIS-based model using spatial filter eigenvectors and Bayesian estimates might be able to generate globally asymptotically stable parameters while quantifying the amplification cycle prior to bridge vector involvement.Bridge vectors are mosquito species that readily feed on humans and mammals [2].
In conclusion, correlation estimates generated from the WV-1 DEM revealed that the terrain covariate elevation was significantly associated with the Cx.quinquefasciatus habitats sampling sites.The spatial hydrological model allowed for hydrologic modeling, view-shed determination, slope/aspect analyses, and 3-d surface visualization of the georeferenced parameters associated to Cx. quinquefasciatus habitat data.Elevation data was represented in the WV-1 DEM by including gridded data where the terrain covariate was estimated for each cell in a regular grid, a triangular irregular network, and contours.The coefficients from the decomposition of the sampled predictor variables were then entered into a Bayesian matrix.This framework allowed flexible model fitting, estimation, and mapping of all high risk habitats based on field-sampled egg-raft count data using a MCMC specification.Residual estimates were then analyzed with an eigenvector spatial filtering algorithm which transformed all the georferenced predictor variables containing spatial dependence into covariates free of spatial dependence by partitioning the original sampled habitat data and the sampled egg-raft abundance count data into two synthetic variates: 1) a spatial filter variate capturing latent spatial dependency and 2) a non-spatial variate that was free of spatial dependence.The geographic distribution of the sampled habitats based on the egg-raft counts exhibited positive autocorrelation in all models tested: log-egg-raft counts aggregated in geo-space.The procedure accounted for 17 % pseudo-replicated information in the datasets.All models were then adjusted based on the residual estimates form the eigenfunction decomposition algorithm and then kriged.The sill indicated that the semivariance values had been reached (i.e., the value of maximum variance was equivalent to the variance of a 0.5m pixel value), while a non-zero intercept value (i.e., nugget variance) of the varigram model was indicative of the variability of the field and remotesampled Cx. quinquefasciatus data.A simple quantitative measure of the interpolation performed was determined by generating RMSE values for the models and by constructing Voroni polygons.Optimizing the RMSE by minimizing the spatial structure in a Cx.quinquefasciatus habitat model generated a pure nugget variogram, of which the level of nugget variance represented noise characteristics in the predictor variables.Interpolation accuracy tests measured by the natural logarithm of the mean squared interpolation error and the Voroni polygons revealed all error effects of the parameter estimates in the kriged-based model and quantified the covariate interaction terms.Newer GIS software, WV-I data and spatial statistics can implement larval control strategies by targeting prolific habitats of Cx. quinquefasciatus based on field-sampled count data in Trinidad.
In the future high quality spatial autoregressive and surface models should be constructed using spatiotemporal-sampled Cx. quinquefasciatus and other mosquito habitat data and Worldview-2 (WV-2) imagery for improving sensitivity in LULC, elevation and kriged-based interpolation parameters in Trinidad.On October 8, 2009 Digitalglobe launched WV-2 which is the first high resolution commercial satellite to offer 8-band capability with high accuracy, agility and spectral diversity (www.digitalglobe.com).This allows the creation of accurate Cx. quinuqefasciatus habitat maps without the need of GCPs.A wide range of features can be extracted using the traditional four spectral bands but WV-2 offers the addition of the Coastal Yellow, Red edge and shorter wave NIR bands which may provide a wider range of potential derivatives for targeting prolific Cx.Quinque- fasciatus and other mosquito habitats.For example, the Coastal-blue bands can be used for bathymetry, benthic, wetland and atmospheric modeling which can assist in detecting subtle differences in land cover and coastal features associated with prolific mosquito habitats.Probably the best feature of WV-2 for geomapping Cx. quinquefasciatus and other mosquito habitats, however, is the ability to derive algorithms and functions from new band combinations for feature extraction and improved classification.For example, the Red-edge NIR band can identify plant health and age through chlorophyll production while the yellow band in conjunction with the Red Edge band can provide agricultural applications and iron oxide mineralogy which can model moisture content.Thus, WV-2 can provide detailed imagery for precise map generation, LULC change detection, construction of DEMs and autoregressive predictive models.WV-2 data can further help implement a remote-based habitat surveillance in Trinidad and in other regions.

Figure 1 .Figure 2 .
Figure 1.Gridded area of urban land cover in the Trinidad study site using WV-1 with georeferenced Cx. quinquefasciatus habitats.
written as SST = SSM + SSE, where SS was notation for sum of squares and T, M, and E were the notation for total, model and error, respectively.The square of the sample correlation was equal to the ratio of the estimates while the sum of squares was related to the total sum of squares: r² = SSM/SST.This formalized the interpretation of r² as explaining the fraction of variability in the sampled egg-raft Cx. quinquefasciatus habitat population in the Trinidad study site explained by the regression-based parameters.The sample variance s y ² was equal to equal to the SST/DF (degrees of freedom), the total sum of squares divided by the total DF.A linear regression equation was constructed using the mean square model (MSM) = also equal to SSE/DF and also the estimate of the variance about the regression line (i.e.,σ²).The MSE is an esti- e., the normalization constant).The Bayesian predictive distribution ( )n p x Xwas provided by averaging the model over the posterior distribution as follows,

.
This led to the functional[ ]F q being minimized under the constraint and then the variation posteriors, ( )


B and ψ were positive definite matrices, and p τ (•) was the multivariate gamma function.In this research the distribution of the inverse of a Wishart-distributed matrix was generated from , which was of size of the matrix constructed from the field and remote-sampled Cx. quinquefasciatus habitat parameters.Under these circumstances 1 p τ (•) was the multivariate gamma function of the sampled data.

*
was decomposed into a white-noise component, ε , and a set of unspecified and/or misspecified models that had the structure * used as a measure of model fitness based on the values of the predictor variables of the habitats sampled.The effective number of parameters included in the Cx.quinquefasciatus habitat model was computed as ( ) was the expectation of θ .The DIC was calculated as: DIC value for the final model was 930.3.The DIC value for the model was 933.4.

Figure 3 .
Figure 3. Observed Cx. quinquefasciatus habitats egg-raft count data in the Trinidad study site using a a) non-spatially adjusted and b) in a spatially adjusted model.thesemivariance) was produced by assigning the values z(x + h) to the value z(x), where h (i.e., the lag) was the inter-sample distance between the georeferenced Culex habitats.Semivariance is an average of the squared deviations of values that are less than the mean[25].In this research, plotting the semivariance versus the lag h gave a semivariogram.The summation was carried out over all distances h.The semivariance increased as the distance increases until it reached a maximum at a certain distant away from a sampled Cx. quinquefasciatus habitat.When the semivariance equals the variance around the averaged values, and no longer increases, this causes a flat region to occur on the semivariogram, (i.e., the sill)[28].The sill indicated that the semivariance values had been reached (i.e., the value of maximum variance was equivalent to the variance of the WV-1 0.5 m pixel value).A non-zero intercept value (i.e., nugget variance) of the variogram model was also generated which indi-

Figure 4 .
Figure 4.An ordinary kriged-based interpolation of egg-raft abundance of Cx. quinquefasciatus habitats in the trinidad study site.

Figure 5 .
Figure 5. Voroni tessellations generated from the predicted spatiotemporal-sampled Cx. quinquefasciatus habitats in the Trinidad study site.

Table 1 . Rational polynomial coefficient results for the WV- 1 digital elevation model.
GCP: Ground Control Points, CP: Control Point, RMS: Root Mean Square