# Assessing the Longitudinal Impact of Environmental Chemical Mixtures on
Children’s Neurodevelopment: A Bayesian Approach

Wei Jia
Division of Biostatistics and Bioinformatics
University of Cincinnati
<EMAIL_ADDRESS>Roman A. Jandarov
Division of Biostatistics and Bioinformatics
University of Cincinnati
<EMAIL_ADDRESS>
##  1 Abstract

This manuscript presents a novel Bayesian varying coefficient quantile
regression (BVCQR) model designed to assess the longitudinal effects of
chemical exposure mixtures on children’s neurodevelopment. Recognizing the
complexity and high-dimensionality of environmental exposures, the proposed
approach addresses critical gaps in existing research by offering a method
that can manage the sparsity of data and provide interpretable results. The
proposed BVCQR model estimates the effects of mixtures on neurodevelopmental
outcomes at specific ages, leveraging a horseshoe prior for sparsity and
utilizing a Bayesian method for uncertainty quantification. Our simulations
demonstrate the model’s robustness and effectiveness in handling high-
dimensional data, offering significant improvements over traditional models.
The model’s application to the Health Outcomes and Measures of the Environment
(HOME) Study further illustrates its utility in identifying significant
chemical exposures affecting children’s growth and development. The findings
underscore the potential of BVCQR in environmental health research, providing
a sophisticated tool for analyzing the longitudinal impact of complex chemical
mixtures, with implications for future studies aimed at understanding and
mitigating environmental risks to child health.

##  2 Introduction

Environmental chemical exposures can affect children’s health such as
neurological development. One important characteristic in environmental
exposure research is that humans are not exposed to an isolated single
chemical. As such, modeling and understanding environmental mixtures are a
priority area of the National Institute of Environmental Health Science
(NIEHS) [6][18]. Many methods have been proposed and applied to address the
mixture effect of chemical exposures to children’s neurodevelopment but few
have analyzed longitudinal effects of mixtures. Specifically, we want to not
only estimate a static effect of chemical mixture on outcome measured at a
single time point, but also estimate effect of chemical mixture on outcome at
a given age.
Numerous approaches have been proposed to estimate mixture effects. Two
popular model classes are response-surface methods and exposure-index methods
[15]. Response-surface approaches model complex mixtures and their
relationship with outcomes non-parametricaly. One representative approach is
Bayesian Kernal Machine Regression (BKMR) [3] which estimates a multi-
dimensional mixture-response surface with a kernel approach. On the other
hand, exposure-index methods estimates mixture effect with linear exposure
mixture index. The most popular method is Weighted Quantile Sum Regression
(WQS) [7]. Both classes of methods have achieved state-of-art results.
However, they are designed to solve time-invariant problems.
Previous models have embedded methods for modeling chemical mixtures with
time-varying models in order to solve longitudinal problems. Lagged kernel
machine regression [14, 13] combined kernel machine with lagged windows
approach to estimate effects of chemical exposures measured at multiple time
points on outcome measured at a single time point. Bayesian Varying
coefficient kernel machine regression (BVCKMR) [12] combined kernel machine
regression with a mixed effect model to estimate effects of exposures measured
at a single time point on outcome measure at multiple time points. It
estimates mixture effect for both baseline age and longitudinal trajectory.
However, kernel methods are known to perform bad when the dimensionality of
data gets high[1]. Also, the interpretation of BVCKMR results can be difficult
since the interpretation mainly depends on posterior analysis. Additionally,
one known limitation of BVCKMR is while modeling high dimensional chemical
exposure data, sparsity of the mixture effect plays an important role in the
analysis. However, methods we mentioned above do not perform any shrinkage or
variable selection to handle the sparsity.

In this paper, we propose a Bayesian varying coefficient quantile regression
(BVCQR) which can (1) characterize longitudinal effects of chemical exposure
mixtures, (2) provide succinct estimations when the dimensionality of data
gets high and preserve the convenience of interpretations, (3) apply a
horseshoe prior to ensure the sparsity of chemical effects estimation, (4)
correct uncertainty quantification using Bayesian method (5) we illustrated
the robustness of BVCQR on estimating exposure effects with a simulations
study and a real-life data ”HOME” cohort.

##  3 Method

In this section, we will start by introducing some related previous
approaches, weighted quantile sum regression (WQS), Quantile G-Computation
(QGC), mixed effect model, and Bayesian varying coefficient kernel machine
regression (BVCKMR). In section 2.3, we introduce our BVCWQS, and details of
how we handle sparsity with horseshoe prior are given in section 2.4.

###  3.1 Weighted Quantile Sum Regression and Quantile G-Computation

For subject i, let $y_{i}$ be the corresponding outcome,
$\\{z_{i1},\ldots,z_{iM}\\}$ be chemical exposures, and $x_{i}$ be a vector of
covariates for the $ith$ subject where $i=1,\ldots,N$. Exposure index models
such as WQS and QGC fit a linear model as:

$y_{i}=\alpha_{0}+\beta(\theta_{1}q_{i1}+\ldots{}+\theta_{M}q_{iM})+\bm{x}_{i}^{T}\gamma+\epsilon_{i}$
(1)

where $q_{ip}$ represent quantile indices of $z_{ip}$ such as $0,1,2,3$ if
$z_{ip}$ is in the 0-25 percentile, 25-50 percentile, 50-75 percentile and
75-100 percentile. $\beta$ represents the overall association between the
weighted sum and outcome. Different strategies and constrains have been
applied to obtain weights and overall association for WQS and QGC. For WQS,
the constrain has been set on weights is $\sum_{m}^{M}\theta_{m}=1$ and
$\theta_{m}>0$. Alternatively, QGC relaxes this single directional constrain
and estimates coefficients for each $q_{ij}$. Component weights are calculated
by $\theta_{j}=\frac{\theta_{j}}{\sum_{j}^{M}\theta_{M}}$. These methods ease
the interpretation because they characterize a single parameter overall
mixture effect and a set of weights to illustrate how each chemical contribute
to the overall mixture effect.

###  3.2 Mixed Effect Model

A mixed effect model is generalization of standard linear regression model
that can be used to analyze and model data which generated from multiple
source of variation [11, 16]. It relates one response variable (outcome) with
some explanatory variables. What sets the mixed-effects model apart is its
incorporation of both fixed and random factors. Fixed effects estimate average
responses across groups, akin to a standard regression model, while random
effects (e.g., individual-specific effects) accommodate natural variability in
responses among different entities, facilitating response estimation for each
entity in the analysis. Given that measurements within the same entity are
often correlated, this correlation must be considered during modeling. The
nature of correlation among repeated responses can vary, leading to diverse
covariance or correlation structures. The model supports various covariance
structures and allows for the estimation of effects and variance parameters.
The dataset’s structure can range from balanced (equal observations per
entity) to unbalanced (unequal observations per entity).

###  3.3 Bayesian varying coefficient kernel machine regression

Let $y_{ij}$ be the outcome of interest where $i=1\ldots n$ for each subject
and measured at age $j=1\ldots J_{i}$. The outcome is related with M
components of environmental chemicals $z_{i}=(z_{1i}\ldots z_{Mi})^{T}$.
$x_{i}=(x_{1i}\ldots x_{pi})^{T}$ are $p$ numbers of clinical covariates.
BVCKMR is an extension of BKMR [3] for analyzing longitudinal exposure data
which achieved by employing random effect model. The model can be writen as:

$y_{ij}=\gamma_{1}+\gamma_{2}*age_{ij}+h_{1}(z_{1i},\ldots,z_{Mi})+h_{2}(z_{1i},\ldots,z_{Mi})*age_{ij}+\bm{x_{i}^{T}}\bm{\beta}+\bm{u_{ij}^{T}b_{i}}+\epsilon_{ij}$
(2)

where $h_{q}(.)$, $q=1,2$ are mixture effects which estimated by using kernel
method. $h_{1}$ represents the mixture effect at the outcome baseline, and
$h_{2}$ represents the mixture effect over outcome trajectory.
$\gamma_{1}+\gamma_{2}*age_{ij}$ represents a global trend line. The model can
be simplified as

$\bm{Y}_{i}=\bm{X}_{i}\bm{\beta}+\bm{W}_{i}\bm{h}_{i}+\bm{U}_{i}\bm{b}_{i}+\bm{\epsilon}_{i}$
(3)

where $\bm{h}$ is the vector of mixture effects for each individual. $\bm{b}$
is the vector which contains random effects of each individual. $\bm{X}_{i}$is
the covariates matrix for subject $i$. $\bm{U}$ and $\bm{W}$ are matrices
represent ages of subjects. $\bm{Y}_{i}$ is a vector of continues outcome for
subject i, where $\bm{Y}_{i}=(Y_{i1},\ldots,Y_{iJ_{i}})^{T}$ and
$\bm{\epsilon}_{i}\thicksim N_{J_{i}}(0,\sigma^{2}\bm{I}_{J_{i}})$. $\bm{h}$
can be hard to estimate by kernel function when M becomes large due to the
curse of dimensionality. Also, kernel based method is know as its high
computational cost.

###  3.4 Bayesian varying coefficient quantile regression

Let $y_{ij}$ be the outcome of interest where $i=1\ldots n$ for each subject
and measured at age $j=1\ldots J_{i}$. The outcome is related with M
components of environmental chemicals $z_{i}=(z_{1i}\ldots z_{Mi})^{T}$.
$x_{i}=(x_{1i}\ldots x_{pi})^{T}$ are $p$ numbers of clinical covariates. Our
data has exposures measured at a single time point and outcome is repeatedly
measured on multiple time points. Consider a random effect model from [12]
which can be written as

$y_{ij}=\gamma_{1}+\gamma_{2}*age_{ij}+h_{1}(.)+h_{2}(.)*age_{ij}+\bm{x_{i}^{T}}\bm{\beta}+\bm{u_{ij}^{T}b_{i}}+\epsilon_{ij}$
(4)

where $h_{l}(.)$, $l=1,2$ are mixture effects and we defined it as:

$h_{l}(.)=\theta_{l1}q_{1i}+\ldots+\theta_{lM}q_{Mi}+\phi_{l}$ (5)

where $\bm{q}_{i}=(q_{1i}\ldots q_{Mi})^{T}$ is a vector of quantile
transformed $\bm{z}_{i}$ such as $q_{1i}=0,1,2,3$ if $z_{1i}$ is in the 0-25
percentile, 25-50 percentile, 50-75 percentile or 75-100 percentile,
$\bm{\theta_{l}}=(\theta_{l1},\ldots,\theta_{lM})$ is estimated coefficients
for chemicals, $\phi_{l}$ is the variance of the mixture effect. $h_{1}$
represents the mixture effect at the baseline age, and $h_{2}$ represents the
mixture effect over time. $\gamma_{1}+\gamma_{2}*age_{ij}$ represents a global
trend line corresponding to fixed effects. Following [12], we can simplify our
model as:

$\bm{Y}_{i}=\bm{X}_{i}\bm{\beta}+\bm{W}_{i}\bm{h}_{i}+\bm{U}_{i}\bm{b}_{i}+\bm{\epsilon}_{i}$
(6)

where $\bm{h}$ is the vector of mixture effects for each individual. $\bm{b}$
is the vector which contains random effects of each individual. $\bm{X}_{i}$is
the covariates matrix for subject $i$. $\bm{U}$ and $\bm{W}$ are matrices
represent ages of subjects. $\bm{Y}_{i}$ is a vector of continues outcome for
subject i, where $\bm{Y}_{i}=(Y_{i1},\ldots,Y_{iJ_{i}})^{T}$ and
$\bm{\epsilon}_{i}\thicksim N_{J_{i}}(0,\sigma^{2}\bm{I}_{J_{i}})$. The global
trend line can be included by adding columns of age and intercept in to
$\bm{X}_{i}$.

$\bm{Y}=\begin{pmatrix}\bm{Y}_{1}\\\ \bm{Y}_{2}\\\ \vdots\\\
\bm{Y}_{n}\end{pmatrix}=\begin{pmatrix}Y_{11}\\\ \vdots\\\ Y_{1J_{1}}\\\
\vdots\\\ Y_{n1}\\\ \vdots\\\
Y_{nJ_{n}}\end{pmatrix},\bm{\beta}=\begin{pmatrix}\beta_{1}\\\ \vdots\\\
\beta_{n}\end{pmatrix},\bm{X}=\begin{pmatrix}\bm{X}_{1}\\\ \vdots\\\
\bm{X}_{n}\end{pmatrix}$

$\bm{U}=\begin{pmatrix}1&age_{11}&&&&&\\\ \vdots&\vdots&&&&&\\\
&&1&age_{21}&&&\\\ &&\vdots&\vdots&&&\\\ &&1&age_{2J_{1}}&&&\\\
&&&&\ddots&&\\\ &&&&&1&age_{n1}\\\ &&&&&\vdots&\vdots&\\\
&&&&&1&age_{nJ_{n}}\\\ \end{pmatrix}$

$\bm{b}=\begin{pmatrix}b_{1,1}\\\ b_{2,1}\\\ b_{1,2}\\\ b_{2,2}\\\ \vdots\\\
\vdots\\\ b_{1,n}\\\
b_{2,n},\end{pmatrix},\bm{W}=\begin{pmatrix}1&0&\ldots&0&age_{11}&0&\ldots&\\\
\vdots&&&&\vdots&&&\\\ 1&0&\ldots&0&age_{1J_{1}}&0&\ldots&\\\
0&1&0&\ldots&0&age_{21}&0&\ldots&\\\ 0&\vdots&&&&\vdots&&&\\\
0&1&0&\ldots&0&age_{2J_{2}}&0&\ldots\\\ 0&0&1&0&\ldots&0&age_{31}&0\\\
0&0&\vdots&&&&\vdots&\ddots\\\ \end{pmatrix},\bm{h}=\begin{pmatrix}h_{1,1}\\\
h_{1,2}\\\ \vdots\\\ h_{1,n}\\\ h_{2,1}\\\ h_{2,2}\\\ \vdots\\\
h_{2,n}\end{pmatrix}$

The proposed model can be expressed as a hierarchical model:

$\bm{y}|\bm{h},\bm{X},\bm{\beta},\sigma^{2},\bm{b}\thicksim
N(\bm{X\beta}+\bm{Wh}+\bm{Ub},\sigma^{2}\bm{I})$ (7)
$\bm{h}|\bm{\theta}_{1},\bm{\theta_{2}},\phi_{1},\phi_{2}\thicksim
N_{2n}\begin{pmatrix}\begin{pmatrix}\bm{q}^{T}\bm{\theta}_{1}\\\
\bm{q}^{T}\bm{\theta}_{2}\end{pmatrix},\begin{pmatrix}\phi_{1}^{2}\bm{I}_{n}&0\\\
0&\phi_{2}^{2}\bm{I}_{n}\\\ \end{pmatrix}\end{pmatrix}$ (8)

We applied horseshoe prior for $\bm{\theta}_{1},\bm{\theta}_{2}$ to handle the
sparsity[8] and will be illustrated in the next subsection.

$\phi_{1}^{2},\phi_{2}^{2}\thicksim gamma^{-1}(\alpha_{0},\gamma_{0})$ (9)
$\bm{b}|\sigma^{2},\bm{D}\thicksim
N_{2n}(\bm{0},\sigma^{2}(\bm{I}_{n}\bigotimes\bm{D}))$ (10)

where covariance matrix D follows:

$\bm{D}^{-1}\thicksim Wishart_{q}(v_{0},C_{0})$ (11) $\sigma^{2}\thicksim
gamma^{-1}(\alpha,\gamma)$ (12) $\bm{\beta}\varpropto 1$ (13)

###  3.5 Handling Sparsity of Mixture Effects with Horseshoe Prior

Consider the $\bm{\theta}_{1}$ and $\bm{\theta_{2}}$ in equation 8 which are
effects of chemicals at baseline and trajectory. We assume both of
$\bm{\theta}_{1}$ and $\bm{\theta_{2}}$ should be sparse. The horseshoe prior
is set as for each $\theta_{lm}$ in
$\bm{\theta}_{l}^{T}=(\theta_{l1},\ldots,\theta_{lM})$ where $m=1\ldots M$

$\theta_{lm}|\lambda_{lm},\tau_{l}\thicksim N(0,\lambda_{lm}^{2},\tau_{l})$
(14)

$\tau_{l}$ is global parameter which shrinks all $\theta_{lm}$ toward zero;
$\lambda_{lm}$ is local parameter which allow some $\theta_{lm}$ to escape the
shrinkage. The level of shrinkage is controlled by $\tau_{l}$ as large value
of $\tau_{l}$ will reflect little shrinkage and value close to zero will
shrink all $\theta_{lm}$ to zero. Follow the parameterization in [17], we set
half-t priors for $\lambda_{lm}$ and $\tau_{l}$

$\lambda_{lm}\thicksim t^{+}(df_{\lambda},0,1)$ (15) $\tau_{l}\thicksim
t^{+}(df_{\tau},0,a_{0}\phi_{l}^{2})$ (16)

where $df_{\lambda}$ and $df_{\tau}$ are degree of freedoms which we set to 1.
The scale parameter for $\tau_{l}$ are set as scaled variance of the
corresponding $h$ function and $\lambda_{lm}$ has unit scale.

The BVCQR is implemented in Stan and all posterior distributions are computed
by Hamiltonian-Monte Carlo algorithm which has been known as the one of most
efficient algorithms for reducing correlation between sample states and
therefor ensure the stability of the estimates [5, 2]

##  4 Simulations

We conducted a simulation study to evaluate the performance of our proposed
BVCQR model for estimating $\bm{h}_{l},l=1,2$. We simulated 36 chemical
exposures based on our ”HOME” cohort. We first compute the covariance matrix
$\bm{C}_{36\times 36}$ of 36 exposures. Then we simulated exposure data as
$Z_{i}\thicksim N_{36}(\bm{0},\bm{C})$. Then we took the quantized version of
$\bm{Z}$ which is $\bm{q}$ . Outcome $Y$ is generated from the model
$y_{ij}=h_{1}(q_{1i},\ldots,q_{Mi})+h_{2}(q_{1i},\ldots,q_{Mi})*age_{ij}+\bm{x_{i}^{T}}\bm{\beta}+\bm{u_{ij}^{T}b_{i}}+\epsilon_{ij}$
where $\epsilon_{ij}\thicksim N(0,1)$, $\bm{x}_{j}=(x_{1j},x_{2,i})$, and
$x_{1i}\thicksim N(0,1)$ and $x_{2i}\thicksim Bernoulli(0.5)$. We simulated
data on two scenarios:
Scenario 1:

$\bm{h}_{1}(q)=5*(0.1q_{1}+0.3q_{12}+0.4q_{24}+0.2q_{35})$

$\bm{h}_{2}(q)=3*(0.1q_{9}+0.5q_{23}+0.2q_{27}+0.2q_{33})$
Scenario 2:

$\bm{h}_{1}(q)=5*(0.1q_{1}+0.3q_{12}-0.3q_{23}+0.4q_{24}+0.2q_{35})$

$\bm{h}_{2}(q)=3*(-0.2q_{6}+0.1q_{9}+0.5q_{23}+0.2q_{27}+0.2q_{33})$
Scenario 1 is designed to test the robustness of BVCWQS when modeling
exposures with only positive effects on the outcome, scenario 2 is designed
for exposures have both positive and negative effects.

For each scenario, we simulated $n=100$ samples. The performance of BVCQR on
estimating $\bm{h}$ was evaluated by regressing the predicted $\hat{\bm{h}}$
on $\bm{h}$. Intercept, slope and root mean squared error (RMSE) for $h_{1}$
and $h_{2}$ are presented. Close to zero intercept, and close to one slope and
$R^{2}$ indicate good estimations. Results for all scenarios are presented in
table 1. The results indicate a good performance of BVCWQS. For all scenarios,
the intercept ranges from -0.167 to 0.003 ; the slope ranges from 0.998 to
1.023; the $\bm{R}^{2}$ ranges from 0.987 to 0.993; the RMSE ranges from 0.186
to 0.309. In both scenarios, intercepts are close to zero and slope and
$\bm{R}^{2}$ are close to one, and their RMSE keep low, for both $\bm{h}_{1}$
and $\bm{h}_{2}$

Simulation Scenario | h | Intercept | Slope | $R^{2}$ | RMSE
---|---|---|---|---|---
1 | With hs prior | 1 | -0.167 | 1.023 | 0.992 | 0.287
|  | 2 | -0.023 | 1.003 | 0.992 | 0.186
| Without hs prior | 1 | -0.240 | 1.031 | 0.982 | 0.434
|  | 2 | -0.023 | 1.004 | 0.988 | 0.236
2 | With hs prior | 1 | -0.048 | 1.010 | 0.987 | 0.309
|  | 2 | -0.003 | 0.998 | 0.993 | 0.189
| Without hs prior | 1 | -0.092 | 1.016 | 0.974 | 0.434
|  | 2 | -0.008 | 1.001 | 0.989 | 0.235

Table 1: Simulation results, regression of $\hat{\bm{h}}$ on $\bm{h}$

The simulation results are further plot and illustrated in Figure 1 and Figure
2, in which we plot out the exposure effects on the outcome for each specific
chemical (The coefficients for chemicals inside function $\bm{h}$) for
scenario 1 and scenario 2. For each scenario, estimated effects and their
$95\%$ confidence intervals from BVCQR are plotted as black point , true
effects are plotted as red points in figures. In both scenarios, chemicals
which are related with outcome are significantly estimated by BVCQR (95%
confidence intervals do not contain zero) and all other chemicals which are
not related with outcome are not significant in the model (95% confidence
intervals contain zero).

Figure 1: Effects plot of Scenario 1: Upper plot is the baseline effect (h1)
and the lower plot is the effect on age trajectory (h2). X axis shows names of
simulated chemical exposures and Y axis shows values of each exposure effect.
True effects are given as red dots and estimated effects are plot in black
dots as well as their 95% confidence.  Figure 2: Effects plot of Scenario 2:
Upper plot is the baseline effect (h1) and the lower plot is the effect on age
trajectory (h2). X axis shows names of simulated chemical exposures and Y axis
shows values of each exposure effect. True effects are given as red dots and
estimated effects are plot in black dots as well as their 95% confidence.

To illustrate the sparsity handling by horseshoe prior, we conducted ablation
study on simulation 1 and 2. We run BVCQR with uninformative prior specified
for $\bm{\lambda}_{1}$ and $\bm{\lambda}_{2}$. The results of linear
regression model of $\hat{\bm{h}}$ and $\bm{h}$ have been given in Table 1.
Our proposed BVCQR regression with horseshoe prior out performs the
implementation with uninformative prior on both scenario 1 and scenario 2 on
intercept, slope and $R^{2}$ closer to 0 and with lower RMSE. Effects
estimation for scenario 1 and scenario 2 by BVCQR without horseshoe prior are
given in Figure 4 and Figure 5. larger confidence intervals are estimated by
BVCQR without horseshoe prior compare to BVCQR regression with horseshoe
prior. Two unwanted effects, q3 and q15, which are not associated to outcome
in the simulation, are estimated as significant effects at baseline (Figure 4
and Figure 5).

Figure 3: Effects plot of Scenario 1 estimated by BVCQR without horseshoe
prior: Compare to the estimations from BVCQR with horseshoe prior, these
estimations have larger confidence interval. q3 and q15 were not used to
generate outcome in the simulation, however, they are falsely estimated at
baseline. Figure 4: Effects plot of Scenario 2 estimated by BVCQR without
horseshoe prior: Compare to the estimations from BVCQR with horseshoe prior,
these estimations have larger confidence interval. q3 and q15 were not used to
generate outcome in the simulation, however, they are falsely estimated at
baseline.

##  5 Application on HOME data

We applied BVCQR to investigate the association between chemical exposure and
children’s growth and development on the Health Outcomes and Measures of the
environment (HOME) Study, which is a prospective cohort from three hospitals
in greater Cincinnati, Ohio, metropolitan area. environmental samples were
collected during pregnancy and the Bayley PDI score were collected when the
children were at 12, 24 and 36 month. We used a sub-sample of this cohort (n=
147). Follow the preprocess procedure from [4], we excluded chemicals with
less than $20\%$ samples with detectable values and end up with 36 continuous
chemicals variables. We assigned a value of $LOD/\sqrt{2}$ for concentrations
below the limit of detection[10]. We then rescaled these variables by dividing
them by 2 times their SD[9]. We encoded the age 24 months as baseline. The
BVCQR is implemented in Stan and compiled with R. The Hamiltonian Monte Carlo
sampler was run 2000 iterations, with first 1000 used as burn-in.
We fitted the data with a quantile version of linear mixed effect model (LMM)
with random slope and random intercept as a preliminary analysis. We want to
identify the association of chemicals and children’s neurodevelopmental
outcome at baseline (age 24 month) and neurodevelopmental tragectory (age 12
to 36 month) while adjustion sex, gestational age and birth weight. The
baseline effects can be estimated as main exposure effects in the LMM and
effects on neurodevelopmental tragectory can be estimated as the interaction
effects of exposures and age. The results from LMM are given in Figure 3 to 7.
At the baseline, PFNA has positive association with Bayley PDI at 95%
significance level; at the tragectory, PCB28 has negative association at 95%
significance level.

Figure 5: Linear Mixed model estimated Bayley PDI at 24 months and Bayley PDI
trajectory (12-36 months main effects and 95% confidence intervals of phthbpa
group chemicals. Figure 6: Linear Mixed model estimated Bayley PDI at 24
months and Bayley PDI trajectory (12-36 months main effects and 95% confidence
intervals of pcb group chemicals. PCB28 has significant negative association
at tragectory according to LMM estimation. Figure 7: Linear Mixed model
estimated Bayley PDI at 24 months and Bayley PDI trajectory (12-36 months main
effects and 95% confidence intervals of oc group chemicals. Figure 8: Linear
Mixed model estimated Bayley PDI at 24 months and Bayley PDI trajectory (12-36
months main effects and 95% confidence intervals of pbde group chemicals.
Figure 9: Linear Mixed model estimated Bayley PDI at 24 months and Bayley PDI
trajectory (12-36 months main effects and 95% confidence intervals of pfc
group chemicals. PFNA has significant association with Bayley PDI at 24 month.

We then use our proposed BVCQR model to fit the data. The posterior means of
the parameters of the global trend line which corresponding to intercept and
slope of individual-level fixed effects were $\gamma_{1}=1.86$ and
$\gamma_{2}=0.33$. The effects estimations are given in Figure 10 to Figure
14. Compare to estimations from LMM, our proposed BVCQR shrinks most of
insignificant chemicals to exactly zero or very close to zero. Both BVCQR and
LMM detected PFNA’s positive association at 24 month. Contrast to the negative
association of PCB28 at tragectory which detected by LMM, BVCQR detected
positive association of PCB28 at 24 month.

Figure 10: BVCQR estimated Bayley PDI at 24 months and Bayley PDI trajectory
(12-36 months main effects and 95% confidence intervals of pfc group
chemicals. phthbpa has significant association with Bayley PDI at 24 month.
Figure 11: BVCQR estimated Bayley PDI at 24 months and Bayley PDI trajectory
(12-36 months main effects and 95% confidence intervals of pcb group
chemicals. PCB28 has significant positive association at tragectory according
to LMM estimation. Figure 12: BVCQR estimated Bayley PDI at 24 months and
Bayley PDI trajectory (12-36 months main effects and 95% confidence intervals
of oc group chemicals. Figure 13: BVCQR estimated Bayley PDI at 24 months and
Bayley PDI trajectory (12-36 months main effects and 95% confidence intervals
of pbde group chemicals. Figure 14: BVCQR estimated Bayley PDI at 24 months
and Bayley PDI trajectory (12-36 months main effects and 95% confidence
intervals of pfc group chemicals. PFNA has significant association with Bayley
PDI at 24 month.

##  6 Discussion

We developed a Bayesian hierarchical model to analyze the effect of high
dimensional chemical exposures on the outcome longitudinal trajectories. We
handling sparsity of chemical exposures in $h$ function with horseshoe prior.
We showed that our proposed BVCQR robustly model chemical exposure effects as
well as handling sparsity on both outcome baseline and trajectory with
simulations. We also included uncertainty quantification and preserved
interpretbility of BVCQR. Using BVCQR to analyzing home data, we found similar
set of significant chemicals compare to a quantile version linear mixed effect
model, while our BVCQR shrinks effects of all insignificant chemicals to
either zero or very close to zero.
There are also some other factors which can be considered and studied in the
future to improve the function of BVCQR: (1) chemicals are belong to different
groups, chemicals belong to the same group are with similar chemical
structures which can lead to similar effect on the health outcomes. This
information can be used to perform the variable selection. (2) interaction and
non-linear effect can be estimated by incorporating other strategies.
Previously, people have applied Multiple Index model [15] on exposure data to
deal with these two factors. We are exploring to incorporate multiple index
model with BVCQR to solve these two problems under the longitudinal
environment.

## References

  * [1] Yoshua Bengio, Olivier Delalleau, and Nicolas Le Roux.  The curse of dimensionality for local kernel machines.  Techn. Rep, 1258(12):1, 2005.
  * [2] Michael Betancourt.  A conceptual introduction to hamiltonian monte carlo.  arXiv preprint arXiv:1701.02434, 2017.
  * [3] Jennifer F Bobb, Linda Valeri, Birgit Claus Henn, David C Christiani, Robert O Wright, Maitreyi Mazumdar, John J Godleski, and Brent A Coull.  Bayesian kernel machine regression for estimating the health effects of multi-pollutant mixtures.  Biostatistics, 16(3):493–508, 2015.
  * [4] Joseph M Braun, Amy E Kalkbrenner, Allan C Just, Kimberly Yolton, Antonia M Calafat, Andreas Sjödin, Russ Hauser, Glenys M Webster, Aimin Chen, and Bruce P Lanphear.  Gestational exposure to endocrine-disrupting chemicals and reciprocal social, repetitive, and stereotypic behaviors in 4-and 5-year-old children: the home study.  Environmental health perspectives, 122(5):513–520, 2014.
  * [5] Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng.  Handbook of markov chain monte carlo.  CRC press, 2011.
  * [6] Danielle J Carlin, Cynthia V Rider, Rick Woychik, and Linda S Birnbaum.  Unraveling the health effects of environmental mixtures: an niehs priority, 2013.
  * [7] Caroline Carrico, Chris Gennings, David C Wheeler, and Pam Factor-Litvak.  Characterization of weighted quantile sum regression for highly correlated data in a risk analysis setting.  Journal of agricultural, biological, and environmental statistics, 20:100–120, 2015.
  * [8] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  Handling sparsity via the horseshoe.  In Artificial intelligence and statistics, pages 73–80. PMLR, 2009.
  * [9] Andrew Gelman.  Scaling regression inputs by dividing by two standard deviations.  Statistics in medicine, 27(15):2865–2873, 2008.
  * [10] Richard W Hornung and Laurence D Reed.  Estimation of average concentration in the presence of nondetectable values.  Applied occupational and environmental hygiene, 5(1):46–51, 1990.
  * [11] SAS Institute.  SAS/STAT software: changes and enhancements, volume 229.  Sas Institute, 1992.
  * [12] Shelley H Liu, Jennifer F Bobb, Birgit Claus Henn, Chris Gennings, Lourdes Schnaas, Martha Tellez-Rojo, David Bellinger, Manish Arora, Robert O Wright, and Brent A Coull.  Bayesian varying coefficient kernel machine regression to assess neurodevelopmental trajectories associated with exposure to complex mixtures.  Statistics in medicine, 37(30):4680–4694, 2018.
  * [13] Shelley H Liu, Jennifer F Bobb, Birgit Claus Henn, Lourdes Schnaas, Martha M Tellez-Rojo, Chris Gennings, Manish Arora, Robert O Wright, Brent A Coull, and Matt P Wand.  Modeling the health effects of time-varying complex environmental mixtures: Mean field variational bayes for lagged kernel machine regression.  Environmetrics, 29(4):e2504, 2018.
  * [14] Shelley H Liu, Jennifer F Bobb, Kyu Ha Lee, Chris Gennings, Birgit Claus Henn, David Bellinger, Christine Austin, Lourdes Schnaas, Martha M Tellez-Rojo, Howard Hu, et al.  Lagged kernel machine regression for identifying time windows of susceptibility to exposures of complex mixtures.  Biostatistics, 19(3):325–341, 2018.
  * [15] Glen McGee, Ander Wilson, Thomas F Webster, and Brent A Coull.  Bayesian multiple index models for environmental mixtures.  Biometrics, 79(1):462–474, 2023.
  * [16] Chava Peretz, A Goren, T Smid, and H Kromhout.  Application of mixed-effects models for exposure assessment.  Annals of Occupational Hygiene, 46(1):69–77, 2002.
  * [17] Juho Piironen and Aki Vehtari.  On the hyperprior choice for the global shrinkage parameter in the horseshoe prior.  In Artificial Intelligence and Statistics, pages 905–913. PMLR, 2017.
  * [18] Kyla W Taylor, Bonnie R Joubert, Joe M Braun, Caroline Dilworth, Chris Gennings, Russ Hauser, Jerry J Heindel, Cynthia V Rider, Thomas F Webster, and Danielle J Carlin.  Statistical approaches for assessing health effects of environmental chemical mixtures in epidemiology: lessons from an innovative workshop.  Environmental health perspectives, 124(12):A227–A229, 2016.