Gestational exposure to endocrine disrupting chemicals in relation to infant birth weight: a Bayesian analysis of the HOME Study

Background Pregnant women are exposed to a mixture of endocrine disrupting chemicals (EDCs). Gestational EDC exposures may be associated with changes in fetal growth that elevates the risk for poor health later in life, but few studies have examined the health effects of simultaneous exposure to multiple chemicals. This study aimed to examine the association of gestational exposure to five chemical classes of potential EDCs: phthalates and bisphenol A, perfluoroalkyl substances (PFAS), polychlorinated biphenyls (PCBs), polybrominated diphenyl ethers (PBDEs), and organochlorine pesticides (OCPs) with infant birth weight. Methods Using data from the Health Outcomes and Measures of Environment (HOME) Study, we examined 272 pregnant women enrolled between 2003-2006. EDC concentrations were quantified in blood and urine samples collected at 16 and 26 weeks gestation. We used Bayesian Hierarchical Linear Models (BHLM) to examine the associations between newborn birth weight and 53 EDCs, 2 organochlorine pesticides (OPPs) and 2 heavy metals. Results For a 10-fold increase in chemical concentration, the mean differences in birth weights (95% credible intervals (CI)) were 1 g (-20, 23) for phthalates, -11 g (-52, 34) for PFAS, 0.2 g (-9, 10) for PCBs, -4 g (-30, 22) for PBDEs, and 7 g (-25, 40) for OCPs. Conclusion Gestational exposure to phthalates, PFAS, PCBs, PBDEs, OCPs or OPPs had null or small associations with birth weight. Gestational OPP, Pb, and PFAS exposures were most strongly associated with lower birth weight. Electronic supplementary material The online version of this article (10.1186/s12940-017-0332-3) contains supplementary material, which is available to authorized users.

The vulnerability of the developing fetus [15] coupled with the pervasive nature of EDCs has raised concerns about reproductive health effects arising from gestational EDC exposure, such as low birth weight [15]. Low birth weight is a predictor for coronary heart disease [16], diabetes mellitus [17] and infant mortality [17]. Additionally, some previously published studies [14,15] found associations between gestational EDC exposure and birth weight. Understanding the role EDCs may play in affecting birth weight will inform efforts to prevent low birth weight and the development of chronic diseases.
Many investigators have explored the relationship between EDC exposures and birth weight [18]. Lenters et al. examined the relationships between 16 chemicals (6 phthalates, 8 PFAS, 2 PCBs and 1 OCP) and birth weight using Elastic Net Regression analyses; 2 phthalates, 1 PFAS, and 1 OCP were associated with lower birth weight [18]. Using a job exposure matrix, which can differ from exposure measurement, Birks et al. (2016) examined occupational EDC exposures in a meta-analysis of European birth cohorts and found that pregnant women exposed to more than 1 EDC class were more likely to have a low birth weight infant [15]. There is sufficient evidence that increased PFAS (especially PFOA) exposure is associated with low birth weight [18][19][20], but mixed results are reported for other EDCs and birth weight [18,21]. For example, a meta-analysis by Govarts et al. (2016), found an association between PCB 153 and low birth weight [21] whereas Lenters et al. [18] did not find an association. Similar inconsistent associations exist for phthalates and BPA [22][23][24], OCPs [18,25] and PBDEs [26,27].
These inconsistencies may be due to a shortage of statistical tools to explore the effects of multiple chemical exposure mixtures on birth weight. One way to model multi-exposure mixtures is to collapse the exposures into a summary term that is included in regression models for the outcome [28] or use proxy variables [29]. These approaches, however, assign equal weight and importance to all chemicals, irrespective of their predictive value. Multivariable linear regression [5,13,24] with selected predictors [30] is another method, but it can produce unreliable estimates due to the highly correlated nature of some EDCs [31]. LASSO (Least Absolute Shrinkage and Selection Operator) [30] and Elastic Net are more sophisticated multivariable linear modeling methods that have been used for analyzing multiple correlated exposure variables.
Bayesian hierarchical linear models (BHLMs) are one possible solution to analyzing the effect of multiple correlated exposures [32]. BHLM is a statistical method belonging to a class of methods known as shrinkage methods that, by placing constraints on the distribution of the variables in the model, compress the beta coefficients toward a common mean, resulting in more stable estimates and avoiding the collinearity problem. This also reduces the Type I error (false positive) rate because the shrunken beta coefficients are closer to zero [33]. BHLM assumes that in a linear model, individual regression coefficients arise from a normal distribution, called the prior distribution. The prior distribution depends on hyperparameters, including a mean μ β and standard deviation σ β, which are informed by the data and can be quantified when fitted to the data. The quantity μ β is the average association between all regression coefficients and the outcome, whereas σ β captures the heterogeneity of effects.
The objective of this study was to use BHLM to simultaneously examine the association between exposure to 5 different EDC classes (phthalates and BPA, PFAS, PCBs, PBDEs and OCPs), 2 organophosphate pesticides (OPPs) and 2 heavy metals, lead (Pb) and mercury (Hg), with birth weight among women and their children in the Health Outcomes and Measures of Environment (HOME) Study [34].

HOME Study
The HOME Study is a prospective birth cohort study in the Cincinnati, Ohio area. Women were recruited between 13 and 19 weeks of pregnancy from prenatal clinics between 2003 and 2006. Inclusion criteria for the HOME Study included women who were >18 years old, <19 weeks gestation at the time of enrollment, and living in a residence built before 1987. Additional details about study eligibility, recruitment, and follow-up are described elsewhere [34]. The Cincinnati Children's Hospital Medical Center (CCHMC) and participating delivery hospitals' Institutional Review Boards (IRB) approved the HOME Study. After research assistants explained the study protocols, women provided consent for themselves and their children. Demographic information was collected via standardized questions administered by trained research assistants as described in more detail elsewhere [34] (Table 1).

Outcome variables
The outcome variable was birth weight, measured in grams (g). Birth weight was abstracted from birth records and analyzed as a continuous variable.

Biomarkers of gestational EDC exposure
Blood and urine samples were collected from participants at approximately 16 weeks and 26 weeks gestation. We used sensitive and specific liquid or gas chromatography mass spectrometry [34] to measure gestational EDC exposure. A total of nine phthalates, BPA, 5 PFAS, 23 PCBs, 9 PBDEs, 6 OCPs, 2 OPPs and 2 heavy metals were measured ( Table 2). EDCs with measurements below the limit of detection (LOD) were assigned LOD/√2 [10,35,36]. Phthalate metabolites and BPA were creatinine-normalized in μg/g creatinine and measured twice. The metals were also measured twice while remaining biomarkers were only measured once. PFAS were expressed in units of μg/L serum. PCBs, PBDEs, OCPs and OPPs were lipidnormalized in units of ng/g serum lipid. All biomarker samples were log 10 transformed and treated as continuous variables. If a woman provided multiple urine or blood samples at both 16 and 26 weeks gestation, then the sample concentrations were averaged following log 10 transformation and treated as continuous variables [34,37]. Biomarkers with <90% detection level or that were only examined in a subset of the participants were excluded from analysis to preserve sample size ( Table 2). Building on a previous analysis of the HOME Study data [14], we examined OPPs as the molar sum of diethylphosphate (DEP) and dimethylphosphate (DMP) metabolite concentrations. Pb, while not conventionally considered an EDC, is a widespread reproductive [38] toxin that impacts neurodevelopment [6] and birth weight at very low exposure levels [39] and so is included in this study.

Potential Confounders
We drew a directed acyclic graph to select confounders based on the relationships between the HOME Study   [10,34] covariates and EDCs and birth weight (Additional file 1: Figure S1). Covariates in the birth weight BHLM included: maternal race (white vs. black/other), age at delivery, infant sex, maternal education (<high school vs. ≥high school), tobacco exposure (<3ng/ml vs. ≥3ng/ml serum cotinine concentrations), household annual income (<$25,000 per annum vs. >$25,000 per annum), employment, maternal insurance status (public/none vs. private), marital status (single vs. living with a spouse or partner), pre-natal vitamin use (yes vs. no), and maternal BMI (underweight, normal, overweight, obese). To account for the effect of gestation duration on birth weight, we controlled for gestational age (GA) non-parametrically using the 'splines' package in R [40].

Statistical Analyses
We used descriptive analyses to examine the distribution of EDC variables, birth weight, and participant characteristics. We calculated basic descriptive statistics for participant variables (Table 1), EDC concentration geometric means and percentiles (Table 2). We also examined pair wise EDC correlation in a heat map (Additional file 1: Figure S2). In multivariable modeling, each EDC exposure was standardized (mean = 0, standard deviation = 1) and BHLM was fitted using Bayesian software Stan and the 'rstan' package in R [40,41]. The final adjusted model was: Y = β 0 + βx 1 X 1 ....βx p X p +βc 1 C 1 +....+βc k C k +spline(GA)+ error, where Y is the outcome birth weight, β 0 is the y-intercept, βx 1 ...βx p are the exposure variable regression coefficients, X 1 ...X p are the exposure variables, spline(GA) is the nonparametric spline for GA, and C 1 …C k are the 9 covariates. All 53 exposure variables from 5 chemical classes were included in the final BHLM model for analysis. BHLM assumes that βx i are Gaussian random effects with group specific mean μ β and standard deviation σ β . The hyperparameter μ β is the average association between the exposure variables within a chemical class and the outcome, whereas σ β governs the heterogeneity of effects for the class regression coefficients. The hyperparameters for the different chemical classes were assigned uninformative prior distributions: μ β normal with mean 0 and variance 10 5 , and σ β normal with mean 0 with variance 10 5 such that σ β >0. Hamilton Monte Carlo sampling was used to obtain 2 chains of 20,000 iterations with 500 warm-up iterations. The 2 OPPs, Pb and Hg were included in the model as covariates without hierarchical modeling. The final results included both 95% and 50% equitailed posterior credible intervals (95%, 50% CI) for parameter estimates. The 50% CI is the smallest interval that contains 50% posterior probability of containing the true parameter. Using 50% CIs gives higher power to detect small effects in the data and they can serve as a screening tool for identifying weak effects that may warrant further study [33]. Reporting the 50% CIs also decreases the Type II error resulting from BHLM shrinkage. However, a limitation of 50% CIs when used for multiple frequentist hypothesis tests is that they have a higher Type I error rate than using alpha level = 0.05.
Additionally, to better understand the impact of the BHLM shrinkage on the regression coefficients, we reanalyzed the data using both LASSO and Elastic Net. We used the R package 'glmnet' [42], and included the BHLM confounders (Table 3) in the regression as unpenalized variables. Cross-validation was used to determine the penalization parameter. Note that confidence intervals for beta coefficients from LASSO and Elastic Net are not available using standard software, and calculating standard errors for parameter estimates is challenging because the shrinkage estimation produces biased estimates for large coefficients. Following the approach of Lenters et al. [18], the predictors selected by LASSO and Elastic Net were incorporated into an omnibus multiple linear regression model to obtain unpenalized, mutually adjusted estimates.

Results
Of 389 singleton live births, 272 (69.9%) had all covariates and were included in analysis. The study participants were mostly white (61%), married (64%), and university educated (50%) ( Table 1). Several maternal demographic characteristics predicted birth weight, including age, race, marital status, annual household income, public or no insurance and other SES variables. Tobacco smoke exposure also predicted lower birth weight with a 210 g reduction among smoking mothers (>3 ng/ml serum cotinine) compared with non-smoking mothers (<3 ng/ml serum cotinine).
With few exceptions (e.g., PFOA), median concentrations among HOME Study participants tended to be similar to concentrations observed in US women ( Table 2). The GM of phthalate and BPA concentrations ranged from 2.4 ng/ml (MCPP) to 41.3 ng/ml (MECPP). The PFAS exhibited a minimum GM of 1 μg/L (PFNA) and maximum GM of 13.6 μg/L (PFOS). PCBs ranged from 0.06 ng/ml (PCB 49) to 11.4 ng/ml (PCB 153), PBDEs held a minimum GM 0.06 ng/ml (PBDE 66) and a maximum concentration of 19.1 ng/ml (PBDE 47). The lowest OCP GM was DDT (1.1 ng/ml), and DDE had the highest GM (74 ng/ml). Individual chemicals within a class tended to be positively correlated with each other, particularly for persistent compounds like PCBs and PBDEs (Additional file 1: Figure S2 [43]). Phthalates and BPA were less strongly correlated with other chemicals and, with the exception of the di-2-ethylhexylphthalate metabolites, weakly correlated with each other.
The phthalates and BPA, PFAS, PCBs, PBDEs, OCPs, OPPs, Pb and Hg were examined using BLHM ( Table 3). The μ β posterior mean for the phthalates and BPA was 1 g with 95% CI (-20, 23), meaning the average association between the 9 phthalates and BPA and birth weight was 1 g. Another way to interpret this is that for every 10fold increase in exposure to phthalates and BPA as a whole, there was a 1 g increase in birth weight. This increase was not statistically significant at the 0.05 level based on a frequentist hypothesis test using the posterior credible interval. The σ β posterior mean was 10 g, which indicated heterogeneity in the phthalate and BPA effects around this estimate, and that roughly two thirds of the beta coefficients were within 10 g (one standard deviation) of μ β . In contrast, for PFAS exposure, the μ β and σ β posterior means were larger in magnitude and equal to -11 g, 95% CI (-52, 34) and 24 g, respectively. Although the 95% CI for μ β for PFAS did cover zero, the estimate was significantly lower than zero based on the upper limit of the 50% CI (-21, -2). For PCB exposure, the μ β posterior mean was 0.2 g, 95% CI (-9, 10). The PBDE μ β estimate was -4 g, 95% CI (-31, 22). Finally, OCP exposure was associated with a μ β of 7 g, 95% CI (-25, 99) with an estimated σ β of 18 g.
We examined associations between individual EDCs and birth weight to understand which individual EDCs within a class drive the value of μ β , given by the posterior distribution of the BHLM regression coefficients βx 1 ... βx p (see Figs. 1 and 2). In the phthalates and BPA class, there was a 9 g birth weight increase in response to a 10-fold increase exposure to MEOHP, and a 10-fold increase in MEHP exposure was associated with a 7 g decrease in birth weight. None of the individual phthalate metabolites or BPA were statistically significant at the 0.05 or 0.5 level, which may explain why the class as a whole was not significantly different from zero. In contrast, all 5 PFAS were associated with reductions in birth weight, although none of the 95% CIs excluded zero. PFHXS exhibited the largest association in the BHLM with a 10-fold exposure increase associated with a 17 g birth weight decrease, whereas PFNA was the weakest predictor of birth weight. PCBs ranged from an 8 g birth weight increase associated with a 10-fold increase of PCB 105 to a 7 g decrease associated with a 10-fold increase in PCB 199. The large number of 23 PCBs examined in this study meant that the BHLM performed more aggressive shrinkage on the model regression coefficients for this class. Of the 9 PBDEs, a 10-fold increase in PBDE 17 was associated with the greatest birth weight increase (12 g) while PBDE 153 was associated with the greatest decrease (23 g). The largest increase in the OCP class was seen with a 10-fold increase in HCH exposure (5 g), while the largest birth weight decrease was seen in response to a 10-fold increase in Nonachlor exposure (10 g). None of the OCP class 95% CIs excluded zero. The two OPPs and two heavy metals were examined individually as covariates with no hierarchical priors. For every 10-fold increase in DEP exposure, there was a 14 g increase in birth weight 95% CI (-51, 79). DMP exposure was associated with a 57 g reduction in birth weight (95% CI (-123, 10), 50% CI (-79, -34)). Blood lead concentration was associated with a 45 g decrease in birth weight (95% CI (-110, 22), 50% CI (-68, -22)), while mercury exposure was associated with a 5 g, 95%CI (-56, 66) increase in birth weight.
To examine the impact of shrinkage and the choice of prior distributions on the BHLM results, we applied LASSO and Elastic Net to the phthalates and BPA,PFAS, PCBs, PBDEs and OCPs. As with the BHLM, Pb, Hg and OPPs were included as covariates. LASSO selected 7 predictor variables (2 phthalates, 2 PCBs, 2 PBDEs and 1 OCP) out of a possible 53. The LASSO results are given in Additional file 1: Figure S3. The LASSO Abbreviations: SD standard deviation, CI credible interval, SES socioeconomic status, BMI body mass index Significance: Italicized 95% CIs and 50% CIs do not contain zero and therefore are considered significant at p = 0.05 or p = 0.5 respectively regression coefficients tended to be larger in magnitude than BHLM with wider and less precise 95% intervals, which we expected because LASSO does not use shrinkage that borrows information between chemicals within a particular class. Elastic Net selected the same variables as LASSO, as well as an additional 5 PCBs and 1 PBDE. The results are given in Additional file 1: Figure S4. Our Elastic Net analysis identified 1 EDC in common with Lenters et al. [18], the phthalate MEOHP. In Lenters, this EDC was associated with a -0.15 g difference in birth weight, which was attenuated after adjusting for gestational age. In our analysis, MEOHP was associated with a 13 g birth weight increase (Additional file 1: Figure S3). Additionally, Pb was associated with a 58 g decrease in birth weight in the LASSO analysis (Additional file 1: Figure S3) and a 53 g decrease in the Elastic Net analysis (Additional file 1: Figure S4). DEP, one OPP, was associated with a 16 g birth weight increase using LASSO (Additional file 1: Figure S3) and DMP, the other OPP, was associated with a 48 g decrease (Additional file 1: Figure S3). In the Elastic Net analysis, DEP was associated with a 21 g increase and DMP with a 54g decrease in birth weight (Additional file 1: Figure S4). Because LASSO and Elastic Net do not estimate the μ β parameter, this makes Additional file 1: Figures S2 and S3 harder to interpret because they do not estimate the average association between the exposure variables and the outcome within a chemical class. GA is a potential mediating variable on the causal pathway between EDC exposure and birth weight, while sex is a potential modifier. Controlling for mediators may attenuate the total effect of EDCs on birth weight. Accordingly, we conducted sensitivity analyses without controlling for GA (Table 3) and stratified the BHLM results by sex (Additional file 1: Tables S1 and S2). While the results are mostly unchanged, the phthalates and BPA exposures were associated with a decrease in birth weight for male infants compared to the combined analysis (3 g decrease compared with 1 g increase). The decrease in male birth weight associated with Pb exposure was roughly double the associated birth weight decrease in the combined analysis (94 g decrease compared with 45 g decrease). There are some potentially interesting differences in the impact of DEP and DMP on birth weight when the study sample is stratified by infant sex. While DEP was associated with a 14 g increase in birth weight in the combined analysis, it was associated with a 51 g increase in males and a 9 g decrease in females in the sex stratified analysis. DMP in the combined analysis was associated with a 57 g decrease in birth weight; looking at the stratified analysis revealed that this decrease is greater in magnitude for male infants (88 g decrease), and lesser in female infants (27 g decrease). The associations of greater magnitude present in male infants correspond with previously published studies [24].

Discussion
In this study we examined the relationship between gestational exposure to 53 different EDCs and birth weight in a prospective birth cohort of mother-infant pairs. We found that exposure to phthalates and BPA, PFAS, PCBs, PBDEs or OCPs had null or small associations with birth weight. Other investigators have reported similar small or statistically non-significant associations between birth weight and phthalates [37], PCBs [15], PBDEs [13], and OCPs [25]. Several other studies, however, have found associations between increased exposure to phthalates, PCBs, PBDEs or OCPs whereas we did not. A case- control study by Zhang et al. [24] revealed that low birth weight infants in China tended to have much higher meconium phthalate levels than infants who were normal weight. Additionally, a meta-analysis by Govarts et al. [21] reported that increased exposure to PCBs was associated with decreased birth weight. Lignell et al. [27] found a statistically significant inverse association between PBDE exposure and birth weight. Finally, Monteagudo et al. [44] and Lenters et al. [18] both found changes in birth weight associated with increased OCP exposure.
We found some evidence that PFAS, DMP and Pb were associated with small reductions in birth weight. On average, a 10-fold increase in exposure to the PFAS class was associated with an -11 g change in birth weight, 50% CI (-21, -2). However, due to the higher Type I error rate associated with using alpha level = 0.5, this association should be interpreted with caution. Similar findings concerning PFAS exposure have been reported elsewhere [5,19]. While the exact mechanism by which PFAS impact birth weight is unknown, several theories have been proposed. One explanation is that PFAS interact with estrogen receptors and disrupt hormonal balances [20]. Additionally, PFAS alter human [20] and animal [3] serum lipid levels, potentially affecting fetal growth and development. Finally, PFAS affect adipose tissue development and the regulatory systems involved in body weight homeostasis, which may impact fetal growth outcomes [20], including those that occur after birth such as infant mortality or cardiovascular disease which are beyond the scope of this study [16,17]. A 10-fold increase in Pb exposure was associated with a 45 g birth weight reduction (50% CI (-68, -22), while increased exposure to DMP was associated with a -57 g (50% CI(-79, -34)) change in birth weight. Both of these findings, as with the PFAS findings, should be interpreted cautiously because alpha = 0.5. They do, however, agree with previously published studies [14,39,45,46].
The BHLM approach has advantages compared with other analytic techniques including LASSO and Elastic Net. First, each EDC class contains between 4 to 23 chemicals with varying associations between any given chemical and birth weight. Having many exposure variables complicates the interpretation of associations between the 5 EDC classes and birth weight when using traditional multi-variable linear methods. However, BHLM examines the average effect of EDC exposure on birth weight at the class level using the parameter μ β , which is less susceptible to random variation [32]. A second advantage of reporting μ β is that it prevents the investigator from only reporting chemicals that are significantly related to birth weight by virtue of it being an average association between a class and the outcome. For example, rather than sorting through 23 PCBs to identify those that are significant, we simply reported the average association of all 23 PCBs and birth weight. Additionally, Bayesian modeling produces the σ β value that measures the heterogeneity of EDC class effects. A final advantage of Bayesian modeling is that posterior distributions of the regression coefficients are shrunk towards pooled mean μ β , reducing the Type I error (false positives) rate when the interval estimates are used for frequentist hypothesis testing.
There were several limitations to the study. First, the BHLM model did not account for interactions between classes, and we assumed that all chemicals within a class share a similar mode of action and arise from a common prior distribution. Second, the model also excluded nonlinear relationships, including potential threshold effects that could yield important information about the relationships between EDC exposure and birth weight. Third, BHLM results are sensitive to the assigned prior distributions of hyperparameters μ β and σ β . We used priors following recommendations of a recent study of vague priors for random effects models [47]. BHLM, because it is a shrinkage method, is also prone to higher Type II errors. We included the 50% CIs in our study to account for this. Fourth, our study includes both persistent and non-persistent EDCs when examining EDC exposure. Non-persistent chemicals exhibit larger degrees of measurement error. For example, phthalates metabolites and BPA have short half-lives and this may have attenuated regression estimates and reduced statistical power [48]. To reduce exposure misclassification, we pooled the exposures at 16 and 26 weeks [40]. Fifth, we adjusted for gestational age non-parametrically. Other methods to adjust for gestational age include birth weight z-scores [49] and semi-parametric modeling [50]. Lastly, we conducted these analyses under the assumption that structurally similar chemicals would have similar modes of action and therefore should be grouped together. Other grouping schemes, however, may be more informative.

Conclusion
This paper introduces BHLM as a method for tackling multiple correlated exposures in relation to health outcomes and compared this method with LASSO and Elastic Net. LASSO and Elastic Net regression coefficients were larger in size but with less precise 95% intervals compared to BHLM. BHLM uses shrinkage estimation to bias the regression coefficient estimates towards the group mean parameter. Statistical theory demonstrates that the estimators have smaller mean squared error than ordinary least squares estimators [51]. BHLM also has the advantage that it yields estimates of the hyperparameters that aid in interpreting the effects of multiple correlated exposures. We showed that gestational exposure to phthalates and BPA, PCBs, PBDEs or OCPs had small or null associations with differences in birth weight. In contrast, and consistent with previous studies, PFAS, Pb, and OPP exposures were more strongly associated with reduced birth weight [19,39,45,46].