Bayesian Multiple Index Models for Environmental Mixtures
BBiometrics
Bayesian Multiple Index Models for Environmental Mixtures
Glen McGee , ∗ , Ander Wilson , Thomas F. Webster , and Brent A. Coull Department of Statistics and Actuarial Science, University of Waterloo, Waterloo, ON, Canada Department of Statistics, Colorado State University, CO, U.S.A. Department of Environmental Health, Boston University, Boston, MA, U.S.A. Department of Biostatistics, Harvard T. H. Chan School of Public Health, Boston, MA, U.S.A.* email: [email protected]
Summary:
An important goal of environmental health research is to assess the risk posed by mixtures of environmen-tal exposures. Two popular classes of models for mixtures analyses are response-surface methods and exposure-indexmethods. Response-surface methods estimate high-dimensional surfaces and are thus highly flexible but difficult tointerpret. In contrast, exposure-index methods decompose coefficients from a linear model into an overall mixtureeffect and individual index weights; these models yield easily interpretable effect estimates and efficient inferenceswhen model assumptions hold, but, like most parsimonious models, incur bias when these assumptions do nothold. In this paper we propose a Bayesian multiple index model framework that combines the strengths of each,allowing for non-linear and non-additive relationships between exposure indices and a health outcome, while reducingthe dimensionality of the exposure vector and estimating index weights with variable selection. This frameworkcontains response-surface and exposure-index models as special cases, thereby unifying the two analysis strategies.This unification increases the range of models possible for analyzing environmental mixtures and health, allowing oneto select an appropriate analysis from a spectrum of models varying in flexibility and interpretability. In an analysis ofthe association between telomere length and 18 organic pollutants in the National Health and Nutrition ExaminationSurvey (NHANES), the proposed approach fits the data as well as more complex response-surface methods and yieldsmore interpretable results.
Key words:
Environmental health; multiple index models; kernel machine regression; variable selection.
This paper has been submitted for consideration for publication in
Biometrics a r X i v : . [ s t a t . M E ] J a n
1. Introduction
An important goal of environmental health research is to quantify the risk posed by theenvironmental exposures to which humans are exposed. Exposures of interest can includethose from multiple domains, such as chemical stressors (e.g. metals, persistant organicpollutants, or particles) and non-chemical stressors such as psychosocial stress or diet.Research has historically focused on the effects of individual environmental exposures, butthis is unrealistic as humans are not exposed to a single pollutant in isolation. As such,a priority of the National Institute of Environmental Health Sciences (NIEHS), and theenvironmental health sciences in general, is to understand the impact of environmental mixtures (Dominici et al., 2010; Carlin et al., 2013; Taylor et al., 2016).Numerous statistical approaches have been proposed to estimate the health effects ofmixtures. These have been thoroughly reviewed in the literature (Billionnet et al., 2012;Braun et al., 2016; Davalos et al., 2017; Hamra and Buckley, 2018; Gibson et al., 2019;Lazarevic et al., 2019; Gibson et al., 2019; Tanner et al., 2020). When choosing an appropriatemethod, practitioners need to identify the primary scientific question of interest (Gibsonet al., 2019) and then choose the appropriate statistical method for that question. This choiceoften includes deciding between highly flexible methods that are often hard to interpret andmore restrictive methods that yield interpretable results. Two of the most popular mixturemodels in environmental health research, exposure-response surface methodology and indexbased methods, exemplify this trade off.Exposure-response surface methods, such as Bayesian kernel machine regression (BKMR;Bobb et al., 2015, 2018) or Gaussian process regression more generally (Williams and Ras-mussen, 2006; Ferrari and Dunson, 2019), estimate a multi-dimensional exposure-responsesurface non-parametrically. This class of models allows for estimation of non-linear and non-additive relationships between exposures and response. As such, these models can describe a
Biometrics, 000 broad array of exposure-response relationships. In the specific case of BKMR, interpretationis primarily based on posterior analysis of the high-dimensional exposure-response surfaceand interpretation is highly reliant on visualization. When there is a moderate to largenumber of components in the mixture, succinct interpretation can be difficult due to thelarge number of main effect and interaction surfaces. For example, in the analysis of 18exposures considered in this paper, BKMR requires visual inspection of 18 main effectexposure-response functions and 306 pairwise interaction surface plots.In contrast, exposure-index methods analyze the relationship between a response and alinear multi-exposure index. Among the most popular approaches in the environmental healthsciences are weighted quantile sum regression (WQS; Carrico et al., 2015; Renzetti et al.,2019; Colicino et al., 2019) and quantile G-computation (QGC; Keil et al., 2020). Theseapproaches first transform predictors to the quantile scale and then fit a linear model. Theregression coefficients are then decomposed into (1) an “overall index effect” and (2) indexweights indicating the relative contribution of each mixture component to the overall effect.This decomposition eases interpretation because there is a single parameter that characterizesthe overall effect of the mixture, and a set of weights that reflects the relative importanceof each exposure to this overall effect. In addition, model parsimony results in efficientinferences when model assumptions hold. However, these methods are ultimately instancesof linear models and may not be sufficiently flexible to accurately model a more complexexposure-response relationship. Though previous work has incorporated higher order termsor interactions (Keil et al., 2020), one must specify these parametrically, and moreover thiscomes at the cost of the convenient interpretation of index weights.From a statistical perspective, these index based methods can be viewed as a single indexmodel with an assumed linear association between the outcome and an index formed as aweighted average of quantiled exposures. More general formulations of the single index model (SIM) relax the linearity assumptions and model the relationship between the response andthe exposure index non-parametrically (Powell et al., 1989; Hardle et al., 1993; Naik and Tsai,2001; Hristache et al., 2001; Yu and Ruppert, 2002; Xia et al., 2004; Lin and Kulasekera,2007; Kong and Xia, 2007; Xia, 2008; Wang and Yang, 2009; Wang et al., 2010). In situationswhere there are natural groupings of components within a mixture, one can further allow formultiple indices and interactions among them in a multiple index model (MIM) framework(Stoker, 1986; Ichimura and Lee, 1991; Samarov, 1993; Hristache et al., 2001; Horowitz, 2012;Xia, 2008). These approaches have received substantial attention in the econometric andstatistical literature but have not enjoyed the same popularity in environmental health: allof eight recent reviews of methods for analyzing multi-pollutant mixtures cited WQS/QGCand BKMR, but none considered MIMs (Taylor et al., 2016; Braun et al., 2016; Davaloset al., 2017; Hamra and Buckley, 2018; Gibson et al., 2019; Lazarevic et al., 2019; Gibsonet al., 2019; Tanner et al., 2020). Recently Wang et al. (2020) applied a frequentist SIM forexposure mixtures, but we are unaware of any application of MIMs to mixtures analysis.In terms of methodological innovation, while Bayesian methods for fitting SIMs have beendeveloped (Antoniadis et al., 2004; Park et al., 2005; Choi et al., 2011; Gramacy and Lian,2012; Alquier and Biau, 2013; Reich et al., 2011; Poon and Wang, 2013), we are not awareof Bayesian methods for MIMs.In this paper, we propose a Bayesian multiple index model (BMIM), which provides aunified framework for estimating the association between exposure mixtures and a healthoutcome. Specifically, the BMIM represents a class of models that includes both responsesurface and index-based models. This class of models also encompasses a spectrum of modelsbetween these two extremes that vary in terms of model flexibility and interpretability, thusgiving analysts greater control in selecting an appropriate method for any given data analysis.The proposed approach: (1) facilitates formal Bayesian inference on interpretable index
Biometrics, 000 weights; (2) allows for a non-linear relation between an index and the outcome; (3) per-mits non-additive interactions among multi-exposure indices; and (4) incorporates Bayesianvariable selection on mixture components. Moreover, we extend the BMIM framework toaccommodate binary responses as well as cluster-correlated data.In Section 2 we introduce the motivating case study based on data from the National Healthand Nutrition Examination Survey (NHANES). We briefly outline two of the most popularmethods for mixtures analyses in the environmental health literature—QGC and BKMR—in Section 3. In Section 4 we describe the proposed Bayesian multiple index framework. Wethen compare methods via simulation study (Section 5) and in the NHANES case study(Section 6). Finally, we conclude with a discussion in Section 7.
2. NHANES Case Study: Telemere Length and a Mixture of Pollutants
We consider a case study of the association between a mixture of environmental pollutantsand leukocyte telomere length (LTL), an important biomarker of cellular aging (Mitro et al.,2016). To investigate this question, Mitro et al. (2016), and later Gibson et al. (2019),analyzed data from the 2001–2002 National Health and Nutrition Examination Survey(NHANES). From the original cohort of 11,039, both authors analyzed a subset of 1003people who were at least 20 years of age and had complete data for a number of variables.Following Gibson et al. (2019), we consider an exposure mixture including 18 persistentorganic pollutants that have been grouped naturally into sets of pollutants hypothesized toact similarly: Group 1 includes eight non-dioxin-like polychlorinated biphenyls (PCBs 74,99, 138, 153, 170, 180, 187, and 194); Group 2 includes two non-ortho PCBs (PCBs 126 and169); Group 3 includes PCB 118, four furans, and three dioxins. The outcome of interest islog-LTL, which is thought to be susceptible to environmental exposures (Mitro et al., 2016).Further details on exposure and outcome measurement, as well as on inclusion criteria, havebeen reported in depth elsewhere (Mitro et al., 2016; Gibson et al., 2019).
In what follows, we analyze the same sample of N =1003 observations with the goal ofcharacterizing the relationship between log-LTL and the mixture of pollutants as well as theindividual pollutant contributions.
3. Existing Methods
First let y i be the outcome of interest (log-LTL), { x i , · · · , x iP } be the exposures ( P =18pollutants), and z i be a vector of covariates for the i th observation ( i = 1 , · · · , N ).3.1 Quantile G-Computation (QGC)
Exposure index methods such as WQS (Gennings et al., 2010; Carrico et al., 2015) and QGC(Keil et al., 2020) fit a linear model as follows: y i = α + β ( b q i + · · · + b P q iP ) + z i T γ + (cid:15) i , (1)where q ip represent pre-transformed versions of x ip , scored as quantiles (e.g., 0,1,2,3 if x ip is inthe 0-1st, 1-2nd, 2-3rd or 3-4th quartile). The parameter β represents the linear associationbetween the index and the outcome. Rather than constructing an index based on weightsfixed a priori , WQS and QGC regression simultaneously estimate an overall index associationand component weights.To identify both the weights and the overall effects some identifiability constraint is needed.For WQS, (cid:80) j b j = 1 and b j (cid:62) ∀ j (which we do not pursue here). QGC relaxes thepositivity assumptions and allows for both positive and negative associations. Ultimately,the QGC approach estimates a linear model and then decomposes the estimated coefficientsinto an overall effect beta and component weights b p (cid:48) = β p (cid:48) (cid:80) p β p post-hoc. When estimates are ofopposite sign, QGC typically reports positive weights ( b + p (cid:48) = β p (cid:48) β p (cid:48) > / (cid:80) p β p β p (cid:48) > b − p (cid:48) = β p (cid:48) β p (cid:48) < / (cid:80) p β p β p (cid:48) < qgcomp package; Keil et al., 2020). As such, each is interpreted as the “proportion of the positiveassociation,” and analogously so for the negative weights (Keil et al., 2020). Biometrics, 000
Bayesian Kernel Machine Regression (BKMR)
We consider BKMR as a popular response surface method for environmental mixtures.BKMR is a flexible approach to modelling mixtures that allows non-linear associations andnon-additive interactions among exposures. The BKMR model is y i = h P ( x i , . . . , x iP ) + z Ti γ + (cid:15) i , (2)where h P ( · ) : R P → R is an unknown and potentially non-linear function represented via akernel. We assume h P ( · ) exists in H K , defined by a positive semidefinite reproducing kernel K : R P × R P → R .The choice of kernel function K ( · , · ) uniquely determines a set of basis functions (Cristianiniet al., 2000). A common choice is the Gaussian kernel, K ( x , x (cid:48) ) = exp (cid:104) − (cid:80) Pp =1 ρ p ( x p − x (cid:48) p ) (cid:105) , where ρ p (cid:62) x = ( x , . . . , x P ). This corresponds to radial basisfunctions (Liu et al., 2007).Estimation can proceed based on an equivalence with a linear mixed effects model (Liuet al., 2007). Under a Kernel representation for h P ( · ), model (2) can be written as y i | h i ∼ N ( h i + z Ti γ , σ ) , ( h , · · · , h N ) T ∼ N ( , λ − σ K ) , where K is the kernel matrix with elements K ij = K ( x i , x j ), and λ > λ favoring a more flexible model (Liu et al.,2007). Estimation is based on the marginal likelihood of y = ( y , · · · , y N ) T with respect to h = ( h , · · · , h N ) T : y ∼ N (cid:2) Z γ , σ ( I + λ − K ) (cid:3) . (3)The model is completed by specifying priors for { ρ , γ , σ , λ } . In our simulations and dataanalysis we adopted the default priors proposed by Bobb et al. (2015, 2018).Component-wise variable selection is incorporated via spike and slab priors on the ρ p . Toidentify effects of exposures that may be highly correlated, Bobb et al. (2015) also introduced a hierarchical variable selection scheme. This permits only one pollutant among a group ofexposures to enter the model at a time, which may be too restrictive in some situations.
4. Proposed Approach: Bayesian Multiple Index Models
Model Framework
We propose a Bayesian multiple index model (BMIM) framework to combine the flexibilityof response surface methods with the interpretability of more parsimonious index models.Suppose x i , · · · , x iP can be partitioned into M ( M ∈ { , . . . , P } ) mutually exclusive groupsdenoted x im = ( x im , · · · , x imL m ) T for m = 1 , . . . , M . The proposed model can be written y i = h M (cid:0) x Ti θ , · · · , x TiM θ M (cid:1) + z Ti γ + (cid:15) i , (4)where θ m are L m -vectors of index weights subject to the identifiability constraints: L m T θ m (cid:62)
0, where L m is the unit vector of length L m , and θ Tm θ m = 1 for m = 1 , · · · , M. Contrastthese constraints with those of the linear index models: like QGC, this allows for positiveand negative associations, but rather than summing to 1 the weights have norm 1.The key notion is that, in contrast to BKMR, one need only estimate an exposure-responsesurface of dimension M (cid:54) P , which may remain small even when P is large. Moreover, withinthe m th index, we can interpret the contributions of each exposure via index weights θ m .We again employ kernel function K ( · , · ) that is now a function of a vector of the M indices.That is, the Gaussian kernel can be written as K ( E , E (cid:48) ) = exp (cid:34) − M (cid:88) m =1 ρ m { ( x m − x (cid:48) m ) T θ m } (cid:35) , (5)for E = ( x T θ , · · · , x TM θ M ) and E (cid:48) = ( x (cid:48) T θ , · · · , x (cid:48) TM θ M ).To reduce the computational burden of sampling a vector from this constrained space inour MCMC algorithm, we reparameterize θ ∗ m = ρ / m θ m as in Wilson et al. (2020). Thekernel in (5) is then: K ( E , E (cid:48) ) = exp (cid:34) − M (cid:88) m =1 { ( x m − x (cid:48) m ) T θ ∗ m } (cid:35) . (6) Biometrics, 000
The weights θ m can be fully identified and recovered from the posterior sample of θ ∗ m .Specifically, we have ρ m = || θ ∗ m || = θ ∗ mT θ ∗ m and θ m = || θ ∗ m || − θ ∗ m .4.2 Prior Specification and Posterior Inference
We specify a prior directly for the unconstrained θ ∗ m . In particular, we allow for variableselection on the exposure component weights via spike and slab priors: P ( θ ∗ ml | δ ml ) = δ ml N (0 , σ θ ) + (1 − δ ml ) P , for l = 1 , · · · , L m s.t. ( TL m θ ∗ m (cid:62) ,P ( δ ml ) = Bernoulli ( π ) ,P ( π ) = Beta ( a , b ) , where P is a point mass at zero. Note that this applies selection at the component level, butcan exert selection on an index when none of its components are selected into the model.Although we impose equal shrinkage at the component level (i.e., each component has equalprior inclusion probability), one could entertain other specifications to, say, apply equalshrinkage to each index .To complete the model specification, we adopt a flat prior for γ , σ − ∼ Gamma( a σ , b σ )and λ − ∼ Gamma( a λ , b λ ), and we set the hyperparameters to a σ = 0 . , b σ = 0 . , a λ =1 , b λ = 0 . , a = 1, and b = 1.Inference follows from the marginal likelihood of y , which takes the same form as (3) withthe relevant kernel matrix. After marginalizing over π , the posterior can be decomposed as: P ( θ ∗ , · · · , θ ∗ M , δ , γ , σ − , λ − | y ) = P ( γ | θ ∗ , · · · , θ ∗ M , δ , σ − , λ − , y ) × P ( σ − | θ ∗ , · · · , θ ∗ M , δ , λ − , y ) × P ( θ ∗ , · · · , θ ∗ M , δ , λ − | y ) . We draw from the posterior via MCMC (see Appendix for details). To avoid sampling σ − and γ at every iteration, we integrate over them and draw from the marginal poste-rior P ( θ ∗ , · · · , θ ∗ M , δ , λ − | y ) via modified Gibbs with Metropolis steps, iterating between drawing λ − and drawing ( θ ∗ , · · · , θ ∗ M , δ ) jointly (as in Bobb et al., 2015). After thinningand applying burn-in, we draw σ − directly from a Gamma distribution and γ from anormal distribution. Finally, we obtain posterior samples of the interpretable weights θ m by decomposing draws of θ ∗ m as described above.We summarize the posterior for the weights as follows. We first report the posteriorinclusion probability (PIP) for an entire index—i.e. the posterior probability that ρ m (cid:54) =0.When a whole index is selected out of the model ( ρ m = 0), the component weights θ ml areundefined. In this case, we subsequently describe the conditional posterior for θ ml | ρ (cid:54) = 0via PIP, mean and 95% credible intervals. To summarize component-wise variable selection,we recover marginal PIPs for the component weights by multipling the PIP for ρ m by theconditional PIP for θ ml , as P ( θ ml (cid:54) = 0) = P ( θ ml (cid:54) = 0 | ρ m (cid:54) = 0) P ( ρ m (cid:54) = 0).4.3 Estimating the Exposure-Response Surface
Posterior draws for h i at the observed exposure levels can be obtained following Bobb et al.(2015). Details can be found in the supplementary material. More often, we are interestedin predicting h new on a grid of G new exposure levels. Let E new = ( E new T , . . . , E newG T ) T be the G × M matrix of new levels. Recall that K is the N × N matrix with elements K ij = K ( E i , E j ). Let K nn be the G × G kernel matrix for the new levels, with elements K nnij = K ( E newi , E newj ). Similarly let K no be the G × N kernel matrix comparing new exposure levelsto observed ones, with elements K noij = K ( E newi , E j ). The posterior predictive distribution is h new ∼ M V N (cid:2) λ − K no ( I + λ − K ) − ( y − Z γ ) , σ λ − { K nn − λ − K no ( I + λ − K ) − K noT } (cid:3) . Indexwise Curves and Other Associations of Interest
As with BKMR, we can predict the response surface at arbitrary exposure levels. As such,BMIM estimation yields analogous component-wise response curves component-wise re-sponse curves by varying a single exposure along a grid and holding all others fixed. Thisapproach also allows one to report an “overall mixture effect” analogous to that of QGC Biometrics, 000 by simultaneously increasing all exposures by a quantile, or an “overall index effect” bysimultaneously increasing all exposures within an index.The BMIM structure lends itself to reporting index-wise response curves. Consider the m th index, and set E newgm (cid:48) to some fixed quantile of posterior means of (cid:0) x Tim (cid:48) θ m (cid:48) (cid:1) for each m (cid:48) (cid:54) = m and g = 1 , . . . , G . Then set ( E new m , . . . , E newGm ) T (i.e., the m th column of E new ) toa grid of constants—e.g., equally spaced values between 5th and 95th percentiles across N observations of the posterior means for each x imT θ m . This is a convenient choice inthat it naturally captures how exposures vary jointly. In contrast, increasing exposuressimultaneously by a quantile would not necessarily capture their joint variability unlessexposures are highly correlated. Note that this implicitly takes the weights as fixed and thusisolates uncertainty in the shape of the index-response curve. That is, it ignores uncertaintyin the weights that comprise the index of interest. This is akin to making inferences about β in model (1), or constructing indices via (fixed) toxic equivalency factors as is common intoxicology (e.g., Mitro et al., 2016).The parsimonious set of index-wise exposure-response functions simplifies the interpreta-tion of a fitted model, as one can plot M index-wise curves rather than P component-wisecurves. The contrast is even starker for interactions, as one could present M × ( M − P × ( P − Relation to Existing Methods and a Spectrum of Models
The proposed framework contains two special cases worth highlighting. First, by setting M = 1, the BMIM reduces to a single index model. In particular, if one were to adopt apolynomial kernel: K ( E , E (cid:48) ) = (cid:104) (cid:80) Mm =1 ρ m ( x Tm θ m )( x (cid:48) Tm θ m ) (cid:105) d , of degree d = 1 (and pre- transform exposures accordingly) one could recover the index model in (1)—with the benefitof uncertainty quantification on the weights θ ml . Even in more flexible single-index specifi-cations, say with a higher order polynomial or a Gaussian kernel, one can still estimate well-defined and interpretable weights. Hence, the class of index models is contained in BMIM,including standard linear models and more flexible models with a nonlinear association.At the other extreme, setting M = P —that is, a BMIM model containing P single-exposureindices—corresponds to the ususal BKMR. This shows that the two common approaches toanalyzing environmental mixtures in fact lie on a opposite ends of a spectrum of modelsthat vary in their flexibility. Rather than forcing analysts to choose between one of thesetwo extremes, the BMIM framework permits one to specify the level of flexibility andinterpretability most appropriate for a given analysis.4.6 Software and Extensions
In addition to the methods described here—including both Gaussian and polynomial kernels—we have extended the BMIM to several common scenarios. When outcomes are cluster-correlated, the BMIM can incorporate random intercepts. It can also accommodate binaryoutcomes via a probit-latent variable formulation. Details of these extensions can be found inthe supplemental material. Associated software is available at github.com/glenmcgee/BMIM .
5. Simulations
Setup
We conducted a series of simulations to compare the behaviour of Bayesian single- andmultiple-index models to the more flexible BKMR and the parsimonious QGC. Using thereal exposures ( P =18) and covariates in the NHANES sample, we generated outcomes from: y i ∼ N ( h i + z Ti γ , σ ) , Biometrics, 000 where z i included age (standardized), age , male (0,1), and indicators of BMI (25—29.9;30+). We set γ = ( − . , . , − . , . , . T to loosely reflect effect sizes in the dataapplication, and we set σ = 0 .
5. Finally, we assumed two different structures for h i = h ∗ ( x i ):(A) Scenario A follows from a single index model. First generate x ∗ i = x Ti w where w =( w T , (0 , T , w T ) T such that w = (8 , , , , , , , T and w = − × , standardizedto have mean 0 and standard deviation 1. Then h i = h A ( x ∗ i ), where h A ( · ) is S-shaped.(B) Scenario B follows from a multiple index model. Generate x ∗ = ( x , . . . , x ) w , x ∗ =( x , x ) w and x ∗ = ( x , . . . , x ) w , where w = (1 , T , standardized as above, and w and w are as defined above. Then h i = h ∗ B ( x ∗ i ) + h ∗ B ( x ∗ i ) + h ∗ B ( x ∗ i ) + 0 . h ∗ B ( x ∗ i ) × h ∗ B ( x ∗ i )where h B is a unimodal function, h B is flat (null), and h B is sigmoidal.Details on the exposure response curves { h A , h B , h B , h B } can be found in the Appendix.In each setting, we generated 100 datasets of 500 observations. We then split each sampleinto a training set of N = 300 observations, and held out the remaining 200 observations onwhich to evaluate model fit. To each training set, we then fit QGC (with q =10, i.e. deciles),a Bayesian single index model (BSIM), a 3-index BMIM using the same indices as above(BMIM-3), full BKMR, and BKMR using hierarchical variable selection using the threegroups of exposures (BKMR-H), each with Gaussian kernel.Intuitively, in scenario (A), we would expect the BSIM to perform best as it is correctlyspecified (with the shape of the response curve unknown). We expect BMIM (with M = 3)and BKMR to be flexible enough to still perform well here, although with more variability dueto less model structure. In scenario (B), the BSIM is mis-specified and hence should performpoorly. Here BMIM is correctly specified (with unknown response curve) and should performbetter than the more flexible BKMR due to its increased structure. QGC is misspecified inboth scenarios, as it assumes linearity on the quantile scale. We evaluated model fit via mean squared error (MSE) for the unknown h i in the testset, as well as MSE for the outcomes in the test set and in 4-fold cross validation (CV).For MSE, we report the mean as well as the standard deviation across datasets in order tocontextualize differences between models. We also report 95% interval coverage (Cvg) andaverage standard errors (SE) for h i in the test set.5.2 Simulation Results
Under the single index data-generating mechanism (A), the BSIM had the lowest MSE. Thenext lowest MSE was for BMIM-3 followed by BKMR. Under the three-index data-generatingmechanism (B), BMIM-3 had the lowest MSE followed by BKMR. The incorrectly specifiedBSIM performed worse than both. In both scenarios, QGC exhibited very high MSE relativeto the other approaches, because it was misspecified. As described previously, one couldincorporate higher order terms into the QGC approach, but we found that incorporatingquadratic terms for all exposures still did not perform well (see supplementary material).Cross validation MSE followed the same pattern, suggesting this can be used for modelselection in practice.Results for interval coverage behaved similarly, with the simplest correctly specified modelachieving near nominal coverage. Interestingly, under the single-index scenario (A), the moreflexible alternatives (BMIM-3 and BKMR) had progressively lower coverage, potentiallybecause of bias due to applying more shrinkage; nevertheless both greatly outperformedQGC. In scenario (B), BMIM-3 and BKMR both achieved near nominal coverage whereasthe misspecified BSIM and QGC had very low coverage. Naturally, the more structuredBSIM had smaller SEs than BMIM-3, which had lower SEs than BKMR in (A) and (B).Interestingly, when we repeated the simulations in low-signal settings (setting σ =1.0 or2.0), the simpler models were no longer the clear winners. Instead, the BSIM, BMIM-3 andBKMR all performed about the same in terms of MSE, interval coverage, and even average Biometrics, 000
SEs. Indeed, cross validation MSEs were effectively equal in low-signal settings. These modelsfit about equally well, and it would be difficult to choose one that best fits the data. Facingsuch a setting in practice, one would likely favor the model with simpler interpretations.Note that we report here results for BKMR using the same priors for ρ m as in the BMIMand BSIM. In the supplementary material we report results for BKMR using the defaultpriors for ρ m from the bkmr package Bobb et al. (2018), which performed very similarly toour implementation when σ was small but worse when σ was high. BKMR with hierarchicalvariable selection generally performed worse in terms of MSE and interval coverage butachieved the lowest average SEs (see supplementary material), because it imposes the mostshrinkage by limiting a maximum of one exposure per group to feature in the model at atime. Therefore it is misspecified for the data generation mechanisms considered here.We also repeated scenario A with a linear response curve, and the results are similar to thosereported here (see Appendix). Even here, QGC remains misspecified, as it assumes linearityon the quantile scale. As such it still performed poorly relative to the kernel approaches.Unsurprisingly, BSIM with a polynomial kernel of degree one performed best.[Table 1 about here.]
6. Analysis of NHANES Case Study
Models and Analyses
We analyzed the NHANES data using the methods described here in order to compare andcontrast the different conclusions that can be drawn under each approach. We considereda broad set of analytic goals: characterizing the overall outcome association of the entiremixture, quantifying the relative contributions of each mixture component, characterizingthe outcome associations of each index, and investigating two-way interactions.We conducted a QGC analysis (with exposures binned into deciles; Q = 10), full BKMR (with componentwise as well as hierarchical variable selection via spike-and-slab priors), aBSIM analysis, and a BMIM analysis with three indices (BMIM-3; M = 3) correspondingto the three exposure groups identified by Gibson et al. (2019), within which exposures areexpected to act similarly. In the kernel-based methods, exposures were first standardized tohave mean 0 and standard deviation 1, and we adopted a Gaussian kernel. All models wereadjusted for age (linear and quadratic), sex, and BMI ( <
25, 25 − (cid:62) ResultsOverall Mixture Association.
All of the methods can estimate an “overall” association be-tween mixture and health. However, the interpretations differ somewhat. QGC quantifies theoverall association through its β parameter. This parameter represents the mean differencein log(LTL) comparing a population to another with all exposures one decile group higher.We estimate the overall effect to be 0.069 (95% confidence interval: [0.030, 0.108]). Thekernel methods can estimate arbitrary contrasts, so we estimated an analogous “overall”mixture-health association, this time comparing a population with all exposures at their60th percentile versus all exposures set to their 50th percentile. Using this approach, thethree kernel methods all produced similar estimates of the overall effect: the BSIM estimateof this was 0.034 (95% credible interval [-0.006, 0.075]), the BMIM-3 estimate was 0.037(95%CI [-0.004, 0.077]), and the BKMR estimate was 0.037 (95% CI [0.005, 0.069]). Thereare two likely reason that the estimates of the overall effect vary between QGC and thekernel-based methods. First, because of the non-linearity in the kernel approaches (as theestimate may depend on the specific quantiles being compared), but also because—despitetheir similar interpretations—the estimands are slightly different, as QGC compares twoquantile categories rather than two specific exposure values. Individual Component Contributions.
All of the methods provide a measure of variable orcomponent importance. The interpretation and inference varies between methods. Biometrics, 000
QGC provides two sets of weights—positive and negative weights—and each represents theproportion of the positive or negative effect that can be attributed to that exposure. Table 2presents estimates of these weights. Among the ten exposures with positive weights, Furan1contributed the most, with a weight 0.18, and PCB99 and Furan4 each had weights of 0.14.Among the negative weights, PCB180 contributed the most, with a weight of 0.29, and theothers were no more than 0.16. A key limitation is that the current software provides noestimates of uncertainty for the weights. Hence, inferential statements about the relative sizeof weights—or even the sign of weights—cannot be readily made.BKMR does not estimate exposure weights. Instead, we can get a sense of variable im-portance via posterior inclusion probabilities (PIPs) for the exposures (Table 2). This isthe posterior probability that ρ m (cid:54) = 0. The PIP identifies the probability that a particularcomponents contributes to the exposure-response function but does not provide a measureof effect size. This is in sharp contrast to the weights in QGC that identify only the relativesize of the association but not statistical certainty. BKMR also identified Furan1 as havingby far the strongest signal, with a PIP of 0.88; all other exposures had PIPs below 0.30.PCB180, which had the largest component weight in the QGC analysis, had a PIP of 0.13indicating weak support for an association with the outcome.The BSIM and BMIM-3 estimate PIPs as well as component weights, along with 95%credible intervals for the weights. The BSIM estimates of the weights for components inthe first group (non-dioxin-like PCBs) were all close to zero (between -0.04 and 0.06; PIPsbetween 0.14 and 0.17). In contrast, two exposures in this group—PCB99 and PCB180—had among the strongest weights under QGC. Among the 10 other exposures, the signs ofthe weights matched those of QGC for all but Furan3. Overall, the results were much moresimilar to those of BKMR, in that Furan1 had by far the strongest weight, with a PIP of θ =0.98, 95% CI: [-0.12, 1.00]). In index 1, the estimated weight for PCB74(0.39; 95% CI [-0.49, 1.00]) was nearly twice that of PCB187 (0.21; 95% CI [-0.50, 0.99]),indicating the association was twice as strong for PCB74 in the direction of the responsecurve for index 1—however, index 1 as a whole was very weakly associated with the outcome.[Table 2 about here.]As Furan1 was identified as the strongest mixture component by each of the kernel meth-ods, we plot the corresponding estimated exposure-response curves (with all other exposuresset to their medians) under each of these approaches in Figure 1. The three plots are nearlyidentical, indicating that BKMR, the BSIM and the BMIM-3 lead to roughly the sameconclusions for the association between Furan1 and the outcome.[Figure 1 about here.] Indexwise Exposure-Response Curves.
One of the advantages of the proposed index modelsis that they naturally facilitate studying response curves at the index level, rather thanreporting P = 18 individual exposure-response curves, or increasing every exposure simul-taneously. Under the BSIM, we visualize the entire index-wise estimated response curve inFigure 2 (a), which appears just slightly sub-linear. Under the BMIM-3, the three estimatedindex-wise exposure response curves are not radically different from that of the BSIM (panelb), with each increasing sublinearly. Specifically, the estimated response curve for index 1is close to null (flat) and has the widest credible intervals; that of index 2 corresponds toa slightly stronger association, and that of index 3 is the strongest. Note in particular thesimilarity between the response curve for index 3 and the component-wise curve for Furan1. Biometrics, 000
Matching this observation, the index-level PIPs (for ρ m ) were correspondingly low for thefirst two indices (0.60 and 0.44), whereas that of index 3 was 0.97.[Figure 2 about here.] Interactions.
A feature of the BMIM and BKMR approaches is that we can considerinteractions, although we consider a different set of interactions in each. Under BKMR,we compare fitted exposure response curves for x j when x j (cid:48) is set to its 10 th , 50 th or 90 th percentile (and all others to their medians). We display these for all pairs j (cid:54) = j (cid:48) in Figure3 (a), although it is not straightforward to interpret the 306 ( P × [ P − indices , so that rather than scanning 306 plots, we need only investigate 6 ( M × [ M − Model Fit.
We also conducted 5-fold cross validation in order to compare model fit via MSE,as in the simulations. Of the four models, QGC had the highest MSE (0.790), and BKMRhad the lowest (0.774). From a statistical perspective, the MSEs for BSIM and BMIM-3were effectively equal to that of BKMR (less than 0.3% higher). As the results of the kernelapproaches identified one strong association, we also fit BKMR with hierarchical variableselection, but it fit no better than the standard BKMR (0.774). Practically speaking, thesimpler index models were able to achieve virtually the same fit as the more complex BKMR,all while simplifying the burden of interpretation substantially.
7. Discussion
The proposed BMIM framework represents a spectrum of models, with linear index modelsand BKMR occupying either end of the spectrum. At present it is common in environmentalmixture studies to adopt one of these two extremes. A key goal of this paper is to increasethe range of options available to an analyst and unify the two seemingly unrelated methods.By bridging the gap between the two, this framework gives analysts the freedom to select ananalysis strategy with an appropriate balance of flexibility and interpretability. At its mostparsimonious, the BMIM extends the standard linear index models to allow for variableselection and non-linear relations. When
M >
1, and in the absence of interactions, it canalso allow for additive index models, which have not been explored in the environmentalhealth literature. Even when a single multi-exposure index is believed to be appropriate,the BMIM allows one to examine interactions between this multi-exposure index and oneor more individual covariates like age or a measure of socioeconomic status (SES), e.g., y i = h (cid:0) x i T θ , age i , SES i (cid:1) + z Ti γ + (cid:15) i .Unlike fully non-parametric approaches, the BMIM framework allows one to incorporateand investigate biologically plausible mechanisms. If, say, a set of compounds is hypothesizedto bind to the same receptor, an index structure might allow one to easily encode this priorknowledge. In toxicology, it is common to construct multi-pollutant indices based on relativepotency factors (e.g. toxic equivalency factors) derived from laboratory experiments (Mitroet al., 2016). The BMIM allows one to incorporate such indices while allowing for a non-linearrelationship and allowing one to investigate possible interactions among indices without fullyspecifying a parametric relationship.In the NHANES case study, some commonalities and a few interesting differences emergedacross methods. All methods weighted the Furan1 association highly. The BSIM, BMIM-3,and BKMR all indicated it had the strongest association with the outcome; QGC estimated Biometrics, 000 it to have the largest positive weight. An interesting difference is that QGC reported aneven larger negative weight for PCB180, whereas the response-surface and Bayesian indexmethods did not weight this exposure highly (based on either posterior inclusion probabilitiesor component weights). This occurs because while the estimated coefficient for this exposurein the underlying linear model is the largest in magnitude, it is not large relative to theuncertainty of the estimate. The Bayesian methods reflect this uncertainty, and therefore thestrength of evidence of an association with a particular exposure, whereas the conventionfor existing index models in the literature has been to not report these uncertainties. Ourresults suggest this is worth doing even if one opts for an existing frequentist index method.A feature of the proposed framework is the ability to estimate indexwise associations. Froma substantive perspective, this more naturally captures how exposures within an index varyjointly—rather than artificially setting them all to some fixed quantile simultaneously. Onemight even estimate an overall mixture effect in the same way, say by contrasting all indicesset to two quantiles. From a statistical perspective, this can also help alleviate sparsity issues.For example it is typical to summarize interactions with BKMR by plotting an estimatedexposure response curve at some fixed level of another exposure and at median levels of allothers. This might be reasonable for a subset of the exposures, but one quickly encounterssparsity when examining the 153 different pairwise comparisons in the NHANES data. Thisis exacerbated when one considers higher order interactions. By instead basing inference oncomposites of multiple exposures, the BMIM is less reliant on such fine-grained comparisons.
Acknowledgements
This research was supported by NIH grants ES028800, ES028811, ES028688, and ES000002.This research was also supported by USEPA grants RD-835872-01 and RD-839278. Its con-tents are solely the responsibility of the grantee and do not necessarily represent the officialviews of the USEPA. Further, USEPA does not endorse the purchase of any commercialproducts or services mentioned in the publication. Supporting Information
Supplementary material is available online. Software is available at github.com/glenmcgee/BMIM . References
Alquier, P. and Biau, G. (2013). Sparse single-index model.
Journal of Machine LearningResearch
Statistica Sinica pages 1147–1164.Billionnet, C., Sherrill, D., Annesi-Maesano, I., et al. (2012). Estimating the health effectsof exposure to multi-pollutant mixture.
Annals of epidemiology
Environmental Health
Biostatistics
Environmentalhealth perspectives
A6–A9.Carlin, D. J., Rider, C. V., Woychik, R., and Birnbaum, L. S. (2013). Unraveling the healtheffects of environmental mixtures: an niehs priority.Carrico, C., Gennings, C., Wheeler, D. C., and Factor-Litvak, P. (2015). Characterizationof weighted quantile sum regression for highly correlated data in a risk analysis setting.
Journal of agricultural, biological, and environmental statistics
Journal of Nonparametric Statistics Biometrics, 000
Colicino, E., Pedretti, N. F., Busgang, S., and Gennings, C. (2019). Per-and poly-fluoroalkylsubstances and bone mineral density: results from the bayesian weighted quantile sumregression. medRxiv page 19010710.Cristianini, N., Shawe-Taylor, J., et al. (2000).
An introduction to support vector machinesand other kernel-based learning methods . Cambridge university press.Davalos, A. D., Luben, T. J., Herring, A. H., and Sacks, J. D. (2017). Current approaches usedin epidemiologic studies to examine short-term multipollutant air pollution exposures.
Annals of epidemiology
Epidemiology(Cambridge, Mass.) arXiv preprint arXiv:1911.01910 .Gennings, C., Sabo, R., and Carney, E. (2010). Identifying subsets of complex mixtures mostassociated with complex diseases: polychlorinated biphenyls and endometriosis as a casestudy.
Epidemiology pages S77–S84.Gibson, E. A., Goldsmith, J., and Kioumourtzoglou, M.-A. (2019). Complex mixtures,complex analyses: an emphasis on interpretable results.
Current Environmental HealthReports Environmental Health computer experiments. Technometrics
Current epidemiology reports The annals of Statistics pages 157–178.Horowitz, J. L. (2012).
Semiparametric methods in econometrics , volume 131. SpringerScience & Business Media.Hristache, M., Juditsky, A., Polzehl, J., Spokoiny, V., et al. (2001). Structure adaptiveapproach for dimension reduction.
The Annals of Statistics
Annals of Statistics pages 595–623.Ichimura, H. and Lee, L.-F. (1991). Semiparametric least squares estimation of multipleindex models: single equation estimation.
Nonparametric and semiparametric methodsin econometrics and statistics pages 3–49.Keil, A. P., Buckley, J. P., O’Brien, K. M., Ferguson, K. K., Zhao, S., and White, A. J.(2020). A quantile-based g-computation approach to addressing the effects of exposuremixtures.
Environmental health perspectives
Biometrika
Environmental health perspectives
Biometrika Biometrics, 000
Liu, D., Lin, X., and Ghosh, D. (2007). Semiparametric regression of multidimensional ge-netic pathway data: Least-squares kernel machines and linear mixed models.
Biometrics
Environmental health perspectives
Biometrika
Journal of Computational and Graphical Statistics
Computational statistics & data analysis
Econometrica: Journal of the Econometric Society pages 1403–1430.Reich, B. J., Bondell, H. D., and Li, L. (2011). Sufficient dimension reduction via bayesianmixture modeling.
Biometrics
Journal of Statistical Software .Samarov, A. M. (1993). Exploring regression structure using nonparametric functionalestimation.
Journal of the American Statistical Association
Econometrica: Journal ofthe Econometric Society pages 1461–1481.Tanner, E., Lee, A., and Colicino, E. (2020). Environmental mixtures and children’s health:identifying appropriate statistical approaches.
Current Opinion in Pediatrics Taylor, K. W., Joubert, B. R., Braun, J. M., Dilworth, C., Gennings, C., Hauser, R., Heindel,J. J., Rider, C. V., Webster, T. F., and Carlin, D. J. (2016). Statistical approaches forassessing health effects of environmental chemical mixtures in epidemiology: lessons froman innovative workshop.
Environmental health perspectives
A227–A229.Wang, J.-L., Xue, L., Zhu, L., Chong, Y. S., et al. (2010). Estimation for a partial-linearsingle-index model.
The Annals of statistics
Statistica Sinica pages 765–783.Wang, Y., Wu, Y., Jacobson, M., Lee, M., Jin, P., Trasande, L., and Liu, M. (2020). A familyof partial-linear single-index models for analyzing complex environmental exposures withcontinuous, categorical, time-to-event, and longitudinal health outcomes.
EnvironmentalHealth .Williams, C. K. and Rasmussen, C. E. (2006).
Gaussian processes for machine learning ,volume 2. MIT press Cambridge, MA.Wilson, A., Hsu, H.-H. L., Mathilda Chiu, Y.-H., Wright, R. O., Wright, R. J., and Coull,B. A. (2020). Kernel machine and distributed lag models for assessing windows ofsusceptibility to mixtures of time-varying environmental exposures in children’s healthstudies. arXiv preprint arXiv:1904.12417 .Xia, Y. (2008). A multiple-index model and dimension reduction.
Journal of the AmericanStatistical Association
Statistica Sinica pages 1–28.Yu, Y. and Ruppert, D. (2002). Penalized spline estimation for partially linear single-indexmodels.
Journal of the American Statistical Association Biometrics, 000 −0.50−0.250.000.250.500.75 −1 0 1 2
Exposure Component E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Component 15A −0.50−0.250.000.250.500.75 −1 0 1 2
Exposure Component E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Index 1 , Component 15B −0.50−0.250.000.250.500.75 −1 0 1 2
Exposure Component E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Index 3 , Component 5C
Figure 1 : Exposure response curve for Furan1 in NHANES analysis. Panel A shows fittedcurve using the bkmr package in R; B is based on single index model; C is based on 3-indexmodel. Exposure E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Index 1 (a) Bayesian single index model ( M =1; L =18) −0.250.000.250.50 −1 0 1 2 Exposure E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Index 1 −0.250.000.250.50 −1 0 1 2
Exposure E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Index 2 −0.250.000.250.50 0 1
Exposure E s t i m a t ed e x po s u r e − r e s pon s e ( h ) Index 3 (b) Bayesian multiple index model ( M =3; L =8, L =2, L =8) Figure 2 : Estimated index-wise response curves from NHANES analysis. The top panelshows the exposure-response function for the BSIM. The bottom panel shows the exposure-response functions for each of the three indices in the BMIM-3. For the BMIM-3, indexexposure-response function is plotted with the other indices fixed at the median value. Theplot also shows 95% credible intervals. Biometrics, 000 D i o x i n1 D i o x i n2 D i o x i n3 F u r an1 F u r an2 F u r an3 F u r an4 P C B P C B P C B P C B P C B P C B P C B P C B P C B P C B P C B Dioxin1Dioxin2Dioxin3Furan1Furan2Furan3Furan4PCB074PCB099PCB118PCB126PCB138PCB153PCB169PCB170PCB180PCB187PCB194 − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . − . − . . . . Q uan t il e . . . h ( e x po s | quan t il e s o f e x po s ) ( a ) B K M R ( M = P = ) I nde x I nde x I nde x Index 1Index 2Index 3 − − . . . . . . . . . . . . . . . Q uan t il e . . . I nde x w i s e I n t e r a c t i on s ( b ) B a y e s i a n m u l t i p l e i nd e x m o d e l ( M = ; L = , L = , L = ) F i g u r e : V i s u a li z i n g t w o - w a y i n t e r a c t i o n s und e r B K M R ( l e f t) a nd t h e B M I M - (r i g h t) i n t h e NHAN E S a n a l y s i s . E a c hp a n e l s h o w s t h ee x p o s u r e - r e s p o n s e r e l a t i o n f o r o n ee x p o s u r e ( p a n e l[ a ] ) o r o n e i nd e x ( p a n e l[ b ] ) w i t h a ll o t h e r e x p o s u r e s o r i nd i ce s fi x e d a t ag i v e np e r ce n t il e ( . , . , o r . ) . P a r a ll e lli n e s i nd i c a t e n o i n t e r a c t i o nb e t w ee n c o m p o n e n t s o r i nd i ce s w h il e d e v i a t i o n i nd i c a t e s i n t e r a c t i o n s . C o l o r v e r s i o n c a nb e f o und o n li n e i n t h ee l ec tr o n i c v e r s i o n o f t h e a rt i c l e . Table 1: Simulation Results across 100 datasets. The table shows mean squared error (MSE),average standard errors (SE) and 95% interval coverage (Cvg) on the estimated h functionin the test dataset, as well as MSE for y evaluated on a test dataset (MSE( y )) and via crossvalidation (CV-MSE( y )). BKMR-H is BKMR with hierarchical variable selection.MSE( h ) SE( h ) Cvg( h ) MSE( y ) CV-MSE( y ) σ M Model Mean SD Mean SD Mean SD0.5 1 QGC 0.150 0.015 0.19 0.63 0.41 0.04 0.44 0.04BSIM 0.031 0.010 0.17 0.94 0.28 0.03 0.30 0.03BMIM-3 0.051 0.011 0.21 0.91 0.30 0.03 0.34 0.03BKMR 0.076 0.012 0.22 0.83 0.32 0.04 0.37 0.03BKMR-H 0.085 0.014 0.14 0.60 0.33 0.04 0.38 0.033 QGC 0.188 0.024 0.19 0.61 0.41 0.04 0.42 0.04BSIM 0.122 0.027 0.17 0.54 0.34 0.04 0.35 0.04BMIM-3 0.035 0.012 0.19 0.95 0.28 0.03 0.29 0.03BKMR 0.048 0.013 0.23 0.94 0.29 0.03 0.30 0.03BKMR-H 0.085 0.021 0.14 0.62 0.33 0.04 0.35 0.031.0 1 QGC 0.226 0.048 0.33 0.81 1.21 0.13 1.27 0.11BSIM 0.086 0.038 0.30 0.94 1.07 0.12 1.13 0.10BMIM-3 0.113 0.033 0.31 0.90 1.10 0.12 1.17 0.10BKMR 0.125 0.034 0.30 0.87 1.11 0.12 1.19 0.10BKMR-H 0.142 0.035 0.22 0.65 1.13 0.13 1.17 0.093 QGC 0.263 0.059 0.33 0.77 1.21 0.12 1.27 0.12BSIM 0.157 0.042 0.30 0.81 1.13 0.11 1.15 0.10BMIM-3 0.090 0.035 0.33 0.95 1.07 0.11 1.11 0.11BKMR 0.094 0.035 0.34 0.95 1.08 0.11 1.11 0.10BKMR-H 0.141 0.046 0.23 0.71 1.12 0.12 1.15 0.102.0 1 QGC 0.528 0.177 0.62 0.91 4.43 0.47 4.63 0.42BSIM 0.227 0.113 0.48 0.92 4.17 0.45 4.28 0.36BMIM-3 0.229 0.103 0.48 0.92 4.18 0.44 4.33 0.37BKMR 0.232 0.104 0.48 0.92 4.18 0.45 4.34 0.34BKMR-H 0.291 0.120 0.37 0.75 4.24 0.46 4.29 0.353 QGC 0.565 0.186 0.62 0.89 4.43 0.45 4.63 0.44BSIM 0.228 0.096 0.50 0.94 4.18 0.42 4.28 0.39BMIM-3 0.200 0.097 0.52 0.96 4.15 0.42 4.29 0.40BKMR 0.194 0.101 0.52 0.96 4.15 0.42 4.29 0.40BKMR-H* 0.277 0.118 0.38 0.80 4.22 0.44 4.28 0.37 *–based on 99 datasets, due to computational instability. Biometrics, 000 T a b l e : Su mm a r i z i n g e x p o s u r e w e i g h t s i n B K M R a nd t h e M I M . F o r B K M R w e r e p o rt p o s t e r i o r i n c l u s i o np r o b a b ili t i e s ( P I P s ) f o r e a c h e x p o s u r e . F o r M I M , w e r e p o rtt h e P I P f o rt h ee n t i r e i nd e xv i a ρ ; w e a l s o s u mm a r i ze t h e d i s tr i bu t i o n o f w e i g h t s θ c o nd i t i o n a l o n ˆ ρ (cid:54) = ( o t h e r w i s e i t i s n o t w e ll d e fin e d ) . E s t i s t h e p o s t e r i o r m e a n s t a nd a r d i ze d t o s a t i s f y t h ec o n s tr a i n t s ; C I i s % c r e d i b l e i n t e r v a l. Q G C B K M R B S I M ( M = ) B M I M - ( M = ) W e i g h t s ˆ r p ˆ ρ ˆ θ l | ˆ ρ (cid:54) = ρ m ˆ θ m l | ˆ ρ m (cid:54) = G r o up E x p o s u r e P o s N e g P I PP I PP I PE s t C I P I PP I PE s t C I P C B . . . . . ( - . , . ) . . . ( - . , . ) P C B . . . . ( - . , . ) . . ( - . , . ) P C B . . . . ( - . , . ) . . ( - . , . ) P C B . . . . ( - . , . ) . . ( - . , . ) P C B . . . - . ( - . , . ) . . ( - . , . ) P C B . . . - . ( - . , . ) . . ( - . , . ) P C B . . . . ( - . , . ) . . ( - . , . ) P C B . . . - . ( - . , . ) . . ( - . , . ) P C B . . . . ( - . , . ) . . . ( - . , . ) P C B . . . . ( - . , . ) . . ( . , . ) P C B . . . . ( - . , . ) . . . ( - . , . ) D i o x i n . . . . ( - . , . ) . - . ( - . , . ) D i o x i n . . . - . ( - . , . ) . - . ( - . , . ) D i o x i n . . . - . ( - . , . ) . . ( - . , . ) F u r a n . . . . ( . , . ) . . ( - . , . ) F u r a n . . . . ( - . , . ) . . ( - . , . ) F u r a n . . . . ( - . , . ) . . ( - . , . ) F u r a n . . . . ( - . , . ) . . ( - . , .94