Environmental and State-Level Regulatory Factors Affect the Incidence of Autism and Intellectual Disability

Many factors affect the risks for neurodevelopmental maladies such as autism spectrum disorders (ASD) and intellectual disability (ID). To compare environmental, phenotypic, socioeconomic and state-policy factors in a unified geospatial framework, we analyzed the spatial incidence patterns of ASD and ID using an insurance claims dataset covering nearly one third of the US population. Following epidemiologic evidence, we used the rate of congenital malformations of the reproductive system as a surrogate for environmental exposure of parents to unmeasured developmental risk factors, including toxins. Adjusted for gender, ethnic, socioeconomic, and geopolitical factors, the ASD incidence rates were strongly linked to population-normalized rates of congenital malformations of the reproductive system in males (an increase in ASD incidence by 283% for every percent increase in incidence of malformations, 95% CI: [91%, 576%], p<6×10−5). Such congenital malformations were barely significant for ID (94% increase, 95% CI: [1%, 250%], p = 0.0384). Other congenital malformations in males (excluding those affecting the reproductive system) appeared to significantly affect both phenotypes: 31.8% ASD rate increase (CI: [12%, 52%], p<6×10−5), and 43% ID rate increase (CI: [23%, 67%], p<6×10−5). Furthermore, the state-mandated rigor of diagnosis of ASD by a pediatrician or clinician for consideration in the special education system was predictive of a considerable decrease in ASD and ID incidence rates (98.6%, CI: [28%, 99.99%], p = 0.02475 and 99% CI: [68%, 99.99%], p = 0.00637 respectively). Thus, the observed spatial variability of both ID and ASD rates is associated with environmental and state-level regulatory factors; the magnitude of influence of compound environmental predictors was approximately three times greater than that of state-level incentives. The estimated county-level random effects exhibited marked spatial clustering, strongly indicating existence of as yet unidentified localized factors driving apparent disease incidence. Finally, we found that the rates of ASD and ID at the county level were weakly but significantly correlated (Pearson product-moment correlation 0.0589, p = 0.00101), while for females the correlation was much stronger (0.197, p<2.26×10−16).


Introduction
Autism spectrum disorders (ASD) are a collection of chronic, complex neuropsychiatric diseases with well-characterized comorbidities and increasing apparent prevalence [1]. With few and limited effective treatments and considerable financial burden, its etiology remains a scientific puzzle. Evidence suggests that autism is highly heritable and clustered within families; consequently, much scientific attention has been dedicated to the discovery of predisposing genetic factors [2][3][4][5][6][7]. There is also evidence for environmental influences, such as prenatal exposure to pesticides or valproate, but it is challenging to account systematically for these factors because they are mostly undocumented. In addition, there are numerous factors that could affect or distort the observed variation in temporal and spatial disease prevalence: evolving diagnostic criteria, socioeconomic, legal, and cultural incentives for diagnosis [8], changing environmental exposures, and the accumulation of genetic burden in the growing human population. However, the relative importance of all these putative causal factors and confounders on ASD prevalence, the nature of interactions between contributing factors, and the underlying biological mechanisms, remain unclear.
Along these lines, geospatial clustering of ASD has been observed in California [9,10], Texas [11,12], North Carolina [13] and Utah [14]. Clustering could indicate the existence of localized risk factors, such as environmental toxins [15] or maternal education [13]. However, studies to date have focused primarily on within-state patterns and socioeconomic predictors, such as the level of parental education, the controversial financial incentives induced by state policies for special-education services, and broad environmental indicators. Now, the increasing availability of large administrative clinical datasets, with national coverage and fine spatial granularity, along with data on possible causal and confounding factors, provides an opportunity to compare the magnitudes of these and other factors within a unified mathematical framework.
Here, we report a mixed-effect Poisson regression analysis of the spatial incidence patterns of ASD and, for comparison, intellectual disability (ID). The data was derived from a very large insurance claims database containing nearly 100 million patients in the United States, which was augmented with census data to introduce additional county-level covariates that captured socioeconomic, demographic, and environmental effects. We present strong statistical evidence for environmental and legal factors driving the apparent spatial heterogeneity of both phenotypes, while documented socioeconomic factors and population structure have much weaker effects.

Results
We analyzed the strength of disparate factors on the apparent incidence rates of ASD and ID by computationally interrogating insurance claims for approximately one third of the US population, using a bivariate-response, three-level, mixed-effects Poisson regression model with 50 free parameters, 44 of which correspond to the fixed effects of known factors while the remaining 6 account for the variance and covariance among the random effects (see Methods). The bivariate outcomes modeled by these parameters were the incidence counts for the ASD and ID phenotypes, tabulated separately for males and females in 3,111 counties, nested within 50 states (plus the District of Columbia) and adjusted for population size.
The results are summarized in Figures 1-5 and Table 1: We observed clear spatial clusters for both ASD and ID. The raw data analyzed prior to complex modeling (see Figure 1) indicated that putative environmental variables were strongly predictive of rates of ASD and ID across the US. This trend persisted after the analysis was corrected for confounding variables using mixedeffect Poisson regression, see  We found that ASD in males (normalized by county population) has a county-level mean rate of 0.1% per male of any age. The distribution of rates across counties is skewed: the median is 0.023% while maximum observed value is 5.2%. Similarly, for ID, the average rate over the whole country is 0.024% per any male, 0.025% per any female. The maximum per county per person rate reached 0.9% and 0.58% for males and females, respectively. Furthermore, for males, the rates of ASD and ID at the county level were weakly but significantly correlated (Pearson product-moment correlation 0.0589, p = 0.00101), while for females the correlation was much stronger (0.197, p,2.26610 216 ). The estimated Poisson rates of incidence produced by model-based inference are shown in Figure 2. When direct comparison was possible, we compared our conclusions with those of prior studies ( [9][10][11][12][13][14], see Table 2) and found that they were consistent.
Accumulating evidence [16][17][18][19][20][21][22][23][24][25][26][27][28] suggests that the rate of birth malformations, especially of those affecting the reproductive system in newborn boys [17], adjusted for population size and structure, could serve as an indicator of average parental exposure to toxins within a geographic unit. After controlling for ethnicity, gender, and socioeconomic factors, the strongest predictor of ASD was the rate of male congenital malformations of the reproductive system, used as an approximate measurement for exposure to teratogens, based on extensive epidemiological evidence ( [16,17], see Figure 3 and Discussion). Every additional percent incidence of male congenital malformations of the reproductive system was predictive of a 283% increase in the rate of the ASD incidence (95% confidence interval, CI: [91%, 576%], p,6610 25 ). Similarly, non-reproductive congenital male malformations accounted for a 31.8% ASD rate increase (CI: [12%, 52%], p,6610 25 ). In contrast, male congenital malformations of the reproductive system were barely significantly predictive for ID (94%, CI: [1%, 250%], p = 0.0383). However, the effect of non-reproductive congenital malformations in males on ID incidence was statistically significant and strong: an increase of 43% (CI: [23%, 67%], p,6610 25 ). Another variable significantly affecting both ASD and ID was population-adjusted incidence of viral infections in males (Table 1, Figure 3, and Discussion). Moreover, comorbidity analysis demonstrated that male children with ASD are 5.53 times more likely to have congenital genital malformations than unaffected males (odds ratio 95% CI [5.22, 5.87], p,2.2610 216 , Fisher's exact test).
Male congenital malformations of the reproductive system are subdivided in the ICD9 taxonomy [29], which was used to encode the data in this analysis, into unspecified malformations, hypospadias (abnormally placed external urethral orifice), epispadias (the urethra does not develop into a full tube), micropenis, congenital chordees and undescended testicles. The US average incidence rate for all male congenital malformations of reproductive system was 0.2687% per male of any age group; of these 18% were hypospadias (0.049% rate), 6% congenital chordees (0.0161% rate), 1% micropenis (0.00275% rate), and 0.083% epispadias (0.00227% rate). The rest of the malformations were in the unspecified category; undescended testicles were not encoded explicitly. Per-county rates of hypospadias and congenital chordees were significantly correlated with each other (Pearson's r = 0.34, 95% CI [0.31, 0.372], p,2.2610 216 ), as were hypospadias and epispadias (r = 0.066, 95% CI [0.031, 0.101], p = 0.00022). The highest per county rates of malformations were 2.4% (all male genital malformations), 1.1% (hypospadias), 0.91% (chordees), 1.1% (epispadias), and 0.23% (micropenis). There were counties with no reported malformations (zero apparent rate). All discussed groups of malformations were significantly correlated with autism rates. Female birth malformations of reproductive system, variable ConGenF, showed very similar disease-specific predictive behavior to the congenital male malformations of reproductive system, variable ConGenM. The female malformations were predictive of an increase in both ASD and ID, but the magnitude of the statistical effect associated with this factor was much smaller than ConGenM, although highly correlated.
Both ASD and ID showed significant gender-specific incidence effects, with males affected more frequently than females; this was

Author Summary
Disease clusters are defined as geographically compact areas where a particular disease, such as a cancer, shows a significantly increased rate. It is presently unclear how common such clusters are for neurodevelopmental maladies, such as autism spectrum disorders (ASD) and intellectual disability (ID). In this study, examining data for one third of the whole US population, the authors show that (1) ASD and ID display strong clustering across US counties; (2) counties with high ASD rates also appear to have high ID rates, and (3) the spatial variation of both phenotypes appears to be driven by environmental, and, to a lesser extent, economic incentives at the state level. more extreme for ASD (Table 1 and Figure 3). Using ethnicity variables to account for genetic heterogeneity of the US population, corrected for socioeconomic factors, such as the mean county-specific income, we found the incidence of the two diseases significantly varied across ethnic groups (Table 1), with Pacific Islanders, for example, having significantly lower risk for both diseases. The per capita income of the county was weakly positively correlated with the incidence rates for both diseases: the income variable was associated with 3.2% rate increase per every additional $1,000 of income above the country average for ASD (CI: [2.3%, 4.2%], p,6610 25 ) and a 2.7% rate increase for ID (CI: [1.8%, 3.7%], p,6610 25 ). Other important socioeconomic predictors included the percentage of urban population in a county; a one percent increase in urbanization predicted about a 3% increase in ASD and ID incidence (Table 1). Our analysis also indicated that state-specific laws [30,31] had a large but only marginally significant effect on the incidence rates of ASD and ID (Table 1, Figure 1). The strictest form of diagnostic evaluation (variable Eval, the state-mandated diagnosis of autism or autism spectrum disorders by a pediatrician or clinician for consideration in the special education system) was predictive of a considerable decrease in ASD and ID incidence rates, 98.6% (CI: [28%,

Discussion
By analyzing the spatial incidence patterns of autism and intellectual disability drawn from insurance claims for nearly one third of the total US population, we found strong statistical evidence that environmental factors drive the apparent spatial heterogeneity of both phenotypes while economic incentives and population structure appear to have relatively large albeit weaker effects. The strongest predictors for autism were associated with the environment: congenital malformations of the reproductive system in males (an increase in ASD incidence by 283% for every per cent of increase in the incidence of malformations), nonreproductive congenital malformations (31.8% ASD rate increase), and viral infections in males (19% ASD rate increase). For ID we observed similar but weaker effects: 93% increase of ID rate for every per cent of increase in congenital malformations of the reproductive system in males, 43% increase for per cent of nonreproductive congenital malformations, and 23% for viral infections in males.
We highlight the role of male congenital genitourinary malformations as an approximate measurement of environmental exposure to unmeasured developmental risk factors, including toxins. Some infants are born with congenital malformations with unknown genetic etiology-not explained by known Mendelian mutations or detectable chromosomal aberrations. At least a fraction of such birth defects may be due to parental exposure to environmental insults. The environmental factors implicated so far include pesticides [18,19], environmental lead [20], sex hormone analogs [21,22], medications [23], plasticizers [24], and other synthetic molecules [25]. More generally, the risk of congenital birth defects is statistically linked to parental occupation [26][27][28]. There is a statistically significant increase in birth defects associated with some maternal occupations (janitor, maid, landscaper), and significant decrease associated with others (nonpreschool teacher) [27,28]. It is very likely that the list of environmental factors potentially affecting development of human embryo is large and yet predominantly undocumented; correspondingly, detailed statistics on these factors do not exist.
It is known that some birth malformations are caused by de novo genetic events, such as large copy number variants that have been found to increase the risk for ASD by approximately 400% [32]. Single-gene deletions, for example, involving CHD7 are known to cause CHARGE syndrome [33,34] associated with genital abnormalities and putatively associated with ASD [35]. However, these genetic events may have currently poorly identified environmental triggers, and 70 to 80% of male congenital malformations of the reproductive system have no clear genetic causes [36]. Instead, they appear to be driven by specific environmental insults that were not serious enough to lead to more serious adverse events during pregnancy, such as spontaneous abortion. Therefore, in this study, we used the rate of birth malformations as a surrogate measure for environmental burden.
The hypospadias of the male urethra can arise during early embryonic development, specifically weeks 9-12 (p. 206 in [36]). This window corresponds to the time when cell division and migration takes place in brain development. Furthermore, maternal exposure to estrogen and estrogen analogs in animal models affects both brain and genital development in male progeny (p. 206 in [36]), and small physical malformations appear enriched in autistic children compared to healthy children [37].
Following similar logic, in addition to causing birth defects, environmental toxins, such as pesticides [38,39] can substantially weaken the human immune system, especially in men, which results in more frequent infections. (The rates of female viral infections were highly correlated with male viral infections; these can serve a somewhat weaker fixed effect predictor, data not shown.) This suggests that per capita rate of viral infection, when socioeconomic and other biological factors have been controlled for, may serve as another environmental indicator, although specifics of the causal, biological mechanisms remain unresolved. In our analysis we found that the rate of viral infection in males was significant for both ASD and ID, see Table 1.
Importantly, the effect of state-level regulations involving ASD appeared relatively large in magnitude (over 98% ASD and ID rate decrease) but with a wide confidence interval and inconsistent effects across states, resulting in only marginal significance. Furthermore, our estimates of random effects at the state and county levels, see Figure 2, suggest that additional yet unknown confounder factors exist at both state and county levels, as is evident from the clear state boundaries seen in Figure 2.
As with other statistical analyses (see Table 3), significant associations are not necessarily causal. Identified predictor variables may reflect underlying mechanisms, or may be correlated with unmeasured causal factors. However, we have  included variables, such as the rate of birth malformations and male viral infections, that have well-documented environmental causes. Overall, this increases our confidence in their scientific relevance. Furthermore, because we have controlled for many county-level socioeconomic variables, strong state-specific effects are almost certainly rooted in legal and regulatory differences that exist at this level.
Our results have implications for the ongoing scientific quest for the etiology of neurodevelopmental disorders. We provide evidence that routinely expanding the scope of inquiry to include environmental, demographic and socioeconomic factors, and governmental policies at a broad scale in a unified geospatial framework. It appears that detailed documentation of environmental factors should be recorded and used in genetic analyses of ASDs and failure to do so risks omitting important information about possibly strong confounders.

Ethics statement
Our analysis involved de-identified patient data and was approved by the University of Chicago Institutional Review Board.
Our multi-level, mixed-effects model predicted the incidence of ASD and ID conditional on several individual-level, county-level, and state-level covariates. For the regression analysis described below, we used county level variables to predict disease rate. In the analysis of the comorbidity between congenital malformations and ASD or ID, we used patient-level data.

Data
We used the Truven Health Analytics MarketScan Commercial Claims and Encounters Database to provide geocoded diagnosis In the figures we color-coded the Empirical Bayes estimates of the state-level random effects, separately for ASD and ID. While county-specific random effects are directly comparable within the same state, comparison of these effects across different states is not meaningful, because each state-specific random effect determines the baseline disease rate for each county in the corresponding state, and these baseline rates vary across states. doi:10.1371/journal.pcbi.1003518.g005 Claims were de-identified, that is, all patient-level personal information was redacted, and geocoded at the county level by Truven, and thus, this did not require any additional processing on the authors' part. The MarketScan insurance claims dataset is not a random sample of the USA population. This is because compilation of this dataset required reaching agreements between Truven and numerous individual insurance providers to share data, and the insurance providers inherently had uneven and non-random coverage of geographic areas. It is possible, therefore, that the dataset carries traces of hidden correlations imposed by the data collection method. Furthermore, while the entire USA is well represented in the data, it is possible that coverage across geographic areas is not perfectly proportional to population density.

Statistical analysis
We framed our analysis as a mixed-effect regression model for Poisson-distributed count data [40], independently implemented in SuperMix [41], lme4 [42,43], MCMCglmm [44], and GLLAMM [45,46]. The choice of the Poisson model was motivated by the countable nature of data and the rarity of the disease incidence events.
The model parameters were estimated using a joint statistical inference as follows. Most of the parameters (44 out of 50) were real-valued coefficients representing regression weights of individual factors, such as average income in county, percentage of ethnic groups, per cent of urban and poor population, see Table 1. The factors that are a priori suspected to be relevant to disease incidence and are deliberately included into the model, are referred to in the mixed effect model formalism as the fixed effects. In addition, the model included zero-centered and normally distributed random effects, uncorrelated for the same disease between a county and the encapsulating state, but geographically correlated between two diseases, see equations below and Figures 4 and 5. The six parameters associated with the random effects included the variances and covariances for the state-and county-level random effects (see Table 2 and equations below). Note that the random effects themselves were not parameters, but Gaussian zerocentered random variables.   Fixed effects by design have stochastic but predictable influence on data, while random effects describe zero-centered random influence, not captured by the fixed effects. Below we present more formal formulation of the model.

Fixed-effect variables
The Truven database was augmented with US census data [30] consisting of county-level measurements for a variety of socioeconomic, demographic, and geospatial factors. Our fixedeffect county-level covariates were gender, Gender, average per capita income, Income, percent ethnicity (separately for American Indians, AmInd, Asians, Asian, White Hispanics, WHisp, White non-Hispanics, W, Black Hispanics, BHisp, Black non-Hispanics, B, and Pacific Islanders, Pacific), and the proportions of various socioeconomic groups (poor, Poor, urban, Urban, insured, Insured). Our county-level environmental indicators used as fixed-effect covariates (normalized by county population size) comprised congenital malformations excluding malformations of the genitals (separately for females and males, CongMrepF and CongMrepM, respectively), congenital malformations of the genitals (separately for females and males, ConGenF and ConGenM, respectively), viral infections (separately for females and males, Viral_F and Viral_M, respectively), ectopic pregnancy (ect_pr), abnormal conception (abnormal_concept), spontaneous abortion(spont_abort), and multiple gestations (mult_gest). The county-level environmental indicators were extracted from the Truven database and normalized by county population, separately for males and females.
To account for variation in policies for special education eligibility and reimbursement, we used four variables derived by hand-coding the state policies for eligibility in special education programs under the Individuals with Disabilities Education Act [31] with categorical variables: (i) CFR, to indicate whether state criteria met (21) (21, 22, 1, 2). Details of the coding are described in the abbreviations. All predictor variables were mean-centered.

Assumptions of the model
The assumptions of the Poisson regression model [40] were as follows. First, we assumed that the data, corresponding to the observed counts of people within each county diagnosed with either ASD or ID, were generated by a Poisson process, with rate (l ijkl ) varying over counties, where: h is a vector of all 50 model parameters, (b, S). The observed counts of disease incidence (the response variable y ijkl ) was defined as the number of disease cases per county for ASD (k = 1) and ID (k = 2). Subscripts i and ij were used to indicate a state and a county nested within that state, and subscript l to indicate gender. The second assumption was that the logarithm of Poisson rate (l ijkl ) was expressed as a linear combination of fixed and random effects.
Here matrix X k is the design matrix for the fixed effects associated with disease k; b k is the corresponding vector of unknown regression weights; z k is a design matrix for random effects; v k is the vector of random effects; v ik , v ijk , are i-stateand ij-countyspecific random effects for disease k. The fixed-effect design matrix is simply a matrix of county-specific zero-centered properties, such as the mean income, or proportions of ethnic groups. The design matrix z has a very simple form: entries of 1 for random effects of a given county and corresponding state, and zeros in all other cells. N ijkl is a state-, county-, disease-and gender-specific offset-the total number of people with a specified gender living within a given county. The third assumption was that data was hierarchical: the zerocentered random effects were independently introduced at i-state and ij-county levels. However, the random effects associated with the two diseases (1 and 2) were geographically correlated,  where s 2 ik , and s 2 ijk , are variances of state-level and county-level random effects for disease k and r i , and r ij are correlation coefficients for random effects for diseases 1 and 2 at the state and county levels, correspondingly.
Together these assumptions define the following likelihood for the given ij-county, i-state, l-gender and k-disease: The full log-likelihood is obtained by summing the individual loglikelihoods specific to each y ijkl over all possible indices. While estimates varied slightly across different implementations of the model and estimation approaches, the major trends were identical across all. Here we present the results of the Markov chain Monte Carlo/Empirical Bayes analysis. The estimation methods and starting parameter values varied considerably across the implementations. For example, SuperMix started with finding an analytical solution of the fixed-effect part of the equations and then estimated parameters for the full model involving random effects. The Markov chain Monte Carlo-based GLLAMM started with a random set of parameter guesses and then rather quickly discovered the high-probability area of parameter values.

Confounders
We tested several putative confounding variables, such as county-specific median mother's age at childbirth, and the proportion of county population in the childbearing age. While these variables produced statistically significant fixed-effect coefficients, they did not affect the relationship between the outcome variable and the compound environmental predictor variables as would be expected if these added variables were true confounders.