The National Birth Defects Prevention Study (NBDPS) is an ongoing multi-site population based case-control study. The goal of the study is to identify environmental and genetic risk factors for birth defects. Over 35 structural birth defects were ascertained, including gastroschisis. The study included 10 study centers; however, for the present analysis only data from Arkansas, California, and Utah were used owing to the unavailability of geocoded data from the other centers. Cases were ascertained from the birth defect surveillance systems; prenatally diagnosed and electively terminated cases were included. Live-born infants with no major structural malformations were chosen as controls to represent the underlying population that gave rise to the cases. Controls were selected from birth certificates (Arkansas [2000–2007] and Utah) and birth hospitals (Arkansas [1997–1999], California).(Canfield and others, 2009; Robitaille and others, 2009; Yoon and others, 2001) In Arkansas and California, the present analysis was comprised of cases and controls born from 1997 through 2007. Utah joined the NBDPS in 2003; therefore, cases and controls were selected among births from 2003 through 2007. The NBDPS catchment area for Arkansas and Utah include the entire state. In California the catchment area is located in the San Joaquin Valley and is comprised of Fresno, Kern, Kings, Madera, Merced, San Joaquin, Stanislaus, and Tulare counties; these counties represent 13% of the birth population in California.

Mothers of case and control infants were contacted and interviewed in English or Spanish within 6 weeks to two years after the estimated due date (EDD) to ascertain demographic data and exposures during pregnancy, including illness history, medication use, occupational history, and residential history. Participation for the interview was 67% for gastroschisis case mothers and 66% for control mothers. Medical records for cases were reviewed by clinical geneticists and cases with known or strongly suspected single gene conditions or chromosome abnormalities were excluded.

In the interview mothers reported their residential addresses from 3 months before pregnancy until delivery. Since multiple residences could be reported, the residence at the midpoint of the first trimester was chosen for this analysis. When reporting addresses, mothers were asked what month and year they started and stopped residing at an address. Since only the month and year were specified, we assigned the first of the month as the start and stop day. If the start or stop month was missing, but the year was available, two assumptions were made in an attempt to fill in the missing month. If the month the mother started residing at an address was missing and the year of residence occurred prior to the pregnancy, then the month of December was assigned to the missing start month. By assigning December to the missing month we were conservatively assuming a mother only lived at that address for one month out of the year. If the date the mother stopped living at an address was missing and the year of residence occurred after the year of her first trimester of pregnancy, then the month of January was assigned to the missing stop month, consequently only assuming 1 month of residence at that address. The addresses were assigned a latitude and longitude (i.e., geocoded) by the Agency for Toxic Substances and Disease Registry. If a mother did not report an address at the mid-point of the first trimester or resided out of state during this time period, she was excluded from the main analysis.

We used the following generalized additive model (GAM) to assess spatial clustering: logit [p(x)] = α + γ’z + S(x1, x2). The left side of the equation is the log of the disease odds, α is the intercept, z is a vector of the covariates, and (S(x1, x2)) is a non-parametric smoothing function, where x1 and x2 are the latitude and longitude of the maternal residence. Without the smoothing function the model reduces to an ordinary logistic regression model. A loess smoother was used, as it adjusts to changes in data density and gives more weight to nearby points and less to those further away. The percentage of data used in the smoothing process was determined by the span size. For instance, when a span size of 0.20 was used, the smoothing process would use 20% of the data around the location. In general, the use of a small span size results in increased variability (e.g., the detection of random patterns) but reduced bias, while the use of a large span size produces a smoother surfaces resulting in less variability but increased bias. To determine the optimal span size and balance the tradeoff between bias and variance, various span sizes were tested and the value that minimized the Akaike’s Information Criterion was chosen as the optimal amount of smoothing.

To create the final maps, a grid was placed over the study area and grid points that fell outside the area were removed. Additionally, areas where people cannot live (e.g., large lakes) or with sparse populations were clipped from the grid. The log odds were predicted at each grid point using the GAM. In the adjusted models, the covariates were held constant and the adjusted log odds were predicted at each point. Using the entire study population as the reference, the log odds were converted to odds ratios (ORs) by dividing the log odds at each grid point by the log odds from the model without the smoothing term. The smoothing term serves as a measure of location; therefore, omitting the smoothing term results in computing the log odds for the entire study population. All modeling was performed with the MapGAM package in R and the maps were visualized using ArcGIS.(ESRI, 2011; R Development Core Team, 2009)

We conducted a global test with a null hypothesis that case status was not dependent on location. The global test compared the deviance from the model with the smoothing term to the model without the smoothing term; since the smoothing term serves as a measure of location, the comparison with and without the smoothing term serves to test the significance of location. We compared the deviance statistic to a distribution of the statistics generated under the null hypothesis. The distribution was obtained by randomly assigning a new location to each mother, under the null hypothesis that case status was not associated with location. The dataset was permutated in this manner 999 times and each time the models were re-run and the deviance statistic was calculated to obtain a distribution of deviance. A p-value cut-off of 0.05 was used to identify significant associations.

If the global test indicated that location was significant, a local test was conducted to identify areas of significantly increased or decreased odds on the map. A distribution of log odds at each grid point was calculated using the permutated data from the global test. To determine if our results were due to chance, we compared the log odds at each grid point to the distribution of log odds. All points from the main analysis that fell in the upper or lower 2.5% of the distribution were considered statistically significant.

The following covariates from the interview were assessed for confounding: maternal age, maternal race/ethnicity, maternal education, income, body mass index, first trimester alcohol use, first trimester smoking, illicit drug use anytime during pregnancy, gestational diabetes, nativity, season of conception, and total number of residences during pregnancy. The covariates were added to the model and the maps were visually inspected. If the variable changed the appearance of the map surface, it was included in the final model. In addition, if the optimal span size changed with the addition of a covariate, the variable was also considered for inclusion in the final model.

We conducted a sensitivity analysis to assess if the results of the main analysis changed when excluded mothers were added to the analysis. To identify a residence for mothers missing an address at the mid-point of the first trimester, all reported addresses were examined and the address closest in time to the mid-point of the first trimester was chosen for the sensitivity analysis.

A total of 198 (16.4%) mothers were excluded from the analysis due to not reporting an address at the mid-point of the first trimester (27 cases and 141 controls) or residing out of state (7 cases and 23 controls), resulting in a total of 100 cases and 907 controls used in the spatial analysis. The cases and controls represented births from January 1998 through December 2007. In Arkansas the NBDPS catchment area encompasses the entire state and the spatial distribution of mothers’ residences can be seen in Figure 1A. To preserve confidentiality the location of the maternal residences were altered in the point data maps but in all analyses the exact locations were used.

The spatial analysis revealed elevated crude ORs in the southwest corner of the state, with the highest ORs equaling 2.2 (p-value: 0.48; span size: 0.95) (Figure 1B). After adjustment for maternal age and race/ethnicity, predicted ORs for a non-Hispanic white woman ≥25 years of age in the southwest corner were attenuated to 2.0 and the global deviance remained non-statistically significant (p-value: 0.57; span size: 0.95) (Figure 1C). For the sensitivity analysis, alternate addresses were identified for 45 mothers that were missing addresses at the midpoint of the first trimester. These mothers were added to the sensitivity analysis, along with the 30 mothers that resided out of state. The addition of the 75 mothers resulted in an attenuation of the predicted ORs in the southwest corner of the state (maximum OR: 1.2; p-value: 0.07; span size: 0.95) (Figure 1D).

In California, 158 (12.7%) mothers were excluded from the analysis due to not reporting an address at the mid-point of the first trimester (20 cases and 112 controls) or residing out of state (2 cases and 24 controls), resulting in a total of 193 cases and 892 controls. Due to the geography of this area and sparse population, the prediction grid was clipped to exclude the mountains in the east and western region of the catchment area. The distribution of mothers’ residences can be seen in Figure 2A and the counties included in the catchment area are outlined. To determine the optimal span size for the models, only data from the 8 counties were used. Once the optimal span size was determined, all the data were used in the final models (i.e. California mothers that resided outside the catchment area in the first trimester were included in the modeling).

The results of the crude spatial analysis indicated a slight elevation in risk along the eastern and western borders of the study area (Figure 2B), with the maximum predicted ORs at 1.3 (p-value: 0.27; span-size: 0.95). After adjustment for maternal age and race/ethnicity the maximum predicted ORs did not change (maximum OR: 1.3) and the pattern of risk did not change substantially (p-value: 0.34; span-size: 0.85) (Figure 2C). A total of 83 mothers were included in the sensitivity analysis (57 with alternate addresses and 26 out of state mothers). Seven of the out of state mothers resided on the east coast at the mid-point of the first trimester, due to the large distance between the modeling area and these addresses, the model would not converge; therefore, for these seven mothers in-state addresses were chosen from their residential history. The predicted ORs were attenuated on the western portion of the study area in the sensitivity analysis, however the maximum predicted OR did not change and remained at 1.3 (p-value: 0.36; span-size: 0.85) (Figure 2D).

In this center, 64 (9.0%) mothers were excluded from the analysis due to not reporting an address at the mid-point of the first trimester (6 cases and 35 controls) or residing out of state (2 cases and 21 controls). A total of 78 cases and 562 controls were used in the spatial analysis.

The NBDPS catchment area in Utah included the entire state; however, the spatial analysis was restricted to the geographic area around Salt Lake Valley where the majority of the population resided (Figure 3A). The grid was clipped to exclude the lakes around Salt Lake Valley as well as the southern part of the state and the mountains in the east due to sparse population. Removing these points ensured the model was not predicting ORs for locations where no one resided. Data from the 5 counties around Salt Lake City that contained the clipped grid were utilized to determine the optimal span size for the final models. When running the final models, all the data in the state were used to predict the ORs at each grid point. Prior to adjustment the maximum predicted OR was 2.5 (p-value: 0.61; span-size: 0.30) and after adjustment for maternal age and race/ethnicity there was a slight decrease in the maximum predicted OR for a non-Hispanic white woman ≥25 years of age (maximum OR: 2.4; p-value: 0.34; span-size: 0.45) (Figures 3B and 3C). For the sensitivity analysis, 38 mothers were added to the analysis. Due to issues with the model converging, out-of-state addresses were not used and instead an in-state address was chosen when available. The appearance of the map did not change substantially with the addition of the excluded mothers (maximum OR: 2.7; p-value: 0.20; span size: 0.45) (Figure 3D).