Effects of gestational exposures to chemical mixtures on birth weight using Bayesian factor analysis in the Health Outcome and Measures of Environment (HOME) Study.

Studying the effects of gestational exposures to chemical mixtures on infant birth weight is inconclusive due to several challenges. One of the challenges is which statistical methods to rely on. Bayesian factor analysis (BFA), which has not been utilized for chemical mixtures, has advantages in variance reduction and model interpretation.
METHODS: We analyzed data from a cohort of 384 pregnant women and their newborns using urinary biomarkers of phthalates, phenols, and organophosphate pesticides (OPs) and serum biomarkers of polychlorinated biphenyls (PCBs), polybrominated diphenyl ethers (PBDEs), perfluoroalkyl substances (PFAS), and organochlorine pesticides (OCPs). We examined the association between exposure to chemical mixtures and birth weight using BFA and compared with multiple linear regression (MLR) and Bayesian kernel regression models (BKMR).
RESULTS: For BFA, a 10-fold increase in the concentrations of PCB and PFAS mixtures was associated with an 81 g (95% confidence intervals [CI] = -132 to -31 g) and 57 g (95% CI = -105 to -10 g) reduction in birth weight, respectively. BKMR results confirmed the direction of effect. However, the 95% credible intervals all contained the null. For single-pollutant MLR, a 10-fold increases in the concentrations of multiple chemicals were associated with reduced birth weight, yet the 95% CI all contained the null. Variance inflation from MLR was apparent for models that adjusted for copollutants, resulting in less precise confidence intervals.
CONCLUSION: We demonstrated the merits of BFA on mixture analysis in terms of precision and interpretation compared with MLR and BKMR. We also identified the association between exposure to PCBs and PFAS and lower birth weight.

What this study adds

In the field of environmental epidemiology, analysis and interpretation of mixture effect are statistically challenging due to collinearity among chemicals. We introduced a novel BFA approach and demonstrated that the precision of the estimates for the mixture-effect is improved in BFA compared with BKMR and MLR. BFA also provided the specific contributions of individual chemicals on the mixture, which demonstrated the advantage in providing a clear and easy-to-understand interpretation of the effect estimates.

Introduction

Exposure to chemicals mixtures during pregnancy has been associated with perinatal complications and adverse fetal development, such as preterm birth and low birth weight.[1-14] Most epidemiologic studies, however, are informed by single-pollutant statistical models, particularly linear and logistic regression models, that do not capture the complex exposure profiles in real-life scenarios among pregnant mothers.[2,15] As researchers move beyond the “one chemical at a time” analysis to evaluate mixture effects,[2,16-18] several challenges related to collinearity among individual chemicals and providing easily interpretable analysis results have arisen.[16,19,20]

To combat the challenge of collinearity, several frequentist approaches such as least absolute shrinkage and selection operator and principal component analysis were developed to reduce collinearity by discarding correlated variables that are less impactful.[16,19] In the context of environmental epidemiology, it is difficult to justify the discarding of chemical variables because variables within a class of chemical mixture often share similar biologic pathways. Bayesian methods,[21] on the other hand, are gaining attention in the field of environmental epidemiology as an approach to address the challenges without discarding variables. They provide a more explicit quantification of uncertainty than conventional measures, such as p values, by modeling parameters as probability distributions.[22] Moreover, Bayesian methods have the ability to improve the precision of parameter estimates in the presence of collinearity among variables in mixtures compared with traditional methods.[10,15,23,24] This is typically achieved by combining Bayesian techniques with regularization, shrinkage, and prior information about model parameters. However, Bayesian procedures tend to increase the computation time and complexity of the analysis.

Factor analysis modeling, also known as latent variable modeling, is widely used in the field of psychology to manage collinearity for characteristics that are difficult to directly measure[25-27] and can be extended to a Bayesian framework. Using latent constructs to linearly quantify the combined effects of an unmeasured variable, like a chemical mixture, is an appealing way to address collinearity challenge, although providing interpretable estimates. However, apart from a publication by Ferrari and Dunson, Bayesian factor analysis (BFA) has not been used in this context.[28] Accordingly, we aimed to illustrate the potential benefits of BFA to estimate the association between chemical mixtures and birth weight using data from the Health Outcomes and Measure of the Environment (HOME) Study, a birth cohort from Cincinnati, Ohio, established to study the health impact of various chemical and their mixtures.[29] We also compared our BFA results with two established methods, multiple linear regression (MLR) and Bayesian Kernel Machine Regression (BKMR) to assess collinearity reductions and interpretability.

Methods

Health outcomes and measures of the environment study

The HOME Study is a prospective birth cohort of pregnant mothers and their infants established in 2003 at the Cincinnati Children’s Environmental Health Center, Ohio.[29] The primary goal of the HOME Study is to examine the impact of environmental toxicants on child health. Pregnant mothers who were >18 years old and at 16 ± 3 weeks of gestation and living in a residence built before 1978 were recruited from seven prenatal clinics and hospitals.[29] Out of the 468 women initially enrolled in the study, we excluded 67 women who dropped out before delivery, three stillbirths, nine sets of twins, and five participants missing covariate data. Therefore, 384 mothers who delivered singleton live births, provided biologic samples and had complete sociodemographic information were included in our analysis.

Biomarkers of environmental chemical mixtures

We collected blood and urine samples from participants at approximately 16- and 26-weeks gestation.[29] The Centers for Disease Control and Prevention Environmental Health Laboratories used gas and liquid chromatography-mass spectrometry to measure the concentrations of environmental chemical biomarkers in serum and urine samples as previously described.[29]

With the specific goal of estimating the effect of exposure to environmental chemical mixtures on infant birth weight, we consulted existing literature on birth outcomes to identify potential environmental chemical mixtures to investigate.[10,13,14,30-32] A total of seven classes of chemical mixtures were identified: polychlorinated biphenyls (PCBs), polybrominated diphenyl ethers (PBDEs), phthalates, organochlorine pesticides (OCPs), organophosphate pesticides (OPs), phenols, and perfluoroalkyl and polyfluoroalkyl substances (PFAS). PCBs, PBDEs, and OCPs are lipophilic and were lipid standardized. Phthalate metabolites, phenols, and OPs were creatinine standardized to account for urine dilution. In addition, to preserve the sample size of our analysis, we further restricted our analysis to biomarkers that are widely detected in the population (>80% detected above the limit of detection). Furthermore, to keep the modeling approaches consistent between all chemical classes and mitigate issues with the excessive dimensionality[33] in regression analysis, we selected a total of 35 biomarkers to be included in our final analysis (Table 1). For PCBs, PBDEs, OCPs, and PFAS, we used samples measured at 16 weeks to maintain consistency across measures. For phthalates, OPs, and phenol biomarkers, we averaged concentrations in samples collected at 16 and 26 weeks to represent the overall concentrations. For all biomarkers, measurements below the limit of detection were replaced using single imputation according to Lubin et al.[34] The imputed values were sampled from a truncated lognormal distribution with the mean and standard deviation of the concentration of the chemical variables. The detection limit of each specific chemical was set as the upper bound value for imputation. The concentrations of these biomarkers were log10 transformed to reduce the effects of right skewness in the distribution and to assist with the interpretation of the results. The regression coefficients are interpreted as the change in birth weight for every 10-fold increase in the chemical concentrations.

Table 1.

Names and abbreviation of environmental chemical mixtures and the associated individual chemical biomarkers from pregnant women for HOME study, 2003–2006, Cincinnati, OH, n = 384

Mixture group	Individual chemical biomarkers
PCBs	PCB 118PCB 138PCB 153PCB 170PCB 180
PBDEs	PBDE 28PBDE 47PBDE 99PBDE 100PBDE 153
OCPs	DDEDDTT_NONAOXYCHLORHCB
OPs	DMDTPDETPDEPDMTPDMPDEDTP
Phthalates	ΣDEHP*MBzPMnBPMiBPMEP
Phenols	BPAMPBBP3PPBTCS
PFAS	PFHXSPFNAPFOAPFOS

*Weighted molar sum of the DEHP metabolites calculated from: Mono-(2-ethylhexyl) phthalate (MEHP), Mono(2-ethyl-5-oxohexyl) phthalate (MEHHP), Mono(2-ethyl-5-oxohexyl) phthalate (MEOHP), Mono-(2-ethyl-5-carboxypentyl) phthalate (MECPP), expressed in units of ng/mL of MECPP (308 g/mol).

BP3indicates benzophenone-3; DDT, dichlorodiphenyltrichloroethane; DETP, diethylthiophosphate; DMDTP, dimethyldithiophosphate; DMP, dimethyl phosphate; DMTP, dimethyl thiophosphate; HCB, hexachlorobenzene; MBzP, mono-benzyl phthalate; MiBP, mono-iso-butyl phthalate; MnBP, mono-n-butyl phthalate; OXYCHLOR, oxychlordane; T_NONA, trans-nonachlor.

Outcome variable

Infant birth weight, measured in grams (g), was abstracted from the medical records and examined as a continuous variable. To examine fetal growth, we adjusted for gestational age and measured in weeks.[14] An alternate way would be to use gestational age-specific birth weight z-scores. However, its interpretation is not straightforward as it reports effect measures in the unit of standard deviation, which results in different absolute amounts of weight across the gestational age spectrum.[35]

Covariates

A direct acyclic graph was drawn to select confounders based on the relationship among potential covariates, the selected five classes of environmental chemical mixtures and birth weight (eFigure 1; http://links.lww.com/EE/A141). Exposures to lead and tobacco smoke have been documented to have effects on infant birth weight.[30,36] Therefore, we included the biomarker measurements of lead and cotinine as covariates. Additional covariates in the statistical models included maternal age at delivery, infant sex, race, marital status, maternal education, maternal BMI, and annual household income. We excluded maternal BMI from the covariates in the analysis of lipophilic chemicals to avoid duplicate adjustment since the concentrations were already lipid-adjusted, which is directly related to BMI.[37] The effect of gestational duration on birth weight has been documented to be nonlinear,[14] therefore, we used the cubic spline approach for the adjustment of gestational age using the “splines” package in R.[38]

Analytic approach

The primary analytic approach for this study was BFA to estimate the mixture effect of environmental chemicals on birth weight. We also compared our results with BKMR, and additionally, two bridging methods that can be viewed as the intermediate step bridging MLR and BFA. Individual pollutant MLR model was used as a sensitivity analysis because it is the most commonly used method in the literature of environmental epidemiology for continuous outcomes. It also serves as a benchmark for comparison with the more advanced methods.[39] BKMR is a combinatorial method of Bayesian approach and nonlinear approximation methods that is gaining attention in the field of environmental epidemiology.[23,40,41]

Approach 1—Bayesian factor analysis

We used BFA to assess the association between each class of mixture and birth weight, although adjusting for covariates. A total of seven BFA models for the seven specific chemical classes were examined. We did not include a BFA model with all chemical classes simultaneously because dimensionality increases dramatically when more variables are included in the model and result in low statistical power with small sample size. We used the confirmatory factor analysis approach for the purpose of generating regression coefficients that can be interpreted with respect to a specific class of chemical mixture.[42,43] Bayesian techniques were used for regularization and easier interpretation of the parameter estimates. Each class of mixture is represented by a latent variable illustrated by the following set of equations:

for i=1, ..., k, where k is the number of chemicals in the mixture, and where , are the Y-intercepts, is the regression coefficient for the latent variable for each mixture class of chemicals, Z is the latent variable representing each mixture class of chemicals, is the vector of regression coefficients of covariates, C is the vector of confounders such as age and household income, is the ith individual chemical within the class of mixture . is the factor loading score of the ith chemical on mixture , and and are the normally distributed random errors. For parameters and , we assigned uninformative normal priors with variance equal to 1,000. For y intercepts , , and random errors , we used the default priors in R package “blavaan.”[44] The variance of Z was set to be 1.0 for identifiability.[44]

Markov Chain Monte Carlo (MCMC) sampling was accomplished with the R package “blavaan”[44] and “rstan”[45] to generate samples from a posterior to estimate parameters of the interest. For each BFA model, the number of iterations were determined experimentally to achieve convergence assessed by the measure of the potential scale reduction factor. As a result, a total of 40,000 iterations were run for samples with 2,000 burn-in iterations.

Approach 2—Bridging methods between BFA and MLR

To examine the mathematical relationship between MLR and BFA, we included two additional bridging methods to obtain regression estimates that can be viewed as the intermediate steps bridging MLR and BFA. As previously shown in Equations 1 and 2, BFA can be conceptually broken down into three hierarchical steps. The first step estimates the latent variable Z for each study participant, which is denoted by the symbol . The second step computes the parameter estimates for the effect of the latent variable. The third step applies Bayesian prior distributions on the parameter estimates. Therefore, the first bridging method was factor analysis (FA) using the R package[44] “lavaan.” The second bridging method was “MLR with extracted factor score,” which is MLR incorporating the estimated latent variable Z obtained from the FA model in the absence of individual chemicals, although adjusted for covariates as the following equation:

where is the Y intercept, is the estimated factor score, is the regression coefficient for the estimated factor score, is the vector of regression coefficients of covariates, C is the vector of confounders, and is the random error. By comparing the different regression estimates from the sensitivity analyses, we can observe the gradual differences in the precision of estimates.

Approach 3—Bayesian Kernel machine regression

The BKMR model can be described using the following equation:

where is the Y-intercept, is the vector of exposure-response functions for each individual chemical within the specific class of mixture, is the vector of regression coefficients of covariates, C is the vector of confounders, and is the random error in the model. The exposure-response functions for each individual chemical were determined by employing Gaussian kernel functions nonparametrically based on the available data structure.[23] Since the exposure-response functions were determined based on the data, the priors for the parameters of each individual chemical were also specified differently according to the exposure-response functions with details explained by Bobb et al.[46] Similar to BFA, Markov Chain Monte Carlo (MCMC) sampling was also used in BKMR to generate samples from posterior distributions for the estimation of the parameters and the dose-response curves illustrated by the cross-section views of the exposure-surface functions. The posterior samples were sampled from a total of 20,000 iterations, which is determined experimentally to achieve convergence assessed by the measure of the potential scale reduction factor. The variable selection feature of “bkmr” R package[46] was not used in our primary analysis since we intended to retain all chemicals in the class in the BKMR model to compare with other methods. A separate analysis of BKMR using the variable selection feature was also used to assess the impact of such feature on the analysis results.

Sensitivity analyses—MLR

We used MLR as our sensitivity analyses to assess the association between each individual chemical and birth weight, although adjusted for covariates. We also used MLR to assess the association between each individual chemical and birth weight, although adjusted for both covariates and copollutants within that class of chemicals.

Results

Descriptive statistics

The study participants consisted of 384 mother-singleton newborn pairs. Due to various degrees of the missingness that could not be imputed (missingness due to incomplete biospecimen collection or insufficient volumes for chemical assays instead of measurements below limit of detection), the sample sizes for our analysis were 360 for OPs mixture, 366 for phthalates mixture, 284 for PBDEs mixtures, 237 for OCPs mixtures, 310 for PCBs mixtures, 296 for phenols mixtures, and 307 for PFAS mixtures. Mothers who participated in the study were mostly White (62.5%), married (65.6%), and had at least a bachelor’s degree (60.1%). The mean infant birth weight was 3,352 g with a standard deviation of 632 g. The infant sex ratio was roughly 1.18 to 1 (54.2% female to 45.8% male). Sociodemographic characteristics that were associated with birth weight included maternal age, household income, and maternal BMI (Table 2). Infant birth weight tended to decrease with increasing maternal age and increased with increasing household income and maternal BMI (eTable 1; http://links.lww.com/EE/A141).

Table 2.

Distribution of birth weight in relation to participant characteristics among women in the HOME study, 2003–2006, Cincinnati, OH

	n (%)	Birth weight (g)mean ± SD
All participants	384 (100%)	3,352 ± 632
Maternal age
<25	89 (23.2%)	3,066 ± 608
25–29	109 (28.0%)	3,447 ± 623
30–34	123 (32.0%)	3,454 ± 579
35+	63 (16.4%)	3,391 ± 674
Education
Bachelor’s degree or higher	233 (60.1%)	3,403 ± 664
Some college or 2-y degree	94 (24.5%)	3,272 ± 573
High school diploma or some high school	57 (14.8%)	3,270 ± 577
Race
White	240 (62.5%)	3,484 ± 635
Black	117 (30.5%)	3,128 ± 538
Other	27 (7.0%)	3,148 ± 680
Marital status
Married, living with partner	252 (65.6%)	3,454 ± 637
Not married, living with partner	53 (13.8%)	3,200 ± 551
Not living with partner	79 (20.6%)	3,128 ± 596
Household income
<$25,000	98 (25.5%)	3,123 ± 515
>$25,000 and <$50,000	83 (21.6%)	3,392 ± 685
>$50,000 and <$100,000	139 (36.2%)	3,472 ± 670
>$100,000	64 (16.7%)	3,390 ± 559
Infant sex
Male	176 (45.8%)	3,473 ± 686
Female	208 (54.2%)	3,249 ± 565
Maternal BMI
Underweight or normal	161 (41.9%)	3,309 ± 582
Overweight	130 (33.9%)	3,381 ± 629
Obese	93 (24.2%)	3,385 ± 718

BMI indicates body mass index; CI, confidence intervals.

A high degree of correlation was detected among chemicals within the same class (Figure 1). For example, all PCB congeners displayed correlation coefficients in the range of 0.51 (PCB 118 and PCB 180) to 0.99 (PCB 170 and PCB 180) with each other.

Figure 1.

Pearson correlation coefficients between environmental chemical biomarkers. The color intensity of shaded circles indicates the magnitude of the correlation. Blue indicates a positive correlation, although red indicates a negative correlation.

BFA analysis results

We ran seven BFA models for the seven different classes of chemical mixtures. The class-specific regression coefficients of each class of the mixture and loading of the individual congeners on the mixture were evaluated by BFA. PCBs and PFAS displayed associations with birth weight reduction and every 10-fold increase in the concentration of the mixture (Figure 2). Specifically, the regression coefficients were −81 g (95% confidence intervals [CI] = −132 to −31 g) for PCBs and −57 g (95% CI = −105 to −10 g) for PFAS.

Figure 2.

The associations between every 10-fold increase of the latent mixture of PCBs (A), the latent mixture of PFAS (B) and birth weight (represented by coefficient β) and the factor loadings of the individual congeners onto the latent mixture (represented by coefficient γ) among mother-child birth pairs in the HOME Study estimated by BFA adjusted for covariates including cubic-spline gestational age, maternal age, maternal education, race, marital status, household income, infant sex.

In addition to the mixture-specific effect estimates, we also estimated the loading coefficient of individual chemicals within the class of chemical mixtures, denoted by the quantity from Equation 2 (Figure 2). This provides a relative measure of importance for these individual chemicals because it measures how much influence each individual chemical variable contributes to the overall latent mixture variable.[43] For the PCB mixture, we observed that PCB 170 and PCB 180 had a stronger impact on the overall latent mixture than PCB 153, PCB 118, and PCB 138. For the PFAS mixture, we observed that perfluorohexanesulfonic acid (PFHXS) and perfluorooctanoic acid (PFOS) had a stronger impact on the overall latent variable compared with perfluorononanoic acid (PFNA) and perfluorooctanoic acid (PFOA). The BFA results of the other chemicals are represented in eFigure 2; http://links.lww.com/EE/A141.

We also observed a slight decrease in birth weight with every 10-fold increase in the concentration of OCPs mixture at −16 g (95% CI = −66 to 34 g) and phenols mixture at −21 g (95% CI = −71 to 28 g). A slight increase in birth weight was associated with every 10-fold increase in the concentration of PBDEs mixture at 25 g (−18 g, 69 g), phthalates mixture at 49 g (95% CI = −3 to 99 g), and OPs mixture at 8 g (95% CI = −40 to 54 g). However, all the credible intervals of OCPs, phenols, PBDEs, phthalates, and OPs were imprecise and contained the null value of zero.

Bridging methods results comparing BFA with MLR

We provide detailed comparisons of regression estimates for PCBs and PFAS and their 95% confidence interval across different methods (Figure 3). We observed that when collinearity between the individual chemicals is present, copollutant adjustments in MLR results in less reliable parameter estimates and poorer precision for both PCBs and PFAS. Small sample size also played an important role. For PCB 153, for example, we observed that for single pollutant MLR, the precision interval of the estimators was −99 g (95% CI = −143 to 147 g). After copollutant adjustments, the precision interval of the estimators inflated to −418 g (95% CI = −1645 to 808 g). When using the bridging methods for PCBs, the precision interval obtained from MLR with extracted factor score was −24 g (95% CI = −89 to 42 g), which represents a range of 131 g for 95% CI. The precision interval obtained from FA was −43 g (95% CI = −110 to 10 g), which represents an absolute range of 120 g for 95% CI. And finally, the precision interval obtained from BFA was −81 g (95% CI: −132 to −31 g), which represents an absolute range of 101 g. For PFAs, the precision interval obtained from MLR with extracted factor score was −58 g (95% CI = −117 to −13 g), which represents a range of 104 g for 95% CI. The precision interval obtained from FA was −57 g (95% CI = −115 to −16 g), which represents an absolute range of 99 g for 95% CI. And finally, the precision interval obtained from BFA was −57 g (95% CI = −105 to −10 g), which represents an absolute range of 95 g.

Figure 3.

Comparisons of the association between PCBs (A), PFAS (B), and change in birth weight (g) according to different methods: The red bars represent the regression estimates β with 95% CI for the single pollutant MLR model adjusted for covariates and copollutants. The green bars represent the regression estimates β with 95% CI for the single pollutant model adjusted for covariates but not copollutants. The blue bars represent regression estimates β with 95% CI for three different mixture-specific models related to the factor analysis model outlined: (1) MLR with the extracted latent variable, (2) FA, and (3) BFA.

Therefore, Figure 3 illustrates that the precision of the regression estimates increased when more conceptual steps of BFA were performed. This demonstrates BFA can provide a more precise measure of mixture effects when multiple correlated copollutants are in the model.

BKMR analysis results

The dose-response function between individual PCBs and PFAS chemicals and the change in birth weight using BKMR are shown (Figure 4), as is the overall association between the chemical mixtures and birth weight (Figure 5). Exposure to PCB congeners and PFAS congeners both displayed inverse associations with birth weight. It is apparent that as the concentration quantile of PCB congeners and PFAS congeners increases, then the mean estimate of birth weight decreases. The confidence intervals of the estimates at the extreme ends were the widest due to the smaller sample sizes. The BKMR results for the rest of the chemicals are in eFigures 3 and 4; http://links.lww.com/EE/A141. All regression estimates obtained from BKMR had 95% interval estimates that were imprecise and contained the null. A separate BKMR analysis using the variable selection feature was also conducted to examine the influence of different variables in the model (eTable 2; http://links.lww.com/EE/A141). According to the results, all individual variables from the seven different classes of chemicals were included in each BKMR model, giving the same results as the BKMR when variable selection is not used. This is reasonable because each class of chemical mixture is selected based on its similar chemical structures. All individual chemicals within the mixtures have high correlation with one other. Therefore, the model selection did not drop any variables to improve the parameter estimates. The posterior inclusion probability for each chemicals are given in eTable 2; http://links.lww.com/EE/A141.

Figure 4.

Dose-response function (95% credible intervals) between every 10-fold increase in concentrations of selected PCB congeners (A) and birth weight while fixing other PCB congener concentrations at median values and PFAS congeners (B) and birth weight while fixing other PFAS congener concentrations at median values estimated by BKMR adjusted for covariates including cubic-spline gestational age, maternal age, maternal education, race, marital status, household income, infant sex.

Figure 5.

Difference in birth weight (95% credible intervals) for different percentiles of the concentrations of all PCB congeners (A) and all PFAS congeners (B) while centering the effect at median concentrations at zero estimated by BKMR adjusted for covariates including cubic-spline gestational age, maternal age, maternal education, race, marital status, household income, infant sex.

Sensitivity analysis results

When the biomarkers were analyzed one at a time in MLR although controlling for covariates, then PCB 170, PCB 180, and PCB 153 from the PCBs mixture, PBDE 153 from the PBDE mixture, dichlorodiphenyldichloroethylene (DDE) from the OCPs mixture, mono-ethyl phthalate (MEP) from the phthalate mixture, diethyl phosphate (DEP), dimethyl phosphate (DMP), and diethyl dithiophosphate (DEDTP) from the OPs mixture, bisphenol A (BPA), methylparaben (MPB), and triclosan (TCS) from the phenol mixture and all biomarkers from the PFAS mixture displayed negative associations with birth weight (Table 3). All the associations, however, contained the null value of zero.

Table 3.

Regression coefficients for the association between individual environmental chemical biomarkers (10-fold increases) and mean birth weight among women in the HOME study, 2003–2006, Cincinnati, OH, using MLR

	β adjusted for covariates (95% CI)	β adjusted for covariates and other chemicals within the mixture class (95% CI)
PCBs (n = 310)
PCB 118	78 g (−143 to 286 g)	3 g (−330 to 336 g)
PCB 138	37 g (−202 to 263 g)	505 g (−235 to 1245 g)
PCB 153	−99 g (−344 to 147 g)	−418 g (−1645 to 808 g)
PCB 170	−118 g (−348 to 111 g)	1267 g (97 to 2437 g)
PCB 180	−194 g (−423 to 35 g)	−1461 g (−2696 to −227 g)
PBDEs (n= 284)
PBDE 28	45 g (−101 to 191 g)	−11 g (−381 to 360 g)
PBDE 47	65 g (−66 to 195 g)	−32 g (−703 to 638 g)
PBDE 99	86 g (−38 to 210 g)	180 g (−186 to 547 g)
PBDE 100	17 g (−105 to 138 g)	24 g (−497 to 545 g)
PBDE 153	−74 g (−189 to 42 g)	−153 g (−392 to 86 g)
OCPs (n= 237)
DDE	−111 g (−358 to 137 g)	−312 g (−627 to 3 g)
DDT	102 g (−66 to 271 g)	179 g (−20 to 377 g)
OXYCHLOR	64 g (−233 to 361 g)	183 g (−522 to 889 g)
HCB	101 g (−284 to 487 g)	132 g (−358 to 622 g)
T_NONA	37 g (−208 to 282 g)	−96 g (−640 to 448 g)
Phthalates (n=366)
ΣDEHP	65 g (−26 to 157 g)	48 g (−55 to 150 g)
MEP	−4 g (−94 to 86 g)	−36 g (−134 to 62 g)
MiBP	66 g (−47 to 181 g)	22 g (−140 to 184 g)
MnBP	73 g (−41 to 187 g)	29 g (−134 to 192 g)
MBzP	62 g (−38 to 162 g)	29 g (−107 to 165 g)
OPs (n=360)
DMDTP	20 g (−26 to 66 g)	34 g (−25 to 94 g)
DETP	38 g (−23 to 98 g)	34 g (−34 to 102 g)
DEP	−29 g (−81 to 24 g)	−19 g (−81 to 43 g)
DMTP	17 g (−58 to 89 g)	35 g (−69 to 140 g)
DMP	−50 g (−106 to 7 g)	−81 g (−157 to −5 g)
DEDTP	−5 g (−54 to 44 g)	8 g (−43 to 59 g)
Phenols (n=296)
BPA	−35 g (−177 to 108 g)	−32 g (−177 to 113 g)
MPB	−29 g (−136 to 78 g)	−69 g (−222 to 84 g)
BP3	40 g (−32 to 111 g)	49 g (−26 to 123 g)
PPB	0 g (−86 to 86 g)	36 g (−86 to 156 g)
TCS	−11 g (−103 to 82 g)	−14 g (−110 to 83 g)
PFAS (n=307)
PFHXS	−109 g (−282 to 63 g)	−41 g (−261 to 179 g)
PFNA	−251 g (−564 to 63 g)	−160 g (−557 to 237 g)
PFOA	−114 g (−339 to 112 g)	22 g (−265 to 310 g)
PFOS	−194 g (−429 to 42 g)	−103 g (−469 to 264 g)

Total sample size for this analysis was reduced to exclude samples with missing values in one or more of the chemical concentrations after the imputation process. The regression coefficients refer to the association with every two-fold increase in the chemical concentration.

Adjusted for all covariates including cubic-spline gestational age, maternal age, maternal education, race, marital status, household income, infant sex).

BMI indicates body mass index; BP3, benzophenone-3; CI, confidence intervals; DDT, dichlorodiphenyltrichloroethane; DETP, diethylthiophosphate; DMDTP, dimethyldithiophosphate; DMP, dimethyl phosphate; DMTP, dimethyl thiophosphate; HCB, hexachlorobenzene; MBzP, mono-benzyl phthalate; MiBP, mono-iso-butyl phthalate; MnBP, mono-n-butyl phthalate; OXYCHLOR, oxychlordane; T_NONA, trans-nonachlor.

The regression coefficients and the variances associated with the regression coefficients were both inflated in magnitude after adjusting for copollutants within the same mixture class (Table 3). Some individual chemicals even showed a reversal in the direction of effect estimates. For example, in the single pollutant model, a 10-fold increase in the concentrations of PCB 170 was associated with a change in birth weight of −118 g (95% CI = −348 to 111 g). In the model adjusted for copollutants, a 10-fold in the concentrations of PCB 170 is associated with a change in birth weight of 1,267 g (95% CI = 97 to 2,437 g). These results show that with the presence of collinearity, MLR is inadequate for mixture analysis because variances associated with parameter estimates would inflate to extreme values, resulting in unreliable and imprecise estimates.

Discussion

Most of the previous mixture analysis of the HOME Study on perinatal outcomes focused on reducing collinearity among individual chemicals. For example, Woods et al[10] used a hierarchical Bayesian approach to reduce collinearity in the data and generated comprehensive estimates of multiple individual chemical congeners. Kalloo et al[32] used both nonparametric (k-means clustering) and parametric approaches (principal component analysis) to generate effect estimates associated with mixtures in conjunction with collinearity reduction. Additionally, numerous papers on the effects of mixtures on other childhood outcomes utilized innovative statistical methods to reduce collinearity.[41,47-53] Yet, few paid attentions to the challenge of interpretability.

We evaluated whether BFA improves precision and interpretability when estimating the health effects of prenatal exposure to chemical mixtures, compared with established methods MLR and BKMR. Among the three methods, BFA produced the most precise effect estimates for the mixture models (Figure 3). The improvement in the precision of the estimate was apparent for both PCBs mixture and PFAs mixtures. Furthermore, the magnitude of the precision improvement was directly related to the degree of correlation among the chemicals. PCBs had higher correlation among each other compared with PFAs and hence had more improvement in precision of the estimates. The improvement of precision is achieved by a combination of latent variable modeling and Bayesian techniques.[26,27] Latent variable modeling alone decreases variance greatly, although Bayesian procedures and prior distributions further stabilize the parameter estimates. Meanwhile, BKMR uses nonlinear smoothing techniques,[23] which resulted in less precise effect estimates compared with BFA. This is due to the additional amount of variance introduced by allowing nonlinearity in the kernel approximating functions. However, both BFA and BKMR performed better than the copollutant adjusted MLR models.

MLR showed poor estimate precision and is, therefore, inadequate for mixture analysis.[15-17] The precision of the MLR regression coefficients is drastically reduced after copollutant adjustment is made (Table 3). This inflation of variances in the regression estimates is directly related to the degree of collinearity present among the exposures.[20] Although the single pollutant regression models could generate more precise estimation, it provided biased estimates because it assumed the absence of the copollutant confounding. Additionally, the effect estimates generated in parallel by single pollutant model cannot be simply added arithmetically for mixture effects.[20]

Among the three methods, BFA had the clearest interpretation of the mixture effect estimates. It achieved this by simultaneously modeling the parameter estimates and error terms as depicted in Equations 1 and 2. The regression coefficients generated by BFA can be directly interpreted as the mixture-specific regression coefficients. For example, a change in birth weight for every 10-fold increase in PCB mixture concentration (consists of PCB 118, 138, 153, 170, and 180) is associated with a birth weight change of −81 g (95% CI = −132 to −31 g). Although a change in birth weight for every 10-fold increase in PFAS mixture concentration (consists of PFHXS, PFNA, PFOA, and PFOS) is associated with a birth weight change of −57 g (95% CI = −105 to −10 g). Additionally, the BFA modeling framework can be explicitly defined by researchers, making replications and comparisons of results across studies possible. BFA can also be used flexibly for different types of research questions including source-specific mixture studies where the latent construct is defined as the exposure source instead of the chemical structure, as was in this study.[54]

The latent variable in BFA is a form of dimensional reduction method that can capture information of all the individual chemicals and produce a single index representing the overall exposure of the mixtures.[40,54] The Bayesian framework applied further restrictions to the parameter estimates to reduce variance. In this study, an uninformative Gaussian prior with mean zero and variance 1,000 were used for the parameter estimates in the BKMR and BFA. This means that minimal prior information was introduced in our model to influence the final estimates and the information of our data dictate our analysis results. Therefore, the assumptions made in our Bayesian analysis were the same as the non-Bayesian analysis, with the exception that BKMR assumes nonlinear relationship while BFA assumes linear relationship. Mathematically, the amount of variance is reduced because only one exposure is modeled instead of all five within the mixture. This simplifies the interpretation and provides a more explicit and specific definition of mixture exposure by capturing all information of the individual chemicals instead of other dimensional reduction methods where a subset of the chemicals is selected.[15,20,55,56]

BKMR provided graphical outputs clearly depicting the dose-response curve. However, the interpretation is challenging because BKMR is an extremely flexible model. During the modeling procedure, different nonlinear transformations were used for each exposure. This is advantageous for detection and characterization of nonlinear effect in the dose-response curves of the chemical mixtures. However, when the actual relationship of the dose-response is linear, it could introduce additional complexities in the model through the flexible approximation functions. Furthermore, these approximation functions are dependent on the specific data structures of the cohort. This makes direct comparisons of analysis results from different cohorts challenging unless the characteristics of the cohorts are generally comparable and the same assumption of the biologic mechanism is made.

In terms of the interpretation of health effects in our analysis, our BFA results were coherent with our BKMR results. Both methods found that prenatal exposure to PCBs and PFAS were associated with reduced birth weight. Our findings are consistent with previous birth weight studies using the HOME Study data, although the magnitude of the effects may vary slightly due to a combination of reasons such as different data transformation or standardization techniques.[31,32] For example, in a single pollutant study, Rauch et al[31] found that every 10-fold increase in the concentration of selected OPs (∑DAP, ∑DEP, ∑DMP) was associated with slight negative but less precise associations with birth weight. In a mixture study, Kalloo et al[31] used principal component and clustering techniques and found that the principal component and cluster with mostly OCPs and phenol compounds were associated with a reduction of birth weight, although the 95% CI contained the null. Although the effect estimate computed by Kalloo et al[31] can be attributed to certain mixtures such as principal components and clusters, it is difficult to interpret the results since the mixture generated were dependent solely on the available data. The mixtures identified can be very different given different make-up of the study population, making direct comparisons of the results difficult.

Although the strength of BFA in estimation precision and interpretability is apparent, BFA has several limitations. BFA assumes a linear relationship among the latent mixture, which is informed by specific hypothetical causal diagrams and prior distributions of the parameters. If the actual relationship among these variables deviates significantly from the assumptions, BFA results may be biased. Another disadvantage of BFA is the relatively longer computation time needed for achieving model convergences when more variables and more sample sizes are supplied. It should also be noted that depending on the internal variance-covariance relationship present in the data, researchers may need to try different priors and to tune parameters such as the number of burn-in iterations and the ratio of adaptive and posterior sample size for the model to converge successfully. Additionally, our study contained limitations that could not be addressed by simply employing different methods. For example, some of the chemicals analyzed, such as phthalates, may vary during pregnancy. Measurement errors may exist in these nonpersistent chemicals because of their short half-lives and that measurements taken at a specific time may not reflect the actual amount of exposure.[57]

In terms of the generalizability for our analysis, the sample size of the HOME Study is relatively modest to examine multiple classes of environmental chemical mixtures. However, the concentrations of environmental chemicals in the HOME Study are similar to those in pregnant women in the US National Health and Nutrition Examination Survey at similar time of enrollment.[58] The estimation of mixture-specific regression coefficients without considering accompanying classes of chemical mixture is a limitation of the study. Ideally, a mixture analysis method includes all classes of chemical exposures, but computational burden is a hindrance. In the future, it is possible to combine different variable selection methods in multiple stages to enhance the estimation of regression coefficients for chemicals from different classes with an increasing sample size.

Conclusion

We examined three different statistical approaches to characterize and quantify the association between birth weight and prenatal exposures to seven classes of environmental chemical mixtures. We found that PCBs and PFAs displayed strong associations with reduced birth weight. We demonstrated the advantages of BFA in estimate precision and interpretability, although BKMR excels at visualizing dose-response relationships. Therefore, BFA and BKMR can complement each other to provide a more comprehensive interpretation of the mixture-specific effect. We also demonstrated the inadequacy of MLR for mixture assessment, especially in the presence of collinearity.