{'fraction_non_alphanumeric': 0.037202344458886136, 'fraction_numerical': 0.01132991637442676, 'mean_word_length': 4.837938439513243, 'pattern_counts': {'":': 0, '<': 6, '<?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': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Cross-validated Ensemble Approach to Robust Hypothesis Testing of Continuous Nonlinear Interactions: Application to Nutrition-Environment Studies

Gene-environment and nutrition-environment studies often involve testing of high-dimensional interactions between two sets of variables, each having potentially complex nonlinear main effects on an outcome. Construction of a valid and powerful hypothesis test for such an interaction is challenging, due to the difficulty in constructing an efficient and unbiased estimator for the complex, nonlinear main effects. In this work we address this problem by proposing a Cross-validated Ensemble of Kernels (CVEK) that learns the space of appropriate functions for the main effects using a cross-validated ensemble approach. With a carefully chosen library of base kernels, CVEK flexibly estimates the form of the main-effect functions from the data, and encourages test power by guarding against over-fitting under the alternative. The method is motivated by a study on the interaction between metal exposures in utero and maternal nutrition on children's neurodevelopment in rural Bangladesh. The proposed tests identified evidence of an interaction between minerals and vitamins intake and arsenic and manganese exposures. Results suggest that the detrimental effects of these metals are most pronounced at low intake levels of the nutrients, suggesting nutritional interventions in pregnant women could mitigate the adverse impacts of in utero metal exposures on children's neurodevelopment.

Investigation of the interplay between multiple lifestyle, biological, and environmental factors contributing to disease risk is a major goal in public health. Classic gene-environment and nutritionenvironment studies focus primarily on the interaction between discrete factors [34,9], or between discrete factors and the linear effect of a few continuous measurements (e.g. [31]). In recent years, however, recognizing the fact that populations are exposed to combinations of continuouslymeasured chemical and non-chemical factors that potentially have a nonlinear effect on outcome, there has been increasing interest in the best ways to statistically quantify the complex interplay of these continuous, nonlinear effects on health.
In this work we analyze data from a birth cohort study on the interaction between in utero exposure to a metal mixture and maternal nutrition intake during pregnancy on children's neurodevelopment in rural Bangladesh [16,23]. Bangladesh has been experiencing unparalleled levels of arsenic (As) other toxic metal poisoning through contaminated groundwater [36]. Bangladesh also has rates of undernutrition that are among the highest in the world [43]. A recent study [44] assessed the relationships between the arsenic (As), manganese (Mn), and lead (Pb) metal pollution mixture and infant neurodevelopment in Bangladesh, and has detected nonlinear, inverted-U shaped exposure-response relationships that differ among population subgroups, suggesting a role of additional cultural/behavioral factors, such as dietary habit, in affecting the impact of this metal mixture on children's health. One possible factor impacting these environmental effects is maternal nutrition during pregnancy. At vulnerable stages of fetal development, mother's overall nutrition intake may exacerbate adverse effects of chemical stressors. Specific nutrients may modify chemical effects because of their influence on the metabolism of the chemicals, on epigenetic programming in response to the chemicals, or through other mechanisms that vary by metal or by outcome. To answer this question, a companion nutritional study [29] was conducted to collect data on mother's nutrition intake during pregnancy, measuring the level of nutrition intake of 27 nutrients grouped in five nutrition categories (macronutrient, minerals, (pro-)vitamin As, vitamin Bs, and other vitamins), thereby providing an unique opportunity for researchers to quantitatively investigate the effect modification between nutrition intake during pregnancy and in utero metal exposures on infant development.
The Bangladesh study posed two challenges that are common in many modern data science applications: (1) high dimensionality of the interaction, as the interaction term contains second-and higher-order interactions between 27 nutrients and 3 metal exposures, and (2) the nonlinearity of the underlying exposure-response relationship, whose mathematical properties are unknown a priori.
In such a scenario, linear-model based methods are known to suffer from misspecification of the main effects model for nutrients (that include nutrient-nutrient interactions) and metals (including metal-metal interactions) even under the null of no nutrient-metal interactions, leading to inflated Type I error and reduced test power [41,9,49]. To boost efficiency and incorporate nonlinearities in the exposure-response relationship, a recent line of research has focused on constructing interaction tests based on kernel machine regression (KMR) [39,37]. Building on the success of previous kernel testing literature [32,48], these tests model the main-effect and interaction-effect functions as elements in reproducing kernel Hilbert spaces (RKHS) generated by pre-specified kernel functions, and build the hypothesis test by re-parametrizing the kernel machine regression as a linear mixed model [32]. In this framework, the interaction term is an additional random effect term controlled by an univariate garrote parameter, on which one can construct a variance-component score test [30] for a test of the null hypothesis of no interaction. Successful applications of such tests include targeted gene effect identification in genetic pathway analysis [33], gene-gene interaction detection in genome-wide association study [27], and also in gene-environment interaction studies with discrete factor such as gender [5] and risk indicators of cardiovascular disease [15].
Applications of interaction tests involving sets of multiple continuous measurements with nonlinear effects, however, remain rare. The key challenge impeding the success of interaction tests in continuous settings lies in designing a proper kernel function for the multi-dimensional, nonlinear main-effect functions of unknown form. The kernel functions for the main effect terms need to generate a RKHS that is rich enough to contain the main-effect functions under the null, while at the same time be sufficiently structured to maintain power for detecting interactions. Earlier work [32,33] approached this problem by selecting the kernel from an assumed parametric family (e.g. the Gaussian radial basis functions (RBF)) through maximum likelihood estimation, risking the specification of overly strong assumptions for these nonlinear functions. More recent approaches alleviate assumptions on the data-generation mechanism by incorporating multiple candidate kernels into the analysis, treating the kernel function as a weighted combination of candidate kernels, and learning kernel weights by maximizing various objective functions such as centered kernel alignment [12] or by an L 1 -regularized model likelihood [40]. However, designed primarily to maximize predictive accuracy, such procedures can be overly flexible under the alternative and potentially result in hypothesis tests with low power [51]. Permutation tests are another popular approach for alleviating the issue of kernel misspecification [8,51]; however, constructing a permutation procedure for an interaction test is usually not possible in observational studies, since the geneenvironment independence condition tends to not hold [7].
In this article, we propose a new approach to test for the interaction effect between groups of continuous features, each having potentially complex main effect functions relating outcome to that set of exposures. Built under the framework of kernel machine regression, we address the issue of kernel misspecification by deploying an ensemble of candidate kernels, and carefully design the ensemble strategy so that it minimizes the generalization error of the overall ensemble [14].
Consequently, the proposed test automatically estimates the form of the kernel under the null from the data and guards against overfitting the interaction effect under the alternative, resulting in a powerful test that is robust under a wide range of data generation mechanisms. As we discuss in Section 4, such a strategy results in an estimator that enjoys oracle property for ensemble selection and good generalization performance in limited samples, thereby achieving a powerful null-model estimator especially suitable for hypothesis testing in epidemiology studies. We term our method the Cross-Validated Ensemble of Kernels (CVEK). In Section 5, we illustrate the robustness of our method by conducting simulation studies that evaluate the finite-sample performance (Type-I error and power) under a range of data-generating scenarios and compare the performance of the proposed approach with other popular interaction tests. Finally, in Section 6, we apply our method to data from the Bangladesh reproductive cohort study [16,23] to investigate the interaction between mother's daily nutrient intake and in-utero exposure to an environmental metal mixture (As, Mn and Pb) on children's neurodevelopment.

Neurodevelopment in Bangladesh Children
The Bangladesh Reproductive Cohort Study (Project Jeebon) was initiated in 2008 to investigate the effects of prenatal and early childhood exposure to As, Mn and Pb on early childhood development.
During 2008-2011, pregnant female participants (with gestational age < 16 weeks) were recruited from two rural health clinics operated by the Dhaka Community Hospital Trust (DCH) in the Sirajdikhan and Pabna Sadar upazilas of Bangladesh. During 2008-2013, data were collected at five time points spanning the entire perinatal and early childhood period, including: initial clinic visit (gestational age < 16 weeks, Visit 1); pre-delivery clinic visit (gestational age = 28 weeks, Visit 2), time of delivery (Visit 3), post-delivery clinic visit (infant age less than 1 month, Visit 4), and a post-natal follow-up visit (infant age between 20-40 weeks, Visit 5). Our central hypothesis is that children born from mother who had lower nutrient intake will be the most susceptible to adverse effects of metal exposures.
Detailed procedures for data collection and measurement protocols have been documented previously [16,23,44]. Briefly, background information on parent's demographic status, including age, education, smoking history and socioeconomic status were collected through structured questionnaires at the two clinic visits during pregnancy (Visits 1-2). Information on infant's biometric measurements, including sex, birth weight, length, head circumference, birth order and gestational age, were recorded at birth. Information on maternal medical history, maternal depression status (in Edinburgh Depression scale), maternal IQ (assessed using the Raven's Progressive Indices [38]) were measured during the pregnancy visits (Visits 1-2), and an infant's early childhood development, medical history, and quality of home environment (in terms of emotional, social, and cognitive stimulation, measured by HOME instrument score [3] were measured during the follow-up visits (Visits 4-5), respectively.
Each infant's exposure to multiple metals As, Mn and Pb (concentrations in µg/dL) during pregnancy were measured using blood samples from infant's umbilical cord venous blood collected at the time of the birth. Mother's overall nutrition intake status during pregnancy were measured for 27 nutrients derived from semi-quantitative Food Frequency Questionnaires (FFQs) specially adapted to Bangladeshi diet [28] at both the pre-and post-delivery visits (Visit 2 and 4). This instrument derives data on these 27 nutrients from measures of the consumption frequency (amount per week) of 42 food items during the 12-month period preceding delivery. The nutrients measured can be grouped into 5 categories including macro-nutrients (5 nutrients: protein, fat, carbohydrate, dietary fiber and ash), minerals (8 nutrients: calcium, iron, magnesium, phosphorus, potassium, sodium, zinc and copper), vitamin A and provitamin As (6 nutrients: vitamin A, retinol, beta-carotene equivalents, alpha-carotene, beta-carotene, and cryptoxanthin), vitamin B (5 nutrients: thiamin (B1), riboflavin (B2), niacin (B3), vitamin B6 and folate (B9)), and other vitamins (3 nutrients: vitamin C (i.e. L-ascorbic acid), vitamin D, and vitamin E). Finally, infant's neurodevelopmental outcomes were assessed at 20-40 months of age (Visit 5) using a translated and culturally-adapted version of the Bayley Scales of Infant and Toddler Development, Third Edition (BSID-III) including five cognitive domains: cognitive, receptive language, expressive language, fine motor and gross motor.

Model and Inference
Assume we observe data from n independent subjects. For the i th subject, let y i be a continuous response, x i be the set of p baseline covariates that can be entered into the model linearly, and z i be the set of q continuous covariates that have a nonlinear effect on y i . Furthermore, we assume that there exists a grouping structure among the z i covariates such that z i = {z 1,i , z 2,i }, where the m th group z m,i ∈ R q m contains q m covariates, m = 1, 2. We discuss the generalization to the case of more than two groups in z i in Section 3.2.
We assume that the outcome y i depends on covariates x i , z i through the model: where β is a p × 1 vector of unknown coefficients for background covariates, h(z i ) : R q → R is an unknown continuous function describing the effect of z i , and ε i is random noise that is independently and identically distributed as N(0, σ 2 ). For identifiability purpose, h is assumed to be square-integrable and subject to the constraint R q h(z)dz = 0.
Our main objective in this work is to test for the interaction between two chosen sets of covariates in z i = {z 1,i , z 2,i }, while accounting for interactions within each covariate set. Without loss of generality, consider testing for the interaction between z 1,i and z 2,i . Then our hypothesis is: where H 12 is the space of "pure interaction" functions that contain only the interaction effect between (z 1,i , z 2,i ). That is, under the null hypothesis, h(z) may depend on the individual main effects of z 1,i , z 2,i , but does not depend on the interaction effect of the set pair (z 1,i , z 2,i ).
We take the penalized likelihood approach to estimate parameters (β , h). Namely, we first specify H the candidate space and λ the penalty parameter, then estimate parametersθ = (β ,ĥ) by minimizing the penalized negative log likelihood: We model H using Kernel Machine Regression (KMR) [39]. Specifically, we assume H to be a Reproducing Kernel Hilbert Space (RKHS) generated by a positive-definite kernel function k(z i , z i ), such that any h ∈ H can be expressed in terms the kernel function as f ( Then by the Representer theorem [6], if we define y n×1 = [y 1 , . . . , . . , α n ] T and also denote K n×n the kernel matrix with its (i, j) th element to be K i, j = k(z i , z j ), then (3) can be re-written as Furthermore, if we define τ = σ 2 λ , h * can arise exactly from a linear mixed model (LMM) [32]

A Variance Component Test for Kernel Interaction
Under the LMM formulation of Kernel Machine regression in (5), Maity and Lin ([33]) built a general test for the hypothesis H 0 : h ∈ H 0 by assuming that h lies in a RKHS generated by a garrote kernel function k δ (z, z ), which is constructed by attaching an extra garrote parameter δ to a regular kernel function. When δ = 0, the garrote kernel function k 0 (z, z ) = k δ (z, z ) δ =0 generates exactly H 0 the space of functions under the null hypothesis. The authors further proposed a REML-based variance component score test for H 0 . In order to adapt the above approach to the hypothesis for interaction H 0 : h ∈ H ⊥ 12 , we construct the garrote kernel function k δ (z, z ) by building its corresponding RKHS for the main-effect and interaction space using the tensor-product construction [18]. Briefly, let 1 = { f | f ∝ 1} be the RKHS of constant functions with kernel function k(z, z ) = 1, and let H m be the RKHS of centered functions (i.e. f (z m ) dz m = 0) with domain on the m th covariate set z m , then H adopts below orthogonal decomposition: where H 0 = H 1 ⊕ H 2 is the space of main-effect functions that does not contain the (z 1,i , z 2,i ) interaction, and H ⊥ 12 = H 1 ⊗ H 2 is "pure interaction" space whose elements contain only the interaction effect between (z 1,i , z 2,i ). Correspondingly, the garrote kernel function for interaction can be constructed as Under the above form of the garrote kernel function, the element of the null derivative kernel matrix i.e. the null derivative kernel matrix ∂ K 0 is simply the kernel matrix K 12 that corresponds to the interaction space. Therefore the test statistic is: The null distribution ofT 0 is a mixture of chi-squares that can be approximated using a scaled chi-square distribution κ χ 2 ν using either Satterthwaite-Welch method [50] or other higher-moment approximations [4].

Generalization to Multiple Groups with Nuisance Interaction
Our description so far assumes there exists no nuisance interaction terms in the model y = x T β + h(z) + ε. However, in more realistic scenario, z usually exhibits complex hierarchical structure subsuming multiple groups, and it is often of interest to test only for the interaction between two small subgroups of z, leaving other interactions as nuisance effect to be accounted for by the null model. For example, consider the case of nutrition-environment interaction in Bangladesh birth cohort, z i is the 30 × 1 vector of during-pregnancy exposure to 27 nutrients and 3 metal pollutants, corresponding the grouping structure z i = {z metal , z nutr }, where z nutr is further divided into Therefore, when testing for the interaction between metal mixture exposures and a specific nutrient group of interest, care should be given to formulate h(z i ) such that it not only explicitly characterizes the interaction of interest, but also account for all nuisance interactions among other z i subgroups.
More specifically, assume z i = {z 1 , z 2 , z 3 }, when testing for the interaction between z 1 and z 2 , the nuisance interactions between z 1 and z 3 , as well as between z 2 and z 3 , should also be included in the null model. To this end, following the tensor-product construction shown in Section 3.1, we adopt the following orthogonal decomposition of h(z): Under such construction, the null hypothesis of no interaction corresponds to h 12 and h 123 equaling zero, i.e. 23 and k a = k 12 + k 123 . Here the k 1 , k 2 , and k 3 are the reproducing kernels for the main-effect space of z 1 , z 1 , z 3 , respectively, and the higher-order interaction kernels are constructed as ∀(i, j), k i j = k i * k j and k 123 = k 1 * k 2 * k 3 similar to Section 3.1. Consequently, denoting K i as the kernel matrix corresponding to k i , the null kernel matrix K 0 and the interaction kernel matrix K 12 are and the test statistic can be constructed as in (7).

Robust Effect Estimation using Cross-validated Ensemble
We motivate the importance of robust null model estimation by considering the possible impact of a misspecified null kernel function k 0 on the performance of the resulting hypothesis test. Specifically, we express the test statisticT 0 in (7) in terms of the model residualε = y −μ −ĥ: where we have used the fact [20]. Therefore, the test statisticT 0 is a scaled quadratic-form statistic that is a function of the model residual. If k 0 is too restrictive, model estimates will underfit the data under the null hypothesis, introducing extraneous correlation among theε i 's that yield inflatedT 0 values and deflated p-values under the null. Therefore, this approach will yield an invalid test having inflated Type I error. On the other hand, if k 0 is too flexible, model estimates will likely overfit the data in small samples, producing underestimated residuals, which leads to underestimated test statistics and overestimated p-values. Accordingly, the resulting test will have low power.
The above observations motivate a kernel estimation strategy that is flexible in that it does not underfit under the null, yet stable so that it does not overfit under the alternative. To this end, we propose estimating h using the convex ensemble of a library of fixed base kernels {k d } D d=1 : whereĥ d is the kernel predictor generated by d th base kernel k d . In order to maximize model stability, the ensemble weights u are estimated to minimize the overall cross-validation error ofĥ.
We term this method the Cross-Validated Kernel Ensemble (CVKE). The exact algorithm proceeds in three stages as follows (see Algorithm 1 for summary): Stage 3: Ensemble Kernel Matrix Using the ensemble hat matrixÂ, estimate the ensemble kernel matrixK by solvingK(K + λ I) −1 = A. Specifically, if we denote U A and {δ A,k } n k=1 as the eigenvector and eigenvalues of A, respectively, then the ensemble kernel matrix K adopts the form : where we recommended setting We remind readers that the CVEK's ensemble form (Step 2) belongs to the general class of model aggregation method known as convex aggregation [42], whose oracle property in model selection has been established both asymptotically and in finite-sample [45,26]. It can be also considered as a special case of ensemble of kernel predictors (EKP) [11], whose generalization behavior is well characterized in terms of the rate of eigenvalue decay of the base kernels. Consequently, under the null hypothesis, with a diverse set of base kernels, the CVEK ensemble converges in O( 1 n ) rate to the "oracle ensemble" made by an oracle that has access to infinite amount of validation data, thereby resulting in correct Type I error by mitigating null model misspecification. Under the alternative, by setting the diverse kernel library to be a mix of parametric kernels (linear, polynomial) and smooth kernels of exponential eigendecay rate (e.g. a collection of Gaussian RBF kernel with different fixed spatial smoothness parameters), CVEK converges to its asymptotic counterpart in O( 1 n ) rate if the data-generation function is indeed parametric, and in the "near-parametric" rate of O( log(n) n ) if the data-generation function is complex and nonlinear, thereby encouraging good test power by not overfitting the interaction effect due to fast generalization rate. The resulting ensemble kernel is therefore a strong candidate for a null model estimator that is suitable for hypothesis testing. We refer readers to Supplementary Section A.2 for detailed discussion.

Numeric Studies
We evaluate the finite-sample performance of the proposed interaction test in a simulation study that mimics a small-sample nutrition-environment interaction study. We generate two groups of input features (z i,1 , z i,2 ) ∈ R p 1 × R p 2 independently from a multivariate Gaussian distribution, reflecting each subject's level of exposure to p 1 environmental pollutants and the levels of a subject's intake of p 2 nutrients during the study. Throughout the simulation scenarios, we keep n = 200, and We generate the outcome y i as: where h 1 , h 2 , h 12 are sampled from RKHSs H 1 , H 2 and H 1 ⊗ H 2 , generated using a ground-truth kernel k true . We standardize all sampled functions to have unit norm, so that δ represents the strength of interaction relative to the main effect. For each simulation scenario, we first generate data using δ and k true , and then use a k model to estimate the null model and obtain p-value using the proposed test. We repeat each scenario 1000 times, and evaluated the test performance using the empirical probability P(p ≤ 0.05).
In this study, we vary k true to produce data-generating functions h δ (z i,1 , z i,2 ) with different smoothness and complexity properties, and varied k model to reflect different common modeling strategies for the null model in addition to using CVEK. We then evaluate how these two aspects impact the hypothesis test's Type I error and power. More specifically, we sample the data-generating function using k true from Matérn kernel family [37]: with two non-negative hyperparameters (ν, σ ). For a function h sampled using a Matérn kernel, ν determines the function's smoothness (i.e. degree of mean-square differentiability), and σ determines the function's complexity in terms of spectral frequency [37]. In this work, we vary ν ∈ { 3 2 , 5 2 , ∞} to generate once-, twice, and infinitely-differentiable functions, and vary σ ∈ {0.5, 1, 1.5} to generate functions with varying degree of complexity.
We consider 12 k model 's that are grouped into five model families (See Table 1 for a complete summary): (1) Polynomial Kernels that is equivalent to polynomial ridge regression. In this work, we use the linear kernel k linear (x, x |p) = x T x and quadratic kernel k quad (x, x |p) = (1 + x T x ) 2 .
(2) Gaussian RBF Kernels: The simulation results are presented graphically in Figure 1 and documented in detail in the Supplementary Section A.3. We first observe that for reasonably specified values of k model , the proposed hypothesis test with CVEK estimator always has the correct Type I error and reasonable power. We also observe that the complexity of the data-generating function h δ (11) plays a role in test performance, in the sense that the power of the hypothesis tests increases as the Matérn k true 's complex parameter σ becomes larger, which corresponds to functions that put more weight on the simpler, slow-varying eigenfunctions in Bochner's spectral decomposition [37].
There exist clear differences in test performance between different model families. In general, polynomial models (linear and quadratic) appear to be too restrictive and underfit the data under both the null and the alternative, producing inflated Type I error and diminished power. On the other hand, lower-order Matérn kernels (Matérn 1/2 and Matérn 3/2, dark blue lines) appear to be too flexible, due to their slow eigenvalue decay. Whenever data are generated from similarly or smoother k true , Matérn 1/2 and 3/2 overfit the data and produce deflated Type I error and severely diminished power, even if the hyperparameter σ is fixed at its true value. Comparatively, Gaussian RBF works well for a wider range of k true 's, but only if the hyperparameter σ is selected carefully. Specifically, RBF-Median (black dashed line) works generally well, despite being slightly conservative (i.e. lower power) when the data-generation function is smooth and of low complexity. RBF-MLE (black solid line), on the other hand, tends to underfit the data under the null and exhibits inflated Type I error, possibly because of the fact that σ is not strongly identified when the sample size is modest [47]. The situation becomes more severe as h δ becomes rougher and more complex. In the more extreme case of non-differentiable h with σ = 1.5, the Type I error is inflated to as high as 0.24. Neural Network kernels also perform well for a wide range of k true , and with the Type I error more robust to the specification of hyperparameters.
Finally, the two ensemble estimators CVEK-RBF and CVEK-NN perform as well or better than the non-ensemble approaches for all k true 's, despite being slightly conservative under the null. As compared to CVEK-NN, CVEK-RBF appears to be slightly more powerful, validating recommendations in Section 4.

Data Analysis
We use the proposed methods to test for nutrition-environment interactions in the Bangladesh birth cohort study (see Section 2 for complete study description). Our aim is to detect whether mother's nutrient intake during pregnancy modifies the effect of metal mixture exposures on children's early-stage fine motor BSID-III scores in the district of Pabna (n = 351). Specifically, our interest concentrates on the interaction between the mixture of As, Mn and Pb, and five major nutrient groups: macro-nutrient, mineral, vitamin A, vitamin B and the other vitamins (denoted as vitamin, other). For each of the five nutrient groups, we test for the overall interaction between the selected group and the joint effect of the As, Pb, Mn mixture. We adjust for parent's demographic status (age, education, smoking history), infant's biometric measurements at birth (sex, birth weight, length, head circumference, birth order and gestational age), and quality of early-childhood home environment (HOME score, maternal depression scale, maternal IQ).
We compare our method with three existing approaches for testing high-dimensional interaction. can be considered as a special case of kernel machine regression with the kernel matrix constructed adaptively using spline basis functions [24]. In order to visualize the identified interaction and thereby provide interpretable findings, we graphically summarize the multivariate interaction effects by examining the joint exposure-response surface between the principal components of the pollutant mixture and those for each nutrient group. Table 2 presents p−values for the interaction between the overall metal mixture and each of the five nutrient groups. We conducted the proposed test using two types of CVEK models: an ensemble of seven neural network kernels with prior variance set between log(σ 2 ) ∈ {−3, −2, −1, 0, 1, 2, 3}

Nutrient -Mixture Interactions
(denoted as CVEK-NN), and an ensemble of seven RBF kernels with bandwidth parameter set between log(σ ) ∈ {−3, −2, −1, 0, 1, 2, 3} (denoted as CVEK-RBF). We compared the results of each to those generated by iSKAT, GKM, and GE-Spline. As shown in the table, most tests yielded strong evidence of interaction (p < 0.05) for vitamin A and the vitamin, other groups, as well as moderate evidence of interaction (p <≈ 0.1) for the mineral and the vitamin B groups. There was no evidence of an interaction between metal exposures and macro-nutrients.
Comparing the performance across different tests, we observed similar patterns for p-values for CVEK-NN and CVEK-RBF, suggesting robustness in test performance with respect to the choice of the family of the base kernels. We also observed similar patterns for p-values between iSKAT (linear kernel) and GKM, the latter of which used a single RBF kernel with REML-based hyperparameter tuning. The results from these two tests are similar to those from the CVEK tests for vitamin A and vitamin, other groups of nutrients, but are less powerful in detecting the interaction for the mineral and the vitamin B groups. This is consistent with our observation in Section 5 that, when the true effect is nonlinear and exhibits a moderate level of smoothness and low complexity (a scenario that is likely to hold for the effect of environmental exposures, see Figure 1 (a)), the hypothesis test based on GKM performs similarly to that based on the iSKAT but is less powerful than the CVEK-based test. This reduction in power is possibly due to the overly strong smoothness assumption imposed by these two models. Finally, we notice that the performance of the test from the GE-spline model appears sub-optimal when compared to that of the other methods. GE-spline produced a much smaller p−value for the interaction for the A vitamins (p = 0.0084), but much higher p−values for the other nutrient groups, failing to detect the interaction for the B vitamins and the other vitamins.
We hypothesize that the observed instability of GE-spline is likely caused by the lack of fit of the null model, due to the difficulty in estimating multivariate splines in high dimensions.

Visualization of Exposure-Response Surface
To better understand the nature of the multivariate interactions between the environmental exposures and nutrition, in Figure 2 Figure 2 (a)) in the overall interaction is expectedly much stronger than that between the 3 rd PCs (e.g. Figure 2 (i)). This explains the lack of significant evidence of overall interaction for the macronutrients in Table 2, since the joint effect between the PCs accounting for more variance (e.g.  Table 2. Finally, we observe that across all nutrient groups, the nutrient PCs interacts the most often with the 1 st pollutant PC, which is strongly associated with As, and also with the 3 rd pollutant PC, which is strongly associated with Mn, suggesting that As and Mn are the two main pollutants driving the overall interaction. Furthermore, the pattern of interaction between nutrient and pollutant PCs are observed to be similar across nutrient groups: at lower levels of nutrients (x-axis), higher levels of metal exposure (y-axis) is associated with lower neurodevelopment scores. At intermediate or high levels of the nutrient, however, this negative association either disappears (see, e.g. Figure 3 (c), (d), (h)) or even becomes positive (see, e.g.

Discussion
Under the framework of kernel machine regression, we have developed a hypothesis testing procedure for detecting nonlinear interactions between groups of continuous covariates. In this context,  [13]. Furthermore, the ensemble weights {u d } D d=1 (see (9)) in CVEK were estimated to maximize the estimator's cross-validation stability. The optimality of such method in terms of the power of the hypothesis test has not been fully investigated. It is desirable to develop an optimal estimation procedure for the ensemble weights {u d } D d=1 that maximizes the power of the hypothesis test, in a manner similar to [17]. Given such a procedure, it is also of theoretical interest to compare the difference between the ensemble weights generated by maximizing crossvalidation stability to those generated by maximizing the power of the test in both finite samples and asymptotically. Finally, due to the abundance of non-Gaussian and non-continuous outcomes in the epidemiology and medical literature, extension of the kernel ensemble approach to discrete or censored continuous outcomes is also of interest.

Kernel Function
Model Name Parameter Value

A.1 Derivation for Ensemble Kernel Matrix
Given the ensemble hat matrix A in Section 4, we consider how to identify the ensemble kernel matrix K by solving: ) the eigenvector and eigenvalues of A and K, respectively, then the above system reduces to: . Therefore the ensemble kernel matrix K is estimated as: Choice of ensemble tunning parameter λ K Notice that we have left the "ensemble tunning parameter" λ K unspecified. In practice, λ K serves only as a constant scaling factor for the kernel matrix K, whose exact value does not impact either the prediction or the p-value calculation, since both procedures are scale invariant with respect to the kernel matrix. Therefore it can be set to a value of our choice. One common choice for λ K is to set ) −1 such that tr( K) ≤ 1, this is because the Rademacher complexity of the overall ensemble can be upper-bounded as a function of tr( K) [25]. Another interesting choice ) . Intuitively, this means the tunning parameter for ensemble kernel matrix K should grow in the same rate as the tunning parameter for the best-peforming base kernel, as the ensemble kernel matrix is expected to perform as well or better than the best-performing base kernel. Further, such choice provides guarantee on the generalization performance of the ensemble by bounding the ensemble kernel matrix's decay rate of the tail sum of the eigenvalues [35].

A.2.1 Oracle Selection
We remind readers that the CVEK's ensemble strategy (Step 2) is similar to that of Jackknife Model Averaging (JMA) [19] in the case of leave-one-out cross validation and the Super Learner [46] in the case of K-fold cross validation. Its oracle property in model selection has been established both asymptotically and in finite-sample [45]. That is, on average, given the same set of estimated kernel

A.2.2 Generalization
In limited samples, the oracle property does not ensure a powerful test, and the estimator's generalization property, i.e. how fast the estimator's behavior approaches its asymptotic counterpart, must also be taken into consideration to guarantee good power under the alternative. This is because even under a oracle selector, an ensemble of flexible estimators in finite samples will still overfit the To this end, we draw upon a classical result from the RKHS literature that the generalization error of a kernel-based estimator is bounded above by its RKHS's local Rademacher complexity [2], a measure of the richness of the class of candidate functions that is characterized by the rate of eigenvalue decay of its kernel function [35]. Consequently, since the learned ensemble estimator h = ∑ D d=1 u d h d lies in convex combination of the RKHSs generated by the base kernels, the generalization property of the CVEK estimator is explicitly characterized by the generalization property of the "selected" base estimators (i.e. theĥ d assigned non-zero weight by the ensemble) in terms of their respective rates of eigenvalue decay. Moreover, it can be shown that for special classes of kernels, the generalization error rate can indeed be improved upon: finite-rank kernels (e.g. linear and polynomial kernels) are able to achieve O( 1 n ), while kernel families with exponential rate eigenvalue decay (e.g. Gaussian RBF kernel) can achieve a rate of O( log(n) n ) [10].

A.2.3 Choice for Base Kernels
In combination of the CVEK's oracle selection property, the above result suggests that, if there exists sets of parametric kernels in the library, CVEK behaves as a parametric model by achieving a rate of O( 1 n ) if the data-generation function is indeed parametric. For more complex data-generation mechanism, CVEK is able to achieve a rate of O( log(n) n ) by including sets of Gaussian RBF kernels in the library. Consequently, we recommend practitioners to construct the kernel library with a mix of parametric kernels (linear, polynomial) and smooth kernels with exponential rate in eigendecay (e.g. a collection of Gaussian RBF kernel with different fixed spatial smoothness parameters). We notice that despite the extreme smoothness of these base kernels, the resulting ensemble is in fact very flexible and hence does not risk underfitting the data generation function. This is because by Bochner's theorem, fitting the model with an ensemble of Gaussian RBF kernels is equivalent to approximating the spectral density of the function with convex mixtures of Gaussian densities, which is shown to be capable of approximating arbitrary continuous density with compact support in R d [1]. More generally, practitioners should be careful about using an ensemble containing only flexible kernels, since in order to represent a large function space, these kernel functions usually have heavy-tailed spectral densities and therefore slow eigenvalue decay [35]. Finally, we notice that there also exists interesting kernel families that lie outside the usual scope of theoretical analysis, but tend to do well in practice. One such example is the Neural Network kernel [37]. We investigate the performance of these kernels and compare to the performance of the Gaussian RBF ensemble in Section 5.

A.3 Detailed Simulation Results
In this section we document the value of estimated P(p < 0.05) from the simulation presented in Section 6 (Simulation Experiment) of the paper. Recall that the simulation data is generated from below mechanism: where h i 's are functions with unit norm sampled from the reproducing kernel Hilbert spaces (RKHSs) generated by k true , and the data is then fitted using Gaussian process with k model .
Each table documents the P(p < 0.05) resulted from fixing k true to a Matérn kernel with specific value of smoothness parameter ν and complexity parameter σ , and then varying the strength of the interaction δ ∈ [0, 1] and the model kernel k model .
Our general observations are: 1. The value of test power increases as the value k true 's complexity parameter σ becomes larger. This is possibly caused by the fact that the interaction becomes easier to detect as the pure interaction function h 12 ∈ H 12 becomes more complex as in it varies more quickly.
2. Given the data-generation mechanism: (a) Polynomial kernels (Linear and Quadratic kernels) exhibits underfit, and result in inflated Type I error but also low power. (d) Neural Network kernels also work well in general. Their performance is also impacted by the hyperparameters. However not as sensitive as Gaussian RBF.