# Semiparametric Discovery and Estimation of Interaction in Mixed Exposures
using Stochastic Interventions

David B. McCoy
Division of Biostatistics, University of California, Berkeley
and
Alan Hubbard
Division of Biostatistics, University of California, Berkeley
and
Mark van der Laan
Division of Biostatistics, University of California, Berkeley
and
Alejandro Schuler
Division of Biostatistics, University of California, Berkeley
The authors gratefully acknowledge funding for Core E of the NIEHS Superfund
Center at Berkeley funded by NIH grant P42ES004705

###### Abstract

This study introduces a nonparametric definition of interaction and provides
an approach to both interaction discovery and efficient estimation of this
parameter. Using stochastic shift interventions and ensemble machine learning,
our approach identifies and quantifies interaction effects through a model-
independent target parameter, estimated via targeted maximum likelihood and
cross-validation. This method contrasts the expected outcomes of joint
interventions with those of individual interventions. Validation through
simulation and application to the National Institute of Environmental Health
Sciences Mixtures Workshop data demonstrate the efficacy of our method in
detecting true interaction directions and its consistency in identifying
significant impacts of furan exposure on leukocyte telomere length. Our
method, called SuperNOVA, advances the ability to analyze multiexposure
interactions within high-dimensional data, offering significant methodological
improvements to understand complex exposure dynamics in health research. We
provide peer-reviewed open-source software that employs or proposed
methodology in the SuperNOVA R package.

Keywords: Targeted Maximum Likelihood Estimation, Mixtures, Interactions,
Decision Trees, Ensemble Learning

##  1 Introduction

In environmental health research, understanding the impact of simultaneous
multiple exposures on health outcomes is crucial yet challenging. Traditional
epidemiological studies frequently rely on generalized linear models (GLM) to
assess main effects, often neglecting complex interactions within exposure
mixtures, potentially oversimplifying the intricate dynamics at play wang2020,
García-Villarino et al. [2022]. While methods such as Principal Component
Analysis (PCA) and penalized regressions, such as LASSO and Ridge, address
issues of high dimensionality and multicollinearity, they fall short in
interpretability, especially in regulatory and policy contexts where specific
guidance on exposure reduction is needed Roberts and Martin [2006],
XinWangBhramarMukherjee2018, McEligot et al. [2020]. The limitations of
existing approaches, including quantile sum g-computation and linear models
without interaction terms, risk simplifying the nuanced effects of
interactions Keil et al. [2019]. Bayesian Kernel Machine Regression (BKMR)
offers a more nuanced framework employing Bayesian hierarchical models and
kernel machine learning to evaluate high-dimensional exposure effects,
providing uncertainty quantification through posterior distributions Bobb et
al. [2014]. Despite BKMR’s advantages, practical challenges emerge,
particularly with defining interactions in high-dimensional settings (which
two or three-way interactions to examine), which lead to combinatorial
complexity and prioritization challenges for researchers. This backdrop
underscores the necessity for our proposed method, which aims to accurately
identify and estimate interactions in the context of environmental health
research, thereby addressing the limitations of current methodologies.

Given that interactions in a mixed-exposure setting are not known a priori,
data-driven approaches for interaction identification are critical. The
functional ANOVA (fANOVA) technique, for example, developed to detect
interactions in flexible models, such as tree-based or additive models, is one
of such approaches. Lengerich et al. [2020]. However, fANOVA, like traditional
ANOVA, primarily quantifies interaction effects in terms of the variance
explained, which may not directly translate into actionable insights for
epidemiologists or policy makers. Similarly, methods such as Narisetty et al.
[2019] use a forward selection of expanded basis functions to find
interactions. In both cases, although modeling informs what interaction
improves model fit, there is no measure of causal association, which is more
pertinent to inform policy decisions. Important insights have been provided on
how to define interactions in a causal inference framework VanderWeele [2009],
VanderWeele and Knol [2014], VanderWeele and Robins [2007], Rothman et al.
[1980]. Much of this work has done so in the context of binary explanatory
variables of interest. However, the wide range of exposures that characterize
environmental investigations often render standard binary causal models
inadequate. This accentuates the need for a versatile, semiparametric
interaction definition that can accommodate a spectrum of data types, a
current gap in methodologies.

Our proposed approach, termed SuperNOVA, synthesizes two critical phases of
analysis, discovery and estimation, using the same data. Initially, we used an
ANOVA-style decomposition to suggest interacting exposures in a semiparametric
model via a sieve approach. Subsequently, we estimate the counterfactual
effect of changing all of these exposures by a small amount versus the
combined effects of changing them individually. The interaction discovery
component of SuperNOVA is similar to the fANOVA approach for the discernment
of interactions. However, our approach finds the best finite sample
approximation in a model using ensemble learning of semiparametric maximum
likelihood estimation regressions van der Laan Mark et al. [2007]. This leads
to ANOVA over basis functions, which, when aggregated, offers holistic
F-statistics sets for pivotal exposure sets used in the best estimator. These
exposure sets are then the target for our nonparametrically defined
interaction parameter. Our target parameter evaluates the effect of feasible
stochastic shift interventions on exposure distributions, as conceptualized in
the literature Iván Díaz Muñoz and Mark van der Laan* [2012]. Such
interventions are highly relevant to public policy, as they often mirror the
realistic impact of regulatory measures that aim to slightly reduce, rather
than completely eliminate or uniformly standardize, exposures across a
population. For example, we might consider the population-level effect of
reducing everyone’s exposure to a persistent organic pollutant by a small,
relative amount. Our parameter defines an interaction effect by comparing the
expected outcome under a joint shift of exposures (like simultaneously
reducing levels of a dioxin and PCB) against the summed expected outcomes of
individual shifts. This approach offers a pragmatic and policy-relevant
interpretation of interactions, grounded in causal assumptions, and is
estimable from observational data using an asymptotically efficient targeted
maximum likelihood estimator that facilitates the construction of confidence
intervals and hypothesis testing.

The concurrent requirements of identifying (initially unknown) interacting
variables and estimating their effects require a two-prong approach to prevent
obvious overfitting. To address this, we split the full sample data: one
segment identifies interactions and the other employs an efficient estimator
for the data-adaptive target parameters. Our approach, termed SuperNOVA
integrates estimates across different data splits with cross-validated
targeted maximum likelihood (CV-TMLE) Hubbard et al. [2016], Zheng and van der
Laan [2010]. SuperNOVA stands distinctively due to several key advantages.
Firstly, it avoids bias of assuming linearity unlike techniques such as PCA or
lasso. Although fANOVA and BKMR also embrace nonparametric approaches,
SuperNOVA sets itself apart through its data-driven discovery of interactions.
Contrary to approaches like fANOVA and BKMR, which estimate the joint impact
of all mixture variables, SuperNOVA offers a more practical solution to
analyze high-dimensional mixtures in finite samples. It integrates a variable
selection mechanism, focusing exclusively on those variables within the
mixture that demonstrate statistically significant associations with the
outcome of interest. This targeted approach in SuperNOVA not only improves
efficiency, but also provides clearer insights into the key variables driving
interactions in complex exposure scenarios. Additionally, its policy-relevant
definition of interaction clearly marks an improvement over the aforementioned
methods. The method’s scalability and speed of inference further bolster its
appeal. By leveraging standard ML/TMLE tools, SuperNOVA avoids the pitfalls of
bespoke optimization algorithms, such as those in fANOVA or time-consuming
MCMC as seen in BKMR. In SuperNOVA, the use of data-adaptive estimators in
both interaction identification and nuisance estimation is key. Moreover, the
application of CV-TMLE ensures robust inference using the full dataset,
thereby retaining power even when data is used for both parameter definition
and estimation. However, in the case where the same interactions are not found
in every fold of the CV procedure, the inferred null must be incorporated into
a pooled result. We offer an approach to estimate pooled TMLE results for
inconsistent interactions, which has higher power compared to simply reporting
fold-specific results.

The structure of this paper is as follows. Section 2.1 provides a background
on the semiparametric methodology that underpins stochastic interventions.
Sections 2.2 and 2.3 delineate the interaction target parameter and the
essential assumptions for our causal interpretations. Section 3. discusses the
TMLE estimator. In Section 4. we discuss inference and elucidate our approach
to determining variable sets. Section 5. presents simulation results,
underscoring the asymptotic unbiasedness of our estimates. Section 6
demonstrates SuperNOVA’s applicability, juxtaposing it against existing
techniques and its relevance in real-world scenarios using synthetic data from
the National Institute of Environmental Health Sciences (NIEHS). We provide
real-world results when applying SuperNOVA to National Heath and Nutrition
Examination Survey, where we investigate the effects of 18 POPs on telomere
length. In Section 7 we wrap up with a comprehensive discussion of strengths
and weaknesses and provide our software.

##  2 Data and Parameter of Interest

###  2.1 Overview

At a high level, consider that an analyst has data on a mixture of 15
chemicals (e.g. 15 POPs). Our approach first splits the entire data into
discovery and estimation samples. In the discovery sample, we employ versatile
machine learning to identify interactions. For example, consider that we find
an interaction between a POP-1 and POP-6. We then use the estimation portion
of the data to estimate our interaction target parameter. Our interaction
parameter is based on comparing a joint stochastic shift, that is, shifting
both exposure simultaneously and observing the expected outcome and comparing
this to the summed expected outcomes given individual shifts of each POP
independently. We cycle through the data, for example, 10 times, and average
estimates across the 10 folds using cross-validated targeted maximum
likelihood estimation (CV-TMLE). Because our estimand depends on the data, if
the full data were used to both discover and estimate, our statistical
inference suffers from overfitting; thus, this data-splitting is necessary.
With this example in mind, subsequent sections will delve deeper, elucidating
each component of our parameter for enhanced clarity.

###  2.2 Notation and Framework

Building upon the concepts presented in Iván Díaz Muñoz and Mark van der Laan*
[2012] on stochastic interventions, we generalize this work to a target
parameter that compares the expected outcome under a joint shift (here we
focus on two exposures for brevity). As such, we first describe our target
parameter for a stochastic shift of multiple exposures (which deviates little
from estimating the impact of a stochastic shift of an individual variable)
and then show how we can build our target parameter for interactions from this
joint-shift target parameter. Consider an experiment in which a vector of
exposure variables $\mathbf{A}$, a continuous or binary outcome $Y$, and a set
of covariates $W$ are measured for $n$ randomly sampled subjects. Let
$O=(W,\mathbf{A},Y)$ represent a random variable with distribution $P_{0}$,
and $O_{1},...,O_{n}$ represent $n$ i.i.d. observations of $O$. The true
distribution $P_{0}$ of $O$ can be factorized as
$P_{0}(O)=P_{0}(Y|\mathbf{A},W)P_{0}(\mathbf{A}|W)P_{0}(W)$. We use $p(\cdot)$
to denote the corresponding densities, e.g. $p(Y|A,W)$. The conditional joint
density of exposures given covariates $p(\mathbf{A}|W)$ we will refer to as
$g_{0}(\mathbf{A}|W)$ to distinguish it. Similarly we use
$\bar{Q}_{0}(\mathbf{A},W)\equiv E_{0}(Y|\mathbf{A},W)$. We denote a “shifted”
exposure density using the notation
$g_{\delta}(\mathbf{A}|W)=g_{0}(\mathbf{A}-\boldsymbol{\delta}|W)$ and we
denote the full “shifted” density over all variables as
$p_{\delta}(O)=p(Y|\mathbf{A},W)g_{\delta}(\mathbf{A}|W)p(W)$. In this
section, while we theoretically describe the shifts $\boldsymbol{\delta}$ as
fixed quantities (in particular, which variables are shifted at all and in
what amounts), practical applications require them to be identified in a data-
adaptive manner. We first describe the parameter that we will estimate
treating exposure shifts as fixed, and then we describe how to identify what
exposures to shift and by how much.

For example, consider $A=(A_{1},A_{2})$ (e.g., $A_{1}$ = dioxin, $A_{2}$ =
furan, both POPs), here we are interested in the probability of the outcome
given joint exposure. We are then interested in reducing this exposure profile
by $\boldsymbol{\delta}=(\delta_{1},\delta_{2})$ which can be static (say some
practical unit or standard deviation) or this can be a function
$\boldsymbol{\delta}(W)$, that is, a change amount based on the history of the
covariates. We use $P$ to denote a generic probability distribution within a
nonparametric statistical model $\mathcal{M}$. The symbol $P_{n}$ represents
the empirical distribution of the observed data, attributing a probability of
$\frac{1}{n}$ to each observation for the sample size $n$. The mapping of
target parameters is defined as $P\rightarrow\Psi(P)$, with $\Psi(P_{0})$ as
the true target parameter under $P_{0}$, and $\hat{\Psi}(P_{n})$ as an
estimator mapping the empirical distribution $P_{n}$ to its estimated value.
The true value of the target parameter is $\psi_{0}=\Psi(P_{0})$, and its
estimate based on $P_{n}$ is represented as $\psi_{n}$. In causal inference,
we often seek to contrast observed data with potential outcomes under
different intervention scenarios (the full causal model). We begin by first
clarifying the distinction between statistical and causal parameters for
stochastic shift interventions for potentially multiple exposures, then we
move on to build our parameter from this identification. We first define our
causal parameter using counterfactuals Rubin [1974, 2005]. For each unit in
our population, we have an infinite family of potential outcomes $Y(\cdot)$
corresponding to every possible value of $\mathbf{A}$. This family,
$Y(\cdot)$, is considered a stochastic process. In practice, for each unit, we
observe only one of these potential outcomes, specifically $Y(\mathbf{A})$,
the outcome corresponding to the observed treatment. However, we can define a
true causal distribution $P^{*}_{0}$ as a distribution over the variables
$Y(\cdot),\mathbf{A},W$. The condition $Y=Y(\mathbf{A})$ defines a mapping
$\mathcal{O}$ from causal distributions $P^{*}$ to observable distributions
$P$. Our causal goal is to define potential outcomes via intervention on
$\mathbf{A}$ to a new value, $\mathbf{A}-\boldsymbol{\delta}$.In this context,
our causal parameter is defined with respect to the causal distribution
$P^{*}$,
$\Psi^{*}(P^{*})=\mathbb{E}_{P^{*}}[Y(\mathbf{A}-\boldsymbol{\delta})]$.

However, because $P^{*}$ includes unobserved potential outcomes, we cannot
estimate $\Psi^{*}(P^{*})$ directly. Instead, we rely on identification
strategies that allow us to express the causal parameter $\Psi^{*}(P^{*})$ in
terms of the observed data distribution, thereby allowing us to make a
practical estimation from the data we have. The assumptions needed for this
are conditional ignorability: $Y(a)\perp A|W$ for all $a$ and positivity:
$0<\frac{g_{\delta}(A|W)}{g(A|W)}<M$ with probability 1. With these
assumptions, we can identify our causal parameter via the following
statistical parameter
$\Psi(P)=E_{P_{\delta}}[Y]=\int_{\mathbf{A}}\int_{W}\bar{Q}(\mathbf{a},w)g_{\delta}(\mathbf{a}|w)\,p(w)\,d\mathbf{a}\,dw$.
This captures the distribution of outcomes $Y$ if we were to simply change the
observed exposure $\mathbf{A}$ by some amount $\boldsymbol{\delta}$, while
keeping all other factors unchanged. Throughout, let $E_{\delta}[\cdot]$
abbreviate $E_{P_{\delta}}[\cdot]$. Note that despite being defined in terms
of $P_{\delta}$, $\Psi(P)$ is indeed a function of $P$ because $P_{\delta}$ is
defined in terms of $P$. Thus, upon invoking these
assumptions,$\Psi^{*}(P^{*})=\Psi(\mathcal{O}(P^{*}))=\Psi(P)$. From now on,
we focus on the target parameter $\Psi$ (defined on the observable
distribution) and our observable estimand $\psi_{0}$ (the value $\Psi$ takes
at our true, unknown distribution $P_{0}$). In practice, identifying causal
effects requires thorough consideration of these assumptions. The conditional
ignorability assumption, in particular, asserts that there are no unmeasured
confounders influencing both the exposure and the outcome, which is a
substantial claim, especially in observational data where mixed exposures are
usually observed. The positivity assumption ensures that there is enough
variability in our data to make meaningful causal inferences. When these
assumptions are reasonable and satisfied, the statistical parameters derived
from the observed data can be insightful proxies for the true causal
parameters. Alternatively, the resulting quantities can be interpreted as pure
statistical parameters of interest, a particular type of association of $A$
with $Y$ controlling for $W$, van der Laan [2006]. In this way, we can
interpret “causal parameters ” as a variable importance measure.

###  2.3 Interaction Target Parameter

Our target parameter above defines the causal parameter for an arbitrary joint
shift of multiple exposures under the stochastic intervention framework.
Realistically, given the usual size and complexity of public health data, no
more than two-way interactions can normally be reliably determined. As such,
we focus on describing the interaction parameter for two exposures that builds
off the joint stochastic shift parameter previously discussed; it is
straightforward to generalize the above to more than two exposures. The symbol
$\delta$ in bold represents the shift of both exposures by two shifted
amounts, $E_{\delta_{1}}[Y]$ represents the change in $A_{1}$ by $\delta_{1}$
leaving $A_{2}$ at the observed levels, and $E_{\delta_{2}}[Y]$ represents a
shift in $A_{2}$ by $\delta_{2}$ leaving $A_{1}$ at the observed levels. As
such, we define our bivariate interaction parameter
$E_{\boldsymbol{\delta}}[Y]-\left(E_{\delta_{1}}[Y]+E_{\delta_{2}}[Y]\right)+E[Y]$.
Which is the expected outcome under the joint shift of $\mathbf{A}$ compared
to the expected additive outcome under each individual shift in $\mathbf{A}$.
The interpretation of this parameter is similar to that of an interaction term
in a linear model.

##  3 Estimation and Interpretation

###  3.1 Efficient Influence Function

In expanding the scope of Diaz and van der Laan’s work on TMLE for these shift
intervention parameters Iván Díaz Muñoz and Mark van der Laan* [2012] to
consider exposures as vectors rather than scalar entities, we simply modify
the notation to possible sets of exposures accordingly. The EIF for
$E_{\delta}[Y]$ is given by
$D_{\delta}(P_{0})(o)=H_{\delta}(\mathbf{a},w)({y-\overline{Q}(\mathbf{a},w)})+\overline{Q}(\mathbf{a}-\boldsymbol{\delta},w)-E_{\delta}[Y]$.
Where $H_{\delta}$ is defined as
$H_{\delta}(\mathbf{a},w)=\frac{g_{0}(\mathbf{a}-\boldsymbol{\delta}\mid
w)}{g_{0}(\mathbf{a}\mid w)}$. Note that vectors $\mathbf{a}$ and
$\boldsymbol{\delta}$ emphasize the element-wise operations and comparisons.
Extension to a vectorial representation allows us to use the original EIF for
multiple exposures. The auxiliary covariate $H_{\delta}(\mathbf{a},w)$ serves
as a ratio of conditional (possibly joint) densities $p(\mathbf{a}|W)$.
Specifically, $g_{0}(\mathbf{a}\mid w)$ denotes the conditional density under
observed values, while $g_{0}(\mathbf{a}-\delta\mid w)$ is the density after a
decrease of $\delta$ to $\mathbf{a}$. This ratio captures the density’s change
upon a $\delta$ shift. In our approach, we use TMLE to remove the plug-in bias
of our $\bar{Q}$ estimators. Given the interaction parameter defined as
$\Psi(P)=E_{\boldsymbol{\delta}}[Y]-\left(E_{\delta_{1}}[Y]+E_{\delta_{2}}[Y]\right)+E[Y]$,
the EIF for this interaction can be succinctly represented in terms of the
EIFs for the individual components
$D(P)=D_{\boldsymbol{\delta}}(P)-\left(D_{\delta_{1}}(P)+D_{\delta_{2}}(P)\right)$,
where $D_{\boldsymbol{\delta}}(P)$ is the EIF of $E_{\boldsymbol{\delta}}[Y]$
and $D_{\delta_{1}}(P)$ and $D_{\delta_{2}}(P)$ are the EIFs of
$E_{\delta_{1}}[Y]$ and $E_{\delta_{2}}[Y]$, respectively.

###  3.2 TMLE

We chose to estimate our parameter using TMLE (rather than estimating
equations, inverse propensity of treatment weights (IPTW) and augmented IPTW )
chiefly due to its favorable finite sample properties due to being a plug-in
estimator. Briefly, to construct the TMLE estimator, we: 1. Use data-adaptive
regression for the initial estimation of $g_{0}(\mathbf{a},W)$ and
$\overline{Q}_{0}(\mathbf{a},W)$, 2. Evaluate $H(\mathbf{a}_{i},w_{i})$ for
each observation, 3. Using the derived estimates, we implement a logistic
regression to refine predictions under the shift
$\text{logit}\overline{Q}_{\epsilon,n}(\mathbf{a}-\delta,w)=\text{logit}\overline{Q}_{n}(\mathbf{a}-\delta,w)+\epsilon
H_{n}(\mathbf{a}-\delta,w)$ then 4. From this model, once $\epsilon$ is
determined, we update the predicted counterfactuals, using these updated
counterfactuals, we then compute the plugin. This is a standard TMLE procedure
for stochastic interventions; we simply shift multiple exposures
simultaneously in the initial estimation and update. We estimate a standard
error by taking the empirical variance of the estimated EIF and use this
construct Wald-type confidence intervals and p-values. We use TMLE to
ascertain unbiased estimates for each part of our interaction target
parameter. That is, in $E_{\boldsymbol{\delta}}[Y]$ we update our
counterfactual outcome under a joint shift by the joint density ratio as
$H(\boldsymbol{a},w)$. Similarly, $E_{\delta_{1}}[Y]$ and $E_{\delta_{2}}[Y]$
are simply updated based on $H(a,w)$ for univariate shifts. In general, we use
TMLE to target each of the three parts of our parameter. Given that our
interaction is a continuous function (linear combination) of these
asymptotically efficient estimates, we use the delta method for variance
derivation.

###  3.3 Joint Density Estimation

In the context of extending the clever covariate to account for joint shifts,
a challenge emerges when confronting the estimation of the joint density,
especially in scenarios with high-dimensional $W$. Direct estimation of joint
conditional densities using conventional techniques scales poorly to large
data Chen et al. [2018]. To avoid the complications inherent in the direct
estimate of joint density, we take an alternative approach. Instead of seeking
an explicit joint density estimator, the methodology takes advantage of the
combination of conditional and marginal densities. In the bivariate scenario
involving exposures $A1$ and $A2$, the joint density can be calculated more
conveniently as $g_{0}(\mathbf{A}|W)=g_{0}(A1|W)\times g_{0}(A2|A1,W)$. Within
this framework, two ensemble-based machine learning density estimators are
deployed. The first is used to approximate the density of $A1$ conditioned on
$W$, while the second estimates the conditional density of $A2$ given both
$A1$ and $W$. The product of these constructed densities provides an effective
estimate of the joint density, which avoids some of the inherent challenges of
direct joint density estimation.

##  4 Exposure Set Identification

The above presentation of the estimation of the marginal, joint and
interaction shift parameters assumes that the subsets of the mixture that are
important with respect to the outcome, as well as the pairs that have evidence
of statistical interaction, are known a priori. However, in finite samples
with a large number of mixture variables, $\mathbf{A}$, including all
components of the mixture as potentially impactful exposures and all possible
interactions will result in unacceptable sample variability (and potential
nonnormality). Therefore, we propose using a data-adaptive target parameter
approach Hubbard et al. [2016] to use sample splitting to “discover” the
relevant components and interactions among the mixture variables and estimate
the resulting data-adaptively defined parameters.

We propose a K-fold cross-validation in which the sample data are split into
K-folds of roughly equal size. For example, consider $K$ = 10, and define
$k=1$. The $K-1/K$ portion of the sample is used to define the parameter of
interest (here 90%) and the $1/K$ (10%) is used to estimate the parameters.
Let $T_{n,k}\subset\\{1,\ldots,n\\}$ denote the observations indices in the
k-th training sample and $V_{n,k}$ denote the observations indices in the k-th
validation sample for $k=1,\ldots,K$. The following paragraphs describe how
through cross-validation the $T_{n,k}$ (training) portion of the data is used
to identify important variable sets in the mixture and $V_{n,k}$ (validation)
is used to estimate our stochastic shift intervention parameters.

###  4.1 Identifying Interactions in the Mixture Using Training Data

The discovery of variables and interactions of mixture variables that are
important for the prediction of the outcome is done using an ensemble method:
the discrete Super Learner (SL) van der Laan Mark et al. [2007]. This
algorithm employs a cross-validation mechanism, optimizing over a library of
candidate algorithms to select the estimator with the best empirical fit. We
use a discrete SL to select the best fitting estimator from a library of
different sieve-based estimators that use basis splines and their tensor
products, allowing the model to capture intricate nonlinear relationships as
linear combinations of lower-dimensional basis functions while remaining
interpretable. These include packages like earth Milborrow. Derived from
mda:mars by T. Hastie and R. Tibshirani. [2011] , polySpline Ripley and
Venables [2021], and hal9001 Coyle et al. [2022]. In our current
implementation, our candidate algorithms used in the library are all based on
earth. Because earth has many hyperparameters such as the number of
interactions, how knot points are searched, etc. we include earth models with
various parameters and select the best one based on cross-validation. This
still falls within the framework of Super Learning. We employ this estimator
for $E[Y|\boldsymbol{A},W]$ used in the data-adaptive discovery procedure
using the training data. For example, for $k=1$ in a 10 fold CV, this includes
all folds not equal to 1. In the context of high-dimensional mixture analysis,
traditional ANOVA techniques are likely inadequate due to their limitations in
handling nonlinear relationships. To address this, we introduce a non-
parametric extension, Non-parametric Analysis of Variance (NP-ANOVA),
following the selection of the best-fitting estimator through the Super
Learner algorithm. The NP-ANOVA approach is designed to decompose the variance
of the response variable by attributing it to contributions from individual
basis functions utilized in the chosen model. This method is versatile and is
applicable across various model structures such as MARS, HAL, among others.
The central idea is to assess the relative importance of each basis function
in explaining the variance of the response variable. The procedure can be
summarized as follows. 1. From the model selected by the Super Learner, we
extract the basis functions employed by the model, 2. Utilize these basis
functions to build a GLM, allowing for a structured analysis of the
contributions of bases, 3. Perform a Type 3 ANOVA on the GLM. This step
involves systematically evaluating the change in explained variance when each
basis function is omitted, thus assessing their individual contributions, 4.
For each basis function, calculate an F-statistic, quantifying its impact on
the explanatory power of the model, 5. Group and sum the F-statistics
according to the variable sets corresponding to each basis, such as individual
variables (e.g., $A_{1}$) or interactions (e.g., $A_{3}-A_{4}$), 6. Employ
quantile-based thresholding on aggregated F-statistics to discern the most
influential variable sets. For example, selecting the top 25% of variable sets
or including all sets by setting the threshold at 0.

###  4.2 Estimation of Shift Interventions on Discovered Variable Sets

Treating these variable sets as fixed, we then estimate the resulting data-
adaptively defined parameters using TMLE for marginal and joint stochastic
shift interventions. Using the training data, the nuisance functions $g_{n}$
and $Q_{n}$ are trained. The TMLE update step is performed using these
training estimates to estimate $\epsilon$. With the validation data, we obtain
estimates for $g_{n}$ and $Q_{n}$ and apply $\epsilon$ to debias our initial
plug-in for our estimated stochastic shift target parameter using the
validation data. In our running example, the validation data is data that
correspond to $k=1$. This procedure describes CV-TMLE Hubbard et al. [2018,
2016] which is the infrastructure we use to estimate our target parameter.
Readers interested in the details of this procedure are directed to these
resources.

With the estimates for our parameter complete for the validation for $k=1$, we
would then move on to $k=2$ and repeat the sections 2-4. Upon completion of
the cross-estimation procedure, a pooled TMLE update provides a summary
measure of the target parameter across k-folds. Specifically, we stack the
predictions for each nuisance parameter from the estimation samples and run a
pooled TMLE update on the cumulative initial estimates. The resulting average
is then used for parameter estimation. The pooled estimate across the folds
can be succinctly represented as
$\theta=\frac{1}{K}\sum_{k}\Psi_{F_{T_{n,k}}}\left(V_{n,k}\right)(P_{0})$.
Where $\Psi_{F_{T_{n,k}}}$ identifies optimal exposure sets and subsequently
formulates a plug-in estimator using training data $T_{n,k}$. $V_{n,k}$ refers
to our data’s estimation sample for deriving our estimates.

It is possible that the same interactions are not found in every fold. Here,
our approach deviates from the standard data-adaptive target parameter method.
In defining parameters normally, one would define, for example, the
interaction with maximum impact defined by some measure (like our F-statistic
or another measure) or perhaps create ranks. However, this may result in many
different interactions being included in the top rank, which makes
interpretation difficult. In this setting, the analyst is then required to
interpret only the fold-specific results, which are statistically valid but
have low power, and perhaps report the proportion of folds the interaction was
found. In the proposed framework, we gain power by pooling in folds with
consistent estimates, but the pooled result is difficult to interpret as well
because these estimates do not incorporate information from folds with null
results. As such, in refining our approach to pooled TMLE estimation, our goal
is to integrate both the non-null estimates from the estimation samples and
the null values. To address this, we propose a weighted average approach that
combines our original pooled estimate with the implied null value, which we
define as 0.

Specifically, let $n_{0}$ denote the number of folds where the null is
implied, and let $n_{1}$ represent the number of folds with non-null
estimates. The variance of the "null" estimate, $var(\Psi_{\text{null}})$, is
estimated using the variance of the original pooled estimates,
$var(\hat{\Psi}_{\text{pooled}})$, as follows
$var(\Psi_{\text{null}})=\left(var(\hat{\Psi}_{\text{pooled}})\right)\cdot\frac{n_{1}}{n_{0}}$.
We then calculate the weighted average of the original pooled estimate and the
null value
$\hat{\Psi}_{\text{pooled}}^{*}=\frac{w_{1}\cdot\hat{\Psi}_{\text{pooled}}+w_{0}\cdot\Psi_{\text{null}}}{w_{0}+w_{1}}$,
where $w_{1}=\frac{1}{var(\hat{\Psi}_{\text{pooled}})}$ and
$w_{0}=\frac{1}{var(\Psi_{\text{null}})}$. The variance of this new pooled
estimate is given by
$var(\hat{\Psi}_{\text{pooled}}^{*})=\frac{1}{w_{1}+w_{0}}$. This approach
ensures that our pooled TMLE estimate reflects both the explicit estimates
from folds where the variable set was selected and the implicit null values
from the remaining folds. Of course, in settings where there is strong
"signal" and interactions are consistently found across the folds this inverse
variance approach is not needed and only the standard pooled TMLE is used.
Fold-specific TMLE updates are also provided, allowing researchers to gauge
consistency by comparing estimates across individual folds. Each fold-specific
update applies a similar logic as the pooled update but remains confined to
the respective fold, utilizing its specific $Q_{k}$ and $h_{k}$ parameters.
Consequently, the EIF for each fold is calculated using only validation data
within the fold, providing distinct variance estimates.

###  4.3 Pooled Estimates under Data-Adaptive Delta

Stochastic interventions, particularly joint stochastic shifts that involve
products of conditional densities, may be vulnerable to positivity violations.
Such violations can introduce bias and inflate variance in exposure effect
estimation after shift, even with sophisticated estimators like TMLE.
Addressing this necessitates a restrained $\mathbf{\delta}$, diminishing
positivity violations, and thereby reducing bias and variance. We propose
treating $\mathbf{\delta}$ as data-adaptive. It is determined in the training
sample, adhering to positivity criteria, and subsequently applied to the
validation sample. Define $H(a_{\delta},w)_{i}$ as the density ratio for
observation $i$ after a shift $\delta$. To maintain ratios below a threshold
$\lambda$, we iteratively diminish $\delta$ by $\epsilon$ until: $\forall
i,\;H(a_{\delta},w)=\frac{g_{n_{-k}}(a_{n_{-k}}-\delta\mid
w)}{g_{n_{-k}}(a_{n_{-k}}\mid w)}\leq\lambda$. Here, $\lambda$ is preset
(default is 50 in our implementation). If any predicted conditional
probabilities exceed this threshold under observed exposure, $\delta$
decreases by $\epsilon$ (defaulted to 10% in our package). This adjustment
employs training data with estimators trained therein. Should $\delta$ remain
constant across folds, the pooling involves averaging the estimates of the
target parameters. For a data-adaptive $\delta$, where $\delta$ may be
different in different folds, we simply take the mean estimate for a pooled
$\delta$ across the folds that is paired with the pooled TMLE estimate.

##  5 Simulations

In this section, we demonstrate using simulations that our approach identifies
the correct exposures and interactions used to generate the outcome and
correctly estimates the counterfactual mean difference in outcome under joint
shift intervention, marginal shifts, and subsequently the comparison, our
interaction parameter. We construct a data-generating process (DGP) where $Y$
is generated from a linear combination of different effects including marginal
impacts and interactions. This DGP was constructed to represent the complexity
of a mixed exposure in which mixture variables are correlated but only some
have an impact on the outcome.

This DGP has the following characteristics, $O=W,A,Y$. Let
$W=(W_{1},W_{2},W_{3})$ denote a random vector of covariates, where $W_{1}$
and $W_{2}$ are bivariate normal with mean vector $\boldsymbol{\mu}=(6,7)$ and
covariance matrix $\boldsymbol{\Sigma}_{W}=\begin{pmatrix}1&0.4\\\
0.4&1\end{pmatrix}$ and $W_{3}$ is a binary variable with Bernoulli
distribution with probability $0.5$. The mean values for the components of the
distribution of mixture $A$ are calculated as $\mu_{1}==e^{W_{1}/2}$,
$\mu_{2}=W_{2}/2$, and $\mu_{3}=5$. That is, exposure 3 is not associated with
the covariates. A exposure vector $A=(A_{1},A_{2},A_{3})$ is generated from a
trivariate normal distribution with a mean vector
$\boldsymbol{\mu}=(\mu_{1},\mu_{2},\mu_{3})$ and a covariance matrix:
$\boldsymbol{\Sigma}_{A}=\begin{pmatrix}1&0.5&0.8\\\ 0.5&1&0.7\\\
0.8&0.7&1\end{pmatrix}$. A fourth mixture variable was created via
$A_{4}\sim\mathcal{N}(4,2)$. The outcome $Y$ is then generated by:
$Y=1.3A_{4}+0.4A_{4}A_{1}+0.1W_{1}+0.3W_{2}+\epsilon+(A_{3}^{2}\mathbf{1}_{\\{W_{3}=1\\}}+A_{3}\mathbf{1}_{\\{W_{3}=0\\}})$.
Where the third component of $A$ depends on $W_{3}$, $A_{3}$ is squared if
$W_{3}$ is equal to 1, and is left unchanged if it is equal to 0. In
generating the outcome, this variable represents effect modification or the
differential impact of $A_{3}$ depending on $W_{3}$. Overall, built into this
outcome-generating process there is a marginal impact of $A_{4}$, a
multiplicative interaction between $A_{1}$ and $A_{4}$ and effect modification
of $A_{3}$ based on $W_{3}$. There is no effect of $A_{2}$ on the outcome.
This DGP was generated to reflect, although more simply, a real-world exposure
situation that includes both interaction, effect modification, and no effects.

To calculate the ground-truth counterfactual mean under a shift intervention,
we simply generate a very large data set from this DGP and apply various
shifts to the exposures. We applied a shift of 1 unit increase to each
exposure. In the first simulation, we deterministically set the variable sets
to evaluate estimates of the counterfactual outcomes under these shifts. This
was done with the var_sets parameter in SuperNOVA which, if is not null,
bypasses the data-adaptive determination of the variable sets in the basis
functions. In the second simulation we include the data-adaptive discovery of
the variable sets to provide estimates for evaluating the variable
relationship identification through the basis function procedure. In our case,
given this DGP, we are interested in assessing the counterfactual mean
differences under shifts for $A_{4}$ individually, $A_{1}$ individually,
$A_{4}$ and $A_{1}$ jointly and comparing this joint shift with the individual
shifts.

###  5.1 Evaluating Performance

We assessed the asymptotic convergence to the true exposure relationships used
in the DGP, as well as the convergence to the true counterfactual differences
for these exposures, in each simulation. To do so, we followed the following
steps 1. From this sample, we generated a random sample of size $n$ from the
DGP, 2. At each iteration, we used the training sample to define the
exposure(s) and create the necessary estimators for the target parameter
dependent on the variable sets. We then used the validation sample to obtain
the updated causal parameter estimate using TMLE. We repeated this process for
all the folds, 3. At each iteration, we output the stochastic shift estimates
given the pooled TMLE.

To evaluate the performance of our approach, we calculated several metrics for
each iteration, including bias, variance, MSE, confidence interval (CI)
coverage, and the proportion of instances in which the true variable
relationships were identified. To visually inspect if the rate of convergence
was at least as fast as $\sqrt{n}$, we show projections of a $\sqrt{n}$
consistent estimator starting from the initial bias. We then calculated the
variance for each estimate and used it to calculate the mean squared error
(MSE) as $\text{MSE}=\text{bias}^{2}+\text{variance}$. To account for
different sample sizes, we multiply the MSE estimates by $n$. For each
estimate, we also calculated the CI coverage of the true stochastic shift
parameter given the exposure determined adaptively from the data. We
calculated these performance metrics at each iteration, performing 50
iterations for each sample size $n=$ (250, 500, 1000, 1500, 2000, 2500, 3000,
5000). We used SuperNOVA with 5-fold cross-validation (to speed up
calculations in the simulations) and default learner stacks (see below) for
each nuisance parameter and data-adaptive parameter. Additionally, the
quantile threshold to filter the basis functions based on F-statistic was set
to 0 to include all basis functions used in the final best-fitting model. For
each sample size, we calculate the proportion of iterations in which the
correct main effects and interaction effects were discovered.

To use SuperNOVA, we need estimators for $\bar{Q}=E(Y|A,W)$ and $g_{n}=p(A|W)$
(conditional density). SuperNOVA provides default algorithms to be used in a
Super Learner van der Laan Mark et al. [2007] that are both fast and flexible.
For our data-adaptive procedure, we include learners from the packages earth
Milborrow. Derived from mda:mars by T. Hastie and R. Tibshirani. [2011],
polspline Ripley and Venables [2021], and hal9001 Coyle et al. [2022]. The
results from each of these packages can be formed into a model matrix, on
which we can fit an ANOVA to obtain the resulting linear model of basis
functions. To estimate $\bar{Q}$, we include estimators for the Super Learner
from glm, elastic net Tay et al. [2023], random forest Wright and Ziegler
[2017], and xgboost Chen and Guestrin [2016]. For the semiparametric density
estimator $g_{n}$, we create estimators based on homoscedastic errors (HOSE)
and heteroscedastic errors (HESE) from the same estimators used in $\bar{Q}$.

###  5.2 Results

SuperNOVA consistently discovers interactions. The interaction effect was
identified in 100% of the folds in all sample sizes. This indicates that, at
least in this simple simulation, the default learners are flexible enough to
capture the interaction effect. Target parameters estimated by SuperNOVA
converge to truth at a $1/\sqrt{n}$ rate, Figure 1 A shows the absolute bias,
CI coverage, mean squared error, and standard deviation as the sample size
increases to 5000. For bias, MSE and standard deviation, there is a converge
to zero at sample size 5000. For coverage, the average coverage for each
target parameter was: individual shift: 98%; joint shift: 96%; interaction:
96%.

Figure 1: A. Bias, CI, Standard Deviation and MSE for Each Parameter Across
Sample Sizes: These metrics represent results from CV-TMLE, that is bias is
w.r.t. the data adaptive target parameter where we use CV to identify the
exposure and then estimate. B. Bias Standardized by Standard Error Compared to
Ground-Truth Outcome Under Shift Interventions

To establish SuperNOVA’s estimator validity, it is essential to demonstrate
that its sampling distribution approximates a normal distribution centered at
0, with its variance decreasing as sample size grows. This is evaluated
through the sampling distribution of each estimate. Figure 1 B illustrates the
probability density of the standardized bias relative to true estimates in 50
iterations for varying sample sizes, indicating convergence to a mean 0 normal
distribution as sample size increases. Specifically, Figure 1 B. Individual
Effect displays the convergence of the standardized bias for a single variable
shift towards a normal distribution with mean 0 and standard deviation 1.
Likewise, Figure 1 B. Joint Effect and Figure 1 B. Interaction depict the
sampling distributions for dual variable shifts and interaction parameter
estimates, respectively. In all scenarios, the distributions align with a
standard normal distribution as the sample size increases, ensuring the
reliability of our inferential approach.

##  6 Applications

###  6.1 Analysis of NIEHS Synthetic Mixtures Data

The NIEHS synthetic mixtures data is a commonly used data set to evaluate the
performance of statistical methods for mixtures. These synthetic data can be
considered the results of a prospective cohort study, where the outcome cannot
cause exposures, and correlations between exposure variables can be thought of
as caused by common sources or modes of exposure. The variable $W$ can be
assumed to be a potential confounder and not a collider. The data set has 7
exposures ($A_{1}-A_{7}$) with a complex dependency structure based on
endocrine disruption. Two groups of exposure ($A_{1},A_{2},A_{3}$ and
$A_{5},A_{6}$) lead to high correlations within each group.
$A_{1},A_{2},A_{7}$ positively contribute to the outcome, $A_{4},A_{5}$ have
negative contributions, while $A_{3}$ and $A_{6}$ have no impact on the
outcome. Rejecting $A_{3}$ and $A_{6}$ is difficult due to their correlations
with the members of the cluster group. This correlation and effects structure
is biologically plausible, as different congeners of a group of compounds may
be highly correlated but have different biological effects. Exposures have
various agonistic and antagonistic interactions, a breakdown of the variables
sets and their relationships are: $A_{1},A_{2}$ both increase Y by
concetration addition, $A_{1}-A_{4},A_{2}-A_{4},A_{1}-A_{4},A_{2}-A_{4}$
competetive antagonism, $A_{1}-A_{7},A_{2}-A_{7}$ supra-additive or
synergistic, $A_{4},A_{5}$ both addivively decrease Y,
$A_{4}-A_{7},A_{5}-A_{7}$ unusual antagonism. Synthetic data and the key for
dataset 1 are available on GitHub. This resource shows the interactions and
marginal dose relationships built into the data. Given these toxicological
interactions, we expect these sets of variables to be determined in SuperNOVA.
For example, we might expect a positive counterfactual result for A1, A2, A7
and negative results for A4 and A5. Likewise, in the case for antagonistic
relationships such as A1-A5 (third row in Table 1), we would expect a joint
shift to come closer to the null, since A5 antagonizes the positive effects of
A1. For A1 and A2, we would expect the joint shift to be close to the sum of
individual shifts (not much interaction), but for A1 and A7 we expect to
estimate more than an additive effect (some interaction). The NIEHS data set
has 500 observations and 9 variables. $W$ is a binary confounder. Of course,
in this data there is no ground truth, like in the above simulations, but we
can gauge SuperNOVA’s performance by determining if the correct variable sets
are used in the interactions and if the correct variables are rejected.

We apply SuperNOVA to this NIEHS synthetic data using a 20-fold CV. For
computational simplicity, we considered a small learner library consisting of
a main terms linear regression and random forest models for estimation $g$ and
$Q$. For the discovery of the exposure variable sets, we built 13 earth models
with various hyperparameters. We parallelize over the cross-validation to test
computational run-time on a newer personal machine an analyst might be using.
We use a delta of 1 for all the exopsures in the mixture. We present the
pooled TMLE results and the pooled TMLE results with inverse variance for null
folds. Due to space limitations we do not show fold specific results but these
are provided to the user.

###  6.2 NIEHS Synthetic Data Results

SuperNOVA accurately rejects exposures $A_{3},A_{6}$, finding the correct
marginal effects for $A_{1},A_{2},A_{4},A_{5},A_{7}$ consistently in all
folds. Specifically, A1: 100%, A2: 85%, A4: 95%, A5: 100%, A7 100%. The
interactions were harder to identify consistently, and the strongest
interaction between A5 and A7 was detected in 70% of the folds. Specifically,
A1-A2: 10%, A1-A5: 45%, A1-A7: 20%, A2-A7: 50%, A4-A5: 15%, A4-A7: 5%, and
A5-A7: 70%. Our analysis of the NIEHS Mixtures Data reveals the marginal
effects of various exposures on endocrine disruption outcomes, measured as the
change in outcome per unit increase in relative exposure. Each exposure’s
effect, along with its confidence interval, p-value, and the proportion of
folds in which the effect was observed, are reported as follows. A1
demonstrated a significant increase in the outcome by 5.69 units, with a 95%
CI of 4.81 to 6.57, identified consistently across all folds (proportion =
1.0), and a p-value < 0.01. A5 Showed a significant decrease in the outcome by
-3.41 units, with a 95% CI of -4.28 to -2.54, consistently observed across all
folds (proportion = 1.0), and a p-value < 0.01. A7 indicated a significant
increase in the outcome by 3.40 units, with a 95% CI of 2.47 to 4.33, found in
all folds (proportion = 1.0), and a p-value < 0.01. A4 had a non-significant
effect of -0.36 units (inverse variance pooled), with a 95% CI of -1.69 to
0.22 and observed in 95% of folds, with p-value of 0.29. A2 exhibited a
significant increase in the outcome by 1.82 units (inverse variance pooled),
with 95% CIs of 0.93 to 2.72, found in 85% of folds, and both with p-values <
0.01. These findings align with the direction and relative strengths of
association built into this synethetic data DGP.

We found three significant interaction effects defined as deviations when
considering joint shifts in exposures compared to the sum of individual
shifts. These interactions were for exposures A1-A5, A1-A7, and A2-A7, which
aligns with the interactions built into this DGP. For the interaction between
A1 and A5, we observe that a unit shift in A1 individually leads to an
increase of 10.40 (9.08-11.73) in the outcome. On the contrary, a unit change
in A5 individually shows a decrease of -2.86 (-4.10 - -1.62), suggesting an
antagonistic relationship between these exposures. When both A1 and A5 are
simultaneously shifted by one unit, the expected outcome under joint shift is
0.49 (-0.66 - 1.63), thus generally nullifying each exposure effect, making
the interaction effect -7.26 (-8.59 - -5.92), p < 0.001. This effect was found
in 45% of folds. The interaction between A1 and A7 shows that individual
changes in A1 and A7 lead to increases of 14.74 (13.14 - 16.34) and 2.42 (0.62
- 4.22), respectively. However, the joint shift in X1 and X7 only increases
endrocine function by 5.31 (3.65 - 6.98), resulting in an interaction effect
of -11.76 (-13.4 - -10.12). This effect indicates that the joint shift is
significantly less than expected under additive assumptions. This interaction
was found in 20 % folds. The shift in A2 led to a 5.56 (4.19 - 6.93) and A7 a
2.83 (1.65 - 4.02), the joint shift a 4.89 (3.75 - 6.04), which led to an
interaction effect of -3.86 (-5.24 - -2.48), p < 0.001, proportion of folds =
50%. Individual shifts in A5 and A7 result in changes of -2.74 (-3.70 - -1.78)
and 2.61 (1.60 - 3.62), respectively. However, the interaction effect of their
joint shift is -0.02 (-1.03 - 0.98), which leads to an interaction effect of
-0.05 (-1.02 - 0.93), found in 70% of the folds. This suggests that these two
exposures have an additive antagonism, which is built into the DGP. Due to
space limitations, we do not provide fold-specific results and more details,
but these are provided on github package page. Although our detection of
interactions was not found in every fold (due potentially to a low sample
size), we identified the majority of interactions built into these data.

###  6.3 Comparison to Quantile G-Computation

Quantile g-computation is a favored approach in environmental epidemiology’s
mixture analyses. Estimates the joint effect of increasing all exposures by
one quantile, based on linear assumptions, as:
$Y_{i}=\beta_{0}+\sum_{j=1}^{d}\beta_{j}X_{ji}^{q}+\beta Z_{i}+\epsilon_{i}$.
Here, $X^{q}$ represents the quantized components of the mixture. The
methodology consists of transforming the components of the mixture into
quantiles and then adding the coefficients of the linear model to produce a
summary measure of the mixture, $\Psi$. However, several limitations emerge:
1. Quantiles may not accurately reflect the exposure-response relationship,
which can be nonmonotonic in substances like endocrine disruptors, 2. For a
comprehensible mixture estimate $\Psi$, it presumes an additive relationship
of quantiles, 3. Given the intricate nature of mixtures, interactions are
vital. Including interactions in quantile g-computation complicates the
interpretation, as the exposure’s proportional contribution becomes variable
depending on other factors’ levels, 3. This method lacks flexibility in
covariate control. Our application to the NIEHS dataset with 4 quantiles,
excluding interactions, yielded a scaled positive effect size of 6.28
(encompassing $A_{1},A_{2},A_{3},A_{7}$) and a negative effect size of -3.68
(involving $A_{4},A_{5},A_{6}$). Discrepancies arise when compared to the
NIEHS baseline, incorrectly including $A_{3},A_{6}$. Incorporating
interactions compounds the challenge. Our approach with all exposure
interactions showed a joint impact of main effects of 0.02 (-3.16 - 3.20), p =
0.99 and interaction effects of 0.59 (-0.71 - 1.90), p = 0.37. Considering the
parameter count and a sample size of $n=500$, it is not surprising that these
estimates are not significant. However, even attempts to rigorously analyze
2-way or 3-way interactions introduce an array of complications, underscoring
mixtures as an inherently data-adaptive issue and the insufficiencies of
conventional methods on even straightforward synthetic datasets.

###  6.4 Methods Applied to NHANES Dataset

The National Health and Nutrition Examination Survey (NHANES) 2001-2002 cycle
serves as the foundation of our analysis. This data source, known for its
credibility in the public health domain, boasts interviews with 11,039
individuals. Of this subset, 4,260 provided blood samples and willingly
consented to DNA analysis. The data set used aligns with that of Mitro et al.
Mitro et al. [2016] focusing on the correlation between persistent exposure to
organic pollutants (POPs), specifically those binding to the aryl hydrocarbon
receptor (AhR) and extended leukocyte telomere length (LTL). However, this
subset was further refined to ensure complete exposure data, yielding 1007
participants for our analysis, compared to 1003 in the initial study by Mitro
et al. In alignment with protocols detailed by Mitro et al. Mitro et al.
[2016], exposure was quantified, focusing on 18 congeners. These include 8
non-dioxin-like PBCs, 2 non-ortho PCBs, 1 mono-ortho PCB, 4 Dioxins, and 4
Furans. All congeners underwent lipid serum adjustments using an enzymatic
summation methodology. The telomere length, was analyzed using the
quantitative polymerase chain reaction (qPCR) methodology Mitro et al. [2016].
This technique measures the T/S ratio by comparing the length of the telomere
to a standardized reference DNA. To enhance accuracy, samples underwent triple
assays in duplicate, generating an average from six data points. The CDC
conducted a blinded quality control assessment. Our modeling accounted for
several covariates, including demographic factors such as age, sex, race /
ethnicity and level of education, as well as health indicators such as BMI,
serum cotinine, and blood cell distribution and count. Categorizations for
race/ethnicity, education, and BMI are consistent with previous studies.
Furthermore, this comprehensive dataset is integrated into the SuperNOVA
package for replicable analysis.

We configured our methodology using default learners and employed a 20-fold CV
without a quantile threshold for the F-statistic of the basis function. The
built-in algorithms of the Super Learner were used to identify variable
associations in the discovery procedure, using linear terms and two-way
interactions of earth models and various hyperparameters (13 models in total).
Given the data range, a $\delta$ value of 2 was fixed for all exposures,
denoting a focus on counterfactual changes in telomere length for an increase
in 2 ng / g or pg / g across exposures that predict telomere length depending
on the scale of exposure.

###  6.5 Comparing NHANES Furan Results to Previous Work

In our analysis of the 18 POPs, furan 2,3,4,7,8-pncdf emerged consistently
across all 20 folds. None of the other chemicals were featured in any fold.
The pooled TMLE estimated a 0.02 (0.03 - 0.001) decrease in telomere length
after a 2 pg / g increase in furan 2,3,4,7,8-pncd which was significant at p =
0.04. Results in individual folds showed similar effects. The Mitro et al.
[2016] summarize the results of the workshop. They state that clustering
methods identified high, medium, and low POP exposure groups with longer log-
LTL observed in the high exposure group. Principal component analysis (PCA)
and exploratory factor analysis (EFA) revealed positive associations between
overall POP exposure and specific POPs with log-LTL. Penalized regression
methods identified three congeners (PCB 126, PCB 118, and furan
2,3,4,7,8-pncdf) as potentially toxic agents. WQS identified six POPs (furans
1,2,3,4,6,7,8-hxcdf, 2,3,4,7,8-pncdf, and 1,2,3,6,7,8-hxcdf, and PCBs 99, 126,
169) as potentially toxic agents with a positive overall effect of the POP
mixture. BKMR found a positive linear association with furan 2,3,4,7,8-pncdf,
suggestive evidence of linear associations with PCBs 126 and 169, a positive
overall effect of the mixture, but no interactions among congeners. These
results (in the supervised methods) controlled for the same covariates.
Interestingly, although we corroborate the finding of furan 2,3,4,7,8-pncdf,
specifically in the BKMR method (the most flexible method used), our results
show a negative association. Additionally, we do not find associations of any
other POP used in any fold, meaning that in the best-fitting spline model, no
basis functions for any of the other POPs were used. The inverse association
we find compared to these other methods could be a result of our $Q$ model
being a Super Learner, which can model relationships perhaps missed in these
other methods. Through this NHANES example, we underscore the ability of
SuperNOVA in discerning pertinent variable relationships in complex mixture
data, subsequently producing estimates. The consolidated findings allow for a
clear understanding of the potential changes in outcomes in different exposure
scenarios, improving the interpretability of the results. Interestingly,
despite the heterogeneity of the findings of other methods, we identify only
one chemical possibly related to a reduced telomere length, furan
2,3,4,7,8-pncdf, which warrants further investigation.

##  7 Discussion

In our study, we developed a semiparametric statistical model to estimate
interaction effects in mixed exposures, focusing on identifying and
quantifying exposure sets that influence outcomes. Using data-adaptive target
parameters and using targeted maximum likelihood estimation (TMLE) within a
cross-validation framework, we achieved asymptotically unbiased estimates of
interaction effects. Our approach offers several advantages, including
robustness to the number of exposures and the complexity of the data-
generating process, and it provides interpretable results that improve
understanding of drug synergies and environmental health impacts. Unlike
traditional interaction analysis methods, our model minimizes assumptions,
yielding more reliable and interpretable estimates of interaction effects and
focuses on causal quantities. Our methodology faces two primary limitations.
First, the computational intensity of the density estimation step, though
mitigated by parallel processing, remains a challenge. Alternatives such as
recasting density estimation as a classification problem could offer
efficiency improvements, but would require iterative refitting of the
classifier. Second, the identification of exposure sets varies between cross-
validation folds, introducing variability that demands careful and consistent
reporting to ensure reproducibility and reliability of results. Future
enhancements, including the incorporation of algorithms like glmnet and the
highly adaptive lasso (HAL), are expected to stabilize identification
processes. It is crucial for findings to be reported only if they demonstrate
consistency across multiple cross-validation folds, with transparent
communication about the variability of these findings.

In future research, we aim to address the limitations and expand the
capabilities of our approach. Key directions include improving density
estimation. We plan to explore the transformation of density estimation into a
classification problem using Bayesian probability. This shift promises to
improve computational efficiency and incorporate a wider range of machine
learning algorithms. Refining pooling strategies and testing is also required;
we will develop and evaluate new pooling strategies, such as kernel estimators
to estimate the sampling distribution of our estimates with null effects and
sampling from null distributions to refine pooled estimates incorporating the
null, hopefully improving the robustness of our inverse variance strategy. We
can also focus beyond two-way interactions to include higher-order
complexities. We also intend to investigate alternative TMLE strategies that
may offer better performance, especially in smaller samples. This includes
direct computation of influence functions for interaction parameters and novel
estimation techniques for clever covariates. To facilitate adoption, we have
developed the SuperNOVA R package, available on GitHub, offering researchers
software to apply our proposed methodology for the analysis of mixtures, this
package has been peer reviewed by the Journal of Open-Source Software McCoy et
al. [2023].

## References

  * Bobb et al. [2014] Jennifer F. Bobb, Linda Valeri, Birgit Claus Henn, David C. Christiani, Robert O. Wright, Maitreyi Mazumdar, John J. Godleski, and Brent A. Coull.  Bayesian kernel machine regression for estimating the health effects of multi-pollutant mixtures.  _Biostatistics_ , 16(3):493–508, 2014.  ISSN 14684357.  doi:10.1093/biostatistics/kxu058.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi:10.1145/2939672.2939785.  URL http://doi.acm.org/10.1145/2939672.2939785.
  * Chen et al. [2018] Wei-Chia Chen, Ammar Tareen, and Justin B Kinney.  Density estimation on small data sets.  _Phys. Rev. Lett._ , 121(16):160605, October 2018.
  * Coyle et al. [2022] Jeremy R Coyle, Nima S Hejazi, Rachael V Phillips, Lars WP van der Laan, and Mark J van der Laan.  _hal9001: The scalable highly adaptive lasso_ , 2022.  URL https://github.com/tlverse/hal9001.  R package version 0.4.3.
  * García-Villarino et al. [2022] Miguel García-Villarino, Antonio J. Signes-Pastor, Margaret R. Karagas, Isolina Riaño-Galán, Cristina Rodríguez-Dehli, Joan O. Grimalt, Eva Junqué, Ana Fernández-Somoano, and Adonina Tardón.  Exposure to metal mixture and growth indicators at 4–5 years. A study in the INMA-Asturias cohort.  _Environmental Research_ , 204, 2022.  ISSN 10960953.  doi:10.1016/j.envres.2021.112375.
  * Hubbard et al. [2018] Alan Hubbard, Chris Kennedy, and Mark Laan.  _Data-Adaptive Target Parameters_ , pages 125–142.  01 2018.  ISBN 978-3-319-65303-7.  doi:10.1007/978-3-319-65304-4_9.
  * Hubbard et al. [2016] Alan E. Hubbard, Sara Kherad-Pajouh, and Mark J. Van Der Laan.  Statistical Inference for Data Adaptive Target Parameters.  _International Journal of Biostatistics_ , 12(1):3–19, 2016.  ISSN 15574679.  doi:10.1515/ijb-2015-0013.
  * Iván Díaz Muñoz and Mark van der Laan* [2012] Iván Díaz Muñoz and Mark van der Laan*.  Population Intervention Causal Effects Based on Stochastic Interventions.  _Biometrics._ , 68(2):541–549, 2012.  ISSN 15378276.  doi:10.1111/j.1541-0420.2011.01685.x.Population.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3624763/pdf/nihms412728.pdf.
  * Keil et al. [2019] Alexander P. Keil, Jessie P. Buckley, Katie M. O’Brien, Kelly K. Ferguson, Shanshan Zhao, and Alexandra J. White.  A quantile-based g-computation approach to addressing the effects of exposure mixtures.  _arXiv_ , 128(April):1–10, 2019.  ISSN 23318422.  doi:10.1097/01.ee9.0000606120.58494.9d.
  * Lengerich et al. [2020] Benjamin Lengerich, Sarah Tan, Chun Hao Chang, Giles Hooker, and Rich Caruana.  Purifying Interaction Effects with the Functional ANOVA: An Efficient Algorithm for Recovering Identifiable Additive Models.  _Proceedings of Machine Learning Research_ , 108:2402–2412, 2020.  ISSN 26403498.
  * McCoy et al. [2023] David McCoy, Alejandro Schuler, Alan Hubbard, and Mark van der Laan.  Supernova: Semi-parametric identification and estimation of interaction and effect modification in mixed exposures using stochastic interventions in r.  _Journal of Open Source Software_ , 8(91):5422, 2023.  doi:10.21105/joss.05422.  URL https://doi.org/10.21105/joss.05422.
  * McEligot et al. [2020] Archana J. McEligot, Valerie Poynor, Rishabh Sharma, and Anand Panangadan.  Logistic lasso regression for dietary intakes and breast cancer.  _Nutrients_ , 12(9):1–14, 2020.  ISSN 20726643.  doi:10.3390/nu12092652.
  * Milborrow. Derived from mda:mars by T. Hastie and R. Tibshirani. [2011] S. Milborrow. Derived from mda:mars by T. Hastie and R. Tibshirani.  _earth: Multivariate Adaptive Regression Splines_ , 2011.  URL http://CRAN.R-project.org/package=earth.  R package.
  * Mitro et al. [2016] Susanna D Mitro, Linda S Birnbaum, Belinda L Needham, and Ami R Zota.  Cross-sectional associations between exposure to persistent organic pollutants and leukocyte telomere length among us adults in nhanes, 2001–2002.  _Environmental health perspectives_ , 124(5):651–658, 2016.
  * Narisetty et al. [2019] Naveen N. Narisetty, Bhramar Mukherjee, Yi-Hau Chen, Richard Gonzalez, and John D. Meeker.  Selection of nonlinear interactions by a forward stepwise algorithm: Application to identifying environmental chemical mixtures affecting health outcomes.  _Statistical Medicine_ , 38(9):1582–1600, Apr 2019.  doi:10.1002/sim.8059.
  * Ripley and Venables [2021] B. D. Ripley and W. Venables.  _polspline: Polynomial Spline Routines_ , 2021.  URL https://CRAN.R-project.org/package=polspline.  R package version 1.1.26.
  * Roberts and Martin [2006] Steven Roberts and Michael A. Martin.  Using supervised principal components analysis to assess multiple pollutant effects.  _Environmental Health Perspectives_ , 114(12):1877–1882, 2006.  ISSN 00916765.  doi:10.1289/ehp.9226.
  * Rothman et al. [1980] Kenneth J Rothman, Sander Greenland, and Alexander M Walker.  Concepts of interaction.  _Am J Epidemiol_ , 112(4):467–470, 1980.
  * Rubin [1974] Donald B Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _J. Educ. Psychol._ , 66(5):688–701, October 1974.
  * Rubin [2005] Donald B Rubin.  Causal inference using potential outcomes.  _J. Am. Stat. Assoc._ , 100(469):322–331, March 2005.
  * Tay et al. [2023] J. Kenneth Tay, Balasubramanian Narasimhan, and Trevor Hastie.  Elastic net regularization paths for all generalized linear models.  _Journal of Statistical Software_ , 106(1):1–31, 2023.  doi:10.18637/jss.v106.i01.
  * van der Laan [2006] Mark J van der Laan.  Statistical inference for variable importance.  _Int. J. Biostat._ , 2(1), February 2006.
  * van der Laan Mark et al. [2007] J. van der Laan Mark, Polley Eric C, and Hubbard Alan E.  Super learner.  _Statistical Applications in Genetics and Molecular Biology_ , 6(1):1–23, 2007.  URL https://EconPapers.repec.org/RePEc:bpj:sagmbi:v:6:y:2007:i:1:n:25.
  * VanderWeele [2009] Tyler J VanderWeele.  On the distinction between interaction and effect modification.  _Epidemiology_ , 20(6):863–871, 2009.
  * VanderWeele and Knol [2014] Tyler J VanderWeele and Mirjam J Knol.  A tutorial on interaction.  _Epidemiologic Methods_ , 3(1):33–72, 2014.
  * VanderWeele and Robins [2007] Tyler J VanderWeele and James M Robins.  Directed acyclic graphs, sufficient causes, and the properties of conditioning on a common effect.  _Am J Epidemiol_ , 166(10):1096–1104, 2007.
  * Wright and Ziegler [2017] Marvin N. Wright and Andreas Ziegler.  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  _Journal of Statistical Software_ , 77(1):1–17, 2017.  doi:10.18637/jss.v077.i01.
  * Zheng and van der Laan [2010] Wenjing Zheng and MJ van der Laan.  Asymptotic theory for cross-validated targeted maximum likelihood estimation.  _U.C. Berkeley Division of Biostatistics Working Paper Series_ , (273), 2010.  URL http://biostats.bepress.com/ucbbiostat/paper273/.