Prenatal Phthalate Exposures and Adiposity Outcomes Trajectories: A Multivariate Bayesian Factor Regression Approach

TL;DR

Prenatal exposure to phthalates may influence childhood obesity, with effects varying by age and sex.

Contribution

A novel Bayesian multivariate factor regression approach is used to model phthalate mixtures and their effects on adiposity outcomes.

Findings

In boys, prenatal phthalate exposure is linked to lower adiposity at younger ages but higher adiposity after age 7.

No significant associations were found between phthalate exposure and adiposity outcomes in girls.

The Bayesian method effectively models non-linear time-varying effects and quantifies uncertainty.

Abstract

Experimental animal evidence and a growing body of observational studies suggest that prenatal exposure to phthalates may be a risk factor for childhood obesity. Using data from the Mount Sinai Children’s Environmental Health Study (MSCEHS), which measured urinary phthalate metabolites (including MEP, MnBP, MiBP, MCPP, MBzP, MEHP, MEHHP, MEOHP, and MECPP) during the third trimester of pregnancy (between 25 and 40 weeks) of 382 mothers, we examined adiposity outcomes—body mass index (BMI), fat mass percentage, waist-to-hip ratio, and waist circumference—of 180 children between ages 4 and 9. Our aim was to assess the effects of prenatal exposure to phthalates on these adiposity outcomes, with potential time-varying and sex-specific effects. We applied a novel Bayesian multivariate factor regression (BMFR) that (1) represents phthalate mixtures as latent factors—a DEHP and a non-DEHP factor, (2) borrows information across highly correlated adiposity outcomes to improve estimation precision, (3) models potentially non-linear time-varying effects of the latent factors on adiposity outcomes, and (4) fully quantifies uncertainty using state-of-the-art prior specifications. The results show that in boys, at younger ages (4–6), all phthalate components are associated with lower adiposity outcomes; however, after age 7, they are associated with higher outcomes. In girls, there is no evidence of associations between phthalate factors and adiposity outcomes.

Figures

Full text

1. Introduction

The increasing prevalence of childhood obesity is a relevant public health concern. The 2017–2018 National Health and Nutrition Examination Surveys (NHANES) estimated that 19.3% of children aged 2–19 years have obesity and another 16.1% are overweight [1]. The prevalence of obesity among children and adolescents has increased more than three times in the last 30 years [2]. Obese children are at increased risk for several health conditions, including obesity in adulthood, type 2 diabetes, heart disease, arthritis, and various cancers, as well as shorter life. Studies suggest that approximately 70% of obese children face a significant risk of heart disease in adulthood [3]. Some estimates suggest that one-third of children born today (and half of Latino and black children) are expected to develop type 2 diabetes at some point in their lives [4]. Childhood obesity also has negative implications for mental health, as it is associated with an increased risk of being bullied [3].

Environmental chemical exposures, particularly exposure to phthalates, have been identified as risk factors for childhood obesity [5,6], potentially by interfering with the body’s endocrine system [7,8]. Phthalates are a group of chemicals widely used in consumer products, such as toys, food packaging, or cosmetics, that are known to have endocrine disrupting properties [9,10]. In vivo and in vitro animal studies suggest that fetal development is a potential critical window of vulnerability to phthalate exposures, which may promote obesity, and the effect may depend on gender [8,11,12]. Observational studies in humans, including birth cohort studies, further support these findings by showing associations between prenatal phthalate exposure and childhood adiposity outcomes [13,14,15,16], some reporting sex-specific effects [14,17,18].

In this paper, we investigate the time-varying health effects of prenatal phthalate exposures on childhood obesity using data from the Mount Sinai Children’s Environmental Health Study (MSCEHS) [13]. The study measured the urinary phthalate metabolites of mothers during the third trimester of pregnancy and followed the health indices of children, including adiposity outcomes, between ages 4 and 9. We would like to identify groups of phthalate metabolites critical to these adiposity outcomes and estimate their time-varying effects and any sex-specific effects using a new modeling framework: the Bayesian multivariate factor regression model (BMFR). Our approach addresses the challenge of fully quantifying the uncertainty in estimating the time-varying effects of the mixture and provides several advantages tailored to our cohort data.

Many studies have examined the time-varying effects of exposure to phthalates on adiposity outcomes or growth trajectories, but the results remain inconsistent [14,18,19,20,21,22,23,24,25]. A key limitation is that traditional statistical models, or even mixture methods not designed to estimate time-varying effects, lack full uncertainty quantification. In traditional statistical models, key modeling choices, such as the functional form of time-varying effects, are chosen or estimated first, then treated as fixed in the subsequent outcome model. For example, linear mixed-effects models or generalized estimating equations typically rely on interaction terms between a single exposure and age to model time-varying effects [21,22,23,25]. These approaches cannot consider all phthalate metabolites simultaneously due to multicollinearity, and assumptions about the functional form of time-varying relationships (e.g., linear or polynomial) are treated as fixed. Fixing these choices underestimates uncertainty, leading to confidence intervals that are too narrow and p-values that are artificially small. More advanced methods, such as growth mixture models, latent class growth models, or functional principal component analysis, can capture nonlinear growth trajectories, but they require a second-stage model to estimate associations between mixture exposures and outcome trajectories [18,24]. In such multistage analyses, outputs from earlier stages (e.g., a number of selected trajectories or factors) are treated as fixed in later stages, again leading to underestimated uncertainty. Other advanced mixture methods, such as Bayesian kernel machine regression [26], quantile g-computation [27], and weighted quantile sum regression [28], do not model nonlinear time-varying effects. Bayesian varying-coefficient kernel machine regression (BVCKMR) is another advanced mixture method that can model nonlinear effects across exposure levels on outcome trajectories [29]. However, at fixed exposure levels, BVCKMR’s estimated effects over time are limited to a functional form.

To address these limitations, BMFR applies state-of-the-art prior specifications to infer modeling decisions and fully quantify their uncertainty. The model includes variable selection priors for covariates, avoiding fixed subsets [30]; Gaussian process priors for flexible modeling of time-varying health effects without assuming a specific functional form [31]; and half-t priors for robust between-subject variances [32]. To model the exposure mixture, BMFR assumes that structured variations in correlated exposures can be attributed to a small number of latent factors that also explain part of the variations in the outcomes, while fully quantifying the uncertainty in the number of factors through the multiplicative gamma process (MGP) prior [33]. This contrasts with principal component analysis (PCA) [20], which maximizes variation in exposures without regard to outcome relevance, and structural equation models [34], which require a fixed number of factors.

Our implementation of BMFR also has custom features specific to our cohort. BMFR jointly models multiple adiposity outcomes, including BMI z-scores, waist circumference, waist-to-hip ration, and fat mass percentage, to improve estimation precision. Previous studies have generally examined these outcomes individually [14,16,17,20,21,23,35], hypothesized that mixtures affect all outcomes similarly [20], or searched for consistent results across outcomes [16,21,23]. Because these outcomes are highly positively correlated and likely reflect shared mechanisms through which phthalates influence obesity, modeling them jointly is expected to improve estimation precision. Finally, BMFR quantifies uncertainty in imputing missing data and measurements below the limit of detection (LOD). Instead of relying on standard fixed-value imputations (e.g., mean or one-half of the LOD), BMFR handles imputation within the Markov Chain Monte Carlo (MCMC) process, propagating uncertainty in imputed values into the credible intervals for the effects of interest. An R package, optimized in C++, for BMFR is freely available at https://github.com/phuchonguyen/famr (accessed on 20 August 2023). Section 2 describes in detail the data from the MSCEHS. Section 3 describes our proposed model BMFR and its prior specifications. Section 4 validates BMFR’s utility through simulation studies. Section 5 describes the analysis of time-varying health effects of prenatal phthalate exposures on childhood obesity using the MSCEHS data.

2. Data

2.1. Study Population

Between 1998 and 2002, MSCEHS recruited 479 first-time mothers with singleton pregnancies from the Mount Sinai Diagnostic and Treatment Center and two adjacent private practices in New York City. Among these women, 75 were excluded due to medical complications (n = 3), infant or fetal death (n = 2), very premature birth (before 32 weeks of gestation or <1500 g; n = 5), miscarriage (n = 1), delivery of an infant with genetic abnormalities or malformations (n = 5), inability to obtain biologic samples before delivery (n = 12), relocation or transfer to a hospital outside of New York City (n = 28), or loss to follow-up (n = 19) [13]. The final cohort consisted of 404 babies with birth data recorded. The children were invited back for three follow-up visits at ages 4–5.5, 6, and 7–9. Of the 404 babies, 382 had their mothers’ prenatal concentrations of phthalate metabolites measured in urine. Two additional observations were excluded because they had very dilute urine (<10 mg/dL creatinine) that can produce inaccurate biomarker measurements [14]. Of the 382 babies, only 180 came back for at least one follow-up visit. This results in 362 total visits for fat mass percentage, 363 total visits for body mass index (BMI), 364 total visits for waist-to-hip ratio, and 364 total visits for waist circumference. Figure 1 shows the pattern of loss to follow-up of these 382 babies. All observed outcome measurements were included in our analysis. Our analysis was a secondary analysis of the de-identified data from the MSCEHS study.

2.2. Phthalate Exposures

Mothers who were pregnant between 25 and 40 weeks (with a mean of 31.5 weeks) provided a urine sample that was analyzed by the CDC laboratory for various phthalate metabolites, including MEP, MnBP, MiBP, MCPP, MBzP, MEHP, MEHHP, MEOHP, and MECPP. The DEHP group consists of MEHP, MEHHP, MEOHP, and MECPP. To account for inaccuracies in analytical standards, correction factors of 0.72 and 0.66 were applied to MBzP and MEP concentrations and limits of detection (LOD), respectively [36]. Urinary concentration was measured using creatinine. To adjust for the dilution of urine samples, we standardize the metabolites’ concentrations by a Cratio, as well as include creatinine concentration as a covariate in the analysis as suggested by [37]. The Cratio is calculated as the ratio between predicted creatinine conditional on observed covariates of the mother including the mother’s age, mother’s BMI, mother’s gestational weight gain adequacy category, mother’s smoking status, mother’s education, mother’s race, and observed creatinine [37]. As seen in Figure 2 (Left), phthalate metabolites are highly and positively correlated with each other, especially those within the DEHP group and those within the non-DEHP group. There are four MiBP, one MEP, one MBzP, four MCPP, one MECPP, one MEHHP, one MEOHP, and 15 MEHP measurements under their respective LOD of detections in total. We impute these within the MCMC sampler as discussed later.

2.3. Adiposity Outcomes

The MSCEHS assessed the weight and body composition of the infants through bioelectrical impedance analysis using a pediatric Tanita scale at three follow-up visits scheduled at approximately 4 years, 6 years, and 7 years, though the age the children at actual visits ranges from 4 to 10 years (or 48 to 122 months). We consider the following four outcomes in our multivariate analysis: fat mass percentage (FMP), BMI z-score (BMIz), waist-to-hip ratio (WHR), and waist circumference (WC). FMP is based on fat mass estimates reported by the Tanita scale (model TBF-300; Tanita Corporation of America) and calculated as (fat mass/weight) $[eqn]$ . BMI is calculated as weight (in kilograms)/height (in meters)^2^. It is then standardized by age and sex using a CDCSAS (version 9.3) macro [38] to produce BMIz. There are 362 total visits for FMP, 363 for BMIz, 364 for WHR, and 364 for WC from 180 subjects across all their follow-up visits. Figure 2 (Right) shows that these outcomes (except WHR) are highly and positively correlated with each other.

2.4. Covariates

Mothers were interviewed for 2 h during enrollment to gather covariate data. The computerized perinatal database at Mount Sinai Hospital was used to obtain pregnancy and delivery information. Gestational weight gain adequacy is calculated by dividing the observed gestational weight gain (last pregnancy weight minus self-reported pre-pregnancy weight) by the expected gestational weight gain based on the 2009 Institute of Medicine guidelines times 100 [39]. We categorize gestational weight gain as inadequate if the ratio is <86%, adequate if 86–120%, and excessive if >120%. The final baseline covariates include the mother’s age, mother’s BMI, mother’s gestational weight gain adequacy category, mother’s smoking status, mother’s education, mother’s race, whether mother breastfed, child’s sex, child’s birth weight, and creatinine concentration. The children’s age in months was also recorded at each follow-up visit and included as a covariate. Table 1 summarizes the baseline covariates of children with at least one follow-up visit included in our analysis. Most covariates have no missing values or only negligible amounts, except for maternal gestational weight gain. The baseline characteristics of the male and female subsamples are similar, with only minor differences in race/ethnicity distribution.

3. Bayesian Multivariate Factor Regression for Time-Varying Effects

3.1. Model Correlated Chemical Mixtures with a Latent Factor Model

Let $[eqn]$ be a p-vector of correlated metabolite concentrations measured during the third trimester of pregnancy for subject i. We assume that variation in $[eqn]$ can be attributed to $[eqn]$ latentvariables:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is a K-vector of unobserved latent factors of subject i, $[eqn]$ is the factor loadings matrix, and $[eqn]$ are idiosyncratic noise variances. We assume the exposures have been mean-centered and remove the intercepts. Independent priors on each $[eqn]$ are chosen to be those often used in factor analysis. We use the multiplicative gamma process (MGP) prior on the factor loadings to learn sparse loadings structure and infer the number of factors K [33]:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

3.2. Model Correlated Outcomes as a Function of Latent Factors

Let $[eqn]$ be a vector of q outcomes at follow-up time $[eqn]$ , where $[eqn]$ is the number of follow-ups with at least one measured outcome for subject i. We assume each outcome at each follow-up has been mean-centered and remove the intercepts. We also assume the variation in the outcomes $[eqn]$ can be decomposed into the variation explained by the latent factors of $[eqn]$ , the variation due to unobserved factors and idiosyncratic noise:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ has a half-t prior with a small degree of freedom for an uninformative prior that still behaves well in the case that between-subject variance $[eqn]$ is close to zero, as suggested by [32]. Random variables $[eqn]$ are subject-level random intercepts. $[eqn]$ is the residual covariance, which describes variances and covariances in the outcome due to unmeasured factors as well as random noise. We can interpret $[eqn]$ as the proportion of total residual variation explained by between-subject variation, and $[eqn]$ as the proportion of residual variation explained by within-subject variation.

3.3. Model Health Effects as Flexible Functions of Time

B is a $[eqn]$ matrix of regression functions, which is of primary interest in our analysis. Element $[eqn]$ in B models the effect of the $[eqn]$ latent factor on the $[eqn]$ outcome that can vary smoothly and flexibly over follow-up times. We consider Gaussian processes as priors to learn these smooth regression functions on a discrete-time grid $[eqn]$ , where T is the total number of unique ages at which children had follow-up visits. At the same time, we want to incorporate our belief that effects across adiposity outcomes, which are driven by similar mechanisms through which phthalates interfere with the body’s hormones, should be correlated. As a result, instead of placing independent Gaussian process priors on elements of B, we adopt the following factorization:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is a $[eqn]$ matrix, with $[eqn]$ , that linearly combines H independent basis functions into elements of B. A similar factorization was used by [40], but their work focused on covariance regression. We also use the MGP prior on the basis functions loadings to learn sparse loadings structure and help infer the number of basis functions H [33]. U is a matrix of independent basis functions. We choose the Gaussian kernel $[eqn]$ for all elements of U to ensure the time-varying effects are smooth functions of time with the same wiggliness encoded in a shared length scale $[eqn]$ . Since the input is on a grid, we place a uniform prior on a grid of plausible values for $[eqn]$ . Conditional on $[eqn]$ , the $[eqn]$ column of B has a separable Gaussian process prior:

$$[eqn]$$

where $[eqn]$ describes the covariance in the regression functions of factor k across outcomes, and $[eqn]$ . Thus, $[eqn]$ , so we set the amplitude of the kernel to 1 for identifiability.

3.4. Model Linear Effects and Interactions of Covariates

Let $[eqn]$ be a L-vector of covariates, including both baseline covariates and those collected at follow-up time t. For covariates where the linear relationships are reasonable, we can add a linear effect term to Equation (4) as follows:

$$[eqn]$$

We endow the regression coefficient matrix $[eqn]$ with a global-local shrinkage prior on matrix normal parameters of [30]:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Shrinkage parameter $[eqn]$ provides variable selection to determine if predictor l is important to all outcomes, which fits our application. Outcome-specific effects within column l can additionally shrink toward zero. The authors in [30] suggested setting the global shrinkage parameter $[eqn]$ to satisfy sufficient conditions for posterior consistency [30]. When $[eqn]$ , this is the horseshoe prior [30,41]. A similar setup can be used for any linear interactions between the latent factors and covariates.

3.5. Imputation of Censored and Missing Data

We have very few missing outcomes at recorded follow-up visits (two missing FMP and one missing BMIz measurement out of 364 observations). In the case that it is reasonable to assume outcomes are missing at random conditional on observed covariates and birth weight data [42], we impute them during MCMC. We sample $[eqn]$ given $[eqn]$ from a conditional multivariate normal, where $[eqn]$ are all unknown parameters:    

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ are parameters corresponding to the indices of the missing values, $[eqn]$ are parameters at observed indices, $[eqn]$ is the covariances between missing and observed indices, $[eqn]$ is the covariance matrix of missing indices, and $[eqn]$ is the covariance matrix of observed indices.

Moreover, we often observe censored metabolite concentrations that are below the limit of detection (LOD). The LOD is defined as the lowest concentration of an analyte in a sample that can be reliably distinguished from the highest concentration of the same analyte in a sample with no such analyte [43]. We can impute metabolite concentrations under the LOD by sampling from a conditional truncated normal at each MCMC iteration:

$$[eqn]$$

where $[eqn]$ is the LOD of the $[eqn]$ chemical, $[eqn]$ is the $[eqn]$ row of loading matrix $[eqn]$ , and $[eqn]$ is a truncated normal distribution with mean m, variance v, and support  $[eqn]$ .

3.6. Posterior Computation

See Appendices for full conditional and adaptive Metropolis-within-Gibbs updates.

4. Simulations

We compare our proposed method BMFR with the following approaches: (1) two-stage univariate regressions, (2) BVCKMR [29], and (3) a baseline mean model. In the two-stage approach, we reduce the dimension of the exposures using PCA, keeping the first few principal components (PCs), and then fit LMM for each outcome separately. The baseline mean model returns the mean of each outcome. We simulate data from the following three scenarios to demonstrate our method’s utility compared to existing approaches. For all scenarios, we generate data for $[eqn]$ subjects. We generate the exposures according to a factor model, creating correlated exposures with group structures, similar to the observed phthalates metabolites:

$$[eqn]$$ $$[eqn]$$

We generate sparse $[eqn]$ so that every five metabolites load onto one factor for $[eqn]$ and $[eqn]$ latent factors. The non-zero entries of $[eqn]$ are sampled from $[eqn]$ . We generate $[eqn]$ outcomes measured at $[eqn]$ time points for all subjects:

$$[eqn]$$ $$[eqn]$$

We use different exposure-response functions $[eqn]$ including time-varying effects, different distributions for the random intercept $[eqn]$ , and random error $[eqn]$ for each scenario. Below is the description of each simulation scenario:

• Scenario 1:Linear exposure-response function where the first two factors are important, independent responses:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the latent factor-response function. The induced effects of X range from 0 to 1. Under this setting, all assumptions for PCA-LMM are satisfied, though it does not propagate the uncertainty from the first stage.

• Scenario 2:Non-linear time-varying exposure-response function where the first two factors are important, positively correlated responses:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the standard normal pdf at t. Covariances $[eqn]$ and $[eqn]$ have a composite symmetry structure with a high correlation of 0.7. Under this setting, all assumptions for our BMFR model are satisfied.

• Scenario 3:Quadratic exposure-response function where three metabolites are important and responses are independent, the most favorable scenario for BVCKMR:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

We choose $[eqn]$ for PCA and our method, and $[eqn]$ for our method, to resemble analyses where $[eqn]$ are close to but not exactly the true latent dimension. For predictive performance evaluation, we calculate the mean predictive square error (MPSE) on a test set of 200 subjects at 10 time points. Additionally, to evaluate how well the methods measure the relative importance of each chemical, we calculate the Spearman correlation between the true relative importance rank and the inferred rank of chemicals in X. The rank is based on the absolute value of the effect at each time point, summed over all time points. We calculate the rank for our model by first calculating the induced effects in the original predictors X at each time point as shown in [44]. Similarly, we can calculate the effects from the PCA-LMM analysis in X as $[eqn]$ , where $[eqn]$ is the first K left singular vectors, and $[eqn]$ is a vector of regression coefficients for the PCs. Table 2 shows the MPSE results. Note that PCA-LMM performed worst in all scenarios, even when all its assumptions were met. Our method performs best when the data are generated according to its model. When the data are more favorable to the other models, our method still performs better than PCA-LMM. Table 3 shows the Spearman correlation results. Here, our method performs best in the first two scenarios and is very close to the best in the third scenario.

5. Analysis of MSCEHS Cohort Data

5.1. Data Preprocessing

We log-transformed WC so that its marginal is more approximately normal. We mean-centered all outcomes. We also used the logarithm of phthalate metabolites as exposures. We corrected for urinary dilution by dividing the metabolite concentrations (not on the log scale) by the Cratio as discussed in Section 2. We standardize other continuous covariates and create dummy variables for categorical covariates. We used R package mice (version 3.18.0) for multiple imputation of missing values in the covariates (using predictive mean matching for continuous, logistic regression for binary, and proportional odds model for ordered categorical covariates). We created a time variable that is age in years based on age in months of the children at follow-up visits. Age ranges from 4 to 10 years old (Figure 3).

5.2. Preliminary Analysis

As a preliminary analysis, we performed two-staged PCA-LMM analyses of each of the four outcomes. We applied PCA to the standardized log chemical exposures. We used the first three principal components (PCs) for the LMM stage because they explained over 90% of the variations in the exposures. The first three PCs can be interpreted as the non-DEHP (excluding MEP) factor, DEHP factor, and MEP factor, respectively. The PCA factor loading matrix is available in Appendix B.

We fitted four independent LMMs for four outcomes with random intercepts and interactions between the PCs and the child’s gender, controlling for all baseline covariates. We fitted a second set of LMMs with interactions between the PCs, child’s gender, and child’s age, and a third set with interactions between the PCs, child’s gender, and polynomials of degree 2 of child’s age. We used AIC, BIC, likelihood ratio tests, and 6-fold cross-validated MPSE for model comparison within each outcome. After Bonferroni correction for multiple testing, we still saw evidence from the likelihood ratio tests that models with linear interactions with age were the best fits for FMP, WHR, and WC. Models with linear interactions with age had the best cross-validated MPSE for BMIz, FMP, and WC. We saw insufficient evidence of the quadratic interactions in age being useful. Full summary tables of AIC, BIC, p-values, and MPSEs can be found in Appendix B.

5.3. Main Analysis via BMFR

In the main analysis, we fitted our proposed Bayesian multivariate factor regression with time-varying effects to all four outcomes simultaneously. We fitted two models, one for male and one for female children, to assess any sex-specific effects. We controlled for linear main effects of covariates and included random intercepts.

We used the following prior specifications. Since we standardized the log chemical exposures to have unit variances, we set hyperparameters for inverse gamma priors on idiosyncratic noise variances $[eqn]$ so that they are less than 1 with 99% probability. For hyperparameters on the MGP prior for loadings $[eqn]$ and $[eqn]$ , we used $[eqn]$ and $[eqn]$ as suggested in the note by [45]. For the IW prior on $[eqn]$ , we set $[eqn]$ to the sample covariance of the four outcomes, and a small $[eqn]$ so that prior is loosely centered around the sample covariance. For a weakly informative prior that $[eqn]$ should be below 100, we use a half-Cauchy with a scale of 25 as in [32].

We selected the number of factors K and the number of basis functions H using grid search. We chose the combination of K and H from a set of options ( $[eqn]$ , and $[eqn]$ ) that had the smallest 6-fold cross-validated MPSE, with $[eqn]$ only. We considered small values for H because of our prior belief that the effects of the exposures on all outcomes were similar. The final model for males had $[eqn]$ and the one for females had $[eqn]$ . We fixed the length scale parameter $[eqn]$ , since we observed from our preliminary analysis that there were likely no highly variable effects. Sensitivity analysis was also performed for models that infer the length scale parameter during MCMC.

5.4. Results of Analysis via BMFR

For out-of-sample predictive performance comparison, we calculated MPSE of a baseline mean model that returned outcome-specific means. The 6-fold cross-validated MPSE of our proposed method is 0.93 compared to 0.98 of the baseline mean model. For the female model, our proposed method’s MPSE is 0.82 compared to 1.04 of the baseline mean model. On the whole data set (combining male and female), our MPSE is 0.88, compared to 0.92 of PCA-LMMs with linear interactions with age, and 1.03 of the baseline mean model.

For interpretation of the results, we resolved rotational and label-switching ambiguity in the factors using the MatchAlign algorithm proposed by [46]. Figure 4 shows the post-processed factor loading matrix and time-varying effects of each latent factor from the model fitted for male children. Time-varying effects were transformed to represent the effects of one unit increase in the latent factors. We identified two main latent factors corresponding to the non-DEHP group and DEHP group of chemicals. This is consistent with previous studies on phthalate [20]. Though both groups of chemicals seem to be associated with lower values in adiposity outcomes at younger ages, their promotive effects of obesity seem to increase over time. The latent factors identified in the model for female also includes non-DEHP and DEHP groups, though all effects seem to remain null over time (Figure 5). The sex-specific effects here could be related to phthalates being anti-androgens [47]. Though the C.I.s from the female model are not statistically significant, the effects of latent factors on WHR seem to be different from their effects on the other three outcomes. This could be related to the fact that WHR is less correlated to the other outcomes than they are to each other.

We included results from sensitivity analysis of inferring the length scale parameter during MCMC in Appendix B. Overall, the sensitivity analysis results were similar to the results in Figure 4 and Figure 5.

6. Discussion

This paper assesses the time-varying health effects of prenatal phthalate exposures on adiposity outcomes measured in children from ages 4 to 9 from the MSCEHS cohort study using the BMFR approach we propose. BMFR represents phthalate mixtures as latent factors—a DEHP and a non-DEHP factor—and borrows information across highly correlated adiposity outcomes to improve estimation precision, models potentially non-linear time-varying effects of the latent factors on adiposity outcomes, and fully quantifies uncertainty using state-of-the-art prior specifications. We find that in boys, at younger ages (4–6 years), all phthalate latent factors (DEHP and non-DEHP) show negative associations with adiposity outcomes. After age 7, these associations begin to become positive. In girls, there is no evidence of associations between phthalate components and outcomes. We also find these time-varying effects to be similar across all adiposity outcomes (BMIz, fat mass percentage, waist-to-hip ratio, and waist circumference). Our introduction of a new Bayesian mixture method for estimating time-varying effects with full uncertainty quantification and our finding of sex-specific time-varying associations of prenatal phthalate exposures with childhood obesity between age 4 to 9 are novel.

We were able to estimate sex-specific time-varying effects not previously identified in analyses of the MSCEHS cohort because BMFR is customized to this research question and includes several custom specifications for the data. BMFR applies prior specifications to infer modeling decisions and avoid the artificially narrow confidence intervals that are an unintended consequence of fixing choices in traditional statistical models or multistage analyses. The model includes variable selection priors for covariates, avoiding fixed subsets; Gaussian process priors for flexible modeling of time-varying effects without assuming a functional form; and half-t priors for robust between-subject variances. BMFR represents highly correlated exposures with a small number of latent factors that are independent of each other but correlated with the outcomes, while fully quantifying uncertainty in the number of components through the MGP prior. BMFR also jointly models multiple outcomes of adiposity—BMIz, waist circumference, waist-hip ratio, and fat mass percentage—available in our cohort to improve estimation precision. Finally, BMFR quantifies the uncertainty in imputing missing data and values below the limit of detection (LOD) within the MCMC sampling, propagating this uncertainty into the credible intervals for the effects of interest. An R package for BMFR is freely available at https://github.com/phuchonguyen/famr (accessed on 20 August 2023). The package implements the MCMC algorithm in Appendix A in C++ for optimal computational speed.

Our finding of time-varying effects is similar to results from [25], which reported that exposure to phthalate in the first trimester was associated with lower BMI six months after birth but higher BMI in older ages, although they did not observe specific sex effects. Other studies have also found a somewhat similar time-varying effect: higher maternal urinary phthalate concentrations associated with lower fetal growth and birth weight, followed by higher growth trajectories later on [48,49], with sex-specific effects [48]. Our results also align with [20], which found non-DEHP metabolites associated with lower BMI, fat mass, and waist circumference in boys aged 5 and 7, and with [50], which reported DEHP associated with higher outcomes among boys between ages 8 and 10, and provide a potential explanation for seemingly inconsistent results between the two studies. The time-varying effect we identify in boys may appear null when aggregating across all ages, which is consistent with a previous analysis of this cohort [42] that estimated average effects over time and did not observe associations with percent fat mass or sex-specific modification. However, there are other studies inconsistent with ours, including reports of DEHP associated with a lower BMI in girls [14] and a higher weight gain in early childhood that stabilized during puberty in girls [24]. Our findings may have important implications for pregnancy care guidelines and child health. The time-varying effects of prenatal exposure suggest that these exposures may influence childhood obesity years after the time of exposure. They further suggest that the third trimester of pregnancy may be a vulnerable window of exposure, making interventions to reduce phthalate exposure during pregnancy important.

This study also has limitations and opportunities for future improvements. BMFR does not model nonlinear dose-response relationships. As a result, we did not investigate different effects at different exposure levels, though this could be performed by stratifying the analysis by tertiles of exposure. This may be important, as mixture effects could vary in nonlinear ways across exposure levels [26,29]. In this cohort, exposures were measured from a single urine sample, with collection times ranging from 25 to 40 weeks of gestation. Given the short half-lives of urinary phthalate metabolites [51] and the likely episodic nature of exposures [42], this limits the precision of exposure measurement. Future research should consider multiple exposure measurements and aim to identify the most critical window of vulnerability during the fetal period to growth trajectories. Data on adiposity outcomes before age 4 and after age 9 are not available in this cohort. Stratification by sex also reduces sample size and power to detect small effects in our analysis. Our findings would be strengthened by replication across a longer time window, from birth through puberty, and with data from a larger cohort study. Finally, residual confounding may remain, as we could not account for child calorie intake, or maternal consumer product preferences, which could influence exposure levels [42].

7. Conclusions

This paper presents a Bayesian multivariate factor regression approach to assessing the time-varying health effects of prenatal phthalate exposures as measured in maternal urine sample during the third trimester of pregnancy on adiposity outcomes measured in young children from age 4 to 9 using data from the MSCEHS cohort study. BMFR addresses challenges in analyzing mixture exposures by representing them as latent factors that predict the outcomes. It also allows non-linear time-varying effects of exposure mixture to be estimated with full uncertainty quantification, while improving estimate precision by borrowing information across correlated adiposity outcomes. The results show that in boys, at younger ages (4–6), all phthalate components show an association with lower adiposity outcomes; however, after age 7, they begin to show an association with higher outcomes. In girls, there is no evidence of associations between phthalate components and adiposity outcomes.