Chagas disease vector control and Taylor's law

Background Large spatial and temporal fluctuations in the population density of living organisms have profound consequences for biodiversity conservation, food production, pest control and disease control, especially vector-borne disease control. Chagas disease vector control based on insecticide spraying could benefit from improved concepts and methods to deal with spatial variations in vector population density. Methodology/Principal findings We show that Taylor's law (TL) of fluctuation scaling describes accurately the mean and variance over space of relative abundance, by habitat, of four insect vectors of Chagas disease (Triatoma infestans, Triatoma guasayana, Triatoma garciabesi and Triatoma sordida) in 33,908 searches of people's dwellings and associated habitats in 79 field surveys in four districts in the Argentine Chaco region, before and after insecticide spraying. As TL predicts, the logarithm of the sample variance of bug relative abundance closely approximates a linear function of the logarithm of the sample mean of abundance in different habitats. Slopes of TL indicate spatial aggregation or variation in habitat suitability. Predictions of new mathematical models of the effect of vector control measures on TL agree overall with field data before and after community-wide spraying of insecticide. Conclusions/Significance A spatial Taylor's law identifies key habitats with high average infestation and spatially highly variable infestation, providing a new instrument for the control and elimination of the vectors of a major human disease.


Methodology/Principal findings
We show that Taylor's law (TL) of fluctuation scaling describes accurately the mean and variance over space of relative abundance, by habitat, of four insect vectors of Chagas disease (Triatoma infestans, Triatoma guasayana, Triatoma garciabesi and Triatoma sordida) in 33,908 searches of people's dwellings and associated habitats in 79 field surveys in four districts in the Argentine Chaco region, before and after insecticide spraying. As TL predicts, the logarithm of the sample variance of bug relative abundance closely approximates a linear function of the logarithm of the sample mean of abundance in different habitats. Slopes of TL indicate spatial aggregation or variation in habitat suitability. Predictions of new mathematical models of the effect of vector control measures on TL agree overall with field data before and after community-wide spraying of insecticide.

Conclusions/Significance
A spatial Taylor's law identifies key habitats with high average infestation and spatially highly variable infestation, providing a new instrument for the control and elimination of the vectors of a major human disease. PLOS

Chagas disease
Vector-borne pathogens contribute to 17% of the global human disease burden [1]. Chagas disease or American trypanosomiasis, one of the World Health Organization's "neglected tropical diseases," is caused by the protozoan Trypanosoma cruzi. It is transmitted mainly by diverse triatomine bug species associated with selected wild, peridomestic and domestic habitats in the Americas. The major vectors of human Chagas disease thrive in human dwellings and peridomestic structures housing domestic animals. (Peri)domestic populations of the major vector Triatoma infestans differ widely depending on the specific local habitat and host species [2]. Here "(peri)domestic" refers to structures that are domestic or peridomestic. We show that Taylor's law (TL), which we describe below, describes well the average and variance of habitat-specific relative population sizes of T. infestans and three other vector species of T. cruzi in all (peri)domestic habitats. The data result from 33,908 habitat searches for triatomine bugs in four areas of Argentina from 1993 to 2010 before and after the large disturbance caused by community-wide insecticide spraying directed to suppress (peri)domestic infestations with T. infestans. One area, well described by TL, had moderate insecticide resistance that caused vector control failures. We determine the effect of insecticide spraying or the history of chemical control interventions on the values of the parameters of TL and describe some implications of TL for vector control and surveillance. The present paper may be the first to demonstrate the connection of TL with any aspect of Chagas disease, and in particular with the population densities of the insect vectors of the disease.

Taylor's law
Large spatial and temporal fluctuations in the population density of living organisms have profound consequences for biodiversity conservation, food production, pest control and disease control, especially vector-borne disease control. In empirical studies of insects and many other species, the sample variance v and the sample mean m of population counts or other measures of population density or abundance approximate well a linear relationship on log-log coordinates, log 10 v % a + b Ã log 10 m [3], which is mathematically equivalent to the power law v % 10 a m b . From a mathematical point of view, the slope b of TL is the proportional (or percentage) rate of increase in the variance for a given infinitesimal proportional (or percentage) increase in the mean. The slope b has been interpreted as an index of spatial aggregation because purely random (Poisson) distributions of individuals have a variance equal to the mean and therefore would be expected to generate TL with b = 1, while some distributions in which the variance grows faster than in proportion to the mean would be expected to generate TL with b > 1. Variations in habitat suitability and other ecological mechanisms could also generate TL with b > 1. (See Future research in the Discussion.) Although Taylor was not the first to publish empirical examples of the above linear relationship, he made it widely known [4,5], and it is usually called Taylor's law (TL) among ecologists, or fluctuation scaling or large-number scaling among physicists [6]. More than 1000 papers have been published on TL and its applications to hundreds of species and many fields besides ecology [6], including weekly cases of measles in 366 communities in England and Wales preand post-vaccination [7], the aggregation of parasite individuals within host individuals (not including any parasites, vectors, or hosts related to the transmission of Chagas disease) [8,9], human population densities [10], crop yields [11], prime numbers [12] and tornado outbreaks [13]. TL can be generated by many different models (e.g., [6,7,[14][15][16]).
TL has important applications in the management of agricultural pests and fisheries. When TL is valid, TL can be used to design more efficient sampling schemes to estimate pest density and decide whether to spray pesticides or release natural enemies in a timely fashion [17,18]. TL provides a stopping rule for fixed precision sampling of fisheries, permitting reduced sampling effort [19].
The uses of TL to identify unusual variability in crop yields [11] and plan more efficient control measures by recognizing the heteroskedasticity of population densities at different mean densities are potentially valid for controlling the insect vectors of major human infectious diseases, including malaria, dengue, Chagas disease, sleeping sickness and the leishmaniases. However, literature searches in Pubmed and Google Scholar (October 16, 2017) using "Taylor's law" (or "Taylor's power law") combined with "malaria mosquito (or Anopheles)", or "dengue mosquito (or Aedes aegypti)", or "tsetse fly (or Glossina fly)", or "Chagas vector (or Triatoma)", or "Leishmaniasis sandfly (or Lutzomyia or Phlebotomus)", identified no paper on Chagas disease vectors and TL and only a few papers which mainly used TL for sample size determination of malaria and dengue mosquitoes [20][21][22][23][24][25][26][27][28][29].
Two widely tested forms of TL are a temporal TL and a spatial TL. In a temporal TL, n populations labeled i = 1, . . ., n are followed over time, and the sample mean size (averaged over time) m i of population i and the sample variance of population size (over time) v i of population i are calculated separately for each population i. Each population is represented by one dot associated with population i on a plot of log 10 v i (vertical axis) as a function of log 10 m i (horizontal axis). If the dots fall approximately along a straight line, the data support a temporal TL.
In a spatial TL, which we pursue here, different populations of a species are grouped into different categories. In this article, each category will be a different habitat in which Chagas vectors may be found, such as a chicken coop or a goat corral. Habitats are labeled h = 1, . . ., H, where H is the number of different habitats. The mean m h and the variance v h of population sizes over all sites of habitat h (e.g., over all chicken coops in a community) are calculated and log 10 v h is plotted as a function of log 10 m h , with one data point for each habitat h. If the H dots fall approximately along a straight line, the data support a spatial TL. efforts were conducted. We use "habitat" for a category of individual places that were surveyed for bugs. For example, chicken coops are one habitat, goat corrals are another habitat, and cow corrals are a third habitat. We use "site" for a particular exemplar of a habitat, such as a particular chicken coop, or a particular goat corral. A "house compound" consists of a domicile for people and near-by buildings for human use and corrals for animals. Each such domicile and building is one site. As indicated above, "(peri)domestic" habitats include all such structures.
Here we demonstrate that TL describes the spatial distribution (in different sites of a habitat) of four of the vector species of Chagas disease. When the means and variances of the number of each species of vector are computed over sites separately for each habitat in a community of house compounds, TL is confirmed with high accuracy and consistency over time and under diverse control procedures. The slope b of TL does not deviate significantly from the range 1 < b < 2. We develop simple mathematical models to help interpret and extend this primary empirical finding. We suggest some practical consequences and potential uses of TL in Chagas disease vector control. The full implications of TL for Chagas disease vector control remain to be worked out in future research and practice and are not the primary objective of this paper. Finally, we suggest some future research.

Methods and materials Field data
The data come from four large research projects in the Argentine Chaco region where Chagas disease was endemic. These projects aimed primarily to control the major vector Triatoma infestans, but also included observations of other local triatomines not considered as the main control targets. The surveys were conducted in well-defined rural areas of Olta (municipalities of General Belgrano and Chamical, of the province of La Rioja, western Argentina), Figueroa and Amamá (Figueroa and Moreno departments, respectively, of the province of Santiago del Estero, northwestern Argentina), and Pampa del Indio (General San Martín department, of the province of Chaco, northeastern Argentina). The studies were organized spatially in a hierarchy with five levels: Argentine Chaco region; four study areas within the region; villages within each area; house compounds (defined above under "Purposes") within each village; and sites within each house compound. The details of each area are described extensively in Detailed Methods in S1 Text. Fig 1(A) maps the areas of these studies and Fig 1B, 1C and 1D illustrates the hierarchy of villages, house compounds, and sites. Table 1 summarizes the quantities of the data collected.
Data sets. S1 Table summarizes, for T. infestans only, the sample means and sample variances of relative bug population density, by study, by survey, and by habitat. S2 Table gives the raw data for all species. From these data, we extracted 83 sets of data from which we estimated a slope and intercept of TL (Results, Table 2), as described below. From these 83 regressions, we selected 79 with sufficient observations to support analysis of the variability of parameter estimates of individual species.

Statistical analysis
Fitting and testing Taylor's law: Summary. We fitted TL and tested its adequacy as a description of the data in three steps. First, for every habitat h in an area (e.g., domiciles, kitchens, chicken coops, pig corrals), we computed the sample mean m h and the sample variance v h of the number of bugs in each site of that habitat per standardized search effort (as described in Detailed Methods) and executed ordinary least-squares linear regression of log 10 v h on log 10 m h , h = 1, . . ., H, across all habitats in the area, separately for each triatomine species, survey, and area. Second, we tested for curvature in the relation of log 10 v h to log 10 m h by fitting a quadratic regression log 10 v h = a + b Ã log 10 m h + c Ã (log 10 m h ) 2 by least squares. If the confidence interval of the coefficient c did not include 0, the data rejected TL because there was statistically significant evidence of curvature. In most cases, there was no statistically significant evidence of curvature, and we examined the residuals of the linear regression models for heteroskedasticity, normality, skewness and kurtosis. Third, when the analyses in steps 1 and 2 did not reject TL as a description of the data, we used analysis of covariance (ANCOVA) to test for differences in the parameters of TL fitted to different subsets of the data. The following subsection gives the details of each step and the standard software used.
Fitting and testing Taylor's law: Details. The fit of the data to TL was tested by ordinary least-squares linear regression of log 10 v h on log 10 m h , h = 1, . . ., H, separately for each triatomine species, survey, and area. Regression analyses were performed with Stata 14.2 [37]. The mean and variance of relative bug abundance for a given triatomine species, habitat, survey, and area were calculated over all identified sites that had been examined for infestation at a given point in time. Analyses included all habitats with mean abundance greater than 0, possibly including some individual sites with abundance equal to 0. Taylor et al. ( [38], p. 721) suggested that at least 15 observations (here sites) should be available to calculate each mean and variance (here, for a given habitat) and that the linear regression (here, for a given species, survey and area) should include at least 5 paired data of v h and m h . The data used here nearly always complied with these suggestions. For example, the Amamá longitudinal data on T. guasayana from the core postintervention surveys of October 1993 and May 1999 are included among the 83 regression estimates in Results, Table 2, but are omitted from further analyses   Dataset: study area. For Amamá , "comparative" compares "core" villages, which had sustained vector control, with "periphery," outlying villages with pulsed vector control. "Longitudinal" presents successive surveys of the core. Triatoma species: "T." = "Triatoma." Intervention: date of survey and whether the survey preceded or followed community- because only two or three habitats had sufficient information to estimate means and variances, giving zero or one df for estimation of TL. For Amamá core and periphery, we included a regression in Table 2 for the pooled data from both zones for each of the triatomine species. To pool the data, we collapsed all the raw data into one file that did not distinguish between zones. For Amamá, each habitat appeared exactly once, even if some habitats may have appeared in one zone and not in the other.
For Olta before and after intervention, Table 2 includes a regression for the pooled data from all of the triatomine species combined. To pool the data, we collapsed all the raw data into one file that did not distinguish between species. However, we excluded these two regressions from further analysis because the regressions were not species-specific. Our further analyses rest on the remaining 79 = 83-2-2 regressions used to estimate and test TL.
In a second step, we tested for curvature in the relation of log 10 v h to log 10 m h by fitting a quadratic generalization of TL originally suggested by Taylor et al. [39, p. 388, their equation (14)] and widely used since: log 10 v h = a + b log 10 m h + c(log 10 m h ) 2 . If the confidence interval of the coefficient c did not include 0, the data rejected TL because there was statistically significant curvature.
Residuals of the linear regression models were tested for normality, skewness and kurtosis using the commands swilk, estat hottest and estat imtest. The program swilk carries out the Shapiro-Wilk test; estat hottest performs three versions of the Breusch-Pagan [40] and Cook-Weisberg [41] test for heteroskedasticity vs homoskedasticity, and estat imtest performs an information matrix test for the regression model and an orthogonal decomposition into tests for heteroskedasticity, skewness and kurtosis [42]. Each residual of the TL linear regression measured the stability of population abundance in the corresponding habitat following [11]. In most cases of fitting the quadratic generalization of TL, the quadratic term was not significant, so these tests of the normality, skewness and kurtosis of the residuals were not performed for the quadratic regressions.
In a third step, when the analyses in steps 1 and 2 did not reject TL, we used analysis of covariance (ANCOVA, implemented in the anova command) to test for differences in the parameters of TL fitted to different subsets of the data. For example, when TL described acceptably the relationship between log sample mean and log sample variance of relative population density among different habitats of two or more species separately, we used ANCOVA to examine whether one or both of the parameters (slope and intercept) of the species-specific TLs differed between species. This ANCOVA treated "species" as a categorical variable and asked whether "species" or the interaction term "species × log sample mean" significantly influenced log sample variance. If "species" influenced log sample variance but not the interaction term, then the intercept of TL differed between species. If the interaction term influenced log sample variance, then "species" affected the slope. If both "species" and the interaction term influenced log sample variance, then both the intercept and the slope of TL depended on the species. We also used ANCOVA to test whether the parameters of TL for a given species differed before and after spraying of insecticides, or according to the history of control measures (sustained versus pulsed).
To compare estimates of slope under two conditions, we used Welch's t-test for two quantities with unequal variances [43]. These calculations used Matlab Version 9.2.0.556344 [44].

Theory
Suppose that Taylor's law (TL) describes well the relation between the mean and the variance of relative population density of a single vector species in the habitats of a study area before the house compounds (including all (peri)domestic structures) are sprayed with insecticides to kill the vectors. What would we expect to be the effect of spraying? Specifically, would we expect TL to hold after spraying? If so, what if any connection should we expect between the intercept and slope of TL before and after spraying? Here we propose two simple models to answer these questions. In the Results, we will compare some of the predictions of these models with observations. It is not necessary to follow the mathematical details to understand the models' predictions or the empirical results. S1 Text gives mathematical proofs.
Both models use the same general notation. Suppose there are H > 2 habitats, such as chicken coop; open shed; oven; piled materials; cow corral; latrine/bathroom; etc. These habitats are labeled h = 1, 2, . . ., H. Let B(h) be a random variable representing the number of vectors of a single species (not all Triatoma species combined) in the various sites of habitat h in the study area Before spraying, and let A(h) be a random variable representing the number of vectors in the various sites of habitat h in the study area After spraying. Because of the gap in time between the survey before spraying and the survey after spraying, the set of sites of a given habitat before spraying may differ from the set of sites of that habitat after spraying.
The population mean (or expectation) of B(h) will be written E(B(h)) and the population variance, Var(B(h)); likewise, the population mean E(A(h)) and population variance Var(A(h)) of A(h). We assume the population mean and population variance exist and are positive. For the vectors before spraying, the log-log form of TL using the population mean E(B(h)) and the population variance Var(B(h)) instead of the corresponding sample mean m and sample variance v is log 10 Var(B(h)) = a + b log 10 E(B(h)). This linear form is mathematically equivalent to the power-law form of TL, namely, The value of the slope b in the log-log form of TL is identical to the value of the exponent b in the power-law form, hence we use the same notation b and we refer to b interchangeably as the slope or the exponent. The value of the intercept a in the log-log form is related to the value of the coefficient C in the power-law form by 10 a = C or log 10 C = a, hence we use different words and symbols (intercept a versus coefficient C).
We assume that spraying reduces the relative population density of the vector, and does not increase it.
Model 1: Constant survival proportion after spraying. Suppose that a fraction s, where 0 < s < 1, of vectors survive spraying. The letter s was chosen as a mnemonic for Survive Spraying. Model 1 assumes that this fraction s is the same for every site of a given habitat (e.g., for every chicken coop in the study area) and for all habitats (e.g., all chicken coops, cow corrals, etc.).
Then for every habitat h = 1, 2, . . ., . Now if the relative population density of vectors satisfies TL before spraying, namely,

then substituting TL into the prior equation and using E(B(h))
TL holds exactly for the relative population density of vectors after spraying with the same exponent b but the coefficient C before spraying changes to s 2-b C after spraying.
In many, but not all, prior empirical studies of insect populations, TL has been confirmed with b < 2. If b < 2, then 2 -b > 0 and hence s 2-b < 1. Then the coefficient s 2-b C of TL after spraying should be smaller than the coefficient C of TL before spraying. On the other hand, if b > 2, then 2 -b < 0 and hence s 2-b > 1. The coefficient s 2-b C of TL after spraying should then be larger than the coefficient C of TL before spraying. If b = 2, the coefficients before and after spraying should be identical.
Model 1 gives six testable predictions.
(1). For every habitat h = 1, 2, . . ., H, we have E(A(h)) = sE(B(h)). This equation relates population means before and after spraying. If we plot the sample mean relative population density of vectors in habitat h before spraying on the horizontal axis and the sample mean relative population density of vectors in habitat h after spraying on the vertical axis, the data points should fall approximately along a straight line through the origin with slope s, apart from sampling variability in both the horizontal and the vertical coordinates of each point.
(2). For every habitat h = 1, 2, . . ., H, we have Var(A(h)) = s 2 Var(B(h)). This equation relates population variances before and after spraying. If we plot the sample variance of relative population density of vectors in habitat h before spraying on the horizontal axis and the sample variance of relative population density of vectors in habitat h after spraying on the vertical axis, the data points should fall along a straight line through the origin with slope s 2 , apart from sampling variability in both the horizontal and the vertical coordinates of each point.
(3). Because we assumed 0 < s < 1, it follows that s 2 < s, so the slope of the line for the sample variances of relative population density before and after spraying should be smaller (lower) than the slope of the previous line for sample means of relative population densities before and after spraying.
(4). The sample means and the sample variances of relative population density before and after spraying should both obey TL if either one does.
(5). The slope or exponent b should remain unchanged before and after spraying, apart from sampling variability in the estimates of b. (6). The coefficient of the power-law form of TL should change from C before spraying to s 2-b C after spraying. Apart from sampling variability in the estimates of s and the parameters of TL, this coefficient after spraying will be smaller than, equal to, or larger than the coefficient of TL before spraying according as b < 2, b = 2, and b > 2.
If model 1 is supported by data, then s provides a valuable statistical summary of the overall effectiveness of spraying. The smaller s is, the more effective the spraying was.
In testing prediction (1) above, if most of the data points (E(B(h)), E(A(h))), h = 1, 2, . . ., H, fall along a straight line E(A(h)) = sE(B(h)) but one or two points fall far above the line, then the control for the habitats corresponding to these outliers needs to be strengthened. If most of the data points (E(B(h)), E(A(h))), h = 1, 2, . . ., H, fall along a straight line but one or two points fall far below the line, then the control for the habitats corresponding to these outliers is unusually effective or the vector populations in those habitats are unusually vulnerable, and some control effort for these habitats might be directed to other more difficult habitats.
In testing prediction (2) above, if most of the data points (Var(B(h)), Var(A(h))), h = 1, 2, . . ., H, fall along a straight line Var(A(h)) = s 2 Var(B(h)) but one or two points fall far above the line, then spraying is not being uniformly applied to all sites of the habitats corresponding to these outliers or the vector populations are variably vulnerable to spraying for reasons that need to be determined (for example, because of differences in the physical complexity of sites considered to belong to the same habitat). If most of the data points (Var(B(h)), Var(A(h))), h = 1, 2, . . ., H, fall along a straight line but one or two points fall far below the line, then these habitats may lack established bug colonies. Further investigation would be required to determine whether such habitats have few adult bugs and the offspring cannot develop a new colony, leading possibly to Poisson variation among sites. This scenario is likely for the sylvatic triatomine species in certain peridomestic habitats.
Model 2: Random survival proportion after spraying. Suppose that the fraction of vectors that survive spraying at a particular site of habitat h is a random variable S(h) where 0 < S(h) < 1. Model 2 assumes that this survival fraction may differ among sites of each habitat h (e.g., may be different for every chicken coop in the study area) and has a distribution that depends on the habitat h (e.g., the distribution for chicken coops may differ from the distribution for goat corrals). Also assume that S(h) and B(h) are independent for all habitats h = 1, 2, . . ., H.
Then for every habitat h = 1, 2, . . ., H, the number of vectors that survive spraying at a particular site of habitat h is A(h) = S(h)B(h). Both sides of this equation are random variables with values that vary among the sites of habitat h. By the assumed independence between S(h) and B(h), the "product rule" holds: The We also prove in the Detailed Methods in S1 Text that in model 2, if b = 2 and the coefficient of variation of S(h) (across sites of habitat h) is the same for every habitat h, then Var(A(h)) and E(A(h)) satisfy TL with b = 2. The assumption that b = 2 in TL before spraying means that the coefficient of variation (across sites) of relative population density is the same for all h.
Because the probability distribution (including its mean and variance) of the spraying survival fraction S(h) is assumed to vary from habitat to habitat, model 2 gives fewer testable predictions than model 1. Assume the means and the variances of relative population density before spraying obey TL.

Results
With m the sample mean of the relative bug population density of each habitat with at least one infestation detected and v the corresponding sample variance, Table 2 gives the estimated intercept a and slope b of Taylor's law (TL) log 10 v % a + b × log 10 m, the standard errors of a and b, the adjusted R 2 and several associated statistical tests of the linear regression, and the minimum, maximum, and range (over all habitats with at least one infestation detected) of log 10 m, for each study, vector species, intervention status, and survey (83 cases).

Amamá
Testing Taylor's law. The data from Amamá confirmed TL. The log-mean relative abundance of T. infestans in (peri)domestic habitats was highly significantly correlated with the logvariance of relative abundance over the pooled core and periphery (b = 1.633 point estimate ± 0.089 standard error, a = 1.310 ± 0.082) (Fig 2A, Table 2). Residuals showed no significant deviations from normality, homoskedasticity and normal kurtosis. Adding a quadratic term to each of the equations did not significantly improve the fit of the models. Similar results held for T. guasayana and T. garciabesi in the data from the pooled core and periphery (Fig 2B  and 2C, Table 2), except for marginally significant results for T. guasayana in the tests for heteroskedasticity and the information matrix, which could easily have been type 1 errors.
Comparing interventions. The Amamá area included two groups of rural communities that had different histories of insecticide control against T. infestans. The core had sustained control. The periphery had pulsed control. S1 Text describes the differences between the core and periphery in greater detail. The data from Amamá do not argue for or against Model 1 because neither region was before spraying or after spraying. Nevertheless, it is informative to compare the core and periphery by the same approach that we will use to compare the distribution of vectors before and after interventions in Olta, Figueroa, and Pampa del Indio.
ANCOVA found no significant differences in the slopes and intercepts of TL for T. infestans between the core (b = 1.610 ± 0.112, a = 1.070 ± 0.105) and the periphery (b = 1.576 ± 0.106, a = 1.080 ± 0.087), indicating no significant effects of the history of insecticide control on TL for this species (Table 2). However, different triatomine species had highly significantly different parameters of TL (Fig 2A-2C).
For T. garciabesi (Fig 2C), the mean and variance of abundance were larger under sustained control (core, open circle) than under pulsed control (periphery, filled circle), as illustrated for the main habitat of T. garciabesi, trees with chickens (ckt). The core had twice as many trees with chickens as the periphery. Moreover, inter-house distances between houses were lower and the proportion of houses raising chickens was higher in Amamá village than they were in rural villages of the periphery, deeper in the forest. Apparently these differences favored the apparently limited dispersal capacity of T. garciabesi in the core.
For Amamá, Fig A(A) in S1 Text gives, for each habitat (individual data points), the total number of T. infestans individuals and the mean (Fig A(B) in S1 Text) and the variance (Fig  A(C) in S1 Text) of the number of T. infestans individuals per site of each habitat, in the core (horizontal axis) and in the periphery (vertical axis). We give two least-squares regressions for each summary statistic, and focus here on the regression through the origin, since in each panel the horizontal and vertical axes measure the same quantity under different conditions.
On the average across habitats, the total number of T. infestans in all sites of a habitat in the periphery was about 23% more than the total number of T. infestans in all sites of a habitat in the core (Fig A(A) in S1 Text). It cannot be concluded from this comparison that the pulsed surveillance in the periphery is inferior, as a method of reducing T. infestans populations, to the sustained surveillance in the core, because these data do not include a baseline of the initial T. infestans populations in the core and periphery prior to any control. Thus we cannot compare bug population sizes before and after control in the two regions.
The mean number of T. infestans per site of each habitat in the periphery was very loosely (adj. R 2 = 0.37) related to the mean number of T. infestans per site of each habitat in the core (Fig A(B) in S1 Text). If there was any relation at all, the mean in the periphery was about twice as great as in the core. Likewise, the variance of the number of T. infestans per site of each habitat in the periphery was very loosely, or hardly at all, related to the variance of the number of T. infestans per site of each habitat in the core (Fig A in S1(C) Text). Storerooms were the outlier in the upper left corner of the plot. If there was any relation at all, the variance in the periphery was about 64% greater than in the core.
The regressions above included only habitats that had a mean greater than 0 in both core versus periphery, to be consistent with the following analyses in Olta, Figueroa, and Pampa del Indio, which included only habitats that had a mean greater than 0 both before and after interventions. This constraint reduced the number of data points in the above regression lines. Including all 15 data points (regardless of whether one of the means was 0) reduced the nointercept regression slope of the mean abundance in the periphery to the mean in the core from 2.0060 (± 0.7969) to 1.5825 (± 0.5928, adj. R 2 = 0.290) and the slope of the variance from 1.6387 (± 2.9043) to 1.1610 (± 1.8592, adj. R 2 = -0.042).
Core area longitudinal surveys. In the 13 surveys from October 1993 to October 2002 of the Amamá core under sustained vector surveillance and control, the mean and variance of the relative abundance of T. infestans bugs in each habitat obeyed TL (Fig 3A and 3B, Table 2). Eleven of the 13 tests rejected the null hypothesis that a = 0, and all point estimates of a were positive. Of the 52 = 4 × 13 tests of various assumptions of the linear model, only 3 had P < 0.05, and 3/52 = 0.058 was close to the fraction 0.05 expected by chance alone. The linear model of TL was not rejected. For 10 of the 13 surveys, the adjusted R 2 of TL exceeded 0.9. In general, TL was a plausible model. Results were similar for the two other species in the Amamá core ( Table 2).
The intercept a and the slope b of TL fluctuated to some extent but displayed no systematic trend over time (Fig 3C, Table 2). The point estimate of the slope reached a peak of b = 2.002 ± 0.720 at survey 7 (May 1997). The greater the range (maximum log 10 mean minus minimum log 10 mean) of the mean numbers of bugs by habitat, the smaller the fluctuations of b around the central value of roughly 1.5 ( Fig 3D, Table 2). Of the 13 estimates of slope, 10 fell between 1.4 and 1.65. Seasonal and secular changes during this decade had little effect on the form of TL or its parameters.
For the less abundant bug species, T. guasayana and T. garciabesi, the patterns were similar (Table 2), though with slightly more fluctuations.

Olta
Testing Taylor's law. The data from Olta confirmed TL. Before community-wide insecticide application, the log-mean relative abundance of T. infestans, T. guasayana and T. garciabesi in peridomestic habitats was highly significantly correlated with the log-variance of bug abundance when the three species were taken together (adj. R 2 = 0.981), with a slope (b = 1.504 ± 0.045) suggesting significant insect aggregation or differences in habitat suitability across habitats (Table 2). Approximately the same patterns were recorded when each species was taken separately before interventions: T. infestans (b = 1.233 ± 0.150; Fig 4A), T. guasayana (b = 1.428 ± 0.167; Fig 4B), and T. garciabesi (b = 1.558 ± 0.190; Fig 4C). All intercepts were significantly different from 0 whereas the quadratic term was not statistically significantly different from 0 for each triatomine species both before and after insecticide spraying ( Table 2). Residuals showed weak deviations from normality, homoskedasticity and normal kurtosis for the three species.
Community-wide insecticide application did not destroy the linear relationship between the log-variance and log-mean of bug abundance for the three triatomine species taken together (b = 1.646 ± 0.240, a = 1.037 ± 0.083) ( Table 2). The slopes of log-variance to logmean bug abundance did not differ significantly among triatomine species either before or after insecticide spraying, nor when each species was compared separately before versus after spraying. Intercept values before and after insecticide spraying were not significantly different.
Effects of spraying: Testing predictions 1, 2, 3. For Olta, Fig 5 gives, for each habitat (individual data points), the total (A), the mean (B), and the variance (C) of the number of T. infestans individuals per site of each habitat, before spraying (horizontal axis) and after spraying (vertical axis). The number of sites of each habitat was identical before and after spraying (S1 Table). We focus here on the regression through the origin, as above.
Both the mean (Fig 5B) and the variance (Fig 5C) after spraying were smaller than their values before spraying, as predicted by Model 1. They were loosely related to the mean and the  variance (adj. R 2 = 0.550, 0.259), respectively, before spraying, so predictions 1 and 2 of model 1 were weakly confirmed. The mean after spraying was approximately 23% of the mean before spraying, averaged over habitats by the linear regression. The variance after spraying was approximately 19% of the variance before spraying, averaged over habitats by the linear regression. Only for chicken trees was the variance after spraying greater than the variance before spraying (S1 Table). That the slope of the regression for the sample variances of relative population density before and after spraying was smaller than the slope of the regression for sample means of relative population densities before and after spraying confirms prediction 3 qualitatively. However, the data did not confirm the quantitative prediction that the slope of the variance regression should be the square of the slope of the mean regression. The slope of the   (Table 2) confirmed prediction 4 (the means and the variances of relative population density before and after spraying should both obey TL if either one does) of model 1. For T. infestans and the other two vector species each considered separately and for all three species combined, both pre-and postintervention, there was no strong (P < 0.01) evidence to reject TL.
The data from Olta (Table 2) also confirmed prediction 5 (the slope or exponent b should remain unchanged before and after spraying). For T. infestans, the preintervention slope and standard error were b = 1.233, SE(b) = 0.150. The postintervention slope and standard error were b = 1.767, SE(b) = 0.347. The Welch test gave little or no evidence (P % 0.098) that the difference of slopes was significantly different from 0 (S3 Table). This conclusion is evident from inspection, since the difference between the pre-and postintervention slopes, 1.233-1.767 = -0.534 is not even twice the postintervention SE(b). The same holds for the other two vector species and for all three species combined.
The data from Olta (Table 2) also confirmed qualitatively prediction 6: the coefficient of the power-law form of TL should change from C before spraying to Cs 2-b after spraying. With b < 2, the coefficient of TL after spraying was smaller than the coefficient of TL before spraying. In Olta, for T. infestans, b < 2 both pre-and postintervention. As predicted, the preintervention a = 1.100 was (very slightly) larger than the postintervention a = 1.095. The same direction of difference held for the other two vector species individually. However, for all three species combined, the preintervention a = 0.995 was slightly (but not significantly) smaller than the postintervention a = 1.037.

Figueroa
Testing Taylor's law. The data from Figueroa confirmed TL. The preintervention logmean relative abundance of T. infestans in (peri)domestic habitats of Figueroa rural houses was highly significantly correlated with the log-variance of bug abundance (adj. R 2 = 0.919). The slope b = 1.401 ± 0.146 suggested significant insect aggregation or differences in the variation of suitability across habitats ( Table 2, Fig 6A). Residuals showed no significant deviations from normality, homoskedasticity and normal kurtosis. Mean bug abundance ranged over nearly two orders of magnitude. Bugs were very rare in latrines at the extreme left, which had very few bugs collected among the 115 sites inspected for infestation and only exceptionally had a bloodmeal host (chickens) inside or leaning against latrine walls. Latrines do not appear to harbor established bug populations. When both data points for latrines were suppressed, the point estimate of the slope (±standard error) rose to 1.494 ± 0.413 and adj. R 2 fell to 0.634. Regardless of whether latrines were included or not, the confidence interval for the coefficient c of the quadratic term in the generalized TL included 0, indicating that the linear form of TL was not rejected by the data.
Community-wide insecticide spraying did not significantly affect TL for T. infestans except possibly at 5 months post-spraying (March 2004) (Fig 6A). Slopes for T. infestans increased from 1.401 ± 0.146 before spraying to 2.092 ± 0.225 at 5 months post-spraying (March 2004) (S3 Table gives  Residuals showed no significant deviations from normality, homoskedasticity and normal kurtosis after spraying. For T. garciabesi and T. guasayana, the preintervention slopes of each did not differ significantly from the respective postintervention slopes (P > 0.10) (Fig 6B and 6C).
Effects of spraying: Testing predictions 1, 2, 3. For Figueroa, Fig B in S1 Text gives, for each habitat (individual data points), the total, the mean and the variance of the number of T. infestans individuals per site of each habitat, before spraying (October 2003) (horizontal axis) and after spraying (October 2004) (vertical axis). The number of sites of each habitat was substantially different before and after spraying (S1 Table). We eliminated three habitats in which the mean and variance were 0 before or after spraying (chicken trees, chicken houses, and open sheds), leaving eight habitats. We focus here on the regression through the origin, as above.
Both the mean (Fig B(B) in S1 Text) and the variance (Fig B(C) in S1 Text) after spraying are usually smaller than their values before spraying, as assumed by model 1. They are only loosely related to the mean and the variance, respectively, before spraying, so predictions 1 and 2 of model 1 are at best weakly confirmed. The mean after spraying was much smaller than that before spraying, approximately 9% of the mean before spraying. The variance after spraying was approximately 2% of the variance before spraying. That the slope of the regression for the sample variances of relative population density before and after spraying was smaller than the slope of the regression for sample means of relative population densities before and after spraying confirmed prediction 3 qualitatively. The quantitative prediction that the slope of the variance regression should be the square of the slope of the mean regression was (at best)  roughly, or perhaps not, confirmed, as the slope of the variance regression, 0.0160, differed from the squared slope of the mean regression, 0.0928 2 = 0.00861 by a factor of 2.
TL and spraying: Testing predictions 4, 5, 6. All three species observed in Figueroa (T. garciabesi, T. infestans and T. guasayana) satisfied TL before spraying (in the survey of October 2003), confirming prediction 4. Only for T. infestans do we have enough data to test TL in the first survey after spraying (March 2004), and for this species and survey, TL held (further confirming prediction 4).
For T. infestans, the preintervention slope b = 1.401 ± 0.146 was marginally significantly smaller than the immediate (March 2004) postintervention b = 2.092 ± 0.225, as noted above (P = 0.016) according to the Welch test (S3 Table, Table) from the preintervention b, again in agreement with model 1's prediction 5, but the quadratic term c differed significantly from zero (P = 0.006, Table 2 Table 2).

Pampa del Indio
Testing Taylor's law. The data from Pampa del Indio confirmed Taylor's law. The logmean relative bug abundance in (peri)domestic habitats was highly significantly correlated with the log-variance of bug abundance before insecticide applications for T. infestans (b = 1.601 ± 0.112, a = 1.384 ± 0.105) and T. sordida (b = 1.671 ± 0.115, a = 1.249 ± 0.135) ( Table 2, Fig 7A and 7B). The slopes suggest a similar degree of aggregation or diversity of habitat suitability across habitats and triatomine species. Intercepts were highly significantly different from 0. Residuals deviated significantly from homoskedasticity only for T. sordida according to the Breusch-Pagan/Cook-Weisberg test but not by Cameron and Trivedi's test. Adding a quadratic term to each of the equations did not significantly improve the fit of the models.
The slopes of TL postintervention varied more widely for T. infestans (range, 1.356-2.145) than for T. sordida (range, 1.316-1.748), possibly because the former had a very much reduced number of infested habitats over the last three surveys (Table 2). Intercepts postintervention varied from 1.224 to 2.460 for T. infestans, and from 0.778 to 1.542 for T. sordida. For T. infestans, the extreme high values of a and b occurred at survey 6 (October 2009) (Fig 7C) and were associated with the smallest range (0.564) between the maximal log 10 mean of relative bug abundance and the minimal log 10 mean of relative bug abundance (Fig 7D). These results illustrate the general pattern that estimates of the slope of any linear regression are unstable when the range of x-axis values is small. A larger range was generally associated with less variable point estimates of the slope. The eight surveys of Pampa del Indio made possible 64 = 8 × 8 pairwise comparisons between the two triatomine species of the slopes of TL. According to the Welch test (S3 Table), the slopes differed significantly between species only once with P < 0.01 and in only 7 of 64 comparisons with P < 0.05, giving no compelling evidence of difference in slopes between species. This detailed evidence for Pampa del Indio is consistent with the lack of strong evidence for a difference in slopes between these two species in Table 2. Effects of spraying: Testing predictions 1, 2, 3. For Pampa del Indio, Fig C in S1 Text gives, for each habitat (each individual data point), the total, the mean and the variance of the number of T. infestans individuals per site of each habitat, before spraying (survey 1) (horizontal axis) and after spraying (survey 3) (vertical axis). The number of sites of each habitat was substantially different before and after spraying (S2 Table). We included here all habitats that had bugs before or after spraying, though we did not use all of these habitats to estimate TL (which requires nonzero means and nonzero variances). We focus here on the regression through the origin, as above.
Both the mean (Fig C(B) in S1 Text) and the variance (Fig C(C) in S1 Text) after spraying were smaller than their values before spraying, with the slight exception of latrines (S1 Table), and were visually linearly related (for the mean before and after, adj. R 2 = 0.631, for the variance before and after, adj. R 2 = 0.244), as predicted by Model 1. According to the linear regressions through the origin, the mean after spraying was approximately 13% of the mean before spraying. The variance after spraying was approximately 7% of the variance before spraying. That the slope of the regression for the sample variances of relative population density before and after spraying was smaller than the slope of the regression for sample means of relative population densities before and after spraying confirmed prediction 3 qualitatively. However, the quantitative prediction that the slope of the variance regression should be the square of the slope of the mean regression was at best weakly, or not, confirmed. The slope of the variance regression, 0.0664, exceeded the squared slope of the mean regression, 0.1253 2 = 0.0157, by more than a factor of four.
For both T. infestans ( Fig 8A) and T. sordida (Fig 8B), TL described well the relation of the log spatial sample variance to the log spatial sample mean of vector relative abundance across all habitats at each of eight periodic surveys conducted before (October 2007) and after (April 2008 to October 2010) community-wide insecticide application. During the three years after community-wide spraying, assessments of postintervention bug infestations were coupled with selective sprays of the residual foci detected except for April 2008. The linear relationship between the log-variance and log-mean of bug abundance for all sites and dates was highly significant for T. infestans (b = 1.299 ± 0.063, a = 1.592 ± 0.092) and T. sordida (b = 1.810 ± 0.267, a = 1.979 ± 0.275) ( Table 2, Fig 8). Residuals showed either no or weak deviations from normality, homoskedasticity and normal kurtosis for both species. The slopes of TL differed marginally between triatomine species by the Welch test (P = 0.0559, S3 Table).
TL and spraying: Testing predictions 4, 5, 6. In Pampa del Indio, T. infestans and T. sordida satisfied TL before and after spraying (prediction 4). Though there was some heteroskedasticity preintervention in both studies, there was no statistically significant evidence of quadratic curvature ( Table 2, Fig 7A and 7B).
For T. infestans, the preintervention (October 2007) slope b = 1.601 ± 0.112 did not differ significantly (P = 0.337, S3 Table) from the postintervention (April 2008) slope b = 1.738 ± 0.297, in agreement with model 1's prediction 5. The preintervention a = 1.384 ± 0.105 was smaller than, but did not differ significantly from, the postintervention a = 1.515 ± 0.276. The latter finding would be consistent with prediction 6 if spraying survival s were close to 1 in the subset of habitats with enough sites after spraying to support calculations of the mean and variance.
For T. sordida, the preintervention slope b = 1.671 ± 0.115 was larger, but not significantly larger, than the postintervention slope b = 1.508 ± 0.141, P = 0.190 (S3 Table), consistent with prediction 5. The preintervention intercept a = 1.249 ± 0.135 exceeded, but not significantly, the postintervention intercept a = 1.060 ± 0.200. The direction of the difference is qualitatively in accord with prediction 6 of model 1.
Is the slope b of Taylor's law a fixed attribute unique to each vector species?
Taylor [3,4] estimated that different insect species had different slopes b. He proposed that the slope b was a characteristic unique and specific to each species, meaning that a given species always had the same value of b and no other species had that value of b. Others observed that Taylor   have identical or statistically indistinguishable point estimates of b due to sampling fluctuations alone. Moreover, Yamamura [45] and Cohen et al. [46] showed by examples that the physical scale on which population density was measured in a spatial TL could substantially affect the value of b, so that b was not necessarily an invariant species-specific characteristic.
The question of whether the slope b is a species characteristic (under some range of circumstances) has practical significance for Chagas disease vector control, because it implies further questions. For a given species, must b be measured anew in every new field situation? Must b be measured independently for different vector species?
S3 Table gives the Welch test statistic, df, and P for testing the null hypothesis of no significant difference between every two estimates of the slope b in Table 2 (excluding 4 of the 83 estimates, as we now describe). All these studies shared a similar ecology in the Gran Chaco and similar sampling methods. S3 Table and the following analyses omit two regression estimates in Table 2 from Amamá for T. guasayana in the core postintervention surveys of October 1993 and May 1999, for which df = 0 and 1, respectively, and the two estimates for Olta with three species combined, pre-and postintervention, leaving 79 = 83 -4 species-specific estimates of b. There are 3,081 = 79×78/2 distinct pairs of estimates of b. Of these 3,081 comparisons, 206, or 6.7%, were less than 0.01, nearly seven times more than the approximately 31 = 3,081×0.01 P values less than 0.01 that would be expected on average by chance alone under the null hypothesis that all underlying values of b are identical.
Because each P value describes the difference of two slopes b and the same b occurs in multiple comparisons with all other values of b, we refrain from attempting to assign a probability to the chance that 3,081 dependent trials with marginal probability of success equal to 0.01 would yield 206 cases with P < 0.01 in the Welch test that two b values are equal. Instead, we simply infer that at least some of the differences in b were probably not due to sampling variation alone, and we focus attention on where the frequency of significant differences in b occurred most frequently. Table 3 summarizes the number of cases where the Welch test's P < 0.01 for each pairing of each of the four species of Triatoma with itself and with each of the other three species. According to the diagonal elements of Table 3C, the percentage of P values that were less than 0.01 exceeded 1% by a factor of 2.8-4.7 when the slope for each of the four species was compared with the slope for the same species in another study. Each diagonal element was less than the percentage directly above it or directly to its right, meaning that there were fewer intraspecific than interspecific significant differences in slope. Each species resembled itself in TL slope b more than it resembled any other species, when we compared pairs of distinct studies. The largest percentages of P values that were less than 0.01 occurred in the comparison of b values of T. guasayana with those of the other three species. Overall, T. guasayana had b values that differed the most from those of the three other species, and overall the slope b of TL varied less within any of the four species than interspecifically. The b values from T. sordida and T. garciabesi hardly differed significantly more often interspecifically (3.2%) than they varied intraspecifically (2.8% for both species).

Discussion
We review our most important empirical findings, discuss the findings in light of our mathematical models, suggest potential practical applications, and indicate some future research.
Taylor's law describes the spatial distribution of Chagas disease vectors among habitats With remarkable precision, Taylor's law (TL) described the relationship, among multiple habitats, of the variance of relative bug population density in sites of a given habitat to the mean of relative bug population density in the same sites of the given habitat. The log variance was well approximated as a linear function of the log mean in four studies that differed in geographic location, methods of bug control, time since spraying, degree of insecticide resistance, and principal bug species. The adjusted R 2 had median 0.94 and lower and upper quartiles 0.91, 0.96 in 79 single-species field surveys with df > 1 in Table 2. The slope b of the linear relationship of log variance to log mean differed most between Triatoma guasayana and the other three species, T. infestans, T. sordida, and T. garciabesi. Interspecific differences exceeded intraspecific differences in slope b.
Slopes of TL generally, but not always, rejected Poisson-distributed vectors with different means in different habitats and were consistent with substantial spatial aggregation or differences in habitat suitability: the median slope b was 1.48 and the lower and upper quartiles were 1.35, 1.63 in these 79 field surveys ( Table 2). Only four of these field surveys had b > 2, and none of these slopes was significantly greater than 2. Only four of these field surveys had b < 1, and none of these slopes was significantly less than 1. Thus none of the 79 TL field surveys significantly rejected the general pattern that 1 < b < 2, as has been widely found in other studies. The first inequality, 1 < b, implies that, when different habitats were compared, the variance in their bug populations increased faster than in proportion to their mean bug populations. The second inequality, b < 2, implies that, when different habitats were compared, the coefficient of variation (standard deviation divided by mean) decreased as the mean number of bugs per habitat increased.

Models of control
Model 1 usefully highlights deviations from the predictions that follow from its very simple assumptions. The increased slope of TL for T. infestans observed after spraying in Olta,  Table 2, summarizing results in S3 Table). (A) "Cell count" is the number of pairwise comparisons. For example, there were 9 values of b for T. sordida and 31 values of b for T. infestans, so there were 9×8/2 = 36 distinct intraspecific comparisons of b for T. sordida, 31×30/2 = 465 intraspecific comparisons of b for T. infestans, and 9×31 = 279 interspecific comparisons of b for T. sordida versus T. infestans. (B) "P<0.01 count" is the number of these comparisons that had P < 0.01 according to the Welch's test (S3 Table). (C) "% P<0.01" is the percentage of comparisons with P < 0.01. For example, for the comparisons of b for T. sordida versus T. infestans, 3.2% = 9/279. If differences in b were due to sampling fluctuations alone and there were no intra-or inter-specific differences in the underlying values of b, then "%P<0.01" = (P<0.01 count)/(Cell count) should approximate 1.0%. Figueroa, and Pampa del Indio indicates that spraying lowered mean vector abundance in some habitats more than in others, and hence was not uniformly effective across habitats. This heterogeneity in the effectiveness of spraying is also evident in the plots of log-mean-abundance-after-spraying as a function of log-mean-abundance-before-spraying (Fig 5, Fig B in S1 Text, Fig C in S1 Text). This evidence rejects the classic assumption made by vector-control agencies that insecticide spraying is uniformly effective. Heterogeneity in the effectiveness of spraying calls for targeted, improved vector control, possibly including appropriate environmental management measures. This issue is of paramount importance for the goal of largescale Chagas disease vector elimination under the aegis of regional intergovernmental programs such as the Southern Cone and Central America Initiatives created in the 1990s and ongoing [47,48]. Model 1 has relevance beyond attempts to control Chagas disease vectors. As part of a study of the ecological consequences of sudden oak death in northeastern United States temperate forests, oak trees were censused and measured in 12 plots in 2007 and again in 2010 in Black Rock Forest, Cornwall, New York. Cohen et al. [49] fitted TL successfully to the counts of oak trees in 2007 before a major intervention in 2008 (killing of some oak trees by girdling) and again in 2010. The girdling of the oaks was analogous to spraying the bugs. The acceptable fits of TL before and after intervention were consistent with prediction 4 of model 1, and the absence of statistically significant evidence of a change in slope as a result of the intervention was expected from prediction 5 of model 1. Cohen et al. [49] offered no explanation of why the intervention did not significantly change the slope or destroy the fit of TL. Model 1 offers an explanation, or at least a phenomenological description.

General implications and potential practical applications
Taylor's law identifies key habitats with high mean and variance of infestation before or after insecticide spraying. TL has not previously been used to discover highly variable infestations by any insect vector, to the best of our knowledge. In general, a habitat with exceptionally high variance, given its mean, would likely have a high likelihood of a vector (here bug) outbreak, and might deserve special attention for control by spraying or environmental alteration. A habitat with exceptionally low variance, given its mean, would be a habitat with a relatively stable endemic bug population. For example (Fig 2A), the spatial variance of T. infestans in chicken houses (ckh) in the periphery of Amamá is at or very near the lower limit of the 95% CI given the rather large spatial mean abundance of bugs in chicken houses in the periphery. Even with no spraying of insecticides, one might expect chicken houses with relatively stable chicken populations to harbor spatially consistent bug populations.
TL is useful for practical control efforts also for habitats that are not outliers from the loglog regression. For example, in the graphs of TL based on T. infestans abundance by habitat in Amamá, the data points at the upper right of the graph were granaries and sheds (i.e., open sheds with a thatched roof, frequently used as a storage area, with wide variation in stored contents and usage between households, sometimes with chickens nesting there). Overlooking or failing to inspect or treat adequately granaries and sheds could have large effects on bug persistence and subsequent house reinfestation rates. People rarely permit insecticide spraying of granaries with stored corn (except when the granaries are empty). The preintervention mean bug population size of granaries was large in Figueroa and a high proportion of granaries was infested [2]. Granaries were again at the top (upper right) of the log-log regression before and 5 months after spraying in Figueroa. Open sheds do not get much attention as potentially important habitats for bug control, perhaps because of the wide variability in bug counts between individual sheds. Granaries and open sheds had greater means and variances of bug abundance than chicken coops, storerooms, and domiciles. The latter three habitats usually take all the attention of bug control personnel. Domiciles, storerooms and kitchens tend to appear on the upper end of the log variance-log mean regressions. Because they concentrate competent hosts, these three habitats concentrate nearly all of the T. cruzi-infected bugs in any study area in the Gran Chaco, especially domiciles where all human-vector contacts occur (e.g., [50]).
The same argument applies to other habitats that receive less attention because they are located outside of domestic compounds in areas less disturbed by human activities. Examples include orchard fences where cavies and other rodents thrive and chickens sometimes nest [51], and abandoned constructions covered by vegetation and used by free-ranging goats, sheep or other domestic or sylvatic mammals [52]. Such areas are very rarely inspected for bugs or sprayed by bug control staff on the assumption that they are rarely infested (i.e., implicitly assuming that such habitats have an exceptionally large variance of bug population density) and because substantial time, effort and insecticide would be needed to spray them rigorously. Over the years, we collected several examples of such highly infested, rarely detected, habitats that would fall in the upper right corner of the graph of TL and would severely hamper any serious elimination campaign.
In the practical control of Chagas disease vectors, researchers and bug control people historically have shown little interest in levels of statistical precision or in using sampling theory. In the case of T. infestans, perhaps one reason was that the main initial program goal was to eliminate the vector completely from most of its range, not to control it or estimate its relative population density, and to treat all infested houses completely and homogeneously. Another possible reason was the crude method of detecting or sampling bugs, which is still used. Bug control programs declare house compounds infested or not. Though they may count the number of bugs per unit of search effort at the level of the house compound (distinguishing domestic from peridomestic habitats), they do not use this information other than for saying "many bugs or few bugs" at the community-wide level. Failure to use sampling theory in past practice does not diminish the present and future need to consider TL and its implications for certifying the interruption of transmission of T. cruzi to humans at a district-or state-wide scale, for example. Figs 2-4, 6 and 7 make very clear which habitats have the largest means and variance of bug populations and how habitats respond to insecticide applications. Sometimes a habitat stays at the extreme right (with large mean and large variance of relative population size) after interventions while other habitats jump to the left extreme or remain there before and after interventions. That a habitat stays at the upper extreme after insecticide applications demonstrates the vulnerability of vector control through insecticides and the need for better chemical control or alternative interventions.
The habitat "other" appears several times with high mean abundance and high variance. Because "other" habitats are hard to classify, they may not be identified as important for insecticide applications. This problem is particularly acute for Triatoma dimidiata and other species where selective control is frequent (i.e., only domiciles, or only domiciles and chicken coops, are sprayed, for example). Neglected "other" habitats may be sources of recurrent infestations.
This collection of surveys and trials shows that insecticide-based vector suppression in several study areas across the Argentine Chaco has been far from uniformly successful over the 28 years from 1993 to 2010. The graphical display of means and variances by habitat may provide useful quantitative guidance for improving insecticide-based vector suppression.

Future research
The relative population densities of Chagas vectors reported here strongly confirm TL but do not identify the mechanisms that may generate TL, such as insect behavior versus habitat suitability (measured by bug birth rates and death rates in different habitats). How important are selective dispersal and migration (insect behavior) versus differences between habitats in their average suitability for bugs and in their variation (across sites) of suitability? Control measures that transitorily lower the suitability of a habitat, such as residual insecticide spraying, could in principle stimulate adaptive behavioral responses by insects (e.g., excito-repellency), which may reduce the effectiveness of spraying at the level of the house compound or the village. Thus understanding the mechanisms that produce TL is important.
In these studies, postintervention surveys occurred from 4-6 months to 12 months after insecticide spraying. Hence the observed relative bug population density postintervention included immigration, local recruitment (new births of bugs), and local survival postintervention. Models 1 and 2 of habitat-specific persistence represent only the last of these processes, and therefore cannot be expected to account precisely for our field observations. Models 1 and 2 assume that insecticide spraying does not reduce to zero the average or the variance of relative population density in sites of a habitat. They are models of successful partial control but not models of successful local elimination of bugs from all sites of any habitat. When the ultimate program goal of insecticide spraying is local elimination (suppression) of bugs from all sites of a habitat rather than control, these models may be helpful in identifying less vulnerable habitats using TL.
It would be useful to supplement these models with an empirical summary and theoretical model of the fraction of sites of each habitat from which spraying eliminated the vectors. If the number of bugs per site of a habitat were described by the negative binomial distribution, then the mean and variance would imply a fraction of sites with zero bugs. Efforts have been made to link TL with the mean and the variance of the negative binomial distribution [8,[52][53][54][55][56][57], but it has recently been recognized that TL cannot hold with constant parameters at the same time that the negative binomial distribution holds with a constant scale parameter and a changing probability parameter ( [9], p. E50) (see S1 Text for further details). A practically important topic for further research is finding a useful way to estimate the fraction of sites of a habitat with no bugs, when the mean and variance of bug abundance obey TL, as here.
The units of analysis in this paper are habitats such as chicken coops, kitchens, or granaries. By contrast, vector control programs use house compounds as units to calculate village-wide infestation rates. Additional empirical and theoretical analyses are needed to link habitats with house compounds and village-level summaries of infestation. Moreover, since the goal of vector control is to reduce or eliminate human T. cruzi infections, it would also be highly desirable to explore the potential of TL to shed light on the distribution of T. cruzi infection in vectors, humans and other animal hosts.
Supporting information S1 Text. Contains detailed methods regarding data, statistical analyses, and models (including mathematical proofs of two variance formulas for model 2), detailed results, and four supplementary Figs A-D. (DOCX) S1 Table. A key to habitat abbreviations and, for Triatoma infestans only, the sample means and sample variances of relative bug population density, by study, by survey, and by habitat. S1 Table includes Table. Comparisons of TL slope b for every pair of values in 79 studies, by area, Triatoma species, and intervention. Sheet "WelchT" gives the t-statistic of the Welch test. Sheet "Welchdf" gives the degrees of freedom (df) of the Welch's test. Sheet "WelchP" gives the Pvalue associated with the corresponding t and df. (XLSX)