{'fraction_non_alphanumeric': 0.025580658305790736, 'fraction_numerical': 0.011454701860823108, 'mean_word_length': 4.59361809893635, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimating the health effects of environmental mixtures using Bayesian semiparametric regression and sparsity inducing priors

Humans are routinely exposed to mixtures of chemical and other environmental factors, making the quantification of health effects associated with environmental mixtures a critical goal for establishing environmental policy sufficiently protective of human health. The quantification of the effects of exposure to an environmental mixture poses several statistical challenges. It is often the case that exposure to multiple pollutants interact with each other to affect an outcome. Further, the exposure-response relationship between an outcome and some exposures, such as some metals, can exhibit complex, nonlinear forms, since some exposures can be beneficial and detrimental at different ranges of exposure. To estimate the health effects of complex mixtures we propose a flexible Bayesian approach that allows exposures to interact with each other and have nonlinear relationships with the outcome. We induce sparsity using multivariate spike and slab priors to determine which exposures are associated with the outcome, and which exposures interact with each other. The proposed approach is interpretable, as we can use the posterior probabilities of inclusion into the model to identify pollutants that interact with each other. We illustrate our approach's ability to estimate complex functions using simulated data, and apply our method to two studies to determine which environmental pollutants adversely affect health.


Introduction
The study of the health impacts of environmental mixtures, including both chemical and non-chemical stressors, on human health is a research priority in the environmental health sciences (Carlin et al. , 2013;Taylor et al. , 2016;Braun et al. , 2016). In the past, the majority of health effect studies upon which environmental regulation is based estimate the independent association between an outcome and a given exposure, either marginally or while controlling for co-exposures. However, exposure to mixtures of chemicals is a real world scenario, and it has been argued that chemical mixture effects are possible even when each chemical exposure is below its no observable adverse effect level (Kortenkamp et al. , 2007).
Accordingly, there has recently been a renewed focus on the development of statistical methods for the analysis of the effects of chemical mixtures, or more generally environmental mixtures defined to include both chemical and non-chemical (such as diet or psychosocial factors) stressors, on health. Recent work on this topic has recognized that investigations on these effects may focus on a variety of research questions of interest, including the identification of individual exposures driving this effect, identification of interactions among chemicals on health, and the estimation of the cumulative health effects of exposure to multiple agents. Recent reviews of relevant statistical methods useful for analyzing the resulting data have highlighted that the choice of statistical method will depend on the primary objective of the study , and recent work in this area has developed both Bayesian and frequentist methods that implement forms of regularization, variable selection, dimension reduction, statistical learning, or smoothing in various combinations to address these questions in the presence of correlated exposures (see Lazarevic et al. (2019) for an excellent review). Notably, there is a gap in this literature in terms of methods that directly quantify the strength of evidence supporting interactions among individual chemicals on a health outcome in the data. For instance, Bayesian kernel machine regression (BKMR; Bobb et al. (2014Bobb et al. ( , 2018) can estimate cumulative health effects of exposure to multiple agents, provide variable importance scores for each individual chemical, and estimate selected summary statistics that reflect interactions among some chemicals at given levels of others. However, it does not yield overall variable importance scores for interactions among chemicals within the larger mixture of exposures.
Other approaches allow for interactions between exposures, perhaps with some variable selection technique applied, but these approaches are often based on models assuming simple, linear associations between exposures or other strong parametric assumptions. It has been recently shown that in some settings, such as in studies on the impacts of exposure to metal mixtures on neurodevelopment in children, the exposureresponse relationship between the multivariate exposure of interest and a health outcome can be a complex, non-linear and non-additive function (Henn et al. , 2010Bobb et al. , 2014). Therefore there is a need for statistical approaches that can handle complex interactions among a large number of exposures. Even in settings with a relatively small number of exposures, the number of potential interactions can become large, making some form of variable selection or dimension reduction necessary. In this paper we provide a Bayesian approach to this problem that allows for nonlinear effects on the outcome, can identify important exposures among a large mixture, and can identify important interactions among the elements of a mixture.

Review of statistical literature
There exists a vast statistical literature on high-dimensional regression models, including approaches to reduce the dimension of the covariate space. The LASSO (Tibshirani, 1996) utilizes an L1 penalty on the absolute value of regression coefficients in a regression model, which effectively selects which variables should be included in the model by setting some coefficients equal to zero. Many approaches have built upon this same methodology (Zou & Hastie, 2005;Zou, 2006); however these all typically rely on relatively simple, parametric models. More recently, a number of papers (Zhao et al. , 2009;Bien et al. , 2013;Hao et al. , 2014;Lim & Hastie, 2015) have adopted these techniques to identify interactions, although these approaches are only applicable to linear models with pairwise interactions. While similar ideas for nonlinear models were used by Radchenko & James (2010), these analyses were restricted to pairwise interactions. Moreover, quantification of estimation uncertainty can be quite challenging in high dimensions.
A number of approaches have been introduced that embed variable selection within nonlinear regression models. Shively et al. (1999) utilized variable selection in Bayesian nonparametric regression based on an integrated Wiener process for each covariate in the model, and then extended these ideas in Wood et al. (2002) to non additive models that allowed for interactions between covariates. They utilized the BIC approximation to the marginal likelihood in order to calculate posterior model probabilities, which requires enumerating all potential models. Even when the number of covariates in the model is not prohibitively large, inclusion of just two-way interactions can lead to a massive number of possible models. Reich et al. (2009) utilized Gaussian process priors for all main effects and two-way interactions within a regression model and performed variable selection in order to assess which covariates are important predictors of the outcome. They do not allow for explicit identification of higher-order interactions, as all remaining interactions are absorbed into a function of the covariates that accounts for higher-order interactions. Yau et al. (2003) and Scheipl et al. (2012) utilized variable selection in an additive regression model framework with interaction terms within a Bayesian framework, though one must specify the terms that enter the model a priori and there can be a large number of potential terms to consider. Qamar & Tokdar (2014) used Gaussian processes to separately model subsets of the covariate space and used stochastic search variable selection to identify which covariates should enter into a given Gaussian process function, implying that they interact with each other. While this approach can be used to capture and identify complex interactions, it relies on Gaussian processes and therefore does not scale to larger data sets.
Complex models that allow for nonlinearities and interactions are becoming popular in the machine learning literature, with Gaussian process (or kernel machine) regression (O'Hagan & Kingman, 1978), support vector machines (Cristianini & Shawe-Taylor, 2000), neural networks (Bengio et al. , 2003), Bayesian additive regression trees (BART, Chipman et al. (2010)), and random forests (Breiman, 2001) being just a few examples. Bobb et al. (2014) used Bayesian kernel machine regression to quantify the health effects of an environmental mixture. This approach imposed sparsity within the Gaussian process to flexibly model multiple exposures. While these approaches and their many extensions can provide accurate predictions in the presence of complex relationships among exposures, they typically only provide a variable importance measure for each feature included in the model, and do not allow for formal inference on higher-order effects. An important question remains how to estimate the joint effect of a set of exposures that directly provides inference on interactions, while scaling to larger data sets without requiring a search over such a large space of potential models involving higher-order interactions. Further, full Markov Chain Monte Carlo (MCMC) inference in kernel machine regression does not scale with the sample size, limiting its applicability in large sample settings. While one can potentially use approximate Bayesian methods, such as variational inference to simplify computation, these methods are only approximate and can yield biased estimators for some parameters (Guo et al. , 2016;Miller et al. , 2016).

Our contribution
In this paper we introduce a Bayesian semiparametric regression approach that addresses several of the aforementioned issues as it allows for nonlinear interactions between variables, identifies interactions between variables, and scales to larger data sets than those handled by traditional Gaussian process regression. We will introduce sparsity via spike and slab priors (Mitchell & Beauchamp, 1988;George & McCulloch, 1993) within a semiparametric regression framework that reduces the dimension of the parameter space, while yielding inference on interactions among individual exposures in the mixture. Section 2 introduces the motivating data set looking at exposure to metals in Bangladesh. Section 3 introduces our model formulation and discusses various estimation issues. Section 4 discusses some results about our prior distribution. Section 5 illustrates the proposed approach through simulated data and compares it with existing approaches. Section 6 estimates the health effects of environmental exposures on neurodevelopment in Bangladeshi children. We will conclude with some final remarks in Section 7.

Chemical mixtures and neurodevelopment in Bangladesh
To study the health impact of chemical mixtures on neurodevelopment, we examine a study of 375 children in Bangladesh whose mothers were potentially exposed to a variety of chemical pollutants during pregnancy. The outcome of interest is the z-score of the motor composite score (MCS), a summary measure of a child's psychomotor development derived from the Bayley Scales of Infant and Toddler Development, Third Edition (BSID-III) (Bayley, 2006). The exposure measures are log-transformed values of Arsenic (As), Manganese (Mn), and Lead (Pb) measured in umbilical cord blood, which is taken to represent prenatal exposure to these metals. As potential confounders of the neurodevelopment -exposure associations, we will also control for the following covariates: gender, age at the time of neurodevelopment assessment, mother's education, mother's IQ, an indicator for which clinic the child visited, and the HOME score, which is a proxy for socioeconomic status. Previous studies have examined these data (Bobb et al. , 2014;Valeri et al. , 2017) and found harmful associations between the exposures and the MCS score. Further, these studies have suggested a possible interaction between Arsenic and Manganese, though it was unclear whether this was simply random variation or a true signal. The assessment of interactions within environmental mixture data is an important problem  as the identification of exposures that interact with one another can help identify groups of exposures, or exposure profiles, that have particularly strong effects on health. The ability to identify them and assign measures of uncertainty around the existence of interactions is important, and will represent an important aspect of our proposed approach.

Preliminary analyses
Before examining any associations with the outcome, we first look at the Pearson correlation among the log exposure values. Manganese and lead are uncorrelated, manganese and arsenic have a correlation of 0.58, and arsenic and lead have a correlation of -0.37. Next, we evaluate marginal associations between the exposures and outcome, unadjusted for any additional covariates or co-exposures. Table 1 shows these results. We see that all three exposures have significant marginal correlations with the MCS score. We also evaluate whether a linear or nonlinear relationship between the MCS score and the exposures was favored by the data by looking at AIC values. The nonlinear AIC value is obtained by regressing the MCS score against a 3 degree of freedom spline basis representation of the exposures. We see that for Arsenic, the nonlinear model is favored, while the linear model is favored for the other two exposures. Arsenic concentrations are positively associated with MCS score, which could possibly be due to unmeasured confounding by, say, nutrients, due to the fact that a significant exposure pathway of Arsenic exposure in this population is diet.

Exposure
Correlation Next, to highlight the importance of the new methodology we propose to address the presence of interactions between environmental exposures, we analyze the data using current state of the art approaches for the analysis of chemical mixtures . Specifically, we use the Bayesian Kernel Machine Regression approach (BKMR, Bobb et al. (2014)), which is now widely used within the environmental epidemiological literature (Braun, 2017;Harley et al. , 2017). BKMR allows for nonlinear relationships, interactions among exposures, and adjusts for additional covariates measured in the population. We use the bkmr package in R, generating 25,000 Markov chain Monte Carlo (MCMC) scans and the default prior specification in the package.  These measures quantify the amount of uncertainty regarding the importance of each of the three expo-sures in predicting MCS score, but do not reflect whether or not the three exposures interact in their effects on the outcome. One could get posterior inclusion probabilities such as these if all three exposures have distinct, important main effects, or if there is a three-way interaction between them. To begin to address this question, one can visualize the estimated exposure response curves between these exposures. We visually inspect the relationship that Manganese and Arsenic have on the MCS score. In particular, Figure 1 shows the relationship between Manganese and the MCS score while fixing Arsenic at different quantile levels (0.1, 0.5, and 0.9). If there is an interaction between these two variables then these curves would differ across the different quantiles of Arsenic, whereas if there were no interaction, then the shape of the three curves would not differ. We see that there is some variability in the shape of the curves, particularly for the 0.5 quantile, which has an inverted U shape. This suggests the possibility of an interaction between Manganese and Arsenic; however, this can not be concluded definitively given the combination of the posterior inclusion probabilities and the exposure-response curves. This shows the importance of new statistical methodology to identify these interactions

Model formulation
Throughout, we assume that we observe where n is the sample size of the observed data, Y i is the outcome, X i is a p-dimensional vector of exposure variables, and C i is an mdimensional vector of covariates for subject i. Building on the construction of Wood (2006), we assume that we have low rank bases that represent smooth functions for each of the exposures. For instance, the main effect of X 1 and an interaction effect between X 1 and X 2 could be modeled as follows: f (X 1 ) = X 1 β 1 f (X 1 , X 2 ) = X 1 β 1 + X 2 β 2 + X 12 β 12 , respectively, where X 1 = [g 1 (X 1 ), . . . , g d (X 1 )] and X 2 = [q 1 (X 2 ), . . . , q d (X 2 )] represent d−dimensional basis function expansions of X 1 and X 2 , respectively. X 12 = [g 1 (X 1 )q 1 (X 2 ), g 1 (X 1 )q 2 (X 2 ), . . . , g d (X 1 )q d (X 2 )] represents a d 2 −dimensional basis expansion of the interaction between X 1 and X 2 . In this work, we specify g() and q() to be natural spline basis functions of X 1 and X 2 , respectively. For simplicity we let the degrees of freedom for each exposure all be set to d, though this can easily be extended to allow the degrees of freedom to vary across exposures. Importantly, the construction above enforces that main effect terms are a subset of the terms included in the interaction terms. As we include higher-order interactions, we impose the constraint that all lower level interactions between the exposures in that interaction term are included in the basis expansion. To ease the presentation of our model, we use subscripts to denote the order of the interaction to which a coefficient applies. For instance, β 12 refers to parameters associated with two-way interaction terms between X 1 and X 2 (excluding main effect terms), while β 1 refers to the main effect parameters for X 1 . This process could be trivially extended to any order interaction up to a p−way interaction, which is the maximum possible interaction in our model. We will use the following model formulation: where the dots in the final line indicate that the summation extends all the way through p−way interactions. Our main goal is the estimation of f (X i ) and the identification of important exposures within X i . We are representing f (X i ) as the sum of k different functions f (h) (X i ), where k is a number chosen to be sufficiently large so that all exposure effects can be captured by the model. This model is both over parameterized and not identifiable without any additional constraints. This over parameterization occurs because the model allows for all possible interactions of any order, which includes far too many parameters to reasonably estimate. It is not identifiable because each of the k functions has the same functional form and the regression coefficients for a particular function f (h) (X i ) are only identifiable up to a constant. To reduce the dimension of the parameter space and to identify parameters of interest, we will adopt the following prior specification for the regression coefficients: Here ζ = {ζ jh } is a matrix of binary indicators of whether the exposure j is included in the h th function in the model. One of our main goals is to identify important exposure effects and interactions between exposures, and these can be seen immediately from {ζ jh }, which tells us both which exposures and which interactions are included in a particular function. We use spike and slab priors to eliminate exposures from the model, and in this paper we consider the spike to be a point mass at 0, while the slab, ψ 1 (), is a multivariate normal distribution centered at 0 with covariance Σ β . We let Σ β be a diagonal matrix with σ 2 σ 2 β on the diagonals. The prior distribution for ζ jh is a Bernoulli distribution with parameter τ h . This prior is accompanied by an indicator function, 1(A h ⊂ A m ∀m or A h = {}), which reflects the fact that each function either contains information from a unique set of exposures or is removed from the model completely. If a function only contained exposures that were all included in another function, then it would be providing redundant information and it is not needed, so we preclude this from happening by introducing this indicator function. More details on this indicator function and its implications for sampling can be found in the following section. Full details of posterior sampling from this model can be found in the Appendix.

Identifiability of variable inclusion parameters
Two scientific goals in the study of environmental mixtures are to identify which exposures are associated with the outcome and to identify interactions among exposures. To this end, our prior encodes the condition that, for any two sets of variables A h and A m included in functions h and m, either A h ⊂ A m and A m ⊂ A h or one of the two is an empty set. An illustrative example of why this condition is necessary is a scenario in which all elements of ζ are zero with the exception of ζ 11 and ζ 12 . In this scenario, the first exposure has a nonzero effect in both the first and second functions of the model, i.e. β (2) 1 = 0, inference regarding ζ would differ across these two models. In the first scenario it would appear that the first exposure is associated with the outcome, even if in truth it is not. Our conditions avoid this scenario by ensuring that there are no functions that simply contain subsets of the exposures found in another function, and therefore ζ will provide a meaningful and interpretable measure of the importance of each exposure and potential interaction.
Another issue affecting our model selection indicators is that low-order interactions are subsets of the models defined by higher-order interactions. For example, if the true model is such that there exists a main effect of both X 1 and X 2 , but no interaction between them, we could mistakenly estimate that there is an interaction between exposure 1 and 2 simply due to the main effects. (Wood, 2006;Reich et al. , 2009) addressed this issue in similar contexts by utilizing constraints that force the two-way interaction between X 1 and X 2 to only include information that is orthogonal to the main effects of the individual exposures. In our context, this is not so much an issue of identifiability, but rather an issue with MCMC chains getting stuck in local modes. When exploring the model space via MCMC we could accept a move that includes the two-way interaction, which is a local mode, and never be able to leave this model even if the simpler model is the true model.
To avoid this scenario we simply need to search the model space in a particular manner. We will illustrate this in the case of a two-way interaction as follows: Imagine that at a given iteration of an MCMC, ζ = ζ p1 and we look to update ζ 21 , which indicates whether exposure 2 enters into function 1. If we do this naively, we simply could compare models ζ p1 and ζ p2 , which would result in a choice of ζ p2 even if there were only main effects of X 1 and X 2 and no interaction between them. If, however, we update ζ 22 at the same time, we also evaluate model ζ p3 . If the true relationship is a main effect of the two exposures, then model ζ p3 will be favored by the data and we will estimate ζ accurately. If the truth is a two-way interaction between X 1 and X 2 then ζ p2 will be favored by the data. This setting generalizes to higher-order interactions and illustrates how we can avoid misleading conclusions about interactions. Whenever we propose a move to a higher-order interaction, we must also evaluate all lower-order models that could lead to incorrectly choosing the higher-order model. While proposing moves in this manner keeps the MCMC from getting stuck at incorrect models, it needs to be carefully constructed to maintain a valid sampler. In the Appendix, we prove that this strategy can be tailored in such a way that the sampler is reversible and therefore maintains the correct target distribution.

Selection of spline degrees of freedom
An important tuning parameter in the proposed model is the number of degrees of freedom, d, used to flexibly model the effects of the p exposures. Throughout we will make the simplifying assumption that the degrees of freedom is the same for all exposures, though this can be relaxed to allow for differing amounts of flexibility across the marginal functions. We propose the use of the WAIC model selection criterion (Watanabe, 2010;Gelman et al. , 2014) to select d. If we let s = 1, . . . , S index samples from the posterior distribution of the model parameters and let θ represent all parameters in the model, the WAIC is defined as The WAIC approximates leave one out cross validation, and therefore provides an assessment of model fit based on out of sample prediction performance. We first draw samples from the posterior distributions of the model parameters under a range of values for d, and then select the model that minimizes WAIC. While this is not a fully Bayesian approach to choosing d, we will see in Section 5 that this approach still leads to credible intervals with good frequentist properties.

Prior choice
There are a number of hyperparameters that need to be chosen to complete the model specification. First, one must choose an appropriate value for Σ β , which controls the amount of variability and shrinkage of β for those coefficients originating from the slab component of the prior. We will make the simplifying assumption that Σ β is a diagonal matrix with σ 2 σ 2 β on the diagonals. It is well known that variable selection is sensitive to the choice of this parameter (Reich et al. , 2009), so it is crucially important to select an appropriate value. We propose the use of an empirical Bayes strategy where we estimate σ 2 β and then obtain posterior samples conditional on this value. A description of the algorithm for finding the empirical Bayes estimate is available in the Appendix. One could alternatively put a hyper prior on σ 2 β for a fully Bayesian analysis, though we have found that the empirical Bayes strategy works better in practice, so we restrict attention to that case.
We also must select a value of τ h , the prior probability of an exposure being included in function h. We will assume that τ h ∼ Beta(M, γ), and set M and γ to fixed values. One can choose values to induce a desired level of sparsity in the model, however in practice, the appropriate degree of sparsity is typically not known and therefore, in order to accommodate higher-dimensional data, we choose values that scale with the number of exposures. We will set γ = p, which allows the prior amount of sparsity to increase as the number of exposures increases. This strategy closely resembles other existing approaches in high-dimensional Bayesian models, such as those proposed by Zhou et al. (2014) and Ročková & George (2016). We will discuss the implications of scaling the prior with p in Section 4.
The only remaining prior distributions to specify are for β C and σ 2 . For both parameters we specify diffuse, conjugate priors. In the case of β C this means independent normal prior distributions with very large variances. For the residual variance σ 2 we specify an inverse gamma distribution with both parameters set to 0.001.

Lower bound on slab variance
One issue with the aforementioned empirical Bayes estimation of σ 2 β is that it will be estimated to be very small if there is little to no association between the exposures and the outcome. This result can lead to problems in the estimation of ζ. If the prior variance for the slab is very small, then it is almost indistinguishable from the spike, since in this case it approximates a point mass at zero. In this situation, updating ζ jh becomes essentially a coinflip for all j and h, which will lead to high posterior inclusion probabilities even when there is no true association. To avoid this, we estimate a lower bound, above which we believe that any high posterior inclusion probabilities reflect true signals in the data.
To estimate this lower bound, we first permute the rows of Y , thereby breaking any associations between the exposures and the outcome. We run our MCMC algorithm for a small number of iterations on the permuted data and keep track of the probability that ζ jh = 1 when τ h = 0.5 for all h. We do this for a decreasing grid of potential values for σ 2 β and we take as our lower bound for σ 2 β the smallest value of σ 2 β such that a pre-chosen threshold is met. Examples of such criteria include the smallest value of σ 2 β such that the probability of including a main effect for any of the exposures is less than 0.25 or the smallest value such that the probability of any two-way interaction is less than 0.05. Alternatively, because one can calculate the average false discovery rate for each value of σ 2 β , one can use this strategy to implement Bayesian counterparts to traditional multiple testing adjustments that control frequentist error rates, such as the false discovery rate. It is important to note that this false discovery rate control is conditional on a fixed value of τ h , which we choose to be 0.5. In practice, we let τ h be random and estimated from the data. Nonetheless, this procedure should provide a reasonable lower bound for σ 2 β . If strict control of false discovery rate is required then one could specify a truncated prior for τ h such that it is never greater than the value of τ h used for the calculation of the lower bound (i.e., 0.5). It is important to note that this procedure is very fast computationally, generally taking only a small fraction of the computation time as the full MCMC for one data set.
To illustrate how this approach works for a given analysis, we simulated data with no associations between the exposure and outcome. In this case, the empirical Bayes estimate of σ 2 β will be very small as it is roughly estimated to be the average squared value of the nonzero coefficients in β. Figure 2 shows the posterior inclusion probabilities averaged over all exposures in this example as a function of σ 2 β . The gray region represents the area where false positives start to arise because σ 2 β is too small, and we see that the empirical Bayes estimate is contained in this region. If we proceed with the empirical Bayes estimate, then we would conclude that all of the exposures had posterior inclusion probabilities near 0.5, incorrectly indicating they are important for predicting the outcome. Our permutation strategy to identifying a lower bound is able to find a value of σ 2 β that avoids high posterior inclusion probabilities for exposures with null associations with the outcome. In practice, we estimate both the lower bound and empirical Bayes estimate and proceed with the maximum of these two values, which in this simulated case is the lower bound.

Properties of Prior Distributions
A primary goal in many analyses of the health effects of an environmental mixture is to identify interactions among exposures on risk. Therefore, it is important to understand the extent to which the proposed prior specification induces shrinkage on the interaction terms. Both τ h and Σ β dictate the amount of shrinkage. We utilize empirical Bayes to select Σ β . Therefore, here we focus on the role of τ h , which directly controls the probability that a j th order interaction is nonzero in the model. The prior for this parameter can be controlled through selection of M and γ. We examine the probability that a j th order interaction is zero by deriving the probability that the first j covariates are included in an interaction term, as this probability will be the same for any other interaction of the same order.
Lemma 3.1: Assuming Model 1, the prior probability that there does not exist a function, k in the model such that ζ 1k = ζ 2k = · · · = ζ jk = 1 and ζ j+1k = · · · = ζ pk = 0 is The Appendix contains a proof of this result, which is useful because it expresses the probability of an interaction of any order being zero as a function of M and γ. This result can guide the user in specifying values of M and γ that best represents prior beliefs. It has been noted in the literature that a desirable feature of the prior is to place more shrinkage on higher-order interactions (Gelman et al. , 2008;Zhou et al. , 2014). Therefore we use (2) to derive conditions that accomplish this type of shrinkage. Specifically, Theorem 3.2 provides a default choice of γ that specifies the amount of shrinkage on an interaction term to be an increasing function of the order of that term. The Appendix contains a proof of this result, which shows that higher-order interactions are more aggressively forced out of the model under either the standard Beta(1, p) prior, which is often used in highdimensional Bayesian models (Zhou et al. , 2014;Ročková & George, 2016), or if γ ≥ p + M − 1, which shrinks higher-order interactions even more aggressively than this beta prior. In practice, a common prior choice for τ h is Beta(M, p), with M being some constant that does not change with the number of exposures. Nonetheless, this result shows that one can achieve stronger penalization of higher-order interaction terms as long as we scale the prior with p.

Simulation study
Here we present results from a simulation study designed to assess the operating characteristics of the proposed approach in a variety of settings. In all cases we take the maximum of the empirical Bayes and lower bound estimates of σ 2 β . For each simulated data set, we fit the model for values of d ∈ {1, 2, 3, 4, 5, 6}, and run a Gibbs sampler for 10,000 iterations with two chains, keeping every 8th sample and discarding the first 2,000 as burn-in. We present results for each value of d, the value that is selected based on the WAIC, the spike and slab GAM (ssGAM) approach of Scheipl et al. (2012), and for the Bayesian kernel machine regression (BKMR) approach of Bobb et al. (2014). We assess the performance of the various models by calculating the mean squared error (MSE) and 95% credible interval coverages of the true relationship between the exposures and the outcome, f (X 1 , ..., X p ), at all observed values of the exposure. We also assess the ability of the proposed model to identify nonzero exposure effects and interaction terms.

Polynomial interactions
First, we generate data with a sample size of n = 200, and generate m = 10 covariates, C, from independent standard normal distributions. We generate p = 10 exposures, X, from a multivariate normal distribution with mean µ X . We will let the mean of X be a linear function of the covariates, i.e. µ X = Cα. We let all pairwise correlations between the exposures be 0.6, which is the average pairwise correlation among the exposures in the NHANES data analyzed in Appendix A. The marginal variances are set to 1 for all exposures. We set the residual variance to σ 2 = 1. For the outcome, we use a polynomial model that contains both a linear and a nonlinear interaction term as follows: For exact values of α and β C see Appendix B. Figure 3 presents the results of the simulation study. The top left panel of Figure 3 shows the matrix of posterior probabilities that a pair of exposures interact with each other. Specifically, cell (i, j) of this heatmap shows the posterior probability that exposures i and j interact with each other in the model averaged over all simulations. Results show that the model correctly identifies the pairs (X 2 , X 3 ) and (X 4 , X 5 ) as interacting with each other. The top right panel shows the average 95% credible interval coverage of f (X i ) for each approach considered. Our model based on WAIC achieves nearly the nominal coverage indicating that we are accounting for all sources of uncertainty. The bottom left panel highlights the ability of our procedure to estimate the marginal effect of X 4 , which has a quadratic effect on the outcome, while holding constant the values of other exposures. The grey lines, which represent individual simulation estimates, reflect the true relationship, denoted by the red line. Finally, the bottom right panel of Figure 3 presents the mean squared error of the global effect estimates of f (X i ). Results show that well chosen values of d, including the automated choice using WAIC, lead to the best performance among the methods considered.

High-dimensional data with three-way interaction
Now we present a scenario with n = 10, 000 and p = 100 to highlight that our approach can handle large sample sizes and large exposure dimensions. We will generate the exposures from independent uniform random variables. We generate one additional covariate, C, from a standard normal distribution, and we set the residual variance to σ 2 = 1. The outcome model is closely related to one that was explored in Qamar & Tokdar (2014) and corresponds to a model that includes a three-way, nonlinear interaction as follows: Y = 2.5 sin(πX 1 X 2 ) + 1.5 cos(π(X 3 X 4 + X 5 )) + 2(X 6 − 0.5) + 2.5X 7 + It is clear that the relationship between the exposures is nonlinear and non-additive. We do not show results for BKMR due to the sample size limitation from using Gaussian processes. Further, there are a large number of exposures considered here, as there are 4,950 possible two-way interactions and 161,700 possible three-way interactions. For this reason, we do not include ssGAM in this simulation as it requires specifying the form of the linear predictor, and there are too many terms to consider. Our approach considers interactions of any order, though it does not require pre-specification of which terms to consider, and aims to find the correct model among all possible interactions. Figure 4 shows results from this scenario. The top left panel shows the posterior probability of a three-way interaction between X 3 , X 4 and X 5 . This probability is 1 for all models except for d = 6, which does not capture this interaction. The top right panel shows the 95% credible interval coverages of f (X i ). We see that the coverage is low for d = 1 and d = 2, which are not sufficiently flexible to capture the true function, while the coverage is low for d = 6 since it does not capture the three-way interaction. Importantly, the model chosen by WAIC achieves the nominal coverage level. The bottom left panel of Figure 4 shows the estimated marginal effect of X 4 on the outcome, fixing the other exposures. Results suggest that in almost all cases the method accurately estimates the true effect, which is close to linear in the range of the data we are looking at. The bottom right panel of Figure 4 shows that the MSE of the model chosen by WAIC is competitive, or better, than any individual choice of d, suggesting that WAIC is an effective method for choosing the degrees of freedom for the spline terms.

Analysis of Bangladesh data
Here we analyze the data on metal mixtures and neurodevelopment introduced in Section 2. We applied our proposed approach to identify the complex relationship between the three metal exposures and neurodevelopment. We ran the Markov chain for 40,000 iterations with M = 3, γ = p, a 0 = 0.001, b 0 = 0.001, and a non-informative, normal prior on β c . We ran the model for values of d ∈ {1, 2, 3, 4, 5}. Convergence of the MCMC was confirmed using trace plots and PSR values, and these results can be seen in the Appendix.

The importance of the degrees of freedom
Here we evaluate the impact that the degrees of freedom d had on the analysis. Table 3 shows the WAIC values for each of the values of d considered, the corresponding main effect posterior inclusion probabilities for each exposure, as well as the posterior inclusion probability for the two-way interaction between Manganese and Arsenic. The WAIC values strongly indicate that linear models are not appropriate for this data, as higher degrees of freedom lead to lower WAIC values and therefore better model fits. One of the most telling takeaways from Table 3 is the correspondence between WAIC values and the posterior inclusion probabilities for the two-way interaction between Manganese and Arsenic. Models with higher degrees of freedom lead to lower WAIC values and higher posterior inclusion probabilities for the two-way interaction. It appears that the lower values of d are not adequate for capturing the interaction between Manganese and Arsenic, and this leads to poorer fits to the data. The marginal inclusion probabilities for each exposure are more stable across d, however, these go up slightly from d = 1 to d = 5 as well. Overall, the WAIC indicates that d = 5 provides the best fit to the data and we will therefore present results from this model from this point onwards.  Table 3: WAIC and posterior inclusion probabilities for differing degrees of freedom.

Posterior inclusion probabilities
One feature of the proposed approach is the ability to examine the posterior distribution of ζ in order to identify whether individual exposures are important for outcome prediction, whether there is evidence of an interaction between two metals or any higher-order interactions. Figure 5 shows the posterior probability of two-way interactions between all metals in the model, and the marginal posterior inclusion probability for each metal exposure. Figure 5 provides strong evidence that each exposure is associated with a child's MCS score. This confirms the preliminary results of Section 2 that showed associations between each of the three exposures and the MCS score. A key distinction of the results here and those seen in Section 2 is that it also provides strong evidence of an interaction between Mn and As, as well as moderate evidence of an interaction between Pb and As. Not shown in the figure is the posterior probability of a three-way interaction among Pb, Mn, and As, which was estimated to be 0.0005.

Visualization of identified interactions
Now that we have identified interactions in the model, we plot posterior predictive distributions to assess the form of these associations. Figure 6 plots the relationship between Mn and MCS for different values of As, while fixing Pb at its median. At the median level of As, we estimate an inverted U-shaped effect of Mn on MCS scores. This matches results seen previously in a kernel machine regression analysis of these data Bobb et al. (2014), which provided graphical evidence of an interaction between Mn and As on the MCS score but didn't provide a measure of the strength of this interaction. We also examine surface plots that show the exposure-response relationship between As, Mn, and MCS across the range of As and Mn values, not just specific quantiles. Figure 7 shows both the posterior mean and pointwise standard deviations of this relationship. We see that the median levels of Mn and As are the most beneficial towards children in terms of their MCS score, which corresponds to the hump in the middle panel of Figure 6. This is a somewhat unexpected result given that As is not beneficial even at low levels, and could be due to residual confounding (Valeri et al. , 2017). The posterior standard deviations of this estimate indicate that the uncertainty around this surface is the smallest in the center, where most of the data is observed, and greater on the fringes of the observed data, where we are effectively extrapolating effects.

Discussion
We have proposed a flexible method for analyzing complex, nonlinear associations between multiple exposures and a continuous outcome. Our method improves upon existing methodology for environmental mixtures in several ways. By not relying on Gaussian processes to introduce flexibility in the exposure-response relationship, our approach can handle larger data sets on the order of tens of thousands of subjects. The spike and slab formulation of our model also allows for formal inference on the presence of interactions among exposures, an important question in environmental epidemiology. Previous studies have used graphical tools and inference on summary statistics characterizing the effects of one pollutant given the levels of others to identify interactions between pollutants. While useful, these approaches do not provide uncertainty assessment of the interaction as a whole. Moreover, it can be difficult to visualize higher-order interactions.
Once we have samples of the ζ matrix, we can easily perform inference on more complex features of the model such as higher-order interactions and whether or not exposures are more likely to enter the model together or separately. We have illustrated the strengths of this approach through simulated data, and a study of the health effects of metal mixtures on child neurodevelopment in Bangladesh. A natural question that arises in the application of our method is how does it perform on higherdimensional mixtures in real-world settings. We evaluated a situation with both a large sample size and number of exposures in Section 5.2 of our simulation study, and found that it performs well at identifying higher-order interactions. To further illustrate this point, we have analyzed a secondary data set from the National Health and Nutrition Examination Survey (NHANES), which has a larger sample size and more exposures than the Bangladesh data. The results of this study can be found in Appendix A.
Throughout the manuscript we have assumed that the relationship between the exposures and outcome is constant across levels of the covariate space C. In some cases, there is interest in understanding if the exposure effect varies across levels of particular covariates such as age or gender, which is frequently referred to as treatment effect heterogeneity. This could easily be explored using our model if the covariates are continuous. The covariates of interest can be placed with the exposures inside of the main f () function, and our approach will be able to identify if any of these covariates modify the exposure effect. Identifying treatment effect heterogeneity in the presence of categorical covariates is not immediately applicable within our modeling framework, but it provides a topic of future research.
In the modeling framework we focused on a high-dimensional X and assumed that the number of covariates in C was relatively small, meaning we do not need to impose any shrinkage on β C . If the number of covariates is large as well, then shrinkage or spike and slab priors would need to be adopted. If interest is in estimation of f (X), this could lead to what is called regularization induced confounding (RIC, Hahn et al. (2017). RIC occurs because shrinkage leads to biased estimates of the coefficients in β C , and this bias leads to confounded estimates of the exposure effects of interest. Future work would look to extend the ideas here to the more difficult setting of a high-dimensional C when RIC could affect inference of the effects of interest. This would require reducing the shrinkage of the elements in C that are most associated with the elements of X in a manner that removes RIC bias.
While our formulation enables the analysis of much larger data sets, there are some disadvantages relative to fully nonparametric techniques such as Gaussian processes. Our model requires us to specify the functional form of the relationships between the predictors and the outcome a priori. In the current implementation, we employed natural splines to flexibly model these relationships. While other formulations are possible, any formulation would require a subjective choice. Further, one must select the number of degrees of freedom of the splines, and we have restricted this number to be the same for all exposures across all features. This can be a limitation if some of the functions to be estimated are quite smooth while others are more complicated. Our approach is also not fully Bayesian, as we use the WAIC to select the degrees of freedom. However, our simulations show the proposed estimation procedure yields uncertainty estimates that are accurate in all of the scenarios we considered. Future extensions could consider a data-adaptive approach to the selection of the degrees of freedom for each function, or the utilization of other flexible structures that are both computationally fast and do not require pre-specification of the functional forms of the exposure-response functions.

Software
An R package implementing the proposed methodology is available at github.com/jantonelli111/NLinteraction Acknowledgements This publication was made possible by USEPA grant RD-835872-01. Its contents are solely the responsibility of the grantee and do not necessarily represent the official views of the USEPA. Further, USEPA does not endorse the purchase of any commercial products or services mentioned in the publication. This work was also supported by grants from the National Institutes of Health (ES028811, ES000002, ES024332, ES007142, ES016454).

A Effects of pollutants on telomere length in NHANES
We will utilize data from the 2001-2002 cycle of the National Health and Nutrition Examination Survey (NHANES). The data is a representative sample of the United States and consists of 1003 adults. The outcome of interest is leukocyte telomere length (LTL). Telomeres are segments of DNA at the ends of chromosomes that are used to help protect chromosomes, and their lengths generally decrease with increasing age. In particular, LTL has been used as a proxy for overall telomere length, and has been shown to be associated with a number of adverse health outcomes (Mitro et al. , 2015). In the NHANES data, we have measurements from 18 persistent organic pollutants, which consist of 11 polychlorinated biphenyls (PCBs), 3 dioxins, and 4 Furans. It is believed that these environmental exposures could impact telomere length, though it is unclear how telomere length is associated with this particular mixture of exposures. Also in the data are additional covariates for which we need to adjust including age, sex, BMI, education status, race, lymphocyte count, monocyte count, cotinine level, basophil count, eosinophil count, and neutrophil count. We will analyze these data to assess whether any of the aforementioned environmental exposures are associated with LTL, and to search for any interactions between the environmental exposures. This quickly becomes a high-dimensional statistical issue even when looking at two-way interactions only as there are 153 two-way interactions alone. Therefore, this data requires a statistical approach to estimate high-dimensional models that allow for interactions and nonlinearities in the associations between the environmental mixture and LTL.
We ran the Markov chain for 40,000 iterations with M = 3, γ = p, a 0 = 0.001, b 0 = 0.001, and a noninformative, normal prior on β c . We applied our approach with degrees of freedom, d ∈ {1, 2, 3, 4, 5, 6, 7} and we will present the results for d = 3 as it was the model that minimized the WAIC. We examined the convergence of our MCMC using both trace plots and the potential scale reduction factor (PSR, Gelman et al. (2014)) in Appendix H, and find that our model has converged. Figure 8 shows the marginal inclusion probabilities for each of the 18 exposures. We see that all of the exposures have a posterior inclusion probability at or near zero with the exception of Furan1, which was included in the model 99% of the time. We do not detect any higher-order interactions, as the posterior inclusion probabilities were exactly zero for all two-way and higher-order interactions. This provides strong evidence of a relationship between Furan1 and telomere length.  Figure 9 shows the effect of Furan1 on telomere levels while fixing other exposures at their medians. It is important to note that fixing other exposures at their medians does not change the shape of the curve in this instance, because there are no interactions between Furan1 and other exposures. We estimate a positive association between Furan1 exposure and telomere length. The curve is somewhat nonlinear as seen in both figure 9 and the fact that the WAIC chose three degrees of freedom for the basis functions.

B Simulation details and additional simulation results
Here we present more details regarding the simulation study in Section 5.1 of the manuscript, an additional simulation study with increased correlation between the exposures and covariates, as well as an additional simulation study that does not include any interactions, but has nonlinear main effect terms.
B.1 Details of data generation for simulation 5.1 In Section 5.1 of the manuscript, we generated exposures from a multivariate normal distribution with mean µ X = Cα, and then generated the outcome from the following model: We can also look at the correlation matrix between X and C to understand the impact of these parameter choices. We can see in Figure 10 that there is a very high amount of correlation between the 10 elements of X denoted by the lower left block in the correlation matrix. Looking at the correlation between X and C we can see that each element of X has certain elements of C that are correlated with it.

B.2 Increasing correlation between X and C
Now we will increase the correlation between X and C to see how inference is impacted by this correlation. We will now let the mean of each exposure be a function of three covariates (instead of 2 as in simulation 5.1) Figure 10: Correlation matrix among all exposures and covariates for simulation 5.1 and we will increase the magnitude of this correlation. The correlation matrix for this new setup can be found in Figure 11, and we can see there is substantially more correlation between X and C. We will compare results from two situations: the results from simulation 5.1 of the manuscript, and this new simulation with increased correlation between the exposures and covariates. The data generating models for the outcome are the same between both these scenarios with the main difference being the correlation matrices for X and C.
The mean squared error from this additional scenario can be seen in Figure 12. We only highlight our approach based on WAIC to choose the degrees of freedom and the BKMR approach as these two were the best performing approaches in the manuscript. We see that the performance of both approaches degrades when the correlation is increased between the covariates and exposure. This increase in MSE is fairly similar between the two approaches, however, and the approach based on WAIC performs the best between the two approaches. Not shown in the figure is that the approach based on WAIC also achieves interval coverages of 91% indicating that performance of the approach is still good even under high levels of correlation between the covariates and exposures.

B.3 Additional simulation study
Now we present a simulation study that is similar to the one presented in Section 5.1 of the manuscript, except the effect of the exposures has been changed. The covariates and exposures are generated from the same multivariate normal distribution as before. Now the outcome is generated from the following model: The results can be seen in Figure 13. The top left panel shows that our model is correctly identifying the four important exposures while excluding the others. Not shown in the figure is that all two-way interactions had posterior inclusion probabilities less than 0.01, indicating that our model correctly removes all higherorder interactions. The top right panel shows that none of the approaches are able to achieve the nominal coverage in this instance, though again ssGAM has confidence intervals that are far too wide leading to a coverage probability of one. Of the other estimators though, our approach with degrees of freedom chosen by WAIC performs best in terms of coverage probability with 86%. The bottom left panel shows that we are correctly capturing the nonlinear relationship between X 4 and Y . Lastly, the bottom right panel shows that our approach obtains the lowest MSE of all approaches considered.  1. Sample σ 2 from the following full conditional: Sample τ h for h = 1, . . . , k from the following full conditional: To update ζ let us first introduce some additional notation. Let ζ B be the subset of ζ we are updating and β B be the corresponding regression coefficients. Further, let Λ be the set of all parameters in our model excluding ζ B and β B . To update ζ B we need to calculate p(ζ B |D, Λ), which we can do via the following:  Figure 12: Mean squared error in estimating f (X i ) for the two approaches considered under both simulation scenarios. The first scenario is that of simulation 5.1 of the manuscript, and the second simulation is with the increased correlation between the covariates and exposures.
where Φ() represents the multivariate normal density function. M and V represent the conditional posterior mean and variance for the elements of β corresponding to the elements of ζ B that are equal to one and are defined as where X * is a design matrix with the appropriate main effects and interactions induced by the model Intuitively, the numerator in (3) is the prior probability that all elements of β corresponding to features in ζ B are zero and ζ B = z B . The denominator in (3) represents the conditional posterior probability that all elements of β corresponding to features in ζ B are zero.

Update β
(h) S from a normal distribution with mean M and variance V as defined above.
5. Update β c from the following distribution:

D Derivation of empirical Bayes variance
As seen in Casella (2001), the Monte Carlo EM algorithm treats the parameters as missing data and iterates between finding the expectation of the "complete data" log likelihood, where expectations are calculated from draws of a gibbs sampler, and maximizing this expression as a function of the hyperparameter values. After removing the terms not involving σ β , we write the "complete data" log likelihood as X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 1 0 95% interval coverage of f(X) where µ β S is the prior mean for coefficient β S , and |β| 0 is the number of nonzero regression parameters. The estimate of σ 2 β at iteration m uses the posterior samples from iteration m − 1, which were sampled with σ 2 β = σ 2 β (m−1) . Adopting the same notation that is typically used with the EM algorithm, we compute the expectation of this expression, where expectations are taken as the averages over the previous iterations posterior samples: We then take the derivative of this expression with respect to σ 2 β to find that the maximum occurs at Intuitively, to obtain the empirical Bayes estimate, one must run an MCMC chain and update σ 2 β every T MCMC iterations. This process is continued until the estimate of σ 2 β converges, and then the MCMC can be run conditional on this estimated value. E Proofs of lemma 3.1 and theorem 3.2 Here we investigate the prior probability that a particular interaction term is equal to zero. Our goal is to penalize interactions of higher order more strongly, and the probability that we force them to zero is one component of this strategy. Since the prior probabilities of inclusion into function k for each covariate are independent, we will investigate, without loss of generality, the probability that ζ 1k = ζ 2k = · · · = ζ jk = 1 and ζ j+1k = · · · = ζ pk = 0. The probability would remain constant if we changed which of the p exposures we chose to be included or excluded from the interaction term, so we will look at the case where the first j terms are included for simplicity.

E.1 Probability increases as j increases if γ ≥ p
It is of interest how the variable inclusion proabilities vary as a function of j, the order of each interaction. Ideally we would want to show that the probability increases as a function of j. We don't expect this to happen in general, but we can find conditions at which it does hold. To simplify expressions with the gamma function we restrict γ to be an integer, though we expect the results to extend to all real values of γ.
We want this ratio to be greater than 1 for all j ∈ {1, . . . , p − 1} for the function to be increasing. Clearly this ratio is getting smaller as j gets larger so we can look at the largest value of j, which is p − 1.
So clearly if we set γ ≥ p + M − 1, we have the desired result.

F Reversibility of updates for variable inclusion matrix
Here we will show a condition under which our proposed approach to sampling will maintain the desired target distribution. Let π() be the conditional density of ζ given the other model parameters and the data. Next let K(j|i) be the transition probability of going to ζ = j from ζ = i. Our goal is to show that reversibility holds, i.e. that π(i)K(j|i) = π(j)K(i|j) Our strategy for updating ζ is to randomly draw one of the pk elements in the matrix. if setting this element equal to 1 leads to an interaction of any order, then we must also evaluate whether the data could favor this move simply due to there being a lower order interaction. If we are evaluating a two-way interaction, then we must also evaluate whether there are simply main effects of each of the two exposures, otherwise we could incorrectly accept a move to an interaction model and bias our MCMC towards higher-order models. If including this element leads to a three way interaction, then we must evaluate whether in truth there is only the original two-way interaction, and then some additional terms (whether they be main effects or two-way interactions) that include the element we have chosen to update. This strategy amounts to creating a list of models to evaluate, and then the probability of moving to any one of them is proportional to their conditional posterior probability, π(). With this in mind, let us first look at K(j|i): where M ij is the set of all possible updates among the pk possibilities that contain a possible move to state j from state i. Further, ζ (m,i) denotes the set of all ζ matrices that are considered if we update m and we are currently at state i. Now let's put this back into the equation for reversibility: π(i)K(j|i) = π(j)K(i|j) ↔π(i) π(j) pk m∈Mij 1 z∈ζ (m,i) π(z) = π(j) π(i) pk m∈Mji 1 z∈ζ (m,j) π(z) . This equality must hold if two conditions hold. First M ij = M ji , i.e. the updates that can take you from state i to state j are the same updates that can take you in the other direction. Second, we need ζ (m,i) = ζ (m,j) , i.e. the space of models that are being considered if we update element m and are currently at state i is the same as those if we are at state j and we update element m. Thankfully we are in full control of the models that we choose to evaluate if we choose to update element m, so we can ensure this condition holds. Our first priority is to ensure that we consider all models that contain lower-order interactions as described in the manuscript so that we don't incorrectly move to a higher-order interaction when the truth is a lower-order model. Next, we will do a check to ensure that the required conditions hold, and if not, we will add any other models to evaluate that will make the conditions satisfied. If the model is complex and includes a large number of interactions, then this search can be computationally intensive, so we detail a strategy below using Metropolis Hastings that does not suffer from the same computational challenges, yet also satisfies reversibility.

F.1 Metropolis Hastings updates for variable inclusion matrix
Here we describe a strategy for proposing models and moving around the model space using Metropolis Hastings. As above, we will again choose one of the pk elements of the matrix randomly. If setting this element equal to 1 leads to an interaction of any order, then we will create a list of possible models that could lead to us choosing this model that includes all lower-order interactions. We then choose one of these models randomly and use it as a proposed new variable inclusion matrix. For a Metropolis Hastings step, we simply need the density π() and the proposal probabilities, Q(j|i) of going from state i to j. The proposal probability is straightforward to calculate here as it is simply one divided by the number of models in the list. These proposal probabilities can easily be calculated in the reverse direction as well to find the probability of going from the proposed model to the current model, and then the standard acceptance probability min 1, π(j)Q(i|j) π(i)Q(j|i) can be used. By being a Metropolis Hastings step this is guaranteed to satisfy reversibility, and it should not get stuck in modes with unnecessary higher-order interactions because it considers all models that could lead to the data favoring a higher-order move.

G MCMC convergence in the Bangladesh data
To assess convergence in the Bangladesh data set we kept track of the potential scale reduction factor (PSR, Gelman et al. (2014)) for all identifiable parameters in our model. Therefore we examined f (X i ) for i = 1, . . . , n, all elements of β C , and σ 2 . We found that the PSR value was below 1.005 for all parameters in our model, which leads us to believe that our MCMC has converged. A plot of the PSR values for f (X i ) can be found in Figure 14. Subsequently, we show the trace plots for all parameters in β C and σ 2 in Figure  15 and find that our chain appears to be converging.

H MCMC convergence in the NHANES data
To assess convergence in the NHANES data set we kept track of the PSR for all identifiable parameters in our model. Therefore we examined f (X i ) for i = 1, . . . , n, all elements of β C , and σ 2 . We found that the PSR value was below 1.005 for all parameters in our model, which leads us to believe that our MCMC has converged. A plot of the PSR values for f (X i ) can be found in Figure 16. Subsequently, we show the trace plots for a subset of parameters in β C and σ 2 in Figure 17 and find that our chain appears to be converging.