Assessing the health effects of environmental chemical mixtures is an important challenge in environmental epidemiology. It is quite difficult to find the relationship between exposure variables and a health outcome when the effects may be non-linear, interaction effects are present and may be subject to detection limits. In the past decade, there have been numerous approaches for analyzing this type of data. Hwang et al. [12] and Zhang et al. [21] proposed a latent class approach that links the exposure profile with disease severity with latent variables. Bobb et al. [3] employed Bayesian kernel machine regression that allows for flexible non-linear estimation of chemical effects on disease outcomes. Further, Carrico et al. [4] proposed a weighted quantile sum approach that generalizes using cumulative chemical exposure for predicting disease outcomes. Herring et al. [11] proposed a Bayesian regression approach that estimates linear main and interactive effects using a Dirichlet process prior and incorporates detection limits for the chemical exposures.

In this paper, we propose a Bayesian approach that incorporates non-linear exposure, interaction effects, and detection limits in a flexible way that does not require parametric assumptions on the exposure distributions. When the number of parameters are moderately large, the inclusion of all pairwise interaction terms results in sparsity that in turn results in poor estimation with maximum-likelihood estimation. We propose a shrinkage prior approach for estimating interactions. First, we investigate a shrinkage prior where interactions are treated the same as main effects, and no relationship between the two is incorporated. In a second approach, we propose a shrinkage prior that incorporates a relationship between the interaction and main effects to increase performance in sparse data situations. The hierarchical principle in linear models specifies that interactions will only be examined in situations where the corresponding main effects are sizable [15]. Models with interactions without main effects place restrictions on the parameter space that are not natural. We will show how incorporating this additional structure will provide efficiency gain when studying the interactions in chemical mixtures. Additionally, we show how to extend the shrinkage prior approach to situations where multiple parameters are used to model the effect of each chemical component on disease risk. This later extension can be used to model non-linear exposure effects as well as to provide a flexible approach for dealing with detection limits in mixtures. We propose a shrinkage prior approach that incorporates the hierarchical principle for the estimation of interactions for linear and nonlinear exposures as well as for exposures subject to detection limits.

In our paper, we describe our methodology in detail in Section 2. Next in Section 3, we discuss our prior specifications and posterior computations. In Section 4 we describe how to extend our model from linear where the relationship of an exposure has a single parameter to the situation where multiple parameters are associated with that exposure (e.g., nonlinear exposures or detection limits). In Section 5 we show the efficiency of our proposed method by simulation studies and results. Finally, we applied our methodology to the NCI-SEER NHL study and described our findings in Section 6.

Let Y = ( Y 1 , Y 2 , … , Y N ) ′ denotes the binary health response for N individuals and X i = ( X i 1 , X i 2 , … , X i p ) ′ be the corresponding p-dimensional vector of continuous chemical exposures. For p chemicals we have p ( p − 1 ) / 2 two-way interactions. We consider a logistic regression model with linear effects in their corresponding interactions of the following form: (2.1) logit P ( Y i = 1 | X i ) = U i ′ α + ∑ j = 1 p X i j β j ∗ + ∑ j = 1 p ∑ k = j + 1 p − 1 X i j X i k γ j k ∗ , i = 1 , 2 , … , N , where logit a = log a 1 − a , U i denotes q-dimensional covariate vector which includes an intercept, α is the corresponding regression coefficient vector, β j ∗ denotes the main effect regression coefficient of the jth chemical, and γ j k ∗ denotes the interaction effect regression coefficient of the jth and kth chemicals. We consider a latent variable approach [1] and approximate equation (2.1) using a robit link [14]. Let’s consider ω = ( ω 1 , ω 2 , … , ω N ) ′ be N-dimensional latent vector such that (2.2) Y i = 1 ω i > 0 , 0 otherwise , where, ω i = U i ′ α + ∑ j = 1 p X i j β j ∗ + ∑ j = 1 p ∑ k = j + 1 p − 1 X i j X i k γ j k ∗ + ϵ i . If ϵ i ∼ F t v , where F t v is a cumulative distribution function of student t-distribution with v degrees of freedom, it is called robit(v) regression [13], i.e. (2.3) P ( Y i = 1 | α , β ∗ , γ ∗ ) = 1 − P ( Y i = 0 | α , β ∗ , γ ∗ ) = F t v ( U i ′ α + ∑ j = 1 p X i j β j ∗ + ∑ j = 1 p ∑ k = j + 1 p − 1 X i j X i k γ j k ∗ ) , where β ∗ = ( β 1 ∗ , β 2 ∗ , … , β p ∗ ) and γ ∗ = ( γ 11 ∗ , γ 12 ∗ , … , γ p ( p − 1 ) / 2 ∗ ). As v → ∞, the robit(v) model becomes the probit regression model. Liu et al. [14] suggested that the robit link with v = 7 degrees of freedom closely approximates the logit link with β j = β j ∗ / 1.5484 and γ j k = γ j k ∗ / 1.5484. Moreover, we use the fact that the t-distribution can be represented as a scale mixture of normal distribution by introducing a mixing variable λ i , such that ϵ i | λ i ∼ N ( 0 , 1 λ i ) and λ i ∼ G ( v 2 , v 2 ), where N ( μ , σ 2 ) denotes a normal distribution with mean μ and variance σ 2 and G ( c 1 , c 2 ) denotes the gamma distribution with mean c 1 / c 2 and variance c 1 / c 2 2 to formulate the likelihood. For simplicity, let consider for the ith individual the interaction term between two exposure variables X i j and X i k is defined by Z i j k = X i j X i k and Z i = ( Z i 11 , Z i 12 , … , Z i p ( p − 1 ) / 2 ) ′ . Hence ω i | λ i ∼ N ( U i ′ α + X i ′ β + Z i ′ γ , 1 λ i ) and λ i ∼ G ( v 2 , v 2 ) where β = ( β 1 , β 2 , … , β p ) and γ = ( γ 11 , γ 12 , … , γ p ( p − 1 ) / 2 ). From equations (2.1) and (2.2) the complete data likelihood can be written as: (2.4) L ( Y | X , w ) = ∏ i = 1 N [ y i 1 { w i > 0 } + ( 1 − y i ) 1 { w i ≤ 0 } ] × ( 2 π ) − 1 2 λ i 1 2 exp ( − λ i 2 ( w i − U i ′ α − X i ′ β − Z i ′ γ ) 2 ) × ( ν 2 ) ν 2 Γ ( ν 2 ) λ i v 2 − 1 exp ( − λ i v 2 ) .

For linear models with interactions, the hierarchical principle implies that interactions should only be included if the corresponding main effects are non-zero [7, 9]. Hence, we consider a dependence structure between the main and interaction effects such that the inclusion of interaction effects depends on the inclusion of the corresponding main effects. To this end, we consider the following prior distributions: (3.1) β j ∼ N ( 0 , 1 a η j ) , γ j k ∼ N ( 0 , 1 b η j η k θ j k ) , η j ∼ G ( a 1 , b 1 ) , θ j k ∼ G ( a 2 , b 2 ) .

In the prior distribution in equation (3.1), a controls the global shrinkage towards the origin for the main effect regression coefficient β j , j = 1 , 2 , … , p, and η = ( η 1 , η 2 , … , η p ) ′ is the predictor specific local shrinkage parameter that allow deviations in the degree of shrinkage between predictors. a 1 , a 2 defines the shape parameter and b 1 , b 2 defines the scale parameter of the gamma distributions. We consider a G ( a 1 = 1 , b 1 = 1 ) distribution as a prior choice for η j with mean and variance 1 to induce some variability between the η j ’s. In this formulation, larger values of η j ’s will induce more shrinkage towards zero, while smaller values of η j result in minimum shrinkage to zero. This specification is based on global-local shrinkage framework [17], where typical recommendation is to consider a heavy tail distribution for the local shrinkage parameter to avoid over-shrinking large signals, and the global shrinkage parameter should have substantial mass near zero. Similarly, for the interaction effect regression coefficients γ j k , the global shrinkage parameter b shrinks all parameter towards zero. In contrast, the predictor specific local shrinkage parameter θ = ( θ 1 , θ 2 , … , θ p ( p − 1 ) / 2 ) ′ captures the interaction specific shrinkage effects. We consider a G ( a 2 = 1 , b 2 = 1 ) prior for θ j k . Note that, to share the information between main and interaction effects, the prior variance of γ j k is also dominated by the term η j η k . The parameter γ j k is shrunk to zero if at least one of η j , η k or θ j k is large. Consistent with the hierarchical principle, an interaction term will tend to be small if either β j or β k is small or θ j k is large. Furthermore, we also consider a G ( 1 , 1 ) prior on both global shrinkage parameter a and b, and a vague prior N ( 0 , 10 2 ) for the regression coefficient α j . The posterior distribution based on the complete data is given by (3.2) π ( α , β , γ , w | X , Y ) ∝ ∏ i = 1 N [ y i 1 { w i > 0 } + ( 1 − y i ) 1 { w i < = 0 } ] × λ i 1 2 exp ( − λ i 2 ( w i − U i ′ α − X i ′ β − Z i ′ γ ) 2 ) × λ i v 2 − 1 × exp ( − λ i v 2 ) × | Λ | − 1 2 exp ( − 1 2 β ′ Λ − 1 β ) × | Ω | − 1 2 × exp ( − 1 2 γ ′ Ω − 1 γ ) × | Ψ | − 1 2 exp ( − 1 2 α ′ Ψ − 1 α ) × ( ∏ j = 1 p η j a 1 − 1 exp [ − η j b 1 ] ) × ( ∏ j = 1 p − 1 ∏ k = j + 1 p θ j k a 2 − 1 exp [ − θ j k b 2 ] ) × ( a a 3 − 1 exp [ − a b 3 ] ) ( b a 4 − 1 exp [ − b b 4 ] ) . In equation (3.2), we define β = ( β 1 , β 2 , … , β p ) ′ , γ = ( γ 12 , γ 13 , … , γ p ( p − 1 ) ) ′ , Λ = diag ( 1 / a η 1 , 1 / a η 2 , … , 1 / a η p ), Ω = diag ( 1 / b η 1 η 2 θ 12 , 1 / b η 1 η 3 θ 13 , … , 1 / b η p − 1 η p θ ( p − 1 ) p ) and Ψ = 10 2 I q , where I q represents the q × q order identity matrix.

The proposed methodology can be easily extended to incorporate other prior distributions such as horseshoe, Cauchy, and Dirichlet-Laplace prior.

Previously in Section 2, we described our method in the context of a linear exposure and outcome relationship. In many realistic settings, more than one parameter will be needed to model individual exposures in the mixtures (e.g., quadratic exposures). Hence, in this section, we extend our methodology to capture those non-linear exposure-outcome relationships using the following logistic regression model: (4.1) logit P ( Y i = 1 | X i j ) = U i ′ α + ∑ j = 1 p g j ( X i j ) + ∑ j = 1 p ∑ k = j + 1 p − 1 f j k ( X i j , X i k ) .

For modeling mixtures with nonlinear exposure relationships, we can use a polynomial representation for the effect of each exposure. Polynomial effects can be incorporated in the main and interaction terms by using equation (4.1) with functions g j ( X i j ) = X i j ′ β j and f j k ( X i j , X i k ) = Z j k ′ γ j k and the logistic regression can be written in the following form: (4.2) logit P ( Y i = 1 | X i , Z i ) = U i ′ α + ∑ j = 1 p X i j ′ β j + ∑ j = 1 p ∑ k = j + 1 p − 1 Z j k ′ γ j k , i = 1 , 2 , … , N , where, X i j = ( X i j , X i j 2 ) ′ and Z j k = ( X i j X i k , X i j 2 X i k , X i j X i k 2 , X i j 2 X i k 2 ) ′ and the regression coefficients β j = ( β j 1 , β j 2 ) ′ and γ j k = ( γ j k 1 , γ j k 2 , γ j k 3 , γ j k 4 ) ′ .

Individual components of chemical mixtures may be subject to lower detection limits where the measurements are censored below these limits. Chiou et al. [6] and Ortega et al. [16] proposed a two-component exposure model for lower detection limits where the effect has one component indicating whether the exposure is above the detection limit and the other reflecting the value of the measurement if the exposure is above the detection limit. These authors showed that this parameterization does not make the unverifiable modeling assumptions inherent in treating lower detection limits as left censored in a parametric exposure distribution. We can incorporate detection limits into the mixture analysis by using equation (4.1) with: (4.3) f j k ( X i j , X i k ) = γ j k 1 I ( X i j ≥ C j ) I ( X i k ≥ C k ) + γ j k 2 ( X i j − C j ) I ( X i j ≥ C j ) I ( X i k ≥ C k ) + γ j k 3 ( X i k − C k ) I ( X i j ≥ C j ) I ( X i k ≥ C k ) + γ j k 4 ( X i j − C j ) ( X i j − C k ) I ( X i j ≥ C j ) I ( X i k ≥ C k ) .

In equation (4.3), β j 1 is the log odds of disease at the value of the detection limit relative to the log-odds of disease below the detection limit. Further, β j 2 is the log-odds ratio of disease for a one unit change in exposure above the detection limit. The parameter vector γ j k . measure the interactive effects between the jth and kth chemical. Specifically, γ j k 1 measures the interactive effect of being above the detection limit on both the jth and kth chemical, γ j k 4 measures the interactive effect of increasing X i j and/or X i k when both markers are above the detection limit. Finally, γ j k 2 and γ j k 3 are cross product interaction effects.

The shared shrinkage prior proposed in Section 3 can be extended to the multiple parameter per exposure case using the grouped shrinkage prior [5]. We show this for the non-linear formulation (4.2), but it applies more generally to other settings such as detection limits (4.3). Specifically, (4.4) ω i | λ i ∼ N ( U i ′ α + ∑ j = 1 p X i j ′ β j + ∑ j = 1 p ∑ k = j + 1 p − 1 Z j k ′ γ j k , 1 λ i ) , λ i ∼ Gamma ( v 2 , v 2 ) , β j ∼ MVN k 1 ( 0 , η j a I 2 ) , η j ∼ Gamma ( k 1 + 1 2 , 1 2 ) , γ j k ∼ MVN k 2 ( 0 , η j η k θ j k b I 4 ) , θ j k ∼ Gamma ( k 2 + 1 2 , 1 2 ) , α ∼ N ( 0 , 10 2 ) , a , b ∼ Gamma ( 1 , 0.1 ) .

Here k 1 and k 2 defines the dimension of the parameters β j and γ j k , respectively. In the non-linear case we share the information between main effect and interaction effect similarly as in linear hierarchical model in the equation (3.1). Posterior calculation follows as in Section 3.1.

We perform a series of the simulation studies to investigate the performance of the proposed methodology. We compare the proposed model with following two models. 1.

Independent Vague Prior: We incorporate the vague prior β j , γ j k ∼ N ( 0 , 10 2 ), where no dependence is introduced between the main and the interaction effects. This approach is approximately a maximum likelihood approach.

The NCI-SEER NHL study [18] is a population-based case-control study (508 controls & 672 cases) of non-Hodgkin lymphoma (NHL), was designed to determine the associations between exposures of chemicals/pesticides found in used vacuum cleaner bags and the risk of NHL. Different epidemiological pieces of evidence confirm that exposure to chemicals increases the risk of certain cancers in humans [20]. Often chemicals enter the household from indoor use or drift in from outdoor and may persist for months and years in carpet and cushion furniture without being degraded by sunlight, rain, and extreme temperatures. Hence carpet dust sampling provides a more objective basis for exposure assessment as it contains integrated pesticide exposure over a long period which is potentially more relevant to disease risk than recent or current exposure. In this study, the samples were collected from used vacuum cleaner bags of 672 cases in Detroit, Iowa, Los Angeles, and Seattle and were analyzed for pesticides [8]. Primarily the laboratory measurements contain missing data due to concentrations being below the minimum detection level. The “fill-in” approach [10, 8] was used for imputation where, for each biomarker, measurements below the detection limit were imputed by first estimating maximum-likelihood estimators of a log-normal distribution from the left-censored data and then imputing values below the detection limit based on this distribution.

Particularly for chemicals with a high percentage of values below their detection limits, results may not be robust to misspecification of the parametric assumptions. The median percent of observations below the detection limit was 61% (across chemicals) with a range of (3% to 93%). Some of the pesticides are highly correlated ( > 0.9) and skewed. Figure 1 shows the correlation between the chemicals. In that case, we chose to use only one of them for the analysis and log-transformed the exposure data. In the final data set, we have a total of p = 14 chemical exposures on N = 1180 individuals (508 controls & 672 cases). Figure 2 shows the correlation between the selected chemicals for the analysis. We considered site, sex, education, and age as covariates [8] in all models for our data application. First we estimated a main effect model using a vague prior and an independent shrinkage prior. Figure 3 shows only one chemical, Diazinon, had a 95% HPD interval that excludes zero, for all other chemicals the 95% HPD interval contained zero. Figure 5 shows a random subset of interaction effects estimated under all three prior distributions. The 95% HPD intervals included zero for all estimated interactions for all priors. The intervals for the shared were narrower than the independence shrinkage, and both were narrower than the vague prior.

The proposed multi-parameter exposure model that incorporates detection limits does not make strong assumptions on the distributions and exposure effects for values below the detection limit. Figures 4 and 6 show the results from fitting the multi-parameter model. These figures show a random set of interaction effects corresponding to the slope-slope interactions (corresponding to γ j k 4 parameters). Unlike for the linear exposure models presented, we see strong evidence for interactive effects. Figure 7 shows a contour plot obtained from the estimated terms corresponding to the interaction of D24 and Diazinon. The figure demonstrates a very strong qualitative interaction where exposure effects the risk of NHL only when an individual is exposed to both agents. (note that the log-odds ratios are near zero when either of the agents have exposure at or below their detection limit).

In this paper we have proposed a shrinkage prior that shares the information between the main effects and interaction effects for estimating the complex relationship between chemical mixtures and disease risk. We proposed a prior that shrinks the interaction term closer to zero when there is a little evidence of a corresponding main effect. This approach is consistent with the hierarchical principle that argues that one should only look for interactions when the corresponding main effects are present. Through simulations, this approach showed sizable efficiency gain in using the shared shrinkage prior when the hierarchical principle holds. Interestingly, we saw that even when there was no main effect, the shared shrinkage prior showed efficiency over priors that did not incorporate this shared information. We presume this is due to the flexibility of the shared shrinkage prior. Bien et al. [2] propose a penalized likelihood lasso approach for incorporating the hierarchical principle for parameter estimation. However, there are advantages to using the Bayesian approach over the likelihood approach is that i) there is no need to estimate the penalization constant, and ii) estimating measures of uncertainty are more straightforward using MCMC. Bayesian kernel machine regression (BKMR, [3]) is an alternative approach that can be used to identify two-way interactions that do not impose an additive structure in the model formulation. However, this approach does not incorporate the hierarchical principle in the estimation of interaction between mixture components. Further, BKMR cannot be easily extended to flexibly incorporate detection limits as discussed in Section 4.

A major contribution of this work is introducing the shared shrinkage prior to the setting of chemical mixtures. A major analytical issue in the analysis of chemical mixtures is the presence of lower detection limits among the multiple chemicals. We proposed two parameter exposure models that incorporate these detection limits. Specifically for each chemical, we include one parameter for the effect of being above the detection limit and another parameter for the linear change (slope) when above this limit. At the expense of adding additional parameters, the proposed approach makes no assumptions about the exposure effects and chemical distributions below detection limits.

The analysis of the NHL mixture data provided interesting insights into the methodology. When using a linear exposure model that imputed values below the detection limits, we found no interactive effects using the independence or shared shrinkage prior. These analyses were based on an imputation approach where chemicals were assumed to be log-normally distributed and values were imputed in the tail of this distribution. The percent of values below the detection limits varied across chemical but was generally very high (median 61%). We would expect the inferences using imputation and assuming a linear relationship to be highly sensitive to the imposed parametric assumptions. Ortega et al. [16] showed the lack of robustness of this type of imputation for a single exposure. When we fit the two-parameter exposure model that accounts for detection limits, we found substantially more evidence for interactive effects. Most of these interactions were with regard to D24. For example, we showed a very strong positive interactive effect between Diazinon and D24 (see Figure 7) that was missed with the linear exposure with the imputation approach.

In the simulations, we considered the case where there were 10 chemicals ( p = 10); the example contained 14 chemicals in the final models. The analysis results were insensitive to the chosen prior distributions. We used a Gamma ( 1 , 1 ) for the local and global shrinkage priors, but results were very similar when we used either Gamma ( 1 , 0.5 ) with variance 2 or Gamma ( 1 , 2 ) with variance 0.5 distributions. In the simulation studies, with 10 chemicals we have a large number of interactive effects (45 for linear exposure and 180 for two parameter exposure model). With these large number of interactive effects, the simulations and example showed the importance of using shrinkage priors rather than a vague prior that is closely related to maximum-likelihood. For a larger number of chemicals, even the shrinkage prior approaches may have poor performance due to sparsity. In this case, a Bayesian variable selection method [17] could be used to a smaller number of chemical exposures to include in the model.

We considered models with two-way interactions. Conceptually, we could consider 3 or high-order interactions and extend the shared shrinkage idea to this more complex situation. In the spirit of Stone’s generalized additive model [19], the proposed hierarchical model could be extended to include B-splines rather than polynomial effects. Such an extension would be computationally intensive and is left for future work.

Data will be available upon request and data sharing agreement. For R code, see the github page.